Total Variation Minimization in Compressed Sensing

TL;DR

This paper reviews recovery guarantees for total variation minimization in compressed sensing, highlighting the limitations of synthesis sparse approaches and extending results from Gaussian to subgaussian measurements.

Contribution

It provides a comprehensive overview of total variation minimization in compressed sensing and introduces generalized guarantees for subgaussian measurement scenarios.

Findings

Total variation minimization has specific recovery guarantees in compressed sensing.

Synthesis sparse signal approaches are inadequate for total variation minimization.

Recent results are extended from Gaussian to subgaussian measurement models.

Abstract

This chapter gives an overview over recovery guarantees for total variation minimization in compressed sensing for different measurement scenarios. In addition to summarizing the results in the area, we illustrate why an approach that is common for synthesis sparse signals fails and different techniques are necessary. Lastly, we discuss a generalizations of recent results for Gaussian measurements to the subgaussian case.

Figures

Equations

$$\Vert x{\Vert }_0:=|\{k\in [n]:x_k\ne 0\}|\le s,$$

$${\sigma }_s(x{)}_q:=\underset{\Vert z{\Vert }_0\le s}{\inf }\Vert z-x{\Vert }_q$$

$$\text{minimize }\Vert z{\Vert }_1\text{ such that }\Vert Az-y{\Vert }_2\le \epsilon .$$

$$\Vert x-\hat{x}{\Vert }_2\le c\left(\frac{{\sigma }_s(x{)}_1}{\sqrt{s}}+\epsilon \right),$$

$$\text{minimize }\Vert z{\Vert }_{TV}=\Vert \nabla z{\Vert }_1\text{ such that }\Vert Az-y{\Vert }_2\le \epsilon .$$

$$\text{minimize }\Vert \nabla z{\Vert }_1\text{such that }Az=y.$$

$$\text{mbox}ker({\mathcal{I}}_{I^c}^T\nabla )\cap \text{mbox}ker(A)=\{0\}.$$

$$\forall w\in \ker (A)\backslash \{0\}:\Vert (\nabla w{)}_I{\Vert }_1<\Vert (\nabla w{)}_{I^c}{\Vert }_1.$$

$${\nabla }^{T}v=A^{T}w,v_I=\text{mbox}sign(\nabla x{)}_I,\Vert v_{I^c}{\Vert }_{\infty }<1.$$

$$\Vert ({\mathcal{I}}_{I^c}^T\nabla ({\nabla }^T{\mathcal{I}}_{I^c}{\mathcal{I}}_{I^c}^T\nabla +A^{T}A{)}^{-1}\nabla \text{mbox}sign(\nabla x){)}_I{\Vert }_{\infty }<1$$

$$\mathbb{P}(\forall x_0:\Vert \nabla x_0{\Vert }_0\le s:\hat{x}=x_0)\ge 1-c_3{\mathrm{e}}^{-c_4\sqrt{m}}.$$

$$m\ge \{\begin{array}{l}{c}_{2}s{\log }^{3}n,\text{ if }d=2\\ c_{d}s\log n,\text{ if }d\ge 3,\end{array}$$

$$\text{minimize }\Vert \nabla z{\Vert }_1\text{such that }\Vert Az-y{\Vert }_2\le \epsilon .$$

$$\mathbb{P}\left(\Vert x_0-\hat{x}{\Vert }_2\le c_2\frac{{\min }_{|S|\le s}\Vert (\nabla x_0{)}_{S^c}{\Vert }_1}{\sqrt{n}}+c_3\frac{\epsilon }{\sqrt{n}}\right)\ge 1-c_0{\mathrm{e}}^{-c_{1}n}.$$

$$\mathbb{P}\left(\Vert \hat{x}-x_0\Vert >\frac{2\epsilon }{c\sqrt[4]{ns}(\sqrt{\log (2n)}-1)}\right)\le {\mathrm{e}}^{-\tilde{c}\sqrt{ns}}.$$

$$(1-\delta )\Vert x{\Vert }_2^2\le \Vert \Phi x{\Vert }_2^2\le (1-\delta )\Vert x{\Vert }_2^2.$$

$$\hat{x}={\text{argmin}}_z\Vert z{\Vert }_{T{V}_2}\text{ such that }\Vert \mathcal{M}(z)-y{\Vert }_2\le \epsilon$$

$$m\gtrsim k(\log (n)+\log ({\eta }^{-1})).$$

$$m\gtrsim s{\log }^3(s){\log }^5(N).$$

$$\mathbb{P}\big{[}(\omega_{1}^{j},\omega_{2}^{j})=(k_{1},k_{2})\big{]}=C_{N}\min\left(C,\frac{1}{k_{1}^{2}+k_{2}^{2}}\right)=:\eta(k_{1},k_{2}),\quad-N/2+1\leq k_{1},k_{2}\leq N/2,$$

$$f^{\#}={\text{argmin}}_{g\in {\mathbb{C}}^{N\times N}}\Vert g{\Vert }_{TV}\text{such that}\Vert \rho \circ ({\mathcal{F}}_{\Omega }g-y){\Vert }_2\le \epsilon \sqrt{m},$$

$$\Vert f-f^{\#}{\Vert }_2\lesssim \frac{\Vert \nabla f-(\nabla f{)}_s{\Vert }_1}{\sqrt{s}}+\epsilon .$$

$$\Vert z{\Vert }_2\lesssim \gamma +\Vert z{\Vert }_{TV}\text{mbox}with\Vert Az{\Vert }_2\le \gamma .$$

$$\Vert \nabla z-\nabla x{\Vert }_2\lesssim \left(\delta \sqrt{k}+C_1\frac{\Vert P\nabla z{\Vert }_1}{\sqrt{k}}\right),\frac{\Vert z-x{\Vert }_2}{\sqrt{n}}\lesssim C_2\left(\frac{\delta }{\sqrt{s}}+C_1\frac{\Vert P\nabla z{\Vert }_1}{k}\right),$$

$$C_1=\log (k){\log }^{1/2}(m)\text{mbox},and{C}_2={\log }^2(k)\log (n)\log (m).$$

$$\Vert \nabla z-\nabla x{\Vert }_2\lesssim \left(\delta \sqrt{k}+C_3\frac{\Vert P\nabla z{\Vert }_1}{\sqrt{k}}\right),\Vert z-x{\Vert }_2\lesssim C_2\left(\delta +C_3\frac{\Vert P\nabla z{\Vert }_1}{k}\right),$$

$$C_3=\log (k)\log (n^2/k){\log }^{1/2}(n){\log }^{1/2}(m).$$

$$\Vert \nabla z-\nabla x{\Vert }_2\lesssim \left(\delta \sqrt{k}+C\frac{\Vert P\nabla z{\Vert }_1}{\sqrt{k}}\right),\frac{\Vert z-x{\Vert }_2}{\sqrt{n}}\lesssim (\delta \sqrt{s}+C\Vert P\nabla z{\Vert }_1),$$

$$w(K):=\mathbb{E}\underset{x\in K-K}{\sup }⟨g,x⟩,$$

$$w(\mathrm{ch}(K))=w(K).$$

Full text

1 Introduction

The central aim of Compressed Sensing (CS) [CRT06, Don06] is the recovery of an unknown vector from very few linear measurements. Put formally, we would like to recover $x\in\mathbb{R}^{n}$ from $y=Ax+e\in\mathbb{R}^{m}$ with $m\ll n$, where $e$ denotes additive noise.

For general $x$, recovery is certainly not possible, hence additional structural assumptions are necessary in order to be able to guarantee recovery. A common assumption used in CS is that the signal is sparse. Here for $x$ we assume

[TABLE]

that is, there are only very few nonzero entries of $x$. And say that $x$ is $s$-sparse for some given sparsity level $s\ll n$. We call a vector compressible, if it can be approximated well by a sparse vector. To quantify the quality of approximation, we let

[TABLE]

denote the error of the best $s$-sparse approximation of $x$.

In most cases, the vector $x$ is not sparse in the standard basis, but there is a basis $\Psi$, such that $x=\Psi z$ and $z$ is sparse. This is also known as synthesis sparsity of $x$. To find an (approximately) synthesis sparse vector, we can instead solve the problem of recovering $z$ from $y=A\Psi z$. A common strategy in CS is to solve a basis pursuit program in order to recover the original vector. For a fixed noise level $\varepsilon$, it is given by

[TABLE]

While this and related approaches of convex regularization have been studied in the inverse problems and statistics literature long before the field of compressed sensing developed, these works typically assumed the measurement setup was given. The new paradigm arising in the context of compressed sensing was to attempt to use the remaining degrees of freedom of the measurement system to reduce the ill-posedness of the system as much as possible. In many measurement systems, the most powerful known strategies will be based on randomization, i.e., the free parameters are chosen at random.

Given an appropriate amount of randomness (i.e., for various classes of random matrices $A$, including some with structure imposed by underlying applications), one can show that the minimizer $\hat{x}$ of (1) recovers the original vector $x$ with error

[TABLE]

see, e.g., [BDDW08] for an elementary proof in the case of subgaussian matrices without structure, and [KR14] for an overview, including many references, of corresponding results for random measurement systems with additional structure imposed by applications. Note that (2) entails that if $x$ is $s$-sparse and the measurements are noiseless, the recovery is exact.

For many applications, however, the signal model of sparsity in an orthonormal basis has proven somewhat restrictive. Two main lines of generalization have been proposed. The first line of work, initiated by [RSV08] is the study of sparsity in redundant representation systems, at first under incoherence assumptions on the dictionary. More recently, also systems without such assumptions have been analyzed [CENR10, KNW15]. The main idea of these works is that even when one cannot recover the coefficients correctly due to conditioning problems, one may still hope for a good approximation of the signal.

The second line of work focuses on signals that are sparse after the application of some transform, one speaks of cosparsity or analysis sparsity [NDEG13], see, e.g., [KR15] for an analysis of the Gaussian measurement setup in this framework. A special case of particular importance, especially for imaging applications, is that of sparse gradients. Namely, as it turns out, natural images often admit very sparse approximations in the gradient domain, see, e.g., Figure 1. Here the discrete gradient at location $i=(i_{1},\dots,i_{n})$ is defined as the vector with its $n$ entries given by \big{(}(\nabla z)_{i}\big{)}_{j}=z_{i+e_{j}}-z_{i}, $j=1,\dots,n$, where $e_{j}$ is the $j$-th standard basis vector.

A first attempt to recover a gradient sparse signal is to formulate a compressed sensing problem in terms of the sparse gradient. When this is possible (for instance in the example of Fourier measurements [CRT06]), applying (1) will correspond to minimizing $\|\nabla z\|_{1}=:\|z\|_{TV}$, the total variation seminorm. Then (under some additional assumptions) compressed sensing recovery guarantees of the form (2) can apply. This proof strategy, however, only allows for showing that the gradient can be approximately recovered, not the signal. When no noise is present and the gradient is exactly sparse (which is not very realistic), this allows for signal recovery via integrating the gradient, but in case of noisy measurements, this procedure is highly unstable.

Nevertheless, the success motivates to minimize the total variation seminorm if one attempts to recover the signal directly, not the gradient. In analogy with (1), this yields the following minimization problem.

[TABLE]

For $A$ the identity (i.e., not reducing the dimension), this relates to the famous Rudin-Osher-Fatemi functional, a classical approach for signal and image denoising [rudin1992nonlinear]. Due to its high relevance for image processing, this special case of analysis sparsity has received a lot of attention recently also in the compressed sensing framework where $A$ is dimension reducing. The purpose of this chapter is to give an overview of recovery results for total variation minimization in this context of compressed sensing (Section 2) and to provide some geometric intuition by discussing the one-dimensional case under Gaussian or subgaussian measurements (to our knowledge, a generalization to the latter case does not appear yet in the literature) with a focus on the interaction between the high-dimensional geometry and spectral properties of the gradient operator (Section 3).

2 An overview over TV recovery results

In this section, we will give an overview of the state of the art guarantees for the recovery of gradient sparse signals via total variation minimization. We start by discussing in Section 2.1 sufficient conditions for the success of TV minimization.

Subsequently, we focus on recovery results for random measurements. Interestingly, the results in one dimension differ severely from the ones in higher dimensions. Instead of obtaining a required number of measurements roughly on the order of the sparsity level $s$, we need $\sqrt{sn}$ measurements for recovery. We will see this already in Subsection 2.2, where we present the results of Cai and Xu [CX15] for recovery from Gaussian measurements. In Section 3, we will use their results to obtain refined results for noisy measurements as well as guarantees for subgaussian measurements, combined with an argument of Tropp [Tro15]. In Subsection 2.3 we will present results by Ward and Needell for dimensions larger or equal than two showing that recovery can be achieved from Haar incoherent measurements.

2.1 Sufficient Recovery Conditions

Given linear measurements $Ax=y$ for an arbitrary $A\in\mathbb{R}^{m\times n}$ and a signal $x$ with $\|\nabla x\|_{0}\leq s$, a natural way to recover $x$ is by solving

[TABLE]

For $I\subset[n]$ we denote $A_{I}$ as the columns of $A$ indexed by $I$, and for a consecutive notation we denote $\mathcal{I}_{I}^{T}\nabla$ as the rows of $\nabla$ indexed by $I$ and $\mathcal{I}$ as the identity matrix. The following results can also be easily applied to analysis $\ell_{1}$-minimization, where any arbitrary matrix $D\in\mathbb{R}^{p\times n}$ replaces $\nabla$ in (3), as well as to any real Hilbert space setting [Kru15].

In many applications it is important to verify whether there is exactly one solution of (3). Since $\nabla$ is not injective here, we cannot easily use the well-known recovery results in compressed sensing [FR13] for the matrix $A\nabla^{\dagger}$. However, a necessary conditon can be given since $x$ can only satisfy $Ax=y$ and $(\nabla x)_{I^{c}}=0$ if

[TABLE]

If $\nabla$ is replaced by the identity, this is equivalent to $A_{I}$ being injective. Since this injectivity condition is unavoidable, we assume for the rest of this section that it is satisfied.

The paper [NDEG13] provides sufficient and necessary conditons for uniform recovery via (3). The conditions rely on the null space of the measurements and are hard to verify similar to the classical compressed sensing setup [TP14]. The following result is a corollary of these conditions. It no longer provides a necessary condition, but is more manageable.

Corollary** 2.1****.**

[NDEG13] For all $x\in\mathbb{R}^{n}$ with $s:=\|\nabla x\|_{0}$, the solution of (3) with $y=Ax$ is unique and equal to $x$ if for all $I\subset[n]$ with $|I|\leq s$ it holds that

[TABLE]

To consider measurements for specific applications, where it is difficult to prove whether uniform recovery is guaranteed, one can empirically examine whether specific elements $x$ solve (3) uniquely. For computed tomography measurements, a Monte Carlo Experiment is considered in [JKL15] to approximate the fraction of all gradient $s$-sparse vectors to uniquely solve (3). The results prompt that there is a sharp transition between the case that every vector with a certain gradient sparsity is uniquely recoverable and the case that TV-minimization will find a different solution than the desired vector. This behavior empirically agrees with the phase transition in the classical compressed sensing setup with Gaussian measurements [Don04].

To efficiently check whether many specific vectors $x$ can be uniquely recovered via (3), one needs to establish characteristics of $x$ which must be easily verifiable. Such a non-uniform recovery condition is given in the following theorem.

Theorem** 2.1****.**

[JKL15] It holds that $x\in\mathbb{R}^{n}$ is a unique solution of (3) if and only if there exists $w\in\mathbb{R}^{m}$ and $v\in\mathbb{R}^{n-1}$ such that

[TABLE]

The basic idea of the proof is to use the optimality condition for convex optimization problems [Roc72]. Equivalent formulations of the latter theorem can be found in [ZMY16, KR15] where the problem is considered from a geometric perspective. However, verifying the conditions in Theorem 2.1 still requires solving a linear program where an optimal $v$ for (4) needs to be found. In classical compressed sensing, the Fuchs Condition [Fuc04] is known as a weaker result as it suggests a particular $w$ in (4) and avoids solving the consequential linear program. The following result generalizes this result to general analysis $\ell_{1}$-minimization.

Corollary** 2.2****.**

If $x\in\mathbb{R}^{n}$ satisfies

[TABLE]

then $x$ is the unique solution of (3).

2.2 Recovery from Gaussian measurements

As discussed above, to date no deterministic constructions of compressed sensing matrices are known that get anywhere near an optimal number of measurements. Also for the variation of aiming to recover approximately gradient sparse measurements, the only near-optimal recovery guarantees have been established for random measurement models. Both under (approximate) sparsity and gradient sparsity assumptions, an important benchmark is that of a measurement matrix with independent standard Gaussian entries. Even though such measurements are hard to realize in practice, they can be interpreted as the scenario with maximal randomness, which often has particularly good recovery properties. For this reason, the recovery properties of total variation minimization have been analyzed in detail for such measurements. Interestingly, as shown by the following theorem, recovery properties in the one-dimensional case are significantly worse than for synthesis sparse signals and also for higher dimensional cases. That is why we focus on this case in Section 3, providing a geometric viewpoint and generalizing the results to subgaussian measurements.

Theorem** 2.2****.**

[CX15] Let the entries of $A\in\mathbb{R}^{m\times n}$ be i.i.d. standard Gaussian random variables and let $\hat{x}$ be a solution of (3) with input data $y=Ax_{0}$. Then

There exist constants $c_{1},c_{2},c_{3},c_{4}>0$, such that for $m\geq c_{1}\sqrt{sn}(\log n+c_{2})$

[TABLE] 2. 2.

For any $\eta\in(0,1)$, there are constants $\tilde{c}_{1},\tilde{c}_{2}>0$ and a universal constant $c_{2}>0$, such that for $s\geq\tilde{c}_{0}$ and $(s+1)<\frac{n}{4}$. If $m\leq\tilde{c}_{1}\sqrt{sn}-\tilde{c}_{2}$, there exist infinitely many $x_{0}\in\mathbb{R}^{n}$ with $\|\nabla x_{0}\|_{0}\leq s$, such that $\mathbb{P}(\hat{x}\neq x_{0})\geq 1-\eta$.

This scaling is notably different from what is typically obtained for synthesis sparsity, where the number of measurements scales linearly with $s$ up to $\log$ factors. Such a scaling is only obtained for higher dimensional signals, e.g., images. Indeed, in [CX15], it is shown that for dimensions at least two the number of Gaussian measurements sufficient for recovery is

[TABLE]

where the constant $c_{d}$ depends on the dimension.

Furthermore, as we can see in Theorem 2.5 below, this is also the scaling one obtains for dimensions larger than $1$ and Haar incoherent measurements. Thus the scaling of $\sqrt{sn}$ is a unique feature of the $1$-dimensional case. Also note that the square-root factor in the upper bound makes the result meaningless for a sparsity level on the order of the dimension. This has been addressed in [KRZ15], showing that a dimension reduction is also possible if the sparsity level is a (small) constant multiple of the dimension.

The proof of Theorem 2.2 uses Gordon’s escape through the mesh Theorem [Gor88]. We will elaborate on this topic in Section 3.

In case we are given noisy measurements $y=Ax_{0}+e$ with $\|e\|_{2}\leq\varepsilon$, we can instead of solving (3) consider

[TABLE]

If $\nabla x_{0}$ is not exactly, but approximately sparse, and our measurements are corrupted with noise, the following result can be established.

Theorem** 2.3****.**

[CX15] Let the entries of $A\in\mathbb{R}^{m\times n}$ be i.i.d. standard Gaussian random variables and let $\hat{x}$ be a solution of (5) with input data $y$ satisfying $\|Ax_{0}-y\|_{2}\leq\varepsilon$. Then for any $\alpha\in(0,1)$, there are positive constants $\delta,c_{0},c_{1},c_{2},c_{3}$, such that for $m=\alpha n$ and $s=\delta n$

[TABLE]

This looks remarkably similar to the recovery guarantees obtained for compressed sensing, note however that the number of measurements needs to be proportional to $n$, which is not desirable. We will present a similar result with improved number of measurements in Section 3.5.

Theorem** 2.4****.**

(Corollary of Theorem 3.4) Let $x_{0}\in\mathbb{R}^{n}$ be such that $\|\nabla x_{0}\|\leq s$ for $s>0$ and $A\in\mathbb{R}^{m\times n}$ with $m\geq C\sqrt{ns}\log(2n)$ be a standard Gaussian matrix. Furthermore, set $y=Ax_{0}+e$, where $\|e\|\leq\varepsilon$ denotes the (bounded) error of the measurement and for some absolute constants $c,\tilde{c}>0$ the solution $\hat{x}$ of (12) satisfies

[TABLE]

Note, however that in contrast to theorem 2.3, this theorem does not cover the case of gradient compressible vectors, but on the other hand Theorem 3.4 also incorporates the case of special subgaussian measurement ensembles. Also, if we set $s=\delta n$, we reach a similar conclusion as in Theorem 2.3.

2.3 Recovery from Haar-incoherent measurements

For dimensions $d\geq 2$, Needell and Ward [NW13a, NW13b] derived recovery results for measurement matrices having the restricted isometry property (RIP) when composed with the Haar wavelet transform. Here we say that a matrix $\Phi$ has the RIP of order $k$ and level $\delta$ if for every $k$-sparse vector $x$ it holds that

[TABLE]

The results of [NW13a, NW13b] build upon a connection between a signal’s wavelet representation and its total variation seminorm first noted by Cohen, Dahmen, Daubechies and DeVore [CDDD03].

Their theorems yield stable recovery via TV minimization for $N^{d}$ dimensional signals. For $d=2$, notably these recovery results concern images of size $N\times N$.

Several definitions are necessary in order to be able to state the theorem. The $d$ dimensional discrete gradient is defined via $\nabla\colon\mathbb{R}^{C^{d}}\to\mathbb{C}^{N^{d}\times d}$ and maps $x\in\mathbb{C}^{N^{d}}$ to its discrete derivative which, for each $\alpha\in[N]^{d}$ is a vector $(\nabla x)_{\alpha}\in\mathbb{C}^{d}$ composed of the derivatives in all $d$ directions. Up to now, we have always used the anisotropic version of the TV seminorm, which can be seen as taking the $\ell_{1}$ norm of the discrete gradient. The isotropic TV seminorm is defined via a combination of $\ell_{2}$ and $\ell_{1}$ norms. It is given by $\|z\|_{TV_{2}}:=\sum_{\alpha\in[N]^{d}}\|(\nabla z)_{\alpha}\|_{2}$. The result in [NW13a] is given in terms of the isotropic TV seminorm but can also be formulated for the anisotropic version.

Furthermore, we will need to concatenate several measurement matrices in order to be able to state the theorem. This will be done via the concatenation operator $\oplus\colon\mathrm{Lin}(\mathbb{C}^{n},\mathbb{C}^{k_{1}})\times\mathrm{Lin}(\mathbb{C}^{n},\mathbb{C}^{k_{2}})\to\mathrm{Lin}(\mathbb{C}^{n},\mathbb{C}^{k_{1}+k_{2}})$, which ’stacks’ two linear maps.

Finally, we need the notion of shifted operators. For an operator $\mathcal{B}\colon\mathbb{C}^{N^{l-1}\times(N-1)\times N^{d-l}}\to\mathbb{C}^{q}$, these are defined as the operators $\mathcal{B}_{0_{l}}\colon\mathbb{C}^{Nˆd}\to\mathbb{C}^{q}$ and $\mathcal{B}^{0_{l}}\colon\mathbb{C}^{Nˆd}\to\mathbb{C}^{q}$ concatenating a column of zeros to the end or beginning of the $l$-th component, respectively.

Theorem** 2.5**** ([NW13a]).**

Let $N=2^{n}$ and fix integers $p$ and $q$. Let $\mathcal{A}\colon\mathbb{C}^{N^{d}}\to\mathbb{C}^{p}$ be a map that has the restricted isometry property of order $2ds$ and level $\delta<1$ if it is composed with the orthonormal Haar wavelet transform. Furthermore let $\mathcal{B}_{1},\ldots,\mathcal{B}_{d}$ with $\mathcal{B}_{j}\colon\mathbb{C}^{(N-1)N^{d-1}}\to\mathbb{C}^{q}$ be such that $\mathcal{B}=\mathcal{B}_{1}\oplus\mathcal{B}_{2}\oplus\cdots\oplus\mathcal{B}_{d}$ has the restricted isometry property of order $5ds$ and level $\delta<\frac{1}{3}$. Consider the linear operator $\mathcal{M}=\mathcal{A}\oplus[\mathcal{B}_{1}]_{0_{1}}\oplus[\mathcal{B}_{1}]^{0_{1}}\oplus\cdots\oplus[\mathcal{B}_{d}]_{0_{d}}\oplus[\mathcal{B}_{d}]^{0_{d}}$. Then $\mathcal{M}\colon\mathbb{C}^{N^{d}}\to\mathbb{C}^{m}$ with $m=2dq+p$ and for all $x\in\mathbb{C}^{N^{d}}$ we have the following. Suppose we have noisy measurements $y=\mathcal{M}(x)+e$ with $\|e\|_{2}\leq\varepsilon$, then the solution to

[TABLE]

satisfies

$\|\nabla(x-\hat{x})\|_{2}\leq c_{1}\left(\frac{\|\nabla x-(\nabla x)_{S}\|_{1,2}}{\sqrt{s}}+\sqrt{d}\varepsilon\right)$, 2. 2.

$\|x-\hat{x}\|_{TV_{2}}\leq c_{2}\left(\|\nabla x-(\nabla x)_{S}\|_{1,2}+\sqrt{sd}\varepsilon\right)$, 3. 3.

$\|x-\hat{x}\|_{2}\leq c_{3}d\log N\left(\frac{\|\nabla x-(\nabla x)_{S}\|_{1,2}}{\sqrt{s}}+\sqrt{d}\varepsilon\right),$

for some absolute constants $c_{1},c_{2},c_{3}$.

From the last point of the previous theorem, we see that for noiseless measurements and gradient sparse vectors $x$, perfect recovery can be achieved provided the RIP assumption holds. Subgaussian measurement matrices, for example, will have the RIP, also when composed with the Haar wavelet transform $H$ (this is a direct consequence of rotation invariance). Moreover, as shown in [KW11], randomizing the column signs of an RIP matrix will, with high probability, also yield a matrix that has the RIP when composed with $H$. An important example is a subsampled Fourier matrix with random column signs, which relates to spread spectrum MRI (cf. [PMG*+*12]).

2.4 Recovery from subsampled Fourier measurements

Fourier measurements are widely used in many applications. Especially in medical applications as parallel-beam tomography and magnetic resonance imaging it is desirable to reduce the number of samples to spare patients burden. In Section 2.1, this is a motivation for introducing algorithmic checks for unique solutions of (3). In this section, we consider a probabilistic approach where an incomplete measurement matrix $A\in\mathbb{C}^{m\times n}$ chosen from the discrete Fourier transform on $\mathbb{C}^{N}$ is considered. Therefore we consider a subset $\Omega$ of the index set $\{-\lfloor n/2\rfloor+1,...,\lceil n/2\rceil\}$, where $\Omega$ consists of $m$ integers chosen uniformly at random and, additionally, $0\in\Omega$. Hence, we want to recover a signal, sparse in the gradient domain, with a measurement matrix $A=(e^{2\pi ikj/n})_{k\in\Omega,j\in[n]}$. In [CRT06] the optimal sampling cardinality for $s$-sparse signals in the gradient domain was given and enables to recover one-dimensional signals signals from $\mathcal{O}(k\log(n))$ Fourier samples. It naturally extends to two dimensions.

Theorem** 2.6****.**

[CRT06] With probability exceeding $1-\eta$, a signal $z$, which is $k$-sparse in the gradient domain is the unique solution of (3) if

[TABLE]

As already discussed in the introduction, the proof of this result proceeds via recovering the gradient and then using that the discrete gradient (with periodic boundary conditions) is injective. Due to the poor conditioning of the gradient, however, this injectivity results do not directly generalize to recovery guarantees for noisy measurements. For two (and more) dimensions, such results can be obtained via the techniques discussed in the previous subsection.

These techniques, however, do not apply directly. Namely, the Fourier (measurement) basis is not incoherent to the Haar wavelet basis; in fact, the constant vector is contained in both, which makes them maximally coherent. As observed in [PVW11], this incoherence phenomenon only occurs for low frequencies, the high frequency Fourier basis vectors exhibit small inner products to the Haar wavelet basis. This can be taken into account using a variable density sampling scheme with sampling density that is larger for low frequencies and smaller for high frequencies. For such a sampling density, one can establish the restricted isometry for the corresponding randomly subsampled discrete Fourier matrix combined with the Haar wavelet transform with appropriately rescaled rows [KW14]. This yields the following recovery guarantee.

Theorem** 2.7****.**

[KW14]

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

[TABLE]

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

[TABLE]

where $C$ is an absolute constant and $C_{N}$ is chosen such that $\eta$ is a probability distribution.

Consider the weight vector $\rho=(\rho_{j})_{j=1}^{m}$ with $\rho_{j}=(1/\eta(\omega^{j}_{1},\omega^{j}_{2}))^{1/2}$, and assume that the noise vector $\xi=(\xi_{j})_{j=1}^{m}$ satisfies $\|\rho\hskip 1.42262pt\circ\hskip 1.42262pt\xi\|_{2}\leq\varepsilon\sqrt{m}$, for some $\epsilon>0$. Then with probability exceeding $1-N^{-C\log^{3}(s)}$, the following holds for all images $f\in\mathbb{C}^{N\times N}$:

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

[TABLE]

where $\circ$ denotes the Hadamard product, approximates $f$ up to the noise level and best $s$-term approximation error of its gradient:

[TABLE]

A similar optimality result is given in [Poo15], also for noisy data and inexact sparsity. In contrast to the previous result, this result includes the one-dimensional case. The key to obtaining such a result is showing that the stable gradient recover implies the stable signal recovery, i.e.,

[TABLE]

Again the sampling distribution is chosen as a combination of the uniform distribution and a decaying distribution. The main idea is to use this sampling to establish (10) via the RIP. We skip technicalities for achieving the optimality in the following theorem and refer to the original article for more details.

Theorem** 2.8****.**

[Poo15] Let $z\in\mathbb{C}^{n}$ be fixed and $x$ be a minimizer of (5) with $\varepsilon=\sqrt{m}\delta$ for some $\delta>0$, $m\gtrsim k\log(n)(1+\log(\eta^{-1}))$, and an appropriate sampling distribution. Then with probability exceeding $1-\eta$, it holds that

[TABLE]

where $P$ is the orthogonal projection onto a $k$-dimensional subspace,

[TABLE]

In the two-dimensional setting the result changes to

[TABLE]

with remaining $C_{2}$ and

[TABLE]

These results are optimal since the best error one can archive [NW13b] is $\|z-x\|_{2}\lesssim\|P\nabla z\|_{1}k^{-1/2}$.

The optimality in the latter theorems is achieved by considering a combination of uniform random samling and variable density sampling. Uniform sampling on its own can achieve robust and stable recovery. However, the following theorem shows that the signal error is no longer optimal but the bound on the gradient error is still optimal up to log factors. Here (10) is obtained by using the Poincaré inequality.

Theorem** 2.9****.**

[Poo15] Let $z\in\mathbb{C}^{n}$ be fix and $x$ be a minimizer of (5) with $\varepsilon=\sqrt{m}\delta$ for some $\delta>0$ and $m\gtrsim k\log(n)(1+\log(\eta^{-1}))$ with random uniform sampling. Then with probability exceeding $1-\eta$, it holds that

[TABLE]

where $P$ is the orthogonal projection onto a $k$-dimensional subspace and $C=\log(k)\log^{1/2}(m)$.

3 TV-recovery from subgaussian measurements in 1D

In this section, we will apply the geometric viewpoint discussed in [Ver15] to the problem, which will eventually allow us to show the TV recovery results for noisy subgaussian measurements mentioned in Section 2.2.

As in the original proof of the 1D recovery guarantees for Gaussian measurements [CX15], the Gaussian mean width will play an important role in our considerations.

Definition** 3.1****.**

The (Gaussian) mean width of a bounded subset $K$ of $\mathbb{R}^{n}$ is defined as

[TABLE]

where $g\in\mathbb{R}^{n}$ is a vector of i.i.d. $\mathcal{N}(0,1)$ random variables.

In [CX15], the mean width appears in the context of the Gordon’s escape through the mesh approach [Gor88] (see Section 3.4 below), but as we will see, it will also be a crucial ingredient in applying the Mendelson small ball method [KM15, Men14].

The mean width has some nice (and important) properties, it is for example invariant under taking the convex hull, i.e.,

[TABLE]

Furthermore, it is also invariant under translations of $K$, as $(K-x_{0})-(K-x_{0})=K-K$. Due to the rotational invariance of Gaussian random variables, that is $Ug\sim g$, we also have that $w(UK)=w(K)$. Also, it satisfies the inequalities

[TABLE]

which are equalities if $K$ is symmetric about [math], because then $K=-K$ and hence $K-K=2K$.

3.1 $M^{*}$ bounds and recovery

In order to highlight the importance of the Gaussian mean width in signal recovery, we present some arguments from [Ver15]. Thus in this section we present a classical result, the $M^{*}$ bound, which connects the mean width to recovery problems, cf. [Ver15]. Namely, recall that due to rotational invariance, the kernel of a Gaussian random matrix $A\in\mathbb{R}^{m\times n}$ is a random subspace distributed according to the uniform distribution (the Haar measure) on the Grassmannian

[TABLE]

Consequently, the set of all vectors that yield the same measurements directly correspond to such a random subspace.

The average size of the intersection of this subspace with a set reflecting the minimization objective now gives us an average bound on the worst case error.

Theorem** 3.1**** ($M^{*}$ bound, Theorem 3.12 in [Ver15]).**

Let $K$ be a bounded subset of $\mathbb{R}^{n}$ and $E$ be a random subspace of $\mathbb{R}^{n}$ of drawn from the Grassmanian $G_{n,n-m}$ according to the Haar measure. Then

[TABLE]

where $C$ is absolute constant.

Given the $M^{*}$-bound it is now straightforward to derive bounds on reconstructions from linear observations. We first look at feasibility programs - which in turn can be used to obtain recovery results for optimization problems. For that, let $K\subset\mathbb{R}^{n}$ be bounded and $x\in K$ be the vector we seek to reconstruct from measurements $Ax=y$ with a Gaussian matrix $A\in\mathbb{R}^{m\times n}$.

Corollary** 3.1****.**

[MPTJ07] Choose $\hat{x}\in\mathbb{R}^{n}$, such that

[TABLE]

then one has, for an absolute constant $C^{\prime}$,

[TABLE]

This corollary directly follows by choosing $C^{\prime}=2C$, observing that $\hat{x}-x\in K-K$, and that the side constraint enforces $A(\hat{x}-x)=0$.

Via a standard construction in functional analysis, the so called Minkowski functional, one can now cast an optimization problem as a feasiblity program so that Corollary 3.1 applies.

Definition** 3.2****.**

The Minkowski functional of a bounded, symmetric set $K\subset\mathbb{R}^{n}$ is given by

[TABLE]

So the Minkowski functional tells us, how much we have to ’inflate’ our given set $K$ in order to capture the vector $x$. Clearly, from the definition we have that if $K$ is closed

[TABLE]

If a convex set $K$ is closed and symmetric, then $\|\cdot\|_{K}$ defines a norm on $\mathbb{R}^{n}$.

Recall that a set $K$ is star shaped, if there exists a point $x_{0}\in K$, which satisfies that for all $x\in K$ we have $\{tx_{0}+(1-t)x\colon t\in[0,1]\}\subset K$. It is easy to see that convex sets are star shaped, but for example unions of subspaces are not convex, but star shaped.

For bounded, star shaped $K$, the notion of $\|\cdot\|_{K}$ now allows to establish a direct correspondence between norm minimization problems and feasibility problems. With this observation, Corollary 3.1 translates to the following result.

Corollary** 3.2****.**

For $K$ bounded, symmetric and star-shaped, let $x\in K$ and $y=Ax$. Choose $\hat{x}\in\mathbb{R}^{n}$, such that it solves

[TABLE]

then

[TABLE]

Here $\hat{x}\in K$ is due to the fact that the minimum satisfies $\|\hat{x}\|_{K}\leq\|x\|_{K}\leq 1$, as $x\in K$ by assumption.

This result directly relates recovery guarantees to the mean width, it thus remains to calculate the mean width for the sets under consideration. In the following subsections, we will discuss two cases. The first one directly corresponds to the desired signal model, namely gradient sparse vectors. These considerations are mainly of theoretical interest, as the associated minimization problem closely relates to support size minimization, which is known to be NP hard in general. The second case considers the TV minimization problem introduced above, which then also yields guarantees for the (larger) set of vectors with bounded total variation.

Note, however, that the $M^{*}$-bound only gives a bound for the expected error. We can relate this result to a statement about tail probabilities using Markov’s inequality, namely

[TABLE]

In the next section we compute the mean width for the set of gradient sparse vectors, that is we now specify the set $K$ in Corollary 3.1 to be the set of all vectors with energy bounded by one that only have a small number of jumps.

3.2 The mean width of gradient sparse vectors in 1d

Here [PV13] served as an inspiration, as the computation is very similar for the set of sparse vectors.

Definition** 3.3****.**

The jump support of a vector $x$ is given via

[TABLE]

The jump support captures the positions, in which a vector $x$ changes its values. With this, we now define the set

[TABLE]

The set $K_{0}^{s}$ consists of all $s$-gradient sparse vectors, which have $2$-norm smaller than one. We will now calculate the mean width of $K_{0}^{s}$ in order to apply Corrolary 3.1 or 3.2.

Note that we can decompose the set $K_{0}^{s}$ into smaller sets $K_{J}\cap B_{2}^{n}$ with $K_{J}=\{x\colon\mathrm{Jsupp}(x)\subset J\}$, $|J|=s$ and $B_{2}^{n}=\{x\in\mathbb{R}^{n}\colon\|x\|_{2}\leq 1\}$. As we can’t add any jumps within the set $K_{J}$, it is a subspace of $\mathbb{R}^{n}$. We can even quite easily find an orthonormal basis for it, if we define

[TABLE]

As we can align all elements of $J=\{j_{1},j_{2},\ldots,j_{s}\}$ with $1\leq j_{1}<j_{2}<\ldots<j_{s}=n$, we see that $\{e_{[1,j_{1}]},e_{[j_{1}+1,j_{2}]},e_{[j_{2}+1,j_{3}]},\ldots,e_{[j_{s-1}+1,j_{s}]}\}$ forms an ONB of $K_{J}$. Now, we can write all elements $x\in K_{J}\cap B_{2}^{n}$ as $x=\sum_{i=1}^{s}\alpha_{i}e_{[j_{i-1}+1,j_{i}]}$ by setting $j_{0}:=0$. The property that $x\in B_{2}^{n}$ now enforces (ONB) that $\|\alpha\|_{2}\leq 1$. Now, note that $K_{0}^{s}=-K_{0}^{s}$, so we have

[TABLE]

Using the decomposition $K_{0}^{s}=\bigcup_{|J|=s}\left(K_{J}\cap B_{2}^{n}\right)$, we get

[TABLE]

Now

[TABLE]

Note that $G_{i}^{J}$ is again a Gaussian random variable with mean [math] and variance $1$. Furthermore, the supremum over $\alpha$ is attained, if $\alpha$ is parallel to $G^{J}$, so we have $\sup_{x\in K_{J}\cap B_{2}^{n}}\langle g,x\rangle=\|G^{J}\|_{2}$. Also note that $G^{J}$ has i.i.d. entries, but for different $J_{1},J_{2}$, the random vectors $G^{J_{1}}$ and $G^{J_{2}}$ may be dependent. Our task is now to calculate $\mathbb{E}\sup_{|J|=s}\|G^{J}\|_{2}$. As it has been shown for example in [FR13], we have that

[TABLE]

and from standard results for Gaussian concentration (cf. [PV13]), we get

[TABLE]

By noting that $|\{J\subset[n]\colon|J|=s\}|={n\choose s}$, we see by a union bound that

[TABLE]

For the following calculation, set $X:=\sup_{|J|=s}\|G^{J}\|_{2}$. By Jensen’s inequality and rewriting the expectation, we have that

[TABLE]

Now, the previous consideration showed, that

[TABLE]

Computing the resulting integrals yields

[TABLE]

Using a standard bound for the binomial coefficients, namely ${n\choose s}\leq\mathrm{e}^{s\log(\mathrm{e}n/s)}$, we see

[TABLE]

or equivalently

[TABLE]

By setting $\lambda=\sqrt{s\log(\mathrm{e}n/s)}$ and assuming (reasonably) large $n$, we thus get

[TABLE]

From this, we see that

[TABLE]

It follows that the Gaussian mean width of the set of gradient sparse vectors is the same as the mean width of sparse vectors due to the similar structure. If we want to obtain accuracy $\delta$ for our reconstruction, according to Theorem 3.1, we need to take

[TABLE]

measurements.

In Compressed Sensing, the squared mean width of the set of $s$-sparse vectors (its so called statistical dimension) already determines the number of required measurements in order to recover a sparse signal with basis pursuit. This is the case because the convex hull of the set of sparse vectors can be embedded into the $\ell_{1}$-ball inflated by a constant factor. In the case of TV minimization, as we will see in the following section, this embedding yields a (rather large) constant depending on the dimension.

3.3 The extension to gradient compressible vectors needs a new approach

In the previous subsection, we considered exactly gradient sparse vectors. However searching all such vectors $x$ that satisfy $Ax=y$ is certainly not a feasible task. Instead, we want to solve the convex program

[TABLE]

with $\|z\|_{TV}=\|\nabla z\|_{1}$ the total variation seminorm. Now if we have that $x\in K_{0}^{s}$, we get that

[TABLE]

with $\alpha$ as in section 3.2, so $K_{0}^{s}\subset K_{TV}^{2\sqrt{s}}:=\{x\in B_{2}^{n}\colon\|x\|_{TV}\leq 2\sqrt{s}\}$. As $K_{TV}^{2\sqrt{s}}$ is convex, we even have $\mathrm{ch}(K_{0}^{s})\subset K_{TV}^{2\sqrt{s}}$. We can think of the set $K_{TV}^{2\sqrt{s}}$ as ’gradient- compressible’ vectors.

In the proof of Theorem 3.3 in [CX15], the Gaussian width of the set $K_{TV}^{4\sqrt{s}}$ has been calculated via a wavelet-based argument. One obtains that $w(K_{TV}^{2\sqrt{s}})\leq C\sqrt{\sqrt{ns}\log(2n})$ with $C\leq 20$ being an absolute constant. In this section we illustrate, that proof techniques different from the ones used in the case of synthesis sparsity are indeed necessary in order to obtain useful results. In the synthesis case, the $1$-norm ball of radius $\sqrt{s}$ is contained in the set of $s$-sparse vectors inflated by a constant factor. This in turn implies that the mean width of the compressible vectors is bounded by a constant times the mean width of the $s$-sparse vectors.

We will attempt a similar computation, that is to find a constant, such that the set $K_{TV}^{2\sqrt{s}}$ is contained in the ’inflated’ set $c_{n,s}\mathrm{ch}(K_{0}^{s})$. Then $w(K_{TV}^{2\sqrt{s}})\leq c_{n,s}w(K_{0}^{s})$. Although this technique works well for sparse recovery, where $c_{n,s}=2$, it pityably fails in the case of TV recovery as we will see below.

Let us start with $x\in K_{TV}^{2\sqrt{s}}$. Now we can decompose $J:=\mathrm{Jsupp}(x)=J_{1}\uplus J_{2}\uplus\ldots J_{p}$ with $|J_{k}|\leq s$ in an ascending manner, i.e., for all $k\in J_{i},l\in J_{i+1}$, we have that $\alpha_{k}<\alpha_{l}$. Note that the number $p$ of such sets satisfies $p\leq\frac{n}{s}$. Similarly as above, we now write $x=\sum_{i=1}^{|J|}\alpha_{i}e_{[j_{i-1}+1,j_{i}]}=\sum_{k=1}^{p}\sum_{i\in J_{k}}\alpha_{i}e_{[j_{i-1}+1,j_{i}]}$. From this, we see that

[TABLE]

The necessary factor $c_{n,s}$ can be found by bounding the size of $\|\alpha_{J_{k}}\|_{2}$, namely

[TABLE]

From this, we see that $K_{TV}^{2\sqrt{s}}\subset\sqrt{\frac{n}{s}}\mathrm{ch}(K_{0}^{s})$. To see that this embedding constant is optimal, we construct a vector, for which it is needed.

To simplify the discussion, suppose that $n$ and $s$ are even and $s|n$. For even $n$, the vector $x_{1}=(\sqrt{\frac{1-(-1)^{k}\varepsilon}{n}})_{k}$ has unity norm, lies in $K_{TV}^{2\sqrt{s}}$ for $\varepsilon<\frac{2\sqrt{s}}{n}$ and has jump support on all of $[n]$!

For a vector $x\in\mathbb{R}^{n}$ and an index set $I\subset[n]$, we define the restriction of $x$ to $I$ by

[TABLE]

By splitting $\mathrm{Jsupp}(x_{1})$ into sets $J_{1},\ldots,J_{n/s}$ and setting $a_{k}=\sqrt{\frac{n}{s}}x_{1}|_{J_{k}}\in K_{0}^{s}$, we see that $x_{1}=\sum_{k=1}^{n/s}\sqrt{\frac{s}{n}}a_{k}$ and in order for this to be elements of $c_{n,s}\mathrm{ch}(K_{0}^{s})$, we have to set $c_{n,s}=\sqrt{\frac{n}{s}}$. This follows from

[TABLE]

and no smaller inflation factor than $\sqrt{\frac{n}{s}}$ can suffice.

So from the previous discussion, we get

Lemma** 3.1****.**

We have the series of inclusions

[TABLE]

In view of the results obtainable for sparse vectors and the $\ell_{1}$-ball, this is very disappointing, because Lemma 3.1 now implies that the width of $K_{TV}^{2\sqrt{s}}$ satisfies

[TABLE]

which is highly suboptimal.

Luckily, the results in [CX15] suggest, that the factor $n$ in the previous equation can be replaced by $\sqrt{sn}$. However, they have to resort to a direct calculation of the Gaussian width of $K_{TV}^{2\sqrt{s}}$. The intuition why the Gaussian mean width can be significantly smaller than the bound given in Lemma 3.1 stems from the fact, that in order to obtain an inclusion we need to capure all ’outliers’ of the set - no matter how small their measure is.

3.4 Exact recovery

For exact recovery, the $M^{*}$-bound is not suitable anymore and, as suggested in [Ver15], we will use ’Gordon’s escape through the mesh’ in order to find conditions on exact recovery. Exact recovery for TV minimization via this approach has first been considered in [CX15].

Suppose, we want to recover $x\in K_{0}^{s}$ from Gaussian measurements $Ax=y$. Given, that we want our estimator $\hat{x}$ to lie in a set $K$, exact recovery is achieved, if $K\cap\{z\colon Az=y\}=\{x\}$. This is equivalent to requiring

[TABLE]

With the descent cone $D(K,x)=\{t(z-x)\colon t\geq 0,z\in K\}$, we can rewrite this condition as

[TABLE]

by introducing the set $S(K,x)=D(K,x)\cap B_{2}^{n}$, we see that if

[TABLE]

we get exact recovery. The question, when a section of a subset of the sphere with a random hyperplane is empty is answered by Gordon’s escape through a mesh.

Theorem** 3.2**** ([Gor88]).**

Let $S\subset\mathbb{S}^{n-1}$ be fixed and $E\in G_{n,n-m}$ be drawn at random according to the Haar measure. Assume that $\hat{w}(S)=\mathbb{E}\sup_{u\in S}\langle g,u\rangle<\sqrt{m}$, then $S\cap E=\emptyset$ with probability exceeding

[TABLE]

So we get exact recovery with high probability from a program given in Theorem 3.1 or 3.2, provided that $m>\hat{w}(S(K,x_{0}))^{2}$.

Let’s see how this applies to TV minimization. Suppose, we are given $x\in K_{0}^{s}$ and Gaussian measurements $Ax=y$. Solving

[TABLE]

amounts to using the Minkowski functional of the set $K=\{z\in\mathbb{R}^{n}\colon\|z\|_{TV}\leq\|x\|_{TV}\}$, which is a scaled TV-Ball.

In [CX15], the null space property for TV minimization given in Corollary 2.1 has been used in order to obtain recovery guarantees.

They consider the set, where this condition is not met

[TABLE]

and apply Gordon’s escape through the mesh to see that with high probability, its intersection with the kernel of $A$ is empty, thus proving exact recovery with high probability. Their estimate to the mean width of the set $\mathcal{S}$,

[TABLE]

with $c<19$ is essentially optimal (up to logarithmic factors), as they also show that $w(\mathcal{S})\geq C\sqrt[4]{ns}$. So uniform exact recovery can only be expected for $m=\mathcal{O}(\sqrt{sn}\log n)$ measurements.

Let us examine some connections to the previous discussion about the descent cone.

Lemma** 3.2****.**

We have that for $K=\{z\in\mathbb{R}^{n}\colon\|z\|_{TV}\leq\|x\|_{TV}\}$ defined as above and $x\in K_{0}^{s}$, it holds that $S(K,x)\subset\mathcal{S}$.

Proof.

Let $y\in S(K,x)$. Then there exists a $x\neq z\in K$, such that $y=\frac{z-x}{\|z-x\|_{2}}$. Set $J=\mathrm{Jsupp}(x)$, then, as $z\in K$, we have that $\|z\|_{TV}\leq\|x\|_{TV}$, or

[TABLE]

Now, by the triangle inequality and this observation, we have

[TABLE]

The last equality follows from the fact that $\nabla x$ is zero outside of the gradient support of $x$. Multiplying both sides with $\frac{1}{\|z-x\|_{2}}$ gives the desired result

[TABLE]

∎

The previous lemma shows that the recovery guarantees derived from the null space property and via the descent cone are actually connected in a very simple way.

Clearly, now if we do not intersect the set $\mathcal{S}$, we also do not intersect the set $S(K,x)$, which yields exact recovery for example with the same upper bounds on $m$ as for $\mathcal{S}$. Even more specifically, in the calculation of $\hat{w}(\mathcal{S})$ given in [CX15], an embedding into a slightly larger set $\tilde{\mathcal{S}}=\{x\in B_{2}^{n}\colon\|x\|_{TV}\leq 4\sqrt{s}\}$ is made. This embedding can also quite easily be done if we note that $\|x\|_{TV}\leq 2\sqrt{s}$, as we showed above and $\|z\|_{TV}\leq\|x\|_{TV}$.

Note that the same discussion also holds for higher dimensional signals, such that the improved numbers of measurements as given in Section 2.2 can be applied.

3.5 Subgaussian measurements

Up to this point, all our measurement matrices have been assumed to consist of i.i.d. Gaussian random variables. We will reduce this requirement in this section to be able to incorporate also subgaussian measurement matrices into our framework.

Definition** 3.4****.**

A real valued random variable $X$ is called subgaussian, if there exists a number $t>0$, such that $\mathbb{E}\mathrm{e}^{tX^{2}}<\infty$. A real valued random vector is called subgaussian, if all of its one dimensional marginals are subgaussian.

An obvious example of subgaussian random variables are Gaussian random variables, as the expectation in Definition 3.4 exists for all $t<1$. Also, all bounded random variables are subgaussian.

Here, we rely on results given by Tropp in [Tro15] using the results of Mendelson [KM15, Men14]. We will consider problems of the form

[TABLE]

where $A$ is supposed to be a matrix with independent subgaussian rows. Furthermore, we denote the exact solution by $x_{0}$, i.e., $Ax_{0}=y$. We pose the following assumptions on the distribution of the rows of $A$.

1. (M1)

$\mathbb{E}A_{i}=0$, 2. (M2)

There exists $\alpha>0$, such that for all $u\in\mathbb{S}^{n-1}$ it holds that $\mathbb{E}|\langle A_{i},u\rangle|\geq\alpha$, 3. (M3)

There is a $\sigma>0$, such that for all $u\in\mathbb{S}^{n-1}$ it holds that $\mathbb{P}(|\langle A_{i},u\rangle|\geq t)\leq 2\exp(-t^{2}/(2\sigma^{2}))$, 4. (M4)

The constant $\rho:=\frac{\sigma}{\alpha}$ is small.

Then the small ball methods yields the following recovery guarantee (we present the version of [Tro15]).

Theorem** 3.3****.**

Let $x_{0}\in\mathbb{R}^{n}$ and $A\in\mathbb{R}^{m\times n}$ be a subgaussian matrix satisfying (M1)-(M4) above. Furthermore, set $y=Ax_{0}+e$, where $\|e\|\leq\varepsilon$ denotes the (bounded) error of the measurement. Then the solution $\hat{x}$ of (12) satisfies

[TABLE]

with probability exceeding $1-\mathrm{e}^{-ct^{2}}$. $D(K,x_{0})$ denotes the descent cone of the set $K$ at $x_{0}$, as defined in the previous section.

From this we see that, provided

[TABLE]

we obtain stable reconstruction of our original vector from (12). Note that the theorem is only meaningful for $t=\mathcal{O}(\sqrt{m})$, as otherwise the denominator vanishes.

In the previous section, we have shown the inclusion $S(K,x_{0})\subset\mathcal{S}$ for $x_{0}\in K_{s}^{0}$ and hence we have that

[TABLE]

So we see that for $m\geq\tilde{C}\rho^{6}\sqrt{ns}\log(2n)$, we obtain the bound

[TABLE]

with high probability. We conclude that, given the absolute constants $c,C$, we need to set $\tilde{C}\geq\tfrac{C^{2}}{c^{2}}$ in order to obtain a meaningful result. Combining all our previous discussions with Theorem 3.3, we get

Theorem** 3.4****.**

Let $x_{0}\in\mathbb{R}^{n}$, $m\geq\tilde{C}\rho^{6}\sqrt{ns}\log(2n)$ and $A\in\mathbb{R}^{m\times n}$ be a subgaussian matrix satisfying (M1)-(M4). Furthermore, set $y=Ax_{0}+e$, where $\|e\|\leq\varepsilon$ denotes the (bounded) error of the measurement, constants $c,C,\tilde{C}>0$ as above and $t\leq\frac{\sigma(c\sqrt{\tilde{C}}-C)\sqrt[4]{ns}\sqrt{\log(2n)}}{\alpha}$. Then the solution $\hat{x}$ of (12) satisfies

[TABLE]

We can for example set $t=\rho(c\sqrt{\tilde{C}}-C)\sqrt[4]{ns}$ (for $n\geq 2$) to obtain the bound

[TABLE]

For example for i.i.d. standard Gaussian measurements, the constant $\rho=\sqrt{\frac{2}{\pi}}$.

Note that in the case of noisefree measurements $\varepsilon=0$, Theorem 3.4 gives an exact recovery result for a wider class of measurement ensembles with high probability. Furthermore with a detailed computation of $w(S(K,x_{0}))$ one may be able to improve the number of measurements for nonuniform recovery. It also remains open, whether the lower bounds of Cai and Xu for the case of Gaussian measurements can be generalized to the subgaussian case. In fact, our numerical experiments summarized in Figure 2 suggest a better scaling in the ambient dimension, around $N^{1/4}$, in the average case. We consider it an interesting problem for future work to explore whether this is due to a difference between Rademacher and Gaussian matrix entries, between uniform and nonuniform recovery, or between the average and the worst case. Also, it is not clear whether the scaling is in fact $N^{1/4}$ or if the observed slope is just a linearization of, say, a logarithmic dependence.

4 Discussion and open problems

As the considerations in the previous sections illustrate, the mathematical properties of total variation minimization differ significantly from algorithms based on synthesis sparsity, especially in one dimension. For this reason, there are a number of questions that have been answered for synthesis sparsity, but which are still open for the framework of total variation minimization. For example, the analysis provided in [RRT12, KMR14] for deterministically subsampled partial random circulant matrices, as they are used to model measurement setups appearing in remote sensing or coded aperture imaging, could not be generalized to total variation minimization. The difficulty in this setup is that the randomness is encoded by the convolution filter, so it is not clear what the analogy of variable density sampling would be.

Another case of practical interest is that of sparse $0/1$ measurement matrices. Recently it has been suggested that such meausurements increase efficiency in photoacoustic tomography, while at the same time, the signals to be recovered (after a suitable temporal transform) are approximately gradient sparse. This suggests the use of total variation minimization for recovery, and indeed empirically, this approaches yields good recovery results [SKB*+*15]. Theoretical guarantees, however, (as they are known for synthesis sparse signals via an expander graph construction [BGI*+*08]) are not available to date for this setup.

Acknowledgements

FK and MS acknowledge support by the Hausdorff Institute for Mathematics (HIM), where part of this work was completed in the context of the HIM trimester program ”Mathematics of Signal Processing”, FK and CK acknowledge support by the German Science Foundation in the context of the Emmy Noether Junior Research Group “Randomized Sensing and Quantization of Signals and Images” (KR 4512/1-1) and by the German Ministry of Research and Education in the context of the joint research initiative ZeMat. MS has been supported by the Austrian Science Fund (FWF) under Grant no. Y760 and the DFG SFB/TRR 109 ”Discretization in Geometry and Dynamics”.