Randomized Numerical Linear Algebra

Matrices Too Large to Hold Still

Suppose a scientist has a table with a billion observations and thousands of measurements per observation. Classical matrix algorithms can be exact and elegant, yet repeatedly scanning or decomposing that table may dominate the work. Randomized numerical linear algebra replaces the original matrix with a smaller random view that preserves the subspace, product, or low-rank structure needed by the downstream calculation.

The everyday analogy is a carefully designed survey. A naive survey samples people uniformly and may miss a small but influential group. A leverage-aware survey gives higher probability to observations that span unusual directions. A projection takes another route: it mixes all rows or columns into fewer randomized combinations, so every original observation can influence the reduced problem.

For a tiny example, let

$$ A=\begin{bmatrix}1&0\\0&1\\1&1\end{bmatrix},\qquad b=\begin{bmatrix}1\\2\\3\end{bmatrix}. $$

Here $b=A(1,2)^\mathsf T$, so the least-squares residual is zero. Keeping only the first two rows still recovers the solution, but keeping two identical or nearly redundant directions need not. The point is not that row sampling always works; it works when the sample or sketch preserves the column space relevant to every candidate residual.

Matrix multiplication sketches approximate $A^\mathsf TB$ without forming every contribution. Leverage scores identify influential rows or columns. Subspace embeddings preserve norms throughout a fixed subspace. Randomized SVD finds an approximate range before computing a small factorization. Frequent Directions is deterministic: it repeatedly shrinks singular values of a compact matrix and gives covariance-error guarantees in a streaming setting.

Subspace Views of Low-Rank Structure

Let $A\in\mathbb R^{m\times n}$ have singular value decomposition $A=U\Sigma V^\mathsf T$. The best rank-$k$ approximation under both spectral and Frobenius norms is $A_k=U_k\Sigma_kV_k^\mathsf T$. Its errors are

$$ \lVert A-A_k\rVert_2=\sigma_{k+1},\qquad \lVert A-A_k\rVert_F^2=\sum_{j>k}\sigma_j^2. $$

A randomized range finder draws a test matrix $\Omega\in\mathbb R^{n\times(k+p)}$, forms $Y=A\Omega$, and computes an orthonormal basis $Q$ for the columns of $Y$. It then works with the small matrix $B=Q^\mathsf TA$. If $Q$ captures the leading left singular subspace, $QB$ approximates $A$. Oversampling $p$ reduces the risk that random test vectors poorly excite an important direction. Power iterations use $(AA^\mathsf T)^qA\Omega$ to amplify spectral gaps, at the cost of passes and possible numerical stabilization needs.

For a full-column-rank matrix $A$, the leverage scores of its rows are the diagonal entries of the projector $UU^\mathsf T$, where $U$ spans $\operatorname{col}(A)$. They satisfy $0\le\ell_i\le1$ and $\sum_i\ell_i=\operatorname{rank}(A)$. Sampling rows proportionally to approximate leverage scores protects rare but geometrically important directions better than uniform sampling when coherence is high.

A $(1\pm\varepsilon)$ subspace embedding $S$ for $\operatorname{col}(A)$ satisfies

$$ (1-\varepsilon)\lVert Ax\rVert_2^2\le\lVert SAx\rVert_2^2\le(1+\varepsilon)\lVert Ax\rVert_2^2 $$

for every $x$. For least squares, the relevant subspace often includes both $A$ and $b$, because residuals are $Ax-b$. Preserving only columns of $A$ individually does not preserve every linear combination.

Frequent Directions processes rows into a fixed-size matrix $B$. When full, it computes an SVD, subtracts a common squared singular-value threshold, and discards the directions shrunk to zero. It never overstates directional energy: for every unit vector $x$, $0\le\lVert Ax\rVert_2^2-\lVert Bx\rVert_2^2$ is bounded in terms of the tail energy and sketch size under the chosen theorem.

flowchart TD
    A[Large matrix A] --> B{Target task}
    B -->|Low rank| C[Random test matrix Omega]
    C --> D[Range sample Y equals A Omega]
    D --> E[Orthonormal basis Q]
    E --> F[Small matrix Q transpose A]
    F --> G[Approximate SVD]
    B -->|Regression| H[Embed span of A and b]
    H --> I[Solve sketched least squares]
    B -->|Row sampling| J[Estimate leverage scores]
    J --> K[Sample and rescale influential rows]
    B -->|Streaming covariance| L[Frequent Directions]
    L --> M[Deterministic shrinkage]
    G --> N[Low-rank error guarantee]
    I --> O[Residual guarantee]
    K --> O
    M --> P[Directional covariance guarantee]
    classDef navy fill:#17365D,color:#fff,stroke:#17365D;
    classDef teal fill:#147D84,color:#fff,stroke:#0D565B;
    classDef crimson fill:#C6284D,color:#fff,stroke:#861B35;
    classDef purple fill:#6C3BAA,color:#fff,stroke:#472675;
    class A,B navy;
    class C,D,E,H,J,L purple;
    class F,I,K,M teal;
    class G,N,O,P crimson;

Sampling, Projection and Spectral Error

For approximate matrix multiplication, write $A^\mathsf TB$ as a sum of outer products of corresponding rows. Sample row index $i$ with probability $p_i$, and use the rescaled outer product $a_i^\mathsf Tb_i/p_i$. Its expectation is the full product. Repeating and averaging reduces variance. Probabilities proportional to suitable row-norm products or leverage-like quantities reduce the influence of any one rescaled sample. Unbiasedness alone does not control spectral or Frobenius error.

For regression, suppose $S$ preserves the span of $A$ and $b$. Every residual $Ax-b$ lies in that subspace, so

$$ (1-\varepsilon)\lVert Ax-b\rVert_2^2\le\lVert S(Ax-b)\rVert_2^2\le(1+\varepsilon)\lVert Ax-b\rVert_2^2. $$

If $\widetilde x$ minimizes the sketched residual and $x_\star$ minimizes the original, comparing the two inequalities yields an approximate original objective. The proof uses preservation for all candidate residuals, not only for $x_\star$.

Randomized low-rank analysis decomposes error into the optimal tail $A-A_k$ and failure to capture the leading range. Spectral error focuses on the worst direction and is sensitive to a single missed singular vector; Frobenius error aggregates energy and may remain small when one direction is moderately wrong. Always match the norm to the application.

Frequent Directions offers a deterministic alternative. When a sketch $B$ fills, compute $B=U\Sigma V^\mathsf T$, choose a threshold $\delta$ from a trailing squared singular value, and replace singular values by $\sqrt{\max(\sigma_i^2-\delta,0)}$. The removed energy is charged uniformly across directions, yielding a telescoping covariance-error bound. Changing the threshold or omitting rescaling changes the theorem.

flowchart TD
    A[State target: product, regression, or rank k] --> B[Choose sketch family and size]
    B --> C{Sampling or projection?}
    C -->|Sampling| D[Probabilities cover influential rows or columns]
    D --> E[Rescale sampled contributions]
    E --> F[Unbiased estimator plus concentration]
    C -->|Projection| G[Embed the relevant fixed subspace]
    G --> H[Preserve every candidate residual]
    H --> I[Compare sketched optimum with original optimum]
    B --> J{Streaming deterministic option?}
    J -->|Yes| K[Frequent Directions shrinkage]
    K --> L[Charge discarded energy across directions]
    F --> M[Norm-specific product bound]
    I --> N[Regression objective bound]
    L --> O[Covariance and low-rank bound]
    D -. uniform sampling misses high leverage .-> X[Failure]
    G -. embed only individual columns .-> Y[Failure]
    classDef navy fill:#17365D,color:#fff,stroke:#17365D;
    classDef teal fill:#147D84,color:#fff,stroke:#0D565B;
    classDef crimson fill:#C6284D,color:#fff,stroke:#861B35;
    classDef purple fill:#6C3BAA,color:#fff,stroke:#472675;
    class A,B,C,J navy;
    class D,G,K purple;
    class E,F,H,I,L teal;
    class M,N,O teal;
    class X,Y crimson;

PCA, Regression and Scientific Data

For an approachable PCA task, randomized range finding can compress a tall matrix of student measurements before an SVD. The leading components summarize variance, but interpreting them as causal traits remains unjustified. Centering, scaling, missing values, and categorical encodings define the matrix and can matter more than the sketch.

For large regression, a subspace embedding reduces billions of rows to a manageable problem. It can preserve the residual objective while coefficients remain sensitive if the design is ill-conditioned. Ridge regularization changes the relevant augmented system and must be included in the sketching analysis.

For scientific data, low-rank structure can reveal latent modes in telescope, climate, genomics, or simulation matrices. A weak spectral gap may make individual singular vectors unstable even when the low-rank subspace and reconstruction error are stable. Domain conclusions should reflect subspace uncertainty, not overinterpret one basis orientation.

Adversarial or misleading cases include rare high-leverage observations, drift after sampling probabilities are estimated, corrupted rows with huge norms, and adaptive workloads that respond to a fixed public sketch. A global Frobenius guarantee can hide one scientifically critical direction; a small regression residual can coexist with biased data collection.

Design a Randomized Matrix Workflow

Start from the output guarantee, not the fashionable sketch. Specify the matrix shape and access pattern, target task, rank or subspace dimension, norm, relative or additive error, failure probability, passes, precision, and adversary model. Then select sampling, projection, randomized range finding, or deterministic shrinkage and validate the result on the original objective.