Data Interpolation by Near-Optimal Splines with Free Knots Using Linear Programming

Abstract

The problem of obtaining an optimal spline with free knots is tantamount to minimizing derivatives of a nonlinear differentiable function over a Banach space on a compact set. While the problem of data interpolation by quadratic splines has been accomplished, interpolation by splines of higher orders is far more challenging. In this paper, to overcome difficulties associated with the complexity of the interpolation problem, the interval over which data points are defined is discretized and continuous derivatives are replaced by their discrete counterparts. The $l_{\infty }$-norm used for maximum rth order curvature (a derivative of order r) is then linearized, and the problem to obtain a near-optimal spline becomes a linear programming (LP) problem, which is solved in polynomial time by using LP methods, e.g., by using the Simplex method implemented in modern software such as CPLEX. It is shown that, as the mesh of the discretization approaches zero, a resulting near-optimal spline approaches an optimal spline. Splines with the desired accuracy can be obtained by choosing an appropriately fine mesh of the discretization. By using cubic splines as an example, numerical results demonstrate that the linear programming (LP) formulation, resulting from the discretization of the interpolation problem, can be solved by linear solvers with high computational efficiency and the resulting spline provides a good approximation to the sought-for optimal spline.

1. Introduction

When an exact functional form in a model is unknown, it is often described by measurements represented by available data points and some expected assumptions about the function. Traditionally, data points were interpolated by polynomials more frequently than by spline functions. In a basic way, spline theory deals with a simple but important question of how one uses a lower degree piecewise polynomial function rather than a single larger degree polynomial to represent given data and attendant properties. Spline theory attempts to obtain appropriate problem formulations to obtain the spline, its existence and characterization of solutions, and properties and computation of such functions. While polynomials have many desirable properties, polynomial interpolation suffers from Runge’s Phenomena [1]—a polynomial oscillation at the edges of an interval. With large data, the derivatives of a single polynomial at each of the data points tend to grow, thereby resulting in multiple local extrema. In particular, it has been shown that polynomials, extensively used, for example, in thermometry [2], are not capable of recovering local behavior with good quantifiable accuracy. Spline interpolation, on the other hand, tends to greatly reduce oscillation between data points due to lower degree of polynomials. Cubic splines, with a multitude of applications, in particular, have been shown to be able to reproduce local characteristics well [2]. In addition, computationally, the problem of finding a polynomial interpolating a large number of data points is equivalent to inverting a Vandermonde matrix with a large condition number, thereby rendering the interpolation problem intractable, inefficient, and unstable in computer implementations [3]. On the other hand, the great flexibility of splines, resulting from the use of several lower degree polynomials smoothly connected at optimally chosen points (usually referred to as “knots”), resolves the major problems mentioned above effectively. In spline applications, the following features are desired: 1. representation of the function with lower degree polynomials, which leads to efficient computations as well as a good approximation of the original function, and 2. continuous or essentially bounded derivatives (often second and third derivatives in most applications, for quadratic and cubic splines, respectively). The boundedness of derivatives in the interpolation by splines is natural since one can always have a function interpolating the given point with arbitrarily high derivatives, for example, the second derivative in the context of quadratic splines. In function approximation, therefore, to estimate the most favorable possible error in function approximation, one would like to determine the smallest curvature (the second derivative) possible of a function that interpolates the points by finding the optimal quadratic spline [4].

Many important applications require interpolation by splines. The recent and ongoing COVID-19 pandemic has spurred advances in spline-related research; specifically, splines have been shown to be useful for modeling and forecasting the spread of COVID-19 [5,6,7,8], and COVID-19 trends [9]. In [10], splines were also used to analyze the relationship between the ambient temperature and the transmission of COVID-19. Moreover, splines were shown to be useful for capturing the nonlinear relationship between fasting blood glucose levels and risk of intensive care unit admission for COVID-19 patients [11] in comparison to other conventional linear, dichotomous, or categorical methods. Optimal spline with free knots as well as more recent application areas of splines will be briefly discussed in Section 2. Section 3 provides our linear programming formulation to obtain near-optimal splines with free knots interpolating a given set of data points efficiently and with the desired level of accuracy. In Section 4, numerical results will be provided to demonstrate computational efficiency of the method developed in terms of the CPU time and accuracy, including Example 1 that shows each step and all computations involved in a five-point cubic spline problem for full clarity of each step. In Section 5, concluding remarks are given.

2. Literature Review

Spline Applications. Splines have come in common use in multiple and diverse disciplines such as robotics and statistics. Cubic splines are used in the problem of coordination and planning of the robotic motion [12]. The trajectories for robot motion are frequently described by smooth cubic splines [13]. Cubic splines ($r=3$) are used to characterize the dependency of the Markov model statistical parameters on the conditional parameters [14]. Many filtering techniques require the use of splines. Splines are used to extract curves to model approximation to the shapes of segmented images that result from smoothing the lines of an object in the presence of its shape uncertainty [15]. In digital signal processing, splines are used to generate new sequences with a higher sampling rate (upsampling) [16]. Other recent spline applications include (1) operations research [17,18,19], (2) dynamic programming [20,21], (3) mathematical programming [22,23,24], (4) statistics [25,26,27,28,29,30,31], (5) control theory [31,32,33,34,35,36], (6) computer graphics [37,38], (7) path planning [31], and (8) portfolio selection under fluctuating interest rates [20,21].

General Overview of Optimal Splines with Free Knots. Spline functions were introduced as early as in the nineteenth century by Lobachevsky ([39], p. 165) and later by Schoenberg [40]. Since then, the subject has grown vigorously and many excellent texts are available on the subject: [25,40,41,42,43,44,45].

Studies have shown that the approximations by splines can be significantly improved if knots are allowed to be free rather than at a priori fixed locations (e.g., data points) [46]. When the knots of a spline are fixed a priori (for example, knots coincide with data points), it is called a spline with fixed knots [47,48]. When the knots are not fixed but are determined optimally, it is called a spline with free or variable knots. Generally, as one would expect, the analysis for fixed knots tends to be simpler than for free knots [25] (p. 190). A spline that is required to pass through the data points exactly when data are fairly accurate is called an interpolatory spline; otherwise, it is an approximating spline. The methods given in this section focus on optimal interpolatory splines, with free knots with cubic splines being a special case. A comparative analysis showing advantages of free-knot splines over data-smoothing techniques has been given in [49]. Typically, the degree of smoothness is controlled by both the number and the position of the knots.

It has been shown in [26] that there always exists an optimal spline with $\parallel f^{\left(r\right)}{\parallel }_{\infty }=k^{*}$ for each polynomial comprising the spline and $(n-r-1)$ knots, where n is the number of data points and r is the degree of the spline that interpolates them. It is noted that, for a given value $k^{*}$, an optimal spline with $\parallel f^{\left(r\right)}{\parallel }_{\infty }=k^{*}$ is not necessarily unique (since some of the polynomials comprising the spline (but not all) that interpolate given data points may have $\parallel f^{\left(r\right)}{\parallel }_{\infty }$ that is less than the optimal value $k^{*}$). The apparatus of the proof of the above result is complex, but a simple proof for the quadratic case $r=2$ is given in [50].

3. Problem Formulation

3.1. General Problem Formulation

For a given set of data points ${\left\{x_i,y_i\right\}}_{i=1}^n$ such that A = x${}_1$< … < x${}_n$ = B, a spline function $S\left(x\right)$ of degree r with the knots at $x_{kno{t}_1},\dots ,x_{kno{t}_K}$ is a real valued function defined on the interval [A,B] such that

(i)

$S\left(x\right)$ passes through all the data points, i.e., $S\left(x_i\right)=y_i$;

(ii)

In each interval $\left(A,x_{kno{t}_1}\right),\left(x_{kno{t}_1},x_{kno{t}_2}\right),\dots ,$ $\left(x_{kno{t}_K},B\right)$, where K is the number of knots, $S\left(x\right)$ is given by a polynomial of degree r or less;

(iii)

$S\left(x\right)$ and its derivatives of orders $1,2,\dots ,r-1$ are continuous in the interval $(A,B)$ and

(iv)

The rth derivative of $S\left(x\right)$ can be discontinuous at a finite number of points.

Thus, a spline function of degree r is a piecewise polynomial function with $(r-1)$ continuous derivatives. We could also say that it is a member of $C^{r-1}$ for which the rth derivative is a step function. In other words, it is an rth order indefinite integral of a step function. Polynomial pieces of a spline are “spliced” appropriately at the knots where the rth derivative jumps from one value to another (often with opposite signs), enabling it to turn or change its shape, while all derivatives up to ${(r-1)}^{th}$ order remain continuous, imparting smoothness to the spline. It should be noted that x-coordinates of given data points ${\left\{x_i,y_i\right\}}_{i=1}^n$ are generally not the knot points.

The problem of finding the optimal spline of rth degree is

$$\underset{f\in C^{r-1}}{min}\parallel f^{\left(r\right)}{\parallel }_{\infty }=\underset{f\in C^{r-1}}{min}max|f^{\left(r\right)}|$$

(1)

$$s.t.f\left(x_i\right)=y_i,$$

(2)

where ${\left\{x_i,y_i\right\}}_{i=1}^n$ are given data points.

The formulation (1) and (2), though simple and elegant, is unfortunately not easy to solve. This formulation amounts to the optimization over a Banach space [51] of continuous functions on a compact set $[A,B]$, which is well-known to be notoriously difficult ([31], p. 3). Within (1) and (2), the decision variable is a function f that minimizes $\parallel f^{\left(r\right)}{\parallel }_{\infty }$. While the problem of data interpolation by quadratic splines has been accomplished [4], interpolation by splines of higher orders is far more challenging. Our key contributions to resolve this difficulty are to consider obtaining a near-optimal solution in a computationally efficient manner, and to prove that the near-optimal splines approach optimal splines as discretization mesh is refined. This will be achieved by discretizing the entire space into segments over which the optimization is performed by replacing continuous function with their discrete counterparts and by converting the resulting problem into a linear programming problem that can be solved in polynomial time by using Simplex method [52,53,54], as will be explained next.

3.2. Linear Optimization Formulation for Near-Optimal Solution to (1) and (2)

This section presents a numerical approach for finding near-optimal splines in maximum norm space (typically referred to as L${}_{\infty }$ space) since we minimize $\parallel f^{\left(r\right)}{\parallel }_{\infty }=max\left|f^{\left(r\right)}\right|$. The main idea of the approach is to formulate a linear optimization problem, the solution of which approximates $\parallel f^{\left(r\right)}{\parallel }_{\infty }$, with certain tolerance, of the optimal spline interpolating given n data points ${\left\{x_i,y_i\right\}}_{i=1}^n$.

The idea is presented by using n equidistant data points (note that in the following development of the formulation, there are no requirements of the point equidistance and the data points are not restricted to be equidistant) such that A = x${}_1$, B = x${}_n$ and $x_i-x_{i-1}=H$, $H>0$, $i=2,\dots ,n$, and $\left[A,B\right]={\bigcup }_{i=1}^{n-1}\left[x_i,x_{i+1}\right]$. In order to discretize the problem, we chose to divide each interval $\left[x_i,x_{i+1}\right]$, $i=1,\dots ,n-1$ in m subintervals of equal width. Thus, we partition interval [A, B] into l evenly spaced subintervals, where $l=\left(n-1\right)m$. Note that the number of points defining the l subintervals is $l+1$ and that the points are denoted by ${\left\{{\widehat{x}}_j\right\}}_{j=1}^{(n-1)m+1=l+1}$, such that the relation ${\widehat{x}}_{(i-1)m+1}=x_i$, $i=1,\dots ,n-1$ holds between original data points $x_i$ and data points ${\widehat{x}}_j$ created by subdivisions.

For example, ${\widehat{x}}_1=x_1=A,$ ${\widehat{x}}_{11}=x_2,$ ${\widehat{x}}_{21}=x_3,\dots ,{\widehat{x}}_{(n-1)m+1=l+1}=x_n=B$. It should also be pointed out that the general formula for ${\widehat{x}}_j$ is ${\widehat{x}}_j=B\frac{\left(j-1\right)}{l}+A\frac{\left(1-j+l\right)}{l}$, so that ${\widehat{x}}_1=A,$ ${\widehat{x}}_2=\frac{B}{l}+A\frac{\left(1+l\right)}{l}$, ${\widehat{x}}_3=\frac{2B}{l}+A\frac{\left(2+l\right)}{l}$,…,${\widehat{x}}_{l+1}=B.$ Additionally, while H denotes the width of each interval $\left[x_i,x_{i+1}\right]$ specified in the data, h denotes the width of each smaller interval $\left[{\widehat{x}}_j,{\widehat{x}}_{j+1}\right]$ created after m subdivisions of each interval. In other words, $H=mh$.

After the interval [A, B] is discretized, the derivatives are replaced by their discrete counterparts. Derivatives of a function are often approximated by finite differences. A discrete analog of the rth derivative of a continuous function is the rth order difference divided by h${}^r$. A general formula of the rth order forward difference is based on r consecutive (sub)intervals and is represented as follows [55,56,57]

$${\Delta }_h^r\left[f\right]\left(x\right)=\sum _{i=0}^r{\left(-1\right)}^i\left(\begin{array}{c}r\\ i\end{array}\right)f\left(x+\left(r-i\right)h\right).$$

Central differences ${\delta }_h^r\left[f\right]\left(x\right)={\sum }_{i=0}^r{(-1)}^i\left(\begin{array}{c}r\hfill \\ i\hfill \end{array}\right)f\left(x+\left(\frac{r}{2}-i\right)h\right)$ provide more accurate approximations as compared to forward (or backward) differences: $\frac{{\Delta }_h^r\left[f\right]\left(x\right)}{h^r}=f^{\left(r\right)}\left(x\right)+O\left(h\right)$ and $\frac{{\delta }_h^r\left[f\right]\left(x\right)}{h^r}=f^{\left(r\right)}\left(x\right)+O\left(h^2\right)$; thus, errors in forward differences are of the order h, while in central differences, they are of the order h${}^2$. Essentially, this is due to considering the values of the function at $h/2$ distance away both to the left and to the right of the function in central differences rather than the entire h distance only on one side of the function in forward differences, resulting in a better approximation of the function derivative by central differences.

Briefly, the first-order difference $\frac{f\left(x+\frac{h}{2}\right)-f\left(x-\frac{h}{2}\right)}{h}$ divided by h can be expanded into Taylor series in the following way:

$$\frac{f\left(x+\frac{h}{2}\right)-f\left(x-\frac{h}{2}\right)}{h}\approx$$

$$\left(\frac{f\left(x\right)}{h}+\frac{f^{{}^{\prime }}\left(x\right)}{2}+\frac{f^{{}^{″}}\left(x\right)}{2}h+\frac{1}{48}{f}^{\left(3\right)}\left(x\right)h^2\right)-\left(\frac{f\left(x\right)}{h}-\frac{f^{{}^{\prime }}\left(x\right)}{2}+\frac{f^{{}^{″}}\left(x\right)}{2}h-\frac{1}{48}{f}^{\left(3\right)}\left(x\right)h^2\right)=$$

$$f^{{}^{\prime }}\left(x\right)+\frac{1}{24}{f}^{\left(3\right)}\left(x\right)h^2$$

Additionally,

$$\frac{f\left(x+h\right)-f\left(x\right)}{h}\approx \frac{f\left(x\right)}{h}+f^{{}^{\prime }}\left(x\right)+\frac{1}{2}{f}^{{}^{″}}\left(x\right)h-\frac{f\left(x\right)}{h}=f^{{}^{\prime }}\left(x\right)+\frac{1}{2}{f}^{{}^{″}}\left(x\right)h,$$

from which one can observe that $\frac{f\left(x+\frac{h}{2}\right)-f\left(x-\frac{h}{2}\right)}{h}-f^{{}^{\prime }}\left(x\right)=O\left(h^2\right)$ while $\frac{f\left(x+h\right)-f\left(x\right)}{h}-f^{{}^{\prime }}\left(x\right)=O\left(h\right)$. Similar arguments hold for higher order finite differences.

Even though the forward difference approximates the rth derivative at ${\widehat{x}}_j$ with the error of the order of $O\left(h\right)$, it approximates the rth derivative at midpoint $\frac{{\widehat{x}}_{j+1}+{\widehat{x}}_{j+2}}{2}$ with the error of the order of $O\left(h^2\right)$, as in the discussion of central differences above. Indeed, if we denote ${\widehat{z}}_j=\frac{{\widehat{x}}_{j+1}+{\widehat{x}}_{j+2}}{2}$, then

$$\frac{f\left({\widehat{x}}_j\right)-rf\left({\widehat{x}}_{j+1}\right)+\frac{r(r-1)}{2}f\left({\widehat{x}}_{j+2}\right)+\cdots +{(-1)}^{r}f\left({\widehat{x}}_{j+3}\right)}{h^r}=$$

$$\frac{f\left({\widehat{z}}_j-\frac{3h}{2}\right)-rf\left({\widehat{z}}_j-\frac{h}{2}\right)+\frac{r(r-1)}{2}f\left({\widehat{z}}_j+\frac{h}{2}\right)+\cdots +{(-1)}^{r}f\left({\widehat{z}}_j+\frac{3h}{2}\right)}{h^r},$$

which is the rth-order central difference, and therefore, the overall accuracy of the approximation is of the order of $O\left(h^2\right)$.

Our approach to find a near-optimal solution of (1) and (2) is to minimize the maximum of the rth-order differences directly:

$$\left.\begin{array}{c}\underset{f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_2\right),\dots ,f\left({\widehat{x}}_{(m-1)n+1}\right)}{min}max\left|\frac{{\Delta }_h^r\left[f\right]\left({\widehat{x}}_j\right)}{h^r}\right|\hfill \\ s.t.,f\left({\widehat{x}}_{(m-1)i+1}\right)=y_i,i=1,\dots n\hfill \end{array}\right\}$$

(3)

Note that in this discretized version (and for the rest of the paper), points denoted by $f\left({\widehat{x}}_j\right)$ are decision variables, where ${\widehat{x}}_j,j=1,\dots ,(n-1)m+1$ are endpoints of subintervals.

Since the objective function in formulation (3) is clearly nonlinear due to $max$ and absolute value operators, in order to deal with this nonlinearity, the problem is converted into an equivalent linear counterpart (similar to ([58], p. 28)):

$$\left.\begin{array}{c}\underset{f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_2\right),\dots ,f\left({\widehat{x}}_{(m-1)n+1}\right)}{min}k\hfill \\ s.t.,\frac{{\Delta }_h^r\left[f\right]\left({\widehat{x}}_j\right)}{h^r}\le k,j=1,\dots ,l-r\hfill \\ \frac{{\Delta }_h^r\left[f\right]\left({\widehat{x}}_j\right)}{h^r}\ge -k,j=1,\dots ,l-r\hfill \\ f\left({\widehat{x}}_{(m-1)i+1}\right)=y_i,i=1,\dots ,n\hfill \end{array}\right\}$$

(4)

Since $f\left({\widehat{x}}_{(m-1)i+1}\right)=y_i$ are given as data, the actual decision variables to determine are $f\left({\widehat{x}}_j\right)$, with the exception of $f\left({\widehat{x}}_{(m-1)i+1}\right)$, $i=1,\dots ,n$, although for more convenient implementation, $f\left({\widehat{x}}_{(m-1)i+1}\right)$ can still be treated as decision variables. Additionally, due to linearity, (4) can be solved with much higher efficiency than (3); regardless of the problem instance, formulation (4) remains the LP problem and the LP optimization methods converge within the polynomial time, as discussed above. The reduction of complexity is generally the goal of all linearization processes used to deal with nonlinearity from a computational perspective.

Theorem 1.

The optimal value of formulation (4) approaches $k^{*}$, which is the optimal value of the formulation (1) and (2), as the number of subintervals approaches infinity $\left(m\to \infty \right)$.

Proof. 

As established above, formulation (4)

$$\left.\begin{array}{c}\underset{f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_2\right),\dots ,f\left({\widehat{x}}_{(m-1)n+1}\right)}{min}k\hfill \\ s.t.,\frac{{\Delta }_h^r\left[f\right]\left({\widehat{x}}_j\right)}{h^r}\le k,j=1,\dots ,l-r\hfill \\ \frac{{\Delta }_h^r\left[f\right]\left({\widehat{x}}_j\right)}{h^r}\ge -k,j=1,\dots ,l-r\hfill \\ f\left({\widehat{x}}_{(m-1)i+1}\right) = y_i,i=1,\dots n\hfill \end{array}\right\}$$

is equivalent to formulation (3):

$$\left.\begin{array}{c}\underset{f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_2\right),\dots ,f\left({\widehat{x}}_{(m-1)n+1}\right)}{min}max\frac{\left|{\Delta }_h^r\left[f\right]\left({\widehat{x}}_j\right)\right|}{h^r}\hfill \\ s.t.,f\left({\widehat{x}}_{(m-1)i+1}\right) = y_i,i=1,\dots n\hfill \end{array}\right\}$$

As the number of subintervals approaches infinity $\left(m\to \infty \right)$, the third-order difference approaches the third derivative, so formulation (4) becomes

$$\left.\begin{array}{c}\underset{f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_2\right),\dots }{min}max\left|f^{\left(3\right)}\left({\widehat{x}}_j\right)\right|\hfill \\ s.t.,f\left(x_i\right) = y_i,i=1,\dots n\hfill \end{array}\right\}$$

(5)

Despite having the same objective function, formulation (5) is theoretically not equivalent to the formulation of the original problem (1) and (2). Any countable (even infinitely countable) set of decision variables $\left\{f\left({\widehat{x}}_j\right)\right\}$ is a set of measure zero (a set of points that can be enclosed by intervals with arbitrarily small width), whereas the space over which problem (1) and (2) is optimized is continuous, namely, the entire interval $\left[A,B\right]$. In other words, formulations (1) and (2), and (5), though they have the same objective function, are optimized over different decision variable spaces, namely, discrete and continuous, respectively.

Notwithstanding the differences between formulations (1) and (2), and (5) noted above, their respective optimal values differ by an infinitesimally small value as the number of subintervals approaches infinity. Indeed, for an arbitrarily small $\epsilon >0$, we can always find a large enough m, the number of subintervals (of each interval), such that $\left|{\widehat{x}}_j-{\widehat{x}}_{j+1}\right|<\epsilon$. Since the third derivative of f exists almost everywhere for a cubic spline by definition, we can define central differences at the midpoints of each subinterval $\left({\widehat{x}}_j,{\widehat{x}}_{j+1}\right)$ as $\frac{{\delta }_h^r\left[f\right]\left(\frac{{\widehat{x}}_j+{\widehat{x}}_{j+1}}{2}\right)}{h^r}=f^{\left(3\right)}\left(x\right)+O\left(h^2\right)$. Now, as explained before, note that

$$\frac{{\delta }_h^r\left[f\right]\left(\frac{{\widehat{x}}_j+{\widehat{x}}_{j+1}}{2}\right)}{h^r}=\frac{{\Delta }_h^r\left[f\right]\left({\widehat{x}}_{j-1}\right)}{h^r}$$

Therefore, the optimal solution to (5) differs from the optimal solutions by a quantity of the order of $O\left(h^2\right)$. In other words, for any two arbitrarily close decision variables $f\left({\widehat{x}}_j\right)$ and $f\left({\widehat{x}}_{j+1}\right)$, the respective third derivative of f evaluated at any point between them differs from the third derivative evaluated at either of these points by a quantity that is arbitrary small except possibly at the points of discontinuity of $f^{\left(r\right)}$, which constitute a set of measure zero since $f^{\left(r\right)}$ exists almost everywhere. Therefore, respective optimal values of formulations (5)

$$\left.\begin{array}{c}\underset{f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_2\right),\dots }{min}max|f^{\left(r\right)}\left({\widehat{x}}_j\right)|\hfill \\ s.t.f\left(x_i\right) = y_i,i=1,\dots ,n\hfill \end{array}\right\}$$

and (1) and (2)

$$\left.\begin{array}{c}\underset{f\in C^2}{min}\parallel f^{\left(r\right)}{\parallel }_{\infty }\hfill \\ s.t.f\left(x_i\right)=y_i\hfill \end{array}\right\}$$

differ by a quantity of the order of $O\left(h^2\right)$. □

4. Numerical Testing

The algorithms described above were implemented in CPLEX on an Intel Core i7 CPU Q 720 at 1.60 GHz and 4.00 GB RAM (see the CPLEX code in Appendix A).

4.1. Numerical Examples

In this section, five examples are presented. Example 1 provides a step by step walk-through of the algorithm. Examples 2–4 illustrate the performance of the algorithm in terms of the CPU time for different test functions. Scalability as well as accuracy results are presented in Example 5. In addition, Example 5 gives computational evidence that the optimal value of formulation (4) approaches the optimal value of (1) and (2) as $O\left(h^2\right)$.

Example 1.

Walk-Through Example. Consider a problem with five data points: ($x_1$ = 1, $y_1$ = 3), ($x_2$ = 2, $y_2$ = 5), ($x_3$ = 3, $y_3$ = 4.5), ($x_4$ = 4, $y_4$ = 6.6), and ($x_5$ = 5, $y_5$ = 2.6) and let $r=3$. We want to find a near-optimal cubic spline while illustrating all steps of the formulation very clearly and providing all basic details.

We divide each interval into two subintervals, so $n=5$, $m=2, l=(n-1)m=8$ and $h=0.5$ with the resulting subintervals [${\widehat{x}}_1,{\widehat{x}}_2$] = [$1,1.5$], [${\widehat{x}}_2,{\widehat{x}}_3$] = [$1.5,2$], [${\widehat{x}}_3,{\widehat{x}}_4$] = [$2,2.5$], [${\widehat{x}}_4,{\widehat{x}}_5$] = [$2.5,3$], [${\widehat{x}}_5,{\widehat{x}}_6$] = [$3,3.5$], [${\widehat{x}}_6,{\widehat{x}}_7$] = [$3.5,4$], [${\widehat{x}}_7,{\widehat{x}}_8$] = [$4,4.5$], and [${\widehat{x}}_8,{\widehat{x}}_9$] = [$4.5,5$]. Note that some of the ${\widehat{x}}_j,j=1,\dots ,9$ are problem data values $x_i,i=1,\dots ,5$: ${\widehat{x}}_1=x_1$, ${\widehat{x}}_3=x_2$, ${\widehat{x}}_5=x_3$, ${\widehat{x}}_7=x_4$, and ${\widehat{x}}_9=x_5$. Formulation (4) takes the following form:

$$\begin{array}{c}\underset{\begin{array}{c}f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_2\right),f\left({\widehat{x}}_3\right),f\left({\widehat{x}}_4\right),f\left({\widehat{x}}_5\right),f\left({\widehat{x}}_6\right),f\left({\widehat{x}}_7\right),f\left({\widehat{x}}_8\right),f\left({\widehat{x}}_9\right)\hfill \end{array}}{min}k\\ s.t.-k\le \frac{f\left({\widehat{x}}_1\right)-3f\left({\widehat{x}}_2\right)+3f\left({\widehat{x}}_3\right)-f\left({\widehat{x}}_4\right)}{{0.5}^3}\le k;-k\le \frac{f\left({\widehat{x}}_2\right)-3f\left({\widehat{x}}_3\right)+3f\left({\widehat{x}}_4\right)-f\left({\widehat{x}}_5\right)}{{0.5}^3}\le k;\\ -k\le \frac{f\left({\widehat{x}}_3\right)-3f\left({\widehat{x}}_4\right)+3f\left({\widehat{x}}_5\right)-f\left({\widehat{x}}_6\right)}{{0.5}^3}\le k;-k\le \frac{f\left({\widehat{x}}_4\right)-3f\left({\widehat{x}}_5\right)+3f\left({\widehat{x}}_6\right)-f\left({\widehat{x}}_7\right)}{{0.5}^3}\le k;\\ -k\le \frac{f\left({\widehat{x}}_5\right)-3f\left({\widehat{x}}_6\right)+3f\left({\widehat{x}}_7\right)-f\left({\widehat{x}}_8\right)}{{0.5}^3}\le k;-k\le \frac{f\left({\widehat{x}}_6\right)-3f\left({\widehat{x}}_7\right)+3f\left({\widehat{x}}_8\right)-f\left({\widehat{x}}_9\right)}{{0.5}^3}\le k;\\ f\left({\widehat{x}}_1\right) = y_1=3;f\left({\widehat{x}}_3\right) = y_2=5;f\left({\widehat{x}}_5\right) = y_3=4.5;f\left({\widehat{x}}_7\right) = y_4=6.6;f\left({\widehat{x}}_9\right)=y_5=2.6;\end{array}$$

where ${\widehat{x}}_1=x_1,{\widehat{x}}_2=\frac{x_1+x_2}{2},{\widehat{x}}_3=x_2,{\widehat{x}}_4=\frac{x_2+x_3}{2},{\widehat{x}}_5=x_3,{\widehat{x}}_6=\frac{x_3+x_4}{2},{\widehat{x}}_7=x_4,{\widehat{x}}_8=\frac{x_4+x_5}{2}$, and ${\widehat{x}}_9=x_5.$ The decision variables are: $f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_2\right),f\left({\widehat{x}}_3\right),f\left({\widehat{x}}_4\right),f\left({\widehat{x}}_5\right),f\left({\widehat{x}}_6\right),f\left({\widehat{x}}_7\right),f\left({\widehat{x}}_8\right)$, and $f\left({\widehat{x}}_9\right)$, but $f\left({\widehat{x}}_1\right),f\left({\widehat{x}}_3\right),f\left({\widehat{x}}_5\right),f\left({\widehat{x}}_7\right)$ and $f\left({\widehat{x}}_9\right)$ are known to be $y_1,y_2,y_3,y_4$ and $y_5$, respectively, so the actual decision variables are only $f\left({\widehat{x}}_2\right),f\left({\widehat{x}}_4\right),f\left({\widehat{x}}_6\right)$, and $f\left({\widehat{x}}_8\right)$. The solution is $f\left({\widehat{x}}_2\right)=4.9625,f\left({\widehat{x}}_4\right)=4.4125,f\left({\widehat{x}}_6\right)=5.6625,f\left({\widehat{x}}_8\right)=6.0125$, and $k=10.4.$

Figure 1 and Figure 2 present a near-optimal cubic spline, and the cubic polynomials comprising the spline. The splines are represented by discrete points, interpolating the above data points by using $m=10$.

Example 2.

Exponential Function. Consider a problem with n = 20 and m = 50 and a “test” function $f\left(x\right)=e^{x/10}$.

Example 3.

Trigonometric Function. Consider a problem with n = 20 and m = 50 and a “test” function $f\left(x\right)=sin\left(x\right)$.

Example 4.

Polynomial Function. Consider a problem with n = 20 and m = 50 and a “test” function:

$$f\left(x\right)=\frac{\mathrm{98,442}}{\mathrm{34,295}}+\frac{2725}{\mathrm{157,757}}x+\frac{\mathrm{97,012}}{\mathrm{788,785}}{x}^2-\frac{8448}{\mathrm{788,785}}{x}^3,x\le 8.125$$

$$f\left(x\right)=-\frac{\mathrm{6,798,459}}{\mathrm{34,295}}+\frac{\mathrm{671,965}}{\mathrm{157,757}}x-\frac{\mathrm{314,828}}{\mathrm{788,785}}{x}^2+\frac{8448}{\mathrm{788,785}}{x}^3,x>8.125$$

After being discretized, functions f within Examples 2–4 can be represented as 20 points ${\left\{x_i,y_i\right\}}_{i=1}^{20}$, which are shown in

Table 1. Data for Examples 2–4.

x	Example 2	Example 3	Example 4
$x_1=1$	$y_1=$ 1.1051	$y_1=$ 0.8415	$y_1=$ 3
$x_2=2$	$y_2=$ 1.2214	$y_2=$ 0.9093	$y_2=$ 3.3113
$x_3=3$	$y_3=$ 1.3498	$y_3=$ 0.1411	$y_3=$ 3.74
$x_4=4$	$y_4=$ 1.4918	$y_4=$ −0.757	$y_4=$ 4.2219
$x_5=5$	$y_5=$ 1.6487	$y_5=$ −0.959	$y_5=$ 4.6928
$x_6=6$	$y_6=$ 1.8221	$y_6=$ −0.279	$y_6=$ 5.0883
$x_7=7$	$y_7=$ 2.0137	$y_7=$ 0.6570	$y_7=$ 5.3443
$x_8=8$	$y_8=$ 2.2255	$y_8=$ 0.9894	$y_8=$ 5.3963
$x_9=9$	$y_9=$ 2.4596	$y_9=$ 0.4121	$y_9=$ 5.1947
$x_{10}=10$	$y_{10}=$ 2.7183	$y_{10}=$ −0.5440	$y_{10}=$ 4.7732
$x_{11}=11$	$y_{11}=$ 3.0042	$y_{11}=$ −1	$y_{11}=$ 4.196
$x_{12}=12$	$y_{12}=$ 3.3201	$y_{12}=$ −0.537	$y_{12}=$ 3.5274
$x_{13}=13$	$y_{13}=$ 3.6693	$y_{13}=$ 0.4202	$y_{13}=$ 2.8317
$x_{14}=14$	$y_{14}=$ 4.0552	$y_{14}=$ 0.9906	$y_{14}=$ 2.1731
$x_{15}=15$	$y_{15}=$ 4.4817	$y_{15}=$ 0.6503	$y_{15}=$ 1.6159
$x_{16}=16$	$y_{16}=$ 4.9530	$y_{16}=$ −0.288	$y_{16}=$ 1.2244
$x_{17}=17$	$y_{17}=$ 5.4739	$y_{17}=$ −0.961	$y_{17}=$ 1.0628
$x_{18}=18$	$y_{18}=$ 6.0496	$y_{18}=$ −0.751	$y_{18}=$ 1.1953
$x_{19}=19$	$y_{19}=$ 6.6859	$y_{19}=$ 0.1499	$y_{19}=$ 1.6863
$x_{20}=20$	$y_{20}=$ 7.3890	$y_{20}=$ 0.9129	$y_{20}=$ 2.6

. Note that in most practical cases, the exact functional form in which the cubic spline is used to approximate is not necessarily polynomial and is unknown, so the optimal k${}^{*}$ of optimal splines interpolating given points may differ from the exact k of given test functions (see

Table 2. Computed optimal values and CPU time for Examples 2–4.

	CPU Time (s)	$\widehat{\mathit{k}}$	Actual k
Example 2 (Exp.)	0.25	0.0062	0.0073891
Example 3 (Sine)	0.24	0.8947	1
Example 4 (Polyn.)	0.2	0.0645	0.0642609

). When the test function is polynomial, the optimal k${}^{*}$ is very close to the exact k (see Example 5).

4.2. Scalability Results

To test scalability of our approach, formulation (1) and (2) is used to test problems with varying number of subintervals. The computation times are summarized in

Table 3. Scalability results for Example 5.

No. of Sub-Intervals (m)	$\widehat{\mathit{k}}$	$\frac{{\mathit{k}}^{*}-\widehat{\mathit{k}}}{{\mathit{k}}^{*}}$	CPU Time (s)
2	5.75	4.17%	0.05
4	5.93548	1.08%	0.05
8	5.98373	0.27%	0.03
16	5.99592	0.07%	0.03
32	5.99898	0.02%	0.04
64	5.99974	0.00%	0.08
128	5.99993	0.00%	0.22
256	5.99988	0.00%	0.38
512	5.99999	0.00%	0.27

.

Example 5.

Scalability and Accuracy Testing. Consider five data points (1, 4); (2, 15); (3, 39.75); (4, 78.25); and (5, 124.75) chosen in such a way that the optimal spline has the known value of a third derivative equal to 6.

As one can see in Table 3, the algorithm scales well since increase in the computation time is almost linear for the number of subintervals. Moreover, as the number of subintervals doubles, the accuracy increases (the gap, the relative difference between $\widehat{k}$ and $k^{*}$, decreases) by roughly four times.

4.3. Comparison with Other Methods

Our method is compared with the methods implemented in MATLAB R2021a by using data for Example 1. With a mesh of discretization of 0.01, the default spline interpolation [59] computes splines within 0.0156 seconds and obtains the $l_{\infty }$-norm (max norm) of the third derivative of 12.15, while our method computed splines within 0.071 seconds and obtains the third derivative of 10.767, which is 12.84% lower than the value obtained by MATLAB. The spline interpolation with “optimal” knots [60] obtains the knot location at x = 3, which is one of the data points, and the $l_{\infty }$-norm of the third derivative is the same as that obtained within [59]. The difference is explained by the fact that results in [60] are motivated by works as in [47,48]; therefore, optimality is with respect to knots selected from an a priori given points, while within our method, the position of knots is not restricted as being fixed, thereby leading to a lower objective function.

5. Conclusions

Used in very diverse areas of applications, classical data interpolation by a general spline with free knots is formulated as a linear programming problem to minimize spline $l_{\infty }$-norm (max norm) of the derivative of order r, for reduced complexity, and the problem is efficiently solved by using linear programming solvers. Theoretical convergence to the true optimal splines is proven, indicating that the interpolatory data points obtained by solving the LP problem are arbitrarily close to the points on the true optimal splines. Testing results based on cubic interpolatory splines provide good approximations to optimal splines within a fast computation time. The main benefit of the new method is the significant reduction in complexity to obtain near-optimal splines with near-optimal curvature. Another advantage is that the method has a potential to be generalizable for 2D dimensional (and higher) cases, whereby “interpolatory surfaces” can be obtained by using 2D finite differences. Moreover, as indicated in the Introduction section, since the splines are expected to recover local behavior with good quantifiable accuracy, the proposed method is suitable for finding local as well as global minima of the implied underlying function quickly. The complexity involved in finding a global optimum of the resulting spline based on “merge sort” is “subquadratic” $O(n·log(n\left)\right)$. A limitation of this study is that the near-optimal position of knots cannot be efficiently obtained in online implementations since the method only outputs points belonging to the near-optimal splines, although near-optimal positions of knots can be estimated offline based on the points where the third derivative (for cubic splines) changes its sign. Future research will be to obtain an exact functional representation of the optimal splines as well as the exact positions of the knots.