Numerical Linear Algebra
where the algebra stops being the arithmetic
Randomised, and the guarantee that changes kind

Built from products alone

A 512-square hierarchical representation, at a relative error of 4·10⁻⁷, from 256 applications of an operator that is never assembled. The compression route reads 262,144 entries; this one reads none, and pays for it with a factor of seven against the representation the entries would have given.
34 min read 24 figures Measured, not assertedTwo routes to a numberAssertions that reject

Worth reading first: The dimension does not appear · An operator with no entries · A block nobody can call sparse.

an-operator-with-no-entries ends on a sentence this essay contradicts: an operator that can only be applied is an operator that cannot be factorised, and every solver in that essay is an iteration for that reason.

A hierarchical representation is the exception, and the reason is structural rather than clever. Every block of it is a submatrix, and a submatrix can be sampled.

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. 1 Applications of an operator that is never assembled, against the entries of the matrix it stands for. One line adds a constant per doubling and the other doubles twice.

The sampling

Take a block A_{τσ} — rows in cluster τ, columns in cluster σ. Build a random matrix Ω that is nonzero only on the rows σ, and apply A. Then

(A Ω)|τ = A{τσ} Ω_σ ,

which is the block applied to a random matrix, which is exactly what the randomised field’s range finder turns into an orthonormal basis for the block’s range.

Once a basis Q for the range is in hand, the second factor needs one more product. Put Q’s columns into n-vectors supported on τ, apply A, and read the rows σ:

(A Q)|σ = A{στ} Q = A_{τσ}ᵀ Q (for symmetric A),

so A_{τσ} ≈ Q (A_{τσ}ᵀQ)ᵀ. Two batches per block and no entry read.

The symmetry is used once and deliberately. On a non-symmetric operator the second step needs a product with Aᵀ instead, the count doubles, and nothing else changes — which is worth knowing and is why it is said rather than assumed.

The count

The blocks at one level of the tree have disjoint column supports, so one batch serves all of them at once. Two batches a level for the two triangles, two more for the projections, and one batch of leaf products recovers every diagonal block in the matrix together, because the leaves are disjoint too. So

products = leaf + 2(2k + p)·levels ,

with k the rank asked of each block and p the oversampling. At leaf 16, k 10, p 4 and four levels that is 16 + 2(24)(4) = 208, and the counter inside the operator says 208.

n levels products entries products per unknown error
64 2 112 4,096 1.75 8.2·10⁻¹¹
128 3 160 16,384 1.25 5.6·10⁻⁹
256 4 208 65,536 0.81 6.7·10⁻⁸
512 5 256 262,144 0.50 4.0·10⁻⁷

Forty-eight more products per doubling, which is 2(2k + p) exactly, and the entries column squares. Products per unknown falls at every doubling: the construction gets cheaper per degree of freedom as the problem gets bigger, which is the same shape the storage claim has and for the same reason.

The subtraction that makes it correct

There is one step in the construction whose omission produces an answer that looks like an answer, and it is worth the section because it cost this build a day.

One batch samples every block at a level at once, because their column supports are disjoint. But reading rows τ of the result picks up the contribution of every cluster in the batch’s support, not only the one paired with τ — and those extra contributions are the blocks at coarser levels, which have already been recovered.

So each batch has the coarser levels peeled off it before the range is taken. Without that, the representation comes out at a relative error of 0.25 with every intermediate quantity looking correct: the bases are orthonormal to 8·10⁻¹⁶, the projections are projections, the counts are right, and the answer is wrong.

There is a second version of the same mistake one level finer, and it is the one that survives a first fix. A block recovered by the other triangle of the same level is the transpose of the one being recovered now — and the projection step reads exactly it — so peeling against it subtracts the answer. The subtraction has to be against completed levels only, which is one line and is the difference between 0.25 and 4·10⁻¹³.

Both are instances of a shape this collection keeps meeting: an error whose symptom is a plausible number. Every gate this site has asks whether a quantity is finite, whether a label fits, whether an assertion held. None of them asks whether the thing being sampled is the thing intended, and the only check that caught this was a comparison against the representation the entries give.

What each round of sketch-and-deflate actually sees, for a sketch kept, half redrawn, and redrawnA randomised range finder forms Y = AΩ, takes Q spanning it, and deflates: A ← A − QQᵀA. With the same Ω the next sketch is A₂Ω = AΩ − QQᵀ(AΩ) = Y − QQᵀY, which is zero exactly, because Q spans Y. Measured, the kept sketch sees 2.62 at the first round and 9.03·10⁻¹⁵ at the second — the unit roundoff, which is to say nothing at all, and the subspace it produces after that is whatever the rounding error happened to point at. A fresh sketch stays at about 3.3, and redrawing half the columns sees about half as much. After 5 rounds of identical arithmetic the kept run is at 0.2845 and the redrawn one at 0.1126.012345610⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹round‖AΩ‖ ⁄ ‖A‖ at that roundredrawn each roundhalf redrawnkept: the sketch is now the zero matrixthe randomness is spentkept, round 12.6kept, round 29·10⁻¹⁵redrawn, round 23.1kept, final error0.28redrawn, final error0.11a random matrix used twiceis not random the second time
Fig. 2 The closest relative of this mistake on the site, from the randomised field: a sketch applied to a matrix it has already been used on, which sees the zero matrix and reports a subspace anyway.

What it costs

The construction is not free and the price is on the same figure. Against the representation the same partition would give if every block were decomposed at the same rank:

n from products from entries excess
64 8.21·10⁻¹¹ 2.91·10⁻¹¹ 2.8×
128 5.61·10⁻⁹ 9.59·10⁻¹⁰ 5.8×
256 6.65·10⁻⁸ 1.01·10⁻⁸ 6.6×
512 4.03·10⁻⁷ 5.49·10⁻⁸ 7.3×

A factor of three to seven, widening slowly with the depth. The widening is the peeling: every level’s data has the coarser levels’ approximations subtracted from it rather than the coarser levels themselves, so a coarse error propagates into every finer level.

That makes this the one accumulation in this whole phase that accumulates. The essay next door measures ninety-eight truncations inside a factorisation costing nothing measurable; this measures five levels of peeling costing a factor of seven, and the difference is that one of them is a projection onto a fixed target and the other is a subtraction of an approximation from a measurement.

Seven is a factor rather than an order of magnitude, which is what keeps the route usable. It is also the thing to watch on a deeper tree, and it is the argument for the variants of this construction that re-sample rather than peel.

What the black-box construction costs: the rank has to be chosen before anything is knownThe compression route is handed an accuracy and returns whatever rank that costs. This one is handed a rank and returns whatever accuracy that buys, because the random matrix has to be drawn before a single entry of the block has been seen. The curve is a straight line at -0.74 decades a column — from 0.0407 at k = 2 to 1.72·10⁻⁹ at k = 12 — so the guess is a guess about a number of digits, and being two columns short costs about a decade and a half. The lower line is the best representation of the same rank, so the vertical gap is what never seeing the matrix costs: 2.9, 4.5, 4.1, 4.4, 6.6, 8.1 times, widening slowly as the rank grows, because a sample of a block whose spectrum has already fallen off a cliff is sampling noise.0246810121410⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank asked of every block‖A − A_H‖ ⁄ ‖A‖built from productsthe best of the same rankan accuracy for a rankerror at k = 20.041error at k = 121.7·10⁻⁹decades a column0.74products at the top240excess over the optimum8.1one route asks for digitsand the other asks for columns
Fig. 3 The other half of the price, and the more awkward one: the rank has to be chosen before anything about the block is known.

The trade the compression route does not make

compress is handed an accuracy and returns whatever rank it costs. This is handed a rank and returns whatever accuracy it buys, because the random matrix has to be drawn before a single entry of the block has been seen.

k error products
2 4.07·10⁻² 80
4 2.13·10⁻³ 112
6 6.19·10⁻⁵ 144
8 1.83·10⁻⁶ 176
10 6.65·10⁻⁸ 208
12 1.72·10⁻⁹ 240

A straight line at 0.73 decades a column. So the guess is a guess about a number of digits, being two columns short costs about a decade and a half, and there is no cliff — which is the same flatness the hierarchy field’s second essay measures from the other direction and is what makes guessing tolerable.

It is the same trade the randomised field already has an essay about. the-sketch-that-is-not-the-answer measures what a sketch buys and what it gives up, and the answer there is that the sample size is chosen against a target the sample cannot see. Here the quantity being guessed is a storage decision rather than an answer, which is a milder failure: guess low and the representation is inaccurate, which is visible; guess high and it is expensive, which is also visible. Neither is silent.

A randomised rank-k solve, 4 seeds a rankRelative error against the rank kept, on a logarithmic vertical axis, with a vertical bar at each rank spanning the seeds. The band is 1.68 wide at rank 8, where the method is at its worst, and 1.004 wide at rank 20, where it is at its best. The same computation on the same data returns a different answer each time.0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.7spread at rank 201best median error0.15widest where the method is worstand the bound does not say so
Fig. 4 The randomised field’s version of the same guess, from the essay that measures it: a target rank chosen before the spectrum is known.

The oversampling, and why the figure is a band

Every number in this construction is a draw, so a line would be a claim about one run.

p mean error spread excess over the optimum
0 9.02·10⁻⁷ 44% 11.6×
1 6.34·10⁻⁷ 17% 8.2×
2 5.33·10⁻⁷ 21% 6.9×
4 3.54·10⁻⁷ 13% 4.6×
8 2.31·10⁻⁷ 12% 3.0×
16 1.96·10⁻⁷ 7% 2.5×

Both columns matter and only one of them is visible in a single run. Oversampling buys a better answer and a more predictable one, and the second is the half a line would hide — at p = 0 the spread across five seeds is 44 per cent of the mean, which means a code that ran once and reported its error has reported a number that could have been half again as large.

And no draw, at any p, beats the flat line. That is the honesty check: the flat line is the best representation of that rank, a sample cannot beat the decomposition it approximates, and an experiment in which one did would be an experiment with a bug in it.

What extra columns buy a construction that never sees the matrix, over five seedsThe band is the range across five seeds and the line is their mean; the flat lower line is the best representation of rank 8 there is, which a decomposition of every block would find. With no oversampling the construction is 11.6 times behind it and the spread across seeds is 44 per cent of the mean. Eight extra columns bring it to 3.0× and the spread to 12 per cent. Both halves matter and only one of them is visible in a single run: what oversampling buys is a better answer and a more predictable one, and a figure drawn from one seed would have shown the first and hidden the second. No draw is ever better than the flat line, which is the check that the comparison is honest — a sample cannot beat the decomposition it is approximating.048121610⁻⁸10⁻⁷10⁻⁶10⁻⁵extra columns in the sample, p‖A − A_H‖ ⁄ ‖A‖the best rank-8 representationfive draws, mean and rangea band, not a lineexcess at p = 012excess at p = 83spread at p = 00.44spread at p = 80.12the optimum of this rank7.8·10⁻⁸one seed shows the meanand five show the risk
Fig. 5 The band, over five seeds, with the optimum of the same rank as a flat line beneath it.
What extra columns buy a construction that never sees the matrix, over five seedsThe band is the range across five seeds and the line is their mean; the flat lower line is the best representation of rank 12 there is, which a decomposition of every block would find. With no oversampling the construction is 50.0 times behind it and the spread across seeds is 146 per cent of the mean. Eight extra columns bring it to 3.4× and the spread to 16 per cent. Both halves matter and only one of them is visible in a single run: what oversampling buys is a better answer and a more predictable one, and a figure drawn from one seed would have shown the first and hidden the second. No draw is ever better than the flat line, which is the check that the comparison is honest — a sample cannot beat the decomposition it is approximating.048121610⁻¹²10⁻¹¹10⁻¹⁰10⁻⁹10⁻⁸extra columns in the sample, p‖A − A_H‖ ⁄ ‖A‖the best rank-12 representationfive draws, mean and rangea band, not a lineexcess at p = 050excess at p = 83.4spread at p = 01.5spread at p = 80.16the optimum of this rank9.5·10⁻¹²one seed shows the meanand five show the risk
Fig. 6 At rank twelve, where the whole picture has slid down and the shape of it has not changed.

Why the blocks can be sampled at all

The step that makes this work is easy to skip past and it is the one that distinguishes a hierarchical matrix from every other structured format on this site.

Every block of the partition is a submatrix of A. Not a transform of one, not a projection of one — the literal set of entries A_{ij} for i in τ and j in σ. So a product with A, restricted to the right rows and columns, is a product with the block, and the randomised field’s machinery applies without modification.

Nothing else in this collection has that property. A circulant’s structure lives in a transform, and it cannot sample the diagonal of F A F* by applying A to a vector supported on part of the index set, because the transform mixes everything. A sparse factorisation’s structure lives in an elimination order, and the factors are not submatrices of anything. A Krylov method’s structure lives in a subspace that does not exist until the products have been done.

That is why this essay exists in this field and not in one of those. The randomised range finder is a general tool for finding the column space of a matrix that can be applied; what a hierarchical partition supplies is a set of matrices that can be applied by restricting products with A, which is the one ingredient the tool needs and almost nothing has.

There is one place the property is nearly lost and it is worth flagging. The peeling above works because a block at a coarse level is also a submatrix, so its action can be subtracted exactly. In the strongly admissible partition — where a level’s off-diagonal part is a mosaic of blocks at several sizes rather than one clean pair — the same subtraction is possible and considerably more bookkeeping, which is a third reason, after storage and after the solve, that the weakly admissible partition is the one codes build.

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. 7 The partition this construction samples, from the hierarchy field: three objects at the top level, and the two dark squares are the same picture one level down.
The diagonals of T, and what each circulant puts in their place (n = 32, ρ = 0.9)Entry value against the diagonal index. The Toeplitz matrix's own diagonals fall geometrically as ρ^k. The wrapped circulant follows them to the halfway point and then jumps to the far diagonal, discarding what it passed. The averaged one blends the two throughout, weighted by how many entries each contributes — and the result is nearer to T in the Frobenius norm, 0.2765 against 0.3123.048121620242800.250.50.751diagonal indexentryT's diagonalswrappedaveragedtwo routes to one circulantformula against the fit8.9·10⁻¹⁶‖C − T‖/‖T‖, averaged0.28‖C − T‖/‖T‖, wrapped0.31one discards a diagonalthe other weights both
Fig. 8 And a structured format whose blocks are not submatrices, from the structure field, where no amount of applying the operator to restricted vectors recovers the structure.

What this makes possible

The construction is the piece that turns everything else in this phase from an observation into a method, and it is worth being explicit about why.

The hierarchy field’s storage and solve results are all measured on representations built by decomposing every off-diagonal block. That reads every entry, costs O(n³) in the worst case, and on the problems this format exists for is the one thing nobody can afford. A field whose central object can only be built that way has described a compression scheme for matrices that have already been formed.

With this, the input is a routine. A boundary element code that can apply its operator by any means — a fast multipole method, a coarse quadrature, another hierarchical matrix — can hand that routine over and get back a representation it can factorise. That is the step from an iteration to a direct method, and it is why this construction is the one part of the phase with a whole essay in a different field.

Error of a difference quotient for J(x)v against ε, on the Bratu problem at n = 64The Jacobian of this problem is a formula, so the error of each quotient is measured against a derivative that is exact rather than against a better quotient. The forward difference falls with a slope of 1.01 — first order — reaches 1.28·10⁻¹⁰ at ε = 10^-6, and rises again with a slope of -1.00 as cancellation takes over. The central difference falls with a slope of 2.00 and bottoms at 1.11·10⁻¹² for two residual evaluations instead of one. Neither gets near the unit roundoff at any ε.10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice
Fig. 9 The field this essay is answering, from the essay that named the problem: an operator that is a subroutine, and a derivative that is accurate to ten digits instead of sixteen.
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. 10 What the representation is for once it exists, from the cost field: a direct solve whose exponent falls towards one.

The three quantities a code has to choose

Nothing in this construction is adaptive, and a reader deciding whether to use it needs the three numbers it asks for and what each of them costs to get wrong.

The rank k. Chosen before anything is known, at 0.73 decades a column. Too low and the representation is inaccurate by a knowable amount; too high and it is expensive by a knowable amount. Neither failure is silent, which is the property that makes the guess acceptable, and an adaptive version — sample a few columns, check whether the residual has fallen enough, sample more — exists and costs an extra batch per block per round.

The oversampling p. Chosen for the spread rather than for the mean. Four is the conventional answer and this measurement supports it: p = 4 costs 24 more products at the largest size and takes the spread across seeds from 44 per cent to 13.

The leaf size. It sets the levels, and the levels appear in the product count linearly — but it also appears additively, as one batch of leaf products for the diagonal blocks. At 512 unknowns with a leaf of 16 that additive term is 16 of the 256 products; at a leaf of 64 it would be 64 of 208. So there is an optimum in the count, it is shallow, and it sits near the leaf sizes a code would choose for other reasons anyway.

The thing none of the three controls is the factor of seven against the decomposition. That is the peeling, it grows with the depth, and no setting of any of the three parameters removes it — which is the honest form of what never seeing the matrix costs.

What extra columns buy a construction that never sees the matrix, over five seedsThe band is the range across five seeds and the line is their mean; the flat lower line is the best representation of rank 4 there is, which a decomposition of every block would find. With no oversampling the construction is 15.0 times behind it and the spread across seeds is 50 per cent of the mean. Eight extra columns bring it to 2.9× and the spread to 24 per cent. Both halves matter and only one of them is visible in a single run: what oversampling buys is a better answer and a more predictable one, and a figure drawn from one seed would have shown the first and hidden the second. No draw is ever better than the flat line, which is the check that the comparison is honest — a sample cannot beat the decomposition it is approximating.048121610⁻⁴10⁻³10⁻²10⁻¹extra columns in the sample, p‖A − A_H‖ ⁄ ‖A‖the best rank-4 representationfive draws, mean and rangea band, not a lineexcess at p = 015excess at p = 82.9spread at p = 00.5spread at p = 80.24the optimum of this rank2.6·10⁻⁴one seed shows the meanand five show the risk
Fig. 11 The second of the three at the narrowest rank measured, where the band is at its widest and the case for oversampling is at its strongest.
160 single-probe estimates of one 40×40 trace, from the two standard probe distributionsTwo clouds of points, each one estimate of the trace from one random vector. Both are centred on the true value of 171.14 — the estimator has no bias at all — and the dashed bands are the standard deviation each distribution is known in closed form to have: 59.14 for the ±1 probe and 72.01 for the normal one. Measured over these draws they come out 56.68 and 68.57.027548110813544119.306194.613269.919345.226420.532drawzᵀAz from one probenormal±1two routes to one spreadthe trace171±1 spread, predicted59±1 spread, measured57normal ÷ ±11.2no bias in either cloudand one of them is narrower for a reason
Fig. 12 And the same decision in the trace field, from the essay that measured it: how many probes, chosen against a variance rather than against an answer.

The refusal

The claim under test is the sentence an-operator-with-no-entries ends on and this essay contradicts: that a matrix that can only be applied is a matrix that cannot be factorised.

The refusal fed to the assertions is the opposite mistake, and it is the one worth guarding against here. A construction that produced a better representation than the decomposition of the same rank would be a construction with a bug in it — a sample cannot beat the optimum it is sampling — so the assertion that the sampled error is below the decomposition’s is fed 6.65·10⁻⁸ against 1.01·10⁻⁸, and it fails.

The claim above is refuted by the table rather than by an assertion, which is worth saying plainly. Two hundred and fifty-six products, a representation at 4·10⁻⁷, and no entry of the matrix ever formed. What the assertions guard is the direction of the comparison, which is where an experiment of this kind goes wrong.

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 80, 112, 144, 176 at n = 64, 128, 256, 512 — 32 more per doubling, which is a logarithm. The other line is n², the entries the compression route reads. At n = 512 that is 176 products against 262,144 entries, and the representation it produces is within 8.0× 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 = 512176entries at n = 5122.6·10⁵per doubling32‖A − A_H‖ ⁄ ‖A‖2.7·10⁻⁴excess over the compression8no entry of the matrixwas ever read
Fig. 13 At a rank of six, where the product count is lower and the error is four decades worse.
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 144, 208, 272, 336 at n = 64, 128, 256, 512 — 64 more per doubling, which is a logarithm. The other line is n², the entries the compression route reads. At n = 512 that is 336 products against 262,144 entries, and the representation it produces is within 6.5× 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 = 512336entries at n = 5122.6·10⁵per doubling64‖A − A_H‖ ⁄ ‖A‖4·10⁻¹⁰excess over the compression6.5no entry of the matrixwas ever read
Fig. 14 And at fourteen, where it is the tightest representation this construction reaches and still costs under 400 products at the largest size.
What the black-box construction costs: the rank has to be chosen before anything is knownThe compression route is handed an accuracy and returns whatever rank that costs. This one is handed a rank and returns whatever accuracy that buys, because the random matrix has to be drawn before a single entry of the block has been seen. The curve is a straight line at -0.90 decades a column — from 0.0328 at k = 2 to 3.43·10⁻¹¹ at k = 12 — so the guess is a guess about a number of digits, and being two columns short costs about a decade and a half. The lower line is the best representation of the same rank, so the vertical gap is what never seeing the matrix costs: 2.9, 3.9, 5.3, 4.7, 5.8, 3.6 times, widening slowly as the rank grows, because a sample of a block whose spectrum has already fallen off a cliff is sampling noise.0246810121410⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹rank asked of every block‖A − A_H‖ ⁄ ‖A‖built from productsthe best of the same rankan accuracy for a rankerror at k = 20.033error at k = 123.4·10⁻¹¹decades a column0.9products at the top184excess over the optimum3.6one route asks for digitsand the other asks for columns
Fig. 15 The guess, at half the size, where the same slope holds and the whole line sits a decade lower.

What the construction is placed against

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 96, 136, 176, 216 at n = 64, 128, 256, 512 — 40 more per doubling, which is a logarithm. The other line is n², the entries the compression route reads. At n = 512 that is 216 products against 262,144 entries, and the representation it produces is within 7.1× 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 = 512216entries at n = 5122.6·10⁵per doubling40‖A − A_H‖ ⁄ ‖A‖10⁻⁵excess over the compression7.1no entry of the matrixwas ever read
Fig. 16 At a rank of eight, where the product count is lower and the representation is two decades worse.
What extra columns buy a construction that never sees the matrix, over five seedsThe band is the range across five seeds and the line is their mean; the flat lower line is the best representation of rank 6 there is, which a decomposition of every block would find. With no oversampling the construction is 13.8 times behind it and the spread across seeds is 16 per cent of the mean. Eight extra columns bring it to 2.7× and the spread to 23 per cent. Both halves matter and only one of them is visible in a single run: what oversampling buys is a better answer and a more predictable one, and a figure drawn from one seed would have shown the first and hidden the second. No draw is ever better than the flat line, which is the check that the comparison is honest — a sample cannot beat the decomposition it is approximating.048121610⁻⁶10⁻⁵10⁻⁴10⁻³extra columns in the sample, p‖A − A_H‖ ⁄ ‖A‖the best rank-6 representationfive draws, mean and rangea band, not a lineexcess at p = 014excess at p = 82.7spread at p = 00.16spread at p = 80.23the optimum of this rank5·10⁻⁶one seed shows the meanand five show the risk
Fig. 17 The band at rank six, where the whole picture is higher and its shape is unchanged.
What extra columns buy a construction that never sees the matrix, over five seedsThe band is the range across five seeds and the line is their mean; the flat lower line is the best representation of rank 10 there is, which a decomposition of every block would find. With no oversampling the construction is 15.7 times behind it and the spread across seeds is 85 per cent of the mean. Eight extra columns bring it to 3.0× and the spread to 5 per cent. Both halves matter and only one of them is visible in a single run: what oversampling buys is a better answer and a more predictable one, and a figure drawn from one seed would have shown the first and hidden the second. No draw is ever better than the flat line, which is the check that the comparison is honest — a sample cannot beat the decomposition it is approximating.048121610⁻¹⁰10⁻⁹10⁻⁸10⁻⁷extra columns in the sample, p‖A − A_H‖ ⁄ ‖A‖the best rank-10 representationfive draws, mean and rangea band, not a lineexcess at p = 016excess at p = 83spread at p = 00.85spread at p = 80.049the optimum of this rank9.6·10⁻¹⁰one seed shows the meanand five show the risk
Fig. 18 And at ten, where it is lower and still narrows in the same way.
What the black-box construction costs: the rank has to be chosen before anything is knownThe compression route is handed an accuracy and returns whatever rank that costs. This one is handed a rank and returns whatever accuracy that buys, because the random matrix has to be drawn before a single entry of the block has been seen. The curve is a straight line at -1.08 decades a column — from 0.0177 at k = 2 to 2.8·10⁻¹³ at k = 12 — so the guess is a guess about a number of digits, and being two columns short costs about a decade and a half. The lower line is the best representation of the same rank, so the vertical gap is what never seeing the matrix costs: 2.2, 2.2, 4.1, 3.5, 2.8, 3.2 times, widening slowly as the rank grows, because a sample of a block whose spectrum has already fallen off a cliff is sampling noise.0246810121410⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²rank asked of every block‖A − A_H‖ ⁄ ‖A‖built from productsthe best of the same rankan accuracy for a rankerror at k = 20.018error at k = 122.8·10⁻¹³decades a column1.1products at the top128excess over the optimum3.2one route asks for digitsand the other asks for columns
Fig. 19 The guess at the smallest size measured, where the same slope holds two decades down.
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. 20 The partition this construction samples, whose blocks are submatrices — which is the property that makes any of it possible.
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. 21 What the representation costs once it exists, from the hierarchy field.
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. 22 And what its accuracy is, which is a backward error whether the representation was compressed or sampled.
‖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 The operation that makes a built representation worth having: a factorisation performed inside the format.
The randomised SVD against the optimum it cannot beatA semi-logarithmic plot of approximation error against target rank. A shaded band shows the spread across seeds, a solid line the optimal error from the exact singular values, and a dashed line the published probabilistic bound well above both.04812162010⁻¹10⁻⁰.⁵1target rank k‖A − A_k‖₂published boundrandomisedσ_{k+1}, optimalhow far apart the three areworst seed spread1.6bound / median at k = 125.9median / optimum at k = 121.960×60, 6 seeds, oversampling p = 5band is best to worst
Fig. 24 And the honest half of every claim in this field, from the essay that established it: a method that exploits structure returns nothing on a matrix with none.

Where this goes next

10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice An operator with no entries At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all. 10²10².³10².⁵⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁶10⁻¹10⁻⁰.⁵1rows in the sketchworst relative distortiondimension 64dimension 2565 seeds per point, band is best to worstthe dimension does not appear The dimension does not appear A random projection preserves the lengths of a set of vectors to within a distortion that depends on how many vectors there are and not on how many coordinates each one has. That is the fact the whole field rests on, and it is genuinely surprising. 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 The sketch that is not the answer Sketch-and-solve throws away the original problem and keeps the small one's answer, which is why its answer moves with the seed. Use the same sketch as a preconditioner instead and the condition number the iteration sees is the same number at every κ from a hundred to ten billion — identically the same, to nine digits, because the spectrum cancels out of it.

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.

AdmissibilityBlack box constructionHierarchical matrixMatrix-freeOff diagonal rankOversamplingRandom projectionRandomised SVDRange finderSubspace embedding