The points that ruin the fit

A Taylor polynomial spends everything it has on being exactly right at one point. An interpolating polynomial spreads the same budget over several: given $n$ points of a curve, there is exactly one polynomial of degree $n-1$ through all of them, and it is a reasonable expectation that using more points gives a better approximation.

The polynomial is correct. It is the unique one of its degree through those eleven points, it was computed in a numerically sound way, and it agrees with the function exactly wherever it was asked to. Between the samples near the ends of the interval it is wrong by $0.25$, on a function whose entire range is $[0, 1]$.

Adding points makes it worse #

The natural response is to take more samples, and it is exactly wrong.

The measured sequence is not subtle. With $5$, $9$, $11$, $15$, $21$, $31$ and $41$ evenly spaced samples the worst errors are $0.236$, $0.224$, $0.249$, $0.337$, $0.604$, $1.87$ and $6.33$. It falls slightly, turns, and then grows without bound: the interpolating polynomials through equally spaced points do not converge to this function at all, and the failure is not slow.

This is Runge’s phenomenon, published by Carl Runge in 1901, and it is worth being precise about what it does and does not resemble. It is not the Taylor series failing outside its radius — that was a series, this is a sequence of interpolants, and no term of either is a term of the other. The two failures happen on the same function and their causes are related without being the same.

The error is a product of distances #

The formula for interpolation error explains everything above, and it has one factor in it that nobody controls and one that is entirely a matter of choice.

For a polynomial through $n$ nodes $x_1, \ldots, x_n$,

$$f(x) - p(x) = \frac{f^{(n)}(\xi)}{n!}\,\prod_{i=1}^{n}(x - x_i)$$

for some $\xi$ in the interval — the same shape as Lagrange’s remainder, with the single power $x^{n}$ replaced by a product of distances to the nodes.

That product is the free parameter. At a point of the interval it is the product of the distances to every sample, and its size depends entirely on where the samples were put. Put them evenly and the product is small in the middle — where a point is close to several nodes — and large near the ends, where every node is on one side and the distances multiply up. The imbalance is enormous: for twenty-one evenly spaced nodes across an interval, the product near the ends exceeds its value in the middle by a factor of thousands.

The other factor, $f^{(n)}(\xi)/n!$, is a property of the function and grows for $1/(1+x^2)$ at just the rate needed to defeat the factorial — which is the same statement as the radius of convergence being $1$, the same two singularities at $\pm i$ doing the work, and the reason the phenomenon is often said to have a complex-analytic explanation.

So a function whose derivatives grow slowly is safe whatever the nodes, and one whose derivatives grow fast is at the mercy of the product. The sine is the first kind.

Moving the points #

The product of distances is the half that is chosen rather than given, so the question is where to put the nodes to make its worst value as small as possible. That question has an exact answer.

Take equally spaced angles on a semicircle and project them onto the diameter — the construction that turns a circle into a wave, used here on a set of points rather than on a moving one. A point near the top of the circle moves a long way horizontally for a small change of angle, and one near the side barely moves at all, so the projections bunch up at the two ends of the interval. Those are the Chebyshev points, and the product of distances over them is as flat as it can possibly be made — its worst value is smaller than for any other choice of $n$ nodes, by a theorem rather than by experiment.

The measured sequence for Chebyshev nodes at $5, 9, 11, 15, 21, 31, 41$ points is $0.189$, $3.32\times10^{-2}$, $1.67\times10^{-2}$, $3.66\times10^{-3}$, $3.51\times10^{-4}$, $7.14\times10^{-6}$, $1.45\times10^{-7}$. It converges, geometrically, on the function the evenly spaced version diverges from.

Same function, same degrees, same arithmetic, opposite conclusions, and the only difference is where fifty points were placed on a line.

The amplification, which is one number #

There is a way of stating the difference that does not mention any particular function, and it is the sharper statement.

Interpolation is linear in the data: change a sample value and the polynomial changes by that amount times a fixed basis polynomial. So the worst factor by which an error in the data can be magnified anywhere in the interval is a single number depending only on the nodes — the Lebesgue constant.

The evenly spaced constant grows like $2^n/(n\log n)$ — exponentially — and the Chebyshev constant like $\log n$, which is the slowest growth this collection has a name for and is slow enough that doubling the node count adds about $0.22$ to the factor. At twenty-one nodes the two are four orders of magnitude apart.

That number settles what is really wrong. An interpolation scheme whose amplification is $30\,000$ is not usable on measured data whatever function is behind it, because a rounding error in the fifth digit becomes an error in the first. The Runge function is a demonstration rather than the disease: the evenly spaced rule was broken before any particular function was chosen, and the choice of function only decided how visibly.

The constant also puts a bound on the whole method, and the bound is the reason it is the right quantity to compute. Whatever the function, the interpolant’s error is at most $(1 + \Lambda)$ times the error of the best polynomial of that degree — the one chosen by an oracle to minimise the worst discrepancy, which interpolation has no way of finding. So a scheme with $\Lambda = 2.9$ is within a factor of four of optimal on every function there is, and a scheme with $\Lambda = 30\,000$ carries no useful guarantee at all.

That inequality is what turns a measurement into a theorem. The seven error values quoted above are about one function; the amplification bound is about all of them, and it says the Chebyshev result was never going to be an accident of the example. A quantity that depends on the nodes alone is worth far more than a table of errors that depends on the function too, which is the general reason to look for a constant that separates the two.

It also explains the shape of the failures. The amplification is worst near the ends of the interval, so that is where the oscillation appears, and it is why every equally spaced picture on this page is excellent in the middle third and hopeless outside it.

Where the two failures differ #

The Taylor series and the interpolants both fail on $1/(1+x^2)$, and it is worth separating them because the fixes are different.

The series fails outside $|x| < 1$ and cannot be repaired by taking more terms, since more terms is exactly what makes it worse there. It can be repaired by moving the centre, which moves the disc.

The interpolants fail near the ends of whatever interval is chosen, and the failure gets worse with more points. It cannot be repaired by moving anything, because there is no centre; it is repaired by redistributing the points, and after redistribution the convergence is geometric on the whole interval — including well outside $|x| < 1$, where the series has nothing to offer at all.

So Chebyshev interpolation succeeds where the Taylor series cannot go. That is the useful comparison and it is the reverse of the usual reading, in which Runge’s phenomenon is presented as a cautionary tale about polynomials. Polynomials are fine. Degree $n$ can approximate this function to $10^{-7}$ across an interval five wide, and the only question was ever which $n+1$ points to ask about.

What is still not fixed #

Chebyshev nodes are only optimal for one measure. They minimise the worst value of the product of distances, which is the right objective when the derivative factor is unknown and could be anything. If more is known about the function a different placement can beat them, and if the objective is a least-squares error rather than a worst case the optimal points are different again.

The nodes are dictated, and the data may not be. All of the above assumes the function can be sampled wherever wanted. Measured data arrives where it arrives, usually evenly spaced because that is how instruments work, and none of this is available: interpolating such data by a high-degree polynomial is exactly the situation the Lebesgue constant condemns, and the answer there is not better nodes but a different method — a spline, which is a low-degree polynomial on each piece, or a least-squares fit of a low degree through all of them.

And convergence on a nice function is not convergence on every function. Chebyshev interpolation converges for functions with a modest amount of smoothness, and there are continuous functions for which no fixed scheme of nodes converges — a theorem of Faber, from 1914. What is guaranteed is that some polynomial of each degree comes close, by Weierstrass’ approximation theorem; what Faber’s result denies is that interpolation at any pre-chosen set of points will find it.

That gap between “a good polynomial exists” and “this procedure finds it” is the whole subject in one sentence, and it is the same gap the least-squares fit closes in its own setting by asking for the nearest point rather than for an exact match.

What it costs, and the form nobody should use #

The polynomial through $n$ points can be written down three ways, and they differ enormously in what a machine makes of them.

Solve for the coefficients. Write the polynomial as $c_0 + c_1x + \cdots$, demand that it pass through every point, and solve the resulting square system. The matrix is the Vandermonde matrix — its columns are the powers of the nodes — and its columns are nearly the same vector for exactly the reason the monomials are a terrible basis: on an interval, $x^{10}$ and $x^{12}$ are almost indistinguishable. The system is solvable and the solution loses most of its digits, and the number of digits lost grows with $n$ faster than the accuracy being sought.

Use the Lagrange basis directly. Each basis polynomial is a product of $n-1$ factors, evaluating one costs $n$ multiplications, and evaluating the interpolant costs $n^2$. It never forms the Vandermonde matrix and is stable, and the cost per evaluation is quadratic, which for a curve drawn at eight hundred points is eight hundred times $n^2$.

Use the barycentric form, which is what every figure here does. Precompute one weight per node, at a cost of $n^2$ once, and every subsequent evaluation costs $2n$ operations — a division and two running sums. It is exact at the nodes by construction, so the assertion that the curve passes through every dot is a check on the implementation rather than on the arithmetic, and it is stable for the Chebyshev nodes in a sense that has been proved rather than observed.

The three compute the same polynomial and the first should never be written. That is worth stating plainly because it is the form a reader derives first, being the one that follows directly from the definition, and it is the form that makes high-degree interpolation look numerically hopeless when the hopelessness is entirely an artefact of the representation. The object is not ill-conditioned; one way of writing it down is, and the same distinction runs through the normal equations one field away.

Runge’s paper, and what it was about #

Runge was not investigating interpolation for its own sake. He was working on numerical methods for differential equations — the family of methods that carries his name comes from the same period — and the function $1/(1+x^2)$ appears in his 1901 paper as an example, chosen because its derivatives are easy to write down and grow.

The paper’s result was read as a warning against high-degree polynomials, and that reading held for most of a century: the standard advice in numerical analysis texts was to keep interpolation degrees low, and splines were developed partly in response. The rehabilitation is recent, and the argument is the one made above — the trouble was the node placement, the fix has been known since Chebyshev’s work in the 1850s, and the two facts took a long time to be stated together.

There is a small lesson in that delay. A counterexample is evidence about the hypothesis it violates, and reading it as evidence against the whole family is a mistake that can persist for a century. Runge’s function shows that evenly spaced interpolation diverges; it says nothing about polynomials, and it was read as though it did.

What the pictures cannot show #

Each figure draws one polynomial at one sample count, and the claim is about a sequence of them. The divergence is a statement about what happens as $n$ grows, and no drawing of finitely many curves establishes a limit; the numbers quoted in the prose are seven terms of a sequence, which is evidence and not proof.

The amplification figure plots a maximum taken over 2001 sampled points of the interval, which is a lower bound on the true worst case. For these node sets the function being maximised is smooth and the sampling cannot miss much, and the figure would still be drawing an estimate if it could.

And the interpolants are drawn on a clipped vertical range. At forty-one evenly spaced nodes the polynomial leaves the frame entirely near the ends, and a figure that showed the whole of it would be a figure of two vertical strokes with a function invisible between them — which is honest and useless, and is why the frame is cut instead.

The ladder from here #

Rungs above: minimax approximation, where the polynomial is chosen to make the worst error as small as possible rather than to pass through anything, and the equioscillation theorem that characterises it. Chebyshev series, which expand a function in the polynomials these nodes belong to and behave like a Fourier series in disguise — the projection from the semicircle is exactly the change of variable that turns one into the other. Splines, which give up a single high degree in exchange for low degree on many pieces and are what actually gets used on measured data. Barycentric interpolation, which is how the polynomials on this page were evaluated and is stable where the textbook formula is not. And Faber’s theorem, which says no fixed scheme of nodes works for every continuous function.

The rule that was not the rule #

The habit is about how a counterexample gets read.

The natural summary of the hero figure is high-degree polynomial interpolation does not work. That summary is memorable, it is what most readers take away, and it is false — the same degree on the same function with the points moved is accurate to seven figures. The true summary is equally spaced interpolation does not work, and the difference between the two is a hypothesis nobody noticed was doing any work, because evenly spaced points are what a person draws when asked for some points.

The question to ask of any counterexample is which of its features is load-bearing, and the way to find out is to vary them one at a time. Here there were three candidates — the degree, the function, and the node placement — and varying each in turn settles it in three pictures: raising the degree makes it worse, changing the function to a sine makes it go away, and moving the nodes makes it go away while the degree and the function stay put.