From Blind deconvolution to Blind Super-Resolution through
convex programming
Augustin Cosse

arXiv:1709.09279v1 [] 26 Sep 2017

Courant Institute of Mathematical Sciences and Center for Data Science, NYU, NYC.
Département de Mathématiques et Applications, Ecole Normale Supérieure, Ulm, Paris,
Université Paris Sciences et Lettres.
September 28, 2017

Abstract
This paper discusses the recovery of an unknown signal x ∈ RL through the result of its convolution
with an unknown filter h ∈ RL . This problem, also known as blind deconvolution, has been studied
extensively by the signal processing and applied mathematics communities, leading to a diversity of
proofs and algorithms based on various assumptions on the filter and its input. Sparsity of this filter, or
in contrast, non vanishing of its Fourier transform are instances of such assumptions. The main result
of this paper shows that blind deconvolution can be solved through nuclear norm relaxation in the case
of a fully unknown channel, as soon as this channel is probed through a few N & µ2m K 1/2 input signals
xn = Cn mn , n = 1, . . . , N, that are living in known K-dimensional subspaces Cn of RL . This result
holds with high probability on the genericity of the subspaces Cn as soon as L & K 3/2 and N & K 1/2 up
to log factors. Our proof system relies on the construction of a certificate of optimality for the underlying
convex program. This certificate expands as a Neumann series and is shown to satisfy the conditions for
the recovery of the matrix encoding the unknowns by controlling the terms in this series. We apply specific
concentration bounds to the first two terms in the series, in order to reduce the sample complexities, and
bound the remaining terms through a more general argument. The first term in the series is bounded
through the subexponential Bernstein inequality. The second term is decomposed into two contributions,
each corresponding to a fourth order gaussian chaos. The first contribution, containing the univariate
fourth order monomials in the random vectors defining the subspaces Cn , is bounded through a matrix
version of the Rosenthal-Pinelis inequality. The second contribution, containing the cross terms, is
bounded by using a decoupling argument for U-Statistics. An incidental consequence of the result of
this paper, following from the lack of assumptions on the filter, is that nuclear norm relaxation can be
extended from blind deconvolution to blind super-resolution, as soon as the unknown ideal low pass filter
has a sufficiently large support compared to the ambient dimension L. Numerical experiments supporting
the theory as well as its application to blind super-resolution are provided.

Acknowledgement. AC was supported by the FNRS, FSMP, BAEF and Francqui Foundations. AC thanks MIT Math, Harvard IACS and The University of Chicago, for hosting him
during this work. AC is grateful to Laurent Demanet and Ali Ahmed for interesting discussions
as well as Joel Tropp for pointing out the matrix version of the Rosenthal-Pinelis inequality.

1

Introduction

In Blind deconvolution, a single unknown signal x[n] is convolved with an unknown filter h[n], resulting in
the output signal y[n] = h[n] ∗ x[n]. An additional additive noise e[n] is also sometimes considered. The
applications of this problem range from image processing, including medical and astronomical imaging as
1

well as super-resolution [44], to communication. The principal motivation behind this problem lies in the
cost of high fidelity imaging and transmission devices. When images obtained through a lower quality device
can be deblurred successfully, this is often more interesting than acquiring images through higher quality
devices. Mathematically, the discrete problem reads
x, h ∈ RL

find
subject to

y[`] = (h ∗ x)[`] + e[`],

` = 1, . . . , L

(1)

Where h ∗ x is used to denote the discrete convolution. Problem (1) is ill-posed in the general case [7]. As
an illutration of the ill-posedness of the problem, note that both x and h are defined up to a scaling α. A
common approach to fix this uncertainty will be to assume that both kxk and khk are unitary. Another
improvement can be obtained by increasing the number of probing signals, thus replacing x with the input
matrix X = [x1 , . . . xN ] and by requiring the probing signals xi to ”live” in lower dimensional subspaces.
When multiple input signals are considered, the discrete convolution reads

yn [`] = (h ∗ xn )[`] =

L−1
X

h[`0 ]xn [(` − `0 )mod L].

(2)

`0 =0

Sensitivity to the noise also seems to remain a serious problem in blind deconvolution [11]. Several approaches
have been introduced in order to first make this problem well-posed and then solve it efficiently. It seems
that the oldest references to blind deconvolution go back to the 1970’s [15, 45]. Most approaches from
engineering and signal processing assume a kernel h with small support which is reasonable when dealing
with applications in communications where the filter is usually representing a blur. When dealing with Green
operators such as in acoustics or inverse scattering, the support of the perturbation m1 (x) is not necessarily
small anymore and those assumptions will not hold. Some of the most relevant attempts at solving blind
deconvolution are listed in section 1.2 below. We essentially skim through those approaches. The list is
non-exhaustive. For a more detailed summary we suggest to turn to [11] and references therein.

1.1

An interesting experiment

When considering a single input, ill-posedness of the blind deconvolution problem manifests itself through
an erroneous decomposition of the recovered image into a contribution to the filter and a weaker (partially
recovered) image. Such a wrong recovery is illustrated in Figs. 1 and 2 below. When considering multiple
snapshots that are sufficiently distinct from each other yet sufficiently compressible in some appropriate
domain; such as multiple slices of a same image volume, or multiple snapshots generated from a time
sequence; the ill-posedness of the problem is mitigated and the inacurate splitting of the recovered image
into a perturbation to the filter and some incomplete recovery of the snapshot itself is not permitted anymore.
Indeed such error terms would accumulate in the subsequent slices and affect the other measurements of the
sequence. In other words, if the first image is inaccurately recovered, then the filter is affected as well and
leads to errors in the second frame as well. Such a wrong reconstruction of the second frame will then imply
a modification of the filter to match the measurements which might not match the modification implied by
the erroneous reconstruction of the first frame. The error is thus gradually reduced throughout the sequence
of inputs. In practice, only a few (N & 1) sufficiently distinct snapshots seem to be enough to prevent
incorrect reconstruction.

1.2

Connections with existing work

Numerous approaches have been introduced to tackle the blind deconvolution problem. Without entering
excessively into the details, we review the most important ones below. An important number of results
rely on alternating methods that estimate each of the unknown signals sequentially rather then jointly. A
2

Figure 1: Blind deconvolution performed through the nuclear norm minimization program (9) with a single
frame (i.e n = 1). (Top left) Original image, (top right) point spread function corresponding to the ideal low
pass filter, (bottom left) Recovered image, (bottom right) recovered low pass filter. Note the ghost image
that appears in the recovered filter. This phenomenon is further highlighted in Fig. 2 below. When only
one probing image is considered, the algorithm is not able to separate the low pass filter from the image.
In the recovered image shown at the bottom, some of the finer details can also be seen to have completely
disappeared. This image should be compared to Fig. 17 and 18 below where multiple slices of a same volume
are used.
Recovered Image

Recovered filter

Figure 2: Further illustration of the overlap between the recovered low pass filter and the recovered image
resulting from an insufficient number of samples. When considering a single image, nuclear norm minimization is unable to split the input image from the filter. The use of a few (sufficiently distinct) compressible
images makes the problem better posed. This phenomenon is also highlighted by the phase diagrams given
in Fig. 3 below.

3

common trend among those approaches is to use a statistical description of the input signal x which enables
the use of the Bayesian framework and associated maximum a posteriori (MAP) estimators based on either
joint or marginal probability distributions.
• In [6] Ayers and Dainty compute the Fourier transform of either of the unknown signals before dividing
the output of the convolution by this Fourier transform to get an estimate of the Fourier transform of
the second unknown. The method then proceeds iteratively.
• Introductions to Bayesian methods in blind deconvolution can be found in [11, 34]. In those papers,
Levin et al. compare the joint MAP estimator against the MAP estimator for the filter only. They show
that the former, even when combined with sparsity prior on the input signal gradient will generally favor
recovery of a blurred input instead of a blurring filter thus leading to an inaccurate reconstruction.
When using the MAP estimator for the filter only, however, for a gaussian prior on the input, it
is possible to recover the true filter. They suggest to solve the resulting formulation through an
expectation maximation algorithm. A similar approach is developed in [33] by Lagendijk et al. For
both estimators, penalty terms favoring sparsity in the derivatives are considered. Other use of penalty
functions include [55] (anisotropic diffusion) and [16] (TV regularization).
• In [43], Reeves and Mersereau model the image as an autoregressive (AR) process and the blurr as
a moving average (MA) process thus reducing problem (2) to the recovery of H and A satisfying
y = H(I − A)−1 n1 + n2 where n1 and n2 are noise models and H and A are respectively the filter
and image (AR) parameters. The approach is somewhat equivalent to the Wiener filter. The original
output y but one pixel is used to get an estimate of the filter and input image parameters and this
estimate is then used to compute the misfit at the remaining pixel. The objective is the average over
all pixels of this one pixel difference.
• Another interesting recent line of work uses the bilinearity of the multiplication in Fourier space to
rewrite the blind deconvolution/amplitude calibration problem as a linear problem in the inverse of
the filter h on the one hand and the input signals on the other. Early references along that line
include [8] as well as [26], where the authors provide numerical experiments in the case of inputs
defined as S-sparse combinations of known basis vectors. The authors consider a filter with unknown
phase, amplitude or both. The experiments show that the recovery is improved when the number of
input signals is increased. This line of work extends the early approach from signal processing (see the
aforementioned [6]). Recovery guarantees along that line were recently provided in [39] as well as [17]
where the subspace assumption is removed. Despite its interest in terms of the lack of incoherence
condition on the filter, the principal disadvantage of the least squares approach is its lack of stability
and, as a consequence, its inability to handle functions with vanishing or very small entries.
The approach developed in this paper follows the recent line of work on semidefinite programming relaxations
of non linear and non convex problems in order to derive robustness and unconditional convergence garantees.
Some of the most relevant results along that particular line of work include the following
• In [2] Ahmed et al. certify recovery of both h and x up to a scaling through minimization of the nuclear
norm of the matrix X used as a proxy for the rank one matrix hx∗ . They consider a single input signal
and construct their certificate of optimality using the golfing scheme [27]. The main difference with the
result of this paper is that their result requires both the impulse response and the input signal to ”live”
in lower dimensional subspaces. They certify recovery through the construction of a dual certificate
whenever L & K1 + K2 up to log factors. K1 and K2 denote respectively the dimensions of the filter
and input subspaces. This paper is the last of a line of work by Ahmed and Romberg including [3]
and [5] where the authors study mixing and recovery of signals generated from a few basis elements
and [4] where the authors consider optimal sampling and recovery of a similar ensemble. In each case
the recovery of the matrix encoding the ensemble is carried out through nuclear norm minimization.
• The idea from [2] was further refined in [1] where the authors remove the subspace assumption on the
filter, while maintaining a sparsity constraint together with sufficient incoherence. In this paper, a
4

series of inputs xn has to be recovered from convolutions of the form y = w ∗ xn where w is assumed
to be S-sparse. As in this paper, the inputs are assumed to live in known K-dimensional subspaces.
Recovery of both the inputs xn and the filter w is guaranteed as soon as K + S log2 S . L/ log4 (LN )
and N & log2 (LN ) for inputs xn that are sufficiently distinct from each other. The relation of [2] to
the present paper is further discussed in section 4.
• In [38], Ling et al. study calibration problems of the form y = DAx for unknown (S-sparse) signals
x and diagonal calibration matrix D ∈ RL×L defined as D = diag(Bh) with B ∈ RL×M such that
B ∗ B = I and h ∈ RM . They show that both x and D can be recovered through minimization of the
1-norm of the matrix X used as a proxy for the rank one matrix X0 = h0 xT0 encoding the unknowns
of the original problem. Exact recovery is certified for A ∈ RL×K (L < K) random gaussian whenever
L & SM up to log factors and for A random Fourier for comparable scalings. Corresponding stability
results are provided.
• In [37] the same authors study the recovery of input signals xi and filters hi when measurements yi
are given by the sum of the outputs of the filters. They prove recovery as soon as L & R max(K, µ2h N )
where R is the number of channels, K is a bound on the support of the filters and N is the size of the
inputs.
• In [29], Kech et al. study the general bilinear inverse problem with sparsity or subspace constraints
on both h and x and certify injectivity of the bilinear map for L & S1 + S2 (when considering sparse
vectors) or L & N1 + N2 (when considering signals living in known subspaces). S1 , S2 and N1 , N2
respectively denote the sparsity and subspace dimension of each of the signals involved in the bilinear
product.
Besides convexifying and linearizing approaches, a recent trend from the statistics and optimization communities studies initialization and convergence guarantees of non convex optimization algorithms. The first
results along that trend were obtained in [13]. Similar ideas have been applied to blind deconvolution by
Cambareri et al. [10] as well as Li et al. [35].
The proof of our main result follows the general idea developed in [14] by Candès et al., although, transposing
this idea to the blind deconvolution framework.
Finally, an incidental result following from the lack of assumptions on the filter in the main result of this
paper, is the extension of nuclear norm minimization from blind deconvolution to blind super-resolution.
This idea is further discussed in section 4. The use of convex programming in blind super-resolution was
recently discussed by Yang et al. in [53, 54]. In these papers, the authors derive recovery guarantees for
the atomic norm, in the case of the modulation of an unknown spike train satisfying a minimum separation
condition and modulated by random waveforms generated from a random low dimensional subspace whose
basis vectors satisfy some incoherence condition.
Our paper derives recovery guarantees for a multiple inputs formulation similar to the one in [8] by using
the nuclear norm relaxation framework introduced in [2]. In the Fourier domain, problem (2) can read
as an affine rank minimization problem which can in turn be solved through nuclear norm minimization
(see [23, 24, 41]). Certifying recovery through the nuclear norm minimization program then relies on the
construction of a so-called dual certificate. The main conclusion of this√paper is that the restrictions on the
impulse response h in [2] can be lifted by considering a few N & µ2m K (up to log factors) input signals
rather than a single one. The recovery also requires a sufficient number of measurements (L & K 3/2 µ2h ).
The rest of the paper is organized as follows: Section 1.3 derives the mathematical formulation of the problem
and introduces the underlying notations. Section 1.4 states the main result of the paper. Section 2 outlines
the main argument of the proof. The proof relies on several lemmas. Most of them are only stated in
section 2 and their proofs are detailled in section 3. Numerical experiments supporting the statement of
Theorem 1 are provided in section (4) in which applications to super-resolution and medical imaging are
discussed. Finally, the result of Theorem 1 is commented in section 5 together with a few related open
questions.
5

1.3

Problem Formulation

As explained above, assuming both an unknown impulse response and unknown probing signals is ill-posed
in the general case. For this reason this chapter considers the problem of estimating the channel response h
and the input signals whenever those input signals belong to K-dimensional generic subspaces Cn of RL . In
other words we reduce the number of unknowns from L × N to K × N + L for LN measurements. Intuitively,
it should now be clear that by taking L sufficiently large, the problem should become better posed. Each of
the generic subspace Cn will be represented by a corresponding basis matrix Cn with gaussian i.i.d entries
such that Cn [k, `] ∼ N (0, L1 ), thus with variance σ 2 = 1/L. Consequently, the input signals can thus read
xn = Cn mn , n = 1, . . . , N for some coefficients mn ∈ RK . For notational convenience, we let m denote
the vector concatenating each of the coefficient vectors mn as m = [m1 , . . . , mN ]. We will work in the
Fourier domain since it turns the convolution (2) into a Hadamard product. Let us introduce the DFT
matrix F whose rows are defined as
1
f` [k] = √ e−2πi(k−1)`/L , `, k ∈ {1, . . . , L}.
(3)
L
Note that the DFT matrix F satisfies F F ∗ = F ∗ F = I where F ∗ is used to define the inverse discrete
Fourier tranform. We will use the notation ŷ to denote the Fourier transform of y and y ∗ to denote the
conjugate transpose of y. Let [L] = {1, . . . , L} and [N ] = {1, . . . , N }. We use ĉ`,n to represent the vector
defined from the basis CN by putting the `th coefficient of the Fourier transform of the k th column of Cn
at position (n − 1)K + k for each 1 ≤ k ≤ K and zeros everywhere else. Using the canonical vector e` ,
c`,n can thus read as ĉ`,n [k] = en ⊗ (f`∗ Cn ). Given this definition, and with the subspace decomposition
xn = Cn mn , in the Fourier domain, problem (2) can be written as
√
(`, n) ∈ [L] × [N ],
ŷn [`] = Lĥ[`]x̂n [`],
√
(4)
(`, n) ∈ [L] × [N ],
= L(F h)[`](F Cn mn )[`],
√
= Lhf` , hihĉ`,n , mi,
(`, n) ∈ [L] × [N ].
Now using the Frobenius inner product,
√
(5)
ŷn [`] = Lhf` ĉ∗`,n , hm∗ i,
(`, n) ∈ [L] × [N ].
√
The scaling L comes from equality between the Fourier transform of the convolution y and the product
of the Fourier
√ transforms of x and h. From now on we will include the scaling in the definition of ĉ`,n ,
ĉ`,n := ĉ`,n L. We let the rank one matrices A`,n encode the outer product f` ĉ∗`,n so that the last line
of (5) can be written compactly as ŷn [`] = hA`,n , hm∗ i. Following definition of those matrices, problem (5)
can now be recast as the following affine rank minimization problem,
minimize

rank(X)

subject to

ŷ`,n = hA`,n , Xi.

X

(6a)
(`, n) ∈ [L] × [N ].

(6b)

The matrix X is used as a proxy for the rank one matrix hm∗ encoding the impulse response and the
coefficient vector. If the original problem is well posed, then the unique solution to problem (6) is given by
hm∗ . Affine rank minimization is hard in the general case. A common approach that has revealed efficient
so far when dealing with affine rank minimization is to replace the hard minimal rank constraint by the
minimization of the nuclear norm (see [23,41]). Following this trend from convex optimization, we thus relax
problem (6) into
minimize

kXk∗

subject to

ŷ[`, n] = hA`,n , Xi

X

(7a)
(`, n) ∈ [L] × [N ].

(7b)

We will sometimes write problem (7) compactly by introducing the linear map A : CL×KN 7→ CLN and
defined from the matrices A`,n as
{A(X)}(`,n) ≡ hA`,n , Xi = A`,n (X).
6

(`, n) ∈ [L] × [N ].

(8)

Using definition (8), formulation (7) can now read compactly as
minimize

kXk∗

(9a)

subject to

A(X) = y.

(9b)

X

The next section summarizes the strategy that will be used to certify recovery of the rank one matrix
X0 = hm∗ through the convex relaxation (7). We assume without loss of generality that khk = 1 and
kmk = 1 which implies kX0 kF = khm∗ k = 1. The next section gives the main result of the paper.

1.4

Main Result

The main result of this paper shows that exact recovery through the nuclear norm relaxation (9) occurs
with high probability (on the genericity of the subspaces) as soon as the dimension K of the subspaces,
the ambient dimension L and the number of input signals N obey L & K 3/2 µ2h and N & K 1/2 µ2m up to
log factors. The coherence factors µ2h and µ2m measure the respective spreading of h and m and are defined as
N · kmn k2
L · |ĥ[`]|
,
(10)
µ2m = sup
(11)We
µ2h = sup
kmk2
khk2
n
`
can now state the main result of this paper.
Theorem 1. Let Cn be random matrices of size L × K with gaussian independent and identically distributed
(i.i.d) entries, i.e, Cn [`, k] ∼ N (0, 1/L) representing corresponding K-dimensional subspaces. Let mn denote
K-dimensional vectors representing the L-dimensional input signals xn in the bases Cn so that xn = Cn mn .
The nuclear norm relaxation of (9) recovers the rank one matrix X0 = hm∗ r
where m = [m1 , . . . , mN ],
 2

Kµh
µ2m
from the measurements yn = h ∗ xn , with probability at least 1 − c1 (LN )−β −
L + N , where c1 , β
are known positive constants, as soon as
L & K 3/2 µ2h ,

N & K 1/2 µ2m .

Numerical experiments are provided in section 4 to support the claim of Theorem 1. The proof of this theorem
relies on the construction of a dual certificate and its analysis by means of a Neumann series and appropriate
concentration bounds. The exponents arising in the sample complexities are due to the proof technique used
to control the terms appearing in the expansion defining the certificate. Explicit concentration bounds are
used for the first two terms, leading to sample complexities of K 1/2 and K 3/2 . It is likely that extending
those specific bounds to higher order terms in the Neumann series can further reduce those exponents.

2

Proof of Theorem 1

This section outlines the argument that will be used to certify exact recovery of the matrix hm∗ through
the nuclear norm minimization program (7). It is organized as follows. Section 2.1 derives uniqueness and
optimality conditions for problem (7). Proving optimality and uniqueness of the solution hm∗ for the nuclear
norm minimization program (7) is equivalent to exhibiting a dual vector Y defining a hyperplane separating
the nuclear norm ball and the affine subspace S ≡ {X : A(X) = b}. Proving Theorem 1 can thus be done
by first suggesting a candidate for Y and then, proving that Y is indeed corresponding to a separating
hyperplane (i.e that it satisfies the conditions of section 2.1). Section 2.2 recalls the important concepts
that are needed to bound subexponential and subgaussian random variables. Such variables will appear
extensively in the construction of Y . Section 2.3 introduces a candidate for the normal vector Y and then
explains how this particular ansatz can be proved to satisfy the optimality conditions from section 2.1 and
thus to correspond to a valid separating hyperplane for the underlying nuclear norm and affine subspace.

7

2.1

Optimality and Uniqueness

Proving optimality of the solution X0 = hm∗ is equivalent to exhibiting a (dual) vector Y that would be
simultaneously normal to one of the supporting hyperplanes to the nuclear norm ball kXk∗ ≤ 1 at hm∗ and
to the affine subspace A(X) = y. This is also equivalent to finding a subgradient to the nuclear norm that is
normal to the affine subspace. Such a vector implies that the affine subspace is locally tangent to the nuclear
norm ball (or equivalently to the nuclear norm) at hm∗ and so, that any further reduction in the nuclear
norm would require to leave the affine subspace. Another way to understand the geometric meaning of such
a normal vector is that it defines a separating hyperplane between the nuclear norm ball (resp. nuclear
norm) and the affine subspace. We now derive the mathematical conditions defining the normal vector. The
tangent space T of the nuclear norm at hm∗ can be defined as

T = hx∗ + ym∗ , x ∈ RKN , y ∈ RL
(12)
The projection PT (Y ) = (Y )T onto this tangent space can be defined as (see for example [14])
PT (Y ) = hh∗ Y + Y mm∗ − hh∗ Y mm∗ .

(13)

As a consequence, the corresponding projector onto the orthogonal complement of T , PT⊥ (Y ) is defined
⊥
simply as (Y )⊥
T = PT (Y ) = Y − PT (Y ). From the definition of those two projectors, for a general matrix
X = U ΣV , the subdifferential of the nuclear norm is known to be defined as [42]

∂kXk∗ = Y | PT (Y ) = U V , kPT⊥ (Y )k∞ ≤ 1
(14)
Here kXk∞ denotes the operator norm of the matrix X. From definition (14) and the discussion above,
proving optimality of X0 thus reduces to exhibiting a vector Y in the range of A∗ such that
(Y )T = hm∗ ,

k(Y )T ⊥ k∞ ≤ 1.

(15)

Finding a vector Y that satisfies the conditions (15) above might be hard. In particular, the strict condition
(Y )T = hm∗ might be difficult to satisfy. Some recent results [2, 21, 40] actually show that such an
textitexact dual vector is not always required and that one can instead certify recovery through the construction of a so-called inexact dual certificate by a careful analysis of the properties of the linear map A. Let Z
be such that A(Z) = 0. One can always define two matrices H⊥ and M⊥ such that [h, H⊥ ] and [m, M⊥ ]
are unitary and such that H⊥ and M⊥ satisfy hH⊥ M⊥ , PT⊥ (Z)i = kPT⊥ (Z)k∗ (it suffices to take M⊥ and
H⊥ to be the left and right unitary matrices in the singular value decomposition of PT⊥ (Z) multiplied by
the sign matrix from the eigenvalues. For M = hm∗ , we then have [40],
∗
kM + Zk∗ ≥ hhm∗ + H⊥ M⊥
, M + Zi
∗

≥ kM k∗ + hhm +

∗
H⊥ M⊥

(16)
− Y , Zi

∗

T
≥ kM k∗ + hhm − PT (Y ), PT (Z)i + hH⊥ M⊥
− PT⊥ (Y ), PT⊥ (Z)i

≥ kM k∗ − khm∗ − PT (Y )kF kPT (Z)kF + (1 − kPT⊥ (Y )k)kPT⊥ (Z)k∗
In the first inequality, we use the duality between the nuclear norm and the operator norm,
kAk∗ =

sup hA, Xi.

(17)

kXk∞ ≤1

In the second line, we use the fact that since Z is in the kernel of the linear map A, any dual vector in the
range of A∗ implies hY , Zi = 0. In the last line, we use Hölder’s inequality. The last line of (16) shows that
kM + Zk∗ > kM k∗ as soon as
(1 − kPT⊥ (Y )k)kPT⊥ (Z)k∗ − khm∗ − PT (Y )kF kPT (Z)kF > 0

(18)

In particular, this condition thus implies uniqueness of the solution X0 = hm∗ on top of optimality which
followed from (15).
8

From the condition (18) also follows the notion of inexact dual certificate. The idea of an inexact dual
certificate relies on relaxing the difficult constraint PT (Y ) = hm∗ at the expense of a slight strengthening of
the operator norm constraint kPT⊥ (Y )k∞ ≤ 1. Formally, this requires finding the relation between kPT⊥ (Z)k∗
and kPT (Z)kF . This is where the norm of A will be needed. To relate the norms of PT (Z) and PT⊥ (Z) we
start by using the fact that Z is in the nullspace of A, so that
kA(PT (Z))k − kA(PT⊥ (Z))k ≤ kA(PT (Z) + PT⊥ (Z))k = 0

(19)

which implies kA(PT (Z))k = kA(PT⊥ (Z))k. We then have
kA(PT (Z))k ≤ kAkkPT⊥ (Z)k.

(20)

The proof of lemma 1 below can be used to show that
kPT A∗ APT − PT k2 ≤

1
,
2

w.h.p.

(21)

which can then be used to derive
kAPT (Z)k22 − kPT (Z)k22 ≤

1
kPT (Z)k22
2

(22)

Equation (22) in turn implies
1
kPT (Z)k22 ≤ kAPT (Z)k22 .
2

(23)

Substituting expression (23) into (20), one gets,
1
√ kPT (Z)k2 ≤ kAkkPT⊥ (Z)kF .
2

(24)

Finally, using the equivalence of the norms yields
1
√ kPT (Z)kF ≤ kAkkPT⊥ (Z)k∗
2

(25)

We can now use this last relation together with the result of (16) to obtain a final lower bound on kM + Zk∗
as,
 √

kM + Zk∗ ≥ kM k∗ + − 2khm∗ − PT (Y )kF kAk + 1 − kPT⊥ (Y )k kPT⊥ (Z)k∗
(26)
Certifying recovery of X0 = hm∗ thus reduces to finding a dual vector Y satisfying the following relation,
for all Z ∈ null(A), γ = kAk,


√
1 − kPT⊥ (Y )k − 2γkhm∗ − PT (Y )kF kPT⊥ (Z)k∗ > 0.
(27)
From (25), if PT⊥ (Z) = 0 then PT (Z) vanishes as well, which implies Z = 0 and we can thus assume that
kPT⊥ (Z)k > 0. With this in mind and using (26), we can thus move on and assume recovery as soon as there
exists a Y satisfying


√
1 − kPT⊥ (Y )k − 2γkhm∗ − PT (Y )kF > 0.
(28)
In the rest of the chapter, we will construct a Y ∈ RanA∗ satisfying the following two conditions,
khm∗ − PT (Y )kF ≤ √

1
,
2γc1

kPT⊥ (Y )k∞ < 1 −

1
c1

(29)

for
p some constant c1 . Computing the operator norm of−βA can be done through proposition 1 and gives
K log(LN ) + β log(L) with probability at least 1 − L . The next section introduces some results on
sub-gaussian and sub-exponential random variables as well as the corresponding Bernstein inequality which
will be used extensively to derive the recovery guarantees.
9

2.2

Sub-gaussian and sub-exponential tails

Concentration through the Bernstein inequality requires the terms in the sum to be bounded which is not
possible with the the linear map A given the infinite tail of the gaussian basis matrices Cn . However, the
gaussian distribution is known to be bounded with respect to the probability measure and so is its moments
generating function. For this reason, we will consider the following ensemble of (Orlicz) norms introduced
for a general measure µ and a non decreasing, convex function Ψ(x) such that Ψ(0) = 0 (see [49]) as,


Z
|X|
dµ ≤ 1}.
(30)
kXkΨ = inf{k > 0 |
Ψ
k
q

For the probability measure and the functions Ψq (x) = ex − 1, we have the following concentration result
which will be used extensively throughout the proofs of Theorem 1 to construct the dual certificate (see [30,
31, 46])
Proposition 1 (Bernstein concentration (Orlicz version)). Let Z1 , . . . , ZP be i.i.d. random matrices of size
m × n with E{Zi } = 0. Define


1/2
1/2
P
P

 1 X
X
1
(31)
E (Zi Zi∗ )
,
E (Zi∗ Zi )
σZ = max

 P
P i=1
i=1
Let X = kZk and assume that kXkΨq ≤ U (q) < ∞ for some q ≥ 1. Further let M = m + n. Then there
exists a constant C > 0 such that ∀t > 0, the following bound holds with probability 1 − e−t .
( r
)

1/q
t + log(M ) (q)
U (q)
t + log(M )
kZ1 + . . . + ZP k ≤ CP max σ
,U
log
(32)
P
σ
P
In the rest of this section we will show how to derive bounds on the Orlicz norm of random variables
with sub-gaussian or sub-exponential tails. Those are two important classes of random variables which
will appear in the proof of Theorem 1. We start by introducing a formal characterization of those random
2
variables. Those two classes are respectively bounded in the kXkΨ2 and kXkΨ1 norms for Ψ2 (x) = ex − 1
x
and Ψ1 (x) = e − 1 as we will see. We now define the notion of sub-gaussian and sub-exponential random
variables. Propositions 2 and 3 can be found, for example in [52].
Proposition 2 (Equivalent characterization of Sub-Gaussian random variables ). For a zero mean random
variable, the following properties are equivalent
i) There is a constant σ such that Eeλ(X−µ) ≤ e

λ2 σ 2
2

for all λ ∈ R.

ii) There is a constant c ≥ 1 and a Gaussian variable Z ∼ N (0, τ 2 ) such that
P (|X| ≥ A) ≤ cP (|Z| ≥ A) ,
λX 2

iii) For all λ ∈ [0, 1), Ee 2σ2 ≤ √

for all A ≥ 0.

(33)

1
.
1−λ

The first condition is also called Laplace transform condition. The second one is known as sub-gaussian tail
estimate. A similar ensemble of equivalent definitions can be derived for sub-exponential random variables,
Proposition 3 (Equivalent characterization of Sub-exponential random variables ). For a zero mean random
variable, the following properties are equivalent
i) There are non-negative numbers (ν, b) such that
E{eλ(X−µ) } ≤ e

ν 2 λ2
2

10

for all |λ| <

1
b

(34)

ii) There is a positive number c0 > 0 such that E{eλX } < ∞ for all |λ| ≤ c0
iii) There are constants c1 , c2 > 0 such that
P (|X| ≥ t) ≤ c1 e−c2 t

for all t > 0.

(35)

From propositions 2 and 3 above and in particular from (i), it should be clear that sub-gaussian random
variables are sub-exponential. We now show that subexponential random variables are always bounded
in the Ψ1 norm whether subgaussian variables are always bounded in the Ψ2 norm. For sub-exponential
random variables, the bound on the expectation of the moment generating function (34) provides a direct
bound on the Ψ1 -norm. Indeed note that by using definition (34) together with (30) for Ψ(x) = Ψ1 (x)
(subexponential), we have,
inf{u > 0 : E exp(X/u) ≤ 2} ≤ inf{u > b : exp(ν 2 /2u2 ) ≤ 2}
ν
).
. max(b, √
2 log 2
. max(b, ν).
The subexponential parameters (b, ν) thus provide a bound on the Orlicz 1-norm of subexponential variables.
For sub-gaussian variables, we can use the following result from [50]:
Proposition 4 (Equivalence of the subexponential and subgaussian norms). A random variable X is subgaussian if and only if X 2 is sub-exponential. Moreover,
kXk2Ψ2 ≤ kX 2 kΨ1 ≤ 2kXk2Ψ2 .

(36)

From proposition 4, to derive abound on the Ψ2 norm of a subgaussian variable X, one can thus simply
use the bound on the moment generating function of the corresponding sub-exponential variable X 2 and
then take the square root of this bound. Two other important results which will be assumed throughout the
paper are kX1 X2 kψ1 ≤ kX1 kψ2 kX2 kψ2 whenever X1 and X2 are subgaussians with kX1 kψ2 , kX2 kψ2 < ∞,
and kX − EXkψq ≤ 2kXkψq (see for example [1] and lemma 7 in [5]).
An important class of variables that will be used in the proof of Theorem 1 are chi-squared variables, Z = X 2
with X ∼ N (0, 1). Let us first show that those variables are subexponential. The expectation of the moment
generating function eλX reads,
Z
2
2
2
1
e−X /2σ eλ(X −1) dX
Eeλ(Z−1) = √
2π
Z
2
1
= √ e−λ e−(1/2−λ)X dX
2π
Z
2
e−λ
1
=√ √
e−X dX
π 1 − 2λ
e−λ
1
.
(37)
=√ √
π 1 − 2λ
The last equality (37) holds whenever λ < 1/2. In order to derive a bound of the form (34) we first
require |λ|√< 1/4 to bound the denominator of (37). For this interval, we want a ν satisfying Eeλ(Z−1) =
2
e−λ e− log (1−2λ)π ≤ eλ ν . One can check that this holds for ν = 2. Z is thus sub-exponential with
parameters (2, 4). The
Pproof of Theorem 1 will also repeatedly require us to bound the Orlicz norm of a
sum of the form Z = k Xk2 for independent Xk ∼ N (0, 1). Since the Xk are independent, the expectation
of the moment generating function of Z reads,
Eeλ

Pn

k=1

Xk2

=

n
Y
k=1

2

EeλXk ≤

n
Y

eλ

2 2

ν

k=1

2 2

≤ enλ

ν

,

for all |λ| < 1/b.

(38)

Pn
√
This shows that χ = k=1 Xk2 is subexponential with parameters (b, ν n). The next section introduces a
candidate for the dual vector Y .
11

2.3

Ansatz

In order to satisfy the conditions (29), whenever the map A∗ A concentrates to the identity, the most obvious
choice for Y would be to consider the certificate Y1 defined as,
Y1 = A∗ A(hm∗ ).

(39)

However, this particular construction doesn’t match the sample complexity observed empirically (see section 4.1 and the discussion therein). In particular, at sample complexities observed empirically, it fails to
satisfy the condition PT (Y1 ) ≈ hm∗ . For this reason, we consider instead the better ansatz,
Y = A∗ (APT )(PT A∗ APT )−1 hm∗ ,

(40)

where, as in [14], the notation (PT A∗ APT )−1 hm∗ really means the element F in T obeying (PT A∗ APT )F =
hm∗ . In particular, the candidate certificate Y is well defined as soon as (PT A∗ APT ) is a one-to-one
mapping from T onto T . As soon as the ansatz (40) is well defined, by construction we immediately have
PT (Y ) = hm∗ , and the first condition in (29) is immediately satisfied. Injectivity of this map is the point
of lemma (1) below. This lemma is proved in section (2.4).
Lemma 1. Let T be defined as in (12) with corresponding projector PT defined as in (13). Further let
µ2m and µ2h be defined as in (10) and (11) respectively. Then for any constant δ1 , β1 , as soon as L &
β1 (1/δ1 )Kµ2h log(LN ) and N & β1 (1/δ1 )µ2m log(LN ),
kPT A∗ APT − PT k . δ1

(41)

with probability at least 1 − (L ∨ N )−β1 .
Lemma 1 shows that the mapping PT A∗ APT − PT is a contraction. For a sufficiently small δ1 (see for
example [32]), anf for hm∗ ∈ T , one can thus express the element (PT A∗ APT )−1 hm∗ as
−1

(PT − (PT − PT A∗ APT ))

hm∗ = (PT A∗ APT )

−1

hm∗

∗

(42)

∗

∗

∗

2

∗

= hm + (PT − PT A APT ) hm + (PT − PT A APT ) hm + . . .
(43)
Adding the second factor from (40) gives the following expansion for the ansatz (40)
−1

A∗ APT (PT A∗ APT )

hm∗ = A∗ APT (PT A∗ APT )
∗

∗

−1

hm∗

∗

= A APT hm + A APT (PT − PT A∗ APT ) hm∗ + . . .

(44)
(45)

In order to show the second condition in (15), we are thus left with showing that the projection onto T ⊥ of
each of the terms in (45) can be controlled. Following the approach in [14] we will only use explicit concentration results to bound the first two terms in the series in order to reduce the general sample complexity
and then use a more general argument for the remaining terms. We start with the first term
kPT⊥ (A∗ APT hm∗ ) k = kPT⊥ (A∗ APT hm∗ − EA∗ APT hm∗ ) k,

(46)

noting that EA∗ APT hm∗ ∈ T . This first term can be bounded through lemma 2 below whose proof is given
in section 2.5.
√
Lemma 2. Let ĉ`,n be defined as in section 1.3 (including the scaling L), where c`,n [k] ∼ N (0, 1/L) are
i.i.d. gaussian,. and A`,n = f` ĉ∗`,n . Let the operator A : CL×KN 7→ CL×N be defined from the matrices
A`,n as in (8). The coherences µ2m and µ2h are defined as in (10) and (11). Then for any constants δ2 , β2 ,
as soon as L & β2 (1/δ2 )Kµ2h and N & β2 (1/δ2 )µ2m
kPT⊥ (A∗ APT hm∗ ) k . δ2
−β2

with probability at least 1 − (LN )

.
12

(47)

For the second term in (45), we want to show that the event E3 defined as

E3 ≡ PT⊥ A∗ A (PT A∗ APT − PT ) hm∗ . δ3

(48)

holds with sufficient probability, for a sufficiently small constant δ3 . We prove this result through lemma 3
below. The proof of this lemma is given in section 2.5.
√
Lemma 3. Let ĉ`,n be defined as in section 1.3 (including the scaling L) where c`,n [k] ∼ N (0, 1/L). Let
A`,n = f` ĉ∗`,n and let A denote the linear map defined from the matrices A`,n as in (8). Then for any
constant δ3 , for coherences µ2m and µ2h defined as in (10) and (11),
PT⊥ A∗ A (PT A∗ APT − PT ) hm∗ . δ3
r
Kµ2h
with probability at least 1 − (1/δ3 )
L +

µ2m
N



(49)

− c3 (LN )−β as soon as L & β(1/δ3 )Kµ2h , N & β(1/δ3 )µ2m .

The last lemma below concludes the proof by bounding the remaining terms in the series (45). The proof of
this lemma, which is almost identical to the proof given in [14] is recalled for clarity in section 3.8.
Lemma 4. Under the assumptions of Theorem 1, there exists a constant Ck0 such that
∞
X

k

kPT⊥ A∗ APT (PT − PT A∗ APT ) hm∗ k ≤ Ck0

(50)

k=k0

with probability at least 1 − (LN )−β . In particular, the constant Ck0 can be made smaller than δ as soon as
1
L & (1/δ)K 1+ k0 µ2h , N & (1/δ)K 1/k0 .
As we only need to bound the terms corresponding to k > 1, we can just take k0 = 2 which gives the sample
complexities of Theorem 1.
4

Combining the results of lemmas 1 to 4, and choosing δ such that (1/δ)
i )}i=1 as well as β ≥
√ ≥ {(1/δ
3/2 2
2
{β1 , β2 , β3 , β4 }, we have that as soon r
as L & (1/δ)βK µh , N & (1/δ)β Kµm for any β1 , β2 < 1, with

 2
P
Kµh
µ2m
probability at least 1 − (LN )− βi −
L + N , we have
kPT⊥ (Y )k < δ1 + δ2 + δ3 + Ck0

(51)

which can be made sufficiently smaller than 1 for sufficiently small constants δi , Ck0 . This concludes the
proof of Theorem 1. The remaining sections proceed with the proofs of each of the lemmas mentioned above.

2.4

The injectivity property

In this section, we prove injectivity of the normal operator A∗ A on T . This condition certifies that the
ansatz (40) is well defined. We start by recalling lemma 1 below.
Lemma 1. Let T be defined as in (12) with corresponding projector PT defined as in (13). Further let
µ2m and µ2h be defined as in (10) and (11) respectively. Then for any constant δ1 , β1 , as soon as L &
β1 (1/δ1 )Kµ2h log(LN ) and N & β1 (1/δ1 )µ2m log(LN ),
kPT A∗ APT − PT k . δ1
with probability at least 1 − (L ∨ N )−β1 .

13

(41)

Proof. To bound the operator norm of PT A∗p Ap PT we will use proposition (1). We start by bounding the
variance (31). We then derive a corresponding bound on the Orlicz norm of each of the terms appearing
within the norm (41).
The operator PT A∗ APT : X 7→ PT A∗ APT (X) expands as
PT A∗ APT (X) =

L X
N
X

PT (A`,n )hA`,n , PT (X)i

`=1 n=1

=

L X
N
X

PT (A`,n )hPT (A`,n ), Xi

(52)

`=1 n=1

=

L X
N
X

vec(PT (A`,n )) ⊗ vec(PT (A`,n ))∗ vec(X).

`=1 n=1

Now taking the operator norm, we get,
L X
N
X

PT A∗`,n A`,n PT − E{PT A∗`,n A`,n PT }

`=1 n=1

=

N
L X
X

PT A∗`,n A`,n PT − PT

`=1 n=1

PL PN
where we used E `=1 n=1 A∗`,n A`,n = I. We start by deriving the bound for the variance (31). Let us
use Z`,n to denote the operators defined as

Z`,n ≡ PT A∗`,n A`,n PT − E PT A∗`,n A`,n PT
Recall that the variance is defined as

L X
N
 X
∗
σ 2 = max
EZ`,n
Z`,n

`=1 n=1

= max σ12 , σ22

2

,

L X
N
X

2
∗
EZ`,n Z`,n




(53)



`=1 n=1

(54)

Since the variables Z`,n are symmetric, the two bounds σ12 and σ22 in (54) are exactly the same in this case
and we can thus focus on either of them.
2

σ =

L X
N
X

2

EPT A∗`,n A`,n PT kPT (A`,n )k2F

−

2
EPT A∗`,n A`,n PT

`=1 n=1

=

L X
N
X

2
∗

2

E (vec(PT (A`,n )) ⊗ vec(PT (A`,n ))) kvec(PT (A`,n ))k −

2
EPT A∗`,n A`,n PT

.

(55)

`=1 n=1

The Frobenius norm kPT A`,n k can be bounded from the definition of the projector PT (13) and from the
definition of the matrices A`,n as
kPT (A`,n )k2F = khh∗ f` c∗`,n + f` c∗`,n mm∗ − hh∗ f` c∗`,n mm∗ k2F
= khk2 |ĥ[`]|2 kc`,n k2 + kmk2 |hc`,n , mi|2 kf` k2 + khk2 |ĥ[`]|2 |hm, c`,n i|2
+ 2hf` c∗`,n mm∗ , hh∗ f` c`,n i − 2hhh∗ f` c∗`,n mm∗ , hh∗ f` i
− 2hhh∗ f` c∗`,n mm∗ , f` c`,n mm∗ i
= |ĥ[`]|2 kc`,n k2 + |hc`,n , mi|2 − |ĥ[`]|2 |hm, c`,n i|2
≤ |ĥ[`]|2 kc`,n k2 + |hc`,n , mi|2 .
14

(56)

In (56) we use the fact that khk = kmk = 1 as well as kf` k = 1. Now note that for any projection operator
P, we always have kPAk ≤ kAk and kAPk ≤ kAk. The norm (55) can thus simplify to
2

σ =

L X
N
X

2
2

E (vec(PT (A`,n )) ⊗ vec(PT (A`,n ))) kvec(PT (A`,n ))k −

2
EPT A∗`,n A`,n PT

(57)

`=1 n=1

≤ PT

L X
N
X

2

!
EA∗`,n A`,n kPT (A`,n )k2F

PT

(58)

`=1 n=1

≤

L X
N
X

2

EA∗`,n A`,n kPT (A`,n )k2F

(59)

`=1 n=1

In the second line, we use the positive semidefiniteness of the variance and the fact that for matrices A and
B with A − B  0, kA − Bk ≤ kAk. The operator A∗`,n A`,n can be written in matrix form as
A∗`,n A`,n (X) = (vec(A`,n )vec(A`,n )∗ ) vec(X)

= (en e∗n ) ⊗ c`,n c∗`,n ⊗ f` f`∗ vec(X).

(60)
(61)

Substituting the last line of (56) together with (61) into (59), and using the expression for the moments of
multivariate gaussian random variables, one can write,
L X
N
X

Evec(A`,n ) ⊗ vec(A`,n )kPT (A`,n )k2F ≤

`=1 n=1

L X
N
X



(en e∗n ) ⊗ KIK ⊗ |ĥ[`]|2 f` f`∗

`=1 n=1

+

L X
N
X

2(en e∗n ) ⊗ (mn m∗n ⊗ f` f`∗ )

(62)

`=1 n=1

Now using the definitions of µ2h and µ2m from (10) and (11), and noting that
terms of expression (62) can be upper bounded respectively as

P

`

f` f`∗ = F F ∗ = I, the two

L X
N
X
`=1



Kµ2h
(en e∗n ) ⊗ KIK ⊗ |ĥ[`]|2 f` f`∗ ≤
,
L
n=1

(63)

and
L X
N
X

2(en e∗n ) ⊗ (mn m∗n ⊗ f` f`∗ ) ≤ 2

`=1 n=1

µ2m
.
N

(64)

P
For
bound, we use the fact that k n (en en ) ⊗ An k ≤ supn kAn k for any given matrix An and
P the first
k ` f` f`∗ k = 1. For the second bound we use kmm∗ ⊗ Ik ≤ kmm∗ k. Combining (63) and (64) into (62)
gives
X

Evec(A`,n ) ⊗ vec(A`,n )

∗

kPT A`,n k2F

`,n


≤

µ2
µ2h K
+2 m
L
N


(65)

∗
∗
As explained earlier, the exact same result holds for the Z`,n
Z`,n since Z`,n = Z`,n
. We now derive a bound
on the Orlicz norm of the Z`,n . We use the norms derived from the functions Ψ2 (z) = exp(z 2 ) − 1 and
Ψ1 (z) = exp(x) − 1 respectively for sub-gaussian and sub-exponential random variables. For a gaussian
random vector c`,n , we have shown above that kc`,n k2 follows a χ2 distribution with K degrees of freedom
and is therefore subexponential.

15

For two matrices A, B, with kAk ≤ kBk, we also have kAkΨ1 ≤ kBkΨ1 and so from (56),
kZ`,n kΨ1 ≤ kkvec(PT (A`,n ) ⊗ vec(PT (A`,n )kkΨ1
≤ kkPT (A`,n )k2F kΨ1
≤ k|ĥ[`]|2 kc`,n k2 + |hc`,n , mi|2 kΨ1
≤ |ĥ[`]|2 kkc`,n k2 kΨ1 + k|hc`,n , mi|2 kΨ1

(66)

The first term in (66) is a (subexponential) chi-squared distribution with K degrees of freedom. Following
the discussion in section 2.2, the first term of (66) can be bounded by
|ĥ[`]|2 kkc`,n k2 kΨ1 .

Kµ2h
L

(67)

To bound the second term, we use the fact that hc`,n , mn i is a sum of zero-mean gaussian random variables
(i.e gaussian mixture) for which the variance is simply given by the mixture of the variances. The square
of this mixture is thus a chi-squared. Applying an argument similar to (37), one can show that the subexponential parameters of this chi-squared are given by (b, ν) = (kmn k2 , 2) which gives the following bound,
k|hc`,n , mn i|2 kΨ1 . kmn k2 .

µ2m
.
N

(68)

The bound on the Ψ1 -norm of Z`,n is thus finally given by combining (67) and (68) into
kZ`,n kΨ1 .

Kµ2h
µ2
+ m
L
N

(69)

We can now apply proposition 1 which gives
N
L X
X

PT A∗`,n A`,n PT − PT .

`=1 n=1

(s
max

)
 2


µ2h K
µ2m p
µh K
µ2m
+2
t + log(LKN ),
+2
log(LKN )(t + log(LKN ))
L
N
L
N

with probability 1 − e−t . Taking t & β1 log LN for a constant β1 , and L & (1/δ1 )Kµ2h and N & (1/δ1 )β1 µ2m
up to log factors concludes the proof.

2.5

The size property

In this section, we prove the second condition from (15). To do so, we will bound each of the terms in (45).
We start with the first one. This term is bounded through lemma (2),
√
Lemma 2. Let ĉ`,n be defined as in section 1.3 (including the scaling L), where c`,n [k] ∼ N (0, 1/L) are
i.i.d. gaussian,. and A`,n = f` ĉ∗`,n . Let the operator A : CL×KN 7→ CL×N be defined from the matrices
A`,n as in (8). The coherences µ2m and µ2h are defined as in (10) and (11). Then for any constants δ2 , β2 ,
as soon as L & β2 (1/δ2 )Kµ2h and N & β2 (1/δ2 )µ2m
kPT⊥ (A∗ APT hm∗ ) k . δ2
with probability at least 1 − (LN )

−β2

(47)

.

Proof. The proof of lemma 2 relies, once again on proposition 1. We want to bound the quantity
kPT⊥ (A∗ A)hm∗ k = kPT⊥ (A∗ A(hm∗ ) − EA∗ A(hm∗ )) k
≤ kA∗ A(hm∗ ) − EA∗ A(hm∗ )k
16

(70)
(71)

Where we use the fact that EA∗ A = I. Again we introduce variables Z`,n to denote each of the terms in
the sum (71).
Z`,n ≡ PT⊥ (A∗`,n A`,n )hm∗ − EPT⊥ (A∗`,n A`,n )hm∗

(72)

The norm (71) expands as
L X
N
X

Z`,n =

L X
N
X

A`,n hA`,n , hm∗ i − EA`,n hA`,n , hm∗ i

(73)

`=1 n=1

`=1 n=1

In order to apply proposition 1, we again start by deriving the bound for the variance. The Z`,n are not
hermitian anymore as is shown by A`,n A∗`,n = kc`,n k2 f` f`∗ and A∗`,n A`,n = kf` k2 c`,n c∗`,n and we need to
consider the two variance bounds in (31) separately. For the first bound, we have
L X
N
X

∗
=
EZ`,n Z`,n

L X
N
X

2

EA`,n A∗`,n |hA`,n , hm∗ i|2 − |EA`,n hA`,n , hm∗ i|

`=1 n=1

`=1 n=1

≤

L X
N
X

EA`,n A∗`,n |hA`,n , hm∗ i|2

`=1 n=1

≤

L X
N
X

Ef` f`∗ kc`,n k2 |hc`,n , ĥ[`]m∗n i|2

`=1 n=1

≤K

L X
N
X

f` f`∗ kĥ[`]m∗n k2

(74)

`=1 n=1

Now using the fact that

PL

`=1

f` f`∗ = I as well as kmk = 1, the last line of (74) can be reduced to
L X
N
X

∗
EZ`,n Z`,n
≤K·

`=1 n=1

µ2h
L

(75)

For the second bound, a similar argument gives,
E

L X
N
X

∗
Z`,n
Z`,n = E

L X
N
X

kf` k2 c`,n c∗`,n |hf` c`,n , hm∗ i|2 ,

`=1 n=1

`=1 n=1

≤ E

N
L X
X

kf` k2 c`,n c∗`,n |hc`,n , ĥ[`]m∗n i|2

`=1 n=1

.

L X
N
X

X

kĥ[`]m∗n k2 (en e∗n ) ⊗ IK +

`=1 n=1

ĥ[`]mn ĥ[`]∗ m∗n

(`,n)∈[L]×[N ]

µ2
. m.
N
Combining (75) with (76), we get the following bound on the variance σ,
 2

µh K
µ2m
σ.
+
L
N

(76)

(77)

We now bound the Orlicz 1-norm of each of the variables Z`,n in (73). We start by establishing the distribution of kA`,n hA`,n , hm∗ i − EA`,n hA`,n , hm∗ ik
kA`,n hA`,n , hm∗ ik2 = kf` c∗`,n hc`,n , ĥ[`]mn i − ĥ[`]f` m∗n k2
µ2h ∗
kc hc`,n , ĥ[`]mn i − m∗n k
L `,n
µ2
≤ h kc`,n k2 |hc`,n , mn i|2
L
≤

17

(78)
(79)
(80)

Both kc`,n k2 as well as |hc`,n , ĥ[`]mn i|2 are subexponential (chi-squared) variables for which the Orlicz-1
norm can be bounded by using the discussion in section (2.2) (apply (37) for a general gaussian X with
mean 0 and variance σ 2 or alternatively use√ lemma 7 in [2] together with lemma 2.2.1 in [22]). Then
using proposition 4, note that kkc`,n kkψ2 . K and k|hc`,n , mn i|kψ2 . kmn k. Finally note that for two
subgaussian random variables X and Y with kXkψ2 , kY kψ2 < ∞, we have kXY kψ1 ≤ kXkψ2 kY kψ2 (see
lemma 7 in [5]).
√
µh µm K
(81)
kkc`,n k|hc`,n , mn i|kΨ1 . √
LN

∗

kkA`,n khA`,n , hm ikkΨ1

√
µh µm K
≤ √
.
LN

(82)

Substituting (82) and (77) into proposition 1, we have
L X
N
X

A∗`,n A`,n hm∗ − hm∗ .

`=1 n=1

(s 
max

)
√

µ2h K
µ2m p
µh µm K
log(LKN )(t + log(LKN ))
+
t + log(LKN ), √
L
N
LN

with probability at least 1 − e−t . Taking t = β2 log(LN ) gives the desired result.
We now bound the second term in the series (45). The general idea behind this second bound is summarized
through lemma 3 which we recall below.
√
Lemma 3. Let ĉ`,n be defined as in section 1.3 (including the scaling L) where c`,n [k] ∼ N (0, 1/L). Let
A`,n = f` ĉ∗`,n and let A denote the linear map defined from the matrices A`,n as in (8). Then for any
constant δ3 , for coherences µ2m and µ2h defined as in (10) and (11),
PT⊥ A∗ A (PT A∗ APT − PT ) hm∗ . δ3
r
Kµ2h
with probability at least 1 − (1/δ3 )
L +

µ2m
N



(49)

− c3 (LN )−β as soon as L & β(1/δ3 )Kµ2h , N & β(1/δ3 )µ2m .

The second term in the Neumann series reads as
A∗ APT (PT A∗ APT − PT ) hm∗

(83)

This term is a 4th order gaussian chaos involving products of dependent gaussians. It is not sub-exponential
anymore and we will thus need to turn to a generalization of the Bernstein concentration bound (1). To
derive a tail bound on the event E3 in (48), we proceed as follows. Let E`,n denote the expectation,
E`,n ≡ EPT (A`,n )hPT (A`,n ), hm∗ i
2

= h|ĥ[`]|
so that clearly,

P

`,n

m∗n

2

(84)
∗

2

2

∗

+ f` |ĥ[`]|kmn k m − h|ĥ[`]| kmn k m

(85)

E`,n = hm∗ . We start by decomposing the norm in (48) as a sum of four contributions

18

contributions.
PT⊥ A∗ A (PT A∗ APT (hm∗ ) − PT (hm∗ ))


*
= PT⊥

X

A`,n

A`,n , 

+
X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

`0 ,n0

(`,n)



*
X

−PT⊥ E(`,n)

A`,n

A`,n , E`0 ,n0 

A`,n

A`,n , 

+
X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )=(`,n)

(`,n)



*
X

−PT⊥ E(`,n)

A`,n

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )=(`,n)



*
+ PT⊥

(88)

+
X

A`,n , E`0 ,n0 

(`,n)

X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 



*
≤ PT⊥

(87)

+
X
`0 ,n0

(`,n)

X

(86)

A`,n

+
X

A`,n , 

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)



*
X

−PT⊥ E(`,n)

A`,n

(89)

+
X

A`,n , E`0 ,n0 

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)

PT⊥ A∗ A (PT A∗ APT (hm∗ ) − PT (hm∗ ))


*
≤ PT⊥

X

A`,n

A`,n , 

+
X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i

(`0 ,n0 )=(`,n)

(`,n)



*
−PT⊥ E(`,n)

X

A`,n



A`,n

A`,n , 

+
X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)



*
−PT⊥ E(`,n)

X

A`,n

A`,n , E`0 ,n0 

X

(92)

+
X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)

+

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i

(`0 ,n0 )=(`,n)

*
+ PT⊥

(91)

+
X

A`,n , E`0 ,n0 

(`,n)

X

(90)

A`,n hA`,n , E`,n i − E

`,n

X

A`,n hA`,n , E`,n i

(93)

`,n

Now let Cov(A∗ A, PT A∗ APT hm∗ ) be defined as
Cov(A∗ A, PT A∗ APT hm∗ ) ≡ PT⊥ E`,n

X

A`,n hA`,n , (PT (A`,n )hPT (A`,n ), hm∗ i)i

(`,n)

*
− PT⊥ E(`,n)

X
(`,n)

A`,n



+
X

A`,n , 

(`0 ,n0 )=(`,n)

19

E`0 ,n0 PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i

(94)

From (94), we can expand (88) into
≤ PT⊥

X

A`,n hA`,n , (PT (A`,n )hPT (A`,n ), hm∗ i)i

(`,n)

−PT⊥ E(`,n)

(95)
X

A`,n hA`,n , (PT (A`,n )hPT (A`,n ), hm∗ i)i

(`,n)

+ Cov(A∗ A, PT A∗ APT hm∗ )


*
+ PT⊥

X

A`,n

+
X

A`,n , 
*

X

A`,n



X

+
X

A`,n , 

(97)

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)

+

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)

−PT⊥ E(`,n)

(96)

A`,n hA`,n , E`,n i − E

`,n

X

A`,n hA`,n , E`,n i

(98)

`,n

We will successively bound each of terms (95), (96), (97) and (98) through corresponding lemmas 5, (6),
(7), and (8) below. We start with the last term. This term is a sum of sub-exponential random variables
and can be bounded through proposition 1. This idea is summarized by lemma 5 below which is proved in
section 3.1,
√
Lemma 5. Let ĉ`,n be defined as in section 1.3 (including the scaling L) with c`,n [k] ∼ N (0, 1/L) i.i.d.
gaussian and A`,n = f` ĉ∗`,n . Let E`,n be defined as in (85). Then as soon as L & (1/δ5 )β5 Kµ2h , N &
(1/δ5 )β5 µ2m for coherences µ2m and µ2h defined as in (10) and (11),
X

A`,n hA`,n , E`,n i − E

`,n

X

A`,n hA`,n , E`,n i . δ5

(99)

`,n

with probability at least 1 − (LN )−β5 .
The covariance term (96) is a purely deterministic term. It can be controlled through lemma (6) below which
is proved in section 3.2.
Lemma 6. Let the covariance Cov(A∗ A, PT A∗ APT hm∗ ) be defined as in (94). Then for any constant δ6 ,
as soon as L & (1/δ6 )Kµ2h and N & (1/δ6 )µ2m ,
Cov(A∗ A, PT A∗ APT hm∗ ) . δ6

(100)

The first term (95) is a sum of LN independent random variables, each defined from fourth order monomials
in the gaussian vectors c`,n , 1 ≤ ` ≤ L, 1 ≤ n ≤ N . We will bound this first term through the RosenthalPinelis inequality (see (3.1) in [25] as well as [47], Theorem 1, for the matrix version.) which extends
proposition 1. This inequality bounds the expectation of the operator norm of a sum of independent random
matrices from the knowledge of the variance and a bound on the expectation of the largest operator norm
among those matrices. It is recalled through proposition 5 below,
Proposition 5 (Rosenthal-Pinelis [47]). Let Z1 , . . . , Zn be i.i.d. random matrices of dimension d1 × d2 with
E{Zi } = 0. Let
)
( n
n
X
X
∗
∗
σZ = max
E (Zi Zi ) ,
E (Zi Zi )
,
(101)
i=1

i=1

20

as well as

1/2
B = E max kZi k2
,

and

i

Then the expectation of the norm of

C(d1 , d2 ) = 4 (1 + 2dlog(d1 + d2 )e) .

(102)

Pn

Zi can be bounded as
2
p
C(d1 , d2 )σ + C(d1 , d2 )B .
EkZk2 ≤
i=1

(103)

Let Q denote the norm
Q ≡ PT⊥

X

A`,n hA`,n , (PT (A`,n )hPT (A`,n ), hm∗ i)i

(104)

X

(105)

(`,n)

−PT⊥ E(`,n)

A`,n hA`,n , (PT (A`,n )hPT (A`,n ), hm∗ i)i

(`,n)

As soon as one can control the expectation of the operator norm, one can then use standard concentration
tools such as Markov’s inequality to derive a tail bound on the event E5 defined as E5 ≡ {Q ≥ δ} as
Pr (Q ≥ t) ≤

EQ
t

(106)

Lemma 7 below precisely derives such a bound on the expectation EQ. This lemma is proved in section (3.3).
√
Lemma 7 (Fourth order dependence). Let ĉ`,n be defined as in section 1.3 (including the scaling L), where
Cn [`, k] ∼ N (0, 1/L) are i.i.d. gaussian. Let A`,n = f` c∗`,n and let the operators A`,n : CL×KN 7→ CL×N
be defined from those matrices A`,n as in (8). The coherences µ2m and µ2h are defined as in (10) and (11).
Then,
E

L X
N
X

A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗

−

EA∗`,n A`,n (PT A∗`,n A`,n PT )hm∗


. max

`,=1 n=1

µ2L µ2m
,
L N


(107)

Finally, the third term (97) is of the form
X

h(vec(A`,n ) ⊗ vec(A`,n )∗ , vec(A`0 ,n0 ) ⊗ vec(A`0 ,n0 )∗ )

(108)

(`,n)6=(`0 ,n0 )

for appropriate measurable functions h(X, Y ). This is a special case of a U -statistics (see [19], chapter 3.
as well as [28]). We bound this last term by using a decoupling argument from de la Penã et al. [20]. This
result is summarized by the following lemma which is proved in section 3.4 below.
√
Lemma 8. Let ĉ`,n be defined as in section 1.3 (including the scaling L) where c`,n [k] ∼ N (0, 1/L) are
i.i.d gaussian. Let A`,n = f` c∗`,n . And A`,n be defined from the matrices A`,n as the corresponding operators
A`,n (X) = hA`,n , Xi. The coherences µ2m and µ2h are defined as in (10) and (11). Then for any constant
δ8 , β8 , as soon as L & β8 (1/δ8 )Kµ2h and N & (1/δ8 )β8 µ2m

+
*
X
X
PT⊥
A`,n A`,n , 
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 
(109)
(`0 ,n0 )6=(`,n)

(`,n)

*
−PT⊥ E

X

A`,n



A`,n , 

(`,n)

+
X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

with probability at least 1 − (LN )−β8
21

. δ8

(110)

Combining the results of lemmas 5 to 8 with the discussion above, with (1/δ3 ) & max{(1/δi )}8i=4 , β ≥
max{βi }, we have that as soon as L & β(1/δ3 )β3 Kµ2h and N & β3 (1/δ3 )µ2m , the bound of lemma 3 is
satisfied with probability at least 1 − c(LN )−β3 where c is a constant. The next section proceeds with the
proofs of each lemma.

3
3.1

Proofs of Auxiliary lemmas
Proof of lemma 5

√
Lemma 5. Let ĉ`,n be defined as in section 1.3 (including the scaling L) with c`,n [k] ∼ N (0, 1/L) i.i.d.
gaussian and A`,n = f` ĉ∗`,n . Let E`,n be defined as in (85). Then as soon as L & (1/δ5 )β5 Kµ2h , N &
(1/δ5 )β5 µ2m for coherences µ2m and µ2h defined as in (10) and (11),
X

A`,n hA`,n , E`,n i − E

`,n

X

A`,n hA`,n , E`,n i . δ5

(99)

`,n

with probability at least 1 − (LN )−β5 .
The norm on the LHS of (99) is the norm of a sum of subexponential random variables and we can thus
use proposition 1. We start by deriving the bound on the variance. Recall that E`,n is defined as E`,n =
EPT (A`,n )hPT (A`,n ), hm∗ i. Let Z`,n be defined as
X
X
Z`,n ≡
A`,n hA`,n , E`,n i − E
A`,n hA`,n , E`,n i
(111)
`,n

`,n

Using (85), |hA`,n , E`,n i|2 expands as
|hA`,n , E`,n i|2 = ĥ[`]3 hmn , c`,n i + kf` k2 ĥ[`]kmn k2 hmn , c`,n i − ĥ[`]3 kmn k2 hmn , c`,n i

2

. |ĥ[`]|6 |hmn , c`,n i|2 + |ĥ[`]|2 kmn k4 |hmn , c`,n i|2 + |ĥ[`]|6 kmn k4 |hmn , c`,n i|2
Following the proofs of lemmas 1 and 2, each of the variance bounds in (31) can be expressed as


 X

X
∗
∗
σ ≤ max
E
Z`,n
Z`,n , E
Z`,n Z`,n


`,n
`,n


 X

X
≤ max
E
f` f`∗ kc`,n k2 |hA`,n , E`,n i|2 , E
c`,n c∗`,n kf` k2 |hA`,n , E`,n i|2


`,n

(112)
(113)

(114)

(115)

`,n

The first term can be bounded as
E

X
`,n

f` f`∗ kc`,n k2 |hA`,n , E`,n i|2 . sup E
`

X

f` f`∗ kc`,n k2 |hA`,n , E`,n i|2

µ6h
µ2
µ6
kmn k2 + K h kmn k6 + K h3 kmn k6
3
L 
L
L
µ2h µ4h µ2m
µ6m
µ4h µ6m
.K
+ 3+ 2 3
L L2 N
N
L N
.K

22

(116)

n

(117)
(118)

For the second term, we can similarly write,
E

X

c`,n c∗`,n kf` k2 |hA`,n , E`,n i|2

.

X  µ4

h
|ĥ[`]|2
L2

`,n

`,n


µ4h
2
4
+ |ĥ[`]| kmn k + 2 |ĥ[`]| kmn k mn m∗n
L
2

4


4
X
µ4h
2 µh
2
4
2
4
|ĥ[`]| 2 + |ĥ[`]| kmn k + 2 |ĥ[`]| kmn k kmn k2 I
+
L
L

(119)

(120)

`,n

Equations (119) and (120) can further be reduced to

E

X


µ4h
µ4h
4
4
. sup
+ kmn k + 2 kmn k kmn k2
L2
L
n

 4
µ4h
µh
4
4
+
km
k
+
km
k
kmn k2
+
n
n
L2
L2
µ4
µ4
. h2 + m2
L
N


c`,n c∗`,n kf` k2 |hA`,n , E`,n i|2

`,n

(121)

(122)
(123)

Finally the Orlicz norm is bounded by noting that the variables are subexponentials and by developing
kZ`,n kΨ1 . Following the same reasoning as in the proof of lemmas 2 and 1, we can write,


kf` c`,n ĥ[`]3 |hc`,n , mn i + |ĥ[`]|kmn k2 hmn , c`,n i − ĥ[`]3 kmn k2 hmn , c`,n i k
(124)

 3
µh µ2m
µh
√
+
(125)
. kkc`,n kkψ2 k|hc`,n , mn i|kψ2
L3/2
L N
 3

√ µm
µh
µh µ2m
√
. K√
+
(126)
N L3/2
L N
√

µ3h µm
µh µ3m
√
√
+
L3/2 N
L N 3/2
√
µh µm K
. √
LN

kkf` c`,n kkΨ1 .





K

(127)
(128)

Using the fact that kc`,n k and |hc`,n , mn i| are sub-gaussians. We can now apply proposition 1 from which
the conclusion follows,

3.2

Proof of lemma 6

The lemma below first shows that the covariance (94) can be made arbitrarily small.
Lemma 6. Let the covariance Cov(A∗ A, PT A∗ APT hm∗ ) be defined as in (94). Then for any constant δ6 ,
as soon as L & (1/δ6 )Kµ2h and N & (1/δ6 )µ2m ,
Cov(A∗ A, PT A∗ APT hm∗ ) . δ6

(100)

Proof. The covariance in (94) expands as
Cov(A∗ A, (PT − PT A∗ APT ) hm∗ ) =PT⊥ E

X

f` c∗`,n hf` c∗`,n , PT (f` c∗`,n )hPT (f` c∗`,n ), hm∗ ii

`,n

− PT⊥

X

Ef` c∗`,n hf` c∗`,n , EPT (f` c∗`,n )hPT (f` c∗`,n ), hm∗ ii

`,n

23

(129)

As explained above, the second term vanishes. Developing the first term in (129), we have,
Cov(A∗ A, (PT A∗ APT ))
L X
N


X
=E
f` c∗`,n |ĥ[`]|2 kc`,n k2 + |hc`,n , mn i|2 − |ĥ[`]|2 |hc`,n , mn i|2 ĥ[`]hc`,n , mn i

(130)

`=1 n=1

=E

L X
N
X



f` c∗`,n |ĥ[`]|2 ĥ[`]kc`,n k2 hc`,n , mn i + |hc`,n , mn i|2 hc`,n , mn iĥ[`]

(131)



f` c∗`,n |ĥ[`]|2 ĥ[`]|hc`,n , mn i|2 hc`,n , mn i

(132)

`=1 n=1

−E

L X
N
X
`=1 n=1

=

L X
N
X

K|ĥ[`]|2 ĥ[`]f` m∗n + Ef` c∗`,n |hc`,n , mn i|2 hc`,n , mn iĥ[`]

(133)

`=1 n=1

−E

L X
N
X

f` c∗`,n |ĥ[`]|2 ĥ[`]|hc`,n , mn i|2 hc`,n , mn i

(134)

`=1 n=1

In the last line we use E(f` c∗`,n )kc`,n k2 hc`,n , mn i = Kf` m∗n . When considering the remaining terms,
|hc`,n , mn i|2 hc`,n , mn i, those terms can be decomposed as
|hc`,n , mn i|2 hc`,n , mn i =

K
X

(c`,n [k])3 m3n [k] +

3c2`,n [k]m2n [k]

X

K
X

c`,n [j]mn [j]

j=1

k=1

k=1

+6

K
X

(135)

c`,n [i]c`,n [j]c`,n [k]mn [i]mn [j]mn [k].

i6=j6=k

When multiplying (135) by f` c∗`,n , the expectation of the sum reduces to


K
K
K
X
X
X
c`,n [j]mn [j] (136)
3c2`,n [k]m2n [k]
Ef` c∗`,n |hc`,n , mn i|2 hc`,n , mn i = Ef` c`,n  (c`,n [k])3 m3n [k] +
+ Ef` c∗`,n 6

j=1

k=1

k=1

X

c`,n [i]c`,n [j]c`,n [k]mn [i]cn [j]cn [k]

(137)

i6=j6=k

=

K
X

f` mn [k]3 (en ⊗ ek )∗ + 3kmn k2 f` m∗n

(138)

k=1

Plugging (138) into (134), gives
Cov(A∗ A, (PT A∗ APT )) =

L X
N
X

K|ĥ[`]|2 ĥ[`]f` m∗n

(139)

`=1 n=1

+

K
L X
N
X
X
`=1 n=1

!
∗

3

f` mn [k] (e`,n ⊗ ek ) +

f` m∗n kmn k2




ĥ[`] + |ĥ[`]|2 ĥ[`] (140)

k=1

(141)
Each of the two terms in (139) and (140) have operator norms respectively bounded as
L X
N
X

2

K|ĥ[`]|

ĥ[`]f` m∗n

`=1 n=1

L
µ2h X
.K
f` ĥ[`] kmk
L

(142)

`=1

.

24

Kµ2h
L

(143)

P

where we use

`

ĥ[`]f` = h as well as khk = kmk = 1.

L X
N
K
X
X
`=1 n=1

!
∗

3

f` mn [k] (en ⊗ ek ) +

f` m∗n kmn k2

k=1



µ2
ĥ[`] + |ĥ[`]|2 ĥ[`] . m
N

(144)

PL
P P
P P
3
∗
where we use `=1 ĥ[`]f` = h and
≤ supn,k |mn [k]|2 k n k mn [k](e`,n ⊗ ek )∗ k
n
k mn [k] (en ⊗ ek )
as kmn k ≤ 1. Finally, (140) can thus be bounded as


X X
X
µ2 µ2

f` mn [k]3 (e`,n ⊗ ek )∗ +
(145)
f` m∗n kmn k2  |ĥ[`]|2 ĥ[`] . h m
L N
`,n

k

`,n

Note that each of the bounds (143), (144) and (145) can be made smaller than δ for any constant δ as soon
as L & (1/δ)µ2h K and N & µ2m (1/δ). This concludes the proof of lemma 6.

3.3

Proof of lemma 7

Before giving the proof of the lemma, we recall the statement for clarity.
√
Lemma 7 (Fourth order dependence). Let ĉ`,n be defined as in section 1.3 (including the scaling L), where
Cn [`, k] ∼ N (0, 1/L) are i.i.d. gaussian. Let A`,n = f` c∗`,n and let the operators A`,n : CL×KN 7→ CL×N
be defined from those matrices A`,n as in (8). The coherences µ2m and µ2h are defined as in (10) and (11).
Then,

E

L X
N
X

A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗

−

EA∗`,n A`,n (PT A∗`,n A`,n PT )hm∗


. max

`,=1 n=1

Proof. Let us first develop the sum
L X
N
X

PL

`=1

PN

n=1

µ2L µ2m
,
L N


(107)

A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗ . We have

A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗

(146)

`=1 n=1

=

L X
N
X



f` c∗`,n hf` c∗`,n , hĥ[`]c∗`,n + f` hc`,n , mn im∗ − hĥ[`]hc`,n , mn im∗ i ĥ[`]hc`,n , mn i

(147)



f` c∗`,n |ĥ[`]|2 kc`,n k2 + |hc`,n , mn i|2 − |ĥ[`]|2 |hc`,n , mn i|2 ĥ[`]hc`,n , mn i

(148)

`=1 n=1

=

L X
N
X
`=1 n=1

To prove (107) through proposition (5), we need to bound the variance and the expected operator norm of
each of the variables within the norm. We first compute the variance. Deriving the bound on the variance is
the point of section 3.3.1 below. Section 3.3.2 then derives a bound on the expectation of any of the variable
operator norm.

3.3.1

Bound on the variance

Squaring the weights of each of the terms in (148), we get
2

|ĥ[`]|2 kc`,n k2 + |hc`,n , mn i|2 − |ĥ[`]|2 |hc`,n , mn i|2 |ĥ[`]|2 |hc`,n , mn i|2

(149)

. |ĥ[`]|6 kc`,n k4 |hc`,n , mn i|2 + |ĥ[`]|2 |hc`,n , mn i|6 + |ĥ[`]|6 |hc`,n , mn i|6 .

(150)

25

Let us, once again, use Z`,n to denote each of the variables
Z`,n ≡ A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗ − EA∗`,n A`,n (PT A∗`,n A`,n PT )hm∗ .

(151)

Recall that the variance bound is defined as
( L N
)
L X
N
XX
X
∗
∗
σ = max E
Z`,n Z`,n , E
Z`,n Z`,n
`=1 n=1

(152)

`=1 n=1

Using (150), the first term in (152) can be expressed as
E

L X
N
X

A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗



A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗

∗

`=1 n=1

=E

L X
N
X



f` f`∗ kc`,n k2 |ĥ[`]|6 kc`,n k4 |hc`,n , mn i|2 + |ĥ[`]|2 |hc`,n , mn i|6 + |ĥ[`]|6 |hc`,n , mn i|6

(153)

`=1 n=1

=E

L X
N
X

f` f`∗ kc`,n k6 |ĥ[`]|6 |hc`,n , mn i|2 + E

`=1 n=1

+E

L X
N
X

L X
N
X

f` f`∗ kc`,n k2 |hc`,n , mn i|6 |ĥ[`]|2

(154)

`=1 n=1

f` f`∗ kc`,n k2 |ĥ[`]|6 |hc`,n , mn i|6

(155)

`=1 n=1

We now bound the operator norm of each of the terms in (154) and (155). For the first term, first note that
kc`,n k6 =

K
X

!3
|c`,n [k]|2

(156)

k=1

=

K
X

|c`,n [k]|6 + 3

X
k6=k0

k=1

X

|c`,n [k]|4 |c`,n [k 0 ]|2 + 6

|c`,n [k1 ]|2 |c`,n [k2 ]|2 |c`,n [k3 ]|2

(157)

k1 6=k2 6=k3

as well as
|hc`,n , mn i|2 =

K
X

|c`,n [k]|2 m2n [k] + 2<e


X


k=1

c`,n [k]c∗`,n [k 0 ]mn [k]mn [k 0 ]




.

(158)

K 3 µ6h
L3

(159)



k,k0

When multiplying (157) by (158) and taking the expectation, we thus get
E

L X
N
X

f` f`∗ kc`,n k6 |hc`,n , mn i|2 |ĥ[`]|6 .

`=1 n=1

L X
N
X

f` f`∗ K 3 kmn k2

`=1 n=1

 µ6h
L3

.

Noting that for X ∼ N (0, σ), E {X n } = σ(n − 1)!! = c, where c is a constant. For the second term, we can
expand the inner product |hc`,n , mn i|6 in a similar way,
Ekf` k2 kc`,n k2 |ĥ[`]|4 |hc`,n , mn i|6 |ĥ[`]|2
≤

µ6
. h3
L
.

(160)

µ6
Ekc`,n k2 h3 |hc`,n , mn i|6
L

µ6h
L3

K
X

(161)

!
|c`,n [k]|2

|hc`,n , mn i|6

(162)

k=1
K
X
k=1

3
!
K
X
X
|c`,n [k]|2 
c`,n [i]c`,n [j]mn [i]mn [j] +
|c`,n [i]|2 mn [i]2 
i,j

i=1

26

(163)

Developing the right factor in the expression above gives
Ekf` k2 kc`,n k2 |ĥ[`]|4 |hc`,n , mn i|6 |ĥ[`]|2
! K
!3
K
X
µ6h X
2
2
2
.E 3
|c`,n [k]|
|c`,n [i]| mn [i]
L
i=1
k=1
3
! K
K
6
X
X
µ
+ h3 E
|c`,n [k]|2 
c`,n [i]c`,n [j]mn [i]mn [j]
L
i,j

(164)
(165)

(166)

k=1

K
X
µ6
+ h3 E
|c`,n [k]|2
L

!

K
X
µ6
+ h3 E
|c`,n [k]|2
L

!

K
X



2
c`,n [i]c`,n [j]mn [i]mn [j]

i,j

k=1

K
X



!
|c`,n [i]|2 mn [i]2

(167)

i=1


c`,n [i]c`,n [j]mn [i]mn [j]

i,j

k=1

K
X

K
X

!2
|c`,n [i]|2 mn [i]2

(168)

i=1

The first term (165) above can read as
K
X

µ6
E h3
L

!

K
X

|c`,n [k]|2

K
X

!
|c`,n [k]|2

K
X



(169)

.

|c`,n [i]|6 mn [i]6 +

i=1

k=1

µ6
+ E h3
L

|c`,n [i]|2 mn [i]2

i=1

k=1

µ6
E h3
L

!3

K
X

K
X


|c`,n [i]|4 c`,n [j]|2 mn [i]4 mn [j]2 

!
|c`,n [k]|2

(170)

i,j



6

X

|c`,n [i1 ]|2 mn [i1 ]2 |c`,n [i2 ]|2 mn [i2 ]2 |c`,n [i3 ]|2 mn [i3 ]2 

(171)

i1 ,i2 ,i3

k=1



K
K
X
X
µ6h X
mn [i]6 +
mn [i]4 mn [j]2 +
mn [i1 ]2 mn [i2 ]2 mn [i3 ]2 
.K 3
L
i=1
i,j
i1 ,i2 ,i3
!


3
µ6 µ6m
µ2m
µ4 µ2
µ6
µ6
. K h3
+ m2 m +
. m3 K h3
3
L
N
N N
N
N
L

(172)

(173)

For (166), we get,
K
X
µ6h
E
|c`,n [k]|2
L3

3
! K
X

c`,n [i]c`,n [j]mn [i]mn [j] .

k=1

µ6
E h3
L

K
X


! K
X
|c`,n [k]|2 
c3`,n [i]c3`,n [j]m3n [i]m3n [j]

µ6
+ E h3
L

K
X

!
|c`,n [k]|2

k=1

X

6

3
Y


c`,n [ik ]c`,n [jk ]mn [ik ]mn [jk ]

(176)

i1 ,j1 ,i2 ,j2 ,i3 ,j3 k=1

k=1
K
X

(175)

i,j

k=1

µ6
+ E h3
L

(174)

i,j

!
|c`,n [k]|2

K
X




|c`,n [i1 ]|2 |c`,n [j1 ]|2 mn [i1 ]2 mn [j2 ]2 c`,n [i2 ]c`,n [j2 ]mn [i2 ]mn [j2 ]

i1 ,j1 ,i2 ,j2

27

(177)

The first and last terms always vanish and (176) contributes to the expectation only for chains of the form
i1 , i1 , i2 , i2 , i3 , i3 for some i1 , i2 , i3 . We can thus write
3
! K
K
X
X
µ6h
E
|c`,n [k]|2 
c`,n [i]c`,n [j]mn [i]mn [j]
(178)
L3
i,j
k=1

!
K
X
µ6 µ6
µ6h X
(179)
|c`,n [k]|2 
mn [i1 ]2 mn [i2 ]2 mn [i3 ]2  . K h3 m3
. 3
L
L N
i ,i ,i
k=1

1

2

3

For the last two terms (167) and (168), we can respectively write,
2
! K
!
K
K
X
X
X
µ6h
2 
2
2

E
|c`,n [k]|
c`,n [i]c`,n [j]mn [i]mn [j]
|c`,n [i]| mn [i]
L3
i,j
i=1
k=1

! K
!
K
K
X
X
X
µ6h
2 
2
2 2
2
2
2

. 3E
|c`,n [i]| |c`,n [j]| mn [i]mn [j]
|c`,n [k]|
|c`,n [i]| mn [i]
L
i,j
i=1
k=1

!
! K
2
K
K
X
X
X Y
µ6h
2
2
2 

c`,n [ik ]c`,n [jk ]mn [ik ]mn [jk ]
|c`,n [i]| mn [i]
|c`,n [k]|
+ 2 3E
L
i=1
i ,j ,i ,j
k=1

.K

1

1

2

2

(180)

(181)

(182)

k=1


µ6h µ2m
µ6 µ6
kmn k4 . K h3 m3
3
L N
L N

(183)

as well as
K
X
µ6h
E
|c`,n [k]|2
L3


! K
!2
K
X
X
2
2

c`,n [i]c`,n [j]mn [i]mn [j]
|c`,n [i]| mn [i]
i,j

k=1

K
X
µ6
|c`,n [k]|2
. h3 E
L


! K
!
K
X
X
4
4


c`,n [i]c`,n [j]mn [i]mn [j]
|c`,n [i]| mn [i]
i,j

k=1

K
X
µ6
+ h3 E
|c`,n [k]|2
L
k=1

!

(184)

i=1

K
X



(185)

i=1



K
X
c`,n [i]c`,n [j]mn [i]mn [j] 
|c`,n [i]|2 mn [i]2 |c`,n [j]|2 mn [j]2 

i,j

(186)

i,j

Both of these lines vanish when taking the expectation because of the second factor. The total bound
on (160) is thus given by
Ekc`,n k2 |ĥ[`]|6 |hc`,n , mn i|6 . K

µ6h µ6m
L3 N 3

(187)

For the second term in (226), one can derive a similar bound by dividing (187) by µ4h /L2 as,
Ekc`,n k2 |ĥ[`]|2 |hc`,n , mn i|6
.

µ2h
L

µ6
K m3
N

(188)
(189)

All three bounds (159), (187) and (189) can be made sufficiently small as soon as L & Kµ2h and N & µ2m .
We now bound the second term in (152). For this term we have to bound the sum,
L X
N
X

∗
Z`,n
Z`,n

`=1 n=1

=

L X
N
X



kf` k2 c`,n c∗`,n |ĥ[`]|6 kc`,n k4 |hc`,n , mn i|2 + |ĥ[`]|2 |hc`,n , mn i|6 + |ĥ[`]|6 |hc`,n , mn i|6

`=1 n=1

28

(190)

The sum (190) is a block diagonal matrix. When taking the operator norm, one can thus only take the
supremum over the n indices as
E

L X
N
X

∗
Z`,n
Z`,n

`=1 n=1

. sup E
n

L
X



kf` k2 c`,n c∗`,n |ĥ[`]|6 kc`,n k4 |hc`,n , mn i|2 + |ĥ[`]|2 |hc`,n , mn i|6 + |ĥ[`]|6 |hc`,n , mn i|6

(191)

`=1

Now, as before, we split the sum into the diagonal contribution and the off-diagonal terms. For the diagonal,
each of the three terms in (191) read as
sup E
n

L
X



kf` k2 diag(|c`,n [k]|2 ) |ĥ[`]|6 kc`,n k4 |hc`,n , mn i|2

(192)



kf` k2 diag(|c`,n [k]|2 ) |ĥ[`]|2 |hc`,n , mn i|6

(193)



kf` k2 diag(|c`,n [k]|2 ) |ĥ[`]|6 |hc`,n , mn i|6

(194)

`=1

sup E
n

L
X
`=1

sup E
n

L
X
`=1

In (192), (193) and (194), we use diag(|c`,n [k]|2 ) to denote the diagonal matrix built from the diagonal of
c`,n c∗`,n . Following the same reasoning as above, we get
E sup
n

L
X
`=1



K 2 µ4h
kmn k2
kf` k2 diag(|c`,n [k]|2 ) |ĥ[`]|6 kc`,n k4 |hc`,n , mn i|2 .
L2

(195)

K 2 µ4h µ2m
.
L2 N
E sup
n

L
X
`=1

E sup
n

L
X
`=1



µ2
kf` k2 diag(|c`,n [k]|2 ) |ĥ[`]|2 |hc`,n , mn i|6 . m
N

(196)



µ4 µ2
kf` k2 diag(|c`,n [k]|2 ) |ĥ[`]|6 |hc`,n , mn i|6 . h2 m
L N

(197)

P
In (195) we use ` |ĥ[`]|2 kf` k2 = khk2 = 1. For (196) and (197) we use the same reasoning as the one used
to derive (187) except that we now have
E sup
n

≤

L
X
`=1

L
X



kf` k2 diag(|c`,n [k]|2 ) |ĥ[`]|2 |hc`,n , mn i|6

(198)

`=1

kf` k2 |ĥ[`]|2 sup E sup diag(|c`,n [k]|2 ) |hc`,n , mn i|6
`,k



(199)

n

Since the inner product |hc`,n , mn i|6 is now multiplied with a single element |c`,n [k]|2 from the diagonal
diag(c`,n [k]2 ), the bounds in (187), (189) also get divided by K. For the off-diagonal contributions in (191),
only the terms exhibiting an even power in the c`,n [k] monomials corresponding to the off-diagonal entry
will remain when taking the expectation. Focusing on those off-diagonal entries, the three contributions

29

of (192), (193) (192) can now read as,
E

N X
L
X



kf` k2 (c`,n c∗`,n − diag(|c`,n [k]|2 )) |ĥ[`]|6 kc`,n k4 |hc`,n , mn i|2

(200)



kf` k2 (c`,n c∗`,n − diag(|c`,n [k]|2 )) |ĥ[`]|2 |hc`,n , mn i|6

(201)



kf` k2 (c`,n c∗`,n − diag(|c`,n [k]|2 )) |ĥ[`]|6 |hc`,n , mn i|6

(202)

n=1 `=1

E

N X
L
X
n=1 `=1

E

N X
L
X
n=1 `=1



In (200), all the even powers generated from the weight |ĥ[`]|6 kc`,n k4 |hc`,n , mn i|2 will vanish. The only
remaining term can thus be bounded as
sup
n

L
X

kf` k2 mn m∗n |ĥ[`]|6 .

`=1

µ4h K 2 µ2m
L2 N

(203)

For the last two terms, we again use (160) to (186). Retaining only (165) to (168), we can thus write
E

N X
L
X

kf` k2 (c`,n c∗`,n − diag(|c`,n [k]|2 ))|hc`,n , mn i|6

(204)

n=1 `=1

=

N X
L
X

kf` k2

L
N X
X

m3n [k]m3n [k 0 ](en ⊗ ek )(en ⊗ ek0 )∗

kf` k2 (mn m∗n )

n=1 `=1

+

N X
L
X

(205)

k,k0

n=1 `=1

+

X

X

mn [k1 ]2 mn [k2 ]2

(206)

k1 ,k2

kf` k2 kmn k4 mn m∗n



(207)

n=1 `=1

+

N X
L
X

kf` k2

n=1 `=1

+

L
N X
X
n=1 `=1

K
X

!
(mn m∗n ) kmn k2

|mn [k]|2

(208)

k=1
2

kf` k

K
X

!2
|mn [k]|

mn m∗n

2

(209)

k=1

For (206), whenever all three factors are different, the expectation reduces to zero. When at least two of
them are the same, such as in (ci c∗j )(cj c∗i )(ck c` ) or in chains of the form (ci c∗j )(cj c∗k )(ck c∗` ). Those chains
will contribute to the corresponding entries of c`,n c∗`,n . In other words, we have the freedom to choose the
first four indices arbitrarily and the last two are fixed. As a consequence, the only part of (176) that will
contribute to the variance bound when multiplying by (c`,n c∗`,n − diag(|c`,n [k]|2 )) and taking the expectation
is of the form
X
c`,n [i1 ]|c`,n [i2 ]|2 |c`,n [i2 ]|2 |c`,n [i3 ]|2 c`,n [i4 ]mn [i1 ]|mn [i2 ]|2 |mn [i2 ]|2 |mn [i3 ]|2 mn [i4 ].
(210)
i1 ,i2 ,i3 ,i4

When multiplying this sum by (c`,n c∗`,n − diag(|c`,n [k]|2 )), each of the (i1 , i4 ) terms contribute to one offdiagonal entry in the matrix, and we thus get
E(c`,n c∗`,n − diag(|c`,n [k]|2 ))

3
X Y
k1 6=k10
k2 6=k20
k3 6=k30



(c`,n [ki ]c`,n [ki0 ]mn [ki ]mn [ki0 ])

X
k2 ,k3

30

(211)

i=1

mn m∗n |mn [k2 ]|2 |mn [k3 ]|2 .

(212)

From this, we get the following bound on the operator norm of (206),
L X
N
X

X

kf` k2 |ĥ[`]|2 E(c`,n c∗`,n − diag(|c`,n [k]|2 ))

3
Y

(c`,n [ki ]c`,n [ki0 ]mn [ki ]mn [ki0 ]) .

k1 6=k10 ,k2 6=k20 i=1
k3 6=k30

`=1 n=1

µ6m
(213)
N3

Similarly, For (208), simply note that the only possibility for (182) to contribute to the off-diagonal entries,
when taking the expectation, is for the second factor in this expression to exhibit a chain of the form
c`,n [i]c`,n [j]c∗`,n [j]c∗`,n [k] = c`,n [i]|c`,n [j]|2 c`,n [k]∗ . In this case, the whole second factor contributes to each
off-diagonal entry and we can write,


K
3
X
X Y
E(c`,n c∗`,n − diag(|c`,n [k]|2 ))kmn k2 
|c`,n [k1 ]|2 mn [k1 ]2
c`,n [kj ]mn [kj ]
(214)
k2 ,k3 j=2

k1 =1

!
= kmn k2

X

|mn [k]|2

mn m∗n

(215)

k

which gives
X

|ĥ[`]|2 kf` k2 kmn k4 mn m∗n .

`,n

µ6m
N3

(216)

Multiplying by |ĥ[`]|2 and taking the operator norm for the terms in (204) to (209), gives
N X
L
X

E

n=1 `=1



µ6
kf` k2 (c`,n c∗`,n − diag(|c`,n [k]|2 )) |ĥ[`]|2 |hc`,n , mn i|6 . m3 .
N

(217)

This first contribution can thus be made sufficiently small as soon as N & µ2m . The same bound applies
to (202) up to multiplication by |ĥ[`]|4 from which we immediately get,

E

N X
L
X
n=1 `=1



µ6 µ4
kf` k2 (c`,n c∗`,n − diag(|c`,n [k]|2 )) |ĥ[`]|6 |hc`,n , mn i|6 . m3 h2 .
N L

(218)

Grouping (218) together with (217) and (203) and adding the diagonal contributions (195) to (197), we get
L X
N
X

∗
Z`,n
Z`,n


.

`=1 n=1

.

µ6m
µ4h K 2 µ2m
µ6m µ4h
+
+
N 3 L2
N3
L2 N

µ2m
µ4 K 2
+ h2
N
L




+

K 2 µ4h µ2m
µ2m
µ4h µ2m
+
+
L2 N
N
L2 N


(219)
(220)

Combining the bound (220) with (159), (187) and (189), we get the final variance bound for the fourth order
terms,


 2
µm
µ4 K 2
µ6
µ4 K 3 µ6h
µ2h µ4m
σ . max
+ h 2 , K 3 h3 + m2
+
K
(221)
N
L
L
N L3
L N2
This bound can be made less than δ by taking L & (1/δ)Kµ2h and N & (1/δ)µ2m . In order to use proposition (5) we still need to bound the expected maximal operator norm among each of the terms in (107). This
is the point of section 3.3.2 below.

31

3.3.2

Expected maximal operator norm

We now derive a bound on the expected largest operator norm (102) for each of the terms in (107). The
operator norm of any of the term in (107) reads as,
kZ`,n k = f` c∗`,n hf` c∗`,n , PT (f` c∗`,n )hPT (f` c∗`,n , hm∗ ii

(222)

− Ef` c∗`,n hf` c∗`,n , PT (f` c∗`,n )hPT (f` c∗`,n , hm∗ ii
The first term can be bounded as follows,
2

E f` c∗`,n hf` c∗`,n , PT (f` c∗`,n )hPT (f` c∗`,n ), hm∗ ii
2

≤ Ekf` k2 kc`,n k2 |ĥ[`]|2 |hc`,n , mn i|2 |ĥ[`]|2 kc`,n k2 + kf` k2 |hmn , c`,n i|2 + |ĥ[`]|2 |hmn , c`,n i|2


. Ekf` k2 kc`,n k2 |ĥ[`]|2 |hc`,n , mn i|2 |ĥ[`]|4 kc`,n k4 + kf` k4 |hmn , c`,n i|4 + |ĥ[`]|4 |hmn , c`,n i|4


. Ekf` k2 |ĥ[`]|6 kc`,n k6 |hmn , c`,n i|2 + kf` k4 |ĥ[`]|2 kc`,n k2 |hmn , c`,n i|6


+ Ekf` k2 kc`,n k2 |ĥ[`]|6 |hmn , c`,n i|6

(223)
(224)
(225)
(226)

Note that the three terms above are in fact very similar to the ones ( (159) (187) and (189)) that appeared
in the derivation of the variance. Using those results, we can bound (226) as


µ4 K 3 µ6h
µ2h µ4m
µ6
µ2
2
E f` c∗`,n hf` c∗`,n , PT (f` c∗`,n )hPT (f` c∗`,n ), hm∗ ii . m K 3 h3 + m2
+
K
(227)
N
L
N L3
L N2
where we again used the fact that
|hf` c∗`,n , PT (f` c∗`,n )i| ≤ |ĥ[`]|2 kc`,n k2 + kf` k2 |hmn , c`,n i|2 + |ĥ[`]|2 |hmn , c`,n i|2

(228)

|hPT (f` c∗`,n ), hm∗ i|2 = |hf` c∗`,n , hm∗ i|2 ≤ |ĥ[`]|2 |hc`,n , mn i|2

(229)

as well as

For the expectation, simply note that for a random variable X, Var(X) = hX − EX, X − EXi = EkXk2 −
kE(X)k2 ≥ 0. Again, (227) can be made less than δ for any constant δ as soon as L & (1/δ)Kµ2h and
N & (1/δ)µ2m .
3.3.3

Conclusion through Rosenthal-Pinelis

We now use the results of section 3.3.1 and 3.3.2 to bound the norm (107). Using the bound on the operator
norm (227) together with the bound on the variance (221), as well as proposition (5), as soon as L & Kµ2h
and N & µ2m , we have that

E

L X
N
X

2

A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗

−

EA∗`,n A`,n (PT A∗`,n A`,n PT )hm∗


.

`,=1 n=1

Kµ2h
µ2
+ m
L
N


(230)

As explained above, we can then use Markov’s inequality,
L X
N
X

A∗`,n A`,n (PT A∗`,n A`,n PT )hm∗ − EA∗`,n A`,n (PT A∗`,n A`,n PT )hm∗ . t

`,=1 n=1

with probability at least 1 −

1
t



Kµ2h
L

+

µ2m
N

1/2

. Taking t = δ gives the desired result.
32

(231)

3.4

Proof of lemma 8

The proof of of lemma 8 relies on an argument developed in the proof of Theorem 3.4.1 in [19], we recall this
argument, as well as its proof, in section 3.4.1. The idea is that ”just” splitting the norm into independent
components is not enough to get a sufficient level of concentration. Fortunately, the sum inside the norm (110)
is a special case of a U -statistics for which it is possible to derive efficient tail bounds through a decoupling
argument. Let us start by recalling the statement of lemma 8.
√
Lemma 8. Let ĉ`,n be defined as in section 1.3 (including the scaling L) where c`,n [k] ∼ N (0, 1/L) are
i.i.d gaussian. Let A`,n = f` c∗`,n . And A`,n be defined from the matrices A`,n as the corresponding operators
A`,n (X) = hA`,n , Xi. The coherences µ2m and µ2h are defined as in (10) and (11). Then for any constant
δ8 , β8 , as soon as L & β8 (1/δ8 )Kµ2h and N & (1/δ8 )β8 µ2m

+
*
X
X
PT⊥
A`,n A`,n , 
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 
(109)
(`0 ,n0 )6=(`,n)

(`,n)

*
−PT⊥ E

X

A`,n



A`,n , 

+
X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

. δ8

(110)

(`0 ,n0 )6=(`,n)

(`,n)

with probability at least 1 − (LN )−β8
Proof. The proof of lemma 8 relies on a 3 steps argument. The first step uses the decoupling inequality due
to de la Penã and Montgomery-Smith [20] which is summarized by proposition 6. We provide an adapted
version of the proof of this proposition in section 3.4.1 for completeness although this proof is essentially the
same as the proof of Theorem 3.4.1 in [19]. Proposition 6 basically shows that the the probability of success
for the event E8 (A`,n ) corresponding to the U -statistics and defined as


+
*

X
X
E8 (A`,n ) ≡
P⊥
A`,n A`,n , 
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 
 T
0
0
(`,n)
(` ,n )6=(`,n)


+
(232)
*

X
X
−PT⊥ E
A`,n A`,n , 
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 
≥δ

0
0
(`,n)

(` ,n )6=(`,n)

can be upper bounded by the probability of success of the decoupled event E8d (A`,n , Ã`,n ) where the Ã`,n
are independent copies of the A`,n and defined as


+
*

X
X
E8d (A`,n , Ã`,n ) ≡
P⊥
A`,n A`,n , 
PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 
 T
0
0
(`,n)
(` ,n )6=(`,n)


+
(233)
*

X
X
−PT⊥ E
A`,n A`,n , 
PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 
≥δ .

0
0
(`,n)

(` ,n )6=(`,n)

Once we know that Pr(E8 ) ≤ Pr(E8d ), we can use the following argument. First note that,

+
*
X
X
PT⊥ E
A`,n A`,n , 
PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0  = 0.
(`,n)

(234)

(`0 ,n0 )6=(`,n)

as the variables PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 are centered,


X
E
PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0  = 0
(`0 ,n0 )6=(`,n)

33

(235)

In particular, we thus have


*
PT⊥

X

A`,n

A`,n , 


*
X

A`,n

(`0 ,n0 )6=(`,n)

A`,n

+
X

A`,n , 

PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)



*
X

−PT⊥ EA`,n

A`,n

PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 

(237)

(`0 ,n0 )6=(`,n)



*
X

+
X

A`,n , 

(`,n)

+ PT⊥ EA`,n

PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 



*
≤ PT⊥

(236)

+
X

A`,n , 

(`,n)

X

PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)

−PT⊥ E

+
X

A`,n

+
X

A`,n , 

PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 

(238)

(`0 ,n0 )6=(`,n)

(`,n)

Where the expectation in (237) is now taken with respect to the outer matrices A`,n only. That is the second
factor in (238) remains a random variable. Moreover, as we have

+
*
X
X
⊥
∗
PT EA`,n
A`,n A`,n , 
PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm i − E`0 ,n0 
= 0,
(239)
(`0 ,n0 )

(`,n)

one can further bound the norm (238) as

*
X
PT⊥ EA`,n
A`,n A`,n , 


*
X

A`,n

A`,n , 

X

(240)

+
X

PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 

(241)

(`0 ,n0 )

(`,n)

− PT⊥ EA`,n

PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0 

(`0 ,n0 )6=(`,n)

(`,n)

= PT⊥ EA`,n

+
X

D

E
A`,n A`,n , PT (Ã`,n )hPT (Ã`,n ), hm∗ i − E`,n

(242)

(`,n)

≤ PT⊥ EA`,n

X

D

E
A`,n A`,n , PT (Ã`,n )hPT (Ã`,n ), hm∗ i − E`,n

(243)

(`,n)

Given those comments. One can proceed with the rest of the proof. The second norm (243) is bounded
through lemma 9 below. This lemma follows from a direct application of proposition 1 as the outer matrices
are averaged out. It is proved in section 3.7.
√
Lemma 9. Let ĉ`,n be defined as in section 1.3 (including the scaling L), with Cn [`, k] ∼ N (0, 1/L) i.i.d.
gaussian. Let A`,n = f` ĉ∗`,n and let the operator A : CL×KN 7→ CL×N be defined from those matrices A`,n
as in (8). The coherences µ2m and µ2h are defined as in (10) and (11). Then for any constants δ9 , β9 , as
soon as L & β9 (1/δ9 )Kµ2h and N & β9 (1/δ9 )µ2m
PT⊥ EA`,n

X

D

E
A`,n A`,n , PT (Ã`,n )hPT (Ã`,n ), hm∗ i − E`,n
≤ δ9

(`,n)

with probability at least 1 − (LN )

−β9

.
34

(244)

The first norm in (237) will be bounded by lemmas 10 and 11. The idea behind these two lemmas is as follows.
Since the variables A`,n are independent from the variables Ã`,n , we can now condition the probability of
having a small norm (237) on the Ã`,n and therefore consider those variables as fixed over a first phase. In
particular, if we let


X
H`,n ≡ 
PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ i − E`0 ,n0  .
(245)
(`0 ,n0 )6=(`,n)
2
Let us now define the coherences µ2H`,n , ρ2H`,n and νH
, for given matrices H`,n as
`,n

µ2H`,n = L · sup

N
X

kH`,n [`2 , ∼ n2 ]k2 ,

(246)

kH`,n [`2 , ∼ n2 ]k2

(247)

`2 n =1
2

ρ2H`,n

= N · sup
n2

L
X
`2 =1

2
νH
= LN · sup kH`,n [`2 , ∼ n2 ]k2
`,n

(248)

`2 ,n2

and let the event E9 be defined as
n
E9 ≡ µ2H`,n . µ2h , ρ2H`,n . µ2m ,

2
νH
. µ2h µ2m
`,n

o

(249)

one can now bound the probability


X
X
Pr  PT⊥
A`,n hA`,n , H`,n i − PT⊥ EA`,n
A`,n hA`,n , H`,n i ≥ δ 
(`,n)

(250)

(`,n)

as




Pr  PT⊥ EA`,n

X

A`,n hA`,n , H`,n i − PT⊥ EA`,n

(`,n)

X

A`,n hA`,n , H`,n i ≥ δ 

(`,n)


≤ Pr  PT⊥ EA`,n


X

A`,n hA`,n , H`,n i − PT⊥ EA`,n

(`,n)

X

(251)

A`,n hA`,n , H`,n i ≥ δ E9  + Pr (E9c )

(`,n)

Because of the decoupling argument, the event E9 is now independent of the matrices A`,n and can be fixed
while bounding the norm. Bounding the two probabilities on the RHS of (251) is the point of lemmas 10
and 11 below which are respectively proved in sections 3.5 and 3.6.
Lemma 10. Let the linear map A be defined as in (8) with A`,n (X) = hA`,n , Xi = hf` c∗`,n , Xi. For any
fixed matrices X`,n independent of the A`,n . Let the coherences of the X`,n be defined as,
µ2X`,n = L · sup

N
X

kX`,n [`2 , ∼ n2 ]k2 ,

(252)

kX`,n [`2 , ∼ n2 ]k2

(253)

`2 n =1
2

ρ2X`,n = N · sup
n2

L
X
`2 =1

2
νX
= LN · sup kX`,n [`2 , ∼ n2 ]k2
`,n

(254)

`2 ,n2

Correspondingly, let µ2X , ρ2X and ν 2X denote the supremas over (`, n) of these quantities,
µ2X = sup µ2X`,n ,
`,n

ρ2X = sup ρ2X`,n ,
`,n

35

and

2
ν 2X = sup νX
.
`,n
`,n

(255)

One can write,
L X
N
X

A∗`,n A`,n (X`,n ) − E

`=1 n=1

L X
N
X

A∗`,n A`,n (X`,n ) .

`=1 n=1

(s
β10 max

)

µ2X K
ρ2X p
ν 2X
+
log(LKN ),
K log(LKN ))
(256)
L
N
LN

2
with probability at least 1 − (LN )−β10 . In particular, letting X`,n = H`,n , assuming µ2H`,n . µ2h , νH
.
`,n
2 2
2
2
2
2
µh µm and ρH`,n . µm , and taking L & β10 (1/δ10 )Kµh as well as N & µm , one can make the norm on
the LHS of (256) less than δ10 . Consequently, the first probability on the RHS of (251) can be bounded by
(LN )−β10 .

Lemma 11. Let A`,n be defined as in (8), and let X`,n (or equivalently H`,n ) be defined as
X
X
X`,n =
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i
(`0 ,n0 )6=(`,n)

(257)

(`0 ,n0 )6=(`,n)

2
where the matrices A`0 ,n0 are defined as in (5). The coherences µ2X`,n , ρ2X`,n and νX
are defined as
`,n
in (246), (247) and (248). Then we have

µ2X`,n . µ2h ,

ρ2X`,n . µ2m ,

and

2
νX
. µ2m µ2h ,
`,n

∀(`, n) ∈ [L] × [N ],

(258)

where each event holds with probability at least 1 − (LN )−β11 for any constant β11 as soon as L & β11 Kµ2h ,
N & β11 µ2m .
Before going through the decoupling argument underlying the bound Pr(E8 ) ≤ Pr(E8d ), we combine all
previous results and conclude the proof of lemma 8. Combining the decoupling argument relating E8 and E8c
together with the sequence of bounds (236) to (238), the sequence (240) to (243), the fact that Pr(kA + Bk ≥
δ) ≤ Pr(kAk ≥ δ/2) + Pr(kBk ≥ δ/2), as well as the results of lemmas 10 and 11 and 9, we get,

+
*
X
X
PT⊥
A`,n A`,n , 
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 
(259)
(`0 ,n0 )6=(`,n)

(`,n)

*
−PT⊥ E

X

A`,n



A`,n , 

+
X

PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E`0 ,n0 

.δ

(260)

(`0 ,n0 )6=(`,n)

(`,n)

with probability at least 1 − c(LN )−β where c, β > 0 are positive constants. This concludes the proof of
lemma 8.
3.4.1

Decoupling strategy

In this section, we explain the decoupling argument used in the proof of lemma 8 to relate E3 and E4 . This
argument is derived from Theorem 3.1.1 and 3.4.1 in [19, 20]. A similar result is given for the expectation
in [51] (Theorem 1.). We provide the proof of this result as well for completeness although the proof is
essentially the same as the proof given in [19].
Proposition 6 (Decoupling inequality for U-statistics (matrix version)). Let us recall the norm,
E = PT⊥

X

X

A`,n hA`,n , PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ ii

(261)

(`,n) (`0 ,n0 )6=(`,n)

−PT⊥ E

X

X

A`,n hA`,n , PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ ii

(`,n)6=(`0 ,n0 ) (`0 ,n0 )6=(`,n)

36

(262)

and let us consider the following norm where the sequence of variables A`0 ,n0 have been replaced by an
independent copy of the sequence {A`,n }(`,n)∈[L]×[N ] ,
Ẽ = PT⊥

X

X

(`,n)

(`0 ,n0 )6=(`,n)

X

−PT⊥ E

A`,n hA`,n , PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ ii

X

A`,n hA`,n , PT (Ã`0 ,n0 )hPT (Ã`0 ,n0 ), hm∗ ii

(263)

(264)

(`,n)6=(`0 ,n0 ) (`0 ,n0 )6=(`,n)



Then Pr(E ≥ t) ≤ CPr Ẽ ≥ t/C where C is a constant.
Proof. The proof is, in essence, the same as the proof given by de la Penã and Giné in [19]. The only
difference lies in its transposition to matrix valued random variables. Let X`,n and Y`0 ,n0 be defined as
follows,
X`,n = (vec(A`,n ) ⊗ vec(A`,n )∗ )

(265)

Y`0 ,n0 = (vec(PT (A`0 ,n0 ))hPT (A`0 ,n0 ), hm∗ i)

(266)

The norm (264) can read compactly as
X
X
X
PT⊥ (X`,n Y`0 ,n0 ) − E
(`,n) (`0 ,n0 )6=(`,n)

=
=

X

X

(`,n)

(`0 ,n0 )6=(`,n)

X

X

X

PT⊥ (X`,n Y`0 ,n0 )

(`,n) (`0 ,n0 )6=(`,n)

(vec(A`,n ) ⊗ vec(A`,n )) (vec(PT (A`0 ,n0 ))hPT (A`0 ,n0 ), hm∗ i) − C`,n,`0 ,n0

(267)

(vec(A`,n ) ⊗ vec(A`,n )) PT (vec(A`0 ,n0 ) ⊗ vec(A`0 ,n0 )∗ ) vec(hm∗ ) − C`,n,`0 ,n0

(268)

(`,n) (`0 ,n0 )6=(`,n)





= h`,n,`0 ,n0 vec(A`,n ) ⊗ vec(A`,n ), vec(A`0 ,n0 ) ⊗ vec(A`0 ,n0 )

(269)

= h`,n,`0 ,n0 (X`,n , X`0 ,n0 )

(270)

where C`,n,`0 ,n0 are constant (i.e. deterministic) matrices and h`,n,`0 ,n0 (X, Y ) is a function defined as
h`,n,`0 ,n0 (X`,n , X`0 ,n0 ) = PT⊥ X`,n PT X`0 ,n0 vec(hm∗ ) − C`,n,`0 ,n0

(271)

PT is the matrix formulation of the projector PT defined in (13). Let {X̃`,n }`,n denote an independent copy
of the sequence {X`,n }`,n . Let {ε`,n } denote a sequence of independent Rademacher random variables and
consider the accompanying sequences of matrices {Z`,n } and {Z̃`,n } defined as


X`,n if ε`,n = 1
X̃`,n if ε`,n = 1
Z`,n =
,
Z̃`,n =
(272)
X̃`,n if ε`,n = −1
X`,n if ε`,n = −1
The law D of (Z1,1 , . . . , ZL,N ) is the same as the law of (A(1,1) , . . . , A(L,N ) ) and similarly,
D(Z1,1 , . . . , ZL,N , Z̃1,1 , . . . , Z̃L,N ) = D(A1,1 , . . . , AL,N , Ã1,1 , . . . , ÃL,N )
Both of these laws are given by (P1,1 , . . . , PL,N ) = P LN and (P1,1 , . . . , PL,N )2 = P 2LN where P is the law
of each A`,n . If we let X denote the sigma algebra σ(A`,n , Ã`,n ) generated by the sequences A`,n and Ã`,n ,
we have
h(Z`,n , Z`0 ,n0 |X ) = h(X`,n , X`0 ,n0 )δ(ε`,n = 1, ε`0 ,n0 = 1)
+ h(X`,n , X̃`0 ,n0 )δ(ε`,n = 1, ε`0 ,n0 = −1)
+ h(X̃`,n , X`0 ,n0 )δ(ε`,n = −1, ε`0 ,n0 = 1)
+ h(X̃`,n , X̃`0 ,n0 )δ(ε`,n = −1, ε`0 ,n0 = −1)
37

(273)

which, when averaging over the Rademacher sequence, implies

1
h(X`,n , X`0 ,n0 ) + h(X`,n , X̃`0 ,n0 ) + h(X̃`,n , X`0 ,n0 ) + h(X̃`,n , X̃`0 ,n0 )
Eε h(Z`,n , Z`0 ,n0 |X ) =
4

(274)

Equivalently, we will use the following relation later in the proof,
4h(Z`,n , Z̃`0 ,n0 ) = (1 + ε`,n )(1 + ε`0 ,n0 )h(X`,n , X̃`0 ,n0 )
+ (1 + ε`,n )(1 − ε`0 ,n0 )h(X`,n , X`0 ,n0 )
+ (1 − ε`,n )(1 + ε`0 ,n0 )h(X̃`,n , X̃`0 ,n0 )

(275)

+ (1 − ε`,n )(1 − ε`0 ,n0 )h(X̃`,n , X`0 ,n0 )
Now we can follow the proof of Theorem 3.4.1 in [19] (Theorem 1 in [20]). We will need the following result
which follows as a special case of Theorem 1.1.3. in [19]. For X, Y and Z i.i.d., we have
Pr (kXk > t) = Pr (k(X + Y ) + (X + Z) − (Y + Z)k > 2t)






2t
2t
2t
+ Pr kX + Zk >
+ Pr kY + Zk >
≤ Pr k(X + Y )k >
3
3
3


2t
≤ 3Pr kX + Y k >
3
Using (278) with the independent sequence {X̃`,n }`,n , one can write


X
X
Pr 
h(X`,n , X`0 ,n0 ) > t

(276)
(277)
(278)

(279)

`,n (`0 ,n0 )6=(`,n)


≤ 3Pr 

X

X

`,n (`0 ,n0 )6=(`,n)


2t 
h(X`,n , X`0 ,n0 ) + h(X̃`,n , X̃`0 ,n0 ) >
3

(280)

Splitting the norm,
X

X

h(X`,n , X`0 ,n0 ) + h(X̃`,n , X̃`0 ,n0 )

(281)

`,n (`0 ,n0 )6=(`,n)

≤

X

X

h(X`,n , X`0 ,n0 ) + h(X̃`,n , X̃`0 ,n0 ) − h(X̃`,n , X`0 ,n0 ) − h(X`,n , X̃`0 ,n0 )

(282)

h(X̃`,n , X`0 ,n0 ) − h(X`,n , X̃`0 ,n0 )

(283)

`,n (`0 ,n0 )6=(`,n)

+

X

X

`,n (`0 ,n0 )6=(`,n)

Susbtituting into (280) yields,


X
X
Pr 
h(X`,n , X`0 ,n0 ) > t

(284)

`,n (`0 ,n0 )6=(`,n)


≤ 3Pr 

X

X

`,n (`0 ,n0 )6=(`,n)


t
h(X`,n , X`0 ,n0 ) + h(X̃`,n , X̃`0 ,n0 ) − h(X̃`,n , X`0 ,n0 ) − h(X`,n , X̃`0 ,n0 ) >
3
(285)


+ 3Pr 


X

X

h(X̃`,n , X`0 ,n0 ) − h(X`,n , X̃`0 ,n0 ) >

`,n (`0 ,n0 )6=(`,n)

38

t
3

(286)




Pr 

X

X

`,n

(`0 ,n0 )6=(`,n)

h(X`,n , X`0 ,n0 ) > t


≤ 3Pr 

X

X

`,n (`0 ,n0 )6=(`,n)

(287)


t
h(X`,n , X`0 ,n0 ) + h(X̃`,n , X̃`0 ,n0 ) − h(X̃`,n , X`0 ,n0 ) − h(X`,n , X̃`0 ,n0 ) >
3
(288)




t
(289)
6
`,n (`0 ,n0 )6=(`,n)


X
X
t
h(X`,n , X̃`0 ,n0 ) > 
+ 3Pr 
(290)
6
0
0
`,n (` ,n )6=(`,n)


X
X
t
h(X`,n , X`0 ,n0 ) + h(X̃`,n , X̃`0 ,n0 ) − h(X̃`,n , X`0 ,n0 ) − h(X`,n , X̃`0 ,n0 ) > 
≤ 3Pr 
3
0
0

+ 3Pr 

X

X

h(X̃`,n , X`0 ,n0 ) >

`,n (` ,n )6=(`,n)

(291)

+ 6Pr 


X

X

h(X̃`,n , X`0 ,n0 ) >

`,n (`0 ,n0 )6=(`,n)

t
6

(292)

In (292) we use the fact that h(X`,n , X̃`0 ,n0 ) and h(X̃`,n , X`0 ,n0 ) have the same distribution. We are thus
left with bounding the U -statistics (291). Let us use U to denote the sum,
X
X
U =4
h(Z`,n , Z̃`0 ,n0 )
(293)
`,n (`0 ,n0 )6=(`,n)

n
o
and let U = h(X`,n , X`0 ,n0 ), h(X`,n , X̃`0 ,n0 ), h(X̃`,n , X`0 ,n0 ), h(X̃`,n , X̃`0 ,n0 ) . Using (275), one can see that
U conditionned on U is a second order Rademacher chaos. finally, to conclude, we use a conditional version
of Jensen’s inequality following from Theorem 3.3.6 [19]. Let us use T to denote the U-statistics,
i
X X h
T =
h(X`,n , X`0 ,n0 ) + h(X`,n , X̃`0 ,n0 ) + h(X̃`,n , X`0 ,n0 ) + h(X̃`,n , X̃`0 ,n0 )
(294)
`,n `0 ,n0 6=(`,n)

Clearly, (273) and (275) show that we have E {U |U} = T . Following the proof of [20], we now use the
following Theorem from [19],
Theorem 2 (Theorem 3.3.6 in [19]). Let F be a normed linear space and let n ∈ N. Let {εi }1≤i≤n denote
a sequence of independent Rademacher variables. Let X = {X0 , X1 , . . . , Xi , . . . , Xi1 ,...,id }, 1 ≤ i, ik ≤ n
denote a sequence of random variables taking values in a normed vector space F. Let ε and X be defined
on different factors of a product probability space. Then for every p there exists a constant Cpd such that the
conditional Rademacher chaos
X = X0 +

n
X
i=1

X

Xi ξi + . . . +

Xi1 ,...,id ε1 . . . εid

(295)

1≤i1 6=i2 6=...6=id ≤n

satisfies
Pr∗ (E {kXkp |X} ≥ 2tp ) ≤ Cp,d Pr∗ (kXk > t) .
Here the notation Pr∗ (X) means the outer probability of the event X.

39

(296)

Using Theorem 2 for p = 1, one can write,
Pr (kT k > t) = Pr (kE{U |U}k > t)

(297)

≤ cPr (2kU k > t)


XX
h(Z`,n , Z̃`0 ,n0 )k
= cPr 2k4

(298)
(299)

(`,n) `0 ,n0





= cPr 2k4

XX
(`,n)

h(X`,n , X̃`0 ,n0 )k

(300)

`0 ,n0

In the last line, we use the fact that (Z`,n , Z̃`0 ,n0 ) and (X`,n , X̃`0 ,n0 ) have the same distribution. This
concludes the proof of proposition 6.

3.5

Proof of lemma 10

Building upon the result of proposition 6 and the discussion of section 3.4, we now prove the remaining two
lemmas 10 and 11
Lemma 10. Let the linear map A be defined as in (8) with A`,n (X) = hA`,n , Xi = hf` c∗`,n , Xi. For any
fixed matrices X`,n independent of the A`,n . Let the coherences of the X`,n be defined as,
µ2X`,n = L · sup

N
X

kX`,n [`2 , ∼ n2 ]k2 ,

(252)

kX`,n [`2 , ∼ n2 ]k2

(253)

`2 n =1
2

ρ2X`,n = N · sup
n2

L
X
`2 =1

2
νX
= LN · sup kX`,n [`2 , ∼ n2 ]k2
`,n

(254)

`2 ,n2

Correspondingly, let µ2X , ρ2X and ν 2X denote the supremas over (`, n) of these quantities,
µ2X = sup µ2X`,n ,
`,n

ρ2X = sup ρ2X`,n ,
`,n

and

2
ν 2X = sup νX
.
`,n

(255)

`,n

One can write,
L X
N
X
`=1 n=1

A∗`,n A`,n (X`,n ) − E

L X
N
X

A∗`,n A`,n (X`,n ) .

`=1 n=1

(s
β10 max

)

µ2X K
ρ2X p
ν 2X
+
K log(LKN ))
(256)
log(LKN ),
L
N
LN

2
.
with probability at least 1 − (LN )−β10 . In particular, letting X`,n = H`,n , assuming µ2H`,n . µ2h , νH
`,n
2 2
2
2
2
2
µh µm and ρH`,n . µm , and taking L & β10 (1/δ10 )Kµh as well as N & µm , one can make the norm on
the LHS of (256) less than δ10 . Consequently, the first probability on the RHS of (251) can be bounded by
(LN )−β10 .

Proof. Let X`,n denote arbitrary matrices of size L × N K and let X[`, ∼ n] denote the submatrix given by
considering the `th row of X and the columns (n − 1)K + 1 ≤ k ≤ nK. I.e
X̂[`, ∼ n] ≡

K
X

f`∗ X[:, (n − 1)K + k](en ⊗ ek ).

k=1

40

(301)

ek is a K-dimensional zero vector with its k th entry set to 1 and en is the N -dimensional zero vector with
its nth entry set to 1. The norm on the LHS of (256) can expand as
X

Z`,n =

X

A`,n hA`,n , X`,n i − EA`,n hA`,n , X`,n i

(302)

A`,n hA`,n , X`,n i − EA`,n hA`,n , X`,n i

(303)

(`,n)

(`,n)

=

X
(`,n)

We will use proposition (1) to derive the tail bound of lemma 10. To use this proposition, we start by deriving
the bound for the variance. Note that the Z`,n are not hermitian as we have A`,n A∗`,n = kc`,n k2 f` f`∗ and
A∗`,n A∗`,n = kf` k2 c`,n c∗`,n . The first variance bound gives
X

X

∗
=
EZ`,n Z`,n

2

EA`,n A∗`,n |hA`,n , X`,n i|2 − |EA`,n hA`,n , X`,n i|

(`,n)

(`,n)

X

≤

EA`,n A∗`,n |hA`,n , X`,n i|2

(`,n)

X

≤

Ef` f`∗ kc`,n k2 |hc`,n , X̂`,n [`, ∼ n]i|2

(`,n)

≤K

X

f` f`∗ kX̂`,n [`, ∼ n]k2

(304)

(`,n)

Now using the coherences defined in (252) to (254) and their supremas defined in (255), the last line of (304)
can be reduced to
X

∗
EZ`,n Z`,n
≤K·

(`,n)∈Γp

µ2X
L

(305)

For the second term, a similar argument gives,

E

X

X

∗
Z`,n
Z`,n = E

(`,n)

kf` k2 c`,n c∗`,n |hf` c`,n , X`,n i|2 ,

(`,n)∈Γp

≤ E

X

kf` k2 c`,n c∗`,n |hc`,n , X̂`,n [`, ∼ n]i|2

(`,n)

.

X

kX̂`,n [`, ∼ n]k2 (en e∗n ) ⊗ IK +

(`,n)


.

ρ2X
N

X

T
X̂`,n
[`, ∼ n]X̂ `,n [`, ∼ n]

(`,n)


.

(306)

The last line comes from orthogonality of the X[`, ∼ n1 ] for X[`, ∼ n2 ] for distinct n1 , n2 . Combining (305)
with (306), we get the following bound on σ,

 2
µX K
ρ2X
σ.
+
(307)
L
N
41

We now compute the bound on the Ψ1 norm. Following the same argument as in previous lemmas, we get
kkA`,n hA`,n , X`,n ikkψ1 ≤ kf` kkkc`,n kkψ2 k|hc`,n , X̂`,n [`, ∼ n]i|kψ2
√
≤ KkX`,n [`, ∼ n]k
√
νX
≤ K√
LN

(308)
(309)

The last line follows from the discussion at the end of section 2.2. Substituting those bounds into proposition 1, we have
X

A∗`,n A`,n X`,n − X`,n .

(`,n)

(s 
max

)

µ2X K
ρ2X p
νX √
+
t + log(LKN ), √
K log(LKN )(t + log(LKN ))
L
N
LN

with probability at least 1 − e−t . Taking t = β log(LN ) concludes the proof of lemma 10.

3.6

Proof of lemma 11

Again, we recall the statement of lemma 11 for clarity before proving this lemma.
Lemma 11. Let A`,n be defined as in (8), and let X`,n (or equivalently H`,n ) be defined as
X
X
X`,n =
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i − E
PT (A`0 ,n0 )hPT (A`0 ,n0 ), hm∗ i
(`0 ,n0 )6=(`,n)

(257)

(`0 ,n0 )6=(`,n)

2
where the matrices A`0 ,n0 are defined as in (5). The coherences µ2X`,n , ρ2X`,n and νX
are defined as
`,n
in (246), (247) and (248). Then we have

µ2X`,n . µ2h ,

ρ2X`,n . µ2m ,

and

2
νX
. µ2m µ2h ,
`,n

∀(`, n) ∈ [L] × [N ],

(258)

where each event holds with probability at least 1 − (LN )−β11 for any constant β11 as soon as L & β11 Kµ2h ,
N & β11 µ2m .
2
Proof. The proof is similar for any of the three relations on either νX
, µ2X`,n or ρ2X`,n . We thus only show
`,n
2
the relation νX
≤ µ2h µ2m . From the definition of νX`,n , we have,
`,n
2
νX
`

1 ,n1

/LN = sup kX̂`1 ,n1 [`2 , ∼ n2 ]k2F

(310)

`2 ,n2

2

= sup
`2 ,n2

X

X

∗

F PT (A`,n )hA`,n , hm i − E

k2 ∼n2 (`,n)6=(`1 ,n1 )

X

∗

F PT (A`,n )hA`,n , hm i

(`,n)6=(`1 ,n1 )

(311)
`2 ,k2

The notation above really means the modulus squared of the entry (`2 , k2 ) of the underlying matrix. We
start by bounding the sum of terms inside the modulus. Let us denote each of the terms of this sum by z`,n .
We have
z`,n = (F PT (A`,n )hA`,n , hm∗ i − EF PT (A`,n )hA`,n , hm∗ i)(`2 ,k2 ) .

42

(312)

for fixed (`, n) 6= (`1 , n1 ), `2 , n2 and k2 . We will use proposition 1 to bound the sum
start by deriving a bound on the variance of the z`,n .

P

(`,n)6=(`1 ,n1 ) z`,n .

We

2

2

E|z`,n |2 = E |F PT (A`,n )hA`,n , hm∗ i|`2 ,k2 − |EF PT (A`,n )hA`,n , Xi|`2 ,k2
2

≤ E |F PT (A`,n )hA`,n , Xi|`2 ,k2 .
Let us use P̂T to denote the Fourier transform of the projector onto the tangent space, P̂T = F PT . From
definition (13), we have
P̂T (A`,n )hA`,n , hm∗ i

2

2
`2 ,k2

≤ P̂T (A`,n )

(k2 ,`2 )

|hc`,n , ĥ[`]mn i|2
2

≤ |ĥ[`]|h[`2 ]c`,n [k2 ] + e` [`2 ]hc`,n , mn im∗ [k2 ] − ĥ[`2 ]ĥ[`]hc`,n , mim∗ [k2 ] |hc`,n , ĥ[`]m∗n i|2


. |ĥ[`]|2 |ĥ(`2 )|2 |c`,n [k2 ]|2 |hc`,n , ĥ[`]m∗n i|2

+ |e` [`2 ]|2 |hc`,n , mn i|2 |mn2 [k2 ]|2 |hc`,n , ĥ[`]m∗n i|2


+ |ĥ[`2 ]|2 |ĥ[`]|2 |hc`,n , mi|2 |mn2 [k2 ]|2 |hc`,n , ĥ[`]m∗n i|2 .
Taking the expectation gives,
E|z`,n |2 . |ĥ[`]|2 |ĥ[`2 ]|2 δ(n, n2 )|ĥ[`]m∗n [k2 ]|2 + δ(`2 , `)kmn k2 |mn [k2 ]|2 kĥ[`]m∗n k2
+ |ĥ[`2 ]|2 |ĥ[`]|2 kmn k2 |mn [k2 ]|2 kĥ[`]m∗n k2 ,

 µ2 µ2
h m
. |ĥ[`]|2 |ĥ[`2 ]|2 δ(n, n2 ) + δ(`2 , `)kmn k2 |mn [k2 ]|2 + |ĥ[`2 ]|2 |ĥ[`]|2 kmn k2 |mn [k2 ]|2
LN
Summing over the measurements, we get
σ2 =

X

E|z`,n |2 . (|ĥ[`2 ]|2 + |m[k2 ]|2 )

(`,n)6=(`1 ,n1 )

µ2h µ2m
LN

(313)

In order to apply proposition 1, we are left with computing the bound on the Ψ1 norm of the z`,n . Again,
we can write
|F PT (A`,n )hA`,n , hm∗ i|`2 ,k2


= ĥ[`]ĥ[`2 ]c`,n [k2 ] + e` [`2 ]hc`,n , mim[k2 ] − ĥ[`2 ]ĥ[`]hc`,n , mim[k2 ] hc`,n , ĥ[`]mn i
. ĥ[`]ĥ[`2 ]c`,n [k2 ]hc`,n , ĥ[`]mn i
+ e` [`2 ]hc`,n , mim[k2 ]hc`,n , ĥ[`]mn i
+ ĥ[`2 ]ĥ[`]hc`,n , mim[k2 ]hc`,n , ĥ[`]mn i
Each of the terms above can be bounded by again using the fact that the product of two subgaussian random
variables follows a subexponential distribution and that the Orlicz one norm of this product can be upper
bounded by the product of the Orlicz-2 norms of each of the subgaussian random variables. In other words,
kXY kΨ1 ≤ kXkΨ2 kY kΨ2 . For each of the three terms in the expression of z`,n − Ez`,n , we can use the fact
that kc`,n [k 0 ]kΨ2 ≤ c, hc`,n , mn iis a mixture of centered gaussians. Note that the corresponding variable
|hc`,n , mn i|2 is a chi-squared that is sub-exponential with parameters (2, 4kmn k2 ) as we saw in the proof of
lemma 1 and using proposition 4, we get,

µ µ
h m
(314)
k |PT (A`,n )hA`,n , hm∗ i|`2 ,n2 kΨ1 . |ĥ[`]ĥ[`2 ]| + e` [`2 ]m[k2 ]kmn k √
LN
Finally for the expectation, simply recall that kX − EXkψq ≤ 2kXkψq
43

Using the bounds (314) and (313), one can now apply proposition 1 to derive the tail bound guaranteeing
that each each term within the modulus of (311) remains sufficiently small. Applying proposition 1 gives,
X

F PT (A`,n )hA`,n , hm∗ i − EF PT (A`,n )hA`,n , hm∗ i

(`,n)6=(`1 ,n1 )

` ,k

2 2
q


 µ µ
√h m
≤ β max
(|ĥ[`2 ]|2 + |m[k2 ]|2 ) , log(N ∨ L) |ĥ[`]ĥ[`2 ]| + e` [`2 ]m[k2 ]kmk
LN

with probability at least 1 − c3 (LN )−β for some constant c3 . We then square the modulus and take the sum
over the indices k2 ∼ n2 , which gives
X
X
sup
|
z`,n |2`2 ,k2 .
(315)
`2 ,n2

k2 ∼n2 (`,n)6=(`1 ,n1 )

2  µ2 µ2
h m
.
max
+ kmn2 k , K log (N ∨ L)
+ |m[k2 ]|kmn k
LN
o µ2 µ2
n
h m
. max Kkĥk2∞ + kmn2 k2 , Kkhk4∞ + kmn k2 kmn2 k2
LN
PK
With probability at least 1 − 2c3 (LN )−β . In (316), we used k2 =1 |mn2 [k2 ]|2 = kmn2 k2


Kkĥk2∞

2

2



kĥk2∞

(316)
(317)

The bound (317) can be made smaller than any constant from the coherences µ2h /L = kĥk2∞ and µ2m /N =
maxn kmn k2 as soon as L & Kµ2h and N & µ2m up to log factors. This concludes the proof of lemma 11.

3.7

Proof of lemma 9

√
Lemma 9. Let ĉ`,n be defined as in section 1.3 (including the scaling L), with Cn [`, k] ∼ N (0, 1/L) i.i.d.
gaussian. Let A`,n = f` ĉ∗`,n and let the operator A : CL×KN 7→ CL×N be defined from those matrices A`,n
as in (8). The coherences µ2m and µ2h are defined as in (10) and (11). Then for any constants δ9 , β9 , as
soon as L & β9 (1/δ9 )Kµ2h and N & β9 (1/δ9 )µ2m
PT⊥ EA`,n

X

D

E
A`,n A`,n , PT (Ã`,n )hPT (Ã`,n ), hm∗ i − E`,n
≤ δ9

(244)

(`,n)

with probability at least 1 − (LN )

−β9

.

Proof. We use proposition 1. The terms within the norm (244) expand as

X 
f` |ĥ[`]|2 c̃∗`,n + kf` k2 hmn , c̃∗`,n im∗n − |ĥ[`]|2 hmn , c̃∗`,n im∗n |ĥ[`]|hc̃`,n , mn i

(318)

`,n



− Ef` |ĥ[`]|2 c̃∗`,n + kf` k2 hmn , c̃∗`,n im∗n − |ĥ[`]|2 hmn , c̃∗`,n im∗n |ĥ[`]|hc̃`,n , mn i

(319)

We will bound each of the first two term through proposition 1. The bound on the last term follows from
the first two. For each of those two terms, letting aside the expectation, the first and second terms in the
variance bound can be expressed as
E

X

f` f`∗ |ĥ[`]|4 kc̃`,n k2 |ĥ[`]|2 |hc̃`,n , mn i|2 .

`,n

E

X

Kµ6h µ2m
L3 N

kf` k2 c`,n c∗`,n |ĥ[`]|4 |ĥ[`]|2 |hc̃`,n , mn i|2 .

`,n

44

µ4h
(kmn k2 + 1)
L2

(320)

(321)

for the first term, and

E

X

f` f`∗ kf` k2 |hmn , c̃∗`,n i|4 kmn k2 |ĥ[`]|2 .

`,n

E

X

kf` k4 |hmn , c̃∗`,n i|4 mn m∗n |ĥ[`]|2 .

`,n

For the Orlicz norm, simply note that


f` |ĥ[`]|2 c̃∗`,n |ĥ[`]|hc̃`,n , mn i

ψ1

µ2h
kmn k4
L

µ4h
(kmn k2 + 1)
L2



f` |ĥ[`]|2 c̃∗`,n |ĥ[`]|

.

(322)

ψ2

(323)

khc̃`,n , mn ikψ2

√ µh µ2
h
. kmn k K √
L L

(324)
(325)

as well as

f` kf` k2 hmn , c̃∗`,n im∗n |ĥ[`]|hc̃`,n , mn i

∗

f` kf` k2 hmn , c̃∗`,n imn

.

ψ2

ψ1

|ĥ[`]|hc̃`,n , mn i

(326)
ψ2

3 µh
. kmn k √
L

(327)

All of these terms can be made less than δ whenever L & (1/δ)βµ2h K and N & µ2m . The conclusion follows
from applying proposition 1 with t = β log(LN ).

3.8

Proof of lemma 4

The proof of lemma 4 follows the approach in [14]. simply note that we have,
k

k

kPT⊥ (A∗ APT ) (PT − PT A∗ APT ) hm∗ k ≤ kPT⊥ (A∗ APT ) k (PT − PT A∗ APT ) hm∗
k

≤ kPT⊥ (A∗ APT ) k kPT − PT A∗ APT k khm∗ k
k

≤ kA∗ kkPT (A)k kPT − PT A∗ APT k khm∗ k

(328)
(329)
(330)

Now taking the sum and using the fact that kA∗ Ak ≤ kAk2 . K log(L ∨ N ) as well as
|kAPT k2 − kPT k| = |kPT (A∗ A)PT k − kPT k| . δ
gives,
∞
X

k

kPT⊥ (A∗ APT ) (PT − PT A∗ APT ) hm∗ k ≤

k=k0

p
kPT − PT A∗ APT kk0
K log(L ∨ N )
khm∗ k
1 − kPT − PT A∗ APT k

(331)

p
Now we can use the fact that lemma 1 gives kPT −PT A∗ APT k . µ2h K/L + µ2m /N . δ as soon as L & Kµ2h
and N & µ2m , the denominator is always larger than 1 − δ. To conclude, it suffices to again apply lemma 1
on the numerator to get
∞
X

kPT⊥

k

p
(A APT ) (PT − PT A APT ) hm k ≤ K log(L ∨ N )
∗

∗

∗

k=k0

45



µ2h K
µ2
+ m
L
N

k0 /2
(332)

In particular, as soon as L & K 1+1/k0 and N & K 1/k0 , the bound (333) reduces to
∞
X

k

kPT⊥ (A∗ APT ) (PT − PT A∗ APT ) hm∗ k ≤ C

(333)

k=k0

for a constant C which can be taken arbitrarly small so as to satisfy the bound on kPT⊥ (Y )k. In our case,
as k0 = 2, we get the sample complexities L & K 3/2 µ2h and N & K 1/2 µ2m .

4

Numerical simulations

In this section we provide additional numerical experiments and study the resulting phase diagrams in order
to quantify the probability of success for the formulation (7). To conduct the numerical experiments, we
first write problem (7) as a trace minimization problem (see [24]). Introducing V as a proxy for the rank
one matrix V0 ,




V11 V12
hhH hmH
V =
≈ V0 =
,
V21 V22
mhH mmH
one can write problem (7) as
min

Tr(V11 ) + Tr(V22 )

s.t A(V12 ) = {yn },

V  0.

(334)

To handle reasonably large instances of (9), we follow the approach in [9] and replace the matrix V by a
low rank factorization,


H
R1
R1
V =
= RRH
(335)
R2
R2
where R1 ∈ CL×r and R2 ∈ CKN ×r for some rank r  L + KN . Introducing multipliers λ`,n for each one
of the constraints {hA`,n , Xi = yn [`]}`,n and σ > 0, we then minimize the augmented Lagrangian L(R)
associated to (334) with respect to R1 and R2 ,
L(R) =

X
Xσ
1
1
kR1 k2F + kR2 k2F +
λ`,n (hA`,n , R1 R2H i − yn [`]) +
(hA`,n , R1 R2H i − yn [`])2 .
2
2
2
`,n

(336)

`,n

In addition to dealing with fewer variables, the factorization introduced in (335) is also a very convenient
way to circumvent the constraint V  0 since V  0 ⇔ V = Y Y H for some Y . In particular, taking
r = L + KN is equivalent to solving the original semidefinite program (9). As explained in [9], one can thus
minimize (336) for a very small rank r  L + KN , check optimality with respect to the original problem and
then increase the rank in case the rank-r factorization doesn’t lead to the solution of the original problem.
In practice we don’t even need to do so and setting r = 4 is enough to achieve convergence to the minimizer
of (9).

4.1

Phase transition

The success rates for different values of the parameters L, N and K are shown in Fig. 3. The first set of
numerical experiments, corresponding to the first diagram of Fig. 3, shows the rate of success as a function
of both the size of the input subspace K and the size of the ambient space L, the number of input signals
N being set to 40. For each of the values of the pairs (K, L), 100 experiments were run by taking gaussian
i.i.d sensing matrices Cn with Cn [`, k] ∼ N (0, 1), as well as gaussian i.i.d vectors h and m. For each of the
100 experiments, we ran 40 iterations of the limited memory BFGS on the augmented Lagrangian (336).
An experiment was classified as a success (white/1) as opposed to failure (black/0) whenever the relative

46

difference between the obtained matrix X and the optimal solution X0 = hmH was less than 2%. In other
words,
kX − hmH kF
< .02
(337)
khmH kF
The second diagram of Figure 3 shows the rate of success as a function of K and N for the same experimental
framework. Here L is fixed to 800.
The third diagram shows the rate of success for various values of L and N for a fixed K = 40. It is worth
noting that as long as there is a sufficient number of columns (> 50), the recovery mainly depends on the
parameters K and L.
Finally, we provide a phase diagram highlighting the independence of the recovery vis a vis the sparsity
of the impulse response h (Sparsity was studied as an additional condition on the recovery for example
in [1, 36, 38]). Figure (4) shows that for fixed L and N , the sparsity has no influence on the recovery which
is essentially driven by the subspace dimension. This observation further motivates the extension to blind
super-resolution discussed below.

4.2

Convergence of the certificate

In this section, we provide numerial evidence supporting the choice of the ansatz (40). As explained in
section 2.3, the most natural certificate construction as soon as the normal map A∗ A concentrates to the
identity sufficiently fast, would be to use Y1 = A∗ Ahm∗ . In practice, however, such a certificate performs
poorly and does not achieve the sample complexity observed empirically in Fig. 3. In particular, it requires
poor ratios L/K and N/K in order for the condition PT (Y1 ) to be met in (29).
As soon as the map PT A∗ APT can be shown to be injective, the second ansatz Y2 = A∗ A(PT A∗ APT )−1 hm∗
provides a much better candidate to satisfy the first condition in (29). In this section, we compare the two
ansatz and provide numerical evidence that the conditions (29) are met for both, albeit in different sample
complexity regimes.
Fig. 5 first compares the projection of the two ansatz on the tangent space T for the particular choice of
dimensions given by K = 8, L = 80 and N = 20, and in a framework where the minimization program (9)
is known to recover the exact solution. In this regime, the direct ansatz fails to meet the condition in (29).
The use of the inverse however enables the second ansatz to satisfy the first condition exactly, in passing
highlighting the fact that in practice, the constant in the first condition of (29) can be taken arbitrarily
large. The norms corresponding to the first condition are respectively given by kPT (Y1 ) − hm∗ k ≈ 0.48 and
kPT (Y2 ) − hm∗ k ≈ 1e−13 .
Fig. 6 then compares the projection of the two ansatz onto the orthogonal complement T ⊥ . The two ansatz
exhibit comparable magnitude on the orthogonal complement T ⊥ and both of them satisfy the second
condition in (29) although the direct certificate seems to perform slightly better for the particular dimensions
considered when considering this second condition. The norms corresponding to the second condition in (29)
are given respectively by kPT⊥ (Y1 )k ≈ 0.67 and kPT⊥ (Y2 )k ≈ 0.72, thus showing that for the particular
choice K = 8, L = 80 and N = 20, the first candidate certificate fails to satisfy the optimality conditions
corresponding to a zero subgradient whether the second candidate certificate meets both of these conditions.

4.3

Applications

In this section, we discuss some of the application of the result of this paper. Among the many applications
of blind deconvolution, one should mention astronomical imaging, movie deblurring, seismic data processing,
super-resolution and medical imaging. In this section we discuss two of these applications in greater details:
Communication and blind super-resolution for medical imaging.

47

1

1

20

50

0.8
0.6

100
K

0.8

40
60

0.6

K 80
0.4

150
0.2

0.4

100
120

0.2

140
200

0

500

1000
L

1500

0

2000

50

100
N

150

200

1

100
0.8

200
300
L

0.6

400
500

0.4

600

0.2

700
0

800
50

100
N

150

200

Figure 3: Recovery of hmH from nuclear norm minimization for (top left) various values of K and L, (top
right) various values of K and N , and (bottom) various values of L and N . White (1) is used to indicate
success while black (0) is used to indicate failure. From the first figure (top left), the empirical recovery
rate is seen to increase when L increases and when K decreases, with an expected phase transition when L
is a multiple of K. From the remaining two figures, we see that a minimum (small) value of N is needed
for recovery, but past this threshold, K and L are the only parameters affecting recovery. All these phase
transitions can be explained by noting that the number of unknowns in the original problem is given by
KN + L whether the number of measurements is given by LN .
4.3.1

Communication and Rayleigh fading

When transmitting signals within a densely built environment such as in wireless communication, the multipath nature of the transmission arising from the many reflections that the signal will face will result in a
channel impulse response that can be considered completely arbitrary. A common assumption in this case
is to view each of the entries in the channel transfer matrix are identically and i.i.d., following a gaussian
distribution [35, 48]. This assumption is known as i.i.d. Rayleigh fading. In such a regime, it is clear that
• No sparsity or non vanishing assumption holds on the channel impulse response, and
• when some of the entries from the channel transfer matrix are very small, or vanish, linearized approaches such as [26, 39] do not hold.

48

80
100
120

K
140
160
180
200
50

100

150

S
Figure 4: Recovery of hmH from nuclear norm minimization for various values of K and S-sparse vectors h
with randomized supports. White (1) is used to indicate success while black (0) is used to indicate failure.
The phase transition here illustrates the fact that recovery through the nuclear norm minimization (9) is
essentially driven by the dimension of the inputs subspace as well as the coherences (11) and (10) more than
the sparsity S of the filter.
This paper is precisely interested in this regime as it removes the need for sparsity or non vanishing assumptions on the filter h. For more details on Rayleigh fading, see chapters 2, 7 in [48].
4.3.2

Blind super-resolution and medical imaging

To evaluate the interest of the nuclear norm minimization (9) for the super-resolution of signals, we consider
3 datasets. Before discussing each of these datasets and provide corresponding results of the nuclear norm
reconstruction for each of these datasets, we briefly recall the framework of blind super-resolution. In superresolution, one is interested in recovering a signal x from the result of its convolution with a known (ideal)
low pass filter h(x, y). In Fourier space, the convolution with the low pass filter thus reads as
ŷ(ω1 , ω2 ) = x̂(ω1 , ω2 )ĥ(ω1 , ω2 )

(338)

p
with ĥ(ω1 , ω2 ) = 0 for |ω| := ω12 + ω22 > Ωc (see for example [12]). The super-resolution framework
naturally extends to the problem of recovering a signal from the result of its convolution with an unknown
low pass filter. This extension, known as blind super-resolution, arises in many engineering applications
such as medical imaging, microscopy, astronomy or even seismic imaging. In these applications, acquisition
at lower resolution is often either performed by adding a blurring filter in order to reduce the influence
of the side lobes arising from the convolution with the sinc, or the ideal low pass filter is corrupted by
spurious reflections, noise or phase ambiguities, and is therefore only partially known. An illustration of the
convolution with an ideal low pass filter with and without gaussian blurr is given in Fig. 16. As explained in
section 1.2, recent developments on the subject include [54] in which Yang et al. discuss the problem in the
case where the output is given by a sum of spikes modulated by unknown (random) waveforms. The spikes
have to satisfy a minimum separation condition and the waveforms are assumed to be generated as random
combinations of random vectors satisfying an incoherence condition.
Although this paper primarily focuses on blind deconvolution, it also applies in the case of blind super49

Projection of the certificate on T

·10−2
3

hm∗
PT (Y1 )

2
1
0

20

40

60

80

100

120

140

160

180

200

220

240

Projection of the certificate on T

·10−2

hm∗
PT (Y2 )

2

1

0

20

40

60

80

100

120

140

160

180

200

220

240

Figure 5: Convergence of the dual certificates on T for (Top) the direct construction Y1 = A∗ A(hm∗ ) and
(bottom) the construction involving the inverse, Y2 = A∗ APT (PT A∗ APT )−1 hm∗ . Here we take K = 8,
N = 20 and L = 80. In this framework, the nuclear norm minimization program (9) is known to recover the
exact solution. The certificates have been downsampled by a factor 50 for clarity.

Comparison of the certificates Y2 = A∗ APT (PT A∗ APT )−1 hm∗ and Y1 = A∗ Ahm∗ on T ⊥
0.1

PT⊥ (Y2 )
PT⊥ (Y1 )

5 · 10−2

0

20

40

60

80

100

120

140

160

180

200

220

240

Figure 6: Convergence of the dual certificates on T ⊥ for the direct construction Y1 = A∗ A(hm∗ ) and the
construction involving the inverse, Y2 = A∗ APT (PT A∗ APT )−1 hm∗ . Here we take K = 8, N = 20 and
L = 80. In this framework, the nuclear norm minimization program (9) is known to recover the exact
solution. The certificates have been downsampled by a factor 50 for clarity.

50

resolution as is shown below. Using a gaussian ideal low pass filter, we now provide three examples of blind
super-resolution through the minimization program (9).
• In the first one, we consider the recovery of one dimensional wavelet trains from their convolution with
a one dimensional gaussian ideal low pass filter. The results are shown in Figs. 13 and 14. In these
examples, nuclear norm minimization is shown to recover the exact wavelet train and to remove the
spurious oscillations arising from the convolution with the sinc.
• In the second example, we consider the 3D Shepp Logan phantom of Schabel1 . For each of the xy slices from this phantom, we use as subspaces Cn , the wavelets corresponding to the K largest
coefficients from the discrete wavelet transform of each of the frames. We consider data generated
from the convolution of the slices with a gaussian (ideal) low pass filter similar to the one shown in
Fig. 16. The low resolution images together with the recovered (super-resolved) ones are shown in
Fig. 17, 18. Figs 7 and 8 illustrate the evolution of the recovery with the cutoff frequency ωc .
• Super-resolution is particularly interesting in medical imaging and Magnetic Resonance Imaging (MRI)
where a reduction in the number of samples acquired in k-space, the equivalent of the two-dimensional
Fourier space, leads to a reduction in the acquisition time. For an expensive imaging modality such
as MRI, a reduction in the acquisition time means a larger number of faster examinations and thus
a reduction of the waiting time for a fixed number of magnets. As a third example, we consider the
MRI dataset distributed by the Stanford data archive2 , originally generated at the University of North
Carolina. This MRI dataset consists of 109 slices of the skull, each being of size 256 by 256. We start
by applying a low pass gaussian filter to each slice in the volume. We then study reconstruction of these
slices through the nuclear norm minimization program (9). The results are shown in Figs 9, 10, 11 and
12. The subspaces Cn are defined from the wavelets corresponding to the largest K coefficients of the
discrete wavelet transform obtained for each slice. Figs 9 and 10 compare the original image to the
low pass image and recovered image as well as the correponding original low pass filter and recovered
filter. Figs 11 and 12 provide further highlights on the details that can be recovered through nuclear
norm minimization.

5

Conclusions and perspectives

In this paper we have considered a version of the blind deconvolution problem where the paradigm is shifted
towards an arbitrary transfer function at the expense of requiring a small number of sufficiently distinct
input signals to guarantee the recovery. Theory shows that whenever the number of inputs N , the ambient
dimension L and the dimension of the inputs subspace K satisfy L & K 3/2 µ2h and N & K 1/2 µ2m up to
log factors and for coherences defined as in (11) and (10), both the unknown filter as well as the unknown
input signals can be recovered from the outputs to the filter by means of the nuclear norm minimization
program (9). When trying to recover an unknown filter of length L whose Fourier transform is sufficiently
”spread out” from its outputs, having a collection of a few outputs is therefore helpful. Such a framework
finds applications in medical imaging, astronomy or microscopy where multiple compressible slices of a same
volume are usually acquired.
The importance of the coherences (10) and (11) in the formulation of Theorem 1 illustrates the fact that
the measured outputs have to carry a sufficient amount of information on both the impulse response h and
the input signals. An intuition for this is that if the energy of the filter is concentrated at a single position
in Fourier space, most of the measurements will be zero and wont therefore carry information on the input
signals. An equivalent statement holds whenever a number of input signals are zero. In other words, what
really matters is the number of ”meaningful” input signals and the fact that the probing of the filter through
those signals should provide a sufficient amount of information on this filter. The result of this paper also
1 https://www.mathworks.com/matlabcentral/fileexchange/9416-3d-shepp-logan-phantom
2 https://graphics.stanford.edu/data/voldata/

51

Figure 7: Evolution of the recovery in the blind super-resolution framework, for different sizes of the filter
support (Part I). The original (low resolution) slice from the 3D Shepp Logan volume is shown on the right
and the result obtained through the nuclear norm relaxation (9) is shown on the left. Here K ≈ L/20,
L = 2002 , |supp(h)| = L/10 and the subspaces Cn are defined from the discrete wavelet transform with
Daubechies wavelets by taking the wavelets corresponding to the K largest coefficients of the DWT.

Figure 8: Evolution of the recovery in the blind super-resolution framework, for different sizes of the filter
support (Part I). The original (low resolution) slice from the 3D Shepp Logan volume is shown on the right
and the result obtained through the nuclear norm relaxation (9) is shown on the left. Here K ≈ L/20,
L = 2002 , supp(h) = L/16 and the subspaces Cn are defined from the discrete wavelet transform with
Daubechies wavelets by taking the wavelets corresponding to the K largest coefficients of the DWT.
shows that the notion of coherence alone is sufficient to describe the recovery. This is in contrast to previous
work where both sparsity and sufficient incoherence were required (see for example [1] or [38]). A direct
consequence of this conclusion is that as soon as the support of the Fourier transform of the filter is on the
order of the ambient dimension, it remains possible to recover both the filter and input signals including
when the support is stricly smaller than this dimension. This observation implies that the nuclear norm
minimization program (9) can be used in the framework of blind super-resolution as well.

52

Figure 9: Ideal low pass filtering and blind super-resolution of MR images through nuclear norm minimization. From top to botton, left to right. (Top left) Original ideal low pass filtering of the Magentic resonance
image shown in Fig. 10. (Top right) recovery through semidefinite programming. (Bottom left) Original low
pass (sinc) filter and (bottom right) recovery of the low pass filter. Additional comparison of the quality
of the super-resolved image vs low pass image is given in Fig. 11. The change in contrast appearing in
the low pass image comes from the blurring of the forehead which has the highest intensity in the original
image. Such blurring homogeneize the forehead with the background then resulting in a reduced contrast
and apparent higher intensity of the global image.
To certify recovery through nuclear norm minimization, we construct a certificate of optimality. Exact
recovery can then be shown by proving that the certificate satisfies the conditions derived from requiring the
subgradient of the Lagrangian to vanish at the solution X0 . In this case, the candidate certificate expands as
a Neumann series, and certifying these conditions can be done by controlling each of the terms in the series.
We applied ad-hoc concentration results to the first two terms only as it is assumed that deriving bounds on
a constant number of higher order terms should follow the same approach and is likely to become gradually
heavier in terms of the derivations needed. It is likely that the sample complexity will benefit from further
applications of the decoupling inequality of de la Penã and Giné [20]. It is not clear, however whether such
an inequality can be used to achieve the sample complexity observed empirically (L & K, N & 1) as the
constant appearing in this inequality scales badly with the order of the U -statistics involved.
This work naturally raises an important open question: How far can we go in the complexity of the filter while
still being able to certify the recovery? So far we have assumed that the filter remained constant spatially. In
many applications, however, the point spread function varies with the position and it is not clear whether the
proof techniques used in this paper can translate easily to that framework. Similar open problems include
the more general field of blind linear system identification or even non linear system identification.
Certainly equally interesting is the question of the efficiency of nuclear norm minimization for blind deconvolution. Why does the nuclear norm perform so well on the blind deconvolution problem? The probabilistic
53

Figure 10: From left to right, Slice from the original MR volume, approximation obtained throug ideal low
pass filtering through a truncated gaussian filter similar to the one depicted in Fig 16 and recovered (superresolved) image obtained through the nuclear norm minimization program (9). An additional background
subtraction step has been added to facilitate the comparison. The original (unprocessed) images are shown
in Fig. 9.

Figure 11: . Additional illustration of the use of the nuclear norm minimization program (9) in the framework
of the blind super-resolution of MR images from a same volume. (Part I)

Figure 12: . Additional illustration of the use of the nuclear norm minimization program (9) in the framework
of the blind super-resolution of MR images from a same volume (Part 2).
method reveals an elegant tool to derive recovery guarantees yet it does not make use of the particular structure of the problem, and as a consequence, is unable to explain the efficiency of nuclear norm minimization
for that particular problem structure. Blind deconvolution however seems a natural candidate for a better

54

Original 1

Original 2

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

0

50

100

150

200

−0.4

0

50

100

150

200

150

200

150

200

Low resolution 2

Low resolution 1
0.2

0.2
0.1

0
0
−0.1
−0.2

−0.2
0

50

100

150

200

0

50

Recovered 1

Recovered 2

0.4

0.4

0.2

0.2

0

0

−0.2

−0.2

−0.4

0

50

100

100

150

200

−0.4

0

50

100

Figure 13: Blind super-resolution of a wavelet train from its convolution with an ideal gaussian low pass
filter. The side lobes arising from the convolution with sinc induced by the windowing in Fourier space are
clearly visible on the right figure. The subspaces Cn are defined from the wavelets corresponding to the K
largest coefficients of the discrete wavelet transform of each train. Elimination of the spurious high frequency
oscillations following from the low pass nature of the filter, through the minimization program (9) are further
highlighted in Fig. 14 where the recovered train is superimposed on its low pass approximation.

55

Original + result of blind super-resolution

Original + result of blind super-resolution
0.2

0.2

0

0

−0.2
−0.2
80

100

120

140

160

180

40

60

Original + low pass

80

100

120

140

160

Original + low pass
0.2

0.2

0

0

−0.2
−0.2
80

100

120

140

160

180

20

40

60

80

100

120

140

160

Figure 14: Further illustration of the blind super-resolution of a wavelet train from its convolution with
an ideal gaussian low pass filter. The recovered (super-resolved) wavelet train (red) is superimposed on its
low pass approximation. Note the side lobes arising from the sinc and resulting from the multiplication
in Fourier space with the ideal low pass filter of Fig. 16. The subspaces Cn are defined from the wavelets
corresponding to the K largest coefficients of the discrete wavelet transform of each train.

56

Figure 15: Low resolution images obtained through ideal low pass filtering with (right) and without (left)
the addition of a blurring kernel. In the left image, which is obtained without the blurring kernel, the effect
of the side lobes arising from the sinc is clearly visible whether on the right image, the blurring step helped
mitigrate those effects. Both images result from a windowed acquisition in frequency.

5
4
3
2
1
250
200

250

150

200
150

100
100

50

50

Figure 16: Representation of a gaussian ideal low pass filter such as the one used in the super-resolution
examples of section 4.3.2. All the coefficients above and below the cutoff frequencies ±ω1 and ±ω2 have
been set to zero.

57

Figure 17: Reconstruction of the 3D Shepp-Logan phantom from blurred, low pass sample images. The
original image before undersampling together with the filter and its corresponding transfer function are
shown in Fig. 1 (top row). The filter used here is gaussian and multiplied by a rectangular window zeroing
out all frequencies above ±(kxc , kyc ) such as shown in Fig. 16. The subspaces Cn are defined by considering
the wavelets corresponding to the K largest coefficients obtained from the discrete wavelet transform of each
of the volume slices.

58

Figure 18: Reconstruction of the 3D Shepp-Logan phantom from blurred, low pass samples (continued).
(Right) Original, blurred and low pass images. (Left) Recovery through the nuclear norm minimization
program (9). The gaussian ideal low pass filter has the same structure as shown in Fig. 16 and the subspaces
Cn are defined by considering the wavelets corresponding to the K largest coeffcients obtained from the
discrete wavelet transform of each of the volume slices.

59

understanding of the propagation of information in semidefinite relaxations such as discussed in [18], in the
framework of matrix completion.

References
[1] A. Ahmed and L. Demanet. Leveraging diversity and sparsity in blind deconvolution. arXiv preprint
arXiv:1610.06098, 2016.
[2] A. Ahmed, B. Recht, and J. Romberg. Blind deconvolution using convex programming. Information
Theory, IEEE Transactions on, 60(3):1711–1732, 2014.
[3] A. Ahmed and J. Romberg.
arXiv:1308.5146, 2013.

Compressive multiplexing of correlated signals.

arXiv preprint

[4] A. Ahmed and J. Romberg. Compressive multiplexing of correlated signals. IEEE Trans. Inform.
Theory, 1:479–498, 2015.
[5] A. Ahmed and J. Romberg. Compressive sampling of ensembles of correlated signals. arXiv preprint
arXiv:1501.06654, 2015.
[6] G. Ayers and J. C. Dainty. Iterative blind deconvolution method and its applications. Optics letters,
13(7):547–549, 1988.
[7] M. Bertero and P. Boccacci. Introduction to inverse problems in imaging. CRC press, 1998.
[8] C. Bilen, G. Puy, R. Gribonval, and L. Daudet. Convex optimization approaches for blind sensor
calibration using sparsity. Signal Processing, IEEE Transactions on, 62(18):4847–4856, 2014.
[9] S. Burer and R. Monteiro. A nonlinear programming algorithm for solving semidefinite programs via
low-rank factorization. Mathematical Programming, 95(2):329–357, 2003.
[10] V. Cambareri and L. Jacques. A non-convex blind calibration method for randomised sensing strategies.
arXiv preprint arXiv:1605.02615, 2016.
[11] P. Campisi and K. Egiazarian. Blind image deconvolution: theory and applications. CRC press, 2007.
[12] E. J. Candès and C. Fernandez-Granda. Towards a mathematical theory of super-resolution. Communications on Pure and Applied Mathematics, 67(6):906–956, 2014.
[13] E. J. Candes, X. Li, and M. Soltanolkotabi. Phase retrieval via wirtinger flow: Theory and algorithms.
IEEE Trans. Inform. Theory, 61(4):1985–2007, 2015.
[14] E. J. Candès and B. Recht. Exact matrix completion via convex optimization. Foundations of Computational mathematics, 9(6):717–772, 2009.
[15] M. Cannon. Blind deconvolution of spatially invariant image blurs with phase. Acoustics, Speech and
Signal Processing, IEEE Transactions on, 24(1):58–63, 1976.
[16] T. F. Chan and C.-K. Wong. Total variation blind deconvolution. Image Processing, IEEE Transactions
on, 7(3):370–375, 1998.
[17] A. Cosse. A note on the blind deconvolution of multiple sparse signals from unknown subspaces.
In Wavelets and Sparsity XVII, volume 10394, page 103941N. International Society for Optics and
Photonics, 2017.
[18] A. Cosse and L. Demanet. Rank-one matrix completion is solved by the sum-of-squares relaxation
of order two. In Proceedings of the 6th IEEE International Workshop on Computational Advances in
Multi-Sensor Adaptive Processing (CAMSAP’15). IEEE, 2015.

60

[19] V. De la Pena and E. Giné. Decoupling: from dependence to independence. Springer Science &amp;
Business Media, 2012.
[20] V. H. de la Peña and S. J. Montgomery-Smith. Decoupling inequalities for the tail probabilities of
multivariate u-statistics. The Annals of Probability, pages 806–816, 1995.
[21] L. Demanet and P. Hand. Stable optimizationless recovery from phaseless linear measurements. J.
Fourier Anal. Applications, 20(1):199–221, 2014.
[22] A. V. der Vaart and J. Wellner. Weak Convergence and Empirical Processes. Springer, 1996.
[23] M. Fazel. Matrix rank minimization with applications. PhD thesis, Stanford University, March 2002.
[24] M. Fazel, H. Hindi, and S. P. Boyd. A rank minimization heuristic with application to minimum order
system approximation. In American Control Conference, 2001. Proceedings of the 2001, volume 6, pages
4734–4739. IEEE, 2001.
[25] E. Giné, D. M. Mason, and J. A. Wellner. High Dimensional Probability II, volume 47. Springer Science
&amp; Business Media, 2012.
[26] R. Gribonval, G. Chardon, and L. Daudet. Blind calibration for compressed sensing by convex optimization. In Acoustics, Speech and Signal Processing (ICASSP), 2012 IEEE International Conference
on, pages 2713–2716. IEEE, 2012.
[27] D. Gross. Recovering low-rank matrices from few coefficients in any basis. Information Theory, IEEE
Transactions on, 57(3):1548–1566, 2011.
[28] W. Hoeffding. A class of statistics with asymptotically normal distribution. The annals of mathematical
statistics, pages 293–325, 1948.
[29] M. Kech and F. Krahmer. Optimal injectivity conditions for bilinear inverse problems with applications
to identifiability of deconvolution problems. arXiv preprint arXiv:1603.07316, 2016.
[30] V. Koltchinskii et al. A remark on low rank matrix recovery and noncommutative bernstein type
inequalities. In From Probability to Statistics and Back: High-Dimensional Models and Processes–A
Festschrift in Honor of Jon A. Wellner, pages 213–226. Institute of Mathematical Statistics, 2013.
[31] V. Koltchinskii, K. Lounici, and A. Tsybakov. Nuclear-norm penalization and optimal rates for noisy
low-rank matrix completion. Ann. Stat., 39(5):2302–2329, 2011.
[32] C. S. Kubrusly. Spectral theory of operators on Hilbert spaces. Springer Science &amp; Business Media,
2012.
[33] R. L. Lagendijk, J. Biemond, and D. E. Boekee. Blur identification using the expectation-maximization
algorithm. In Acoustics, Speech, and Signal Processing, 1989. ICASSP-89., 1989 International Conference on, pages 1397–1400. IEEE, 1989.
[34] A. Levin, Y. Weiss, F. Durand, and W. T. Freeman. Understanding and evaluating blind deconvolution
algorithms. In Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on,
pages 1964–1971. IEEE, 2009.
[35] X. Li, S. Ling, T. Strohmer, and K. Wei. Rapid, robust, and reliable blind deconvolution via nonconvex
optimization. arXiv preprint arXiv:1606.04933, 2016.
[36] Y. Li, K. Lee, and Y. Bresler. Identifiability in blind deconvolution with subspace or sparsity constraints.
IEEE Transactions on Information Theory, 62(7):4266–4275, 2016.
[37] S. Ling and T. Strohmer. Blind deconvolution meets blind demixing: Algorithms and performance
bounds. arXiv preprint arXiv:1512.07730, 2015.

61

[38] S. Ling and T. Strohmer.
arXiv:1501.06864, 2015.

Self-calibration and biconvex compressive sensing.

arXiv preprint

[39] S. Ling and T. Strohmer. Self-calibration via linear least squares. arXiv preprint arXiv:1611.04196,
2016.
[40] B. Recht. A simpler approach to matrix completion. J. Mach. Learn. Res., 12(12):3413–3430, December
2011.
[41] B. Recht, M. Fazel, and P. Parrilo. Guaranteed minimum-rank solutions of linear matrix equations via
nuclear norm minimization. SIAM Review, 52(3):471–501, 2010.
[42] B. Recht, M. Fazel, and P. A. Parrilo. Guaranteed minimum-rank solutions of linear matrix equations
via nuclear norm minimization. SIAM review, 52(3):471–501, 2010.
[43] S. J. Reeves and R. M. Mersereau. Blur identification by the method of generalized cross-validation.
Image Processing, IEEE Transactions on, 1(3):301–311, 1992.
[44] C. A. Segall, R. Molina, and A. K. Katsaggelos. High-resolution images from low-resolution compressed
video. Signal Processing Magazine, IEEE, 20(3):37–48, 2003.
[45] T. G. Stockham Jr, T. M. Cannon, and R. B. Ingebretsen. Blind deconvolution through digital signal
processing. Proceedings of the IEEE, 63(4):678–692, 1975.
[46] J. A. Tropp. An introduction to matrix concentration inequalities. arXiv preprint arXiv:1501.01571,
2015.
[47] J. A. Tropp. The expected norm of a sum of independent random matrices: An elementary approach.
In High Dimensional Probability VII, pages 173–202. Springer, 2016.
[48] D. Tse and P. Viswanath. Fundamentals of wireless communication. Cambridge university press, 2005.
[49] A. W. Van Der Vaart and J. A. Wellner. Weak convergence. In Weak Convergence and Empirical
Processes, pages 16–28. Springer, 1996.
[50] R. Vershynin. Introduction to the non-asymptotic analysis of random matrices.
arXiv:1011.3027, 2010.

arXiv preprint

[51] R. Vershynin. A simple decoupling inequality in probability theory. preprint, 2011.
[52] J. Wainwright. High-dimensional statistics: A non-asymptotic viewpoint. preparation. University of
California, Berkeley, 2015.
[53] D. Yang, G. Tang, and M. B. Wakin. Non-stationary blind super-resolution. In Acoustics, Speech and
Signal Processing (ICASSP), 2016 IEEE International Conference on, pages 4727–4731. IEEE, 2016.
[54] D. Yang, G. Tang, and M. B. Wakin. Super-resolution of complex exponentials from modulations with
unknown waveforms. IEEE Transactions on Information Theory, 62(10):5809–5830, 2016.
[55] Y.-L. You and M. Kaveh. Blind image restoration by anisotropic regularization. Image Processing,
IEEE Transactions on, 8(3):396–407, 1999.

62

