Fast Length-Squared Sampling for Positive-Semidefinite Matrices

Summary

This paper presents a simple rejection-sampling-based algorithm for performing length-squared sampling on an n × n positive-semidefinite (psd) matrix, i.e., sampling a column with probability proportional to its squared l2-norm. The main result is stated in Theorem 1: "Let A ∈ Rn×n be nonzero and psd. There exists an algorithm (Algorithm 1) that outputs a random index I ∈ [n] sampled from the length-squared distribution Pr[I = i] = ∥A∗,i∥22 / ∥A∥2F for all i ∈ [n]. The algorithm reads O(n) entries of A and performs O(n) arithmetic operations in expectation."

Algorithm 1 works as follows: it reads the diagonal of A and defines the diagonal sampling distribution D(i) = Aii / tr(A). It then samples two indices I and J independently from D, and accepts the proposal (I, J) with probability A2IJ / (AII AJJ). This is a valid probability because of Fact 3, which states that for psd matrices, A2ij ≤ Aii Ajj for all i, j. Conditional on acceptance, (I, J) is drawn from the squared-entry distribution E(i, j) = A2ij / ∥A∥2F, and both I and J individually follow the length-squared distribution. The key proof step is showing the acceptance probability is Ω(1/n), giving O(n) expected runtime. The algorithm requires 2n entry queries and O(n) arithmetic operations in expectation, and with probability 1−δ terminates using at most n + ⌈n log(1/δ)⌉ entry accesses.

The runtime is proven optimal in Theorem 5: "Any algorithm that, given any nonzero psd matrix A ∈ Rn×n outputs m indices i1,..., im ∈ [n] sampled independently from the length-squared distribution (1.1) must read Ω(min mn, n2) entries of A in expectation in the worst case." This holds even for sampling multiple indices, showing one cannot do better than repeating the algorithm m times.

The paper also presents a Frobenius norm approximation algorithm (Algorithm 2) in Theorem 2: "Let A ∈ Rn×n be a psd matrix. For any ϵ, δ ∈ (0, 1) there exists an algorithm (Algorithm 2) that outputs an estimate F̂ satisfying (1−ϵ)∥A∥2F ≤ F̂ ≤ (1+ϵ)∥A∥2F with probability ≥ 1−δ. The algorithm reads O(n log(1/δ)/ϵ2) entries of A and performs O(n log(1/δ)/ϵ2) arithmetic operations." The algorithm samples I, J independently from the diagonal distribution and averages the unbiased estimator X:= A2IJ / (AII AJJ) over s = ⌈3n log(2/δ)/ϵ2⌉ trials, then rescales by tr(A)2. Theorem 7 shows this is asymptotically optimal: any algorithm estimating the Frobenius norm to (1±ϵ) relative error with constant probability must read Ω(min n/ϵ2, n2) entries in expectation.

Additionally, Theorem 17 shows a small modification of the algorithm can sample columns according to the lp column-norm distribution for any p ≥ 0, running in O(n) expected time.

The paper discusses several applications:

1. Eigenvalue approximation (Theorem 9): For psd matrices, the sublinear time eigenvalue approximation algorithm of [BDD+24, SW25] can be implemented using entry queries alone, without assuming access to column norms. The algorithm returns an (ϵ∥A∥F)-additive approximation to the spectrum of A with probability at least 2/3, querying O((n/ϵ2) log4 n log(21/ϵ)) entries in expectation.

2. Additive-error low-rank approximation (Theorem 12): The algorithm outputs factors X, Y ∈ Rn×k of a rank-k approximation XY⊤ satisfying ∥A − XY⊤∥2F ≤ ∥A − Ak∥2F + ϵ∥A∥2F with probability at least 2/3. It uses O(nk/ϵ2) entry queries and O(n(k/ϵ2)ω−1) arithmetic operations in expectation, where ω is the exponent of fast matrix multiplication. This is based on combining length-squared sampling (Lemma 10 from [DKM06b]) with approximate matrix multiplication (Lemma 11 from [DKM06a]).

3. Robust low-rank approximation (Theorem 16): For a corrupted matrix B = A + N where A is psd and N is arbitrary noise with ∥N∥2F ≤ η∥A∥2F, the algorithm produces rank-k factors satisfying ∥A − XY⊤∥2F ≤ ∥A − Ak∥2F + (ϵ + 9√η)∥A∥2F with probability at least 2/3. It queries O(φ2max nk/ϵ2) entries of B in expectation, where φmax ≥ max 1, maxi∈[n] Aii/Bii is a given corruption parameter. This matches the error bound of [BCW20] but removes their assumption (A2) on row-wise noise spread, at the cost of an extra 1/ϵ factor in query complexity.

The paper notes that for general (non-psd) matrices, length-squared sampling requires Ω(n2) entry accesses, so the psd structure is essential. The optimality proofs use hard instances involving random permutation matrices and diagonal-plus-block structures. The paper concludes that the algorithm's simplicity and minimal assumptions make it widely applicable across different linear algebraic tasks.

Improvement 1: Sublinear-Time Matrix Analysis Without Precomputed Norms

The AI system can now perform length-squared sampling on positive-semidefinite (PSD) matrices using only O(n) entry reads, eliminating the need for precomputed column norms or full matrix access. This enables real-time spectral analysis, low-rank approximation, and Frobenius norm estimation on massive PSD matrices (e.g., covariance matrices, graph Laplacians) stored in distributed or streaming systems, where full column norm computation is infeasible.

Improvement 2: Optimal Query Complexity for Sampling and Estimation

The system can guarantee Ω(n) entry reads for length-squared sampling and Ω(n/ε2) for Frobenius norm estimation, matching theoretical lower bounds. This allows AI pipelines to allocate computational budgets precisely, knowing that no faster exact or (1±ε)-approximate algorithm exists, improving resource efficiency in large-scale machine learning (e.g., kernel matrix approximations, feature selection).

Improvement 3: Unified Framework for lp-Norm Column Sampling

By extending the rejection-sampling trick to any p ≥ 0, the AI can sample columns by lp norms (e.g., l1 for sparsity-focused tasks, l∞ for outlier detection) in O(n) expected time. This enables adaptive sampling strategies for heterogeneous data distributions, improving robustness in clustering, anomaly detection, and randomized numerical linear algebra.

Improvement 4: Noise-Robust Low-Rank Approximation Without Distributional Assumptions

The system can compute rank-k approximations of corrupted PSD matrices (e.g., sensor data with adversarial noise) with error ∥A−A k∥2 F + (ε + 9√η)∥A∥2 F, using only O(nk/ε2) entry queries. Unlike prior work, it does not require row-wise noise spread assumptions, making it applicable to arbitrary noise patterns—critical for AI systems operating on real-world, partially corrupted data (e.g., missing or malicious entries in recommendation systems).

Improvement 5: Entry-Query-Only Eigenvalue Approximation

The AI can approximate the entire spectrum of a PSD matrix to additive error ε∥A∥ F using only O((n/ε2) log4 n) entry reads, without column norm oracle access. This enables fast spectral analysis of implicit matrices (e.g., graph kernels, Hessians in deep learning) where only entries are accessible, improving model interpretability and optimization diagnostics.

Improved AI System Capabilities

• Real-time monitoring: Track spectral properties or low-rank structure of dynamically changing PSD matrices (e.g., evolving social networks, streaming sensor data) with sublinear per-update cost.

• Scalable feature engineering: Perform length-squared sampling on billion-dimensional PSD kernels for efficient kernel-based learning (e.g., Gaussian processes) without materializing the full matrix.

• Robust recommendation systems: Generate low-rank approximations of user-item interaction matrices corrupted by noise, with provable error guarantees and minimal data access.

• Privacy-preserving analytics: Estimate Frobenius norms or spectra of sensitive PSD matrices (e.g., covariance of medical records) by querying only a fraction of entries, reducing exposure of raw data.

• Automated algorithm selection: Given a PSD matrix, the system can decide between exact and approximate methods based on query budget, using the proven optimality bounds to guarantee performance.