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

A rank that is a number of digits

Ask a kernel block for two digits and it costs two columns; ask for fourteen and it costs nine. The curve is a straight line at 0.55 columns a decade, and the bound the geometry gives is a straight line too — at 3.32, which is the same shape and six times the price.
37 min read 27 figures Measured, not assertedTwo routes to a numberPrecision as the knob

Worth reading first: A block nobody can call sparse · Rank is a decision · The condition number is an amplifier.

The previous essay ended on a sentence that reads like a hedge and is not one: the numerical rank of a block is a property of a matrix and an accuracy together. This essay is what that costs, and the answer is a straight line.

How many columns a decade of accuracy costs, measured and predicted, at q = 0.500The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 1. The measured curve is a straight line at 0.55 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 3.32 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 5.3× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121501122334455digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.55bound, a decade3.3rank at 10⁻⁸5bound at 10⁻⁸28q0.5the shape is rightand the constant is not
Fig. 1 Columns kept against digits asked for, on a 128 × 128 block of the kernel 1/r between two intervals separated by their own width. Two straight lines with different slopes.

The measurement

One block. Two clusters of 128 points, on [0, 1] and on [2, 3], so the separation ratio q = (r_x + r_y)/d is exactly 1/2. Count the singular values above ε·σ₁ for a sequence of ε:

ε columns kept what the geometry promises
10⁻² 2 8
10⁻⁴ 3 15
10⁻⁶ 4 21
10⁻⁸ 5 28
10⁻¹⁰ 6 35
10⁻¹² 7 41
10⁻¹⁴ 9 48

The measured column is a straight line against the number of digits, at 0.554 columns a decade. Two digits cost one column. Fourteen digits cost nine.

That is a different kind of statement from this block is low rank. It says the exchange rate is constant: there is no accuracy at which the price per digit jumps, no cliff, and no regime where buying one more digit is suddenly unaffordable. Everything a code does with this format is a consequence of that flatness, because it means an accuracy can be chosen from what the answer needs rather than from what the storage will bear.

The other line

The right-hand column of the table came from four numbers and no entry of the matrix.

The expansion in the previous essay converges at rate q and its truncation at order p has relative error at most q^(p+1)/(1 − q). Solved for the smallest p reaching ε, that is a predicted rank, computable from the two intervals’ endpoints before anything is evaluated. At q = 1/2 it comes to 3.32 columns a decade — because log₁₀ 2 is 0.301, and one over that is 3.32.

So both quantities are linear in the number of digits, which is the agreement, and the constants differ by a factor of six, which is the disagreement. The bound is an upper bound at every point measured, and it has to be: the decomposition is optimal and the expansion is one particular rank-p matrix.

Three routes to the error of a rank-k approximation of one block, at q = 0.5The upper line is q^(p+1)/(1 − q), which used four numbers about two intervals and no entry of the matrix. The middle line is the truncated expansion actually evaluated — a rank p + 1 matrix written down as a table of powers, with no decomposition anywhere in it — and it falls at -0.36 decades a column, close to the log₁₀ q = -0.30 the bound predicts. The lower line is the decomposition, which is optimal by construction, and it falls at -1.74 — 4.8 times as fast. Two of these curves share no arithmetic. What they agree on is that the error is geometric in the rank; what they disagree on is the base, by a factor in the exponent rather than a constant, and the disagreement is why a partition allocated from the bound is safe and wasteful at the same time.0246810121410⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹rank of the approximationrelative errorthe boundthe expansion, evaluatedthe decompositiontwo routes, one shapebound at rank 50.063expansion at rank 50.0048decomposition at rank 52.8·10⁻⁹written, decades a column0.36best, decades a column1.7one route used the matrixand the other used four numbers
Fig. 2 Where the six comes from, drawn as three error curves. The middle one is the expansion evaluated rather than bounded, and it falls at 0.365 decades a column against the bound’s 0.301 and the decomposition’s 1.74.

Which of the two a code uses, and it is not the smaller one

The obvious reading is that the bound is bad and the measurement is what matters. That is wrong in a way worth spelling out, because it is the reading that produces a program which cannot run.

A hierarchical format has to decide how much memory to allocate for a block before it compresses it, and on the problems this is for it has to decide without forming the block. The bound is the only one of the two numbers available at that moment. It comes from the geometry, it costs four subtractions, and it is an upper bound — which is the safe direction, because a rank allocated too generously wastes memory and a rank allocated too meanly loses accuracy silently.

So the six is not an error in the bound; it is the price of the bound being computable in advance. What the measurement buys is the knowledge that the price is a factor and not an order of magnitude, which is what makes the bound usable at all. A bound loose by a thousand would be a bound nobody allocates from.

The proved bound against the error it bounds, at 24 significand bitsTwo curves against the Hilbert size on a logarithmic vertical axis: the error the computation actually made, which is known because the exact answer is the integers by construction, and the bound proved for it in interval arithmetic. The bound is above the error at every size, by 9.949 to 23.12 times. Past n = 5 the method refuses to return a bound at all.345678910111210⁻⁵10⁻³10⁻¹Hilbert size nerrorrefusedrefusedrefusedrefusedrefusedthe proved boundthe actual errorbound ÷ errorslack at n = 313slack at n = 423slack at n = 59.9the error is known, not estimatedthe exact answer is the integers 1 … n
Fig. 3 The same relationship in the field that has an essay about it: a bound that holds at every point and is attained at none, and what the gap between the two costs whoever budgets from it.

The gap moves both lines, and not by the same amount

The one geometric quantity in the whole story is q, and it is a knob.

Move the two intervals apart and q falls. At a gap of a quarter of an interval width q is 0.8 and the block needs 7 columns at eight digits, against a bound of 90. At a gap of eight it is 0.111, the block needs 3, and the bound says 9.

gap q measured at 10⁻⁸ bound
0 1 12 —
0.25 0.800 7 90
0.5 0.667 6 49
1 0.500 5 28
2 0.333 4 18
4 0.200 4 12
8 0.111 3 9

The bound is nearly thirteen times too pessimistic at the tightest separation and three times at the widest — and the tightest separation is exactly where a partition puts most of its blocks, because a partition that only compressed well-separated pairs would compress almost nothing. The bound is worst where the decision is made.

At a gap of zero the bound does not exist at all. q is 1, the expansion does not converge, and the entry in the table is a dash. That case is the subject of the fifth essay in this field, and it is the reason a partition has a test in it rather than a rule.

How many columns a decade of accuracy costs, measured and predicted, at q = 0.800The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 0.25. The measured curve is a straight line at 0.84 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 10.34 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 12.7× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121503570105140digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.84bound, a decade10rank at 10⁻⁸7bound at 10⁻⁸90q0.8the shape is rightand the constant is not
Fig. 4 The two lines at a separation of a quarter of an interval width, where the bound is nine times the measurement and both are still straight.
How many columns a decade of accuracy costs, measured and predicted, at q = 0.111The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 8. The measured curve is a straight line at 0.32 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 1.00 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 3.0× about the constant, which is the safe direction for a quantity you have to allocate storage from.03691215051015digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.32bound, a decade1rank at 10⁻⁸3bound at 10⁻⁸9q0.11the shape is rightand the constant is not
Fig. 5 And at eight, where they nearly meet. The bound is right about the shape everywhere and about the constant only where the geometry is generous.

What a straight line here means for a program

Three consequences, and the third is the one this field turns on.

A storage budget is a linear function of the digits. Doubling the digits does not double the memory; it adds a fixed number of columns per block. Storage against accuracy is therefore a gentle trade, not a cliff, and a code that runs out of memory can give back digits one at a time.

There is no natural stopping point. Nothing in the curve says stop here. A sparsity pattern says what to keep; a numerical rank does not, so the choice has to come from outside the matrix. This collection has a whole field about problems whose answer is not determined by the data, and this is the same shape of question about a representation rather than about a solution.

And what to set it from is not the matrix. Since the price per digit is constant, the sensible question is how many digits the answer needs — which is a question about the condition number of the problem and about what the result is for. The essay on the compression as a backward error is where that closes, and the sentence it arrives at is: decide the digits the answer needs, divide by κ, and compress to that.

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. 6 Where the choice gets made, from later in this field: the accuracy knob turns out to be a backward error, so the digits it costs and the digits it buys are related by the condition number and by nothing else.

What sets the slope, which is not the kernel

The 0.554 is worth taking apart, because a constant that appears in a table without a mechanism is a number a reader has to take on trust.

The exchange rate is one over the logarithm of the ratio between consecutive singular values. This block’s ratio is about 1/52 — the values run 1, 2.40·10⁻², 4.59·10⁻⁴, 8.47·10⁻⁶ — so log₁₀ 52 is 1.72, and one over it is 0.58. The fitted 0.554 differs from that by the small non-geometric part of the head of the spectrum, which is what a fit over seven points ought to differ by.

So the slope is set by the decay ratio, and the decay ratio is set by the geometry. What it is not set by is which smooth kernel was chosen. Replace 1/r by log r on the identical points and the spectrum is 1, 2.09·10⁻², 4.24·10⁻⁴, and the count at eight digits is the same 5. That is not a coincidence and it is not a general law either: the two kernels are related by a derivative, which preserves the expansion’s structure, and a kernel unrelated to either — one that oscillates, say — behaves completely differently. The fourth essay in this field is about that case and it is the one where all of this stops working.

The reading to carry is narrower than smooth kernels are compressible. It is that the decay rate belongs to the pair of clusters, and the kernel decides only whether an expansion of that kind exists at all. Everything in the table above would be reproduced by any kernel with a convergent separable expansion on those two intervals, and none of it survives a kernel without one.

The singular values of one off-diagonal block, for four kernels on the same 64 pointsTwo intervals that do not touch — [0, 1] and [2, 3] — and the 64 × 64 block between them, for four kernels. 1/r falls a factor of 52 a column and is below 10⁻⁸ after 5; log r behaves the same way, for the reason the expansion makes obvious. An independent draw per entry gives 64 singular values above 10⁻⁸ out of 64, on the same size and the same density, which is what makes the other curves a measurement rather than a property of sorted numbers. The block has full algebraic rank in every case; what differs is where the numbers stop mattering, and that is a decision rather than a fact about the matrix.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, in orderσ ⁄ σ₁eight digitsan independent draw per entry1 ⁄ ra cliff, and a control1/r rank at 10⁻⁸5log r rank at 10⁻⁸5noise rank at 10⁻⁸64σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns
Fig. 7 Two smooth kernels on one pair of clusters, with a control. The two upper curves are 1/r and log r and they sit almost on top of each other; the flat one is what a matrix with no expansion looks like.
Two kernels on identical geometry, as the picture is scaled up at a fixed wavelengthTwo parallel segments of length L a distance 2L apart, so q = 0.50 at every point on this sweep and nothing about the arrangement changes. 1/r needs 6, 6, 6, 6, 6 columns — one number. cos(κr)/r at κ = 40 needs 12, 16, 22, 33, 53, growing like a power of the size with no size at which it stops, because what decides it is the number of wavelengths across the pair — 3.2, 6.4, 12.7, 25.5, 50.9 — and no ratio of lengths can see that. The segments are parallel rather than collinear on purpose: in one dimension cos(κ|x − y|) obeys an addition formula and is exactly rank two, which makes the collinear version of this look like a confirmation of the smooth case and is an artefact of the arrangement.-1012301020304050log₂ of the segment lengthcolumns above 10⁻⁸cos(κr) ⁄ r at κ = 401 ⁄ r, same pointsthe case with no answerq, held fixed0.51/r at L = 0.561/r at L = 86cos(κr)/r at L = 0.512cos(κr)/r at L = 853the geometry did not moveand the rank did
Fig. 8 And the kernel that has one and does not benefit from it, from the essay two ahead: the same geometry, a fixed wavelength, and a rank that grows with the size of the picture.

The one place the line bends

Everything above is measured on an admissible pair — two clusters that do not touch. There is one setting of the gap where the whole account changes character and it is the leftmost row of the second table.

At a gap of zero the two intervals share an endpoint. q is exactly 1, the expansion does not converge, the bound is undefined, and the entry is a dash. The block is still compressible — 12 columns at eight digits out of 128 — but nothing above predicts that number, and the third essay in this field measures what makes it different: the 12 becomes 13 when the block is sampled twice as finely, and the 5 of an admissible pair does not become anything.

That is the distinction a partition is built around. It is not that a touching pair is incompressible; it is that its cost is a function of the discretisation as well as of the geometry, so a code cannot budget for it from the geometry alone. Given the choice between a block whose price is known before the problem is discretised and a block whose price has to be measured afterwards, a format takes the first and refuses the second — which is what the admissibility test in the fifth essay does, in one comparison of four numbers.

How many columns a decade of accuracy costs, measured and predicted, at q = 0.667The number of singular values above ε, against the number of digits ε asks for, on a 128 × 128 block between two intervals separated by a gap of 0.5. The measured curve is a straight line at 0.70 columns a decade: a digit costs the same handful of columns wherever you buy it, which is why the accuracy is a knob and not a cliff. The upper line is what the geometry alone promises — q^(p+1)/(1 − q) below ε, solved for p, using four numbers and no entry of the matrix — at 5.66 columns a decade. Both are straight and they are not the same straight: the bound is right about the shape and loose by 8.3× about the constant, which is the safe direction for a quantity you have to allocate storage from.0369121501938577695digits asked for, −log₁₀ εcolumns keptwhat the geometry promiseswhat the matrix costsa rank is a number of digitscolumns a decade0.7bound, a decade5.7rank at 10⁻⁸6bound at 10⁻⁸49q0.67the shape is rightand the constant is not
Fig. 9 Halfway to the boundary, where q is two thirds and the bound has grown to forty-nine columns for what the matrix does in six.
sep(A, B) and the smallest |λᵢ(A) + μⱼ(B)| as the departure from normality grows, n = 8A and B are upper triangular with diagonals 1, 2, …, 8 and μ above them, so every eigenvalue sum is at least 2 at every point on the axis and the flat line is exact rather than nearly flat. sep — the smallest ‖AX + XB‖_F over ‖X‖_F = 1, and the reciprocal of the amplification a perturbation of C receives — starts equal to the gap at μ = 0, where the matrices are normal, and falls to 4.02·10⁻⁴ at μ = 8. The number a reader is invited to consult is the one that does not move.01234567810⁻⁵10⁻⁴10⁻³10⁻²10⁻¹110¹μ, the entry above the diagonalsep, and the eigenvalue gapmin |λᵢ + μⱼ|sep(A, B)the spectra never moveeigenvalue gap, throughout2sep at μ = 02sep at μ = 84·10⁻⁴amplification there2488solvability is the eigenvaluesand conditioning is not
Fig. 10 A separation against a gap in the spectra field, where the same shape of quantity decides how well an invariant subspace is determined. The two are unrelated in subject and identical in form.

Two decades below the roundoff, and the tail

The last row of the first table is worth reading carefully, because it is the only one that is not on the line.

At 10⁻¹⁴ the count is 9, where the line predicts about 8. The reason is not that the decay changed: it is that the block’s tenth singular value is 2.8·10⁻¹⁶, which is the unit roundoff, and the ones below it are numerical noise rather than singular values of the matrix. A ratio of 10⁻¹⁴ sits close enough to that floor for the count to pick up a column that is not carrying information.

That is a small effect here and it is a large one in a code, because a formatted arithmetic asked for a tolerance below the working precision will happily allocate rank for it. The essay on the site about the zeros nobody meets and everybody writes is about exactly this class of parameter, and the rule it arrives at applies unchanged: a tolerance below the unit roundoff is not a tolerance, it is a request for every column there is.

Where the Arnoldi recurrence stops, against the tolerance it calls a zero, n = 27The spectrum has 3 coarse clusters of width 10^-4 and 9 fine ones inside them, so this one matrix has three near-invariant subspaces and all three are real. The tolerance decides which the recurrence stops at: it stops at 3 at the loose end, at 9 in the middle, and at 27 at the tight end. The residual it leaves runs from 1.51·10⁻⁵ to 5.48·10⁻¹⁷. Neither end is a mistake and neither is a trade: stopping too early throws away an answer that was available, and stopping too late means dividing by rounding error.10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²0612182430tolerance called a zerostep the recurrence stops atsubspace of 3subspace of 9subspace of 27too late: dividing by roundingtoo early: an answer thrown awayone matrix, three answerssteps at the tightest tolerance27at the loosest3residual left, tightest5.5·10⁻¹⁷residual left, loosest1.5·10⁻⁵both directions are failuresand the constant decides which
Fig. 11 The floor a tolerance runs into, from the field that measured it. Below the working precision a threshold stops selecting and starts counting rounding error.
The singular values of one off-diagonal block, for four kernels on the same 128 pointsTwo intervals that do not touch — [0, 1] and [2, 3] — and the 128 × 128 block between them, for four kernels. 1/r falls a factor of 52 a column and is below 10⁻⁸ after 5; log r behaves the same way, for the reason the expansion makes obvious. An independent draw per entry gives 128 singular values above 10⁻⁸ out of 128, on the same size and the same density, which is what makes the other curves a measurement rather than a property of sorted numbers. The block has full algebraic rank in every case; what differs is where the numbers stop mattering, and that is a decision rather than a fact about the matrix.051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, in orderσ ⁄ σ₁eight digitsan independent draw per entry1 ⁄ ra cliff, and a control1/r rank at 10⁻⁸5log r rank at 10⁻⁸5noise rank at 10⁻⁸128σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns
Fig. 12 The spectrum the count is taken from, at the size this essay measures. The last four values on it are the arithmetic rather than the matrix.

What the straight line is not

Two readings this figure invites and neither of them survives contact with the next essay.

It is not a claim that the block is cheap in absolute terms. Five columns out of 128 is cheap; five columns out of 8 would not be. Everything here is a statement about a ratio, and the ratio only becomes a saving when the block is large — which is a statement about what a partition does with the matrix rather than about any single block, and the essay on storage against size is where it is settled.

And it is not a claim about the whole matrix. The line above is one block. A hierarchical representation is a hundred blocks of a dozen different shapes, each with its own q and therefore its own line, and the accuracy asked of it is applied to each block against that block’s own largest singular value. What that does to the accuracy of the assembled matrix is a separate question with a surprising answer — measured, the whole matrix comes out about thirty times more accurate than the per-block tolerance asked for, and the factor belongs to the kernel rather than to the number of blocks.

The general shape of that mistake is one this collection keeps meeting. A quantity controlled locally and reported globally is two quantities, and the ratio between them is nobody’s parameter. the-part-of-a-solver-that-may-be-rounded is the same shape about precision, and a-small-residual-is-not-a-small-error is the same shape about a residual. Here it is about a tolerance, and the honest thing is to measure the ratio rather than assume it is one.

The accuracy asked for, against the accuracy obtained, for two kernels on one partitionε is applied to each block against that block's own largest singular value; the error is then reported against the whole matrix's norm. The two are not the same number and the dashed diagonal is where they would be. For 1/r the obtained error is 29 times smaller than the tolerance asked for; for log r on the identical partition it is 13 times smaller. Both curves are straight and parallel to the diagonal, so the knob does what a knob should — a decade in buys a decade out — and neither of them sits on it. The factor is how much of the matrix's mass lives off the diagonal, which is a property of the kernel; it is not the number of blocks, which is 112 and would put the curves on the other side of the diagonal.10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³ε asked for, per block‖A − A_H‖ ⁄ ‖A‖ obtainedasked = obtainedlog r1 ⁄ rtwo numbers, not one1/r, obtained at 10⁻⁸3.4·10⁻¹⁰log r, obtained at 10⁻⁸7.8·10⁻¹⁰1/r, obtained ⁄ asked0.034log r, obtained ⁄ asked0.078blocks in the partition112the tolerance is per blockand the error is per matrix
Fig. 13 The accuracy asked for against the accuracy obtained, on two kernels sharing one partition. The dashed diagonal is where a per-block tolerance would be a whole-matrix tolerance.
Backward and forward error against the condition numberA log–log plot over twelve decades of condition number. The backward error is a flat line at ten to the minus sixteen; the forward error rises in proportion to the condition number.110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴10⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹condition number κ(A)relative errorforward errorbackward errorpredicted: κ · u8×8, 20 seeds per κ; dashed is the worstthe problem worsens, not the method
Fig. 14 And the collection’s oldest instance of the same mistake, from the error field: a small residual and a small error are two numbers, and their ratio is the condition number rather than one.

The refusal

The claim under test is the one every course states in its second week: that the rank of a matrix is a property of the matrix.

It is true in the algebra and this essay does not dispute it. What is refused is its use — the reading that a matrix has a rank which a program can go and find. The assertion that the same 128 × 128 block needs the same number of columns at 10⁻² and at 10⁻¹⁴ is fed those two counts, 2 and 9, and it fails.

The refusal is deliberately not fed a hard case. It is fed the easiest possible one: one matrix, two readings, twelve decades apart, nothing changed in between. If the count were a property of the matrix the two numbers would be equal, and they are not equal by a factor of four and a half. That is the whole of it, and it is the sentence the rest of this field is written on top of.

Singular values of a rank-4 matrix with noise of relative size 10⁻⁸Ten singular values on a logarithmic axis. The first four sit near one; the rest sit at the noise level, and the vertical distance between the two groups is the evidence for the rank.1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision
Fig. 15 The essay that established the point in the spectra field, where the decision is about whether a matrix is rank-deficient. Here it is about how much of it to keep, which is the same decision put to work rather than reported.
The diagonal of a pivoted R against the singular values, on a matrix of planted rank 4Two curves down the index, lying on top of each other. Both fall by more than five orders of magnitude between index 4 and index 5, so both report rank 4 at any threshold placed in the gap.1611162126313610⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹indexmagnitude|r_kk|σ_kone factorisation, two verdicts‖AP − QR‖/‖A‖1.1·10⁻¹⁵|r_nn|10⁻¹²σ_min10·10⁻¹³column interchanges37|r_nn| is never below σ_minso the cheap verdict errs one way only
Fig. 16 And the case where the decision is hardest to make: a spectrum with no gap in it, where every threshold is as defensible as every other.
The rank of an admissible block and of a touching one, against how finely they are sampledBoth blocks are of the kernel 1/r; both are 32, 64, 128, 256 points a side; both are truncated at 10⁻⁸. The admissible pair — [0, 1] against [2, 3] — needs 5, 5, 5, 5 columns, which is one number. The touching pair — [0, 1] against [1, 2] — needs 9, 11, 12, 13, climbing by about one per doubling, which is a logarithm. Neither of them grows like the block, and only one of them stops. That difference is what the admissibility test in a partition is buying, and it is why the touching pair is kept dense rather than compressed at all.45678903691215log₂ of the points a sidecolumns above 10⁻⁸two intervals that touchtwo that do notthe measurement that is a null resultadmissible, n = 325admissible, n = 2565touching, n = 329touching, n = 25613stored ⁄ dense at largest0.039the rank belongs to the geometryand not to the sampling
Fig. 17 The next question, which is the next essay: whether the five columns grow when the block does.
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. 18 And what the answer to it buys, four essays ahead: numbers stored per unknown, against a dense matrix that stores n of them.

Where the exchange rate is spent

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 60, 94, 131, 170 and 60, 88, 120, 156, adding about 37 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 33 and 30 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.65, 1.48, 1.37, 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 = 5128.7·10⁴weak, n = 5128·10⁴dense, n = 5122.6·10⁵per doubling37‖A − A_H‖ ⁄ ‖A‖7.4·10⁻¹⁴the dense line doublesand the other two add a constant
Fig. 19 The same rate applied to a whole matrix at twelve digits: a straight line against log n, sitting a constant above the eight-digit one.
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 48, 64, 80, 96 and 36, 48, 62, 76, adding about 16 per doubling — a constant per doubling is a logarithm, and it is the whole claim. At n = 512 that is 19 and 15 per cent of the dense matrix, and the share falls at every size. The exponent between consecutive sizes is 1.42, 1.32, 1.26, falling towards one and never arriving, which is what a logarithm looks like from inside.56789100108216324432log₂ nnumbers stored per unknowndenseweak: every off-diagonal blockstrong: only the admissible onesa constant per doublingstrong, n = 5124.9·10⁴weak, n = 5123.9·10⁴dense, n = 5122.6·10⁵per doubling16‖A − A_H‖ ⁄ ‖A‖5.2·10⁻⁶the dense line doublesand the other two add a constant
Fig. 20 And at four, where the whole line drops and its slope does not, which is the separation this essay’s straightness is about.
The largest rank in each partition, and what refusing to compress a touching pair costsThe strong rule's worst block is rank 3 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 5, 6, 7, 7, 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 — 49,152 numbers against 38,912 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.56789100246810log₂ nlargest rank in the partitionweak: touching pairs compressedstrong: touching pairs refusedthe test bounds a rank and costs storagestrong, blocks250weak, blocks94strong, numbers4.9·10⁴weak, numbers3.9·10⁴strong ⁄ weak1.3the better partitionis the more expensive one
Fig. 21 What the ranks in a partition do under the same knob: both sets lift, and only one of them keeps its slope with the size.
Multiplications in a hierarchical solve and in a dense factorisation, and where they crossBoth counted rather than estimated: the recursion carries a counter and reports what it actually did. The dense factorisation is n³/3 and rises at exactly three per doubling. The hierarchical solve rises at 2.18, 1.93, 1.78 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 128 and n = 256: below it the format is the more expensive way to solve the system, at 1.97 times the dense count, and at n = 512 it is 4.4 times cheaper. Every point returns an answer at a backward error of about 1.1·10⁻¹², so the comparison is between two ways of getting the same thing.567891010⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 642ratio at n = 5120.23exponent, first doubling2.2exponent, last doubling1.8backward error1.1·10⁻¹²cheaper is a sizenot a property
Fig. 22 The cost of a solve at ten digits, where the extra columns have walked the crossover to the right.
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. 23 A factorisation’s residual against the same knob, tracking at a slope of one over ten decades — the exchange rate arriving at the object the site’s rule is about.
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. 24 The same trade taken from the other side, from the randomised field: a rank chosen first, and the accuracy it happens to buy.
Columns needed by the block between a separator's two halves, against the length of the separatorThe block is 3×4, 5×6, 7×8, 11×12 across the sweep, so its size more than triples. In the separator's own numbering — the unknowns in the order they sit on the line — it needs 3, 3, 3, 4 columns at 10⁻⁴: about one more per doubling, which is a logarithm. Renumber the same Schur complement by one symmetric permutation and the same block needs 3, 5, 7, 10 — all of them at eight digits and tighter, and a multiple of the geometric count at every tolerance drawn. Nothing about the matrix changed. What changed is that a cluster of a shuffled numbering is scattered along the whole separator, so every block mixes near interactions with far ones and there is no smooth kernel left to compress. The ordering that works here is not chosen: it is the one the separator's geometry already has, which is why nested dissection leaves a compressible fill and a fill-reducing ordering chosen on the graph alone need not.0481216202402468101214unknowns on the separatorcolumns above 10⁻⁴the same matrix, renumberedin the separator's own orderthe ordering the geometry hands overseparator 73separator 234renumbered, largest10the block, largest11share of the square stored0.35the fill is totaland it is not independent
Fig. 25 And in the sparsity field, where the accuracy is a knob on how much of a Schur complement has to be stored.
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. 26 What the knob does not control, from later in this field: the cost of the roundings, which is four per cent at every setting of it.
What a hierarchical preconditioner costs, against how accurate it is, at κ = 1.75e+2The iteration is conjugate gradients on the true matrix and the preconditioner is a solve with the representation, so what is measured is how much of the compression the outer loop repairs and what the repair costs. The count falls monotonically — 16, 14, 7, 4, 3, 2 steps against 72 with no preconditioner — and it keeps falling all the way to the tightest accuracy, which is the number a solver prints and the one everybody watches. The total work, one dense product with A plus one solve with A_H per step, is the other curve, and its minimum is at ε = 0.5 where the preconditioner has rank 1. Past that point every extra digit in the preconditioner is bought and not needed: at the tightest setting the run costs 1.59 times the cheapest, having taken 14 fewer steps to get there.-9-7-5-3-1110⁵10⁶10⁷log₁₀ ε of the preconditionermultiplications for the whole solveno preconditioner: 72 stepsiteration count, rescaledmultiplications spenttwo curves, opposite waysκ175steps, no preconditioner72cheapest ε0.5its rank1tightest ⁄ cheapest1.6the count is what is printedand the cost is what is spent
Fig. 27 And where the knob has an optimum rather than a monotone price, which is the one place on this site it does.

Where this goes next

051015202530354010⁻¹⁷10⁻¹⁴10⁻¹¹10⁻⁸10⁻⁵10⁻²10¹singular value, in orderσ ⁄ σ₁eight digitsan independent draw per entry1 ⁄ ra cliff, and a control1/r rank at 10⁻⁸5log r rank at 10⁻⁸5noise rank at 10⁻⁸96σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns A block nobody can call sparse A 96 × 96 block of a kernel matrix has ninety-six nonzero singular values and five that matter. It has no zero entries, it is not described by fewer numbers than it contains, and neither of the two ways this collection already knows to make a large matrix affordable applies to it. 1234567891010⁻¹⁸10⁻¹⁵10⁻¹²10⁻⁹10⁻⁶10⁻³1indexsingular valuecutoff, σ₁ · 10⁻¹⁰numerical rank 10gap 8.2·10⁶an opiniontrue rank 410×10, built with 4 nonzero valuesrank is a decision Rank is a decision A floating-point matrix does not have a rank. It has a spectrum of singular values, and somewhere in that spectrum is a place where the values stop being signal and start being noise. Deciding where is a judgement, and the evidence for it is a gap. 0816243240110²10⁴10⁶10⁸10¹⁰10¹²10¹⁴matrix size ngrowth factor max|u| / max|a|the 2ⁿ⁻¹ boundworst of 30 randommedian randomWilkinson's matrix sits on the bound30 Gaussian matrices per sizeat n = 40: bound 5.5·10¹¹, worst 4.8 The bound that is never attained Partial pivoting's stability guarantee permits the entries to double at every step — a factor of 5.5·10¹¹ at n = 40. The measured growth on random matrices of that size is about three. The gap is eleven orders of magnitude, and the guarantee is still worth having.

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.

AdmissibilityKernel matrixLow-rank approximationLower boundNumerical rankOff diagonal rankSeparable expansionSpectral decayToleranceTruncated svd