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

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.
34 min read 26 figures Where the algebra stops being the arithmeticMeasured, not assertedAssertions that reject

Worth reading first: Rank is a decision · The factor is not sparse · An operator with no entries.

This collection has two answers to the sentence the matrix is too large to store, and both of them are things the problem hands over rather than choices anybody makes.

The first is sparsity. Most of the entries are zero, so keep the ones that are not, and the whole sparsity field follows from what that costs at elimination. The second is not forming it at all: an operator can be a function, y = A(x), with no entries anywhere, and the essay on that ends where this one begins — an operator that can only be applied is an operator that cannot be factorised, which is why every solver in it is an iteration.

There is a third, and it is the one this field is about. It applies to matrices that are not sparse and not small, it is the only one of the three whose cost is chosen rather than handed over, and the choosing is the whole subject.

The singular values of one off-diagonal block, for four kernels on the same 96 pointsTwo intervals that do not touch — [0, 1] and [2, 3] — and the 96 × 96 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 96 singular values above 10⁻⁸ out of 96, 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⁻⁸96σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns
Fig. 1 The singular values of one off-diagonal block, for four kernels on the same points. Three of the curves fall off a cliff. The fourth is what makes that a measurement.

The matrix

Take 96 points spread evenly over the interval [0, 1] and another 96 over [2, 3], and build the matrix whose (i, j) entry is

1 / |x_i − y_j| .

Every entry is nonzero. The smallest is about 1/3 and the largest about 1, so there is not a single one that could be dropped and no threshold at which the picture becomes sparse. It is not a Toeplitz matrix either — the points are the same spacing, so in this particular arrangement it happens to be one, but move the points off a uniform grid and it stops being, while everything below survives unchanged. The structure field’s whole subject is a matrix described by fewer numbers than it contains, and this is not that.

So neither existing answer applies. There are 9,216 entries and no reason yet to think fewer than 9,216 numbers are needed.

The measurement

Take its singular values and divide by the largest:

k σ_k ⁄ σ₁
1 1
2 2.40·10⁻²
3 4.59·10⁻⁴
4 8.47·10⁻⁶
5 1.54·10⁻⁷
6 2.79·10⁻⁹
7 5.01·10⁻¹¹
8 8.97·10⁻¹³
9 1.60·10⁻¹⁴
10 2.84·10⁻¹⁶

A factor of about fifty a column, held for nine columns, until it hits the unit roundoff and stops. Fitted over the first five it is 1.71 decades a column.

Five of the ninety-six are above 10⁻⁸ of the first. By the theorem this collection already has an essay about — the best rank-k approximation is the one the singular value decomposition writes down, and its error is the next singular value — a rank-five approximation of this block is wrong in the eighth digit. Two factors of 96 × 5 are 960 numbers, against 9,216 entries.

Error of the best rank-k approximation to a 10×10 matrixApproximation error against k on a logarithmic axis, with the measured error and the next singular value drawn as separate curves lying exactly on top of one another.12345678910⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1rank k of the approximation‖A − Aₖ‖measured, 2-normσₖ₊₁, from theorymeasured, Frobeniusthe first two agreeto 4.3·10⁻⁹worst |‖A−Aₖ‖₂ − σₖ₊₁| / σₖ₊₁4.3·10⁻⁹worst Frobenius discrepancy4.3·10⁻⁹κ = 10⁹; 30 random rank-3 matrices, none closerthe error is σₖ₊₁
Fig. 2 The theorem the measurement above is being read through, from the field that established it: the truncation is optimal, and the error it leaves is the singular value you stopped at.

The control, which is the argument

A curve that falls is not evidence of anything. Every matrix has decaying singular values in the sense that they are sorted, and a reader shown one falling line has been shown a convention.

So the figure carries a second block, of exactly the same size, on exactly the same points, whose entries are independent draws from a normal distribution instead of a function of the distance between two points. Its singular values fall too — from 1 to about 10⁻² over ninety-six columns, which is what a random matrix does — and 96 of its 96 are above 10⁻⁸.

That is the whole comparison. Same size, same density, no zero entries in either, and one of them needs five columns while the other needs all ninety-six. Nothing about the shape of the matrix distinguishes them. What distinguishes them is that one of them is a function of the geometry of two sets of points and the other is not.

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. 3 The decision this measurement rests on, from the field that named it: a numerical rank is a cut across a spectrum, and where you put the cut is a choice rather than a property.

Why, and the answer is one line of algebra

The kernel 1/(x − y) can be expanded about the two intervals’ centres. Write x = c_x + u and y = c_y + v, with d = c_x − c_y the distance between the centres. Then

1/(x − y) = 1/(d + u − v) = (1/d) · Σ_k (−1)^k ((u − v)/d)^k ,

which converges whenever

q = (r_x + r_y) / d < 1 ,

with r_x and r_y the half-widths of the two intervals. Truncate the sum at order p and expand (u − v)^k by the binomial theorem, and every term separates into a function of u times a function of v. That is a rank p + 1 matrix, written down without any linear algebra at all — two tables of powers of a shifted coordinate — and its relative error is bounded by q^(p+1)/(1 − q).

For the two intervals here, r_x = r_y = 1/2 and d = 2, so q = 1/2 and the expansion converges.

Two things follow and they are worth separating, because the rest of this field lives in the gap between them.

  • The block is approximately low rank, and the reason is a property of the geometry rather than of the entries. That is the finding.
  • The rank the expansion asks for is not the rank the matrix needs. The expansion loses a factor of two per column and the matrix loses a factor of fifty. Both are geometric; they are not the same geometric.
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. 4 Three routes to the error of a rank-k approximation of this block: the bound from four numbers, the expansion actually evaluated, and the decomposition. Two of them share no arithmetic.

Two routes, which is why the middle curve is on the figure

The expansion above is a second route to a number the decomposition already gives, and this collection’s standing habit is that a number arrived at by one route has been wrong every time it has been published.

So it is evaluated rather than described. The two tables of powers are built, multiplied together, and subtracted from the block, and what comes out is a relative error of 0.224 at order 0, falling to 8.2·10⁻⁶ at order 12 — a slope of 0.365 decades a column, against the log₁₀(1/q) = 0.301 the bound promises. The expansion behaves as its own analysis says it does, and it never beats the decomposition, which it cannot: one of them is optimal by construction.

The three curves are therefore an ordering and a disagreement. The ordering — bound above expansion above decomposition — holds at every rank measured, and it has to. The disagreement is that the constant is wrong by a factor in the exponent, and that is a measurement rather than an error.

The exact solution of a 13×13 Hilbert system beside the computed oneTwo columns of numbers: the exact answer, which is the integers one to thirteen, and the answer double-precision elimination returns, with the number of correct digits beside each.H13 x = b, b formed exactly so that x = (1, 2, …, 13)12345678910111213exact1.00002.00003.00133.97865.18864.993010.46720.048921.2657-2.575319.21478.906013.5113computed6.6 correct digits4.8 correct digits3.4 correct digits2.3 correct digits1.4 correct digit0.8 correct digitno correct digitsno correct digitsno correct digitsno correct digitsno correct digits0.6 correct digit1.4 correct digitbackward error2.2·10⁻¹⁷κ = 1.7·10¹⁸The algorithm solved a neighbouring problem perfectly. That problem's answer is this one.right-hand side built in BigInt rationalsthe truth is known
Fig. 5 The habit this figure is an instance of, from the field that established it: a quantity computed one way is a quantity nobody has checked.

What the third answer actually is

Sparsity is a property of the entries. An entry is zero or it is not, and no decision anybody makes changes how many are.

Matrix-free is a property of the code. Either a routine that applies the operator exists or the code have the entries; the problem decided which before a numerical analyst arrived.

Numerical rank is neither. It is a property of a matrix and an accuracy together, and the accuracy is a number somebody types. Ask for two digits and this block costs two columns; ask for fourteen and it costs nine. That is the sentence this field is organised around, and the next essay is about the shape of that trade — which turns out to be a straight line, so a digit costs the same handful of columns wherever it is bought.

The singular values of one off-diagonal block, for four kernels on the same 32 pointsTwo intervals that do not touch — [0, 1] and [2, 3] — and the 32 × 32 block between them, for four kernels. 1/r falls a factor of 53 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 32 singular values above 10⁻⁸ out of 32, 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.05101520253010⁻¹⁷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⁻⁸32σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.6·10⁻⁹the block has full rankand five useful columns
Fig. 6 The same four kernels sampled a third as finely. The cliff has not moved, which is the subject of the third essay in this field and is worth noticing here.
The singular values of one off-diagonal block, for four kernels on the same 256 pointsTwo intervals that do not touch — [0, 1] and [2, 3] — and the 256 × 256 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 256 singular values above 10⁻⁸ out of 256, 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⁻⁸256σ₂ ⁄ σ₁0.024σ₆ ⁄ σ₁2.8·10⁻⁹the block has full rankand five useful columns
Fig. 7 And nearly three times more finely. The random block’s curve grows longer, because a random matrix’s rank is its size; the smooth one’s does not.

What it does not say

Three things this measurement is not, since each of them is a reading somebody will take.

It is not a statement about the matrix’s rank. The block has rank 96, exactly, in the algebra. The singular values below the tenth are not zero, they are 10⁻¹⁶ — and this collection has an essay about why calling those zero is a decision with two ways to be wrong. Nothing here changes that. What is being claimed is that a rank-five matrix is within 10⁻⁷ of it, and the two sentences are different.

It is not a statement about conditioning. A block whose singular values span fourteen decades is not a well-conditioned object, and a solve with it would be in trouble. Nothing here proposes to. The blocks a hierarchical format compresses are the off-diagonal ones; the diagonal blocks — where the points of the two clusters are the same points — are kept dense precisely because none of this applies to them.

And it is not free. Computing those singular values costs a decomposition of the block, which is cubic. The next four essays in this field are about what to do with the fact; the one after them is about a case where none of it works; and the one on building the representation from products alone is about the fact that the cheapest available route does not compute the decomposition at all.

What a rank-10 approximation can achieve, by spectrumA semi-logarithmic plot of singular value against index for three spectra — geometric decay, algebraic decay, and flat — with the rank-ten approximation error marked on each.1112131415110⁻³10⁻²10⁻¹1index jσⱼσ11geometric, 0.85ʲalgebraic, j^−1.00flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic0.091rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 8 Spectral decay, from the field that measures what it is worth: a randomised method finds what a decomposition would find only when the spectrum has already fallen away.
A randomised rank-k solve, 4 seeds a rankRelative error against the rank kept, on a logarithmic vertical axis, with a vertical bar at each rank spanning the seeds. The band is 1.68 wide at rank 8, where the method is at its worst, and 1.004 wide at rank 20, where it is at its best. The same computation on the same data returns a different answer each time.0481216202428321rank keptrelative errormedianthe answer movesspread at rank 81.7spread at rank 201best median error0.15widest where the method is worstand the bound does not say so
Fig. 9 And what a probabilistic bound on the same quantity looks like. The measurement above is deterministic; the field it borrows its machinery from is not.

Where else this shape has appeared here

The site has met a matrix with hidden structure three times before and named the structure differently each time.

the-matrix-that-is-one-row is a circulant: n numbers, no zero entries, diagonalised exactly by a transform. a-limit-the-matrix-never-reaches is Toeplitz: 2n − 1 numbers, a symbol, and a spectrum that has a limit it never attains. Both of those are matrices described by fewer numbers than they contain, which is a statement that can be checked by counting.

This is not that, and the difference matters. A kernel matrix is described by exactly n² numbers, each of them the value of a function at a pair of points, and there is no shorter description. What it has instead is that large sub-blocks of it are close to matrices with short descriptions, and close is a word with a parameter in it.

The 32 eigenvalues of a circulant, two waysA circulant matrix of size 32 has its eigenvalues in closed form: they are the discrete Fourier transform of its first column. Plotted against an eigensolver's answer for the same matrix, the two curves lie on top of each other to 4.4·10⁻¹⁵ relative. The eigenvectors are the same for every circulant of this size and are known before any entry is looked at.04812162024283210⁻¹1index kλthe transformthe eigensolverC = F* Λ F is a factorisationworst relative disagreement4.4·10⁻¹⁵‖Cx − b‖/‖b‖ from the transform solve9.6·10⁻¹⁶imaginary part of a real spectrum1.9·10⁻¹⁶n = 32, and the whole matrix is 32 numberseigenvectors known in advance
Fig. 10 A matrix genuinely described by fewer numbers than it contains, from the field about those. Its spectrum is exact and its solve is n log n, and no accuracy parameter appears anywhere.
κ of the ρ = 0.5 Toeplitz family, against the limit it never reachesThe condition number of the n×n section of the Kac–Murdock–Szegő matrix ρ^|i−j| at ρ = 0.5, plotted against the size on logarithmic axes, with Szegő's asymptotic value ((1+ρ)/(1−ρ))² = 9 drawn as a horizontal line. The measured curve climbs towards it from below and reaches 99.9% of it at n = 128.10¹10²10¹size ncondition numberlimit 9measureda limit, as a fraction of itselfreached at n = 1281still to go0.0013κ at n = 8, as a fraction0.83every point is below the line and none of them is on itthe limit is not a value
Fig. 11 And the Toeplitz limit, which is the other structure of that kind. Both of them are properties; this field’s is a decision.

The refusal

The claim under test is the one this essay opens by contradicting: that a matrix with no zero entries has to be stored entry by entry.

It is refuted by the measurement, and it would be refuted by any falling curve, which is why the refusal is fed the other block. The assertion rank ≤ 12 at 10⁻⁸ is handed the 96 × 96 block of independent draws on the same points, and it fails — 96 of 96 survive the cut — which is what it must do. Without that, the essay’s first figure is a picture of a sorted list.

The refusal has a second half worth stating. Compressibility is not a property that a matrix either has or lacks; it is a property a matrix has with respect to an ordering of its indices, and the seventh essay in this field takes the same 256 × 256 matrix, renumbers it with one permutation, and finds that it stores 118,208 numbers where the original stored 27,008. Nothing a norm can see moves. The refusal above rules out one wrong reading; that essay rules out the other.

The randomised SVD against the optimum it cannot beatA semi-logarithmic plot of approximation error against target rank. A shaded band shows the spread across seeds, a solid line the optimal error from the exact singular values, and a dashed line the published probabilistic bound well above both.04812162010⁻¹10⁻⁰.⁵1target rank k‖A − A_k‖₂published boundrandomisedσ_{k+1}, optimalhow far apart the three areworst seed spread1.6bound / median at k = 125.9median / optimum at k = 121.960×60, 6 seeds, oversampling p = 5band is best to worst
Fig. 12 The matrix with nothing to find, from the field that named it. A method that exploits structure returns nothing on a matrix that has none, and that is the honest half of every claim in this field.
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. 13 And the measurement the next two essays are about: the rank of an admissible block against how finely it is sampled, beside a touching pair where the same kernel behaves differently.

What it costs to find out

There is an awkwardness in everything above and it is worth naming before the field is four essays further on, because a reader who notices it and is not told will assume it was missed.

The five columns were found by taking a singular value decomposition of the block. That costs O(n³) operations and reads every one of the 9,216 entries — so the result is a 960-number object and the route to it touched every entry of the thing it was replacing. On the problems this format exists for, where the matrix is a boundary integral operator and n is in the millions, forming n² entries is the one thing nobody can afford, and a method whose first step is to form them has answered no question at all.

There are three ways out and this field takes all of them.

The first is that most blocks are never formed. A partition splits the matrix so that the compressed blocks are large and few, and the dense blocks are small and near the diagonal; the total number of entries ever visited is the thing that has to be counted, not n².

The second is the expansion above. It writes down a rank p + 1 approximation from four numbers, so a code that knows its kernel analytically can build the factors and never look at the block. That it is loose by a factor of five in rank is the price of not looking, and for a great many kernels it is worth paying.

The third is the one that took the longest to find and is the most general: sample the block. Apply the operator to a handful of random vectors supported on the block’s columns, read the block’s rows, and what comes back spans the block’s range. The whole representation of a 512-square can be built from 256 applications of an operator that is never assembled — a measurement that belongs to a later essay, and the reason this one is careful not to claim that a decomposition is the only route.

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. 14 The third way out, from later in this field: applications of an operator against entries of the matrix it stands for, on the same axes.
Sketch distortion against sketch width, for 60 vectorsA log-log plot of the worst relative change in vector length against the number of rows in the sketch, for vectors of two dimensions a factor of four apart. The two curves lie almost on top of one another and both fall steadily.10²10².³10².⁵⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁹⁶10⁻¹10⁻⁰.⁵1rows in the sketchworst relative distortiondimension 64dimension 2565 seeds per point, band is best to worstthe dimension does not appear
Fig. 15 And the machinery it borrows, from the randomised field: what a random projection does to lengths, which is the reason a handful of vectors finds a subspace at all.

The number to keep

Five columns, 96 × 96, eight digits. Nine hundred and sixty numbers instead of nine thousand two hundred and sixteen, for a matrix in which every entry is nonzero and no entry is small.

The rest of this field is what follows from that, and each of the next steps is a question the measurement above leaves open. What does an extra digit cost? Does the five grow when the block does? Which pairs of clusters is it true of, and who decides? And — the question that turns a storage format into a numerical method — if the whole matrix is now an approximation, what exactly has been approximated?

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. 16 The first of those, which is the next essay: how many columns a decade of accuracy costs, and the bound the geometry gives for the same quantity.
A 256 × 256 kernel matrix partitioned by the strong rule, with each compressed block's rank112 blocks: 46 kept dense and 66 stored as two thin factors, whose ranks run from 4 to 5. The rule decides one thing — whether a pair of index clusters may be compressed — and it decides it from four numbers about where those clusters sit, before any entry of the matrix is read. The strong rule refuses any pair whose clusters touch, and subdivides instead, so the diagonal is fringed with small dense blocks and no rank on the picture exceeds 5. The whole thing stores 27,008 numbers against 65,536 entries, and reproduces the matrix to 3.41·10⁻¹⁰.545545545455554545545545545454555555454545545545545455554545545545rows and columns, in the order the points arrivedark: kept dense · light: two thin factors, rank printedthe strong partitionblocks112kept dense46largest rank5numbers stored2.7·10⁴‖A − A_H‖ ⁄ ‖A‖3.4·10⁻¹⁰the picture is decidedbefore a number is read
Fig. 17 And where it is going: the whole matrix cut into blocks by a rule that reads four numbers per pair and no entry at all.

The field from here

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. 18 What a digit costs on a tighter pair of clusters, which is the next essay’s subject at the separation a partition actually spends most of its blocks on.
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 7, 7, 7, 7 columns, which is one number. The touching pair — [0, 1] against [1, 2] — needs 13, 15, 17, 19, 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.456789048121620log₂ of the points a sidecolumns above 10⁻¹²two intervals that touchtwo that do notthe measurement that is a null resultadmissible, n = 327admissible, n = 2567touching, n = 3213touching, n = 25619stored ⁄ dense at largest0.055the rank belongs to the geometryand not to the sampling
Fig. 19 The null result at twelve digits rather than eight, where the rank is larger and still does not depend on how finely the block was sampled.
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 κ = 20 needs 9, 12, 16, 22, 33, 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 — 1.6, 3.2, 6.4, 12.7, 25.5 — 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.-101230714212835log₂ of the segment lengthcolumns above 10⁻⁸cos(κr) ⁄ r at κ = 201 ⁄ r, same pointsthe case with no answerq, held fixed0.51/r at L = 0.561/r at L = 86cos(κr)/r at L = 0.59cos(κr)/r at L = 833the geometry did not moveand the rank did
Fig. 20 The case where all of it fails, four essays ahead: identical geometry, a fixed wavelength, and a rank that grows with the size of the picture.
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. 21 What a whole partition’s ranks look like once the observation is applied to a matrix rather than to a block, under the two rules for deciding which pairs qualify.
What one symmetric permutation does to the storage, on a matrix it does not changeThe same 128 × 128 matrix, its rows and columns renumbered by one permutation and its columns by the same one. The condition number is 20.8802 either way, to eight digits; the Frobenius norm is 2155.267 either way, to twelve. Nothing a norm can see has moved. Clustered, the partition stores 10,112 numbers. Shuffled, the admissibility test finds no admissible pair anywhere — every cluster of a shuffled numbering spans the whole interval, so every q is infinite — and the format degenerates to dense storage exactly. The rule with no test to fail does worse than that: it compresses every off-diagonal block regardless, gets ranks up to 57 out of 64, and stores 28,096 numbers — 1.71 times the matrix it was compressing. A rank-119 factorisation of a 128-column block is a more expensive way to write down the block than the block.numbers stored, 128 × 128clustered, strong10,112clustered, weak8,960the dense matrix16,384shuffled, strong16,384shuffled, weak28,096the same matrix, twiceκ, clustered21κ, shuffled21‖A‖_F, clustered2155‖A‖_F, shuffled2155shuffled weak ⁄ dense1.7the compressibility is in the numberingand the numbering is not in the matrix
Fig. 22 And the qualification every sentence above needs: the same matrix under one permutation, where the five columns become a hundred and nineteen.
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 36 times smaller than the tolerance asked for; for log r on the identical partition it is 30 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 46 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⁻⁸2.8·10⁻¹⁰log r, obtained at 10⁻⁸3.3·10⁻¹⁰1/r, obtained ⁄ asked0.028log r, obtained ⁄ asked0.033blocks in the partition46the tolerance is per blockand the error is per matrix
Fig. 23 The gap between the accuracy asked of a block and the accuracy an assembled matrix has, which is a factor of thirty and belongs to the kernel.
The Schur complement left on a 15-unknown separator, shaded by the size of each entryA 15 × 15 grid, its middle column taken as a separator, and both halves eliminated exactly. What is left on the separator is 100 per cent nonzero — the sparsity field's result, unchanged: eliminating a variable couples everything it touched, and by the end everything is coupled to everything. The shading is what that field does not measure. The entries fall away smoothly from the diagonal, because the Schur complement is a discrete Green's function — the operator mapping data on the separator to response on the separator — and away from the diagonal that is an integral operator with a smooth kernel. The outlined block is the 7 × 8 between the separator's two halves: every entry nonzero, and 5 columns describe it to eight digits. The rank structure was not put there by the elimination. It was in the differential operator before anything was discretised.the separator's 15 unknowns, in the order they sit on the lineoutlined: the block between the two halvesdense, and not independententries nonzero1the block56columns it needs5numbers stored150entries in the square225every entry is nonzeroand six columns describe them
Fig. 24 The same cliff arrived at from the sparsity field, where the matrix is what an elimination left behind rather than what a kernel was evaluated at.
Multiplications in a hierarchical solve and in a dense factorisation, and where they crossBoth counted rather than estimated: the recursion carries a counter and reports what it actually did. The dense factorisation is n³/3 and rises at exactly three per doubling. The hierarchical solve rises at 2.13, 1.88, 1.73 — falling, because the cost is n log²n and its exponent is on its way to one. The two cross between n = 64 and n = 64: below it the format is the more expensive way to solve the system, at 0.92 times the dense count, and at n = 512 it is 10.3 times cheaper. Every point returns an answer at a backward error of about 2.3·10⁻⁸, so the comparison is between two ways of getting the same thing.567891010⁴10⁵10⁶10⁷10⁸10⁹log₂ nmultiplicationsdense factorisation, n³⁄3the recursion, countedwhere the format starts payingratio at n = 640.92ratio at n = 5120.097exponent, first doubling2.1exponent, last doubling1.7backward error2.3·10⁻⁸cheaper is a sizenot a property
Fig. 25 What the observation eventually buys: a solve whose exponent falls towards one, against a dense factorisation whose exponent is three at every size.
‖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. 26 And the object the site’s own rule applies to — a factorisation performed inside the format, with the residual that says how far its factors are from the matrix.

Where this goes next

the matrix, lower triangle408 entriesits Cholesky factor1739 entries · 1331 created→‖A − LLᵀ‖/‖A‖1.4·10⁻¹⁶fill, symbolic1331fill, numeric1331n = 144 · density 3.2% · bandwidth 12same matrix, renumberedthe answer is identical to rounding The factor is not sparse A sparse matrix has a factor that is not sparse, and the gap between them is the entire reason iterative methods exist. The entries elimination creates can be counted before any arithmetic runs, from the graph alone. 10⁻¹⁷10⁻¹⁵10⁻¹³10⁻¹¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10⁻¹⁴10⁻¹²10⁻¹⁰10⁻⁸10⁻⁶10⁻⁴10⁻²1εrelative error in J(x)vforwardcentralcancellationtruncationagainst a derivative that is exactforward floor1.3·10⁻¹⁰central floor1.1·10⁻¹²truncation slope, forward1truncation slope, central2no ε reaches the roundoffand the analytic derivative is free of the choice An operator with no entries At the sizes where linear algebra is expensive the matrix does not exist. What exists is a subroutine that returns Av. Every Krylov method survives that unchanged; every algorithm that reads an entry disappears. And the derivative such a code computes is accurate to ten digits instead of sixteen, which turns out to cost nothing at all. 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.

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.

Eckart–YoungKernel matrixLow-rank approximationMatrix-freeNumerical rankOff diagonal rankSingular valuesSparsitySpectral decayTruncated svd