Coherence-Pattern Guided Compressive Sensing with Unresolved Grids

Introduction

Reconstruction of a high-dimensional sparse signal from sparse linear measurements is a fundamental problem relevant to imaging, inverse problems and signal processing.

Consider, for example, the problem of spectral estimation in signal processing. Let the uncontaminated signal $y(t)$ be a linear combinations of $s$ time-harmonic components

where $c_{j}$ are the amplitudes. Suppose that $y(t)$ is contaiminated by noise $n(t)$ and the received signal is

The task is to find out the frequencies $\Omega=\{\omega_{j}\}$ and the amplitudes $\{c_{j}\}$ by sampling $b(t)$ at discrete times.

We cast the spectral estimation problem into the form

Small gridding error requires that the objects are a priori close to the grid points. The gridding error is related to basis mismatch analyzed in .

Sparse reconstruction with $N,s\ll M$, where $s$ is the sparsity of $\mathbf{x}$, has recently attracted a lot of attention in various areas thanks to the breakthroughs in compressive sensing (CS) . The main thrust of CS is the $L^{1}$-minimization principle, Basis Pursuit (BP) and Lasso, for solution characterization. Many $L^{1}$-based algorithms as well as the alternative, greedy algorithms, which are not directly based on global optimization, require either incoherence or Restricted Isometry Property (RIP) to have good performances.

One commonly used characterization of incoherence in CS is in terms of the mutual coherence $\mu$. Let the pairwise coherence between the $k$-th and $j$-th columns be

The mutual coherence of ${\mathbf{A}}$ is the maximum pairwise coherence among all pairs of columns

Without any prior information about the object support, the gridding error for the resolved grid, however, can be as large as the data themselves, creating a unfavorable condition for sparse reconstruction. To reduce the gridding error, it is natural to consider the fractional grid

Figure 2 shows the coherence pattern $[\mu(j,k)]$ of a $100\times 4000$ matrix (4) with $F=20$ (left panel). The bright diagonal band represents a heightened correlation (pairwise coherence) between a column vector and its neighbors on both sides (about 30). The right panel of Figure 2 shows a half cross section of the coherence band across two RLs. Sparse recovery with large $F$ exceeds the capability of currently known algorithms as the condition number of the $100\times 30$ submatrix corresponding to the coherence band in Figure 2 easily exceeds $10^{15}$. The high condition number makes stable recovery impossible.

The difficulty with unresolved grids is not limited to the problem of spectral estimation in signal processing. Indeed, the issue is intrinsic and fundamental to discretization of PDE-based inverse problems such as remote sensing and medical imaging . While Figure 2 is typical of the coherence pattern from discretization of one-dimensional problem. In two or three dimensions, the coherent pattern is more complicated than Figure 2. Nevertheless the coherence band typically reflects proximity in the physical space. The proximity between the object support and its reconstruction can be described by the Bottleneck or the Hausdorff distance . More generally, coherent bands can arise in sparse and redundant representation by overcomplete dictionaries (see Section 6 for an example). Under this circumstance, the Bottleneck or Hausdorff distance may not have a direct physical meaning.

In any case, the hope is that if the objects are sufficiently separated with respect to the coherence band, then the problem of a huge condition number associated with unresolved grids can be somehow circumvented and the object support can be approximately reconstructed.

Under this additional assumption of widely separated objects, we propose in the present work several algorithmic approaches to recovery with unresolved grids and provide some performance guarantee for these algorithms.

The paper is organized as follows. In Section 2 we introduce the technique of band exclusion (BE) to modify the Orthogonal Matching Pursuit (OMP) and obtain a performance guarantee for the improved algorithm, called Band-excluded OMP (BOMP). In Section 3 we introduce the technique of Local Optimization (LO) and propose the algorithms, Locally Optimized OMP (LOOMP) and Band-excluded LOOMP (BLOOMP). In Section 4 we introduce the technique of Band-Excluded Thresholding (BET), which comes in two form, Band-excluded Matched Thresholding (BMT) and Band-excluded Locally Optimized Thresholding (BLOT) and propose the algorithms, Band-excluded, Locally Optimized Subspace Pursuit (BLOSP), Band-excluded, Locally Optimized CoSaMP (BLOCoSaMP), Band-excluded, Locally Optimized Iterative Hard Thresholding (BLOIHT) and BP/Lasso with BLOT (BP/Lasso-BLOT). In Section 5 we present numerical study of the comparative advantages of various algorithms. In Section 6 we compare the performance of our algorithms with the existing algorithms, Spectral Iterative Hard Thresholding (SIHT) and coherent-dictionary-based BP recently proposed in and , respectively. We conclude in Section 7.

Band Exclusion (BE)

The first technique that we introduce to take advantage of the prior information of widely separated objects is called Band Exclusion and can be easily embedded in the greedy algorithm, Orthogonal Matching Pursuit (OMP) .

First let us recall a standard performance guarantee for OMP .

Suppose that the sparsity $s$ of the signal vector $\mathbf{x}$ satisfies

where $x_{\text{min}}=\displaystyle\min_{k}|x_{k}|=|x_{s}|$. Denote by $\hat{\mathbf{x}}$, the output of the OMP reconstruction. Then

where $\text{supp}(\mathbf{x})$ is the support of $\mathbf{x}$.

In the ideal case where ${\mathbf{e}}=0$, (10) reduces to

which is near the threshold of OMP’s capability for exact reconstruction of arbitrary objects of sparsity $s$.

Intuitively speaking, if the objects are not in each other’s coherence band, then it should be possible to localize the objects approximately within their respective coherence bands, no matter how large the mutual coherence is.

Let us first define precisely the notion of coherence band. Let $\eta>0$. Define the $\eta$-coherence band of the index $k$ to be the set

and the $\eta$-coherence band of the index set $S$ to be the set

Due to the symmetry $\mu(i,k)=\mu(k,i),\forall i,k$, $i\in B_{\eta}(k)$ if and only if $k\in B_{\eta}(i)$.

To imbed BE into OMP, we make the following change to the matching step

meaning that the double $\eta$-band of the estimated support in the previous iteration is avoided in the current search. This is natural if the sparsity pattern of the object is such that $B_{\eta}(j),j\in\hbox{supp}(\mathbf{x})$ are pairwise disjoint. We call the modified algorithm the Band-excluded Orthogonal Matching Pursuit (BOMP) which is formally stated in Algorithm 1.

A main theoretical result of the present paper is the following performance guarantee for BOMP.

Let $\mathbf{x}$ be $s$-sparse. Let $\eta>0$ be fixed. Suppose that

Let $\hat{\mathbf{x}}$ be the BOMP reconstruction. Then $\hbox{supp}(\hat{\mathbf{x}})\subseteq B_{\eta}(\hbox{supp}(\mathbf{x}))$ and moreover every nonzero component of $\hat{\mathbf{x}}$ is in the $\eta$-coherence band of a unique nonzero component of $\mathbf{x}$.

Suppose $\hbox{supp}(\mathbf{x})=\{J_{1},\ldots,J_{s}\}$. Let $J_{\rm max}\in\hbox{supp}(\mathbf{x})$ be the index of the largest component in absolute value of $\mathbf{x}$.

by assumption (15). On the other hand, $\forall l\notin B_{\eta}(\hbox{supp}(\mathbf{x}))$,

then the right hand side of (2) is greater than the right hand side of (2) which implies that the first index selected by BOMP must belong to $B_{\eta}(\hbox{supp}(\mathbf{x}))$.

Now suppose without loss of generality that the first $(k-1)$ indices $I_{1},...,I_{k-1}$ selected by BOMP are in $B_{\eta}(J_{i}),J_{i}\in\hbox{supp}(\mathbf{x}),i=1,...,k-1$, respectively. Write the residual as

First, we estimate the coefficients $c_{I_{1}},...,c_{I_{k-1}}$. Since $\left\langle\mathbf{r}^{k-1},{\mathbf{a}}_{I_{1}}\right\rangle=0$,

Let $c_{\text{max}}=\displaystyle\max_{j=1,...,k-1}|c_{I_{j}}|$. Inequality (19) implies that

Moreover, condition (16) implies that $\eta(s-1)<\frac{1}{5}$ and $\frac{1}{1-\eta(k-2)}\leq\frac{5}{4}$. Hence

We claim that $B^{(2)}_{\eta}(S^{k-1})$ and $\{J_{k},...,J_{s}\}$ are disjoint.

If the claim is not true, then there exists $J_{i},$ for some $i\in\{k,\ldots,s\}$, $J_{i}\in B^{(2)}_{\eta}(I_{l})$ for some $l\in\{1,...,k-1\}$. Consequently, $J_{i}\in B^{(3)}_{\eta}(J_{l})$ or equivalently $B_{\eta}(J_{i})\cap B^{(2)}_{\eta}(J_{l})\neq\emptyset$ which is contradictory to the assumption (15).

Now we show that the index selected in the $k$-th step is in $B_{\eta}(\{J_{k},...,J_{s}\})$.

One the other hand, we have that $\forall l\notin B^{(2)}_{\eta}(S^{k-1})\cup B_{\eta}(\{J_{k},...,J_{s}\})$,

in view of $B_{\eta}(\{J_{1},\ldots,J_{k-1}\})\subseteq B^{(2)}_{\eta}(S^{k-1})$ as a result of the induction assumption.

If the right hand side of (2) is greater than the right hand side of (2) or equivalently

then the $k$-th index selected by BOMP must be in $B_{\eta}(\{J_{k},...,J_{s}\})$ because

and because the $k$-th selected index does not belong in $B^{(2)}_{\eta}(S^{k-1})$ according to the band-exclusion rule. Condition (16) implies (23) by setting the maximal $k=s$ in (23) and noting that $\eta(5s-4)<1$ under (16). ∎

In the case of the matrix (4), if every two indices in ${\hbox{supp}(\mathbf{x})}$ is more than one RL apart, then $\eta$ is small for sufficiently large $N$, cf. Figure 2.

When the dynamic range ${x_{\text{max}}}/{x_{\text{min}}}=\mathcal{O}(1)$, Theorem 1 guarantees approximate recovery of $\mathcal{O}(\eta^{-1})$ sparsity pattern by BOMP.

The main difference between Theorem 1 and Proposition 1 lies in the role played by the dynamic range $x_{\text{max}}/x_{\text{min}}$ and condition (15).

First, numerical evidence points to degradation in BOMP’s performance for large dynamic ranges (Figure 3). This is consistent with the prediction of (16).

Secondly, condition (15) means that BOMP can resolve 3 RLs. Numerical experiments show that BOMP can resolve objects separated by close to 1 RL when the dynamic range is close to 1 (Figure 7).

Local Optimization (LO)

As our numerical experiments show, the main shortcoming with BOMP is in its failure to perform even when the dynamic range is only moderate.

To overcome this problem, we now introduce the second technique: the Local Optimization (LO).

LO is a residual-reduction technique applied to the current estimate $S^{k}$ of the object support. To this end, we minimize the residual ${\|{\mathbf{A}}\hat{\mathbf{x}}-\mathbf{b}\|_{2}}$ by varying one location at a time while all other locations held fixed. In each step we consider $\hat{\mathbf{x}}$ whose support differs from $S^{n}$ by at most one index in the coherence band of $S^{n}$ but whose amplitude is chosen to minimize the residual. The search is local in the sense that during the search in the coherence band of one nonzero component the locations of other nonzero components are fixed. The amplitudes of the improved estimate is carried out by solving the least squares problem. Because of the local nature of the LO step, the computation is not expensive.

Embedding LO in BOMP gives rise to the Band-excluded, Locally Optimized Orthogonal Matching Pursuit (BLOOMP).

We now give a condition under which LO does not spoil the BOMP reconstruction of Theorem 1.

Let $\eta>0$ and let $\mathbf{x}$ be a $s$-sparse vector such that (15) holds. Let $S^{0}$ and $S^{k}$ be the input and output, respectively, of the LO algorithm.

and each element of $S^{0}$ is in the $\eta$-coherence band of a unique nonzero component of $\mathbf{x}$, then each element of $S^{k}$ remains in the $\eta$-coherence band of a unique nonzero component of $\mathbf{x}$.

Because the iterative nature of Algorithm 2, it is sufficient to show that each element of $S^{1}$ is in the $\eta$-coherence band of a unique nonzero component of $\mathbf{x}$.

Suppose $J_{1}\in{\hbox{\rm supp}}{(\mathbf{x})}$ and $i_{1}\in B_{\eta}(J_{1})$. Let

We want to show that $r<r^{\prime},\forall j\in B_{\eta}(i_{1})\backslash B_{\eta}(J_{1})$ so that the LO step is certain to pick a new index within the $\eta$-coherence band of $J_{1}$, reducing the residual in the meantime. For the subsequent analysis, we fix $j\in B_{\eta}(i_{1})\backslash B_{\eta}(J_{1})$.

Reset the $J_{1}$ component of $\mathbf{x}$ to zero and denote the resulting vector by $\mathbf{x}^{\prime}$. Hence the sparsity of $\mathbf{x}^{\prime}$ is $s-1$. It follows from the definition of $r$ that

where ${\hbox{\rm supp}}({\mathbf{z}})=\{i_{2},\ldots,i_{k}\}$.

Because of (15), $j,J_{1}\not\in B_{\eta}({\hbox{\rm supp}}(\mathbf{x}^{\prime})\cup\{i_{2},\ldots,i_{k}\})$. By the definition of $\eta$-coherence band, we have

To prove $r<r^{\prime}$, it suffices to show

Considering the worst case scenario, we replace $|x_{J_{1}}|$ by $x_{\rm min}$ and $k$ by $s$ to obtain the condition (24).

Let $\hat{\mathbf{x}}$ be the output of BLOOPM. Under the assumptions of Theorems 1 and 2, $\hbox{supp}(\hat{\mathbf{x}})\subseteq B_{\eta}(\hbox{supp}(\mathbf{x}))$ and moreover every nonzero component of $\hat{\mathbf{x}}$ is in the $\eta$-coherence band of a unique nonzero component of $\mathbf{x}$.

Even though we can not improve the performance guarantee for BLOOMP, in practice the LO technique greatly enhances the success probability of recovery that BLOOMP has the best performance among all the algorithms tested with respect to noise stability and dynamic range (see Section 5). In particular, the LO step greatly enhances the performance of BOMP w.r.t. dynamic range.

Band-Excluded Thresholding (BET)

The BE technique can be extended and applied to selecting $s$ objects all at once in what is called the Band-Excluded Thresholding (BET).

We consider two forms of BET. The first is the Band-excluded Matched Thresholding (BMT) which is the band-exclusion version of the One-Step Thresholding (OST) recently shown to possess compressed-sensing capability under incoherence conditions .

For the purpose of comparison with BOMP, we give a performance guarantee for BMT under similar but weaker conditions than (15)-(16).

Let $\mathbf{x}$ be $s$-sparse. Let $\eta>0$ be fixed. Suppose that

Let $\hat{\mathbf{x}}$ be the BMT reconstruction. Then $\hbox{supp}(\hat{\mathbf{x}})\subseteq B_{\eta}(\hbox{supp}(\mathbf{x}))$ and moreover every nonzero component of $\hat{\mathbf{x}}$ is in the $\eta$-coherence band of a unique nonzero component of $\mathbf{x}$.

Let $\hbox{supp}(\mathbf{x})=\{J_{1},\ldots,J_{s}\}$. Let $J_{\rm max}\in\hbox{supp}(\mathbf{x})$ be the index of the largest component of $\mathbf{x}$ in absolute value.

On the other hand, $\forall l\notin B_{\eta}({\hbox{supp}(\mathbf{x})})$,

Therefore, the condition (27) implies that the right hand side of (4) is greater than the right hand side of (4). This means $|\mathbf{b}^{*}{\mathbf{a}}_{k}|>|\mathbf{b}^{*}{\mathbf{a}}_{l}|,\forall k=1,..,s,\forall l\not\in B_{\eta}({\hbox{supp}(\mathbf{x})})$ and hence the $s$ highest points of $|\mathbf{b}^{*}{\mathbf{a}}_{k}|$ are in $B_{\eta}({\hbox{supp}(\mathbf{x})})$.

From step ii) of BMT and (26) it follows the second half of the statement, namely every nonzero component of $\hat{\mathbf{x}}$ is in the $\eta$-coherence band of a unique nonzero component of $\mathbf{x}$. ∎

Condition (26) roughly means that the objects are separated by at two RLs which is weaker than (15). In numerical simulations, however, BOMP performs far better than BMT. In other words, BMT is not a stand-alone algorithm but should instead be imbedded in other algorithms such as Subspace Pursuit (SP) , the Compressive Sampling Matching Pursuit (CoSaMP) and the Normalized Iterative Hard Thresholding (IHT) . This gives rise to Band-excluded Subspace Pursuit (BSP), Band-excluded Compressive Sampling Matching Pursuit (BCoSaMP) and Band-excluded Normalized Iterative Hard Thresholding (BNIHT) which we demonstrate their performance numerically.

In addition to BMT, the second form of BET, namely the Band-excluded, Locally Optimized Thresholding (BLOT), can further enhance the performance in reconstruction with unresolved grids.

Now we state the algorithm Band-excluded, Locally Optimized Subspace Pursuit (BLOSP). The Band-excluded Locally Optimized CoSaMP (BLOCoSaMP) is similar and omitted here.

Embedding BLOT in NIHT turns out to have a nearly identical performance to embedding BLOT in the Iterative Hard Thresholding (IHT) . Since the latter is simpler to implement and more efficient to compute, we state the resulting algorithm, the Band-excluded, Locally Optimized IHT (BLOIHT), below.

In addition, the technique BLOT can be used to enhance the recovery capability with unresolved grids of the $L^{1}$-minimization principles, Basis Pursuit (BP)

where $\sigma$ is the standard deviation of the each noise component and $\lambda$ is the regularization parameter. In this case, BLOT is applied to the BP and Lasso estimates $\hat{\mathbf{x}}$ to produce a $s$-sparse reconstruction. The resulting algorithms are called BP-BLOT and Lasso-BLOT, respectively.

The thresholded Lasso, the Lasso followed by a hard thresholding, has been considered previously (see and references therein). The novelty of our version lies in the BE and LO steps which greatly enhance the performance in dealing with unresolved grids.

Numerical study

We test our various band-exclusion algorithms on the matrix

where $\xi_{k}$ is uniformly and independently distributed in $(0,1)$. When $F=1$, ${\mathbf{A}}$ is the random partial Fourier matrix analyzed in and, with sufficient number of samples, successful recovery with ${\mathbf{A}}$ is guaranteed with high probability. For large $F$, however, ${\mathbf{A}}$ has a high mutual coherence. Unless otherwise stated, we use $N=100$, $M=4000$ and $F=20$ and the i.i.d. Gaussian noise ${\mathbf{e}}\sim N(0,\sigma^{2}I)$ in our simulations. Recall the decay profile of pairwise coherence as the index separation increases in Figure 2 (right). The $\eta$-coherence band is about $0.7$ RL in half width with $\eta=0.3$.

We compare performance of various algorithms in terms of success probability versus dynamic range, noise level, number of measurements and resolution. Unless otherwise stated, we use in our simulations 10 randomly phased and located objects, separated by at least 3 RLs. A reconstruction is counted as a success if every reconstructed object is within 1 RL of the object support. This is equivalent to the criterion that the Bottleneck distance between the true support and the reconstructed support is less than 1 RL.

For subsets in one dimension, the Bottleneck distance can be calculated easily. Let $A=\{a_{1},\ldots,a_{n}\}$ and $B=\{b_{1},\ldots,b_{n}\}$ be listed in the ascending order. Then

In higher dimensions, however, it is more costly to compute the Bottleneck distance . The Bottleneck distance is a stricter metric than the Hausdorff distance which does not require one-to-one correspondence between the two target sets.

In the first set of experiments, we test various greedy algorithms equipped with the BE step (only). This includes BOMP, Band-excluded Subspace Pursuit (BSP), Band-excluded CoSaMP (BCoSaMP) and Band excluded Normalized Iterative Hard Thresholding (BNIHT). For comparison, we also show the performance of OMP without BE.

As shown in Figure 3, BOMP has the best performance with respect to dynamic range, noise and, in the case of higher dynamic range ($\geq 3$, bottom right panel), number of measurements. In the case of dynamic range equal to 1, BSP is the best performer in terms of number of measurements followed closely by BCoSaMP and BNIHT (bottom panel). In the case of dynamic range equal to 1, BOMP and OMP have almost identical performance with respect to noise (middle left) and number of measurements (bottom left). The performance of BSP and BCoSaMP, however, depends crucially on the BE step without which both SP and CoSaMP fail catastrophically (not shown).

In the next set of experiments, we test BLO-equipped algorithms. For the purpose of comparison, we also test the algorithm, Locally Optimized OMP (LOOMP) which is the same as Algorithm 3 but without BE.

Lasso-BLOT is implemented with the regularization parameter

which is proposed in . Other larger values have been proposed in . Our numerical experiments indicate that for matrix (32) with large $F$ the choice (33) is nearly optimal among all $\lambda/\sqrt{\log{M}}\leq 10$ and relative noise up to $5\%$. The superiority of the choice (33) to (34) (and other choices) manifests clearly across all performance figures involving both of them.

FIgure 4 shows success probability versus dynamic range in the presence of noise. The top performers are LOOMP and BLOOMP both of which can handle large dynamic range. In the noiseless case, the success rate for LOOMP, BLOOMP, BLOSP, BLOOMP and BLOCoSaMP stays near unity for dynamic range up to as high as $10^{14}$. With $1\%$ noise (left panel), BLOSP, BLOCoSaMP and BLOIHT perform better than Lasso-BLOT with either (33) or (34) while with $3\%$ noise (right panel), BLOSP, BLOCoSaMP and BLOIHT performance curves have dropped below that of Lasso-BLOT with (33). But the noise stability of Lasso-BLOT never catches up with that of LOOMP/BLOOMP even as the noise level increases as can be seen in Figure 5.

Figure 5 shows that LOOMP and BLOOMP remain the top performers with respect to noise while Lasso-BLOT with (34) has the worst performance. Lasso-BLOT with (33), however, is a close second in noise stability. As seen in Figures 4 and 5, the performance of Lasso-BLOT depends significantly on the choice of the regularization parameter.

With respect to number of measurements (Figure 6), BP/Lasso-BLOT with (33) is the best performer, followed closely by BLOSP and BLOCoSaMP for dynamic range 1 (left panels) while for dynamic range 10, BLOOMP and LOOMP perform significantly better than the rest (right panels). As clear from the comparison of the top left and right panels of Figure 6, at low level of noise the performance of BLOOMP and LOOMP improves significantly as the dynamic range increases from 1 to 10. In the meantime, the performance of BLOSP, BCoSaMP and BLOIHT improves slightly while the performance of BP/Lasso-BLOT deteriorates. At $10\%$ noise, however, the performance of BLOOMP and LOOMP is roughly unchanged as dynamic range increases while the performance of all other algorithms deteriorates significantly (bottom right).

Next we compare the resolution performance of the various algorithms for 10 randomly phased objects of unit dynamic range in the absence of noise. The 10 objects are consecutively located and separated by equal length varying from $0.1$ to 3 RLs. The whole object support is, however, randomly shifted for each of the 100 trials. For closely spaced objects, it is necessary to modify the band exclusion and local optimization rules: If $h$ is the object spacing, we use $h/2$ to replace 2 RLs of the original BE rule and 1 RL of the original LO rule.

Figure 7 shows the averaged Bottleneck distance between the reconstruction and the true object support (top panels) and the residual (bottom panels) as a function of the object spacing. For this class of objects, BP-BLOT has the best resolution for dynamic range up to 10 followed closely by BLOIHT for dynamic range 1 and by BLOOMP/LOOMP for dynamic range 10. The high precision (i.e. nearly zero Bottleneck distance) resolution ranges from about 1.5 RLs for BP-BLOT to about 1.7 RLs for the rest. Consistent with what we have seen in Figure 6, the resolving power of BLOOMP/LOOMP improves significantly as the dynamic range increases while that of BP-BLOT deteriorates. Note that in the case of unit dynamic range, BOMP recovers the support as well as BLOOMP/LOOMP does. BOMP, however, produces a high level of residual error even when objects are widely separated. There is essentially no difference between the resolving power of BLOOMP and LOOMP.

It is noteworthy that the relative residuals of all tested algorithms peak at separation between 1 and 1.5 RLs and decline to zero as the object separation decreases. In contrast, the average Bottleneck distances increase as the separation decreases except for BP-BLOT when the separation drops below 0.5 RL. When the separation drops below 1 RL, the Bottleneck distance between the objects and the reconstruction indicates that the objects are not well recovered by any of the algorithms (top panels). The vanishing residuals in this regime indicates nonuniqueness of sparse solutions.

Figures 4-7 show negligible difference between the performances of LOOMP and BLOOMP. To investigate their distinction more closely, we test the stability with respect to the gridding error for various $F$’s (cf. Figure 1). We consider randomly phased objects of dynamic range 10 that are randomly located in $[0,1000)$ and separated by at least $3$ RLs. We compute the reconstruction error in the Bottleneck distance and $L^{2}$-norm averaged over 100 trials with $5\%$ external noise and show the result in Figure 8. Evidently the advantage of BLOOMP over LOOMP lies in the cases when the refinement factor $F$ is less than 10 and the gridding error is sufficiently large. When $F\geq 10$, the difference between their performances with respect to gridding error is negligible. On the other hand, for $F=5$ both BLOOMP and LOOMP’s reconstructions would have been considered a failure given the magnitudes of error in the Bottleneck distance.

Comparison with other algorithms in the literature

The present work is inspired by the performance guarantee established in that the MUSIC algorithm aided by BMT produces a support estimate that is within 1 RL of the locations of sufficiently separated objects.

In comparison to other CS literature on coherent and redundant dictionary, our work resembles those of and , albeit with a different perspective. In fact, the algorithms developed in and can not be applied to the spectral estimation problem formulated in the Introduction. The algorithms developed here, however, can be applied to their frame-based setting and this is what we will do below for the purpose of comparison.

Following we consider the following problem

where ${\mathbf{\Phi}}$ is a $N\times R$ i.i.d Gaussian matrix of mean and variance $\sigma^{2}$. The signal to be recovered is given by $\mathbf{y}={\mathbf{\Psi}}\mathbf{x}$ where ${\mathbf{\Psi}}$ is the over-sampled, redundant DFT frame

As before $F$ is the refinement factor. Combining (35) and (36) we have the same form (5) with ${\mathbf{A}}={\mathbf{\Phi}}{\mathbf{\Psi}}$ whose coherence pattern is shown in Figure 9 similar to Figure 2.

The algorithm proposed in , Spectral Iterative Hard Thresholding (SIHT), relies on a measurement matrix ${\mathbf{\Phi}}$ satisfying some form of RIP and so does the frame-adapted BP proposed in

which, in addition, requires the analysis coefficients ${\mathbf{\Psi}}^{*}\mathbf{y}$ to be $s$-sparse or $s$-compressible. The RIP is satisfied by i.i.d. Gaussian matrices of proper sizes. Their approaches are not applicable, however, when the sensing matrix ${\mathbf{A}}$ can not be decomposed into a product ${\mathbf{\Phi}}{\mathbf{\Psi}}$ where ${\mathbf{\Phi}}$ has RIP.

In the language of digital signal processing (37) is a $L^{1}$-analysis method while the standard BP or Lasso and the BLOT-version are $L^{1}$-synthesis methods. Both methods are based on the same principle of representational sparsity. In principle, the synthesis approach (such as the BP, Lasso and all BLO-based algorithms) is more general than the analysis approach since every analysis method can be recast as a synthesis method while many synthesis formulations have no equivalent analysis form. In practice, however, each of them performs the best on different types of signals .

For example, the bottom right panel of Figure 10 shows the absolute values of the component of the vector ${\mathbf{\Psi}}^{*}\mathbf{y}$ in the order of descending magnitude. Clearly ${\mathbf{\Psi}}^{*}\mathbf{y}$ is neither sparse nor compressible. The shape of the curve can be roughly understood as follows. The DFT frame ${\mathbf{\Psi}}$ has a coherence band of roughly 1.5 RLs, corresponding to 30 columns, and since the 10 components in $\mathbf{x}$ are widely separated ${\mathbf{\Psi}}^{*}\mathbf{y}$ has about $300$ significant components. The long tail of the curve is due to the fact that the pairwise coherence of ${\mathbf{\Psi}}$ decays slowly as the separation increases. Therefore the analysis approach (37) would require a far higher number of measurements than 100 for accurate reconstruction.

For the analysis approach, the main quantity of interest is $\mathbf{y}$. So in our comparison experiments, we measure the performance by the relative error $\|\hat{\mathbf{y}}-\mathbf{y}\|_{2}/\|\mathbf{y}\|_{2}$ as in . We compute the averaged relative errors as dynamic range, noise level and the number of measurements vary. In each of the 100 trials, 10 randomly phased and located objects (i.e. $\mathbf{x}$) of dynamic range 1 and i.i.d. Gaussian $\Phi$ are generated.

The results are shown in Figure 10. Consistently across the top left, top right and bottom left panels, the smallest error is achieved by BP-BLOT and Lasso-BLOT with (33) with respect to dynamic range (top left), noise (top right) and number of measurements (bottom left). The latter two plots are for dynamic range 1. BLOOMP and BLOSP perform among the best with respect to dynamic range (top left) and noise (top right) and can achieve the minimum error with increasing number of measurements (bottom left). The SIHT algorithm requires much higher number of measurements to get its error down (bottom left) and produces the highest level of error w.r.t. dynamic range (top left) and noise (top right).

In Figure 10, we include the performance curves of OMP with various sparsities (s, 2s, 5s) as well as the standard BP/Lasso without BLOT. Not surprisingly, the relative error of OMP reconstruction decreases with increasing sparsity.

It is noteworthy that the BLOT technique reduces the BP/Lasso reconstruction errors: BP without BLOT produces $0.5\%$ relative error with respect to dynamic range (top left, not clearly visible) while BP-BLOT, along with BLOSP and BLOOMP, produces $10^{-16}$ relative error. Moreover, Lasso with the optimized parameter (33) but without BLOT produces significantly higher errors than Lasso-BLOT, BLOOMP and BLOSP (top right).

Conclusions and Discussions

We have developed and tested various algorithms for sparse recovery with highly coherent sensing matrices arising in discretization of imaging problems in continuum such as radar and medical imaging when the grid spacing is below the Rayleigh threshold .

We have introduced two essential techniques to deal with unresolved grids: band exclusion and local optimization. We have embedded these techniques in various CS algorithms and performed systematic tests on them. When embedded in OMP, both BE and LO steps manifest their advantage in dealing with larger dynamic range. When embedded in SP, CoSaMP, IHT, BP and Lasso the effects are more dramatic.

We have studied these modified algorithms from four performance metrics: dynamic range, noise stability, sparsity and resolution. With respect to the first two metrics (dynamic range and noise stability), BLOOMP is the best performer. With respect to sparsity, BLOOMP is the best performer for high dynamic range while for dynamic range near unity BP-BLOT and Lasso-BLOT with the optimal regularization parameter have the best performance. BP-BLOT also has the highest resolving power up to certain dynamic range. Lasso-BLOT’s performance, however, is sensitive to the choice of regularization parameter

One of the most surprising attributes of BLOOMP is improved performance with respect to sparsity at larger dynamic range and low noise.

The algorithms BLOSP, BLOCoSaMP and BLOIHT are good alternatives to BLOOMP and BP/Lasso-BLOT: they are faster than both BLOOMP and BP/Lasso-BLOT and shares, to a lesser degree, BLOOMP’s desirable attribute with respect to dynamic range.

Comparisons with existing algorithms (SIHT and frame-adapted BP) demonstrate the superiority of BLO-enhanced algorithms for reconstruction of sparse objects separated above the Rayleigh length.

Finally to add to the debate of analysis versus synthesis , the performance of BLO-based algorithms for sparse, widely separated objects are independent of the refinement factors $F$ representing redundancy, and, since the discretization error decreases with $F$, the reconstruction errors of the BLO-based synthesis methods also decrease with $F$ in stark contrast to the examples presented in which show that the synthesis approach degrades with redundancy.

References