Stable and robust sampling strategies for compressive imaging

Introduction

The measurement process in a wide range of imaging applications such as radar, sonar, astronomy, and computer tomography, can be modeled – after appropriate approximation and discretization – as taking samples from weighted discrete Fourier transforms . Similarly, it is well known in the medical imaging literature that the measurements taken in Magnetic Resonance Imaging (MRI) are well modeled as Fourier coefficients of the desired image. Within all of these scenarios, one seeks strategies for taking frequency domain measurements so as to reduce the number of measurements without degrading the quality of image reconstruction. A central feature of natural images that can be exploited in this process is that they allow for approximately sparse representation in suitable bases or dictionaries.

The theory of compressive sensing, as introduced in , fits comfortably into this set-up: its key observation is that signals which allow for a sparse or approximately sparse representation can be recovered from relatively few linear measurements via convex approximation, provided these measurements are sufficiently incoherent with the basis in which the signal is sparse.

Much work in compressive sensing has focused on the setting of imaging with frequency domain measurements and, in particular, towards accelerating the MRI measurement process . For images that have sparse representation in the canonical basis, the incoherence between Fourier and canonical bases implies that uniformly subsampled discrete Fourier transform measurements can be used to achieve near-optimal oracle reconstruction bounds: up to logarithmic factors in the discretization size, any image can be approximated from $s$ such frequency measurements up to the error that would be incurred if the image were first acquired in full, and then compressed by setting all but the $s$ largest-magnitude pixels to zero . Compressive sensing recovery guarantees hold more generally subject to incoherence between sampling and sparsity transform domains. Unfortunately, natural images are generally not directly sparse in the standard basis, but rather with respect to transform domains more closely resembling wavelet bases. As low-scale wavelets are highly correlated (coherent) with low frequencies, sampling theorems for compressive imaging with partial Fourier transform measurements have remained elusive.

A number of empirical studies, including the very first papers on compressive sensing MRI , suggest that better image restoration is possible by subsampling frequency measurements from variable densities preferring low frequencies to high frequencies. In fact, variable-density MRI sampling had been proposed previously in a number of works outside the context of compressive sensing, although there did not seem to be a consensus on an optimal density . For a thorough empirical comparison of different densities from the compressive sensing perspective, see .

Guided by these observations, the authors of estimated the coherence between each element of the sensing basis with elements of the sparsity basis separately as a means to derive optimal sampling strategies in the context of compressive sensing MRI. In particular, they observe that incoherence-based results in compressive sensing imply exact recovery results for more general systems if one samples a row from the measurement basis proportionally to its squared maximal correlation with the sparsity-inducing basis. For given problem dimensions, they find the optimal distribution as the solution to a convex problem. In , the approach of optimizing the sampling distribution is combined with premodulation by a chirp, which by itself is another measure to reduce the coherence . A similar variable-density analysis already appeared in in the context of sampling strategies and reconstruction guarantees for functions with sparse orthogonal polynomial expansions and will also be the guiding strategy of this paper (cf. Section 5). After the submission of this paper, the idea of variable-density sampling has been extended to the context of block sampling , motivated by practical limitations of MRI hardware.

2 Contributions of this paper

In this paper we derive near-optimal reconstruction bounds for a particular type of variable-density subsampled discrete Fourier transform for both wavelet sparsity and gradient sparsity models. More precisely, up to logarithmic factors in the discretization size, any image can be approximated from $s$ such measurements up to the error that would be incurred if the wavelet transform or the gradient of the image, respectively, were first computed in full, and then compressed by setting all but the $s$ largest-magnitude coefficients to zero. Note that the reconstruction results that have been derived for uniformly-subsampled frequency measurements only provide such guarantees for images which are exactly sparse.

A major role in determining an appropriate sampling density will be played by the local coherence of the sensing basis with respect to the sparsity basis, as introduced in Section 5. Consequently, an important ingredient of our analysis is Theorem 6.2, which provides frequency-dependent bounds on inner products between rows of the orthonormal discrete Fourier transform and rows of the orthonormal discrete Haar wavelet transform. In particular, the maximal correlation between a fixed row in the discrete Fourier transform and any row of the discrete Haar wavelet transform decreases according to an inverse power law of the frequency, and decays sufficiently quickly that the sum of squared maximal correlations scales only logarithmically with the discretization size $N$. This implies, according to the techniques used in , that subsampling rows of the discrete Fourier matrix proportionally to the squared correlation results in a matrix that has the restricted isometry property of near-optimal order subject to appropriate rescaling of the rows.

For reconstruction, total variation minimization will be our algorithm of choice. In the papers and , total variation minimization was shown to provide stable and robust image reconstruction provided that the sensing matrix is incoherent with the Haar wavelet basis. Following the approach of , we prove that from variable density frequency samples, total variation minimization can be used for stable image recovery guarantees.

3 Outline

The remainder of this paper is organized as follows. Preliminary notation is introduced in Section 2. The main results of this paper are contained in Section 3. Section 4 reviews compressive sensing theory and Section 5 presents recent results on sampling strategies for coherent systems. The main results on the coherence between Fourier and Haar wavelet bases are provided in Section 6, and proofs of the main results are contained in Section 7. Section 8 illustrates our results by numerical examples. We conclude with a summary and a discussion of open problems in Section 9.

Preliminaries

Clearly, $\sigma_{s}(f)_{p}=0$ if $f$ is $s$-sparse. Informally, $f$ is called compressible if $\sigma_{s}(f)_{1}$ decays quickly as $s$ increases.

For two nonnegative functions $f(t)$ and $g(t)$ on the real line, we write $f\gtrsim g$ (or $f\lesssim g$) if there exists a constant $C>0$ such that $f(t)\geq Cg(t)$ (or $f(t)\leq Cg(t)$, respectively) for all $t>0$.

Here we note that our definition is the anisotropic version of the total variation semi-norm. The isotropic total variation semi-norm becomes the sum of terms

The isotropic and anisotropic total variation semi-norms are thus equivalent up to a factor of $\sqrt{2}$.

2 Bases for sparse representation and measurements

The Haar wavelet basis is a simple basis which allows for good sparse approximations of natural images. We will work primarily in two dimensions, but first introduce the univariate Haar wavelet basis as it will nevertheless serve as a building block for higher dimensional bases.

We will also work with discrete Fourier measurements.

indexed by discrete frequencies in the range $-N/2+1\leq k_{1},k_{2}\leq N/2.$

We denote by ${\cal F}$ the two-dimensional discrete Fourier transform $f\rightarrow\big(\left\langle f,\varphi_{k_{1},k_{2}}\right\rangle\big)_{k_{1},k_{2}}$ and, again, also the associated unitary matrix. Finally, we denote by ${\cal F}_{\Omega}$ its restriction to a set of frequencies $\Omega\subset[N]^{2}$.

Main results

Fix integers $N=2^{p},m,$ and $s$ such that $s\gtrsim\log(N)$ and

Select $m$ frequencies $\{(\omega_{1}^{j},\omega_{2}^{j})\}_{j=1}^{m}\subset\{-N/2+1,\dots,N/2\}^{2}$ i.i.d. according to

Given noisy partial Fourier measurements $y={\cal F}_{\Omega}f+\xi$, the estimation

approximates $f$ up to the noise level and best $s$-term approximation error of its gradient:

Fix integers $N=2^{p},m,$ and $s$ such that $s\gtrsim\log(N)$ and

approximates $f$ up to the noise level and best $s$-term approximation error in the bivariate Haar basis:

Even though the required number of samples $m$ in Theorem 3.2 is smaller than the number of samples required for the total variation minimization guarantees in Theorem 3.1, one finds that total variation minimization requires fewer measurements empirically. This may be due to the fact that the gradient of a natural image has stronger sparsity than its Haar wavelet representation. For this reason we focus on total variation minimization. Independent of this observation, we strongly suspect that the additional logarithmic factors in the number of measurements stated in Theorem 3.1 are an artifact of the proof, and that it should be possible to strengthen the result to obtain a similar recovery guarantee with the number of measurements as in Theorem 3.2. Moreover, one should be able to reduce the number of necessary log-factors with a RIP-less approach . These are important follow-up questions, as the current number of logarithmic factors may limit the direct applicability of our results to practical problems.

Compressive sensing background

In particular, reconstruction is exact, $x^{\#}=x$, if $x$ is $s$-sparse and $\varepsilon=0$.

There are stronger versions of this result which allow for weaker constraints on the restricted isometry constant . However, our version is a corollary of the following proposition, which appears as Proposition 2 in , and generalizes the results from . This proposition will also play an important role in the proof of our main results.

Suppose further that for a subset $S$ of cardinality $|S|=k$, the signal $u$ satisfies a cone constraint

Indeed, Proposition 4.2 follows from Proposition 4.3 by noting that the minimality of $x^{\#}$ implies a cone constraint for the residual $x-x^{\#}$ over the support of the $s$ largest-magnitude entries of $x$. The proof of Proposition 4.3 can be found in .

2 Bounded orthonormal systems

While the strongest known results on the restricted isometry property concern random matrices with independent entries such as Gaussian or Bernoulli, a scenario that has proven particularly useful for applications is that of structured random matrices with rows chosen from a basis incoherent to the basis inducing sparsity (see below for a detailed discussion on the concept of incoherence). The resulting sampling schemes correspond to bounded orthonormal systems, and such systems have been extensively studied in the compressive sensing literature (see for an expository article including many references).

Consider a set $T$ equipped with probability measure $\nu$.

An orthonormal system is said to be bounded with bound $K$ if $\sup_{j\in[N]}\|\psi_{j}(x)\|_{\infty}\leq K$.

For example, the basis of complex exponentials $\psi_{j}(x)=\exp{(i2\pi jx)}$ forms a bounded orthonormal system with optimally small constant $K=1$ with respect to the uniform measure on $T=\{0,\frac{1}{N},\dots,\frac{N-1}{N}\}$, and $d$-dimensional tensor products of complex exponentials form bounded orthonormal systems with respect to the uniform measure on the set $T^{d}$. A random sample of an orthonormal system is the vector $(\psi_{1}(x),\ldots,\psi_{N}(x))$, where $x$ is a random variable drawn according to the associated distribution $\nu$. Any matrix whose rows are independent random samples of a bounded orthonormal system, such as the uniformly subsampled discrete Fourier matrix, will have the restricted isometry property:

for some $s\gtrsim\log(N)$ For matrices consisting of uniformly subsampled rows of the discrete Fourier matrix, it has been shown in that this constraint is not necessary., then with probability at least $1-N^{-C\log^{3}(s)},$ the restricted isometry constant $\delta_{s}$ of $\frac{1}{\sqrt{m}}\Psi$ satisfies $\delta_{s}\leq\delta$.

Local coherence

The sparse recovery results in Corollary 4.6 based on mutual coherence do not take advantage of the full range of applicability of bounded orthonormal systems. As argued in , Proposition 4.5 implies comparable sparse recovery guarantees for a much wider class of sampling/sparsity bases through preconditioning resampled systems. In the following, we formalize this approach through the notion of local coherence.

and choose $m$ (possibly not distinct) indices $j\in\Omega\subset[N]$ i.i.d. from the probability measure $\nu$ on $[N]$ given by

Note that as the matrix $\Psi$ with rows $\psi_{k}$ is unitary, the vectors $\eta_{j}:=\Psi\phi_{j}$, $j=1,\dots,N$, form an orthonormal system with respect to the uniform measure on $[N]$ as well. We show that the system $\{\widetilde{\eta}_{j}\}=\{d_{j}\eta_{j}\}$ is an orthonormal system with respect to $\nu$ in the sense of Definition 4.4. Indeed,

hence the $\widetilde{\eta}_{j}$ form an orthonormal system with respect to $\nu$. Noting that $|\eta_{j}(k)|=|\langle\varphi_{j},\psi_{k}\rangle|\leq\kappa_{j}$ and hence this system is bounded with bound $\|\kappa\|_{2}$, the result follows from Proposition 4.5. ∎

Note that the local coherence not only appears in the embedding dimension $m$, but also in the sampling measure. Hence a priori, one cannot guarantee the optimal embedding dimension if one only has suboptimal bounds for the local coherence. That is why the sampling measure in Theorem 5.2 is defined via the (known) upper bounds $\kappa$ and $\|\kappa\|_{2}$ rather than the (usually unknown) exact values $\mu_{loc}$ and $\|\mu_{loc}\|_{2}$, showing that suboptimal bounds still lead to meaningful bounds on the embedding dimension.

For $\mu\leq KN^{-1/2}$ (as in Corollary 4.6), one has $\|\mu^{loc}\|_{2}\leq K$ , so Theorem 5.2 is a direct generalization of Corollary 4.6. As one has equality if and only if $\mu^{loc}$ is constant, however, Theorem 5.2 will be stronger in most cases.

Local coherence estimates for frequencies and wavelets

Due to the tensor product structure of both of these bases, the two-dimensional local coherence of the two-dimensional Fourier basis with respect to bivariate Haar wavelets will follow by first bounding the local coherence of the one-dimensional Fourier basis with respect to the set of univariate building block functions of the bivariate Haar basis.

To estimate this expression, we note that

If $0\neq|k|\leq 2^{p-2}$, we bound $|1-e^{2\pi i2^{-p}k}|\geq 2^{-p}|k|$ and apply (6.6) to obtain

For $2^{p-2}<|k|\leq 2^{p-1}$, and hence $2^{-p}\leq\frac{1}{2}|k|^{-1}$, we note that $|1-e^{2\pi i2^{-p}k}|\geq\frac{\sqrt{2}}{2}$ and bound, again using (6.6),

This lemma enables us to derive the following incoherence estimates for the bivariate case.

and one has $\|\kappa\|_{2}\leq\|\kappa^{\prime}\|_{2}\leq 52\sqrt{p}=52\sqrt{\log_{2}(N)}$.

First note that the bivariate Fourier coefficients decompose into the product of univariate Fourier coefficients:

For $k_{i}\neq 0$, the factors can be bounded using Lemma 6.1, which, for $k_{1}\neq 0\neq k_{2}$, yields the bound

In both cases, we obtain $\mu^{loc}_{k_{1},k_{2}}\leq\frac{18\pi}{\max(|k_{1}|,|k_{2}|)}$. The bound $\mu^{loc}_{k_{1},k_{2}}\leq 1$ follows directly from the Cauchy-Schwarz inequality. We conclude $\mu^{loc}_{k_{1},k_{2}}\leq\kappa(k_{1},k_{2})\leq\kappa^{\prime}(k_{1},k_{2})$.

where we used that $p\geq 8$. Taking square root implies the result. ∎

We believe that the factor of $\sqrt{\log_{2}N}$ which appears in the proposition is due to lack of smoothness for the Haar wavelets. Hence we hope this factor can be removed by considering smoother wavelets.

As the infimum of a strictly decreasing function and a strictly increasing function is bounded uniformly by its value at the intersection point of the two functions, Lemma 6.1 also gives frequency-dependent bounds for the local coherence between frequencies and wavelets in the univariate setting.

Recovery guarantees

In this section we present proofs of the main results.

2 Preliminary lemmas for the proof of Theorem 3.1

The proof of Theorem 3.1 proceeds along similar lines to that of Theorem 3.2, but we need a few more preliminary results relating the bivariate Haar transform to the gradient transform. The first result, Proposition 7.1, is derived from a more general statement involving the continuous bivariate Haar system and the bounded variation seminorm from .

See for a derivation of Proposition 7.1 from Theorem $8.1$ of .

We will also need two lemmas about the bivariate Haar system.

3 Proof of Theorem 3.1

By Theorem 5.2 combined with the bivariate incoherence estimates from Theorem 6.2, we know that with high probability ${\cal A}:=\frac{1}{\sqrt{m}}D{\cal F}_{\Omega}{\cal H}^{*}$ has the restricted isometry property of order $s$ and level $\delta$ once

Thus, for the stated number of measurements $m$ with an appropriate hidden constant, we can assume that ${\cal A}$ has the restricted isometry property of order

where the exact value of the constant $\widetilde{C}$ will be determined below. In the remainder of the proof we show that this event implies the result.

Let $u=f-f^{\#}$ denote the residual error of (3.3). Then we have

Cone Constraint on $\nabla u$. Let $S$ denote the support of the best $s$-sparse approximation to $\nabla f$. Since $f^{\#}=f-u$ is the minimizer of (TV) and $f$ is also a feasible solution,

Cone Constraint on $w^{u}={\cal H}u$. Proposition 7.1 allows us to pass from a cone constraint on the gradient to a cone constraint on the Haar transform. More specifically, we obtain

recalling that $w^{u}_{(1)}$ is the coefficient associated to the constant wavelet. Now consider the set $\widetilde{S}$ consisting of the $s$ edges indexed by $S$. By Lemma 7.2, the set $\Lambda$ indexing those wavelets which change sign across edges in $\widetilde{S}$ has cardinality at most $|\Lambda|=6s\log(N)$. Decompose $u$ as

and note that by linearity of the gradient,

Now, by construction of the set $\Lambda$, we have that $(\nabla u_{\Lambda^{c}})_{S}=0$ and so $(\nabla u)_{S}=(\nabla u_{\Lambda})_{S}$. By Lemma 7.3 and the triangle inequality,

Combined with Proposition 7.1 concerning the decay of the wavelet coefficients and the cone constraint (7.1), and letting

this gives rise to a cone constraint on the wavelet coefficients:

Since both $f$ and $f^{\#}$ are in the feasible region of (3.3), we have for $u=f-f^{\#}$,

Using the derived cone and tube constraints on ${\cal H}u$ along with the assumed RIP bound on ${\cal A}$, the proof is complete by applying Proposition 4.3 using $\gamma=\widetilde{C}\log(N^{2}/s)\leq 2\widetilde{C}\log(N)$, $k=6s\log N$, and $\xi=\widetilde{C}\log(N^{2}/s)\|\nabla f-(\nabla f)_{S}\|_{1}$. In fact, this is where we need that the RIP order is $\overline{s}$, to accommodate for the factors $\gamma$ and $k$. ∎

Numerical illustrations

In this section, we will provide numerical examples for our results. As there have been papers entirely devoted to the empirical investigation of optimal sampling strategies , the goal will be to illustrate our results rather than provide a thorough empirical analysis.

First, we consider a $256\times 256$ spine image and visually compare the reconstruction quality for different spine images. While the inferior reconstruction quality for uniform sampling is obvious, the difference between variable density sampling and using only the low frequencies is less apparent, both visually and in the reconstruction error.

One should remark that all of these experiments were performed without preconditioning in the regularization term, while our results contain such a step. Preliminary experiments suggest that this may be an artifact of the proof, and for this reason, our experiments were carried out with the standard noise model in the reconstruction procedure. A more in depth comparison of various noise models and weighting in the reconstruction poses an interesting object of study for future work. Note that weighted noise models similar to the one resulting from our analysis have been explored in .

Summary and outlook

We established reconstruction guarantees for variable-density discrete Fourier measurements in both the wavelet sparsity and gradient sparsity setup. Our results build on local coherence estimates between Fourier and wavelet bases. Although we derive local coherence estimates only for 1D and 2D Fourier/wavelet systems, such estimates can be extended to higher dimensions by induction, using the tensor-product structure of these bases.

Variable density sampling in compressive imaging has often been justified as taking into account the tree-like sparsity structure of natural images in wavelet bases (e.g., in ). We note that our theory does not directly take such signal statistics into account, and depends only on the local incoherence between Fourier and wavelet bases. Incorporating this additional structure to derive stronger reconstruction guarantees, by either improved sampling strategies or improved reconstruction strategies, remains an interesting and important direction of future research.

All the recovery guarantees in this paper are uniform, that is, we seek measurement ensembles which allow for approximate reconstruction of all images. For non-uniform recovery guarantees, we expect that the number of measurements required in our main results can be reduced by several logarithmic factors by following a probabilistic and “RIP-less” approach .

It should also be noted that this paper does not address the important issue of errors arising from discretization of the image and Fourier measurements. In particular, as observed for example in , the use of discrete rather than continuous Fourier representations can be a significant source of error in compressive sensing. The authors of propose to resolve this issue using uneven sections, that is, the number of discretization points in frequency is chosen to be larger than the number of discretization points in time. Nevertheless, the results in are again just formulated for incoherent samples. Recently, it has been proposed to overcome this issue by sampling all of the low frequencies in addition to uniformly sampling the higher frequencies . After the submission of this paper, reconstruction guarantees for such a setup were provided in , also for a generalization to multilevel sampling schemes. In addition to an asymptotic notion of coherence (related to the local coherence we look at in this paper), also considers an asymptotic notion of sparsity, which relates to the additional structure of wavelet expansions mentioned above.

We think that it should be an interesting to study how our approach can be applied to infinite dimensional image models – due to the variable density, it may even be possible to sample from all of the infinite set rather than restricting to a finite subset. Such a generalization would prove challenging for the optimization-based approaches such as in , which will always be specific to the given problem dimension. In this sense, we expect that the additional understanding provided by this paper can eventually lead to optimized sampling schemes. All these questions, however, are left for future work.

Acknowledgments

The authors would like to thank Ben Adcock, Anders Hansen, Deanna Needell, Holger Rauhut, Justin Romberg, Amit Singer, Mark Tygert, Robert Vanderbei, Yves Wiaux, and the anonymous reviewers for helpful comments and suggestions. They are grateful for the stimulating research environment of the Mathematisches Forschungsinstitut Oberwolfach, where part of this work was completed. Rachel Ward was supported in part by an Alfred P Sloan Research Fellowship, a Donald D. Harrington Faculty Fellowship, an NSF CAREER grant, and DOD-Navy grant N00014-12-1-0743.

References