# algorithm for convex $C^2$ interpolation

Let $$x_0 and $$f_0,f_1,\ldots,f_n$$ be real numbers and $$s_i=(f_i-f_{i-1})/(x_i-x_{i-1}),~~~c_i=(s_{i+1}-s_i)/(x_{i+1}-x_{i-1}).$$

If $$f$$ is a convex function defined on $$[x_0,x_n]$$ with $$f(x_i)=f_i$$ for $$i=0,\ldots,n$$ then all $$c_i$$ are nonnegative. Conversely, this condition guarantees that a convex function $$f$$ with this property exists, namely the piecewise linear interpolant. Necessary and sufficient condition for realizing a twice continuously differentiable function $$f$$ are given in

R. Delbourgo, Shape preserving interpolation to convex data by rational functions with quadratic numerator and linear denominator, IMA J. Numer. Anal. 9 (1989), 123-136.

with an algorithm for constructing such a $$C^2$$ function that involves solving a nonlinear system of equations. Is there a simpler algorithm that does not require solving nonlinear equations?

• The premise is insufficient to guarantee the existence of a $C^2$ function. For example, $c_2>c_1=c_3=0$ would only permit piece-wise linear interpolant from $x_0$ to $x_4$. $c_i$ are all positive would suffice though. – Hans Apr 19 at 2:48
• @Hans: Yes, this is already in Delbourgo's paper, as noted in my statement. – Arnold Neumaier Apr 24 at 7:27

If $$c_i$$'s are all positive, there are infinitely many such convex $$C^2$$ functions. As I have pointed out in my comment above, nonnegativity of $$c_i$$'s is insufficient to guarantee the existence of a $$C^2$$ function. One very simple construction via Bezier curve is as follows.

Draw a straight line through each point $$(x_i,f_i)$$ such that all other points lie above the line. In each interval, construct a quartic Bezier curve as follows. Set the control points of the quadratic Bezier curve resulting from the previously drawn straight lines. Then make the midpoint of each line segment a control point (doubling the control number minus one). Draw the quartic Bezier curve from these control points.

The reason for the construction is that the tangent vector (first derivative) with respect to the parameter ($$t$$ in the Wikipedia article) of a Bezier at an end point is the attached line segment while the second derivative of the curve is the difference of the two closest line segments. We are making the tangent vectors on both sides of the end (data) point coincide and the difference vectors vanish thus equal. You get a $$C^2$$ curve with a continuous first derivative and a continuous second derivative vanishing at each data point.

The above algorithm, proving the existence of the convex $$C^2$$ interpolation function, has a vanishing second derivative at each data point. That makes the first derivative run parallel to the $$x$$ axis each time the curve reaches a data point, making the first derivative wiggly. It does not have to Having proved the existence of an $$C^2$$ interpolation, we can make the first derivative of the convex $$C^2$$ interpolating curve smoother by constructing a higher order Bezier curve by connecting neighbouring data points with many small line segments of almost equal lengths each turning almost a constant angle. This will eliminate the horizontal running points from the first derivative and makes it appear smoother.

• Could you please explain why this Bezier construction gives a $C^2$ function? – Arnold Neumaier Apr 24 at 6:28
• @ArnoldNeumaier: I have added more details as well as the motivation of the construction. Let me know if it is clear. – Hans Apr 25 at 0:09
• OK. As I suspected, quadratic Bezier curves as in the previous version do not provide the required smoothness, but quartics do. This solves my problem from a theoretical point of view. Unfortunately, as in the case of BernM's answer, the derivative is very wiggly, hence it is not suitable for visual applications of the derivative. – Arnold Neumaier Apr 25 at 1:46
• @ArnoldNeumaier: This is just one of the easier examples to give a solution and a constructive proof. The wiggle in the derivative comes from the second derivative being zero. This is not necessary. You can construct higher order Bezier curves to smooth out the derivative while preserving $C^2$. – Hans Apr 25 at 2:06
• I agree that this is a nice costructive proof of sufficiency. – Arnold Neumaier Apr 25 at 2:25

This reference Mulansky, Bernd; Schmidt, Jochen W. Constructive methods in convex C2 interpolation using quartic splines. Numer. Algorithms 12 (1996), no. 1-2, 111–124 may be helpful, but there certainly are more recent ones.

• This paper indeed gives a finite recipe without the need for solving an auxiliary nonlinear problem or optimization problem. Hence it solves my problem as posed. However, the explicit choice given is not visually pleasing (strong oscillations of the curvature) unless the freedom in the placement of the intermediate nodes is used to minimize a shape-quality measure via a quadratic program. Thus from a practical point of view it still uses a nontrivial optimization step. – Arnold Neumaier Apr 17 at 14:54
• Of course, you are right. Here is some theoretical background: There always exists a (strictly) convex interpolating POLYNOMIAL to strictly convex data (all $c_i$ positive). Let S be a finite dimensional subspace of $C^1$. Then there exist strictly convex data which do not admit convexity preserving interpolation from S. Reference: Mulansky, Neamtu: Interpolation and approximation from convex sets. J. Approx. Theory 92 (1998), no. 1, 82–100 – BerndM Apr 26 at 7:43

Edit: After some thoughts, I have concluded that the natural cubic spline can not guarantee to be convex for any convex data. The following algorithm does produce a unique natural cubic spline going through all the data points but it is not guaranteed to be convex. The regression/smoothing algorithm guarantees $$C^2$$ and convexity but not going through each data points.

The Bezier curve construction I supply in the other answer however does provide a simple solution.

You can use cubic spline to not only interpolate these discrete data points but also regress them so that the resulting cubic spline is convex. For the regression problem, minimize $$\begin{equation}\label{eq:splineLoss} L[g] = (1-\lambda)\sum_j w_j(g(t_j)-y_j)^2+\lambda\int_a^b g''(t)^2dt \end{equation}$$ for $$g\in C^2[a,b],\quad g''(t)\ge 0, \quad \lambda\in[0,1].$$ It becomes a quadratic programming problem under the constraint that $$g''(t_j)\ge0,\,\forall j$$.

The procedure is even more specifically spelled out in Berwin A. Turlach, Shape constrained smoothing using smoothing splines

Both of these accounts present the regression (smoothing) algorithm which includes the interpolation which is what OP is asking as a special case.

• but regression destroys the interpolation property – Arnold Neumaier Apr 13 at 19:08
• @ArnoldNeumaier: I am saying this is a more powerful method than just interpolation. It is "not only" good for interpolation which is easier "but also" good for regression. – Hans Apr 13 at 19:11
• The problem is that there is no convex cubic spline interpolant to arbitrary strictly convex data, unless one introduces variable intermediate nodes. This turns the proposed optimization into a nonlinear problem. – Arnold Neumaier Apr 13 at 19:23
• @ArnoldNeumaier: That is not true. There is a unique natural cubic spline interpolation for any discrete data. If the data is convex, the interpolating function is convex.The optimization problem is indeed nonlinear but it is, as I have said, a quadratic programming problem. – Hans Apr 13 at 19:57
• Where can I find a proof of this property of the natural cubic spline? – Arnold Neumaier Apr 13 at 20:52

Too long to comment:

Here's another take, restricting to higher order splines. Suppose the polynomials comprised in the spline are given by $$\{p_{0,1}(x),\cdots,p_{n-1,n}(x)\}$$. Let $$p_{i,i+1}(x) = a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x + \cdots + a^{(n)}_{i,i+1}x^{n-1}$$, $$\forall ~i$$, $$n$$ being odd. The variables here would, of course, be the coefficients of these polynomials (degree, we'll talk about it later). The following points are in order now:

(1) Interpolation constraints $$\{x_i,f_i\}$$ imply linear equality constraints on the coefficients. In particular those will be: $$a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x_i + \cdots + a^{(n)}_{i,i+1}x_i^{n-1} = f_i, ~\forall ~i\\ a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x_{i+1} + \cdots + a^{(n)}_{i,i+1}x_{i+1}^{n-1} = f_{i+1},~\forall ~i.$$ $$C^2$$ requirement will imply another bunch of linear constraints. In particular those will be the following set of equations: $$a^{(1)}_{i-1,i} + 2a^{(2)}_{i-1,i}x_i + \cdots + (n-1)a^{(n)}_{i-1,i}x_i^{n-2} = a^{(1)}_{i,i+1} + 2a^{(2)}_{i,i+1}x_i\cdots + (n-1)a^{(n)}_{i,i+1}x_{i}^{n-3},~\forall ~i,$$ and $$2a^{(2)}_{i-1,i} + 6a^{(3)}_{i-1,i}x_i+ \cdots + (n-1)(n-2)a^{(n)}_{i-1,i}x_i^{n-3} = 2a^{(2)}_{i,i+1} + 6a^{(3)}_{i,i+1}+ \cdots + (n-1)(n-2)a^{(n)}_{i,i+1}x_{i}^{n-2},~\forall ~i.$$

(2) Requirement of convexity over $$[x_i,x_{i+1}]$$, implies that the polynomial given by $$p_{i,i+1}''(x)$$ is positive over $$[x_i,x_{i+1}]$$. A uni-variate polynomial over $$[x_i,x_{i+1}]$$ is non-negative if and only if (see Victorial Powers and Bruce Reznick, Polynomials that Are Positive on an Interval for a nice explanation): $$p_{i,i+1}''(x) = g_{i,i+1}(x) + (x_{i+1}-x)(x-x_i)h_{i,i+1}(x),$$ where $$g_{i,i+1}$$ and $$h_{i,i+1}$$ are SOS polynomials of degree at most $$n$$ and $$n-2$$, respectively. Now, $$g_{i,i+1}(x)=z^\top G_{i,i+1}z$$, where $$G\succeq 0$$ and $$z=[1 ~x \cdots x^{(n-1)/2}]^\top$$. Similarly, $$h_{i,i+1}(x)=y^\top H_{i,i+1}y$$, where $$H_{i,i+1}\succeq 0$$ and $$y=[1 ~x \cdots x^{(n-1)/2-1}]^\top$$. The matrices $$G_{i,i+1}$$ and $$H_{i,i+1}$$ are additional variables. Comparing coefficients on either side of the equation given by: $$2a^{(2)}_{i,i+1} + 6a^{(3)}_{i,i+1}+ \cdots + (n-1)(n-2)a^{(n)}_{i,i+1}x_{i}^{n-2} = g(x) + (x_{i+1}-x)(x-x_i)h(x),$$ yields linear affine equations in coefficients of $$p_{i,i+1}(x)$$, $$G$$ and $$H$$. For brevity, let these equations be given as $$\mathcal{L}_{i,i+1}\left(a^{(0)}_{i,i+1},\cdots,a^{(n)}_{i,i+1},G_{i,i+1},H_{i,i+1}\right)=0$$. Note that you get a set of linear equations $$\forall ~i$$.

(3) An SDP solver (CVXPY, or alike) can now be used to find a feasible solution, i.e.: $$\min ~~1~~\mbox{subject to}\\ a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x_i + \cdots + a^{(n)}_{i,i+1}x_i^{n-1} = f_i, ~\forall ~i\\ a^{(0)}_{i,i+1} + a^{(1)}_{i,i+1}x_{i+1} + \cdots + a^{(n)}_{i,i+1}x_{i+1}^{n-1} = f_{i+1},~\forall ~i\\ a^{(1)}_{i-1,i} + \cdots + (n-1)a^{(n)}_{i-1,i}x_i^{n-2} = a^{(1)}_{i,i+1} + \cdots + (n-1)a^{(n)}_{i,i+1}x_{i}^{n-3},~\forall ~i\\ 2a^{(2)}_{i-1,i} + \cdots + (n-1)(n-2)a^{(n)}_{i-1,i}x_i^{n-3} = 2a^{(2)}_{i,i+1} + \cdots + (n-1)(n-2)a^{(n)}_{i,i+1}x_{i}^{n-2},~\forall ~i\\ \mathcal{L}_{i,i+1}\left(a^{(0)}_{i,i+1},\cdots,a^{(n)}_{i,i+1},G_{i,i+1},H_{i,i+1}\right)=0\\ G_{i,i+1},H_{i,i+1} \succeq 0, \forall i.$$ Lack of a feasible point, however, does not mean there does not exist a $$C^2$$ function, of a spline with degree higher than the one chosen or outside the realm of splines.

(4) The choice of the degree of the polynomial is something I don't have a clear answer for. I would assume that a degree greater than 8 should work -- for every polynomial piece, there are 2 interpolation constraints, 4 for $$C^2$$, and 2 LMI constraint.

(5) A disadvantage of this method is that its not practical to work with several thousand data points.

Hope it helps.

• You have not provided an (operable) algorithm. Neither does your scheme guarantee a $C^2$ convex function interpolation. – Hans Apr 21 at 6:41
• @Hans, thanks for the comment. Will type out all the necessary details in the next edit. However, I do not see why it does not guarantee $C^2$ and convexity. To clarify (1)$C^2$ is guaranteed as long as the poly parts and two derivatives are continuous at interpolation points. This translates the linear equalities as I've mentioned in the answer. And (2)a uni-variate function is convex if and only if its second derivative is non-negative at every point over the desired range. This translates to convex LMI, which comes from a SOS poly-rep of uni-variate positive polys over a compact interval. – DSM Apr 21 at 7:44
• I am saying you need to guarantee the continuity of $f, f', f''$ at the nodes $x_i$'s. Not all sums of squares polynomials ensure that. You have not provided any algorithm much less a proof to ensure the existence of such sums of squares polynomials as the second derivatives in the intervals and first derivatives at the nodes. Also, it would be great if you could use the full name of a nomenclature before using its acronym. – Hans Apr 21 at 7:57
• Have added a lot of more details. Hope it helps. For the acronyms part, had to do because of the character limit. Sorry for the inconvenience. – DSM Apr 21 at 10:56
• Thank you for providing the details, particularly the crucial paper. I like the formulation. However, it still lacks a proof of the existence of a feasible $G$ and $H$ for this SDP problem. My answer using the Bezier curve or Bernstein polynomial interestingly provides just such an existence proof. – Hans Apr 21 at 17:32