Cubic spline interpolation

Motivation

In an interpolation problem, we're given a set of points (we'll be using 2D points X,Y throughout this post) and are asked to estimate Y values for Xs not in this original set, specifically for Xs that lie between Xs of the original set (estimation for Xs outside the bounds of the original set is called extrapolation).

As a concrete example, consider the set of points (0, 1), (1, 3), (2, 2); here they are plotted in the usual coordinate system:

Interpolation is estimating the value of Y for Xs between 0 and 2, given just this data set. Obviously, the more complex the underlying function/phenomenon, and the fewer original points we're given, interpolation becomes more difficult to do accurately.

There are many techniques to interpolate between a given set of points. Polynomial interpolation can perfectly fit N points with an N-1 degree polynomial, but this approach can be problematic for large a N; high-degree polynomials tend to overfit their data, and suffer from other numerical issues like Runge's phenomenon.

Instead of interpolating all the points with a single function, a very popular alternative is using Splines, which are piece-wise polynomials. The idea is to fit a low-degree polynomial between every pair of adjacent points in the original data set; for N points, we get N-1 different polynomials. The simplest and best known variant of this technique is linear interpolation:

Linear interpolation has clear benefits: it's very fast, and when N is large it produces reasonable results. However, for small Ns the result isn't great, and the approximation is very crude. Here's the linear spline interpolation of the Sinc function sampled at 7 points:

We can certainly do much better.

How about higher-degree splines? We can try second degree polynomials, but it's better to jump straight to cubic (third degree). Here's why: to make our interpolation realistic and aesthetically pleasing, we want the neighboring polynomials not only to touch at the original points (the linear splines already do this), but to actually look like they're part of the same curve. For this purpose, we want the slope of the polynomials to be continuous, meaning that if two polynomials meet at point P, their first derivatives at this point are equal. Moreover, to ensure smoothness and to minimize needless bending [1], we also want the second derivatives of the two polynomials to be equal at P. The lowest degree of polynomial that gives us this level of control is 3 (since the second derivative of a quadratic polynomial is constant); hence cubic splines.

Here's a cubic spline interpolating between the three points of the original example:

And the Sinc function:

Because of the continuity of first and second derivatives, cubic splines look very natural; on the other hand, since the degree of each polynomial remains at most 3, they don't overfit too much. Hence they're such a popular tool for interpolation and design/graphics.

All the plots in this post have been produced by JavaScript code that implements cubic spline interpolation from scratch. Let's move on to learn how it works.

Setting up equations for cubic spline interpolation

Given a set of N points, we want to produce N-1 cubic polynomials between these points. While these are distinct polynomials, they are connected through mutual constraints on the original points, as we'll see soon.

More formally, we're going to define N-1 polynomials in the inclusive range :

For each polynomial, we have to find 4 coefficients: a, b, c and d; in total, for N-1 polynomials we'll need 4N-4 coefficients. We're going to find these coefficients by expressing the constraints we have as linear equations, and then solving a system of linear equations. We'll need 4N-4 equations to ensure we can find a unique solution for 4N-4 unknowns.

Let's use our sample set of three original points to demonstrate how this calculation works: (0, 1), (1, 3), (2, 2). Since N is 3, we'll be looking for two polynomials and a total of 8 coefficients.

The first set of constraints is obvious - each polynomial has to pass through the two points it's interpolating between. The first polynomial passes through the points (0, 1) and (1, 3), so we can write the equations:

The second polynomial passes through the points (1, 3) and (2, 2), resulting in the equations:

We have 4 equations, and need 4 more.

We constrain the first and second derivatives of the polynomials to be equal at the points where they meet. In our example, there are only two polynomials that meet at a single point, so we'll get two equations: their derivatives are equal at point (1, 3).

Recall that the first and second derivatives of a cubic polynomial are:

The equation we get from equating the first derivatives is:

Or, expressed as a linear equation of all coefficients:

Similarly, the equation we get from equating the second derivatives is:

Expressed as a linear equation of all coefficients:

This brings us to a total of 6 equations. The last two equations will come from boundary conditions. Notice that - so far - we didn't say much about how our interpolating polynomials behave at the end points, except that they pass through them. Boundary conditions are constraints we create to define how our polynomials behave at these end points. There are several approaches to this, but here we'll just discuss the most commonly-used one: a natural spline. Mathematically it says that the first polynomial has a second derivative of 0 at the first original point, and the last polynomial has a second derivative of 0 at the last original point. In our example:

Substituting the second derivative equations:

We have 8 equations now:

To restate the obvious - while our example only uses 2 polynomials, this approach generalizes to any number. For N original points, we'll interpolate with N-1 polynomials, resulting in 4N-4 coefficients. We'll get:

• 2N-2 equations from setting the points these polynomials pass through
• N-2 equations from equating first derivatives at internal points
• N-2 equations from equating second derivatives at internal points
• 2 equations from boundary conditions

For a total of 4N-4 equations.

The code that constructs these equations from a given set of points is available in this file.

Solving the equations

We now have 8 equations with 8 variables. Some of them are trivial, so it's tempting to just solve the system by hand, and indeed one can do it very easily. In the general case, however, it would be quite difficult - imagine interpolating 10 polynomials resulting in 36 equations!

Fortunately, the full power of linear algebra is now at our disposal. We can express this set of linear equations as a matrix multiplication problem , where A is a matrix of coefficients, x is a vector of unknowns and b is the vector of right-hand side constants:

Solving this system is straightforward using Gaussian elimination. Our JavaScript implementation does this in a few steps:

• Performs Gaussian elimination to bring A into row-echelon form, using the algorithm outlined on Wikipedia. This approach tries to preserve numerical stability by selecting the row with the largest (in absolute value) value for each column [2].
• Further transforms the resulting matrix into reduced row-echelon form (a.k.a. Gauss-Jordan elimination)
• Extracts the solution.

In our example, the solution ends up being the vector (-0.75, 0, 2.75, 1, 0.75, -4.5, 7.25, -0.5); therefore, the interpolating polynomials are:

Performing the interpolation itself

Now that we have the interpolating polynomials, we can generate any number of interpolated points. For all x between 0 and 1 we use , and for x between 1 and 2 we use . In our JavaScript code this is done by the doInterpolate function. We've already seen the result:

Code

The complete code sample for this post is available on GitHub. It includes functions for constructing equations for cubic splines from an original set of points, code for solving linear equations with Gauss-Jordan elimination, and a demo HTML page that plots the points and linear/spline interpolations.

The code is readable, heavily-commented JavaScript with no dependencies (except D3 for the plotting).

An additional demo that uses similar functionality is line-plotting; it plots arbitrary mathematical functions with optional interpolation (when the number of sampled points is low).