The roots are not the coefficients

A scalar polynomial is a matrix polynomial of size one. Its companion matrix is the linearisation, the eigenvalues of that matrix are the roots, and everything the polynomial field says about linearising applies with n = 1.

Which means the field’s central claim can be tested at a size anybody can check by hand — and the test is the oldest example in this subject, sixty years old and still the sharpest.

The polynomial #

Take ∏(x − i) for i = 1 to n and expand it. The roots are the integers 1 to n, which is as known as an answer gets — no closed form to evaluate, no reference computation, just twenty integers.

The coefficients are integers too, and enormous: at n = 20 the constant term is 20! = 2.43·10¹⁸ and the coefficient of x¹⁹ is −210. Expanded exactly in BigInt, as they are here, so the only approximation anywhere in the demonstration is writing them down as doubles.

Then hand them to the route every numerical library takes: build the companion matrix, factorise it, read off the eigenvalues.

What comes back #

At degree twenty a root whose true value is an integer comes back wrong in the third significant figure. Not the fifteenth, not the eighth — the third.

The second column is the prediction, and it is made from the coefficients without reference to the computation. For a simple root r, perturbing coefficient aₖ by a relative ε moves r by about ε|aₖ||r|^k/|p′®|, so the relative sensitivity of r is the largest of those over k, divided by |r|. It reaches 1.0·10¹³ at degree twenty, at the root fourteen — not at the largest root and not at the smallest, but in the middle, where |r|^k is large and the derivative of the product is not.

The prediction crosses from over-estimate to under-estimate between degree 10 and 14. Measured against predicted, the ratio runs 0.83, 7.6, 9.1 and 5.0 at degrees 10, 14, 18 and 22 — so κ·u is conservative at degree 10 and optimistic by nearly an order of magnitude at 18.

The coefficient column says where the crossing comes from. 2⁵³ is 9.007·10¹⁵, and the largest coefficient runs 1.28·10⁷, 3.92·10¹¹, 3.4·10¹⁶ and 6.76·10²¹ at those four degrees — so at degree 10 and 14 the polynomial handed to the routine is the exact one, and at 18 and 22 it is not. Past that point two errors are present rather than one: the conditioning amplifying the arithmetic, which κ·u predicts, and a perturbation committed before the routine was called, which it cannot. That is why the prediction is conservative below the threshold and optimistic above it, and it is a sharper account than a factor of nine.

Measured against predicted agrees to within a factor of nine at every degree and a factor of three at the last one. The error is the condition number times the unit roundoff, which is the condition number is an amplifier doing exactly what it says, on a problem where every other quantity is exactly known.

And the computation did nothing wrong #

The temptation is to blame the eigenvalue routine, and the accounting refuses it.

The companion matrix at degree twenty has a norm of 2.4·10¹⁸, dominated by its first row. A real Schur factorisation of it is backward stable in the ordinary sense: the computed eigenvalues are exact for a companion matrix a relative rounding away. That theorem holds here, and it is the same one the form a real matrix can reach rests on.

But a perturbation of the matrix is not a perturbation of the polynomial. Change a single entry of the companion matrix by one unit in the last place and what results is the companion matrix of a polynomial whose coefficient has moved by one ulp — and one ulp of 2.43·10¹⁸ is about 500. The roots of a polynomial whose constant term has moved by 500 are not the integers.

So the answer is exact for a nearby polynomial, and a nearby polynomial has quite different roots. Which is the exact answer to a nearby problem with the second half of the sentence doing the damage, and it is the same structure as a backward-stable answer to a problem nobody asked one field over — there the mismatch is between a matrix and three coefficient matrices, and here between a matrix and a list of coefficients.

The coefficients are not binary64 numbers #

There is a second failure underneath the first, and it arrives at a specific degree.

An integer is exactly representable in binary64 up to 2⁵³ ≈ 9.0·10¹⁵. The largest coefficient of the Wilkinson polynomial passes that at degree nineteen: 18! is 6.4·10¹⁵ and 19! is 1.2·10¹⁷.

So above degree eighteen the polynomial handed to the routine is not the polynomial that was written down. It has been rounded before any algorithm ran, by a relative 10⁻¹⁶ in a coefficient whose condition number is 10¹³, which accounts for three of the digits lost at degree twenty on its own.

That is worth separating from the conditioning argument, because it is a different kind of failure and has a different repair. Conditioning is a property of the question; representability is a property of the format, and it is the same boundary the units that overflow before the answer does measures for a matrix polynomial. Here it arrives at 2⁵³ rather than at 10³⁰⁸ because the quantity that has to be exact is an integer rather than merely finite.

What the site’s own machinery says about it #

This site has an exact-arithmetic library, and it is worth saying plainly what it can and cannot do here.

exact.js expands the polynomial in BigInt, so the coefficients are exact at any degree — that is how the reference above is built. It can also evaluate the polynomial exactly at any rational point, so a root can be bracketed exactly: a sign change between two rationals is a certificate, and bisection on exact evaluations converges to any accuracy wanted with no rounding anywhere.

What it cannot do is make the double-precision route accurate, because the double-precision route’s problem is that the data it is given is already wrong. And what nobody can do is repair the conditioning: a polynomial whose roots move by 10⁻³ when a coefficient moves in its last bit has roots that are not determined by its coefficients to more than three digits, and no algorithm recovers information that the data does not contain.

That is the whole of what rank is a decision says about a different quantity, and it applies here word for word: the question “what are the roots of this polynomial, given these coefficients as doubles” has an answer good to three digits, and any routine claiming more is claiming something about a polynomial it was not given.

The perturbation that does the damage, written out #

Wilkinson’s own account is worth reproducing because the arithmetic is small enough to check and because it makes the size of the effect concrete rather than asymptotic.

Perturb the coefficient of x¹⁹ — which is −210, an integer any format holds — by 2⁻²³, a relative change of about 10⁻⁹. The root at 20 moves by roughly

(2⁻²³ · 20¹⁹) / |p′(20)|,

and 20¹⁹ is 5.2·10²⁴ while p′(20) = 19! = 1.2·10¹⁷. The quotient is about 4·10⁻¹, so a change in the ninth digit of one coefficient moves a root by nearly half a unit. In Wilkinson’s original the pair at 16 and 17 leave the real axis altogether and become a complex conjugate pair — a polynomial with twenty distinct real roots, perturbed in one coefficient’s ninth digit, acquiring complex roots.

The measurement here is at the rounding level rather than at 2⁻²³, so the effect is smaller — and whether it is different in kind is a question the exact library can settle rather than assert.

The change of kind, checked both ways #

The true answer does change kind. Adding 2⁻²³ to the coefficient of x¹⁹, exactly, and evaluating the resulting integer polynomial on a grid of step 1/64 across [1, 25] with nothing rounded anywhere: there are ten sign changes and no more, at

1.016   2.000   3.016   4.000   5.016   6.000   7.000   8.000   9.156   9.516

Ten real roots where there were twenty, all of them below 9.6, and the entire upper half of the spectrum gone from the axis. That is Wilkinson’s result, reproduced with a certificate rather than an eigensolver: a sign change between two rationals is a proof of a root between them, and the absence of one over the whole of (9.6, 25] is evidence of the same strength that there is none there.

The computed answer does not change kind, and that is worse. The companion-matrix route returns maxImag of exactly zero — at every degree from 6 to 22, and at every perturbation of the x¹⁹ coefficient from 10⁻¹² to 10⁻⁶. At a relative perturbation of 10⁻⁸ the worst root is wrong by more than its own size, and the routine still reports twenty real numbers.

So the sentence to write is the opposite of the tempting one. The failure is not a real polynomial returning complex roots; it is a polynomial whose roots are half complex returning twenty reals, with no imaginary part anywhere in the output to indicate that anything of the kind has happened. A caller checking “are the roots real, as the model says they should be” gets a yes at every perturbation measured, including the ones where half of them are not.

That is the essay’s own argument with the last quantity closed. Every number a caller can see — the residual of the eigenproblem, the reality of the spectrum, the plausibility of the values — reports success on a computation whose answer has lost not merely its digits but its shape.

Where the site has met this shape before #

Three times, and collecting them is the point of putting this essay in the error field rather than in the polynomial one.

The number that decides nothing is about a determinant: a scalar that is unrepresentable long before it is uninformative, and wrong sixteen sizes before it is absent. Here the coefficients are that scalar’s relatives — 20! is a determinant-sized number sitting in a coefficient — and they fail for the same arithmetic reason.

Two condition numbers of one matrix is about a quantity that is enormous in one norm and modest in another, so that the choice of measure decides whether a problem looks hard. The root condition number is the same kind of object: it is large because the monomial basis makes it large.

And the units the matrix is measured in is the repair that sometimes works — a scaling of the variable, x → αx, changes every coefficient and can improve the conditioning. Swept here:

   α        0.1      0.2      0.5       1       2        5       10
 error   2.5·10⁻²  5.9·10⁻³ 6.1·10⁻⁴ 7.6·10⁻³ 2.1·10⁻² 1.4·10⁻² 2.8·10⁻²

Worth a factor of twelve at α = 0.5, and not worth more. It is more than “does not save this polynomial” suggests — an order of magnitude is three digits back — and it is nowhere near enough: the best α still leaves a twenty-digit answer good to four. The roots span a factor of twenty and no single α balances them, which is why the curve has a shallow optimum near a half rather than a sharp one anywhere.

The shape of that curve is worth one more sentence, because it is the same shape as every other scaling on this site and it is not the shape a reader expects from “the repair”. A scaling that fixed the problem would have a sharp minimum at a particular α, because there would be a right answer to find. This one is flat to within a factor of five over two decades of α with a shallow dip in the middle, which is what a knob looks like when it is addressing a symptom rather than a cause: it can move the balance between the coefficients and it cannot remove the fact that a degree-twenty monomial basis is a bad description of twenty roots spread over a factor of twenty.

That is the same reading the scaling that buys ten orders gives from the other end of the same idea. There the scaling recovers everything, because the problem was a mismatch of units and a unit change is exactly the inverse of it. Here it recovers a decade, because the problem is a basis and a variable scaling is not a change of basis.

Why anybody still does it this way #

Given all that, it is fair to ask why the companion-matrix route is what every library ships.

Because the alternative is worse, and because the case above is rarer than it looks. A polynomial that arrives as a list of coefficients is unusual outside of textbooks — most polynomials in practice arrive as a characteristic polynomial nobody forms, a product of factors that can be kept factored, or an interpolant in a basis that is not the monomial one. The monomial basis is the bad one: its conditioning degrades exponentially with degree for roots away from the unit circle, and the same polynomial expressed in a Chebyshev basis, with a colleague matrix rather than a companion one, is far better behaved.

Which is the general form of the lesson, and it is the same one the polynomial field spends four essays on: the trouble is in the representation rather than in the algorithm. Change the units and a quadratic eigenvalue problem loses ten orders; change the basis and a scalar polynomial loses ten digits. In both cases the object is unchanged, every routine involved is backward stable, and what moved was the description.

What to do instead, priced #

Three alternatives, and each is a real answer for a different situation.

Keep the factors. If the polynomial arrived as a product — a characteristic polynomial of a matrix, a product of linear factors from a model — never expand it. The roots are already there, and expanding is the step that destroys them. This is the same advice the inverse that is never formed gives about a matrix: the object that is easy to write down is not the object to compute with.

Change the basis. 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.

Or work exactly. For integer or rational coefficients, exact evaluation and bisection give any accuracy wanted, at the cost of BigInt arithmetic whose numbers grow with the degree. On this polynomial that is entirely practical up to degree fifty, and it is how the reference for the measurement above is produced.

What none of the three does is make the double-precision monomial route accurate. That route answers a question about a polynomial that is not quite the one asked about, and the gap between the two is the whole of the error — which is, once more, the site’s oldest sentence with the emphasis on the last three words.

Why this belongs in the error field #

The polynomial field would have taken this essay, and it is filed here instead for a reason worth making explicit.

Every quantity in it is one of the error field’s own. There is a condition number, computed from the data and not from the run. There is a backward error, at the rounding level, belonging to an object that is not the one posed. There is a forward error that is the product of the two to within a factor of nine at every degree. And there is an exactly known answer, which is what makes all three measurable rather than estimable.

The identity that runs through this site —

forward error ⪅ condition number × backward error

— is usually demonstrated on a linear system where all three terms are about the same matrix. Here the three terms are about different objects: the condition number is the polynomial’s, the backward error is the companion matrix’s, and the forward error is the roots’. Getting the attribution right requires noticing that, and the whole failure of intuition this essay records is the failure to notice it.

That is the error field’s subject rather than the polynomial field’s. What the polynomial field adds is that the same three-way mismatch occurs at every degree and every size, and that the repair — change the representation — is the same one.

What the measurement needed in order to exist #

A closing note on method, since this essay’s numbers rest on three separate exactnesses and losing any of them would have produced a plausible wrong answer.

The roots are known. They are the integers, by construction, so the forward error is a difference from a known value rather than from a better computation.

The coefficients are exact. They are expanded in BigInt, so the polynomial being discussed is the polynomial that was meant, and the rounding that occurs when they are cast to doubles is identified as a separate, measurable effect rather than absorbed into the total.

And the prediction is independent. The root condition number is computed from the coefficients and the known roots, without reference to the computed answer, so “measured equals predicted” is a statement about two quantities rather than a restatement of one.

Take away the first and the error is unmeasurable. Take away the second and the polynomial drifts before the algorithm runs, with no way to say by how much. Take away the third and the agreement between measurement and prediction becomes circular. That is a lot of scaffolding for five numbers, and it is the reason those five numbers can be believed.