Compressed Sensing with Coherent and Redundant Dictionaries

Introduction

Compressed sensing is a new data acquisition theory based on the discovery that one can exploit sparsity or compressibility when acquiring signals of general interest, and that one can design nonadaptive sampling techniques that condense the information in a compressible signal into a small amount of data . In a nutshell, reliable, nonadaptive data acquisition, with far fewer measurements than traditionally assumed, is possible. By now, applications of compressed sensing are abundant and range from imaging and error correction to radar and remote sensing, see and references therein.

$A$ is an $m\times n$ sensing matrix with $m$ typically smaller than $n$ by one or several orders of magnitude (indicating some significant undersampling) and $z$ is an error term modeling measurement errors. Sensing is nonadaptive in that $A$ does not depend on $x$. Then the theory asserts that if the unknown signal $x$ is reasonably sparse, or approximately sparse, it is possible to recover $x$, under suitable conditions on the matrix $A$, by convex programming: we simply find the solution to

where $\|x\|_{0}=|\{i:x_{i}\neq 0\}|$. In words, $x_{s}$ is the best $s$-sparse approximation to the vector $x$, where we shall say that a vector is $s$-sparse if it has at most $s$ nonzero entries. Put differently, $x-x_{s}$ is the tail of the signal, consisting of the smallest $n-s$ entries of $x$. In particular, if $x$ is $s$-sparse, $x-x_{s}=0$. Then with this in mind, one of the authors improved on the work of Candès, Romberg and Tao and established that $(L_{1})$ recovers a signal $\hat{x}$ obeying

provided that the $2s$-restricted isometry constant of $A$ obeys $\delta_{2s}<\sqrt{2}-1$. The constants in this result have been further improved, and it is now known to hold when $\delta_{2s}<0.4652$ , see also . In short, the recovery error from $(L_{1})$ is proportional to the measurement error and the tail of the signal. This means that for compressible signals, those whose coefficients obey a power law decay, the approximation error is very small, and for exactly sparse signals it completely vanishes.

The definition of restricted isometries first appeared in where it was shown to yield the error bound (1.3) in the noiseless setting, i. e. when $\varepsilon=0$ and $z=0$.

For an $m\times n$ measurement matrix $A$, the $s$-restricted isometry constant $\delta_{s}$ of $A$ is the smallest quantity such that

With this, the condition underlying (1.3) is fairly natural since it is interpreted as preventing sparse signals from being in the nullspace of the sensing matrix $A$. Further, a matrix having a small restricted isometry constant essentially means that every subset of $s$ or fewer columns is approximately an orthonormal system. It is now well known that many types of random measurement matrices have small restricted isometry constants . For example, matrices with Gaussian or Bernoulli entries have small restricted isometry constants with very high probability whenever the number of measurements $m$ is on the order of $s\log(n/s)$. The fast multiply matrix consisting of randomly chosen rows of the discrete Fourier matrix also has small restricted isometry constants with very high probability with $m$ on the order of $s(\log n)^{4}$.

Although there are countless applications for which the signal of interest is represented by some overcomplete dictionary, the compressed sensing literature is lacking on the subject. Consider the simple case in which the sensing matrix $A$ has Gaussian (standard normal) entries. Then the matrix $AD$ relating the observed data with the assumed (nearly) sparse coefficient sequence $x$ has independent rows but each row is sampled from $\mathcal{N}(0,\Sigma)$, where $\Sigma=D^{*}D$. If $D$ is an orthonormal basis, then these entries are just independent standard normal variables, but if $D$ is not unitary then the entries are correlated, and $AD$ may no longer satisfy the requirements imposed by traditional compressed sensing assumptions. In recovery results are obtained when the sensing matrix $A$ is of the form $\Phi D^{*}$ where $\Phi$ satisfies the restricted isometry property. In this case the sampling matrix must depend on the dictionary $D$ in which the signal is sparse. We look for a universal result which allows the sensing matrix to be independent from the signal and its representation. To be sure, we are not aware of any such results in the literature guaranteeing good recovery properties when the columns may be highly – and even perfectly – correlated.

Before continuing, it might be best to fix ideas to give some examples of applications in which redundant dictionaries are of crucial importance.

The Discrete Fourier Transform (DFT) matrix is an $n\times n$ orthogonal matrix whose $k$th column is given by

with the convention that $0\leq t,k\leq n-1$. Signals which are sparse with respect to the DFT are only those which are superpositions of sinusoids with frequencies appearing in the lattice of those in the DFT. In practice, we of course rarely encounter such signals. To account for this, one can consider the oversampled DFT in which the sampled frequencies are taken over even smaller equally spaced intervals, or at small intervals of varying lengths. This leads to an overcomplete frame whose columns may be highly correlated.

Recall that for a fixed function $g$ and positive time-frequency shift parameters $a$ and $b$, the $k$th column (where $k$ is the double index $k=(k_{1},k_{2})$) of the Gabor frame is given by

Radar and sonar along with other imaging systems appear in many engineering applications, and the goal is to recover pulse trains given by

Due to the time-frequency structure of these applications, Gabor frames are widely used . If one wishes to recover pulse trains from compressive samples by using a Gabor dictionary, standard results do not apply.

Curvelets provide a multiscale decomposition of images, and have geometric features that set them apart from wavelets and the likes. Conceptually, the curvelet transform is a multiscale pyramid with many directions and positions at each length scale, and needle-shaped elements at fine scales . The transform gets its name from the fact that it approximates well the curved singularities in an image. This transform has many properties of an orthonormal basis, but is overcomplete. Written in matrix form $D$, it is a tight frame obeying the Parseval relations

where we let $\{d_{k}\}$ denote the columns of $D$. Although columns of $D$ far apart from one another are very uncorrelated, columns close to one another have high correlation. Thus none of the results in compressed sensing apply for signals represented in the curvelet domain.

In many applications a signal may not be sparse in a single orthonormal basis, but instead is sparse over several orthonormal bases. For example, a linear combination of spikes and sines will be sparse when using a concatenation of the coordinate and Fourier bases. One also benefits by exploiting geometry and pointwise singularities in images by using combinations of tight frame coefficients such as curvelets, wavelets, and brushlets. However, due to the correlation between the columns of these concatenated bases, current compressed sensing technology does not apply.

These and other applications strongly motivate the need for results applicable when the dictionary is redundant and has correlations. This state of affair, however, exposes a large gap in the literature since current compressed sensing theory only applies when the dictionary is an orthonormal basis, or when the dictionary is extremely uncorrelated (see e.g. ).

2 Do we really need incoherence?

Current assumptions in the field of compressed sensing and sparse signal recovery impose that the measurement matrix have uncorrelated columns. To be formal, one defines the coherence of a matrix $M$ as

where $M_{j}$ and $M_{k}$ denote columns of $M$. We say that a dictionary is incoherent if $\mu$ is small. Standard results then require that the measurement matrix satisfy a strict incoherence property , as even the RIP imposes this. If the dictionary $D$ is highly coherent, then the matrix $AD$ will also be coherent in general.

Coherence is in some sense a natural property in the compressed sensing framework, for if two columns are closely correlated, it will be impossible in general to distinguish whether the energy in the signal comes from one or the other. Recall that when the dictionary $D$ is sufficiently incoherent, standard compressed sensing guarantees that we recover $x$ and thus $f=Dx$, provided $x$ is $s$-sparse with $s$ sufficiently small. For example, imagine that we are not undersampling and that $A$ is the identity so that we observe $y=Dx$. Suppose the first two columns are identical, $d_{1}=d_{2}$. Then the measurement $d_{1}$ can be explained by the input vectors $(1,0,\ldots,0)$ or $(0,1,0,\ldots,0)$ or any convex combination. Thus there is no hope of reconstructing a unique sparse signal $x$ from measurements $y=ADx$. However, we are not interested in recovering the coefficient vector $x$, but rather the actual signal $Dx$. The large correlation between columns in $D$ now does not impose a problem because although it makes it impossible to tell apart coefficient vectors, this is not the goal. This simple example suggests that perhaps coherence is not necessary. If $D$ is coherent, then we clearly cannot recover $x$ as in our example, but we may certainly be able to recover the signal $f=Dx$ from measurements $y=Af$ as we shall see next.

3 Gaussian sensing matrices

Our main result is that the solution to $(P_{1})$ is very accurate provided that $D^{*}f$ has rapidly decreasing coefficients. Our result for the Gaussian case is below while the general theorem appears in Section 1.5.

Let $D$ be an arbitrary $n\times d$ tight frame and let $A$ be a $m\times n$ Gaussian matrix with $m$ on the order of $s\log(d/s)$. Then the solution $\hat{f}$ to $(P_{1})$ obeys

for some numerical constants $C_{0}$ and $C_{1}$, and where $(D^{*}f)_{s}$ is the vector consisting of the largest $s$ entries of $D^{*}f$ in magnitude as in (1.2).

4 Implications

In words, if the columns of the Gram matrix are reasonably sparse and if $f$ happens to have a sparse expansion, then the frame coefficient sequence $D^{*}f$ is also sparse. All the transforms discussed above, namely, the Gabor, curvelet, wavelet frame, oversampled Fourier transform all have nearly diagonal Gram matrices – and thus, sparse columns.

We now turn to the implications of our result to the applications we have already mentioned, and instantiate the theorem in the noiseless case due to the optimality of the noise level in the error.

To recover multitone signals, we use an oversampled DFT, which is not orthonormal and may have very large coherence. However, since each “off-grid” tone has a rapidly decaying expansion, $D^{*}f$ will have rapidly decaying coefficients. In practice, one smoothly localizes the data to a time interval by means of a nice window $w$ to eliminate effects having to do with a lack of periodicity. One can then think of the trigonometric exponentials as smoothly vanishing at both ends of the time interval under study. Thus our result implies that the recovery error is negligible when the number of measurements is about the number of tones times a log factor.

An easy example. As above, let $D$ be the $n\times 2n$ dictionary consisting of a concatenation of the identity and the DFT, normalized to ensure $D$ is a tight frame (below, $F$ is the DFT normalized to be an isometry):

We wish to create a sparse signal that uses linearly dependent columns for which there is no local isometry. Assume that $n$ is a perfect square and consider the Dirac comb

5 Axiomization

We now turn to the generalization of the above result, and give broader conditions about the sensing matrix under which the recovery algorithm performs well. We will impose a natural property on the measurement matrix, analogous to the restricted isometry property.

Let $\Sigma_{s}$ be the union of all subspaces spanned by all subsets of $s$ columns of $D$. We say that the measurement matrix $A$ obeys the restricted isometry property adapted to $D$ (abbreviated D-RIP) with constant $\delta_{s}$ if

($\gamma$ is an arbitrary positive numerical constant) will satisfy the D-RIP with overwhelming probability, provided that $m\gtrsim s\log(d/s)$. This can be seen by a standard covering argument (see e.g. the proof of Lemma 2.1 in ). Many types of random matrices satisfy (1.5). It is now well known that matrices with Gaussian, subgaussian, or Bernoulli entries satisfy (1.5) with number of measurements $m$ on the order of $s\log(d/s)$ (see e.g. ). It has also been shown that if the rows of $A$ are independent (scaled) copies of an isotropic $\psi_{2}$ vector, then $A$ also satisfies (1.5). Recall that an isotropic $\psi_{2}$ vector $a$ is one that satisfies for all $v$,

for some constant $\alpha$. See for further details. Finally, it is clear that if $A$ is any of the above random matrices then for any fixed unitary matrix $U$, the matrix $AU$ will also satisfy the condition.

The D-RIP can also be analyzed via the Johnson-Lindenstrauss lemma (see e.g. ). There are many results that show certain types of matrices satisfy this lemma, and these would then satisfy the D-RIP via (1.5). Subsequent to our submission of this manuscript, Ward and Krahmer showed that randomizing the column signs of any matrix that satisfies the standard RIP yields a matrix which satisfies the Johnson-Lindenstrauss lemma . Therefore, nearly all random matrix constructions which satisfy standard RIP compressed sensing requirements will also satisfy the D-RIP. A particularly important consequence is that because the randomly subsampled Fourier matrix is known to satisfy the RIP, this matrix along with a random sign matrix will thus satisfy D-RIP. This gives a fast transform which satisfies the D-RIP. See Section 4.2 for more discussion.

We are now prepared to state our main result.

Let $D$ be an arbitrary tight frame and let $A$ be a measurement matrix satisfying D-RIP with $\delta_{2s}<0.08$. Then the solution $\hat{f}$ to $(P_{1})$ satisfies

where the constants $C_{0}$ and $C_{1}$ may only depend on $\delta_{2s}$.

Remarks. We actually prove that the theorem holds under the weaker condition $\delta_{7s}\leq 0.6$, however we have not tried to optimize the dependence on the values of the restricted isometry constants; refinements analagous to those in the compressed sensing literature are likely to improve the condition. Further, we note that since Gaussian matrices with $m$ on the order of $s\log(d/s)$ obey the D-RIP, Theorem 1.2 is a special case of Theorem 1.4.

6 Organization

The rest of the paper is organized as follows. In Section 2 we prove our main result, Theorem 1.4. Section 3 contains numerical studies highlighting the impact of our main result on some of the applications previously mentioned. In Section 4 we discuss further the implications of our result along with its advantages and challenges. We compare it to other methods proposed in the literature and suggest an additional method to overcome some impediments.

Proof of Main Result

We now begin the proof of Theorem 1.4, which is inspired by that in . The new challenge here is that although we can still take advantage of sparsity, the vector possessing the sparse property is not being multiplied by something that satisfies the RIP, as in the standard compressed sensing case. Rather than bounding the tail of $f-\hat{f}$ by its largest coefficients as in , we bound a portion of $D^{*}h$ in an analagous way. We then utilize the D-RIP and the fact that $D$ is a tight frame to bound the error, $\|f-\hat{f}\|_{2}$.

Let $f$ and $\hat{f}$ be as in the theorem, and let $T_{0}$ denote the set of the largest $s$ coefficients of $D^{*}f$ in magnitude. We will denote by $D_{T}$ the matrix $D$ restricted to the columns indexed by $T$, and write $D_{T}^{*}$ to mean $(D_{T})^{*}$. With $h=f-\hat{f}$, our goal is to bound the norm of $h$. We will do this in a sequence of short lemmas. The first is a simple consequence of the fact that $\hat{f}$ is the minimizer.

The vector $D^{*}h$ obeys the following cone constraint,

Proof. Since both $f$ and $\hat{f}$ are feasible but $\hat{f}$ is the minimizer, we must have $\|D^{*}\hat{f}\|_{1}\leq\|D^{*}f\|_{1}$. We then have that

This implies the desired cone constraint.

We next divide the coordinates $T_{0}^{c}$ into sets of size $M$ (to be chosen later) in order of decreasing magnitude of $D_{T_{0}^{c}}^{*}h$. Call these sets $T_{1},T_{2},\ldots$, and for simplicity of notation set $T_{01}=T_{0}\cup T_{1}$. We then bound the tail of $D^{*}h$.

Setting $\rho=s/M$ and $\eta=2\|D_{T_{0}^{c}}^{*}f\|_{1}/\sqrt{s}$, we have the following bound,

Proof. By construction of the sets $T_{j}$, we have that each coefficient of $D_{T_{j+1}}^{*}h$, written $|D_{T_{j+1}}^{*}h|_{(k)}$, is at most the average of those on $T_{j}$:

This along with the cone constraint in Lemma 2.1 gives

With $\rho=s/M$ and $\eta=2\|D_{T_{0}^{c}}^{*}f\|_{1}/\sqrt{s}$, it follows from Lemma 2.1 and the Cauchy-Schwarz inequality that

Next we observe that by the feasibility of $\hat{f}$, $Ah$ must be small.

Proof. Since $\hat{f}$ is feasible, we have

We will now need the following result which utilizes the fact that $D$ satisfies the D-RIP.

Proof. Since $D$ is a tight frame, $DD^{*}$ is the identity, and this along with the D-RIP and Lemma 2.2 then imply the following:

Since we also have $\|D_{T_{0}}^{*}h\|_{2}\leq\|h\|_{2}$, this yields the desired result.

We now translate these bounds to the bound of the actual error, $\|h\|_{2}$.

The error vector $h$ has norm that satisfies,

Proof. Since $D^{*}$ is an isometry, we have

where the last inequality follows from Lemma 2.2.

We next observe an elementary fact that will be useful. The proof is omitted.

For any values $u$, $v$ and $c>0$, we have

We may now conclude the proof of Theorem 1.4. First we employ Lemma 2.6 twice to the inequality given by Lemma 2.5 (with constants $c_{1}$, $c_{2}$ to be chosen later) and the bound $\|D_{T_{0}}^{*}h\|_{2}\leq\|h\|_{2}$ to get

Using the fact that $\sqrt{u^{2}+v^{2}}\leq u+v$ for $u,v\geq 0$, we can further simply to get our desired lower bound,

It only remains to choose the parameters $c_{1}$, $c_{2}$, and $M$ so that $K_{1}$ is positive. We choose $c_{1}$=1, $M=6s$, and take $c_{2}$ arbitrarily small so that $K_{1}$ is positive when $\delta_{7s}\leq 0.6$. Tighter restrictions on $\delta_{7s}$ will of course force the constants in the error bound to be smaller. For example, if we set $c_{1}=1/2$, $c_{2}=1/10$, and choose $M=6s$, we have that whenever $\delta_{7s}\leq 1/2$ that $(P_{1})$ reconstructs $\hat{f}$ satisfying

Note that if $\delta_{7s}$ is even a little smaller, say $\delta_{7s}\leq 1/4$, the constants in the theorem are just $C_{1}=10.3$ and $C_{2}=7.33$. Note further that by Corollary 3.4 of , $\delta_{7s}\leq 0.6$ is satisfied whenever $\delta_{2s}\leq 0.08$. This completes the proof.

Numerical Results

In these experiments, we test the performance on a simulated real-world signal from the field of radar detection. The test input is a superposition of six radar pulses. Each pulse has a duration of about 200 ns, and each pulse envelope is trapezoidal, with a 20 ns rise and fall time, see Figure 2. For each pulse, the carrier frequency is chosen uniformly at random from the range 50 MHz to 2.5 GHz. The Nyquist interval for such signals is thus 0.2 ns. Lastly, the arrival times are distributed at random in a time interval ranging from $t=0\text{ s}$ to $t\approx 1.64\,\mu\text{s}$; that is, the time interval under study contains $n=8192$ Nyquist intervals. We acquire this signal by taking $400$ measurements only, so that the sensing matrix $A$ is a Gaussian matrix with $400$ rows. The dictionary $D$ is a Gabor dictionary with Gaussian windows, oversampled by a factor of about $60$ so that $d\approx 60\times 8,192=491,520$. The main comment about this setup is that the signal of interest is not exactly sparse in $D$ since each pulse envelope is not Gaussian (the columns of $D$ are pulses with Gaussian shapes) and since both the frequencies and arrival times are sampled from a continuous grid (and thus do not match those in the dictionary).

Because $D$ is massively overcomplete, the Gram matrix $D^{*}D$ is not diagonal. Figure 5 depicts part of the Gram matrix $D^{*}D$ for this dictionary, and shows that this matrix is “thick” off of the diagonal. We can observe visually that the dictionary $D$ is not an orthogonal system or even a matrix with low coherence, and that columns of this dictionary are indeed highly correlated. Having said this, the second plot in Figure 5 shows the rapid decay of the sequence $D^{*}f$ where $f$ is the signal in Figure 2.

Discussion

2 Fast Transforms

These results yield a transform with a fast multiply which satisfies the D-RIP. The number of measurements and the multiply and storage costs of the matrix are of the same magnitude as those that satisfy the RIP. The D-RIP is, therefore, satisfied by matrices with the same benefits as those in standard compressed sensing. This shows that compressed sensing with redundant and coherent dictionaries is viable with completely the same advantages as in the standard setting.

Acknowledgements

This work is partially supported by the ONR grants N00014-10-1-0599 and N00014-08-1-0749, the Waterman Award from NSF, and the NSF DMS EMSW21-VIGRE grant. EJC would like to thank Stephen Becker for valuable help with the simulations.

References