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.

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 ⁄ r, falling to the floora 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ᵢ − yⱼ| .

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 σₖ ⁄ σ₁
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.

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.

Why, and the answer is one line of algebra

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

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

which converges whenever

q = (rx + ry) / d < 1 ,

with rx and ry 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, rx = ry = 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. 2 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 factor of five is one point on an axis

The two slopes above — 0.365 decades a column for the written expansion, 1.71 for the decomposition — divide to 4.7, and the sentence they were used to write said the bound is loose by a factor of five in rank. That sentence has a hidden argument. Both numbers were measured on one pair of intervals, [0, 1] against [2, 3], where q is exactly 1/2, and a factor is only a summary of a relationship if it holds along the axis the relationship lives on.

The axis here is the separation, and it is not an incidental parameter: the whole business of deciding which pairs of clusters are allowed to be compressed is a rule about where a pair sits on it. So the sweep is the thing to do. Hold the block at 128 square and the accuracy at eight digits, slide the second interval from touching to far away, and compare the rank the four-number bound asks for against the rank the block turns out to need.

gap q rank needed rank the bound asks for factor
0.05 0.952 10 440 44.0
0.1 0.909 9 219 24.3
0.25 0.800 7 90 12.9
0.5 0.667 6 49 8.2
1 0.500 5 28 5.6
2 0.333 4 18 4.5
4 0.200 4 12 3.0

The bound remains a bound at every row, which is the half of it that is a theorem and is worth no less for being unsurprising. What it is not is loose by a constant. Across pairs that are all admissible — every q below 1, every expansion convergent — its looseness runs from 3x to 44x, and it runs monotonically. The factor of five quoted above is neither a typical value nor an upper one. It is the value at the single geometry the figures happen to be drawn on.

The mechanism is visible in the two rates rather than in the two ranks. As the intervals close, q rises to 1 and the bound’s decay rate log₁₀(1/q) goes to zero with it — at q = 0.952 the bound is promising a fifth of a decade every ten columns. The block’s own decay does not go to zero at all: it is still 0.88 decades a column there, and the matrix still compresses to ten columns at eight digits. A barely-admissible block is not a barely-compressible one. The bound thinks it is because the quantity the bound is built out of, the ratio of radii to distance, is the wrong quantity for a spectrum near the touching limit.

Those seven rows are seven spectra, and drawn rather than tabulated they say something the rank column cannot. The same 96 × 96 block, the same four kernels, the same eight-digit cut — only the second interval moves.

The singular values of one off-diagonal block, for four kernels on the same 96 pointsTwo intervals that do not touch — [0, 1] and [1.05, 2.05] — and the 96 × 96 block between them, for four kernels. 1/r falls a factor of 7 a column and is below 10⁻⁸ after 10; 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 ⁄ r, falling to the floora cliff, and a control1/r rank at 10⁻⁸10log r rank at 10⁻⁸9noise rank at 10⁻⁸96σ₂ ⁄ σ₁0.19σ₆ ⁄ σ₁8·10⁻⁵the block has full rankand five useful columns
Fig. 3 Barely admissible: the intervals are a twentieth of their own width apart, q = 0.952. The 1/r curve still falls off a cliff, at σ₂/σ₁ = 0.192 — a factor of five a column — and is under eight digits after ten columns.
The singular values of one off-diagonal block, for four kernels on the same 96 pointsTwo intervals that do not touch — [0, 1] and [1.25, 2.25] — and the 96 × 96 block between them, for four kernels. 1/r falls a factor of 15 a column and is below 10⁻⁸ after 7; 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 ⁄ r, falling to the floora cliff, and a control1/r rank at 10⁻⁸7log r rank at 10⁻⁸7noise rank at 10⁻⁸96σ₂ ⁄ σ₁0.085σ₆ ⁄ σ₁1.5·10⁻⁶the block has full rankand five useful columns
Fig. 4 A quarter of a width apart, q = 0.8. The first step is now a factor of twelve, and seven columns buy the same eight digits.

The ranks are 10 and 7 against the sweep’s 10 and 7, measured here at 96 square where the table was taken at 128, and the whole column agrees row for row. What is not in the table is the rate. σ₂/σ₁ runs 0.192, 0.143, 0.0846, 0.0489, 0.024, 0.00981 and 0.0034 across the seven separations — the first column of the cliff steepens from a factor of 5.2 to a factor of 294, which is a factor of 56 in the decay rate against a factor of 2.5 in the rank.

That ratio is the answer to why a barely-admissible block is not a barely-compressible one, and it is arithmetic rather than luck: a rank is a number of digits divided by a rate, and a rate is already a logarithm. Eight digits against the 0.88 decades a column fitted at q = 0.952 is nine or ten columns; the same eight digits against the 2.47 decades of the first step at q = 0.2 is three or four. So a factor of fifty-six in the decay buys a factor of two and a half in the rank — and the bound, which is built out of q directly rather than out of its logarithm, moves by a factor of forty-four over the same range.

The singular values of one off-diagonal block, for four kernels on the same 96 pointsTwo intervals that do not touch — [0, 1] and [3, 4] — and the 96 × 96 block between them, for four kernels. 1/r falls a factor of 127 a column and is below 10⁻⁸ after 4; 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 ⁄ r, falling to the floora cliff, and a control1/r rank at 10⁻⁸4log r rank at 10⁻⁸4noise rank at 10⁻⁸96σ₂ ⁄ σ₁0.0098σ₆ ⁄ σ₁3.2·10⁻¹¹the block has full rankand five useful columns
Fig. 5 Two widths apart, q = 1/3: σ₂/σ₁ = 9.81·10⁻³, and four columns are enough.
The singular values of one off-diagonal block, for four kernels on the same 96 pointsTwo intervals that do not touch — [0, 1] and [5, 6] — and the 96 × 96 block between them, for four kernels. 1/r falls a factor of 368 a column and is below 10⁻⁸ after 4; 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 ⁄ r, falling to the floora cliff, and a control1/r rank at 10⁻⁸4log r rank at 10⁻⁸3noise rank at 10⁻⁸96σ₂ ⁄ σ₁0.0034σ₆ ⁄ σ₁1.6·10⁻¹³the block has full rankand five useful columns
Fig. 6 And four widths, q = 0.2, the far end of the sweep. A factor of 294 a column, three columns of log r, four of 1/r — and the same 96 of 96 for the independent draws.

The control does not move at any of them. The block of independent draws keeps 96 of its 96 singular values above the cut at every separation on the sweep, because there is no geometry in it to separate — which is what makes the four falling curves a measurement of the two intervals rather than of the sorting. A reader who suspects the cliff of being an artefact of the particular pair [0, 1] against [2, 3] can watch the pair move by a factor of eighty in q and watch the control not notice.

That is the direction that matters for a partition. An admissibility rule accepts a pair the moment q drops below its threshold, so the pairs it accepts pile up near the threshold rather than far past it — the comfortably separated blocks are the rare ones. The blocks whose rank a bound is actually asked to predict are the blocks it is worst at predicting, and the essay on what the test costs against what it saves is where that becomes a number rather than an inconvenience.

assertTheBoundsLoosenessIsNotAConstant holds the whole sweep: that the bound never goes under, that its looseness varies by more than a factor of five across the range, that it tightens monotonically with the gap, and that at the tightest pair the rank is under twelve while the prediction is over two hundred. The habit it enforces is the one this collection keeps arriving at from different directions — a quantity measured at one setting and quoted as a property is a quantity nobody has swept. The kernel that has nothing to compress at all is the other end of the same warning, and the next essay’s straight line is only straight because somebody moved the knob.

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 ⁄ r, falling to the floora 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. 7 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.

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.

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 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 essay that numbers the same matrix twice 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.

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.

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?

What links here

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

Reads more easily once this is understood

Essays that name this one as worth reading first.

Shares its objects with

Essays that name at least two of the same things, and that neither author linked.

Named objects

A flat tag is an object no other essay names yet.

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