Numerical Linear Algebra
where the algebra stops being the arithmetic
Reduction, and what a model is for

Where to put the poles of a rational function

Three times in one field the same question arrives from different directions — ADI shifts, rational approximation of a square root, the decay of a Gramian — and it has one answer. Cluster them geometrically towards wherever the function is difficult, and the alternative that looks reasonable costs orders.
35 min read 32 figures Identical algebra, different arithmeticMeasured, not assertedWhere the algebra stops being the arithmetic

Worth reading first: Why a Gramian can be truncated at all · The bound that is known in advance.

Why a Gramian can be truncated at all reduces the decay of a Gramian to a question about rational approximation and leaves the question itself unanswered. This essay answers it, and the answer turns out to be the same one three different problems in this collection arrive at from three directions.

The question is: given k poles to place, where do they go?

The rational function 6 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 6 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.1187. The upper curve uses 6 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.8922 — a factor of 7.52 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts6κ of the spectrum389geometric max0.12equally spaced max0.89the factor between7.5the same k solvesand one choice of where
Fig. 1 The rational function six ADI shifts amount to, chosen two ways. One of them equioscillates and one of them never approaches the small end.

The objective

The ADI iteration’s error after k steps is max|r| over the spectrum of −A, squared, where

r(λ) = ∏ⱼ (λ − pⱼ)/(λ + pⱼ)

and the pⱼ are the shifts. Every factor is a number that is zero at its own shift and approaches one far from it, so r is small near the shifts and large between them. The whole design question is how to space the zeros so that the largest value between any two of them is as small as possible.

That objective has a name in approximation theory — it is the Zolotarev problem — and it has a classical answer. What is worth measuring rather than citing is how much the answer is worth against the choice a working programmer makes without one.

The two choices

The spectrum here runs from 9.9 to 3855, a ratio of 389. Six shifts:

Geometric. pⱼ = a(b/a)^((2j+1)/2k), spaced evenly on a logarithmic scale. The resulting curve touches zero six times and rises to the same height between every consecutive pair — it equioscillates, which is the signature of an optimum. Maximum: 5.2·10⁻².

Even. pⱼ spread linearly across [a, b]. On a spectrum spanning three decades, every one of the six lands in the top decade: the smallest is at 331, and there are no shifts at all below it while the spectrum continues down to 9.9. The curve is small at the top and essentially one at the bottom. Maximum: 8.5·10⁻¹.

Sixteen times worse, in a quantity the iteration pays squared. Measured on the iteration itself after eight steps: 2.0·10⁻³ against 5.9·10⁻¹, a factor of 290, for the same eight solves and the same recurrence, with one line different.

Why even spacing fails, precisely

It is not that the shifts are in the wrong place on average. It is that the objective is a maximum and the spectrum is logarithmic.

The factor |(λ − p)/(λ + p)| for λ far below every p is close to one, and it stays close to one no matter how many shifts are added above. So the value of r at the bottom of the spectrum is decided entirely by whether any shift is near the bottom, and adding shifts to the top decade of a spectrum spanning three does nothing about it. Adding a seventh, eighth and ninth evenly spaced shift moves the maximum from 0.85 to 0.83.

The eigenvalues down there are not incidental. They are the slow modes, and in a Gramian they are exactly the directions with the largest energy — the ones the accuracy that is thrown away would call the part of the answer worth computing accurately — the ones a low-rank factor most needs to represent. So the failure is not a small loss at an unimportant end of the spectrum; it is a total loss at the important end.

This is the same reading coarsening in one direction only gives to an anisotropic multigrid smoother, and the same one a direction the smoother cannot see gives to the error it leaves behind: a method that treats a spectrum uniformly when the spectrum is not uniform leaves the part it never addressed, and the part it never addressed is the part that matters.

Whether geometric is optimal, measured rather than cited

Geometric spacing is not exactly the Zolotarev optimum, which is expressible through elliptic functions. The gap is worth measuring, and the measurement has to be set up carefully.

A descent started at the geometric shifts that finds no improvement proves nothing — it is indistinguishable from a descent that does not work. So the search starts from the evenly spaced shifts, which are known to be bad, and is allowed to walk.

It walks them from 8.9·10⁻¹ to 1.9·10⁻¹, which is a real improvement and establishes that the descent works. And it lands beside the geometric value of 1.2·10⁻¹ rather than beating it — within a factor of 1.6.

Two conclusions. A code without an elliptic-function library loses about sixty per cent by using geometric shifts, which against a factor of 290 for getting the spacing wrong is nothing. And the objective has a wide flat optimum: many sets of shifts are within a factor of two of the best, which is why the spread is the decision and the exact placement is not.

That is a reassuring shape and it is not universal. It is worth contrasting with four knobs and one floor, where four different rules for a regularisation parameter also land within a few per cent of each other — and where the reason is the opposite one. There the rules agree because none of them knows anything; here they agree because the objective is genuinely flat near its optimum, which is a property of the problem rather than of the rules.

The same question, three times in one phase

The reason this essay exists rather than being a paragraph in the previous one is that the identical question and the identical answer arrive three times across this collection, from three directions that share no machinery.

ADI shifts, above: poles of a rational function that has to be small on the spectrum of −A. The answer is geometric clustering across the interval.

The decay of a Gramian’s spectrum, which is the previous essay: the bound is the value of the best such rational function, so the same optimisation appears as an explanation rather than as an algorithm.

Rational approximation of a square root, which is an error committed before the arithmetic in the polynomial field: a nonlinear eigenvalue problem whose function has a branch cut, approximated by a rational function whose poles have to go on the cut. The answer there is geometric clustering towards the branch point, and the measurement is that an evenly spaced set of poles is orders worse on a target set that stands clear of the singularity.

Three problems, one answer, and the answer is always the same sentence: cluster geometrically towards wherever the function is difficult. For ADI the difficulty is the small end of the spectrum, where the factors approach one. For the square root it is the branch point. In both cases a linear spacing puts everything where the function is easy.

That recurrence is worth carrying as a prior. On meeting a construction that asks for a set of poles or shifts or interpolation points, the first guess should be geometric towards the hard end, and the second should be to find out where the hard end is.

The limit of the clustering, which is a real one

There is a place the rule stops, and the polynomial field’s version of this problem is where it is visible.

Clustering the poles closer to the singularity than the target set reaches buys nothing and costs conditioning. Two poles nearer the branch cut than any point where the function is being fitted give two basis functions the fit cannot tell apart, and the least-squares problem that produces the approximation becomes ill conditioned. The measurement in that field: clustering wins by a factor of two to nine where there is room, and by nothing at all on a target set that touches the singularity.

Here the equivalent limit is the interval itself. Shifts outside [a, b] are wasted, and shifts clustered below a are approximating a part of the spectrum that is not there — the same waste a parameter that counts steps measures for an iteration continued past the point where its own quantity has stopped moving. The geometric rule has the interval built into it, which is why it does not run into the limit — but it is worth knowing that the rule is clustering relative to the region of interest rather than clustering towards a point, and the two coincide only when the region reaches the point.

Where the shifts come from when the ends are not known

The rule needs a and b, and on a large problem neither is available exactly. That is not a problem and it is worth saying why, because “estimate the extreme eigenvalues first” sounds like it might be as expensive as the thing it precedes.

It is not. b — the largest — is what a handful of power-method steps returns, and the estimate can be loose in the safe direction: taking b too large adds shifts above the spectrum, which are wasted but harmless. a — the smallest — is harder, and is what an inverse power method or a few Lanczos steps gives. Taking a too small is the dangerous direction, since it spends shifts on a region that is not there and leaves the actual small end thinly covered.

So the practical rule is asymmetric: overestimate b freely and underestimate a cautiously. A code that has no estimate at all can bracket with ‖A‖ and ‖A⁻¹‖⁻¹ — the first from a norm and the second from one solve — which gives an interval containing the spectrum and possibly wider than it. A wider interval means a slower rate and a correct answer, which is the right way for an estimate to fail.

What the rule cannot survive is a spectrum that is not an interval. If A has complex eigenvalues the region is two-dimensional and the shifts have to be chosen in the plane; the Zolotarev problem over a region is a genuinely harder object, the geometric rule has no direct analogue, and the practical answer is a heuristic rather than a theorem. That case is where this field’s remaining difficulty lives and it is worth marking rather than glossing.

What the equioscillation certifies

The signature of the geometric choice on the figure is that every hump between consecutive zeros reaches the same height, and that is not decoration — it is the certificate.

If one hump were lower than the others, the shifts could be moved to trade a little of it for a reduction of the highest, and the maximum would fall. So an optimum has to equioscillate, and any set of shifts whose curve has one hump noticeably below the rest is not optimal and can be improved by inspection. That is a check anyone can run on a figure without knowing anything about elliptic functions.

It is also why the search from the evenly spaced shifts converges to something that looks geometric. The descent is doing, numerically and slowly, what the equioscillation condition does in closed form: pushing shifts towards regions where the curve is high until every region is equally high. Watching it land within a factor of 1.6 of the geometric answer is the strongest evidence available here that the geometric rule is approximating the right object rather than merely being a respectable convention.

The same certificate appears in the best approximation there is in a different guise — there the optimality of a truncated SVD is certified by what the residual looks like rather than by a search — and the habit generalises: an optimum usually leaves a visible signature in the object it produces, and finding out what that signature is worth more than any amount of trust in the derivation.

What it costs to get this right

Nothing. That is the part worth stating flatly.

The shifts are k numbers computed from two: the ends of the spectrum of −A, which can be estimated by a few steps of a power method and its inverse. The iteration is identical either way, the number of solves is identical, and the difference in the answer is a factor of 290.

That is the same accounting the scaling that buys ten orders does for a change of units in the polynomial field: two lines, computed from norms already in hand, worth more than every subsequent algorithmic choice put together. This site keeps finding that the free decisions dominate the expensive ones, and this is the clearest instance of it in this field.

The rule, and the three places it has been checked

Stated once, plainly, since the whole essay is an argument for one sentence.

Given k poles to place so that a rational function is uniformly small over a region, space them geometrically across the region, clustered towards wherever the function being approximated is difficult. For an ADI iteration the region is the spectrum of −A and the difficulty is at the small end; for a rational approximation of a square root it is a target interval and the difficulty is the branch point; for the bound on a Gramian’s decay it is the same interval as the first and the statement is an explanation rather than an algorithm.

The three have been checked separately and they agree. The evenly spaced alternative is sixteen times worse in the first, a factor of two to nine worse in the second where there is room to cluster, and produces the same bound in the third because the third is the optimum rather than a choice.

What is worth noticing about the agreement is that the three problems share no code. The ADI shifts and the Gramian bound share a rational function; the square-root approximation shares neither the function nor the algorithm nor the field. That the same rule wins in all three is evidence that the rule is about rational approximation rather than about any of the settings it appears in — which is why it is stated here as a rule rather than as a property of an iteration.

The corollary is a working habit. On meeting a construction that takes a set of poles, shifts, or interpolation points as an input, assume geometric clustering until measured otherwise, and spend the effort on finding out where the hard end is rather than on refining the spacing. On the evidence here the first decision is worth two orders and the second is worth sixty per cent.

At other settings

The rational function 2 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 2 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.6181. The upper curve uses 2 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.9731 — a factor of 1.57 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts2κ of the spectrum389geometric max0.62equally spaced max0.97the factor between1.6the same k solvesand one choice of where
Fig. 2 Two shifts, where the equioscillation is easiest to read.
The rational function 4 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 4 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.2711. The upper curve uses 4 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.9341 — a factor of 3.45 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts4κ of the spectrum389geometric max0.27equally spaced max0.93the factor between3.4the same k solvesand one choice of where
Fig. 3 Four, and the geometric curve has four zeros and three equal humps.
The rational function 8 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 8 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.05198. The upper curve uses 8 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.8497 — a factor of 16.3 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts8κ of the spectrum389geometric max0.052equally spaced max0.85the factor between16the same k solvesand one choice of where
Fig. 4 Eight, where the evenly spaced curve has hardly moved from its value at two.
The rational function 12 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 12 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.009967. The upper curve uses 12 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.766 — a factor of 76.9 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts12κ of the spectrum389geometric max0.01equally spaced max0.77the factor between77the same k solvesand one choice of where
Fig. 5 And twelve, where the two are five orders apart at the small end of the spectrum.
The ADI error is the rational function, not merely bounded by it‖P − ZZᵀ‖/‖P‖ after k steps of low-rank ADI on a 30-state model, against max|r|² for the same shifts, with P taken from its closed form rather than from another solve. The two curves are parallel: their ratio is 0.756 at k = 2 and 0.746 at k = 10, a spread of 1.0340 while the error itself falls from 0.289 to 3.87·10⁻⁴. A bound that is a constant factor from its quantity over four orders is an identity with a projection in front of it — the constant is how much of P lies where |r| is largest, which depends on the system and not on k.234567891010⁻⁴10⁻²1ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.77spread of the ratio1error at k = 103.9·10⁻⁴its residual1.1·10⁻⁴four orders of errorand one constant
Fig. 6 What the shift choice buys, measured on the iteration rather than on the function.
The ADI error is the rational function, not merely bounded by it‖P − ZZᵀ‖/‖P‖ after k steps of low-rank ADI on a 16-state model, against max|r|² for the same shifts, with P taken from its closed form rather than from another solve. The two curves are parallel: their ratio is 0.759 at k = 2 and 0.748 at k = 10, a spread of 1.0384 while the error itself falls from 0.193 to 5.03·10⁻⁵. A bound that is a constant factor from its quantity over four orders is an identity with a projection in front of it — the constant is how much of P lies where |r| is largest, which depends on the system and not on k.234567891010⁻⁵10⁻³10⁻¹ADI steps krelative error, and max|r|²max|r|²the measured errora bound that is an identitysteps drawn9ratio, worst0.78spread of the ratio1error at k = 105·10⁻⁵its residual1.6·10⁻⁵four orders of errorand one constant
Fig. 7 The same on a narrower spectrum, where every choice is better.
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 30-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 2.48·10⁻⁸ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.86 and 3834, whose ratio is 388.8 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates30κ of the spectrum389λ11 ÷ λ₁2.5·10⁻⁸the bound there5.2·10⁻⁴the cliff everything rests onand the reason for it
Fig. 8 The decay the best shifts are a bound on.
Why a Gramian can be truncated at all: a rational approximation problemλₖ₊₁(P)/λ₁(P) for the controllability Gramian of a 40-state model, against Z_k², where Z_k is the smallest a degree-(k,k) rational function with poles on one side can be made on the other. The Gramian's spectrum falls from 1 to 1.11·10⁻⁷ in 10 steps, and the bound falls with it. Nothing about the model enters the bound except the two ends of the spectrum of −A, 9.86 and 6714, whose ratio is 680.6 — so the decay that makes model reduction possible is predictable from two numbers before the Gramian exists. That is the fact this whole field rests on and it is a fact about rational approximation.1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates40κ of the spectrum681λ11 ÷ λ₁1.1·10⁻⁷the bound there0.001the cliff everything rests onand the reason for it
Fig. 9 And on a wider spectrum, where more shifts are needed for the same accuracy.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from 0.602 to 8.99, and its left end is 1 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 2.22·10⁻⁵; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 1.71·10⁻⁴. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.110¹10⁻⁸10⁻⁶10⁻⁴distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.3left end, from the cut1rational, worst2.2·10⁻⁵polynomial, worst1.7·10⁻⁴linearisation size, both54committed before the solveand invisible to it
Fig. 10 The same question in the polynomial field, with a branch cut instead of a spectrum.
The error that happens before any arithmetic: |g − r| on the target setγ√(λ + c) minus its approximant, plotted against the distance from the branch point at −0.4. The target set runs from -0.2 to 8.99, and its left end is 0.2 from the branch point — a choice made before anything is computed. The rational approximant with 8 poles clustered on the cut reaches 8.46·10⁻⁴; the polynomial of degree 9, which linearises to a matrix of the same size, reaches 0.00676. Both errors are committed at this stage. No residual computed later contains either of them, because a solver's residual is taken against whichever of these two functions it was handed.110¹10⁻⁷10⁻⁵10⁻³distance from the branch point, λ + c|g − r|the eigenvaluesdegree 98 poleschosen, not computedreach0.06left end, from the cut0.2rational, worst8.5·10⁻⁴polynomial, worst0.0068linearisation size, both54committed before the solveand invisible to it
Fig. 11 And a target set that reaches the singularity, where the clustering runs out of room.
What two lines of scaling are worth: the same quadratic in nine systems of unitsThe forward error against the closed form for an overdamped chain of 8 masses, before and after Fan–Lin–Van Dooren scaling. The unscaled curve runs 7.49·10⁻¹⁴ to 0.001258 — every digit gone by the far end — and the scaled one runs 1.26·10⁻¹³ to 8.23·10⁻¹⁴, flat to within a factor of 2.01. The scaling is γ = √(‖K‖/‖M‖) and δ = 2/(‖K‖ + γ‖C‖), computed from three norms and nothing else, and the map back is λ = γμ with no rounding in the statement. Flatness is the half that matters: after scaling every stop of the sweep IS the same problem, so there is nothing left for the change of units to do.0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were
Fig. 12 A free decision worth ten orders, in a third field.
Two condition numbers of one 8×8 system, as its rows are put into different unitsFour curves against the spread of the row units, in decades. κ_∞ of the scaled matrix rises from 9.83 to 1.9·10¹⁰ while the componentwise condition number stays at 6.98 throughout — the same system, the same solution, and one of the two numbers is a fact about the units. Hilbert's two numbers are drawn flat beside them at 3.4·10¹⁰ and 1.2·10¹⁰: a matrix whose sensitivity no scaling repairs.0246810110²10⁴10⁶10⁸10¹⁰10¹²spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread9.8κ_∞ at 10 decades1.9·10¹⁰cond, either end7Hilbert, equilibrated1.3·10¹⁰the solution is the same at every spreadand one of these curves knows it
Fig. 13 And the units that decision is about.
Hankel singular values, and the half of them a product cannot seeσₖ divided by σ₁, for a 20-state model of McMillan degree 12. The lower curve is the square-root route — Cholesky-like factors of the two Gramians and one SVD of RᵀS — and it descends to the level at which the Gramians were computed. The upper curve eigendecomposes the product PQ, agrees for the first 6 values and then stops, flattening at 9.67·10⁻¹⁰. The dashed line is σ₁√u = 2.32·10⁻⁹, which is where a route that squares the conditioning must stop: the eigenvalues of PQ are σ², an absolute error of uσ₁² on them is a relative error of √u on σ. κ(PQ) is 1.8·10³⁶ against 1.34·10¹⁸. Every σ below the line is a term in the error bound, so this is not an academic loss.1357911131510⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1index kσₖ ÷ σ₁σ₁√utwo factorstheir productthe normal equations, againroutes agree to k =6product floors at9.7·10⁻¹⁰σ₁√u2.3·10⁻⁹κ(P)κ(Q)1.8·10³⁶√ of it1.3·10¹⁸do not form the productthe σ below the line are the bound
Fig. 14 The numbers a truncated Gramian eventually produces.
The error a reduced model has, and the bound computed before it existed‖H − Hᵣ‖∞ for balanced truncation of a 20-state model, against 2Σ_{k>r}σₖ — a number available from two Lyapunov solves before any reduced model is formed — and against σᵣ₊₁, which is a lower bound no reduced model of that order can beat. The measured error and the upper bound are the same curve: the ratio is 1.0000 to 1.0000 across the sweep, spread 1.0000. The bound is attained rather than approached, because on a spectrum falling this fast the tail sum is its own first term and removing one state costs exactly twice it. The measured error sits at 2.03 times the lower bound, so the whole question of what order to truncate at is answered by a curve nobody has to compute a model to see.123456710⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹order r kepterror and bounds2Σσ, and the errorσᵣ₊₁a bound that is an equalityorders7bound ÷ error, worst1spread over the sweep1error ÷ σᵣ₊₁, worst2.1computed before the modeland attained by it
Fig. 15 And the bound they are a sum of.
Smoothing factor against anisotropy, at ω = 0.667Three curves of the smoothing factor against the anisotropy parameter on a logarithmic axis. One rises to one as the anisotropy grows; the other two coincide and stay near a third.10⁻⁴10⁻³10⁻²10⁻¹100.250.50.751anisotropy εsmoothing factor μpoint + fully-line + fullpoint + semi-y⅓, the one-dimensional answertwo routes, three curvesgap between the repairs0scan against closed form1.1·10⁻¹⁶the dashed curve lies on the solid one beneath itone repair, written two ways
Fig. 16 A method that treats a spectrum uniformly when it is not.
What one weighted Jacobi sweep does to each mode, at ω = 0.67A curve of the damping factor against mode frequency, falling from one at the left through zero and towards minus one at the right. The right-hand half of the plot is shaded, and two horizontal lines mark the worst damping inside it.00.250.50.751-1-0.500.51mode frequency θ / πdamping factor per sweeppredicted±0.333measureda coarse grid seestwo routes to one factorsmoothing factor, scanned0.33smoothing factor, closed form0.33worst mode disagreement3.3·10⁻¹⁶63 interior points, one sweepthe left-hand end is what the coarse grid is for
Fig. 17 And what it leaves behind.
The fine operator: five pointsA three-by-three arrangement of discs carrying the stencil's coefficients, with the zero positions drawn small, beside a bar chart of stored entries per row at each level of the hierarchy.0-0.250-0.251-0.250-0.250the discretisationstored entries per rowlevel 0 · 31×314.87/rowlevel 1 · 15×158.22/rowlevel 2 · 7×77.37/rowlevel 3 · 3×35.44/rowlevel 4 · 1×11.00/rowstill a stencil, still annihilates a constantentries in an interior row5weight outside the 3×30row sum0the isotropic model problemfive points
Fig. 18 Coarsening in the direction the problem is coupled in.
Two filters on one sum, λ = 0.001The weight each term of the solution is given, against its index. Truncation is a step: one for the first 32 terms and zero after. Tikhonov is σ²/(σ² + λ²), which falls smoothly through the same place. The unregularised solution is the constant one, which is why it divides noise by a σ of 1.7·10⁻¹³.081624324048566400.250.50.751index kfilter factor fₖno regularisation: fₖ = 1truncationTikhonovthe same sum, three weightsTikhonov, relative error0.16truncation, relative error0.13no filter at all5.5·10⁸both filters are one expression with a different weightfₖ = 1 is the catastrophe
Fig. 19 Another rational filter, applied to a spectrum for a different reason.
The filter 12 conjugate gradient steps apply, measured and predictedFilter factors against the singular-value index. The factors measured off the iterate and the polynomial predicted from the recurrence coefficients agree to 10⁻¹³ and are drawn as one curve. It rises above one — 1.120 at its largest — and changes direction several times. Tikhonov's filter at the matching cutoff is monotone and never exceeds one.081624324048566400.250.50.7511.25index kfilter factoronezeroCG, both routesTikhonovone of these is not a weightlargest CG factor1.1measured vs predicted10⁻¹³largest Tikhonov factor0.99two routes to the same curveand a curve that goes above one
Fig. 20 And a polynomial one, which is what a Krylov method builds.
What 3 shifts do to each direction of the starting vectorAmplification of each eigen-component against its eigenvalue, on a logarithmic vertical axis. The points are measured — the component of the vector along each eigenvector, before and after — and the curve is |Π(λ − θ)| with the roots at the 3 discarded Ritz values. They agree to 3.9·10⁻¹¹. The 3 wanted directions are amplified by at least 1 and everything else by at most 0.802.1357910⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1eigenvalue λamplificationkeptdiscardeda filter on the starting vectormeasured against the polynomial3.9·10⁻¹¹worst kept direction1best discarded direction0.8the roots are the discarded Ritz valuesand the vertical lines are where they sit
Fig. 21 Shifts chosen to filter a spectrum rather than to approximate a function.
Two routes to e^A as a 8×8 bidiagonal matrix walks up to a defective oneA_δ is upper bidiagonal with diagonal λ, λ+δ, …, so its eigenvalues are exact and its eigenvectors have a closed form. As δ falls the eigendecomposition route's relative error rises to 2.92·10⁶⁵ while scaling and squaring's falls to 4.09·10⁻¹², which is the distance ‖A_δ − A₀‖ = 7·10⁻¹² to the matrix whose exponential is known in closed form. The two agree at the right of the axis, where the matrix really is far from the limit; the divergence is κ(V), which reaches 3.3·10⁸².10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²10⁻¹⁴10⁻¹10¹²10²⁵10³⁸10⁵¹10⁶⁴δ, the gap between consecutive eigenvaluesrelative error in e^Aan answer with no correct digitsV f(Λ) V⁻¹scaling and squaring‖A_δ − A₀‖exact eigenvalues throughoutκ(V) at the smallest δ3.3·10⁸²eigen route2.9·10⁶⁵scaling and squaring4.1·10⁻¹²distance to the limit7·10⁻¹²the eigenvalues are the diagonaland they are exact at every stop
Fig. 22 Rational approximation of a matrix exponential, one field over.
Relative error of scaling and squaring against the number of squarings, μ = 12Each point computes e^(A/2ˢ) by a Padé approximant and squares it s times, against a closed-form exponential. At s = 0 the approximant is being used far outside its range and the error is 5.55·10⁻¹¹; the best is 9.08·10⁻¹⁶ at s = 2; at s = 15 it has risen again to 3.16·10⁻¹². The right-hand rise is the squaring phase amplifying its own rounding, which is why a library chooses s from ‖A‖ rather than taking as many as it can afford.0246810121410⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸squarings, srelative error against the closed formbest: s = 2the approximant is out of rangethe squarings amplifyboth sides costerror at s = 05.6·10⁻¹¹best, at s = 29.1·10⁻¹⁶error at s = 153.2·10⁻¹²‖A‖ / 2^s at the best3the squarings are exact in the algebraand each one doubles the rounding
Fig. 23 And the cost of getting its argument small enough first.
Four knobs on one problem at 1.0% noiseRelative error against the fraction of each method's own range, on a logarithmic vertical axis. A truncation, a Tikhonov parameter, a conjugate gradient step count and a randomised rank each have an interior minimum, and the four minima are 0.1445, 0.1406, 0.1426, 0.1449 — a spread of 3%. The horizontal line is the best of them.00.250.50.75110⁻¹110¹10²10³fraction of the method's own rangerelative errorfloor 0.141truncation KTikhonov λCGLS steprandomised rankfour methods, one floortruncation K0.14Tikhonov λ0.14CGLS step0.14randomised rank0.14four knobs from four fieldsand one obstruction underneath them
Fig. 24 Four rules landing together for the opposite reason.
The L-curve, and where four rules put λThe norm of the solution against the norm of its residual, on logarithmic axes, as λ sweeps eight decades. The curve has a corner: to the left of it the noise is being amplified and to the right the signal is being thrown away. Four points are marked — the three rules that use only the data, and the oracle, which requires the exact answer and is not a method.10⁻¹110¹10²10³10⁴10⁵‖Ax − b‖‖x‖the oraclediscrepancyL-curvegeneralisedscored against a truth none hasoracle, relative error0.14discrepancy principle, as a multiple1L-curve corner, as a multiple1.3generalised cross-validation, as a multiple1the oracle needs the exact answer and is not a methodit is the reference the others are scored on
Fig. 25 And what each of them is reading.
Two bases for one subspace, at an interpolation point of 20κ₂ of the moment basis (σI − A)⁻¹B, (σI − A)⁻²B, … at a single point σ = 20, against κ₂ of the same number of solves spent at points spread across the model's poles and orthogonalised as they are built. The moment basis runs 10.9 to 7.75·10⁹, multiplying by about 65.8 per vector, while the spread one stays at 1. The reason is the one behind every Krylov basis on this site: repeated application of one operator converges to a direction, so every new moment is nearly the last. The rate is set by the ratio of distances from σ to the poles, which means the user's own choice of interpolation point decides how fast their basis becomes unusable — and spending the solves at different points instead is both a better approximation and a basis that survives.2345678110²10⁴10⁶10⁸10¹⁰vectors in the basisκ₂ of the basisas derivedsolves spread outone subspace, two spanning setseight moments at one point7.7·10⁹eight points, spread1growth per vector661/u4.5·10¹⁵the same subspaceand only one of them usable
Fig. 26 Points placed badly in a different construction in the same field.
A reduced model of order 4, and the 4 places it is exact|H(s) − Hᵣ(s)| ÷ |H(s)| along the real axis, for a rational-Krylov reduction of a 24-state model at the interpolation points 1.5, 5, 16, 50. At each of them the curve falls to 5.3·10⁻¹⁶ — the reduced function passes through the original, and its derivative does too, because the projection is two-sided. Between and beyond them it reaches 8.37·10⁻⁴, and there is no bound on it: the method buys 8 exact conditions for 4 solves and offers nothing anywhere else. That is the trade against balanced truncation, which asks for nothing and bounds everything at a cost of two Lyapunov solves.10⁻¹110¹10²10³10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵s, on the real axis|H − Hᵣ| ÷ |H|exact where askedpoints4conditions bought8worst at a point5.3·10⁻¹⁶worst away from one8.4·10⁻⁴4 points, 8 conditionsand no bound in between
Fig. 27 And placed well.
What a rank-10 approximation can achieve, by spectrumA semi-logarithmic plot of singular value against index for three spectra — geometric decay, algebraic decay, and flat — with the rank-ten approximation error marked on each.1112131415110⁻³10⁻²10⁻¹1index jσⱼσ11geometric, 0.85ʲalgebraic, j^−0.00flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic1rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 28 A spectrum where no placement helps.
Relative error of each singular value of a bidiagonal matrix graded over 11 decades, by four routesThe matrix is 8×8, bidiagonal, and every entry is a power of two times a small integer, so the matrix the float routines are handed and the matrix the exact routine reasons about are the same matrix. The reference is a Sturm bisection in BigInt rationals, which never rounds. σ_min is 1.934·10⁻¹¹. One-sided Jacobi holds every value to 7.1·10⁻¹⁶ relative and the zero-shift bidiagonal sweep to 4.41·10⁻¹⁶; the eigenvalues of BᵀB are wrong by 14.5 — not a rounding error but the whole number — and the damage is entirely at the small end.1234567810⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, largest firstrelative errorone-sided Jacobizero-shift QRshifted QReigenvalues of BᵀBagainst a rational bisectionσ_min, exactly1.9·10⁻¹¹worst, one-sided Jacobi7.1·10⁻¹⁶worst, zero-shift QR4.4·10⁻¹⁶worst, eigenvalues of BᵀB15a relative error is a ratioand the denominator is the answer
Fig. 29 And one spanning ten decades, where the rate is slow and still geometric.
The condition number a sketched preconditioner leaves, against the condition number it was givenTwo curves against κ(A), both axes logarithmic. The matrix's own condition number climbs the diagonal from 100 to 10¹⁰; κ(AR⁻¹), where R comes from a QR of a 4n-row sketch, is 2.6304 at every one of them — the same number to ten digits, not a similar one. The reason is two lines of algebra: with G = SU the preconditioned singular values are those of (GᵀG)⁻¹, which has no spectrum of A in it at all.10²10⁴10⁶10⁸10¹⁰110²10⁴10⁶10⁸10¹⁰condition number of the matrixcondition number seen by the iterationκ(A), unchangedκ(AR⁻¹)a bound with no κ(A) in itκ(AR⁻¹), every κ(A)2.6κ(SU), the other route2.6κ(A), across the sweep10⁸κ(AR⁻¹), across the sweep1the sketch never sees the spectrumand the spectrum cancels out of the answer
Fig. 30 Changing a spectrum on purpose, which is the other way to make κ smaller.
The two 3×3 matrices of a Sylvester equation, and the 9×9 matrix it meansA and B are 3×3 with 5 nonzeros each. The coefficient matrix of the map X ↦ AX + XB is 9×9 with 21 nonzeros, drawn at the same scale so the ratio is visible rather than quoted. Bartels and Stewart's algorithm never forms it: two Schur reductions and a back-substitution give the same X to 5.5·10⁻¹⁷, and that X satisfies AX + XB = C to 8.8·10⁻¹⁷.A3×3B3×3I ⊗ A + Bᵀ ⊗ I9×9one equation, two objectsentries in A and B18entries in the coefficient matrix81two routes, relative gap5.5·10⁻¹⁷‖AX + XB − C‖/‖C‖8.8·10⁻¹⁷the small squares are the problemand the large one is the notation
Fig. 31 The equation being solved.
The rational function 6 ADI shifts amount to, chosen two ways|∏(λ − pⱼ)/(λ + pⱼ)| over the spectrum of −A for a 30-point discretisation, whose ends are 9.86 and 3834 and whose condition number is 388.8. The lower curve uses 6 geometrically spaced shifts and touches zero at each of them, equioscillating between with a maximum of 0.1187. The upper curve uses 6 equally spaced ones: every shift sits in the top decade of the spectrum, the small end is never approached, and the maximum is 0.8922 — a factor of 7.52 worse, which the iteration pays squared. The two runs are the same algorithm, the same number of solves, and one different line where the shifts are chosen.10¹10²10³10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹λ, over the spectrum of −A|r(λ)|geometric maxequally spacedgeometricone line of codeshifts6κ of the spectrum389geometric max0.12equally spaced max0.89the factor between7.5the same k solvesand one choice of where
Fig. 32 The hero again, beside the decay it explains.

Where this goes next

1234567891010⁻¹⁹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹kλₖ₊₁ ÷ λ₁, and the boundZ_k²the Gramianpredicted from two numbersstates30κ of the spectrum389λ11 ÷ λ₁2.5·10⁻⁸the bound there5.2·10⁻⁴the cliff everything rests onand the reason for it Why a Gramian can be truncated at all Every method in this field rests on one fact nobody states the reason for — the eigenvalues of a Gramian fall off a cliff. The equation defining it has a rank-one right-hand side and no low-rank structure anywhere — and the answer's decay is a rational approximation problem with a closed-form rate. 0246810⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹log₁₀ γ, the change of unitsforward erroras writtenafter scalingone change of variableunscaled, worst0.0013scaled, worst1.7·10⁻¹³orders recovered10scaled coefficient spread4.5the answer was never the problemthe units were The scaling that buys ten orders Two lines computed from three norms, a change of variable that is exact in both directions, and the whole of the loss the previous essay measured comes back — flat, at every stop, because after scaling every stop is the same problem. 081624324010⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²conjugate gradient steprelative residualabove, on their own scale: rank asked 5, rank kept 3solid: the budgeted residual · light: no budgetbudget 3rank asked for5rank kept3residual, budgeted3.8·10⁻⁴residual, unbudgeted1.6·10⁻¹⁴numbers stored150the step asks for moreat every step An iterate that must be made smaller Applying a Kronecker-sum operator to a low-rank iterate multiplies its ranks by d and adding two of them adds their ranks, so a solver in a compressed format cannot keep what it produces. Every step is followed by a truncation — and whether that truncation is a floor on the residual depends on the right-hand side rather than on the truncation.

What links here

Computed from the collection, not written here: the essays that point at this one.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

ADI iterationCondition numberEquioscillationGramianLow-rank approximationLyapunov equationRational approximationShift selectionZolotarev number