Interpolating Polynomials

John D'Errico is back today to talk about interpolating polynomials.

Contents

Introduction

I talked in my last blog about polynomial regression models. We saw some utility in polynomial models, and that increasing the order of the polynomial caused the errors to decrease. How far can we go? Can we reduce the errors to zero?

At some point, I might choose to increase the degree of the polynomial to one less than the number of data points. Assuming that my data have no replicated points, this is an interpolating polynomial that fits our data exactly, at least to within the double precision accuracy of our computations. A straight line can pass through any two points, a quadratic passes through three points, a cubic hits four points exactly, etc. If you prefer to think in terms of statistical degrees of freedom, if you have n points, a polynomial of order n-1 (with n coefficients) passes through your data exactly.

I've started this all out by talking about polynomial modeling in a least squares sense, or regression analysis. But really, I wanted to talk about interpolation, as opposed to the approximations provided by polynomial regression. Perhaps this was an indirect approach, but one of the things I feel important is to distinguish interpolation from the more general modeling/curve fitting tools used in mathematics. This difference is often unappreciated by many users of MATLAB.

Let's start by sampling five points from an exponential function.

format short
x = (-2:1:2)';
y = exp(x);
plot(x,y,'o')
grid on

Exact Fit

While it is true that polyfit gives an exact fit with an order n-1 polynomial, a direct method is more efficient.

Build an interpolating polynomial using vander coupled with the use of backslash.

P = vander(x)\y
xev = linspace(-2,2,500)';
yev = polyval(P,xev);

plot(x,y,'ko',xev,yev,'r-')
grid on
title('Quartic polynomial through 5 points on exp(x)')

P =
    0.0492
    0.2127
    0.4939
    0.9625
    1.0000

Many books teach you to use the Lagrange form for interpolation. I won't get into that here because I don't really advise the use of polynomials in general.

Interpolation should yield zero residuals at the data points

Verify that the residuals are essentially zero at the data points.

format long g
pred = polyval(P,x) - y

pred =
    -1.11022302462516e-016
    -1.11022302462516e-016
                         0
                         0
                         0

Does the Data Have Replicates?

Suppose I have two points with the same independent variable value, but distinct values for the dependent variable? Since for any given x, a single-valued function returns only a single value for the prediction, this must cause a problem. In this example, what value should I expect for an interpolant when x = 2? Should f(x) be 3 or 4?

x = [0 1 2 2 3 4]';
y = [1 2 3 4 5 6]';
plot(x,y,'-o')

You can see a hint of the problem if you look at the rank of the Vandermonde matrix. The rank should have been 6, but the replicated point causes a singularity.

rank(vander(x))

ans =
     5

This singularity results in a rather useless polynomial.

P = vander(x)\y

Warning: Matrix is singular to working precision.
P =
   NaN
   NaN
   NaN
  -Inf
   Inf
     1

Average the Replicates

If the data truly is not single-valued, the common approach is to use a parametric relationship. (This is something to save for a possible future blog.)

Much of the time however, I just want to do something simple, like average the replicate points. consolidator, a tool I contributed to the File Exchange, can do just that. If you've never found or looked at the file exchange, please take some time to wander through the great tools there.

High Order Polynomials Not a Sure Bet

Sometimes polynomials don't do so well. It is important to know your data, know your problem. Here the data comes from a function known to be monotone increasing. Polynomials don't like to be monotone, so forcing a high order polynomial to interpolate such a function is a sure way to fail.

x = (-5:1:5)';
y = erf(x);

P = vander(x)\y
xev = linspace(-5,5,1000)';
yev = polyval(P,xev);

plot(x,y,'ko',xev,yev,'r-')
grid on
title('Interpolating polynomial through 11 points on erf(x)')

P =
     3.49245965480804e-020
      2.0991424886728e-005
    -1.58112816820803e-018
      -0.00119783140998794
     2.17562811018059e-017
        0.0239521737329468
    -9.91270557701033e-017
        -0.211403900317105
     1.85037170770859e-016
          1.03132935951897
                         0

Runge's Phenomenon

Runge found this example of a function that causes high order interpolating polynomials to fare very poorly.

x = (-5:1:5)';
y = 1./(1+x.^2);
P = vander(x)\y
xev = linspace(-5,5,1000)';
yev = polyval(P,xev);
ytrue = 1./(1+xev.^2);

plot(x,y,'ko',xev,yev,'r-.',xev,ytrue,'g-')
legend('Data points','Interpolating polynomial','y = 1./(1+x.^2)')
grid on
title('Runge example')

P =
    -2.26244343891403e-005
     1.04238735515992e-019
       0.00126696832579186
    -4.47734652027112e-018
       -0.0244117647058824
     5.61044648135202e-017
          0.19737556561086
    -2.38551402729053e-016
         -0.67420814479638
     1.38777878078145e-016
                         1

Numerical Problems

You also need to worry about numerical issues in high order polynomial interpolants. The matrices you would create are often numerically singular.

For example, the rank of the matrix vander(1:n) should be n. However, this fails when working in finite precision arithmetic.

rank(vander(1:10))
rank(vander(1:20))
rank(vander(1:30))

ans =
    10
ans =
     7
ans =
     6

Centering and scaling your data can improve the numerical performance of a high order polynomial to some extent, but that too has a limit. You are usually better off using other interpolation methods.

What can you do? You want to be able to interpolate your data. Polynomials seem like a good place to look, but they have their issues. High order polynomial interpolation often has problems, either resulting in non-monotonic interpolants or numerical problems. You might consider other families of functions to build your interpolant, for example trig or bessel functions, or orthogonal polynomials. Much depends on the relationship you want to model. KNOW YOUR PROBLEM!

By the way, for those who are interested in the methodology of modeling and curve fitting, I list a few ideas in this thread that should be of interest.

We can still do much with polynomials however. They form the basis for many modeling tools.

A nice trick is to go simple. Why use a single high order polynomial to interpolate our data? Instead, break the problem up and use smaller segments of low order polynomials, pieced together. In my next blog I'll begin talking about these piecewise interpolants. Until then please tell me here where you have found interpolating polynomials of use, or if you have found problems that they cannot solve.

Get the MATLAB code

Published with MATLAB® 7.6