Numerical Linear Algebra
where the algebra stops being the arithmetic
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.
36 min read 25 figures Measured, not assertedAssertions that rejectWhere the algebra stops being the arithmetic

Worth reading first: Where the format starts paying · What a float can hold · 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.

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.

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-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. 3 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.

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

The previous phase of this collection measured three unrelated accumulations and found all three behaving 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.

Three accumulations, each against its own quantity, each divided by its own first pointA left-to-right sum against the number of terms, at a fitted slope of 0.486; a chain of rotations against the number of steps, at 0.554; and a conjugate gradient residual recurrence against the largest iterate, at 0.507. The three share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against, drawn here as the upper line, and each comes out at half of it. Nothing is rescaled except the division by each series' own first point, which is what makes three quantities of different sizes comparable in slope and in nothing else.110¹10²10³110¹10²10³the accumulating quantity, relative to its first valuethe error, relative to its first valuethe bounds: slope 1what all three do: slope ½three mechanisms, one exponenta left-to-right sum0.49a chain of rotations0.55a residual recurrence0.51every bound's slope1spread of the three0.067a bound is a sum of the roundingsand the roundings have signs
Fig. 4 The three walks the previous phase measured, from the essay that established them: three mechanisms sharing no arithmetic, all at about a half, against bounds that are all linear.
Every column's 1-norm on a 24×24 matrix built to stop Hager's walk one column short24 bars, one per column, each the sum of the absolute values in that column — the quantity the estimator maximises, whose largest is the matrix's 1-norm. The walk visits 1 of them and returns 24.00; column 2 has 1-norm 474.24 and is never visited, because its entries alternate in sign and cancel against every sign vector the walk stands on.the estimator maximises this quantity over the columns it visitscolumn 1 ‹visited›24column 2 ‹the answer›474.2column 322.8column 422.8column 522.8column 622.8column 722.8column 822.8column 922.8column 1022.8column 1122.8column 1222.8column 1322.8column 1422.8column 1522.8column 1622.8column 1722.8column 1822.8column 1922.8column 2022.8column 2122.8column 2222.8column 2322.8column 2422.8estimate 24.0a walk that stopped earlythe estimate returned24the true 1-norm474columns visited1products with the matrix5the walk's own stopping test firedand every column it could see was smaller
Fig. 5 And the longest chain of them, from the same field: three thousand rotations, and an error growing like the square root of the count.

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.

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 − A_H‖⁄‖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.

‖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. 6 The measurement, drawn: the factorisation’s residual and the representation’s own error against how many truncations were performed between them.

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.

The residual of a formatted Cholesky, against the accuracy its blocks were compressed atA leaf of 16 means 34 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.027 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.93, 1.83, 0.88, 0.91, 0.91 — 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 factorisationtruncations34residual at 10⁻⁴2.2·10⁻⁵residual at 10⁻¹⁰9.9·10⁻¹²slope1worst ratio to the representation1.8ten decadesand a slope of one
Fig. 7 And the same factorisation across ten decades of tolerance, tracking at a slope of one, with the ratio to the representation’s error printed at every point.
‖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 6.134·10⁻¹⁴ to 1.282·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.95, 0.90, 0.85, 0.85, never above one. A hundred approximate operations contributed nothing measurable to the answer.02040608010010⁻¹⁴10⁻¹³10⁻¹²truncations performed by the factorisationrelative errorthe representation's own error‖A − LLᵀ‖ ⁄ ‖A‖no decomposition without its residualresidual, 0 truncations6.1·10⁻¹⁴residual, 98 truncations1.3·10⁻¹³representation, deepest1.5·10⁻¹³residual ⁄ representation0.85levels, deepest5a hundred approximate stepsand an exact-looking factorisation
Fig. 8 At twelve digits, where the individual truncations are cutting far further down each block’s spectrum and the two curves still do not separate.

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 this measurement does not reach a tree deep enough to have that many. What it reaches is 98, and at 98 there is nothing.

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. 9 The resource the leaf size actually trades against, from the partition field: numbers stored per unknown, most of which at a large leaf is the dense fringe.
Products with A to return all 2 copies, against the block widthOne bar per block width. Widths below 2 return fewer copies than the eigenvalue has, at every step count up to a basis of half the problem's dimension, and are marked rather than drawn. The cheapest width that works is 2, at 32 products with A; the widest drawn costs 84.an eigenvalue repeated 2 times, in a 200-dimensional problemblock of 1never returns them allblock of 232 products, 16 stepsblock of 351 products, 17 stepsblock of 460 products, 15 stepsblock of 570 products, 14 stepsblock of 684 products, 14 stepsthe narrowest that workscheapest width2its products with A32the widest drawn84narrower than the multiplicity is not slowit is absent
Fig. 10 And the other half of the trade, from the cost field: a block size chosen from a memory hierarchy rather than from a matrix.

The order still matters, slightly

Formatted addition is not associative and not commutative, so the order the terms arrive in changes the answer. This collection has an essay about the scalar version of that question 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.

Relative error of three summation algorithms in binary32A log–log plot of relative error against the number of terms for naive, pairwise and compensated summation, each measured against the exactly rounded sum.10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸
Fig. 11 The scalar version, from the arithmetic field, where the ordering is worth orders of magnitude rather than two per cent and the rule runs the other way.

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 and the answer to it is not in this essay.

The residual of a formatted Cholesky, against the accuracy its blocks were compressed atA leaf of 8 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. 12 The same factorisation at the deepest tree measured, where ninety-eight truncations track the tolerance over ten decades.
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. 13 The one term the budget has, from the error field, which is what this essay licenses carrying over to a factorisation.

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 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. 14 At rank two, the narrowest running object measured, where the excess is still under five per cent over thirty-two truncations.
The singular values of a sum of two rank-8 blocks, and the 8 a truncation has to discardA rank-8 block times a vector is a rank-8 block times a vector. A rank-8 block times a rank-8 block is a rank-8 block. A rank-8 block PLUS a rank-8 block is a rank-16 block, exactly, and the 16 bars here are why: the sum of two 8-dimensional spaces is generally 16-dimensional, and none of the 16 singular values is small. Truncating back to 8 costs 54.4 per cent of the block. Below the 16th 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.914σ3, kept0.899σ4, kept0.842σ5, kept0.827σ6, kept0.773σ7, kept0.749σ8, kept0.681σ9, discarded0.655σ10, discarded0.62σ11, discarded0.589σ12, discarded0.562σ13, discarded0.543σ14, discarded0.499σ15, discarded0.461σ16, discarded0.389σ172.96·10⁻¹⁶σ182.77·10⁻¹⁶the operation that is not closedrank of each term8rank of the sum16truncated back to8cost of the truncation0.54the best there is0.54two planesmake a four-space
Fig. 15 And the operation at rank eight, where sixteen singular values arrive and eight have to go.
‖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.459·10⁻⁶ to 2.226·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.95, 0.89, 1.83, 1.78, 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.5·10⁻⁶residual, 98 truncations2.2·10⁻⁵representation, deepest1.3·10⁻⁵residual ⁄ representation1.8levels, deepest5a hundred approximate stepsand an exact-looking factorisation
Fig. 16 The factorisation at four digits, where each truncation is discarding a per cent of its block and the residual still does not notice.

The measurement at other depths

‖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 6.566·10⁻⁸ to 1.462·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.96, 0.92, 0.88, 0.85, 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 truncations6.6·10⁻⁸residual, 98 truncations1.5·10⁻⁷representation, deepest1.7·10⁻⁷residual ⁄ representation0.85levels, deepest5a hundred approximate stepsand an exact-looking factorisation
Fig. 17 At six digits, where each truncation is discarding a thousandth of its block and ninety-eight of them still do not move the residual.
‖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 3.97·10⁻⁴ to 0.0023, 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.95, 0.94, 0.93, 0.83, 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 truncations4·10⁻⁴residual, 98 truncations0.0023representation, deepest0.0028residual ⁄ representation0.83levels, deepest5a hundred approximate stepsand an exact-looking factorisation
Fig. 18 And at two, where each is discarding a per cent.
The residual of a formatted Cholesky, against the accuracy its blocks were compressed atA leaf of 32 means 10 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.028 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.94, 0.89, 0.92, 0.92, 0.93 — one excursion above one, at 0.01, 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 factorisationtruncations10residual at 10⁻⁴1.1·10⁻⁵residual at 10⁻¹⁰9.6·10⁻¹²slope1worst ratio to the representation0.94ten decadesand a slope of one
Fig. 19 The same factorisation at a shallower tree, where ten truncations replace ninety-eight and the line does not move.
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. 20 The synthetic chain at rank six, where the excess over the optimum is smaller still.
The singular values of a sum of two rank-3 blocks, and the 3 a truncation has to discardA rank-3 block times a vector is a rank-3 block times a vector. A rank-3 block times a rank-3 block is a rank-3 block. A rank-3 block PLUS a rank-3 block is a rank-6 block, exactly, and the 6 bars here are why: the sum of two 3-dimensional spaces is generally 6-dimensional, and none of the 6 singular values is small. Truncating back to 3 costs 56.8 per cent of the block. Below the 6th 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.881σ3, kept0.814σ4, discarded0.696σ5, discarded0.639σ6, discarded0.521σ72.17·10⁻¹⁶σ81.89·10⁻¹⁶the operation that is not closedrank of each term3rank of the sum6truncated back to3cost of the truncation0.57the best there is0.57two planesmake a four-space
Fig. 21 The operation at rank three, where six singular values arrive and three have to go.
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. 22 The ranks a formatted arithmetic is performed on, from the partition field, and the reason a rank costs more than its storage.
The backward error of a hierarchical solve, against the error of the representation it usedSolve A_H x = b exactly and the residual against the matrix that was wanted is b − Ax = (A_H − A)x, so ‖b − Ax‖ ⁄ ‖A‖‖x‖ cannot be anything but the representation's own error. Measured across six accuracies spanning ten decades, the two track at a slope of 0.983 and sit a constant 7.4× apart, which is the difference between a norm of a matrix and a norm of that matrix applied to one vector. The reading is the one this field is built on: the accuracy knob is not an accuracy, it is a backward error chosen in advance. Everywhere else on this site a backward error is something an algorithm produced and somebody then measured; here it is a line in the program, and its size is known before the solve starts.10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³‖A − A_H‖ ⁄ ‖A‖, the representation‖b − Ax‖ ⁄ ‖A‖‖x‖, the solveequala backward error, chosenslope0.98representation at 10⁻⁸1.4·10⁻⁹backward error there1.4·10⁻¹⁰representation at 10⁻¹²1.2·10⁻¹³backward error there2.3·10⁻¹⁴the compression is not an approximationit is a perturbation of the problem
Fig. 23 The one term the error budget has left after this measurement, 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. 24 And the one accumulation in this field that does accumulate, from the randomised field: a construction whose every level peels the approximations above it.
Three accumulations, each against its own quantity, each divided by its own first pointA left-to-right sum against the number of terms, at a fitted slope of 0.486; a chain of rotations against the number of steps, at 0.554; and a conjugate gradient residual recurrence against the largest iterate, at 0.507. The three share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against, drawn here as the upper line, and each comes out at half of it. Nothing is rescaled except the division by each series' own first point, which is what makes three quantities of different sizes comparable in slope and in nothing else.110¹10²10³110¹10²10³the accumulating quantity, relative to its first valuethe error, relative to its first valuethe bounds: slope 1what all three do: slope ½three mechanisms, one exponenta left-to-right sum0.49a chain of rotations0.55a residual recurrence0.51every bound's slope1spread of the three0.067a bound is a sum of the roundingsand the roundings have signs
Fig. 25 The three accumulations this one is not, from the sequence field: three mechanisms sharing no arithmetic, all growing like the square root of the count.

Where this goes next

10¹10²10³10⁴10⁵10⁶10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹number of terms addedrelative error against the exact sumin orderin a treecompensatedbinary32 · terms are 1/icompensated: 3·10⁻⁸ The order they are added in Addition is associative in the algebra and is not associative in the arithmetic. The same million numbers, added in a different order, give answers that differ in the third significant figure — and the fix is not a wider float, it is a different order. 110¹10²10³110¹10²10³the accumulating quantity, relative to its first valuethe error, relative to its first valuethe bounds: slope 1what all three do: slope ½three mechanisms, one exponenta left-to-right sum0.49a chain of rotations0.55a residual recurrence0.51every bound's slope1spread of the three0.067a bound is a sum of the roundingsand the roundings have signs Three walks and one bound A left-to-right sum, a chain of three thousand rotations and a conjugate gradient residual recurrence share no arithmetic and no vocabulary. Each has a standard bound that is linear in whatever it accumulates against. All three come out at a half — 0.486, 0.554 and 0.507 — and nothing is rescaled. 10²10³10⁴10⁻¹⁷10⁻¹⁶10⁻¹⁵10⁻¹⁴10⁻¹³10⁻¹²steps taken‖RᵀR − AᵀA‖ ⁄ ‖AᵀA‖the bound, linear in the steps√k · uevery step safe, the chain notdrift after the run3.9·10⁻¹⁴the bound there3.3·10⁻¹³√k · u there6.1·10⁻¹⁵worst single amplification2.7worst leverage met0.69refactorisations1backward stable onceand three thousand times is a different claim Stable once, and three thousand times A sliding window adds a row and removes one at every step and never looks at the data again. No single step of it amplifies by more than 2.72, no downdate fails, and after three thousand steps the triangular factor in memory is 3.9·10⁻¹⁴ from the matrix it is supposed to be a factor of — six hundred times growth from a per-step bound that says nothing about chains.

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