A basis built from the points

The valley with no bottom fitted a polynomial of degree eleven twice — once in monomials, whose design matrix had a condition number ten orders of magnitude larger, and once in Chebyshev polynomials — and drew the two fitted curves lying on top of each other. It concluded that the basis decides the coefficients and not the curve: the data did not determine the coefficients along the valley floor, the arithmetic had nothing to do with it, and the curve every reader cares about is the same.

That conclusion is true of the data the essay fitted, which was a smooth function sampled with no noise, so the least-squares residual was at the level of rounding. It turns out to be a statement about the residual rather than about the basis. With noise in the data the curves separate, the separation grows with the condition number of the basis, and at high degree the monomial curve is wrong in the third decimal place while the orthogonal one is not.

That makes the choice of basis a numerical decision after all, and it reopens the question of which orthogonal basis. Chebyshev polynomials are the textbook answer and they are the right answer on the interval. A least-squares fit does not see an interval. It sees a set of points, and a basis that is orthogonal on the points is a different object that has to be built from them.

The curve on exact data, reproduced #

The first step is to reproduce the valley essay’s observation, so that everything after it is a change of one input rather than a change of experiment. Two hundred equispaced points on [−1, 1], a smooth function sampled exactly, and a fit of degree 40. The monomial design matrix has a condition number of 7.4·10¹⁴, within a factor of twelve of what a double can represent at all.

The monomial curve and the curve from an orthonormal basis agree to 2.7·10⁻¹⁴ at the sample points, and both are within 2·10⁻¹⁵ of the function in root-mean-square. So the observation holds at forty, as it did at eleven: a condition number of 10¹⁴ in the basis costs the fitted curve nothing, on this data.

The reason that is possible is the one the projection and the right angle is built on. The fitted values Ax are the orthogonal projection of the data onto the column space of A, and the column space is the space of polynomials of the chosen degree whichever basis spans it. A backward-stable least-squares routine — Householder QR — returns the exact projection onto the column space of a slightly perturbed matrix A + δA, with ‖δA‖ a modest multiple of u‖A‖. Perturbing the matrix tilts the space, and tilting a space moves the projection of a point by an amount proportional to how far the point is from the space.

That distance is the residual. On exact data the residual is at rounding, the point is in the space to within rounding, and no tilt of the space moves its projection by anything visible. The condition number of the basis decides how much a given ‖δA‖ tilts the space, and it is multiplied by a residual of 10⁻¹⁵.

Add noise, and the curves come apart #

The perturbation theory of least squares puts that argument into a formula: the error in the fitted values from a backward-stable solve is of order κ(A) · u · ‖r‖, on top of the rounding every route pays. With ‖r‖ at rounding the product is invisible at any κ a double can hold. With ‖r‖ at the noise level it is not.

The same function, the same 200 points, and noise of 0.1% added to the data. At degree 40 the monomial curve and the orthogonal-basis curve now differ by 1.7·10⁻⁵ at the sample points, where on exact data they agreed to 2.7·10⁻¹⁴ — a factor of 650 million for a change in nothing but the residual. At degree 48, where the monomial matrix is beyond what a double can represent, the two curves differ by 0.004, in the third decimal place of a fit to data of size one.

The prediction holds to within the looseness of its constant. Measured against κ·u·‖r‖ over every degree where the monomial condition number is representable and the prediction is above the rounding floor, at two noise levels and on two point sets, the ratio of measured to predicted separation lies between 0.097 and 2.9 — inside a factor of ten in both directions at every one of twenty measurements. The factor of ten is required in both directions, because a separation a hundred times smaller than predicted would mean the mechanism was not the one argued.

The proportionality to the residual is the part to check directly, and it is clean. At 10⁻⁶ noise instead of 10⁻³, the monomial separation at degree 40 is 1.66·10⁻⁸ against 1.73·10⁻⁵: a thousandfold less noise, a factor of 1,042 less separation. At degree 36 it is 5.8·10⁻¹⁰ against 6.0·10⁻⁷, a factor of 1,030. The basis contributes the condition number; the data contributes the residual; the product is the error.

At 10% noise, which a scatter plot from a laboratory can easily carry, the monomial fit of degree 48 is off by 0.4. That is not a subtle numerical effect on a plot of the data: it is a different curve. The fit is also a poor model at that degree — forty-eight coefficients chasing 10% noise — and a practitioner would not fit it. The point is not that degree 48 is wise. It is that the claim “the basis changes the coefficients and not the curve” has a hidden condition, and the condition is a small residual, which is the one thing noisy data does not have.

A small residual is not a small error is the collection’s standing warning in one direction. This is the other direction: a large residual — honestly large, because the data is noisy — is what lets a large condition number reach the answer.

Degree by degree, the separation follows the basis #

The measurement behind that factor of ten is worth reading row by row, because the rows show the mechanism working rather than merely agreeing with it. On 200 equispaced points with 0.1% noise, the monomial curve’s separation from the orthogonal-basis curve is 5.1·10⁻¹² at degree 20, 2.8·10⁻¹¹ at 24, 2.0·10⁻¹⁰ at 28, 1.2·10⁻⁸ at 32, 6.0·10⁻⁷ at 36 and 1.7·10⁻⁵ at 40. It rises by between five and sixty times for every four degrees, and so does the monomial condition number that drives it; the prediction κ·u·‖r‖ over the same degrees is 1.7·10⁻¹², 5.7·10⁻¹¹, 1.8·10⁻⁹, 6.1·10⁻⁸, 2.1·10⁻⁶ and 7.0·10⁻⁵. Only the first row sits above its prediction, by 2.9 times, at a separation small enough that the rounding every route pays is part of it. From degree 24 on, the measured separation is a steady fifth to a half of the prediction.

A constant fraction below a bound is what a bound that is attained up to a constant looks like, and it is the reason the check runs in both directions. If the monomial curve were drifting from the orthogonal one for some other reason — the fit’s own ill-posedness at high degree, which the valley with no bottom is about, or a defect in the reference — the ratio would wander with the degree instead of holding. On the gap set it holds too, between 0.10 and 1.04 over degrees 16 to 36.

The monomial basis is bad on every point set, and the point distribution barely helps it. On Chebyshev points, the set most favourable to polynomial approximation, the monomial matrix still reaches 8.3·10¹⁶ at degree 48, against 1.5·10¹⁷ on equispaced points and 1.2·10¹⁷ on the clustered set. A basis whose columns are powers of x is ill-conditioned because the powers look alike, and no arrangement of the points makes x⁴⁰ and x⁴² look different enough.

On the clustered set the Chebyshev curve’s separation at degree 48 is 2.9·10⁻¹³, a hundred times its separation on equispaced points, which is the Chebyshev condition number of 1.03·10⁶ doing to the curve what the monomial one does, a hundred thousand times more gently. That is the whole difference between the two classical bases on data that is not spread like either’s weight: one fails by four digits and the other by eleven, and the gap closes as the degree rises.

What building the basis costs #

A basis made from the points is not free, and the price is easy to state. Each degree adds one column: a multiplication by the points, m operations for m points, and an orthogonalisation against every column built so far, about m times the degree, done twice. Over d degrees that is about m·d² operations — the same order as the Householder factorisation of an m × (d + 1) design matrix, which the fit needs whichever basis is used. The Arnoldi basis therefore costs a constant multiple of the factorisation it makes well conditioned, and after it is built the least-squares solve is a matrix–vector product, because its columns are orthonormal and the coefficients are simply their inner products with the data.

The same arithmetic is what the right-hand side as one more column measured for Gram–Schmidt as a least-squares method: orthogonalising against the columns of the problem, rather than forming a product with a computed Q, is what keeps the answer’s accuracy. The coefficients here are inner products with an orthonormal basis that is orthonormal on the data to 10⁻¹⁰, so that concern does not arise; it would arise for a basis orthogonalised once against columns ill-conditioned enough to lose orthogonality, which is why the process is run twice.

Chebyshev polynomials are orthogonal on the interval, not on the data #

So the basis matters on noisy data, and an orthogonal basis is the repair. Chebyshev polynomials are the standard orthogonal basis on [−1, 1], and on equispaced points their design matrix is well conditioned: 3.04 at degree 8, 3.52 at degree 24, 17.7 at degree 48. On 200 Chebyshev points — clustered towards the ends like the Chebyshev weight — it is 1.41 at every degree.

That last number is the clue to what the basis is orthogonal with respect to. Chebyshev polynomials are orthogonal under the weight 1/√(1 − x²), and a sample whose points are distributed like that weight makes their discrete inner products nearly the continuous ones. A sample distributed differently does not.

On 200 points clustered towards one end of the interval, Chebyshev’s condition number is 5.57 at degree 8, 13.8 at degree 24 and 1.03·10⁶ at degree 48. On two intervals with a gap in the middle it is 9.52 at degree 8, 2,490 at degree 24 and 1.55·10⁷ at degree 48. Its conditioning grows exponentially in the degree on both sets, like everybody else’s, only later.

The gap set is not exotic. Measurements taken on two separate occasions, a sensor that fails for a stretch, a response measured at low and high doses and not in between: data with holes is ordinary, and on it the interval’s orthogonal basis is an ill- conditioned basis for the data’s own problem. The fitted curve inherits the condition number exactly as the monomial one does. At degree 48 on the gap set with 0.1% noise, the Chebyshev curve is off by 3.5·10⁻¹² — small, because 1.55·10⁷ times u times the residual is small — and the monomial curve by 0.0037.

Building the basis from the points #

A basis orthogonal on the data has to be built from the data, and the construction is short. Start with the constant vector on the sample points, normalised. Multiply the latest basis vector, entry by entry, by the sample points x — which is what multiplying a polynomial by x does to its values. Orthogonalise the product against every basis vector built so far, in the ordinary inner product of the sample. Normalise, and repeat.

That is Arnoldi’s process applied to the diagonal matrix of sample points with the constant vector as the starting vector, and Brubeck, Nakatsukasa and Trefethen named it Vandermonde with Arnoldi when they proposed it for exactly this problem; the idea of orthogonal polynomials generated from the data goes back to Forsythe’s fitting routines of the 1950s. The columns span the same space as the monomials of the same degree — each new column is x times the previous one plus a combination of earlier ones — so the fit is the same least-squares problem. What changes is only the matrix it is computed from, and by construction that matrix has orthonormal columns on the data. Its condition number is one.

Measured, it is one on every set: below 1 + 10⁻¹⁰ on equispaced, clustered, gap and Chebyshev points at degrees 8, 24 and 48. The orthogonalisation is done twice, which is two Gram–Schmidts’s lesson applied: a single pass of modified Gram–Schmidt loses orthogonality in proportion to the conditioning of the columns it is handed, and here those columns are x times an orthonormal column, which is well conditioned, so once is nearly enough. The second pass is there so that κ = 1 is a measurement at every degree drawn rather than a hope about the ones not checked. It is the process an orthogonalisation nobody calls one describes inside every Krylov method, used here for a problem with no linear system in it at all.

Evaluating the fit anywhere else #

A basis built from the sample points has a practical objection: a fitted curve is wanted at points that were not sampled, and the basis vectors are only defined at the sample.

The objection is answered by the same recurrence that built them. Arnoldi’s process records, at each step, the coefficients it subtracted — the Hessenberg matrix of the process. Those coefficients are all that is needed to run the recurrence again at any other set of points: multiply by the new x, subtract the recorded combination, divide by the recorded norm. The polynomials that come out are the same polynomials, evaluated elsewhere, and the fitted coefficients in the Arnoldi basis apply to them unchanged.

That is checked as two routes to one curve. On the gap set at degree 24, the Arnoldi fit pushed through its Hessenberg matrix onto a fine grid between the data’s points is compared with the Chebyshev fit evaluated by cosines on the same grid — a degree where the Chebyshev matrix’s condition number is 2,490, small enough for its curve to be trusted. The two agree to 5.6·10⁻¹⁵ on the data’s own intervals.

The recurrence is stable where the sample is informative and not where it is not. Evaluating it far outside the data, or across a gap much wider than the data’s spacing, extrapolates a high-degree polynomial into a region where nothing constrains it, and no basis repairs that; it is a property of the fit, not of its arithmetic. The agreement above is measured on the data’s intervals and claimed nowhere else.

What the valley essay’s conclusion becomes #

The valley essay’s central sentence — the data never contained the digits the coefficients lack — survives intact, and this essay depends on it. What needs a clause is its corollary about the curve. On exact data the basis cannot reach the curve. On noisy data it reaches it by κ·u·‖r‖, and a basis that is orthogonal on some other set of points does not escape that; only a basis orthogonal on the sample does.

That changes the practical advice in one place. Choosing a basis by its behaviour on the interval — Chebyshev for polynomials, Legendre for integrals — is right when the sample is spread like the basis’s weight, and a condition number in the millions when it is not. Building the basis from the points costs one Gram–Schmidt sweep per degree, a matrix of recurrence coefficients to store, and nothing in accuracy, on every set measured here.

An earlier essay names this kind of discovery: the roots are not the coefficients found that a polynomial’s conditioning depends on how it is represented, and the representation that was natural for the mathematics was the one that was ill-conditioned for the arithmetic. The monomial basis is natural for writing a polynomial down. The basis that is natural for a least-squares fit is the one the points define.

When the matrix is wrong too is the case where the perturbation δA is not rounding but measurement error in the matrix itself, and there the same κ·‖r‖ term is a statement about the data rather than the arithmetic. That is a different problem with a different repair, and an orthogonal basis does not remove it: a matrix that is genuinely uncertain has an uncertain column space whatever basis spans it.

Where this goes from here #

Rational functions instead of polynomials. A polynomial basis is the wrong basis for a function with a nearby singularity, and the same authors’ rational fitting methods build their bases from the data in the same way. The points the algorithm chose measures one of them, AAA, as an interpolation method; its least-squares version and the conditioning of its barycentric basis on noisy data are the next measurement.

Several variables. Arnoldi’s recurrence extends to polynomials in two or three variables by multiplying by each coordinate in turn, and the number of basis vectors grows combinatorially. How the condition number of the multivariate monomial basis grows on scattered two-dimensional data, and whether the Arnoldi basis keeps it at one at a price that is still worth paying, is unmeasured.

And the degree. Every separation above is at a fixed degree. Choosing the degree is the regularisation question for polynomial fits — too low loses the function, too high fits the noise — and it is the question the regularisation field answers for deconvolution. Whether cross-validation and the discrepancy principle choose the degree as well as they choose a Tikhonov parameter, on a basis whose conditioning no longer interferes, is where fitting meets regularisation.