Approximation of Irregular Geometric Data by Locally Calculated Univariate Cubic $L^1$ Spline Fits

Abstract

$L^1$ splines have been under development for interpolation and approximation of irregular geometric data. We investigate the advantages in terms of shape preservation and computational efficiency of calculating univariate cubic $L^1$ spline fits using a steepest-descent algorithm to minimize a global data-fitting functional under a constraint implemented by a local analysis-based interpolating-spline algorithm on 5-node windows. Comparison of these locally calculated $L^1$ spline fits with globally calculated $L^1$ spline fits previously reported in the literature indicates that the locally calculated spline fits preserve shape on the average slightly better than the globally calculated spline fits and are computationally more efficient because the locally-calculated-spline-fit algorithm can be parallelized.

1 Introduction

Over the past decade and more, $L^1$ splines have been under development for interpolation and approximation of irregular geometric data, that is, data with abrupt changes in magnitude, direction and/or spacing. $L^1$ splines, which are based on nonlinear programming that implements a geometric principle of shape preservation, have the advantage of being free of the extraneous oscillation common in other interpolants and approximants such as unconstrained and constrained polynomial and rational splines, including those with penalties, a posteriori filtering and fairing. Recent advances in univariate and bivariate $L^1$ interpolating splines (see [1, 2, 5, 7, 8, 12, 14, 16] and [10, 13], resp.) have been accompanied by some but fewer advances in univariate and bivariate $L^1$ approximating splines (see [9, 11] and [3], resp.). This is unfortunate, since approximation is, from a practical point of view, more important than interpolation.

In 2007, Auquiert, Gibaru and Nyiri discovered that $L^1$ interpolating splines calculated on small sets of 5 nodes preserve linearity of 3-node subsets of the 5 nodes [2]. This discovery led to the development of univariate $L^1$ interpolating splines calculated locally on 5-point windows [7, 8, 14, 16]. These locally calculated $L^1$ interpolating splines have the dual advantage of enhanced shape preservation and sharply reduced computing time, the latter of which is due an algorithm in which parallel analysis replaces sequential raw computing. The local-calculation approach that has been successful for interpolating splines has, however, not previously been investigated for approximating splines. This paper is a first step toward developing locally calculated $L^1$ approximating splines to replace the previously globally calculated $L^1$ approximating splines.

Approximation by $L^1$ (and other) splines can be carried out using spline fits or smoothing splines. Spline fits are calculated by minimizing a data fitting functional over a manifold of interpolating splines. Smoothing splines are calculated by minimizing a linear combination of a data fitting functional and an interpolating spline functional over a spline space. Smoothing splines depend on a “balance parameter” that represents the relative weight of fitting the data vs. that of smoothing out local variations in the data. Many researchers have investigated the balance parameter for classical smoothing splines (see [15] and references therein) but practically relevant, optimal balance parameters for large classes of data are not yet available. For $L^1$ smoothing splines, no theoretical information is available at all [9, 11]. For this reason, spline fits, which do not require the user to provide a balance parameter, are of high interest and are the approximating splines that we choose here.

The precise goal of the present paper is to investigate the advantages in terms of shape preservation and computational efficiency of calculating univariate cubic $L^1$ spline fits using a steepest-descent algorithm to minimize a global data-fitting functional under a constraint implemented by the local analysis-based interpolating-spline algorithm of [16] on 5-node windows. We will compare these “locally calculated $L^1$ spline fits” with the globally calculated $L^1$ spline fits presented in [11].

2 Cubic $L^1$ Spline Fits: Definitions and Algorithms

Let the data to be approximated be $(\hat{x}_m,\hat{z}_m)$, $m = 1, 2, \dots ,M$, where all of the $\hat{x}_m$ are in a finite real interval $[x_0,x_I]$. Let there be given strictly monotonically increasing but otherwise arbitrary spline nodes $x_i$, $i = 0, 1, . . . , I$. Let the function values and the first derivatives of the spline fit at these nodes be denoted by $z_i$ and $b_i$, respectively. Let $\mathbf{z}$ denote $(z_0,z_1,\dots ,z_I)$ and $\mathbf{b}$ denote $(b_0,b_1,\dots ,b_I)$. Let $h_i=x_{i+1}-x_i$ and $t=\frac{x-x_i}{h_i}-\frac{1}{2}$ for $-0.5 < t < 0.5$. We consider $C^1$-smooth piecewise cubic polynomials $z(x)$ with nodes $x_i$, $i = 0, 1, \dots , I$, of the (Hermite) form

$$\begin{aligned} \begin{array}{l} z(x) = \left( \frac{1}{2} - \frac{3}{2}t + 2t^3\right) z_i + \left( \frac{1}{2} + \frac{3}{2}t - 2t^3\right) z_{i+1} \\ \quad \quad \quad \quad + \left( \frac{1}{8} - \frac{1}{4} t - \frac{1}{2} t^2 + t^3\right) h_i b_i + \left( - \frac{1}{8} - \frac{1}{4} t + \frac{1}{2} t^2 + t^3\right) h_i b_{i+1} \end{array} \end{aligned}$$

(1)

in each interval $[x_i,x_{i+1}]$, $i=0,1,\dots ,I-1$. A cubic $L^1$ spline fit is a function $z(x)$ of form (1) that minimizes the data fitting functional

$$\begin{aligned} F(\mathbf{z},\mathbf{b}):=\sum _{m=1}^M \left| z(\hat{x}_m) - \hat{z}_m \right| \end{aligned}$$

(2)

over a manifold of cubic $L^1$ interpolating splines. The free parameters in a cubic $L^1$ spline fit are the function values $z_i=z(x_i)$, $i=0,1,\dots ,I$, at the nodes. The first derivatives $b_i=\text{ d }z/\text{ d }x(x_i)$, $i=0,1,\dots ,I$ of the spline fit are dependent variables that are calculated locally or globally from the $z_i$ by minimizing an $L^1$ interpolating spline functional as described in the next two paragraphs.

A globally calculated cubic $L^1$ interpolating spline is a function of form (1) that minimizes the global functional

$$\begin{aligned} \int \limits _{x_0}^{x_I} \left| \frac{\text{ d }^2 z}{\text{ d }x^2} \right| \, \text{ d }x. \end{aligned}$$

(3)

[11] over all functions of form (1). In this minimization process, the $z_i$ are fixed and the free variables are the first derivatives $b_i=\text{ d }z/\text{ d }x(x_i)$, $i=0,1,\dots ,I$ of the spline. Nonuniqueness is handled by choosing the feasible $b_i$ that are closest to 0, a “regularization procedure” that was used for the computational examples in [11] with which we will compare the new results in this present paper. [Other regularizations could be used. In [16], nonuniqueness is handled by choosing the feasible $b_i$ that is closest to $(z_{i+1}-z_{i-1})/(x_{i+1}-x_{i-1})$].

A locally calculated cubic $L^1$ interpolating spline [16] is a function of form (1) in which each of the first derivatives $b_i$, $i=3,4,\dots ,I-3$ for (1) is calculated by minimizing the local 5-node-window spline functional

$$\begin{aligned} \int \limits _{x_{i-2}}^{x_{i+2}} \left| \frac{\text{ d }^2 z}{\text{ d }x^2} \right| \, \text{ d }x. \end{aligned}$$

(4)

over the free variables $b_k$, $k=i-2,i-1,i,i+1,i+2$, with $z_k$, $k=i-2,i-1,i,i+1,i+2$, fixed. The first derivatives $b_0$, $b_1$ and $b_2$ are calculated by minimizing functional (4) with $i=2$ and the first derivatives $b_{I-2}$, $b_{I-1}$ and $b_I$ are calculated by minimizing functional (4) with $i=I-2$. Nonuniqueness is handled by choosing the feasible $b_i$ that is closest to 0.

The coefficients of a globally calculated $L^1$ spline fit are the solution of a two-level or “bilevel” optimization problem, namely, that of minimizing the $\ell ^1$ data-fitting functional (2) for the $z_i$ subject to a constraint that the $b_i$ are functions of the $z_i$ that are determined by minimizing the global $L^1$ interpolating spline functional (3). Many approaches for solving bilevel optimization problems exist, including branch-and-bound, complementary-pivoting, descent, penalty-function and trust-region approaches [6]. However, there are still many challenges, especially for nonlinear bilevel problems, a category that includes calculation of the coefficients of globally calculated $L^1$ spline fits. In [11], a “Lagrange-multiplier primal-affine algorithm” was used to solve this bilevel problem. On each iterative step, this algorithm applies a primal-affine approach to reduce the minimization of the data-fitting functional (2) and the minimization of the global $L^1$ interpolating spline functional (3) to weighted least-squares problems. One formulates the linear systems that correspond to these two problems, links these two linear systems by Lagrange multipliers and computes the next iterate by solving the resulting larger linear system. A theoretical proof of of convergence is not available but computational results indicate good performance. At the same time, the method is relatively computationally expensive due to the global structure of the calculations.

The recent development of locally calculated $L^1$ interpolating splines, which are better at shape preservation and much more computationally efficient than globally calculated $L^1$ interpolating splines, opens up opportunities for approximation by locally calculated $L^1$ spline fits that we wish to begin to explore in this paper. The opportunities are threefold: (1) increased computational efficiency by shifting from minimization of a global $L^1$ interpolating spline functional to parallel minimization of many local $L^1$ interpolating spline functionals, (2) increased shape-preservation capability and (3) potential availability of theoretical proof of convergence due to knowledge of the analytical structure of the $b_i$ when they are calculated locally from the $z_i$. There exists no known analogue of this analytical structure for globally calculated $L^1$ splines. The present paper focuses on the first two of these opportunities.

The algorithm that we propose for locally calculated $L^1$ spline fits is a steepest-descent algorithm for minimizing $F(\mathbf{z},\mathbf{b}(\mathbf{z}))$ with $\mathbf{z}$ as a free, unconstrained vector of variables. The initial iterate from which the steepest-descent process starts is the linear spline fit of the data. The linear spline fit is the function $\zeta (x)$ that minimizes the data fitting functional

$$\begin{aligned} \sum _{m=1}^M \left| \zeta (\hat{x}_m) - \hat{z}_m \right| \end{aligned}$$

(5)

over the manifold of $C^0$-smooth piecewise linear polynomials $\zeta (x)$ with nodes $x_i$, $i = 0, 1, \dots , I$, of the form

$$\begin{aligned} \begin{array}{l} \zeta (x) = \frac{(x_{i+1}-x)}{h_i} z_i + \frac{(x-x_i)}{h_i} z_{i+1} \end{array} \end{aligned}$$

(6)

in each interval $[x_i,x_{i+1}]$. The linear spline fit is geometrically close to the cubic $L^1$ spline fit (an observation for which theoretical proof is not yet available), which is a factor in promoting good performance of the steepest-descent algorithm that we introduce here.

To proceed from an approximant $\mathbf{z}^{(k)}$ on the $k$th step to the next approximant $\mathbf{z}^{(k+1)}$, we first compute a numerical gradient $\bigtriangledown ^{(k)}$ of $F(\mathbf{z},\mathbf{b}(\mathbf{z}))$ with respect to $\mathbf{z}$ at $\mathbf{z}=\mathbf{z}^{(k)}$. The derivative $\text{ d }\mathbf{b}/\text{ d }\mathbf{z}$ that occurs in this gradient is calculated by analytical differentiation of the $\mathbf{b}(\mathbf{z})$ of [7, 16] but with nonuniqueness resolved by choosing the feasible $b_i$ that is closest to 0 instead of by choosing the $b_i$ closest to $(z_{i+1}-z_{i-1})/(x_{i+1}-x_{i-1})$. This formula is not fully applicable at boundaries where the analytical derivative of $\mathbf{b}$ with respect to $\mathbf{z}$ is discontinuous. This aspect of the algorithm will be adjusted in the future. However, since the initial iterate, the linear spline fit, is close to the final solution, that is, to the cubic $L^1$ spline fit, there are not a large number of boundaries with discontinuous derivatives between the initial iterate and the final solution and the current version of the algorithm performs well. After calculating $\bigtriangledown ^{(k)}$, we use golden section line search to numerically calculate a step length $\gamma ^{(k)}$ that minimizes $F(\mathbf{z}^{(k)}-\gamma \bigtriangledown ^{(k)}, \mathbf{b}(\mathbf{z}^{(k)} -\gamma \bigtriangledown ^{(k)}))$ over $0\le \gamma \le 1$. Finally, we set $\mathbf{z}^{(k+1)}=\mathbf{z}^{(k)} -\gamma ^{(k)} \bigtriangledown ^{(k)}$. If $F(\mathbf{z}^{(k+1)}, \mathbf{b}(\mathbf{z}^{(k+1)}))>F(\mathbf{z}^{(k)}, \mathbf{b}(\mathbf{z}^{(k)}))$, we stop and output $\mathbf{z}^{(k)}$. Otherwise, we increase $k$ by 1 and repeat the process. A maximum number of iterations (100 in our computational experiments) is set for potential cases of slow convergence or non-convergence.

3 Computational Results

Computations were carried out on a LENOVO Thinkpad T410 laptop with an Intel Core i5 2.67GHz CPU and 4.00GB RAM. The software environment was Microsoft Windows 7 Professional and MATLAB R2012a. Locally calculated $L^1$ spline fits were computed using a MATLAB code and globally calculated $L^1$ spline fits were computed using a legacy C++ code, both codes implemented sequentially. On the T410 laptop, MATLAB was 112.3 times slower than a C++ code for a task consisting of $10^7$ additions and $10^7$ multiplications. In reporting CPU times below, we have normalized the CPU times required by the MATLAB code to corresponding C++ code CPU times by division by 112.3.

We present in Table 1 computational results for locally and globally calculated cubic $L^1$ spline fits for “data-and-node sets” 1, 2 and 4, resp., of [11]. These cubic $L^1$ spline fits are shown in Figs. 1, 2, 3, 4, 5, 6. In these figures, the data points are represented by the symbol $+$ and the spline nodes are represented by the symbol $\circ$.

Fig. 1

Locally calculated cubic $L^1$ spline fit for data-and-node set 1

Fig. 2

Globally calculated cubic $L^1$ spline fit for data-and-node set 1

Fig. 3

Locally calculated cubic $L^1$ spline fit for data-and-node set 2

Fig. 4

Globally calculated cubic $L^1$ spline fit for data-and-node set 2

Fig. 5

Locally calculated cubic $L^1$ spline fit for data-and-node set 4

Fig. 6

Globally calculated cubic $L^1$ spline fit for data-and-node set 4

Table 1 Computational results

Comparison of the locally calculated $L^1$ spline fits of Figs. 1, 3 and 5 with the corresponding globally calculated $L^1$ spline fits of [11] in Figs. 2, 4 and 6 indicates that the locally calculated spline fits preserve shape in these computational experiments on the average slightly better than the globally calculated spline fits, a situation that is reflected in the $\ell ^1$ errors for data-and-node sets 1 and 4. This observation is analogous to previous results in the literature that indicate that locally calculated $L^1$ interpolating splines preserve shape slightly better than globally calculated $L^1$ interpolating splines. The locally calculated $L^1$ spline fit for data-and-node set 2 preserves shape better than the globally calculated $L^1$ spline fit in the interval [4, 5] but not as well in the intervals [1, 3.5] and [7, 9.5]. The sequential computing times for locally and globally calculated $L^1$ spline fits reported in Table 1 are of the same order of magnitude. (No further conclusion about relative computing time is appropriate at present due to the small amount of data and the fact that the MATLAB computing time for the locally calculated $L^1$ spline fits was artificially normalized to a corresponding C++ computing time.) The sequential structure of the current locally-calculated-$L^1$-spline-fit code will be parallelized in the future, which will lead to a large increase in computational efficiency. The globally-calculated-$L^1$-spline-fit code cannot be parallelized.

4 Conclusion

The evidence that we have presented here is preliminary but is already sufficient to suggest that further investigation of the local-calculation approach for $L^1$ spline fits may well lead to significant advances in shape-preserving approximation. While it is commonly expected that one must trade off shape-preservation capability vs. computational efficiency, the results presented here indicate that the local-calculation approach with its parallel structure may allow both shape-preservation capability and computational efficiency to be enhanced simultaneously. One option that has not yet been investigated but that will be investigated in the future is replacement of the global data-fitting functional (2) by local $\ell ^1$ data-fitting functionals, perhaps on the same 5-node windows that are used for the local $L^1$ interpolating spline functional (4).

Theoretical proof of convergence of the algorithm is hampered by the complexity of the dependence of $\mathbf{b}$ on $\mathbf{z}$. The authors expect that the interplay of the regularization strategy for the $L^1$ interpolating spline functional (that is, how a specific solution is chosen when the minimum occurs over a whole interval) with the resulting smoothness of $\mathbf{b}(\mathbf{z})$ will play a large role in this proof when the authors develop it in the future.

The results about univariate $L^1$ spline fits that we have presented here are not merely a step toward improved univariate approximation but also a basis for bivariate and multivariate approximation. In contrast to the situation with univariate conventional splines, univariate $L^1$ splines generalize in natural ways to higher dimensions, a statement that is valid not merely for interpolating splines but also for spline fits. Further univariate results will pave the way for future investigation of the shape-preservation capability and computational efficiency of bivariate and multivariate $L^1$ spline fits.