Why is Lagrange interpolation numerically unstable?

Here is my understanding of the polynomial interpolation problem:

Interpolating by inverting the Vandermonde matrix is unstable because the Vandermonde matrix is ill-conditioned, so "difficult" to invert. Lagrange interpolation is much more clever as it computes the same polynomial using a different basis where the matrix is diagonal, so easy to invert.

Furthermore, the algorithm for computing Lagrange polynomials is straightforward and I don't see anything that makes it unstable.

So why is it commonly claimed in textbooks that Lagrange interpolation is numerically unstable?

Thanks!

Solution 1:

Instability does not always refer to numerical issues in an algorithm. It can refer to the effect of perturbations in data inputs to the output solution. For example, in the picture below we fit a 4th degree polynomial to two sets of 5 data points. The data sets match except that the data point at the origin has been perturbed slightly (can think of it as measurement error). Although the perturbation is small, the change in the Lagrange polynomial is large.

This behavior is not caused by numerical inaccuracies in the algorithm to compute the coefficients of the polynomial, it is inherent in the interpolation problem itself.

Solution 2:

There are a few issues to deal with in Lagrange interpolation. First and formost is the problem of the [Runge Phenomenon][1](this was a hyperlink originally, but apparently my reputation isnt high enough). This essentially boils down to the fact that on evenly spaced points (and for other point distributions) the interpolation problem is ill posed even in exact arithmetic.

In order to guarantee stability of the interpolation scheme, we must choose our interpolation points judiciously. Chebyshev points (essentially you cluster points towards the endpoints of the interval by projecting down from evenly spaced points on the unit circle) map the local extrema of the interpolating polynomial to the collocation points, essentially insuring that the Runge phenomenon is avoided. There are other such sets of points that avoid the Runge Phenomenon also (Legendre, Gaussian) but for the purposes of this post we'll stick to Chebyshev points.

Now the canonical form of Lagrange interpolation requires O(n^2) multiplications and subtractions to evaluate the interpolation at a point. This cost means that for certain problems the numerical errors in evaluating your interpolant can grow exponentially (I think, numerical instabilities are a weak point of mine) with the amount interpolation points you choose. So even if we interpolate on Chebyshev points, the Lagrange formulation still suffers from numerical instabilities for very large values of n.

This problem has a beautiful resolution in the form of the Barycentric Lagrange formula. This form of the interpolation is the equivalent to the canonical form, but offers us a saving grace from numerical instabilities. The main problem in barycentric form is that in order to calculate the weights w_i, we need to take a large product.

But due to the work of Salzer in 1978, when interpolating on Chebyshev points, the Barycentric weights are known exactly. This removes any problems we had with numerical instability.

Berrut and Trefethen have a paper from 2004 that reviews all of these issues in depth, and shows that Barycentric Lagrange Interpolation on Chebyshev points is damn near optimal both analytically and numerically.

(I'm only an undergrad with a heavy interest in approximation theory, but the above is taken pretty much directly from my numerical analysis class this morning.)