Numerical Linear Algebra
where the algebra stops being the arithmetic
When the problem arrives again

Where the drift lands

The standing rule for when a preconditioner has gone stale is to rebuild it once the matrix has changed by more than some fraction of itself. Two drifts of exactly the same relative size cost 19 iterations and 5 on the same matrix, and the quantity that separates them is not in the rule at all — the perturbation is divided by the eigenvalue it lands on.
28 min read 24 figures Where the algebra stops being the arithmeticWhose fault is itMeasured, not assertedTwo routes to a number

Worth reading first: A factorisation kept past its date · Changing the condition number on purpose · The condition number is an amplifier.

A preconditioner is the expensive object in an iterative solve and the obvious one to keep. Every code that solves a sequence keeps one, every one of them has a rule for when to throw it away, and the rule is almost always the same: rebuild when the matrix has changed by more than some fraction of itself.

Written down, that is a threshold on ‖A_t − A₀‖ / ‖A₀‖. It is a perfectly reasonable-looking quantity — it is relative, it is cheap, and it goes to zero when nothing has happened.

It is also the wrong quantity, and the size of the mistake is the condition number.

Preconditioned iterations against the size of the drift, for a drift on each half of the spectrumThe preconditioner is the exact Cholesky factorisation of the first matrix, so it solves that matrix in one step and every iteration counted here is bought by the drift. Both curves are drifts of the same relative Frobenius norm; the upper one sits on the half of the spectrum with the small eigenvalues and the lower on the half with the large. At a relative drift of 10⁻² they cost 38 and 9 iterations, because the drift is divided by the eigenvalue it lands on: the preconditioned matrix has eigenvalues 1 + w/λ, and λ runs over 10⁴.10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10112233445566relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²38large end at 10⁻²9κ(M⁻¹A), small end36κ(M⁻¹A), large end1.3‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about
Fig. 1 Preconditioned iterations against the size of the drift, for two drifts of identical relative Frobenius norm placed on the two halves of one spectrum.

The experiment

The preconditioner is the exact Cholesky factorisation of the first matrix, which is not what anybody uses and is exactly what this measurement needs: with it, M⁻¹A₀ is the identity to the last bit and conjugate gradients solves the first member in a single step. So every iteration counted afterwards was bought by the drift and by nothing else.

The drift is then built in the spectrum rather than added and measured afterwards. A₀ is constructed with a known eigendecomposition at a prescribed condition number, and the perturbation is

E = s · Σ v_i v_iᵀ

summed over the forty smallest eigenvalues or the forty largest, with s chosen so that the Frobenius norm of E is the same fraction of ‖A₀‖ in both cases. The two drifts are the same size in the norm the rule is written in, to fourteen digits, and that equality is asserted rather than assumed.

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.2284 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.2κ(SU), the other route2.2κ(A), across the sweep10⁸κ(AR⁻¹), across the sweep1the sketch never sees the spectrumand the spectrum cancels out of the answer
Fig. 2 What a preconditioner is worth when it is built for the matrix it is used on, from the field that priced it. Everything here is measured against that baseline of one step.

The measurement

‖E‖ / ‖A₀‖ on the small eigenvalues on the large ones
10⁻⁶ 4 iterations 2
10⁻⁵ 5 3
10⁻⁴ 10 4
10⁻³ 19 5
10⁻² 38 9
10⁻¹ 55 19
1 58 38

At a relative drift of 10⁻³ — a tenth of a per cent, a number most rules would not even react to — the two runs cost 19 iterations and 5. At 10⁻² they cost 38 and 9. The rule cannot distinguish them, because in the rule’s own quantity they are identical.

Preconditioned iterations against the size of the drift, for a drift on each half of the spectrumThe preconditioner is the exact Cholesky factorisation of the first matrix, so it solves that matrix in one step and every iteration counted here is bought by the drift. Both curves are drifts of the same relative Frobenius norm; the upper one sits on the half of the spectrum with the small eigenvalues and the lower on the half with the large. At a relative drift of 10⁻² they cost 82 and 16 iterations, because the drift is divided by the eigenvalue it lands on: the preconditioned matrix has eigenvalues 1 + w/λ, and λ runs over 10⁶.10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹101632486480relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²82large end at 10⁻²16κ(M⁻¹A), small end2911κ(M⁻¹A), large end3.7‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about
Fig. 3 At a condition number of 10⁶ the gap between the two curves opens further, because the gap is the condition number.
Preconditioned iterations against the size of the drift, for a drift on each half of the spectrumThe preconditioner is the exact Cholesky factorisation of the first matrix, so it solves that matrix in one step and every iteration counted here is bought by the drift. Both curves are drifts of the same relative Frobenius norm; the upper one sits on the half of the spectrum with the small eigenvalues and the lower on the half with the large. At a relative drift of 10⁻² they cost 10 and 6 iterations, because the drift is divided by the eigenvalue it lands on: the preconditioned matrix has eigenvalues 1 + w/λ, and λ runs over 100.10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹1061218243036relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²10large end at 10⁻²6κ(M⁻¹A), small end1.5κ(M⁻¹A), large end1‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about
Fig. 4 And at 10² it nearly closes. A rule written in the norm is approximately right on a well-conditioned problem and wrong by orders of magnitude on the problems preconditioners exist for.

Why, which is one line

With M the exact factorisation of A₀ and A_t = A₀ + E,

M⁻¹A_t = I + M⁻¹E

and if E is built from A₀’s own eigenvectors with weight w_i on the ith, the eigenvalues of the preconditioned matrix are exactly

1 + w_i / λ_i .

The drift is divided by the eigenvalue it lands on. A perturbation sitting on λ_min is amplified by 1/λ_min; the same perturbation on λ_max is divided by λ_max; and ‖E‖ is the same number in both cases because a Frobenius norm does not know where in the spectrum its mass is.

So the quantity that predicts the cost is max(1 + w_i/λ_i) over min(1 + w_i/λ_i), and at a relative drift of 10⁻³ it is 4.47 on the small end and 1.0327 on the large one. Both are available before either solve starts, from the construction rather than from an eigensolve of the product.

How much a perturbation of the right-hand side is amplified, κ = 10⁴The cumulative distribution of the amplification factor over two hundred random perturbation directions, with the condition number marked as the upper limit.110¹10²10³10⁴10⁵00.250.50.751amplification of the input perturbationfraction of directions at or belowκ = 10000worst found 76206×6, 200 directionsmedian reaches 0.29 of κ
Fig. 5 The amplifier the collection is organised around, appearing here in a place it is not usually looked for: not amplifying the data’s error into the answer’s, but amplifying a perturbation’s position into a cost.
Two condition numbers of one 12×12 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 15 to 3.3·10⁶ while the componentwise condition number stays at 13.2 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.8·10¹⁶ and 1.1·10¹⁶: a matrix whose sensitivity no scaling repairs.0123456110³10⁶10⁹10¹²10¹⁵10¹⁸spread of the row units (decades)condition numberκ_∞(DA)cond(DA)Hilbert κ_∞Hilbert condone system, two numbersκ_∞ at no spread15κ_∞ at 6 decades3.3·10⁶cond, either end13Hilbert, equilibrated1.3·10¹⁶the solution is the same at every spreadand one of these curves knows it
Fig. 6 And the essay on the units a matrix is measured in, whose point is the same one from the other side: a norm is a summary, and two matrices a norm calls equally far apart need not be.

The predictor orders and over-predicts

The predictor gets the ordering right at every drift measured, which is the finding. What it does not get right is the size, and saying so is the second half of an honest measurement.

The conjugate gradient bound applied to a preconditioned condition number of 3.47·10³ asks for 699 iterations. The run takes 58. The bound over-predicts by a factor of twelve, and the fitted slope of iterations against predicted conditioning is 0.176 rather than the bound’s 0.5.

The reason is the one this collection already has an essay about: a stale preconditioner does not spread the spectrum, it splits it. Half of M⁻¹A_t’s eigenvalues are exactly 1, because the drift was placed on the other half, and a method that minimises a polynomial over the spectrum spends nothing on a cluster it has already annihilated. The condition number is a summary of two numbers and the method is reading the whole distribution.

So: where the drift lands decides the cost, and the standard bound built on the same quantity is loose by an order of magnitude. Both statements are about the same figure and neither is a correction to the other.

The spectrum of C⁻¹T against T's own, n = 64, ρ = 0.8Two sorted spectra on a linear vertical axis. The unpreconditioned Toeplitz matrix spreads its eigenvalues across its whole range; the preconditioned operator puts 97% of them within a tenth of 1, with 0 below zero.08162432404856640123456789index, sortedeigenvalueT aloneC⁻¹Teigenvalues, not singular valuesfraction within 0.1 of 10.97eigenvalues below zero0λ_min of the preconditioner0.11the dashed line is 1, where the cluster formsthe solid line is zero
Fig. 7 Why the bound over-predicts, from the essay that established it: a clustered spectrum is solved far faster than its condition number allows.
GMRES on the Laplacian and on the cyclic shift, both 40×40A semi-logarithmic plot of relative residual against iteration. One curve falls steadily; the other is flat at one for every step until the last, where it drops to zero.071421283510⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1iteration‖r‖ / ‖b‖Laplaciancyclic shiftno progress at allevery eigenvalue of the shift is on the unit circleand it predicts nothing
Fig. 8 And the standing warning that the spectrum does not predict everything either. The predictor here is an ordering, not a rate.

The object a real code holds

The exact Cholesky above isolates the phenomenon. Nobody keeps one. What a code keeps is an incomplete factorisation, and it ages differently in a way worth drawing.

The drift is the one this collection already has a field about: the five-point operator’s two directions stop being equally weighted, ε running from 1 down to 0.033 over nineteen members, with the sparsity pattern unchanged — so the incomplete factorisation from the first member stays applicable for the whole run, which is exactly the situation in which it gets kept.

The kept factorisation goes

18 21 24 26 29 32 34 36 37 38 38 38 40 42 44 46 48 50 52

iterations. A factorisation rebuilt at every member goes

18 18 18 18 18 17 17 17 17 16 16 16 15 15 15 15 14 14 13

in the opposite direction. The problem is getting easier — an anisotropic operator is a friendlier problem for a factorisation that knows about the anisotropy — and the kept preconditioner is getting worse at it. The two start at the same 18 by construction and never meet again.

An incomplete Cholesky kept while the operator turns anisotropic, against one rebuilt at every memberThe five-point operator's two directions stop being equally weighted, ε running from 1 down to 0.0331, with the sparsity pattern unchanged throughout — so the factorisation from the first member stays applicable for the whole run, which is the situation in which it gets kept. The kept one goes from 18 iterations to 52. The rebuilt one goes from 18 to 13, because an anisotropic operator is an easier problem for a factorisation that knows about the anisotropy. The two start at the same point by construction and never meet again.0246810121416182001020304050member of the sequencepreconditioned conjugate gradient iterationskept from the first memberrebuilt every memberone operator, two policieskept, first member18kept, last member52rebuilt, first member18rebuilt, last member13ε at the last member0.033the problem got easierand the kept preconditioner got worse at it
Fig. 9 The two curves. They agree exactly at the member the factorisation was built from, which is the only point on the figure that is not a measurement of drift.
An incomplete Cholesky kept while the operator turns anisotropic, against one rebuilt at every memberThe five-point operator's two directions stop being equally weighted, ε running from 1 down to 0.0331, with the sparsity pattern unchanged throughout — so the factorisation from the first member stays applicable for the whole run, which is the situation in which it gets kept. The kept one goes from 23 iterations to 55. The rebuilt one goes from 23 to 17, because an anisotropic operator is an easier problem for a factorisation that knows about the anisotropy. The two start at the same point by construction and never meet again.024681012141618200102030405060member of the sequencepreconditioned conjugate gradient iterationskept from the first memberrebuilt every memberone operator, two policieskept, first member23kept, last member55rebuilt, first member23rebuilt, last member17ε at the last member0.033the problem got easierand the kept preconditioner got worse at it
Fig. 10 The same comparison on a larger grid, where both curves lift and the separation between them does not change.

The failure mode nobody sees

The incomplete case is the one worth worrying about in practice, and not because it is worse. It is milder, and that is the problem.

A stale factorisation stops converging — a cliff, an error, something in a log. A stale preconditioner never stops converging. It costs 18 iterations, then 30, then 52, and every member returns a correct answer to the tolerance it was asked for. Nothing in the output says that the run is now taking three times as long as it needs to, because nothing in the run knows what it would have taken with a fresh preconditioner.

That is a class of defect this collection keeps meeting: the symptom is absence. There is no wrong number to find, no assertion to trip, and no gate that could be written on the output of a single member. The only way to see it is to run the comparison, and running the comparison means building the fresh preconditioner, which is the thing being avoided.

LSQR steps to a normal residual of 10⁻¹², with and without a sketched preconditionerTwo curves against κ(A). The unpreconditioned solve takes 15, 22, 32, 59, 150 steps as the conditioning worsens, and at the worst one it does not reach the tolerance inside 150. The preconditioned solve takes 10, 10, 10, 10, 10 — flat across eight decades. Both return the same vector, to the accuracy the problem allows: the preconditioned answer and the direct one differ by about κ(A)·u at every κ, 5.7·10⁻⁶ at the worst. The SAME sketch used as an answer instead of as a preconditioner is 4.5 away.10²10⁴10⁶10⁸10¹⁰110¹10²10³condition number of the matrixLSQR steps to 10⁻¹²no preconditionersketchedthe cost, not the answersteps, preconditioned10steps, without150answer gap, preconditioned5.7·10⁻⁶κ(A)·u at the worst κ1.1·10⁻⁶answer gap, sketch-and-solve4.5one sketch, two things to do with itand only one of them moves the answer
Fig. 11 The iteration count as a currency, from the field that established it. It is the only thing that moves when a preconditioner ages, and it moves quietly.
The cost of a drifting sequence of nonlinear solves, against how often the Jacobian is refactorisedTwenty members, each warm-started from the last, each solved by a chord iteration on a factorisation that may be several members old. Refactorising at every member costs 16.42 MFlop; refactorising every 5 costs 9.27. Past a period of 20 the chord iteration stops converging altogether, which is drawn as an open circle on the ceiling rather than omitted. The filled square is the rule that refactorises when the observed contraction ratio exceeds 0.2: 4 factorisations, 9.22 MFlop, and it was never told the drift rate.024681012141618202210⁶10⁷10⁸members served by one factorisationmultiplications for the whole sequencedoes not convergethe contraction ruledrift 0.01 a memberevery member1.6·10⁷every 5 members9.3·10⁶contraction rule9.2·10⁶its factorisations4cliff at a period of20a factorisation has a shelf lifeand the cliff is past the optimum
Fig. 12 The other object, which fails loudly. The two failure modes want different rules and it is worth not assuming otherwise.

Two routes to the same number

The predicted conditioning is computed here from the construction — the weights w_i were chosen, the eigenvalues λ_i are known, and 1 + w_i/λ_i is arithmetic on two lists. That is the cheap route and it is only available because the experiment built the drift on purpose.

It would be worth nothing if it did not agree with the expensive route, so the agreement is the check. Forming M⁻¹A_t explicitly and taking the ratio of its extreme eigenvalues gives the same number, and it must: the two are the same quantity computed from two different objects, one of them never assembled. The collection’s standing habit is that a number arrived at by one route has been wrong every time it has been published, and this is the cheapest available second route.

What the agreement licenses is the sentence in the figure’s badge: the conditioning of the preconditioned matrix is knowable before the solve, from quantities a code that built the preconditioner already has. It is not knowable cheaply for a general drift — nobody has A₀’s eigenvectors — but the experiment is not about what is cheap. It is about which quantity the answer depends on, and the answer depends on that one and not on the norm.

How far a perturbation of size ε moves an eigenvalue, 10×10A log–log plot of eigenvalue movement against perturbation size. The symmetric case lies on a line of slope one; the non-symmetric case lies on a line of slope one eighth, and at a perturbation of ten to the minus sixteen it has already moved by a hundredth.10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹size of the perturbation ‖δA‖how far the eigenvalues moveJordan block, ε^(1/10)symmetric, ≤ ‖δA‖rounding error alone moves it to 10⁻²six seeds per symmetric point; Jordan is closed formsymmetry beats precision
Fig. 13 Eigenvalues moving under a perturbation, from the field that measures how far they move. Here the same perturbation is being read for what it does to a solver rather than to a spectrum.
Where 24 perturbations of size 10⁻⁸ put the eigenvalues of two 6×6 matrices with the same spectrumTwo clouds of eigenvalues in the complex plane. Both matrices have every eigenvalue at 0.8 exactly; both were perturbed by the same 24 random matrices of norm 10⁻⁸. The normal matrix's eigenvalues stay within 8.86·10⁻⁹ of where they were — the size of the perturbation — and the bidiagonal's spread out to 0.0701, a factor of 7.9·10⁶ further on the same data.100.1750.35real partimaginary part0.8bidiagonalnormalone spectrum, two matricesthe perturbation10⁻⁸normal, furthest moved8.9·10⁻⁹bidiagonal, furthest0.07ratio7.9·10⁶the two matrices have identical eigenvaluesand one of them says so under perturbation
Fig. 14 And the same question asked of a matrix with no symmetry, where the movement is far larger. Every statement in this essay is about a symmetric positive definite matrix and does not survive that change.

What the rule should be written in

If not ‖E‖/‖A₀‖, then what?

The honest answer from this measurement is ‖A₀⁻¹E‖, or anything that behaves like it: the drift measured in the units the preconditioner works in rather than the units the matrix is written in. That quantity is 4.47 and 0.033 for the two drifts above, which is a factor of 135 between two things the naive rule calls equal.

Estimating it is not free — it is a norm of a product with an inverse, which is the object the condition-estimation essays are about — but it is not expensive either, and this collection already has the machinery: an estimate of ‖A₀⁻¹E‖ costs a handful of solves with a factorisation that is already in memory.

The cheaper answer, and the one the essay on what a rebuild is worth is about, is to stop trying to predict and read the iteration count instead. That has its own failure and the two essays are best read together.

How close Hager's estimate is to the true κ₁, over 200 seeded 8×8 matricesFive bars. The estimate is exactly the true condition number on 81% of the sample and inside ten per cent on 87%; the worst underestimate in the whole sample returns 38% of the truth. The last bar is the matrix built to defeat it, at 7.7%, well below anything the random sample reached.each bar is a percentage — of the sample, or of the true condition numberexactly right80.5%inside 10%87.0%inside a factor of 287.0%worst in the sample, ×10037.7%the constructed matrix, ×1007.7%usually exactexact share0.81worst of the sample0.38the constructed matrix0.077a routine that is right most of the timeand never wrong in the safe direction
Fig. 15 Estimating a norm that involves an inverse, from the field that measured how well it can be done. The quantity the rule wants is of exactly this kind.
The cheapest rebuild period for one drifting sequence, against what a rebuild costsOne sequence, one drift rate, one preconditioner — and six answers, because the answer is not a property of the sequence. Costed in iteration-equivalents, the optimal period runs from every 1 member at a setup worth 5 iterations to every 12 at a setup worth 200. A dense Cholesky at this size is worth 6.7 iterations, which is at the left of the axis, so for it the answer is always to rebuild. The rule that rebuilds when a solve takes 1.5 times what the last fresh one did beats the best fixed period at 1 of the 6 ratios and loses at the rest, because it cannot see the setup cost at all.10¹10²10³03691215what one rebuild is worth, in preconditioned iterationscheapest number of members between rebuildsa dense Cholesky here is worth 6.7 iterationsringed: where the growth rule beats every fixed periodcheapest period, by setup costsetup worth 51setup worth 102setup worth 203setup worth 506setup worth 1008setup worth 20012how long to keep itis a question about what it cost to build
Fig. 16 And the other approach: not predicting the ageing at all, and reading the cost of a solve that has already happened.

What would have to be true for the norm rule to work

It is worth stating the condition under which the standard rule of thumb is right, because it is not empty and because naming it is more useful than saying the rule is wrong.

A threshold on ‖E‖/‖A₀‖ predicts the cost correctly when the drift is spread evenly over the spectrum — when the perturbation has about as much weight on the small eigenvalues as anywhere else. Then max(1 + w_i/λ_i) is proportional to ‖E‖/λ_min, the condition number enters once as a constant, and a threshold calibrated on one problem transfers to another of similar conditioning.

That is a real situation. A drift that is genuinely random with respect to the operator’s eigenvectors — a perturbation from measurement noise, say, or from a rounding — looks like that, and for it the rule is a reasonable proxy.

What is not like that is any drift with structure, which is most of them: a time step that changes one coefficient, a load that acts on part of a domain, an operator whose anisotropy grows. Those have their weight somewhere in particular, and where they have it is exactly the question the norm discards. So the rule is not a bad approximation to the right quantity; it is the right quantity multiplied by an unknown that runs from 1 to the condition number, and nothing in the norm says which.

The refusal

The claim under test is that a relative drift threshold is a usable rule. It is fed the case where it happens to be right, and the case is chosen to be one where being right means nothing.

At a relative drift of 10⁻⁶, the small end costs 4 iterations and the large end costs 2. The assertion that the small end costs more than two and a half times the large one is fed those numbers and fails — correctly, because at that drift neither placement costs anything worth measuring and there is no distinction to make. The threshold rule is right there, in the sense that it makes no error; it is right the way a rule that always says “do nothing” is right on a problem where nothing needs doing.

That is what the refusal is for: to stop the main measurement from being read as a claim that the norm is always wrong. It is wrong where the answer matters and correct where it does not, which is the least useful pattern a rule can have.

Preconditioned iterations against the size of the drift, for a drift on each half of the spectrumThe preconditioner is the exact Cholesky factorisation of the first matrix, so it solves that matrix in one step and every iteration counted here is bought by the drift. Both curves are drifts of the same relative Frobenius norm; the upper one sits on the half of the spectrum with the small eigenvalues and the lower on the half with the large. At a relative drift of 10⁻² they cost 21 and 7 iterations, because the drift is divided by the eigenvalue it lands on: the preconditioned matrix has eigenvalues 1 + w/λ, and λ runs over 1000.10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10918273645relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²21large end at 10⁻²7κ(M⁻¹A), small end4.9κ(M⁻¹A), large end1.1‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about
Fig. 17 Between the two extremes, where the ordering is already clear and the gap is a factor of three.
Preconditioned iterations against the size of the drift, for a drift on each half of the spectrumThe preconditioner is the exact Cholesky factorisation of the first matrix, so it solves that matrix in one step and every iteration counted here is bought by the drift. Both curves are drifts of the same relative Frobenius norm; the upper one sits on the half of the spectrum with the small eigenvalues and the lower on the half with the large. At a relative drift of 10⁻² they cost 64 and 12 iterations, because the drift is divided by the eigenvalue it lands on: the preconditioned matrix has eigenvalues 1 + w/λ, and λ runs over 10⁵.10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10132639526578relative drift ‖E‖ / ‖A₀‖preconditioned conjugate gradient iterationson the small eigenvalueson the large onessame drift, two placessmall end at 10⁻²64large end at 10⁻²12κ(M⁻¹A), small end315κ(M⁻¹A), large end1.9‖E‖/‖A₀‖, both0.01how much the matrix changedis not what the preconditioner cares about
Fig. 18 And near the other end, where it is a factor of six and the norm still cannot see it.
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.2284 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.2κ(SU), the other route2.2κ(A), across the sweep10⁸κ(AR⁻¹), across the sweep1the sketch never sees the spectrumand the spectrum cancels out of the answer
Fig. 19 The preconditioned condition number when the preconditioner is fresh, which is the number the whole comparison is against.
Conjugate gradients at κ = 104, against the bound κ permitsA semi-logarithmic plot of the relative A-norm error against iteration count. The measured curve falls below a smooth dashed curve showing the classical condition-number bound.0408012016020024010⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹iteration‖e‖_A / ‖e₀‖_Ameasuredκ bound119 steps40×40, spectrum spread evenly in logbound permits 1417
Fig. 20 And the rate that condition number predicts — the bound this measurement finds twelve times too pessimistic once the spectrum is in clusters.
What 12 members of a sequence cost, by what changes between themEvery bar is 12 solves of an 120×120 system, costed in multiplications with a dense Cholesky at n³/3 and a solve at 4n². When nothing changes the answer is reused and the cost is one factorisation. When only the right-hand side changes, one factorisation serves every member — 16.5 per cent of the independent cost, with every solve still at a relative residual of 1.54·10⁻¹³. When the matrix drifts, the factorisation is reusable for a while and the while has to be measured: 4 factorisations for 12 members. When everything changes, nothing carries.nothing8.3%the right-hand side16.5%the matrix, slowly85.9%everything100.0%what can be reused: the answerwhat can be reused: the factorisationwhat can be reused: the factorisation, for a whilewhat can be reused: nothingshare of the cost of a sequence that shares nothingwhat 12 members cost in factorisationssame: factorisations1rhs: factorisations1drift: factorisations4independent: factorisations12what changes between the membersdecides what may be carried
Fig. 21 Where this sits in the field’s ledger: the same drifting bar, priced for a different object.
Chord steps per member of a drifting sequence, with the factorisation rebuilt every 4Each point is one member of a continuation, warm-started from the previous member's answer and solved with a factorisation that is between zero and 3 members old. The lower series is the same sequence with a fresh factorisation at every member, which costs 6 to 5 steps throughout. The kept run costs 6 steps on a fresh factor and 12 at its worst, for 5 factorisations against 20.03691215182103691215member of the sequencechord steps to reach the tolerancea fresh factorisation each memberthe ringed points are the refactorisationsrebuilt every 4factorisations5chord steps160worst member's steps12multiplications1.2·10⁷with a fresh factor each time1.8·10⁷a factor is cheapest the member it was built forand dearest the member before it is replaced
Fig. 22 And the object whose ageing is loud, for contrast with the one whose ageing is silent.

One more thing the norm cannot see

There is a second quantity hiding in the incomplete-factorisation figure, and it is worth naming because it is the same mistake in a different disguise.

The drift there is an anisotropy: the operator’s two directions stop being equally weighted. In the Frobenius norm that is a large change — by the last member the matrix is a substantial distance from the first — and a norm-based rule would have rebuilt many members earlier. It would have been right to, but not for the reason it thought: the rule would have fired because the matrix moved, and the correct reason to fire is that the smoothing directions the incomplete factorisation encodes have stopped being the directions the operator is stiff in.

This collection has a whole field about that distinction. Line relaxation along the strong direction and coarsening only along it are both statements about where in the problem the difficulty is, and neither of them is a statement about how far the matrix has moved. A preconditioner is a claim about the structure of the difficulty, and it goes stale when the structure moves rather than when the entries do.

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. 23 The structure a preconditioner is a claim about, from the field that measures it: what a smoother actually reduces, direction by direction.
The same aggregates, with and without the prolongator smoothing (ε = 0.01)Relative residual against V-cycle on a logarithmic vertical axis. Both hierarchies are built from the identical aggregation of the same matrix. Without the smoothing sweep the convergence factor is 0.7051; with it the factor is 0.1923, reached in 13 cycles. The smoothed hierarchy stores 2.71 times the fine matrix against 1.51.036912151810⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹V-cyclerelative residualpiecewise constantsmoothedone sweep on the columns of Pfactor, unsmoothed0.71factor, smoothed0.19operator complexity, smoothed2.7the same aggregates in bothand one sweep between them
Fig. 24 And the same operator’s difficulty, which has a direction in it. A norm has no directions at all, which is the whole of this essay in one sentence.

Where this goes next

024681012141618202210⁶10⁷10⁸members served by one factorisationmultiplications for the whole sequencedoes not convergethe contraction ruledrift 0.01 a memberevery member1.6·10⁷every 5 members9.3·10⁶contraction rule9.2·10⁶its factorisations4cliff at a period of20a factorisation has a shelf lifeand the cliff is past the optimum A factorisation kept past its date One Cholesky factor can serve five members of a drifting sequence and save 44 per cent of the work. Kept for twenty it does not lose accuracy — it stops converging altogether. The optimum and the cliff are four members apart, both move with the drift, and a rule written in a ratio the iteration has already computed finds them without being told what the drift is. 051015202530354010⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹iteration‖r‖ / ‖b‖plain CGIC(0) CGwhat the preconditioner didκ(A)48κ(L⁻¹AL⁻ᵀ)5.1‖A − LLᵀ‖/‖A‖0.0832D Laplacian, n = 100√κ ratio predicts 3.07× Changing the condition number on purpose Preconditioning is usually introduced as a trick that makes an iteration converge faster. It is not a trick. It is solving a different system with the same solution and a condition number chosen rather than inherited, and the new condition number is computable. 10⁻¹⁴10⁻¹³10⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1relative size of the entrywise perturbationrelative forward errorκ_∞ · εcond(A,x) · εmeasuredboth bounds holdκ_∞(A)1.9·10⁸cond(A, x)4.8ratio of the bounds4·10⁷both curves above the data are boundsand only one of them is a measurement A condition number scaling cannot move Skeel's componentwise condition number is invariant under any row scaling — exactly, before any norm is taken, because two diagonal factors cancel entry by entry. It is never larger than the normwise one and can be arbitrarily smaller, and the ratio between them is a diagnostic for which kind of ill-conditioning a matrix has.

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.

AnisotropyCholesky factorisationCondition numberConjugate gradientsEigenvalueFlop countIncomplete factorisationPreconditioningSpectral clustering