Fast Length-Squared Sampling for Positive-Semidefinite Matrices

Abstract: We describe a simple rejection-sampling-based algorithm to perform length-squared sampling on an $n \times n$ positive-semidefinite (psd) matrix: that is, to sample a column with probability proportional to its squared $\ell_2$-norm. The algorithm runs in just $O(n)$ expected time, which is significantly sublinear in the input matrix size. The runtime is optimal, even when the input is assumed to be diagonal. Our result has several applications. Length-squared sampling is used by a number of sublinear time algorithms for matrix problems, like low-rank approximation and eigenvalue approximation. Often, it is assumed that the algorithm is given access to the matrix column norms, and thus can perform length-squared sampling efficiently. Our result shows that, at least for psd matrices, we can remove this assumption. We also discuss an application to an asymptotically optimal algorithm for estimating the Frobenius norm of a psd matrix to relative error. Finally, we show that our sampling algorithm yields a very simple sublinear time algorithm for the robust psd low-rank approximation problem introduced by Bakshi et al. (FOCS, 2020), which nearly matches the more complex method developed there.

• The paper introduces a rejection-sampling algorithm that samples exactly from a PSD matrix’s length-squared distribution using O(n) expected entry queries, with termination within O(n log(1/δ)) queries at failure probability δ.
• The paper proves these costs are optimal, including Ω(min{mn, n²}) queries for m independent samples and Ω(min{n/ε², n²}) queries for relative-error Frobenius norm estimation.
• The paper removes column-norm access assumptions in eigenvalue and low-rank approximation algorithms, while extending sampling to ℓp distributions and robust settings with explicit accuracy and noise-dependent query bounds.

Overview

Length-squared sampling—drawing a column of a matrix with probability proportional to its squared $\ell_2$-norm—is a foundational primitive in randomized numerical linear algebra, underlying algorithms for low-rank approximation, approximate matrix multiplication, randomized Kaczmarz, and sublinear eigenvalue estimation. In the entry-query model, however, exact length-squared sampling has generally been assumed rather than implemented: computing even a single column norm requires $n$ queries, and computing the full sampling distribution requires reading all $n^2$ entries. The paper under review, by Bhattacharjee, Epperly, Musco, and Tian (2608.12503), shows that for positive-semidefinite (psd) matrices this assumption can be removed entirely: an extremely simple rejection-sampling procedure samples exactly from the length-squared distribution using only $O(n)$ expected entry queries and arithmetic operations. This is asymptotically optimal even among diagonal matrices, and it immediately simplifies several known sublinear-time algorithms that previously required column norms as input.

The sampling algorithm

The core algorithm exploits the psd off-diagonal inequality: for psd $A$, every entry satisfies $A_{ij}^2 \le A_{ii}A_{jj}$, since the $\{i,j\}$-principal submatrix has nonnegative determinant. This inequality is what rules out the $\Omega(n^2)$ hard case for general matrices—a single nonzero entry hidden at an arbitrary location—for on a psd matrix large off-diagonal entries can only occur where both corresponding diagonal entries are large.

The procedure reads the full diagonal once ($n$ queries) to form the diagonal distribution $\mathcal{D}(i) = A_{ii}/\operatorname{tr}(A)$, then repeatedly proposes independent indices $n$0, accepting the pair $n$1 with probability $n$2. By the off-diagonal inequality this acceptance probability lies in $n$3. Conditioning on acceptance, the pair is distributed according to the entry-squared distribution $n$4, so each marginal index follows the length-squared distribution exactly.

The runtime analysis hinges on the per-round acceptance probability

$n$5

where the last step is Cauchy–Schwarz. The number of rounds is geometric with mean $n$6, giving $n$7 expected total cost; with probability $n$8 the algorithm terminates within $n$9 operations, and truncating at that point yields a sampler within $n^2$0 total variation distance of the target distribution. Since many downstream applications require only approximate length-squared sampling (probability at least a constant fraction of the true value), constant $n^2$1 typically suffices.

The same rejection-sampling template extends to $n^2$2 column-norm sampling for any $n^2$3: replacing diagonal weights $n^2$4 by $n^2$5 yields an $n^2$6-time sampler from the $n^2$7th-power entry distribution, with an identical proof.

Optimality

The $n^2$8 bound is tight in two senses. First, even positive diagonal matrices require reading the full diagonal to perform length-squared sampling, so $n^2$9 is unavoidable. Second—and less obvious—the paper proves that amortization across multiple samples does not help: any algorithm outputting $O(n)$0 independent length-squared samples must read $O(n)$1 entries in expectation. The lower bound reduces from a "heavy column" problem on a randomly permuted block matrix with $O(n)$2 heavy columns: identifying a heavy column requires finding one of $O(n)$3 nonzero off-diagonal entries hidden among $O(n)$4 positions, costing $O(n)$5 queries via Yao's minimax principle. Thus repeating the rejection sampler $O(n)$6 times, at $O(n)$7 total cost, is optimal up to constants.

Frobenius norm estimation

The same proposal mechanism yields an unbiased estimator for $O(n)$8: for $O(n)$9 independently, the variable $A$0 satisfies $A$1 and, by the off-diagonal inequality, $A$2. Averaging $A$3 iid copies and rescaling by $A$4 gives a $A$5-relative approximation to $A$6 with probability $A$7, in $A$8 time. The paper also proves a matching lower bound of $A$9 queries, via a permutation-hidden construction where distinguishing two instances requires estimating the mean of a Boolean vector whose means differ by $A_{ij}^2 \le A_{ii}A_{jj}$0, combining standard mean-estimation lower bounds with a certificate-counting argument showing that finding $A_{ij}^2 \le A_{ii}A_{jj}$1 nonzero off-diagonal entries costs $A_{ij}^2 \le A_{ii}A_{jj}$2 queries. Both the upper and lower bounds are therefore optimal.

Applications

Eigenvalue approximation. The length-squared-based eigenvalue approximation algorithms of Bhattacharjee et al. and Swartworth–Woodruff assumed access to column norms and $A_{ij}^2 \le A_{ii}A_{jj}$3. Combining the new sampler with the Frobenius norm estimator removes these assumptions for psd inputs. Two technical adaptations are needed: a $A_{ij}^2 \le A_{ii}A_{jj}$4 approximation to $A_{ij}^2 \le A_{ii}A_{jj}$5 suffices both for the zeroing-out (truncation) step and for approximately rescaling sampled rows—the distortion introduced by the estimated normalizer is absorbed into the accuracy parameter via a clean rescaling argument—and a Poissonization trick converts independent length-squared samples into the row-wise Poisson sampling scheme required by the prior analysis. The resulting algorithm returns an additive $A_{ij}^2 \le A_{ii}A_{jj}$6 approximation to the full spectrum using $A_{ij}^2 \le A_{ii}A_{jj}$7 entry queries in expectation.

Additive-error low-rank approximation. Sampling $A_{ij}^2 \le A_{ii}A_{jj}$8 columns via the new sampler, computing their top-$A_{ij}^2 \le A_{ii}A_{jj}$9 singular vectors $\{i,j\}$0, and forming the right factor $\{i,j\}$1 by leverage-score-based approximate matrix multiplication yields factors $\{i,j\}$2 with $\{i,j\}$3 using $\{i,j\}$4 queries and $\{i,j\}$5 arithmetic operations. This is substantially simpler than prior sublinear psd low-rank approximation algorithms, though strictly weaker than the relative-error guarantee of Bakshi, Chepurko, and Woodruff, who achieve $\{i,j\}$6 with $\{i,j\}$7 queries—matching a known lower bound up to polylogs. The paper explicitly poses as open whether $\{i,j\}$8 length-squared samples suffice for additive-error approximation, which would close this gap.

Robust low-rank approximation. In the corruption model $\{i,j\}$9 with $\Omega(n^2)$0, the paper adapts its sampler to the clipped matrix $\Omega(n^2)$1 obtained by clipping each entry to $\Omega(n^2)$2. Clipping provably does not increase the noise ($\Omega(n^2)$3), and clipping induces a modified off-diagonal inequality $\Omega(n^2)$4 involving the diagonal corruption parameter $\Omega(n^2)$5, enabling rejection sampling from $\Omega(n^2)$6's length-squared distribution in $\Omega(n^2)$7 expected queries per sample. Plugging this into the additive-error LRA algorithm gives the same error guarantee as Bakshi et al.—$\Omega(n^2)$8—while removing their row-wise noise-spreading assumption entirely, at the cost of an extra $\Omega(n^2)$9 factor in query complexity ($n$0 versus $n$1).

Limitations and open questions

The results are confined to psd matrices; for general matrices, length-squared sampling provably requires $n$2 queries, so no extension along those lines is possible without additional structure. The robust LRA result inherits the $n$3 error floor inherent to the model and pays a factor-$n$4 query overhead relative to the sample-optimal algorithm of Bakshi et al.; whether the simpler approach can be brought to $n$5 samples is left open, as noted above. The eigenvalue approximation carries polylogarithmic factors ($n$6) inherited from the prior analyses, which may not be tight. Finally, the optimality result for multiple samples assumes independence of the drawn indices; correlated or adaptive sampling schemes fall outside the lower bound's scope.

Conclusion

This paper establishes that exact length-squared sampling on psd matrices is achievable in optimal $n$7 expected time via a three-line rejection-sampling scheme built on the psd off-diagonal inequality, together with matching lower bounds for both single-sample and multi-sample regimes and for Frobenius norm estimation. Beyond its intrinsic simplicity, the primitive eliminates the column-norm access assumption embedded in several existing sublinear algorithms, yielding streamlined implementations for spectral approximation and (robust) low-rank approximation of psd matrices.