Numerical Linear Algebra
where the algebra stops being the arithmetic
Neither sparse nor dense

The accuracy worth paying for

Used as a preconditioner, a hierarchical representation gets better at every accuracy — the iteration count falls monotonically all the way to the tightest tolerance. The total work does not. Its minimum sits at a rank-one preconditioner on an easy problem and six decades further along on a hard one.
33 min read 25 figures Measured, not assertedWhose fault is itPrecision as the knob

Worth reading first: The exact answer to a nearby problem · Changing the condition number on purpose · Which pairs are allowed to be small.

The previous essay established that compressing a matrix is perturbing the problem. A perturbed problem is what a preconditioner is, so the obvious thing to do with a deliberately bad representation is to put it inside an iteration and let the iteration repair it.

That works. What it costs is the interesting part.

What a hierarchical preconditioner costs, against how accurate it is, at κ = 2.09e+1The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 13, 8, 6, 3, 2, 2 steps against 41 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 0.5 where the preconditioner has rank 1. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 2.42 times the cheapest, having taken 11 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 41 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ21steps, no preconditioner41cheapest ε0.5its rank1tightest ⁄ cheapest2.4the count is what is printedand the cost is what is spent
Fig. 1 Two curves going opposite ways. The dashed one is the iteration count and the solid one is the work, against how accurate the preconditioner is.

The set-up

Conjugate gradients on the true matrix A, preconditioned by a solve with A_H. The iteration never sees the representation as the operator — every step forms one product with the real A — so the answer it converges to is the answer to the real problem, at whatever tolerance the iteration is asked for, and the compression cannot affect it.

What the compression affects is the rate. A_H is a nearby matrix, so M⁻¹A is near the identity, and how near is exactly the backward error the previous essay measures.

At 128 unknowns and κ = 20.9, with 41 steps unpreconditioned:

ε of the preconditioner rank steps multiplications
5·10⁻¹ 1 13 4.81·10⁵
10⁻¹ 2 8 5.23·10⁵
10⁻² 3 6 6.26·10⁵
10⁻⁴ 6 3 7.13·10⁵
10⁻⁶ 8 2 7.43·10⁵
10⁻⁸ 10 2 1.16·10⁶
none — 41 6.72·10⁵

The step column falls monotonically. The work column rises monotonically. And the best entry in the work column is the rank-one preconditioner, which is a representation nobody would accept as an answer to anything.

Why the two columns disagree

Each step of the preconditioned iteration costs one product with A, which is n² and fixed, plus one solve with A_H, which is not fixed. The solve’s cost grows with the ranks, and the ranks grow with the accuracy.

So the total is

steps × (n² + cost of one solve with A_H) ,

and the accuracy pushes the first factor down and the second up. There is no reason for the product to be monotone in either direction, and it is not.

At κ = 20.9 the first factor runs out of room almost immediately. Thirteen steps is already cheap; getting to two saves eleven products with A, and buying them costs a preconditioner of rank ten where rank one would have done, whose solve is several times more expensive. The trade stops paying after the first step of it.

The iteration count is the number a solver prints and the number every convergence plot in this collection is drawn against, and here it is the wrong number. That is not a criticism of iteration counts — they are the right currency when the preconditioner is fixed, which is the situation every other preconditioning essay on this site is in. It is a warning about what happens when the preconditioner acquires a knob.

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. 2 The currency this essay is arguing with, from the field that established it: iteration count as the measure of what a preconditioner is worth, which it is whenever the preconditioner does not have a parameter.

The optimum moves, and with exactly one thing

If the trade is between iteration count and cost per iteration, the balance ought to shift when the iterations get more expensive to avoid — which is to say when the problem gets harder.

The matrix has a shift in it that does nothing to the geometry. Every distance between every pair of clusters is unchanged, so every rank in the partition is unchanged, and the condition number moves over three decades:

shift κ steps, unpreconditioned cheapest ε its rank its cost
1/n 2.09·10¹ 41 5·10⁻¹ 1 4.81·10⁵
4/n 1.75·10² 72 5·10⁻¹ 1 5.92·10⁵
16/n 1.53·10³ 105 10⁻⁶ 6 4.27·10⁵
64/n 1.07·10⁴ 119 10⁻⁸ 5 3.45·10⁵

Six decades of tolerance between the easiest problem and the hardest. On the easy one the best preconditioner in the comparison is rank one; on the hard one it is the tightest accuracy on the sweep.

And the ranks the partition asks for at a fixed accuracy do not move across that table — because the shift is not a geometric quantity, and this field’s ranks are functions of geometry. So the answer is a statement about the conditioning of the problem and about nothing in the matrix’s structure at all.

The accuracy worth paying for in a preconditioner, against the condition number of the problemFour complete cost sweeps, each reduced to the ε that minimised the total work, against the κ of the problem it was run on. The shift that moves κ over 2.7 decades leaves every distance between every pair of clusters exactly where it was, so the ranks the partition asks for stay at 10, 9, 7, 5 throughout and nothing geometric is moving. The optimum walks 0.5 → 0.5 → 10⁻⁶ → 10⁻⁸ — six decades of accuracy, bought because the problem got harder and for no other reason. On the easiest problem the best hierarchical preconditioner in the comparison has rank one; on the hardest it is the tightest one on the sweep. How accurate an approximate inverse should be is a question with an answer, and the answer is not in the matrix's structure.10¹10²10³10⁴10⁵-9-7-5-3-11condition number of the problemlog₁₀ ε at the cheapest totalrank 1rank 1rank 6rank 5the knob answers to κκ = 210.5κ = 1.1·10⁴10⁻⁸rank at the easy end1rank at the hard end5ranks at 10⁻⁸, spread5the geometry did not moveand the right accuracy did
Fig. 3 Four complete sweeps reduced to four points: the accuracy that minimised the work, against the condition number of the problem it was run on.
What a hierarchical preconditioner costs, against how accurate it is, at κ = 1.07e+4The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 15, 15, 14, 5, 3, 2 steps against 119 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 10⁻⁸ where the preconditioner has rank 5. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 1.00 times the cheapest, having taken 0 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 119 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ1.1·10⁴steps, no preconditioner119cheapest ε10⁻⁸its rank5tightest ⁄ cheapest1the count is what is printedand the cost is what is spent
Fig. 4 The hardest of the four sweeps, where the cost curve has its minimum at the far end and the unpreconditioned line is far above everything.

What the rank-one preconditioner is

A representation at ε = ½ deserves a paragraph, because rank one sounds like a joke and the object is not one.

Every admissible block is replaced by its single dominant term: the outer product of the block’s leading left and right singular vectors, scaled by the leading singular value. Every dense block is kept exactly. So the preconditioner is the exact near field plus a rank-one summary of every far-field interaction, and it reproduces the matrix to 7.8 per cent.

Physically that is the most natural approximation there is. The interaction between two well-separated clusters, to leading order, is the total of one cluster acting on the total of the other — a monopole against a monopole — and rank one is exactly that. The second term is the dipole, the third the quadrupole, and the columns-per-decade measurement in this field’s second essay is the statement that each multipole order buys a couple of digits.

So the first row of the table is not a degenerate corner of the sweep. It is the classical approximation the whole subject started from, and on an easy problem it is the best preconditioner in the comparison. What the later rows buy is multipole orders nobody needed.

Two consequences worth carrying. The first is that a hierarchical preconditioner at low accuracy is cheap to build as well as to apply: rank one per block means one dominant singular pair per block, which a few steps of a power iteration find without a decomposition. The second is that the whole first half of this field’s storage story — the partition, the tree, the admissibility test — is doing its work at ε = ½ exactly as much as at 10⁻⁸. The structure is what is expensive to decide and it is decided once.

A 256 × 256 kernel matrix partitioned by the weak rule, with each compressed block's rank46 blocks: 16 kept dense and 30 stored as two thin factors, whose ranks run from 8 to 12. The rule decides one thing — whether a pair of index clusters may be compressed — and it decides it from four numbers about where those clusters sit, before any entry of the matrix is read. The weak rule compresses every off-diagonal block there is, including the two touching halves at the top level, whose rank is the largest number on the picture and is the one that climbs as the matrix grows. The whole thing stores 24,064 numbers against 65,536 entries, and reproduces the matrix to 1.08·10⁻⁹.121210109988889988881010998888998888rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe weak partitionblocks46kept dense16largest rank12numbers stored2.4·10⁴‖A − A_H‖ ⁄ ‖A‖1.1·10⁻⁹the picture is decidedbefore a number is read
Fig. 5 The structure that is decided once, from the partition field. At ε = ½ every light block in this picture carries the number one.
Error of the best rank-k approximation to a 10×10 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.12345678910⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 4.3·10⁻⁹worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁4.3·10⁻⁹worst Frobenius discrepancy4.3·10⁻⁹κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 6 And what one term of an expansion is worth, from the spectra field: the error of the best rank-one approximation is the second singular value, which for an admissible block is a few per cent.

Which is the same shape as four other knobs here

This collection has a field about what happens when a method acquires a parameter and it found that four of the site’s knobs are one knob: a regularisation parameter, a truncation rank, a stopping index and a filter width all control how much of the data to believe, and how much is worth believing belongs to the problem rather than to the method.

This is a fifth of that kind and it belongs in the same sentence with one word changed. How accurate an approximate inverse should be is a question whose answer is in the problem — specifically in κ — and not in the approximate inverse.

What is different, and what makes it worth a separate essay, is that the four knobs in that field all control an accuracy of the answer and this one does not. The answer’s accuracy is set by the iteration’s stopping tolerance and is the same at every row of both tables above. What the knob controls is only how much work the iteration does to get there, so the trade is entirely inside the cost and none of it is inside the result.

That is a rare and comfortable position. Every other parameter on this site has to be chosen against a quantity nobody can measure; this one has to be chosen against a quantity — total work — that the code can measure exactly, on the run it is doing, and adjust.

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. 7 The field that found four knobs to be one, from the essay that put them on one problem: what each of them controls is how much of the data to believe.

Against no preconditioner at all, which is the other comparison

Every row of both tables beats the unpreconditioned run, and by how much is worth reading because it is the part that makes the whole exercise worth doing.

At κ = 20.9 the unpreconditioned iteration takes 41 steps and 6.72·10⁵ multiplications; the best preconditioned run takes 13 and 4.81·10⁵. That is a saving of 28 per cent, which is real and unexciting. At κ = 1.07·10⁴ the unpreconditioned run takes 119 steps and 1.95·10⁶; the best preconditioned run takes 2 and 3.45·10⁵. That is a factor of 5.7.

So the preconditioner earns its keep in proportion to how hard the problem is, which is what a preconditioner is supposed to do and is worth confirming rather than assuming. It also means the two tables are answering two different questions and both are worth asking: is a hierarchical preconditioner worth having has the answer yes and increasingly so, and how accurate should it be has the answer that moves.

There is a third question the tables do not answer and it is the one the next essay takes. At every row above, the representation is being used as a preconditioner — applied inside an iteration that converges to the true answer. It could instead be used as a solver, applied once, returning an answer at the backward error the previous essay measured. On this matrix at 128 unknowns that costs 5.66·10⁵ multiplications, which is worse than four of the six rows in the first table, and the comparison flips completely by 512 unknowns.

The filter 8 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 6.9·10⁻¹⁴ and are drawn as one curve. It rises above one — 1.070 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 predicted6.9·10⁻¹⁴largest Tikhonov factor0.97two routes to the same curveand a curve that goes above one
Fig. 8 What a preconditioner does to the spectrum an iteration is working against, from the field that measures it.
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. 9 And why the counts here are so much better than the standard bound: a preconditioned spectrum in clusters is solved far faster than its condition number allows.

What a code should actually do

The measurement suggests an obvious policy and the obvious policy is nearly right.

Since the total is measurable and the curve has one minimum, a code can start with a cheap preconditioner, solve, and refine the representation only if the iteration count came out high. The refinement is cheap in a way worth noticing: a hierarchical representation compressed at 10⁻² can be re-compressed at 10⁻⁴ from the blocks it already has, if the original block decompositions were kept, so a code that is prepared to spend the memory can walk down the curve without recomputing.

Two things spoil this and both are worth stating.

The first is that the setup cost is not on the table. Every row above prices only the iterations; the representation had to be built, and at a tight tolerance building it costs more. This collection has an essay in the sequence field about exactly that failure — a rule that reads a free measurement and beats every constant, except that the free measurement cannot see what a rebuild cost, and at the ends of the range the setup cost is the whole decision. The same hole is here.

The second is that the right preconditioner accuracy depends on how many right-hand sides there are. One solve amortises the setup over one iteration count; a hundred solves with the same matrix amortise it over a hundred, and then a tighter representation is worth building. The sequence field is the field about that, and this essay’s table is the one-right-hand-side column of a larger one.

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 8 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 0 of the 6 ratios and loses at the rest, because it cannot see the setup cost at all.10¹10²10³0246810what 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 504setup worth 1006setup worth 2008how long to keep itis a question about what it cost to build
Fig. 10 The failure this policy shares, from the sequence field: a rule that reads what the iteration is doing cannot see what the setup cost, and at the ends of the range the setup cost is the decision.
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. 11 And the field where the amortisation happens: an expensive object built once and used many times, and the arithmetic that decides how many is many.

The half of the sweep that is monotone, and the half that is not

There is a detail in the four condition-number sweeps worth pulling out, because it is the thing that makes the policy in the last section workable at all.

The step column is monotone at every condition number. At κ = 20.9 it goes 13, 8, 6, 3, 2, 2; at κ = 1.07·10⁴ it goes 15, 15, 14, 5, 3, 2. It never rises. That is not obvious in advance — a preconditioner is not required to improve when it gets more accurate, and this collection has a whole essay about a stale preconditioner that improves the problem while getting worse at it — but it is what happens, and it means the code never has to worry about overshooting the accuracy in the direction of the iteration count.

The cost column is not monotone at any condition number and its minimum is interior at two of the four. That is what makes the search a search rather than a ratchet.

And there is a shape in the step columns worth naming. At the two harder problems the count is flat for the first two or three tolerances — 15, 15, 14 — and then falls off a cliff. A preconditioner too crude for a hard problem does nothing at all, several decades of refinement do nothing at all, and then one further decade does everything. A code that walked the curve from the cheap end and stopped when refining stopped helping would stop at the first row and never find the minimum.

So the honest policy is not refine while it helps. It is refine past the point where it stops helping, at least once, which is an uncomfortable rule and is what the measurement says.

What a hierarchical preconditioner costs, against how accurate it is, at κ = 1.75e+2The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 16, 14, 7, 4, 3, 2 steps against 72 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 0.5 where the preconditioner has rank 1. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 1.59 times the cheapest, having taken 14 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 72 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ175steps, no preconditioner72cheapest ε0.5its rank1tightest ⁄ cheapest1.6the count is what is printedand the cost is what is spent
Fig. 12 The second sweep, where the count is 16, 14, 7, 4, 3, 2 and the plateau at the cheap end has almost disappeared.
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. 13 The rate the count is against, from the iterative field, and the bound this collection has repeatedly found twelve times too pessimistic once the spectrum is in clusters.

The refusal

The claim under test is one that is stated in nearly every paper on preconditioning and is nearly always harmless: that a preconditioner which takes fewer iterations is a better preconditioner.

It is harmless when the preconditioners being compared cost the same to apply, which is the usual case. It is false as soon as one of them has an accuracy knob, and the refusal is fed the two ends of the easiest sweep: the tightest preconditioner takes 2 steps against the cheapest one’s 13, and costs 1.16 million multiplications against 0.48 million. The assertion that the tighter one costs less is handed those numbers and it fails.

The general form is one this collection keeps meeting. A currency is only a currency while the exchange rate is fixed, and an iteration count is a currency for work exactly as long as an iteration costs a constant. Two other essays here are about the same substitution going wrong — a flop count that stopped predicting the time once the memory hierarchy appeared, and a residual that stopped predicting the error once the condition number did.

What a hierarchical preconditioner costs, against how accurate it is, at κ = 1.53e+3The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 17, 17, 8, 4, 2, 2 steps against 105 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 10⁻⁶ where the preconditioner has rank 6. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 1.34 times the cheapest, having taken 0 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 105 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ1527steps, no preconditioner105cheapest ε10⁻⁶its rank6tightest ⁄ cheapest1.3the count is what is printedand the cost is what is spent
Fig. 14 The sweep where the two readings disagree least, in the middle of the range, where the cost curve is nearly flat across four decades of tolerance and almost any choice is defensible.
Words moved against the block size, n = 96, M = 100A U-shaped curve on logarithmic axes. Too small a block refactors the panel too often; too large a block does not fit in the fast memory and the tiled update evicts what it is about to read. The minimum is at b = 8, which is √M − 2. The three-line count says √(M/3) = 5, which has the right scaling and the wrong constant.110¹10²10⁵block size bwords movedthe count: √(M/3) = 5measured best: b = 8words movedat the best block1.2·10⁵at b = 15.5·10⁵at b = 243.2·10⁵derived from M with no measurement, and scannedthe two agree
Fig. 15 The same substitution failing in the cost field, from the essay that measured it: a count of operations that stopped being a prediction of the time.
Multiplications in a hierarchical solve and in a dense factorisation, and where they crossBoth counted rather than estimated: the recursion carries a counter and reports what it actually did. The dense factorisation is n³/3 and rises at exactly three per doubling. The hierarchical solve rises at 2.13, 1.93, 1.74 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 64 and n = 128: below it the format is the more expensive way to solve the system, at 1.48 times the dense count, and at n = 512 it is 6.2 times cheaper. Every point returns an answer at a backward error of about 1.4·10⁻¹⁰, so the comparison is between two ways of getting the same thing.567891010⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 641.5ratio at n = 5120.16exponent, first doubling2.1exponent, last doubling1.7backward error1.4·10⁻¹⁰cheaper is a sizenot a property
Fig. 16 And the direct alternative to all of this, from the next essay: solving with the representation outright, which on this matrix is cheaper than every preconditioned run in the first table.

The sweep at its other settings

What a hierarchical preconditioner costs, against how accurate it is, at κ = 1.75e+2The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 16, 14, 7, 4, 3, 2 steps against 72 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 0.5 where the preconditioner has rank 1. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 1.59 times the cheapest, having taken 14 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 72 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ175steps, no preconditioner72cheapest ε0.5its rank1tightest ⁄ cheapest1.6the count is what is printedand the cost is what is spent
Fig. 17 The second of the four, where the plateau at the cheap end has almost disappeared and the cost curve is nearly flat.
What a hierarchical preconditioner costs, against how accurate it is, at κ = 1.07e+4The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 15, 15, 14, 5, 3, 2 steps against 119 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 10⁻⁸ where the preconditioner has rank 5. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 1.00 times the cheapest, having taken 0 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 119 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ1.1·10⁴steps, no preconditioner119cheapest ε10⁻⁸its rank5tightest ⁄ cheapest1the count is what is printedand the cost is what is spent
Fig. 18 And the hardest, where the minimum has walked to the far end of the sweep and the unpreconditioned run is a factor of five above everything.
The accuracy worth paying for in a preconditioner, against the condition number of the problemFour complete cost sweeps, each reduced to the ε that minimised the total work, against the κ of the problem it was run on. The shift that moves κ over 2.6 decades leaves every distance between every pair of clusters exactly where it was, so the ranks the partition asks for stay at 9, 8, 6, 4 throughout and nothing geometric is moving. The optimum walks 0.5 → 0.5 → 0.5 → 10⁻⁴ — six decades of accuracy, bought because the problem got harder and for no other reason. On the easiest problem the best hierarchical preconditioner in the comparison has rank one; on the hardest it is the tightest one on the sweep. How accurate an approximate inverse should be is a question with an answer, and the answer is not in the matrix's structure.110¹10²10³10⁴10⁵-9-7-5-3-11condition number of the problemlog₁₀ ε at the cheapest totalrank 1rank 1rank 1rank 2the knob answers to κκ = 170.5κ = 634010⁻⁴rank at the easy end1rank at the hard end2ranks at 10⁻⁸, spread5the geometry did not moveand the right accuracy did
Fig. 19 The staircase at half the size, where the whole cost surface is lower and its minimum sits in the same place.
The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.983 and sit a constant 7.4× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.98representation at 10⁻⁸1.4·10⁻⁹backward error there1.4·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.3·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 20 What the preconditioner’s accuracy knob actually sets, from the error field, and the reason the iteration converges to the true answer whatever it is set to.
Numbers stored per unknown, against the size of the matrix, at 0.01The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 45, 57, 67, 78 and 28, 34, 42, 50, adding about 11 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 15 and 10 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.33, 1.25, 1.21, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5124·10⁴weak, n = 5122.6·10⁴dense, n = 5122.6·10⁵per doubling11‖A − A_H‖ ⁄ ‖A‖6.8·10⁻⁴the dense line doublesand the other two add a constant
Fig. 21 What a rank-one preconditioner costs to store, from the partition field: the dense fringe and almost nothing else.
Multiplications in a hierarchical solve and in a dense factorisation, and where they crossBoth counted rather than estimated: the recursion carries a counter and reports what it actually did. The dense factorisation is n³/3 and rises at exactly three per doubling. The hierarchical solve rises at 1.82, 1.80, 1.62 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 64 and n = 64: below it the format is the more expensive way to solve the system, at 0.28 times the dense count, and at n = 512 it is 47.4 times cheaper. Every point returns an answer at a backward error of about 2.6·10⁻⁴, so the comparison is between two ways of getting the same thing.567891010³10⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 640.28ratio at n = 5120.021exponent, first doubling1.8exponent, last doubling1.6backward error2.6·10⁻⁴cheaper is a sizenot a property
Fig. 22 And the direct alternative at the same accuracy, which at these sizes is the cheaper of the two.
A 256 × 256 kernel matrix partitioned by the strong rule, with each compressed block's rank112 blocks: 46 kept dense and 66 stored as two thin factors, whose ranks run from 4 to 5. The rule decides one thing — whether a pair of index clusters may be compressed — and it decides it from four numbers about where those clusters sit, before any entry of the matrix is read. The strong rule refuses any pair whose clusters touch, and subdivides instead, so the diagonal is fringed with small dense blocks and no rank on the picture exceeds 5. The whole thing stores 27,008 numbers against 65,536 entries, and reproduces the matrix to 3.41·10⁻¹⁰.545545545455554545545545545454555555454545545545545455554545545545rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe strong partitionblocks112kept dense46largest rank5numbers stored2.7·10⁴‖A − A_H‖ ⁄ ‖A‖3.4·10⁻¹⁰the picture is decidedbefore a number is read
Fig. 23 The structure that is decided once and used at every accuracy on the sweep.
Incomplete Cholesky on the 10×10 grid: κ 48.4 → 5.12A semi-logarithmic plot of relative residual against iteration for conjugate gradients with and without an incomplete Cholesky preconditioner, the preconditioned curve falling faster.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×
Fig. 24 What a preconditioner does to a spectrum, from the field that established the currency this essay is arguing with.
Iterative refinement on a factorisation of a matrix that is not the one being solvedThe last member of the sequence, solved with the first member's pivot order. Every other repair in this collection works because the outer loop recomputes the residual from the matrix, and this one does that too — and it does not work. The perturbed pivots are not a rounding, they are a change to the matrix, and the residual that would reveal it is computed with the same badly scaled rows that caused it. Six steps move the backward error from 4.82·10⁻⁹ to 5.46·10⁻¹⁰, a factor of 8.8, and stop. The equilibrated run has nothing to repair: it starts at 6.55·10⁻¹⁷.0123456710⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷refinement stepsbackward error of the solvethe kept order, in the units as giventhe same order, rows equilibratedthe repair that does not repairbefore refinement4.8·10⁻⁹after six steps5.5·10⁻¹⁰the factor gained8.8equilibrated, before6.5·10⁻¹⁷pivots replaced2the residual is recomputedand it cannot see what was done to the matrix
Fig. 25 And the one case on this site where an outer loop does not repair a deliberately wrong inner object, which is the boundary of the reading this essay depends on.

Where this goes next

10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.99representation at 10⁻⁸1.2·10⁻⁹backward error there1.3·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.5·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem An accuracy that is a backward error Every backward error on this site is something an algorithm produced and somebody then measured. This one is a line in the program. Solving with a compressed matrix gives a residual that is the compression's own error, at a slope of 1.000 over ten decades, so the knob that sets the storage sets the backward error directly. 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. 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 Four knobs and one floor A truncation, a Tikhonov parameter, a step count and a randomised rank, on one problem with an answer that is known. Their best errors are 0.1445, 0.1406, 0.1426 and 0.1449 — a spread of 3% across four methods that share no arithmetic.

Reads more easily once this is understood

Essays that name this one as worth reading first.

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.

Backward errorCondition numberConjugate gradientsFlop countHierarchical matrixParameter choicePreconditioningStopping criterionToleranceWoodbury identity