\section{Introduction} \label{ch11:Introduction}
-Monotonicity preserving interpolation and approximation have received substantial attention in the last thirty years because of their numerous applications in computer aided design, statistics and machine learning \cite{Dierckx1995_book,Kvasov2000_book,deboor2001_book}. Constrained splines are particularly popular because of their flexibility in modelling different geometrical shapes, sound theoretical properties and availability of numerically stable algorithms \cite{Dierckx1995_book,Schumaker1981_book,deboor2001_book}.
+Monotonicity preserving interpolation and approximation have received substantial attention in the last thirty years because of their numerous applications in computer aided design, statistics and machine learning \cite{Dierckx1995_book,Kvasov2000_book,deboor2001_book}. Constrained splines \index{spline} \index{constrained splines} \index{monotonicity} are particularly popular because of their flexibility in modelling different geometrical shapes, sound theoretical properties and availability of numerically stable algorithms \cite{Dierckx1995_book,Schumaker1981_book,deboor2001_book}.
% It is surprising though that few parallel spline algorithms are available.
In this work we examine parallelisation and adaptation for GPUs of a few algorithms of monotone spline interpolation and data smoothing, which arose in the context of estimating probability distributions.
Estimating cumulative probability distribution functions (cdf) from data is quite common in data analysis. In our particular case we faced this problem in the context of partitioning univariate data with the purpose of efficient sorting. It was required to partition large data sets into chunks of approximately equal size, so that these chunks could be sorted independently and subsequently concatenated. In order to do that, empirical cdf of the data was used to find the quantiles, which served to partition the data. Cdf was estimated from the data based on a number of pairs $(x_i,y_i), i=1,\ldots,n$, where $y_i$ was the proportion of data no larger than $x_i$. As data could come from a variety of distributions, a distribution-free nonparametric fitting procedure was required to interpolate the above pairs. Needless to say that the whole process was aimed at GPU, and hence the use of CPU for invoking serial algorithms had to be minimised.
-The above mentioned application is one of many examples (e.g. mass spectrography \cite{Kearsley_2006}, global warming data \cite{Yohai} and so on) where univariate data needs to be fitted by monotonicity preserving interpolants. Of course, cdf is a monotone increasing function, whose inverse, called quantile function, can be used to calculate the quantiles. Spline interpolation would be the most suitable nonparametric method to fit the cdf, except that polynomial splines do not preserve monotonicity of the data, as illustrated on Figure \ref{ch11:fig1}.
+The above mentioned application is one of many examples (e.g. mass spectrography \cite{Kearsley_2006}, global warming data \cite{Yohai} and so on) where univariate data needs to be fitted by monotonicity preserving interpolants. Of course, cdf is a monotone increasing function, whose inverse, called quantile function, can be used to calculate the quantiles. Spline interpolation would be the most suitable nonparametric method to fit the cdf, except that polynomial splines do not preserve monotonicity of the data, as illustrated on Figure \ref{ch11:fig1}.
The failure of splines to preserve monotonicity has prompted fundamental research in this area since 1960s. One of the first methods to remedy this problem were splines in tension by Schweikert \cite{Sch}, where a tension parameter controlled the shape of exponential splines \cite{Spath1969}. Later on several monotonicity preserving polynomial spline algorithms were proposed \cite{Schumaker1983,PasRoul1977,AndElf1987,Andersson1991_JAT,McAllister1981_ACM,PasRoul1977}. These algorithms typically rely on introducing additional spline knots between the abscissae of the data. Algorithmic developments are active to this day, see for example \cite{Kvasov2000_book,Abbas2011}.
-When in addition to the pairs $(x_i, y_i)$ the slopes of the function are available, i.e., the data comes in triples $(x_i, y_i, p_i)$, the interpolation problem is called Hermite, and the Hermite splines are used. However, even when the sequence $y_i$ is increasing and the slopes $p_i$ are non-negative, cubic Hermite splines may still fail to be monotone, as illustrated in Figure \ref{ch11:fig2}. Thus monotone Hermite splines are needed \cite{Gregory1982}.
+When in addition to the pairs $(x_i, y_i)$ the slopes of the function are available, i.e., the data comes in triples $(x_i, y_i, p_i)$, the interpolation problem is called Hermite, and the Hermite splines are used. However, even when the sequence $y_i$ is increasing and the slopes $p_i$ are non-negative, cubic Hermite splines may still fail to be monotone, as illustrated in Figure \ref{ch11:fig2}. Thus monotone Hermite splines are needed \cite{Gregory1982}. \index{Hermite splines}
Another issue with monotone approximation is noisy data. In this case, inaccuracies in the data make the input sequence $y_i$ itself non-monotone, and hence monotone spline interpolation algorithms will fail. Monotone spline smoothing algorithms are available, e.g. \cite{Andersson1991_JAT,Elfving1989_NM}. Such algorithms are based on solving a quadratic (or another convex) programming problem numerically, and have not been yet adapted to parallel processing.
In this work we examined several monotone spline fitting algorithms, and selected the ones that we believe are most suitable for parallelisation on GPUs. We paid attention to numerical efficiency in terms of numerical calculations and memory access pattern, and favoured one-pass algorithms. We also looked at smoothing noisy data, and developed a parallel version of the Minimum Lower Sets algorithm for isotonic regression problem \cite{Best1990, Robertson_book}.
+\index{isotone regression}
The rest of the chapter is organised as follows. Section \ref{ch11:splines} discusses monotone spline interpolation methods and presents two parallel algorithms. Section \ref{ch11:smoothing} deals with smoothing problem. It presents isotonic regression problem and discusses the Pool Adjacent Violators (PAV) and Minimum Lower Sets (MLS) algorithms. Combined with monotone spline interpolation, the parallel MLS method makes it possible to build a monotone spline approximation to noisy data entirely on GPU. Section \ref{ch11:conc} concludes.
\begin{figure}[h]
\centering
\includegraphics[angle=0,width=8cm]{Chapters/chapter11/gregory1_plot1.pdf}
-\caption{Cubic spline (solid) and monotone quadratic spline (dashed) interpolating monotone data from \cite{Gregory1982}. Cubic spline fails to preserve monotonicity of the data.}
+\caption[Cubic spline (solid) and monotone quadratic spline (dashed) interpolating monotone data]{Cubic spline (solid) and monotone quadratic spline (dashed) interpolating monotone data from \cite{Gregory1982}. Cubic spline fails to preserve monotonicity of the data.}
\label{ch11:fig1}
\end{figure}
\begin{figure}[h]
\centering
\includegraphics[angle=00,width=8cm]{Chapters/chapter11/gregory1_plot2_b.pdf}
-\caption{Hermite cubic spline (solid) and Hermite rational spline interpolating monotone data from \cite{Gregory1982} with non-negative prescribed slopes. Despite non-negative slopes, Hermite cubic spline is not monotone.}
+\caption[Hermite cubic spline (solid) and Hermite rational spline interpolating monotone data]{Hermite cubic spline (solid) and Hermite rational spline interpolating monotone data from \cite{Gregory1982} with non-negative prescribed slopes. Despite non-negative slopes, Hermite cubic spline is not monotone.}
\label{ch11:fig2}
\end{figure}
\section{Monotone splines} \label{ch11:splines}
+\index{constrained splines} \index{monotonicity}
Splines are piecewise continuous functions very popular in numerical approximation and computer aided design \cite{deboor2001_book,Dierckx1995_book}. An example of a spline is broken line interpolation. Typically polynomial splines are used, and the first (and often second) derivatives of the polynomial pieces are required to match at the knots. The knots of the splines are usually the abscissae of the input data, although this condition is not always required (e.g., splines with free knots \cite{Jupp_1978,Dierckx1995_book,Beliakov2003_amc}).
Polynomial splines are often represented in the B-spline basis, in which case their coefficients are computed from the input data by solving a banded system of linear equations \cite{Lyche1973, Dierckx1995_book, deboor2001_book}. Tridiagonal systems arise in cubic spline interpolation, while pentadiagonal systems arise in cubic spline smoothing \cite{Lyche1973}. Spline possess important extremal properties \cite{Holladay1957,Lyche1973}, in particular splines of degree $2m-1$ are the most ``smooth" functions that interpolate (or approximate, in the least squares sense) the data. The smoothness term is Tihkonov regularisation functional, the $L_2$ norm of the $m$-th derivative of the interpolant \cite{Lyche1973}.
2\delta_{1}-d_2, & \mbox{if } \delta_{1}(2\delta_1-d_2)>0, \\
0 & \mbox{otherwise},
\end{array}
- \right. \;
+ \right.
+ $$
+ $$
d_n=\left\{\begin{array}{ll}
2\delta_{n-1}-d_{n-1}, & \mbox{if } \delta_{n-1}(2\delta_{n-1}-d_{n-1})>0,\\
0 & \mbox{otherwise}.
It is almost straightforward to parallelise this scheme for GPUs, by processing each subinterval $[x_i,x_{i+1}]$ independently in a separate thread. However, it is not known in advance whether an extra knot $t_i$ needs to be inserted or not, and therefore calculation of the position of the knot in the output sequence of knots ${t_i}$ is problematic for parallel implementation (for a sequential algorithm no such issue arises). To avoid serialisation, we decided to insert an additional knot in every interval $[x_i,x_{i+1}]$, but set $t_i=x_i$ when the extra knot is not actually needed. This way we know in advance the position of the output knots and the length of this sequence is $2(n-1)$, and therefore all calculations can now be performed independently. The price we pay is that some of the spline knots can coincide. However, this does not affect spline evaluation, as one of the coinciding knots is simply disregarded, and the spline coefficients are replicated (so for a double knot $t_i=t_{i+1}$, we have $\alpha_i=\alpha_{i+1}$, $\beta_i=\beta_{i+1}$, $\gamma_i=\gamma_{i+1}$). Our implementation is presented in Figures \ref{ch11:algcoef}-\ref{ch11:algcoef1}.
+\lstinputlisting[label=ch11:algcoef1,caption=Calculation of monotone spline knots and coefficients.]{Chapters/chapter11/code2.cu}
+
+
At the spline evaluation stage we need to compute $s(z_k)$ for a sequence of query values ${z_k}, k=1,\ldots,K$. For each $z_k$ we locate the interval $[t_i,t_{i+1}]$ containing $z_k$, using bisection algorithm presented in Figure \ref{ch11:algeval}, and then apply the appropriate coefficients of the quadratic function. This is also done in parallel.
The bisection algorithm could be implemented using texture memory (to cache the array \texttt{z}), but this is not shown in Figure \ref{ch11:algeval}.
-\lstinputlisting[label=ch11:algcoef,caption=Implementation of the kernel for calcuating spline knots and coefficients. Function fmax is used to avoid division by zero for data with coinciding abscissae.]{Chapters/chapter11/code1.cu}
+\pagebreak
+\lstinputlisting[label=ch11:algcoef,caption=Implementation of the kernel for calculating spline knots and coefficients. Function fmax is used to avoid division by zero for data with coinciding abscissae.]{Chapters/chapter11/code1.cu}
%% \begin{figure}[!hp]
%% \end{figure}
-\lstinputlisting[label=ch11:algcoef1,caption=Calculation of monotone spline knots and coefficients.]{Chapters/chapter11/code2.cu}
%% \begin{figure}[!hp]
%% \renewcommand{\baselinestretch}{1}
\lstinputlisting[label=ch11:algeval,caption=Implementation of the spline evaluation algorithm for GPU.]{Chapters/chapter11/code3.cu}
\subsection{Monotone Hermite splines}
-
+\index{Hermite splines} \index{monotonicity}
In this section, in addition to the points $(x_i,y_i)$ we have the slopes $p_i$, and hence we consider monotone Hermite interpolation. In our motivating application of cdf estimation, the values $p_i$ are easily obtained together with $y_i$, and their use may help to build a more accurate interpolant.
Of course, for monotone non-decreasing functions we must have $p_i\geq 0$. However this does not guarantee that the spline interpolant is monotone, as can be seen in Figure \ref{ch11:fig2}. Fritsch and Carlson \cite{Fritsch1980} show that non-negative $p_i$ is not a sufficient condition to guarantee monotonicity, and design a process for modification of derivatives, so that the necessary and sufficient conditions for monotonicity of a piecewise cubic are met. Hence the values $p_i$ are not matched exactly. In contrast, Gregory and Delbourgo \cite{Gregory1982} design piecewise rational quadratic spline, for which the non-negativity of $p_i$ is both necessary and sufficient condition.
\section{Smoothing noisy data via parallel isotone regression} \label{ch11:smoothing}
+
Inaccuracies in the data are common in practice, and need to be accounted for during spline approximation process. Smoothing polynomial splines were presented in \cite{Lyche1973}, where the data are fitted in the least squares sense while also minimising the $L_2$ norm of the $m-$th derivative of the spline. Monotone smoothing splines were dealt with in several works, in particular we mention \cite{Andersson1991_JAT,Elfving1989_NM}. The presented algorithms rely on solving quadratic programming problems. Monotone approximating splines with fixed knots distinct form the data have been presented in \cite{Beliakov2000_ata}, where an instance of a quadrating programming problem is solved as well.
+\index{isotone regression} \index{monotonicity}
Another approach consists in monotonising the data, so that the sequence $y_i$ becomes monotone. This approach is known as isotone regression \cite{Best1990, Robertson_book}. It is different from monotone spline smoothing, as the regularisation term controlling the $L_2$ norm of the $m-$the derivative is not taken into account. Usually the data is monotonised by minimising the squared differences to the inputs. It becomes a quadratic programming problem, usually solved by active sets methods \cite{Best1990}.
A popular PAV algorithm (PAVA) is one method that provides efficient numerical solution.
+\index{PAV algorithm}
PAVA consists of the following steps. The sequence ${y_i}$ is processed form the start. If violation of monotonicity $y_i>y_{i+1}$ is found, both values $y_i$ and $y_{i+1}$ are replaced with their average $y'_i$, and both values form a block. Since the new value $y'_i$ is smaller than $y_i$, monotonicity may become violated with respect to the datum $y_{i-1}$. If this is the case, the $i-1$st, $i$th and $i+1$st data are merged into a block and their values are replaced with their average. We continue back-average as needed to get monotonicity.
Various serial implementations of the PAVA exist. It is noted \cite{Kearsley_2006} that in PAVA, which is based on the ideas from convex analysis, a decomposition theorem holds, namely performing PAVA separately on two contiguous subsets of data, and then performing PAVA on the result produces isotonic regression on the whole data set. Thus isotonic regression is parallelisable, and divide-and-conquer approach, decomposing the original problem into two smaller subproblems, can be implemented on multiple processors. However, to our knowledge, no parallel PAVA for many-core systems such as GPUs exist.
-Another approach to isotonic regression is called the Minimum Lower Sets algorithm (MLS) \cite{Best1990, Robertson_book}. It provides the same solution as the PAVA, but works differently. For each datum (or block), MLS selects the largest contiguous block of subsequent data with the smallest average. If this average is smaller than that of the preceding block, the blocks are merged, and the data in the block are replaced with their average. MLS is also an active set method \cite{Best1990}, but its complexity is $O(n^2)$ as opposed to $O(n)$ of the PAVA, and of another active set algorithm proposed in \cite{Best1990} under the name Algorithm A.
+\index{MLS algorithm} \index{Minimum Lower Sets}
+Another approach to isotonic regression is called the Minimum Lower Sets algorithm (MLS)
+\cite{Best1990, Robertson_book}. It provides the same solution as the PAVA, but works differently. For each datum (or block), MLS selects the largest contiguous block of subsequent data with the smallest average. If this average is smaller than that of the preceding block, the blocks are merged, and the data in the block are replaced with their average. MLS is also an active set method \cite{Best1990}, but its complexity is $O(n^2)$ as opposed to $O(n)$ of the PAVA, and of another active set algorithm proposed in \cite{Best1990} under the name Algorithm A.
In terms of GPU parallelisation, neither PAVA nor Algorithm A appear to be suitable, as the techniques that achieve $O(n)$ complexity are inheritably serial.
In this work we focused on parallelising MLS. First, we precompute the values
From the results in Table \ref{ch11:table1} we conclude that serial PAVA is superior to MLS for $n>10^4$. While it is possible to transfer data from GPU to CPU and run PAVA there, it is warranted only for sufficiently large data $n\geq 5 \times 10^5$ , for otherwise the data transfer overheads will dominate CPU time. For smaller $n$, isotone regression is best performed on GPU.
-We also see that the use of GPU accelerated MLS by a factor of at least 100. The cost of serial MLS is prohibitive for $n>10^6$.
+We also see that the use of GPU accelerated MLS by a factor of at least 100. The cost of serial MLS is prohibitive for $n>10^6$.
We should mention that not all isotone regression problems allow a PAV-like algorithm linear in time. When the data may contain large outliers, monotonizing the data is better done not in the least squares sense, but using other cost functionals, such as by minimizing the sum of absolute deviations \cite{Wang} or using M-estimators \cite{Yohai}, which are less sensitive to outliers. It is interesting than in all such cases the solution to isotone regression problem can be found by solving maximin problem
$$
-u_i=\max_{k\leq i} \min_{l \geq i} \hat y(k,l),
+u_i=\max_{k\leq i} \min_{l \geq i} \hat y(k,l),
$$
with $\hat y(k,l)$ being the unrestricted maximum likelihood estimator of $y_k\ldots,y_l$. For quadratic cost function $\hat y(k,l)$ is the mean, as in PAV and MLS algorithms, for the absolute deviations it becomes the median, and for other cost functions an M-estimator of location. The MLS algorithm can be applied to such isotone regression problems with very little modification, while linear in time algorithm may not be available. Our parallel MLS algorithm will be valuable in such cases.
%% %\renewcommand{\baselinestretch}{1}
-%% \begin{table}[!h]
-%% \begin{center}
-%% \caption{The average CPU time (sec) of the serial PAVA, MLS and parallel MLS algorithms. } \label{ch11:table1}
-%% \begin{tabular}{|r|r|r|r|}
-
-%% Data & PAVA & MLS & GPU MLS \\ \hline
-
-%% monotone increasing $f$ & & & \\
-%% $n=5\times 10^4$ &0.01&5& 0.092\\
-%% $n=10^5$ &0.03&40& 0.35\\
-%% $n=5\times 10^5$ &0.4&1001&8.6 \\
-%% $n=10^6$ &0.8& 5000& 38 \\
-%% $n=2 \times 10^6$ & 1.6 &-- &152 \\
-%% $n=10 \times 10^6$ & 2 &-- & 3500 \\
-%% $n=20 \times 10^6$ & 4.5&-- & --\\
-%% $n=50 \times 10^6$ & 12 &-- & --\\
-%% \hline
-
-%% constant or decreasing $f$ & & & \\
-%% $n=10^6$ &0.2&0.1& 38\\
-%% $n=10 \times 10^6$ &1.9& 1.9& 3500 \\
-%% $n=20 \times 10^6$ &3.5& 4.0&-- \\
-%% $n=50 \times 10^6$ &11& 11& -- \\
-
-%% \end{tabular}
-%% \end{center}
-%% \end{table}
+\begin{table}[!h]
+\begin{center}
+\begin{tabular}{|r|r|r|r|}
+\hline
+Data & PAVA & MLS & GPU MLS \\ \hline
+
+monotone increasing $f$ & & & \\
+$n=5\times 10^4$ &0.01&5& 0.092\\
+$n=10^5$ &0.03&40& 0.35\\
+$n=5\times 10^5$ &0.4&1001&8.6 \\
+$n=10^6$ &0.8& 5000& 38 \\
+$n=2 \times 10^6$ & 1.6 &-- &152 \\
+$n=10 \times 10^6$ & 2 &-- & 3500 \\
+$n=20 \times 10^6$ & 4.5&-- & --\\
+$n=50 \times 10^6$ & 12 &-- & --\\
+\hline
+
+constant or decreasing $f$ & & & \\
+$n=10^6$ &0.2&0.1& 38\\
+$n=10 \times 10^6$ &1.9& 1.9& 3500 \\
+$n=20 \times 10^6$ &3.5& 4.0&-- \\
+$n=50 \times 10^6$ &11& 11& -- \\
+\hline
+\end{tabular}
+\end{center}
+\caption{The average CPU time (sec) of the serial PAVA, MLS and parallel MLS algorithms. }
+\label{ch11:table1}
+\end{table}
%% %\renewcommand{\baselinestretch}{2}
%% \begin{alltt}
%% \begin{center}
%% \begin{minipage}{13cm}\small
-%% template<typename Tx>
+%% template<typename Tx>
%% __device__ Tx Aver(Tx z,int i,int j, Tx *z) \{return (z-z[j+1])/(j-i+1);\}
%% template<typename Tx>
-%% __global__ void monotonizekernel(Tx *y, Tx *z, Tx *u, int *key, int n)
+%% __global__ void monotonizekernel(Tx *y, Tx *z, Tx *u, int *key, int n)
%% \{ int i = threadIdx.x + blockIdx.x * blockDim.x;
%% if(i<n) \{
%% int smallestJ = i;
%% thrust::inclusive_scan(y_reverse_b, y_reverse_end, z_reverse_b+1);
-%% monotonizekernel<<<grid, block>>>(y, thrust::raw_pointer_cast(&z_d[0]),
+%% monotonizekernel<<<grid, block>>>(y, thrust::raw_pointer_cast(&z_d[0]),
%% u, thrust::raw_pointer_cast(&keys_d[0]), n );
%% thrust::sort(keys_d.begin(), keys_d.end());
-%% thrust::inclusive_scan_by_key(keys_d.begin(), keys_d.end(),
+%% thrust::inclusive_scan_by_key(keys_d.begin(), keys_d.end(),
%% u_d, u_d, binary_pred, binary_op2);
%% \}
%% \end{minipage}
\section{Conclusion} \label{ch11:conc}
We presented three GPU-based parallel algorithms for approximating monotone data: monotone quadratic spline, monotone Hermite rational spline and minimum lower sets algorithm for monotonizing noisy data. These tools are valuable in a number of applications that involve large data sets modeled by monotone nonlinear functions.
-The source code of the package monospline is available from \texttt{www.deakin.edu.au/$\sim$ gleb/monospline.html }
+The source code of the package monospline is available from \texttt{www.deakin.edu.au/$\sim$gleb/monospline.html }
+
+
\putbib[Chapters/chapter11/biblio11]
