Randomized SVD and sketching

Intuition Beginner

The singular value decomposition reads every column and every row of a matrix before it reports the matrix's shape. For a matrix with millions of rows this is expensive: you pay full price even when the data is nearly low-rank and only a handful of directions carry the information. Randomised methods bet that a few random combinations of the columns are enough to expose those directions, the way a small sample of voters can expose an election's trend.

The bet works because random combinations wash out the small directions and keep the large ones. Multiply your matrix $A$ by a tall thin table $\Omega$ (Greek capital omega) of random numbers and you get a compressed sketch $Y=A\Omega$ with far fewer columns. The columns of $Y$ still span the dominant part of $A$, because the big directions appear in almost every random combination while the small ones tend to cancel. Orthonormalise $Y$ and you hold a cheap proxy $Q$ for the range of $A$.

Sketching plays the same trick on the least-squares fit. Instead of solving one huge over-determined problem, you compress the rows with a single random map and solve a much smaller fit. The compressed fit costs a fraction of the work and, with the right random map, its answer scores almost as well on the original objective as the exact answer does. Speed and accuracy trade off through one knob: how many random rows you draw.

Visual Beginner

The table contrasts the exact and randomised routes on the question that separates them: how much of the matrix you must touch.

Step	Exact SVD	Randomised SVD
Reads from $A$	all rows and columns	$A\Omega$, a few random combinations
Cost driver	full pass over every entry	a handful of matrix-multiplies
Output	exact dominant directions	approximate dominant directions
Error knob	none (exact)	oversampling $p$ and power iterations $q$

The difference is the rightmost column. The exact method cannot return a partial answer; the randomised method lets you spend more random samples to shrink the approximation, trading a controllable probability of error for a large cut in work.

The geometry is the whole lesson. A random multiplier collapses dimension while leaving the large directions intact; the small directions, which the exact decomposition also has to read, are averaged away.

Worked example Beginner

Take a rank-one matrix whose two columns are parallel:

$$A=\left(\begin{array}{cc}3& 6\\ 1& 2\end{array}\right),\text{column two}=2\times \text{column one}.$$

Draw a single random test vector, say $\Omega =\left(\genfrac{}{}{0ex}{}{1}{1}\right)$, and form the sketch:

$$Y=A\Omega =\left(\genfrac{}{}{0ex}{}{3\cdot 1+6\cdot 1}{1\cdot 1+2\cdot 1}\right)=\left(\genfrac{}{}{0ex}{}{9}{3}\right).$$

Normalise $Y$ to unit length. Its length is $\sqrt{81+9}=\sqrt{90}=3\sqrt{10}$, so the proxy basis vector is $Q=\left(\genfrac{}{}{0ex}{}{3/\sqrt{10}}{1/\sqrt{10}}\right)$.

Now project $A$ onto this proxy: $B=Q^{\mathsf{T}}A=\frac{1}{\sqrt{10}}(3,1)A=\frac{1}{\sqrt{10}}(10,20)=(\sqrt{10},2\sqrt{10})$. The single stretch factor of the small matrix $B$ is its length $\sqrt{10+40}=\sqrt{50}=5\sqrt{2}$, and the recovered factorisation is $A\approx U\Sigma V^{\mathsf{T}}$ with $\Sigma =5\sqrt{2}$, $U=Q$, and $V^{\mathsf{T}}=(1,2)/\sqrt{5}$.

Checking, $U\Sigma V^{\mathsf{T}}=\frac{5\sqrt{2}}{\sqrt{10}\sqrt{5}}\left(\genfrac{}{}{0ex}{}{3}{1}\right)(1,2)=\left(\genfrac{}{}{0ex}{}{3}{1}\right)(1,2)=A$. One random sample recovered the matrix exactly, because the data had only one true direction. What this tells us: when the matrix is low-rank, a few random samples suffice; oversampling handles the noise and the rank you did not assume in advance.

Check your understanding Beginner

Formal definition Intermediate+

Let $A\in {\mathbb{R}}^{m\times n}$ have singular value decomposition $A=U\Sigma V^{\mathsf{T}}$ with singular values ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge 0$, and let $A_k={\sum }_{j=1}^k{\sigma }_j{u}_j{v}_j^{\mathsf{T}}$ be the optimal rank-$k$ truncation established in 01.01.12. The deterministic SVD computes the full decomposition at cost $O(mn\min (m,n))$; the randomised methods below trade a controllable probability of error for reading $A$ through only a small number of matrix-vector products.

Random test matrices. A Gaussian test matrix is an $\Omega \in {\mathbb{R}}^{n\times \ell }$ whose entries are independent $N(0,1)$ variables. Structured alternatives — the subsampled randomised Fourier / Hadamard transform (SRHT) and the sparse sign matrix — apply in $O(mn\log \ell )$ or $O(\mathrm{nnz}(A))$ time while retaining the same accuracy guarantees up to constants ^{[Halko, N., Martinsson, P.-G. & Tropp, J. A. — Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions]}.

The randomised range finder. Given a target rank $k$ and an oversampling parameter $p\ge 1$, set $\ell =k+p$. Draw $\Omega \in {\mathbb{R}}^{n\times \ell }$, form the sketch $Y=A\Omega \in {\mathbb{R}}^{m\times \ell }$, and compute a thin QR factorisation $Y=QR$ with $Q\in {\mathbb{R}}^{m\times \ell }$ orthonormal. The columns of $Q$ are the randomised range finder: they span an $\ell$-dimensional subspace that, with high probability, contains the dominant left singular directions $u_1,\dots ,u_k$. The approximation error is $\Vert (I-Q{Q}^{\mathsf{T}})A\Vert$, measured in the spectral or Frobenius norm.

The fixed-rank randomised SVD. Combining the range finder with one small SVD yields the decomposition:

1. Draw $\Omega \in {\mathbb{R}}^{n\times (k+p)}$; form $Y=A\Omega$; orthonormalise $Y=QR$.
2. Project to a small matrix: $B=Q^{\mathsf{T}}A\in {\mathbb{R}}^{\ell \times n}$.
3. Compute the SVD of the small matrix $B=\tilde{U}\tilde{\Sigma }{\tilde{V}}^{\mathsf{T}}$.
4. Lift back: $U=Q\tilde{U}$, and return $(U,\tilde{\Sigma },\tilde{V})$ truncated to rank $k$.

The expensive step touches $A$ only $\ell =k+p$ times (to form $A\Omega$ and $Q^{\mathsf{T}}A$); the SVD runs on the $\ell \times n$ matrix $B$. The dominant cost is $O(mn\ell )$ rather than $O(m{n}^2)$, the savings over 43.05.08's full SVD route when $k\ll \min (m,n)$.

Subspace embeddings and sketch-and-solve. A matrix $\Pi \in {\mathbb{R}}^{d\times m}$ is an $\epsilon$-subspace embedding for a subspace $S\subseteq {\mathbb{R}}^m$ when $(1-\epsilon )\Vert x{\Vert }_2^2\le \Vert \Pi x{\Vert }_2^2\le (1+\epsilon )\Vert x{\Vert }_2^2$ for every $x\in S$. Given the overdetermined least-squares problem ${\min }_x\Vert Ax-b{\Vert }_2$ of 43.04.01, sketch-and-solve draws $\Pi$ that embeds the column space of $[A\mid b]$ and solves the compressed problem ${\min }_x\Vert \Pi (Ax-b){\Vert }_2$ instead. A Gaussian $\Pi$ with $d=O(n/{\epsilon }^2)$ rows is an $\epsilon$-embedding of any fixed $n$-dimensional subspace; a sparse count-sketch $\Pi$ achieves the same in $O(\mathrm{nnz}(A))$ application time ^{[Woodruff, D. P. — Sketching as a Tool for Numerical Linear Algebra]}.

The interpolative decomposition. A rank-$k$ interpolative decomposition (ID) writes $A\approx CZ$ where $C=A[:,J]$ consists of $k$ actual columns of $A$ (indexed by $J$) and $Z\in {\mathbb{R}}^{k\times n}$ contains a $k\times k$ identity on those columns with entries bounded by $2$ in magnitude. The randomised ID first builds the range finder $Q$, then runs a column-pivoted QR on $Q^{\mathsf{T}}$ to select the columns $J$, producing an interpretable factorisation (real columns of $A$) at the cost of the randomised SVD.

Counterexamples to common slips

• The randomised SVD does not return the exact singular values; it returns estimates whose accuracy is governed by $\sqrt{1+k/(p-1)}\Vert A-A_k{\Vert }_F$ in expectation. Demanding exactness defeats the point: set $p$ to meet a tolerance, do not expect the deterministic answer.

• Sketch-and-solve preserves the least-squares objective, not the minimiser. The sketch-optimal $\hat{x}$ satisfies $\Vert A\hat{x}-b{\Vert }_2\le \frac{1+\epsilon }{1-\epsilon }\Vert A{x}_{\star }-b{\Vert }_2$, but $\Vert \hat{x}-x_{\star }{\Vert }_2$ may be large when $A$ is ill-conditioned, because $\Pi$ preserves the residual geometry, not the solution geometry.

• Oversampling $p$ is not the target rank. The target rank is $k$; $p$ is the slack (typically $5$ to $10$) that drives the failure probability down. Setting $p=0$ discards the randomisation guarantee and leaves a brittle estimate.

• The power method is not free. Forming $(A{A}^{\mathsf{T}}{)}^{q}A\Omega$ directly loses precision when the spectrum spans many orders of magnitude, because the powered singular values overflow and underflow; the stabilised implementation re-orthonormalises between each application of $A$ and $A^{\mathsf{T}}$.

Key theorem with proof Intermediate+

Theorem (Johnson-Lindenstrauss; Gaussian embedding). Let $x_1,\dots ,x_N\in {\mathbb{R}}^n$ and $\epsilon \in (0,1)$. Let $\Pi =\frac{1}{\sqrt{d}}G\in {\mathbb{R}}^{d\times n}$ where the entries of $G$ are independent $N(0,1)$ variables. If

$$d\ge \frac{3\ln N}{{\epsilon }^2/4-{\epsilon }^3/6},$$

then with probability at least $1-1/N$ the random map $\Pi$ satisfies

$$(1-\epsilon )\Vert x_i-x_j{\Vert }_2^2\le \Vert \Pi x_i-\Pi x_j{\Vert }_2^2\le (1+\epsilon )\Vert x_i-x_j{\Vert }_2^2$$

for every pair $i,j$ simultaneously. ^{[Johnson, W. B. & Lindenstrauss, J. — Extensions of Lipschitz mappings into a Hilbert space]}

Proof. Fix a pair $i,j$ and write $z=x_i-x_j\in {\mathbb{R}}^n$. The coordinates of $\Pi z=\frac{1}{\sqrt{d}}Gz$ are independent Gaussian variables with mean zero and variance $\Vert z{\Vert }_2^2/d$, since each row of $G$ dotted with $z$ is $N(0,\Vert z{\Vert }_2^2)$. Therefore $\Vert \Pi z{\Vert }_2^2=\frac{\Vert z{\Vert }_2^2}{d}{\sum }_{r=1}^d{g}_r^2$ where $g_1,\dots ,g_d$ are independent $N(0,1)$, and the ratio $\Vert \Pi z{\Vert }_2^2/\Vert z{\Vert }_2^2=\frac{1}{d}{\chi }_d^2$ is a chi-squared variable with $d$ degrees of freedom, scaled by $1/d$.

The standard Chernoff bound on the chi-squared distribution gives

$$\mathrm{Pr}\left[\frac{1}{d}{\chi }_d^2>1+\epsilon \right]\le \exp (-d({\epsilon }^2/4-{\epsilon }^3/6)),\mathrm{Pr}\left[\frac{1}{d}{\chi }_d^2<1-\epsilon \right]\le \exp (-d({\epsilon }^2/4-{\epsilon }^3/6)),$$

using the moment generating function $E[e^{s{g}^2}]=(1-2s{)}^{-1/2}$ for $s<1/2$ and the inequality $\epsilon -\log (1+\epsilon )\ge {\epsilon }^2/2-{\epsilon }^3/3$ together with its lower-tail analogue. By the union bound over both tails, the probability that the single pair $(i,j)$ is not preserved within the factor $(1\pm \epsilon )$ is at most $2\exp (-d({\epsilon }^2/4-{\epsilon }^3/6))$.

There are $\left(\genfrac{}{}{0ex}{}{N}{2}\right)\le N^2/2<N^2$ pairs. A second union bound over all pairs bounds the probability that any pair fails:

$$\mathrm{Pr}[\exists i,j:\text{pair }(i,j)\text{ fails}]\le N^2\cdot 2\exp (-d({\epsilon }^2/4-{\epsilon }^3/6)).$$

Choosing $d\ge \frac{3\ln N}{{\epsilon }^2/4-{\epsilon }^3/6}$ makes the exponent satisfy $d({\epsilon }^2/4-{\epsilon }^3/6)\ge 3\ln N$, so the right-hand side is at most $2N^2\cdot N^{-3}=2/N$, and a sharper constant tightens this to $1/N$. With the stated complement probability, every pair is preserved simultaneously. $\square$

Bridge. This theorem builds toward the subspace-embedding view of sketching, and it appears again in 43.04.01, where sketch-and-solve compresses the overdetermined least-squares fit through a random map that preserves every vector norm in the row space of $A$. The foundational reason the random projection works is the concentration of the Gaussian quadratic form $\sum g_r^2$, and this is exactly the fact that a chi-squared variable with $d$ degrees of freedom sits within $(1\pm \epsilon )$ of its mean with probability $1-\exp (-\Theta (d{\epsilon }^2))$. Putting these together, the same dimension-reduction mechanism that preserves pairwise point distances also preserves the geometry of an entire subspace, and the bridge is that the randomised range finder of 01.01.12's SVD and the sketch-and-solve map are two instances of one principle: a tall thin random multiplier that collapses dimension while leaving the relevant norms intact.

Exercises Intermediate+

Advanced results Master

Theorem (HMT range-finder expected error). Let $A\in {\mathbb{R}}^{m\times n}$ have singular values ${\sigma }_1\ge {\sigma }_2\ge \cdots$, let $k\ge 1$ be a target rank, and let $p\ge 2$ be Gaussian oversampling. Draw a Gaussian $\Omega \in {\mathbb{R}}^{n\times (k+p)}$, form $Y=A\Omega$, and let $Q$ orthonormalise $Y$. Then the expected Frobenius error of the range finder satisfies

$$(\mathbb{E}\Vert A-Q{Q}^{\mathsf{T}}A{\Vert }_F^2{)}^{1/2}\le \sqrt{1+\frac{k}{p-1}}\cdot \Vert A-A_k{\Vert }_F.$$

The same bound holds for the output of the fixed-rank randomised SVD, since the dominant error term is the range-finder projection. ^{[Halko, N., Martinsson, P.-G. & Tropp, J. A. — Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions]}

The bound reduces to the deterministic Eckart-Young tail $\Vert A-A_k{\Vert }_F$ as $p\to \infty$, recovering exactness in the limit of infinite oversampling; for the modest choice $p=5$ it inflates the optimal tail by at most $\sqrt{1+k/4}$. A matching high-probability tail bound, also due to Halko-Martinsson-Tropp, shows the error exceeds $(1+t)\sqrt{1+k/(p-1)}\Vert A-A_k{\Vert }_F$ with probability at most $\exp (-\Theta (p{t}^2))$, so the failure probability falls exponentially with the oversampling budget. The dominant computational cost $O(mn(k+p))$ is independent of $p$ in its leading term whenever $k$ dominates, which is why $p$ between $5$ and $10$ is the practical sweet spot.

Theorem (power method / subspace iteration accelerates slow decay). Let the range finder be applied not to $A$ but to the powered matrix ${\tilde{A}}_q=(A{A}^{\mathsf{T}}{)}^{q}A$, with the same Gaussian $\Omega$ and oversampling $p\ge 2$. By the identity ${\sigma }_j({\tilde{A}}_q)={\sigma }_j(A{)}^{2q+1}$, the HMT expected-error bound applied to ${\tilde{A}}_q$ controls the error on the powered matrix by

$$(\mathbb{E}\Vert (I-Q_q{Q}_q^{\mathsf{T}}){\tilde{A}}_q{\Vert }_F^2{)}^{1/2}\le \sqrt{1+\frac{k}{p-1}}\cdot (\sum _{j>k}{\sigma }_j^{2(2q+1)}{)}^{1/2},$$

so the only singular values that survive in the tail are those already smaller than $1$ raised to the power $2q+1$; for any gap ${\sigma }_{k+1}<{\sigma }_k$ the ratio $({\sigma }_{k+1}/{\sigma }_k)$ enters the reconstruction error raised to the $(2q+1)$-st power, so a modest $q=1$ or $q=2$ converts a barely-decaying spectrum into a sharply-decaying one.

The stabilised implementation avoids forming $(A{A}^{\mathsf{T}}{)}^{q}A\Omega$ directly, which would overflow or underflow once the powered singular values span many orders of magnitude. Instead it alternates an application of $A$ or $A^{\mathsf{T}}$ with a re-orthonormalising QR, so each iteration operates on an orthonormal basis and the dynamic range stays bounded; the cost grows to $(2q+1)$ passes over $A$, still linear in the sample size ^{[Martinsson, P.-G. & Tropp, J. A. — Randomized numerical linear algebra: Foundations and algorithms]}.

Theorem (subspace embeddings from count-sketch; input-sparsity-time least squares). Let $S\subseteq {\mathbb{R}}^m$ be a fixed $n$-dimensional subspace and $\epsilon ,\delta \in (0,1)$. There is a distribution on sparse matrices $\Pi \in {\mathbb{R}}^{d\times m}$, each column of $\Pi$ having a single nonzero entry, with $d=O(n^2/{\epsilon }^2)$ rows, such that $\Pi$ is an $\epsilon$-subspace embedding for $S$ with probability at least $1-\delta$. Moreover $\Pi A$ is computable in $O(\mathrm{nnz}(A))$ time. ^{[Clarkson, K. L. & Woodruff, D. P. — Low-rank approximation and regression in input sparsity time]}

The count-sketch construction hashes each row of $A$ to one of $d$ buckets with a random sign; the product $\Pi A$ sums signed rows into buckets, so its cost is proportional to the number of nonzeros of $A$, not to $m$. Combined with the sketch-and-solve guarantee of Exercise 7, this yields an $O(\mathrm{nnz}(A))+O(n^3/{\epsilon }^2)$ algorithm for $(1+\epsilon )$-approximate least squares — a qualitative improvement over the $O(m{n}^2)$ cost of forming the normal equations when $m\gg n$ and $A$ is sparse. The Gaussian embedding achieves the same guarantee with $d=O(n/{\epsilon }^2)$ rows but costs $O(mnd)$ to apply; the fast Johnson-Lindenstrauss transform of Ailon and Chazelle interpolates, applying in $O(mn\log n)$ while retaining $d=O(n/{\epsilon }^2)$ ^{[Ailon, N. & Chazelle, B. — The fast Johnson-Lindenstrauss transform and algorithmic applications]}.

Synthesis. Randomised numerical linear algebra trades a vanishing probability of error for a large cut in computational complexity, and the foundational reason it works is that a random projection preserves the geometry of a fixed low-dimensional subspace with exponentially high probability — the same concentration fact that powers the Johnson-Lindenstrauss lemma proved above. This is exactly what lets the randomised range finder of 01.01.12 recover the dominant singular subspace from a handful of random matrix-vector products, and the bound $\sqrt{1+k/(p-1)}\Vert A-A_k{\Vert }_F$ generalises the deterministic Eckart-Young guarantee to a probabilistic one whose failure probability the power method and oversampling drive down. Putting these together, the randomised SVD, the sketch-and-solve least squares of 43.04.01, and the interpolative decomposition are one construction read three ways: sketch the range, solve in the sketch, then read the factorisation off the small compressed matrix. The bridge is that the deterministic low-rank approximation of 43.05.08 and 01.01.12 and its randomised counterpart are the same Eckart-Young optimisation, with randomisation replacing the exhaustive column scan by a random probe whose accuracy the lemma controls; the method appears again in 43.07.02, where block Krylov and randomised range finding merge into the randomised Arnoldi and Lanczos iterations used for very large eigenvalue problems.

Full proof set Master

Proposition (power iteration sharpens the spectrum). For $A=U\Sigma V^{\mathsf{T}}$ and integer $q\ge 0$, the powered matrix $(A{A}^{\mathsf{T}}{)}^{q}A$ has the singular value decomposition $U{\Sigma }^{2q+1}{V}^{\mathsf{T}}$, and in particular its $j$-th singular value is ${\sigma }_j(A{)}^{2q+1}$.

Proof. Compute $A^{\mathsf{T}}=V{\Sigma }^{\mathsf{T}}{U}^{\mathsf{T}}$, so $A{A}^{\mathsf{T}}=U\Sigma V^{\mathsf{T}}V{\Sigma }^{\mathsf{T}}{U}^{\mathsf{T}}=U\Sigma {\Sigma }^{\mathsf{T}}{U}^{\mathsf{T}}=U{\Sigma }^2{U}^{\mathsf{T}}$, using $V^{\mathsf{T}}V=I$. Inductively $(A{A}^{\mathsf{T}}{)}^q=U{\Sigma }^{2q}{U}^{\mathsf{T}}$, since orthogonal $U$ satisfies $U^{\mathsf{T}}U=I$ and powers of the diagonal matrix ${\Sigma }^2$ are taken entrywise. Then $(A{A}^{\mathsf{T}}{)}^{q}A=U{\Sigma }^{2q}{U}^{\mathsf{T}}\cdot U\Sigma V^{\mathsf{T}}=U{\Sigma }^{2q}\Sigma V^{\mathsf{T}}=U{\Sigma }^{2q+1}{V}^{\mathsf{T}}$. The diagonal of ${\Sigma }^{2q+1}$ is ${\sigma }_j^{2q+1}$, and $U,V$ are the singular vectors of $A$, so this is an SVD with the claimed singular values. The range finder applied to $(A{A}^{\mathsf{T}}{)}^{q}A$ therefore confronts the amplified spectrum ${\sigma }_1^{2q+1}\ge \cdots \ge {\sigma }_n^{2q+1}$, in which any ratio ${\sigma }_{k+1}/{\sigma }_k<1$ is raised to the $(2q+1)$-st power; this is the rigorous content of the power-method acceleration. $\square$

Proposition (sketch-and-solve approximation guarantee). Suppose $\Pi$ is an $\epsilon$-subspace embedding for the column space of $[A\mid b]$. Let $\hat{x}$ minimise $\Vert \Pi Ax-\Pi b{\Vert }_2$ and $x_{\star }$ minimise $\Vert Ax-b{\Vert }_2$. Then $\Vert A\hat{x}-b{\Vert }_2\le \sqrt{\frac{1+\epsilon }{1-\epsilon }}\Vert A{x}_{\star }-b{\Vert }_2$.

Proof. Write $r(x)=Ax-b$. For every $x\in {\mathbb{R}}^n$, the residual $r(x)=Ax-b$ lies in $\mathrm{col}([A\mid b])$, the embedded subspace, so the two-sided embedding bound applies: $(1-\epsilon )\Vert r(x){\Vert }_2^2\le \Vert \Pi r(x){\Vert }_2^2\le (1+\epsilon )\Vert r(x){\Vert }_2^2$. Apply the upper bound to $x_{\star }$: $\Vert \Pi r(x_{\star }){\Vert }_2^2\le (1+\epsilon )\Vert r(x_{\star }){\Vert }_2^2$. Because $\hat{x}$ minimises the sketched residual, $\Vert \Pi r(\hat{x}){\Vert }_2^2\le \Vert \Pi r(x_{\star }){\Vert }_2^2$, hence by transitivity $\Vert \Pi r(\hat{x}){\Vert }_2^2\le (1+\epsilon )\Vert r(x_{\star }){\Vert }_2^2$.

Apply the lower bound to $\hat{x}$: $\Vert \Pi r(\hat{x}){\Vert }_2^2\ge (1-\epsilon )\Vert r(\hat{x}){\Vert }_2^2$. Chaining the two inequalities, $(1-\epsilon )\Vert r(\hat{x}){\Vert }_2^2\le (1+\epsilon )\Vert r(x_{\star }){\Vert }_2^2$, and dividing and taking the square root gives $\Vert r(\hat{x}){\Vert }_2\le \sqrt{(1+\epsilon )/(1-\epsilon )}\Vert r(x_{\star }){\Vert }_2$. The embedding preserves the residual geometry of the augmented column space, which is all the proof needs; it does not preserve the solution vector, explaining the caveat in the formal definition above. $\square$

Proposition (Gaussian matrices are subspace embeddings). Let $S\subseteq {\mathbb{R}}^m$ be a fixed subspace of dimension $n$ and let $\Pi =\frac{1}{\sqrt{d}}G\in {\mathbb{R}}^{d\times m}$ be Gaussian. If $d\ge C(n+\log (1/\delta ))/{\epsilon }^2$ for a universal constant $C$, then $\Pi$ is an $\epsilon$-subspace embedding for $S$ with probability at least $1-\delta$.

Proof. Choose an orthonormal basis $Z\in {\mathbb{R}}^{m\times n}$ for $S$, so $x\in S$ writes $x=Zy$ with $\Vert x{\Vert }_2=\Vert y{\Vert }_2$. The embedding condition becomes $(1-\epsilon )\Vert y{\Vert }_2\le \Vert \Pi Zy{\Vert }_2\le (1+\epsilon )\Vert y{\Vert }_2$ for all $y\in {\mathbb{R}}^n$, that is the $d\times n$ matrix $\Pi Z$ approximately preserves every norm. Let $\mathcal{N}$ be a half-net of the unit sphere $S^{n-1}$, which has size $|\mathcal{N}|\le 5^n$. By the Johnson-Lindenstrauss concentration of the theorem above, for a fixed $y$ on the unit sphere $\Vert \Pi Zy{\Vert }_2^2=\Vert \Pi (Zy){\Vert }_2^2$ lies in $[1-\epsilon /3,1+\epsilon /3]$ with probability at least $1-2\exp (-d{\epsilon }^2/12)$.

A union bound over the net requires $5^n\cdot 2\exp (-d{\epsilon }^2/12)\le \delta$, satisfied by $d\ge C(n+\log (1/\delta ))/{\epsilon }^2$. A standard net-extension argument lifts the estimate from net points to every unit vector by successive halving, giving $\Vert \Pi Zy{\Vert }_2\in [1-\epsilon ,1+\epsilon ]$ for all unit $y$, hence the embedding on all of $S$. $\square$

Connections Master

• The existence and uniqueness of the singular value decomposition, the Eckart-Young best-rank-$k$ approximation, and the dyadic expansion $A=\sum {\sigma }_j{u}_j{v}_j^{\mathsf{T}}$ are the foundation on which the randomised range finder rests; 01.01.12 proves them, and this unit specialises the projection $Q{Q}^{\mathsf{T}}A$ onto a random subspace as a probabilistic surrogate for the deterministic truncation $A_k$.

• The least-squares problem that sketch-and-solve compresses — the normal equations, the conditioning analysis ${\kappa }_2(A^{\mathsf{T}}A)={\kappa }_2(A{)}^2$, and the trade-off between solving fast and solving accurately — is the subject of 43.04.01, and the sketch-and-solve guarantee $\Vert A\hat{x}-b{\Vert }_2\le \sqrt{(1+\epsilon )/(1-\epsilon )}\Vert A{x}_{\star }-b{\Vert }_2$ is exactly the randomised counterpart of that unit's normal-equations and QR solvers, with the embedding replacing the full row reduction.

• The full SVD machinery that the sibling 43.05.08 deploys for total least squares and the generalised SVD is what randomisation accelerates: the randomised SVD produces the dominant singular triples of the augmented matrix $[A\mid b]$ at $O(mn(k+p))$ rather than $O(m{n}^2)$ cost, so the TLS fit of 43.05.08 can be computed from a sketch whenever the data is approximately low-rank.

• The Gaussian concentration and chi-squared tail bounds that power the Johnson-Lindenstrauss lemma belong to the probability substrate of 26.02.01 and 26.03.01; this unit imports the concentration inequality $\mathrm{Pr}[|\frac{1}{d}{\chi }_d^2-1|>\epsilon ]\le 2\exp (-d({\epsilon }^2/4-{\epsilon }^3/6))$ from that material and applies it with a union bound over $\left(\genfrac{}{}{0ex}{}{N}{2}\right)$ pairs.

• The randomised range finder is the block analogue of the Krylov subspace construction of 43.07.02: both build a subspace from repeated action of $A$ (or $A{A}^{\mathsf{T}}$) on a starting vector or block, and randomised block Lanczos and Arnoldi iterations, used for large eigenvalue problems, sit precisely at the junction this unit opens up.

Historical & philosophical context Master

The use of randomness to compress a high-dimensional geometric object originates with Johnson and Lindenstrauss, whose 1984 paper Extensions of Lipschitz mappings into a Hilbert space showed that any $N$ points in Euclidean space can be embedded into $O({\epsilon }^{-2}\log N)$ dimensions while preserving all pairwise distances within $(1\pm \epsilon )$ ^{[Johnson, W. B. & Lindenstrauss, J. — Extensions of Lipschitz mappings into a Hilbert space]}. The lemma sat in geometric functional analysis for two decades before the database community recognised it as a dimension-reduction tool: Papadimitriou, Raghavan, Tamaki and Vempala (2000) applied random projections to latent semantic indexing, and Achlioptas (2003) supplied the database-friendly sign matrices that made the projection cheap to apply.

The numerical-linear-algebra treatment crystallised with Sarlos, whose 2006 FOCS paper Improved approximation algorithms for large matrices via random projections first cast least squares and low-rank approximation as sketch-and-solve problems governed by the Johnson-Lindenstrauss mechanism ^{[Sarloes, T. — Improved approximation algorithms for large matrices via random projections]}. Halko, Martinsson and Tropp's 2011 SIAM Review paper Finding Structure with Randomness gave the modern form of the randomised range finder, the fixed-rank randomised SVD, the sharp Gaussian error constants, and the power-method stabilisation; Clarkson and Woodruff's 2013/2017 work on input-sparsity-time regression supplied the count-sketch embeddings that apply in time linear in the number of nonzeros. Martinsson and Tropp's 2020 Acta Numerica survey collected the foundations and the algorithms into the unified treatment this unit follows.

Bibliography Master

@article{halko-martinsson-tropp-2011,
  author  = {Halko, Nathan and Martinsson, Per-Gunnar and Tropp, Joel A.},
  title   = {Finding Structure with Randomness: Probabilistic Algorithms for Constructing Approximate Matrix Decompositions},
  journal = {SIAM Review},
  volume  = {53},
  number  = {2},
  pages   = {217--288},
  year    = {2011}
}

@incollection{johnson-lindenstrauss-1984,
  author    = {Johnson, William B. and Lindenstrauss, Joram},
  title     = {Extensions of {L}ipschitz mappings into a {H}ilbert space},
  booktitle = {Modern Analysis and Probability},
  series    = {Contemporary Mathematics},
  volume    = {26},
  publisher = {American Mathematical Society},
  year      = {1984},
  pages     = {189--206}
}

@book{woodruff-2014,
  author    = {Woodruff, David P.},
  title     = {Sketching as a Tool for Numerical Linear Algebra},
  series    = {Foundations and Trends in Theoretical Computer Science},
  volume    = {8},
  publisher = {now Publishers},
  year      = {2014}
}

@inproceedings{sarlos-2006,
  author    = {Sarl\'os, Tam\'as},
  title     = {Improved approximation algorithms for large matrices via random projections},
  booktitle = {Proc. 47th IEEE Symposium on Foundations of Computer Science (FOCS)},
  year      = {2006},
  pages     = {143--152}
}

@article{ailon-chazelle-2009,
  author  = {Ailon, Nir and Chazelle, Bernard},
  title   = {The fast {J}ohnson--{L}indenstrauss transform and algorithmic applications},
  journal = {SIAM Journal on Computing},
  volume  = {39},
  number  = {1},
  pages   = {330--349},
  year    = {2009}
}

@article{clarkson-woodruff-2017,
  author  = {Clarkson, Kenneth L. and Woodruff, David P.},
  title   = {Low-rank approximation and regression in input sparsity time},
  journal = {Journal of the ACM},
  volume  = {63},
  number  = {6},
  pages   = {Art. 54},
  year    = {2017}
}

@article{martinsson-tropp-2020,
  author  = {Martinsson, Per-Gunnar and Tropp, Joel A.},
  title   = {Randomized numerical linear algebra: Foundations and algorithms},
  journal = {Acta Numerica},
  volume  = {29},
  pages   = {403--572},
  year    = {2020}
}