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


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





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


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…

3 Responses to Cubic Interpolation (Interpolation, part II)

  1. […] To be continued… Share this:StumbleUponDiggRedditTwitterMoreFacebookEmailLike this:LikeBe the first to like this post. […]

  2. […] previous posts, we discussed linear and cubic interpolation. So let us continue where we left the last entry: […]

  3. […] the last four installments of this series, we have seen linear interpolation, cubic interpolation, Hermite […]

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )


Connecting to %s

%d bloggers like this: