Numerical Linear Algebra
where the algebra stops being the arithmetic
Randomised, and the guarantee that changes kind

Randomisation does not create structure

On a matrix whose singular values are all equal, a rank-ten randomised approximation has error 1.0 — and so does the optimal deterministic one. Neither achieved anything, and only one of them is usually sold with the implication that it might.
17 min read 8 figures Assertions that rejectMeasured, not assertedWhose fault is itWhere the algebra stops being the arithmetic

Randomised numerical linear algebra is presented, more often than not, as a way of getting the same answers for less. The framing is not wrong and it hides the condition it depends on, and the condition is the whole thing.

Take a 60×60 matrix whose singular values are all exactly 1 — a flat spectrum, no decay at all — and ask for the best rank-ten approximation. The randomised method, with oversampling and two power iterations, returns an error of 1.0.

That looks like a damning result until the deterministic answer is computed. The optimal rank-ten approximation, from an exact SVD, has error 1.0 as well.

Neither method achieved anything, and the failure belongs to the matrix. Eckart–Young says the best possible rank-k error is σk+1\sigma_{k+1}, and here every singular value is 1, so the best possible error is the whole of σ₁. There is no low-rank structure to find and no algorithm can find it.

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. 1 The general shape of the argument. Every randomised run sits above the line marked σk+1\sigma_{k+1}, which is the deterministic optimum. On a matrix with a decaying spectrum that line is low and there is something to gain; on a flat one it sits at σ₁ and the whole plot collapses onto it.

What randomisation actually buys

Stated precisely, because the imprecise version is what causes the trouble.

It buys passes over the data. It does not buy accuracy.

The error of the best rank-k approximation is σk+1\sigma_{k+1} whatever produces it. A randomised method is therefore never more accurate than the deterministic one — the floor in every figure of this field is the proof, asserted at every seed and every k. What it offers is an approximation of comparable quality after one or two passes over A, where computing a full SVD requires many.

That is a statement about cost, and no statement about cost can rescue a matrix with no structure to find. The two claims get conflated because they are usually true together: the matrices people reach for these methods with do have decaying spectra, so the method is both cheap and effective, and the second gets attributed to the first.

The site’s refusal is exactly this. The claim that a rank-ten sketch approximates a matrix with a flat spectrum is fed the measured error and must throw. It does.

The three spectra

The library builds matrices by choosing the singular values first, because the spectrum is the only property that matters to any claim in this field. Three shapes cover the range.

Geometric decay, σj=0.85j\sigma_j = 0.85^j. The case randomisation is good at: the tail falls off fast, a rank-k approximation captures nearly everything, and a sketch with modest oversampling finds it.

Algebraic decay, σj=(j+1)−p\sigma_j = (j+1)^{-p}. The honest middle. At p = 1 the tail is substantial and the approximation is mediocre — good enough to be useful, poor enough that the choice of k matters. At p = 0.5 it is bad enough that power iteration is required to get anywhere, which is what a bound that holds with probability’s slider measures.

Flat, σ_j = 1. The case where nothing works, included precisely so that it is shown rather than omitted.

Choosing the spectrum rather than the entries is what makes these comparisons honest. A figure built from a random matrix and then described as showing that a method works has confounded the method with whatever spectrum that random matrix happened to have.

Why a flat spectrum is not an exotic case

It would be convenient if the flat spectrum were a constructed pathology, and it is not — it is what noise looks like.

A matrix of independent random entries has a spectrum described by the Marchenko–Pastur law, which for a square matrix spreads the singular values over a bounded range with no decay to speak of. So pure noise has an approximately flat spectrum, and any matrix that is mostly noise inherits it.

That is the practical form of the warning. A data matrix with genuine low-rank structure buried under substantial noise has a spectrum consisting of a few large values — the structure — followed by a broad flat shelf, which is the noise. Rank-k approximation with k inside the structure works well. Rank-k approximation with k reaching into the shelf is approximating noise, and every additional direction costs a pass over the data to capture something that means nothing.

Which is exactly rank is a decision’s subject arriving in this field: the cutoff is a judgement about where signal stops, the evidence for it is the gap, and a randomised method asked for a rank past the gap will return one without complaint. It has no way to know it has been asked for noise, and the error it reports will look reasonable — because approximating noise at rank k does have error σk+1\sigma_{k+1}, and σk+1\sigma_{k+1} on the shelf is not a large number.

Sketch-and-solve, and what it costs

The other application worth measuring, because it is the one where a number can be attached to “worse but cheaper”.

An overdetermined least-squares problem with 400 rows and 12 columns: sketch the rows down to 60 and solve the small problem instead. The saving is direct — a QR of a 60×12 system rather than a 400×12 one — and the question is what the answer is worth.

Across twenty seeds, the residual of the sketched solution relative to the exact least-squares residual has a median of 1.10 and a worst case of 1.27. And the minimum across all twenty seeds is 1.00: no sketch ever beats the exact solution.

All three numbers are asserted, and the third is the important one. The exact least-squares solution minimises the residual over all x, so nothing can do better — a sketched solve returning a smaller residual would mean the exact solve was not exact, and the assertion would catch it. It is the same floor as Eckart–Young’s, in the least-squares setting.

Ten per cent worse for a seventh of the work, with a tail reaching 27%. Whether that is a good trade is not a numerical question, and the point of measuring it is that the trade can be stated rather than guessed at.

It is also worth noticing which quantity is 10% worse. The residual is 10% larger; the coefficients can be much further off than that, and by a factor governed by the conditioning of the design matrix. A small residual is not a small error is the essay for that gap, and it applies here unchanged — a sketched solve reporting a residual within 10% has said nothing about how close its coefficients are to the ones the full problem would have given.

There is a refinement worth naming: sketch-and-precondition. Use the sketch not to produce the answer but to build a preconditioner for an iterative solve of the full problem — the sketched R factor makes AR⁻¹ well conditioned, and a few iterations then give the exact answer rather than a 10%-worse one. It is the same move changing the condition number on purpose describes, with the preconditioner built by sampling, and it converts a fast approximate method into a fast exact one.

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. 2 Why the sketched problem resembles the original at all. The projection preserves lengths to within a measured distortion, so the sketched residual is close to the true one — and the 10% figure above is that distortion arriving in the answer.

Where the flat spectrum shows up in this site’s own machinery

The claim that a flat spectrum defeats approximation has a companion elsewhere on the site, and putting them together makes both sharper.

The valley with no bottom measures a fit whose coefficients can move 39% of their own size along one direction before the residual changes in the sixth digit. That is a small singular value: one direction in which the problem determines almost nothing.

A flat spectrum is the opposite extreme — no direction is more determined than any other — and it produces the opposite pathology. There is nothing to discard, so any truncation loses as much as any other, and the matrix cannot be summarised at all.

Between them is where every useful matrix sits, and the quantity that says where is the shape of the spectrum rather than any single number drawn from it. The condition number reports the ratio of the two ends and says nothing about the middle; a matrix with κ = 10⁶ might have one small singular value and a tight cluster, or a smooth decay across six decades, and those are entirely different situations for every method in this field.

That is a limitation of κ worth carrying back to the fields that use it. The rate the condition number predicts finds the conjugate gradient bound loose by an order of magnitude for exactly this reason — κ discards the middle of the spectrum, and the middle is what the method responds to. Here the same discarding makes κ useless as a predictor of approximability: the flat matrix has κ = 1, the best possible value, and is the least approximable matrix in this essay.

The seed spread, which is what a flat spectrum still has

One observation from the flat-spectrum experiment that is worth keeping, because it is a little counterintuitive.

The method returns error 1.0 on the flat matrix at every seed, so on that matrix it is reproducible in the only sense that matters — uselessly, but consistently. On a matrix with algebraic decay, thirty seeds give rank-k errors spanning a factor of 1.79.

That is where the refusal for this essay lives. The claim that the randomised SVD returns the same error at every seed is fed those thirty values and must throw, and it does.

It was not the first refusal written. The first fed the published bound a small oversampling and expected some seed to exceed it; no seed did, in thirty runs, because the bound is conservative enough that a violation is not something a build can rely on catching. An assertion that depends on observing a rare event is an assertion that will fail on somebody else’s Tuesday, so it was replaced by one that is reliably false — and the replacement happens to state the property that distinguishes this field from every other one on the site.

What is asserted here

On a flat spectrum the randomised error equals the optimal one, to within 5%, and the optimal one is the whole of σ₁, to 10⁻⁶.

No sketch beats the exact least-squares residual, across twenty seeds.

A sketch of 60 rows from 400 stays within 30% of it, at the median.

And it has a tail, with the worst seed above the median — asserted, because a random method whose worst case equalled its median would not be a random method.

The refusals: the claim that sketching finds structure in a flat spectrum must throw, and so must the claim that a randomised method returns one answer. Both do.

When to reach for it, and when not

The decision is not subtle once the two claims are kept apart, and it comes down to three questions.

Does the matrix have a decaying spectrum? If not, nothing here helps, and neither does anything else — the answer is that a low-rank approximation is the wrong tool, not that a better algorithm is needed. Rank is a decision is about how to tell, and the answer is the gap in the spectrum, which is a measurement rather than a hope.

Is a pass over the data the expensive thing? These methods pay for themselves when A is too large for memory, distributed across machines, arriving as a stream, or available only as a routine that applies an operator with no entries to factorise. In every one of those the deterministic method is not slower, it is unavailable. When the matrix fits comfortably in memory — which is the situation the technique is most often reached for — an exact SVD is fast, deterministic, and gives a better answer.

Is an approximate answer acceptable? If the result feeds a decision that tolerates a few per cent, yes. If it feeds something where a 27% tail matters, either run it several times and look at the spread, or use sketch-and-precondition and get the exact answer anyway.

The cost model, spelled out once

The claim “randomisation buys passes over the data” deserves the arithmetic behind it, since the whole argument rests on cost and cost is the thing this essay has so far only asserted.

A full SVD of an m×n matrix costs O(mn²) operations, and — more to the point on a large problem — it requires the matrix to be available for repeated access. A bidiagonalisation touches every entry many times, in an order that is not sequential.

The randomised range finder touches A twice: once for AΩ, once for QᵀA. With q power iterations it is 2q + 2 times. Everything else happens on matrices of width k + p, which is small, and costs O(mk²) rather than O(mn²).

So on a matrix of 10⁶ rows and 10⁴ columns where a rank-50 approximation is wanted, the deterministic route is 10¹⁴ operations against the randomised route’s 10⁹ — and the deterministic route additionally needs the whole matrix resident, which at that size it is not.

That is where the technique earns its place, and it is worth noticing that the argument never mentions accuracy at any point. Both routes are approximating the same object with the same Eckart–Young floor. The randomised one gets close to it in two passes; the exact one attains it in many. Whether two passes is the deciding constraint is a property of the deployment rather than of the mathematics, which is why the same method is indispensable in one setting and pointless in another.

The general form

The reframing this essay is for, stated so it survives outside the subject.

A method that is faster on the problems where the answer is easy has not been shown to be faster. Randomised low-rank approximation works on matrices with decaying spectra, and matrices with decaying spectra are the ones where a rank-k approximation is worth computing at all. The technique’s domain and its success are the same set, which makes the success very hard to interrogate from the outside.

The way to interrogate it is the one this essay takes: find the case where the problem has no answer, and check that the method reports that rather than producing one. The flat spectrum is that case. Both methods return 1.0, both are correct to do so, and the only thing separating an honest account of the field from an advertisement is whether that number appears.

One more consequence of that framing, and it is the least comfortable one. Because the technique’s value depends on the deployment rather than on the matrix, the decision cannot be made by looking at the matrix alone — which is where every other decision on this site is made. Whether to pivot, whether to use Householder rather than Gram–Schmidt, whether the normal equations are safe: all of those are settled by properties of A. Whether to sketch is settled by how A is stored, on what hardware, reached how often, and none of that is in the numbers.

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. 3 The deterministic floor, checked as an equality rather than a bound. The best rank-k error is σk+1\sigma_{k+1} — so the shape of the spectrum decides what any method can achieve, before any method is chosen.
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. 4 And the decision that comes first. Whether a matrix has low-rank structure is a judgement backed by a gap in its singular values, and a randomised approximation of a rank chosen without that evidence is an approximation to a number nobody established.
How far the coefficients can move without changing the fit, degree 9Relative increase in the residual against relative change in the coefficients, along the least determined direction. The residual does not move measurably until the coefficients have changed by more than a factor of one.10⁻⁴10⁻³10⁻²10⁻¹110¹10²10³10⁻¹⁶10⁻¹³10⁻¹⁰10⁻⁷10⁻⁴10⁻¹relative change in the coefficients, along the worst directionrelative increase in the residualcoefficients doubled39% change, fit unmoved in the sixth digit308×: the third digit movesκ(A) = 3.6·10⁶. Exact arithmetic would pick one point on this floor. It would not raise it.24 points, degree 9, monomial basisthe data leaves them free
Fig. 5 The opposite pathology to a flat spectrum. Here one singular value is small and the coefficients along its direction are undetermined; there, none is small and nothing can be discarded. Both are matrices that cannot be summarised, for opposite reasons.
Least squares by QR and by the normal equations in binary32Relative error of the computed coefficients against epsilon, on log axes. The QR route is a flat line near the bottom; the normal-equations route climbs and then stops, at the epsilon where the cross-product matrix becomes exactly singular.10⁻⁸10⁻⁷10⁻⁶10⁻⁵10⁻⁴10⁻³10⁻²10⁻¹10⁻⁹10⁻⁷10⁻⁵10⁻³10⁻¹10¹ε in Läuchli's matrix (smaller ε, larger κ)relative error in the coefficientsAᵀA exactly singularnormal equationsQRκ from 1.7·10⁸ to 17the cliff is at √u = 2.4·10⁻⁴
Fig. 6 The exact route a sketched solve is measured against. The least-squares residual it cannot beat is the one this figure’s better road produces, and the 10% is the distance from it.
Fitting the same degree-11 polynomial in two basesOn the left, the data and two fitted curves that lie on top of each other. On the right, the condition numbers of the two design matrices, ten orders of magnitude apart.00.10.20.30.40.50.60.70.80.910.20.40.60.81xyboth fits, drawn on top of one anothercondition number of the design matrixmonomial1.2·10⁸Chebyshev2.5largest fitted coefficientmonomial113Chebyshev0.51rms residual: 3.9·10⁻⁶ and 2.9·10⁻⁷ — the data is fitted either way.30 points, degree 11, single precisionthe basis is part of the problem
Fig. 7 And why a residual within 10% says little about the coefficients. The conditioning of the design matrix decides how far a small change in residual can move the fit, which is a factor the sketch’s guarantee does not mention.
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^−2.50flatmeasured, and equal to σ₁₁rank-10 error, geometric0.2rank-10 error, algebraic0.0025rank-10 error, flat160×60, spectrum chosen rather than the entriesthe matrix decides, not the method
Fig. 8 The three cases on one axis, with the rank-ten error marked on each. The flat curve is this essay’s subject: its optimum is the whole of σ₁, and the measurement is asserted as an equality against Eckart–Young rather than read off the plot. Nothing to the right of that line is a shortcoming of a method.

Where this goes next

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 A bound that holds with probability Every other guarantee on this site is deterministic. The randomised low-rank approximation offers one that holds with a probability, the seed changes the answer, and the honest figure is a band rather than a line. 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 The dimension does not appear A random projection preserves the lengths of a set of vectors to within a distortion that depends on how many vectors there are and not on how many coordinates each one has. That is the fact the whole field rests on, and it is genuinely surprising. 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 σₖ₊₁ The best approximation there is The error of the best rank-k approximation is not bounded by the next singular value. It is equal to it. That is an unusually sharp theorem, and it makes the theorem itself usable as an independent check on the computation.