`1 -Analysis Minimization and Generalized (Co-)Sparsity:
When Does Recovery Succeed?
Martin Genzel, Gitta Kutyniok, and Maximilian März

arXiv:1710.04952v1 [] 13 Oct 2017

Technische Universität Berlin, Department of Mathematics
Straße des 17. Juni 136, 10623 Berlin, Germany
E-Mail: [genzel,kutyniok,maerz]@math.tu-berlin.de
Abstract. This paper investigates the problem of signal estimation from undersampled noisy sub-Gaussian measurements under the assumption of a cosparse model. Based on generalized notions of sparsity, we derive novel
recovery guarantees for the `1 -analysis basis pursuit, enabling highly accurate predictions of its sample complexity. The corresponding bounds on the number of required measurements do explicitly depend on the Gram matrix
of the analysis operator and therefore particularly account for its mutual coherence structure. Our findings defy
conventional wisdom which promotes the sparsity of analysis coefficients as the crucial quantity to study. In fact,
this common paradigm breaks down completely in many situations of practical interest, for instance, when applying a redundant (multilevel) frame as analysis prior. By extensive numerical experiments, we demonstrate that, in
contrast, our theoretical sampling-rate bounds reliably capture the recovery capability of various examples, such as
redundant Haar wavelets systems, total variation, or random frames. The proofs of our main results build upon
recent achievements in the convex geometry of data mining problems. More precisely, we establish a sophisticated
upper bound on the conic Gaussian mean width that is associated with the underlying `1 -analysis polytope. Due to a
novel localization argument, it turns out that the presented framework naturally extends to stable recovery, allowing
us to incorporate compressible coefficient sequences as well.
Key words. Compressed sensing, cosparse modeling, analysis sparsity, `1 -analysis basis pursuit, stable recovery,
redundant frames, total variation, Gaussian mean width

1 Introduction
Initiated by the pioneering works of Candès, Donoho, Romberg, and Tao [CRT06; CT05; Don06], a considerable amount of research on compressed sensing during the last decade has dramatically changed our
methodology for exploiting structure in many signal processing tasks. The classical setup in this field
considers the problem of reconstructing an unknown sparse signal vector x∗ ∈ Rn from non-adaptive,
linear measurements of the form
y = Ax∗ + e,
(1.1)
where A ∈ Rm×n is a known measurement matrix and e ∈ Rm captures potential distortions, usually
due to noise. The success of compressed sensing is based on the fundamental insight that, by explicitly
incorporating a sparsity prior, this task becomes feasible even if m  n. In particular, there exist
numerous convex and greedy recovery methods that enjoy both efficient implementations and a rich
theoretical foundation. Among them, probably the most popular approach is the basis pursuit [CDS98]:
min k xk1

x ∈Rn

subject to

k Ax − yk2 ≤ η,

(BPη )

where η ≥ 0 is such that kek2 ≤ η. The crucial ingredient of (BPη ) is the `1 objective functional
which promotes sparse solutions of the minimization program. Remarkably, it can be shown that
(BPη ) indeed recovers an S-sparse vector1 x∗ with the optimal sampling rate of m = O(S · log(2S/n))
if the measurement design A is drawn according to an appropriate random distribution.
Although such types of theoretical guarantees are elegant and practically appealing, the traditional
assumption of sparsity is not satisfied in most real-world applications. Fortunately, many signals-ofinterest do at least exhibit a low-complexity representation with respect to a certain transformation, for instance, the wavelet or Fourier transform. In this work, we study the the so-called analysis sparsity model
(also known as cosparse model), which has gained increasing attention within the past years [CENR11;
1 That

means, at most S entries of x∗ are non-zero, or more formally, k x∗ k0 := |supp( x∗ )| ≤ S.

1

2

1 I NTRODUCTION

KR15; NDEG13]. The key idea of this approach is to test (“analyze”) the signal x∗ ∈ Rn with a collection of analysis vectors ψ1 , . . . , ψN ∈ Rn , i.e., one computes
Ψx∗ = (hψ1 , x∗ i, . . . , hψN , x∗ i) ∈ R N ,

(1.2)

where the matrix Ψ := [ψ1 . . . ψN ]T ∈ R N ×n is called the analysis operator. If Ψ is able to reflect the
underlying structure of x∗ , one might expect that these analysis coefficients are dominated by only a few
large entries. This hypothesis of transform sparsity motivates the following generalization of (BPη )
that is typically referred to as the analysis basis pursuit (or `1 -analysis minimization):
min kΨxk1

x ∈Rn

subject to

k Ax − yk2 ≤ η.

(BPΨ
η)

In fact, the adapted strategy of (BPΨ
η ) has turned out to work surprisingly well for numerous problems
settings, such as in total variation minimization [Cha04; ROF92] or for multiscale transforms in classical
signal and image processing tasks [COS10; DPSS14; LDSP08; PM11]. Apart from that, it has become
also relevant to regularizing physics-driven inverse problems [KABG16; NG12] and operator learning
approaches based on test data sets [HKD13; RB11; RB13]. But despite this empirical success, many
theoretical properties of the analysis basis pursuit remain unexplored and “its rigorous understanding
is still in its infancy.” [KR15, p. 174]

1.1 Does (Co-)Sparsity Explain the Success (or Failure) of the Analysis Basis
Pursuit?
As already foreshadowed by (1.2), a central objective of the analysis formulation is to come up with an
operator Ψ ∈ R N ×n that provides a coefficient vector Ψx∗ of “low complexity.” The traditional theory
of compressed sensing—where Ψ is just the identity—would suggest that the sparsity of Ψx∗ is the
key concept to look at. Indeed, a large part of the related literature precisely builds upon this intuition.
Although many of those approaches rely on different proof strategies, e.g., the D-RIP [CENR11] or an
adapted null space property [KR15], they eventually promote results of a very similar type: Recovery
of x∗ ∈ Rn via (BPΨ
η ) succeeds if the number of measurements obeys
m ≥ C · S · PolyLog( 2N
S ),

(1.3)

where S := kΨx∗ k0 and C > 0 is a constant that might depend on Ψ. This bound on the sampling rate
clearly resembles classical guarantees for the basis pursuit (BPη ) and therefore forms a quite natural
extension towards analysis sparsity.
Another important branch of research takes a somewhat contrary perspective, and identifies the
cosparsity L := N − S as a crucial quantity for the success of the analysis model. A remarkable observation of [NDEG13] was that the location of vanishing coefficients in Ψx∗ is the driving force behind the
analysis formulation, rather than the number of non-zero components. This viewpoint naturally leads
to the so-called cosparse signal model, which is typically described by union-of-subspaces [BD09]. Following this terminology, it has turned out that successful recovery via combinatorial searching can be
guaranteed if the number of observations is of the order of the signal’s manifold dimension [NDEG13,
Sec. 3]. However, such a simple relationship does not seem to carry over to tractable methods like
(BPΨ
η ), see [GPV15]. Indeed, many sampling-rate bounds relying on cosparsity can be easily translated
into sparsity-based statements, simply meaning that S is replaced by N − L in (1.3), e.g., see [KR15,
Thm. 1] or [GNEGD14, Thm. 3.8].
Even tough the above approaches are quite appealing due to their simplicity, interpretability, and
consistency, it still remains unclear whether they are sufficient for a sound foundation of (BPΨ
η ) in its
general form. The notions of analysis (co-)sparsity are merely determined by the support of Ψx∗ , which
in turn does not account for the coherence structure of the individual analysis vectors ψ1 , . . . , ψN ; in
other words, the underlying “geometry” of Ψ is ignored.

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

3

To get a first glimpse of this issue, let us consider a simple example: Figure 1 shows the results of
a numerical simulation that reconstructs a block-signal (see Figure 1(a)) using three different analysis
operators. The plot of Figure 1(b) exhibits the phase transition behavior of (BPΨ
η =0 ) for a redundant,
discrete Haar wavelet transform Ψrdwt and the analysis operator Ψirdwt associated with its inverse transform (i.e., a certain dual frame of Ψrdwt ). Although the (co-)sparsity (S = 906 and L = 886) are exactly
the same for both choices, their recovery capability indeed differs dramatically! In conclusion, just investigating the parameters S and L does by far not explain why the transition of Ψirdwt happens much
earlier (m ≈ 85) than the one of Ψrdwt (m ≈ 240). Even more striking, the prediction of (1.3) deviates
from the truth by orders of magnitudes.
The plot of Figure 1(c) reveals another insight: While (1.3) is reliable for an orthonormal Haar wavelet
irdwt
transform Ψdwt ∈ Rn×n , a comparison with the actual recovery rate of (BPΨ
η =0 ) indicates that redundancy can be beneficial in the analysis model. Finally, it is also worth mentioning that a compressibility
argument does “save the day” here: Figure 1(d) demonstrates that the transition takes place far before
the remaining coefficients would be negligibly small.
We emphasize that the above example is not too specific or artificial, but rather illustrates a scenario
that often occurs in applications: Due to linear dependencies within Ψ, the analysis sparsity oftentimes
cannot become arbitrarily small. For instance, if Ψ corresponds to a highly redundant dictionary, one
typically has S  n, whereas the true sample complexity is relatively small (below the space dimension
n). This important phenomenon gives rise to the following fundamental question:
If (co-)sparsity does not entirely explain what is going on, which general principles make the analysis
basis pursuit (BPΨ
η ) succeed or fail?
To be a more precise about this concern, let us formulate three central issues that we will address in
this paper:
(Q1) Sampling rate. How many measurements are needed for an accurate estimate of x∗ via (BPΨ
η )?
What crucial quantities and parameters determine the sample complexity?
(Q2) Compressibility. What if x∗ is only compressible, i.e., its analysis coefficients Ψx∗ are not exactly
sparse but only close to a sparse vector? Is (BPΨ
η ) stable under such model inaccuracies?
(Q3) Interpretability. What practical guidelines can we derive from our recovery results? Which characteristics of Ψ deserve special attention in applications?

1.2 Main Contributions and Overview
A major concern of this work is to shed more light on the problem of `1 -analysis recovery and to
provide a deeper understanding of its underlying mechanisms. Our first main result in Section 2
(Theorem 2.5) focuses on the setup of noiseless Gaussian measurements, i.e., e = 0 in (1.1) and A is a
standard Gaussian matrix. In this prototypical situation, we can even expect an exact reconstruction of
x∗ via (BPΨ
η =0 ), supposed that m is sufficiently large. In fact, Theorem 2.5 states a non-asymptotic and
non-uniform bound on the sampling rate which intends to meet the following desiderata:
(i) Accurate. The number of required measurements should be close to the (optimal) sample complexity.
(ii) Computable. The expressions should be explicit with respect to the the underlying analysis operator Ψ and numerically evaluable.
(iii) Interpretable. The involved parameters should have a clear meaning and resemble well-known
principles.
(iv) Generic. The bound should not be tailored to a specific operator Ψ (e.g., wavelets or total variation), but apply in a general setting.

4

1 I NTRODUCTION

(a)

(b)

(c)

(d)

Figure 1: The phase transition of (BPΨ
η =0 ) for wavelets in 1D with noiseless Gaussian measurements (η = 0).
For a detailed description of this experiment, see Subsection 3.1. (a) Signal vector x∗ ∈ R256 (blocks, see
Ψirdwt
rdwt
[DJ94]). (b) Success rate of exact recovery via (BPΨ
η =0 ) and (BPη =0 ), respectively. (c) Success rate of exact
∗
1
∗
dwt
recovery via (BPΨ
η =0 ). (d) Sorted and ` -normalized magnitudes of the coefficients Ψrdwt x and Ψirdwt x .
The decay behavior is very similar in both cases.

The key ingredients of our bound are three parameters, generalizing the notions of sparsity and cosparsity (Definition 2.3). While the support of the coefficients Ψx∗ still plays a big role, these novel quantities do also take account of the (mutual) coherence structure of Ψ by incorporating the (off-diagonal)
entries of its Gram matrix. Indeed, this refinement will allow us to satisfy the above “wish list” for
many different examples of interest. Besides the well-known case of orthonormal bases, this particularly includes highly redundant analysis operators, which are only poorly understood in that context so
far.
The proof of Theorem 2.5 (provided in Section 6) loosely follows a recent technique proposed by
[ALMT14, Recipe 4.1]. More precisely, we establish a highly non-trivial upper bound on the Gaussian
mean width1 of the descent cone of x 7→ kΨxk1 at x∗ (Theorem 6.8), which is known to characterize
the phase transition behavior of (BPΨ
η =0 ) for Gaussian measurements. Noteworthy, Kabanava, Rauhut,
and Zhang follow a very similar proof strategy in [KRZ15]. It is therefore not surprising that their
1 This

is equivalent to bounding the so-called statistical dimension, introduced in [ALMT14]. But we will keep using the notion of
‘Gaussian mean width’, which is more common in the literature.

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

5

approach bears a certain resemblance to ours, but there are in fact crucial differences as a more detailed
comparison in Section 3 and Subsection 4.4 will show.
Our second guarantee (Theorem 2.8) tackles the problem of stability (Q2) in the setting of noisy subGaussian measurements. It is based on the simple geometric idea of approximating the ground truth
x∗ by another signal vector x̄∗ ∈ Rn of lower complexity, in the sense that its associated Gaussian mean
width is significantly smaller. This novel result is achieved by means of a localized variant of the (conic)
Gaussian mean width [MPT07] which allows us to quantify the mismatch caused by approximating
x∗ . Let us emphasize that this approach goes beyond the “naive” reasoning for compressible vectors
that just suggests to determine the best S-term approximation of the analysis coefficient vector Ψx∗ .
Section 3 presents several applications of the framework developed in Section 2. A special focus
here is on the little studied example of redundant wavelet frames, but the variety of our results will be
also demonstrated for other important choices of Ψ, such as the total variation operator and random
frames. In this course, we assess the predictive quality of our theoretical findings by means of various
numerical experiments. Section 4 is then dedicated to an overview of the related literature, revisiting
the initial concern of Subsection 1.1: Does (co-)sparsity explain the recovery performance of (BPΨ
η )? This discussion is concluded by Section 5, where we address the issues of (Q3) in greater detail: Even though
the main objective of this work is a mathematical foundation, our novel perspective on (co-)sparsity
might have practical implications as well, for instance, when designing or learning an analysis operator for a specific task. Finally, we point out some unresolved challenges left for future works, in
particular, the question of how close to optimal our prediction bounds are.

1.3 General Notation
Let us fix some notations and conventions that will be frequently used in the remainder of this paper:
• For an integer d ∈ N, we set [d] := {1, . . . , d}. If I ⊂ [d], the set complement of I in [d] is
given by I c := [d] \ I . Vectors and matrices are denoted by lower- and uppercase boldface
letters, respectively. Unless states otherwise, their entries are indicated by subscript indices and
0
lowercase letters, e.g., x = ( x1 , . . . , xd ) ∈ Rd for a vector and B = [bk,k0 ] ∈ Rd×d for a matrix.
For I ⊂ [d], the vector xI ∈ R|I| is the restriction of x ∈ Rd to the components of I . Similarly,
0
0
BI ∈ R|I|×d restricts a matrix B ∈ Rd×d to the rows corresponding to I .
• Let x = ( x1 , . . . , xd ) ∈ Rd . The support of x is defined by the set of its non-zero entries supp( x) :=
{k ∈ [d] | xk 6= 0} and the sparsity of x is k xk0 := |supp( x)|. For 1 ≤ p ≤ ∞, we denote the
` p -norm on Rd by k · k p . The associated unit ball is given by Bdp := { x ∈ Rd | k xk p ≤ 1} and the
Euclidean unit sphere is Sd−1 := { x ∈ Rd | k xk2 = 1}. If ΣdS := { x ∈ Rd | k xk0 ≤ S} denotes the
set of all S-sparse vectors in Rd , the best S-term approximation error (with respect to the ` p -norm) of
a vector x0 ∈ Rd is
σS ( x0 ) p := min k x0 − xk p .
x∈ΣdS

• We denote the conic hull of a set K ⊂ Rd by cone(K ). If U ⊂ Rd is a linear subspace, the associated
orthogonal projection onto U is denoted by PU ∈ Rd×d . Then, we have PU ⊥ = Id − PU , where
U ⊥ ⊂ Rd is the orthogonal complement of U and Id ∈ Rd×d is the identity matrix.

• If g is a mean-zero Gaussian random vector in Rd with covariance matrix Σ ∈ Rd×d , we just write
g ∼ N (0, Σ). Moreover, a random vector a in Rd is sub-Gaussian if k akψ2 < ∞, where k · kψ2 is
the sub-Gaussian norm; see [Ver12, Def. 5.22].
• Let Ψ ∈ R N ×n be a matrix with row vectors ψ1 , . . . , ψN ∈ Rn . Then, the collection F :=
{ψk }k∈[ N ] ⊂ Rn forms a frame for Rn with frame bounds 0 < a ≤ b < ∞ if
a · k xk22 ≤ kΨxk22 ≤ b · k xk22

for all x ∈ Rn .

6

2 M AIN R ESULTS

If there is no danger of confusion, we will identify the analysis operator Ψ with F . Finally, F is
said to be tight if a = b. See [CK13] for a detailed introduction to finite-dimensional frame theory.
• For s, t ∈ R, we define the clip function clip(s; t) := sign(s) · min{|s|, t} and the positive part
[s]+ := max{s, 0}.
• The letter C is always reserved for a (generic) constant, whose value could change from time to
time. We refer to C as a numerical constant if its value does not depend on any other involved
parameter. If an (in-)equality holds true up to a numerical constant C, we sometimes simply
write A . B instead of A ≤ C · B.

2 Main Results
In this section, we present the main theoretical results of this work. Starting with some notations and
model assumptions in Subsection 2.1, our novel sparsity parameters are introduced in Subsection 2.2.
Based on these central notions, Subsection 2.3 then establishes an exact recovery result for noiseless
Gaussian observations. Finally, in Subsection 2.4, we show that the analysis basis pursuit (BPΨ
η ) is also
robust against noise and stable under model inaccuracies, meaning that the analysis coefficient vector
is allowed to be compressible. For the sake of readability, all corresponding proofs are postponed to
Section 6.

2.1 Model Setup and Notation
In this part, we state the standing assumptions and model hypothesis of our recovery framework; see
also Table 1 for a summary. The following notations and conventions are supposed to hold true for the
remainder of this section. Let us first recall the linear measurement scheme of (1.1) and set up a formal
random observation model:
Model 2.1 (Noisy Linear Measurements) Let x∗ ∈ Rn be a fixed vector, which is typically referred to
as the signal (or source). The measurement vectors a1 , . . . , am ∈ Rn are assumed to be independent copies
of an isotropic, mean-zero, sub-Gaussian random vector a in Rn with k akψ2 ≤ κ. These vectors form
the rows of the measurement matrix A := [ a1 . . . am ]T ∈ Rm×n . The actual measurements of x∗ are then
given by
y := Ax∗ + e ∈ Rm ,
where e = (e1 , . . . , em ) ∈ Rm models noise, which could be deterministic and systematic. We assume
that the noise is `2 -bounded, i.e., kek2 ≤ η for some η ≥ 0.
Next, we fix our notation for the analysis operator and coefficients. Note that the dimension of the
analysis domain R N does not necessarily have to be larger than the dimension of Rn .
Definition & Notation 2.2 (1) The matrix Ψ = [ψk,j ] ∈ R N ×n is called an analysis operator (or analysis
matrix) if none of its rows equals the zero vector. The rows of Ψ, denoted by ψ1 , . . . , ψN ∈ Rn , are
called the analysis vectors. Moreover, we define the Gram matrix of Ψ as
G = [ gk,k0 ] := ΨΨT = [hψk , ψk0 i] ∈ R N × N .

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

Term
(Ground truth) Signal vector
Measurement vectors
Measurement matrix
Noise variables
Measurement variables
Analysis vectors
Analysis operator
Gram matrix
Analysis coefficients of x ∈ Rn
Analysis support and sparsity
Analysis cosupport and cosparsity
Analysis sign vector

7

Notation
x ∗ ∈ Rn
a 1 , . . . , a m ∈ Rm
A = [ ai,j ] = [ a1 . . . am ]T ∈ Rm×n
e = ( e1 , . . . , e m ) ∈ Rm
y = (y1 , . . . , ym ) = Ax∗ + e ∈ Rm
ψ1 , . . . , ψN ∈ Rn
Ψ = [ψk,j ] = [ψ1 . . . ψN ]T ∈ R N ×n
G = [ gk,k0 ] = ΨΨT = [hψk , ψk0 i] ∈ R N × N
Ψx = (hψ1 , xi, . . . , hψN , xi) ∈ R N
S x = supp(Ψx) ⊂ [ N ] and S = |S x |
S xc = supp(Ψx)c ⊂ [ N ] and L = |S xc | = N − S
σx = (σx1 , . . . , σxN ) = sign(Ψx) ∈ {−1, 0, 1} N

Table 1: A summary of important notations used in this paper.

(2) The analysis coefficients of a vector x ∈ Rn (with respect to Ψ) are given by
Ψx = (hψ1 , xi, . . . , hψN , xi) ∈ R N .
The analysis support of x is denoted by S x := supp(Ψx) ⊂ [ N ], and if S = |S x |, we say that x is Sanalysis-sparse.1 Analogously, we call the complement S xc := supp(Ψx)c ⊂ [ N ] the analysis cosupport of
x and speak of an L-analysis-cosparse vector if L = |S xc | = N − S.

2.2 Generalized (Co-)Sparsity
Before presenting the actual recovery results, we need to introduce three adapted notions of sparsity,
which will form the heart of our sampling-rate bounds:
Definition 2.3 (Generalized (Co-)Sparsity) Let Ψ ∈ R N ×n be an analysis operator and let x ∈ Rn . We
define the generalized sparsity of x (with respect to Ψ) by
SSσxx :=
where

0

∑
0

σxk · σxk · gk,k0 ,

k,k ∈S x

σxk := sign(hψk , xi) ∈ {−1, 0, 1},

k = 1, . . . , N,

and σx := (σx1 , . . . , σxN ) ∈ {−1, 0, 1} N is called the analysis sign vector of x. Moreover, we introduce the
terms
LSxc :=
Lr
S c :=
x

∑
0

√

k,k ∈S xc

∑c

√

2
gk,k
0

gk,k · gk0 ,k0

gk,k ,

,
(2.1)

k ∈S x

which are both referred to as the generalized cosparsity of x.2
Note that, for the sake of readability, we have omitted the dependence on Ψ here. Considering
the canonical choice of an orthonormal basis, it becomes actually clear why we speak of ‘generalized’
1 If

there is no danger of confusion, we may omit the term ‘analysis’ and just speak of ‘coefficients’, ‘support’, ‘sparsity’, etc.
label ‘r ’ indicates that Lr
S c only operates on the diagonal entries of G.

2 The

x

8

2 M AIN R ESULTS

Figure 2: The graph of Φ. This function is strictly monotonically decreasing (with limρ&0 Φ0 (ρ) = −∞), and
at its boundary points, we have limρ&0 Φ(ρ) = 1 and limρ→∞ Φ(ρ) = 0.

sparsity: Since G = ΨΨT = In in this case, one has that SSσxx = S = kΨxk0 and LSxc = Lr
S xc = L = n − S;
hence, the notions of Definition 2.3 precisely coincide with their traditional counterparts.
In general, this correspondence is more complicated. The definition of the generalized sparsity SSσxx
still operates on the analysis support of x, but also involves a weighted sum over all Gram matrix
entries associated with S x . The same holds true for the generalized cosparsity term LSxc , respectively.
Such an incorporation of the off-diagonal entries of G is in fact quite appealing because, to a certain
extent, it respects the mutual coherence structure of the analysis vectors ψ1 , . . . , ψN . Finally, let us state
two basic observations, characterizing when the above sparsity parameters do not vanish (see also
Lemma 6.11):
SSσxx = kΨT σx k22 > 0
LSxc , Lr
S xc

if, and only if,

x 6∈ ker Ψ,

> 0 if, and only if, S xc 6= ∅.

(2.2)

2.3 Sampling-Rate Function and Exact Recovery
To highlight the key features of our approach, we restrict our analysis to the simplified case of noiseless
Gaussian observations in this part, i.e., η = 0 and a ∼ N (0, In ) in Model 2.1. The analysis basis pursuit
then takes the form
minn kΨxk1 subject to Ax = y,
(BPΨ
η =0 )
x ∈R

and we can even hope for an exact retrieval of x∗ ∈ Rn . Indeed, the recent work of [ALMT14] has
made the remarkable observation that a convex program of this type typically undergoes a sharp phase
transition as m varies: Recovery of x∗ fails with overwhelmingly high probability if m is below a certain
threshold, but once m exceeds a small transition region, recovery succeeds with overwhelmingly high
probability; see Figure 1(b) for an example. This minimal number of required measurements (also
depending on the desired probability of success) is often referred to as the sample complexity (or optimal
sampling rate) of an estimation problem. However, computing such a quantity in an explicit way is
usually a challenging task.
Our actual goal is therefore rather to come up with an upper bound on the sample complexity of
(BPΨ
η =0 ) that is tight in many situations. For this purpose, let us introduce the following function,
which essentially determines the sampling rate proposed by Theorem 2.5 below:

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

9

Definition 2.4 Let Ψ ∈ R N ×n be an analysis operator and let x∗ ∈ Rn with x∗ 6∈ ker Ψ. Then, we
define the sampling-rate function of Ψ and x∗ by

σ ∗ !
2

( Lr
SS x ∗

S xc∗ )

x

·Φ
, if S xc∗ 6= ∅,
n−
LS c∗
LS c∗
M(Ψ, x∗ ) :=
x
x



n,
otherwise,
where
Φ(ρ) := erf



h −1 ( ρ )
√
2

with
h : (0, ∞) → (0, ∞), τ 7→

q



,

ρ > 0,

2
2 e−τ /2
π
τ

+ erf( √τ ) − 1.
2

(2.3)

It is not hard to verify that the univariate functions h, h−1 , and Φ —each one mapping from (0, ∞)
to (0, ∞) —are well-defined; see first paragraph of Appendix A.3 for more details. Together with (2.2),
this also implies the well-definedness of M(Ψ, x∗ ). Figure 2 shows the graph of Φ, visualizing how
σ ∗
the ratio SSx∗ /LS c∗ affects the sampling-rate function.
x
x
With this notion at hand, we are now ready to state our first main result:
Theorem 2.5 (Exact Recovery via (BPΨ
η =0 )) Assume that Model 2.1 is satisfied with a ∼ N (0, In ) and η = 0.
N
×
n
Let Ψ ∈ R
be an analysis operator such that x∗ 6∈ ker Ψ. Then, for every u > 0, the following holds true
2
with probability at least 1 − e−u /2 : If the number of measurements obeys
m>

q

M(Ψ, x∗ ) + u

2

+ 1,

(2.4)

∗
then (BPΨ
η =0 ) recovers x exactly.

Theorem 2.5 gives a first answer to (Q1): Roughly speaking, reconstruction succeeds with high
probability as m slightly exceeds the sampling-rate function M(Ψ, x∗ ). Since M(Ψ, x∗ ) is completely
determined by our generalized sparsity parameters from Definition 2.3, we can even conclude that the
bound (2.4) meets the desiderata of (ii), (iii), and (iv) requested in Subsection 1.2.
Unlike many approaches from the literature (cf. (1.3)), the statement of Theorem 2.5 is highly nonuniform. Indeed, the condition of (2.4) does not just involve the analysis sparsity S = kΨx∗ k0 , but
explicitly depends on the support S x∗ as well as on the sign vector σx∗ = sign(Ψx∗ ). We will return to
this point in the course of our experiments in Section 3, demonstrating that the performance of (BPΨ
η =0 )
is oftentimes not entirely explainable by means of S, even if Ψ is fixed.
On the other hand, our wish for an almost tight bound on the sample complexity has not been
clarified yet (see (i) in Subsection 1.2). This issue unfortunately turns out to be very challenging, and in
general, we do not have a quantitative error estimate for our prediction. However, we will give at least
numerical evidence in Section 3 and discuss some theoretical aspects of optimality in Subsection 6.3
(see Remark 6.15).
Remark 2.6 (1) The assumption of x∗ 6∈ ker Ψ in Theorem 2.5 is not very restrictive and merely of
∗
technical nature. If x∗ ∈ ker Ψ, the analysis basis pursuit (BPΨ
η =0 ) uniquely recovers x if, and only if,
ker Ψ ∩ ( x∗ + ker A) = { x∗ }.
In the setting of Gaussian measurements, where ker A can be identified almost surely with a random
subspace of dimension n − m, the latter condition holds true if, and only if, m ≥ dim(ker Ψ). Thus, the
dimension of ker Ψ precisely yields the sample complexity of (BPΨ
η =0 ) in this case.

10

2 M AIN R ESULTS

(2) An important feature of Theorem 2.5 is scaling invariance: Replacing Ψ by a scaled version λΨ with
λ 6= 0 does not affect the sampling-rate function, i.e., M(λΨ, x∗ ) = M(Ψ, x∗ ). This is due to the fact
σ ∗
that the parameters SSx∗ , LS c∗ , and Lr
S c do only appear in terms of appropriate fractions.
x∗

x

x

(3) Since M(Ψ, x∗ ) ≤ n, the condition of (2.4) does not lead to situations where m needs to be much
larger than n in order to achieve successful recovery. Such a requirement would be overly restrictive
because the equation system y = Ax∗ is almost surely uniquely solvable if m ≥ n, no matter what
operator Ψ is applied. In contrast, this simple observation is not always reflected by a naive bound of
the form (1.3), at least when the domain of analysis coefficients R N is much higher dimensional, i.e.,
N  n.
♦
The function Φ from Definition 2.4 is quite easy to understand from an analytical perspective, but it
is however non-elementary. For that reason, let us conclude with a simplification of the bound (2.4) in
Theorem 2.5. This obviously comes along with a certain loss of accuracy but sheds more light on the
asymptotic behavior of our sampling-rate function M(Ψ, x∗ ).
Proposition 2.7 Let Ψ ∈ R N ×n be an analysis operator and let x∗ ∈ Rn with x∗ 6∈ ker Ψ and S xc∗ 6= ∅.
σ ∗
Setting S := SSx∗ , L := LS c∗ , and L̄ := Lr
S c , it holds
x∗

x

x

n

M(Ψ, x∗ ) ≤ min n −

L̄2
L

+ ( LL̄ )2 · [2S · log( S+S L ) + S], n −

2
π

·

L̄2
S+ L

o

.

(2.5)

In particular, if Ψ ∈ Rn×n forms an orthonormal basis, we have S = kΨx∗ k0 and L = L̄ = n − S, so that

M(Ψ, x∗ ) ≤ 2S · log( Sn ) + 2S . S · log( 2n
S ).

(2.6)

The expression of (2.5) is indeed much more explicit than M(Ψ, x∗ ). The left branch particularly
reminds us of a well-known property of traditional bounds, where the sampling rate is supposed to
grow linearly with S up to a logarithmic factor. This intuition is also confirmed by (2.6), showing that
our approach is consistent with the standard setup of compressed sensing. On the contrary, the right
σ ∗
branch of (2.5) intends to mimics the asymptotics of M(Ψ, x∗ ) if the value of SSx∗ is large compared to
x

σ ∗

LS c∗ , i.e., SSx∗ /LS c∗  1.
x

x

x

2.4 Stable and Robust Recovery
In this part, we turn to our wish for stability as stated in (Q2). For illustration, let us recall the above
situation of exact recovery and inspect the sampling-rate function x∗ 7→ M(Ψ, x∗ ) more closely. While
this mapping actually insinuates a strong dependence on x∗ , a brief look at our generalized sparsity
parameters reveals that it is already determined by the sign vector σx∗ . In particular, the latter quantity
does not take account of the magnitudes of the analysis coefficients. The size of M(Ψ, x∗ ) may therefore change dramatically (in a discontinuous manner) if x∗ is slightly varied in such a way that the zero
entries of Ψx∗ become negligibly small but do not vanish anymore. For a full coefficient vector Ψx∗
with S = N, we would even have M(Ψ, x∗ ) = n, which typically occurs in real-world examples, such
as Figure 3. In these situations, the “binary” statement proposed by Theorem 2.5 is too pessimistic,
meaning that one cannot expect exact recovery with m  n. Instead of undergoing a sharp phase
transition, the reconstruction error of (BPΨ
η =0 ) will then decay rather smoothly as m grows.
Our approach to this problem relies on a simple geometric idea: Approximate x∗ ∈ Rn by a vector
∗
x̄ ∈ Rn whose analysis coefficients Ψ x̄∗ are much sparser than Ψx∗ . Such a “surrogate” of x∗ is
supposed to be of lower complexity, in the sense that M(Ψ, x̄∗ ) is small, and if x̄∗ and x∗ are not too
distant, we speak of analysis compressibility. In that case, one may hope that the actual performance
∗
∗
of (BPΨ
η ) is almost the same as if x̄ would be the ground truth signal, i.e., y = A x̄ . This desirable
property of a recovery method is usually referred to as stability.

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

11

A very natural, albeit naive, strategy to come up with a good low-complexity approximation of x∗ is
to solve
minn k x∗ − x̄k2 subject to |S x̄ | = kΨ x̄k0 ≤ S
(2.7)
x̄∈R

for some appropriately chosen sparsity threshold S. Any minimizer is then called a best S-analysis
sparse approximation of x∗ (with respect to Ψ). Since the constraint set of (2.7) is actually a union of
linear subspaces, it is also convenient to consider the following equivalent reformulation:
min k x∗ − PUS x∗ k2 ,

(2.8)

S⊂[ N ]
|S|≤S

where US := ker ΨS c ⊂ Rn . Interestingly, the same optimization program was recently studied in
[GNEGD14] under the name of optimal projections, forming a crucial step in greedy-like algorithms for
combinatorial cosparse modeling. Although (2.8) is tractable in some simple cases, it eventually turned
out to be NP-hard in general [TGP14], due to its combinatorial nature. However, such algorithmic
issues are only of minor importance to us, since we merely aim for a theoretical justification of stability
in `1 -analysis recovery.
For these reasons, we shall state our second main result for Euclidean projections onto an arbitrary
subspace U ⊂ Rn , leaving space for applying less challenging (near-optimal) approximation methods.
Note that, compared to Theorem 2.5, we now also allow for sub-Gaussian measurements as well as the
presence of noise.
N ×n
Theorem 2.8 (Stable and Robust Recovery via (BPΨ
η )) Suppose that Model 2.1 is satisfied and let Ψ ∈ R
∗
n
∗
be an analysis operator with x 6∈ ker Ψ. Let U ⊂ R be a linear subspace and assume that PU x 6∈ ker Ψ.
Then
kΨx∗ k1
· P x ∗ ∈ Rn
(2.9)
x̄∗ :=
kΨPU x∗ k1 U

is well-defined and there exists a numerical constant C > 0 such that, for every R > 0 and u > 0, the following
2
holds true with probability at least 1 − e−u /2 : If the number of measurements obeys1
m > m0 := C · κ 4 ·



R +1
R

·

q

2

M(Ψ, x̄∗ ) + 1 + u + 1,

(2.10)

then any minimizer x̂ of (BPΨ
η ) satisfies
n
2η
√
k x̂ − x∗ k2 ≤ max Rk x∗ − x̄∗ k2 , √m−1−
m

o
0 −1

.

(2.11)

In the Gaussian case of a ∼ N (0, In ), we particularly have C = κ = 1.
The statement of Theorem 2.8 makes the above geometric principle precise: Instead of evaluating the
sampling-rate function M(Ψ, ·) at the ground truth x∗ , we rather consider M(Ψ, x̄∗ ) in (2.10), which
could be significantly smaller. The price to pay is an additional error term in (2.11), governed by the
distance between x∗ and x̄∗ . Thus, one can formulate the following rule-of-thumb:
Recovery accuracy can be exchanged with investigating less complex signals (taking fewer measurements), and vice versa.
The benignity of this trade-off clearly depends on the ansatz space U. We already pointed out above
that selecting U is indeed very delicate and specific about the analysis operator Ψ. Moreover, even
if one could identify an optimal support set S according to (2.8), it is still not clear whether setting
U = US = ker ΨS c leads to the smallest sampling rate M(Ψ, x̄∗ ) among all S-analysis-sparse approximations of x∗ . For these reasons, we defer a general discussion of this issue to future works. A simple
1 One

may also set R = ∞ here, with the convention that

R +1
R

:= 1 and R · 0 := 0.

12

2 M AIN R ESULTS

(a)

(b)

Figure 3: (a) Visualization of the ‘cameraman’ x∗ ∈ R256×256 . (b) Sorted `1 -normalized magnitudes (in
logarithmic scale) of the analysis coefficients Ψx∗ , where Ψ corresponds to a compactly supported shearlet
frame with four scales, provided by Shearlab [KLR16].

greedy strategy is however presented in Appendix B.1, producing satisfactory results for redundant
wavelets frames (see Subsection 3.1.4).
Compared to the setup of the previous subsection, the sub-Gaussianity of A comes along with an
extra factor of C · κ 4 for the sampling rate, whereas the noise level η affects the actual√error estimate.
√
Roughly speaking, (2.11) states that we can expect reliable outcomes as long as kek2 / m ≤ η/ m =
O(1), which is a well-known observation, e.g., see [Tro15, Cor. 3.5].
Finally, we would like to highlight that our approach is based on approximations in the signal domain Rn . Working in the space of analysis coefficients R N instead is actually more common in the
literature (see Subsection 4.2), but suffers from disregarding the coherence structure of Ψ, which is in
turn problematic in highly redundant scenarios. There is however a connection between both regimes,
as the following bound for the approximation error in (2.11) shows:

k x∗ − x̄∗ k2 ≤

kPU x∗ k2
kΨPU x∗ k1

· kΨPU ⊥ x∗ k1 + kPU ⊥ x∗ k2 ,

see Proposition 6.7 for a proof. In fact, if Ψ is an orthonormal basis and U corresponds to the coordinate
space of the S-largest entries of Ψx∗ in magnitude, it is not hard to see that this estimate is closely
related to traditional bounds from compressibility theory (cf. [FR13, Thm. 4.19]).
Remark 2.9 (1) The scaling factor in (2.9) ensures that kΨx∗ k1 = kΨ x̄∗ k1 . Hence, x̄∗ is not precisely equal to the optimal approximation of x∗ in U, unless kΨx∗ k1 = kΨPU x∗ k1 . This technicality
is required for the proof of Theorem 2.8 (see Corollary 6.6), but we believe that it is not necessary in
general.

(2) The parameter R in Theorem 2.8 allows us to fine-tune the trade-off described above: Enlarging
R causes a decrease in accuracy in (2.11) but improves the sampling rate of (2.10) at the same time. In
this sense, the fraction ( R + 1)/R can be also regarded as oversampling factor.
♦

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

13

3 Applications and Numerical Experiments
The purpose of this section is to thoroughly evaluate the predictive quality of the theorems presented
in Section 2. A series of experiments based on redundant Haar wavelet frames (Subsection 3.1), total
variation minimization (Subsection 3.2), and tight random frames in general position (Subsection 3.3)
is conducted, which covers a wide range of analysis operators in one and two spatial dimensions
(Subsection 3.4). By creating phase transition plots, we shall visualize the sampling-rate bound of
Theorem 2.5 and investigate in how far our wish for accuracy is fulfilled. Combined with a simple
greedy approximation strategy for compressible signals, this also allows for a verification of our result
on stable recovery (Theorem 2.8).
While our theoretical framework includes a more general class of sub-Gaussian random matrices,
we do only consider the benchmark case of noiseless Gaussian measurements here, as it is often done
in the compressed sensing literature. Let us point out again that the statement of Theorem 2.5 is nonasymptotic so that it can be directly compared with the true recovery performance of (BPΨ
η =0 ). In contrast, a visualization of asymptotic results, such as stated in [CENR11; NW13], is inappropriate for
this setup due to unspecified numerical constants. The only other non-asymptotic bound, that we are
aware of, was established by Kabanava, Rauhut, and Zhang in [KRZ15], and will therefore often serve
as an object of comparison; see also Subsection 4.4 for a more formal presentation of their findings.
The notation of this part again follows Subsection 2.1, in particular Model 2.1. As it is common in the
literature (cf. [ALMT14]), only the sampling-rate function x∗ 7→ M(Ψ, x∗ ) is reported when analyzing
the quality of Theorem 2.5, whereas the probability parameter u in (2.4) is neglected. Unless stated
otherwise, we are solving the convex program (BPΨ
η =0 ) using the Matlab software package cvx [GB08;
GB14] with the default settings in place and the precision set to best. An outcome x̂ ∈ Rn is considered
to be an ‘exact’ recovery of x∗ ∈ Rn if k x̂ − x∗ k2 < 10−5 . Indeed, this threshold has proven to produce
very stable solutions and seems to reflect the numerical accuracy of cvx.

3.1 (Highly) Redundant Wavelet Frames
A major difficulty in analysis modeling is to select an appropriate operator Ψ ∈ R N ×n that enables
good reconstructions for a whole class of signals. One of the most popular classes is formed by piecewise
constant functions (e.g., the blocks signal of [DJ94]), oftentimes regarded as a prototype system in the
literature. In this context, translation-invariant Haar wavelet frames have turned out to be very useful
for various signal- and image-processing tasks, most prominently denoising [CD95; DJ94]. Depending
on the number of scaling levels, such wavelet transforms are typically highly redundant, i.e., the total
number of wavelet coefficients N may be orders of magnitudes larger than the signal’s dimension n.
However, despite their success and popularity, there exists to the best of our knowledge no sound
numerical or theoretical analysis of this important use case in respect of `1 -analysis recovery.
We emphasize that the situation of high redundancy is not just a flaw of this specific setup, but is
ubiquitous in applied harmonic analysis. For instance, when applying Gabor frames in time-frequency
analysis, or (directional) representation systems like curvelets or shearlets in imaging, one often has to
face problems with N/n ≥ 50 [KLR16]. Similar observations were recently also made for approaches to
analysis operator learning from training data: For example, [HKD13] reported that a redundancy factor
of N/n = 128 yields the best results for many classical imaging tasks. Apart from being considerably
redundant, these systems usually share the following properties:
• Linear dependency. There are (necessarily) linear dependencies among the analysis vectors of the
frame.
• Translation-invariance. The analysis vectors are obtained by translating a family of generators.
• Intrinsic localization. The frame is well-localized in the sense of [FG05; Grö04]. Roughly speaking,
this means that the off-diagonal decay of the Gram matrix G = ΨΨT is fast.

14

3 A PPLICATIONS AND N UMERICAL E XPERIMENTS

3.1.1 Setup and Notation of Wavelet Frames
Inspired by these observations, we build the simulations of this subsection on a combination of piecewise constant signals and Haar wavelet systems, serving as benchmark example of (highly) redundant
frames. Our specific implementation of wavelet transforms relies on the Matlab software package
spot [BF13], which is in turn adapted from the Rice Wavelet Toolbox [B+17]. The redundant wavelet
transform rdwt : Rn → R N of this package is based on the so-called algorithme à trous. Here, the downsampling steps of the discrete wavelet transform are omitted, whereas instead the filter responses
are upsampled by padding with zeros, which reminds us of creating holes (‘trous’ in French) in the
non-zero coefficients. In the following, let Ψrdwt ∈ R N ×n denote the (matricized) analysis operator
associated with this frame. The redundant forward transform of [BF13] comes along with an inverse
wavelet transform irdwt : R N → Rn satisfying
irdwt ◦ rdwt = IdRn .

(3.1)

From the viewpoint of frame theory, the application of irdwt corresponds to the synthesis operation
with a dual frame of Ψrdwt . More precisely, if Ψirdwt ∈ R N ×n denotes the analysis operator of this dual
frame, the rule (3.1) translates into ΨT
irdwt · Ψrdwt = In . It is worth mentioning that Ψrdwt and Ψirdwt
are not merely academic concepts that are inaccessible in practice. Indeed, the analysis and synthesis
operations of both frames are computable via the algorithm à trous with a complexity of O(n log(n)).
Finally, if n is a power of two, the non-redundant discrete wavelet transform dwt : Rn → Rn is orthonormal. Its associated analysis operator is denoted by Ψdwt ∈ Rn×n , forming an orthonormal basis for
Rn . Unless stated otherwise, we will always choose a Haar wavelet filter and 6 decomposition levels
for Ψrdwt , Ψirdwt , and Ψdwt . For a more detailed introduction to (redundant) wavelet transforms, the
interested reader is referred to [Fow05; Mal09].
3.1.2 Recovery of blocks Signal
Our first experiment revisits the simulation of Figure 1 in Subsection 1.1, where the reconstruction of
the classical blocks signal [DJ94] was investigated for different choices of Haar wavelet frames. Technically, we fix the ambient dimension n = 256, initialize x∗ ∈ Rn as blocks, and run Experiment 3.1
below for Ψ ∈ {Ψrdwt , Ψirdwt , Ψdwt } and M = {1, . . . , n}. The plots of Figure 1(b) and Figure 1(c) show
the empirical success rates1 for Ψ ∈ {Ψrdwt , Ψirdwt } and Ψ = Ψdwt , respectively.
Experiment 3.1 (Recovery of a fixed signal)
Input: Fixed ambient dimension n, signal vector x∗ ∈ Rn , analysis operator Ψ ∈ R N ×n , range of
measurements M ⊂ {1, 2, . . . , n}.
Compute: Repeat the following procedure 50 times for every m ∈ M:

I Draw a standard i.i.d. Gaussian random matrix A ∈ Rm×n and determine the measurement
vector y = Ax∗ .
n
I Solve the analysis basis pursuit (BPΨ
η =0 ) to obtain an estimator x̂ ∈ R .

I Compute and store the recovery error k x∗ − x̂k2 .
Figure 1(b) reveals that a perfect reconstruction of x∗ with Ψrdwt requires almost m ≈ n measurements. This observation is remarkable, since using Ψrdwt for denoising via soft thresholding has
proven to be very effective [CD95], whereas it seems to be quite inappropriate for the `1 -analysis formulation. In contrast, the choice of Ψirdwt performs significantly better, succeeding from m ≈ 105

1 As

pointed out above, ‘success’ means that x∗ is recovered up to an accuracy of k x∗ − x̂k2 < 10−5 .

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

15

Figure 4: Predicted error decay of (3.3) proposed by [KR15, Thm. 2].

onward. Interestingly, the same superiority remains true for analysis basis pursuit denoising (i.e., applying (BPΨ
η ) with A = In ), which indicates that the analysis approach is fundamentally different from
classical denoising via soft thresholding.
The subsampled orthonormal basis Ψdwt behaves slightly worse than Ψirdwt and allows for retrieval
with m ≈ 125 onward (see Figure 1(c)). This gives further evidence that the redundancy of a sparsifying transform can help to improve the recovery capability of analysis-based priors.
A closer look at the structure of Ψrdwt and Ψirdwt reveals that they actually correspond to the same
frame up to a reweighting of their analysis vectors. In fact, it holds
ψrdwt,k =

ψirdwt,k
∈ Rn
kψirdwt,k k2

(3.2)

for k = 1, . . . , N, particularly implying that the respective (co-)supports S x∗ and S xc∗ coincide in both
cases. Following traditional approaches based on analysis sparsity [CENR11; KR15] or cosparse modeling [NDEG13], one would therefore conclude that the number of samples required for exact recovery
should be equal as well.
Since this is far from being true, one might be tempted to take the perspective of analysis compressibility, instead of insisting upon perfect sparsity. However, the sharp phase transition in Figure 1(b) as
well as the decay of the analysis coefficients in Figure 1(d) show that this principle does not provide a
valid explanation. Indeed, even the most accurate and explicit bound known from the literature [KR15,
Thm. 2] predicts an error rate of
2σS (Ψx∗ )1
k x∗ − x̂k2 ≤ √
,
(3.3)
a·S
supposed that, roughly speaking, m ≥ 10 · (b/a) · S · log( N/S). Here, a and b are lower and upper
frame bounds of Ψ, respectively, and σS (Ψx∗ )1 denotes the best S-term approximation error of Ψx∗ (cf.
Subsection 1.3). The visualization of (3.3) in Figure 4 indicates that a reliable reconstruction can be only
expected if m exceeds the ambient dimension n by a factor of 3, contradicting the actual performance
of (BPΨ
η =0 ). Moreover, the plots do even certify a superiority of Ψrdwt over Ψirdwt , whereas the complete opposite is true. We finally would like to point out that an argumentation based on the matrix
condition is also misleading in this setup because both matrices Ψrdwt and Ψirdwt actually have the
same condition number (which equals the square root of the frame bound ratio b/a); see Table 2.
Turning towards the predictive quality of Theorem 2.5, we observe in Figure 1(b) that our bound
on the sampling rate very precisely captures the location of the phase transitions for both Ψrdwt and
Ψirdwt . While M(Ψrdwt , x∗ ) almost perfectly hits the center of transition curve (50% success rate),

16

3 A PPLICATIONS AND N UMERICAL E XPERIMENTS
Ana. op. Ψ

Scales

a

b

b/a

N

S

L

SSσ

LS c

Lr
Sc

M(Ψ, x∗ )

δ(Ψ, x∗ )

Ψdwt
Ψrdwt
Ψirdwt
Ψirdwt

6
6
6
3

1
2
1/64
1/8

1
64
1/2
1/2

1
32
32
4

256
1792
1792
1024

41
906
906
306

215
886
886
718

41
28462
33
39

215
5180
197
194

215
886
212
209

114
241
100
109

114
240
84
95

Table 2: Characteristic parameters of the simulation in Figure 1. Here, a and b denote the lower and upper frame bound of Ψ, respectively. The (co-)sparsity of the analysis coefficients Ψx∗ are S = |S x∗ | and
σ∗
L = |S xc∗ | = N − S, whereas SSx∗ , LSxc∗ , and LSrc∗ correspond to the generalized sparsity terms from Defx

x

inition 2.3. The sampling-rate function M(Ψ, x∗ ) is defined according to Definition 2.4. The true sample
complexity is denoted by δ(Ψ, x∗ ), see Appendix B.2. Note that all values are rounded.

M(Ψirdwt , x∗ ) is slightly more pessimistic, located at a rate of 75% success. There seems to be a common pattern in all of our examples: The higher the sample complexity1 δ(Ψ, x∗ ), the more accurately it
is approximated M(Ψ, x∗ ). On the other hand, the prediction for Ψdwt in Figure 1(c) is perfect, which is
due to the fact that our bound is provably tight in the orthonormal case (cf. Proposition 2.7). Tabular 2
lists a number of key quantities for all previously discussed choices of Ψ. For the sake of completeness,
in the last line, we have also documented the same parameters for a version of Ψirdwt with 3 scales
instead of 6. In a nutshell, we can conclude our discussion as follows:
The capacity of the analysis basis pursuit using highly redundant frames is not solely captured by
(co-)sparsity. In contrast, the sampling-rate bound proposed by Theorem 2.5 reliably states when
recovery is possible and when not.

3.1.3 Asymptotic Behavior of the Sampling Rate in n
In the setup of the previous subsection, we have fixed both the ambient dimension n and the ground
truth x∗ ∈ Rn . Moreover, changing n does also require an adaptation of the analysis operator Ψ ∈
R N ×n , so that the value of M(Ψ, x∗ ) is affected in a highly non-trivial way. For these reasons, we
next investigate the accuracy of M(Ψ, x∗ ) if n is varied, while x∗ and Ψ are just adapted in terms of
different resolution levels. Since the theoretical guarantee of [KRZ15, Thm. 3.3] is of a similar type, it
is also reported in our simulation below.
In order to create the plots of Figure 5, we again consider a combination of x∗ ∈ Rn as the blocks signal [DJ94] and Ψirdwt ∈ R N ×n as analysis operator, where n ∈ {64, 128, 256, 512, 1024, 2048, 4096, 8192}.
Although the latter range is restricted to powers of two for technical reason, Figure 5(a) indicates that
both M(Ψirdwt , x∗ ) and M(Ψdwt , x∗ ) = δ(Ψdwt , x∗ ) scale almost logarithmically in n. Compared to
that, the prediction of [KRZ15, Thm. 3.3] rather exhibits a linear growth. To assess the mismatch of our
bound with the true sample complexity, Figure 5(b) visualizes the relative error

M(Ψirdwt , x∗ ) − δ(Ψirdwt , x∗ )
,
n
showing that the error curve drops to ∼ 4% for increasingly higher dimensions.
We emphasize that the setup of Figure 5 considers an asymptotic model: While the “resolution” of Ψ
gets increasingly finer as n → ∞, the signal x∗ ∈ Rn is always generated by discretizing a piecewise
constant function consisting of exactly 11 discontinuities. Hence, one would expect that the complexity
of x∗ is only slightly affected, and indeed, δ(Ψirdwt , x∗ ) and M(Ψdwt , x∗ ) grow both much more slowly
than n. But in certain sense, this experiment pushes our theoretical approach to its limits because it
is rather designed for a non-asymptotic setting, where n remains fixed. Let us briefly summarize our
findings:
1 The

true sample compexity is numerically computed by means of the statistical dimension [ALMT14], see Appendix B.2 for
details.

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

(a)

17

(b)

Figure 5: Subplot (a) visualizes how the predicted sampling-rate bounds scale as the ambient dimension n is increased, where x∗ ∈ Rn is blocks. Subplot (b) reports the relative mismatch produced by
M(Ψirdwt , x∗ ) with respect to n. Note that both plots use a logarithmic scale on the horizontal axis.

The proposed sampling-rate of Theorem 2.5 scales almost logarithmically in the ambient dimension,
resembling the asymptotics of traditional bounds from compressed sensing. The relative deviation
from the true sample complexity appears to remain constant for larger values of n.
3.1.4 Stable Recovery of Compressible Signals
As discussed in the course of Subsection 2.4, the significance of Theorem 2.5 relies on the assumption
that sufficiently many analysis coefficients are zero. However, perfect sparsity is rarely satisfied in
practice, since real-world signals are usually only compressible (cf. Figure 3). Although one cannot
expect exact recovery in these cases, it still makes sense to ask for an approximate reconstruction of
x∗ ∈ Rn if m is reasonably large, corresponding to a smooth error decay. In this part, we intend
to verify this heuristic numerically and demonstrate how the stability result of Theorem 2.8 can be
applied to such situations.
For this purpose, we continue to use Ψirdwt ∈ R N ×n as analysis operator, whereas the blocks signal
clearly needs to be modified in order to achieve a compressible coefficient vector with respect to Ψirdwt .
Our adapted choice of the source signal, denoted by xc∗ ∈ Rn , is displayed in Figure 6(a) for n = 256.
More precisely, it is constructed by replacing one of the piecewise constant sections of blocks with
a smoother segment, such that the number of non-zero coefficients increases from 906 to 1140, i.e.,
kΨirdwt xc∗ k0 = 1140. Figure 6(b) shows that this modification particularly results in a larger region
of rapidly decaying analysis coefficients (between 800 and 1140). For a more appropriate comparison,
we do not consider the original blocks signal as “sparse” counterpart, but define a simplified version
∗ ∈ Rn which does not have any discontinuities within the curvy segment of x∗ , see Figure 6(a) for
xsp
c
∗ and Figure 1(a) for blocks.
xsp
irdwt
The performance of the analysis basis pursuit (BPΨ
η =0 ) is visualized in Figure 6(c): Here, we have
∗
∗ with M = {80, 81, . . . , 160} and report
applied the instructions of Experiment 3.1 to both xc and xsp
∗ drops to zero at
the averaged reconstruction errors. Similar to Figure 1(b), the estimation error of xsp
∗
m ≈ 95. In contrast, xc is not exactly recovered before m ≈ 150, but starting from m ≈ 100, the error
curve however smoothly tends to zero. This underpins our intuition that one cannot expect a perfect,
but still an accurate, outcome for compressible signals. For illustration, we have also plotted two
irdwt
exemplary reconstructions of xc∗ via (BPΨ
η =0 ) in Figure 6(e) and Figure 6(f) for m = 100 and m = 125,
respectively.
To invoke our theoretical framework from Subsection 2.4, we first need to come up with a meaning-

18

3 A PPLICATIONS AND N UMERICAL E XPERIMENTS

(a)

(b)

(c)

(d)

(e)

(f)

∗ . The projection P x∗ yields
Figure 6: (a) Compressible signal xc∗ and its piecewise constant counterpart xsp
US c
∗
1
a low-complexity approximation of xc . (b) Sorted ` -normalized magnitudes of the analysis coefficients
∗ . (c) Averaged recovery error via (BPΨirdwt ). (d) Relative recovery error predicted by
Ψirdwt xc∗ and Ψirdwt xsp
η =0
∗
m
irdwt
Theorem 2.8 (see Experiment 3.2). (e) + (f) Minimizers x̂ of (BPΨ
η =0 ) with input y = Axc ∈ R .

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

19

ful approximation strategy for xc∗ . More precisely, our task is to identify a subspace US ⊂ Rn , such that
the rescaled projection (cf. (2.9))
x̄c∗ :=

kΨirdwt xc∗ k1
· P x ∗ ∈ Rn
kΨirdwt PUS xc∗ k1 US c

is both of lower complexity and the resulting approximation error k xc∗ − x̄c∗ k2 is small. If US is appropriirdwt
ately chosen, Theorem 2.8 certifies accurate recovery via (BPΨ
η =0 ), while the involved sampling-rate
function M(Ψirdwt , x̄c∗ ) is of moderate size. As already pointed out in the course of (2.8), finding a best
S-analysis-sparse approximation of xc∗ with respect to Ψirdwt is not straightforward and may lead to an
NP-hard problem. Hence, we merely propose a greedy method in Appendix B.1: Roughly speaking,
Algorithm B.1 provides a subspace US ⊂ Rn for a given S ∈ [ N ] that tries to meet the above two criteria, even though this choice might be suboptimal. Figure 6(a) shows the “evolution” of our approach
by plotting an intermediate approximation PUS xc∗ with S = kΨirdwt PUS xc∗ k0 = 900. The algorithm indeed first targets the area of highest curvature, approximating xc∗ by piecewise constant sections here,
so that the analysis coefficients of x̄c∗ become sparser. However, we suspect that an optimal projection
would rather replace the smooth segment in a more “zig-zag” like fashion.
To verify the statement of Theorem 2.8, we apply the steps of Experiment 3.2 below. Note that this
template permits any subspace US with kPUS xc∗ k0 ≤ S.
Experiment 3.2 (Approximation of a compressible signal)
Input: Signal vector xc∗ ∈ R256 , analysis operator Ψirdwt ∈ R1792×256 , approximation subspace
US ⊂ Rn for every S ∈ [ N ].
Compute: Repeat the following procedure 50 times for every m ∈ {80, 81, . . . , 160}:

I Starting with S = 1140, decrease S until M(Ψirdwt , PUS xc∗ ) ≤ m is satisfied.
I Determine x̄c∗ :=

kΨirdwt xc∗ k1
kΨirdwt PUS xc∗ k1

· PUS xc∗ according to (2.9).

I Compute and store the approximation error k xc∗ − x̄c∗ k2 .
The results of this simulation are presented in Figure 6(d), plotting the error curve over the selected
range for m. To illustrate the benefit of our greedy method, Algorithm B.1, we have repeated the instructions of Experiment 3.2 using a standard approximation strategy, i.e., one sets US := ker[(Ψirdwt )S c ]
where S just corresponds to the S largest coefficients of Ψirdwt xc∗ in magnitude. Note that the plots of
Figure 6(d) are only suited for a relative comparison between both approaches. In fact, we have disregarded the impact of the tuning parameter R in Theorem 2.8, so that the respective curves do rather
describe the qualitative behavior of the error bound (2.11). This particularly explains why the labels
on the vertical axis are omitted in Figure 6(d).
For m ≥ 100 measurements, the predicted curve based on the greedy choice of US strongly resembles
the true recovery error in Figure 6(c). This stands in contrast to the naive approach (‘largest coeff.’),
which is apparently not able to exploit the compressibility of xc∗ . The unfavorable “jump” at m ≈ 138
is due to the ignorance of linear dependencies within Ψirdwt : For example, when enforcing that

hPUS xc∗ , ψk i = 0,
for a high-scale wavelet frame vector ψk ∈ Rn of Ψirdwt , the same orthogonality relation does automatically hold true for many more coefficients at lower scales. Such a clustering of coefficient eventually
leads to poor `2 -approximations of xc∗ . Our greedy selection procedure avoids this drawback by implicitly respecting the multilevel structure of wavelets; see also wavelet footprints [DV03]. Let us again
draw a brief conclusion:

20

3 A PPLICATIONS AND N UMERICAL E XPERIMENTS

The recovery error of analysis compressible signals typically tends smoothly to zero as m grows.
The stability result of Theorem 2.8, combined with an appropriate greedy approximation scheme,
adequately reflects this decay behavior.
3.1.5 Phase Transition for Piecewise Constant Signals
Up to now, we have only analyzed a specific signal (blocks) with respect to certain Haar wavelet
frames. Although the resolution of x∗ ∈ Rn was adapted to different dimensions n, its geometric
features (the number of piecewise constant segments) remained untouched. On the other hand, it
is also of interest to assess how well our main result, Theorem 2.5, predicts the phase transition of
∗
(BPΨ
η =0 ) if the “complexity” of x is changed. In the classical situation of an orthonormal analysis
operator (e.g., Ψdwt ), the success of recovery is completely determined by the sparsity S = kΨx∗ k0 ,
implying that it is quite natural to create phase transition plots over S, see [ALMT14]. But such a
simple indicator for complexity does often not exist for (highly) redundant systems, which is in fact
one of the key findings of this work. However, recalling that wavelets were specifically developed for
the detection of singularities, one may argue that the number of jumps characterizes a signal much
better when working, for example, with Ψirdwt . The following experiment shows that this heuristic
of total variation sparsity (TV-sparsity)1 indeed serves as an appropriate surrogate of analysis sparsity in
the context of piecewise constant functions.
Our first simulation is generated according to the following template with n = 128, Ψ = Ψirdwt , and
S = {1, 2, . . . , 128}:
Experiment 3.3 (Phase transition for piecewise constant signals)
Input: Fixed ambient dimension n, analysis operator Ψ ∈ R N ×n , range of TV-sparsity
S ⊂ {0, 1, . . . , n − 1}.
Compute: Repeat the following procedure 50 times for every STV ∈ S :

I Select a random set STV ⊂ [n − 1] with |STV | = STV and determine an orthonormal basis B of
∗
c ]. Then draw a standard Gaussian random vector c and set x := Bc.
ker[(ΨTV-1 )STV
I Repeat the following steps 10 times for each m ∈ {1, 2, . . . , n}:
B Draw a standard i.i.d. Gaussian random matrix A ∈ Rm×n and determine the measurement vector y = Ax∗ .
n
B Solve the analysis basis pursuit (BPΨ
η =0 ) to obtain an estimator x̂ ∈ R .

B Declare success if k x∗ − x̂k2 < 10−5 .
Two examples of the signal-generation step in Experiment 3.3 are shown in Figure 7. The actual phase
transition plot of Figure 8(a) is then created by computing the empirical mean of the success rates from
Experiment 3.3. The annotated curves visualize the sampling-rate bounds proposed by Theorem 2.5
and [KRZ15, Thm. 3.3], respectively. For the second run in Figure 8(b), we have invoked Experiment 3.3
with higher dimension n = 512 and Ψ = Ψirdwt , but using a slightly coarser grid S = {1, 5, 9, . . . , 509}
to reduce the computational burden.
Interestingly, the results of Figure 8 do strongly resemble classical phase transitions, e.g., as reported
in [ALMT14]. This observation is somewhat surprising because the (averaged) coefficient sparsity
S = kΨirdwt x∗ k0 , displayed on the top of the plots in Figure 8, appears detached in our setting. Regarding accuracy, we can conclude that M(Ψirdwt , x∗ ) captures the location of the phase transition
fairly well, whereas [KRZ15, Thm. 3.3] provides a much worse prediction. Thus, the key message of
this subsection reads as follows:
TV-sparsity (or gradient sparsity) of x∗ ∈ Rn is given by STV := kΨTV-1 x∗ k0 , where ΨTV-1 ∈ R(n−1)×n is the standard total
variation operator in 1D. For a precise definition, see Subsection 3.2.

1 The

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

(a)

21

(b)

Figure 7: Two signal vectors x∗ ∈ Rn generated by the first step of Experiment 3.3 with different choices of
TV-sparsity.

(a)

(b)

Figure 8: Phase transition plots for piecewise constant signals analyzed with Ψirdwt (STV is the number of
discontinuities of x∗ ∈ Rn ). The blue curve is obtained by computing the empirical mean of the samplingrate function M(Ψirdwt , x∗ ) within each iteration step STV ∈ S of Experiment 3.3. The same procedure was
performed for the orange curve corresponding to [KRZ15, Thm. 3.3].

The number of discontinuities governs the complexity of piecewise constant signals, enabling us to
create appropriate phase transition plots for redundant Haar wavelet frames. The resulting transition curves are very accurately described by our sampling-rate bound in Theorem 2.5.
Remark 3.4 At this point, it is also worth revisiting Subsection 3.1.3: The simulation of Figure 5
investigates the behavior of the sample complexity when the TV-sparsity of x∗ ∈ Rn remains constant
while n grows. Conceptually, this corresponds to invoking Experiment 3.3 for very large values of
n and studying the left end of the resulting phase transition plot. The relative error in Figure 5(b)
therefore particularly reflects the small deviation of our prediction from the truth in Figure 8 if STV is
small. In practice, however, the choice of Ψ is usually adapted to the signal’s resolution level, implying
that one is rather interested in those vertical cross sections of transition plots for which the ratio STV /n
is of moderate size.
♦

22

3 A PPLICATIONS AND N UMERICAL E XPERIMENTS

(a)

(b)

TV-1
Figure 9: Phase transition plots for (BPΨ
η =0 ) using (a) signals built of random piecewise constant segments
and (b) signals built of “dense jumps.” The blue curve is obtained by computing the empirical mean of
M(ΨTV-1 , x∗ ) for every iteration step STV ∈ S of Experiment 3.3, and the orange curve corresponds to
[KRZ15, Thm. 3.1].

3.2 Total Variation
In this section, we consider a fundamentally different, yet classical example of an analysis operator,
namely total variation in 1D (TV-1). It originates from the seminal work of [ROF92], investigating the
problem of signal denoising. Although conceptually quite simple, total variation has proven to be a
highly effective prior in regularizing inverse problems and therefore became one of the most popular
analysis operators used in practice.
Perhaps, the most striking structural difference to the wavelet-based approach of the previous subsection is that the total variation operator


−1 1
0 ... 0
 0 −1 1
0 0


(n−1)×n
ΨTV-1 :=  .
. ∈ R
..
..
 ..
. .. 
.
0

...

0

−1 1

does not constitute a frame for Rn . Note that there exist numerous variants of total variation imposing different boundary conditions. The above choice uses forward differences with von Neumann
boundary conditions and is common in the field of compressed sensing.
To assess the quality of Theorem 2.5, it makes sense to follow the same strategy as in Subsection 3.1.5.
Ψ
Indeed, the TV-based basis pursuit (BPη =TV-1
0 ) promotes piecewise constant output signals, which indicates that TV-sparsity is again the right quantity to study. Our first phase transition plot in Figure 9(a)
is generated according to Experiment 3.3, with n = 300, Ψ = ΨTV-1 and S = {1, 4, . . . , 298}. The second simulation in Figure 9(b) repeats the same procedure, but the first step of Experiment 3.3 is now
replaced by

I Generate x∗ ∈ Rn by setting x ∗j = (−1) j−1 for j = 1, 2, . . . , STV and x ∗j = 0 otherwise.
Even though the respective TV-sparsity coincides in both signal-generation steps, the two transition
curves of Figure 9 look somewhat different. This observation has a remarkable implication: While
TV-sparsity is usually considered to be a heuristic measure of complexity for piecewise constant funcΨ
tions, it does not entirely characterize the capacity of (BPη =TV-1
0 ). Similar to the example of wavelets
in Subsection 3.1, we therefore suspect that a sound analysis of the total variation operator must be

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

23

non-uniform to a certain extent, depending on the specific structure of x∗ . In fact, the sampling-rate
function x∗ 7→ M(ΨTV-1 , x∗ ) seems to meet this desire for non-uniformity, since the shape of the blue
curves in Figure 9 adapts to each of the cases. But we have to clearly confess that our predictions are
only reliable in the regimes of mid- and high-level sparsity. In contrast, the bound of [KRZ15, Thm. 3.1]
is not very significant and particularly remains unchanged for both signal classes. Let us conclude our
discussion as follows:
Total variation minimization exhibits a non-uniform recovery behavior, which is also captured by
the sampling-rate bound of Theorem 2.5. However, if the TV-sparsity is very small, our prediction
is much less accurate than for redundant wavelet frames.
Remark 3.5 (Related Literature) Compared to redundant wavelet frames, the `1 -analysis formulation
of total variation is actually much better theoretically understood. For example, the work of Cai and Xu
√
Ψ
[CX15] proves that the optimal sampling rate of (BPη =TV-1
0 ) is given by m = Θ ( STV · n ) for Gaussian
measurements, where STV = kΨTV-1 x∗ k0 and log-terms are ignored. The corresponding proofs are
also based on estimating the Gaussian mean width, although using very different techniques that are
specifically tailored to ΨTV-1 . The guarantees of [CX15] underpin once again the fundamental role
of TV-sparsity, but as pointed out above, a refined non-asymptotic analysis would eventually rely on
additional geometric properties of the ground truth signal. For a more extensive discussion of total
variation in compressed sensing, see [KKS17] and the references therein.
♦

3.3 Tight Random Frames
Proposition 2.7 confirms the optimality of our sampling-rate bound for orthonormal bases, and indeed,
x∗ 7→ M(Ψ, x∗ ) precisely describes the phase transition curve of (BPΨ
η =0 ) in this case, see [ALMT14,
Prop. 4.5]. However, the situation already becomes more complicated for a closely related class of
analysis operators, namely those generated from standard Gaussian matrices. These constructions
almost surely yield (tight) frames in general position1 and have been widely studied in the literature
as a benchmark example of the analysis formulation [KR15; NDEG13].
In this part, we roughly follow a construction of tight random frames from [NDEG13]: First, draw
an N × n Gaussian random matrix and compute its singular value decomposition UΣV T . If N ≥ n,
we replace Σ by the matrix [ In , 0]T ∈ R N ×n , which yields a tight frame
Ψrand := U ([ In , 0]T )V T ∈ R N ×n .
If N < n, we replace Σ by [ IN , 0] ∈ R N ×n and set Ψrand := U [ IN , 0]V T ∈ R N ×n . But note that this
operator does not form a frame, and in particular, ΨT
rand Ψrand 6 = In .
Our first phase transition plot is created according to Experiment 3.6 below with n = 300, N = 350,
and
S = {(51, N − 51), (54, N − 54), (57, N − 57), . . . , (297, N − 297)}.
The set of tuples is chosen such that for every pair (S, L) ∈ S , it holds S + L = N = 350. Moreover, we
do only consider S > 50, since otherwise ker[(Ψrand )S c ] = {0} due to the general position property of
Ψrand . The experimental result in Figure 10 shows that the sampling-rate bound of Theorem 2.5 almost
perfectly hits the transition curve.
Experiment 3.6 (Phase transition for random analysis operators)
Input: Fixed ambient dimension n, range of sparsity-cosparsity tuples
S ⊂ {0, 1, . . . , N } × {0, 1, . . . , N }.
Compute: Repeat the following procedure 5 times for each pair (S, L) ∈ S :
1A

matrix Ψ ∈ R N ×n with N ≥ n is in general position if every subset of n rows is linearly independent.

24

3 A PPLICATIONS AND N UMERICAL E XPERIMENTS

N ×n . The blue curve is obtained
rand
Figure 10: Phase transition of (BPΨ
η =0 ) for tight random frames Ψrand ∈ R
∗
by computing the empirical mean of M(Ψrand , x ) for every iteration step (S, L) ∈ S of Experiment 3.6.
The labels of the horizontal axis do only display the first component of S .

I Set N = S + L and construct a random operator Ψrand ∈ R N ×n as described above. Select a
random set S ⊂ [ N ] with |S| = S and determine an orthonormal basis B of ker[(Ψrand )S c ].
Then draw a standard Gaussian random vector c and set x∗ := Bc.
I Repeat the following steps 10 times for each m ∈ {1, 2, . . . , n}:
B Draw a standard i.i.d. Gaussian random matrix A ∈ Rm×n and determine the measurement vector y = Ax∗ .
Ψ

n
B Solve the analysis basis pursuit (BPη =rand
0 ) to obtain an estimator x̂ ∈ R .

B Declare success if k x∗ − x̂k2 < 10−5 .
In Figure 10, we have selected the pairs (S, L) such that both n and N remained fixed. The purpose of
the second and third run of Experiment 3.6 is to highlight the non-trivial impact of both parameters by
varying either one of them (and adapting N). Figure 11(a) is obtained by fixing L = 250 and applying
Experiment 3.6 with n = 300 and
S = {(1, L), (4, L), (7, L), . . . , (100, L)}.
In Figure 11(b), we have fixed S = 50 and invoked Experiment 3.6 with n = 300 and
S = {(S, 200), (S, 203), (S, 206), . . . , (S, 299), (S, 300)}.
The transition curves of Figure 11(a) and Figure 11(b) are obviously not constant in S and L, respectively, which verifies once again that neither sparsity nor cosparsity can thoroughly quantify the
recovery performance of `1 -analysis minimization. Note that the plots of Figure 11 do actually consist
of two regions for which the shapes of the curves are somewhat different. This is due to the case distinction between N ≥ n and N < n, i.e., whether Ψrand ∈ R N ×n forms a frame or not. Interestingly,
Figure 11(b) also reveals another drawback of the traditional bound in (1.3): If S remains fixed, the
analysis dimension N = S + L increases with L, so that one would expect a logarithmic growth of
the number of required measurements. But Figure 11(b) shows that the true sample complexity even

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

(a)

25

(b)

N ×n with either L or S fixed. (a) L =
rand
Figure 11: Phase transition of (BPΨ
η =0 ) for random frames Ψrand ∈ R
250 is fixed and the labels of the horizontal axis do only display the sparsity of S . (b) S = 50 is fixed and
the labels of the horizontal axis do only display the cosparsity of S .

decreases for larger values of L. Regarding Theorem 2.5, we can again conclude that the sampling-rate
function M(Ψrand , x∗ ) captures the phase transition almost perfectly in both scenarios. The prediction
of [KRZ15, Thm. 3.3] also reflects the shapes of the transition curves here, but is still slightly worse
than ours.1

3.4 Analysis Operators for 2D Signals
We have only considered one-dimensional signals up to now, in particular, the class of piecewise constant functions in Subsection 3.1 and Subsection 3.2. However, signals in higher dimensions typically
exhibit a richer geometric structure, such as anisotropic features. For example, it has been observed in
the literature that the sample complexity of total variation minimization in 2D scales very differently
compared to its counterpart in 1D [CX15; KKS17; NW13]. Let us therefore also conduct a simple experiment in 2D, using total variation as well as a redundant Haar wavelet frame as analysis operators. In
order to reduce the immense computational burden of creating a phase transition plot, we just restrict
ourselves to a specific image signal here. Figure 12(a) illustrates our choice of x∗ , which is a realistic
brain phantom for magnetic resonance imaging [GLPU12] at a resolution of 100 × 100.
Recovery based on total variation in 2D aims at retrieving images from compressed measurements
by promoting sparse gradients, i.e., piecewise constant signals. As in 1D, there exist multiple variants
of total variation minimization, depending on different boundary conditions and on how the `1 -norm
of the discrete gradient is calculated. To keep the exposition as brief as possible, we focus on the simple
case where an anisotropic total variation operator, based on periodic forward differences, is used and
x∗ corresponds to an ñ × ñ grayscale image. Thus, up to modifications of the boundary values, the
entries of the discrete (forward) gradient in 2D are defined by
"
#

∗
∗
∗) 0 
x
−
x
(∇
x
0
0
1
j,j
j
+
1,j
j,j
:= ∗
(∇ x∗ ) j,j0 :=
.
(∇2 x∗ ) j,j0
x j,j0 +1 − x ∗j,j0
For the sake of convenience, we will now identify the signal domain Rñ×ñ with its canonical column2
wise vectorization Rñ , i.e., the ambient space is of dimension n := ñ2 and x∗ ∈ Rn . The associated 2D
that, although [KRZ15, Thm. 3.3] is only stated for frames, it literally holds true for any choice of Ψ ∈ R N ×n where the
upper frame bound just needs to be replaced by kΨk22→2 .

1 Note

26

3 A PPLICATIONS AND N UMERICAL E XPERIMENTS

(a)

(b)

(c)

2
∼ R100×100 (‘brain phantom’
Figure 12: (a) Visualization of the ground truth image signal x∗ ∈ R100 =
ΨTV-2
Ψirdwt-2
[GLPU12]). (b)+(c) Average recovery error achieved by (BPη =0 ) and (BPη =0 ), respectively.

total variation operator then takes the following form:
ΨTV-2 := [ΨTV-1 ⊗ Iñ , Iñ ⊗ ΨTV-1 ]T ∈ R2(ñ−1)ñ×n .
To define a redundant Haar wavelet frame in 2D, we again rely on the (two-dimensional) inverse
wavelet transform with 2 scales provided by the software package spot [BF13], and denote it by
Ψirdwt-2 ∈ R7n×n . We assess the recovery capability very similarly to Figure 6(c): Experiment 3.1 is
invoked for both Ψ = ΨTV-2 and Ψ = Ψirdwt-2 , where n = ñ2 = 10.000 and M = {1, 26, 51, . . . , 9976};
note only 12 repetitions were performed instead of 50. The plots of Figure 12(b) and 12(c) visualΨ
Ψirdwt-2
ize the averaged error curves for (BPη =TV-2
), respectively. Unfortunately, solving these
0 ) and (BPη =0
minimization problems via cvx with n = 10.000 is not feasible anymore, so that we have applied an
inhouse implementation based on the alternating direction method of multipliers instead; see [BPCPE11]
and references therein.
Similar to the 1D case, it can be observed that the coefficient sparsity kΨirdwt-2 x∗ k0 = 25.897 is
substantially larger than the ambient dimension n = 10.000, even though Ψirdwt-2 is rather mildly
redundant. Nevertheless, there appears to be a sharp phase transition at m ≈ 5.000 measurements in
Figure 12(c). The prediction m ≈ 5.500 of Theorem 2.5 is also quite accurate in this situation, corresponding to a relative error of about 5%. In contrast, our bound is slightly more pessimistic for 2D total
variation (relative error about 10 %) but still quite reliable. This is obviously not the case for [KRZ15],
whose predictions deviate significantly from the true sample complexity. Finally, it is worth mentioning that the above experimental setup may not be suited for a comparison between total variation and

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

27

Haar wavelets in respect of analysis modeling. We suspect that redefining Ψirdwt-2 with more scales
irdwt-2
and a more sophisticated weighting would considerably improve the outcomes of (BPΨ
).
η =0

4 Related Literature
In this part, we survey existing theoretical approaches to the analysis formulation and put them into
context with our findings. Our discussion starts in Subsection 4.1 with a comparison to the synthesis formulation, which is also widely used in compressed sensing and promotes a somewhat different
viewpoint on sparse representations. In Subsection 4.2, we return to our initial concern from Subsection 1.1, presenting several results relying on traditional analysis sparsity (cf. (1.3)). Subsection 4.3 then
points our the importance of cosparse modeling. Finally, more details on the recent work of Kabanava,
Rauhut, and Zhang [KRZ15] are provided in Subsection 4.4.

4.1 Analysis versus Synthesis
One of the cornerstones in the literature on sparsity priors is [EMR07], which was the first work systematically studying the relationship between the predominant synthesis formulation and the analysis
formulation. Instead of solving (BPΨ
η ), the synthesis formulation rather considers the convex program
!
D·

argmin kzk1

subject to

z ∈R N

k ADz − yk2 ≤ η ,

(BPSηD )

where D ∈ Rn× N is a (possibly redundant) dictionary in Rn . The rationale behind this approach is that
the signal vector x∗ ∈ Rn possesses a sparse representation by D, i.e., x∗ = Dz for a coefficient vector
z ∈ R N with kzk0  n.
D
By investigating the respective polytopal geometry of (BPΨ
η ) and (BPSη ), the authors of [EMR07]
develop a theoretical model that describes the differences between both methods. In particular, they
point out the following fundamental issue: While the synthesis approach (BPSηD ) seems to benefit from
a higher redundancy of the dictionary D, it is not clear how (BPΨ
η ) is influenced by the redundancy
degree of Ψ. Theorem 2.5 shows that this structural property indeed has a wide impact, which is
however highly non-trivial. We hope that our results could give rise to further progress in this matter.
The numerical simulations of [EMR07] certify a considerable gap between both strategies and the observed recovery performances indicate that the analysis formulation outperforms its synthesis-based
counterpart in many situations of interest. Although never stated as a general conclusion, the superiority of analysis-based priors is often confirmed in the literature, e.g., see [SF09]. Moreover, depending
on the redundancy of Ψ and D, one may also argue that solving (BPΨ
η ) is computationally more apn
pealing, since the actual optimization takes place in R . In contrast, (BPSηD ) operates on the possibly
much higher dimensional coefficient space of R N . We refer the interested reader to [RSV08] for more
information on synthesis recovery.

4.2 Approaches Based on Analysis Sparsity
The first compressed-sensing-based approach to the analysis formulation was undertaken by the seminal work of Candès, Eldar, Needell, and Randall [CENR11]. While previous results for redundant dictionaries did rather study coefficient recovery in the synthesis formulation, e.g., based on incoherent
frame atoms [RSV08], a major breakthrough of [CENR11] was that such assumptions can be avoided
by switching over to the analysis perspective. The theoretical analysis of [CENR11] relies on the socalled D-RIP, which is an adaptation of the classical restricted isometry property to sparse representations
in dictionaries:

28

4 R ELATED L ITERATURE

Definition 4.1 (D-RIP, [CENR11]) Let D ∈ Rn× N be a tight frame1 for Rn and A ∈ Rm×n . Then, A
satisfies the D-RIP with parameters δ > 0 and S ∈ [ N ] if

(1 − δ)k Dzk22 ≤ k ADzk22 ≤ (1 + δ)k Dzk22
holds true for all S-sparse vectors z ∈ R N .
Supposed that A fulfills the D-RIP with constant δ < 0.08, it was shown in [CENR11] that any
T
minimizer x̂ ∈ Rn of (BPΨ
η ) with Ψ := D obeys

k x̂ − x∗ k2 ≤ C0 · η + C1 ·

σS (Ψx∗ )1
√
,
S

(4.1)

where C0 , C1 > 0 are numerical constants and σS (Ψx∗ )1 denotes the best S-term approximation error
of the coefficient vector Ψx∗ . This result was the starting point of several generalizations and refinements: For instance, [Fou14] has extended the above setup to arbitrary frames and Weibull matrices,
incorporating an adapted robust null space property. The work [LML12] modifies the concept of DRIP in a such way that (BPΨ
η ) can be applied to arbitrary dual operators of frames. Several achievement
towards practical applicability were also made by [KNW15], allowing for structured measurements if
a certain incoherence condition on the dictionary D and the subsampled measurement matrix is satisfied. Interestingly, guarantees of a similar flavor were even proven in an infinite-dimensional setting
by [Poo17], using different proof techniques.
Regarding the measurement model of this work (Model 2.1), it has turned out that a standard Gaussian matrix A ∈ Rm×n fulfills the D-RIP with high probability provided that the number of observations obeys m & S · log( N/S). In this case, we again obtain a uniform error bound of the form (4.1). A
similar sampling rate for Gaussian measurements was also recently achieved by Kabanava and Rauhut
in [KR15]. Their proof is however based on yet another statistical tool, namely a modified version of
Gordon’s Escape Through a Mesh—note that we also make use of this fundamental principle in Theorem A.1. In a nutshell, the first (non-uniform) guarantee of [KR15] reads as follows: Given x∗ ∈ Rn
with kΨx∗ k0 ≤ S, a Gaussian measurement matrix A ∈ Rm×n satisfying
m2
2b
≥
·S·
m+1
a

r
log



e· N
S



q

+

a·log(ε−1 )
b·S

2
(4.2)

enables recovery of x∗ via (BPΨ
η =0 ) with probability at least 1 − ε. It is worth mentioning that, compared
to (4.1), this bound does not involve any unspecified numerical constants. Moreover, by introducing
a generalized null space property, the statement of (4.2) can be even extended to uniform recovery
[KR15, Thm. 9].
Compared to our sampling-rate bound in Theorem 2.5, a condition of the form m & S · log( N/S)
is quite attractive and intuitive from an aesthetic viewpoint, since it mimics traditional results in compressed sensing. This resemblance is not very astonishing, since most proofs in the literature do actually “operate” on the space of analysis coefficients R N and then just pull back the corresponding
estimates to the signal domain Rn . For various examples of less redundant analysis operators, e.g.,
tight random frames, such a strategy indeed seems to provide accurate predictions of the required
number of measurements.
However, as already pointed out in the course of Subsection 1.1, the situation becomes much more
delicate for highly redundant and coherent systems, in which the above approaches do by far not
explain the capacity of (BPΨ
η ). The theoretical framework of Section 2 and the preceding experiments
in Section 3 demonstrate, for instance, that the mutual coherence structure of an operator Ψ deserves
special attention as well. Apart from that, our results do not assume that Ψ forms a spanning system
contrast to our convention in Subsection 1.3, the columns of D form a frame for Rn here, meaning that D is actually the
associated synthesis operator.

1 In

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

29

for Rn , which in turn was an important hypothesis in [CENR11; KR15]. This observation particularly
gives rise to doubts that conditioning-related quantities, such as the ratio b/a in (4.2), need to appear
explicitly in a sampling-rate bound.
Let us finally turn to the issue of stability: The error estimate of (4.1) relies on the best S-term approximation error. Since this parameter is completely determined by the analysis coefficient sequence Ψx∗ ,
the linear dependencies within Ψ are again completely disregarded. The compressibility experiment of
Figure 6 shows that this methodology can be problematic for redundant systems. On the other hand,
our approach in Subsection 2.4 is based on approximations in the signal space Rn , and when combined
with an appropriate greedy scheme (e.g., Algorithm B.1), it yields much more reliable outcomes.

4.3 Approaches Based on Cosparse Modeling
A very influential contribution to analysis-based priors in inverse problem theory is the work of Nam,
Davies, Elad, and Gribonval [NDEG13], which has coined the term cosparse analysis model. Perhaps,
this recent branch of research has grown out of the same observation as we have made above:
“The fact that this representation may contain many non-zeroes (and especially so when
N  n) should be of no consequence to the efficiency of the analysis model.” [NDEG13, p.
35]
In fact, a key insight of this paper was that the viewpoints of analysis and synthesis sparsity are actually based on completely different signal models and therefore also draw their strength from different
features. While the synthesis approaches puts its emphasis on the number of non-zero coefficients,
[NDEG13] argues that the analysis model is rather specified by the number of vanishing coefficients,
that is, the cosparsity. Indeed, the class of L-cosparse signals can be simply written as a union of subspaces
[
WΛ ,
(4.3)
Λ⊂[ N ],
|Λ|= L

where WΛ := ker ΨΛ . A careful comparison of both setups reveals that one should not—as suggested by the approaches presented in Subsection 4.2—regard the analysis formulation from a synthesis perspective, thereby treating “the analysis operator as a poor man’s sparse synthesis representation” [NDEG13, p. 35]. For example, the subspace collection of (4.3) may become extremely complicated if Ψ corresponds to a redundant frame with strong linear dependencies, and it is quite evident
that its geometric arrangement is then not just reflected by the value of S (or L).
Towards a theoretical analysis, [NDEG13, Thm. 7] establishes the following equivalence: Every x∗ ∈
n
R with cosupport Λ = S xc∗ ⊂ [ N ] is a unique minimizer of (BPΨ
η =0 ) if, and only if,
sup |hΨΛc z, sign(ΨΛc x)i| < kΨΛ zk1

for all z ∈ ker( A) \ {0}.

(4.4)

x ∈Rn
Ψ Λ x =0

We note that similar conditions were also studied for stable recovery in [VPDF13]. Moreover, it is
worth mentioning that the vector sign(ΨΛc x) plays a central role in (4.4), which reminds us of our
generalized sparsity term SSσxx from Definition 2.3. In contrast, such an expression is entirely missing in
the concept of D-RIP.
Regarding sample complexity, it has turned out that the quantity
κΨ ( L) := max dim WΛ
|Λ|≥ L

essentially determines the (optimal) number of measurements required for noiseless recovery via a
combinatorial search over all L-cosparse signals [NDEG13, Sec. 3]. This finding is particularly consistent with the philosophy of compressed sensing after which an estimation succeeds as long as the

30

4 R ELATED L ITERATURE

sampling rate slightly exceeds the signal’s manifold dimension. Somewhat surprisingly, the recent
work of Giryes, Plan, and Vershynin [GPV15] has shown that the situation is unfortunately not that
simple in `1 -analysis minimization: When asking for robustness or a tractable algorithm instead of
combinatorial searching, the sample complexity is not appropriately captured by κΨ ( L) anymore. For
that reason, we hope that our theoretical framework sheds more light on this issue, featuring the complex relationship between sparsity and cosparsity in the study of (BPΨ
η ).

4.4 Non-Asymptotic Sampling-Rate Bounds
To the best of our knowledge, the only work that allows for a direct comparison to our results is by
Kabanava, Rauhut, and Zhang [KRZ15], hence often serving as a “competitor” in Section 3. In fact,
it is also motivated by a similar observation as we have made in this paper: Bounds on the number
of needed measurements for frame-based analysis operators (cf. (4.2)) and total variation [CX15] may
render vacuous statements if the sparsity level is too high.
Focusing on the case of Gaussian measurements, Theorem 3.3 in [KRZ15] provides a non-asymptotic
∗
sampling-rate bound when Ψ ∈ R N ×n is a frame: Any minimizer x̂ ∈ Rn of (BPΨ
η ) satisfies k x̂ − x k2 ≤
2η/τ with probability 1 − ε, provided that
s
2
q
2(∑i∈S c kψi k2 )2
m2
x∗
−1 ) + τ
,
(4.5)
≥
n−
+
2
log
(
ε
π ·b· N
m+1
where b is an upper frame bound of Ψ. An analogous statement holds true if Ψ is the total variation
operator in 1D (cf. Subsection 3.2), where (4.5) is replaced by
s 
2

2  q
m2
1
S +1
−
1
≥
n· 1− π 1− N
+ 2 log(ε ) + τ ,
m+1
see Theorem 3.1 in [KRZ15] for details. Ignoring terms of lower order, both bounds involve the dimension n of the ambient space as an additive term. This resemblance to the sampling-rate function
M(Ψ, x∗ ) from Definition 2.4 is actually not very surprising, since the proof strategy of [KRZ15] builds
upon the same principle as in this work, namely estimating the conic mean width.
However, there are several substantial differences between the approach of [KRZ15] and ours:
• The results of [KRZ15] are not consistent with classical compressed sensing theory: If Ψ ∈ Rn×n
is an orthonormal basis, (4.5) basically degenerates into


2S2
m ≥ 1 − π2 · n + 4S
π − π ·n ,
which is clearly suboptimal for small values of S.
• The guarantees of [KRZ15] are specific to Gaussian measurements and do not include stable
recovery that allows for coping with compressible signals.
• In all cases examined in Section 3, our bounds are vastly superior in terms of predicting the phase
transition behavior of (BPΨ
η =0 ); see Figure 8, Figure 9, and Figure 11.
Comparing the sampling-rate bound (2.4) of Theorem 2.5 with (4.5), we observe that the latter does
neither incorporate the cross-correlations between the individual analysis vectors nor the sign vector
σx∗ = sign(Ψx∗ ). On the other hand, both features seem to be essential for achieving sound recovery
results. A careful study of the corresponding proofs reveals the underlying reason for the superiority
of our approach: While Step 1 in the proof of Theorem 6.8 in Subsection 6.3 is rather standard (cf.
[ALMT14]), our argumentation in Step 2 is quite different. Indeed, the dual vector w0 in [KRZ15] is
chosen such that it just maximizes the term T1 in (6.11). Our choice of w0 in (6.10) is more sophisticated,

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

31

particularly taking the variable τ into account. This appears to be a minor issue at first sight but
actually affects the remainder of the proofs drastically. For example, we had to come up with a highly
non-trivial estimate for T2 in Lemma 6.12, eventually leading to the generalized cosparsity term LS c∗
x
(see Definition 2.4). Furthermore, the optimization over τ is very challenging in our case, whereas it is
just a simple quadratic problem in [KRZ15].

5 Discussion and Outlook
One of the key findings of this work was that the widely-used concept of analysis sparsity is oftentimes not capable of explaining the success (and failure) of `1 -analysis recovery. This observation defies conventional wisdom that would justify the superiority of an operator Ψ ∈ R N ×n only by means
of sparsity, i.e., smaller values of S = kΨx∗ k0 lead to better outcomes. The definition of our three
generalized sparsity parameters (see Definition 2.3), in contrast, is somewhat more delicate and particularly involves the Gram matrix G = ΨΨT ∈ R N × N . In that way, the mutual coherence structure
of Ψ is also taken into account, which is in turn a missing feature of most traditional approaches. This
important refinement enabled us to prove much more reliable bounds on the sample complexity of
∗
N
the analysis basis pursuit (BPΨ
η ), even when a coefficient sequence Ψx ∈ R is just compressible (see
Subsection 2.4). Hence, we can conclude that our main results, Theorem 2.5 and Theorem 2.8, provide
fairly general answers to the initial challenges of (Q1) and (Q2), as stated in Subsection 1.1.
A more technical yet crucial role in these statements is played by the sampling-rate function M(Ψ, x∗ ),
which essentially determines the minimal number of measurements needed to invoke the actual error
estimates. Our numerical simulations in Section 3 have demonstrated that this quantity indeed very
accurately localizes the phase transition of many different types of analysis operators, such as redundant Haar wavelet systems, total variation, or random frames. This variety certifies once again that
the proposed bounds meet the desiderata of (i)–(iv) from Subsection 1.2 in many situations of interest.
Thus, let us make the following closing remark:
The functional x∗ 7→ M(Ψ, x∗ ) can be regarded as a surrogate measure of complexity, quantifying
how well Ψ captures the low-dimensional structure of a signal x∗ . Compared to the number of nonzero analysis coefficients, its value has turned out to be an appropriate indicator for the success or
failure of `1 -analysis minimization.

5.1 Practical Scope and Guidelines
The previous sections have rather focused on the predictive power of our theoretical framework,
whereas its practical implications were only marginally addressed so far. Therefore, we shall now
investigate the issues of (Q3) in greater detail. Taking the perspective of a practitioner, these problem
statements are typically associated with the following tasks: How to come up with a good analysis operator
Ψ for my specific application? What properties and design rules are of particular relevance? Unfortunately, we
are not able to answer these questions in full generality here, but our theoretical approach could at least
give rise to several novel solution strategies. The above discussion in fact suggests to use M(Ψ, x∗ ) as
a measure of quality that allows for an assessment and comparison of different choices of Ψ. Let us
a recall a simple example from Subsection 3.1.2: In the course of relation (3.2), it has turned out that,
although Ψrdwt and Ψirdwt correspond to the same frame except from a scale-wise reweighting of their
analysis vectors, the predicted sampling rates differ dramatically, as illustrated in Figure 1(b). This
observation indicates that the recovery performance might be greatly improved by weighting, which
is a principle often referred to as weighted sparsity in the literature [FMSY12; OKH12]. The numerical
results of Section 3 even show that such a rule-of-thumb could be applied to various other structural
properties of analysis operators, for example:
• Intrinsic localization of the Gram matrix. The off-diagonal decay of G = ΨΨT may strongly affect
the recovery capability of Ψ (cf. Remark 6.15).

32

6 P ROOFS OF M AIN R ESULTS

• Conditioning. The impact of the linear dependency structure and frame bounds of Ψ is expected
to be highly non-trivial.
• Redundancy. Highly redundant frames are sometimes superior over orthonormal bases (see Figure 1).
While analyzing M(Ψ, x∗ ) with respect to certain characteristic quantities yields a first approach to
(Q3), a systematic methodology for designing good operators is still missing. Interestingly, the language of statistical learning theory can be very helpful to tackle this challenge. The underlying mathematical formalism is sketched in Appendix C, but a detailed elaboration is however deferred to future
research.

5.2 Open Problems
Beside the emerging challenge of the previous subsection, we would like to state some further open
issues which might be studied in future works as well:
• Sharper bounds and optimality. Our analysis of the total variation operator in Subsection 3.2 has
demonstrated that the sampling-rate bound of Theorem 2.5 could be somewhat inaccurate in situations of low sparsity. While there are indeed several non-sharp estimates in the proof of Theorem 6.8 (e.g., see Proposition 6.9 and Lemma 6.12), the most heuristic step is actually our choice
of the dual vector w0 in (6.10), see also Remark 6.15. This leaves certain space for refinements,
but note that, no matter how w0 is chosen, one eventually has to face a non-trivial minimization
problem over τ in (6.8). Apart from that, it would be interesting to specify conditions on Ψ and
x∗ under which our bounds are already optimal, at least in an asymptotic sense.
• Greedy algorithms. The analysis basis pursuit is just a very popular example of algorithms that
try to exploit analysis sparsity in signal estimation. In fact, there exist numerous alternatives
to solving (BPΨ
η ), for instance, the greedy analysis pursuit [NDEG13] or those of [GNEGD14]. We
believe that our proof techniques can be adapted to many of these methods, ultimately leading
to similar recovery guarantees.
• Structured measurements. Our argumentation in Section 6 does strongly rely on the row-independence
of the measurement matrix A ∈ Rm×n . This assumption is unfortunately not satisfied for many
types of structured measurement schemes, so that our statistical tools from Appendix A.1 do not
apply anymore. In particular, the concept of Gaussian mean width might become inappropriate,
implying that most parts of our proofs break down. For such situations, one probably would
have to come up with a completely new proof strategy.
• Analysis versus synthesis. We have already discussed the crucial difference between analysisand synthesis-based priors in Subsection 4.1, but a fundamental question remains widely open:
Which of the two formulations is more appropriate for a given signal class and measurement
ensemble?

6 Proofs of Main Results
Before presenting all proofs in detail, let us briefly sketch the roadmap of this section: In Subsection 6.1,
we first provide an abstract framework, which even applies to a generalized version of the basis pursuit. The associated recovery results rely on variants of the Gaussian mean width, which forms a key concept in analyzing the sample complexity of many convex estimation programs [ALMT14; CRPW12].
These quantities are however very implicit and hard to compute in general. For that reason, we show
in Subsection 6.2 (Theorem 6.8) that the square of the Gaussian mean width of the descent cone of
x 7→ kΨxk1 at x∗ is at least upper bounded by our sampling-rate function M(Ψ, x∗ ), introduced in
Definition 2.4. Combining this key step with the guarantees from Subsection 6.1, we will immediately

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

33

obtain the statements of Theorem 2.5 and Theorem 2.8. The more technical proof of Theorem 6.8 is
then postponed to Subsection 6.3.

6.1 Structured Signal Recovery and Gaussian Mean Width
In this part, we analyze the generalized basis pursuit
min f ( x)

subject to

x ∈Rn

k Ax − yk2 ≤ η,

f

(Pη )

where f : Rn → R is a convex function. Note that for our purposes, it would be sufficient to consider
the `1 -analysis functional f ( x) = kΨxk1 , but in fact, each of the following arguments does literally hold
f
true for the general case. Conceptually, the optimization program of (Pη ) seeks for the most structured
n
∗
signal x ∈ R that is consistent with the measured data y = Ax + e. The optimal sampling rate of this
approach heavily relies on how well the objective function f captures the complexity of the ground
truth signal x∗ ∈ Rn . In order to formalize this intuition, let us first introduce the notion of (Gaussian)
mean width:
Definition 6.1 Let K ⊂ Rn be a set.
(1) The (global) mean width of K is defined as
w (K ) := E[suph g, hi],
h∈K

where g ∼ N (0, In ) is a standard Gaussian random vector.
(2) The conic mean width of K is given by
w∧ (K ) := w (cone(K ) ∩ Sn−1 ).
(3) The local mean width of K at scale t > 0 is defined as
wt (K ) := w ( 1t K ∩ Sn−1 ).
Remark 6.2 If K is convex and 0 ∈ K, the mapping t 7→ wt (K ) is non-increasing, and we also have
wt (K ) ≤ w∧ (K ), since 1t K ⊂ cone(K ) for all t > 0. Moreover, it holds wt (cone(K )) = w∧ (cone(K )) =
w ∧ ( K ).
♦
A particularly interesting choice of K ⊂ Rn is the set of all descent directions of f at a certain “anchor
point”:
Definition 6.3 Let f : Rn → R be a convex function and let x∗ ∈ Rn . The descent set of f at x∗ is given
by
D ( f , x∗ ) := {h ∈ Rn | f ( x∗ + h) ≤ f ( x∗ )},
and its corresponding descent cone is denoted by D∧ ( f , x∗ ) := cone(D ( f , x∗ )). In this context, we refer
to w∧ (D ( f , x∗ )) as the conic mean width of f at x∗ .
Note that, since f is assumed to be convex, both D ( f , x∗ ) and D∧ ( f , x∗ ) are convex sets. The geof
metric idea behind these definitions is quite simple: Every minimizer of (Pη ) is close to x∗ if, and only
n
if, the intersection between the constraint set { x ∈ R | k Ax − yk2 ≤ η } and x∗ + D ( f , x∗ ) has a small
diameter. It was therefore a remarkable observation of [CRPW12] that, for Gaussian measurements,
the occurrence of this event is completely characterized by a single parameter, namely the conic mean
width of f at x∗ . For this reason, w∧ (D ( f , x∗ )) may be also regarded as a measure of complexity for a
signal x∗ with respect to a given structure-imposing functional f .

34

6 P ROOFS OF M AIN R ESULTS

Our first theoretical guarantee actually extends this fundamental principle by incorporating the notion of local mean width. We will see below that this refinement indeed allows us to establish stable
recovery, which is an essential feature of Theorem 2.8.
Theorem 6.4 Let Model 2.1 be satisfied and let f : Rn → R be a convex function. There exists a numerical
constant C > 0 such that, for every fixed t > 0 and u > 0, the following holds true with probability at least
2
1 − e−u /2 : If
2η
 ,
t ≥ √
(6.1)
2
m − 1 − C · κ · wt (D ( f , x∗ )) + u +
f

then any minimizer x̂ of (Pη ) satisfies k x̂ − x∗ k2 ≤ t. In the Gaussian case, a ∼ N (0, In ), we have C = κ = 1.
A proof of Theorem 6.4 and all its corollaries below can be found in Appendix A.1. Roughly speaking, the error bound of Theorem 6.4 states the following: To achieve a desired reconstruction accuracy
f
t > 0 via (Pη ), increase m until (6.1) is fulfilled. While this procedure always succeeds for sufficiently
large values of m, the condition of (6.1) is still quite implicit, since both sides depend on t. Fortunately,
due to Remark 6.2, one may just estimate the local mean width by its conic counterpart. This leads to
our next result, which is well-known from the literature, e.g., see [Tro15, Cor. 3.5].
Corollary 6.5 Let Model 2.1 be satisfied and let f : Rn → R be convex. There exists a numerical constant
C > 0 (the same as in Theorem 6.4) such that, for every u > 0, the following holds true with probability at least
2
f
1 − e−u /2 : Any minimizer x̂ of (Pη ) satisfies

k x̂ − x∗ k2 ≤ √

2η
 .
m − 1 − C · κ 2 · w∧ (D ( f , x∗ )) + u +

In particular, if η = 0 and
m > C2 · κ 4 · w∧ (D ( f , x∗ )) + u
one has x̂ =

2

+ 1,

(6.2)

x∗ .

As already briefly discussed at the beginning of Subsection 2.3, the condition of (6.2) is not only
sufficient but essentially also necessary to achieve exact recovery in the Gaussian case [ALMT14].1
Hence, the conic mean width is the crucial quantity to study when determining the sample complexity
f
of (Pη =0 ).
While Corollary 6.5 already includes robustness against noise, the wish for stable reconstructions
f
via (Pη ), as requested in Subsection 2.4, remains unfulfilled. A major difficulty is that the mapping
x∗ 7→ w∧ (D ( f , x∗ )) is typically discontinuous,
since it involves forming the conic hull of the descent
√
set D ( f , x∗ ). If w∧ (D ( f , x∗ )) ≈ n, the statement of Corollary 6.5 looses its significance, even when
there exists a vector close to x∗ whose conic mean width is much smaller. For example, such a situation
might occur if x∗ just possesses a compressible coefficient sequence which is not too sparse; see also
Subsection 2.4 as well as [Gen17, Sec. III.D] for more details. Our second corollary of Theorem 6.4
resolves this issue by analyzing the local mean width at a specific scale t > 0. It shows that one may
investigate the complexity of f with respect to a surrogate vector x̄∗ ∈ Rn by sacrificing estimation
accuracy in the order of the approximation error k x∗ − x̄∗ k2 .
Corollary 6.6 Let Model 2.1 be satisfied and let f : Rn → R be convex. There exists a numerical constant
C > 0 (the same as in Theorem 6.4) such that, for every R > 0 and u > 0, the following holds true with
2
f
probability at least 1 − e−u /2 : Let x̄∗ ∈ Rn with f ( x̄∗ ) = f ( x∗ ). Then, any minimizer x̂ of (Pη ) satisfies
n
k x̂ − x∗ k2 ≤ max Rk x∗ − x̄∗ k2 ,
1 The

[

√

2η
√
m−1− m̃0 −1]+

o

,

(6.3)

authors of [ALMT14] are focusing on the so-called statistical dimension rather than on the Gaussian mean width, but these
notions are actually equivalent (see [ALMT14, Prop. 10.2]).

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

where
m̃0 := C2 · κ 4 ·



R +1
R

2


· w∧ (D ( f , x̄∗ )) + 1 + u + 1.

35

(6.4)

In particular, if η = 0, we have k x̂ − x∗ k2 ≤ Rk x∗ − x̄∗ k2 provided that m > m̃0 .
f

In a nutshell, Corollary 6.6 certifies a good performance of (Pη ) if both x̄∗ approximates x∗ sufficiently
well and w∧ (D ( f , x̄∗ )) is small. For a more detailed discussion of this surrogate principle, we again
refer to Subsection 2.4, which considers the special case of f (·) = kΨ(·)k1 . In this course, we do also
demonstrate a general selection strategy for x̄∗ : If there exists a subspace U ⊂ Rn which consists
of low-complexity vectors according to x 7→ w∧ (D ( f , x)), it is natural to set x̄∗ as the renormalized
orthogonal projection of x∗ onto U, such that f ( x̄∗ ) = f ( x∗ ). If f is a semi-norm, one may even invoke
the following upper bound on the approximation error, leading to a more explicit expression in (6.3):
Proposition 6.7 Assume that f = k · k is a semi-norm on Rn . Let x∗ ∈ Rn and let U ⊂ Rn be a linear
k x∗ k
subspace such that kPU x∗ k 6= 0. Setting x̄∗ := kP x∗ k · PU x∗ , we have k x∗ k = k x̄∗ k and
U

k x∗ − x̄∗ k2 ≤

kPU x∗ k2
kPU x∗ k

· kPU ⊥ x∗ k + kPU ⊥ x∗ k2 .

Proof. Since x∗ = PU x∗ + PU ⊥ x∗ , we have

k x∗ − x̄∗ k22 =



1−

= 1−

k x∗ k
kPU x∗ k

k x∗ k 2
kPU x∗ k



PU x∗ + PU ⊥ x∗

· kPU x∗ k22

=

kPU x∗ k−kPU x∗ +PU ⊥ x∗ k 2
kPU x∗ k

≤

kPU ⊥ x∗ k 2
kPU x∗ k

≤

x∗ k
U⊥
kPU x∗ k

 kP

2
2

+ kPU ⊥ x∗ k22
· kPU x∗ k22 + kPU ⊥ x∗ k22

· kPU x∗ k22 + kPU ⊥ x∗ k22
2
· kPU x∗ k2 + kPU ⊥ x∗ k2 .



6.2 Proofs of Theorem 2.5 and Theorem 2.8
Let us now focus on the `1 -analysis (semi-)norm f (·) = kΨ(·)k1 . The abstract results of the previous
part indicate that the conic mean width w∧ (D (kΨ(·)k1 , x∗ )) plays a key role in analyzing (BPΨ
η ). While
this quantity may indeed yield an accurate bound on the sampling rate, our desiderata (ii)–(iv) from
Subsection 1.2 are still hardly met, since the definition of the mean width is quite abstract and noninformative. Therefore, we need to come up with a more explicit expression or at least a meaningful
upper bound on w∧ (D (kΨ(·)k1 , x∗ )). This is precisely the purpose of the following theorem, which
forms the main technical ingredient of our approach. Its proof is postponed to Subsection 6.3.
Theorem 6.8 Let Ψ ∈ R N ×n be an analysis operator and let x∗ ∈ Rn with x∗ 6∈ ker Ψ. Then
w2∧ (D (kΨ(·)k1 , x∗ )) ≤ M(Ψ, x∗ ).
Our two main results are now a direct consequence of this bound and the guarantees from Subsection 6.1:
Proof of Theorem 2.5. By Theorem 6.8 and the condition (2.4), we obtain
m>

q

M(Ψ, x∗ ) + u

2


2
+ 1 ≥ C2 · κ 4 · w∧ (D (kΨ(·)k1 , x∗ )) + u + 1,

36

6 P ROOFS OF M AIN R ESULTS

implying that in the ‘in particular’-part of Corollary 6.5 can be applied with f (·) = kΨ(·)k1 . Note that
C2 · κ 4 = 1 here, due to the Gaussianity of A. This already yields the claim of Theorem 2.5.

Proof of Theorem 2.8. By the definition of x̄∗ ∈ Rn in (2.9), we have x̄∗ 6∈ ker Ψ and therefore Theorem 6.8 gives
w2∧ (D (kΨ(·)k1 , x̄∗ )) ≤ M(Ψ, x̄∗ ).
Moreover, we observe that kΨx∗ k1 = kΨ x̄∗ k1 . Hence, Corollary 6.6 is applicable with f (·) = kΨ(·)k1
and we particularly have m̃0 ≤ m0 < m; note that the constants C in the definitions of m̃0 (in (6.4)) and
m0 (in (2.10)) just differ in taking the square. The claim now follows from the error bound of (6.3). 

6.3 Proof of Theorem 6.8 (Bounding the Conic Mean Width)
To complete the proof, it remains to verify Theorem 6.8. For this, we loosely follow [ALMT14, Recipe 4.1],
but note again that we do not attempt to exactly compute w2∧ (D (kΨ(·)k1 , x∗ )) here, but rather establish
an upper bound.
Step 1: A polar bound on the conic mean width
Our first proof step is based on a well-known polarity argument, providing an upper bound on the
mean width of a descent cone. This result even holds true in the very general setup Subsection 6.1.
Proposition 6.9 ([ALMT14, Prop. 4.1]) Let f : Rn → R be convex and let x∗ ∈ Rn . Then, the subdifferential of f at x∗ is well-defined:
∂ f ( x∗ ) := {z ∈ Rn | f ( x) ≥ f ( x∗ ) + hz, x − x∗ i for all x ∈ Rn }.
If ∂ f ( x∗ ) is non-empty, compact, and does not contain the origin, we have1
h
i
w2∧ (D ( f , x∗ )) ≤ inf E
inf k g − τ · zk22 ,
τ >0

z∈∂ f ( x∗ )

(6.5)

where g ∼ N (0, In ).
Remark 6.10 For several special cases, such as norms [ALMT14, Thm. 4.3] or total variation [ZXCL16],
it is known that (6.5) provides a very accurate bound. For a more detailed discussion of the sharpness
of Proposition 6.9, we refer to [ALMT14, Sec. 4.2].
♦
In order to invoke Proposition 6.9, we need to specify the subdifferential of f (·) = kΨ(·)k1 at x∗ .
For the sake of brevity, we now just write S := S x∗ = supp(Ψx∗ ), S c := S xc∗ as well as σ := σx∗ =
sign(Ψx∗ ).
Lemma 6.11 We have
n
o
∂kΨ(·)k1 ( x∗ ) = ΨT σ + ΨT w | w ∈ R N , supp(w) ⊂ S c , kwk∞ ≤ 1 .
In particular, ∅ 6= ∂kΨ(·)k1 ( x∗ ) ⊂ ΨT [−1, 1] N and the following equivalences hold true:
0 6∈ ∂kΨ(·)k1 ( x∗ )

⇔

x∗ 6∈ ker Ψ

⇔

SSσ = kΨT σ k22 > 0.

Proof. First, we apply the chain rule for subdifferentials (see [Roc70, Thm. 23.9])
∂kΨ(·)k1 ( x∗ ) = ΨT (∂k · k1 (Ψx∗ )) ,
1 Hereafter,

the expectation is always taken with respect to the Gaussian random vector g.

(6.6)

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

37

and together with
n
o
∂k · k1 (Ψx∗ ) = z ∈ R N | zS = sign(Ψx∗ ), kzS c k∞ ≤ 1 ,
| {z }
=σ

we obtain that
n
o
∂kΨ(·)k1 ( x∗ ) = ΨT σ + ΨT w | w ∈ R N , supp(w) ⊂ S c , kwk∞ ≤ 1 .
To verify the ‘in particular’-part, suppose that 0 ∈ ∂kΨ(·)k1 ( x∗ ). Then there exists a w ∈ R N with
supp(w) ⊂ S c and kwk∞ ≤ 1 such that 0 = ΨT σ + ΨT w. Thus,
0 = h x∗ , ΨT σ + ΨT wi = hΨx∗ , σ i + hΨx∗ , wi = kΨx∗ k1 ,
| {z } | {z }
=kΨx∗ k

1

(6.7)

=0

that means x∗ ∈ ker Ψ. On the other hand, if x∗ ∈ ker Ψ, we have σ = sign(Ψx∗ ) = 0 and therefore
ΨT σ = 0. Selecting w := 0 shows that 0 ∈ ∂kΨ(·)k1 ( x∗ ). This proves the first equivalence.
By the definition of SSσ (see Definition 2.3), it holds
SSσ = hσ, Gσ i = hσ, ΨΨT σ i = hΨT σ, ΨT σ i = kΨT σ k22 .
Thus, the second equivalence of (6.6) follows from
SSσ = kΨT σ k22 = 0

⇔

ΨT σ = 0

⇔

σ = sign(Ψx∗ ) = 0

⇔

Ψx∗ = 0,

where we have again made use of (6.7) with w = 0 in the second last equivalence.
In conclusion, if x∗ 6∈ ker Ψ, Proposition 6.9 and Lemma 6.11 imply that
h
i
w2∧ (D (kΨ(·)k1 , x∗ )) ≤ inf E
inf
k g − τΨT σ − τΨT wk22 .
τ >0

w ∈R N ,
wS =0,kwS c k∞ ≤1



(6.8)

Step 2: Selecting the `1 -dual vector w.
Evaluating the right-hand side of (6.8) particularly asks us to study the optimization program inside
of the expected value. It can be reformulated as a least-squares problem with box-constraint:
min k g − τΨT σ − τΨT wk22

w ∈R N

subject to

wk = 0 for k ∈ S
and |wk | ≤ 1 for k ∈ S c .

(6.9)

Unfortunately, (6.9) does not seem to have a closed-form minimizer in general. Hence, instead of
seeking for an exact solution, we rather consider an ansatz vector w0 = (w10 , . . . , w0N ) ∈ R N of the form
(
0,
k ∈ S,
0 :
wk =
(6.10)
−
1
clip(τ αk hψk , g i; 1), k ∈ S c ,
where αk > 0 are flexible tuning parameters that will be specified later on. Note that w0 ∈ [−1, 1] N
depends on both g and τ. In Remark 6.15 below, we will also state the KKT-conditions of (6.9), showing
that the choice of (6.10) is at least plausible (although not always optimal) if the Gram matrix of Ψ is
well-localized.

38

6 P ROOFS OF M AIN R ESULTS

Since wk0 is a symmetric transformation of hψk , g i ∼ N (0, kψk k22 ) for every k ∈ S c , we observe that
E[w0 ] = 0. Now, we can proceed with (6.8) as follows:
h
i
E
inf
k g − τΨT σ − τΨT wk22
w ∈R N ,
wS =0,kwS c k∞ ≤1

h
i
≤ E k g − τΨT σ − τΨT w0 k22
i
i
h
h
= E k g − τΨT w0 k22 + E kτΨT σ k22 − E[2h g, τΨT σ i] + E[2hτΨT w0 , τΨT σ i]
i
h
= E k g − τΨT w0 k22 + τ 2 kΨT σ k22 −2hE[ g ], τΨT σ i + 2hτΨT E[w0 ], τΨT σ i
| {z }
|{z}
| {z }
(6.6) σ
= SS

=0

=0

i
i
h
h
= E k g k22 − 2τE[h g, ΨT w0 i] + τ 2 E kΨT w0 k22 + τ 2 · SSσ

= n − 2τE[hΨg, w0 i] + τ 2 E[hw0 , Gw0 i] +τ 2 SSσ = n − T1 + T2 + τ 2 · SSσ .
{z
}
|
{z
} |

(6.11)

=:T2

=:T1

Hence, it remains to resolve the expected values of T1 and T2 .
Step 3: A general bound on the mean width.
Let us start with the term T1 : Using the definition of w0 , we obtain
T1 = 2τE[hΨg, w0 i] = 2τ

∑ c E[hψk , g i · wk0 ]

k ∈S

= 2τ

∑c E

h

hψk , g i · sign(hψk , g i) · min{τ −1 αk |hψk , g i|, 1}

i

k ∈S

=2

∑

h
i
E |hψk , g i| · min{αk |hψk , g i|, τ } ,

k ∈S c

and since hψk , g i ∼ N (0, gk,k ), an elementary calculation shows that

√τ
.
T1 = 2 ∑ αk · gk,k · erf
αk

k ∈S c

2gk,k

Handling the term T2 turns out to be more complicated: Since
such that the clipped variables are standard Gaussians:

√

gk,k = kψk k2 , we first rewrite T2

T2 = τ 2 E[hw0 , Gw0 i]

= τ2

∑
0

k,k

=

∈S c

gk,k0 · E[wk0 · wk0 0 ]

∑
0

h
i
gk,k0 · E clip(αk hψk , g i; τ ) · clip(αk0 hψk0 , g i; τ )

∑
0

gk,k0 · αk · αk0 ·

k,k ∈S c

=

k,k ∈S c

√

h

ψ
gk,k · gk0 ,k0 · E clip h kψ kk , g i; α
k 2

√τ
k gk,k




i
ψ
· clip h kψ k00k , g i; α 0 √τg 0 0 .
k

2

k

k ,k

(6.12)
In general, there does not exist a closed-form expression for the covariance between clipped Gaussian
variables. However, we may invoke the following bound, which is based on Price’s Theorem and
proven in Appendix A.2:

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

39

Lemma 6.12 Let v1 , v2 ∈ Sn−1 and let g ∼ N (0, In ). Fix numbers β 1 , β 2 ≥ 0 and set βmin := min{ β 1 , β 2 },
βmax := max{ β 1 , β 2 }. Then, we have
h
i
E clip(hv1 , g i; β 1 ) · clip(hv2 , g i; β 2 )
q
h


i
min 
max 
−( βmax )2 /2
≤ |hv1 , v2 i| · erf β√ + 1 − erf β√ − π2 e βmax
· βmin · βmax
2
2
h
i
βmin 
= |hv1 , v2 i| · erf √ − h( βmax ) · βmin · βmax ,
(6.13)
2

where h(τ ) =

q

2

2 e−τ /2
π
τ

+ erf( √τ ) − 1; see (2.3). If v1 = v2 , (6.13) holds true with equality.
2

Let us apply Lemma 6.12 to equation (6.12) by setting
v1 :=
β 1 :=
Observing that

√

gk,k0
gk,k · gk0 ,k0

T2 ≤

ψk
,
kψk k2
√τ
αk gk,k ,

ψk0
, and
kψk0 k2
√τ
αk0 gk0 ,k0 .

v2 :=
β 2 :=

= hv1 , v2 i, we end up with

∑

k,k0 ∈S c

where
βmin
k,k0 := min

τβmin
0 
2
√k,k
− h(τβmax
k,k0 ) · τ
2


2
gk,k
erf
0 · αk · αk0 ·

n
αk

√1
√1
gk,k , αk0 gk0 ,k0

o

and

βmax
k,k0 := max


max
· βmin
·
β
,
0
0
k,k
k,k

n
αk

√1
√1
gk,k , αk0 gk0 ,k0

o

.

(6.14)

Combining our findings for T1 and T2 with the estimates of (6.8) and (6.11) yields the following
interim result:
Proposition 6.13 Let Ψ ∈ R N ×n be an analysis operator and let x∗ ∈ Rn with x∗ 6∈ ker Ψ. For fixed numbers
αk > 0, k ∈ S c , it holds
n
w2∧ (D (kΨ(·)k1 , x∗ )) ≤ inf n + τ 2 · SSσ − 2
τ >0

+

∑
0

k,k ∈S c

∑ c αk · gk,k · erf

k ∈S


2
gk,k
erf
0 · αk · αk0 ·

√τ
αk



2gk,k

τβmin
0 
2
√k,k
− h(τβmax
k,k0 ) · τ
2

max
· βmin
k,k0 · β k,k0

o

,

(6.15)

max
where βmin
k,k0 and β k,k0 are defined according to (6.14).

The bound of Proposition 6.13 is in fact quite general and leaves much space for selecting the parameters αk . However, we still need to optimize over τ, which seems to be a very challenging task in
general. Our final step is therefore to make a specific choice of αk , eventually leading to the explicit
bound of Theorem 6.8.

Step 4: Selecting the tuning parameters αk .
In order to simplify the optimization problem of (6.15), we would like to factor out all terms depending
on τ. Such a separation is easily obtained with αk := √λg for k ∈ S c and some λ > 0 which is fixed
k,k
later on. With this choice of αk , (6.15) reads as follows:

40

6 P ROOFS OF M AIN R ESULTS

w2∧ (D (kΨ(·)k1 , x∗ )) ≤ inf

τ >0
λ >0

n

n + τ 2 · SSσ − 2λ ·



∑c

√



gk,k

· erf

τ
√
λ 2

k ∈S

|

+ λ2 ·

{z

}

∑

√

(2.1) r
= LS c



k,k0 ∈S c

|

2
gk,k
0
gk,k · gk0 ,k0

{z



 
· erf

τ
√



λ 2

− h( λτ ) · ( λτ )2

o

,

}

(2.1)

= LS c

and substituting τ by λ · τ leads to
n


w2∧ (D (kΨ(·)k1 , x∗ )) ≤ inf n + SSσ − LS c h(τ ) τ 2 λ2 + λ2 LS c − 2λLr
c erf
S
τ >0
{z
λ >0 |

√τ
2

o

.

(6.16)

}

=:F (τ,λ)

It is not hard to see that (τ, λ) 7→ F (τ, λ) defines a smooth function that can be continuously extended
to [0, ∞) × [0, ∞). The next lemma shows that F has in fact a unique global minimum. Its technical
proof relies on basic calculus and is deferred to Appendix A.3.
Lemma 6.14 Assuming that SSσ , LS c , Lr
S c > 0, the function F : [0, ∞ ) × [0, ∞ ) → R, defined according to
(6.16), attains its unique global minimum at
SSσ
LS c


(τ 0 , λ0 ) := h−1



,

Lr
Sc
LS c



.

Now, we can easily derive the statement of Theorem 6.8:
∗
Proof of Theorem 6.8. First, suppose that S c 6= ∅. This implies that LS c , Lr
S c > 0, and since x 6 ∈ ker Ψ,
σ
we also have SS > 0 by Lemma 6.11. Lemma 6.14 and (6.16) then yield

w2∧ (D (kΨ(·)k1 , x∗ )) ≤ F (τ 0 , λ0 )


= n + SSσ − LS c · h(τ 0 ) ·(τ 0 )2 · (λ0 )2 + (λ0 )2 · LS c − 2λ0 · Lr
S c · erf
|
{z
}
|
{z
}
=0

= n−

( Lr
)2
Sc
LS c

· erf

τ0
√
2



=−( Lr
)2 /LS c
Sc



√1
2

· h −1

SSσ
LS c



= n−

( Lr
)2
Sc
LS c

·Φ



SSσ
LS c



= M(Ψ, x∗ ).

(6.17)

On the other hand, if S c = ∅, we simply observe that
w2∧ (D (kΨ(·)k1 , x∗ )) ≤ lim F (τ, λ) = F (0, 0) = n = M(Ψ, x∗ ).
τ,λ&0



Remark 6.15 (KKT-Conditions and Optimality) A key step in the above proof was to select the dual
vector w0 in (6.10), which eventually enabled us to derive a closed-form bound on the conic mean width
w2∧ (D (kΨ(·)k1 , x∗ )). Such a choice is however not always optimal. Indeed, any minimizer w ∈ R N of
(6.9) fulfills the following KKT-conditions (cf. [MR10]):


τ 2 ΨS c (ΨS c )T + diag(µS c ) wS c = τΨS c g − τ 2 ΨS c ΨT σ,
µk ≥ 0,

k ∈ [ N ],

µk (1 − wk ) = 0,

k ∈ S c,

w2k ≤ 1,

k ∈ S c,

wk = 0,

k ∈ S,

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

41

1
for some µ = (µ1 , . . . , µ N ) ∈ R N . Supposed that G = ΨΨT = diag(α1−1 , . . . , α−
N ) for certain num0
bers αk > 0, it is not hard to see that w = w from (6.10) precisely satisfies these conditions. This
observation particularly explains why our sampling-rate bound is optimal in the case of orthonormal
bases (cf. Proposition 2.7), where G = In .
On the other hand, if the operator Ψ is redundant, i.e., N > n, the Gram matrix cannot be diagonal
anymore and the solution of (6.9) may take a different form. But as long as G is well-localized, in the
1
sense that G = ΨΨT ≈ diag(α1−1 , . . . , α−
N ), one can expect that (6.10) is at least relatively close to the
true minimizer. While this heuristic actually suggests to set αk := 1/gk,k for k ∈ S c , it is not clear how
to proceed with (6.15) in that situation. Therefore, we have made a slightly different choice of αk in Step
4 above. But as our numerical simulations in Section 3 indicate, this strategy still works surprisingly
well for many different classes of analysis operators.
♦

6.4 Proof of Proposition 2.7 (Simplified Sample Complexity Bound)
We prove that M(Ψ, x∗ ) is upper bounded by each of the two expressions on the right-hand side of
(2.5). Recall that the left branch of the minimum is supposed to mimic the (asymptotic) behavior of
σ ∗
σ ∗
M(Ψ, x∗ ) if the ratio SSx∗ /LS c∗ is small, whereas the right branch is appropriate as SSx∗ /LS c∗ becomes
x
x
x
x
large.
σ ∗

Proof of Proposition 2.7. By Lemma 6.11 and x∗ 6∈ ker Ψ, we know that S = SSx∗ > 0, while the assumpx

tion S xc∗ 6= ∅ implies that L = LS c∗ > 0 and L̄ = Lr
S c > 0. Therefore, by Lemma 6.14 and (6.17), we
x∗

x

have

M(Ψ, x∗ ) = F (τ 0 , λ0 ) ≤ F (τ, λ0 )
√
2
for all τ > 0. Since h(τ ) ≥ 0 and erf(τ/ 2) ≥ 1 − e−τ /2 for all τ > 0 (see [FR13, Prop. 7.5]), one can
bound F (τ, λ0 ) as follows:

F (τ, λ0 ) = n + S − L · h(τ ) · τ 2 · ( L̄L )2 −

≤ n + S · ( L̄L )2 · τ 2 −
= n−
Finally, setting τ = τ 00 :=

q

L̄2
L

L̄2
L

· (1 − e − τ

+ ( L̄L )2 · (S · τ 2 + L · e

L̄2
L
2 /2

· erf

√τ
2



)

−τ 2 /2

).

L
2 log( S+
S ) > 0, we obtain

M(Ψ, x∗ ) ≤ F (τ 00 , λ0 ) = n −

L̄2
L


+ ( L̄L )2 · 2S · log( S+S L ) +

L·S
S+ L



|{z}
≤S

≤ n−
It remains to verify that M(Ψ, x∗ ) ≤ n −
this is equivalent to
Φ
ξ:=h−1 (S/L)

⇔

⇔

 
S
L

= erf

erf( √ξ ) ≥
2

2
π



·

√1
2

L̄2

L̄ 2
L + (L)

L̄2
S+ L




· 2S · log( S+S L ) + S .

holds true as well. By the definition of M(Ψ, x∗ ),


2
L
2
· h−1 ( SL ) ≥ ·
= ·
π S+L
π

2
1
·
,
π h(ξ ) + 1

ξ>0

G (ξ ) := erf( √ξ ) · (h(ξ ) + 1) ≥
2

2
,
π

ξ > 0.

S
L

1
+1

42

A P ROOFS OF A UXILIARY R ESULTS

First, we observe that
lim G (ξ ) = lim

ξ &0

h

ξ &0

erf( √ξ ) ·
2

q

2
2 e−ξ /2
π
ξ

i

=

2
.
π

The claim then follows from the fact that G is monotonically increasing, which is in turn a consequence
of G 0 (ξ ) > 0 for all ξ > 0: indeed, we have
q
q


2
2
0 < G 0 (ξ ) = π2 e−ξ /2 · (ξ 2 − 1) · erf( √ξ ) + ξ · π2 · e−ξ /2
2
q
2
ξ
2
2
⇔ 0 < G̃ (ξ ) := (ξ − 1) · erf( √ ) + ξ · π · e−ξ /2 ,
2

while the latter statement holds true due to G̃ (0) = 0 and
G̃ 0 (ξ ) = 2ξ · erf( √ξ ) > 0,
2

ξ > 0.



A Proofs of Auxiliary Results
A.1 Proofs of Recovery Results From Subsection 6.1
First, let us state an important preliminary result from random matrix theory.
Theorem A.1 (Minimum Spheric Singular Value) Assume that Model 2.1 is satisfied and let K ⊂ Sn−1 .
There exists a numerical constant C > 0 such that, for every u > 0, we have
√
sup k Ahk2 > m − 1 − C · κ 2 · (w (K ) + u)
h∈K

with probability at least 1 − e−u

2 /2

. If a ∼ N (0, In ), we have C = κ = 1.

Proof. This statement is a consequence of [LMPV17, Thm. 1.3]. To obtain the optimal constants in the
Gaussian case, one can use [Tro15, Prop. 3.3] instead.

The Gaussian case of Theorem A.1 is well-known as Gordon’s Escape Through a Mesh Theorem, originating from [Gor88]. The above formulation is typically interpreted as a lower bound on the minimum
singular value of the random matrix A ∈ Rm×n restricted to cone(K ).
With this result at hand, we are now ready to prove Theorem 6.4:
f

Proof of Theorem 6.4. Any minimizer x̂ of (Pη ) satisfies x̂ − x∗ ∈ D ( f , x∗ ) as well as k A x̂ − yk2 ≤ η. Our
goal is therefore to show that every vector x ∈ Rn with k x − x∗ k2 > t cannot be a minimizer, if both
x − x∗ ∈ D ( f , x∗ ) and k Ax − yk2 ≤ η are fulfilled.
First, we apply Theorem A.1 for K = 1t D ( f , x∗ ) ∩ Sn−1 : Let h ∈ tK and assume k Ax − yk2 ≤ η
where x := x∗ + h. Without loss of generality, suppose that
√
m − 1 > C · κ 2 · (wt (D ( f , x∗ )) + u),
since otherwise, the condition (6.1) is not satisfied anyway. On the event of Theorem A.1, we have
δ :=

h√

m − 1 − C · κ 2 · (wt (D ( f , x∗ )) + u)

i
+

=

√

m − 1 − C · κ 2 · (w (K ) + u)

< 1t k Ahk2 ≤ 1t (k Ah − ek2 + kek2 ) = 1t (k Ax − yk2 + kek2 ) ≤

2η
.
t

Rearranging the latter inequality gives t < 2η/δ which contradicts the assumption (6.1). This implies
f
that x = x∗ + h cannot be a minimizer of (Pη ).

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

43

f

Now, assume that x̂ ∈ D ( f , x∗ ) + x∗ is a minimizer of (Pη ) with k x̂ − x∗ k2 > t. Setting ĥ := x̂ − x∗ ,
we have t̂ := k ĥk2 > t. The convexity of D ( f , x∗ ) and 0 = x∗ − x∗ ∈ D ( f , x∗ ) imply that there exists
h ∈ tK = D ( f , x∗ ) ∩ tSn−1 such that ĥ = t̂t h. In other words, h is the intersection point between tSn−1
and the line that connects ĥ with 0. Hence, on the event of Theorem A.1, we also observe that
δ < 1t k Ahk2 = 1t̂ k A ĥk2 ≤

2η
2η
,
<
t
t̂
f

which is again a contradiction to (6.1). From this, we can conclude that any minimizer x̂ of (Pη ) satisfies
k x̂ − x∗ k2 ≤ t.

Proof of Corollary 6.5. First, assume that η > 0. We would like to apply Theorem 6.4 with
t := √

m−1−C

· κ2

2η
 .
· (w∧ (D ( f , x∗ )) + u) +

Since w∧ (D ( f , x∗ )) ≥ wt (D ( f , x∗ )) by Remark 6.2, the assumption (6.1) is satisfied and the claim now
follows from Theorem 6.4.
Now, let us consider the noiseless case of η = 0. Similar to the proof of Theorem 6.4, we apply
Theorem A.1 with K = D∧ ( f , x∗ ) ∩ Sn−1 . Let h ∈ K and x := x∗ + h with k Ax − yk2 = k Ahk2 = 0. On
the event of Theorem A.1, we have
(6.2)

0 <

√

m − 1 − C · κ 2 · (w∧ (D ( f , x∗ )) + u) < k Ahk2 = 0,

(A.1)

which is a contradiction. If h ∈ D∧ ( f , x∗ ) \ {0} with k Ahk2 = 0, we just invoke (A.1) for h/k hk2 .
f
Hence, the only feasible point of (Pη ) is x∗ , so that x̂ = x∗ .

To prove Corollary 6.6, we need the following lemma. It shows how to select a scale t > 0 for the
local mean width such that it can be bounded by the conic mean width of a surrogate signal:
Lemma A.2 Let x∗ ∈ Rn and let K ⊂ Rn be convex. For every t > 0 and x̄∗ ∈ K, we have
wt ( K − x∗ ) ≤
In particular, w Rk x̄∗ − x∗ k2 (K − x∗ ) ≤

R +1
R

t+k x̄∗ − x∗ k2
t

· (w∧ (K − x̄∗ ) + 1).

· (w∧ (K − x̄∗ ) + 1) for every R > 0.

Proof. We set t0 := k x̄∗ − x∗ k2 . For every x ∈ K ∩ ( x∗ + tSn−1 ), we have

k x − x̄∗ k2 ≤ k x − x∗ k2 + k x∗ − x̄∗ k2 ≤ t + t0 ,
and therefore, x ∈ Kx̄∗ ,t+t0 := K ∩ ( x̄∗ + (t + t0 ) B2n ). This implies
w t ( K − x ∗ ) = E[

sup
1
h∈ t (K − x∗ )∩Sn−1

h g, hi] = 1t E[

≤ 1t E[ sup h g, x − x̄∗ i] + 1t E[
x∈K x̄∗ ,t+t0

x∈K ∩( x∗ +tSn−1 )

t+t0
t

·

=

t+t0
t

· E[

1
t + t 0 E[

x∈K ∩( x∗ +tSn−1 )

sup h g, x − x̄∗ i]
x∈K x̄∗ ,t+t0

sup
1
h∈ t+t0 (K − x̄∗ )∩ B2n

h g, hi]

h g, x̄∗ − x∗ i]

sup

|

=

h g, x − x∗ i]

sup

{z

=0

}

44

A P ROOFS OF A UXILIARY R ESULTS
x̄∗ ∈K

≤

CSI

≤

t+t0
t

· E[

t+t0
t


· E[(

(∗)

≤

t+t0
t

h g, hi]

sup
h∈cone(K − x̄∗ )∩ B2n

h g, hi)2 ]

sup

1/2

h∈cone(K − x̄∗ )∩ B2n

· (w2∧ (K − x̄∗ ) + 1)1/2 ≤

t+t0
t

· (w∧ (K − x̄∗ ) + 1),

where (∗) follows from [ALMT14, Prop. 3.1(5) and Prop. 10.2].



Proof of Corollary 6.6. We would like to apply Theorem 6.4 with



2η
i
t = max Rk x∗ − x̄∗ k2 , h√
R
+
1

2

m − 1 − C · κ · ( R · (w∧ (D ( f , x̄∗ )) + 1) + u)




+

.




Thus, it suffices verify that the assumption (6.1) is satisfied. For this, we first observe that
wt (D ( f , x∗ )) = wt (D ( f , x∗ ) + x∗ − x∗ )
t≥ Rk x∗ − x̄∗ k2

≤ w Rk x∗ − x̄∗ k2 ((D ( f , x∗ ) + x∗ ) − x∗ )

Lemma A.2

≤
=

R +1
R
R +1
R

· (w∧ ((D ( f , x∗ ) + x∗ ) − x̄∗ ) + 1)
· (w∧ (D ( f , x̄∗ )) + 1),

where we have used the monotony of the local mean width in the second step, and the last step follows
from
D ( f , x∗ ) + x∗ − x̄∗ = { h ∈ Rn | f ( x∗ + h) ≤ f ( x∗ ) } + x∗ − x̄∗
| {z }
= f ( x̄∗ )

h0 :=h+ x∗ − x̄∗

= {h0 ∈ Rn | f (h0 + x̄∗ ) ≤ f ( x̄∗ )} = D ( f , x̄∗ ).

This implies1
t ≥ h√

2η
m − 1 − C · κ 2 · ( RR+1 · (w∧ (D ( f , x̄∗ )) + 1) + u)

i
+

2η
 ,
≥ √
2
m − 1 − C · κ · (wt (D ( f , x∗ )) + u) +
so that the claim again follows from Theorem 6.4.



A.2 Proof of Lemma 6.12
The main idea behind the proof of Lemma 6.12 is to regard the “clipped covariance”
E[clip(hv1 , g i; β 1 ) · clip(hv2 , g i; β 2 )]
as a function in the correlation between hv1 , g i and hv2 , g i, which is given by E[hv1 , g i · hv2 , g i] =
hv1 , v2 i. In fact, we will even see that it is monotone and (6.13) holds with equality if |hv1 , v2 i| = 1. To
make this argument precise, we need the following
1 This

argument does actually not cover the case of η = 0, but here we can use that t = Rk x∗ − x̄∗ k2 > 0 if x∗ 6= x̄∗ .

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

45

Lemma A.3 Let β 1 , β 2 ≥ 0 and define the function
Θ : [−1, 1] → R, ρ 7→ E

h




( g1 ,g2 )∼N 0,

1ρ
ρ1



i
clip( g1 ; β 1 ) · clip( g2 ; β 2 ) .

Then, Θ is continuous, Θ(0) = 0, and Θ is convex on [0, 1]. In particular, Θ(ρ) ≤ ρ · Θ(1) for all 0 ≤ ρ ≤ 1.
The proof of Lemma A.3 is a consequence Price’s Theorem [McM64; Pri58], which is a powerful
tool to calculate the derivative of non-linearly distorted, correlated Gaussian random variables. The
following version of Price’s Theorem was established by Voigtlaender [Voi17] and even applies to
tempered distributions:1
Theorem A.4 (Price, [Voi17, Cor. 2]) For ρ ∈ (−1, 1), let
φρ : R2 → (0, ∞), ( x1 , x2 ) 7→

2π

√1

1− ρ2


· exp −

x12 + x22 −2ρ· x1 · x2
2(1− ρ2 )



 h i
1ρ
be the density function of a random vector g ∼ N 0, ρ 1 . If f ∈ S 0 (R2 ) is a tempered distribution, then
the function

F : (−1, 1) → C, ρ 7→ h f , φρ iS 0 ,S

is infinitely often differentiable with
D
F (k) ( ρ ) =

∂2k
∂x1k ∂x2k

f , φρ

E
S 0 ,S

for all ρ ∈ (−1, 1) and k ≥ 0.

 h i
1ρ
In particular, if f ∈ L1,loc (R2 ), we have F (ρ) = E[ f ( g )] for g ∼ N 0, ρ 1 .
Proof of Lemma A.3. We consider the function
f β1 ,β2 : R2 → R, ( x1 , x2 ) 7→ clip( x1 ; β 1 ) · clip( x2 ; β 2 ).
It is not hard to see that f β1 ,β2 is 1-Lipschitz and therefore f β1 ,β2 ∈ S 0 (R2 ). The (weak) derivative of
f β1 ,β2 is given by the set indicator function of the box (− β 1 , β 1 ) × (− β 2 , β 2 ) ⊂ R2 :
∂2
∂x1 ∂x2 f β 1 ,β 2

= χ(− β1 ,β1 )×(− β2 ,β2 ) .

For any φ ∈ S (R2 ), Fubini’s theorem and the fundamental theorem of calculus now yield
E
D 4
E
D
∂
∂2
f
=
χ
,
φ
,
φ
β
,β
2
2
(− β 1 ,β 1 )×(− β 2 ,β 2 ) ∂x ∂x
1 2
∂x ∂x
1

2

S 0 ,S

1

=

Z β2 Z β
1
− β2

− β1

2

S 0 ,S

∂2
∂x1 ∂x2 φ ( x1 , x2 ) dx1 dx2

= φ( β 1 , β 2 ) − φ( β 1 , − β 2 ) − φ(− β 1 , β 2 ) + φ(− β 1 , − β 2 ).
Hence, by Theorem A.4, the map Θ is smooth on (−1, 1) and
Θ00 (ρ) = φρ ( β 1 , β 2 ) − φρ ( β 1 , − β 2 ) − φρ (− β 1 , β 2 ) + φρ (− β 1 , − β 2 )

= 2(φρ ( β 1 , β 2 ) − φρ ( β 1 , − β 2 ))

(A.2)

for ρ ∈ (−1, 1), where we have also used the symmetry of φρ . Now, let us show the continuity of Θ on
p
the interval [−1, 1]. For this, let g := ( g, g0 ) ∼ N (0, I2 ) and set gρ := ρ · g + 1 − ρ · g0 for ρ ∈ [−1, 1].
1 For

a detailed introduction to Schwartz functions and tempered distributions, see [Fol99, Chap. 8+9] for example.

46

A P ROOFS OF A UXILIARY R ESULTS

 h i
1ρ
Then ( g, gρ ) ∼ N 0, ρ 1 , so that

|Θ(ρ) − Θ(ρ0 )| ≤ Eg [| f β1 ,β2 ( g, gρ ) − f β1 ,β2 ( g, gρ0 )|]
f β1 ,β2 is 1-Lipschitz

≤ Eg [| gρ − gρ0 |]
p
p
ρ0 →ρ
≤ |ρ − ρ0 | · E[| g|] + | 1 − ρ − 1 − ρ0 | · E[| g0 |] −−−→ 0.

Moreover, we observe that
Θ(0) = Eg [ f β1 ,β2 ( g, g0 )] = E[clip( g; β 1 )] · E[clip( g0 ; β 1 )] = 0.
{z
} |
{z
}
|
=0

=0

Finally, let us verify the convexity of Θ on [0, 1]. It is sufficient to show Θ00 (ρ) ≥ 0 for ρ ∈ (0, 1). Indeed,
by (A.2), we have

⇔

Θ00 (ρ) = 2(φρ ( β 1 , β 2 ) − φρ ( β 1 , − β 2 )) ≥ 0




β2 + β2 −2ρ· β · β
β2 + β2 +2ρ· β · β
exp − 1 2(21−ρ2 ) 1 2 ≥ exp − 1 2(21−ρ2 ) 1 2

⇔

β21 + β22 − 2ρ · β 1 · β 2 ≤ β21 + β22 + 2ρ · β 1 · β 2

⇔

0 ≤ 4ρ · β 1 · β 2 ,

where the latter inequality is obviously fulfilled for every ρ ∈ (0, 1).



Now, using Lemma A.3, we can easily derive the statement of Lemma 6.12:

Proof of Lemma 6.12. First, we observe that
 h
(hv1 , g i, hv2 , g i) ∼ N 0, hv 1,v
1

hv1 ,v2 i
1
2i

i

.

If hv1 , v2 i ≥ 0, Lemma A.3 implies
E[clip(hv1 , g i; β 1 ) · clip(hv2 , g i; β 2 )] = Θ(hv1 , v2 i) ≤ hv1 , v2 i · Θ(1),
and an elementary calculation shows that
Θ(1) = E[clip( g ; β 1 ) · clip( g; β 2 )]
|{z}
∼N (0,1)



= 1 − erf

βmax 
√
−
2

q

max )2 /2
2 e−( β
π
βmax



· βmin · βmax + erf

βmin 
√
.
2

On the other hand, if hv1 , v2 i ≤ 0, we have
E[clip(hv1 , g i; β 1 ) · clip(hv2 , g i; β 2 )] = −E[clip(h−v1 , g i; β 1 ) · clip(hv2 , g i; β 2 )]

= −Θ(h−v1 , v2 i)
h−v1 ,v2 i≥0

≥ −h−v1 , v2 i · Θ(1)
= hv1 , v2 i · Θ(1),

and by taking the absolute value, the claim follows.



G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

47

A.3 Proof of Lemma 6.14
Before proving the actual statement of Lemma 6.14, we recall the definition of the function h defined
in (2.3):
q
2

2 e−τ /2
π
τ

h : (0, ∞) → (0, ∞), τ 7→

+ erf( √τ ) − 1.
2

The well-definedness of h, i.e., h(τ ) > 0 for all τ > 0, is a consequence of standard bounds on the error
function, e.g., see [FR13, Prop. 7.5]. Since limτ &0 h(τ ) = ∞, limτ →∞ h(τ ) = 0, and
h0 (τ ) = −

q

2
2 e−τ /2
π τ2

<0

for all τ > 0, we can conclude that h−1 : (0, ∞) → (0, ∞) is also well-defined.
Proof of Lemma 6.14. Adapting the notation from Proposition 2.7, let us set S := SSσ , L := LS c , L̄ := LSrc ,
and for the sake of convenience, let us also restate the definition of F:



F (τ, λ) = n + S − L · h(τ ) · τ 2 · λ2 + λ2 · L − 2λ · L̄ · erf √τ , τ, λ ≥ 0.
2

First, we calculate the partial derivatives of F for τ, λ > 0:
q

2
∂F
2
(τ, λ) = 2λ · τ · S − L · h(τ ) + 2λ · π2 e−τ /2 · (λ · L − L̄),
∂τ

∂F
(τ, λ) = 2λ · τ 2 · S − L · h(τ ) + 2 erf( √τ ) · (λ · L − L̄).
2
∂λ
Hence, the gradient of F takes the following form:


q


2 −τ 2 /2
2 · τ 2λ ·
2λ
e
S − L · h(τ )
π


∇ F (τ, λ) =
·
.
λ · L − L̄
2λ · τ 2
2 erf( √τ )
2
|
{z
}
=:M

Since S, L, L̄ > 0 and
q


2
det( M ) = 4λ2 · τ · erf( √τ ) − τ · π2 e−τ /2 > 0
2

for all τ, λ > 0, the only critical point of F on (0, ∞) × (0, ∞) is attained if

0 = S − L · h(τ ) = λ · L − L̄ ⇔ (τ, λ) = (τ 0 , λ0 ) = h−1

S
L



,

L̄
L



.

Computing the value of F at this point, we observe that (cf. (6.17))


( L̄)2
< n = F (τ, 0) = F (0, λ)
F (τ 0 , λ0 ) = n − L · erf √1 · h−1 SL
2

for all λ, τ ≥ 0, implying that F cannot attain a global minimum on the boundary of [0, ∞) × [0, ∞).
Consequently, it suffices to seek for a global minimum on (0, ∞) × (0, ∞). For this purpose, let us first
consider the quadratic function λ 7→ Fτ (λ) := F (τ, λ) for λ > 0 and a fixed τ > 0. It is not hard to see
that Fτ has its unique minimum at
λ0 (τ ) :=

L̄ · erf( √τ )
2

L · erf( √τ ) + τ 2 · (S − L · h(τ ))

.

2

√
Note that λ0 (τ ) > 0, since L · (erf(τ/ 2) − τ 2 · h(τ )) > 0 for all τ > 0. Next, we minimize the

48

B D ETAILS ON I MPLEMENTATION

parameterized mapping
τ 7→ F̃ (τ ) := F (τ, λ0 (τ )) = n −

L̄2 · erf( √τ )2
2

L · erf( √τ ) + τ 2 · (S − L · h(τ ))
2

on (0, ∞). An elementary calculation shows that
0 = F̃ 0 (τ )

 q

2
0 = (S − L · h(τ )) · τ · π2 e−τ /2 − erf( √τ )
2
{z
}
|

⇔

<0

0

⇔

τ=τ =

h−1 SL



,

and F̃ 00 (τ 0 ) > 0. Thus, F̃ attains its unique minimum at τ 0 .
Finally, let (τ, λ) ∈ (0, ∞) × (0, ∞) be arbitrary. By the above observations,
F (τ, λ) ≥ F (τ, λ0 (τ )) ≥ F (τ 0 , λ0 (τ 0 )) = F (τ 0 , λ0 ),
| {z }
=λ0

where equality holds if, and only if, τ = τ 0 and λ = λ0 (τ 0 ) = λ0 . Therefore, F has indeed its unique
global minimum at (τ 0 , λ0 ).


B Details on Implementation
B.1 Greedy Approach to Constructing Sparse Representations
In this section, we present a greedy strategy for finding a good approximation space U ⊂ Rn according
to the argumentation of Subsection 2.4. Recall that our basic geometric idea is to select U ⊂ Rn in such
a way that the orthogonal projection of PU x∗ enjoys a small value of M(Ψ, PU x∗ ) while not being too
distant from x∗ at the same time. This is precisely the purpose of Algorithm B.1 below, where we fix
a target analysis sparsity S in advance and then identify a subspace U with kΨPU x∗ k0 ≤ S and small
approximation error k x∗ − PU x∗ k2 .
Algorithm B.1 (Greedy subspace selection)
Input: Analysis operator Ψ ∈ R N ×n , signal vector x∗ ∈ Rn , target level of sparsity S ∈ [ N ],
(machine) precision parameter ε ≥ 0.
Output: Linear approximation space U ⊂ Rn .
Initialize: S0 := N, S0 := [ N ], U := {0}.

I While S0 > S:
(1) Select an index k ∈ S0 such that kPUk x∗ k2 becomes minimal, where Uk := span{ψk }.
(2) Determine all indices k0 ∈ S0 with |hψk0 , PU ⊥ x∗ i| ≤ ε and exclude them from S0 .
k

(3) Set U := span({ψk | k ∈

S0c })⊥ .

(4) Update S0 := kΨPU x∗ k0 .

I Return U.

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

49

Since k x∗ − PU ⊥ x∗ k2 = kPUk x∗ k2 , the actual intention of Step (1) in Algorithm B.1 is to identify an
k
atom ψk ∈ Rn of Ψ that yields a best approximation of x∗ while enforcing hψk , PU ⊥ x∗ i = 0. The coherk
ence structure of Ψ is then incorporated by Step (2), where all (almost) vanishing analysis coefficients
of PU ⊥ x∗ are excluded from the index set S0 . Let us emphasize that this procedure could be suboptik
mal in the sense of best S-analysis-sparse approximations (cf. (2.8)), but it still may allow for accurate
predictions by means of Theorem 2.8, as demonstrated in Subsection 3.1.4. Finally, let us mention that
the statement of Theorem 2.8 actually considers a scaled version of PU x∗ in (2.9), which is merely an
artifact of its proof and therefore ignored here.

B.2 Numerical Estimation of the Sample Complexity
To validate the quality of our sampling-rate predictions via M(Ψ, x∗ ), we have also reported the true
∗
sample complexity δ(Ψ, x∗ ) of (BPΨ
η =0 ) in Subsection 3.1 for several choices of Ψ and x . Technically,
this actually relies on computing the so-called statistical dimension, which is known to capture the sample complexity of many convex programs [ALMT14].
Let us sketch how this quantity can be numerically calculated. Recalling that D (kΨ(·)k1 , x∗ ) denotes
the descent set of the `1 -analysis objective kΨ(·)k1 at x∗ (see Definition 6.3), the statistical dimension
of the associated descent cone D∧ (kΨ(·)k1 , x∗ ) is given by
δ(Ψ, x∗ ) := E

h

sup

h g, hi

2 i

,

h∈D∧ (kΨ(·)k1 ,x∗ )∩ B2n

where g ∼ N (0, In ); see [ALMT14, Prop. 3.1(5)]. Since D (kΨ(·)k1 , x∗ ) is a polytope containing the
origin, it is not hard to see that there exists t0 > 0, such that
t−1 D (kΨ(·)k1 , x∗ ) ∩ B2n = D∧ (kΨ(·)k1 , x∗ ) ∩ B2n
for all 0 < t ≤ t0 , and in particular,
δ(Ψ, x∗ ) = E

h

sup

h g, hi

2 i

.

(B.1)

h∈t−1 D (kΨ(·)k1 ,x∗ )∩ B2n

The inner expression of this expected value can be rewritten as a simple convex program:
sup h g, hi

h ∈Rn

subject to

kΨ( x∗ + th)k1 ≤ kΨx∗ k1 and khk2 ≤ 1,

(B.2)

where we have used that h ∈ t−1 D (kΨ(·)k1 , x∗ ) is equivalent to kΨ( x∗ + th)k1 ≤ kΨx∗ k1 . Thus, to
approximate the expectation in (B.1), we first select a sufficiently small value of t (typically t = 0.01
works well) and then solve program (B.2) for independent samples g1 , . . . , gk ∼ N (0, In ), using the
Matlab software package cvx. Due to the concentration behavior of empirical Gaussian processes, the
arithmetic mean of k = 200 samples already yields a good estimate of δ(Ψ, x∗ ).
Note that the notion of statistical dimension is essentially equivalent to the conic Gaussian mean
width considered in Section 6. More precisely, we have (see [ALMT14, Prop. 10.2])
w2∧ (D (kΨ(·)k1 , x∗ )) ≤ δ(Ψ, x∗ ) ≤ w2∧ (D (kΨ(·)k1 , x∗ )) + 1.
For numerical evaluations, the statistical dimension however seems to be more appealing due to the
convexity of (B.2).

50

C A N E XCURSION TO A NALYSIS O PERATOR L EARNING

C An Excursion to Analysis Operator Learning
According to Theorem 2.5, the value of M(Ψ, x∗ ) determines the number of measurements required to
ensure recovery. The number of measurements, in turn, can be regarded as the “costs” of successfully
applying (BPΨ
η =0 ). In that sense, the mapping

M : R N ×n × Rn → R, (Ψ, x∗ ) 7→ M(Ψ, x∗ )
actually plays the role of a loss function (or cost function). This perspective enables us to formalize the
procedure of Subsection 5.1 by means of a (machine) learning problem:
Problem C.1 (Operator Learning – Expected Loss Minimization) Let X ⊂ Rn be a collection of feasible
signal vectors and assume that µ is a probability measure on X . Then, solve the minimization problem
min

Ψ∈H

Z
X

M(Ψ, x∗ )dµ( x∗ ),

(C.1)

where H ⊂ R N ×n is called hypothesis set (or parameter set), containing all potential choices of analysis
operators.
Intuitively, the purpose of (C.1) is to select an operator Ψ ∈ H which minimizes the average costs
of reconstructing signals from X . We would like to emphasize that the significance of this approach
clearly depends on both the considered signal class X and the hypothesis set H. For instance, if X is
just a singleton, i.e., X = { x∗ }, solving (C.1) can lead to vacuous minimizers, being strongly adapted
to x∗ . Of more interest are larger classes X whose elements share certain (geometric) characteristics.
In fact, we have already studied a prototypical example in Section 3, namely discrete piecewise constant functions; more formally, X might consist of all those vectors x∗ ∈ Rn with a fixed number of
discontinuities (that is, kΨTV-1 x∗ k0 ≤ STV for some STV > 0).
The choice of H is typically driven by the belief in what properties of Ψ are appropriate to capture the
structural features of X . Once again, Subsection 3.1.2 provides an interesting example: In the scenario
of Haar wavelet systems, H could contain all those Ψ ∈ R N ×n arising from a scale-wise weighting of
Ψirdwt . Then, a minimizer of (C.1) would yield an optimally weighted version of Ψirdwt (with respect
the class of piecewise constant signals). It is also worth mentioning that the dimension of the parameter
space R N ×n is considerably reduced in that way, since the degrees of freedom of H would just equal
the number of scales.
The optimization task of (C.1) is however often unfeasible in practice because the underlying probability measure µ is usually unknown, e.g., when the x∗ are supposed to be natural images. Instead,
one has merely access to a collection of independent samples drawn from a random vector in Rn distributed according to µ. This leads to an empirical version of Problem C.1:
Problem C.2 (Operator Learning – Empirical Loss Minimization) Let X ⊂ Rn be the set of feasible signals. Suppose that x1∗ , . . . , x∗L ∈ X are independent samples drawn according to a probability measure
µ on X . Then, solve the minimization problem
1
Ψ∈H L

min

L

∑ M(Ψ, x∗l ),

(C.2)

l =1

where H ⊂ R N ×n is again a hypothesis set.
The idea behind (C.2) is quite simple: By computing the empirical mean instead of an integral, one
may end up with a tractable problem, whose solutions are still very close to minimizers of (C.1). In
fact, solving and analyzing programs of the type (C.2) is one of the key objectives in statistical learning
theory [MRT12; Vap00]. This field of research offers a rich toolbox that allows us to investigate many
related issues, such as regularization, non-convex optimization, iterative algorithms, or sample complexity. We hope that this methodology, combined with our refined sampling-rate bounds, could be

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

51

a starting point of novel design rules and guidelines in analysis-based modeling. On the other hand,
due to the non-convexity and discontinuity of M(Ψ, x∗l ), the specific program of (C.2) is expected to
be (algorithmically) challenging in general. A detailed study would go beyond the scope of this paper
and is therefore deferred to future research. Let us finally point out that similar learning strategies
were recently proposed in the literature as well, but using very different kinds of loss functions, e.g.,
see [HKD13; RB13; SWGK16; YNGD12].

Acknowledgments
The authors thank Claire Boyer, Peter Jung, Jackie Ma, and Pierre Weiss for fruitful discussions. Our special thanks goes to Ali
Hashemi, who referred us to Price’s theorem, and Felix Voigtlaender, who assisted us by working out a fairly general version
of it [Voi17], which is also used in this work. M.G. is supported by the European Commission Project DEDALE (contract no.
665044) within the H2020 Framework Program as well as the Einstein Center for Mathematics Berlin. G.K acknowledges partial
support by the Einstein Foundation Berlin, the DFG Collaborative Research Center TRR 109 “Discretization in Geometry and
Dynamics,” and the European Commission Project DEDALE (contract no. 665044) within the H2020 Framework Program, DFG
Grant KU 1446/18 as well as DFG-SPP 1798 Grants KU 1446/21 and KU 1446/23. M.M. is supported by the DFG Priority
Programme DFG-SPP 1798.

References
[ALMT14]

D. Amelunxen, M. Lotz, M. B. McCoy, and J. A. Tropp. “Living on the edge: phase transitions in convex
programs with random data”. Inf. Inference 3.3 (2014), 224–294.

[B+17]

R. Baraniuk, H. Choi, R. Neelamani, et al. Rice Wavelet Toolbox, Version 3. URL: https://github.com/ricedsp/
rwt. Aug. 2017.

[BD09]

T. Blumensath and M. E. Davies. “Sampling Theorems for Signals From the Union of Finite-Dimensional
Linear Subspaces”. IEEE Trans. Inf. Theory 55.4 (2009), 1872–1882.

[BF13]

E. van den Berg and M. P. Friedlander. Spot – A Linear-Operator Toolbox. URL: http://www.cs.ubc.ca/labs/
scl/spot/index.html. Aug. 2013.

[BPCPE11]

S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein. “Distributed Optimization and Statistical Learning via
the Alternating Direction Method of Multipliers”. Found. Trends Mach. Learn. 3.1 (2011), 1–122.

[CD95]

R. R. Coifman and D. L. Donoho. “Translation-Invariant De-Noising”. Wavelets and Statistics. Ed. by A. Antoniadis and G. Oppenheim. Vol. 103. Lecture Notes in Statistics. Springer New York, 1995, 125–150.

[CDS98]

S. S. Chen, D. L. Donoho, and M. A. Saunders. “Atomic Decomposition by Basis Pursuit”. SIAM J. Sci. Comput.
20.1 (1998), 33–61.

[CENR11]

E. J. Candès, Y. C. Eldar, D. Needell, and P. Randall. “Compressed sensing with coherent and redundant
dictionaries”. Appl. Comput. Harmon. Anal. 31.1 (2011), 59–73.

[Cha04]

A. Chambolle. “An Algorithm for Total Variation Minimization and Applications”. J. Math. Imaging Vision 20.1
(2004), 89–97.

[CK13]

P. G. Casazza and G. Kutyniok, eds. Finite Frames: Theory and Applications. Applied and Numerical Harmonic
Analysis. Birkhäuser, 2013.

[COS10]

J.-F. Cai, S. Osher, and Z. Shen. “Split Bregman Methods and Frame Based Image Restoration”. Multiscale
Model. Simul. 8.2 (2010), 337–369.

[CRPW12]

V. Chandrasekaran, B. Recht, P. A. Parrilo, and A. S. Willsky. “The convex geometry of linear inverse problems”. Found. Comput. Math. 12.6 (2012), 805–849.

[CRT06]

E. J. Candès, J. K. Romberg, and T. Tao. “Stable signal recovery from incomplete and inaccurate measurements”. Comm. Pure Appl. Math. 59.8 (2006), 1207–1223.

[CT05]

E. J. Candès and T. Tao. “Decoding by linear programming”. IEEE Trans. Inf. Theory 51.12 (2005), 4203–4215.

[CX15]

J.-F. Cai and W. Xu. “Guarantees of total variation minimization for signal recovery”. Inf. Inference 4.4 (2015),
328–353.

[DJ94]

D. L. Donoho and J. M. Johnstone. “Ideal spatial adaptation by wavelet shrinkage”. Biometrika 81.3 (1994),
425–455.

[Don06]

D. L. Donoho. “Compressed sensing”. IEEE Trans. Inf. Theory 52.4 (2006), 1289–1306.

[DPSS14]

A. Deniţiu, S. Petra, C. Schnörr, and C. Schnörr. “Phase Transitions and Cosparse Tomographic Recovery of
Compound Solid Bodies from Few Projections”. Fund. Inform. 135.1-2 (2014), 73–102.

52

R EFERENCES

[DV03]

P. L. Dragotti and M. Vetterli. “Wavelet footprints: theory, algorithms, and applications”. IEEE Trans. Signal
Process. 51.5 (2003), 1306–1323.

[EMR07]

M. Elad, P. Milanfar, and R. Rubinstein. “Analysis versus synthesis in signal priors”. Inverse Probl. 23.3 (2007),
947–968.

[FG05]

M. Fornasier and K. Gröchenig. “Intrinsic Localization of Frames”. Constr. Approx. 22.3 (2005), 395–415.

[FMSY12]

M. P. Friedlander, H. Mansour, R. Saab, and O. Yilmaz. “Recovering compressively sampled signals using
partial support information”. IEEE Trans. Inf. Theory 58.2 (2012), 1122–1134.

[Fol99]

G. B. Folland. Real Analysis: Modern Techniques and Their Applications. 2nd Edition. Wiley, 1999.

[Fou14]

S. Foucart. “Stability and robustness of `1 -minimizations with Weibull matrices and redundant dictionaries”.
Linear Algebra Appl. 441 (2014), 4–21.

[Fow05]

J. E. Fowler. “The Redundant Discrete Wavelet Transform and Additive Noise”. IEEE Signal Process. Lett. 12.9
(2005), 629–632.

[FR13]

S. Foucart and H. Rauhut. A Mathematical Introduction to Compressive Sensing. Applied and Numerical Harmonic Analysis. Birkhäuser, 2013.

[GB08]

M. Grant and S. Boyd. “Graph implementations for nonsmooth convex programs”. Recent Advances in Learning
and Control. Ed. by V. Blondel, S. Boyd, and H. Kimura. Vol. 371. Lecture Notes in Control and Information
Sciences. Springer London, 2008, 95–110.

[GB14]

M. Grant and S. Boyd. CVX: Matlab Software for Disciplined Convex Programming, version 2.1. URL: http://
cvxr.com/cvx. Mar. 2014.

[Gen17]

M. Genzel. “High-Dimensional Estimation of Structured Signals From Non-Linear Observations With General
Convex Loss Functions”. IEEE Trans. Inf. Theory 63.3 (2017), 1601–1619.

[GLPU12]

M. Guerquin-Kern, J. Lejeune, K. P. Pruessmann, and M. Unser. “Realistic analytical phantoms for parallel
Magnetic Resonance Imaging”. IEEE Trans. Med. Imaging 31.3 (2012), 626–636.

[GNEGD14]

R. Giryes, S. Nam, M. Elad, R. Gribonval, and M. E. Davies. “Greedy-like algorithms for the cosparse analysis
model”. Linear Algebra Appl. 441 (2014), 22–60.

[Gor88]

Y. Gordon. “On Milman’s inequality and random subspaces which escape through a mesh in Rn ”. Geometric
aspects of functional analysis. Ed. by J. Lindenstrauss and V. D. Milman. Vol. 1317. Lecture Notes in Mathematics.
Springer, 1988, 84–106.

[GPV15]

R. Giryes, Y. Plan, and R. Vershynin. “On the Effective Measure of Dimension in the Analysis Cosparse Model”.
IEEE Trans. Inf. Theory 61.10 (2015), 5745–5753.

[Grö04]

K. Gröchenig. “Localization of Frames, Banach Frames, and the Invertibility of the Frame Operator”. J. Fourier
Anal. Appl. 10.2 (2004), 105–132.

[HKD13]

S. Hawe, M. Kleinsteuber, and K. Diepold. “Analysis Operator Learning and its Application to Image Reconstruction”. IEEE Trans. Image Process. 22.6 (2013), 2138–2150.

[KABG16]

S. Kitić, L. Albera, N. Bertin, and R. Gribonval. “Physics-Driven Inverse Problems Made Tractable With Cosparse
Regularization”. IEEE Trans. Signal Process. 64.2 (2016), 335–348.

[KKS17]

F. Krahmer, C. Kruschel, and M. Sandbichler. “Total Variation Minimization in Compressed Sensing”. arXiv
preprint: 1608.08852. 2017.

[KLR16]

G. Kutyniok, W.-Q. Lim, and R. Reisenhofer. “Shearlab 3D: Faithful digital shearlet transforms based on compactly supported shearlets”. ACM Trans. Math. Softw. 42.1 (2016), Art. 5.

[KNW15]

F. Krahmer, D. Needell, and R. Ward. “Compressive Sensing with Redundant Dictionaries and Structured
Measurements”. SIAM J. Math. Anal. 47.6 (2015), 4606–4629.

[KR15]

M. Kabanava and H. Rauhut. “Analysis `1 -recovery with Frames and Gaussian Measurements”. Acta Appl.
Math. 140.1 (2015), 173–195.

[KRZ15]

M. Kabanva, H. Rauhut, and H. Zhang. “Robust analysis `1 -recovery from Gaussian measurements and total
variation minimization”. Eur. J. Appl. Math. 26.6 (2015), 917–929.

[LDSP08]

M. Lustig, D. L. Donoho, J. M. Santos, and J. M. Pauly. “Compressed Sensing MRI”. IEEE Signal Process. Mag.
25.2 (2008), 72–82.

[LML12]

Y. Liu, T. Mi, and S. Li. “Compressed Sensing With General Frames via Optimal-Dual-Based `1 -Analysis”.
IEEE Trans. Inf. Theory 58.7 (2012), 4201–4214.

[LMPV17]

C. Liaw, A. Mehrabian, Y. Plan, and R. Vershynin. “A Simple Tool for Bounding the Deviation of Random Matrices on Geometric Sets”. Geometric Aspects of Functional Analysis. Ed. by B. Klartag and E. Milman. Vol. 2169.
Lecture Notes in Mathematics. Springer, 2017, 277–299.

[Mal09]

S. Mallat. Cover image A Wavelet Tour of Signal Processing: The Sparse Way. 3rd Edition. Elsevier, 2009.

[McM64]

E. McMahon. “An extension of Price’s theorem (corresp.)” IEEE Trans. Inf. Theory 10.2 (1964), 168–168.

G ENZEL , K UTYNIOK , AND M ÄRZ : `1 -A NALYSIS M INIMIZATION AND G ENERALIZED (C O -)S PARSITY

53

[MPT07]

S. Mendelson, A. Pajor, and N. Tomczak-Jaegermann. “Reconstruction and subgaussian operators in asymptotic geometric analysis”. Geom. Funct. Anal. 17.4 (2007), 1248–1282.

[MR10]

J. L. Mead and R. A. Renaut. “Least squares problems with inequality constraints as quadratic constraints”.
Linear Algebra Appl. 432.8 (2010), 1936–1949.

[MRT12]

M. Mohri, A. Rostamizadeh, and A. Talwalkar. Foundations of Machine Learning. MIT press, 2012.

[NDEG13]

S. Nam, M. E. Davies, M. Elad, and R. Gribonval. “The cosparse analysis model and algorithms”. Appl. Comput.
Harmon. Anal. 34.1 (2013), 30–56.

[NG12]

S. Nam and R. Gribonval. “Physics-driven structured cosparse modeling for source localization”. Proceedings
of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). 2012, 5397–5400.

[NW13]

D. Needell and R. Ward. “Stable Image Reconstruction Using Total Variation Minimization”. SIAM J. Imaging.
Sci. 6.2 (2013), 1035–1058.

[OKH12]

S. Oymak, M. A. Khajehnejad, and B. Hassibi. “Recovery threshold for optimal weight `1 minimization”. 2012
IEEE International Symposium on Information Theory Proceedings (ISIT). 2012, 2032–2036.

[PM11]

G. Plonka and J. Ma. “Curvelet-Wavelet Regularized Split Bregman Iteration for Compressed Sensing”. Int. J.
Wavelets Multiresolution Inf. Process. 09.01 (2011), 79–110.

[Poo17]

C. Poon. “Structure dependent sampling in compressed sensing: Theoretical guarantees for tight frames”.
Appl. Comput. Harmon. Anal. 42.3 (2017), 402–451.

[Pri58]

R. Price. “A useful theorem for nonlinear devices having Gaussian inputs”. IRE Trans. Inf. Theory 4.2 (1958),
69–72.

[RB11]

S. Ravishankar and Y. Bresler. “MR Image Reconstruction From Highly Undersampled k-Space Data by Dictionary Learning”. IEEE Trans. Med. Imaging 30.5 (2011), 1028–1041.

[RB13]

S. Ravishankar and Y. Bresler. “Learning Sparsifying Transforms”. IEEE Trans. Signal Process. 61.5 (2013), 1072–
1089.

[Roc70]

R. T. Rockafellar. Convex Analysis. Princeton University Press, 1970.

[ROF92]

L. Rudin, S. Osher, and E. Fatemi. “Nonlinear total variation based noise removal algorithms”. Physica D
60.1–4 (1992), 259–268.

[RSV08]

H. Rauhut, K. Schnass, and P. Vandergheynst. “Compressed Sensing and Redundant Dictionaries”. IEEE Trans.
Inf. Theory 54.5 (2008), 2210–2219.

[SF09]

I. W. Selesnick and M. A. T. Figueiredo. “Signal restoration with overcomplete wavelet transforms: comparison
of analysis and synthesis priors”. Proceedings of SPIE, Wavelets XIII. Ed. by V. K. Goyal, M. Papadakis, and
D. V. D. Ville. Vol. 7446. 2009.

[SWGK16]

M. Seibert, J. Wörmann, R. Gribonval, and M. Kleinsteuber. “Learning Co-Sparse Analysis Operators With
Separable Structures”. IEEE Trans. Signal Process. 64.1 (2016), 120–130.

[TGP14]

A. M. Tillmann, R. Gribonval, and M. E. Pfetsch. “Projection onto the cosparse set is NP-hard”. Proceedings of
the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). 2014, 7148–7152.

[Tro15]

J. A. Tropp. “Convex Recovery of a Structured Signal from Independent Random Linear Measurements”.
Sampling Theory, a Renaissance. Ed. by G. E. Pfander. Applied and Numerical Harmonic Analysis. Birkhäuser,
2015, 67–101.

[Vap00]

V. Vapnik. The Nature of Statistical Learning Theory. Springer New York, 2000.

[Ver12]

R. Vershynin. “Introduction to the non-asymptotic analysis of random matrices”. Compressed Sensing Theory
and Applications. Ed. by Y. C. Eldar and G. Kutyniok. Cambridge University Press, 2012, 210–268.

[Voi17]

F. Voigtlaender. “A general version of Price’s theorem”. arXiv preprint: 1710.03576. 2017.

[VPDF13]

S. Vaiter, G. Peyre, C. Dossal, and J. Fadili. “Robust Sparse Analysis Regularization”. IEEE Trans. Inf. Theory
59.4 (2013), 2001–2016.

[YNGD12]

M. Yaghoobi, S. Nam, R. Gribonval, and M. E. Davies. “Noise aware analysis operator learning for approximately cosparse signals”. Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). 2012, 5409–5412.

[ZXCL16]

B. Zhang, W. Xu, J.-F. Cai, and L. Lai. “Precise phase transition of total variation minimization”. Proceedings of
the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). 2016, 4518–4522.

