algorithm for convex $C^2$ interpolation Let $x_0<x_1<\ldots<x_n$ 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?
 A: 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.
A: 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.
A: 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$.
These are described and proved in detail in
P.J. Green, Bernard. W. Silverman, Nonparametric Regression and Generalized Linear Models: A roughness penalty approach (Chapman & Hall/CRC Monographs on Statistics & Applied Probability Book 58).
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.
A: 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.       
