Performance Bounds For Co-/Sparse Box Constrained Signal Recovery

TL;DR

This paper establishes performance bounds for recovering co-/sparse signals with box constraints using convex regularization, focusing on TV minimization and its application to compressed sensing and tomography.

Contribution

It introduces dual certificates and statistical dimension analysis to estimate undersampling rates for TV-based recovery, filling a gap in theoretical understanding.

Findings

Derived dual certificates for uniqueness testing.

Estimated undersampling rates via statistical dimension.

Validated bounds empirically for tomographic measurements.

Abstract

The recovery of structured signals from a few linear measurements is a central point in both compressed sensing (CS) and discrete tomography. In CS the signal structure is described by means of a low complexity model e.g. co-/sparsity. The CS theory shows that any signal/image can be undersampled at a rate dependent on its intrinsic complexity. Moreover, in such undersampling regimes, the signal can be recovered by sparsity promoting convex regularization like $\ell_1$- or total variation (TV-) minimization. Precise relations between many low complexity measures and the sufficient number of random measurements are known for many sparsity promoting norms. However, a precise estimate of the undersampling rate for the TV seminorm is still lacking. We address this issue by: a) providing dual certificates testing uniqueness of a given cosparse signal with bounded signal values, b) approximating the undersampling rates via the statistical dimension of the TV descent cone and c) showing empirically that the provided rates also hold for tomographic measurements.

Equations

$$A\overline{x}=b,$$

$$\underset{x}{\min }f(x)\mathrm{subject}\text{leavevmode}\mathrm{to}Ax=A\overline{x},$$

$$f_1(x)$$

$$f_2(x)$$

$$f_3(x)$$

$$f_4(x)$$

$$f_5(x)$$

$$f_6(x)$$

$$\nabla =\left(\begin{array}{c}{I}_n\otimes {\partial }_m\\ {\partial }_n\otimes I_m\end{array}\right)\in {\mathbb{R}}^{p\times n},$$

$$\min f(x)\text{subject to}Ax=b⟺\underset{x\in {\mathbb{R}}^n}{\min }g(x),g(x):=f(x)+{\delta }_{\mathcal{X}}(x),$$

$$\mathcal{X}:=\{x\in {\mathbb{R}}^n:Ax=b\},$$

$$b:=A\overline{x}.$$

$$\mathrm{argmin}g=\{\overline{x}\}\iff 0\in \mathrm{int}\text{leavevmode}\partial g(\overline{x}).$$

$$\underset{x\in {\mathbb{R}}^n}{\min }\Vert Dx{\Vert }_1+{\delta }_{\mathcal{X}}(x)+{\delta }_{[l,u]}(x),$$

$$\partial \left(\Vert D\overline{x}{\Vert }_1+{\delta }_{\mathcal{X}}(\overline{x})+{\delta }_{{[l,u]}^n}(\overline{x})\right)=\partial \Vert D\overline{x}{\Vert }_1+\partial {\delta }_{\mathcal{X}}(\overline{x})+\partial {\delta }_{{[l,u]}^n}(\overline{x}).$$

$$\partial \Vert D\overline{x}{\Vert }_1=\{D^{\top }\alpha :{\alpha }_{{\Lambda }^c}=\mathrm{sign}(D_{{\Lambda }^c,\bullet }\overline{x}),\Vert {\alpha }_{\Lambda }{\Vert }_{\infty }\le 1\},$$

$$\partial {\delta }_{\mathcal{X}}(\overline{x})=N_{\mathcal{X}}(\overline{x})=\mathcal{N}(A{)}^{\perp },$$

$$\partial {\delta }_{[l,u]}(x)=N_{[l,u]}(\overline{x}).$$

$$N_{[l,u]}(\overline{x})=N_{[l_1,u_1]}({\overline{x}}_1)\times \cdots \times N_{[l_n,u_n]}({\overline{x}}_n)={\mathbb{R}}_{-}^{|S_l|}\times \{0{\}}^{|S_{lu}|}\times {\mathbb{R}}_{+}^{|S_u|},$$

$$N_{[l_i,u_i]}({\overline{x}}_i)=\{\begin{array}{ll}0,& {\overline{x}}_i\in (l_i,u_i),\\ (-\infty ,0],& {\overline{x}}_i=l_i,\\ [0,\infty ),& {\overline{x}}_i=u_i,\end{array}$$

$$\partial {\delta }_{{\mathbb{R}}_{+}^n}(\overline{x})=N_{{\mathbb{R}}_{+}^n}(\overline{x})=N_{{\mathbb{R}}_{+}}({\overline{x}}_1)\times \cdots \times N_{{\mathbb{R}}_{+}}({\overline{x}}_n)=\{0{\}}^{|S|}\times {\mathbb{R}}_{-}^{|S^c|},$$

$${\Psi }_{ii}=\{\begin{array}{ll}-1,& {\overline{x}}_i=l_i,\\ 1,& {\overline{x}}_i=u_i,\\ 0,& \text{otherwise.}\end{array}$$

$$0\in \mathrm{int}\left(\partial \left(\Vert D\overline{x}{\Vert }_1+{\delta }_{\mathcal{X}}(\overline{x})+{\delta }_{{[l,u]}^n}(\overline{x})\right)\right).$$

$${\mathbb{R}}^n=\mathrm{aff}(\partial \Vert D\overline{x}{\Vert }_1)+\mathrm{aff}(\partial {\delta }_{\mathcal{X}}(\overline{x}))+\mathrm{aff}(\partial {\delta }_{[l,u]}(\overline{x})).$$

$$\displaystyle\operatorname{aff}(\partial\|D\overline{x}\|_{1})=\$$

$$\displaystyle=\$$

$$\mathrm{aff}(\partial {\delta }_{\mathcal{X}}(\overline{x}))=\mathrm{aff}(N_{\mathcal{X}}(\overline{x}))=\mathrm{aff}(\mathcal{N}(A{)}^{\perp })=\mathcal{N}(A{)}^{\perp }.$$

$$\displaystyle\operatorname{aff}(\partial\delta_{[l,u]}(\overline{x}))=\$$

$$\displaystyle=\$$

$$\mathcal{R}(D_{\Lambda ,\bullet }^{\top }{)}^{\perp }\cap (\mathcal{N}(A{)}^{\perp }{)}^{\perp }\cap \mathcal{R}(\Psi {)}^{\perp }$$

$$\displaystyle=\$$

$$\mathrm{rint}(C_1+C_2)=\mathrm{rint}(C_1)+\mathrm{rint}(C_2),$$

$$0=$$

$$\displaystyle=\$$

Full text

Jan Kuske and Stefania Petra

††footnotetext: † This is an invited paper to the special issue of Analele Stiintifice ale Universitatii Ovidius Constanta (ASUOC) on the occasion of the 11th Workshop on Mathematical Modelling of Environmental and Life Sciences Problems organised in October 12-16, 2016, Constanta, Romania (http://www.ima.ro/mmelsp.htm). It was accepted for publication in 2017 but not yet published.

The authors would like to stress that this is not a recent paper and is not referencing other works that have appeared since 2017.

1. 1 Introduction

2. 1.1 Contribution, Related Work, Organization

3. 1.2 Notation

4. 2 Individual Uniqueness and Dual Certificates

5. 2.1 Uniqueness of a Minimizer of a Polyhedral Function

6. 2.2 Certifying Recovery of Individual Vectors via Dual Certificates

7. 2.3 Case Studies: (1.3a)-(1.3f)

8. 3 Probabilistic Uniqueness

9. 3.1 Phase Transitions in Random Linear Inverse Problems

10. 3.2 Phase Transitions for Case Studies: (1.3a)-(1.3f)

11. 4 Experiments

12. 4.1 1D Empirical Phase Transitions and Theoretical Bounds

13. 4.1.1 Testset

14. 4.1.2 $\ell_{1}$-Minimization

15. 4.1.3 $\operatorname{TV}$-Minimization

16. 4.2 2D Empirical Phase Transitions and Theoretical Bounds

17. 4.2.1 Testset

18. 4.2.2 Tomographic Measuremnts

19. 4.2.3 $\ell_{1}$-Minimization

20. 4.2.4 $\operatorname{TV}$-Minimization

21. 5 Conclusions

1. Introduction

The current work is motivated by the image reconstruction problem from few linear measurements, e.g. discrete tomography [HK99], where typically bounds on the image values are available. This explains why we are concerned with the recovery of box constrained signals. Moreover, we aim to reconstruct an image that exhibits some structure, e.g. sparsity/gradient-sparsity.

Recovering a structured signal/image from few linear measurements is a central point in both compressed sensing (CS) [FR13] and discrete tomography [HK99]. In CS the signal structure is described by means of a low complexity model. For instance, signal sparsity $s:=\|x\|_{0}$, that equals the number of non-zero components in $x$, and the associated $\ell_{0}$-minimization is typically relaxed to $\ell_{1}$-minimization. The theory of CS implies that if the number of measurements $m\geq C\cdot s\log(n/s)$ and the entries of the measurement matrix $A$ follow a normal distribution, i.e. $a_{ij}\sim\mathcal{N}(0,1)$, the solution of $\ell_{0}$-minimization is equal to the solution of $\ell_{1}$-minimization. Thus, a combinatorial problem can be solved by a convex optimization problem.

In the present work we consider the noiseless setting

[TABLE]

where the unknown (structured) signal $\overline{x}\in\mathbb{R}^{n}$ is undersampled either by a Gaussian matrix $A\in\mathbb{R}^{m\times n}$, i.e. $A$ has independent standard normal entries, or by a tomographic projection matrix $A$.

For the reconstruction of the sparse or gradient-sparse box constrained signal of interest we consider some structure enforcing regularizer $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$, that is a proper convex function, and solve

[TABLE]

where $f$ is specialized to one of the following functions

[TABLE]

In the first part of the paper $D\in\mathbb{R}^{p\times n}$ is an arbitrary matrix, while in the latter part $D$ is specialized to the discrete gradient operator, in view of the (anisotropic) $\operatorname{TV}$-regularizer we are interested in.

1.1. Contribution, Related Work, Organization

Our main objectives are:

*Individual Uniqueness: *

We derive practical tests for individual uniqueness: Given a signal $\overline{x}$ and an instance of the optimization problem above, can we efficiently verify whether the signal is the unique solution? Our approach will be to apply the uniqueness condition $0\in\operatorname{int}{\partial f(\overline{x})}$ [Gil17] to several types of polyhedral functions $f$ including $\ell_{1}$- and $\operatorname{TV}$-minimization with and without box constraints, as detailed in Sect. 2. Such testable uniqueness conditions for $\operatorname{TV}$-minimization were also developed in [JKL15, ZMY16], but these use different mathematical tools. Our derivation is shorter and additionally considers box constrained $\operatorname{TV}$-minimization.

*Probabilistic Uniqueness: *

We wish to accurately describe the relation between the signal sparsity and the number of measurements that guarantee uniqueness with high probability. We build on recent CS theory [ALMT14, OH16, DJM13, FM14] that upper bounds the statistical dimension of the descent cone of the structure enforcing regularization. The difficulty lies in the fact that certain CS guarantees are missing for $\operatorname{TV}$-minimization. Moreover, the CS recovery performance bounds based on the statistical dimension (also called phase transitions) do not apply directly to tomography, since tomographic measurements are known to be ill-conditioned on the common structured signal classes from CS. Phase transitions in CS are discussed in Section 3.

*Empirical validation: *

We verify the theoretical probabilistic results for Gaussian matrices along with the theoretical uniqueness tools and compare them to the recovery and uniqueness performance of typical tomographic measurements. We stress that tomographic projections fall short of the common assumptions (e.g. restricted isometry property) underlying compressed sensing, see [PS14], and CS recovery guarantees do not apply for tomography. In view of its practical importance (reducing radiation exposure) we emphasize the need of accurately describing how much undersampling is tolerable in tomographic recovery. We empirically validate our theoretical results in Sect. 4. Hence, our work is closely related to the work in [JKL15], but additionally provides approximations to the theoretical phase transition curves that are validated empirically.

1.2. Notation

For $n\in\mathbb{N}$, we use the shorthands $[n]=\{1,2,\dotsc,n\}$. We define $\mathbb{R}_{+}^{n}:=\{x\in\mathbb{R}^{n}\colon x_{i}\geq 0,\forall i\in{[n]}\}$ and $\mathbb{R}_{++}^{n}:=\{x\in\mathbb{R}^{n}\colon x_{i}>0,\forall i\in{[n]}\}$. Analogously, we define $\mathbb{R}_{-}^{n}$ and $\mathbb{R}_{--}^{n}$. The extended real line is denoted by $\overline{\mathbb{R}}:=\mathbb{R}\cup\{\pm\infty\}$. Vectors $x\in\mathbb{R}^{n}$ are column vectors and indexed by superscripts. $x^{\top}$ denotes the transposed vector $x$ and $\langle x,y\rangle$ or $x^{\top}y$ the Euclidean inner product. $\top$ stands for the transpose. For a vector $x\in{\mathbb{R}}^{n}$, $\operatorname{sign}(x)\in{\mathbb{R}}^{n}$ is the sign vector of $x$ defined component-wise as ${\operatorname{sign}}(x)_{i}=1$ if $x_{i}>0$ and ${\operatorname{sign}}(x)_{i}=-1$ if $x_{i}<0$. For some matrix $A\in\mathbb{R}^{m\times n}$, the nullspace is denoted by $\mathcal{N}(A)$ and its range by $\mathcal{R}(A)$. For some matrix $\Omega\in\mathbb{R}^{p\times n}$ we denote the submatrix indexed by the rows in $\Lambda\subset{[p]}$ by ${\Omega}_{\Lambda,\bullet}$. Similarly, for some matrix $A\in\mathbb{R}^{m\times n}$ the columns indexed by $S\subset{[n]}$ are denoted by ${A}_{\bullet,S}$. The complement of an index set $S\subset{[n]}$ will be denoted by $S^{c}:={[n]}\setminus S$. For a subspace $X$, $X^{\perp}$ will denote its orthogonal complement. The affine hull of a set $X$, see Appendix, is denoted by $\operatorname{aff}(X)$. The conic hull is denoted by $\operatorname{cone}X$. We use the $\ell_{1}$-norm $\|x\|_{1}=\sum_{i\in{[n]}}|x_{i}|$ and the maximum norm $\|x\|_{\infty}:=\max_{i\in{[n]}}|x_{i}|$. The $\ell_{0}$-measure (not a norm!) stands for the cardinality of the support of $x$, i.e. $\|x\|_{0}:=|\operatorname{supp}(x)|$, with $\operatorname{supp}(x)=\{i\in{[n]}\colon x_{i}\neq 0\}$. The cosupport of $x$ is denoted by $\operatorname{cosupp}(x)=\{i\in{[n]}\colon x_{i}=0\}$.

We define the one-dimensional discrete derivative operator ${\partial_{n}}\colon\mathbb{R}^{n}\to\mathbb{R}^{n-1}$, $(\partial_{n})_{ij}=-1$ if $i=j$, $(\partial_{n})_{ij}=1$ if $j=i+1$ and $(\partial_{n})_{ij}=0$ otherwise. For $m\times n$ images the discrete gradient operator is defined as

[TABLE]

where $\otimes$ denotes the Kronecker product and $I_{m},I_{n}$, are identity matrices. The anisotropic discretized total variation (TV) is given by $\operatorname{TV}(x):=\|\nabla x\|_{1}$. The image gradient sparsity is given by $\|\nabla x\|_{0}$.

We further denote the indicator function of a convex set as $\delta_{C}:\mathbb{R}^{n}\to\overline{\mathbb{R}}$. Recall, $\delta_{C}(x)=0$ if $x\in C$ and $\delta_{C}(x)=\infty$, otherwise. The subdifferential of a function $f$ at $x$ will be denoted by $\partial f(x)$, see Appendix for the definition. We will use that the subdifferential of the indicator function equals the normal cone $\partial\delta_{C}(x)=N_{C}(x):=\{v\in\mathbb{R}^{n}\colon v^{\top}(y-x)\leq 0,\forall y\in C\}$. The distance of a vector $x$ to $C$ is denoted by $\operatorname{dist}(x,C)$ and defined as $\operatorname{dist}(x,C)=\min_{y\in C}\|x-y\|$. Finally, $\mathbb{E}(X)$ will denote the expectation of the random variable $X$.

2. Individual Uniqueness and Dual Certificates

In this section we derive testable uniqueness conditions for our six problems of interest. We will first provide some recovery conditions called dual certificates. More precisely, consider that we are given a specific vector $\overline{x}$ and we have to decide whether it is the unique solution of (1.2). We will provide necessary and sufficient conditions which certify the existence and uniqueness of a solution (1.2) for the case of a polyhedral function $f$. Our analysis closely follows [Gil17]. These conditions are formulated in terms of a solution to the dual problem of (1.2). For the cases enumerated in (1.3), we will see that it is possible to test uniqueness by simply solving a linear program.

In this section, it will be useful to recast (1.2) as an unconstrained optimization problem

[TABLE]

where $\delta_{\mathcal{X}}$ denotes the indicator function of the feasible set of (1.2) with

[TABLE]

$f$ the objective in (1.2) and

[TABLE]

Hence, the function $g$ considered in this section is convex polyhedral, meaning that its epigraph is a convex polyhedron.

2.1. Uniqueness of a Minimizer of a Polyhedral Function

We will use that the condition $0\in{\operatorname{int}}\leavevmode\nobreak\ \partial g(\overline{x})$ is equivalent to the uniqueness of the solution $\overline{x}$, in the case of a polyhedral convex function $g$. The following result gives us a necessary and sufficient condition for the uniqueness of a polyhedral function minimizer and was recently emphasized in [Gil17].

Lemma 2.1**.**

[Gil17, Lem 2.2]** Let $g\in\mathcal{F}_{0}(\mathbb{R}^{n})$ be any proper, convex and lsc function that is in addition polyhedral. Then

[TABLE]

We note: the necessity part uses polyhedrality, the sufficiency part is straightforward and does not require polyhedrality. We will see, in Sect. 3, Prop. 3.1, an alternative uniqueness condition for the problem (1.2). The advantage of (2.4) is that it directly allows the derivation of testable uniqueness conditions as detailed next.

2.2. Certifying Recovery of Individual Vectors via Dual Certificates

In this section we show how to decide if a given solution for (1.3) is unique. To this end, we derive uniqueness conditions for the problem

[TABLE]

where $l_{i}\in\mathbb{R}\cup\{-\infty\}$, $u_{i}\in\mathbb{R}\cup\{+\infty\}$, $l_{i}<u_{i}$, $i\in{[n]}$. Hence, (2.5) can be regarded as a generalization of (1.3a)-(1.3f). For example (1.3b) can be obtained from (2.5) by setting $D:=I_{n}$, $l:=0\in\mathbb{R}^{n}$ and $u_{i}:=+\infty$, $i\in{[n]}$. In order to apply the uniqueness criterion of Lem. (2.1) we will need to calculate the subdifferential of the objective function in (2.5). To compute the subdifferentials at a feasible point $\overline{x}$ one can apply the sum rule [RW09, Cor. 10.9]

[TABLE]

Since all individual components are polyhedral functions and $\overline{x}\in\operatorname{dom}\|D\cdot\|_{1}\cap{\mathcal{X}}\cap{{[l,u]}^{n}}$ holds by construction, we can simply add the subdifferentials of the individual terms

[TABLE]

The normal cone corresponding to the box constraints reads

[TABLE]

since

[TABLE]

where $S_{l}=\{i\in{[n]}\colon x_{i}=l_{i}\}$, $S_{u}=\{i\in{[n]}\colon x_{i}=u_{i}\}$ and $S_{lu}=(S_{l}\cup S_{u})^{c}=\{i\in{[n]}\colon x_{i}\in(l_{i},u_{i})\}$. If for example we specialize (2.5) to (1.3b), the normal cone of the nonnegative orthant becomes relevant and specializes to

[TABLE]

since $S_{u}=\emptyset$, $S_{lu}=\{i\in{[n]}\colon\overline{x}_{i}\in(0,+\infty)\}=\operatorname{supp}(\overline{x})=:S$ and $S_{l}=S^{c}$. Based on the considerations above we now derive the main result of this section.

Theorem 2.2** (Dual Certificate).**

A solution $\overline{x}$ to (2.5) is unique, if and only if the following conditions hold

• (i)

$\mathcal{N}(A)\cap\mathcal{N}({D}_{\Lambda,\bullet})\cap\mathcal{N}(\Psi)=\{0\},$**

• (ii)

$\exists\alpha\in\mathbb{R}^{p},\exists\mu\in\mathbb{R}^{n}:$* *

$D^{\top}\alpha+\Psi\mu\in\mathcal{R}(A^{\top}),\alpha_{\Lambda^{c}}=\operatorname{sign}({D}_{\Lambda^{c},\bullet}\overline{x}),\|\alpha_{\Lambda}\|_{\infty}<1,\mu>0,$**

where $\Lambda:=\operatorname{cosupp}(D\overline{x})$ and $\Psi=\Psi(\overline{x})$ is a diagonal matrix that depends on $\overline{x}$ with entries

[TABLE]

Proof.

By Lem. 2.1 the uniqueness of $\overline{x}$ as the unique solution of (2.5) is equivalent to

[TABLE]

We rewrite (2.12) in the form of the two conditions

• (a)

$\mathbb{R}^{n}=\operatorname{aff}\left(\partial\left(\|D\overline{x}\|_{1}+\partial\delta_{\mathcal{X}}(\overline{x})+\partial\delta_{{[l,u]}}(\overline{x})\right)\right)$,

• (b)

$0\in\operatorname{rint}\left(\partial\left(\|D\overline{x}\|_{1}+\partial\delta_{\mathcal{X}}(\overline{x})+\partial\delta_{{[l,u]}}(\overline{x})\right)\right)$.

Condition (a) ensures that the interior of the subdifferential is not-empty, and can be recast as

[TABLE]

By (2.7a) we get

[TABLE]

By (2.7b) we obtain

[TABLE]

And finally, by (2.8) we write

[TABLE]

Hence, (2.13) is equivalent to $\mathbb{R}^{n}=\mathcal{R}({D}_{\Lambda,\bullet}^{\top})+\mathcal{N}(A)^{\perp}+\mathcal{R}(\Psi)$. Taking the orthogonal, condition (a) now becomes

[TABLE]

that yields (i).

Further we reformulate condition (b). We use that for two convex sets $C_{1},C_{2}$ it holds

[TABLE]

see (5.2), and so we get

[TABLE]

This shows (ii) and concludes the proof. ∎

Remark 2.1**.**

By Thm. 2.2 we can show that $\overline{x}$ is the unique solution of (2.5) if and only if $(i)$ and $(ii)$ holds. Condition $(i)$ is tested by evaluating whether $\begin{pmatrix}A^{\top},&{D}_{\Lambda,\bullet}^{\top},&\Psi^{\top}\end{pmatrix}^{\top}$ has full column rank. For verifying $(ii)$ we need to test whether there exist $y\in\mathbb{R}^{m}$, $\alpha\in\mathbb{R}^{p}$, $\mu\in\mathbb{R}^{n}:$ $D^{\top}\alpha+\Psi\mu=A^{\top}y,\alpha_{\Lambda^{c}}=\operatorname{sign}({D}_{\Lambda^{c},\bullet}\overline{x}),\|\alpha_{\Lambda}\|_{\infty}<1,\mu>0.$

The second condition $(ii)$ could be verified in practice by minimizing $\|\alpha_{\Lambda}\|_{\infty}$ w.r.t. $y,\alpha,\mu$ while respecting the equality constraints and converting the strict inequality constraint $\mu>0$ to $\mu>=\varepsilon\mathbb{1}$ for a small $\varepsilon$, e.g. $\varepsilon=10^{-8}$, to ensure that the inequality is satisfied strictly. This results in a linear program (LP),

[TABLE]

For the optimal solution $(t^{\ast},y^{\ast},\alpha^{\ast},\mu^{\ast})$ we have by definition the smallest possible $t:=\|\alpha_{\Lambda}\|_{\infty}$. If $t^{\ast}$ is not smaller than one, then no $y$ exists with smaller $t$. We therefore declare $\overline{x}$ the unique minimizer if $t^{\ast}<1$, and if $t^{\ast}\geq 1$, $\overline{x}$ cannot be the unique minimizer. Numerically, we would test whether $t^{\ast}\leq 1-\varepsilon$. Technically, by applying the above procedure we risk rejecting a unique solution $\overline{x}$, for which $1-\varepsilon<t^{\ast}<1$. Recall that the choice for $\varepsilon$ is ad hoc.

In order to resolve this issue we provide next a theoretically well-founded methodology to deal with the strict feasibility problem

[TABLE]

The above feasibility problem can be recast as a linear system of inequalities, as we will see next. For this purpose we transform the above problem into the form $Mz=q,Pz<p$.

Theorem 2.3**.**

Let $M$ be a matrix so that $\mathcal{N}(M^{\top})=\{0\}$. Then there is a point $\overline{z}$ with $M\overline{z}=q,P\overline{z}<d$ if and only if $v=0$ is the only feasible solution of the problem

[TABLE]

Proof.

We consider the following pair of linear programs

\displaystyle\max\ $\displaystyle 0$

(P)

s.t. $\displaystyle Mz=q$

$\displaystyle Pz\leq d,$

\displaystyle\min\ $\displaystyle q^{\top}u+d^{\top}v$

(DP)

s.t. $\displaystyle M^{\top}u+P^{\top}v=0$

$\displaystyle v\geq 0$

where (P) is the primal problem and (DP) its dual.

First, note that under the assumption $\mathcal{N}(M^{\top})=\{0\}$, vector $v=0$ is the only feasible solution of (2.35) if and only if $v=0$ is the only solution of (DP). Indeed, the common constraint $M^{\top}u+P^{\top}v=0$ implies due to the assumptions $v=0$ and $\mathcal{N}(M^{\top})=\{0\}$ that $(u,v)=(0,0)$.

Hence, we can show the statement of the theorem by showing that (P) is strictly feasible if and only if $v=0$ is the only solution of (DP).

Assume on one hand that (P) is strictly feasible. Since any feasible solution is also a solution of (P) we can deduce from the existence of a primal solution the existence of a dual solution. The optimality conditions yield in particular $(d-P\overline{z})^{\top}v^{\ast}=0$ for such a primal solution $\overline{z}$ of (P) and any dual solution $(u^{\ast},v^{\ast})$ of (DP). Since $P\overline{z}<d$ holds, it follows that $v^{\ast}=0$ is the only solution.

Assume on the other hand that for each dual solution $(u^{\ast},v^{\ast})$ it holds $v^{\ast}=0$. In view of the strict complementarity condition [Van14, Thm 10.7] it exists a pair of primal and dual solutions such that

[TABLE]

Since by assumption $v^{\ast}=0$, (2.36) implies the existence of a primal solution $z^{\ast}$ with $Pz^{\ast}<d$.

This concludes the proof. ∎

Remark 2.2**.**

In view of Thm.  2.3 we can check feasibility of the system (2.31) – (2.34) by linear programming. We first associate to the system (2.31) – (2.34) with strict inequality constraints the primal linear program

[TABLE]

as in the Thm.  2.3. Assumption $\mathcal{N}(M^{\top})=\{0\}$ in Thm.  2.3 now corresponds to the first condition from Thm. 2.2. Assuming that it holds, we can now verify system (2.31) – (2.34) and strict feasibility in (P*<*) and if

[TABLE]

has an optimal objective value equal to zero.

2.3. Case Studies: (1.3a)-(1.3f)

In the previous section we provided verifiable uniqueness conditions for problem (2.5). In this section we specialize these conditions to verifiable uniqueness optimality conditions for our objective functions (1.3a)-(1.3f).

Corollary 2.4**.**

(Case $f_{1}$) A feasible point $\overline{x}$ of problem (1.2) with the objective function defined in (1.3a) is the unique solution if and only if

\begin{array}[]{c}\exists y\in\mathbb{R}^{m}:{A}_{\bullet,S}^{\top}y=\operatorname{sign}(\overline{x}_{S})\wedge\|{A}_{\bullet,S^{c}}^{\top}y\|_{\infty}<1\text{ and}\\ {A}_{\bullet,S}\text{ is injective,}\end{array}**

where $S:=\operatorname{supp}(\overline{x})$ and ${A}_{\bullet,S}$ selects the columns indexed by $S$.

Proof.

We derive the above claim from Thm.  2.2. In view of the definition in (1.3a) we set $D=I$ the identity matrix and $\Psi=0$ the zero matrix. The latter is due to the fact that $\mathbb{R}^{n}$ can be seen as $\delta_{[l,u]}(\overline{x})$ with $l_{i}=-\infty$, $u_{i}=\infty$, for all $i\in{[n]}$. Hence $\overline{x}$ always lies in the interior of $\mathbb{R}^{n}$ and $N_{\mathbb{R}^{n}}(\overline{x})=0$, compare (2.8). These choices specialize the conditions of Thm. 2.2 to the following conditions

• (i)

$\mathcal{N}(A)\cap\mathcal{N}({I}_{S^{c},\bullet})=\{0\}$,

• (ii)

$\exists\alpha\in\mathbb{R}^{n}:\alpha\in\mathcal{R}(A^{\top}),\alpha_{S}=\operatorname{sign}(\overline{x}_{S}),\|\alpha_{S^{c}}\|_{\infty}<1$,

taking into account also that $\Lambda=S^{c}$ and $\mathcal{N}(\Psi)=\mathbb{R}^{n}$. Condition (i) can be further simplified to

[TABLE]

which yields the condition that ${A}_{\bullet,S}$ has to be injective. We conclude the proof by transforming the condition (ii) to

[TABLE]

∎

Corollary 2.5**.**

(Case $f_{2}$) A feasible point $\overline{x}$ of problem (1.2) with the objective function defined in (1.3b) is the unique solution if and only if

\begin{array}[]{c}\exists y\in\mathbb{R}^{m}:{A}_{\bullet,S}^{\top}y=\mathbb{1}_{S}\wedge{A}_{\bullet,S^{c}}^{\top}y<\mathbb{1}_{S^{c}}\text{ and}\\ {A}_{\bullet,S}\text{ is injective,}\end{array}**

with $S:=\operatorname{supp}(\overline{x})$.

Proof.

As in the proof of Cor.  2.4 we apply Thm.  2.2 with $D=I$ and $\Psi_{ii}=0$, for $i\in S$ and $\Psi_{ii}=-1$, for $i\in S^{c}$, since $N_{R^{n}_{+}}(\overline{x})=\{0\}^{|S|}\times\mathbb{R}_{-}^{|S^{c}|}$. Conditions (i) and (ii) in Thm. 2.2 become

• (i)

$\mathcal{N}(A)\cap\mathcal{N}({I}_{S^{c},\bullet})\cap\mathcal{N}(\Psi)=\{0\}$,

• (ii)

$\exists\alpha\in\mathbb{R}^{n},\exists\mu\in\mathbb{R}^{n}:$

$\alpha+\Psi\mu\in\mathcal{R}(A^{\top}),\alpha_{S}=\operatorname{sign}(\overline{x}_{S}),\|\alpha_{S^{c}}\|_{\infty}<1,\mu>0$.

Since $\mathcal{N}({I}_{S^{c},\bullet})=\mathcal{N}(\Psi)$ condition (i) above simplifies as in the proof of Cor.  2.4 to condition: ${A}_{\bullet,S}$ is injective.

Since $\{\Psi\mu\colon\mu\in\mathbb{R}^{n}_{++}\}=\Psi R^{n}_{++}=\{0\}^{|S|}\times\mathbb{R}^{|S^{c}|}_{--}$ and

[TABLE]

condition (ii) above is equivalent to the existence of a vector ${y}\in{\mathbb{R}}^{m}$ such that

[TABLE]

that can be written as

[TABLE]

This completes the proof. ∎

Corollary 2.6**.**

(Case $f_{3}$) A binary feasible point $\overline{x}\in\{0,1\}^{n}$ of problem (1.2) with the objective function defined in (1.3c) is the unique solution if and only if

[TABLE]

with $S:=\operatorname{supp}(\overline{x})$.

Proof.

As above we specialize Thm.  2.2 with $D=I$ and $\Psi_{ii}=-1$, for $i\in S^{c}$ and $\Psi_{ii}=1$, for $i\in S$, since $N_{[0,1]^{n}}(\overline{x})=\mathbb{R}^{|S|}_{+}\times R_{-}^{|S^{c}|}$ for $\overline{x}\in\{0,1\}^{n}$. Hence $\mathcal{N}(\Psi)=\{0\}$. As a consequence condition (i) Thm.  2.2 is automatically fulfilled. We now focus on condition (ii) that now reads

[TABLE]

since for $\overline{x}\in\{0,1\}^{n}$ we have $\alpha_{S}=\operatorname{sign}(\overline{x}_{S})=\mathbb{1}$. In view of

[TABLE]

condition (2.42) is equivalent to the existence of a vector ${y}\in{\mathbb{R}}^{m}$ such that

[TABLE]

that can be written as

[TABLE]

∎

Corollary 2.7**.**

(Case $f_{4}$) A feasible point $\overline{x}\in\mathbb{R}^{n}$ of problem (1.2) with the objective function defined in (1.3d) is the unique solution if and only if

• (i)

$\mathcal{N}(A)\cap\mathcal{N}({D}_{\Lambda,\bullet})=\{0\},$**

• (ii)

$\exists\alpha:D^{\top}\alpha\in\mathcal{R}(A^{\top}),\alpha_{\Lambda^{c}}=\operatorname{sign}({D}_{\Lambda^{c},\bullet}\overline{x}),\|\alpha_{\Lambda}\|_{\infty}<1,$**

where $\Lambda=\operatorname{cosupp}(D\overline{x})$.

Proof.

As in the proof of Cor. 2.4 we conclude that $\Psi=0$ is the zero matrix. Hence, the two conditions of Thm.  2.2 simplify to the two conditions above, that do not involve the variable $\mu$. ∎

Corollary 2.8**.**

(Case $f_{5}$) A feasible point $\overline{x}\in\mathbb{R}^{n}_{+}$ of problem (1.2) with the objective function defined in (1.3e) is the unique solution if and only if

• (i)

$\mathcal{N}(A)\cap\mathcal{N}({D}_{\Lambda,\bullet})\cap\mathcal{N}(\Psi)=\{0\},$**

• (ii)

$\exists\alpha,\exists\mu:D^{\top}\alpha+\Psi\mu\in\mathcal{R}(A^{\top}),\alpha=\operatorname{sign}({D}_{\Lambda^{c},\bullet}\overline{x}),\|\alpha_{\Lambda}\|_{\infty}<1,\mu>0$**

where $\Lambda=\operatorname{cosupp}(D\overline{x})$ and

[TABLE]

Proof.

Immediate from Thm.  2.2 and the definition of $f_{5}$ in (1.3e) and $\Psi$. ∎

Corollary 2.9**.**

(Case $f_{6}$) A binary feasible point $\overline{x}\in\{0,1\}^{n}$ of problem (1.2) with the objective function defined in (1.3f) is the unique solution if and only if

• (i)

$\mathcal{N}(A)\cap\mathcal{N}({D}_{\Lambda,\bullet})\cap\mathcal{N}(\Psi)=\{0\},$**

• (ii)

$\exists\alpha,\exists\mu:D^{\top}\alpha+\Psi\mu\in\mathcal{R}(A^{\top}),\alpha=\operatorname{sign}({D}_{\Lambda^{c},\bullet}\overline{x}),\|\alpha_{\Lambda}\|_{\infty}<1,\mu>0$**

where $\Lambda=\operatorname{cosupp}(D\overline{x})$ and

[TABLE]

Proof.

In view of the definition of $f_{6}$ in (1.3f) and $\Psi$ above $\mathcal{N}(\Psi)=\{0\}$. As a consequence condition (i) Thm.  2.2 is automatically fulfilled. Condition (ii) from Thm.  2.2 is kept unchanged. ∎

3. Probabilistic Uniqueness

In this section we are concerned with probabilistic solution uniqueness when minimizing some structure enforcing regularizer $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ via (1.2) in the noiseless setting of (1.1). In particular we address the question of predicting uniqueness with high probability for a specific number of random linear measurements, when knowing only the solution co-/sparsity.

For $f=\|\cdot\|_{1}$ and a Gaussian matrix $A$ it is well-known that the recovery of a sparse vector $x$ depends only on its sparsity level, and is independent of the locations or values of the nonzero entries. There are precise relations between the sparsity $s$, ambient dimension $n$, and the number of samples $m$ that guarantees success of (1.2). Moreover, several authors have shown that there is a transition from absolute success to absolute failure, and they have accurately characterized the location of the threshold, also called phase transition, in the $\ell_{1}$-case and other norms, e.g. nuclear norm. For the $\operatorname{TV}$-seminorm such a complete analysis is still missing. As mentioned in the introduction, the authors in [ZXCL16] have recently shown that in the case of 1D TV-minimization, case (1.3d), the phase transition can be accurately described by an effective bound for the statistical dimension of a descent cone, which is based on the squared distance of a standard normal vector to the subdifferential of the objective function at the sought solution $\overline{x}$,

[TABLE]

We will explore next if the same relation that describes phase transitions for $\ell_{1}$-minimization and 1D $\operatorname{TV}$-minimization also hold for all our objective functions in (1.3a)-(1.3f).

3.1. Phase Transitions in Random Linear Inverse Problems

We collect here some recent results of convex signal reconstruction with a Gaussian sampling model and briefly explain how classical results for Gaussian processes lead to a sharp bound for the number of Gaussian measurements that suffice for exact recovery.

Definition 3.1**.**

(Statistical dimension, [ALMT14]) The statistical dimension $\delta(K)$ of a closed, convex cone $K\subset\mathbb{R}^{n}$ is the quantity

[TABLE]

where $\Pi_{K}$ is the Euclidean projection onto $K$ and $X$ is a standard normal vector, i.e. $X\sim\mathcal{N}(0,I_{n})$.

Definition 3.2**.**

(Descent cone) The descent cone of a proper convex function $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ at a point $x$ is

[TABLE]

i.e. the conic hull of the directions that do not increase $f$ near $\overline{x}$.

Definition 3.3**.**

(Gaussian width, [CRPW12]) The Gaussian width of $C\subset\mathbb{R}^{n}$ is

[TABLE]

where $X$ is a standard normal vector, i.e. $X\sim\mathcal{N}(0,I_{n})$.

Remark 3.1**.**

One has, see [ALMT14], the relationship

[TABLE]

A classic result for solution uniqueness, equivalent to the exact recovery of individual vectors [FR13, Thm. 4.35], is given next.

Proposition 3.1**.**

The vector $\overline{x}$ is the unique solution of the convex program (1.2) if and only if

[TABLE]

Condition (3.4) holds with high probability for a random matrix $A$ is equal to the probability that a randomly rotated convex cone, here $\mathcal{N}(A)$, shares a ray with a fixed convex cone $\mathcal{D}_{f}(\overline{x})$. This probability is bounded in terms of the statistical dimension $\delta(\mathcal{D}_{f}(\overline{x}))$ as stated next.

Theorem 3.2**.**

[ALMT14, Thm. II]** Fix a tolerance parameter $\varepsilon\in(0,1)$. Let $\overline{x}\in\mathbb{R}^{n}$ be a fixed vector, and let $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ be a proper convex function. Suppose $A\in\mathbb{R}^{m\times n}$ has independent standard normal entries, and let $b=A\overline{x}$. Then

• •

$m\leq\delta(\mathcal{D}_{f}(\overline{x}))-a_{n}\sqrt{n}$* $\Longrightarrow$ (1.2) succeeds with probability $\leq\varepsilon$;*

• •

$m\geq\delta(\mathcal{D}_{f}(\overline{x}))+a_{n}\sqrt{n}$* $\Longrightarrow$ (1.2) succeeds with probability $\geq 1-\varepsilon$,*

with $a_{n}:=\sqrt{8\log(4/\varepsilon)}$.

There is a useful tool for bounding the statistical dimension in terms of $\partial f(\overline{x})$. Interestingly, this bound is tight for some classes of $f$, e.g. norms.

Proposition 3.3**.**

[CRPW12, ALMT14]** Let $f:\mathbb{R}^{n}\to\overline{\mathbb{R}}$ be a proper convex function and let $x\in\mathbb{R}^{n}$. Assume that the subdifferential $\partial f(\overline{x})$ is nonempty, compact, and does not contain the origin. Define the function

[TABLE]

where $X\sim\mathcal{N}(0,I)$. We have

[TABLE]

Furthermore, the function $J$ is strictly convex, continuous at $\tau=0$, and differentiable for $\tau\geq 0$. It achieves its minimum at a unique point.

The above results provide a recipe to upper bound $\delta(\mathcal{D}_{f}(\overline{x}))$. This consists of the following steps.

• •

Compute $\partial f(\overline{x})$;

• •

For each $\tau>0$, compute $J(\tau)=\mathbb{E}[\operatorname{dist}(X,\tau\partial f(\overline{x})]$;

• •

Find the unique solution to $J^{\prime}(\tau)=0$.

We also have the following error bound.

Theorem 3.4**.**

[ALMT14, Thm. 4.3.]** Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be a norm on $\mathbb{R}^{n}$ and let $x\in\mathbb{R}^{n}\setminus\{0\}$. Then

[TABLE]

For the $\ell_{1}$-case the r.h.s. can be made very accurate for large sparsity parameters $s$ providing an accurate estimate of the statistical dimension of the $\ell_{1}$-descent cone. On the other hand, this error estimate is not very accurate when the sparsity $s$ is small. The work in [FM14] contains a bound that improves this result and also extends to more general functions with some additional properties.

Theorem 3.5**.**

[FM14, Prop. 1]** Suppose that, for $\overline{x}\neq 0$, the subdifferential $\partial f(\overline{x})$ satisfies a weak decomposability assumption, i.e.

[TABLE]

Then

[TABLE]

with $J(\tau)=\mathbb{E}[\operatorname{dist}(X,\tau\partial f(\overline{x})]$.

3.2. Phase Transitions for Case Studies: (1.3a)-(1.3f)

Above we discussed probabilistic uniqueness for the general problem

[TABLE]

In this section we discuss how and if we can specialize these results to predict uniqueness for our objective functions (1.3a)-(1.3f).

Remark 3.2**.**

In the case of $f_{1}(x)=\|x\|_{1}$ from (1.3a) and $f_{4}(x)=\|Dx\|_{1}$ from (1.3d) the statistical dimension of the descent cone $\mathcal{D}_{f}(\overline{x})$ can be well approximated in terms of the expected squared distance to the scaled subdifferential. The reason is that the subdifferential $\partial f(\overline{x})$ is compact, see (2.7a), and thus by Prop. 3.3 we know that $J(\cdot):\mathbb{R}_{+}\to\mathbb{R}_{+}$ has a unique minimizer. Moreover, $\min_{\tau\geq 0}J(\tau)$ provides a tight approximation in view of Thm. 3.4 and Thm. 3.5, since $f_{1}$ is a norm and $\partial f_{4}$ is weakly decomposable for the choice $D$ equal to the 1D finite difference operator, as shown in [ZXCL16]. For the 2D finite difference operator the weakly decomposable property has not been previously investigated.

Remark 3.3**.**

The above situation changes for the objectives

[TABLE]

from (1.3b), (1.3c), (1.3e) and (1.3f). It is clear that $\partial f(\overline{x})$ is not compact, for all above choices, compare (2.7a) and (2.7c). Thus, we have no guarantee that, $J$ has even a minimizer. Since $\partial f(\overline{x})$ is not weakly decomposable, even if $\min_{\tau\geq 0}J(\tau)$ exists, we are not guaranteed a tight upper bound on the statistical dimension $\delta(\mathcal{D}_{f}(\overline{x}))$ and consequently on the number of sufficient Gaussian measurements $m$ for exact recovery of $\overline{x}$ with high probability.

Explicit curves for $\min_{\tau}J(\tau)$ can be computed only in the case of $f_{1}$, $f_{2}$, $f_{3}$ along the lines of [ALMT14]. We skip the details here and illustrate this curves in Sect. 4. For $f_{4}$, $f_{5}$ and $f_{6}$ explicit curves for $J$ and for $\min_{\tau}J(\tau)$ are missing. We use an approximation to the upper bound $J(\tau)=\mathbb{E}\left[\operatorname{dist}^{2}(X,\tau\partial f(\overline{x}))\right]$ where $X\sim\mathcal{N}(0,I)$. To approximate the expected value we use a very large $k\in\mathbb{N}$ and calculate

[TABLE]

where each $X_{i}$ is sampled from the normal distribution $\mathcal{N}(0,I)$. Note that for computing each $\operatorname{dist}^{2}(X_{i},\tau\partial f(\overline{x}))$ one is faced with solving a quadratic program, i.e.

[TABLE]

Note that the constraints are linear since $\partial f(\overline{x})$ is described by linear equalities and inequalities, compare (2.7a) and (2.7c).

Based on the empirical results from Sect. 4 we conjecture.

Conjecture 3.1**.**

For the choices of $f$ in Remark 3.3 the statistical dimension is well approximated similar to (3.9) by the squared Euclidean distance to the scaled subdifferential, i.e.

[TABLE]

4. Experiments

Our objective in this section is to empirically verify whether phase transitions occur at

[TABLE]

which is a guaranteed upper bound to the statistical dimension $\delta(\mathcal{D}_{f}(\overline{x}))$ only in some of the considered cases. We note further, that precise expressions for $J(\tau)$ are only available for $\ell_{1}$-minimization and cases (1.3a)–(1.3c). For the $\operatorname{TV}$-minimization and cases (1.3d)–(1.3f) we construct phase diagrams by numerically estimating $\min_{\tau\geq 0}J(\tau)$ and comparing the resulting bounds to the average-case results for both solving the reconstruction problems (1.2) and certifying uniqueness by the procedure derived in Sect. 2.

4.1. 1D Empirical Phase Transitions and Theoretical Bounds

In this section we address $\ell_{1}$-minimization and one dimensional $\operatorname{TV}$-minimization without and with constraints. The 1D TV-regularizer reads $\|\partial_{n}x\|_{1}$, where $\partial_{n}\in\mathbb{R}^{n-1\times n}$ is the one-dimensional discrete derivative operator. Applying $D=\partial_{n}$ to a signal $x$ gives the the entries’ differences for consecutive indices

[TABLE]

4.1.1. Testset

Our testset consists of several randomly generated signals with specified relative sparsity $\rho$ ranging from $0.05$ to $0.95$ with step size $0.05$. For each of these relative sparsities we create a signal so that

[TABLE]

Hence, we cover almost the full range from highly sparse to dense signals. In addition, we create three different classes of signals: with real-valued, non-negative or binary entries, i.e. $x_{i}\in\mathbb{R}$, $x_{i}\geq 0$ are $x_{i}\in\{0,1\}$ for all $i\in{[n]}$.

Creating random real-valued, non-negative or binary signals of a given sparsity is immediate. To create signals that are sparse in a transformed domain we use an idea from of [NDEG13]. Instead of choosing the support of $\overline{x}$ uniformly at random, ones chooses a set $\Lambda\subseteq{[p]}$ in order to create a signal with ${D}_{\Lambda,\bullet}\overline{x}=0$. Having the subset $\Lambda$ and a randomly created vector $v$ with non-zero entries we obtain the desired signal by calculating

[TABLE]

Since $D\overline{x}$ represents differences of signal entries for consecutive indices, taking the absolute value of each signal entry will most likely not change the cosupport of the resulting signal. Therefore, we can use this method to generate non-negative signals. Binary signals can be obtained by using Algorithm 1.

4.1.2. $\ell_{1}$-Minimization

In this section we compare the existing tight upper bounds for the statistical dimension $\delta(\mathcal{D}_{f}(\overline{x}))$ for $f$ defined in (1.3a)–(1.3c) to the empirical phase transition obtained for the above set of signals.

To this end we first set the ambient dimension $n=100$ and generate 10 instances of a sparse 1D signal, see Figure 4.1, and a random Gaussian matrix $A\in\mathbb{R}^{m\times}$. For each pair of relative sparsity and given number $m$ of random linear measurements we verify if the signal is the unique solution of (1.2) with $f$ defined in (1.3a)–(1.3c). Uniqueness is tested as detailed in Sect. 2. We used Mosek111https://www.mosek.com/ to solve the optimization problems. The gray value plots in Figure 4.4, first row, show the empirical probability that a given 1D sparse signal is uniquely reconstructed by the convex relaxation. White means recovered and unique, and black non recovered. The red curves show $\min_{\tau\geq 0}J(\tau)$, plotted with Mathematica222https://www.wolfram.com/mathematica/, and separate these regions accurately.

4.1.3. $\operatorname{TV}$-Minimization

To generate the phase diagrams we proceed as above but use the gradient sparse 1D test signals, see Figure 4.2. Since explicit curves for $J$ are missing we use an approximation to the upper bound $J(\tau)=\mathbb{E}\left[\operatorname{dist}^{2}(X,\tau\partial f(\overline{x}))\right]$, where $X\sim\mathcal{N}(0,I)$. To approximate the expected value we use a large $k=10000$ and calculate

[TABLE]

where each $X_{i}$ is sampled from the normal distribution $\mathcal{N}(0,I)$. Additionally we know from [ALMT14, Lem C.1] that the minimum is unique and lies in the interval $\left[0,\frac{2\|\overline{x}\|_{2}}{b}\right]$ with $b\leq\|w\|_{2}$ for all $w\in\partial f(\overline{x})$, provided $\partial f(\overline{x})$ is bounded. This latter condition is only satisfied for $f_{4}$, but not for $f_{5}$ and $f_{6}$. Nevertheless, we successfully minimize numerically the univariate function $J_{\rm approx}$ above by using the BiSection [BF85] method. Note that for computing each $\operatorname{dist}^{2}(X_{i},\tau\partial f(\overline{x}))$ one is faced with solving a quadratic program. Note that $\partial f(\overline{x})$ is described by linear equalities and inequalities, compare (2.7a) and (2.7c). The empirical probability of uniqueness shown as gray value plots in Figure 4.4, middle row, is plotted along with the curves describing $\min_{\tau\geq 0}J_{\rm approx}(\tau)$. We emphasize the perfect agreement of the phase transition.

4.2. 2D Empirical Phase Transitions and Theoretical Bounds

4.2.1. Testset

Creating gradient sparse images with a given gradient sparsity is more involved than in the one dimensional case. Using random support sets $\Lambda\subseteq[p]$ with $|\Lambda|\geq n$ and the projecting technique used in the 1D case, one would most likely obtain constant 2D images. Thus we use a different approach to construct random gradient sparse images. To this end, we randomly add binary images with homogeneous areas and use the modulo operation to binarize again the result. We show some results in Figure 4.3. Since it is easy to identify all connected components in a binary image, we can assign random real values to the different connected components. The resulting discrete gradient of the new image has the same number of non-zeros, but the image is non binary.

4.2.2. Tomographic Measuremnts

Making our experiments more general, we applied the idea from [RLH14] to reduce the influence of angles in tomographic reconstructions. The basic idea is, that each projection should carry the same amount of information independent of the angle. In Figure 4.6 we illustrate the difference the effect on the tomographic projections when sensing an images embedded in a rectangle shape resp. to a circular shape. We emphasize that our testset for 2D contains only images embedded in a circular shape.

4.2.3. $\ell_{1}$-Minimization

The phase diagrams for 2D images coincide with the ones for 1D signals, shown in Figure 4.4, first row, and are omitted here. The two regions are accurately separated by the $\min_{\tau}J(\tau)$ curve.

4.2.4. $\operatorname{TV}$-Minimization

For fixed ambient dimension $n=64\cdot 64$ we choose a relative image sparsity $\frac{\|D\overline{x}\|_{0}}{p}\in[0,1]$, with $D\in\mathbb{R}^{p\times n}$, corresponding to a gradient sparse test image and also choose the number of measurements $m$. Measurements can be Gaussian or tomographic. We consider three types of tomographic matrices:

• •

binary, as described in [RLH14];

• •

perturbed binary, that is we slightly perturb the nonzero entries above, in order to remove linear dependencies between columns, and

• •

standard tomographic matrices, with nonnegative real entries. We use the MATLAB routine paralleltomo.m from the AIR Tools package [HSH12] that implements such a tomographic matrix for a arbitrary vector of angles. We choose equidistant angles, set N $=64$ the image size and use the default value of p, i.e. the number of parallel rays for each angle p = round(sqrt(2)*N) to obtain a tomographic matrix of size $m\times n$.

Results are presented in Figure 4.4, last row, and Figure 4.5. We show the empirical probability over the 10 repetitions. The success rate of image reconstruction equals the uniqueness result and is displayed by gray values: black $\leftrightarrow 0\%$ uniqueness rate, white $\leftrightarrow 100\%$ uniqueness rate.

All plots display a phase transition and thus exhibit regions where exact image reconstruction has probability equal or close to one. All regions are accurately separated by our approximation to the statistical dimension.

5. Conclusions

The present work was motivated by the need of accurately describing how much undersampling is tolerable to uniquely recover a sparse or gradient sparse box constrained signal from tomographic measurements. Our results show sharp average-case phase transitions from non-uniqueness and non-recovery to uniqueness and exact recovery as guaranteed by compressed sensing for random measurements. Moreover, we show that the phase transition occurs approximately at the statistical dimension of the descent cone of the structure enforcing objective, irrespectively of the employed measuring device. The approximation of the statistical dimension is computed empirically but cross-checked by the verifiable uniqueness test developed in this work.

Appendix

We collect few results from [Roc70, RW09]. The subdifferential of a function $f\colon\mathbb{R}^{n}\to\overline{\mathbb{R}}$ at a point $\overline{x}$ where $f(\overline{x})\in\mathbb{R}$, is defined by

[TABLE]

We set $\partial f(\overline{x})=\emptyset$ if $f(\overline{x})$ is not finite. A function $f\colon\mathbb{R}^{n}\to\overline{\mathbb{R}}$ is said to be subdifferentiable at $\overline{x}$ if $\partial f(\overline{x})\not=\emptyset$. Then the elements of $\partial f(\overline{x})$ are called subgradients of $f$ at $\overline{x}$. For a proper, convex function $f$ and $\overline{x}\in\operatorname{rint}(\operatorname{dom}f)$, the subdifferential $\partial f(\overline{x})$ is nonempty. Furthermore, $\partial f(\overline{x})$ is nonempty and bounded if and only if $\overline{x}\in\operatorname{int}(\operatorname{dom}f)$. Hence only for $f_{1}$ and $f_{4}$ from (1.3a) and (1.3f) $\partial f(\overline{x})$ is compact. For the remaining choices of $f$ in (1.3) $\partial f(\overline{x})$ is unbounded, since in these cases our signals $\overline{x}$ lie on the boundary of $\operatorname{dom}f=\mathbb{R}^{n}_{+}$ or of $\operatorname{dom}f=\{0,1\}^{n}$.

The relative interior of a non-empty convex set $C$ is the interior relative to its affine hull,

[TABLE]

The relative boundary of $C$ is

[TABLE]

We always have [Roc70, Cor. 6.6.2]

[TABLE]

Recall that the affine hull of some set $C\subset\mathbb{R}^{n}$ is the set of all affine combinations of its points,

[TABLE]

while the conic hull equals

[TABLE]

Acknowledgment

JK gratefully acknowledges the support by the German Science Foundation, grant GRK 1653.