OMP-MMV Algorithm for Joint Sparse Recovery

• OMP-MMV is a greedy algorithm that extends classic OMP to recover jointly sparse solutions from multiple measurement vectors or system matrices.
• It employs tailored support selection rules and residual updates to enforce joint or simultaneous sparsity, ensuring accurate recovery under specific RIP conditions.
• Practical applications such as MRI RF pulse design and channel estimation benefit from its low per-iteration complexity, balancing recovery performance with computational efficiency.

The Orthogonal Matching Pursuit Multiple Measurement Vectors (OMP-MMV) algorithm refers to a class of greedy, iterative algorithms designed to recover jointly sparse solutions to sets of underdetermined linear systems where multiple observations (measurement vectors) or, more generally, multiple system matrices are available. In the canonical MMV model, several measurement vectors share the same sparsity profile, while in advanced formulations such as the multiple-system single-output (MSSO) scenario, multiple unknown vectors, each subject to a different system matrix, combine to form a single observation. OMP-MMV extends the classic OMP algorithm to these richer problem settings, with adaptations to the support selection rule and the definition of residuals to enforce either joint sparsity or simultaneous sparsity constraints.

1. Mathematical Framework and Problem Definitions

Let $d \in \mathbb{C}^M$ denote an observation vector, and let $\{F_1, \ldots, F_P\}$ ($F_p \in \mathbb{C}^{M \times N}$) be system matrices. The MSSO model is

$d = F_1 g_1 + F_2 g_2 + \cdots + F_P g_P,$

where $g_p \in \mathbb{C}^N$ are unknowns. The objective is to recover $\{g_p\}$ such that the collection is simultaneously sparse: the support (locations of non-zeros) is shared across $p$.

Equivalently, stacking the $g_p$ as $g_{\text{tot}} \in \mathbb{C}^{PN}$, and defining $F_{\text{tot}} = [F_1 \ F_2 \ \cdots \ F_P]$, yields $d = F_{\text{tot}} g_{\text{tot}}$.

In traditional MMV,

$Y = A X + E,$

with $Y \in \mathbb{C}^{M \times L}$, $A \in \mathbb{C}^{M \times N}$, $X \in \mathbb{C}^{N \times L}$, and $E$ representing noise. $X$ is row-sparse: non-zeros are concentrated in the same rows.

In the generalized MMV (GMMV) framework, measurement matrices may differ,

$$y^{(i)} = A^{(i)} x^{(i)} + e^{(i)}, \qquad i = 1, \dots, d,$$

joint support is still required across $\{x^{(i)}\}$.

2. Core OMP-MMV Algorithmic Procedure

Adaptations of OMP to the MMV and MSSO settings replace traditional vector–atom correlation with a criterion that aggregates evidence across system matrices or measurement channels in a manner that respects joint or simultaneous sparsity.

MSSO Formulation:

• Define $$C_n = [f_{1,n}, f_{2,n}, \ldots, f_{P,n}] \in \mathbb{C}^{M \times P}$$, collecting the $n$th column across all system matrices.
• Group unknown coefficients as $$h_n = [g_1[n], g_2[n], \ldots, g_P[n]]^T \in \mathbb{C}^P$$.
• Observations are modeled as

$d = \sum_{n=1}^N C_n h_n,$

with a constraint that only $T$ of the $h_n$ are nonzero (the simultaneous sparsity pattern).

• At each OMP iteration $k$, select

$$q_k = \arg\max_{n \notin I_{k-1}} \; r_{k-1}^H Q_n r_{k-1}, \qquad Q_n = C_n C_n^\dagger,$$

where $I_{k-1}$ denotes previously selected indices and $r_{k-1}$ is the current residual.

• Update the residual via a least-squares projection onto the span of all selected $C_n$:

$r_k = d - S_k(S_k^\dagger d),$

with $S_k = [C_{q_1}, \ldots, C_{q_k}]$.

MMV/S-OMP Setting:

• At each iteration, select index $j^*$ via

$$j^* = \arg \max_j \sum_{\ell=1}^L | \langle a_j, r^{(\ell)} \rangle |,$$

where $r^{(\ell)}$ is the residual for measurement $\ell$.

• The support set is updated for all vectors simultaneously; then the least-squares problem is solved jointly.

Generalized MMV (MOMP):

• When measurement matrices vary, selection is by maximizing

$$l^* = \arg\max_l \sum_{i=0}^{d-1} \left| \langle a^{(i)}_l, r^{(i)}_{k-1} \rangle \right|^2.$$

• Residuals and supports are maintained for each measurement vector independently but updated for the same index at every iteration.

3. Theoretical Guarantees and Recovery Conditions

The exact recovery guarantees of OMP-MMV depend on matrix properties, sparsity, and noise:

• Restricted Isometry Property (RIP): For matrix $A$ (or generalizations for MSSO), RIP of order $k+1$ with sufficiently small $\delta_{k+1}$ (for example, $\delta_{k+1} < 1/(2\sqrt{k} + 1)$ in the MMV case) ensures that OMP-MMV recovers the exact support in the noiseless setting (Ding et al., 2011).
• Noisy Recovery: Under general measurement and sensing matrix perturbations, support recovery is still guaranteed if the minimal norm of the nonzero rows $t_0$ greatly exceeds the effective noise/perturbation level $\varepsilon_h$, with $\delta_{k+1} < Q(k, t_0/\varepsilon_h)$ for $Q$ a decreasing function (Ding et al., 2011). The relative reconstruction error scales almost linearly with noise levels.
• Instance Optimality: OMP exhibits $(p,q)$ instance optimality for $1 \leq p \leq q \leq 2$ (except the deterministic $(2,2)$ case) under suitable RIP conditions. In the probabilistic case, $(2,2)$ optimality is achieved for random matrices with high probability (Xu, 2010).
• Generalization to Multiple System Matrices: In the MSSO setting, the choice rule and residual update are fundamentally the same as described above. However, the increased block structure complexity results in only intermediate recovery performance compared to more computationally intensive algorithms (e.g., LSMP) or convex approaches (e.g., IRLS, SOCP) (0907.2083).

4. Experimental Behavior and Computational Implications

In the extensive experimental evaluation for the MSSO scenario (0907.2083):

• OMP-MMV consistently yields support recovery performance between that of basic (single-step) matching pursuit and the more sophisticated least-squares matching pursuit.
• In MRI RF pulse design—an application with stringent image fidelity requirements—OMP-MMV does not match the performance of LSMP or convex relaxations.
• OMP-MMV's main operational advantages are low per-iteration complexity (relative to convex programs) and avoidance of repeated index selection. However, it requires pseudoinversion of an expanding submatrix at every iteration, which increases computational burden compared to MP.
• Greedy methods (MP, OMP, LSMP) are dominated in support recovery accuracy by methods based on SOCP or IRLS in both noiseless and noisy cases.

5. Advanced Variants and Extensions

Alternative OMP-MMV algorithms are employed for specific scenarios:

• Weighted OMP-MMV (SOMP-NS): Incorporates noise stabilization by weighting each measurement vector's influence proportional to its estimated inverse variance, improving robustness in the presence of heterogeneous noise (Determe et al., 2015).
• Generalized MMV (GMMV): Measurement matrices allowed to differ across observations; MOMP and related algorithms aggregate support selection across measurements. Recovery probabilities decay exponentially in the number of measurements under mild average isometry conditions—diversity in sensing matrices may improve performance (Heckel et al., 2012).
• Structural MMV OMP (GM-OMP): Enables the recovery of structured sparsity patterns in $X$ beyond strict joint support, using feasible sets specified by graphical connectivity or continuity constraints, facilitating structurally faithful reconstructions in high-dimensional applications (Boßmann, 2017).
• Multidimensional/Sparse Tensor OMP: Extension to high-dimensional dictionaries assembled as the product of smaller ones, with algorithmic modifications to exploit separability and reduce computation/memory requirements (MOMP, SMOMP) for applications such as mmWave channel estimation and localization (Palacios et al., 2022, Palacios et al., 2022).

6. Practical Applications and Limitations

OMP-MMV algorithms, including their structural and multidimensional variants, have been used in magnetic resonance RF pulse design (0907.2083), compressed sensing, direction-of-arrival estimation, neuromagnetic source localization, and mmWave joint channel and position estimation (Palacios et al., 2022). The key strengths are efficient greedy support identification in large-scale problems and adaptability to measurement and system matrix structure. The main trade-offs are:

• Computation: OMP-MMV typically entails solving growing (block) least-squares problems, increasing per-iteration cost compared to MP but remaining manageable versus convex solvers; multidimensional variants (MOMP, SMOMP) provide major complexity reductions in high tensor-dimensional problems.
• Recovery Performance: Performance is reliable under appropriate RIP conditions and with sufficiently strong signals relative to the effective noise floor; however, for applications demanding precise support identification or high-fidelity reconstructions in noise, convex approaches or LSMP can be superior.
• Limitation in Simultaneous Sparsity: While OMP-MMV stabilizes support selection, it may not always maximize recovery fidelity in simultaneous sparsity models, particularly when the mixture structure imposes strong dependencies among system matrices or measurement vectors.

7. Summary of Key Formulas and Algorithmic Steps

MSSO OMP-MMV Key Steps:

1. For $k = 1, \dots, K$:
   • Select

   $$q_k = \arg\max_{n \notin I_{k-1}} \left\{ r_{k-1}^H (C_n C_n^\dagger) r_{k-1} \right\}$$

• Update support $I_k = I_{k-1} \cup \{q_k\}$
• Aggregate $S_k = [C_{q_1}, \ldots, C_{q_k}]$
• Residual update:

  $r_k = d - S_k (S_k^\dagger d)$

1. Terminate upon reaching prescribed sparsity or convergence.

MMV S-OMP Atom Selection: $$j^* = \arg \max_j \sum_{\ell=1}^L | \langle a_j, r^{(\ell)} \rangle |$$

GMMV MOMP Selection: $$l_p = \arg \max_l \sum_{i=0}^{d-1} \left| \langle a_l^{(i)}, r_{p-1}^{(i)} \rangle \right|^2$$

Recovery Condition (MMV, noiseless) (Ding et al., 2011):

If $\Phi$ has RIP constant $\delta_{k+1} < 1/(2\sqrt{k}+1)$, OMP-MMV recovers the exact support in $k$ iterations.

---

In conclusion, OMP-MMV and its variants generalize classical greedy algorithms to hybrid, structurally rich, and multi-measurement settings. Their success depends on appropriate adaptation of support selection and residual projections to enforce joint or simultaneous sparsity, their compliance with matrix-theoretic recovery conditions (RIP/RIC), and careful management of computational complexity, especially as problem dimensionality or system matrix diversity increase. The algorithm is widely used across signal processing domains but is best suited for scenarios balancing the need for efficient, scalable support recovery and relaxed requirements on noise robustness and optimal recovery accuracy (0907.2083, Ding et al., 2011, Heckel et al., 2012, Palacios et al., 2022, Palacios et al., 2022).