Subspace Pursuit Algorithm for Sparse Recovery

• Subspace Pursuit is a greedy iterative algorithm for sparse recovery that leverages RIP conditions to robustly recover K-sparse signals from underdetermined systems.
• It combines support identification, least-squares estimation, and pruning to iteratively refine signal approximations with competitive computational efficiency.
• Extensions to block, decentralized, and collaborative variants broaden its applications to image recovery, sensor networks, and joint-sparse signal processing.

The Subspace Pursuit (SP) algorithm is a greedy iterative method for sparse signal recovery from underdetermined linear systems. Designed for settings where the unknown signal vector is known or hypothesized to be $K$-sparse, SP delivers computational complexity similar to Orthogonal Matching Pursuit (OMP) while attaining recovery guarantees on par with convex relaxation ($\ell_1$ minimization). The algorithm is widely analyzed under the Restricted Isometry Property (RIP) and has seen extensions for block sparsity, noisy measurement models, and decentralized settings.

1. Algorithmic Structure

Subspace Pursuit operates within the canonical linear model $y = A x + e$, where $A \in \mathbb{R}^{m \times n}$ is the measurement matrix (often $m \ll n$), $x \in \mathbb{R}^n$ is a $K$-sparse signal, $y \in \mathbb{R}^m$ the measurements, and $e$ an optional noise vector. The algorithm proceeds as follows (0803.0811, Li et al., 2014, Satpathi et al., 2014):

1. Initialization: Select support $S^0$ as the $K$ indices of largest $|A^Ty|$. Compute the least-squares estimate on $S^0$.
2. Iteration ($n=1,2,\ldots$):
   • Compute residual $r^{n-1} = y - A x^{n-1}$.
   • Identification: Select $\Gamma^n$ as the $K$ indices with largest $|A^T r^{n-1}|$.
   • Augmentation: Form $U^n = S^{n-1} \cup \Gamma^n$.
   • Least-Squares: Solve $$u^n = \arg\min_{z: \operatorname{supp}(z) \subseteq U^n} \|y - Az\|_2$$.
   • Pruning: Update support $S^n = \operatorname{supp}(H_K(u^n))$, where $H_K$ is the hard thresholding to the $K$ largest entries.
   • Compute $$x^n = \arg\min_{z: \operatorname{supp}(z) \subseteq S^n} \|y - Az\|_2$$.
   • Stop on $\|y - Ax^n\|_2 \leq \epsilon$ or if residual does not decrease.

This mechanism combines aggressive support identification (via top correlations), iterative least-squares, and rigorous pruning, thereby repeatedly correcting previous support mistakes (0803.0811).

2. Theoretical Guarantees and RIP Conditions

The recovery properties of SP depend crucially on the Restricted Isometry Property (RIP) of $A$. A matrix $A$ satisfies the RIP of order $T$ with constant $\delta_T$ if

$$(1 - \delta_T)\|v\|_2^2 \leq \|A v\|_2^2 \leq (1 + \delta_T)\|v\|_2^2$$

for all $T$-sparse vectors $v$ (Li et al., 2014).

• Exact Recovery (Noiseless): If $\delta_{3K} < 0.165$ (original), SP recovers any $K$-sparse $x$ exactly (0803.0811). Improved analyses show guarantees for $\delta_{3K} < 0.4859$ (Satpathi et al., 2014), which is not attained by OMP or ROMP under similar constraints.
• Noisy Recovery: For $\delta_{3K}<0.083$, SP achieves $\|\hat{x} - x\|_2 \leq C_K' \|e\|_2$, with $C_K'$ depending on $\delta_{3K}$ (0803.0811). Extensions establish "near-oracle" MSE under random noise and $\delta_{3K}<0.139$, i.e.,

$$\|\hat{x}_{\text{SP}} - x\|_2^2 \leq C_{\text{SP}}^2 \cdot 2(1+a)\ln N \cdot K \cdot \sigma^2$$

with $C_{\text{SP}} \leq 21.4$ for $\delta_{3K}=0.139$ (Giryes et al., 2010).

3. Iteration Complexity and Convergence Bounds

The number of SP iterations required for exact recovery has been substantially refined. The sharpest bound to date (Satpathi et al., 2014) is: $$\#\text{iter}_{\rm SP} \leq \left\lceil \frac{\ln(4/\rho_{3K}^2)}{\ln(1/\rho_{3K}^2)} K \right\rceil, \quad\rho_{3K} = \sqrt{\frac{2 \delta_{3K}^2 (1 + \delta_{3K}^2)}{1 - \delta_{3K}^2}}$$ This result leverages the $\ell_2$-decay of the missed support across iterations and block-partition arguments of the true support. The new bound is strictly smaller than previous bounds (e.g., $\lceil 1.5K/\ln(1/\rho_{3K})\rceil$ from Dai & Milenkovic) except at very small $\delta_{3K}$ (Satpathi et al., 2014).

Numerically, for $\delta_{3K}=0.3$, the new factor $c\approx 2.4$ (old $3.1$); for $\delta_{3K}=0.45$, $c\approx 5.2$ (old $9.0$).

4. Computational Complexity and Implementation

Each iteration of SP, for $A \in \mathbb{R}^{m \times n}$, involves:

• Matrix-vector product $A^T r$: $O(mn)$.
• Support selection (top $K$): $O(n\log K)$ or $O(n)$ using selection algorithms.
• Least-squares solves on $2K$ ($K$) columns: $O(K^3)$.
• Overall per-iteration cost: $O(m n + K^3)$ (Li et al., 2014, 0803.0811).

SP typically converges in $O(K)$ iterations. For $K^2 \ll n$ (very sparse regime), the total complexity is $O(m n \log K)$ (0803.0811).

In contrast with OMP (which never revisits support decisions once made), SP’s backtracking and pruning prevent error propagation; this makes it robust under ill-conditioned dictionaries (0803.0811).

5. Extensions: Block, Decentralized, and Collaborative Variants

Block and Group-Sparse Models: Recent work generalizes SP to group and block sparse recovery. The Group Projected Subspace Pursuit (GPSP) algorithm introduces Subspace Projection Criterion (SPC) for block selection and Response Magnitude Criterion (RMC) for pruning, both theoretically and practically improving support identification in block-sparse contexts under Block-RIP (He et al., 2024). The GPSP converges if the block-RIP constant $\delta_{M,2k} \lesssim 0.1188$ and achieves stable recovery with error bounded as $\|\hat c - c^*\|_2 \leq G_{M,k}\|e\|_2$, where $G_{M,k} = (1+2\delta)/( \delta(1- \delta))$.

Decentralized and Collaborative SP: SP has inspired decentralized and collaborative variants (DCSP/GDCSP) suitable for distributed sensor networks and joint sparsity pattern recovery. Nodes execute local SP steps and share only $K$-length index sets with neighbors, leveraging majority-vote fusions for global support estimation, which minimizes communication overhead. Convergence and accuracy are on par with centralized SP for similar RIP regimes (Li et al., 2014, Li et al., 2014).

6. Relationship to Other Sparse Recovery Algorithms

SP occupies a position midway between pure greedy and convex relaxation approaches:

• Unlike OMP, SP tests and retracts support candidates, supplementing correlation-based identification with projection and pruning.
• Compared to CoSaMP, which utilizes similar candidate generation but slightly different pruning and update mechanisms, SP maintains sharper iteration bounds under practical RIP constants (Satpathi et al., 2014).
• Enhanced variants, such as Subspace Thresholding Pursuit (STP), interleave SP with iterative hard thresholding; such hybrids admit weaker RIP requirements (e.g., $\delta_{3s} < 0.5340$ for STP with $\mu=1$) and deliver improved empirical phase transitions at low measurement rates (Song et al., 2013).
• In the block/group setting, GPSP is distinguished by SPC-based expansion and RMC-based pruning, offering superior identification especially in heterogeneous or noisy regimes (He et al., 2024).

Table: RIP and Guarantee Comparison

Algorithm	Required RIP	Guarantee Type
SP	$\delta_{3K} < 0.4859$	Exact recovery / near-oracle noise MSE
CoSaMP	$\delta_{4K} \leq 0.1$	Similar as SP, higher $\delta_{*}$
OMP	$\delta_{K+1} < c$ (stronger)	Requires smaller $\delta_{*}$
STP	$\delta_{3s} < 0.5340$ (with $\mu=1$)	Improved over SP/CoSaMP
GPSP (blocks)	$\delta_{M,2k}<0.1188$	Block-sparse, BRIP-based

7. Practical Applications and Empirical Observations

SP is widely used in compressive sensing, including signal and image recovery, face recognition, PDE system identification, and decentralized sensor networks (0803.0811, He et al., 2024, Li et al., 2014). Empirical studies demonstrate:

• SP exhibits critical sparsity thresholds matching or exceeding those of $\ell_1$-minimization, with far fewer computations (0803.0811, Song et al., 2013).
• In joint-sparse and block settings, GPSP and group extensions outperform classical block-OMP/CoSaMP, especially for heterogeneous, noisy, or underdetermined signals (He et al., 2024).
• Communication-efficient decentralized variants (DCSP/GDCSP) achieve high support recovery accuracy with $O(K^2)$ to $O(KN)$ per-node messages per iteration in typical sensor networks (Li et al., 2014).

Robustness to noise, absence of coefficient-magnitude-dependent performance, and scalability to distributed and structure-enforcing settings have established SP and its derivatives as cornerstones in modern sparse approximation theory and applications.