Cleve’s Corner: Cleve Moler on Mathematics and Computing

Splines and Pchips 2

Posted by Cleve Moler,

MATLAB has two different functions for piecewise cubic interpolation, spline and pchip. Why are there two? How do they compare?

Contents

Data

Here is the data that I will use in this post.

x = 1:6
y = [16 18 21 17 15 12]
x =

     1     2     3     4     5     6


y =

    16    18    21    17    15    12

Here is a plot of the data.

set(0,'defaultlinelinewidth',2)
clf
plot(x,y,'-o')
axis([0 7 7.5 25.5])
title('plip')

plip

With line type '-o', the MATLAB plot command plots six 'o's at the six data points and draws straight lines between the points. So I added the title plip because this is a graph of the piecewise linear interpolating polynomial. There is a different linear function between each pair of points. Since we want the function to go through the data points, that is interpolate the data, and since two points determine a line, the plip function is unique.

The PCHIP Family

A PCHIP, a Piecewise Cubic Hermite Interpolating Polynomial, is any piecewise cubic polynomial that interpolates the given data, AND has specified derivatives at the interpolation points. Just as two points determine a linear function, two points and two given slopes determine a cubic. The data points are known as "knots". We have the y-values at the knots, so in order to get a particular PCHIP, we have to somehow specify the values of the derivative, y', at the knots.

Consider these two cubic polynomials in $x$ on the interval $1 \le x \le 2$ . These functions are formed by adding cubic terms that vanish at the end points to the linear interpolatant. I'll tell you later where the coefficients of the cubics come from.

$$ s(x) = 16 + 2(x-1) + \textstyle{\frac{49}{18}}(x-1)^2(x-2) - \textstyle{\frac{89}{18}}(x-1)(x-2)^2 $$

$$ p(x) = 16 + 2(x-1) + \textstyle{\frac{2}{5}}(x-1)^2(x-2) - \textstyle{\frac{1}{2}}(x-1)(x-2)^2 $$

These functions interpolate the same values at the ends.

$$ s(1) = 16, \ \ \ s(2) = 18 $$

$$ p(1) = 16, \ \ \ p(2) = 18 $$

But they have different first derivatives at the ends. In particular, $s'(1)$ is negative and $p'(1)$ is positive.

$$ s'(1) = - \textstyle{\frac{53}{18}}, \ s'(2) = \textstyle{\frac{85}{18}} $$

$$ p'(1) = \textstyle{\frac{3}{2}}, \ \ \ p'(2) = \textstyle{\frac{12}{5}} $$

Here's a plot of these two cubic polynomials. The magenta cubic, which is $p(x)$, just climbs steadily from its initial value to its final value. On the other hand, the cyan cubic, which is $s(x)$, starts off heading in the wrong direction, then has to hurry to catch up.

x = 1:1/64:2;
s = 16 + 2*(x-1) + (49/18)*(x-1).^2.*(x-2) - (89/18)*(x-1).*(x-2).^2;
p = 16 + 2*(x-1) + (2/5)*(x-1).^2.*(x-2) - (1/2)*(x-1).*(x-2).^2;

clf
axis([0 3 15 19])
box on
line(x,s,'color',[0 3/4 3/4])
line(x,p,'color',[3/4 0 3/4])
line(x(1),s(1),'marker','o','color',[0 0 3/4])
line(x(end),s(end),'marker','o','color',[0 0 3/4])

If we piece together enough cubics like these to produce a piecewise cubic that interpolates many data points, we have a PCHIP. We could even mix colors and still have a PCHIP. Clearly, we have to be specific when it comes to specifying the slopes.

One possibility that might occur to you briefly is to use the slopes of the lines connecting the end points of each segment. But this choice just produces zeros for the coefficients of the cubics and leads back to the piecewise linear interpolant. After all, a linear function is a degenerate cubic. This illustrates the fact that the PCHIP family includes many functions.

spline

By far, the most famous member of the PCHIP family is the piecewise cubic spline. All PCHIPs are continuous and have a continuous first derivative. A spline is a PCHIP that is exceptionally smooth, in the sense that its second derivative, and consequently its curvature, also varies continuously. The function derives its name from the flexible wood or plastic strip used to draw smooth curves.

Starting about 50 years ago, Carl de Boor developed much of the basic theory of splines. He wrote a widely adopted package of Fortran software, and a widely cited book, for computations involving splines. Later, Carl authored the MATLAB Spline Toolbox. Today, the Spline Toolbox is part of the Curve Fitting Toolbox.

When Carl began the development of splines, he was with General Motors Research in Michigan. GM was just starting to use numerically controlled machine tools. It is essential that automobile parts have smooth edges and surfaces. If the hood of a car, say, does not have continuously varying curvature, you can see wrinkles in the reflections in the show room. In the automobile industry, a discontinuous second derivative is known as a "dent".

The requirement of a continuous second derivative leads to a set of simultaneous linear equations relating the slopes at the interior knots. The two end points need special treatment, and the default treatment has changed over the years. We now choose the coefficients so that the third derivative does not have a jump at the first and last interior knots. Single cubic pieces interpolate the first three, and the last three, data points. This is known as the "not-a-knot" condition. It adds two more equations to set of equations at the interior points. If there are n knots, this gives a well-conditioned, almost symmetric, tridiagonal $n$ -by- $n$ linear system to solve for the slopes. The system can be solved by the sparse backslash operator in MATLAB, or by a custom, non-pivoting tridiagonal solver. (Other end conditions for splines are available in the Curve Fitting Toolbox.)

As you probably realized, the cyan function $s(x)$ introduced above, is one piece of the spline interpolating our sample data. Here is a graph of the entire function, produced by interpgui from NCM, Numerical Computing with MATLAB.

x = 1:6;
y = [16 18 21 17 15 12];
interpgui(x,y,3)

sppchip

I just made up that name, sppchip. It stands for shape preserving piecewise cubic Hermite interpolating polynomial. The actual name of the MATLAB function is just pchip. This function is not as smooth as spline. There may well be jumps in the second derivative. Instead, the function is designed so that it never locally overshoots the data. The slope at each interior point is taken to be a weighted harmonic mean of the slopes of the piecewise linear interpolant. One-sided slope conditions are imposed at the two end points. The pchip slopes can be computed without solving a linear system.

pchip was originally developed by Fred Fritsch and his colleagues at Lawrence Livermore Laboratory around 1980. They described it as "visually pleasing". Dave Kahaner, Steve Nash and I included some of Fred's Fortran subroutines in our 1989 book, Numerical Methods and Software. We made pchip part of MATLAB in the early '90s.

Here is a comparison of spline and pchip on our data. In this case the spline overshoot on the first subinterval is caused by the not-a-knot end condition. But with more data points, or rapidly varying data points, interior overshoots are possible with spline.

interpgui(x,y,3:4)

spline vs. pchip

Here are eight subplots comparing spline and pchip on a slightly larger data set. The first two plots show the functions $s(x)$ and $p(x)$. The difference between the functions on the interior intervals is barely noticeable. The next two plots show the first derivatives. You can see that the first derivative of spline, $s'(x)$, is smooth, while the first derivative of pchip, $p'(x)$, is continuous, but shows "kinks". The third pair of plots are the second derivatives. The spline second derivative $s''(x)$ is continuous, while the pchip second derivative $p''(x)$ has jumps at the knots. The final pair are the third derivatives. Because both functions are piecewise cubics, their third derivatives, $s'''(x)$ and $p'''(x)$, are piecewise constant. The fact that $s'''(x)$ takes on the same values in the first two intervals and the last two intervals reflects the "not-a-knot" spline end conditions.

splinevspchip

Locality

pchip is local. The behavior of pchip on a particular subinterval is determined by only four points, the two data points on either side of that interval. pchip is unaware of the data farther away. spline is global. The behavior of spline on a particular subinterval is determined by all of the data, although the sensitivity to data far away is less than to nearby data. Both behaviors have their advantages and disadvantages.

Here is the response to a unit impulse. You can see that the support of pchip is confined to the two intervals surrounding the impulse, while the support of spline extends over the entire domain. (There is an elegant set of basis functions for cubic splines known as B-splines that do have compact support.)

x = 1:8;
y = zeros(1,8);
y(4) = 1;
interpgui(x,y,3:4)

interp1

The interp1 function in MATLAB, has several method options. The 'linear', 'spline', and 'pchip' options are the same interpolants we have been discussing here. We decided years ago to make the 'cubic' option the same as 'pchip' because we thought the monotonicity property of pchip was generally more desirable than the smoothness property of spline.

The 'v5cubic' option is yet another member of the PCHIP family, which has been retained for compatibility with version 5 of MATLAB. It requires the x's to be equally spaced. The slope of the v5 cubic at point $x_n$ is $(y_{n+1} - y_{n-1})/2$. The resulting piecewise cubic does not have a continuous second derivative and it does not always preserve shape. Because the abscissa are equally spaced, the v5 cubic can be evaluated quickly by a convolution operation.

Here is our example data, modified slightly to exaggerate behavior, and interpgui modified to include the 'v5cubic' option of interp1. The v5 cubic is the black curve between spline and pchip.

x = 1:6;
y = [16 18 21 11 15 12];
interpgui_with_v5cubic(x,y,3:5)

Resources

A extensive collection of tools for curve and surface fitting, by splines and many other functions, is available in the Curve Fitting Toolbox.

doc curvefit

"NCM", Numerical Computing with MATLAB, has more mathematical details. NCM is available online. Here is the interpolation chapter. Here is interpgui. SIAM publishes a print edition.

Here are the script splinevspchip.m and the modified version of interpgui interpgui_with_v5cubic.m that I used in this post.


Get the MATLAB code

Published with MATLAB® 7.14

2 CommentsOldest to Newest

I’ve been using pchip since I first found it in the early 80s. The funny thing is in most cases, most users won’t see those subtle second derivative discontinuities, but they do prefer the good behavior of a shape preserving tool like pchip. For example, given data with a simple sigmoid shape like that of erf. For an interpolant here, especially if the transition is a sharp one that happens over only a few points, pchip is nice to have. Try this as an example:

x = -4:4;
y = erf(1.5*x);
plot(x,y,’o')

It is a curve that is well behaved, but very difficult for a tool like spline to model.

yhats = interp1(x,y,xhat,’spline’);
yhatp = interp1(x,y,xhat,’pchip’);
plot(x,y,’o',xhat,yhats,’-r’,xhat,yhatp,’-b’)

If that monotone behavior is of paramount importance to you, then there is no choice, use pchip. In fact, that is a fundamental curve shape that I used to see all the time when I worked at Eastman Kodak. See that the spline generates ringing (like that of Gibb’s phenomena) that would be objectionable to a client who needed those flat tails, but pchip does a beautiful job.

However, there are some places where the properties of pchip become a problem. I once had a client who was bothered by those tiny second derivative discontinuities. But more importantly, try this curve out:

x = -4:4;
y2 = exp(-(x-.5).^2);
plot(x,y2,’o')

One basic result/artifact of the shape preserving property of pchip is when the curve goes through a maximum or a minimum, the extremum of the interpolant occurs exactly at the extremal data point. A second characteristic of that tool is when there are two consecutive points with the same value, the interpolant is constant over that interval.

yhat2p = interp1(x,y2,xhat,’pchip’);
yhat2s = interp1(x,y2,xhat,’spline’);
plot(x,y2,’o',xhat,yhat2s,’-r’,xhat,yhat2p,’-b’)

Here I don’t terribly like either of the interpolants. See that pchip did what I would want along the baseline, where we see no ringing. But at the peak, where I know that my function actually hits exactly 1.0, the interpolant misses the behavior I want to see. Even if that exponential had been shifted a bit so it would have peaked at 0.3, pchip would still have produced an interpolant that lacked a well shaped peak.

One solution is to employ a different shape preserving interpolant. For many years I used a rational quadratic spline, based on the ideas of Delbourgo and Gregory.

J.A.Gregory,R.Delbourgo, Piecewise rational quadratic interpolation to monotonic data, IMA Journal of Numerical Analysis, 2(1982),123-130.

A rational quadratic interpolant based on that paper is a bit smoother than pchip, in fact generally C2. Its a bit slower to generate than pchip (which can be blazingly fast, even compared to a cubic spline) and often a bit slower to generate than a cubic spline too. As well, you need to trap for singularities over intervals where the function is constant. There are some other shape preserving interpolants that have been generated over the years, but any interpolant with C2 behavior will please some users. However at a flat topped peak like the data above, these shape preserving interpolants still tend to fail. They simply don’t “know” when it is right to generate a nicely shaped peak that falls between a pair of data points, yet still not create a curve that rings like a bell in other places.

A common source of data with a shape like that of the Gaussian I generated above is seen by anyone who works with illuminant spectral power functions. For example, look at the spectral curve for a fluorescent lamp. Some examples are given in this link:

http://www5.konicaminolta.eu/measuring-instruments/color-light-knowledge/light-sources-and-illuminants.html

The last curve, shown as CIE illuminant F11 is a good example of one with sharp peaks, but also valleys that go down to zero. While we won’t wish to see an interpolant for that curve that goes negative in the valleys, we might also want an interpolant that models the peaks well, yet does not generate severe ringing in other places.

One trick that can work on SOME problems is to recognize that a transformation may be appropriate. Thus, do your interpolation in a domain where ringing is not an issue. For example, on my Gaussian curve, I VERY much prefer the curve I get from a log transformation, then a retransformation via an exponential.

yhat2t = exp(interp1(x,log(y2),xhat,’spline’));
plot(x,y2,’o',xhat,yhat2s,’-r’,xhat,yhat2p,’-b’,xhat,yhat2t,’-g’)

I think most people will prefer the green curve. It has the behavior one expects near a zero baseline, yet it hits my peak quite well too.

Some thoughts I forgot to add in my previous response…

First of all, did I say how happy I am to see Cleve writing these blogs? I’ll bet that many others have been reading and enjoying these blogs just like me. I hope to see many more blogs from Cleve in the future.

Next, I would add that my own SLM tools use an idea under the hood that I took from Fritsch and Carlson to generate monotone approximations. They showed that a cubic segment will be monotone on an interval if the coefficients (when transformed properly) fall inside a convex region with an elliptical boundary. SLM approximates that elliptical region with a slightly smaller (6 sided) polygonal domain that still enforces true monotonicity, creating a set of sufficient constraints that can be implemented with a solver like lsqlin.

These postings are the author's and don't necessarily represent the opinions of MathWorks.