Why is the cubic Spline Interpolation wiggly in flat areas?

Dear community

I am using Matlab's griddedInterpolant() to fit 2-D data with non-uniform samples. Please consider the following example:

clc; clear; close all; format shorteng;
Z = [728.8979e-003   636.7566e-003   606.6035e-003   565.2786e-003   550.3203e-003   544.8502e-003   542.6799e-003   541.8001e-003   541.4345e-003
   521.0865e-003   462.6134e-003   468.5421e-003   472.2405e-003   475.2979e-003   476.8306e-003   477.0386e-003   477.5314e-003   477.6013e-003
   478.5540e-003   431.6547e-003   451.0393e-003   462.2832e-003   466.5187e-003   469.1578e-003   470.0577e-003   470.9129e-003   470.9733e-003
   454.4821e-003   422.4798e-003   447.9669e-003   459.6772e-003   465.3374e-003   467.3769e-003   468.7602e-003   469.0878e-003   469.5432e-003
   444.1855e-003   420.1184e-003   446.5284e-003   459.0328e-003   464.4620e-003   466.8222e-003   467.8873e-003   468.6215e-003   468.6755e-003
   441.7790e-003   418.7854e-003   445.8779e-003   458.3203e-003   463.9560e-003   466.7335e-003   467.9822e-003   468.2045e-003   468.2045e-003
   440.6139e-003   418.6750e-003   445.4984e-003   457.8616e-003   463.8469e-003   466.8397e-003   467.4772e-003   467.8873e-003   468.6773e-003
   440.1721e-003   418.5457e-003   445.2773e-003   458.1930e-003   463.8622e-003   466.3328e-003   467.3769e-003   468.3087e-003   468.4529e-003
   439.5239e-003   418.0573e-003   445.3255e-003   457.9355e-003   463.7742e-003   466.3328e-003   467.3769e-003   468.3087e-003   468.4529e-003];
[X,Y] = ndgrid(logspace(log10(210e-9),log10(100e-6),9));

surf(X,Y,Z,'FaceAlpha',0.5)  
set(gca,'XScale','log')
set(gca,'YScale','log')
title('Original')

F = griddedInterpolant(X,Y,Z,'spline');
[XX,YY] = ndgrid(logspace(log10(210e-9),log10(100e-6),3*9));
Z_approx = F(XX,YY);
figure
surf(XX,YY,Z_approx,'EdgeColor','black')
set(gca,'XScale','log')
set(gca,'YScale','log')
title('Approximation')

As you can see, the original data is mostly flat except for small X and Y values. However, when I compute the cubic Splines and use a finer grid to sample data points from the interpolated curve, then this flat area becomes wiggly.

a) Why? Shouldn't the cubic Spline Interpolation be able to come up with lower degree polynomials (constants) at those points? I mean, with degree zero polynomials (constants) it should also be possible to fulfil the C2 smoothness condition.

b) The documentation states that the 'spline' argument uses not-a-knot end conditions. Does Matlab also support other end conditions?

c) Any suggestions on how to get rid of that wiggly surface? Are there better methods? I also tried with multi-variate least-squares but the results aren't really better. Also, due to the nature of where the data comes from I am unable to have uniformly sampled data.

Answer 1

[John D'Errico]

Should not the cubic spline be able to be ABSOLUTELY flat in some areas? NO. Sorry, but that will not be the case. A cubic spline is twice differentiable everywhere. And that forces the spline to oscillate into regions where the data is otherwise flat. A simple example in 1-d will help to illustrate what happens.

x = 1:10;
y = [ones(1,5),zeros(1,5)];
spl = spline(x,y);
fnplt(spl,'r')
hold on
plot(x,y,'go')
grid on
hold off

So the spline oscillates everywhere, even though you might hope that in the intervals [1,5] and [6,10], the function could arguably be absoutely flat, based on the data we have.

This is NOT a question of end conditions. If I used different end conditions, this issue would still arise.

What happens is in order for the function to be twice differentiable EVERYWHERE, (and three times differentiable at all points that are not breakpoints) we get a function that is forced to oscillate.

There are alternatives to a simple spline though. pchip and makima, for example are subtly different interpolants that are not twice differentiable at the data points. They are only once differentiable at those locations. For the particular data I have made up, they will probably produce the same results, but that will not always be the case.

spl2 = makima(x,y);

fnplt(spl2,'b')
hold on
plot(x,y,'go')
grid on
hold off

Now look at the derivatives of these functions.

spldiff = fnder(spl,2);
fnplt(spldiff,'r')
hold on
spl2diff = fnder(spl2,2);
fnplt(spl2diff,'b')
grid on
hold off

Do you see that at the locations 5 and 6, spl2 (in blue) has a second derivative discontinuity? At the same time, the second derivative was always purely continuous for the spline.

You may want for a spline to behave however you wish, but mathematics is cruel. Mathematics follows its own hard set of rules.

Again, end conditions are completely, utterly, absolutely irrelevant here.

If you want to understand what is happening, consider the curve as produced by makima. Just slightly to the left of x==6, and just slightly to the right of x==5, the curve must obviously have a maximal amount of curvature. Correct? Heavily curved fuctinos at those locations cannot suddenly turn into absoutely flat functions right next door, at least not and be twice continuously differentiable.

And while the data I have used here is far more extreme that that which you have, it is also only a 1-d probnlem. Things are always more difficult in 2-d.

Are there better methods available? The simplest answer is to just use what is known as a bilinear interpolant. That is, essentially a lower order spline. In this case, a LINEAR spline in two dimensions. Interp2 produces this, when the 'linear' option is chosen. Would you still complain? Possibly. I can show you cases where such a bilinear interpolant produces artifacts that will confuse most people. Is the bilinear interpolant produced by interp2 better? In some respects, yes, in others, no. You will no longer see a smooth surface, but more of a faceted result. But where your data is flat, the surface produced will now still be flat.

[Dominik Hiltbrunner]

Great explaination! Thank you.