Neither sparse nor dense

The rounding that was not the problem

A rank-k block plus a rank-k block is a rank-2k block, exactly, so every arithmetic in this format truncates after every addition. A Cholesky performed inside it does ninety-eight of those and its residual is 1.14·10⁻⁹ against a representation error of 1.40·10⁻⁹ — the roundings cost nothing measurable.

Worth reading first: Where the format starts paying · The same program, twice · A factorisation with nothing to pivot for.

Everything this field has done so far reads a matrix and writes a representation. A representation that can only be read is a storage format; a representation that can be computed with is a numerical method, and the difference between them is one operation.

The singular values of a sum of two rank-4 blocks, and the 4 a truncation has to discardA rank-4 block times a vector is a rank-4 block times a vector. A rank-4 block times a rank-4 block is a rank-4 block. A rank-4 block PLUS a rank-4 block is a rank-8 block, exactly, and the 8 bars here are why: the sum of two 4-dimensional spaces is generally 8-dimensional, and none of the 8 singular values is small. Truncating back to 4 costs 61.9 per cent of the block. Below the 8th the values are the unit roundoff, which is the check that the doubling is exact rather than approximate. Every product, every factorisation and every Schur complement inside this format is a chain of these, and there is nothing else to do: without the truncation the ranks double at every level and the format is dense by the bottom.σ ⁄ σ₁ of the sum, 64 × 64σ1, kept1σ2, kept0.922σ3, kept0.812σ4, kept0.785σ5, discarded0.778σ6, discarded0.758σ7, discarded0.658σ8, discarded0.571σ92.63·10⁻¹⁶σ102.3·10⁻¹⁶the operation that is not closedrank of each term4rank of the sum8truncated back to4cost of the truncation0.62the best there is0.62two planesmake a four-space
Fig. 1 The operation. Two rank-four blocks, their sum, and the four singular values a truncation has to discard — none of which is small.

The operation the format is not closed under

A rank-k block applied to a vector is a rank-k block applied to a vector. A rank-k block times a rank-k block is a rank-k block. Both of those are closures and neither is in doubt.

A rank-k block plus a rank-k block is a rank-2k block, and it is exact: concatenate the two U factors side by side, concatenate the two V factors, and the product of the wide pair is the sum. There is nothing approximate about it and nothing to be done about it either, because the sum of two k-dimensional column spaces is generally 2k-dimensional.

So every arithmetic in this format carries a truncation after every addition. Not as a refinement — as the only thing that keeps the format a format, since without it a formatted factorisation doubles its ranks at every level and is dense by the bottom.

The truncation itself is done through the factors and never on the block: a thin orthogonalisation of each factor, a decomposition of the small 2k × 2k core, and the cut applied there. That it gives exactly the same answer as truncating the assembled block is asserted rather than assumed — measured at 1.4·10⁻¹⁶ apart — and it is the only reason the operation is affordable at all. The cut itself is the best approximation there is: the discarded tail is the smallest error any rank-k object could have made, which is what makes the comparison below one against an optimum rather than against another code.

What the standard warning says

The received reading of formatted arithmetic is that these truncations accumulate. Each one commits an error of the size of the discarded tail; a chain of m of them commits m such errors; and a formatted factorisation doing hundreds is therefore an object whose accuracy nobody can bound tightly — the same shape as the bound on three thousand rotations, which counts every step and finds the count is not what happens.

That reading has a bound behind it and the bound is correct. It is also, measured, about a factor of thirty too pessimistic on the objects this format actually produces.

The chain, measured

Add m rank-four terms, truncating back to rank four after each, and compare against the best rank-four approximation of the exact sum — the answer no code has and none could beat.

Two regimes, because the warning is written for one of them and formatted arithmetic runs in the other.

Drifting. Every term is a rank-four object whose column space has rotated a little from the last one’s. That is what a sequence of Schur updates to one block looks like: related objects, arriving one at a time, whose span moves slowly.

Independent. Terms with nothing in common, whose sum leaves the rank-four manifold at the second term and stays off it.

m drifting, excess independent, excess
2 1.000 1.000
4 1.000 1.013
8 1.001 1.031
16 1.007 1.036
32 1.034 1.036

The worst point on either curve is 1.04. A bound linear in the number of truncations would say 32.

What a chain of truncations costs, against the best rank-4 answer, in two regimesBoth curves are the error of a running rank-4 object divided by the error of the best rank-4 approximation of the exact sum — the answer nothing in a real code has and none could beat. The standard reading of formatted arithmetic is that a truncation after every operation must accumulate, so a chain of 32 of them is dangerous. Measured, the worst point on either curve is 1.040×. A naive bound multiplies the per-step error by the number of steps and would put the line at 32. The roundings are doing something — the energy they have discarded, which a truncation can only remove, grows from 1.8·10⁻⁷ to 0.03 across the sweep — they are just not doing it to the answer. What does limit the answer is whether the exact sum is still nearly rank 4, and the two regimes disagree about that by a factor of 1022 at the second term.0816243211.11.21.31.4terms added, each followed by a truncationerror ⁄ best rank-k erroroptimalterms with nothing in commona subspace that driftsthe rounding nobody should have feareddrifting, worst excess1independent, worst excess1a linear bound would say32energy discarded, first1.8·10⁻⁷energy discarded, last0.03thirty-two roundingsand four per cent
Fig. 2 The two regimes on one figure, with the optimum as a horizontal line. Thirty-two truncations, four per cent.

The roundings are doing something

A claim that a quantity is small is not evidence unless something could have made it large, so the same sweep carries a second measurement: the deficit, which is how much of the exact sum’s energy the running object has lost. A truncation can only remove, never add, so the deficit can only grow.

It grows from 1.8·10⁻⁷ at the second term to 3.0·10⁻² at the thirty-second — five orders of magnitude. The truncations are removing a great deal.

What they are not doing is removing anything the best rank-four answer would have kept. The exact sum of thirty-two drifting rank-four terms is not itself rank four — its own best rank-four approximation is 0.199 away — and the running object is 1.034 times that. Almost all of the error is the format, and almost none of it is the rounding.

That distinction is the whole essay and it is the one the standard warning collapses. Successive truncation of a nearly-low-rank object tracks the optimum, and what limits the answer is whether the object was ever nearly low rank.

What a chain of truncations costs, against the best rank-2 answer, in two regimesBoth curves are the error of a running rank-2 object divided by the error of the best rank-2 approximation of the exact sum — the answer nothing in a real code has and none could beat. The standard reading of formatted arithmetic is that a truncation after every operation must accumulate, so a chain of 32 of them is dangerous. Measured, the worst point on either curve is 1.034×. A naive bound multiplies the per-step error by the number of steps and would put the line at 32. The roundings are doing something — the energy they have discarded, which a truncation can only remove, grows from 1.8·10⁻⁷ to 0.031 across the sweep — they are just not doing it to the answer. What does limit the answer is whether the exact sum is still nearly rank 2, and the two regimes disagree about that by a factor of 1109 at the second term.0816243211.11.21.31.4terms added, each followed by a truncationerror ⁄ best rank-k erroroptimalterms with nothing in commona subspace that driftsthe rounding nobody should have feareddrifting, worst excess1independent, worst excess1a linear bound would say32energy discarded, first1.8·10⁻⁷energy discarded, last0.031thirty-two roundingsand four per cent
Fig. 3 The narrowest blocks the slider draws. The worst excess over the optimal truncation is 1.0342× for the running object and 1.0337× for independently truncated ones — the two are the same number to three figures.
What a chain of truncations costs, against the best rank-3 answer, in two regimesBoth curves are the error of a running rank-3 object divided by the error of the best rank-3 approximation of the exact sum — the answer nothing in a real code has and none could beat. The standard reading of formatted arithmetic is that a truncation after every operation must accumulate, so a chain of 32 of them is dangerous. Measured, the worst point on either curve is 1.039×. A naive bound multiplies the per-step error by the number of steps and would put the line at 32. The roundings are doing something — the energy they have discarded, which a truncation can only remove, grows from 1.9·10⁻⁷ to 0.042 across the sweep — they are just not doing it to the answer. What does limit the answer is whether the exact sum is still nearly rank 3, and the two regimes disagree about that by a factor of 1000 at the second term.0816243211.11.21.31.4terms added, each followed by a truncationerror ⁄ best rank-k erroroptimalterms with nothing in commona subspace that driftsthe rounding nobody should have feareddrifting, worst excess1independent, worst excess1a linear bound would say32energy discarded, first1.9·10⁻⁷energy discarded, last0.042thirty-two roundingsand four per cent
Fig. 4 Rank three: 1.0385× and 1.0343×.

The two columns stay within a few per cent of each other and of one, and the order between them is about to swap, which is the part that needs the slider rather than a frame.

What a chain of truncations costs, against the best rank-6 answer, in two regimesBoth curves are the error of a running rank-6 object divided by the error of the best rank-6 approximation of the exact sum — the answer nothing in a real code has and none could beat. The standard reading of formatted arithmetic is that a truncation after every operation must accumulate, so a chain of 32 of them is dangerous. Measured, the worst point on either curve is 1.049×. A naive bound multiplies the per-step error by the number of steps and would put the line at 32. The roundings are doing something — the energy they have discarded, which a truncation can only remove, grows from 1.8·10⁻⁷ to 0.033 across the sweep — they are just not doing it to the answer. What does limit the answer is whether the exact sum is still nearly rank 6, and the two regimes disagree about that by a factor of 1005 at the second term.0816243211.11.21.31.4terms added, each followed by a truncationerror ⁄ best rank-k erroroptimalterms with nothing in commona subspace that driftsthe rounding nobody should have feareddrifting, worst excess1independent, worst excess1a linear bound would say32energy discarded, first1.8·10⁻⁷energy discarded, last0.033thirty-two roundingsand four per cent
Fig. 5 Rank six: 1.0372× drifting against 1.0489× independent. The running object is now the better of the two.
What a chain of truncations costs, against the best rank-8 answer, in two regimesBoth curves are the error of a running rank-8 object divided by the error of the best rank-8 approximation of the exact sum — the answer nothing in a real code has and none could beat. The standard reading of formatted arithmetic is that a truncation after every operation must accumulate, so a chain of 32 of them is dangerous. Measured, the worst point on either curve is 1.058×. A naive bound multiplies the per-step error by the number of steps and would put the line at 32. The roundings are doing something — the energy they have discarded, which a truncation can only remove, grows from 1.7·10⁻⁷ to 0.036 across the sweep — they are just not doing it to the answer. What does limit the answer is whether the exact sum is still nearly rank 8, and the two regimes disagree about that by a factor of 1005 at the second term.0816243211.11.21.31.4terms added, each followed by a truncationerror ⁄ best rank-k erroroptimalterms with nothing in commona subspace that driftsthe rounding nobody should have feareddrifting, worst excess1independent, worst excess1.1a linear bound would say32energy discarded, first1.7·10⁻⁷energy discarded, last0.036thirty-two roundingsand four per cent
Fig. 6 At rank eight, where both curves are flatter still. A wider running object is closer to the exact sum in absolute terms and no closer to optimal in relative ones — 1.0400× against 1.0582×.
k drifting independent difference
2 1.0342× 1.0337× drifting worse by 0.0005
3 1.0385× 1.0343× drifting worse by 0.0042
4 1.0343× 1.0401× drifting better by 0.0058
6 1.0372× 1.0489× drifting better by 0.0117
8 1.0400× 1.0582× drifting better by 0.0182

The worst excess is between 3.4% and 5.8% at every rank, against a linear bound that says 32×. The bound overstates by a factor of about thirty at every stop on the slider — 32 against 1.04 and 32 against 1.058 — which is the essay’s headline turned from one measurement into six.

And the drifting column is the better one from rank four upward, by a margin that widens: 0.0058, 0.0117, 0.0182. That is the opposite of the intuition the linear bound encodes. Truncations do not compound here; carrying a running low-rank object through a sequence of additions ends up closer to the optimal rank-k approximation of the exact sum than truncating each addend independently and adding those. Each truncation discards what the current running object cannot represent, which is information about where the sum is actually going — and that is better-aimed than discarding what each addend cannot represent in isolation.

The rank does not decide the size of the effect. Across a factor of four in k the worst excess moves by two percentage points and the bound does not move at all.

The accumulation that is not a walk, and is not a sum either

Three unrelated accumulations in this collection were measured and all three behaved as random walks: a left-to-right sum at slope 0.486, three thousand Givens and hyperbolic rotations at 0.554, a conjugate gradient residual recurrence at 0.507 — against bounds that are all linear.

The obvious question is which of those this is, and the answer is neither, which is why it needed its own measurement.

A truncation error has a sign in one sense and not in another. In energy it has one: the running object is always smaller than the exact sum, and the deficit is monotone, which is a systematic bias and grows linearly. In direction it does not: the discarded subspace at step i is unrelated to the one at step j, so the errors do not line up and the total does not grow like their sum.

But the total error does not grow like a walk either, and it is worth saying why not. In both walks and sums the accumulated error is a function of the number of steps. Here it is not: the running object at step m is close to the best rank-k approximation of the sum of m terms, and how large that is depends on how far the sum has left the manifold — which for the drifting family goes like m² and for the independent family is nearly constant from the second term.

So the right sentence is not about accumulation at all. A running truncation is a projection onto a moving target, and it tracks the target. Whether the answer is good is a question about the target.

The real object

A synthetic chain of additions is a model of a formatted arithmetic and not one. The object worth measuring is a factorisation, so the field’s own rule applies: no decomposition is drawn without its residual printed, and here the residual is the whole measurement — on a hierarchy with no grid behind it, where the partition is the only thing the accuracy depends on.

A Cholesky performed entirely inside the format. The recursion is the textbook one and every step of it stays in the format: the off-diagonal factor L₂₁ is exact, because a triangular solve against one factor of a rank-k block leaves a rank-k block, and the Schur complement A₂₂ − L₂₁L₂₁ᵀ is a rank-k addition to a hierarchical matrix — which pushes down the tree, hits every admissible block in the subtree, and truncates each of them.

So the number of truncations is a count, not a parameter. It is one per admissible block in the right subtree, at every level, and the leaf size sets it:

leaf levels truncations ‖A − AH‖⁄‖A‖ ‖A − LLᵀ‖⁄‖A‖ ratio
128 1 0 2.12·10⁻¹⁰ 2.12·10⁻¹⁰ 1.00
64 2 2 9.75·10⁻¹⁰ 9.07·10⁻¹⁰ 0.93
32 3 10 1.08·10⁻⁹ 9.84·10⁻¹⁰ 0.92
16 4 34 1.08·10⁻⁹ 9.89·10⁻¹⁰ 0.91
8 5 98 1.40·10⁻⁹ 1.14·10⁻⁹ 0.81

The first row is the control: a tree one level deep performs no truncation at all, so its residual is the representation’s error exactly. The last row performs ninety-eight, and its residual is 0.81 times the representation’s error.

Not larger. Smaller, at every depth, and the ratio falls as the depth grows.

Why the ratio falls rather than rises

The residual of the factorisation is not the sum of the representation’s error and the truncations’. It is a norm of a difference of matrix products, and the two contributions are not aligned.

More usefully: the deeper tree compresses more of the matrix — that is why its representation error is 1.40·10⁻⁹ rather than 2.12·10⁻¹⁰ — and the additional compression is of blocks near the diagonal that carry a large share of the norm. The factorisation of that more-compressed matrix has a residual in proportion. The truncations inside the recursion contribute something, and it is below the level at which this measurement can see it.

The right conclusion is not that formatted arithmetic is exact. It is that on this object, at these depths, the approximate arithmetic contributed less than the approximation it was performed on — and since the approximation it was performed on is a deliberate parameter, the truncations are inside the budget rather than beside it.

What the depth knob is actually for

The leaf size appears in this essay as the thing that sets the truncation count, which is a convenient accident. It is worth saying what it is really for, because the null result above changes the answer.

A larger leaf means a shallower tree: fewer levels, fewer blocks, more of the matrix stored densely. That costs storage — the dense fringe is O(n·leaf) — and buys simplicity and better locality, since a 128 × 128 dense solve at a leaf is a well-behaved kernel and a 8 × 8 one is not.

A smaller leaf means the opposite: more of the matrix compressed, less storage, more levels, worse locality, and a great many more truncations in any formatted operation.

Before this measurement the last of those was a reason to be cautious about small leaves. A code that believed the standard warning would keep its leaves large to keep the truncation count down, and would pay for that in storage on every problem.

After it, the leaf size is a pure storage-against-locality trade with no accuracy term in it at all. That is a genuinely simpler decision and it is the practical thing this essay is for: the depth of the tree can be chosen from the machine.

There is one caveat and it belongs on the record. The five rows above span leaf sizes from 128 down to 8, and the ratio column falls monotonically from 1.000 to 0.812 — it does not rise, but it does move. Extrapolating the count rather than the ratio is what would be dangerous: 98 truncations is not 980, and a single matrix size does not reach a tree deep enough to have that many. The next section goes as far as this measurement can.

Two and a half times as far

The caveat above is an invitation and it is cheap to accept: the truncation count is set by the tree, and the tree has two knobs rather than one. Holding the leaf at 8 and doubling the matrix from 256 to 512 adds a level, and a level near the bottom is where most of the blocks are.

n     leaf   levels   truncations   rep. error    residual    ratio
256   128      1           0        2.12·10⁻¹⁰   2.12·10⁻¹⁰   1.000
256     8      5          98        1.40·10⁻⁹    1.14·10⁻⁹    0.812
512   128      2           2        4.09·10⁻¹⁰   4.03·10⁻¹⁰   0.985
512    32      4          34        1.13·10⁻⁹    9.99·10⁻¹⁰   0.887
512     8      6         258        1.44·10⁻⁹    1.14·10⁻⁹    0.793

Two hundred and fifty-eight truncations, and the residual is still below the representation’s own error — 0.793 of it, which continues the fall rather than turning. Across all ten configurations measured the ratio never exceeds 1.000. A bound linear in the count would read 259 at the deepest row against a measured 0.793, so it is out by a factor of three hundred and twenty-six.

The wider sweep also answers a question the single size could not, which is whether the count is the variable at all. Four truncation counts — 2, 10, 34 and 98 — occur at both matrix sizes, and the ratio at a fixed count moves by up to 9% when n doubles, in both directions. Against that, the whole 129-fold range of counts moves it by 21%. The count is the stronger of the two influences and it is a weak one, which is the null result read from the other side: a quantity that a hundred-fold change in the number of approximate operations moves by a fifth, and that a doubling of the matrix at a fixed count moves by a tenth, is not being driven by the approximate operations.

And it does not depend on where the truncations cut. At the deep tree of 512, across six decades of tolerance:

eps        truncations    ratio at leaf 8    ratio at leaf 128
10⁻⁶           258            0.834               0.969
10⁻⁸           258            0.793               0.985
10⁻¹²          258            0.831               0.962

Six per cent of movement over six decades, with the deep tree below the shallow one at every tolerance. That is the reason the result can be stated as a property of the tree rather than of the cut: change the size of what each truncation discards by a factor of a million and the ratio does not notice.

What remains genuinely unmeasured is one order further out. 258 is not 2,580, and the dense reference this sweep is checked against is what stops it — an answer that is known has to be affordable for the measurement to mean anything, and above n = 512 it is not. The honest form of the caveat is therefore narrower than it was: the null result holds over the range where an exact reference exists, which is two decades of truncation count and six of tolerance, and anything past it is an extrapolation.

The order still matters, slightly

Formatted addition is not associative and not commutative, so the order the terms arrive in changes the answer. The scalar version of that question is elsewhere in this collection and the answer there is smallest-first.

Measured here on twenty-four drifting terms: small-first gives 0.11444, the given order 0.11223, and large-first 0.11159, against a best-possible 0.11020 that is identical in all three because addition is commutative and only the rounding is not.

So the ordering is worth about two per cent, and it goes the opposite way from the scalar rule — large-first is best. That is not a contradiction, because what a truncation costs is not set by the size of the term. It is set by how much of the running object’s energy the new term is orthogonal to, and adding the large terms first builds a running subspace the small ones mostly lie inside.

Two per cent is not a result anybody should act on. It is on the record because the scalar question is already on this site with a different answer, and a reader carrying that answer here would be carrying it in the wrong direction.

Where this leaves the format as a numerical method

Three things follow, and together they are what turns a storage scheme into something a program can compute with.

The arithmetic is closed enough. Addition leaves the format and a truncation returns it, at a cost of four per cent over the best possible. Multiplication, triangular solve and factorisation are all built out of additions, so all of them inherit that.

The error budget has one term in it. Before this measurement a formatted factorisation had two sources of error — the representation and the arithmetic — and no way to say which dominated. After it, at these depths, there is one: the representation, which is a parameter. That is what makes the previous field’s sentence hold for a factorisation and not only for a solve: the compression is a backward error, and the arithmetic performed inside it does not add a second one worth naming.

And the residual is still the site’s rule. The last two columns of the depth table are ‖A − LLᵀ‖⁄‖A‖ for a factorisation that never assembled A, never assembled L, and performed ninety-eight approximate operations on the way. It is printed on every figure in this field for the same reason it is printed everywhere else here, and in this case it is not decoration — it is the entire result.

What is not established is anything about a formatted arithmetic that has to be stable rather than merely accurate. Nothing here pivots, nothing here is indefinite, and the matrix is positive definite by construction with a condition number near twenty. A formatted LU with pivoting on an ill-conditioned matrix is a different question — one where the pivot reads the units before it decides anything — and the answer to it is not in this essay.

The residual of a formatted Cholesky on 256 unknowns, against the accuracy its blocks were compressed atA leaf of 8 on 256 unknowns means 98 truncations inside the factorisation, and the ranks of its blocks run to 4, 7, 9, 12, 14 across the sweep. The residual tracks the tolerance at a slope of 1.043 over ten decades: the knob does what a knob should, and nothing in the approximate arithmetic bends it. The ratio to the representation's own error stays at 0.83, 1.78, 0.85, 0.81, 0.86 — one excursion above one, at 10⁻⁴, where the tolerance sits between two singular values of one block and the rank it takes is a rounding of a decision rather than a measurement.10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²ε the blocks were compressed at‖A − LLᵀ‖ ⁄ ‖A‖the factorisationthe knob, on a factorisationtruncations98residual at 10⁻⁴2.2·10⁻⁵residual at 10⁻¹⁰10⁻¹¹slope1worst ratio to the representation1.8ten decadesand a slope of one
Fig. 7 The same factorisation at the deepest tree measured, where ninety-eight truncations track the tolerance over ten decades.

The refusal

The claim under test is the standard warning: that truncating after every operation makes a long chain of them dangerous.

The assertion that thirty-two truncations cost a factor of two over the best rank-four answer is fed 1.034 and 1.040 — the worst points of the two regimes — and it fails. Fed a bound linear in the count it would pass by a factor of sixteen.

The second refusal is the one that makes it about the real object rather than the model. The assertion that a hundred truncations cost twice what none do is fed the first and last rows of the depth table, 1.000 and 0.812, and it fails in the direction nobody would have predicted.

What is left standing after both is a narrower and more useful claim than the warning, and it is the sentence this essay exists for: the rounding is not what limits a formatted arithmetic. The format is. A code that spends effort on a cleverer truncation is optimising the four per cent; a code that spends it on a partition whose blocks are genuinely low rank is optimising the rest.

What links here

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

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.

Cholesky factorisationEckart–YoungError accumulationFormatted arithmeticHierarchical matrixLow-rank approximationRandom walkRecompressionResidualTruncated SVD