Fast approximation of matrix coherence and statistical leverage

Petros Drineas; Malik Magdon-Ismail; Michael W. Mahoney; David P. Woodruff

Source-linked AI summary

Fast approximation of matrix coherence and statistical leverage

Petros Drineas, Malik Magdon-Ismail, Michael W. Mahoney, David P. Woodruff

arXiv:1109.3843v2Sep 18, 2011cs.DScs.DMcs.LG

TL;DR

Exact statistical leverage-score computation costs Θ(nd^2), motivating faster algorithms for large matrices. The paper uses randomized sketches to approximate all scores and extends the approach to cross-leverage scores, general matrices, and streaming settings.

Problem
Exact leverage-score computation requires Θ(nd^2) time, motivating algorithms that run in o(nd^2) time.
Method
The algorithm constructs a randomized sketch of the input matrix and computes the Euclidean norms of its rows to estimate leverage scores.
Results
O(nd log n/ε^2) time yields relative-error approximations to every leverage score, versus Θ(nd^2) for the naïve algorithm.
Takeaways & Limitations
The sketching framework also supports approximating large cross-leverage scores and extends to general fat matrices and streaming environments.
Takeaways & Limitations
Empirical performance depends heavily on numerical implementation and whether matrices are tall or more general.

Abstract

from arXiv · show

The statistical leverage scores of a matrix $A$ are the squared row-norms of the matrix containing its (top) left singular vectors and the coherence is the largest leverage score. These quantities are of interest in recently-popular problems such as matrix completion and Nyström-based low-rank matrix approximation as well as in large-scale statistical data analysis applications more generally; moreover, they are of interest since they define the key structural nonuniformity that must be dealt with in developing fast randomized matrix algorithms. Our main result is a randomized algorithm that takes as input an arbitrary $n \times d$ matrix $A$, with $n \gg d$, and that returns as output relative-error approximations to all $n$ of the statistical leverage scores. The proposed algorithm runs (under assumptions on the precise values of $n$ and $d$) in $O(n d \log n)$ time, as opposed to the $O(nd^2)$ time required by the naïve algorithm that involves computing an orthogonal basis for the range of $A$. Our analysis may be viewed in terms of computing a relative-error approximation to an underconstrained least-squares approximation problem, or, relatedly, it may be viewed as an application of Johnson-Lindenstrauss type ideas. Several practically-important extensions of our basic result are also described, including the approximation of so-called cross-leverage scores, the extension of these ideas to matrices with $n \approx d$, and the extension to streaming environments.

1 Introduction

Statistical leverage scores and coherence quantify nonuniformity in matrix singular vectors and support randomized matrix algorithms. The paper introduces randomized approximations that are faster than orthogonal-basis computation, with extensions to cross-leverage scores and broader matrix settings.

Definitions: Statistical leverage scores are squared row norms of the top left singular-vector matrix, while coherence is the largest leverage score.They are also diagonal and maximum-diagonal quantities of the projection matrix onto the column span.
Motivation: These quantities are used in large-scale data analysis, matrix completion, Nyström-based low-rank approximation, and randomized matrix algorithms.Leverage scores also support outlier detection in regression diagnostics and importance sampling.
Main result: O(nd log n/ε^2) time computes relative-error approximations to all leverage scores, versus Θ(nd^2) for the naïve orthogonal-basis approach, under stated n and d assumptions.Theorem 1 gives success probability at least 0.8 for full-rank n×d matrices with n ≫ d and ε ∈ (0,1/2].
Scope and caveats: The theoretical treatment assumes full column rank, while numerical treatment of rank-deficient matrices is left for future work.Empirical performance is reported as dependent on numerical implementation and matrix shape, with broader empirical analysis deferred.
Main result: A randomized sketch based on a Fast Johnson–Lindenstrauss Transform enables leverage-score estimation by computing row norms of the sketch.The same sketch supports cross-leverage estimation, while identifying large inner products is additionally required for that task.
Extensions: Under d ln d = o(n) and ln n = o(d), the algorithm runs in o(nd^2) time for constant ε and approximates coherence to relative error.The paper also describes extensions to general fat matrices, streaming environments, and large cross-leverage scores.

2 Preliminaries on linear algebra and fast random projections

The preliminaries define the matrix and norm notation, statistical projection concepts, and the randomized projections used to accelerate computations.

Linear algebra notation: A thin SVD represents A as UΣV^T, with U and V spanning the relevant row and column spaces.The Moore–Penrose pseudoinverse is A† = VΣ^-1U^T.
Random projections: An ε-JLT projects points while approximately preserving their geometry with high probability.The construction uses sparse random entries and a dimension determined by the accuracy and failure probability.
Random projections: An ε-FJLT preserves the orthogonality of U while enabling fast multiplication by a projection matrix.For any X, ΠX can be computed in O(nd ln r) time.
Random projections: The SRHT combines random sign flips, a Hadamard transform, and uniform row sampling to form an efficient projection.The transform spreads vector energy and supports fast computation of selected transformed entries.

3 Our main algorithmic results

The paper constructs randomized sketches that approximate leverage and large cross-leverage scores without explicitly computing a full orthogonal basis or all pairwise products.

Basic approach: The core sketch A(Π1A)†Π2 combines an ε-FJLT and an ε-JLT to approximate leverage-related quantities.The same sketch supports both diagonal leverage-score and large cross-leverage-score approximations.
Basic approach: Computing an orthogonal basis, pseudoinverse, and matrix products naïvely requires O(nd^2) time for arbitrary tall matrices.These operations are the main bottlenecks targeted by the randomized construction.
Leverage scores: Uniformly sampling far fewer rows can fail to preserve rank, so the algorithm first applies an ε-FJLT before computing a smaller pseudoinverse.The resulting reduced matrix Π1A is used to construct an approximate basis for A.
Leverage scores: Algorithm 1 returns estimates ˜ℓi for every row after constructing a randomized sketch and computing squared row norms.A QR factorization may replace the SVD when constructing the orthogonalizer.
Cross-leverage scores: Pairwise dot products in the sketch approximate leverage and cross-leverage scores, while heavy-hitter processing avoids enumerating all n^2 pairs.The full pairwise approach has running time T(Ω) + O(n^2d), with r2 = O(ε^-2 ln n).
Cross-leverage scores: Algorithms 2 and 3 identify large cross-leverage pairs and return their estimates using the sketch from Algorithm 1.Algorithm 2 adjusts the heavy-pair threshold before returning the detected pairs.

4 Proofs of our main theorems

The proofs establish that the two randomized projections preserve the geometry needed for leverage estimates and enable efficient identification of large cross-leverage pairs.

Proof strategy: The analysis conditions on Π1 being an ε-FJLT and Π2 being an ε-JLT, yielding theorem success probability at least 0.8.The estimates are formed from squared norms of rows of the sketched matrix.
Runtime: The sketch-based approach improves over the exact basis-plus-pairwise-inner-products cost by a factor of d/ln n.The exact procedure requires O(nd^2 + n^2d) time.
Proof strategy: Π2 preserves pairwise inner products because it is an ε-JLT for n^2 points, improving runtime without changing the relevant geometry substantially.The proof derives inner-product preservation from norm preservation of vectors and pairwise sums.
Leverage-score proof: The FJLT makes Π1U approximately orthogonal, which supports approximating the leverage scores through the smaller matrix Π1A.The pseudoinverse can be computed in O(nd ln r1 + r1d^2) time.
Approximation guarantees: The diagonal case of the inner-product approximation yields relative-error leverage estimates, while off-diagonal terms yield additive-error cross-leverage estimates.The proof combines the approximation bounds for the intermediate sketch and the second projection.
Runtime: The sketch construction has total running time O(nd ln r1 + ndr2 + r1d^2 + r2d^2).With r2 = O(ε^-2 ln n), the resulting asymptotic cost is determined by the reduced dimensions and projection operations.
Heavy-pair detection: Algorithm 3 finds all heavy pairs of a matrix in O(nr + κr^2 + n ln n) time.The method sorts row norms and searches only norm-heavy pairs.
Heavy-pair detection: The cross-leverage proof combines preserved inner products with fast heavy-pair identification to relate detected pairs to the actual scores.Returned pairs satisfy a lower bound involving d/κ and an additive error term.

5 Extending our algorithm to general matrices

The paper extends leverage-score approximation to general matrices by relaxing the ill-posed target of the best rank-k approximation to any sufficiently accurate rank-k approximator. Randomized spectral- and Frobenius-norm algorithms then estimate normalized leverage scores for such approximating matrices.

Problem formulation: The best rank-k approximation target can be ill-posed because its top-k singular space may be nonunique, as for the identity matrix.When singular values near the kth position are indistinguishable, leverage scores for the associated singular vectors may vary substantially.
Problem formulation: The paper defines S_ε as rank-k matrices whose spectral- or Frobenius-norm error is at most (1 + ε) times the best rank-k error.Approximations may target leverage scores of any X ∈ S_ε rather than the potentially nonunique A_k.
Problem formulation: β-approximations are outputs that approximate the normalized leverage scores of some X ∈ S_ε, whose left singular vectors determine those scores.This relaxation removes the ill-posedness while retaining matrices sufficiently close to A_k for the approximation task.
Spectral norm: Algorithm 4 computes normalized leverage-score approximations for general matrices with probability at least 0.7 under its stated rank and accuracy conditions.Its outputs are interpreted as approximations to the normalized leverage scores of a suitable matrix in S_ε.
Spectral norm: The spectral-norm algorithm uses a Gaussian sketch Π ∈ R^{d×2k} and represents a suitable rank-k matrix as X = BY.The sketch matrix B can be computed in O(ndk^q) time, and the construction yields X ∈ S_ε with probability at least 0.9.
Frobenius norm: For Frobenius-norm approximators, the algorithm returns normalized leverage scores exactly equal to those of a closed-form rank-k matrix X with constant probability.The resulting scores form a 1-approximation to the normalized leverage scores associated with the rank-k task.

6 Discussion

The discussion connects leverage-score estimation to prior estimators, under-constrained least squares, and streaming algorithms. It highlights faster approximate least-squares computation and one-pass streaming estimates of leverage-score statistics.

Relation to prior estimators: A prior estimator achieved an O(log^2 n) approximation in roughly O(nd^2 / log n) time, while the paper’s approach maintains positivity in its row-norm estimates.The prior method uses random projections, truncation, and renormalization; the paper contrasts this with a sketch of AA† that preserves positivity.
Under-constrained least squares: Under-constrained least squares can be solved by sampling columns of A with probabilities proportional to the leverage scores of A^T.Algorithm 6 uses the sampled columns to construct an estimator of the minimal ℓ2-norm solution.
Under-constrained least squares: Theorem 5 gives a high-probability approximation guarantee for Algorithm 6 when A has full row rank, n ≪ d, and sampling probabilities satisfy the stated condition.The algorithm and its running-time analysis depend on the supplied accuracy, failure-probability, and sampling parameters.
Under-constrained least squares: Assuming ε, β, and δ are constants and n log n = o(d), Algorithm 6 runs in o(n^2d) time.The paper notes that leverage scores can be approximated using Theorem 1 and Algorithm 1 instead of computed exactly.
Streaming environments: In streaming settings, leverage-score sampling can use O(sdε^-2 log^4 n) bits of space and one pass, with a second pass optionally retrieving sampled rows.The method returns row identities according to their leverage-score distribution.
Streaming environments: The entropy of the leverage-score distribution can be estimated within additive ε using O(dε^-4 log^6 n log^14(1/ε)) bits of space and one pass.The construction maintains sketches of TA and ΠA, then computes the estimate after the pass.