Numerical Linear Algebra
where the algebra stops being the arithmetic
Where the flop count stopped predicting the time

Where the format starts paying

A hierarchical solve costs 1.48 times a dense factorisation at 64 unknowns and 0.16 times it at 512. The crossover is between 64 and 128, it walks right when the accuracy is tightened, and the exponent between consecutive sizes is 2.13, 1.93, 1.74 — falling towards one and never arriving.
33 min read 25 figures Measured, not assertedWhere the algebra stops being the arithmeticAssertions that reject

Worth reading first: The accuracy worth paying for · The same arithmetic at a different price.

Every asymptotic claim in numerical linear algebra is a claim about sizes nobody has. This essay is about the size somebody does have, which is a different question with a different answer.

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. 1 Multiplications in a hierarchical solve and in a dense factorisation, both counted rather than estimated, and the size at which they cross.

What the recursion is

At the top level, a hierarchical matrix in the weakly admissible partition is

A = [ A₁₁ U₁V₂ᵀ ] = D + Ũ Ṽᵀ , [ U₂V₁ᵀ A₂₂ ]

with D block diagonal and the off-diagonal part of rank k₁ + k₂. The Woodbury identity says

A⁻¹ = D⁻¹ − D⁻¹Ũ (I + ṼᵀD⁻¹Ũ)⁻¹ ṼᵀD⁻¹ ,

so a solve with A is a solve with D against k₁ + k₂ + 1 right-hand sides, plus one dense solve of size k₁ + k₂ — and a solve with D is two solves one level down. The recursion bottoms out on dense leaf blocks. Nothing anywhere in it forms an n × n matrix, and nothing in it is approximate: the answer is the exact solution of A_H x = b, and how far that is from the answer to Ax = b is the previous field’s measurement.

Counting what it costs is not a modelling exercise. The recursion carries a counter, every multiplication it performs increments it, and the numbers below are what the code did.

The measurement

n hierarchical dense, n³⁄3 ratio backward error
64 1.29·10⁵ 8.74·10⁴ 1.48 1.2·10⁻¹⁰
128 5.66·10⁵ 6.99·10⁵ 0.81 1.4·10⁻¹⁰
256 2.16·10⁶ 5.59·10⁶ 0.39 1.3·10⁻¹⁰
512 7.23·10⁶ 4.47·10⁷ 0.16 1.4·10⁻¹⁰

The last column is what makes the comparison legitimate. Every one of those runs returns an answer at the same backward error, so the two methods are being asked for the same thing.

The ratio column has a crossover in it, between 64 and 128 unknowns. Below it the format is the more expensive way to solve the system. Above it, and increasingly, it is not.

The exponent, which is the interesting number

A cost of n log²n has no exponent. What it has is an exponent that falls towards one and never arrives, and the way to see that is to take it between consecutive sizes rather than fitting it over a sweep:

doubling exponent
64 → 128 2.13
128 → 256 1.93
256 → 512 1.74

Against the dense factorisation’s 3.00 at every step, exactly, because it is a cube.

Three things are worth reading off that column. The first is that it is well under three throughout, which is why the ratio improves. The second is that it is falling, which is what a logarithm looks like from inside — at n = 512 there are five levels, so log₂ n contributes about a fifth of the answer, and it will contribute a tenth at n = 2¹⁰ and a twentieth much later. The third is that at n = 64 it is above two: at the small end this is not merely slower than the dense factorisation, it is growing faster than a quadratic, because the recursion’s overhead is concentrated in the small dense solves at the leaves and there are more of them per unknown when the tree is shallow.

Reporting a fitted single exponent over the whole sweep would give 1.43 and would say none of that. It would also be a number with no meaning, since the quantity it purports to describe has no exponent.

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.25, 1.95, 1.78 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 128 and n = 256: below it the format is the more expensive way to solve the system, at 2.36 times the dense count, and at n = 512 it is 3.4 times cheaper. Every point returns an answer at a backward error of about 2.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 = 642.4ratio at n = 5120.29exponent, first doubling2.2exponent, last doubling1.8backward error2.4·10⁻¹⁴cheaper is a sizenot a property
Fig. 2 At twelve digits, where the ranks are larger and the crossover has walked right — the dense line cannot move, because a factorisation has no accuracy parameter.
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. 3 And at two, where the ranks are tiny and the format is winning by n = 128.

What the accuracy does to the crossover

The dense line is fixed. A factorisation has no tolerance in it, so n³/3 is n³/3 at every accuracy, and the whole trade lives in the other curve.

Tightening ε raises the ranks, raises the cost of every level of the recursion, and walks the crossover right. That is the economics of the format in one sentence: what the accuracy knob buys is a size at which the format starts paying, and asking for more digits than the problem needs pushes that size out of reach.

The previous field’s essay is what closes this. The accuracy is a backward error and the digits it should be set to are the digits the answer needs divided by the condition number — so a well-conditioned problem needs a loose tolerance, gets small ranks, and crosses over early. An ill-conditioned one needs a tight tolerance, gets larger ranks, and crosses over late, which is a second and quite separate way in which conditioning costs money.

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.994 and sit a constant 6.7× 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.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
Fig. 4 The knob this essay’s curve is a function of, from the error field: what the compression tolerance actually sets.

Where the multiplications go

The counter reports one number and the number is a sum of four quite different things, which is worth separating because three of them behave well and one does not.

The leaves. At the bottom of the recursion sit n/leaf dense solves of size leaf, each performed against every right-hand side that reaches it. That is where the small-n overhead lives: at n = 64 with a leaf of 16 there are four of them and the tree is two levels deep, so the fixed cost of the recursion is spread over very little work.

The block products. Each level applies the thin factors Ũ and Ṽᵀ, which costs O(k·n) per right-hand side per level. Linear in the rank, linear in the size, and there are log n levels — this is the term the whole format is built around.

The small dense solves. At each level the Woodbury correction requires factorising (I + ṼᵀD⁻¹Ũ), which is (k₁ + k₂) square. That is a cube in the rank and a constant in n, which is why the rank matters more here than the storage tables suggest.

And the right-hand sides multiply. This is the term that surprises. Solving with A against one vector requires solving with D against k₁ + k₂ + 1 vectors, and solving with those requires solving one level down against (k₁ + k₂ + 1) times the next level’s rank plus one. The number of right-hand sides at depth ℓ grows like the product of the ranks above it. For small ranks that is fine; for the ranks a touching block acquires at 512 unknowns it is the reason the exponent is 1.74 rather than closer to one.

That last term is the strongest argument in this whole phase for the admissibility test the hierarchy field measures as more expensive. A rank that climbs from 9 to 13 costs 44 per cent in storage and considerably more than that here, because it appears in a product along the depth of the tree.

The largest rank in each partition, and what refusing to compress a touching pair costsThe strong rule's worst block is rank 5 at every size — one number across a factor of eight — because it never compresses a pair of clusters that touch. The weak rule compresses them and its worst rank climbs 9, 10, 12, 13, at about one per doubling, which is the touching block's logarithm arriving inside a whole partition. That is what the test buys. What it costs is on the badge: at every size measured, the finer partition stores MORE — 67,968 numbers against 61,440 at n = 512 — because it pays in blocks what it saves in rank, and the blocks near the diagonal are dense. The bounded rank is an asymptotic argument and the sizes here are not asymptotic.567891003691215log₂ nlargest rank in the partitionweak: touching pairs compressedstrong: touching pairs refusedthe test bounds a rank and costs storagestrong, blocks250weak, blocks94strong, numbers6.8·10⁴weak, numbers6.1·10⁴strong ⁄ weak1.1the better partitionis the more expensive one
Fig. 5 The rank that appears in that product, from the hierarchy field: one partition holds it at five and the other lets it climb.
Two-level convergence factor for three interpolations, at ε = 0.001Three bars per operator: the classical formula, the unconstrained energy minimiser, and the minimiser constrained to reproduce a constant. On the isotropic and aligned operators the constraint costs a factor of 4.8; on the rotated one it buys a factor of 1.33.lower is better; each bar is the factor the error falls by per two-level cycleisotropic, classical0.0609isotropic, minimiser0.0609isotropic, constrained0.0848aligned, classical0.0612aligned, minimiser0.0615aligned, constrained0.2933rotated 45°, classical0.3464rotated 45°, minimiser0.2843rotated 45°, constrained0.2131what the constraint is worthconstraint at isotropic1.4constraint at aligned4.8constraint at rotated 45°0.75the constraint costs where the method worksand buys where it does not
Fig. 6 A cost that is cubic in a small dimension and constant in a large one, from the least squares field, which is the same shape as the Woodbury correction here.

What a multiplication count is not

This is the cost field, and the cost field’s standing result is that an operation count stopped predicting the time some decades ago. It is worth applying that to the figure above rather than letting the figure stand as a runtime claim.

The dense factorisation is the best-behaved kernel in scientific computing. It is three nested loops over contiguous memory, it blocks perfectly, and a tuned implementation reaches a large fraction of a machine’s peak rate. Its 4.47·10⁷ multiplications at n = 512 are cheap multiplications.

The hierarchical solve is not that. It is a recursion over a tree, its leaves are small dense solves that do not block, its intermediate objects are thin factors with poor locality, and it spends a significant fraction of its time on bookkeeping that this counter does not count at all. Its 7.23·10⁶ multiplications are expensive multiplications.

So the ratio of 0.16 at n = 512 is a ratio of two counts and not a ratio of two runtimes, and the runtime ratio is certainly worse. How much worse is a property of the machine, which is this field’s oldest observation: the-same-arithmetic-at-a-different-price measures two orderings of one factorisation with identical flop counts and a memory traffic differing by a factor that decides which is usable.

The honest form of the claim is therefore about the exponent rather than the constant. A method whose exponent is falling towards one beats a method whose exponent is three, at some size, whatever the constants are — and where that size is, in seconds rather than in multiplications, is a question no count answers.

A 6.25-bit block format against 8-bit E4M3, read two waysFour curves of relative error against the spread of the data. Two, for the block format, rise steeply when read as a median and gently when read as a norm; the other two are nearly flat.048121610⁻²10⁻¹1octaves of spread within a blockrelative errorblock, medianblock, 2-normE4M3, medianE4M3, 2-normthe same data, two readingsblock format, bits a value6.3norm across 16 octaves1.8median across 16 octaves62entries deleted at the wide end383640 values, blocks of 32a norm is dominated by what a block format keeps
Fig. 7 The field’s oldest result, from the essay that established it: two computations with identical operation counts and different costs.
Communication for a 512×12 factorisation on 16 processorsTwo counts for the same factorisation. Rounds on the critical path: 48 for the column sweep, 4 for the reduction tree, 4 for Cholesky QR — a factor of n between the first and the other two. Words sent: 1170, 1170 and 2160 — the sweep and the tree send the same number, and the method with the fewest rounds sends the most.rounds on the critical pathHouseholder sweep48reduction tree4Cholesky QR4words sentHouseholder sweep1170reduction tree1170Cholesky QR2160two counts, two rankingsrounds, sweep ÷ tree12words, Cholesky ÷ tree1.8arithmetic, tree ÷ sweep1.5the rounds separate the threeand the words do not
Fig. 8 And the currency that replaced the operation count, from the essay that counted it: words moved rather than arithmetic performed.

Why the exponent is not one

The literature’s headline for this family of methods is linear complexity, and the exponent column above says 1.74 at the largest size measured. Both are correct and the gap between them is worth closing.

The n log n storage argument in this field’s fifth essay needs the rank of every block to be a constant. The solve needs more than that: it needs the number of right-hand sides propagating down the recursion to stay bounded, and that number is a product of ranks along a path from the root to a leaf. With a bounded rank k and L levels it is O(k^L) in the worst case and O(kL) in practice, because the correction at each level acts on a subspace the level below has already reduced.

What the measurements here have is a weakly admissible partition whose ranks are not constant — 9, 10, 12, 13 across the sweep, as this phase’s third essay establishes — so the product along the path is growing at both ends. The exponent is therefore an honest 1.74 for this partition at these sizes, and it would be lower for the strongly admissible one whose ranks are five throughout, if that partition had a comparably simple solve, which it does not.

So linear is a statement about a limit, about a constant rank, and about a formulation with recompression at every level of the recursion rather than the plain Woodbury one used here. Each of those is a real thing and none of them is what the figure measures. What the figure measures is the exponent of a specific, simple, exactly-arithmetic recursion, and 1.74 falling is what that is.

The rank of an admissible block and of a touching one, against how finely they are sampledBoth blocks are of the kernel 1/r; both are 32, 64, 128, 256 points a side; both are truncated at 10⁻⁸. The admissible pair — [0, 1] against [2, 3] — needs 5, 5, 5, 5 columns, which is one number. The touching pair — [0, 1] against [1, 2] — needs 9, 11, 12, 13, climbing by about one per doubling, which is a logarithm. Neither of them grows like the block, and only one of them stops. That difference is what the admissibility test in a partition is buying, and it is why the touching pair is kept dense rather than compressed at all.45678903691215log₂ of the points a sidecolumns above 10⁻⁸two intervals that touchtwo that do notthe measurement that is a null resultadmissible, n = 325admissible, n = 2565touching, n = 329touching, n = 25613stored ⁄ dense at largest0.039the rank belongs to the geometryand not to the sampling
Fig. 9 The rank that is not constant, from the hierarchy field’s third essay: the touching block’s, which is the one the weakly admissible partition compresses.

Against the other things available

The comparison on the figure is against a dense factorisation, which is the fair one — both are direct methods returning an answer at a stated backward error. Two other options are worth pricing.

An iterative solve on the dense matrix. Conjugate gradients at n = 512 on this matrix takes about 50 steps at n² multiplications each, which is 1.3·10⁷ — twice the hierarchical solve and three times cheaper than the dense factorisation. That is a reminder that a well-conditioned matrix rarely needs a factorisation at all, and the whole comparison above is between two methods that earn their keep when there are many right-hand sides or when the conditioning is bad.

A hierarchical solve used as a preconditioner. The previous essay measures that: at 128 unknowns the best preconditioned run costs 4.81·10⁵ multiplications against the direct hierarchical solve’s 5.66·10⁵, so the preconditioner wins slightly at that size. At 512 it does not, because the direct solve’s exponent is 1.74 and the preconditioned iteration’s is 2 — each of its steps is a dense product with A.

Which is a pleasing symmetry to end on. The preconditioned iteration is cheaper on small problems and the direct hierarchical solve is cheaper on large ones, for the same reason that a dense factorisation is cheaper than either on very small ones: at every size the winner is the method whose overhead matches the problem, and the overheads sort in the opposite order from the exponents.

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. 10 The third option priced, from the previous essay: the same representation used inside an iteration rather than as the solve.
Conjugate gradients at κ = 102, 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κ bound58 steps40×40, spectrum spread evenly in logbound permits 142
Fig. 11 And the second, from the iterative field: what a well-conditioned matrix costs an iteration that never factorises anything.

The comparison the figure cannot make

There is a size argument missing from all of this and naming it is more useful than pretending otherwise.

The dense factorisation at n = 512 needs 512² = 262,144 numbers of memory and it needs them all at once. The hierarchical representation at the same size needs 67,968. On this machine both are trivial; on the problems this format exists for, neither is, and the memory is what decides.

That is why the storage figure in the hierarchy field crosses over at a much smaller size than this one. Storage is already in the format’s favour at n = 64 — 3,456 numbers against 4,096 — while the arithmetic is still 48 per cent against it. So there is a band of sizes, roughly 64 to 128 here, where the format costs more arithmetic and less memory, and which of those a code should prefer is not a question about the matrix at all.

A second thing the figure cannot show is what happens with many right-hand sides. Both methods amortise a setup: the dense one factorises once and back-substitutes per vector at n² each; the hierarchical one builds a representation once and applies the recursion per vector at roughly k·n·log n each. So for many vectors the comparison is between n² and k n log n per solve, which crosses over far earlier than the figure’s crossover and in the format’s favour. Everything on the figure is the one-right-hand-side column.

And the third is the setup itself. The hierarchical representation on the figure was built by decomposing every off-diagonal block, which costs more than the dense factorisation it is being compared with — a fact this phase’s essay on building it from products alone exists to fix, and which the figure quietly assumes away by pricing only the solve.

Applications of an operator that is never assembled, against the entries of the matrix it stands forThe blocks at one level of the tree have disjoint column supports, so one batch of random vectors samples all of them at once: two batches a level for the ranges, two more for the projections, and one batch of leaf products for every diagonal block in the matrix together. That is leaf + 2(2k + p)·levels products, and the counter inside the operator says 112, 160, 208, 256 at n = 64, 128, 256, 512 — 48 more per doubling, which is a logarithm. The other line is n², the entries the compression route reads. At n = 512 that is 256 products against 262,144 entries, and the representation it produces is within 7.3× of the one that read them all.567891010¹10²10³10⁴10⁵10⁶log₂ nnumbers the construction touchedentries, n²products with the operatoran operator, applied a few hundred timesproducts at n = 512256entries at n = 5122.6·10⁵per doubling48‖A − A_H‖ ⁄ ‖A‖4·10⁻⁷excess over the compression7.3no entry of the matrixwas ever read
Fig. 12 The setup this figure does not price, and the route that makes it affordable: a few hundred applications of an operator that is never assembled.
Numbers stored per unknown, against the size of the matrix, at 10⁻⁸The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 54, 79, 106, 133 and 50, 70, 94, 120, adding about 26 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 26 and 23 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.55, 1.42, 1.33, 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 = 5126.8·10⁴weak, n = 5126.1·10⁴dense, n = 5122.6·10⁵per doubling26‖A − A_H‖ ⁄ ‖A‖3.8·10⁻¹⁰the dense line doublesand the other two add a constant
Fig. 13 And the resource that crosses over first, from the hierarchy field.

The refusal

The claim under test is the one every abstract about this subject makes: that a hierarchical solver is faster than a dense factorisation.

The verdict is asymptotic rather than false, because the claim is true and the word missing from it is eventually. The refusal is fed the ratio column: the assertion that the hierarchical solve costs less at every size on the sweep is handed 1.48, 0.81, 0.39, 0.16, and it fails on the first entry.

At 64 unknowns the format performs half again as many multiplications as the thing it replaces, for an answer that is no better, using a representation that took a decomposition of every off-diagonal block to build. It is comprehensively the wrong choice, and it is the wrong choice at a size that is perfectly ordinary — a 64 × 64 dense system is the sort of thing that appears inside larger algorithms by the thousand.

That is why the crossover is worth measuring rather than asserting, and why a figure of this kind belongs in the cost field rather than in the field that owns the method. A method’s own field is the last place a comparison it loses is likely to be drawn.

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.88, 1.73 — 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.92 times the dense count, and at n = 512 it is 10.3 times cheaper. Every point returns an answer at a backward error of about 2.3·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 = 640.92ratio at n = 5120.097exponent, first doubling2.1exponent, last doubling1.7backward error2.3·10⁻⁸cheaper is a sizenot a property
Fig. 14 The sweep at six digits, where the crossover sits between the first two points and the largest size is five times cheaper than the dense factorisation.
Numbers stored per unknown, against the size of the matrix, at 10⁻⁸The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 54, 79, 106, 133 and 50, 70, 94, 120, adding about 26 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 26 and 23 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.55, 1.42, 1.33, 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 = 5126.8·10⁴weak, n = 5126.1·10⁴dense, n = 5122.6·10⁵per doubling26‖A − A_H‖ ⁄ ‖A‖3.8·10⁻¹⁰the dense line doublesand the other two add a constant
Fig. 15 The other resource, from the hierarchy field, where the crossover happens earlier: storage is already in the format’s favour at n = 64.
‖AX̂ − I‖ and ‖X̂A − I‖ for a computed Hilbert inverse, to n = 12Two residuals of the same computed inverse, each divided by ‖A‖‖X̂‖ so both are dimensionless. The exact inverse satisfies both to zero. The computed one satisfies ‖AX̂ − I‖ at 3.5·10⁻¹⁸ — the side its columns were solved along — and ‖X̂A − I‖ at 6.8·10⁻¹⁶, 195 times larger, at n = 12. Each column of X̂ is the exact solution of a slightly perturbed system, but a different perturbation for each column, so there is no single nearby matrix whose inverse X̂ is.468101210⁻²⁰10⁻¹⁹10⁻¹⁸10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³nresidual / (‖A‖ ‖X̂‖)‖X̂A − I‖‖AX̂ − I‖the two sides, comparedn = 6, ratio6.4n = 8, ratio8.3n = 10, ratio82n = 12, ratio195the exact inverse satisfies bothand the computed one satisfies the side it was computed along
Fig. 16 And the collection’s standing warning about solving through an object rather than with one, from the field that measured what forming an inverse costs.

The crossover under the knob

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.18, 1.93, 1.78 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 128 and n = 256: below it the format is the more expensive way to solve the system, at 1.97 times the dense count, and at n = 512 it is 4.4 times cheaper. Every point returns an answer at a backward error of about 1.1·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 = 642ratio at n = 5120.23exponent, first doubling2.2exponent, last doubling1.8backward error1.1·10⁻¹²cheaper is a sizenot a property
Fig. 17 At ten digits, where the ranks are larger and the crossover has walked between the second and third points.
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.06, 1.87, 1.66 — 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.61 times the dense count, and at n = 512 it is 17.5 times cheaper. Every point returns an answer at a backward error of about 1.5·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 = 640.61ratio at n = 5120.057exponent, first doubling2.1exponent, last doubling1.7backward error1.5·10⁻⁶cheaper is a sizenot a property
Fig. 18 And at four, where it has walked back to the left of the figure entirely.
Numbers stored per unknown, against the size of the matrix, at 10⁻⁸The dense matrix stores n numbers per unknown, which is the straight line through the origin and doubles whenever n does. The hierarchical representations store 54, 79, 106, 133 and 50, 70, 94, 120, adding about 26 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 26 and 23 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.55, 1.42, 1.33, 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 = 5126.8·10⁴weak, n = 5126.1·10⁴dense, n = 5122.6·10⁵per doubling26‖A − A_H‖ ⁄ ‖A‖3.8·10⁻¹⁰the dense line doublesand the other two add a constant
Fig. 19 The resource that crosses over first, and by a wide margin: storage is in the format’s favour at every size measured.
The largest rank in each partition, and what refusing to compress a touching pair costsThe strong rule's worst block is rank 5 at every size — one number across a factor of eight — because it never compresses a pair of clusters that touch. The weak rule compresses them and its worst rank climbs 9, 10, 12, 13, at about one per doubling, which is the touching block's logarithm arriving inside a whole partition. That is what the test buys. What it costs is on the badge: at every size measured, the finer partition stores MORE — 67,968 numbers against 61,440 at n = 512 — because it pays in blocks what it saves in rank, and the blocks near the diagonal are dense. The bounded rank is an asymptotic argument and the sizes here are not asymptotic.567891003691215log₂ nlargest rank in the partitionweak: touching pairs compressedstrong: touching pairs refusedthe test bounds a rank and costs storagestrong, blocks250weak, blocks94strong, numbers6.8·10⁴weak, numbers6.1·10⁴strong ⁄ weak1.1the better partitionis the more expensive one
Fig. 20 The rank that appears in a product along the depth of the recursion, which is why the exponent here is 1.74 rather than closer to one.
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.994 and sit a constant 6.7× 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.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
Fig. 21 The knob this whole curve is a function of, from the error field.
Applications of an operator that is never assembled, against the entries of the matrix it stands forThe blocks at one level of the tree have disjoint column supports, so one batch of random vectors samples all of them at once: two batches a level for the ranges, two more for the projections, and one batch of leaf products for every diagonal block in the matrix together. That is leaf + 2(2k + p)·levels products, and the counter inside the operator says 112, 160, 208, 256 at n = 64, 128, 256, 512 — 48 more per doubling, which is a logarithm. The other line is n², the entries the compression route reads. At n = 512 that is 256 products against 262,144 entries, and the representation it produces is within 7.3× of the one that read them all.567891010¹10²10³10⁴10⁵10⁶log₂ nnumbers the construction touchedentries, n²products with the operatoran operator, applied a few hundred timesproducts at n = 512256entries at n = 5122.6·10⁵per doubling48‖A − A_H‖ ⁄ ‖A‖4·10⁻⁷excess over the compression7.3no entry of the matrixwas ever read
Fig. 22 And the setup this figure does not price, from the randomised field, where it costs a few hundred products rather than a decomposition of every block.
‖A − LLᵀ‖ ⁄ ‖A‖ for a Cholesky computed in the format, against how many truncations it performedThe recursion is the textbook one and every step stays in the format: the off-diagonal factor is exact, because a triangular solve against one factor of a rank-k block leaves a rank-k block, and the Schur complement is a rank-k addition to a hierarchical matrix, which is the only approximate step there is. So the number of truncations is a count — one per admissible block in the subtree being updated, at every level — and it runs from 0 at a leaf of 128 to 98 at a leaf of 8. The residual goes from 2.119·10⁻¹⁰ to 1.138·10⁻⁹, and all of that movement is the representation, whose own error is the upper curve and moves by the same factor: the ratio between them is 1.00, 0.93, 0.92, 0.91, 0.81, never above one. A hundred approximate operations contributed nothing measurable to the answer.02040608010010⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸truncations performed by the factorisationrelative errorthe representation's own error‖A − LLᵀ‖ ⁄ ‖A‖no decomposition without its residualresidual, 0 truncations2.1·10⁻¹⁰residual, 98 truncations1.1·10⁻⁹representation, deepest1.4·10⁻⁹residual ⁄ representation0.81levels, deepest5a hundred approximate stepsand an exact-looking factorisation
Fig. 23 What a formatted factorisation of the same object costs in accuracy, which is nothing measurable.
Four factorisations of one 512×12 matrix, κ = 7151A ranked bar chart of the departure from orthogonality, on a logarithmic scale. Classical Gram–Schmidt loses 4.6·10⁻¹⁰, modified 1.5·10⁻¹², Householder 2·10⁻¹⁴. The reduction tree, which touches the matrix once instead of 12 times, is of the same size as Householder's — and its error does not grow with the depth of the tree.classical Gram–Schmidt4.62·10⁻¹⁰modified Gram–Schmidt1.49·10⁻¹²Householder, one sweep2.03·10⁻¹⁴reduction tree, 16 leaves1.48·10⁻¹⁵departure from orthogonality, logarithmicthe tree, at four depths‖AᵀA − RᵀR‖/‖AᵀA‖, depth 13.4·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 24.3·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 31.7·10⁻¹⁵‖AᵀA − RᵀR‖/‖AᵀA‖, depth 41.5·10⁻¹⁵the same algebra, four timestwo of them are products of reflections
Fig. 24 A rearrangement that changes the currency rather than the arithmetic, from this field’s own essays, which is the discipline this comparison is borrowing.
Stored entries per row at each level, 400 unknowns and 1192 edgesA rising curve of stored entries per row against level in the hierarchy, beside a flat one for a geometric hierarchy on a grid.01234567110¹10²levelstored entries per rowalgebraicgeometrictwo complexities, one hierarchygrid complexity3.1operator complexity18geometric, for comparison1.5a ring plus random chords — nothing is anywherelevel 6 is 100% dense
Fig. 25 And the multigrid field’s arithmetic for the same shape of claim: a hierarchy whose levels each cost a fraction of the one below.

Where this goes next

10²10³10⁴10⁵matrix size ncountoperations, bothwords, unblockedwords, blocked (b = 6)the answer does not move‖PA − LU‖/‖A‖, unblocked2.8·10⁻¹⁶‖PA − LU‖/‖A‖, blocked2.8·10⁻¹⁶difference between them0the dashed curve is both orderings' operation countthe solid pair is what they cost The same arithmetic at a different price A blocked and an unblocked elimination perform 72,568 operations each — the same operations, associated differently — choose the same pivots, and return a factorisation identical to the last bit: ‖PA − LU‖/‖A‖ = 4.487946226420872·10⁻¹⁶ in both. One of them moves 41,332 words between fast and slow memory and the other moves 19,476. 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. 10²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1κ₂(A)relative errorforward, A⁻¹bforward, A\bbackward, A⁻¹bbackward, A\bat κ = 10¹⁴η, LU solve2.2·10⁻¹⁷η, via the inverse4.5·10⁻⁵forward, LU solve2.8·10⁻⁴forward, via the inverse0.015one factorisation, two ways to use itand one of them forfeits the backward error The inverse that is never formed x = A⁻¹b is how the solution of a linear system is written and it is not how it is computed. The usual reason given is cost — three times the arithmetic. The real reason is that one of the two routes is backward stable and the other is not, and at κ = 10¹⁴ they differ by twelve orders of magnitude in the number that says whose fault a wrong answer is.

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.

Asymptotic analysisBackward errorBlock methodsData movementFlop countHierarchical matrixLU factorisationOptimal complexityTriangular solveWoodbury identity