Iterative Sampling Algorithm

• Iterative sampling is a method that alternates between reducing data dimensionality via random projections and recovering refined sampling probabilities to preserve matrix structures.
• The algorithm uses leverage scores and generalized stretch to efficiently approximate tall-and-skinny matrices while maintaining a (1 ± ε) norm guarantee.
• It achieves state-of-the-art computational efficiency for large-scale regression and graph sparsification, balancing accuracy with reduced sample sizes.

An iterative sampling algorithm is a class of randomized algorithm that progressively constructs a high-fidelity sample or summary from data or a distribution by alternating repeated rounds of coarse approximation and refinement. In computational mathematics and large-scale data analysis, such algorithms are crucial for reducing problem dimensionality, controlling sample quality, and achieving resource efficiency—particularly in settings where direct methods are computationally prohibitive or when the data exhibits highly nonuniform “importance.” Recent advances have unified concepts from randomized numerical linear algebra, matrix sketching, graph sparsification, and leverage score sampling under unified iterative sampling schemes.

1. Iterative Reduction and Recovery Framework

The archetypal iterative sampling algorithm for tall-and-skinny matrices (where $n \gg d$) operates as a two-phase process (Li et al., 2012):

• Reduction Phase: The algorithm repeatedly compresses the input matrix $A$ (or its approximation at level $\ell$, denoted $A^{(\ell)}$) by partitioning rows into blocks (e.g., of size $R$) and mapping these blocks to lower-dimensional spaces via random projections (e.g., multiplying with a random Gaussian matrix $U$). Each reduction approximately preserves the column space structure, and after $L$ reductions, the algorithm obtains a geometrically smaller instance $A^{(L)}$.
• Recovery Phase (Backward Pass): Starting from the highly compressed $A^{(L)}$, the procedure propagates improved approximations of the row sampling probabilities—quantified as leverage scores or generalized “stretch”—up through the sequence of reduced matrices. At every level, these estimates are tightened and “lifted” towards the original matrix $A$, using the small approximants constructed in reduction.

The process ensures that, at every iteration, the sampled matrix $B^{(\ell)}$ is a $(1 \pm \epsilon)$-approximation for $A^{(\ell)}$ in the sense of norm preservation.

Invariant: For all $x \in \mathbb{R}^d$, $$(1 - \epsilon)\|A x\|_2 \leq \|B x\|_2 \leq (1 + \epsilon)\|A x\|_2$$.

2. Leverage Scores and Generalized Stretch

Leverage scores are central to iterative sampling and quantify the influence of each row in the column space:

$\tau_{(i)} = a_i (A^\top A)^{+} a_i^\top,$

where $(A^\top A)^+$ is the Moore–Penrose pseudoinverse. These sum to $\mathrm{rank}(A) \leq d$.

The algorithm generalizes this via the stretch relative to a reference matrix $B$: $\mathrm{STR}_B(a_i) = a_i (B^\top B)^+ a_i^\top,$ and global stretch: $$\mathrm{STR}_B(A) = \| (B B^\top)^{-1/2} A^\top \|_F^2.$$ Coarse approximations to these scores, as obtained during reduction, are robust guides for sampling and are successively refined during the recovery phase.

Key insight: Even loose upper bounds on leverage scores suffice to preserve norm structure in subsampling, and these can be iteratively improved without full (costly) recomputation.

3. Algorithmic Complexity and Theoretical Guarantees

The iterative algorithm in (Li et al., 2012) achieves, for a given $\epsilon > 0$, with high probability (failure probability ${\leq} d^{-c}$ for any constant $c$):

• Output: A matrix $B$ composed of appropriately rescaled rows of $A$ with $O(d \log d \, \epsilon^{-2})$ rows.
• Guarantee: For all $x \in \mathbb{R}^d$,

$$(1 - \epsilon)\|A x\|_2 \leq \|B x\|_2 \leq (1 + \epsilon)\|A x\|_2.$$

• Time Complexity:

$$O(\mathrm{nnz}(A) + d^{\omega + \theta} \epsilon^{-2}),$$

where $\mathrm{nnz}(A)$ is the number of non-zeros in $A$, $\omega$ is the matrix multiplication exponent (currently $\sim$2.3727), and $\theta > 0$ is arbitrarily small.

This matches or improves upon “one-shot” random projection approaches, especially regarding the dependence of sample size on $d$ (moving from quadratic to nearly linear), and offers sharply defined trade-offs between computational cost and approximation quality.

4. Mathematical Structure and Formulation

The central properties maintained during the iterative process are matrix inequalities: $$(1 - \epsilon)A^\top A \preceq B^\top B \preceq (1 + \epsilon)A^\top A,$$ where $\preceq$ denotes the Loewner partial order for semidefinite matrices.

Additionally, the use of upper bounds for the sum of leverage scores ($\sum_i \tau_{(i)} \leq d$) and the connection to stretch/frobenius norms underpins the estimation and refinement strategy.

5. Application Domains and Data Reduction

Regression and Sampling for Optimization

The algorithm is specifically constructed to address large-scale least-squares ($\ell_2$) and $\ell_p$ regression: $\min_x \|A x - b\|_p,$ where direct manipulation of $A$ is prohibitive for $n \gg d$. Substituting $A$ by the succinct $B$ from iterative sampling allows reduction of the problem to $O(d \log d \epsilon^{-2})$ constraints, with guarantees that solutions carry over up to $(1 \pm \epsilon)$ distortion.

Preserving Data Structure: Because each row in $B$ is an exact (rescaled) copy of a row in $A$, the procedure is “structure-preserving.” This is critical for downstream machine learning or signal processing applications where data provenance is essential.

Streaming and Large-Scale Environments

The iterative approach is especially suited for environments with restricted access models (e.g., streaming), as each phase processes only a manageable, summary-sized sketch.

6. Connections to Graph Sparsification and Robustness

Iterative sampling as presented in (Li et al., 2012) is conceptually and technically linked to graph sparsification. In that domain, the goal is to approximate the Laplacian quadratic form of a graph via a sparse subgraph, often by sampling edges according to their effective resistance—a direct analog of leverage scores for matrices. The iterative method draws on these ideas: concentration bounds (e.g., matrix Chernoff inequalities), combinatorial preconditioning, and alternation between coarse (spanner-like) reductions and finer recovery.

Robustness Mechanism: Even if the first round of approximations is rough, subsequent iterations improve the quality, analogous to how a rough sparse graph can be incrementally improved to respect quadratic forms.

7. Implications and Impact

The iterative sampling paradigm enables:

• Tighter theoretical sample complexity for matrix approximation in regression.
• Algorithms with input-sparsity running time (scaling with $\mathrm{nnz}(A)$) and minimal expensive matrix operations.
• Robustness to errors in importance estimation, due to backward refinement.
• Direct applicability to graph algorithms, randomized linear algebra, and large-scale data analysis where preserving the inherent structure of the underlying matrix or graph is desirable.

By unifying random projection-based sketching, leverage score estimation, and graph sparsification, iterative sampling algorithms offer an extensible framework for scalable linear algebra and optimization in modern data-intensive applications (Li et al., 2012).