Algorithm Power Hour: Stagewise Orthogonal Matching Pursuit (StOMP)

This algorithm power hour is devoted to Stagewise Orthogonal Matching Pursuit (StOMP), presented in the paper: D. L. Donoho, Y. Tsaig, I. Drori, and J.-L. Starck, “Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit,” IEEE Trans. on Info. Theory, 2006, but apparently not accepted or published elsewhere. Report here.

Given an $$s$$-sparse vector $$\vx \in \mathcal{R}^N$$,
we sense it by a matrix $$\mathbf{\Phi} \in \mathcal{R}^{m\times N}$$ with $$N > m$$,
which produces the measurement vector $$\vu = \mathbf{\Phi}\vx$$.
Define the index of the columns of $$\mathbf{\Phi}$$ as $$\Omega = \{1, 2, \ldots, N\}$$.
StOMP proceeds in the following way.
Given the index set $$\Omega_k \subset \Omega$$ at iteration $$k$$,
and defining $$\mathbf{\Phi}_{\Omega_k}$$ as the columns from $$\mathbf{\Phi}$$ indexed by $$\Omega_k$$,
StOMP augments this set $$\Omega_{k+1} = \Omega_k \cup J^*$$
where $$J^* \subset \Omega$$ is a set of indices that collectively satisfy an adaptive threshold — just like ROMP.
StOMP defines this set as
$$J^* := \{i \in \Omega : |[\mathbf{\Phi}^*\vr_k]_i| > t_k\sigma_k\}$$
where $$t_k\sigma_k$$ is a threshold parameter, and $$\vr_k$$ is the orthogonal projection of $$\vx$$ onto the range of $$\mathbf{\Phi}_{\Omega_k}$$.
In their work, Donoho et al. define $$\sigma_k = ||\vr_k||_2/\sqrt{m}$$, and $$2\le t_k \le 3$$, motivated from either avoiding false alarms (false alarm rate, which requires knowing the sparsity of the signal) or missed detection (false discovery rate) using concepts from multiple access interference in digital communications. (Very nice! Maleki and Donoho use the same approach to tune their recommended two-stage thresholding algorithm.)
StOMP then orthogonally projects the signal onto the expanded set of atoms to form the new residual, and repeats the procedure above until some stopping condition.
Donoho et al. suggests running the above procedure 10 times.

In their code included in SparseLab, I see a threshold parameter, which I assume is $$t_k$$, but this parameter is defaulted to $$0.5$$. I wonder if that is $$1/t_k$$? I find that when I make it 2, nothing works. When I make it $$0.1$$, results look awful. And when I make it $$0.9$$, or $$0.33$$, it doesn’t look nearly as good as with $$0.5$$. So I will assume $$0.5$$ is where it’s at, technically speaking.