Univariate Cubic L 1 Interpolating Splines: Spline Functional, Window Size and Analysis-based Algorithm

: We compare univariate L 1 interpolating splines calculated on 5-point windows, on 7-point windows and on global data sets using four different spline functionals, namely, ones based on the second derivative, the ﬁrst derivative, the function value and the antiderivative. Computational results indicate that second-derivative-based 5-point-window L 1 splines preserve shape as well as or better than the other types of L 1 splines. To calculate second-derivative-based 5-point-window L 1 splines, we introduce an analysis-based, parallelizable algorithm. This algorithm is orders of magnitude faster than the previously widely used primal afﬁne algorithm.

In the literature, there are a few indications of limitations of the primal affine and primal-dual algorithms for bivariate L 1 splines for large data sets [9,12].There is also unpublished computational experience of the authors and others who have noticed issues of incomplete convergence or not completely correct convergence of the active set, primal affine and primal-dual algorithms for both univariate and bivariate L 1 splines.It is in this context that we seek in this paper a new algorithmic approach for calculating L 1 splines.Auquiert, Gibaru and Nyiri [16] have developed a subdifferential-based procedure for calculating second-derivative-based L 1 splines on 5-point windows.We propose here an algorithm for second-derivative-based 5-point-window L 1 splines based on the analysis in Section 2 of [17], which links, via analytical properties of the spline functional, local geometric properties of 5-point windows of the data set with geometric properties of the L 1 spline interpolant.
In considering second-derivative-based 5-point-window L 1 splines, three types of information are needed for full "situational awareness," namely, 1) information about whether use of 5-point windows is superior to use of windows of other sizes and to use of global data sets, 2) information about whether use of the second derivative in the spline functional is superior to use of the first derivative, function value or antiderivative and 3) information about whether the new, analysis-based algorithm mentioned above can achieve computational results superior to those of the primal affine algorithm, which was the previously most widely used algorithm for calculating L 1 splines.The third item here has not yet been considered in the literature and the information in the literature on the first and second items is sketchy at best.The computational and analytical results in [2,16,17] do suggest that 5-point windows have advantages vs. global calculations.However, it is not yet known whether windows of other sizes might also have advantages.Results presented in [10,11] indicate that first-derivative-based and function-value-based L 1 splines may have advantages over standard second-derivative-based L 1 splines.However, these two publications considered only the behavior of L 1 splines on the global scale and did not consider behavior on the fine, interval-to-interval scale.
The present paper addresses these needs.In Section 2, we give a brief description of the primal affine algorithm that has previously been widely used to calculate L 1 splines.In Section 3, we compare L 1 splines calculated by minimizing, on 5-point windows, on 7-point windows and on global data sets, four different spline functionals, namely, ones based on the second derivative, the first derivative, the function value and the antiderivative.These L 1 splines are calculated by the primal affine algorithm.The results of this section provide motivation for the development of a new, analysis-based algorithm for calculating 5-point-window, second-derivative-based L 1 splines.In Section 4, we present this new algorithm, which is based on the analysis in Section 2 of [17].In Section 5, we show that, while the results of both the new algorithm and the primal affine algorithm look good on the macro level, there are differences on the micro level.Specifically, the results of the new algorithm are accurate on the micro level while those of the primal affine algorithm are occasionally only approximate.Finally, in Section 6, we summarize the results presented in the previous sections and point out the potential for future algorithms for locally calculated univariate L 1 approximating splines and locally calculated bivariate L 1 interpolating and approximating splines.
All of the quantities in this paper are real quantities.The nodes x i , i = 0, 1, . . ., I, are a strictly monotonic but otherwise arbitrary partition of the finite interval [x 0 , x I ].Let h i = x i+1 − x i .At each node x i , the function value z i is given, i = 0, 1, . . ., I. The slope of the line segment connecting (x i , z i ) and (x i+1 , z i+1 ) is The local and global cubic L 1 splines discussed in this paper are cubic polynomials in each interval (x i , x i+1 ), i = 0, 1, . . ., I − 1, and are C 1 continuous at the nodes.The first derivative of the spline at node x i , i = 0, 1, . . ., I, is denoted by b i (to be determined by minimization of the L 1 spline functional).We use δ i to denote the slope of the chord between neighboring points: ( We use ζ to denote the linear spline: For the interpolation problem under consideration in the present paper, the nodes and the function values at the nodes are given.Minimization of a spline functional means "determination of the first derivatives that yield the minimum."The first derivatives determined from the minimization along with the fixed nodal and functional values yield the piecewise cubic L 1 spline by the standard Hermite interpolation formula.For reference, we present here the spline functionals that define traditional, globally calculated L 1 splines ("global L 1 splines").Second-derivative-based, first-derivative-based, function-value-based and antiderivative-based cubic L 1 splines are calculated by minimizing and respectively, over the finite-dimensional spline space of C 1 piecewise cubic polynomials z that interpolate the data.The L 1 splines that minimize these functionals can be nonunique, a situation that will be handled by adding "regularization terms" to the functionals when they are minimized by the primal affine algorithm or, in the new algorithm proposed in this paper, by applying a "choice procedure" as described below in Section 4. Second-derivative-based L 1 splines are the L 1 splines commonly encountered in the literature.First-derivative-based L 1 splines have been investigated in [10,11].Function-value-based L 1 splines have been treated in [11].Antiderivative-based L 1 splines are newly introduced in this present paper to provide additional insight into how the order of the derivative in the spline functional affects the geometric shape preservation properties of the L 1 spline.
In the present paper, computations of different algorithms are based on the following challenging data set consisting of 56 irregular data points that lie on flat, linear, quadratic, cubic and oscillatory functions and on protuberances, patched together with discontinuities of function values and first derivatives and with extreme irregular spacing-with lengths of neighboring intervals differing by up to a factor of 100: This data set was used in [11].The data (x i , z i ), i = 27, 28, 29, 30, 31, 32, 33, 34, in the region from x = 27 to x = 35 lie on the the quadratic function 44 − 2.75(x − 31) 2 .The data (x i , z i ), i = 36, 37, 38, 39, 40, 41, 42, 43, in the region from x = 37 to x = 45 lie on the cubic function −16(x − 41) + (x − 41) 3 .Note the large gaps in the intervals [x 30 , x 31 ] = [27.3, 34.7] and [x 39 , x 40 ] = [37.3, 44.7].In the figures, the data will be represented by dots "•".The linear spline for these data is given in Figure 1.

Primal Affine and Other Previously Available Algorithms for L 1 Splines
Minimization of the global L 1 spline functionals (4), ( 5), ( 6) and ( 7) is a nonlinear programming problem.Direct minimization of (4) has been accomplished by an active set method [3].Active set algorithms for minimizing (5), ( 6) and ( 7) have not been developed, but these functionals and functional (4) have been minimized by primal affine algorithms and primal-dual algorithms [12,14].Primal affine and primal-dual algorithms are linear (not nonlinear) programming procedures.To create a linear program suitable for application of these algorithms, the integrals in the L 1 spline functionals need to be discretized.For the primal affine algorithm used in the present paper and in [5,6,8,10,11], the spline functionals were discretized by the midpoint rule with K equal subintervals in each interval (x i , x i+1 ).In this paper, K = 100.For a detailed description of the primal affine algorithm, see [7].For calculation of L 1 splines by the primal affine method using the global functionals ( 4), ( 5), ( 6) and ( 7), we add to these functionals the regularization terms When ε is sufficiently small, the L 1 spline that minimizes the functional with the regularization term is unique.For the computational experiments of the present paper, ε = 10 −4 .
In addition to calculating L 1 splines by minimizing the global functionals (4), ( 5), ( 6) and ( 7) by the primal affine method, we will calculate L 1 splines using these functionals on 5-point and 7-point windows.To calculate the derivative b i of a "windowed" L 1 spline at node x i , we minimize functionals that are the same as ( 4), ( 5), ( 6) and (7) except that the integral is over a local set of intervals ("window") rather than over the global domain.For 5-point-window L 1 splines, the window [x î−2 , x î+2 ] is used for nodes x î, î = 2, 3, . . ., I − 2, the window [x 0 , x 4 ] is used for nodes x î, î = 0, 1 and the window [x I−4 , x I ] is used for nodes x î, î = I − 1, I. For 7-point-window L 1 splines, the window [x î−3 , x î+3 ] is used for nodes x î, î = 3, 4, . . ., I − 3, the window [x 0 , x 6 ] is used for nodes x î, î = 0, 1, 2 and the window [x I−6 , x I ] is used for nodes x î, î = I − 2, I − 1, I. When minimizing a "windowed" L 1 spline functional for node x î by the primal affine algorithm, we add to the spline functional the regularization term For the computational experiments of the present paper, ε = 10 −4 .

Computational Results for Windowed and Global L 1 Splines
A window with an odd number of points is desirable, since it has a "middle point" at which one can use a locally calculated derivative as the derivative of the global interpolant.Conversely, a window with an even number of points is not desirable, because it lacks a middle point.Using the primal affine algorithm described in the previous section, we generated computational results for 5-point-window and 7-point-window L 1 splines as well as for global L 1 splines.
In the literature, there have been reports [10,11] about L 1 splines based on spline functionals involving first derivatives and function values, that is, derivatives of degree lower than the standard second degree.Results in [10,11] suggest that, for global L 1 splines, spline functionals based on the first derivative or function value (that is, on functional ( 5) or ( 6)) result in improved shape preservation.Those results emphasized overall global preservation of shape on the macro scale but did not treat local preservation of shape on the fine, interval-to-interval scale.Following up on these prior investigations, we revisit here the comparison of L 1 splines based on derivatives of degree 2, 1 and 0, add to the comparison L 1 splines based on derivatives of degree −1 (antiderivatives) and add the new dimension of considering the window size (5 points, 7 points or global).For this comparison, we use the primal affine algorithm described in Section 2.
Computational results are presented in  3, 34.7], even though they do approximate the cubic function in [37.3, 44.7].In classical approximation theory, the ability to approximate or even reproduce quadratic, cubic and higher-degree functions is a goal.In geometric modeling of irregular data, however, reproducing any function of degree higher than 1 is generally considered disadvantageous because it leads to extraneous oscillation, overshoot or undershoot.It is not yet known how to construct L 1 splines that generically avoid approximating functions of degree higher than 1.In the present article, we acknowledge this issue but we do not take it into account when assessing the shape-preservation capabilities of the various types of L 1 splines.We do not claim here that the results presented in Figures 2-13 should yet be broadly interpreted as general descriptions of the shape-preservation capabilities of the various types of L 1 splines.However, these results do provide insight into the shape-preservation capabilities of these types of L 1 splines as well as a piece of evidence supporting the use of second-derivative-based 5-point-window L 1 splines.In these computational results, the performance of second-derivative-based L 1 splines, the hitherto most widely used variant of L 1 splines is improved by using 5-point windows.The other two "top performers" of Figures 2-13 Elementary theory has been published for first-derivative-based L 1 splines but only for the global case [15].Based on these considerations, we choose the second-derivative-based 5-point-window L 1 spline as the L 1 spline to be used in the remainder of this paper.

Algorithm for Minimization of Second-derivative-based 5-point-window Spline Functional
We propose here an algorithm for generation of second-derivative-based 5-point-window L 1 splines based on the analytical results of Section 2 of [17].Recall that, in the interpolation situation of interest in this paper, the task of calculating an L 1 spline is the task of calculating the first derivatives b i by minimization of a spline functional with given nodes x i and given function values z i at the nodes.In the ith 5-point window, In [17], minimization of where and G 1 (b i ) and G 2 (b i ) have similar structure, which leads to the representation where G(q; c) is a function that is defined and analyzed in [17].
The analysis of [17] is based on the signs of △z i−2 − △z i−1 , △z i−1 − △z i and △z i − △z i+1 .There are 27 cases, as shown in Table 1.For the simple cases, the optimal b * i are listed in the last column of the table.For cases 14,15,17,18,23,24,26 and 27, the b * i will be discussed in Algorithm 1.
Table 1.27 cases used in the 5-point window algorithm (see [17]) Whenever the optimal b i is non-unique at a node x i , 2 ≤ i ≤ I − 2, we choose the solution b * i to be the real number in the optimal set closest to δ i , that is, median{b u i , b l i , δ i }.In all 27 cases, the solution is either uniquely determined by a simple expression or lies inside the interval (△z i−1 , △z i ) and can be found by a line search.The expression When line search is required, the solution of ( 10) is obtained by solving (The explicit form of this expression is provided in [17].)For the results generated in this paper, we used the secant method (previously used in an analogous way in [3]).Calculation of b * i for i = 0, 1, I − 1 and I is carried out differently, as described in Step 3 of Algorithm 2.
The complete algorithm for calculating second-derivative-based 5-point-window L 1 splines consists of a subroutine (Algorithm 1) for the local window calculation embedded in an "outer loop" (Algorithm 2) as described in the remainder of this section.Recall that the quantities δ i are defined in (2).
and △z i − △z i+1 and determine the case to which these quantities correspond.STEP 2 Calculate b * i using information from [17] contained in Table 1 or described below in this step.Here, c 1 = △z i−2 − △z i−1 and c 2 = △z i+1 − △z i .
-In Cases 1,2,3,5,6,8,9,11,12,20 • Subcase 14-2: If Use a line search to find b * i such that • Subcase 14-3: Use a line search to find b * i such that Return b * i .-Case 15: Use a line search to find b * i such that Use a line search to find b * i such that This transforms Case 23 into Case 15.Obtain the optimal solution b ′ i of the transformed problem.Return b This transforms Case 24 into Case 17. Obtain the optimal solution b ′ i of the transformed problem.Return b -In Case 27, let This transforms Case 27 into Case 14. Obtain the optimal solution b ′ i of the transformed problem.Return b * i = −b ′ i .STEP 3 Stop.
The "outer loop" (Algorithm 2) consists of repeated or parallel application of Algorithm 1.
Algorithm 2 (Analysis-based Algorithm).STEP 1 Given a set of points (x i , z i ), i = 0, 1, . . ., I, calculate △z i , i = 0, . . ., I − 1, and δ i , i = 2, 3, . . ., I − 2, by formulas ( 1) and (2).In this section, we compare computational results generated by the analysis-based algorithm introduced in Section 4 and the widely used primal affine algorithm described in Section 2. These results differ significantly (by more than 10 −4 ) at a number of nodes as shown in Table 2.It was confirmed by hand calculation that the results for the analysis-based algorithm are accurate for the original nondiscretized L 1 spline functional (a confirmation of both the analysis-based algorithm and the code).The differences in the results produced by the primal affine algorithm and the analysis-based algorithm are due mainly to the discretization (midpoint rule with 100 subintervals, as described in Section 2).
Remark On 5-point windows, the primal affine algorithm converges well.However, for global L 1 splines on large data sets, the primal affine algorithm can converge slowly or not at all, especially in bivariate situations, as indicated in Section 6 of [9], a conclusion supported by additional unpublished computational experience of the authors.In such situations, the L 1 spline functional may result in a nearly degenerate nonlinear program and this near degeneracy is retained in the discretized version.The analysis-based algorithm introduced in this paper does not require discretization of the L 1 spline functional and does not require regularization terms to be added to this functional.The algorithm is simple, accurate and inherently parallelizable.Moreover, it is computationally much cheaper than the primal affine algorithm.Both of the algorithms were implemented in C++ 6.0.Computational results for the data set of Section 1 were generated on an IBM laptop running under the Windows XP operating system at 1.66 GHz CPU with 1.50 GB RAM and were presented in Section 3. The computing times of the primal affine algorithm were • 177.7 milliseconds for a 5-point window (sequential calculations), • 335.2 milliseconds for a 7-point window (sequential calculations), • 94.9 milliseconds for global for the second-derivative-based L 1 splines of Figures 2, 3 and 4, respectively.The computing time of the analysis-based algorithm was • 0.0531 milliseconds for a 5-point window (sequential calculations) for the second-derivative-based L 1 spline that corresponds to Figure 2. In these computational experiments, the sequential analysis-based algorithm is thus 177.7/0.0531= 3347 times faster than the sequential primal affine algorithm.The speed-up for the parallel versions of these algorithms would be roughly the same, since each parallel version would be faster than the corresponding sequential version by a factor roughly equal to the number of points in the data set (56 in this case).The sequential analysis-based algorithm is faster than the global primal affine algorithm by a factor of 94.9/0.0531= 1787.The speed-up of the parallel analysis-based algorithm vs. the global primal affine algorithm would, therefore, be roughly a factor of 1787.2 × 56 ≈ 10 5 , an impressive amount.For larger data sets, the speed-up is correspondingly larger.
Remark The computing time for the primal affine algorithm could be reduced by decreasing the number of subintervals K for discretization (via the midpoint rule) of the spline functional from 100, as was used for the results generated for this paper, to, say, 10, or less.However, since the degree to which a discretized spline functional approximates the continuum spline functional decreases as the number of subintervals is decreased, using only a few subintervals is generally not an advantageous choice when accuracy is an issue, as it is here.

Conclusion
The results presented here indicate that the new, analysis-based algorithm for calculating univariate second-derivative-based L 1 interpolating splines on 5-point windows is a highly computationally efficient algorithm for generating a type of L 1 spline that has good shape-preservation properties on the micro scale as well as (generally) the macro scale.The analysis-based algorithm produces computational results that are accurate, in contrast to the widely used primal affine algorithm, which produces computational results that are only approximate.The low computing time and the inherently parallelizable nature of the calculations in the analysis-based algorithm are strong advantages.
The extension of the results in the present paper to algorithms for locally calculated univariate L 1 approximating splines and for locally calculated bivariate L 1 interpolating and approximating splines is a topic for future research.For approximation, either smoothing splines [4,6,8,9,12] or spline fits [8] could be used.In the bivariate case, analytical results for windows consisting of contiguous triangles or rectangles will have to be developed and, on the basis of those analytical results, algorithms analogous to the univariate algorithm introduced in the present paper will need to be created.The success of the algorithm for univariate interpolation of the present paper indicates that further extension in the directions stated here may be fruitful.

Figure 1 .
Figure 1.Data set and linear spline.
Figures 2-13.In each figure, we highlight in dashed boxes the intervals [25, 27], [35, 37] and [45, 55], where differences among the various splines are most prominent.The second-derivative-based 5-point-window spline of Figure 2 and the first-derivative-based and function-value-based global splines of Figures 7 and 10 preserve linearity in the intervals [25, 27], [35, 37] and [45, 48] and avoid extraneous overshoot, extraneous undershoot and extraneous oscillation in the interval [48.1, 55].None of the other nine splines of Figures 2-13 is able to preserve linearity and avoid extraneous overshoot, undershoot and oscillation in all of these intervals.In particular, the splines of Figures 3, 4, 8, 11, 12 and 13 do not preserve linearity well on the intervals [25, 27], [35, 37] and/or [45, 48].The splines of Figures 3, 5, 8 and 11 have extraneous oscillation on the interval [48.1, 50.9].The splines of Figures 6, 9 and 12 have undershoot in [48.1, 50.9].The differences in the large intervals [27.3, 34.7] and [37.3, 44.7] is a separate issue that is discussed in the remark below.Remark The data points immediately to the right and left of the large intervals [27.3, 34.7] and [37.3, 44.7] lie on a quadratic and a cubic function, respectively.The second-and first-derivative-based splines of Figures 2-7 approximate both the quadratic and the cubic function.In contrast, the function-value-based and antiderivative-based splines of Figures 8-13 avoid approximating the quadratic function in [27.

This transforms Case 18 into
Case 15.Obtain the optimal solution b ′ i of the transformed problem.Return b * i = b ′ i .-In Case 23, let This transforms Case 26 into Case 15.Obtain the optimal solution b ′ i of the transformed problem.Return b * i = −b ′ i .

Table 2 .
b * i generated by primal affine algorithm and analysis-based algorithm.