A basis that has to know where the roots are

The roots are not the coefficients took the oldest demonstration in the subject and measured it. Expand $\prod_{i=1}^{n}(x - i)$ exactly, round its coefficients to doubles, and compute the roots as the eigenvalues of the companion matrix: at degree twenty a root whose true value is an integer comes back wrong in the third significant figure, and the root condition number times the unit roundoff predicts it. The computation is backward stable. The polynomial it answers for is a nearby one, and a nearby polynomial has quite different roots.

The essay closed with three ways out, priced. One was a sentence long and was not measured: “A polynomial expressed in Chebyshev coefficients on an interval containing its roots, with the corresponding colleague matrix, has a conditioning that does not degrade exponentially with degree. The routine is the same — build a matrix, factorise it — and the representation is the thing that changed.”

The routine and the representation are both easy to build, so the sentence can be measured on its own polynomial. The change of basis is worth ten orders at degree twenty. It is not what the sentence said it was.

The same polynomial, in another basis #

To make the comparison fair the Chebyshev coefficients have to be exact before they are rounded, exactly as the monomial ones were. On $[1, n]$ write $t = (2x - n - 1)/(n - 1)$, so each factor $x - i$ becomes a multiple of $t - t_i$, and multiply the factors in one at a time in the Chebyshev basis using $2t\,T_k = T_{k+1} + T_{k-1}$. In integer arithmetic over a common denominator that is exact at any degree. The coefficients are then rounded to doubles once, and the roots are the eigenvalues of the colleague matrix — the Chebyshev basis’s companion matrix, tridiagonal with halves on the off-diagonals and the coefficients in its last row, which two approximants and one matrix size wrote out as the comrade matrix’s simplest case — computed by the same Francis iteration the monomial route uses, reduced to the form the form a real matrix can reach describes, with its convergence checked.

Beside the integers, a control: a polynomial of the same degree whose roots are placed at the Chebyshev points of $[1, n]$, the clustering towards the ends that the Chebyshev basis is built around. Both polynomials are handed to both routes.

The figure at the top of the page is the measurement: worst relative root error against degree, for both routes on both root sets. Three of its four lines climb and one does not, and the four lines say three different things.

The monomial route loses the integers’ roots as the earlier essay found — $1.7\cdot10^{-13}$ at degree six, $7.6\cdot10^{-3}$ at twenty, everything by forty — and it loses the Chebyshev points’ roots just as fast: $2.5\cdot10^{-13}$ at six, $4.2\cdot10^{-3}$ at twenty. The monomial basis does not care where the roots are, and its own condition numbers say so before any eigenvalue is computed: the componentwise condition number of the worst root, times $u$, is $3.5\cdot10^{-13}$ for the integers and $3.3\cdot10^{-13}$ for the Chebyshev points at degree six, $2.2\cdot10^{-7}$ and $1.2\cdot10^{-7}$ at fourteen, $6.0\cdot10^{-3}$ and $2.7\cdot10^{-3}$ at twenty — within a factor of 2.2 at every degree measured. On $[1, n]$ the functions $x^k$ of high degree are all small except near the right end, where they are all large and nearly proportional to each other, so a coefficient vector describes the polynomial well there and poorly everywhere else, whichever roots it has.

The Chebyshev route on the Chebyshev points does not degrade at all. Its worst error is between $4\cdot10^{-15}$ and $9\cdot10^{-14}$ at every degree from six to forty, which is the eigensolver’s own rounding.

And the Chebyshev route on the integers is the one the sentence was about. At degree twenty its worst error is $4.4\cdot10^{-13}$ — seventeen billion times better than the monomial route’s, ten orders. But it climbs: $1.3\cdot10^{-15}$ at six, $5.8\cdot10^{-11}$ at twenty-four, $6.8\cdot10^{-8}$ at thirty-two, $2.8\cdot10^{-4}$ at forty. Fitted from twelve to thirty-two it loses 0.35 decades a degree, against the monomial route’s 0.76 from six to twenty. That is a slower exponential, not the absence of one.

The middle of an equispaced polynomial is small #

Why the integers and not the Chebyshev points is visible in the polynomials themselves, before any basis is chosen.

With roots at the integers the polynomial’s humps between neighbouring roots are of very different heights: the ones near the ends of the interval reach its largest value, and the ones in the middle are $3.8\cdot10^{-5}$ of it at degree twenty. The polynomial is tall at the edges and nearly flat in the middle. With roots at the Chebyshev points every hump reaches the same height — the polynomial equioscillates, which is what Chebyshev points are for.

A root’s sensitivity to the coefficients, in any basis whose functions stay between $-1$ and $1$ on the interval, is roughly the polynomial’s size divided by its slope at the root. For an integer root in the middle of the interval the slope is set by the small humps around it and the size by the large ones at the ends, and the ratio of the two is the hump heights’ range: $10^{4.4}$ at degree twenty, and growing by a fixed factor with every degree, because each extra root deepens the trough in the middle. Drag the dial: at degree forty the middle humps are below $10^{-9}$ of the end ones. It is the same shape that makes interpolation at equispaced points fail — the Runge phenomenon seen from the roots’ side instead of the interpolant’s.

Each root’s error confirms it. At degree thirty the Chebyshev route keeps the outermost roots to between $10^{-13}$ and $10^{-15}$ and loses the middle ones to $3\cdot10^{-8}$ — five to seven orders between the ends and the centre. The monomial route at degree sixteen has the same profile, skewed toward the large roots because its basis functions $x^k$ are largest there. In both bases the hard roots are the middle ones, and the reason is the polynomial, not the basis.

Two factors, both exponential #

The size of the Chebyshev route’s error can be accounted for, and it takes two numbers rather than one.

The first is the roots’ sensitivity to the coefficients measured the way the colleague route perturbs them — normwise, relative to the whole coefficient vector, $\|c\|_1 / |p'(r)|$ — times the unit roundoff. That is the conditioning, and by the argument above it grows exponentially for equispaced roots: $7\cdot10^{-17}$ at degree six, $1.3\cdot10^{-13}$ at twenty, $4\cdot10^{-8}$ at forty.

The second belongs to the route rather than the problem. The colleague matrix’s last row is the coefficients divided by twice the leading one, $c_k/(2c_n)$, so a backward error of one rounding in that matrix corresponds to a perturbation of the coefficients that is larger by $\|c\|/|c_n|$.

For the integers that ratio — the largest coefficient over the leading one — is 1 at degree six, 46 at twenty, 1,330 at thirty and 44,000 at forty: the polynomial’s Chebyshev series has most of its weight in its low terms and a small leading coefficient, because the equispaced polynomial is far from the Chebyshev polynomial of its degree. For the Chebyshev points the ratio is one at every degree — the polynomial is, to within the grid its roots were rounded to, exactly $2^{1-n}T_n$.

Multiplied, the two factors bound the measured error from degree fourteen on, with the error between 0.08 and 0.46 of the product at every degree measured — a factor-of-six band across a quantity that changes by eleven orders. Below fourteen both are under $10^{-14}$ and the error is the eigensolver’s rounding, which neither factor models. So the Chebyshev route on the integers loses exponentially for two reasons that compound: the roots are exponentially sensitive to the coefficients, and the colleague matrix amplifies its own rounding by an exponentially growing ratio before handing it to them.

How much is the description and how much the solve #

A computed root can be wrong for two reasons, and the earlier essay separated them for the monomial route at degree twenty: the coefficients handed over were already rounded, so the polynomial being solved was not quite the one written down; and the eigenvalue computation then made its own error, exact for a nearby companion matrix. With the exact coefficients in hand both can be measured. The rounding error of every coefficient is computable exactly in integer arithmetic, and what it does to each root, to first order, is that error weighted by the basis function at the root over the polynomial’s slope there — the exact answer to a nearby problem with the nearby problem written out coefficient by coefficient.

The two routes lose their digits in different places. In the monomial basis with roots at Chebyshev points, where no coefficient is an integer and every one is rounded, the rounding alone accounts for the measured error to within a factor of five at every degree from ten to twenty: $10^{-3}$ of the $4.2\cdot10^{-3}$ at twenty. The monomial description is the damage; the eigensolver adds little. On the integers the monomial coefficients are exact integers up to degree eighteen — below $2^{53}$, the boundary what a float can hold draws — so there the description is perfect and every lost digit is the solve’s.

In the Chebyshev basis on the integers it is the reverse of the first case. The rounded Chebyshev coefficients are an almost perfect description of the roots: they move them by $5.9\cdot10^{-15}$ at degree twenty and $1.9\cdot10^{-9}$ at forty, where the measured errors are $4.4\cdot10^{-13}$ and $2.8\cdot10^{-4}$. From degree eighteen on the rounding is less than a thirtieth of the error, and from twenty-eight less than a thousandth. Nearly all of what the Chebyshev route loses, it loses in the colleague matrix — which is the second factor above, the last row’s $c_k/(2c_n)$, seen from the other side. The change of basis fixed the description almost completely, and left the solver with a matrix it handles badly.

The rate follows the roots #

If the basis has to match the roots, then the rate should change continuously as the roots move from one distribution to the other.

Move each root a fraction $\alpha$ of the way from its integer to the corresponding Chebyshev point, and fit the Chebyshev route’s loss from degree twelve to thirty-two. At $\alpha = 0$ it is 0.35 decades a degree; at a quarter of the way, 0.23; at half-way, 0.08; from three quarters on, 0.05 to 0.06, which is the rounding floor drifting with the size of the matrix rather than any loss of conditioning. The monomial route at degree twenty is between $4\cdot10^{-3}$ and $10^{-2}$ at every one of these mixes.

So the two bases answer the question about where the roots are in opposite ways. The monomial basis is indifferent: its basis functions are concentrated at the end of the interval far from the origin, and any polynomial of high degree is badly represented in it. The Chebyshev basis is all about it: it represents a polynomial well exactly when that polynomial looks like a Chebyshev polynomial — equioscillating, with its roots clustered at the ends — and it degrades smoothly as the roots spread out from that pattern toward equal spacing.

What the earlier sentence should have said #

The earlier essay’s sentence is right about the size of the improvement and wrong about its kind. At degree twenty the Chebyshev route is ten orders better on Wilkinson’s own polynomial, which is far more than any other choice in that essay buys. But the conditioning in the Chebyshev basis does degrade exponentially for these roots — slower than in the monomial basis, by a factor of about two in the exponent — and the claim that it does not is true only for roots spread like the Chebyshev points themselves.

The same clustering turns up where nobody designed it. The points the algorithm chose found an adaptive rational approximant placing its support points in a geometric cluster at a branch point — the distribution its own residual asked for, as the Chebyshev points are the distribution a polynomial on an interval asks for. A representation is good where its points or its roots sit the way the representation expects them.

The general statement is the one two approximants and one matrix size reached for rational and polynomial approximants of a matrix function, and a basis built from the points reached for polynomial fitting: a basis is good for a problem when it is adapted to where the problem’s information sits. For least-squares fitting at data points, Arnoldi’s recurrence builds a basis orthonormal on the points themselves. For roots, the analogue would be a basis adapted to where the roots are — which is a basis that needs the answer to be chosen.

That is not a reason to keep the monomial basis. It is a reason to say what the repair buys: for roots that are not too far from the ends of an interval, everything; for equispaced roots, a slower exponential, ten orders at degree twenty and three and a half at forty. And it is a reason to keep the condition number is an amplifier in view, since what the change of basis reduced was the condition number — not to one, but by as much as the polynomial’s resemblance to a Chebyshev polynomial allowed.

What forty degrees do not show #

Every polynomial here is monic with real, simple roots, on an interval chosen to contain them exactly, and every root set is symmetric about the interval’s centre, so half the Chebyshev coefficients are zero. A polynomial whose roots lie outside the interval, or are complex, is represented less well in the Chebyshev basis on that interval, and nothing here measures by how much. The Chebyshev points are rounded to a grid of 1/1024 so that the exact expansion stays rational; that rounding is what makes their leading coefficient differ from $2^{1-n}$ in the third digit, and it is far below anything the root errors depend on. The colleague matrix is solved by the same nonsymmetric Francis iteration as the companion matrix; a solver that exploited its tridiagonal-plus-rank-one structure would be cheaper, and whether it would be as accurate on these polynomials is not measured here.

Still open: balancing the colleague matrix, an adapted basis, and exact bisection priced #

Balancing. The route’s second factor, $\|c\|/|c_n|$, is the size of the colleague matrix’s last row relative to the rest of it. Diagonal balancing of the companion matrix is a standard step before its eigenvalues are computed. The prediction with a sign is that balancing the colleague matrix removes most of the second factor on the integers — bringing the error to within a factor of ten of the normwise conditioning alone at every degree from fourteen to forty — which would leave the conditioning as the only exponential.

A basis adapted to equispaced roots. Orthogonal polynomials on the equispaced points $1, \dots, n$ — the discrete Chebyshev, or Gram, polynomials — are to equispaced roots what Chebyshev polynomials are to Chebyshev points. The prediction is that Wilkinson’s polynomial in the Gram basis, solved by the corresponding comrade matrix, keeps its roots to within a hundred times the unit roundoff to degree forty, and that the Chebyshev-point polynomial in the same basis degrades — the two distributions trading places.

Exact bisection, priced. The earlier essay’s third way out, exact evaluation and bisection, has no rounding at all and was called “entirely practical up to degree fifty”. The prediction is that its cost in BigInt operations to isolate every root of Wilkinson’s polynomial to a relative $10^{-15}$ grows like $n^3 \log n$, and that it overtakes the colleague route in accuracy-per-millisecond from degree twenty-four on.