## Cubic Interpolation (Interpolation, part II)

In a previous entry, we had a look at linear interpolation and concluded that we should prefer some kind of smooth interpolation, maybe a polynomial.

However, we must use a polynomial of sufficient degree so that neighboring patches do not exhibit very different slopes on either side of known points. This pretty much rules out quadratic polynomials, because polynomials of the form $ax^2+bx+c$ are only capable of expressing (strictly) convex (or concave) patches. A quadratic piece-wise function would look something like:

Here, the solid lines represent the part of the patch used for interpolation, while the dashed line show the extent of the quadratic patch. We see that this yields a saw-tooth pattern, clearly something we do not want for “smooth” interpolation.

Since a cubic polynomial is capable of (slightly more) expressiveness, namely it can have two bumps (in opposing directions), and uses four points (instead of just three as with the quadratic patches), it captures possibly better the oscillations in the (unknown) underlying function. On the same points, maybe cubic patches would yield the following interpolation:

Maybe the slope doesn’t quite match on either side of the points, but at least, it doesn’t change as wildly than the quadratic patches.

Let again $\{(x_i,y_i)\}_{i=1}^n$ be a series of known points between which we wish to interpolate. So between any two points $(x_i,y_i)$ and $(x_{i+1},y_{i+1})$, we will use a cubic polynomial to interpolate. However, a degree $m$ polynomial needs $m+1$ points to be uniquely defined (that’s the unisolvence theorem), and we will need, in addition to $(x_i,y_i)$ and $(x_{i+1},y_{i+1})$, the points $(x_{i-1},y_{i-1})$ and $(x_{i+2},y_{i+2})$. We need to find the parameters to the equation

$ax^3+bx^2+cx+d=y$

that passes through all four points. That gives us four equations, four unknowns that we must solve:

$ax_{i-1}^3+bx_{i-1}^2+cx_{i-1}+d=y_{i-1}$

$ax_{i}^3+bx_{i}^2+cx_{i}+d=y_{i}$

$ax_{i+1}^3+bx_{i+1}^2+cx_{i+1}+d=y_{i+1}$

$ax_{i+2}^3+bx_{i+2}^2+cx_{i+2}+d=y_{i+2}$

But writing equations this way is cumbersome, and not propitious to fast and efficient solution. Noticing that

$ax^3+bx^2+cx+d=y$

can be rewritten as

$\hat{y}(x)=\left[~x^3~x^2~x~1~\right]~\left[\begin{array}{c} a\\ b\\ c\\ d\end{array}\right]$

That is, as a dot product. We can rewrite the four equations as a matrix/vector system of the form $Ma=y$:

$\left[\begin{array}{cccc} x_{i-1}^3 & x_{i-1}^2 & x_{i-1} & 1\\ x_{i}^3 & x_{i}^2 & x_{i} & 1\\ x_{i+1}^3 & x_{i+1}^2 & x_{i+1} & 1\\ x_{i+2}^3 & x_{i+2}^2 & x_{i+2} & 1 \end{array}\right]~\left[\begin{array}{c} a\\ b\\ c\\ d\end{array}\right]=\left[\begin{array}{c} y_{i-1}\\ y_{i}\\ y_{i+1}\\ y_{i+2}\end{array}\right]$

and we solve for $\left[~a~b~c~d~\right]$. Fortunately, it’s not too hard:

$\left[\begin{array}{c} a\\ b\\ c\\ d\end{array}\right]= \left[\begin{array}{cccc} x_{i-1}^3 & x_{i-1}^2 & x_{i-1} & 1\\ x_{i}^3 & x_{i}^2 & x_{i} & 1\\ x_{i+1}^3 & x_{i+1}^2 & x_{i+1} & 1\\ x_{i+2}^3 & x_{i+2}^2 & x_{i+2} & 1 \end{array}\right]^{-1} \left[\begin{array}{c} y_{i-1}\\ y_{i}\\ y_{i+1}\\ y_{i+2} \end{array}\right]$

Matrix inverses are really painful to compute in general (not even when using pen-and-paper computation) but if we have the special case where $x_{i-1}=-1$, $x_{i}=0$, $x_{i+1}=1$, and $x_{i+1}=2$, then the matrix simplifies to

$M=\left[\begin{array}{cccc} -1 & 1 & -1 & 1\\ 0 & 0 & 0 & 1\\ 1 & 1 & 1 & 1\\ 8 & 4 & 2 & 1 \end{array}\right]$

with the inverse

$M^{-1}=\displaystyle\frac{1}{6}\left[\begin{array}{cccc} -1 & 3 & -3 & 1\\ 3 & -6 & 3 & 0\\ -2 & -3 & 6 & -1\\ 0 & 6 & 0 & 0 \end{array}\right]$

and solving for the parameters are obtained by

$\left[\begin{array}{c} a\\ b\\ c\\ d \end{array}\right]=M^{-1}~\left[\begin{array}{c} y_{i-1}\\ y_{i}\\ y_{i+1}\\ y_{i+2}\end{array}\right]$

*
* *

We can compute the product $M^{-1}y$ in much less than $4\times{4}$ (scalar) products. First, some entries are zero, are are $\pm{1}$ and therefore reduce to simple addition/subtraction; then we notice that some are similar except for sign. With constant folding, we can get something efficient.

*
* *

So we have a rather efficient way of fitting a cubic through four points, but the cubic patches still do not quite fit together as nicely as they should. In particular, the derivative around $(x_i,y_i)$ is not the same whether we consider the patch interpolating between $(x_{i-1},y_{i-1})$ and $(x_i,y_i)$ and the patch interpolating between $(x_i,y_i)$ and $(x_{i+1},y_{i+1})$. So how do we force them to match?

Well, we can add the constraints explicitly in the equation system and solve for different parameters. One simple way of doing so is to use the Hermite splines.

To be continued…