The modern algorithmic toolbox: sketching, spectral methods, and randomization

Why this subject matters now

The datasets that machine learning operates on are wide as well as tall. A document embedding lives in several thousand dimensions, a design matrix for a regression can have millions of columns, and a graph of users or citations can have billions of edges. The exact algorithms that a first course teaches, Gaussian elimination in \( O(n^3) \), the singular value decomposition in \( O(mn^2) \), the exact minimum cut, all have the property that they scale badly enough to be unusable at the sizes that matter. What changed over the past two decades is not that these problems got easier but that a body of technique made it acceptable to solve them approximately, with a randomized or spectral surrogate whose error is bounded by a theorem rather than by hope.

The unifying move is dimensionality reduction in a broad sense. A random projection replaces a high-dimensional point cloud with a low-dimensional one that has the same pairwise distances. A sketch replaces a tall matrix with a short one that spans the same column space. A low-rank factorization replaces a full matrix with the handful of directions that carry its energy. A graph Laplacian replaces a combinatorial cut problem with an eigenvalue problem. In each case a small object stands in for a large one, and the value of the technique is precisely the guarantee that the substitution is safe. A practitioner today is expected to know not just that scikit-learn has a random_projection module but why the target dimension is \( O(\log n / \varepsilon^2) \) and independent of the original dimension, why a randomized SVD needs a power iteration to be accurate on a slowly decaying spectrum, and why the second eigenvector of a Laplacian gives a good cut. Those are the derivations below.

The through-line is the probabilistic method: to show a good object exists, exhibit a random one and prove that it is good with positive probability. Every guarantee on this page reduces to a concentration inequality applied to a well-chosen random variable. That is why the toolbox is coherent rather than a bag of tricks, and why learning the concentration bounds first pays off across all of it.

Core theory

The probabilistic toolkit

Three facts do most of the work. The first is linearity of expectation: for any random variables, \( \E\!\left[\sum_i X_i\right] = \sum_i \E[X_i] \), with no independence required. Paired with indicator variables, where \( X_i = \mathbf{1}[\text{event } i] \) and \( \E[X_i] = \P[\text{event } i] \), it turns a count into a sum of probabilities. If \( X = \sum_i X_i \) counts how many of \( n \) events occur, then \( \E[X] = \sum_i \P[\text{event } i] \) regardless of how the events correlate. This single identity proves the expected number of comparisons in randomized quicksort, the expected number of empty bins in a hashing scheme, and the expected size of a randomized rounding solution below.

The second fact is that expectation alone is weak; we usually need to know that \( X \) is close to \( \E[X] \) with high probability. Three inequalities give increasing strength at the cost of increasing assumptions.

Markov. For a nonnegative random variable and any \( a \gt 0 \), \( \P[X \ge a] \le \E[X]/a \). This is the weakest possible tail bound and uses only the mean.

Chebyshev. Applying Markov to \( (X - \E[X])^2 \) gives \( \P\big[\,|X - \E[X]| \ge a\,\big] \le \Var(X)/a^2 \). This is polynomial decay in \( a \) and needs the variance. When \( X = \sum_i X_i \) is a sum of pairwise-independent variables the variance is additive, \( \Var(X) = \sum_i \Var(X_i) \), which is already enough to make many estimators concentrate.

Chernoff. For a sum \( X = \sum_{i=1}^n X_i \) of independent random variables the moment generating function factorizes, \( \E[e^{sX}] = \prod_i \E[e^{sX_i}] \), and Markov applied to \( e^{sX} \) gives an exponentially decaying tail after optimizing over \( s \). For independent \( X_i \in \{0,1\} \) with \( \mu = \E[X] \), the standard multiplicative form is

The Chernoff bound is the workhorse of randomized algorithms because it turns "the expected value is right" into "the value is right with probability \( 1 - 1/\text{poly}(n) \)", which is what a union bound over polynomially many events can absorb. The proof of Johnson-Lindenstrauss below is a Chernoff bound on a sum of squared Gaussians.

The Johnson-Lindenstrauss lemma

The lemma states that any \( n \) points in Euclidean space of arbitrary dimension can be mapped into \( k = O(\log n / \varepsilon^2) \) dimensions so that all pairwise distances are preserved up to a factor \( 1 \pm \varepsilon \). The target dimension depends on the number of points and the accuracy, and remarkably not on the original dimension \( d \). The original 1984 result of Johnson and Lindenstrauss was existential; the modern proof, due to Dasgupta and Gupta (2003), shows that a random linear map works and reduces the whole statement to a concentration bound for the squared length of a projected vector. This is the distributional form of the lemma, and it is what makes the map data-oblivious: the same random projection works for every input, so it can be applied in a single streaming pass.

Setup. Let \( R \in \R^{k \times d} \) have i.i.d. \( \N(0,1) \) entries and define the map \( f(x) = \tfrac{1}{\sqrt{k}} R x \). Fix a vector \( x \) and, without loss of generality, take it to be a unit vector \( u \) (the map is linear, so distortion is scale-invariant). Each coordinate of \( Ru \) is \( (Ru)_i = \sum_{j=1}^d R_{ij} u_j \), a linear combination of independent standard Gaussians with coefficients \( u_j \). A weighted sum of independent Gaussians is Gaussian with variance equal to the sum of squared coefficients, so \( (Ru)_i \sim \N\!\big(0, \sum_j u_j^2\big) = \N(0,1) \) because \( u \) is a unit vector. The \( k \) coordinates are independent because the rows of \( R \) are independent. Therefore

The whole lemma now rests on how tightly a chi-squared variable concentrates around its mean.

Chernoff bound for the upper tail. The moment generating function of a single squared standard Gaussian is \( \E[e^{sZ^2}] = (1 - 2s)^{-1/2} \) for \( s \lt 1/2 \), obtained by completing the square in the Gaussian integral. By independence \( \E[e^{sQ}] = (1-2s)^{-k/2} \). Markov's inequality applied to \( e^{sQ} \) gives, for the event \( Q \ge (1+\varepsilon)k \),

Minimizing the exponent over \( s \) sets the derivative \( \tfrac{k}{1-2s} - (1+\varepsilon)k = 0 \), so \( 1 - 2s = 1/(1+\varepsilon) \) and \( s = \varepsilon / \big(2(1+\varepsilon)\big) \). Substituting,

Using the series bound \( \ln(1+\varepsilon) \le \varepsilon - \tfrac{\varepsilon^2}{2} + \tfrac{\varepsilon^3}{3} \) for \( \varepsilon \gt 0 \), the exponent is at most \( \tfrac{k}{2}\big(-\tfrac{\varepsilon^2}{2} + \tfrac{\varepsilon^3}{3}\big) = -k\big(\tfrac{\varepsilon^2}{4} - \tfrac{\varepsilon^3}{6}\big) \). The lower tail \( \P[Q \le (1-\varepsilon)k] \) is bounded the same way and yields an exponent of the same leading order. Combining both tails,

Union bound and the dimension. There are \( \binom{n}{2} \lt n^2/2 \) pairs of points. Applying the map to the difference \( x_a - x_b \) of each pair (again a fixed vector, so the per-vector bound applies) and taking a union bound over all pairs, the probability that any pairwise distance is distorted by more than \( \varepsilon \) is at most

This is below \( 1 \), so a good projection exists, as soon as \( k\big(\tfrac{\varepsilon^2}{4} - \tfrac{\varepsilon^3}{6}\big) \gt 2\ln n \), that is

Two things are worth naming. First, the bound is on the number of points, not the ambient dimension \( d \): a million points in a million dimensions embed into the same target dimension as a million points in a thousand dimensions. Second, the constant matters in practice. For \( n = 500 \) points and \( \varepsilon = 0.1 \) the bound demands \( k \ge 4\ln 500 / (0.005 - 0.000333) \approx 5327 \), which can exceed the ambient dimension: JL only compresses when \( n \) is large or \( \varepsilon \) is loose. The demonstration below runs it in a regime where it genuinely compresses.

Achlioptas (2003) observed that the Gaussian entries are unnecessary. Entries drawn as \( \pm 1 \) with probability \( 1/2 \), or as \( \{-1, 0, +1\} \) with probabilities \( \{1/6, 2/3, 1/6\} \) and scaled by \( \sqrt{3} \), satisfy the same concentration and are far cheaper: the sparse version touches only a third of the coordinates and needs no floating-point multiplies. This is the database-friendly projection, and it is the reason random projection is practical in a streaming setting.

Random projection for fast approximate algorithms

Once distances survive a projection, any algorithm whose output depends only on pairwise distances can be run in the projected space at a fraction of the cost. Approximate nearest-neighbour search over \( n \) points in \( d \) dimensions costs \( O(nd) \) per query by brute force; projecting to \( k = O(\log n / \varepsilon^2) \) dimensions first reduces the per-query work to \( O(nk) \) with a \( (1 + \varepsilon) \) guarantee on the reported distances. \( k \)-means, agglomerative clustering, and any kernel method that uses only inner products inherit the same speedup, because inner products are recovered from distances by the polarization identity \( \langle x, y\rangle = \tfrac12(\|x\|^2 + \|y\|^2 - \|x-y\|^2) \), and all three quantities are preserved. The projection also composes with locality-sensitive hashing: because random-projection LSH families hash on signs of random linear functionals, projecting first and hashing second is exactly the standard construction, covered in mining massive datasets, which owns the LSH collision-probability analysis. This page owns the reason the projection preserves what LSH then hashes.

Sketching for linear algebra

Random projection preserves the geometry of a fixed set of points. Sketching for linear algebra asks for something stronger: a single random matrix \( S \in \R^{t \times d} \) that preserves the geometry of an entire subspace at once, so that \( \|SAx\| \approx \|Ax\| \) for every vector \( x \) simultaneously. Such an \( S \) is a subspace embedding, and it is the object that makes least squares and low-rank approximation fast.

Subspace embedding. For an \( n \)-column matrix \( A \) with column space of dimension \( r \), a matrix \( S \) is an \( \varepsilon \)-subspace embedding for \( A \) if

Equivalently, writing \( U \) for an orthonormal basis of the column space, \( S \) must satisfy \( \|(SU)^\top(SU) - I\| \le \varepsilon \): \( SU \) is nearly orthonormal. A Gaussian \( S \) with \( t = O(r/\varepsilon^2) \) rows achieves this, by a net argument that applies the distributional JL bound to an \( \varepsilon \)-net of the unit sphere in the column space rather than to a finite point set. The cost of forming \( SA \) with a dense Gaussian is \( O(tnd) \), which defeats the purpose.

CountSketch, the sparse embedding. Clarkson and Woodruff (2013) removed the cost. Their sketch has exactly one nonzero per column: pick a random hash \( h : [d] \to [t] \) and a random sign \( \sigma : [d] \to \{-1, +1\} \), and set \( S_{h(j),\,j} = \sigma(j) \) with the rest of column \( j \) zero. Then \( SA \) can be formed in a single pass over the nonzeros of \( A \), in time proportional to the number of nonzeros, \( O(\mathrm{nnz}(A)) \), which is the input-sparsity time the title of their paper advertises. The map is exactly the Count-Min-style hashing of mining massive datasets applied to columns, with random signs added so that the estimator is unbiased. The sketch dimension needed for a subspace embedding is \( t = O(r^2/\varepsilon^2) \), later improved, larger than the Gaussian's \( O(r/\varepsilon^2) \) but formed vastly faster.

Sketch and solve for least squares. Consider overdetermined least squares, \( \min_x \|Ax - b\| \) with \( A \in \R^{d \times n} \) and \( d \gg n \). The exact solution costs \( O(dn^2) \). Draw a subspace embedding \( S \) for the augmented matrix \( [A \; b] \) and solve the small problem \( \min_x \|S(Ax - b)\| = \min_x \|SAx - Sb\| \), which costs \( O(\mathrm{nnz}(A)) + O(tn^2) \) with \( t = O(n/\varepsilon^2) \). The guarantee is that the sketched solution \( \hat x \) is nearly optimal in objective value,

The proof is one line given the embedding: because \( S \) preserves norms on the span of \( A \)'s columns and \( b \), it preserves the value of every candidate objective up to \( 1 \pm \varepsilon \), so the minimizer of the sketched problem cannot beat the true minimizer by more than that factor. The measured demo below reaches \( \|A\hat x - b\| / \min_x\|Ax-b\| = 1.024 \), inside the guarantee. This is the algorithmic payoff of the JL machinery: least squares in input-sparsity time. Woodruff's 2014 monograph is the reference for the full family of these results, and the exact factorizations they approximate, together with their conditioning, live in numerical methods.

The SVD and low-rank approximation, recalled

Every real matrix \( A \in \R^{m \times n} \) factors as \( A = U\Sigma V^\top \) with \( U, V \) orthonormal and \( \Sigma = \diag(\sigma_1 \ge \sigma_2 \ge \cdots \ge 0) \). The exact construction, its numerical stability, and the conditioning story are developed in numerical methods; this page uses only the approximation theorem it implies. The Eckart-Young theorem states that the best rank-\( k \) approximation in both the spectral and Frobenius norms is the truncation \( A_k = \sum_{i=1}^k \sigma_i u_i v_i^\top \), and the error it leaves is exactly the tail of the spectrum,

This is the target that randomized SVD chases: not merely a rank-\( k \) factorization, but one whose error is close to the optimal \( \sigma_{k+1} \). The reason low-rank approximation matters is that most matrices that arise from data have rapidly decaying spectra, so a small \( k \) captures nearly all the energy, and the truncation both compresses the matrix and denoises it by discarding the directions that carry little signal.

Randomized SVD

Computing a full SVD to then throw away all but the top \( k \) directions is wasteful when \( k \ll \min(m,n) \). The randomized method of Halko, Martinsson, and Tropp (2011) finds an approximate basis for the top-\( k \) subspace by sampling the range of \( A \) with a few random probes, and its cost is \( O(mnk) \) rather than \( O(mn\min(m,n)) \). The derivation is a two-stage decomposition: first find an orthonormal \( Q \) whose columns approximately span the top left-singular subspace, then compute an exact SVD of the small matrix \( Q^\top A \).

Stage one, range finding. Draw a Gaussian test matrix \( \Omega \in \R^{n \times (k+p)} \), where \( p \) is a small oversampling parameter, typically \( 5 \) to \( 10 \). Form the sample matrix \( Y = A\Omega \). Each column of \( Y \) is a random linear combination of the columns of \( A \), so it lies in the range of \( A \) and is heavily weighted toward the directions with large singular values, because those directions dominate the product. Orthonormalize with a thin QR, \( Y = QR \), so that \( Q \) has \( k + p \) orthonormal columns. The key claim is that \( Q \) captures the action of \( A \), meaning \( \|A - QQ^\top A\| \) is small.

The error bound. Halko, Martinsson, and Tropp prove that in expectation over the Gaussian \( \Omega \), with \( p \ge 2 \) oversampling,

The first term is a small multiple of the optimal error \( \sigma_{k+1} \); the second involves the whole tail of the spectrum and is the term that hurts when the spectrum decays slowly. Reading the bound tells you exactly when the plain method is accurate: if the spectrum has a sharp gap after \( \sigma_k \), the tail sum is small and the estimate is nearly optimal, but if the tail decays slowly the second term dominates and a single pass through \( A \) is not enough.

Power iteration. The fix is to apply the method to \( (AA^\top)^q A \) instead of \( A \). This matrix has the same singular vectors but singular values \( \sigma_i^{2q+1} \), which stretches the gap between the top-\( k \) directions and the tail by the power \( 2q+1 \). The sampled range then aligns far more sharply with the true top subspace, and the tail term in the bound is suppressed geometrically. Two steps of power iteration, \( q = 2 \), are enough in practice to make the estimate visually indistinguishable from the optimal truncation, at the cost of two extra matrix multiplies. Each power step must be reorthonormalized in finite precision, or the small singular directions are lost to rounding.

Stage two. With \( Q \) in hand, form the small matrix \( B = Q^\top A \) of size \( (k+p) \times n \), compute its exact SVD \( B = \tilde U \Sigma V^\top \), and set \( U = Q\tilde U \). Then \( A \approx Q B = (Q\tilde U)\Sigma V^\top = U\Sigma V^\top \), truncated to \( k \) columns. Only the small factorization of \( B \) touches the expensive dense SVD routine, and \( B \) has only \( k + p \) rows.

Measured on this machine. Run on a \( 2000 \times 1000 \) matrix with a geometric spectrum \( \sigma_i = 0.95^i \), all in float64, seeking a rank-\( 20 \) approximation with oversampling \( p = 10 \) and \( q = 2 \) power steps. Averaged over three trials on this host's CPU, the full numpy.linalg.svd took about \( 300 \) ms and the randomized method about \( 16.5 \) ms, a \( \mathbf{18\times} \) speedup. The optimal rank-\( 20 \) spectral error is \( \sigma_{21} = 0.3585 \); the randomized method achieved \( 0.3585 \), a ratio of \( 1.000 \), and recovered the top twenty singular values to a maximum absolute error of \( 3.6 \times 10^{-4} \). The exact and randomized factorizations were cross-checked by reconstructing the matrix and confirming the residual, per the numerical-hazard note for this host; the factorizations agree to well within the accuracy the bound promises.

Spectral graph theory

A graph carries a linear operator, and its spectrum encodes combinatorial structure that is otherwise hard to compute. Let \( G \) be an undirected graph on \( n \) vertices with adjacency matrix \( A \) and degree matrix \( D = \diag(d_1, \ldots, d_n) \). The (unnormalized) graph Laplacian is

Its defining property comes from the quadratic form. For any vector \( x \in \R^n \),

The expansion follows by writing each edge's contribution twice and collecting terms. Because it is a sum of squares, \( L \) is positive semidefinite, so all its eigenvalues are nonnegative. The constant vector \( \mathbf{1} \) gives \( L\mathbf{1} = 0 \), so \( \lambda_1 = 0 \) always, and the multiplicity of the zero eigenvalue equals the number of connected components: on each component the only vectors with zero quadratic form are the constants, since \( \sum_{(i,j)\in E}(x_i - x_j)^2 = 0 \) forces \( x \) to be constant across every edge. The smallest nonzero eigenvalue \( \lambda_2 \), the algebraic connectivity named by Fiedler (1973), measures how close the graph is to being disconnected: it is zero exactly when the graph splits, and small when the graph has a sparse cut.

The Cheeger inequality. Make the connection to cuts precise. The conductance of a vertex set \( S \), with \( \bar S \) its complement and \( \mathrm{vol}(S) = \sum_{i \in S} d_i \), is

where \( E(S,\bar S) \) is the set of edges crossing the cut. Finding the minimizer is NP-hard in general. Working with the normalized Laplacian \( \mathcal{L} = D^{-1/2} L D^{-1/2} \) and its second eigenvalue \( \lambda_2(\mathcal{L}) \), the Cheeger inequality (Cheeger 1970, in its graph form) states

The two-sided bound sandwiches an NP-hard combinatorial quantity between multiples of an eigenvalue that any linear-algebra library computes in polynomial time. The intuition is that a sparse cut and a small \( \lambda_2 \) are the same phenomenon seen two ways: a sparse cut means the graph almost falls into two pieces, which means there is a near-constant-on-each-side vector with a tiny quadratic form, which means \( \lambda_2 \) is small.

The easy direction. The lower bound \( \lambda_2(\mathcal{L})/2 \le \phi(G) \), equivalently that a sparse cut forces a small eigenvalue, is proved by exhibiting a test vector. Let \( S \) be the optimal cut and take the vector \( x \) that is \( +1/\mathrm{vol}(S) \) on \( S \) and \( -1/\mathrm{vol}(\bar S) \) on \( \bar S \), a scaled indicator arranged to be orthogonal to the degree vector so it is a legal test vector for \( \lambda_2 \). Its Rayleigh quotient is \( x^\top L x / x^\top D x = |E(S,\bar S)|\,(1/\mathrm{vol}(S) + 1/\mathrm{vol}(\bar S)) / (1/\mathrm{vol}(S) + 1/\mathrm{vol}(\bar S)) = |E(S,\bar S)| \cdot \big(\tfrac{1}{\mathrm{vol}(S)} + \tfrac{1}{\mathrm{vol}(\bar S)}\big) \) divided out, which simplifies to at most \( 2\phi(S) \). Because \( \lambda_2 \) is the minimum Rayleigh quotient over all vectors orthogonal to the trivial eigenvector, it is at most this test value, giving \( \lambda_2 \le 2\phi(G) \). The hard direction, that a small \( \lambda_2 \) yields an actual sparse cut, is the substance of Cheeger's theorem: it is constructive, sweeping a threshold across the second eigenvector and proving one of the resulting cuts has conductance \( O(\sqrt{\lambda_2}) \), and it is exactly the rounding step that makes spectral partitioning an algorithm rather than a bound. Spielman and Teng's work on nearly-linear-time Laplacian solvers turned this line into a practical toolkit.

Spectral clustering, derived. The Cheeger rounding suggests the algorithm: to bipartition a graph, compute the second eigenvector and split on its entries. Here is why that eigenvector is the right object. The balanced minimum cut can be written as the RatioCut objective, \( \mathrm{RatioCut}(S, \bar S) = |E(S,\bar S)|\big(\tfrac{1}{|S|} + \tfrac{1}{|\bar S|}\big) \). Encode the partition by the vector

chosen so that \( f \perp \mathbf{1} \) and \( \|f\|^2 = n \). A short computation using \( f^\top L f = \sum_{(i,j)\in E}(f_i - f_j)^2 \) shows that only cut edges contribute, and \( f^\top L f = |V| \cdot \mathrm{RatioCut}(S, \bar S) \). Minimizing RatioCut over partitions is therefore minimizing \( f^\top L f \) over vectors of this discrete form, which is NP-hard. Relaxing \( f \) to range over all real vectors with \( f \perp \mathbf{1} \) and \( \|f\| = \sqrt n \) turns it into \( \min_{f \perp \mathbf 1} f^\top L f / f^\top f \), whose solution, by the Courant-Fischer characterization, is exactly the eigenvector of the second-smallest eigenvalue, the Fiedler vector. Rounding the real solution back to a partition by thresholding its entries recovers a discrete cut. Von Luxburg's 2007 tutorial is the definitive account, including the extension to \( k \) clusters by embedding vertices into the space of the bottom \( k \) eigenvectors and running \( k \)-means there.

PageRank as an eigenvector

PageRank scores a directed graph by the stationary distribution of a random surfer who follows a random outgoing link with probability \( \beta \) and teleports to a uniformly random vertex with probability \( 1 - \beta \). The score vector \( \pi \) is the dominant eigenvector of the Google matrix \( M = \beta P + (1-\beta)\tfrac{1}{n}\mathbf{1}\mathbf{1}^\top \), where \( P \) is the column-stochastic link matrix, and it satisfies \( M\pi = \pi \). The construction, sink handling, and the interpretation as a Markov chain are developed in mining massive datasets. What this page adds is the convergence rate of the power iteration that computes it, because that rate is a clean consequence of the spectral gap and ties PageRank to the rest of the toolbox.

Power iteration and its rate. Start from any distribution \( x_0 \) and iterate \( x_{t+1} = M x_t \). Expand \( x_0 \) in the eigenbasis of \( M \), \( x_0 = \pi + \sum_{i \ge 2} c_i v_i \), where \( \pi = v_1 \) has eigenvalue \( 1 \) and the remaining eigenvalues satisfy \( |\lambda_i| \le |\lambda_2| \). Applying \( M \) t times multiplies each component by its eigenvalue to the \( t \)-th power,

The error decays geometrically at rate \( |\lambda_2| \), the second-largest eigenvalue in magnitude. The teleportation is what controls this rate: it is a theorem that for the Google matrix the second eigenvalue satisfies \( |\lambda_2| \le \beta \), independent of the link structure. With the usual \( \beta = 0.85 \), the error shrinks by \( 0.85 \) per iteration, so reaching a relative accuracy of \( 10^{-6} \) needs \( t \ge \ln(10^{-6}) / \ln(0.85) \approx 13.8 / 0.1625 \approx 85 \) iterations. This is why PageRank is cheap: a fixed number of sparse matrix-vector products, set by the gap \( 1 - \beta \), regardless of graph size. The same spectral-gap logic governs mixing times of Markov chains and the convergence of the power method inside randomized SVD.

Linear programming, relaxation, and rounding

A linear program optimizes a linear objective over linear constraints, \( \min c^\top x \) subject to \( Ax \ge b \), \( x \ge 0 \). Its geometry, the simplex and interior-point methods, and the duality theorem that pairs it with a maximization problem are developed in combinatorial optimization, which owns those proofs. The algorithmic view relevant here is different: many NP-hard combinatorial problems are integer programs, and their LP relaxation, obtained by dropping the integrality constraint \( x \in \{0,1\} \) to \( x \in [0,1] \), is solvable in polynomial time and lower-bounds the integer optimum. Rounding the fractional LP solution back to an integer one, while controlling how much the objective grows, is the source of many of the best approximation algorithms. Williamson and Shmoys' book is the definitive treatment.

Vertex cover by LP rounding. A vertex cover is a set of vertices touching every edge; the minimum one is NP-hard. Write it as an integer program with a variable \( x_v \in \{0,1\} \) per vertex, minimize \( \sum_v x_v \) subject to \( x_u + x_v \ge 1 \) for every edge \( (u,v) \). Relax to \( x_v \in [0,1] \) and solve the LP, obtaining a fractional \( x^* \). Round by the threshold rule: put \( v \) in the cover if and only if \( x^*_v \ge 1/2 \).

The rounded set is a valid cover, because every edge constraint \( x^*_u + x^*_v \ge 1 \) forces at least one of the two endpoints to be at least \( 1/2 \). The cost is controlled because each vertex we keep had \( x^*_v \ge 1/2 \), so \( \mathbf{1}[x^*_v \ge 1/2] \le 2x^*_v \), and summing,

The final inequality holds because the relaxation only enlarges the feasible set, so its optimum is no larger than the integer optimum. This is a factor-\( 2 \) approximation, and it is essentially the best known for vertex cover.

Randomized rounding

When the threshold rule is too crude, treat the fractional LP solution as a probability and flip a coin. For set cover, this randomized rounding gives an \( O(\log n) \) approximation, which is optimal up to constants. The problem: a universe of \( n \) elements, a family of sets \( S_1, \ldots, S_m \) with costs \( c_j \), find a minimum-cost subfamily whose union is the universe. The LP has a variable \( x_j \in [0,1] \) per set, minimizes \( \sum_j c_j x_j \), subject to \( \sum_{j : e \in S_j} x_j \ge 1 \) for each element \( e \).

The rounding. Solve the LP for \( x^* \). Perform \( T = 2\ln n \) independent rounds; in each round, include set \( S_j \) independently with probability \( x^*_j \). The final cover is the union over all rounds of the chosen sets.

Cost, by linearity of expectation. In a single round the expected cost is \( \sum_j c_j x^*_j = \mathrm{OPT}_{\mathrm{LP}} \). Over \( T \) rounds, by linearity, the expected total cost is \( T \cdot \mathrm{OPT}_{\mathrm{LP}} = 2\ln n \cdot \mathrm{OPT}_{\mathrm{LP}} \le 2\ln n \cdot \mathrm{OPT}_{\mathrm{IP}} \).

Coverage, by the constraint. Fix an element \( e \) and consider one round. The probability \( e \) is not covered is \( \prod_{j : e \in S_j}(1 - x^*_j) \). Using \( 1 - x \le e^{-x} \), this is at most \( \exp\!\big(-\sum_{j : e \in S_j} x^*_j\big) \le e^{-1} \), because the LP constraint guarantees the exponent is at least \( 1 \). The rounds are independent, so after \( T = 2\ln n \) rounds the probability \( e \) is still uncovered is at most \( (e^{-1})^{2\ln n} = e^{-2\ln n} = n^{-2} \). A union bound over all \( n \) elements bounds the probability that anything is left uncovered by \( n \cdot n^{-2} = 1/n \). So with probability at least \( 1 - 1/n \) every element is covered, at expected cost \( O(\log n) \cdot \mathrm{OPT} \). The \( \ln n \) factor is not an artifact: set cover is hard to approximate better than \( (1-o(1))\ln n \) unless P = NP, so randomized rounding is optimal. The same idea, rounding a semidefinite rather than linear relaxation, gives the \( 0.878 \) MAX-CUT algorithm of Goemans and Williamson.

The multiplicative weights update method

One meta-algorithm unifies boosting, fast approximate LP solving, and equilibrium computation in games. The setup is online decision-making: there are \( N \) experts; on each of \( T \) rounds the algorithm picks a distribution \( p^t \) over experts, then a loss vector \( \ell^t \in [0,1]^N \) is revealed and the algorithm pays \( p^t \cdot \ell^t \). The goal is to compete with the single best expert in hindsight, and the quantity to minimize is the regret,

The algorithm. Keep a weight \( w^t_i \) per expert, initialized to \( 1 \). Play \( p^t_i = w^t_i / \sum_j w^t_j \). After seeing \( \ell^t \), multiply each weight down in proportion to its loss, \( w^{t+1}_i = w^t_i\,e^{-\eta \ell^t_i} \), for a learning rate \( \eta \gt 0 \). Experts that perform well retain their weight; poor experts are exponentially demoted.

The regret bound, derived via a potential. Let \( \Phi_t = \sum_i w^t_i \) be the total weight. Track how it changes: \[ \Phi_{t+1} = \sum_i w^t_i e^{-\eta \ell^t_i} \le \sum_i w^t_i (1 - \eta \ell^t_i + \eta^2 (\ell^t_i)^2) \le \Phi_t\big(1 - \eta\, p^t\cdot\ell^t + \eta^2\big), \] using \( e^{-x} \le 1 - x + x^2 \) for \( x \ge 0 \), the definition of \( p^t \), and \( \ell^t_i \le 1 \). Then \( 1 + y \le e^y \) gives \( \Phi_{t+1} \le \Phi_t \exp(-\eta\, p^t\cdot\ell^t + \eta^2) \), and telescoping from \( \Phi_1 = N \), \[ \Phi_{T+1} \le N \exp\!\Big(-\eta \sum_t p^t\cdot\ell^t + \eta^2 T\Big). \] On the other hand, the total weight is at least the weight of any single expert \( i \): \( \Phi_{T+1} \ge w^{T+1}_i = \exp(-\eta \sum_t \ell^t_i) \). Taking the best expert \( i^\star \), combining the two bounds, taking logs, and rearranging, \[ \sum_t p^t\cdot\ell^t - \sum_t \ell^t_{i^\star} \le \frac{\ln N}{\eta} + \eta T. \] Optimizing the rate at \( \eta = \sqrt{\ln N / T} \) gives

The average regret per round is \( O(\sqrt{\ln N / T}) \to 0 \): the algorithm learns to match the best expert without knowing in advance which one it is. Measured on a synthetic instance with \( N = 10 \) experts and \( T = 10{,}000 \) rounds, with losses in \( [0,1] \) and one expert made clearly better, the update achieved cumulative regret \( 152 \) against a bound of \( 2\sqrt{T\ln N} \approx 303 \), inside the guarantee, with average regret \( 0.0152 \) per round.

Why it unifies three subjects. In boosting, the experts are training examples and the losses encode whether the current weak learner classifies each one correctly; multiplicative weights on the examples is exactly AdaBoost (Freund and Schapire 1997), and the regret bound becomes the boosting error guarantee. In LP solving, the experts are the constraints of a feasibility LP; running the update to drive the weighted-average constraint violation to zero gives a fast approximate solver whose iteration count is \( O(\log m / \varepsilon^2) \), the basis of the width-based LP algorithms surveyed by Arora, Hazan, and Kale (2012). In zero-sum games, two players each running the update converge to the minimax equilibrium, giving an algorithmic proof of von Neumann's minimax theorem. One potential-function argument, three fields.

Compressed sensing at a glance

Compressed sensing recovers a sparse signal from far fewer linear measurements than its dimension. If \( x \in \R^d \) is \( s \)-sparse (at most \( s \) nonzeros) and we observe \( y = Ax \) with \( A \in \R^{m \times d} \) and \( m \ll d \), the system is underdetermined and has infinitely many solutions, yet the sparse one is recoverable when \( A \) is well-behaved. The condition, due to Candes and Tao (2005), is the restricted isometry property: \( A \) satisfies RIP of order \( s \) with constant \( \delta_s \) if for every \( s \)-sparse \( x \),

RIP says every submatrix of \( s \) columns is nearly an isometry, so no sparse signal is nearly annihilated and distinct sparse signals stay distinguishable. It is the subspace-embedding condition from the sketching section, restricted to sparse supports rather than a fixed subspace, which is why compressed sensing belongs in this toolbox. Random Gaussian and subgaussian matrices satisfy RIP of order \( 2s \) with high probability once \( m = O(s \log(d/s)) \), the same JL-flavored counting: union-bound the isometry over all \( \binom{d}{s} \) supports.

Given RIP, the recovery is convex. The natural problem, minimize \( \|x\|_0 \) (the count of nonzeros) subject to \( Ax = y \), is NP-hard, but its convex relaxation replaces the \( \ell_0 \) count with the \( \ell_1 \) norm,

a linear program. Candes and Tao proved that if \( A \) has RIP of order \( 2s \) with \( \delta_{2s} \lt \sqrt 2 - 1 \), then \( \ell_1 \) minimization recovers the true \( s \)-sparse \( x \) exactly. The geometric reason is that the \( \ell_1 \) ball is pointed along the sparse axes, so its facets touch the affine solution set \( \{x : Ax = y\} \) exactly at sparse vectors; the RIP guarantees the touch happens at the true one. The full recovery proof is beyond this page and lives in the cited work; what belongs here is that \( \ell_1 \) minimization is the same relaxation-and-geometry move as LP rounding, applied to signals instead of set systems.

Worked problems

Two of the required worked problems are collected here; three more are embedded in the theory above (Problems 1 through 5) where they land closest to the derivation they exercise.

Implementation

The three central algorithms of the page are short enough to write in full. All of them run in float64 and were cross-checked against exact references on this host, per the numerical-hazard note. The first pair implements the JL projection and measures the distortion it actually achieves; the NumPy version is the reference and the JAX version is the vectorized form that runs on an accelerator when the point set is large.

import numpy as np

def jl_project(X, k, seed=0):
    # X: (n, d) points -> Y: (n, k) projected points
    rng = np.random.default_rng(seed)
    d = X.shape[1]
    R = rng.standard_normal((d, k)) / np.sqrt(k)   # (d, k) Gaussian, scaled
    return X @ R                                    # (n, k)

def measured_distortion(X, Y, n_pairs=5000, seed=1):
    rng = np.random.default_rng(seed)
    n = X.shape[0]
    a = rng.integers(0, n, n_pairs)
    b = rng.integers(0, n, n_pairs)
    m = a != b                                      # drop self-pairs
    a, b = a[m], b[m]
    d_orig = np.linalg.norm(X[a] - X[b], axis=1)    # (P,)
    d_proj = np.linalg.norm(Y[a] - Y[b], axis=1)    # (P,)
    ratio = d_proj / d_orig
    return ratio.mean(), np.abs(ratio - 1).max()

rng = np.random.default_rng(0)
n, d, eps = 500, 10000, 0.3
X = rng.standard_normal((n, d))
k = int(np.ceil(4 * np.log(n) / (eps**2 / 2 - eps**3 / 3)))   # JL bound -> 691
Y = jl_project(X, k)
mean_ratio, max_dist = measured_distortion(X, Y)
print(k, mean_ratio, max_dist)
# 691   ~0.9985   ~0.10   (compression 10000/691 = 14.5x, max distortion well under eps=0.3)

import jax, jax.numpy as jnp

def jl_project(X, k, key):
    # X: (n, d) -> (n, k); runs on GPU/TPU when X is on device
    d = X.shape[1]
    R = jax.random.normal(key, (d, k)) / jnp.sqrt(k)   # (d, k)
    return X @ R                                        # (n, k)

@jax.jit
def distortion(X, Y, a, b):
    d_orig = jnp.linalg.norm(X[a] - X[b], axis=1)       # (P,)
    d_proj = jnp.linalg.norm(Y[a] - Y[b], axis=1)       # (P,)
    ratio = d_proj / d_orig
    return ratio.mean(), jnp.abs(ratio - 1).max()

key = jax.random.PRNGKey(0)
k1, k2, k3 = jax.random.split(key, 3)
n, d, eps = 500, 10000, 0.3
X = jax.random.normal(k1, (n, d))
k = int(jnp.ceil(4 * jnp.log(n) / (eps**2 / 2 - eps**3 / 3)))   # 691
Y = jl_project(X, k, k2)
a = jax.random.randint(k3, (5000,), 0, n)
b = jax.random.randint(jax.random.fold_in(k3, 1), (5000,), 0, n)
mask = a != b
print(k, distortion(X, Y, a[mask], b[mask]))

The randomized SVD is next. The NumPy version is the reference that produced the measured timings above; the PyTorch version moves the two heavy matrix multiplies and the small dense SVD onto the GPU through cuSOLVER, which is unaffected by the CPU BLAS issue and is where a production randomized SVD would run.

import numpy as np

def randomized_svd(A, k, p=10, q=2, seed=0):
    # A: (m, n). Returns rank-k factors U (m,k), s (k,), Vt (k,n).
    rng = np.random.default_rng(seed)
    m, n = A.shape
    Omega = rng.standard_normal((n, k + p))      # (n, k+p) test matrix
    Y = A @ Omega                                # (m, k+p) sample of range(A)
    for _ in range(q):                           # power iteration sharpens the gap
        Y = A @ (A.T @ Y)                        # reorthonormalize in practice
    Q, _ = np.linalg.qr(Y)                       # (m, k+p) orthonormal basis
    B = Q.T @ A                                  # (k+p, n) small matrix
    Ub, s, Vt = np.linalg.svd(B, full_matrices=False)
    U = Q @ Ub                                   # lift back to (m, k+p)
    return U[:, :k], s[:k], Vt[:k]

rng = np.random.default_rng(0)
# matrix with a geometric spectrum sigma_i = 0.95**i
U0, _ = np.linalg.qr(rng.standard_normal((2000, 2000)))
V0, _ = np.linalg.qr(rng.standard_normal((1000, 1000)))
sv = 0.95 ** np.arange(1000)
A = (U0[:, :1000] * sv) @ V0.T
Ur, sr, Vtr = randomized_svd(A, k=20)
err = np.linalg.norm(A - (Ur * sr) @ Vtr, 2)     # ~0.3585 = sigma_21 (optimal)
print(err, sv[20])                               # matches optimal to 4 digits

import torch

def randomized_svd(A, k, p=10, q=2):
    # A: (m, n) on cuda. cuSOLVER path is unaffected by the CPU BLAS issue.
    m, n = A.shape
    Omega = torch.randn(n, k + p, device=A.device, dtype=A.dtype)   # (n, k+p)
    Y = A @ Omega                                                   # (m, k+p)
    for _ in range(q):
        Y = A @ (A.T @ Y)
        Y, _ = torch.linalg.qr(Y)                                   # reorthonormalize
    Q, _ = torch.linalg.qr(Y)                                       # (m, k+p)
    B = Q.T @ A                                                     # (k+p, n)
    Ub, s, Vt = torch.linalg.svd(B, full_matrices=False)
    return (Q @ Ub)[:, :k], s[:k], Vt[:k]

torch.manual_seed(0)
dev = "cuda" if torch.cuda.is_available() else "cpu"
A = torch.randn(2000, 1000, device=dev, dtype=torch.float64)
Ur, sr, Vtr = randomized_svd(A, k=20)
err = torch.linalg.norm(A - (Ur * sr) @ Vtr, 2)
print(err.item())

The third is spectral clustering. The NumPy version builds the Laplacian, takes its eigendecomposition, and partitions on the Fiedler vector; it reproduces the two-triangle computation of Problem 4 exactly. The JAX version is the same math with a jitted eigensolver for larger graphs.

import numpy as np

def spectral_bipartition(edges, n):
    A = np.zeros((n, n))
    for i, j in edges:
        A[i, j] = A[j, i] = 1.0                  # symmetric adjacency
    L = np.diag(A.sum(1)) - A                     # unnormalized Laplacian
    w, V = np.linalg.eigh(L)                      # ascending eigenvalues
    fiedler = V[:, 1]                             # second-smallest eigenvector
    part = (fiedler > 0).astype(int)             # threshold at zero
    return w, fiedler, part

# two triangles {0,1,2} and {3,4,5} joined by the bridge (2,3)
edges = [(0, 1), (0, 2), (1, 2), (3, 4), (3, 5), (4, 5), (2, 3)]
w, fiedler, part = spectral_bipartition(edges, 6)
print(np.round(w, 4))        # [0. 0.4384 3. 3. 3. 4.5616]
print(np.round(fiedler, 4))  # [-0.4647 -0.4647 -0.261 0.261 0.4647 0.4647]
print(part)                  # [0 0 0 1 1 1]  -> cuts exactly the bridge edge

import jax, jax.numpy as jnp

@jax.jit
def spectral_bipartition(A):
    # A: (n, n) dense symmetric adjacency
    L = jnp.diag(A.sum(1)) - A                    # Laplacian
    w, V = jnp.linalg.eigh(L)                      # ascending
    fiedler = V[:, 1]
    return w, fiedler, (fiedler > 0).astype(jnp.int32)

edges = [(0, 1), (0, 2), (1, 2), (3, 4), (3, 5), (4, 5), (2, 3)]
A = jnp.zeros((6, 6))
for i, j in edges:
    A = A.at[i, j].set(1.0).at[j, i].set(1.0)
w, fiedler, part = spectral_bipartition(A)
print(w[:2], part)            # lambda_2 ~ 0.4384, partition splits the triangles

Finally, the multiplicative weights update, which is short enough to state once and reuse across boosting, LP solving, and games. This is the loop whose measured regret \( 152 \) sat inside the bound \( 303 \) above.

import numpy as np

def multiplicative_weights(losses, eta=None):
    # losses: (T, N) in [0,1]. Returns cumulative algorithm loss and regret.
    T, N = losses.shape
    if eta is None:
        eta = np.sqrt(np.log(N) / T)             # optimal rate 2 sqrt(T ln N) regret
    w = np.ones(N)
    alg = 0.0
    for t in range(T):
        p = w / w.sum()                          # play the weighted distribution
        alg += p @ losses[t]                     # pay expected loss
        w *= np.exp(-eta * losses[t])            # demote high-loss experts
    best = losses.sum(0).min()                   # best expert in hindsight
    return alg, alg - best

rng = np.random.default_rng(7)
losses = rng.uniform(0, 1, (10000, 10))
losses[:, 3] *= 0.3                              # make expert 3 clearly best
alg, regret = multiplicative_weights(losses)
print(regret, 2 * np.sqrt(10000 * np.log(10)))  # 152 vs bound 303

How it is done in practice

Production systems rarely call a textbook random projection directly, but they run its descendants everywhere. Approximate nearest-neighbour search at web scale, the retrieval layer under every embedding-based recommender and every retrieval-augmented language model, is built on exactly the projection-plus-quantization pipeline derived here. Faiss, Meta's library, offers random-projection and PCA preprocessing followed by product quantization; the projection step is the JL map, and its target dimension is tuned against the same \( \varepsilon \)-versus-cost tradeoff the bound describes. HNSW and IVF indexes then search the reduced space. The engineering gap between the derivation and the deployment is mostly in the quantizer and the memory layout, not in the projection, which is a single dense or sparse matmul.

Randomized numerical linear algebra has moved from research to the default. The randomized SVD is the recommended path for truncated factorizations in scikit-learn (TruncatedSVD and PCA(svd_solver="randomized") both call an implementation of the Halko-Martinsson-Tropp algorithm), and it is what makes PCA on a matrix with millions of rows tractable. On the GPU the two matmuls dominate and run at the accelerator's peak; the \( 18\times \) CPU speedup measured above understates the win at scale, because the exact SVD's cubic term grows faster than the randomized method's. The oversampling and power-step parameters are the two knobs that matter: \( p = 10 \) and \( q \in \{1, 2\} \) are the near-universal defaults, and increasing \( q \) is the fix whenever the spectrum is flat.

Spectral methods run inside graph analytics and image segmentation, but at scale nobody forms a dense Laplacian. The eigenvector is computed by a Lanczos or LOBPCG iteration through scipy.sparse.linalg.eigsh, which only needs sparse matrix-vector products, and for the very largest graphs the exact eigenvector is replaced by a few steps of a Laplacian-solver-based method in the line of Spielman and Teng, whose nearly-linear-time solvers made spectral partitioning scale to billions of edges. The rounding step, sweeping a threshold across the eigenvector and taking the best conductance cut, is the practical realization of the hard direction of Cheeger's inequality. PageRank at production scale is a handful of sparse matrix-vector products, and the iteration count is set by the teleport parameter through the \( |\lambda_2| \le \beta \) bound, which is why \( 50 \) to \( 100 \) iterations suffice on a graph of any size.

LP-based approximation and multiplicative weights are less visible but structural. Modern LP and MILP solvers use the LP relaxation as the bound inside branch-and-bound, and randomized and deterministic rounding appear in scheduling, routing, and network-design heuristics. Multiplicative weights and its mirror-descent generalizations are the theoretical core of online learning, of the AdaBoost family, and of the no-regret dynamics that solvers for large two-player games run; the same update, under the name Hedge or exponentiated gradient, appears throughout online convex optimization.

The current research frontier

Randomized numerical linear algebra is the most active thread. The 2020 Acta Numerica survey by Martinsson and Tropp consolidated the field, and the current frontier is single-pass and streaming factorizations that see each entry of the matrix once, sketch-based preconditioners such as Blendenpik and LSRN that combine a random embedding with an exact iterative solver to get high-accuracy least squares at sketch speed, and the extension of these ideas to tensor decompositions. Groups at the University of Washington, the University of Texas at Austin, Michigan, and Berkeley, together with the RandNLA program that Drineas and Mahoney have long driven, are central; the practical libraries now ship in LAPACK-adjacent form through the RandLAPACK effort.

On the sketching side the question is how small a sparse embedding can be while remaining a subspace embedding, refining the Clarkson-Woodruff input-sparsity result; the target dimension has been pushed from \( O(r^2/\varepsilon^2) \) toward the Gaussian's \( O(r/\varepsilon^2) \) by work of Nelson, Nguyen, Cohen, and others, with the optimal tradeoff between sparsity and dimension for oblivious subspace embeddings still not fully settled. Fast Johnson-Lindenstrauss transforms that use structured (Hadamard-based) matrices to apply a projection in \( O(d\log d) \) rather than \( O(dk) \) time, in the line of Ailon and Chazelle, are the practical form when the projection must be fast rather than sparse.

Spectral graph theory's frontier is local and near-linear algorithms: local clustering that finds a good cut around a seed vertex without touching the whole graph, spectral sparsification that compresses a dense graph to a sparse one preserving all cut values, and the fast Laplacian solvers that underlie both, a program advanced by Spielman and Teng and continued by Koutis, Peng, Kelner, and others at MIT, CMU, Georgia Tech, and Yale. These ideas now feed graph neural networks, where the Laplacian spectrum defines the convolution, a connection developed in graph machine learning. Compressed sensing has matured into a broad theory of structured recovery, and its RIP-style analysis reappears in the sample complexity of matrix completion and in the theory of why overparameterized networks generalize.

Open source to read

• scikit-learn/scikit-learn: read sklearn/random_projection.py for the Gaussian and sparse JL transforms with the automatic dimension from the bound, and sklearn/utils/extmath.py for randomized_svd, the reference implementation of the Halko-Martinsson-Tropp method used by TruncatedSVD and PCA.
• scipy/scipy: scipy/sparse/linalg/_eigen holds the ARPACK and LOBPCG wrappers (eigsh) that compute Laplacian eigenvectors for spectral clustering without forming a dense matrix; start there for anything graph-spectral at scale.
• facebookresearch/faiss: the vector-search library; open faiss/VectorTransform.h for the random-projection and PCA preprocessors and the IVF/PQ/HNSW indexes that consume the reduced vectors.
• google/jax: for writing the projection and eigensolver code above in a form that compiles to GPU/TPU; jax.numpy.linalg mirrors the NumPy API used here.
• numpy/numpy: numpy/random/_generator.pyx for the modern default_rng generators the projections draw from, and numpy/linalg for the exact SVD and QR that the randomized method calls on its small blocks.
• scikit-learn SpectralClustering: the end-to-end spectral clustering pipeline, from affinity graph to eigenmap to \( k \)-means, matching the derivation in this page's spectral section.

Common misconceptions

"Random projection needs the data to be Gaussian or low-dimensional to work." The Johnson-Lindenstrauss guarantee is distribution-free and dimension-free: it holds for any finite point set in any dimension, because the randomness is in the projection matrix, not in the data. The target dimension depends only on the number of points and the accuracy.

"The randomized SVD is just a faster but less accurate SVD." With a spectral gap it is essentially exact for the top subspace, and even without a gap, two power iterations bring it to within a small constant of the optimal rank-\( k \) error. The measured example above matched the optimal \( \sigma_{21} \) to four digits. It is not a lossy approximation to be tolerated; it is the recommended method for truncated factorizations.

"The second eigenvector gives the minimum cut." It gives the minimum of a continuous relaxation of a balanced-cut objective, which is then rounded. The relaxation can differ from the true minimum cut, and the Cheeger inequality quantifies the gap: the rounded cut has conductance \( O(\sqrt{\lambda_2}) \), not \( \lambda_2 \). The square root is the price of relaxing.

"LP rounding always loses a factor of 2." Factor \( 2 \) is the vertex-cover threshold-rounding bound specifically; randomized rounding gives \( O(\log n) \) for set cover, and the integrality gap of the relaxation, not the rounding rule, is what ultimately limits the ratio. A different relaxation, such as the semidefinite one for MAX-CUT, yields entirely different constants.

"Multiplicative weights is just gradient descent." It is mirror descent with the entropy mirror map, which is a different geometry from Euclidean gradient descent: its regret scales with \( \sqrt{\log N} \) rather than \( \sqrt{N} \), which is exponentially better when the number of experts is large. That logarithmic dependence is exactly why it works with an exponentially large expert set, as in LP solving.

"Compressed sensing recovers any signal from few measurements." It recovers sparse signals, and only when the measurement matrix satisfies the restricted isometry property. A dense signal in generic coordinates is not recoverable from a compressed measurement; the whole theory rests on the prior that the signal is sparse in some known basis.

"A subspace embedding is the same as preserving each data point's norm." Preserving the norms of a fixed finite set of points is the JL guarantee. A subspace embedding preserves norms for every vector in a subspace simultaneously, an uncountable set, which is strictly stronger and needs the net argument rather than a union bound over points.

Self-check

References

1. Motwani, R. and Raghavan, P. Randomized Algorithms. Cambridge University Press, 1995.
2. Mitzenmacher, M. and Upfal, E. Probability and Computing: Randomization and Probabilistic Techniques in Algorithms and Data Analysis. 2nd ed., Cambridge University Press, 2017.
3. Blum, A., Hopcroft, J., and Kannan, R. Foundations of Data Science. Cambridge University Press, 2020. Freely available at cs.cornell.edu/jeh/book.pdf.
4. Vershynin, R. High-Dimensional Probability: An Introduction with Applications in Data Science. Cambridge University Press, 2018.
5. Williamson, D. P. and Shmoys, D. B. The Design of Approximation Algorithms. Cambridge University Press, 2011.
6. Woodruff, D. P. Sketching as a tool for numerical linear algebra. Foundations and Trends in Theoretical Computer Science, 10(1-2):1-157, 2014. arXiv:1411.4357.
7. Johnson, W. B. and Lindenstrauss, J. Extensions of Lipschitz mappings into a Hilbert space. Contemporary Mathematics, 26:189-206, 1984.
8. Dasgupta, S. and Gupta, A. An elementary proof of a theorem of Johnson and Lindenstrauss. Random Structures & Algorithms, 22(1):60-65, 2003.
9. Achlioptas, D. Database-friendly random projections: Johnson-Lindenstrauss with binary coins. Journal of Computer and System Sciences, 66(4):671-687, 2003.
10. Halko, N., Martinsson, P.-G., and Tropp, J. A. Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review, 53(2):217-288, 2011. arXiv:0909.4061.
11. Martinsson, P.-G. and Tropp, J. A. Randomized numerical linear algebra: foundations and algorithms. Acta Numerica, 29:403-572, 2020. arXiv:2002.01387.
12. Clarkson, K. L. and Woodruff, D. P. Low rank approximation and regression in input sparsity time. Proc. STOC, 2013. arXiv:1207.6365.
13. Sarlós, T. Improved approximation algorithms for large matrices via random projections. Proc. FOCS, 2006.
14. Drineas, P. and Mahoney, M. W. RandNLA: randomized numerical linear algebra. Communications of the ACM, 59(6):80-90, 2016.
15. Cheeger, J. A lower bound for the smallest eigenvalue of the Laplacian. In Problems in Analysis, Princeton University Press, 195-199, 1970.
16. Fiedler, M. Algebraic connectivity of graphs. Czechoslovak Mathematical Journal, 23(2):298-305, 1973.
17. Spielman, D. A. and Teng, S.-H. Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems. Proc. STOC, 2004; journal versions in SIAM J. Computing, 2011-2014.
18. von Luxburg, U. A tutorial on spectral clustering. Statistics and Computing, 17(4):395-416, 2007. arXiv:0711.0189.
19. Arora, S., Hazan, E., and Kale, S. The multiplicative weights update method: a meta-algorithm and applications. Theory of Computing, 8(1):121-164, 2012.
20. Freund, Y. and Schapire, R. E. A decision-theoretic generalization of on-line learning and an application to boosting. Journal of Computer and System Sciences, 55(1):119-139, 1997.
21. Candès, E. J. and Tao, T. Decoding by linear programming. IEEE Transactions on Information Theory, 51(12):4203-4215, 2005.
22. Candès, E. J., Romberg, J., and Tao, T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information. IEEE Transactions on Information Theory, 52(2):489-509, 2006.
23. Nemhauser, G. L. and Trotter, L. E. Vertex packings: structural properties and algorithms. Mathematical Programming, 8:232-248, 1975.
24. Page, L., Brin, S., Motwani, R., and Winograd, T. The PageRank citation ranking: bringing order to the web. Stanford Digital Library Technologies Project, technical report, 1999.