Two approximants and one matrix size

An error committed before the arithmetic settles the region and the pole placement. What is left is the choice everybody thinks of first: rational or polynomial.

The answer is that it depends on the region, and the measurement that shows so has to be set up carefully, because the two approximants are not naturally comparable.

Making the contest fair #

A rational approximant with m poles linearises to a matrix of n(m+1) rows. A polynomial of degree d linearises to nd. So degree is not the quantity to fix — a polynomial of degree twenty against a rational with ten poles is a comparison of a 120 × 120 matrix with a 66 × 66 one, and the larger matrix wins for reasons that have nothing to do with the basis.

The eigensolve is what the whole chain is charged for, and it is O(size³). So the fair contest fixes the size, which means d = m + 1, and everything else — the sample points, the fitting procedure, the error measure — is held identical. Same geometric grid, same least-squares QR, same maximum over the same interval.

The result #

Ratio of the polynomial’s error to the rational’s, at equal size:

Two different pictures.

On a target set that stands clear of the branch point, the ratios flip either side of one as the size grows, which looks like a comparison between two equally good methods. Read at more sizes it is not that at all, and what it is instead is the subject of a section below.

On a target set that reaches towards it, the rational pulls away monotonically, to two orders by seventy-eight rows and still widening.

The mechanism is the one the previous essay set out. Away from the singularity √ is analytic on a neighbourhood of the target set and a polynomial converges to it geometrically — as fast as a rational, because there is nothing for the poles to do. Approaching the singularity the polynomial’s rate degrades towards algebraic while the rational’s stays geometric, because the poles can sit on the cut and the polynomial has nothing to sit anywhere.

Neither method is better, which is the same verdict six routes to one spectrum reaches about six linearisations of a quadratic and for a related reason: what varies between them is small beside what varies with the decision made before either was chosen. The question is different, and which question it is was decided by whoever said where the eigenvalues might be.

The easy row is one shape, not a flat one #

Five sizes are too few to tell a flat line from a rise followed by a collapse, and the clear target set is the second. At eleven sizes, with everything else held:

The ratio rises monotonically from 0.18 to 4.45 across the first seven sizes — a factor of twenty-four — which is the same shape the hard row has and not the flat one the five-point table suggested. So the rational is pulling away on the clear target set too; it is pulling away more slowly, which is what the mechanism predicts, and the earlier reading of that row as two equally good methods was reading a curve through its crossing point.

And sweeping the target set rather than the size shows the same crossing from the other side.

The bottom of the range crosses one between reach 0.3 and 0.1, and the top peaks in the middle. Across reaches of 1, 0.5, 0.3, 0.1 and 0.02 the ratio’s range runs 0.45–0.0, 0.71–27.2, 0.87–89.4, 1.14–83.0 and 1.33–49.2. The lower end climbs monotonically — 0.45, 0.71, 0.87, 1.14, 1.33 — so the reach at which the rational wins even at the smallest size is between 0.3 and 0.1.

The upper end is not monotone, and that is the interior optimum showing through. It runs 0.0, 27.2, 89.4, 83.0 and 49.2 — peaking at reach 0.3 and falling by nearly half by 0.02. The rational approximant’s advantage is therefore largest on a moderately hard target set rather than on the hardest one, because at reach 0.02 the rational is running into the same floor the size sweep found at m = 8 while the polynomial is still descending. So neither the reach nor the degree drives the comparison monotonically, and the two non-monotonicities are the same one seen along two different axes.

What ends the rise is not the polynomial improving on merit. Look at the first row: the rational’s error falls by six orders up to m = 8 and then stops — 6.76·10⁻⁹, 7.28·10⁻⁹, 1.21·10⁻⁸, 5.96·10⁻⁹, with four more degrees buying nothing. That is exactly the interior optimum the section below describes, and it is where the clear row’s flip comes from: the rational hits its own floor while the polynomial is still falling by an order every two degrees.

So the two phenomena this essay describes separately are one phenomenon. There is a single shape — the rational pulls away at a rate the target set sets, until its fitting basis runs out of independent directions — and the difference between the rows is only where in the range of sizes anyone would use that happens. On the clear set it happens at m ≈ 10; on the set that reaches the cut the rational is still improving at m = 14, from 9.9·10⁻² down to 3.0·10⁻⁶, so the ratio there rises the whole way to 234.

One conclusion below has to be strengthened rather than softened by this. On a problem whose eigenvalues are known to lie well away from the singularity, the polynomial is the right answer for reasons of simplicity rather than of accuracy — measured past the rational’s floor, the polynomial is right for accuracy as well, by a factor of 23 at seventy-eight rows and 216 at ninety. It is not a tie there and it is not a wash; the polynomial is simply better, because the thing it is being compared against has stopped improving.

There is a rule of thumb in it that the five-point table could not have produced. The rational’s floor arrives when the number of poles exceeds the number of independent directions the region supports, and the height of the floor is what a fit at that many poles achieves — 4·10⁻⁸ here. So a caller who knows their target set can ask the question the other way round: what accuracy does the rational reach before it stalls? If that is enough, the rational is the right basis and the size at which it stalls is the size to use. If it is not, no number of poles will supply the missing digits and the polynomial is the only one of the two that can still be pushed. Neither of those decisions needs an eigensolve, and both are settled by the fit alone.

assertTheEasySetsTieIsTheRationalsOwnFloor holds the monotone rise on both target sets, the rational’s stall on the clear one, the polynomial’s continued fall over the same degrees, and the fact that the hard set’s rational is still improving at every degree measured — which is what makes the two rows one shape read at two places rather than two findings.

The polynomial’s linearisation, which is not a companion matrix #

The polynomial approximant is fitted in the Chebyshev basis, and it linearises to a comrade matrix rather than a companion one. The two constructions are identical except for the three-term recurrence in the corner: a companion matrix has the monomial recurrence xᵏ⁺¹ = x·xᵏ and a comrade matrix has 2xTₖ = Tₖ₊₁ + Tₖ₋₁.

Which one a code uses is decided by which basis the coefficients arrived in, and the difference between the two bases is the roots are not the coefficients: a polynomial of degree twenty in the monomial basis has coefficients spanning twenty orders and one in the Chebyshev basis does not.

Using the monomial basis here would confound the comparison badly. The polynomial approximant would lose, and it would lose for a reason about coefficient representation rather than about approximation power — which is exactly the sort of thing this site insists on separating. So the polynomial is given its best available form and still loses where the rational wins.

The interior optimum #

Adding poles improves the approximation and enlarges the matrix. Both of those cost accuracy in opposite directions, so the total error has a minimum somewhere in between.

It is visible in the measurements. On a target set clear of the cut, the rational approximant’s error runs 1.9·10⁻⁴, 2.9·10⁻⁶, 4.3·10⁻⁸, 6.8·10⁻⁹ — and then 1.2·10⁻⁸ at twelve poles, which is worse than at ten. The fit’s basis has become ill conditioned: twelve basis functions 1/(λ − ξⱼ) on a region that supports about ten independent directions is a least-squares problem whose coefficients are no longer determined.

The same effect appears at the other end of the chain. The residual against the approximant — the free one from the problem the solver was actually given — stays at 10⁻¹⁵ at two poles and rises to 1.5·10⁻¹² at seven, because the linearisation is larger and worse scaled. So the arithmetic term grows while the approximation term falls, and the crossing is where more poles stop buying anything.

That is a familiar shape from the last digit is the cheapest inverted: there each additional digit costs less than the last, and here each additional pole buys less than the last and eventually buys negative. What is unusual is that both curves are computable in advance — ‖g − r‖ from the fit, and the linearisation’s conditioning from its own norms — so the optimum can be located without running the eigensolve at every degree.

What the choice costs when it is made wrong #

Both directions are worth pricing.

Rational on an easy target set costs nothing in accuracy and costs the machinery: poles have to be placed, a rule for placing them has to be chosen, and the fit’s conditioning has to be watched. A polynomial of the same size does as well with none of that. So on a problem whose eigenvalues are known to lie well away from the singularity, the polynomial is the right answer for reasons of simplicity rather than of accuracy.

Polynomial on a hard target set costs two orders at seventy-eight rows, and the loss grows with the size. Doubling the matrix to recover them is an eightfold cost in the eigensolve, which is the expensive step. So this direction is the one worth being careful about, and it is the one a code defaults into: a polynomial approximant needs no decisions and is therefore what gets written first.

The asymmetry says which way to lean. When the region is known to be clear of the singularity, the simpler approximant. When it is not known, the rational one — because the cost of being wrong is much larger in one direction than the other, and not knowing where the eigenvalues are is the ordinary situation — which is the same asymmetry the swap that is not optional prices for pivoting, where the cheap precaution is taken always because the case it guards against cannot be detected in advance.

The comrade matrix, written out #

It is worth writing the construction down once, because it is short and because the way it is built from a recurrence rather than from a shift is the whole difference between the two bases.

For a scalar polynomial p(t) = Σ cₖ Tₖ(t), the comrade matrix is d × d with a symmetric tridiagonal skeleton — a 1 in the first row’s second entry, and ½ on both off-diagonals below — and the last row corrected by −cⱼ/(2c_d) for j = 0 … d−1. Its eigenvalues are the roots of p. The companion matrix has the same shape with the tridiagonal replaced by a subdiagonal of ones, which is the monomial recurrence.

For a matrix polynomial the same construction goes through blockwise, with each scalar coefficient replaced by its coefficient matrix. Here every coefficient is a multiple of the identity except the constant term, which carries A, and the linear term, which carries the −λ from the eigenvalue problem itself. The result is nd × nd and its eigenvalues, mapped back through the affine change of variable that put the target set on [−1, 1], are the eigenvalues of the approximated problem.

Two details are load-bearing and easy to get wrong. The affine map has to be folded into the coefficients before the comrade matrix is formed — q(t) = p(t) − αt − β, where λ = αt + β — rather than applied to the eigenvalues afterwards, because the eigenvalue problem’s own λ is not in the Chebyshev variable. And the leading coefficient divides the correction row, so a fit whose top Chebyshev coefficient has been driven to rounding produces a comrade matrix with an enormous last row: a degree chosen higher than the fit can support fails here rather than in the approximation.

That second point is the interior optimum arriving from a third direction. The approximation error falls with degree, the linearisation’s size grows with degree, and the construction of the linearisation degrades with degree once the fit stops determining its own top coefficient.

What neither of them fixes #

Both approximants have the same relationship to the spurious eigenvalues that come back with the answer. A rational approximant with m poles produces n(m+1) eigenvalues of which n are wanted; a polynomial of degree d produces nd of which n are wanted. Neither construction returns fewer, neither labels them, and the ratio of noise to answer grows with the size in both cases. The eigenvalues that are answers to nothing is that measurement, and it applies to both bases equally.

And neither of them changes what the residual against the approximant means. It is a backward error for whichever object was linearised, it is at rounding in both cases, and it is silent about the approximation in both cases. The choice of basis moves one term in the error decomposition and leaves the reporting problem exactly as it was.

Why this comparison is worth making at all #

A reader might reasonably ask why anyone would run both, given that the answer depends on a region whose position is exactly what is not known.

The reason is that the comparison is cheap and the decision is not otherwise available. Fitting both approximants at one size costs two least-squares problems on a few hundred sample points — call it a millisecond — against an eigensolve of an n(m+1) matrix, which is the whole cost of the method. So a code can fit both, compare the two maxima, and take the better one, at a cost that does not appear in any accounting of the method.

That is a better answer than choosing by rule, because the rule requires knowing where the region sits relative to the singularity, and the singularity’s position is sometimes as unknown as the eigenvalues’. Fitting both measures the thing the rule would have had to be told.

It also gives a diagnostic. If the two maxima are within a factor of two of each other, the target set is clear of whatever difficulty the function has and the problem is easy in the sense that matters. If they differ by an order, the region is close to a singularity, and that is worth knowing independently of which approximant gets used — because it is also the signal that the target set should be narrowed, which an error committed before the arithmetic shows is worth more than either choice.

This is the pattern this site keeps arriving at and it is worth naming once more: two routes to the same quantity, run together, cost almost nothing and answer a question neither of them was asked. Here the question is how hard the problem is, and the answer is a ratio of two numbers that were going to be computed anyway.

The size that is the real cost #

One accounting note, because “equal linearisation size” has been doing a lot of work and its consequences are worth spelling out.

An eigensolve of an N × N matrix is O(N³), and N = n(m+1) here. So doubling the number of poles is roughly an eightfold cost, and the whole method is dominated by that one step: the fit is milliseconds, the linearisation is O(N²) to assemble, and the factorisation is everything else.

That changes what “a better approximation” means in practice. Going from ten poles to twelve on the easy target set buys nothing in accuracy — the interior optimum is already passed — and costs a factor of about 1.7 in the eigensolve. Going from four poles to ten on the hard one buys four orders and costs a factor of about eight. The second is obviously worth it and the first obviously is not, and the point is that both are decidable from the fit alone, before the expensive step.

For a large sparse problem the accounting is different again, because nobody forms an n(m+1) matrix and factorises it — the linearisation is solved by a Krylov method, whose cost is linear in N per iteration and whose iteration count depends on the linearisation’s conditioning. Then the growth term is not the cube but the degradation of the conditioning with m, which is the same quantity the rising residual measures and which this essay has priced only indirectly. That is the honest limit of what is measured here: the dense accounting is exact and the sparse one is inferred.

The one line #

Fix the matrix size, not the degree. Fit both bases, compare the two maxima, and take the better one — it costs a millisecond against an eigensolve. If they agree, the region is easy and the simpler approximant is right; if they differ by an order, the region is close to a singularity and that is worth knowing for its own sake.