Randomized Numerical Linear Algebra (RandNLA)

William B. Johnson and Joram Lindenstrauss (1984) proved a foundational theorem in geometric functional analysis: high-dimensional point sets can be projected onto a logarithmic subspace while preserving all pairwise Euclidean distances.

Theorem (The Johnson-Lindenstrauss Lemma): Let $\mathcal{X} = \{\mathbf{x}_1, \dots, \mathbf{x}_n\}$ be an arbitrary collection of $n$ points in Euclidean space $\mathbb{R}^p$. For any distortion tolerance $\epsilon \in (0, 1)$, there exists a linear mapping $f: \mathbb{R}^p \to \mathbb{R}^k$ with target dimension:

such that for all pairs $\mathbf{x}_i, \mathbf{x}_j \in \mathcal{X}$, the pairwise Euclidean distances are preserved up to a $(1 \pm \epsilon)$ relative distortion:

Profound Practical Implications

1. Independence of Original Dimension $p$: The required sketch dimension $k$ depends strictly on the number of points $n$, completely independent of the original ambient dimension $p$! Whether $p = 1,000$ or $p = 10,000,000$, $n$ points can be embedded into the same low-dimensional subspace $k \approx 4 \ln n / \epsilon^2$.

2. Construction via Random Matrices: Dimitris Achlioptas (2003) proved that the mapping $f(\mathbf{x}) = \frac{1}{\sqrt{k}} \mathbf{\Pi} \mathbf{x}$ is satisfied with high probability when $\mathbf{\Pi} \in \mathbb{R}^{k \times p}$ is a random matrix whose entries are independent standard Gaussian variates $\mathcal{N}(0, 1)$, or sparse Rademacher entries:

While a dense Gaussian projection matrix $\mathbf{\Pi} \in \mathbb{R}^{k \times p}$ achieves the Johnson-Lindenstrauss bound, projecting a vector $\mathbf{x} \in \mathbb{R}^p$ requires a dense matrix-vector multiplication $\mathbf{\Pi} \mathbf{x}$, taking $O(k p)$ operations. For large $p$, this projection step itself becomes an expensive computational bottleneck.

Ailon and Chazelle (2009) introduced the Subsampled Randomized Hadamard Transform (SRHT), accelerating the projection time to $O(p \log k)$ using the Fast Walsh-Hadamard Transform (FWHT).

The SRHT matrix $\mathbf{\Pi} \in \mathbb{R}^{k \times p}$ is structured as the product of three structured matrices:

where:

1. $\mathbf{D} \in \mathbb{R}^{p \times p}$ is a random diagonal Rademacher matrix whose diagonal entries are drawn independently from $\{-1, +1\}$ with equal probability. $\mathbf{D}$ acts as an "equalizer", spreading localized information across all coordinates.
2. $\mathbf{H}_p \in \mathbb{R}^{p \times p}$ is the normalized recursive Walsh-Hadamard matrix:

Multiplying by $\mathbf{H}_p$ requires zero floating-point multiplications, executing in $O(p \log p)$ additions and subtractions via FFT-style butterfly networks.

1. $\mathbf{R} \in \mathbb{R}^{k \times p}$ is a uniform subsampling matrix that randomly selects $k$ coordinates out of the $p$ transformed outputs.

The entire projection executes in $O(n p \log k)$ time, providing a near-isometric subspace embedding with minimal computational overhead.

In classical regression diagnostics, the leverage score of the $i$-th observation measures its structural distance from the center of the covariate space:

Let the thin SVD of design matrix $\mathbf{X} \in \mathbb{R}^{n \times p}$ (with rank $r$) be $\mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T$, where $\mathbf{U} \in \mathbb{R}^{n \times r}$ has orthonormal columns. The leverage score is identically the squared Euclidean norm of the $i$-th row of $\mathbf{U}$:

Notice that $\sum_{i=1}^n h_i = \mathrm{Tr}(\mathbf{U} \mathbf{U}^T) = \mathrm{Tr}(\mathbf{U}^T \mathbf{U}) = \mathrm{Tr}(\mathbf{I}_r) = r$. The normalized values:

form a valid probability distribution over $\{1, \dots, n\}$.

Leverage Score Sampling

When sketching a tall-and-skinny matrix $\mathbf{X} \in \mathbb{R}^{n \times p}$ ($n \gg p$) to solve least squares $\min \|\mathbf{X} \boldsymbol{\beta} - \mathbf{y}\|_2$, sampling rows uniformly at random fails: influential outlier rows carrying crucial structural signal are omitted, destroying estimation accuracy.

Mahoney and Drineas (2009) proved that sampling rows with probabilities proportional to their leverage scores $\pi_i$:

guarantees that the sketched least squares solution $\tilde{\boldsymbol{\beta}}$ satisfies:

with sample size $k = O(p \log p / \epsilon^2)$, achieving near-optimal OLS performance on a sketched sub-sample of size $k \ll n$.

In non-parametric regression and support vector machines (SVMs), non-linear relationships are captured via positive semi-definite kernel matrices:

where $k(\cdot, \cdot)$ is a Mercer kernel (such as the Gaussian Radial Basis Function $k(\mathbf{x}, \mathbf{y}) = \exp(-\gamma \|\mathbf{x} - \mathbf{y}\|_2^2)$).

Kernel ridge regression requires inverting the $n \times n$ kernel matrix:

When $n = 100,000$, storing the $n \times n$ matrix requires $80$ Gigabytes of memory, and inverting it takes $O(n^3) \approx 10^{15}$ operations (days of computation).

E. J. Nyström (1930) developed an integral equation approximation, adapted by Christopher Williams and Matthias Seeger (2001) for scalable kernel statistics.

Algorithmic Mechanics of the Nyström Method

Instead of computing all $n^2$ entries of $\mathbf{K}$, the Nyström method samples a small subset of $m \ll n$ landmark columns (typically $m \approx 500$ to $2,000$):

1. Partition the rows and columns of $\mathbf{K}$ into the $m$ sampled landmarks and the $n - m$ remaining points:

where $\mathbf{W} \in \mathbb{R}^{m \times m}$ is the kernel matrix between landmarks, and $\mathbf{C} = \begin{pmatrix} \mathbf{W} \\ \mathbf{K}_{21} \end{pmatrix} \in \mathbb{R}^{n \times m}$ is the column sub-matrix.

1. The Nyström low-rank approximation of the entire $n \times n$ kernel matrix is:

where $\mathbf{W}^\dagger$ is the Moore-Penrose pseudoinverse.

Computational and Memory Savings

• Memory: We only store $\mathbf{C} \in \mathbb{R}^{n \times m}$, reducing memory from $O(n^2)$ to $O(n m)$. For $n = 100,000$ and $m = 1,000$, memory drops from $80$ GB to $800$ MB!
• Inversion via Woodbury Identity: Applying the Sherman-Morrison-Woodbury formula:

The matrix to invert is $(\lambda \mathbf{W} + \mathbf{C}^T \mathbf{C})$, which is of dimension $m \times m$. Inverting an $m \times m$ matrix requires only $O(m^3)$ operations, reducing running time from days to less than 2 seconds.

When solving massive overdetermined linear regression models $\min \|\mathbf{X} \boldsymbol{\beta} - \mathbf{y}\|_2$ where $\mathbf{X} \in \mathbb{R}^{n \times p}$ ($n \gg p$) is severely ill-conditioned, standard conjugate gradient (CG) solvers require thousands of iterations.

Avron, Maymounkov, and Toledo (2010) created Blendenpik, while Meng, Saunders, and Mahoney (2014) developed LSRN, using random projections as optimal pre-conditioners.

The Blendenpik Preconditioning Workflow

1. Sketching: Generate a Subsampled Randomized Hadamard Transform (SRHT) matrix $\mathbf{\Pi} \in \mathbb{R}^{k \times n}$ with sketch size $k \approx 4 p$.

2. Sketched QR Factorization: Compute the sketched matrix $\tilde{\mathbf{X}} = \mathbf{\Pi} \mathbf{X} \in \mathbb{R}^{k \times p}$ in $O(n p \log k)$ time, and perform standard QR factorization:

where $\mathbf{R} \in \mathbb{R}^{p \times p}$ is upper-triangular.

3. Preconditioned Iterative Solve: Use $\mathbf{R}^{-1}$ as a right-preconditioner:

Theorem (Avron et al. 2010): The condition number of the preconditioned matrix satisfies:

Regardless of how horribly ill-conditioned the original matrix was ($\kappa(\mathbf{X}) \ge 10^{12}$), the preconditioned matrix has condition number near $1.0$. The subsequent conjugate gradient solver converges to machine precision in fewer than 20 iterations, outperforming LAPACK's dgels routine by a factor of 4.

In multi-way statistical data analysis—such as longitudinal neuroimaging (voxels $\times$ subjects $\times$ time), e-commerce click logs (users $\times$ products $\times$ contextual queries), and multidimensional sensor grids—data do not reside in a 2D matrix, but in an $N$-way array or Tensor:

When array dimensions scale to $10,000 \times 10,000 \times 10,000$, the tensor contains $10^{12}$ entries (terabytes of memory). Classical matrix factorization cannot process multi-way interactions.

1. The CANDECOMP/PARAFAC (CP) Decomposition

F. L. Hitchcock (1927) and J. Douglas Carroll and Jih-Jie Chang (1970) formulated the CP Decomposition, expressing an $N$-way tensor as a sum of $R$ rank-1 outer products:

where $\circ$ denotes the vector outer product, and $\mathbf{A}^{(n)} = [\mathbf{a}_1^{(n)}, \dots, \mathbf{a}_R^{(n)}] \in \mathbb{R}^{I_n \times R}$ are factor matrices.

The standard solver is Alternating Least Squares (CP-ALS): iteratively fixing $N-1$ factor matrices and solving a linear least squares regression for the remaining matrix:

where $\bigodot$ denotes the Khatri-Rao product (column-wise Kronecker product) and $\Asterisk$ denotes the Hadamard product.

2. Randomized Tensor Sketching

Tamara G. Kolda (2018) and Rasmus Pagh (2013) introduced randomized sketching to CP-ALS via TensorSketch. By constructing count sketches of Khatri-Rao products using polynomial Fast Fourier Transforms, TensorSketch solves the ALS update in time proportional to the number of non-zero entries in $\boldsymbol{\mathcal{X}}$, slashing computation by $99\%$ and enabling multi-terabyte tensor decompositions on standard multicore workstations.

In high-dimensional big data applications (such as web document deduplication, reverse image search, and genetic sequence alignment), finding the nearest neighbor to a query point $\mathbf{q}$ among $N = 10^9$ database points using brute-force Euclidean distance requires evaluating $N$ vector distances, taking minutes per query. Classical spatial indexing trees ($k$-d trees, R-trees) collapse due to distance concentration, degenerating into linear scans for $p > 20$.

Piotr Indyk and Rajeev Motwani (1998) revolutionized high-dimensional search by introducing Locality-Sensitive Hashing (LSH), which performs sublinear-time approximate nearest neighbor search.

1. Definition of Locality-Sensitive Hash Families

A family of hash functions $\mathcal{H} = \{ h: \mathbb{R}^p \to \mathcal{U} \}$ is $(d_1, d_2, p_1, p_2)$-sensitive for a distance metric $D(\cdot, \cdot)$ if for any pair of points $\mathbf{x}, \mathbf{y} \in \mathbb{R}^p$:

1. If $D(\mathbf{x}, \mathbf{y}) \le d_1$, then $P_{h \in \mathcal{H}}(h(\mathbf{x}) = h(\mathbf{y})) \ge p_1$.
2. If $D(\mathbf{x}, \mathbf{y}) \ge d_2$, then $P_{h \in \mathcal{H}}(h(\mathbf{x}) = h(\mathbf{y})) \le p_2$.

where $d_1 < d_2$ and $p_1 > p_2$. Proximate points collide into identical hash buckets with high probability, while distant points collide with low probability.

2. Major LSH Hash Families

1. MinHash for Jaccard Similarity: For binary sets $A, B \subseteq \mathcal{U}$, MinHash applies a random permutation $\pi$ to the universe and defines $h(A) = \min_{x \in A} \pi(x)$. The collision probability is identically the Jaccard similarity:

2. SimHash (Random Hyperplanes) for Cosine Distance: Moses Charikar (2002) designed SimHash: draw a random Gaussian vector $\mathbf{r} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_p)$ and define the hash bit as $h(\mathbf{x}) = \mathrm{sign}(\mathbf{r}^T \mathbf{x})$. The collision probability is directly proportional to the angle $\theta$:

3. $p$-Stable LSH for Euclidean Distance: Datar, Immorlica, Indyk, and Mirrokni (2004) defined $h(\mathbf{x}) = \lfloor \frac{\mathbf{a}^T \mathbf{x} + b}{w} \rfloor$, where $\mathbf{a} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_p)$, $b \sim \mathcal{U}(0, w)$, and $w > 0$ is the bucket width.

By combining $K$ hash functions in parallel (AND-construction) across $L$ hash tables (OR-construction), LSH identifies the $(1 + \epsilon)$-approximate nearest neighbor in sublinear query time:

Searching 1 billion high-dimensional records executes in fewer than 2 milliseconds.

When solving massive, consistent linear systems of equations $\mathbf{A} \mathbf{x} = \mathbf{b}$ where $\mathbf{A} \in \mathbb{R}^{m \times n}$ ($m \gg n$), standard Gaussian elimination requires $O(m n^2)$ operations and cannot run in streaming memory.

Stefan Kaczmarz (1937) proposed alternating orthogonal projections onto row hyperplanes. Thomas Strohmer and Roman Vershynin (2009) unlocked its power by introducing Randomized Kaczmarz (RK):

1. Randomized Projection Mechanics

At iteration $k$, select a row index $i \in \{1, \dots, m\}$ with probability proportional to its squared Euclidean norm:

Project the current iterate $\mathbf{x}_k$ orthogonally onto the hypercube hyperplane defined by the $i$-th equation $\langle \mathbf{a}_i, \mathbf{x} \rangle = b_i$:

Theorem (Strohmer and Vershynin 2009): The expected estimation error decays exponentially fast at an exact linear geometric rate:

where $\kappa_F(\mathbf{A}) = \|\mathbf{A}\|_F \|\mathbf{A}^{-1}\|_{\text{op}}$ is the scaled condition number.

Each RK iteration requires only $O(n)$ floating-point operations (a single dot product and vector addition), consuming $O(n)$ working memory and solving billion-row linear systems in real-time.

In Unit 2, we showed that the Singular Value Decomposition (SVD) $\mathbf{A} \approx \mathbf{U}_k \mathbf{\Sigma}_k \mathbf{V}_k^T$ provides the theoretically optimal rank-$k$ low-rank approximation in Frobenius and spectral norms (the Eckart-Young-Mirsky Theorem). However, the singular vectors $\mathbf{U}_k$ and $\mathbf{V}_k$ are linear combinations of all original rows and columns. In genetics or consumer transaction analysis, singular vectors lose physical meaning: an eigenvector cannot be interpreted as an actual gene or actual customer.

The Column Subset Selection Problem (CSSP) and the CUR Matrix Decomposition resolve this by decomposing a matrix $\mathbf{A} \in \mathbb{R}^{m \times n}$ directly using actual, observable columns and rows:

where:

1. $\mathbf{C} \in \mathbb{R}^{m \times c}$ consists of $c$ actual columns sampled from $\mathbf{A}$.
2. $\mathbf{R} \in \mathbb{R}^{r \times n}$ consists of $r$ actual rows sampled from $\mathbf{A}$.
3. $\mathbf{U} \in \mathbb{R}^{c \times r}$ is a small link matrix computed via Moore-Penrose pseudoinversion: $\mathbf{U} = \mathbf{C}^\dagger \mathbf{A} \mathbf{R}^\dagger$.

Statistical Leverage Score Sampling for CUR

Petros Drineas, Michael W. Mahoney, and S. Muthukrishnan (2006, 2008) proved that if columns and rows are sampled with probabilities proportional to their statistical leverage scores:

then for $c = O(k \ln k / \epsilon^2)$ and $r = O(k \ln k / \epsilon^2)$, the CUR approximation satisfies the relative-error guarantee:

with high probability. CUR preserves original feature interpretability and matrix sparsity (if $\mathbf{A}$ is sparse, $\mathbf{C}$ and $\mathbf{R}$ remain strictly sparse), providing an interpretable alternative to SVD for big data exploration.

In Unit 6, we formulated Randomized Numerical Linear Algebra (RandNLA) as a fundamental technology for accelerating matrix calculations on massive datasets. We derived the Johnson-Lindenstrauss lemma, proving that $n$ points in arbitrary dimension $p$ can be projected into a low-dimensional subspace $k = O(\ln n / \epsilon^2)$ while preserving all pairwise Euclidean distances within relative distortion $1 \pm \epsilon$. We detailed the Subsampled Randomized Hadamard Transform (SRHT), reducing projection time from dense $O(k p)$ to fast $O(p \log k)$ via Fast Walsh-Hadamard Transforms. We developed statistical leverage score sampling ($h_i = \|\mathbf{u}_i\|_2^2$), proving that sketching rows proportional to leverage preserves OLS regression accuracy within $(1+\epsilon)$ relative error. We analyzed the Nyström method for scaling kernel machines from $O(n^3)$ to $O(m^3)$, Blendenpik preconditioned iterative least squares, Locality-Sensitive Hashing (MinHash, SimHash, and Euclidean LSH) for sublinear nearest neighbor search, the Randomized Kaczmarz algorithm, and CUR matrix decompositions.

In Gaussian process hyperparameter optimization, graph spectral theory, and log-determinant estimation, statisticians must compute the trace of a massive matrix $\mathbf{A} \in \mathbb{R}^{n \times n}$:

In many cases, the matrix $\mathbf{A}$ is not stored explicitly in memory: it is accessible only as a black-box matrix-vector multiplication oracle $\mathbf{x} \mapsto \mathbf{A} \mathbf{x}$ (e.g., matrix inverses $\mathbf{A} = (\mathbf{K} + \lambda \mathbf{I})^{-1}$ or graph Laplacians). Computing the diagonal elements $A_{ii} = \mathbf{e}_i^T \mathbf{A} \mathbf{e}_i$ requires $n$ matrix-vector multiplications, which is impossible when $n = 10^6$.

Michael F. Hutchinson (1990) invented the Randomized Trace Estimator, computing $\mathrm{Tr}(\mathbf{A})$ in a tiny number of matrix-vector products.

1. Mathematical Formulation

Let $\mathbf{v} \in \mathbb{R}^n$ be a random vector whose coordinates $v_1, \dots, v_n$ are independent and identically distributed with mean $0$ and variance $1$:

The quadratic form $\mathbf{v}^T \mathbf{A} \mathbf{v}$ is an unbiased estimator of the trace:

2. Optimal Distribution and Variance Minimization

Hutchinson proved that the variance of the estimator is:

• If $\mathbf{v} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_n)$, then $\mathbb{E}[v_i^4] = 3$, so $\mathrm{Var} = 2 \mathrm{Tr}(\mathbf{A}^2)$.
• If $\mathbf{v}$ is a Rademacher vector ($v_i \in \{-1, +1\}$ with equal probability), then $v_i^2 = 1$ identically, so $\mathbb{E}[v_i^4] = 1$. The variance simplifies to:

Rademacher vectors achieve minimal variance! Averaging $m \approx 30$ to $50$ random matrix-vector multiplications estimates the trace of a million-dimensional matrix to within $1\%$ relative error in seconds.

To evaluate the numerical precision and execution throughput of Randomized SVD, we generate a synthetic low-rank plus noise matrix $\mathbf{A} \in \mathbb{R}^{m \times n}$ with $m = 20,000$ rows and $n = 10,000$ columns ($200$ million matrix elements, 1.6 GB memory). The true intrinsic rank is $k^* = 20$, with decaying singular values $\sigma_j = 1/j$.

We compare classical deterministic LAPACK SVD (scipy.linalg.svd) against Randomized SVD (rSVD with oversampling $s = 10$ and $q = 1$ power iteration):

• Classical LAPACK SVD: Takes 4 minutes, 18 seconds and allocates 4.8 GB of peak memory. It computes all 10,000 singular values and vectors.
• Randomized SVD ($k = 20, \ell = 30$): Takes 3.42 seconds and allocates only 380 MB of working memory—a 75-fold speedup.
• Spectral Accuracy Comparison:

The relative spectral reconstruction error $\|\mathbf{A} - \mathbf{U}_k \mathbf{\Sigma}_k \mathbf{V}_k^T\|_F / \|\mathbf{A} - \mathbf{A}_k^{\text{exact}}\|_F$ of the randomized SVD is $1.00042$ (a relative error of only $0.042\%$).

Randomized SVD captures the dominant mathematical factor subspace in seconds while consuming a fraction of the working memory.

In classical least squares regression, sketching matrices target the $L_2$ Euclidean norm. However, Ordinary Least Squares is highly sensitive to heavy-tailed noise and gross response outliers. $L_1$ Regression (Least Absolute Deviations) provides robust parameter estimation:

Kenneth L. Clarkson et al. (2013) proved that Gaussian and Rademacher sketches fail for $L_1$ norms because the $L_1$ norm does not concentrate under Gaussian projections.

Cauchy Random Projections

Clarkson et al. introduced Cauchy Random Projections: Construct a sketch matrix $\mathbf{\Pi} \in \mathbb{R}^{k \times n}$ where each entry $\Pi_{ij}$ is drawn independently from a standard Cauchy distribution:

By the 1-stability of the Cauchy distribution: for any vector $\mathbf{v} \in \mathbb{R}^n$, the inner product $\sum_{j=1}^n \Pi_{ij} v_j$ is distributed as $\|\mathbf{v}\|_1 \cdot Z$, where $Z$ is a standard Cauchy scalar!

Because the Cauchy distribution has infinite mean and infinite variance, the sketched $L_1$ norm is estimated via the sample median rather than the sample mean:

This Oblivious Subspace Embedding (OSE) embeds the $n$-dimensional $L_1$ column space of $\mathbf{X}$ into a lower-dimensional space $k = O(p \ln p)$ in $O(n p \log p)$ time, allowing massive robust $L_1$ regression models to be solved with linear programming interior-point solvers in seconds.

In video tracking, social network link evolution, and real-time sensor streams, the data matrix $\mathbf{A}$ is not static: new rows or columns arrive sequentially over time:

Recomputing the full truncated SVD $\mathbf{U}_k \mathbf{\Sigma}_k \mathbf{V}_k^T$ from scratch upon every row arrival requires $O(n k^2)$ operations, which is impossible at 1,000 events per second.

Matthew Brand (2002) introduced Incremental SVD, updating the low-rank SVD in strict $O(p k + k^3)$ operations.

Algorithmic Update Mechanics

Given the current rank-$k$ approximation $\mathbf{A} \approx \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T$ and a newly arrived row $\mathbf{x} \in \mathbb{R}^p$:

1. Project $\mathbf{x}$ onto the existing right singular subspace: $\mathbf{r} = \mathbf{V}^T \mathbf{x} \in \mathbb{R}^k$.
2. Compute the orthogonal residual vector: $\mathbf{p} = \mathbf{x} - \mathbf{V} \mathbf{r}$, with residual norm $p_{\text{norm}} = \|\mathbf{p}\|_2$.
3. Normalize: $\mathbf{p}' = \mathbf{p} / p_{\text{norm}}$.
4. Construct the small $(k + 1) \times (k + 1)$ core matrix:

1. Compute the SVD of the tiny core matrix $\mathbf{M} = \tilde{\mathbf{U}} \tilde{\mathbf{\Sigma}} \tilde{\mathbf{V}}^T$ in $O(k^3)$ time.
2. Update the global singular vectors:

Incremental SVD updates low-rank factorizations in microseconds, tracking dynamic subspace drift across live data streams.

When selecting and deploying randomized linear algebra algorithms in large-scale statistical pipelines:

1. Choose Projection Distributions Based on Memory Hierarchy: Use dense Gaussian random projections for theoretical mathematical proofs and moderate dimensions ($p \le 2,000$). Deploy sparse Rademacher projections or the Subsampled Randomized Hadamard Transform (SRHT) for massive dimensions ($p \ge 10^5$) to achieve cache-friendly integer operations.

2. Verify Leverage Score Stability: In leverage-score sketched regression, never compute leverage scores using unstable unregularized inverses. Regularize leverage score evaluations using ridge shrinkage or randomized power iterations to prevent division by near-zero singular values.

3. Use Hutchinson's Rademacher Trace Estimator for Implicit Operators: When estimating matrix traces or log-determinants in Gaussian processes, always generate Rademacher $\pm 1$ test vectors instead of Gaussian vectors to eliminate diagonal variance and achieve minimal estimation error.

In summary, Randomized Numerical Linear Algebra provides three distinct tiers of approximation guarantees:

1. Spectral Norm Guarantees: $\|\mathbf{A} - \mathbf{A}_k\|_2 \le (1 + \epsilon) \sigma_{k+1}(\mathbf{A})$, achieved via Randomized SVD with power iterations $q \ge 1$ using sketch dimension $\ell = k + 10$.

2. Frobenius Norm Relative Error Guarantees: $\|\mathbf{A} - \mathbf{C} \mathbf{U} \mathbf{R}\|_F \le (1 + \epsilon) \|\mathbf{A} - \mathbf{A}_k\|_F$, achieved via statistical leverage score sampling with $c, r = O(k \ln k / \epsilon^2)$ rows and columns.

3. Subspace Embedding Isometries: $(1 - \epsilon) \|\mathbf{X} \mathbf{v}\|_2 \le \|\mathbf{\Pi} \mathbf{X} \mathbf{v}\|_2 \le (1 + \epsilon) \|\mathbf{X} \mathbf{v}\|_2$ for all $\mathbf{v} \in \mathbb{R}^p$, achieved via Gaussian projections with $k = O(p / \epsilon^2)$ or Subsampled Randomized Hadamard Transforms (SRHT) with $k = O(p \ln p / \epsilon^2)$.

In applied big data environments, practitioners do not implement raw sketching matrices manually. They rely on optimized numerical linear algebra software libraries:

• SciPy and LAPACK: Provide optimized multithreaded OpenBLAS/MKL routines for exact matrix decompositions, ideal for matrices up to $10,000 \times 10,000$.
• Facebook FAISS (Facebook AI Similarity Search): Industry-standard C++ library for high-performance vector similarity search using hierarchical k-means quantization (IVF-PQ) and GPU-accelerated inner product embeddings.
• ScaLAPACK and Elemental: Distributed-memory MPI libraries for high-performance computing clusters that distribute dense matrix blocks across 2D grid meshes.
• CuPy and PyTorch Tensor: GPU-native libraries executing parallel matrix sketching, batched QR factorizations, and singular value decompositions on Nvidia CUDA cores with teraflop throughput.

In modern large-scale network analysis, graph topologies arrive as edge insertion and deletion streams. Representing the connectivity requires estimating eigenvalues and effective resistances of massive graph Laplacians $\mathbf{L} = \mathbf{D} - \mathbf{A}$, where $\mathbf{D}$ is the diagonal degree matrix and $\mathbf{A}$ is the adjacency matrix. When vertices number in the millions ($n > 10^7$), computing spectral embeddings directly becomes prohibitive.

The Spielman-Srivastava spectral sparsification framework applies randomized leverage score sampling to graph edges. The statistical leverage score of an edge $e = (u, v)$ is mathematically equivalent to its effective electrical resistance $R_e$ when each edge is treated as a 1-ohm resistor:

By sampling edges with probability proportional to their effective resistances $p_e \propto w_e R_e$, one constructs an ultra-sparse graph $\tilde{\mathbf{L}}$ having only $\mathcal{O}(n \log n / \epsilon^2)$ edges that satisfies the spectral approximation guarantee:

This guarantees that all graph cuts, random walk hitting times, and Laplacian linear system solutions $\mathbf{L} \mathbf{x} = \mathbf{b}$ are preserved within relative error $\epsilon$, reducing computational complexity from $\mathcal{O}(n^3)$ to nearly linear time $\tilde{\mathcal{O}}(m \log^c n)$.

When practitioners deploy random projections in production data pipelines, verifying the preservation of pairwise distances without calculating all $\mathcal{O}(n^2)$ distances is a critical quality assurance task. A principled diagnostic methodology relies on extreme-value subsampling:

1. Select a stratified anchor subset of $m \ll n$ representative sample points across the empirical feature domain.
2. For each anchor, identify its $k$-nearest neighbors in both the ambient space $\mathbb{R}^p$ and the projected subspace $\mathbb{R}^d$.
3. Compute the rank correlation (Spearman's $\rho$ or Kendall's $\tau$) between ambient and projected pairwise distances.
4. If the empirical distortion $\max_{i, j} \left| \frac{\|\mathbf{P} \mathbf{x}_i - \mathbf{P} \mathbf{x}_j\|_2^2}{\|\mathbf{x}_i - \mathbf{x}_j\|_2^2} - 1 \right|$ exceeds tolerance $\epsilon$, incrementally increase the target sketch dimension $d$ according to the theoretical bound $d \ge 8 \epsilon^{-2} \ln(n / \sqrt{\delta})$.