The best nodes have no formula

The points that ruin the fit showed that fitting a polynomial through equally spaced samples of a perfectly smooth function can go badly wrong, and that moving the sample points to Chebyshev’s positions — crowded towards the ends of the interval — repairs it. It measured the damage with one number, the Lebesgue constant, and an ellipse, not a disc ended by observing that Chebyshev’s points, although nearly the best, are not exactly the best. At every degree there are points with a slightly smaller constant.

This essay finds those points. They are characterised by a theorem from 1978, they are easy to compute, and nobody has a formula for them. The essay draws them, measures how much better they are than Chebyshev’s, and measures how far they sit from the simple points they so nearly equal.

What the Lebesgue function measures #

Interpolation takes values at $n + 1$ points $x_0, \ldots, x_n$ in the interval from $-1$ to $1$ and returns the unique polynomial of degree $n$ through them. It is linear: the polynomial through the sum of two data sets is the sum of their polynomials. So it can be written as $p(x) = \sum_j y_j \ell_j(x)$, where $\ell_j$ is the polynomial that is 1 at $x_j$ and 0 at every other node — the cardinal function of the $j$-th node.

If every data value is wrong by at most $\varepsilon$, the polynomial is wrong at $x$ by at most $\varepsilon \sum_j |\ell_j(x)|$, and that sum is the Lebesgue function of the nodes. It is 1 at each node, where only one cardinal function is non-zero, and rises between nodes as the cardinal functions overshoot. Its largest value over the interval is the Lebesgue constant $\Lambda$, the worst factor by which interpolation at those nodes can magnify an error in the data.

The first figure shows the three shapes. For equally spaced nodes the bumps grow enormously towards the ends, reaching 29.9 at eleven nodes; this is the instability behind Runge’s phenomenon. For Chebyshev’s nodes — the projections onto the interval of equally spaced points on a semicircle — the bumps are all below 2.5. For the optimal nodes they are all exactly 2.05, and that is the whole characterisation.

Why the ends need more points #

Chebyshev’s nodes are crowded towards the ends of the interval, and the crowding is the whole point. Place $n + 1$ points equally around a semicircle and drop them straight down onto its diameter: they land close together near the ends and far apart in the middle, at density proportional to $1/\sqrt{1 - x^2}$. The same construction, with the points of a regular polygon in place of the semicircle, is the one the polygon an equation forces drew for roots of unity, and that density is the same arcsine law that half the time is the rarest answer found governing how long a fair game stays on one side — it is the distribution that equalises influence on an interval, in the sense that a unit charge spread with that density produces the same electrical potential at every point.

That is why it is the right density for interpolation. A polynomial’s size at a point is controlled by the product of its distances to all the nodes, and that product is evenly balanced across the interval exactly when the nodes are spread with the equilibrium density. Equally spaced nodes have too few points near the ends, so the products there are exponentially larger than in the middle, and the Lebesgue constant grows exponentially with the degree. Any family of nodes whose density approaches the arcsine law has Lebesgue constant growing more slowly than any exponential; to get the optimal logarithmic growth the nodes must match the law closely, not merely in the limit.

The three logarithmic families in the figures all have exactly that density, and they differ only at the ends of the interval, by amounts that shrink as the degree grows. The equal-bump condition is the final adjustment, deciding precisely how far each node should sit from where the arcsine law puts it, and that adjustment is what has no formula.

A small error, magnified #

The constant is not a technicality. Every measured value has an error, and every computed value has rounding error, so the data of an interpolation problem are always slightly wrong. The figure takes a smooth function, perturbs each of 21 samples by 0.002 — far too little to see — and fits the polynomial. Through Chebyshev’s nodes the result stays within 0.005 of the true function. Through equally spaced nodes the same tiny errors grow to more than 14 near the ends of the interval, because the Lebesgue constant of 21 equally spaced nodes is about eleven thousand.

In both cases the damage stays inside what the constant allows, which the computation checks: the deviation is at most $\Lambda \varepsilon$. The constant is a guarantee, and also, for equally spaced nodes, a warning that turns out to be nearly attained. That is why Chebyshev’s points are the default in numerical practice, and why the question of whether they can be improved is worth asking even though the improvement is small.

Four kinds of node, as the degree grows #

The constants for equally spaced nodes grow exponentially, like $2^n/(e\,n\log n)$: past a thousand at degree 16 and past ten thousand at degree 20. For the other three kinds the growth is logarithmic, and by a theorem of Paul Erdős and of Faber before him, no choice of nodes can do better than logarithmic growth: every set of $n + 1$ nodes has a Lebesgue constant of at least about $(2/\pi)\log n$. Chebyshev’s nodes achieve that rate. Their constant is $(2/\pi)(\log(n+1) + \gamma + \log(8/\pi))$ to within a few hundredths — where $\gamma$ is Euler’s constant — and the computed values follow the formula at every degree tested.

The three logarithmic curves differ only in their constant terms. Chebyshev’s nodes leave the ends of the interval uncovered: the outermost nodes sit a little inside $\pm 1$, and the Lebesgue function rises in the uncovered end pieces. Stretching the nodes by a fixed factor so that the outermost ones land exactly on $\pm 1$ removes that excess at a stroke. The stretched constant is $(2/\pi)\log 2 \approx 0.44$ lower than Chebyshev’s for large degrees, and it is only about two hundredths above the optimal value at every degree computed, from 4 to 40.

Every bump the same height #

The optimal nodes are found from a characterisation that was conjectured by Sergei Bernstein and Paul Erdős in the 1930s and proved by Thomas Kilgore, and independently by Carl de Boor and Allan Pinkus, in 1978: among all sets of nodes containing the endpoints, the ones with the smallest Lebesgue constant are those for which the local maxima of the Lebesgue function between consecutive nodes are all equal.

The figure shows the characterisation at work for thirteen nodes. Chebyshev’s bumps are tallest at the ends, by almost half a unit. The stretched nodes have nearly equal bumps — within two hundredths — but the middle ones are slightly lower than the ones near the ends. The optimal nodes adjust every gap slightly so that every bump is exactly 2.1571. The computation finds them by that rule: starting from the stretched nodes, it shortens the gaps whose bumps are too tall and lengthens the ones whose bumps are too short, and repeats until the heights agree to ten decimal places, which takes a few dozen rounds.

Why equal bumps should be optimal can be seen by trying to improve on them. Moving one interior node slightly changes every bump, but not equally: the two bumps on either side of the moved node change in opposite directions, one rising and one falling, while the distant ones barely move. If the bumps are unequal, there is a direction of movement that lowers the tallest bump at the cost of raising shorter ones, and the constant — which is the tallest bump — goes down. Only when every bump is the same height is there no such move, because lowering any one of them raises a neighbour that is already as tall as it. The hard part of the theorem is showing that this local condition has exactly one solution, so that the equal-bump nodes are the global optimum and not merely a place where small moves stop helping.

The principle is an old friend. The error that keeps coming back to its worst showed that the best approximation to a function has an error curve whose extreme values are all equal in size, alternating in sign — Chebyshev’s equioscillation theorem. The optimal nodes satisfy the same kind of condition one level up: it is not the error of an approximation that equioscillates, but the worst-case magnification of errors, as a function of where the nodes are put.

A condition number for a problem, not a matrix #

The Lebesgue constant is the condition number of interpolation: the factor by which the problem itself, before any method is chosen, can amplify an error in its input. What a map does to a circle drew the corresponding number for a matrix — the ratio of its largest stretch to its smallest — and the two play the same role. A problem with a large condition number cannot be solved accurately from slightly wrong data by any method, however careful, because the answer genuinely depends sensitively on the data.

Writing the interpolating polynomial in different ways does not change its condition number, only how much extra error a method adds on top. The same map in a better basis made that point about matrices: the right basis can make a computation stable, but it cannot make an ill-conditioned problem well-conditioned. For interpolation, the barycentric formula used in these figures is the stable way to evaluate the polynomial, and with it the computed values are as accurate as the Lebesgue constant allows. Through equally spaced nodes no formula rescues the computation; through Chebyshev’s or the optimal nodes, the problem is well-conditioned at every degree anyone uses.

So choosing the nodes is choosing the problem. That is why the optimum matters in principle even though its practical advantage over Chebyshev’s nodes is a fraction of a per cent: it is the best-conditioned version of interpolation that exists, and its constant is the floor for every interpolation scheme with that many points.

How far the best nodes sit from Chebyshev’s #

The optimal nodes are very close to the stretched Chebyshev nodes. At degree 10 no node moves more than 0.0006; at degree 30 none moves more than 0.00006. The shifts are symmetric about the middle and vary smoothly along the interval: nodes in the outer parts move slightly outwards, those near the middle hardly move at all. The pattern is so regular that a formula seems to be just out of sight.

None is known. The stretched Chebyshev nodes have one of the simplest descriptions in numerical analysis — cosines of equally spaced angles, multiplied by a constant, the zeros of the Chebyshev polynomial $T_{n+1}$ that a solvable chaos of every degree iterated as a map — and the optimal nodes, a ten-thousandth away, have none. They are defined implicitly by the equal-bump condition, a system of $n$ equations in $n$ unknowns, and they can be computed to any accuracy, but no expression in terms of familiar functions produces them, and it is not known whether one exists. That is the question the previous essay left, and it remains exactly where it was.

How much the optimum is worth #

In practice the optimal nodes are almost never used. The improvement over stretched Chebyshev nodes is about two hundredths in the Lebesgue constant, which means that the worst-case magnification of an error is about one per cent smaller. Against that, Chebyshev’s nodes have a formula, are the same for every function, and connect interpolation to the cosine series of an ellipse, not a disc, which allows the interpolating polynomial to be computed by a fast transform instead of solving equations. The optimal nodes have none of those advantages.

What the optimal nodes do settle is the limit of the subject. They show that the constant $(2/\pi)\log n$ is not merely the growth rate of good nodes but the best possible growth, with a known constant term: Péter Vértesi established that the optimal Lebesgue constant is $(2/\pi)(\log n + \gamma + \log(4/\pi)) + o(1)$, and the computed optimal values approach that formula from above as the degree grows. Chebyshev’s nodes lose $(2/\pi)\log 2$ to the optimum because of their uncovered ends, the stretched nodes lose almost nothing, and no other choice can do better.

Nodes in higher dimensions #

The question changes character on a square or a triangle instead of an interval. Interpolation by polynomials in two variables requires the nodes to be placed so that no polynomial of the given degree vanishes on all of them — a condition that equally spaced grids often fail — and the Lebesgue constant then depends on a two-dimensional arrangement. Good node sets for the square are known, the Padua points among them, with Lebesgue constants growing like the square of the logarithm of the degree. For the triangle, the best known sets come from numerical optimisation, and their constants are known only by computation.

The pattern of the interval repeats there: simple explicit points that are nearly optimal, and optimal points that can be computed but not described. In two dimensions even the rate of the optimal constant is not settled in general — whether some arrangement achieves logarithmic growth, as on the interval, or whether the square of the logarithm is forced — and the one-dimensional characterisation by equal bumps has no proved analogue.

What the computation does not show #

The figures compute Lebesgue functions from barycentric formulas, locate each bump by sampling and refining, and find the optimal nodes by equalising the bumps until they agree to ten decimal places. That they are optimal rests on the theorem of Kilgore and de Boor–Pinkus, not on the computation, which only finds nodes satisfying the theorem’s condition. The degrees computed run to 40, and the observation that the stretched nodes stay about two hundredths above the optimum is a measurement over that range; the asymptotic behaviour of the gap is a separate question.

The noise figure uses one function and one pattern of errors. The Lebesgue constant bounds the damage for every function and every pattern; a particular pattern may cause less, and the figure shows the bound being approached for equally spaced nodes rather than proving it is attained.

There is also the question of uniqueness, which the computation takes for granted. The equalising procedure starts from the stretched Chebyshev nodes and moves towards equal bumps, and it converges every time it is run; that it converges to the same nodes from any reasonable start, and that those are the only nodes with equal bumps, is part of what the 1978 theorem proves. A computation that found two different sets of nodes with equal bumps would contradict the theorem, and none of the runs here does — but the runs start from one place each, so they test the theorem only lightly.

Still open: a description of the best nodes #

No formula is known for the nodes that minimise the Lebesgue constant at a given degree. They are characterised completely, unique up to the choice of interval, computable to any accuracy, and a ten-thousandth away from points with one of the simplest formulas in the subject — and that combination is what makes the absence of a description a question rather than a curiosity. It is not even known whether the shift from the stretched Chebyshev nodes, scaled by the degree, converges to a definite function as the degree grows, which would be a first step towards describing it.

The corresponding questions in two and more dimensions are wider open: the optimal growth rate of the Lebesgue constant on the square and the triangle is not known, nor is an analogue of the equal-bump characterisation that would identify the optimal nodes when they are found.