A bound that holds with probability

Every claim this site has made so far is deterministic. A backward error is below a bound or it is not. ‖QᵀQ − I‖ is 10⁻¹⁵ or it is 1. The Eckart–Young error equals $\sigmaₖ₊₁$, and the best approximation there is checks it as an equality rather than as a bound. Where a claim is about typical behaviour it is asserted across seeds, and the point of acrossSeeds is that it holds at all of them.

The methods in this field do not offer that, and pretending otherwise would be the one genuinely dishonest thing this site could do with them.

A randomised low-rank approximation has an error bound that holds with a probability. Run it enough times and some run, somewhere, exceeds the bound — not because the implementation is wrong, but because that is what the theorem says. So the interesting question is not whether the bound holds. It is how often, by how much, and what the spread looks like, and all three are measurable.

The method, in four lines #

Multiply A by a random matrix Ω with k + p columns to get a sample of its range. Orthonormalise the result to get Q. Project: B = QᵀA. Factorise the small thing that results.

That is the whole construction. Its cost is one pass over A for the multiplication, a QR of something narrow, and an SVD of something small — against a full SVD’s cost of touching the matrix many times.

Two parameters, and both matter more than the description suggests.

p is the oversampling: the few extra columns beyond k. Its job is to make the sample of the range robust, and it is what turns a bound that usually fails into one that usually holds. Five is the standard choice and the improvement from zero to five is enormous.

q is the number of power iterations: instead of AΩ, sample with $(AA^{\mathsf T})^q A\Omega$. Each pass replaces the spectrum by its cube, which pulls the sampled directions towards the leading ones. It is the tool for a spectrum that decays slowly, and it costs a full pass over the data each time — the cost the method exists to avoid, which is why nobody uses more than two.

Two errors, and confusing them is the easy mistake #

This cost the site a failed assertion, and the failure was correct.

The first version of the check asserted that no randomised run beats $\sigmaₖ₊₁$. It failed immediately: a run at k = 10 came back at 0.192 against a floor of 0.197.

Eckart–Young was not violated. The quantity being measured was the wrong one.

‖A − QQᵀA‖ is what the range finder left behind. Q has k + p = 15 columns, so this is the error of a rank-15 approximation, and the Eckart–Young floor for rank 15 is σ₁₆ — comfortably below σ₁₁. Every published bound on the range finder is stated against $\sigmaₖ₊₁$ for exactly this reason: the oversampling is what buys the slack the bound needs, and a bound that could not go below $\sigmaₖ₊₁$ would be a very weak statement.

‖A − Aₖ‖ is the error of the rank-k matrix the method finally returns, and that one cannot beat $\sigmaₖ₊₁$, ever, by any algorithm. It is what a user of the method receives.

Both are now computed and both are asserted, in the directions they actually hold: the rank-k error is never below the floor, at any of twenty seeds, and the projection error goes below it in 3 of 20 runs. Asserting both ways means the two can never be silently swapped again.

What the bound does, measured #

The Halko–Martinsson–Tropp average-case bound is quoted in the form the theorem gives it, because a bound restated in a convenient form is a bound that no longer bounds anything.

Then it is measured, and the result is not the one the check was written to find.

Over forty seeds at p = 5 and forty more at p = 20, the bound was violated zero times. The check that was written first — count the failures — counts none, and that result is reported rather than tuned away, because it is the true shape of the thing: on a matrix with an ordinary spectrum the average-case bound is conservative and not close to being attained.

So the assertions were rewritten to say what is informative.

The bound is loose, by a factor of 5.43 against the median error. Anyone provisioning from it is provisioning for something that will not happen — the same practical failure the conjugate gradient bound produces in the rate the condition number predicts, arriving here for a different reason.

And the error moves with the seed, spanning a factor of 1.60 across the forty. That is the assertion with no analogue anywhere else on this site. Every other quantity measured across seeds is asserted to be the same at all of them; this one is asserted to differ, and a run that came back identical at every seed would fail the build.

What a clean sweep does not prove #

Worth stating plainly, because zero failures in eighty runs invites the wrong conclusion.

A bound that is never violated in eighty runs is not thereby a deterministic bound. It is a probabilistic bound, and the distinction survives however many seeds come back clean.

The next sentence is easier to get wrong, and the obvious one — its failure rate is too small to see at this sample size — names the sample size as the binding constraint. It is not. Fifteen hundred seeds at each oversampling, 120×60, k = 10:

Zero violations in fifteen hundred, at every oversampling — and the bound sits twenty-three to thirty-nine standard deviations above the mean of the observed distribution. That is not a tail a larger sample reaches. The whole distribution sits a factor of three below the bound and its own width is a factor of two, so an excursion to the bound would have to be about one and a half times the entire observed range.

So the conclusion above stands exactly and its diagnosis changes. No violation was observed, and that does not make the bound deterministic — but the reason none was observed is that the bound is enormously loose, not that the sample was small. A reader who takes “too small to see at this sample size” at face value will go looking for more seeds, and more seeds is not the missing ingredient.

The looseness is the useful number in its own right, and it tightens exactly where the method is used: 4.56 times the worst observed error at p = 2, and 2.62 at p = 5. A bound a factor of three above the worst of fifteen hundred runs is a planning figure with a large margin — which is what a user wants from it, and is not the same thing as a description of the error. The distinction is the one a bound that is proved draws between what a bound guarantees and what it predicts, and the bound that is never attained is the collection’s standing case of the gap.

This is measured, not asserted producing an uncomfortable answer rather than a satisfying one. Usually a measured quantity either confirms a claim or contradicts it. Here it does neither: the measurement is consistent with the theorem and does not confirm it, and the honest report of eighty clean runs is “no violation was observed at this sample size” rather than “the bound holds”.

The refusal is written accordingly. An earlier version fed the bound a small oversampling and expected some seed to exceed it; no seed did, in thirty runs. An assertion that depends on catching a rare event is an assertion that will fail on somebody else’s Tuesday, so the refusal was moved to something reliably false — the claim that the method returns the same answer at every seed — which is the property that actually distinguishes this field from every other one on the site — and which randomisation does not create structure turns into the field’s own warning.

Why oversampling works at all #

Five extra columns is a small change and it produces a large one, which is worth an explanation rather than a table entry.

The range finder is trying to capture the space spanned by the top k left singular vectors. A random Ω with exactly k columns gives a sample whose projection onto that space is a random k×k matrix, and a random square matrix is invertible with probability one and badly conditioned with uncomfortable probability. The failure mode is not that the sample misses the space; it is that the sample spans it so unevenly that recovering the space amplifies whatever error is present.

That is the same fact the condition number is an amplifier measures in a different setting: the smallest singular value of a random square matrix is small far more often than intuition suggests, and its reciprocal is what multiplies the error.

Adding p extra columns makes the relevant matrix k×(k+p) rather than square, and a rectangular random matrix is far better conditioned than a square one — the smallest singular value of a k×(k+p) Gaussian is bounded away from zero with a probability that improves rapidly in p. The bound’s leading factor is (1 + √(k/(p−1))), which is where the p − 1 in the denominator comes from and why p = 1 is useless while p = 5 is enough.

So the oversampling is not a safety margin added to a working method. It is the ingredient that makes the method work, and the difference between p = 0 and p = 5 is the difference between an algorithm with no useful guarantee and one with a good one.

Drawn at four oversamplings, the band does what the argument says it should, and the bound does something the argument did not mention.

Read the two numbers together across p = 2, 5, 8 and 20 and they move the same way: the seed-to-seed spread at k = 12 runs 1.65, 1.56, 1.28, 1.10 and the bound’s looseness runs 11.2, 5.93, 5.33, 4.05. Oversampling buys a narrower answer and a tighter guarantee, which is what makes it the cheap knob — the extra columns cost one thin matrix product and the passes over A do not change.

The figure shows the other half of it: raising p from 5 to 20 lowers the median error from 0.157 to 0.084 and narrows the spread across seeds. Both are asserted, the second because a method that got better on average while becoming more erratic would be a worse method, and the median alone would not say so.

Power iteration, and diminishing returns #

The slider is the number of power iterations, and what it shows is a decision every implementation has already made.

On a matrix whose singular values decay slowly — algebraic decay with exponent 0.5, which is about as unhelpful as a spectrum can be while still being approximable — the median rank-k error at k = 10 is:

The first pass buys 0.175. The second buys 0.016. A third would buy less again, and each costs a full pass over the matrix. Two is what almost every implementation defaults to, and the table is why.

The build asserts all of it: each iteration improves the result, none of them passes the optimum, and the returns diminish — the second gain is smaller than the first. That last assertion is the one that makes the table an argument rather than three numbers.

What the table does not show is what the iterations do to the band, and it is a larger effect than what they do to the median.

Across q = 0, 1, 2 and 3 the spread at k = 12 runs 1.56, 1.085, 1.016 and 1.008. So a power iteration is not mainly an accuracy device; it is a determinism device. One pass takes the answer from something that moves by half its own size between seeds to something that moves by eight per cent, and the second takes it to a figure a reader could quote. The section below argues that a randomised computation run once has produced a sample rather than an answer — that is true at q = 0 and very nearly false at q = 2, and the difference is one pass over the matrix.

And the bound goes the other way, which is the part worth keeping. Its looseness against the median at k = 12 runs 5.93, 10.6, 11.0 and 11.0 as q goes 0, 1, 2, 3. The published bound is a statement about the un-iterated range finder, so it does not know that the answer improved; the method gets better, the guarantee does not, and the gap between what is delivered and what is promised doubles. That is the opposite direction from oversampling, which tightened both — and the two knobs being opposite on this one axis is the reason to report the looseness beside the error rather than either alone.

There is a detail in the implementation worth naming because omitting it is a classic error. Between power iterations the sample must be re-orthonormalised. Without it, the columns of Y collapse onto the leading singular vector within about three passes — every column converging to the same direction, exactly as a power method is supposed to — and the sketch loses the rest of the range it was built to capture. The rounding that destroys it is the same rounding two Gram–Schmidts is about, and the fix is the same one.

What is asserted here #

No seed beats $\sigmaₖ₊₁$, at every k drawn and every seed — the Eckart–Young floor, checked as a floor.

The projection error is not floored by $\sigmaₖ₊₁$, since it is a rank-(k+p) quantity, with 3 of 20 runs below it. Asserted so the two errors cannot be confused again.

The band is ordered — worst above median above best — which would catch a plotting error that no amount of correct arithmetic would.

And the seed changes the answer, spanning more than 2%.

The refusals: an error claimed below the Eckart–Young floor must throw, and so must a claim that the method returns one answer. Both do.

What the band is for #

Every figure in this field draws a band between the best and worst seed with the median through it, and that is a deliberate break from the rest of the site.

Elsewhere, a seeded generator produces the figure. The seed exists to make the picture reproducible on every build, and any seed would serve because the claim being drawn holds at all of them — acrossSeeds is used to prove exactly that. A single line is the honest representation because there is nothing else to represent.

Here a single line would be a lie by omission. The method returns a different answer at every seed, by a factor of 1.6 in these runs, and a reader shown only the median has been shown the one number that would be reported identically by a deterministic method. The band is the result; the median is a summary of it.

The practical consequence is a habit worth carrying into any use of these methods. A randomised computation run once has produced a sample, not an answer. If the quantity matters, run it several times and look at the spread — which costs a small multiple of a method chosen for being cheap, and is the only way to know whether the number in hand is typical or is the tail.

Where the guarantee actually helps #

The obvious objection to a probabilistic guarantee is that a deterministic one is available: compute the SVD and take the first k columns, and the error is exactly $\sigmaₖ₊₁$ with certainty.

The answer is entirely about cost, and it is worth being precise since the field is often oversold.

Randomisation buys passes over the data, not accuracy. The error of the best rank-k approximation is $\sigmaₖ₊₁$ whatever produces it, so a randomised method is never more accurate — the figure’s floor is the proof. What it offers is an approximation of comparable quality after one or two passes over A, where a full SVD needs many.

That matters exactly when a pass over A is the expensive thing: a matrix too large for memory, one distributed across machines, one available only as a stream, or one that exists as a function rather than an array — where A is a routine that applies an operator and there are no entries to factorise. In every one of those the deterministic method is not slower, it is unavailable.

And it matters not at all when the matrix fits comfortably in memory, which is the case in which the technique is most often reached for.

Reproducibility, which the seed only half solves #

One consequence of the band deserves its own note, because it is where randomised methods meet software practice badly.

Every figure on this site is byte-identical on every build, and in this field that is achieved by fixing the seed — the generator is written out in lib/random.js precisely so that a figure does not depend on a platform’s random number source. So the pictures are reproducible.

The method is not, and fixing a seed does not make it so. A randomised computation with a fixed seed is reproducible on the same machine with the same library version and the same number of threads. Change any of those and the sequence of random numbers, or the order they are consumed in, can change — and the answer changes with it, by the factor the band shows.

That is a genuinely awkward property for a numerical library to have. A deterministic solver returns the same answer on any conforming implementation, to within rounding, and a difference is evidence of a bug. A randomised one returns a different answer legitimately, so the usual test — run it twice, compare — cannot distinguish a bug from the method working as designed. That is the same difficulty two machines, one certificate meets from the reproducibility side, arriving here from the algorithm’s.

The defence is to test the distribution rather than the value: assert that the error is below a threshold at every seed in a fixed set, that the spread is within a range, that the floor is respected. Which is exactly what the assertions in this field do, and it is why they took three attempts to write.