This algorithm power hour is devoted to Subspace Pursuit, presented in the paper: W. Dai and O. Milenkovic, “Subspace pursuit for compressive sensing signal reconstruction,” IEEE Trans. Inf. Theory, vol. 55, pp. 2230-2249, May 2009.

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\}\).

SP 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\),

SP redefines this set \(\Omega_{k+1} = J_s^*\)

where \(J_s^* \subset \Omega\) is a set of at most \(s\) indices found in the following way.

First, SP finds

$$

J := \textrm{Supp}\{T_1 \left ( |\mathbf{\Phi}^*\vr_k| ; \tau_{k,1} \right ) \}

$$

where \(T_1(\vy)\) is a thresholding function applied to each element of \(\vy\)

according to the threshold \(\tau_{k,1} \ge 0\),

and \(\vr_k\) is the orthogonal projection of \(\vx\) onto the range of \(\mathbf{\Phi}_{\Omega_k}\).

(Wherever an element of \(\vy\) is below \(\tau_{k,1}\) it is set to zero.)

SP sets \(\tau_{k,1}\) such that \(|J| \le s\).

Then with \(\Omega_{k}\cup J\),

SP computes the weights such that \(||\vu – \mathbf{\Phi}\MI_{\Omega_{k}\cup J}\vb||_2\) is minimized,

where \(\MI_{\Omega_{k}\cup J}\) is a \(N\times N\) diagonal matrix

with ones at \((j,j)\), \(j \in \Omega_{k}\cup J\), and zeros elsewhere.

SP then thresholds \(\vb\) to find the new optimal support

$$

J_s^* = \textrm{Supp}\{T_2\left ( |\vb|; \tau_{k,2} \right )\}

$$

where \(T_2(\vy)\) is another thresholding function,

and \(\tau_{k,2} \ge 0\) is set such that a maximum of \(s\) elements of \(\vb\) are retained.

SP repeats the procedure above until some stopping condition,

for instance, \(||\vr_k|| >||\vr_{k-1}|| \), in which case \(\Omega_{k+1} \leftarrow \Omega_{k}\).

The similarity of SP to CoSaMP is unmistakable. In fact, it is so easy to take the correct code for CoSaMP, and make it SP: subspacepursuit.m. *(Note, this code has been corrected here.)*

Below we see a comparison of SP and CoSaMP.

The red lines are CoSaMP.

I am surprised that such a large difference between the two comes from selecting \(s\) elements instead of \(2s\) in the first threshold stage.

(By the way, that nasty hump in CoSaMP probably comes from the fact that I am averaging only 50 realizations. :)

Anyhow, these results agree with those of Maleki and Donoho

where they find SP outperforms CoSaMP.