Infinite dimensional compressed sensing from anisotropic measurements and applications to inverse problems in PDE

Infinite dimensional compressed sensing from anisotropic measurements and applications to inverse problems in PDE

Abstract

We consider a compressed sensing problem in which both the measurement and the sparsifying systems are assumed to be frames (not necessarily tight) of the underlying Hilbert space of signals, which may be finite or infinite dimensional. The main result gives explicit bounds on the number of measurements in order to achieve stable recovery, which depends on the mutual coherence of the two systems. As a simple corollary, we prove the efficiency of nonuniform sampling strategies in cases when the two systems are not incoherent, but only asymptotically incoherent, as with the recovery of wavelet coefficients from Fourier samples. This general framework finds applications to inverse problems in partial differential equations, where the standard assumptions of compressed sensing are often not satisfied. Several examples are discussed, with a special focus on electrical impedance tomography.

alberti@dima.unige.it

matteo.santacesaria@helsinki.fi

1Introduction

The recovery of a sparse signal from a small number of samples is the fundamental question of compressed sensing (CS). A signal is called sparse if it can be expressed as a linear combination of a small number of known vectors. The seminal papers [23] have triggered an impressive amount of research in the last decade, from real world applications (MRI, X-ray tomography, etc.) to theoretical generalizations in broader mathematical frameworks [33].

In the finite dimensional case, the general CS problem can be stated as follows. Given an unknown sparse vector and a measurement operator represented by a matrix , we want to reconstruct from samples of the form , for . This is done by solving the convex optimization problem

where is the projection matrix on the entries indexed by . It is natural to ask under what conditions the solution of the minimization problem coincides with . These can be formulated as a lower bound on the number of measurements , which depends on the sparsity of the signal , the dimension of the ambient space, and the matrix . An interesting feature is that the lower bound on does not guarantee exact recovery for all set of indices with , but only for most of them.

One of the first contributions [23] considered the case where is the discrete Fourier transform: exact recovery is guaranteed with high probability provided that is selected uniformly at random with . If is a general unitary transformation, the problem has been addressed for the first time in [19], introducing the coherence . In this case, the bound becomes .

Similar results have been recently obtained in the infinite dimensional setting, where one considers signals belonging to a separable Hilbert space and the measurement operator is represented as a bounded linear map . (Note that may be expressed by scalar products with a family of vectors , namely .) The sparsity of a signal is characterized by the sparsity of , where , , is the analysis operator associated to a family of vectors . The first results in this framework were presented in [4], in the case where both and are unitary operators, i.e. correspond to orthonormal bases; the orthonormality of is a standard assumption taken in virtually all works on CS with deterministic measurements. These results were further extended in [5], introducing the more advanced concepts of asymptotic incoherence, local coherence, and local sparsity. An additional improvement was given in [63], which deals with the case where is a Parseval frame (see also [57]). Allowing the system to be a frame, i.e. is not necessarily invertible, is useful whenever we wish to use a redundant representation to sparsify the signals in (e.g. redundant wavelets [34], curvelets [21], ridgelets [20] and shearlets [52]).

In a large number of inverse problems, where one does not have complete freedom in the measurement process, the assumption on being unitary is not verified, thereby preventing the application of CS to many domains. As a result, two large research areas as inverse problems in partial differential equations and sparse signal recovery have been almost completely separated so far. The purpose of this paper is to provide a solid foundation that is expected to allow a fruitful interaction between these two domains.

In order to do so, in this work we present a very general CS result that deals with any bounded and injective linear operators and , defined on any separable complex Hilbert space (finite or infinite dimensional). Equivalently, the families and are simply required to be frames of (not necessarily tight). Since we do not need the measurement operator to be unitary, our results cover the case of anisotropic measurements. These have already been studied in the finite dimensional case using random and not deterministic measurements in [49]. As far as we know, our result is new also in the finite dimensional case.

Another generalization is related to the sampling strategy. Recently, it has been observed in several works [46] that, when precise bounds for the mutual coherence are available, uniform sampling strategies do not give sharp estimates for the minimum number of measurements. Our techniques are also able to cover this case, also known as structured sampling, just as a simple corollary of the main result for the uniform sampling. To our knowledge, this is the first infinite dimensional result under asymptotic incoherence assumptions, where there is no need to use multi-level sampling strategies and local coherence. For instance, we analyze the recovery of wavelet coefficients from Fourier samples, and justify the use of a nonuniform sampling scheme corresponding to the so called log sampling. This represents only a first step, and we believe that many other interesting estimates may be derived as corollaries of the main general result.

As mentioned above, our main motivation in dealing with the infinite dimensional anisotropic framework comes from inverse problems arising from partial differential equations [43]. These inverse problems are intrinsically infinite dimensional, and often the measurement operator cannot be chosen as a unitary transformation. Moreover, in order to obtain a solution to these problems, an infinite number of measurements is often needed, even when the signal to be recovered belongs to a known finite dimensional subspace. CS can thus provide a rigorous, explicit and numerically viable way to find solutions to such problems when only a finite number of measurements is available. We explore applications of our main result to the problems of nonuniform Fourier sampling, (linearized) electrical impedance tomography (EIT) and photoacoustic tomography. In particular, we show for the first time that an electrical conductivity may be recovered from a number of linearized EIT measurements proportional to the sparsity of the signal with respect to a wavelet basis, up to log factors. Many other inverse problems can be tackled with a similar approach and will be the subject of future work.

The plan of the paper is the following. In Section 2 we introduce the mathematical framework of infinite dimensional CS using the language of frames. We define the mutual coherence for general frames as well as the balancing property. The main result is presented in Section 3, which contains also the main corollary about structured sampling and asymptotic incoherence. Section 4 is devoted to applications to three inverse problems, while Section 5 contains the proof of the main result.

2Main assumptions

Let denote the set of all positive natural numbers. Let be a separable complex Hilbert space, representing our signal space, which may be either finite or infinite dimensional. The problem we study in this paper is the recovery of an unknown signal from partial measurements of the form , under a sparsity assumption on with respect to a suitable family of vectors . The main assumption of this paper is the following: these families of vectors are required to be frames of [24].

The measurements and the sparsity are expressed by the analysis operators and , defined by

By construction, the dual operators are given by and , where is the canonical basis of . By Hypothesis ?, since and , we have that and are bounded and the operator norms satisfy

The recovery problem can then be stated as follows: given noisy partial measurements of , namely for some (finite) set , recover the signal , under the assumption that is sparse. Here we have used the notation for the orthogonal projection onto (if we simply write ). The classical way to solve this problem is via minimization, namely

where is the noise level.

Equivalently, one may adapt a more abstract point of view, starting from a bounded operator with bounded inverse. It is immediate to verify that gives rise to a frame, as in Hypothesis ?. The formulation with allows to consider any linear inverse problem of the form

The only requirement is that, with full data, the inverse problem should be uniquely and stably solvable. In particular, any linear invertible operator may be considered, and not necessarily isometries as in the standard compressed sensing setting (see Section 4).

Given the generality of our setting, we need to consider the dual frames of and . By classical frame theory (see [27]), the frame operators and are invertible, and we can consider the dual frames

which have frame constants and , respectively. Equivalently, we may write and , where and are the Moore–Penrose pseudoinverses of and , respectively, defined as follows:

Note that they are left inverses of and , respectively. Therefore, and are the analysis operators of the dual frames, and so arguing as in we obtain

With an abuse of notation, we have denoted and by and , respectively. It can be immediately checked that they are right inverses of and , i.e. and . For later use, set and , so that by and we obtain

The frames and are the canonical dual frames, but in general many other choices are possible. These are in correspondence with all possible bounded left inverses of and , and it is possible to characterize all dual frames [27].

We need to measure the coherence between the sensing system and the representation system or, equivalently, between the measurement operator and the representation operator .

Let us now discuss a particular case.

The partial measurements are indexed by , where will be chosen uniformly at random in . In the infinite dimensional case, the upper bound has to be chosen big enough, depending on the sparsity of . This is quantified by the balancing property, introduced in [4] and extended here to general frames.

For , , we use the notation

which measures the compressibility of the signal by means of -sparse signals . As in [46], we introduce the localization factor

where

and . Note that when is an orthonormal basis or, equivalently, when is a unitary operator.

Following [63], for we denote

where we have used the notation

for an operator .

In the other extreme case, it may happen that , even in finite dimension with a Parseval frame, as the following example shows.

For , let be the smallest integer such that and

where .

3Main results

3.1Finite and infinite dimensional recovery

We now state the main result of this work. Recall that is any separable Hilbert space: we deal with the finite and infinite dimensional case simultaneously.

Theorem ? directly generalizes Theorems 6.1, 6.3 and 6.4 of [4] to the case of anisotropic measurements. It also extends the results of [5] to the case of general frames and . Our result can also be seen as an infinite dimensional generalization of the finite dimentional result in [49] for anisotropic random measurements.

3.2Asymptotic incoherence and virtual frames

The above result shows that with random sampling one needs a number of measurements proportional to the sparsity of the signal (up to logarithmic factors), provided that the coherence is small enough, namely, . While this happens in finite dimension () in some particular situations, e.g. with signals that are sparse with respect to the Dirac basis and with Fourier measurements, in many cases of practical interest the above result becomes almost meaningless since the coherence is of order one. For instance, this happens when is the discrete Fourier transform and the discrete wavelet transform. As it was shown in [5], this is always the case in infinite dimension.

Since the early stages of compressed sensing, it was realized that this issue may be solved by using variable density random sampling [68]. For instance, in the Fourier-Wavelet case, one needs to sample lower frequencies with higher probability than the higher frequencies. We now give a result that deals with this situation; in particular, it takes into account a priori estimates on the coherence and nonuniform sampling. As it is clear from the proof, it follows as a simple corollary of Theorem ?, thanks to the flexibility of its assumptions. More precisely, Theorem ? is applied to a virtual frame obtained from by artificially repeating its elements. More complicated transformations, also involving , may be considered (taking into account, for instance, asymptotic sparsity [5]): we leave these investigations to future work, and we limit ourselves to an example to show the potential of this framework.

For , , set

Let for . We want to apply Theorem ? to and , where the virtual frame , with associated operator , is given as follows. For normalize by and repeat it times, namely

The new index set coincides with if and is finite if is finite (more precisely, we have ). Note that has the same frame bounds of , i.e. , by construction, since

We now want to prove that satisfies the balancing property with respect to , , and . We first notice that , since

In passing, we remark that this identity tells us that the dual frame of the virtual frame coincides with the virtual frame of the dual frame , which is constructed as in . Arguing in the same way, and terminating the above sums to and , respectively, we readily derive

This immediately yields properties and , since and . Similarly, follows from the identities

We have the following straightforward upper bound for :

The factor associated to and , which we denote by , verifies . Indeed, from the definition of , we only need to check that , which follows by :

The factors and do not change since they do not depend on but only on , which is left unchanged.

Let us calculate the new coherence

For there exists with and , so that by and we obtain

since . Therefore .

Now, the factor in the estimate for given in Theorem ? applied to and becomes , so that the estimate on in Theorem ? transforms into . Finally, note that selecting elements uniformly at random from corresponds to the variable density sampling scheme of the statement (with ). Further, setting , we have

This concludes the proof.

3.3Recovery of wavelet coefficients from Fourier samples

For , let be the signal space. Let be the Fourier basis of , namely

Consider a nondecreasing ordering of , namely a bijective map , , such that

where is any norm of . Set for . Let be a separable wavelet basis of (ordered according to the wavelet scales). Note that both systems and are orthonormal bases, so that and . Under certain decay conditions on the scaling function (which may be relaxed to a condition satisfied by all Daubechies wavelets if one considers a different ordering of the frequencies ), it was shown in [44] that

for some . In other words, the wavelet basis and the Fourier basis are asymptotically incoherent. Thanks to the first of these inequalities, we have that assumption of the corollary is satisfied with , so that . Further, by the second of these inequalities and Remark ? we have As a consequence, estimate for the number of measurements becomes

for some constant depending only on . The number of measurements required for the success of the recovery using minimization is directly proportional to the sparsity of the signals, up to log factors, and so, up to log factors, estimate is optimal. It is worth observing that one log factor may be removed by using multilevel sampling, asymptotic sparsity and finer properties of wavelets [5].

According to our result, the measurements must be chosen at random from with probabilities

This nonuniform sampling scheme corresponds to a uniform sampling scheme for the virtual frame , in which each is repeated times (and suitably normalized, see the proof of Corollary ?). It is then natural to wonder how the frequencies associated to are arranged. Consider for simplicity the one-dimensional case with only the positive frequencies , so that the ordering is simply the identity map. By construction, the frequency satisfies , and so , which gives

This is the so called log sampling scheme [2] (up to a suitable scaling of the parameters), which is a 1D model for higher dimensional spiral trajectories; these are common in Magnetic Resonance Imaging. A comparison of the two sampling schemes is shown in Figure 1; the plot contains the negative frequencies as well. As expected, the log sampling scheme yields a smooth approximation of the sampling scheme . As far as the authors are aware, this is the first time that theoretical support is given to the use of the log sampling scheme in CS. It is worth observing that the normalization of the elements by the square root of the number of repetitions corresponds to the standard normalization of weighted Fourier frames by the measures of Voronoi regions [3].

Figure 1: A comparison between the log sampling scheme and the scheme associated to the virtual frame \{\hat\psi_{\hat l}\}_{\hat l}.
Figure 1: A comparison between the log sampling scheme and the scheme associated to the virtual frame .

4Applications

4.1Nonuniform Fourier sampling

The most classical compressed sensing problem formulated in the continuous setting is the recovery of a function from Fourier samples

where is a finite set of frequencies where the measurements are taken. Here is the discrete Fourier transform given by scalar products with the sinusoids , which form an orthonormal basis of . If is sparse with respect to a suitable orthonormal basis with analysis operator , this reconstruction problem fits in the framework discussed in the previous section, and may be recovered by minimization, provided that enough random measurements are taken. The standard theory of compressed sensing may be applied in this case, since both and are unitary operators (see Example ?).

In several applications (such as Magnetic Resonance Imaging, Computed Tomography, geophysical imaging, seismology and electron microscopy), nonuniform Fourier sampling arises naturally, i.e. the frequencies are not taken uniformly in . In this case, the operator fails to be an isometry, since the corresponding family of sinusoids may be only a frame, and not an orthonormal basis. The results discussed in the previous section may be directly applied to this case too.

Let us now give a quick overview of nonuniform Fourier frames; we follow [3]. For additional details, the reader is referred to [26] for the one-dimensional case, and to [13] for the multi-dimensional case.

Let be the space of square-integrable functions with support contained in a compact, convex and symmetric set , i.e. . For , we consider measurements of the form

namely scalar products with the sinusoids

Instead of considering the case when is a cartesian grid of , which gives rise to uniform Fourier sampling, we wish to give more general conditions on the set so that is a frame of .

The first of these conditions requires that the samples are fine enough to capture all the frequency information in a given direction.

The second condition limits the concentration of samples, in order to avoid large energies in small frequency regions.

Under these conditions, the family of sinusoids with frequencies in forms a frame for .

Now, assuming that is a frame for , we can apply Theorem ? and Corollary ? to this setting. This would provide, to our knowledge, the first result about recovery of a sparse signal from nonuniform Fourier measurements via minimization. Even if the measurement frame is generally not tight, we can provide explicit bounds for the recovery of the wavelet coefficients from nonuniform Fourier samples.

Some numerical simulations related to this framework are presented in [36].

4.2Electrical impedance tomography

EIT is an imaging technique in which one wants to determine the electrical conductivity inside a body from boundary voltage and current measurements. It is a non-linear inverse boundary value problem whose mathematical formulation was presented for the first time by Calderón [18].

Let , , be an open bounded domain with Lipschitz boundary and , for almost every , be the electrical conductivity. Given a voltage on the boundary of the domain, the associated potential is the unique solution of the following Dirichlet problem for the conductivity equation:

where , , are the classical Sobolev spaces. The boundary current associated to the voltage is represented by the trace of the normal derivative of the potential on . More precisely, we define the Dirichlet-to-Neumann (DN) map as

where is the unit outer normal to and is the unique solution of .

Calderón’s inverse conductivity problem asks if it is possible to determine a conductivity from the knowledge of its associated DN map . Positive answers to this question have been given since 1987 [66].

If is sufficiently smooth, the problem can be reduced to the so-called Gel’fand-Calderón problem for the Schrödinger equation,

via the change of variables in . This inverse problem consists in the reconstruction of the potential from the knowledge of the DN map

One of the biggest open questions concerning inverse boundary value problems such as Calderón’s or Gel’fand-Calderón’s is the determination of a conductivity/potential from a finite number of boundary measurements. A priori assumptions on the unknown are needed in this case, and to the best of our knowledge the only result concerns piecewise constant coefficients with discontinuities on a single convex polygon [35]. For conductivities/potentials belonging to some finite dimensional subspaces, an infinite number of measurements have always been a fundamental requirement to guarantee uniqueness and reconstruction. For instance, several works have studied the general piecewise constant case with infinitely many measurements [8]. In what follows, we will consider finitely many measurements, and present a first result in this direction for the linearized Gel’fand-Calderón problem, using the theory developed in this paper.

In order to linearize the problem, we assume that where is known and is small. Given two boundary voltages we have Alessandrini’s identity [7]:

where (resp. ) solves the Schrödinger equation with potential (resp. ) and Dirichlet data (resp. ). The quantity on the left of this identity is known since is known and is the boundary measurement corresponding to the chosen potential ( should be seen as a test function). Since for we have , the linearization consists in assuming that we can measure the quantity for given . Focusing on the solutions themselves instead of on their boundary values, this inverse problem may be rephrased as follows.

Without loss of generality, we can assume that , where . Extend by zero to and assume that for some (the letter will be used also for the sparsity, but there is never room for ambiguity). In the rest of this subsection, with an abuse of notation, several different positive constants depending only on , and will be denoted by the same letter . By a classical uniqueness result for the Calderón problem [66] we have that for every and we can construct solutions of in of the form

where are such that and

These solutions are known as exponentially growing solutions, Faddeev-type solutions [32] or complex geometrical optics (CGO) solutions.

We need to consider an ordering of , namely a bijective map , . For each fix and define the measurement operator by

We call the family a CGO frame (see Lemma ? below). Using the same ordering of , we define the discrete Fourier transform by

where .

We can now state the following consequence of Corollary ? for the linearized EIT problem.

Corollary ? provides a first general recipe to recover or approximate a sparse or compressible conductivity from a small number of linearized EIT measurements. The assumption that the sparsifying basis and the Fourier basis must be asymptotically incoherent is not restrictive: as already mentioned in Section 3.3, a large class of wavelet bases satisfy [44]. Note that in the 1D Fourier-Wavelet case we have (Remark ?).

Note that for the incoherence and balancing property we made assumptions only on the Fourier basis and not directly on the CGO frame . This is possible thanks to the following lemma, which also shows that is an invertible operator with bounded inverse, provided that the s are chosen big enough.

Since , setting we have

where we used estimate and the Sobolev embedding for . This implies , so that (note that the series is finite as ). Choosing immediately yields . The first part of follows. Hence , since is an isometry. Moreover, we have the Neumann series expansion

and so , as desired. In addition, the latter identity readily gives the second part of , since . The last estimate can be proven as follows. We readily derive

Since , the series on the right is convergent. Hence, by definition of the Sobolev norm we obtain

provided that . This finishes the proof of the lemma.

We are now in a position to prove Corollary ?.

We want to apply Corollary ? to and . Set , where is given by Lemma ?, so that .

First note that by and since is an orthonormal basis. We begin by showing that satisfies the balancing property (with the original right hand side in the definition) with respect to , , and . Consider condition ( and can be shown by using the same argument): by the triangle inequality we have