Technical Summary: Fast One-Pass Sparse Approximation of the Top Eigenvectors of Huge Approximately Low-Rank Matrices

Problem Statement

The paper addresses the challenge of computing the top eigenvectors of extremely massive matrices $A \in \mathbb{C}^{N \times N}$ that are too large to store in memory (e.g., matrices with $\sim 10^{16}$ entries). Traditional methods like Power Iteration or Lanczos often struggle with high computational costs and slow convergence in this regime. Furthermore, standard streaming algorithms for eigenvector approximation often require multiple passes over the data or assume specific structural properties that limit their applicability to sparse recovery.

The specific goal is to develop a one-pass algorithm that:

1. Processes the entries of $A$ only once.
2. Utilizes a memory footprint proportional to the size of the desired sparse eigenvector approximations (rather than the full matrix size).
3. Produces accurate sparse approximations of the leading eigenvectors, assuming the matrix is approximately low-rank and the eigenvectors are sparse or compressible.

Methodology: The MAM* Approach

The proposed method, termed MAM*, combines classical linear algebra, randomized sketching, and compressive sensing (CS) into a two-stage pipeline.

1. One-Pass Sketching (The MAM* Matrix)

Instead of storing $A$, the algorithm computes a compact sketch matrix $\tilde{A} = MAM^* \in \mathbb{C}^{m \times m}$, where $m \ll N$.

• Process: The algorithm iterates through the entries (or columns) of $A$ exactly once.
• Measurement Matrix: $M \in \mathbb{C}^{m \times N}$ is a structured, memory-efficient matrix (e.g., based on Hadamard transforms or combinatorial designs) that acts as a compressive sensing operator.
• Computation: The sketch is formed as $\tilde{A} = \sum_{j,k} a_{j,k} M e_j e_k^* M^* = M A M^*$. This is computed incrementally without ever storing $A$ or the full product $MA$.

2. Eigenvector Recovery via Compressive Sensing

Once $\tilde{A}$ is computed, the algorithm proceeds as follows:

1. Local Eigendecomposition: Compute the top $\ell$ eigenvectors $\{\tilde{u}_j\}$ of the small matrix $\tilde{A}$.
2. Global Recovery: Use these local eigenvectors as "measurements" in a fast compressive sensing reconstruction algorithm to recover the global sparse eigenvectors $\{u_j\}$ of the original matrix $A$.
   • The core theoretical insight is that if $A$ is approximately low-rank and $M$ satisfies specific properties (like the Restricted Isometry Property or Johnson-Lindenstrauss embedding), then the eigenvectors of $\tilde{A}$ are close to the projected eigenvectors of $A$ (i.e., $\tilde{u}_j \approx M u_j$).
   • Standard CS algorithms (e.g., CoSaMP) or sublinear-time CS algorithms can then invert these measurements to recover $u_j$.

The paper outlines two specific implementations (Algorithms 1 and 2):

• Linear-Time Variant: Uses standard CS algorithms (like CoSaMP) requiring $O(N \log N)$ time per eigenvector.
• Sublinear-Time Variant: Uses specialized combinatorial CS algorithms (based on bit-testing matrices and Khatri-Rao products) that run in time sublinear in $N$, depending primarily on the sparsity $s$ and sketch size $m$.

Key Contributions and Theoretical Results

1. Generalized Eigenvector Perturbation Theory

The paper establishes new perturbation bounds relating the eigenvectors of $A$ to those of $MAM^*$. Unlike previous work focusing on Gaussian sketches, this theory applies to structured matrices $M$ that admit fast matrix-vector multiplication.

• Theorem 3.3: Proves that for Hermitian, Positive Semi-Definite (PSD) matrices with exponentially decaying eigenvalues ($\lambda_j = cq^j$), the error between the sketch eigenvector $\tilde{u}_j$ and the projection $M u_j$ is bounded by $O(\sqrt{\epsilon} \cdot q^{1-j})$. This holds even when $A$ is only approximately low-rank.

2. Two Main Theorems on Approximation Accuracy

The authors prove two main theorems (Theorem 1.1 and Theorem 1.3) guaranteeing the accuracy of the recovered sparse eigenvectors:

• Theorem 1.1 (Linear-Time Recovery):

  • Input: Sketch $\tilde{A}$ and a linear-time CS algorithm.
  • Guarantee: With high probability, the recovered vector $u'_j$ satisfies:
    $$\min_{\phi} \| e^{i\phi} u_j - u'_j \|_2 < C \cdot \max \left( \eta, \frac{1}{\sqrt{s}}\|u_j - (u_j)_s\|_1 + \sqrt{\epsilon} \cdot q^{1-j} \right)$$
  • Complexity: $O(m^3 + \ell N \log N \cdot \log(1/\eta))$ time; $O(N)$ memory.

• Theorem 1.3 (Sublinear-Time Recovery):

  • Input: Sketch $\tilde{A}$ and a sublinear-time CS algorithm.
  • Guarantee: Similar error bounds, with the error term involving a constant $\beta_M^A$ related to the "flatness" of the noise vector.
  • Complexity: $O(m^3)$ time (sublinear in $N$); $O(N)$ memory (or $O(m)$ if $D$ is not stored).
  • Trade-off: The measurement requirement $m$ scales quadratically with $r$ and $s$ (worse than Theorem 1.1), but the runtime is significantly faster for huge $N$.

3. Experimental Validation

The paper presents preliminary experiments on matrices with $N \approx 2^{27} \approx 1.3 \times 10^8$ (total entries $\sim 10^{16}$).

• Setup: Random symmetric PSD matrices with rank $r \in [20, 50]$ and sparsity $s \in [30, 100]$.
• Results:
  • The measurement error (distance between $\tilde{u}_j$ and $M u_j$) remains low and stable across varying ranks and sparsities.
  • The CS recovery error increases as sparsity decreases, consistent with standard CS theory, but the sketching process itself does not appear to be the bottleneck.
  • The method successfully recovers eigenvectors even when the spectral decay is linear rather than exponential, and in the presence of dominant rank-1 Gaussian noise (up to a threshold).
  • The constant $\beta_M^A$ in the sublinear theorem was empirically observed to be $O(1)$ (specifically $<8$) in all tests.

Significance and Claims

The paper claims that MAM* is the first known one-pass algorithm capable of providing sparse approximations for the top eigenvectors of extremely massive matrices. Its significance lies in:

1. Memory Efficiency: It reduces the memory footprint from $O(N^2)$ (storing $A$) or $O(N)$ (storing columns) to $O(N)$ for the sketch generation and $O(m)$ for the final recovery, where $m$ is small.
2. One-Pass Constraint: It operates strictly in a single pass over the data, making it suitable for streaming settings where re-accessing $A$ is prohibitive.
3. Sublinear Runtime: The sublinear-time variant (Theorem 1.3) allows for the recovery of eigenvectors in time that does not scale linearly with the matrix dimension $N$, depending instead on the sparsity and sketch size.
4. Broad Applicability: While the theory assumes Hermitian PSD matrices with specific spectral decay, the authors note that the framework is flexible enough to handle generalized eigenvalue problems and potentially non-Hermitian matrices via asymmetric sketches (future work).

The authors position this work as a bridge between randomized numerical linear algebra (sketching) and compressive sensing, enabling the analysis of "stupendously large" datasets (e.g., for sparse PCA) that were previously intractable for standard eigensolvers.