arXiv:1611.04196v3 [] 21 Nov 2017

Self-Calibration and Bilinear Inverse Problems
via Linear Least Squares∗
Shuyang Ling†, Thomas Strohmer‡
November 22, 2017

Abstract
Whenever we use devices to take measurements, calibration is indispensable. While the
purpose of calibration is to reduce bias and uncertainty in the measurements, it can be quite
difficult, expensive, and sometimes even impossible to implement. We study a challenging
problem called self-calibration, i.e., the task of designing an algorithm for devices so that the
algorithm is able to perform calibration automatically. More precisely, we consider the setup
y = A(d)x+ε where only partial information about the sensing matrix A(d) is known and where
A(d) linearly depends on d. The goal is to estimate the calibration parameter d (resolve the
uncertainty in the sensing process) and the signal/object of interests x simultaneously. For three
different models of practical relevance, we show how such a bilinear inverse problem, including
blind deconvolution as an important example, can be solved via a simple linear least squares
approach. As a consequence, the proposed algorithms are numerically extremely efficient, thus
potentially allowing for real-time deployment. We also present a variation of the least squares
approach, which leads to a spectral method, where the solution to the bilinear inverse problem
can be found by computing the singular vector associated with the smallest singular value of a
certain matrix derived from the bilinear system. Explicit theoretical guarantees and stability
theory are derived for both techniques; and the number of sampling complexity is nearly optimal
(up to a poly-log factor). Applications in imaging sciences and signal processing are discussed
and numerical simulations are presented to demonstrate the effectiveness and efficiency of our
approach.

1

Introduction

Calibration is ubiquitous in all fields of science and engineering. It is an essential step to guarantee
that the devices measure accurately what scientists and engineers want. If sensor devices are not
properly calibrated, their measurements are likely of little use to the application. While calibration
is mostly done by specialists, it often can be expensive, time-consuming and sometimes even impossible to do in practice. Hence, one may wonder whether it is possible to enable machines to calibrate
themselves automatically with a smart algorithm and give the desired measurements. This leads
to the challenging field of self-calibration (or blind calibration). It has a long history in imaging
sciences, such as camera self-calibration [33, 22], blind image deconvolution [12], self-calibration in
∗

This research was supported by the NSF via Award Nr. dtra-dms 1322393 and Award Nr. DMS 1620455.
Courant Institute of Mathematical Sciences and the Center for Data Science, New York University, NY 10003
(Email: sling@cims.nyu.edu)
‡
Department of Mathematics, University of California Davis, CA 95616 (Email: strohmer@math.ucdavis.edu).
†

1

medical imaging [35], and the well-known phase retrieval problem (phase calibration) [16]. It also
plays an important role in signal processing [18] and wireless communications [38, 32].
Self-calibration is not only a challenging problem for engineers, but also for mathematicians. It
means that one needs to estimate the calibration parameter of the devices to adjust the measurements as well as recover the signal of interests. More precisely, many self-calibration problems are
expressed in the following mathematical form,
y = A(d)x + ε,

(1.1)

where y is the observation, A(d) is a partially unknown sensing matrix, which depends on an unknown parameter d and x is the desired signal. An uncalibrated sensor/device directly corresponds
to “imperfect sensing”, i.e., uncertainty exists within the sensing procedure and we do not know
everything about A(d) due to the lack of calibration. The purpose of self-calibration is to resolve
the uncertainty i.e., to estimate d in A(d) and to recover the signal x at the same time.
The general model (1.1) is too hard to get meaningful solutions without any further assumption
since there are many variants of the general model under different settings. In (1.1), A(d) may
depend on d in a nonlinear way, e.g., d can be the unknown orientation of a protein molecule
and x is the desired object [44]; in phase retrieval, d is the unknown phase information of the
Fourier transform of the object [16]; in direction-of-arrival estimation d represents unknown offset,
gain, and phase of the sensors [47]. Hence, it is impossible to resolve every issue in this field,
but we want to understand several scenarios of self-calibration which have great potential in real
world applications. Among all the cases of interest, we assume that A(d) linearly depends on the
unknown d and will explore three different types of self-calibration models that are of considerable
practical relevance. However, even for linear dependence, the problem is already quite challenging,
since in fact we are dealing with bilinear (nonlinear) inverse problems. All those three models have
wide applications in imaging sciences, signal processing, wireless communications, etc., which will
be addressed later. Common to these applications is the desire or need for fast algorithms, which
ideally should be accompanied by theoretical performance guarantees. We will show under certain
cases, these bilinear problems can be solved by linear least squares exactly and efficiently if no
noise exists, which is guaranteed by rigorous mathematical proofs. Moreover, we prove that the
solution is also robust to noise with tools from random matrix theory. Furthermore, we show that
a variation of our approach leads to a spectral method, where the solution to the bilinear problem
can be found by computing the singular vector associated with the smallest singular matrix of a
certain matrix derived from the bilinear system.

1.1

State of the art

By assuming that A(d) linearly depends on d, (1.1) becomes a bilinear inverse problem, i.e., we
want to estimate d and x from y, where y is the output of a bilinear map from (d, x). Bilinear
inverse problems, due to its importance, are getting more and more attentions over the last few
years. On the other hand, they are also notoriously difficult to solve in general. Bilinear inverse
problems are closely related to low-rank matrix recovery, see [15] for a comprehensive review. There
exists extensive literature on this topic and it is not possible do justice to all these contributions.
Instead we will only highlight some of the works which have inspired us.
Blind deconvolution might be one of the most important examples of bilinear inverse problems [12], i.e., recovering f and g from y = f ∗ g, where “∗” stands for convolution. If both
f and g are inside known low-dimensional subspaces, the blind deconvolution can be rewritten
2

as F(y) = diag(Bd)Ax, where F(f ) = Bd, F(g) = Ax and “F” denotes the Fourier transform. In the inspiring work [4], Ahmed, Romberg and Recht apply the “lifting” techniques [13]
and convert the problem into estimation of the rank-one matrix dx∗ . It is shown that solving a
convex relaxation enables recovery of dx∗ under certain choices of B and A. Following a similar
spirit, [30] uses “lifting” combined with a convex approach to solve the scenarios with sparse x
and [29] studies the so called “blind deconvolution and blind demixing” problem. The other line of
blind deconvolution follows a nonconvex optimization approach [3, 26, 25]. In [3], Ahmed, Romberg
and Krahmer, using tools from generic chaining, obtain local convergence of a sparse power factorization algorithm to solve this blind deconvolution problem when h and x are sparse and B and A
are Gaussian random matrices. Under the same setting as [3], Lee et al. [25] propose a projected
gradient descent algorithm based on matrix factorizations and provide a convergence analysis to
recover sparse signals from subsampled convolution. However, this projection step can be hard to
implement. As an alternative, the expensive projection step is replaced by a heuristic approximate
projection, but then the global convergence is not fully guaranteed. Both [3, 25] achieve nearly
optimal sampling complexity. [26] proves global convergence of a gradient descent type algorithm
when B is a deterministic Fourier type matrix and A is Gaussian. Results about identifiability
issue of bilinear inverse problems can be found in [27, 23, 28].
Another example of self-calibration focuses on the setup yl = DAl x, where D = diag(d).
The difference from the previous model consists in replacing the subspace assumption by multiple
measurements. There are two main applications of this model. One application deals with blind
deconvolution in an imaging system which uses randomly coded masks [5, 37]. The measurements
are obtained by (subsampled) convolution of an unknown blurring function D with several random
binary modulations of one image. Both [5] and [2] developed convex relaxing approaches (nuclear
norm minimization) to achieve exact recovery of the signals and the blurring function. The other
application is concerned with calibration of the unknown gains and phases D and recovery of the
signal x, see e.g. [10, 9]. Cambareri and Jacques propose a gradient descent type algorithm in [10, 11]
and show convergence of the iterates by first constructing a proper initial guess. An empirical
study is given in [9] when x is sparse by applying an alternating hard thresholding algorithm.
Recently, [1, 2] study the blind deconvolution when inputs are changing. More precisely, the authors
consider yl = f ∗ gl where each gl belongs to a different known subspace, i.e., gl = Cl xl . They
employ a similar convex approach as in [4] to achieve exact recovery with number of measurements
close to the information theoretic limit.
An even more difficult, and from a practical viewpoint highly relevant, scenario focuses on selfcalibration from multiple snapshots [47]. Here, one wishes to recover the unknown gains/phases
D = diag(d) and a signal matrix X = [x1 , · · · , xp ] from Y = DAX. For this model, the
sensing matrix A is fixed throughout the sensing process and one measures output under different
snapshots {xl }pl=1 . One wants to understand under what conditions we can identify D and {xl }pl=1
jointly. If A is a Fourier type matrix, this model has applications in both image restoration from
multiple filters [21] and also network calibration [6, 31]. We especially benefitted from work by
Gribonval and coauthors [20, 8], as well as by Balzano and Nowak [6, 7]. The papers [6, 7] study
the noiseless version of the problem by solving a linear system and [31] takes a total least squares
approach in order to obtain empirically robust recovery in the presence of noise. If each xl is
sparse, this model becomes more difficult and Gribonval et al. [8, 20] give a thorough numerical
study. Very recently, [43] gave a theoretic result under certain conditions. This calibration problem
is viewed as a special case of the dictionary learning problem where the underlying dictionary DA

3

possesses some additional structure. The idea of transforming a blind deconvolution problem into a
linear problem can also be found in [32], where the authors analyze a certain non-coherent wireless
communication scenario.

1.2

Our contributions

In our work, we consider three different models of self-calibration, namely, yl = DAl x, yl = DAl xl
and yl = DAxl . Detailed descriptions of these models are given in the next Section. We do not
impose any sparsity constraints on x or xl . We want to find out xl (or x) and D when yl (or y) and
Al (or A) are given. Roughly, they correspond to the models in [5, 1, 8] respectively. Though all of
the three models belong to the class of bilinear inverse problems, we will prove that simply solving
linear least squares will give solutions to all those models exactly and robustly for invertible D and
for several useful choices of A and Al . Moreover, the sampling complexity is nearly optimal (up to
poly-log factors) with respect to the information theoretic limit (degree of freedom of unknowns).
As mentioned before, our approach is largely inspired by [8] and [6, 7]; there the authors convert
a bilinear inverse problem into a linear problem via a proper transformation. We follow a similar
approach in our paper. The paper [8] provides an extensive empirical study, but no theoretical analysis. Nowak and Balzano, in [6, 7] provide numerical simulations as well as theoretical conditions
on the number of measurements required to solve the noiseless case.
Our paper goes an important step further: On the one hand we consider more general selfcalibration settings. And on the other hand we provide a rigorous theoretical analysis for recoverability and, perhaps most importantly, stability theory in the presence of measurement errors.
Owing to the simplicity of our approach and the structural properties of the underlying matrices, our framework yields self-calibration algorithms that are numerically extremely efficient, thus
potentially allowing for deployment in applications where real-time self-calibration is needed.

1.3

Notation and Outline

We introduce notation which will be used throughout the paper. Matrices are denoted in boldface
or a calligraphic font such as Z and Z; vectors are denoted by boldface lower case letters, e.g. z.
The individual entries of a matrix or a vector are denoted in normal font such as Zij or zi . For any
matrix Z, kZk denotes its operator
qP norm, i.e., the largest singular value, and kZkF denotes its the
2
Frobenius norm, i.e., kZkF =
ij |Zij | . For any vector z, kzk denotes its Euclidean norm. For

both matrices and vectors, Z T and z T stand for the transpose of Z and z respectively while Z ∗
and z ∗ denote their complex conjugate transpose. For any real number z, we let z+ = 12 (z + |z|).
We equip the matrix space CK×N with the inner product defined by hU , V i := Tr(U ∗ V ). A special
case is the inner product of two vectors, i.e., hu, vi = Tr(u∗ v) = u∗ v. We define the correlation
u∗ v
. For a given vector v, diag(v) represents the
between two vectors u and v as Corr(u, v) = kukkvk
diagonal matrix whose diagonal entries are given by the vector v.
C is an absolute constant and Cγ is a constant which depends linearly on γ, but on no other
parameters. In and 1n always denote the n × n identity matrix and a column vector of “1” in
Rn respectively. And {ei }m
el }pl=1 stand for the standard orthonormal basis in Rm and Rp
i=1 and {e
respectively. “∗” is the circular convolution and “⊗” is the Kronecker product.

The paper is organized as follows. The more detailed discussion of the models under consideration and the proposed method will be given in Section 2. Section 3 presents the main results
4

of our paper and we will give numerical simulations in Section 4. Section 5 contains the proof for
each scenario. We collect some useful auxiliary results in the Appendix.

2

Problem setup: Three self-calibration models

This section is devoted to describing three different models for self-calibration in detail. We will also
explain how those bilinear inverse problems are reformulated and solved via linear least squares.

2.1

Three special models of self-calibration

Self-calibration via repeated measurements Suppose we are seeking for information with
respect to an unknown signal x0 with several randomized linear sensing designs. Throughout this
procedure, the calibration parameter D remains the same for each sensing procedure. How can
we recover the signal x0 and D simultaneously? Let us make it more concrete by introducing the
following model,
yl = DAl x0 + εl , 1 ≤ l ≤ p
(2.1)

where D = diag(d) ∈ Cm×m is a diagonal matrix and each Al ∈ Cm×n is a measurement matrix.
Here yl and Al are given while D and x0 are unknown. For simplicity, we refer to the setup (2.1)
as “self-calibration from repeated measurements”. This model has various applications in selfcalibration for imaging systems [10, 11], networks [6], as well as in blind deconvolution from random
masks [5, 37].
Blind deconvolution via diverse inputs Suppose that one sends several different signals
through the same unknown channel, and each signal is encoded differently. Namely, we are considering
yl = f ∗ Cl xl + εl , 1 ≤ l ≤ p.

How can one estimate the channel and each signal jointly? In the frequency domain, this “blind
deconvolution via diverse inputs” [1, 2] problem can be written as (with a bit abuse of notation),
yl = DAl xl + εl ,

1≤l≤p

(2.2)

where D = diag(d) ∈ Cm×m and Al ∈ Cm×n are the Fourier transform of f and Cl respectively.
We aim to recover {xl }pl=1 and D from {yl , Al }pl=1 .
Self-calibration from multiple snapshots Suppose we take measurements of several signals
{xl }pl=1 with the same set of design matrix DA (i.e., each sensor corresponds one row of A and has
an unknown complex-valued calibration term di ). When and how can we recover D and {xl }pl=1
simultaneously? More precisely, we consider the following model of self-calibration model from
multiple snapshots:
yl = DAxl + εl , 1 ≤ l ≤ p.
(2.3)

Here D = diag(d) is an unknown diagonal matrix, A ∈ Cm×n is a sensing matrix, {xl }pl=1 are
n × 1 unknown signals and {yl }pl=1 are their corresponding observations. This multiple snapshots
model has been used in image restoration from multiple filters [21] and self-calibration model for
sensors [47, 6, 31, 20, 8].
5

2.2

Linear least squares approach

Throughout our discussion, we assume that D is invertible, and we let S := diag(s) = D−1 . Here,
D = diag(d) stands for the calibration factors of the sensor(s) [6, 8] and hence it is reasonable to
assume invertibility of D. For, if a sensor’s gain were equal to zero, then it would not contribute
any measurements to the observable y, in which case the associated entry of y would be zero.
But then we could simply discard that entry and consider the correspondingly reduced system of
equations, for which the associated D is now invertible.
One simple solution is to minimize a nonlinear least squares objective function. Let us take (2.1)
as an example (the others (2.2) and (2.3) have quite similar formulations),
min
D,x

p
X
l=1

kDAl x − yl k2 .

(2.4)

The obvious difficulty lies in the biconvexity of (2.4), i.e., if either D or x is fixed, minimizing over
the other variable is a convex program. In general, there is no way to guarantee that any gradient
descent algorithm/alternating minimization will give the global minimum. However, for the three
models described above, there is one shortcut towards the exact and robust recovery of the solution
via linear least squares if D is invertible.
We continue with (2.1) when εl = 0, i.e.,
diag(yl )s = Al x,

1≤l≤p

(2.5)

where Syl = diag(yl )s with si = d−1
and S is defined as diag(s). The original measurement
i
equation turns out to be a linear system with unknown s and x. The same idea of linearization
can be also found in [20, 8, 43, 6]. In this way, the ground truth z0 := (s0 , x0 ) lies actually inside
the null space of this linear system.
Two issues arise immediately: One the one hand, we need to make sure that (s0 , x0 ) spans
the whole null space of this linear system. This is equivalent to the identifiability issue of bilinear
problems of the form (2.4), because if the pair (α−1 d0 , αx0 ) for some α 6= 0 is (up to the scalar
α) unique solution to (2.1), then (αs0 , αx0 ) spans the null space of (2.5), see also [27, 23, 28]. On
the other hand, we also need to avoid the trivial scenario (s0 , x0 ) = (0, 0), since it has no physical
meaning. To resolve the latter issue, we add the extra linear constraint (see also [8, 6])
  
s
w,
= c,
(2.6)
x
where the scalar c can be any nonzero number (we note that w should of course not be orthogonal
to the solution). Therefore, we hope that in the noiseless case it suffices to solve the following linear
system to recover (d, x0 ) up to a scalar, i.e.,


diag(y1 ) −A1
 
 

..
.. 
0
s


.
.
(2.7)
=


c (mp+1)×1
x (m+n)×1
diag(yp ) −Ap 
{z
} |
{z
}
|
w∗
z
(mp+1)×(m+n)
b
{z
}
|
Aw

6

In the presence of additive noise, we replace the linear system above by a linear least squares
problem
p
X
minimize
k diag(yl )s − Al xk2 + |w ∗ z − c|2
l=1

with respect to s and x, or equivalently,

minimize kAw z − bk2
z

(2.8)

 
 
s
0
where z =
,b=
, and Aw is the matrix on the left hand side of (2.7). Following the same
x
c
idea, (2.2) and (2.3) can also be reformulated into linear systems and be solved via linear least
squares. The matrix Aw and the vector z take a slightly different form for those cases, see (3.1)
and (3.3), respectively.
Remark 2.1. Note that solving (2.8) may not be the optimal choice to recover the unknowns
from the perspective of statistics since the noisy perturbation actually enters into Aw instead of b.
More precisely, the noisy perturbation δA to the left hand side of the corresponding linear system
for (2.1), (2.2) and (2.3), is always in the form of


diag(ε1 ) 0

.. 
..

.
.
(2.9)
δA := 
.
diag(εp ) 0
0
0
The size of δA depends on the models. Hence total least squares [31] could be a better alternative
while it is more difficult to analyze and significantly more costly to compute. Since computational
efficiency is essential for many practical applications, a straightforward implementation of total least
squares is of limited use. Instead one should keep in mind that the actual perturbation enters only
into diag(yl ), while the other matrix blocks remain unperturbed. Constructing a total least squares
solution that obeys these constraints, doing so in a numerically efficient manner and providing
theoretical error bounds for it, is a rather challenging task, which we plan to address in our future
work.

Remark 2.2. Numerical simulations imply that the performance under noisy measurements depends on the choice of w, especially how much w and z0 are correlated. One extreme case is that
hw, z0 i = 0, in which case we cannot avoid the solution z = 0. It might be better to add a constraint
like kwk = 1. However, this will lead to a nonlinear problem which may not be solved efficiently
and not come with rigorous recovery guarantees. Therefore, we present an alternative approach in
the next subsection.

2.3

Spectral method

In this subsection, we discuss a method for solving the self-calibration problem, whose performance
does not depend on the choice of w as it avoids the need of w in the first place. Let S be the

7

matrix Aw excluding the last row (the one which contains w). We decompose S into S = S0 + δS
where S0 is the noiseless part of S and δS is the noisy part1 .
We start with the noise-free scenario: if δS = 0, there holds S = S0 and the right singular
vector of S corresponding to the smallest singular value is actually z0 . Therefore we can recover z0
by solving the following optimization problem:
min kSzk.

kzk=1

(2.10)

Obviously, its solution is equivalent to the smallest singular value of S and its corresponding singular
vector.
If noise exists, the performance will depend on how large the second smallest singular value of
S0 is and on the amount of noise, given by kδSk. We will discuss the corresponding theory and
algorithms in Section 3.4, and the proof in Section 5.4.

3

Theoretical results

We present our theoretical findings for the three models (2.1), (2.2) and (2.3) respectively for
different choices of Al or A. In one of our choices the Al are Gaussian random matrices. The
rationale for this choice is that, while a Gaussian random matrix is not useful or feasible in most
applications, it often provides a benchmark for theoretical guarantees and numerical performance.
Our other choices for the sensing matrices are structured random matrices, such as e.g. the product
of a deterministic partial (or a randomly subsampled) Fourier matrix or a Hadamard matrix2 with a
diagonal binary random matrix3 . These matrices bring us closer to what we encounter in real world
applications. Indeed, structured random matrices of this type have been deployed for instance in
imaging and wireless communications, see e.g. [19, 41].
By solving simple variations (for different models) of (2.8), we can guarantee that the ground
truth is recovered exactly up to a scalar if no noise exists and robustly if noise is present. The
number of measurements required for exact and robust recovery is nearly optimal, i.e., close to the
information-theoretic limit up to a poly-log factor. However, the error bound for robust recovery
is not optimal. It is worth mentioning that once the signals and calibration parameter D are
identifiable, we are able to recover both of them exactly in absence of noise by simply solving a
linear system. However, identifiability alone cannot guarantee robustness.
Throughout this section, we let dmax := max1≤i≤m |di,0 | and dmin := min1≤i≤m |di,0 | where
{di,0 }m
i=1 are the entries of the ground truth d0 . We also define Aw,0 as the noiseless part of Aw
for each individual model.

3.1

Self-calibration via repeated measurements

For model (2.1) we will focus on three cases:
1

This is a slight abuse of notation, since the δA defined in (2.9) has an additional row with zeros. However, it will
be clear from the context which δA we refer to, and more importantly, in our estimates we mainly care about kδAk
which coincides for both choices of δA.
2
A randomly subsampled Fourier matrix is one, where we randomly choose a certain number of rows or columns
of the Discrete Fourier Transform matrix, and analogously for a randomly subsampled Hadamard matrix
3
At this point, we are not able to prove competitive results for fully deterministic sensing matrices.

8

(a) Al is an m × n complex Gaussian random matrix, i.e., each entry in Al is given by
√i N (0, 1).
2

√1 N (0, 1) +
2

(b) Al is an m × n “tall” random DFT/Hadamard matrix with m ≥ n, i.e., Al := HMl where H
consists of the first n columns of an m × m DFT/Hadamard matrix and each Ml := diag(ml )
is a diagonal matrix with entries taking on the value ±1 with equal probability. In particular,
there holds,
A∗l Al = Ml∗ H ∗ HMl = mIn .
(c) Al is an m × n “fat” random partial DFT matrix with m < n, i.e., Al := HMl where H
consists of m columns of an n × n DFT/Hadamard matrix and each Ml := diag(ml ) is a
diagonal matrix, which is defined the same as case (b),
Al A∗l = HH ∗ = nIm .
Our main findings are summarized as follows:
Theorem 3.1. Consider the self-calibration model given in (2.1), where Aw is as in (2.7) and
Aw,0 is the noiseless part of Aw . Then, for the solution ẑ of (2.8) and α = w∗cz0 , there holds


2
kẑ − αz0 k
≤ κ(Aw,0 )η 1 +
kαz0 k
1 − κ(Aw,0 )η
if κ(Aw,0 )η < 1 where η =

2kδAk
√
mp .

κ(Aw,0 ) ≤

s

The condition number of Aw satisfies

max{d2max kxk2 , m}
6(mp + kwk2 )
,
min{mp, kwk2 | Corr(w, z0 )|2 } min{d2min kxk2 , m}

√
mp
s
√
2 3
max{d2max kxk2 , m}
κ(Aw,0 ) ≤
,
| Corr(w, z0 )| min{d2min kxk2 , m}
 n
log2 (m + n);
(a) with probability 1 − (m + n)−γ if Al is Gaussian and p ≥ c0 γ max 1, m
where Corr(w, z0 ) =

w ∗ z0
kwkkz0 k ,

and for kwk =

(b) with probability 1 − (m + n)−γ − 2(mp)−γ+1 if each Al is a “tall” (m × n, m ≥ n) random
Hadamard/DFT matrix and p ≥ c0 γ 2 log(m + n) log(mp);

(c) with probability 1 − (m + n)−γ − 2(mp)−γ+1 if each Al is a “fat” (m × n, m ≤ n) random
Hadamard/DFT matrix and mp ≥ c0 γ 2 n log(m + n) log(mp).
Remark 3.2. Our result is nearly optimal in terms of required number of measurements, because
the number of constraints mp is required to be slightly greater than n + m, the number of unknowns.
Theorem 3.1 can be regarded as a generalized result of [10, 11], in which D is assumed to be positive
and Al is Gaussian. In our result, we only need D to be an invertible complex diagonal matrix and
Al can be a Gaussian or random Fourier type matrix. The approaches are quite different, i.e., [10]
essentially uses nonconvex optimization by first constructing a good initial guess and then applying
gradient descent to recover D and x. Our result also provides a provable fast alternative algorithm
to “the blind deconvolution via random masks” in [5, 37] where a SDP-based approach is proposed.
9

3.2

Blind deconvolution via diverse inputs

We now analyze model (2.2). Following similar steps that led us from (2.1) to (2.7), it is easy to
see that the linear system associated with (2.2) is given by

 

diag(y1 ) −A1
0
···
0
s
diag(y2 )



0
−A2 · · ·
0 
 

 x1 
0




..
..
..
.
.
x
..
.. 
=
.
(3.1)

 2
.
.
.
c (mp+1)×1

 .. 

diag(yp )



.
0
0
· · · −Ap
x
w∗
p (m+n)×1
(mp+1)×(np+m)
{z
}|
{z
}
|
z

Aw

We consider two scenarios:

(a) Al is an m × n complex Gaussian random matrix, i.e., each entry in Al yields
√i N (0, 1).
2

√1 N (0, 1)
2

+

(b) Al is of the form
Al = Hl Ml ,

(3.2)

where Hl ∈ Cm×n is a random partial Hadamard/Fourier matrix, i.e., the columns of Hl
are uniformly sampled without replacement from an m × m DFT/Hadamard matrix; Ml :=
diag(ml ) = diag(ml,1 , · · · , ml,n ) is a diagonal matrix with {ml,i }ni=1 being i.i.d. Bernoulli
random variables.
Theorem 3.3. Consider the self-calibration model given in (2.2), where Aw is as in (3.1). Let
xmin := min1≤l≤p kxl k and xmax := max1≤l≤p kxl k. Then, for the solution ẑ of (2.8) and α = w∗cz0 ,
there holds


2
kẑ − αz0 k
≤ κ(Aw,0 )η 1 +
kαz0 k
1 − κ(Aw,0 )η
if κ(Aw,0 )η < 1 where η =

κ(Aw,0
and for kwk =

√

2kδAk
√
.
m

v
u
u
)≤t

The condition number of Aw,0 obeys

max{pd2max , x2m }
6x2max (m + kwk2 )
min
,
x2min min{m, kwk2 | Corr(w, z0 )|2 } min{pd2min , x2m }
max

m, there holds
v
u
√
u max{pd2max , x2m }
2 3xmax
min
t
κ(Aw,0 ) ≤
,
xmin | Corr(w, z0 )| min{pd2min , x2m }
max

(a) with probability at least 1 − (np + m)−γ if Al is an m × n (m > n) complex Gaussian random
matrix and


n
1
1
+
(γ + 1) log 2 (np + m) ≤ .
C0
p m
4

10

(b) with probability at least 1 − (np + m)−γ if Al yields (3.2) and


n−1
1
1
+
γ 3 log4 (np + m) ≤ ,
C0
p m−1
4

m ≥ 2.

Remark 3.4. Note that if δA = 0, i.e., in the noiseless case, we have z = αz0 if mp ≥ (np +
m)poly(log(np + m)). Here mp is the number of constraints and np + m is the degree of freedom.
Therefore, our result is nearly optimal in terms of information theoretic limit. Compared with a
similar setup in [1], we have a more efficient algorithm since [1] uses nuclear norm minimization
to achieve exact recovery. However, the assumptions are slightly different, i.e., we assume that D
is invertible and hence the result depends on D while [1] imposes “incoherence” on d by requiring
kF (d)k∞
relatively small, where F denotes Fourier transform.
kdk

3.3

Self-calibration from multiple snapshots

We again follow a by now familiar procedure to derive the linear system associated with (2.3),
which turns out to be
 


s
diag(y1 ) −A 0 · · ·
0
 x1 
diag(y2 ) 0 −A · · ·

0
 
 


0
 x2 

..
..
..
.. 
.
.
=
.
(3.3)
 

.
.
.
.
. 
c (mp+1)×1
 .. 




diag(yp ) 0

.
0 · · · −A
x
w∗
p (m+np)×1
(mp+1)×(m+np)
{z
}|
{z
}
|
z

Aw

For this scenario we only consider the case when A is a complex Gaussian random matrix.

Theorem 3.5. Consider the self-calibration model given in (2.3), where Aw is as in (3.3) and Aw,0
corresponds to the noiseless part of Aw . Let xmin := min1≤l≤p kxl k and xmax := max1≤l≤p kxl k.
Then, for the solution ẑ of (2.8) and α = w∗cz0 , there holds


2
kẑ − αz0 k
≤ κ(Aw,0 )η 1 +
kαz0 k
1 − κ(Aw,0 )η
if κ(Aw,0 )η < 1 where η =

κ(Aw,0
and for kwk =

√

2kδAk
√
.
m

v
u
u
)≤t

Here the upper bound of κ(Aw,0 ) obeys

max{pd2max , x2m }
6x2max (m + kwk2 )
min
.,
x2min min{m, kwk2 | Corr(w, z0 )|2 } min{pd2min , x2m }
max

m, there holds
v
u
√
u max{pd2max , x2m }
2 3xmax
min
t
κ(Aw,0 ) ≤
,
xmin | Corr(w, z0 )| min{pd2min , x2m }
max

with probability at least 1 − 2m(np + m)−γ if A is a complex Gaussian random matrix and




n
1
kGk kGk2F
,
(3.4)
+
log2 (np + m) ≤
C0 max
2
p
p
m
16(γ + 1)
11

where G is the Gram matrix of

n

C0

xl
kxl k



op

l=1

n
1
+
p m

. In particular, if G = Ip and kGkF =



log2 (np + m) ≤

√

p, (3.4) becomes

1
.
16(γ + 1)

Remark 3.6. When kδAk = 0, Theorem 3.5qsays that the solution to (2.3) is uniquely determined
up to a scalar if mp = O(max{np + mkGk, mnp2 + m2 kGk2F }), which involves the norm of the
Gram matrix G. This makes sense if we consider two extreme cases: if {vl }pl=1 are all identical,
then we have kGk = p and kGkF = p and if the {vl }pl=1 are mutually orthogonal, then G = Ip .
Remark 3.7. Balzano and Nowak [6] show exact recovery of this model when A is a deterministic
DFT matrix (discrete Fourier matrix) and {xl }pl=1 are generic signals drawn from a probability
distribution, but their results do not include stability theory in the presence of noise.
Remark 3.8. For Theorem 3.3 and Theorem 3.5 it does not come as a surprise that the error
bound depends on the norm of {xl }pl=1 and D as well as on how much z0 and w are correlated.
We cannot expect a relatively good condition number for Aw if kxl k varies greatly over 1 ≤ l ≤ p.
Concerning the correlation between z0 and w, one extreme case is hz0 , wi = 0, which does not rule
out the possibility of z = 0. Hence, the quantity hz0 , wi affects the condition number.

3.4

Theoretical results for the spectral method

Let S be the matrix Aw excluding the last row and consider S = S0 + δS where S0 is the noiseless
part of S and δS is the noise term. The performance under noise depends on the second smallest
singular value of S0 and the noise strength kδSk.
Theorem 3.9. Denote ẑ as the solution to (2.10), i.e., the right singular vector of S w.r.t. the
smallest singular value and kẑk = 1. Then there holds,
kδSk
(I − ẑ ẑ ∗ )z0
kα0 ẑ − z0 k
≤
=
kz0 k
kz0 k
[σ2 (S0 ) − kδSk]+
α0 ∈C
min

where σ2 (S0 ) is the second smallest singular value of S0 , z0 satisfies S0 z0 = 0, and ẑ is the right
singular vector with respect to the smallest singular value of S, i.e., the solution to (2.10).
Moreover, the lower bound of σ2 (S0 ) satisfies
q
√
(a) σ2 (S0 ) ≥ p2 min{ m, dmin kxk} for model (2.1) under the assumption of Theorem 3.1;
(b) σ2 (S0 ) ≥

√1 xmin min
2

(c) σ2 (S0 ) ≥

√1 xmin min
2

n√
√ o
m
pdmin , xmax
for model (2.2) under the assumption of Theorem 3.3;
n√
√ o
m
for model (2.3) under the assumption of Theorem 3.5.
pdmin , xmax

Remark 3.10. Note that finding z which minimizes (2.10) is equivalent to finding the eigenvector
with respect to the smallest eigenvalue of S ∗ S. If we have a good approximate upper bound λ of
kS ∗ Sk, then it suffices to find the leading eigenvector of λI − S ∗ S, which can be done efficiently by
using power iteration.
How to choose λ properly in practice? We do not want to choose λ too large since this would
imply slow convergence of the power iteration. For each case, it is easy to get a good upper bound
of kSk2 based on the measurements and sensing matrices,
12

(a) for (2.1), λ = k

Pp

∗
l=1 diag(yl ) diag(yl )k

+k

Pp

∗
l=1 Al Al k;

(b) for (2.2), λ = k diag(yl ) diag(yl∗ )k + max1≤l≤p kAl k2 ;
(c) for (2.3), λ = k diag(yl ) diag(yl∗ )k + kAk2 .
Those choices of λ are used in our numerical simulations.

4

Numerical simulations

This section is devoted to numerical simulations. Four experiments for both synthetic and real data
will be presented to address the effectiveness, efficiency and robustness of the proposed approach.
For all three models presented (2.1), (2.2) and (2.3), the corresponding linear systems have
simple block structures which allow for fast implementation via the conjugate gradient method for
non-hermitian matrices [34]. In our simulations, we do not need to set up Aw explicitly to carry
out the matrix-vector multiplications arising in the conjugate gradient method. Moreover, applying
preconditioning via rescaling all the columns to be of similar norms can give rise to an even faster
convergence rate. Therefore, we are able to deal with medium- or large-scale problems from image
processing in a computationally efficient manner.
In our simulations the iteration stops if either the number of iterations reaches at most 2000 or
the residual of the corresponding normal equation is smaller than 10−8 . Throughout our discussion,
the SNR (signal-to-noise ratio) in the scale of dB is defined as
 Pp
2
l=1 kyl k
P
.
SNR := 10 log10
p
2
l=1 kεl k
We measure the performance with RelError (in dB) := 20 log 10 RelError where
(
)
kα1 x̂ − x0 k
kα2 dˆ − d0 k
RelError = max min
, min
.
kx0 k
kd0 k
α1 ∈C
α2 ∈C

Here d0 and x0 are the ground truth. Although RelError does not match the error bound in our
theoretic analysis, certain equivalent relations hold if one assumes all |di,0 | are bounded away from
0 because there holds ŝ ≈ s0 if and only if dˆ ≈ d0 .
For the imaging examples, we only measure the relative error with respect to the recovered
x̂−x0 k
.
image x̂, i.e., minα1 ∈C kα1kx
0k

4.1

Self-calibration from repeated measurements

Suppose we have a target image x and try to estimate x through multiple measurements. However,
the sensing process is not perfect because of the missing calibration of the sensors. In order to
estimate both the unknown gains and phases as well as the target signal, a randomized sensing
procedure is used by employing several random binary masks.
We assume that yl = DHMl x0 + εl where H is a “tall” low-frequency DFT matrix, Ml is a
diagonal ±1-random matrix and x0 is an image of size 512×512. We set m = 10242 , n = 5122 , p = 8
2
pm
= 6.4. We
and D = diag(d) ∈ 10242 × 10242 with d ∈ C1024 ; the oversampling ratio is m+n
compare two cases: (i) d is a sequence distributed uniformly over [0.5, 1.5] with w = 1m+n , and (ii)
13


0m
. We
d is a Steinhaus sequence (uniformly distributed over the complex unit circle) with w =
1n
pick those choices of w because we know that the image we try to reconstruct has only non-negative
values. Thus, by choosing w to be non-negative, there are fewer cancellation in the expression w ∗ z0 ,
which in turn leads to a smaller condition number and better robustness. The corresponding results
of our simulations are shown in Figure 1 and Figure 2, respectively. In both cases, we only measure
the relative error of the recovered image.


Figure 1: Here m = 10242 , n = 5122 , p = 8, SNR=5dB, D = diag(d) where di is uniformly
distributed over [0.5, 1.5]. Left: original image; Middle: uncalibrated image, RelError = −13.85dB;
Right: calibrated image, RelError = −20.23dB

Figure 2: Here m = 10242 , n = 5122 , p = 8, SNR=5dB, D = diag(d) where d is a Steinhauss
sequence. Left: original image; Middle: uncalibrated image; Right: calibrated image, RelError =
−10.02dB
In Figure 1, we can see that both, the uncalibrated and the calibrated image are quite good.
Here the uncalibrated image is obtained by first applying the inverse Fourier transform and the
inverse of the mask to each yi and then taking the average of p samples. We explain briefly why
the uncalibrated image still looks good. Note that
x̂uncali

p

p

l=1

l=1

1 X −1 †
1 X −1 †
Ml H (DHMl x0 ) = x0 +
Ml H (D − I)HMl x0
=
p
p

1
H ∗ is the pseudo inverse of H. Here D−I is actually a diagonal matrix with random
where H † = m
1
entries ± 2 . As a result, each Ml−1 H † (D −I)HMl x0 is the sum of m rank-1 matrices with random
± 12 coefficients and is relatively small due to many cancellations. Moreover, [14] showed that most
2-D signals can be reconstructed within a scale factor from only knowing the phase of its Fourier
transform, which applies to the case when d is positive.
However, when the unknown calibration parameters are complex variables (i.e., we do not know
much about the phase information), Figure 2 shows that the uncalibrated recovered image is totally

14

meaningless. Our approach still gives a quite satisfactory result even at a relatively low SNR of
5dB.

4.2

Blind deconvolution in random mask imaging

The second experiment is about blind deconvolution in random mask imaging [5, 37]. Suppose we
observe the convolution of two components,
y l = h ∗ M l x 0 + εl ,

1≤l≤p

where both, the filter h and the signal of interests x0 are unknown. Each Ml is a random ±1-mask.
The blind deconvolution problem is to recover (h, x0 ). Moreover, here we assume that the filter is
actually a low-pass filter, i.e., F(h) is compactly supported in an interval around the origin, where
F is the Fourier transform. After taking the Fourier transform on both sides, the model actually
ends up being of the form (2.1) with Al = HMl where H is a “fat” partial DFT matrix and d is
the nonzero part of F(h). In our experiment, we let x0 be a 128 × 128 image and d = F(h) be a
2-D Gaussian filter of size 45 × 45 as shown in Figure 3.
In those experiments, we choose w = 1m+n since both d and x0 are nonnegative. Figure 4
shows the recovered image from p = 32 sets of noiseless measurements and the performance is quite
pm
≈ 3.52. We can see from Figure 5 that the blurring
satisfactory. Here the oversampling ratio is m+n
effect has been removed while the noise still exists. That is partially because we did not impose any
denoising procedure after the deconvolution. A natural way to improve this reconstruction further
would be to combine the blind deconvolution method with a total variation minimization denoising
step.

4.3

Blind deconvolution via diverse inputs

We choose Al to be random Hadamard matrices with m = 256 and n = 64 and D = diag(d0 )
with d0 being a positive/Steinhaus sequence, as we did previously. Each xl is sampled from
1m
standard Gaussian distribution. We choose w =
if d0 is uniformly distributed over [0.5, 1.5]
0np×1

√
me1
for Steinhaus d0 . 10 simulations are performed for each level of SNR. The
and w =
0np×1
pm
is 2, 2.67 and 3 for
test is also given under different choices of p. The oversampling ratio pn+m
p = 4, 8, 12 respectively. From Figure 6, we can see that the error scales linearly with SNR in
dB. The performance of Steinhaus d0 is not as good as that of positive d0 for SNR ≤ 10 when
∗
we use
 squares method. That is because w z0 is quite small when d0 is complex and
√the least

me1
. Note that the error between ẑ and z0 is bounded by κ(Aw,0 )η 1 + 1−κ(A2 w,0 )η .
w =
0np×1
Therefore, the error bound does not depend on kδAk linearly if κ(Aw,0 )η is close to 1. This may
explain the nonlinear behavior of the relative error in the low SNR regime.
We also apply spectral method to this model with complex gains d0 and Al chosen as either a
random Hadamard matrix or a Gaussian matrix. Compared to the linear least squares approach,
the spectral method is much more robust to noise, as suggested in Figure 7, especially in the low
SNR regime.

15

Figure 3: Left: Original image; Right: Gaussian filter in Fourier domain. The support of the filter
is 45 × 45 and hence m = 452 = 2025, n = 1282 , p = 32.

Figure 4: Left: Blurred image without noise, Right: Recovered image, RelError = -45.47dB

Figure 5: Left: Blurred image with SNR = 5dB, Right: Recovered image, RelError = -5.84dB

4.4

Multiple snapshots

We make a comparison of performances between the linear least squares approach and the spectral
method when d0 is a Steinhaus sequence and A is a Gaussian random matrix A. Each xl is sampled
from the standard Gaussian distribution and hence the underlying Gram matrix G is quite close
to Ip (this closeness could be easily made more precise, but we refrain doing so here). The choice
of w and oversampling ratio are the same as those in Section 4.3. From Figure 8, we see that
the performance is not satisfactory for the Steinhaus case using the linear least squares approach,
especially in the lower SNR regime
√ (SNR
 ≤ 10). The reason is the low correlation between w and
me1
z0 if d0 is Steinhaus and w =
. The difficulty of choosing w is avoided by the spectral
0np×1
method. As we can see in Figure 8, the relative error given by spectral method is approximately
7dB smaller than that given by linear least squares approach when d0 is a complex vector and the
SNR is smaller than 20dB.
16

m = 256, n = 64, d:Positive, LS method, Hadamard

m = 256, n = 64, d:Steinhaus, LS method, Hadamard

0

0
p=4
p=8
p = 12

−20

−30

−40

−50

−60

−70

−80

p=4
p=8
p = 12

−10

Average of sin(θ) of 10 Samples(dB)

Average of sin(θ) of 10 Samples(dB)

−10

−20

−30

−40

−50

−60

−70

0

10

20

30

40

50

60

−80

70

0

10

20

30

SNR(dB)

40

50

60

70

SNR(dB)

Figure 6: Performance of linear least squares approach: RelError (in dB) vs. SNR for yl =
DAl xl + εl , 1 ≤ l ≤ p where m = 256, n = 64, D = diag(d0 ) and each Al is a random Hadamard
matrix. d0 is a Steinhaus sequence.
m = 256, n = 64, d:Steinhaus, spectral method, Gaussian

m = 256, n = 64, d:Steinhaus, spectral method, Hadamard

0

0
p=4
p=8
p = 12

−20

−30

−40

−50

−60

−70

−80

p=4
p=8
p = 12

−10

Average of sin(θ) of 10 Samples(dB)

Average of sin(θ) of 10 Samples(dB)

−10

−20

−30

−40

−50

−60

−70

0

10

20

30

40

50

60

−80

70

0

10

20

30

SNR(dB)

40

50

60

70

SNR(dB)

Figure 7: Performance of spectral method: RelError (in dB) vs. SNR for yl = DAl xl +εl , 1 ≤ l ≤ p
where m = 256, n = 64, D = diag(d0 ). Here Al is either a random Hadamard matrix or a Gaussian
matrix. d0 is a Steinhaus sequence.

5

Proofs

For each subsection, we will first give the result of noiseless measurements. We then prove the
stability theory by using the result below. The proof of spectral method can be found in Section 5.4.
Proposition 5.1. [46] Suppose that Au0 = b is a consistent and overdetermined system. Denote
û as the least squares solution to k(A + δA)u − bk2 with kδAk ≤ ηkAk. If κ(A)η < 1, there holds,


kû − u0 k
2
≤ κ(A)η 1 +
.
ku0 k
1 − κ(A)η
To apply the proposition above, it suffices to bound κ(A) and η.

17

m = 256, n = 64, d:Steinhaus, LS method

m = 256, n = 64, d:Steinhaus, spectral method

0

0
p=4
p=8
p = 12

−20

−30

−40

−50

−60

−70

−80

p=4
p=8
p = 12

−10

Average of sin(θ) of 10 Samples(dB)

Average of sin(θ) of 10 Samples(dB)

−10

−20

−30

−40

−50

−60

−70

0

10

20

30

40

50

60

−80

70

0

10

SNR(dB)

20

30

40

50

60

70

SNR(dB)

Figure 8: Comparison between linear least squares approach and spectral method: RelError (in
dB) vs. SNR for yl = DAxl + εl , 1 ≤ l ≤ p where m = 256, n = 64, D = diag(d0 ) and A is a
Gaussian random matrix. The gain d0 is a random vector with each entry uniformly distributed
over unit circle.

5.1

Self-calibration from repeated measurements

Let us start with (2.7) when εl = 0 and denote Λl := diag(Al v) and v :=

x
kxk

∈ Cn ,

  ∗


Λ1 − √1m A1 
diag(y1 ) −A1

 Dkxk √ 0

..
..  =  ..
..
A0 := 
.

.
.   .
.
0
mIn
1
∗
Λp − √m Ap
diag(yp ) −Ap

Then we rewrite A∗w Aw as

A∗w Aw = A∗0 A0 + ww ∗ =
where
Zl :=

"

Λl

− √1m A∗l

By definition,

#

∈C

(m+n)×m

,

p
X

P Zl Zl∗ P ∗ + ww ∗

l=1




D∗ kxk
0
√
P :=
∈ C(m+n)×(m+n) .
0
mIn

(5.1)

#
∗
√1 Λl Al
−
Λ
Λ
l
l
m
∈ C(m+n)×(m+n) .
Zl Zl∗ =
1
∗A
− √1m A∗l Λ∗l
A
l
m l
"

Our goal is to find out P
the smallest and the largest eigenvalue of Aw . Actually it suffices to under∗
stand the spectrum of " pl=1 Z
# l Zl . Obviously, its smallest eigenvalue is zero and the corresponding
√1
2

1m
√
m

. Let al,i be the i-th column of A∗l and we have E(al,i a∗l,i ) = In under
"
#
√1 1m v ∗
I
−
m
m
all the three settings in Section 3.1. Hence, C := E(Zl Zl∗ ) =
.
− √1m v1∗m
In
It is easy to see that rank(C) = m + n − 1 and the null space of C is spanned by u1 . C has an
eigenvalue with value 1 of multiplicity m + n − 2 and an eigenvalue with value 2 of multiplicity 1.
eigenvector is u1 :=

v

18

More importantly, the following proposition holds and combined with Proposition 5.1, we are able
to prove Theorem 3.1.
Proposition 5.2. There holds

p
X
l=1

(a) with probability 1 − (m +

n)−γ

Zl Zl∗ − pC ≤

p
2

 n
if Al is Gaussian and p ≥ c0 γ max 1, m
log2 (m + n);

(b) with probability 1 − (m + n)−γ − 2(mp)−γ+1 if each Al is a “tall” (m × n, m ≥ n) random
Hadamard/DFT matrix and p ≥ c0 γ 2 log(m + n) log(mp);
(c) with probability 1 − (m + n)−γ − 2(mp)−γ+1 if each Al is a “fat” (m × n, m ≤ n) random
Hadamard/DFT matrix and mp ≥ c0 γ 2 n log(m + n) log(mp).
Remark 5.3. Proposition 5.2 actually addresses the identifiability issue of the model (2.1) in
absence of noise. More precisely, the invertibility of P is guaranteed by that
P of D. By Weyl’s
theorem for singular value perturbation in [36], mP
+ n − 1 eigenvalues of pl=1 Zl Zl∗ are greater
than p2 . Hence, the rank of A0 is equal to rank( pl=1 Zl Zl∗ ) = m + n − 1 if p is close to the
n
). In other words, the
information theoretic limit under the conditions given above, i.e., p ≥ O( m
null space of A0 is completely spanned by z0 := (s0 , x0 ).
5.1.1

Proof of Theorem 3.1

Note that Proposition 5.2 gives the result if εl = 0. The noisy counterpart is obtained by applying
perturbation theory for linear least squares.
Let Aw := Aw,0 + δA where Aw,0 is the noiseless part and δA is defined in (2.9).
  Note
0
by the
that αz0 with α = w∗cz0 is actually a solution to the overdetermined system Aw,0 z =
c
definition of Aw,0 . Proposition 5.1 implies that it suffices to estimate the condition number κ(Aw,0 )
of Aw,0 and η such that kδAk ≤ ηkAw,0 k holds. Note that
!
p
X
A∗w,0 Aw,0 = P
(5.2)
Zl Zl∗ P ∗ + ww ∗
Proof:

= P

l=1
p
X

Zl Zl∗

l=1

e :=
where w

P −1 w.

ew
e
+w

∗

!

e ∗
P ∗ =: P CP

(5.3)

From Proposition 5.2 and Theorem 1 in [36], we know that
!
!
p
p
X
X
p
5p
Zl Zl∗ ≥ , λn+m
λ2
Zl Zl∗ ≤
2
2
l=1
l=1
Pp
where λ1 ≤ · · · ≤ λn+m and λ1 ( l=1 Zl Zl∗ ) = 0. Following from (5.2), we have
kA∗w,0 Aw,0 k

2

≤ kP k


p
X
l=1

Zl Zl∗ + kwk2 ≤ 3p max{d2max kxk2 , m} + kwk2

kwk2
≤ 3 p+
m



λ2max (P ).
19

(5.4)

On the other hand,
kA∗w,0 Aw,0 k ≥

p
X
l=1

A∗l Al ≥

mp
2

(5.5)

follows from Proposition 5.2. In other words, we have found the lower and upper bounds for
kAw,0 k or equivalently, λmax (A∗w,0 Aw,0 ). Now we proceed to the estimation of λmin (A∗w,0 Aw,0 ).
#
"
1m
√
P
Pm+n
m+n
1
2
m
with m+n
Let u := j=1 αj uj be a unit vector, where u1 = √2
j=1 |αj | = 1 and {uj }j=1
v
P
e defined in (5.3) has a lower
are the eigenvectors of pl=1 Zl Zl∗ . Then the smallest eigenvalue of C
bound as follows:
!
p
X
e
ew
e ∗u
Zl Z ∗ u + u∗ w
u∗ Cu
= u∗
l

l=1

≥

≥
e ≥
which implies λmin (C)
e ∗,
P CP

1
2

m+n
X
j=2

m+n
X
j=2

ew
e ∗ u1 u∗1 u
λj |αj |2 + u∗ u1 u∗1 w

λj |αj |2 +

o
n
∗ z |2
0
≥
min p, |w
mkxk2

1
2

|w ∗ z0 |2
|α1 |2
2mkxk2

(5.6)

n
o
2
2
0 )|
min p, kwk | Corr(w,z
. Combined with A∗w,0 Aw,0 =
m


e 2 (P ) ≥ 1 min mp, kwk2 | Corr(w, z0 )|2 λ2 (P ).
λmin (A∗w,0 Aw,0 ) ≥ λmin (C)λ
min
min
2m

Therefore, with (5.4), the condition number of A∗w,0 Aw,0 is bounded by
κ(A∗w,0 Aw,0 ) ≤
From (5.5), we set η =

2kδAk
√
mp

≥

6(mp + kwk2 )
κ2 (P ).
min{mp, kwk2 | Corr(w, z0 )|2 }

kδAk
kAw,0 k .

Applying Proposition 5.1 gives the following upper bound of the estimation error


κ(Aw,0 )η 1 + 1−κ(A2 w,0 )η where α = w∗cz0 and
κ(Aw,0 ) ≤
If kwk =

√

s

6(mp + kwk2 )
κ(P ).
min{mp, kwk2 | Corr(w, z0 )|2 }

mp, the upper bound of κ(Aw,0 ) satisfies
κ(Aw,0 ) ≤

s

s
√
2 3
12
max{d2max kxk2 , m}
κ(P
)
≤
.
| Corr(w, z0 )|2
|Corr(w, z0 )| min{d2min kxk2 , m}

20

kẑ−αz0 k
kαz0 k

≤

5.1.2

Proof of Proposition 5.2(a)

Proof: [Proof of Proposition 5.2(a)] From now on, we assume al,i ∈ Cn , i.e., the i-th column
of A∗l , obeys a complex Gaussian distribution, √12 N (0, In )+ √i2 N (0, In ). Let zl,i be the i-th column
of Zl ; it can be written in explicit form as
"
#
"
#
| hal,i , vi |2 ei e∗i − √1m ei v ∗ al,i a∗l,i
hal,i , viei
∗
, zl,i zl,i =
zl,i =
1
∗
− √1m al,i
− √1m al,i a∗l,i ve∗i
m al,i al,i .
∗ − E(z z ∗ ), we obtain
Denoting Zl,i := zl,i zl,i
l,i l,i
"
#
(| hal,i , vi |2 − 1)ei e∗i
− √1m ei v ∗ (al,i a∗l,i − In )
Zl,i :=
.
1
∗
− √1m (al,i a∗l,i − In )ve∗i
m (al,i al,i − In )

P
Obviously each Zl,i is independent. In order to apply Theorem 5.12 to estimate k l,i Zl,i k, we
Pp Pm
∗
need to bound maxl,i kZl,i kψ1 and
i=1 E(Zl,i Zl,i ) . Due to the semi-definite positivity of
l=1

∗ , we have kZ k ≤ max{kz k2 , k E(z z ∗ )k} ≤ max | ha , vi |2 + 1 ka k2 , 2 and hence
zl,i zl,i
l,i
l,i
l,i
l,i l,i
l,i
m

n
1
kZl,i kψ1 ≤ (| hal,i , vi |2 )ψ1 + (kal,i k2 )ψ1 ≤ C 1 +
m
m

n
which follows from Lemm 5.14 and k·kψ1 is a norm. This implies R := maxl,i kZl,i kψ1 ≤ C 1 + m
.
P
P
p
m
∗
∗
∗
Now we consider σ02 =
i=1 E(Zl,i Zl,i ) by computing (Zl,i Zl,i )1,1 and (Zl,i Zl,i )2,2 , i.e.,
l=1
∗ ,
the (1, 1)-th and (2, 2)-th block of Zl,i Zl,i


1 ∗
2
2
2
∗
∗
(Zl,i Zl,i )1,1 = (| hal,i , vi | − 1) + v (al,i al,i − In ) v ei e∗i ,
m
1
1
∗
(al,i a∗l,i − In )vv ∗ (al,i a∗l,i − In ) + 2 (al,i a∗l,i − In )2 .
(Zl,i Zl,i
)2,2 =
m
m
Following from (5.36), (5.37), (5.38) and Lemma 5.10, there holds

p X
p X
m 
m
∗ )
X
X
0
E(Zl,i Zl,i
1,1
∗
2
E(Zl,i Zl,i ) ≤ 2
σ0 =
∗ )
0
E(Zl,i Zl,i
2,2
l=1 i=1
l=1 i=1



p X
m
n
X
1+ m
ei e∗i
0 
= 2
n
1
0
m + m2 In
l=1 i=1




n
n
Im
0 
1+ m
.
= 2p 1 +
= 2p
n
0
1 + m In
m
By applying the matrix Bernstein inequality (see Theorem 5.13) we obtain
p X
m
nr 
X
n p
≤ C0 max
Zl,i
p 1+
t + log(m + n),
m
l=1 i=1
o p

n
(t + log(m + n)) log(m + n) ≤
1+
m
2

 n
log2 (m+
with probability 1−e−t . In particular, by choosing t = γ log(m+n), i.e, p ≥ c0 γ max 1, m
n), the inequality above holds with probability 1 − (m + n)−γ .
21

5.1.3

Proof of Proposition 5.2(b)

Proof: [Proof of Proposition 5.2(b)] Each Zl is independent by its definition in (5.1) if
Al := HMl where H is an m × n partial DFT/Hadamard matrix with m ≥ n and H ∗ H = mIn
and Ml = diag(ml ) is a diagonal random binary ±1 matrix. Let Zl := Zl Zl∗ − C ∈ C(m+n)×(m+n) ;
in explicit form
#
"
− √1m (Λl Al − 1m v ∗ )
Λl Λ∗l − Im
.
Zl =
− √1m (A∗l Λ∗l − v1∗m )
0
where A∗l Al = mIn follows from the assumption. First we take a look at the upper bound of kZl k.
It suffices to bound kZl Zl∗ k since kZl k = kZl Zl∗ − Ck ≤ max {kZl Zl∗ k, kCk} and C is positive
semi-definite. On the other hand, due to Lemma 5.10, we have kZl Zl∗ k ≤ 2 max{kΛl Λ∗l k, 1} and
hence we only need to bound kΛl k. For kΛl k, there holds
max kΛl k = max kAl vk∞ =

1≤l≤p

1≤l≤p

max

1≤l≤p,1≤i≤m

| hal,i , vi |.

Also for any pair of (l, i), hal,i , vi can be rewritten as
| hal,i , vi | = | hdiag(ml )hi , vi | = | ml , diag(h̄i )v |
where hi is the i-th column of H ∗ and k diag(h̄i )vk = kvk = 1. Then there holds
P



max kΛl k ≥

1≤l≤p

p


p X
m


X
p
2γ log(mp)
≤
P | hal,i , vi | ≥ 2γ log(mp)
l=1 i=1



p
≤ mpP | hal,i , vi | ≥ 2γ log(mp)

≤ 2mp · e−γ log(mp) ≤ 2(mp)−γ+1 ,

(5.7)

where the third inequality follows from Lemma 5.11. Applying Lemma 5.10 to Zl ,
R := max kZl k ≤ 2 max {kΛl Λ∗l k, 1} ≤ 4γ log(mp),
1≤l≤p

1≤l≤p

with probability at least 1 − 2(mp)−γ+1
P. Denote the event {max1≤l≤p kZl k ≤ 4γ log(mp)} by E1 .
Now we try to understand σ02 = k pl=1 E(Zl Zl∗ )k. The (1, 1)-th and (2, 2)-th block of Zl Zl∗ are
given by
(Zl Zl∗ )1,1 = (Λl Λ∗l − Im )2 +
(Zl Zl∗ )2,2 =

1
(Λl Al − 1m v ∗ )(Λl Al − 1m v ∗ )∗ ,
m

1
(A∗ Λ∗ − v1∗m )(Λl Al − 1m v ∗ ).
m l l

By using (5.41), (5.42) and Al A∗l = HH ∗  mIm , we have
∗

E((Λl Λ∗l − Im )2 ) =
∗ ∗

E(Λl Al − 1m v )(Λl Al − 1m v )
E(A∗l Λ∗l

−

v1∗m )(ΛAl

=

∗

− 1m v ) =
=

22

E(Λl Λ∗l )2 − Im  2Im ,

E(Λl Al A∗l Λ∗l ) − 1m 1∗m  mIm ,
E(A∗l Λ∗l Λl Al ) − mvv ∗
m
X
| hal,i , vi |2 al,i a∗l,i − mvv ∗ 
i=1

(5.8)
(5.9)

3mIn .

(5.10)

Combining (5.8), (5.9), (5.10) and Lemma 5.10,



p 
X
0
E(Zl Zl∗ )1,1
2Im + Im 0
2
≤ 2p
≤ 6p.
σ0 ≤ 2
0
E(Zl Zl∗ )2,2
0
3In
l=1

By applying (5.32) with t = γ log(m + n) and R ≤ 4γ log(mp) over event E1 , we have
p
X
p
√ p
(Zl Zl∗ − C) ≤ C0 max{ p (γ + 1) log(m + n), γ(γ + 1) log(mp) log(m + n)} ≤
2
l=1

with probability 1 − (m + n)−γ − 2(mp)−γ+1 if p ≥ c0 γ 2 log(m + n) log(mp).
5.1.4

Proof of Proposition 5.2(c)

Proof: [Proof of Proposition 5.2(c)] Each Zl is independent due to (5.1). Let Zl := Zl Zl∗ −
C ∈ C(m+n)×(m+n) ; in explicit form
#
"
Λl Λ∗l − Im
− √1m (Λl Al − 1m v ∗ )
.
Zl =
1
∗
− √1m (A∗l Λl − v1∗m )
m Al Al − In
Here Al = HMl where H is a “fat” m × n, (m ≤ n) partial DFT/Hadamard matrix satisfying
HH ∗ = nIm and Ml is a diagonal ±1-random matrix. There holds A∗l Al = Ml∗ H ∗ HMl ∈ Cn×n
where H ∈ Cm×n and E(A∗l Al ) = mIn . For each l, kA∗l Al k = kH ∗ Hk = kHH ∗ k = n. Hence,
there holds,


nn
o
1
kA∗l Al k, kΛl k2 , 1 ≤ 2 max
, γ log(mp)
kZl k ≤ max {kZl Zl∗ k, kCkk ≤ 2 max
m
m

with probability at least 1 − 2(mp)−γ+1 , which
Pfollows exactly from (5.7) and Lemma 5.10.
Now we give an upper bound for σ02 := k pl=1 E(Zl Zl∗ )k. The (1, 1)-th and (2, 2)-th block of
Zl Zl∗ are given by
1
(Λl Al − 1m v ∗ )(ΛAl − 1m v ∗ )∗ ,
m

2
1 ∗
1
(A∗l Λ∗l − v1∗m )(Λl Al − 1m v ∗ ) +
Al Al − In .
m
m

(Zl Zl∗ )1,1 = (Λl Λ∗l − Im )2 +
(Zl Zl∗ )2,2 =

By using (5.41), (5.42) and Al A∗l = HH ∗ = nIm , we have
∗

E((Λl Λ∗l − Im )2 ) =
∗ ∗

E(Λl Al − 1m v )(Λl Al − 1m v )
E(A∗l Λ∗l

−

v1∗m )(Λl Al

=

∗

− 1m v ) =
=

For E

E(Λl Λ∗l )2 − Im  2Im ,

E(Λl Al A∗l Λ∗l ) − 1m 1∗m  nIm ,
E(A∗l Λ∗l Λl Al ) − mvv ∗
m
X
| hal,i , vi |2 al,i a∗l,i − mvv ∗ 
i=1

(5.11)
(5.12)

3mIn . (5.13)

2
− In , we have

2
1 ∗
n − 2m ∗
2
1
Al Al − In = 2 A∗l Al A∗l Al − A∗l Al + In =
Al Al + In
m
m
m
m2

1
∗
m Al Al

23

where Al A∗l = nIm . Note that E(A∗l Al ) = mIn , and there holds,

2
1 ∗
n−m
A Al − In =
In .
E
m l
m

(5.14)

Combining (5.11), (5.12), (5.13), (5.14) and Lemma 5.10,



n
6np
Im
0 
2+ m
2
σ0 ≤ 2p
≤
.
n
0
2 + m In
m

By applying (5.32) with t = γ log(m + n), we have
r
p

X
p
np p
n
∗
log(m + n)} ≤
(γ + 1) log(m + n), (γ + 1) γ log(mp) +
(Zl Zl − C) ≤ C0 max{
m
m
2
l=1

with probability 1 − (m + n)−γ − 2(mp)−γ+1 if mp ≥ c0 γ 2 n log(m + n) log(mp).

5.2

Blind deconvolution via diverse inputs

We start with (3.1) by setting εl = 0. In this way, we can factorize the matrix Aw
last row) into
√

pD √ 0
  diag(A

A1
1 v1 )
√
√
−
·
·
·
0
mIn
kx1 kIm · · ·
0

p
m
 0

kx1 k


.
.
.
.
.
.
.
..
..
..
..  
..
..
..
A0 := 

..
  ...

.

Ap
diag(Ap vp )
0
· · · kxp kIm
√
0
· · · − √m
p
|
{z
}|
0
0
{z
}
Q
{z
|
mp×(np+m)
Z∈C

P

(excluding the
···
···
..
.
···

0
0
..
.
√

mIn
kxp k








}

(5.15)
where vl =
is the normalized xl , 1 ≤ l ≤ p. We will show that the matrix Z is of rank
 
v1
 .. 
np + m − 1 to guarantee that the solution is unique (up to a scalar). Denote v :=  .  ∈ Cnp×1
vp
P
√
and pl=1 eel ⊗ vl = v with kvk = p where {e
el }pl=1 is a standard orthonormal basis in Rp .
xl
kxl k

5.2.1

Proof of Theorem 3.3

The proof of Theorem 3.3 relies on the following proposition. We defer the proof of Proposition 5.4
to the Sections 5.2.2 and 5.2.3.
Proposition 5.4. There holds,
1
kZ ∗ Z − Ck ≤ ,
2

"

1
Im
− √mp
1m v ∗
C := E(Z ∗ Z) =
1
Inp
v1∗m
− √mp

#

(5.16)

(a) with probability at least 1 − (np + m)−γ with γ ≥ 1 if Al is an m × n (m > n) complex Gaussian
random matrix and


n
1
1
+
(γ + 1) log 2 (np + m) ≤ ;
C0
p m
4
24

(b) with probability at least 1 − (np + m)−γ − 2(mp)−γ+1 with γ ≥ 1 if Al yields (3.2) and


n−1
1
1
+
γ 3 log4 (np + m) ≤ .
C0
p m−1
4
Remark 5.5. Note that C has one eigenvalue equal to 0 and all the other eigenvalues are at
least 1. Hence the rank of C is np + m − 1. Similar to Remark 5.3, Proposition 5.4 shows
that the solution (s0 , {xl }pl=1 ) to (3.1) is uniquely identifiable with high probability when mp ≥
(np + m) · poly(log(np + m)) and kZ ∗ Z − Ck ≤ 12 .
Proof: [Proof of Theorem 3.3] From (3.1), we let Aw = Aw,0 + δA where Aw,0 is the noiseless
part of Aw . By definition of Aw,0 , we know that αz0 is the solution to the overdetermined system
0
Aw,0 z =
where α = w∗cz0 . Now, (5.15) gives
c
A∗w,0 Aw,0 = A∗0 A0 + ww ∗ = P ∗ Z ∗ Q∗ QZP + ww ∗ .

Define xmax := max1≤l≤p kxl k and xmin := min1≤l≤p kxl k. From Proposition 5.4 and Theorem 1
in [36], we know that the eigenvalues {λj }1≤j≤np+m of Z ∗ Z fulfill λ1 = 0 and λj ≥ 12 for j ≥ 2
since kZ ∗ Z − Ck ≤ 12 ; and the eigenvalues of C are 0, 1 and 2 with multiplicities 1, np + m − 2, 1
respectively.
The key is to obtain a bound for κ(Aw,0 ). From (5.15),
λmax (A∗w,0 Aw,0 ) ≤ kP k2 kQk2 kZ ∗ Zk + kwk2 ≤ 3x2max λmax (P ∗ P ) + kwk2


kwk2
2
λmax (P ∗ P )
≤ 3xmax 1 +
m

(5.17)

where x2max λmax (P ∗ P ) ≥ x2max x2m ≥ m. On the other hand, (5.16) gives
min

λmax (A∗w,0 Aw,0 ) ≥ max {kA∗l Al k} ≥
1≤l≤p

since

1
∗
m Al Al

m
2

− In ≤ 12 . For λmin (A∗w,0 Aw,0 ),

e .
A∗w,0 Aw,0  λmin (Q∗ Q)P ∗ Z ∗ ZP + ww ∗  x2min P ∗ Z ∗ ZP + ww ∗ =: P ∗ CP
#
" 1
√ 1m
m
e = x2 Z ∗ Z + w
ew
e ∗ where w
e = P −1 w. By
Denote u1 := √12 √
such that Zu1 = 0 and C
1
min
v
p
using the same procedure as (5.6),

where u :=
Since

Pnp+m

j=1
|u∗1 (P −1 )∗ w|2

e ≥ x2
u∗ Cu
min

αj uj with
=

P

1
∗
2mp |w z0

np+m
X
j=2

λj |αj |2 + |α1 |2 |u∗1 P −1 w|2

1
2
j |αj | = 1 and λj ≥ 2 for
|2 , the smallest eigenvalue of

j ≥ 2 follows from Proposition 5.4.
e satisfies
C



|w ∗ z0 |2
x2min
min 1,
2
mpx2min


1
x2min
2
2
min 1, kwk | Corr(w, z0 )|
≥
2
m

e ≥
λmin (C)

25

∗

2

∗

2

z0 |
z0 |
≥ |w
≥ kwk2 | Corr(w, z0 )|2 follows from kz0 k2 ≥ px2min .
where |w
kz0 k2
px2min
Therefore, the smallest eigenvalue of A∗w,0 Aw,0 satisfies

e min (P ∗ P )
λmin (A∗w,0 Aw,0 ) ≥ λmin (C)λ

x2
≥ min min m, kwk2 | Corr(w, z0 )|2 λmin (P ∗ P )
2m

(5.18)

Combining (5.17) and (5.18) leads to
κ(A∗w,0 Aw,0 ) ≤
Applying Proposition 5.1 and η =

6x2max (m + kwk2 )
κ(P ∗ P ).
x2min min{m, kwk2 | Corr(w, z0 )|2 }

2kδAk
√
m

≥

kδAk
kAw,0 k ,

we have



2
kẑ − αz0 k
≤ κ(Aw,0 )η 1 +
,
kαz0 k
1 − κ(Aw,0 )η

α=

c
w ∗ z0

if κ(Aw,0 )η < 1 where κ(Aw,0 ) obeys
v
u
max{pd2max , x2m }
u
6x2max (m + kwk2 )
min
t
.
κ(Aw,0 ) ≤
x2min min{m, kwk2 | Corr(w, z0 )|2 } min{pd2min , x2m }
max

In particular, if kwk =

√

m, then κ(Aw,0 ) has the following simpler upper bound:
v
u
√
u max{pd2max , x2m }
2 3xmax
min
t
κ(Aw,0 ) ≤
|Corr(w, z0 )|xmin min{pd2min , x2m }
max

which finishes the proof of Theorem 3.3.
5.2.2

Proof of Proposition 5.4(a)

In this section, we will prove that Proposition 5.4(a) if al,i ∼ √12 N (0, In ) + √i2 N (0, In ) where
al,i ∈ Cn is the i-th column of A∗l . Before moving to the proof, we compute a few quantities which
will be used later. Define zl,i as the ((l − 1)m + i)-th column of Z ∗ ,
"
#
√1 hal,i , vl iei
p
zl,i :=
, 1 ≤ l ≤ p, 1 ≤ i ≤ m
− √1m eel ⊗ al,i
(np+m)×1

m
where {ei }m
el }pl=1 ∈ Rp are standard orthonormal P
basis P
in Rm and Rp respectively;
i=1 ∈ R and {e
p
∗
“⊗” denotes Kronecker product. By definition, we have Z ∗ Z = l=1 m
i=1 zl,i zl,i and all zl,i are
independent from one another.
#
"
1
1
2
∗
− √mp
hvl , al,i i ei (e
e∗l ⊗ a∗l,i )
p | hal,i , vl i | ei ei
∗
zl,i zl,i =
1
1
el ee∗l ⊗ al,i a∗l,i
hal,i , vl i (e
el ⊗ al,i )e∗i
− √mp
me

26

and its expectation is equal to
"

∗
)=
E(zl,i zl,i

Pp

It is easy to verify that C =

l=1

1
∗
p ei ei

1
− √mp
(e
el ⊗ vl )e∗i

Pm

∗
i=1 E(zl,i zl,i ).

#
1
− √mp
ei (e
e∗l ⊗ vl∗ )
.
1
el ee∗l ⊗ In
me

Proof: [Proof of Proposition 5.4(a)] The
is to use apply the matrix Bernstein inPpkeyPhere
m
∗ . Let Z := z z ∗ − E(z z ∗ ) and
∗
equality in Theorem 5.13. Note that Z Z = l=1 i=1 zl,i zl,i
l,i l,i
l,i
l,i l,i
we have


1
1 1
1
2
∗
2
2
kZl,i k ≤ kzl,i k + k E(zl,i zl,i )k ≤ | hal,i , vl i | + kal,i k + 2 max
,
p
m
p m
∗ )k ≤ 2 max{ 1 , 1 } follows from Lemma 5.10. Therefore, the exponential norm of
since k E(zl,i zl,i
p m
kZl,i k is bounded by






1
1
n
1
1 1
2
2
kZl,i kψ1 ≤ 2
(| hal,i , vl i | )ψ1 + (kal,i k )ψ1 + 2 max
,
+
≤C
p
m
p m
p m


n
.
which follows from Lemma 5.14 and as a result R := maxl,i kZl,i kψ1 ≤ C p1 + m
P
P
p
m
∗
Now we proceed by estimating the variance σ02 :=
i=1 Zl,i Zl,i . We express Zl,i as
l=1
follows:
#
"
1
1
2
∗
− √mp
ei (e
e∗l ⊗ (vl∗ (al,i a∗l,i − In )))
p (| hal,i , vl i | − 1)ei ei
.
Zl,i =
1
1
el ee∗l ⊗ (al,i a∗l,i − In )
(e
el ⊗ ((al,i a∗l,i − In )vl ))e∗i
− √mp
me
∗ are
The (1, 1)-th and the (2, 2)-th block of Zl,i Zl,i
∗
)1,1 =
(Zl,i Zl,i

and

2
1
1 ∗
| hal,i , vl i |2 − 1 ei e∗i +
v (al,i a∗l,i − In )2 vl ei e∗i ,
2
p
mp l

∗
)2,2 = eel ee∗l ⊗
(Zl,i Zl,i




1
1
(al,i a∗l,i − In )vl vl∗ (al,i a∗l,i − In ) + 2 (al,i a∗l,i − In )2 .
mp
m

Following from (5.36), (5.37) and (5.38), we have




1
1
n
n
∗
∗
∗
∗
E(Zl,i Zl,i )1,1 =
+
In + 2 In .
ei ei , E(Zl,i Zl,i )2,2 = eel eel ⊗
p2 mp
mp
m

Due to Lemma 5.10, there holds,
σ02 :=

p X
m
X
l=1 i=1

=



2 

1
p

∗
Zl,i Zl,i
≤2

+

n
m

0



Im



p X
m 
X
E(Zl,i Z ∗ )1,1
l,i

l=1 i=1

0
1
p

+

n
m



27

Inp

0



 ≤2



0
∗ )
E(Zl,i Zl,i
2,2

n
1
+
p m



.



P P
∗
Note that σ02 ≥ k pl=1 m
i=1 E(Zl,i Zl,i )1,1 k =
matrix. By applying (5.34),
∗

kZ Z − Ck =

p X
m
X
l=1 i=1

1
p

+

Zl,i =

nr 1

n
m

∗ ) is a positive semi-definite
since E(Zl,i Zl,i

p X
m
X
l=1 i=1

∗
)−C
E(zl,i zl,i

np
≤ C0 max
+
t + log(np + m),
p m


o 1
1
n
+
(t + log(np + m)) log(np + m) ≤ .
p m
2

∗
With
 t = γlog(mp + n), there holds kZ Z − Ck ≤
n
(γ + 1) log2 (np + m) ≤ 14 .
C0 1p + m

5.2.3

1
2

with probability at least 1 − (np + m)−γ if

Proof of Proposition 5.4(b)

We prove Proposition 5.4 based on assumption (3.2). Denote al,i and hl,i as the i-th column of A∗l
and Hl∗ and obviously we have al,i = Ml hl,i . Denote Λl = diag(Al vl ) and let Zl be the l-th block
of Z ∗ in (5.15), i.e.,
#
"
√1 Λl
p
∈ C(np+m)×m .
Zl =
− √1m eel ⊗ A∗l
With A∗l Al = mIn , we have
"
Zl Zl∗ =

1
∗
p Λl Λl

1
eel ⊗ (Λl Al )∗
− √mp

#
1
− √mp
ee∗l ⊗ (Λl Al )
∈ C(np+m)×(np+m)
(e
el ee∗l ) ⊗ In

where the expectation of i-th row of Λl Al yields E(Λl Al )i = E(vl∗ al,i a∗l,i ) = vl∗ . Hence E(Λl Al ) =
1m vl∗ ∈ Cm×n . Its expectation equals
"
#
1
∗ ⊗ (1 v ∗ )
√1 e
e
I
−
m
m
l
p
mp l
E(Zl Zl∗ ) =
.
1
− √mp
eel ⊗ (vl 1∗m )
eel ee∗l ⊗ In

Proof: [Proof of Proposition 5.4(b)] Note that each block Zl is independent and we want
to apply the matrix Bernstein inequality to achieve the desired result. Let Zl := Zl Zl∗ − E(Zl Zl∗ )
and by definition, we have
#
"
1
1
∗
− √mp
(Λl Al − 1m vl∗ )
p (Λl Λl − Im )
kZl k =
1
− √mp
(Λl Al − 1m vl∗ )∗
0
≤

1
1
kΛl Λ∗l − Im k + √ kΛl Al − 1m vl∗ k.
p
mp

Note that

√
√
kΛl Al − 1m vl∗ k ≤ kΛl kkAl k + k1m vl∗ k ≤ mkΛl k + m.


p
Since (5.7) implies that P max1≤l≤p kΛl k ≥ 2γ log(mp) ≤ 2(mp)−γ+1 , we have
p
2γ log(mp) + 1
2γ log(mp) + 1
kΛl k2 + 1 kΛl k + 1
≤
+
+
R := max kZl k ≤
√
√
1≤l≤p
p
p
p
p
28

with probability at least 1 − 2(mp)−γ+1 . We proceed with estimation of σ02 :=
looking at the (1, 1)-th and (2, 2)-th block of Zl Zl∗ , i.e.,
(Zl Zl∗ )1,1 =
(Zl Zl∗ )2,2 =

Pp

∗
l=1 E(Zl Zl )

by

1
1
(Λl Λ∗l − Im )2 +
(Λl Al − 1m vl∗ )(Λl Al − 1m vl∗ )∗ ,
2
p
mp
1
(e
el ee∗l ) ⊗ ((Λl Al − 1m vl∗ )∗ (Λl Al − 1m vl∗ )).
mp

Note that E(Λl Λ∗l − Im )2 = E((Λl Λ∗l )2 ) − Im . The i-th diagonal entry of (Λl Λ∗l )2 is | hal,i , vl i |4 =
| ml , diag(h̄l,i )vl |4 where al,i = Ml hl,i = diag(hl,i )ml and (5.42) implies E(| hal,i , vl i |4 ) ≤ 3
since diag(h̄l,i )vl is still a unit vector (note that diag(h̄l,i ) is unitary since hl,i is a column of Hl∗ ).
Therefore,
E(Λl Λ∗l − Im )2 = E((Λl Λ∗l )2 ) − Im  3Im − Im  2Im .
(5.19)
By using Lemma 5.18, we have
E(Λl Al − 1m vl∗ )(Λl Al − 1m vl∗ )∗

= E(Λl Al A∗l Λ∗l ) − 1m 1∗m
m(n − 1)Im
(n − 1)(mIm − 1m 1∗m )

.
=
m−1
m−1

(5.20)

With al,i = Ml hl,i and independence between {hl,i }m
i=1 and Ml , we have
E((Λl Al − 1m vl∗ )∗ (Λl Al − 1m vl∗ ))

= E(A∗l Λ∗l Λl Al ) − mvl vl∗
=E

m
X
i=1

=E

m
X
i=1

| hal,i , vl i |2 al,i a∗l,i

!

− mvl vl∗

∗

!

diag(hl,i ) | ml , diag(h̄l,i )vl |2 ml ml diag(h̄l,i )

 3mIn − mvl vl∗  3mIn

− mvl vl∗
(5.21)

where E | ml , diag(h̄l,i )vl |2 ml m∗l  3In follows from (5.41) and diag(hl,i ) diag(h̄l,i ) = In . By
using (5.19), (5.20), (5.21) and Lemma 5.10, σ02 is bounded above by
!

p 
p
X
X
E(Zl Zl∗ )1,1
0
∗
2
≤2
Zl Zl
σ0 := E
0
E(Zl Zl∗ )2,2
l=1
l=1

#
"
p
2
n−1
X
I
0
+
m
2
(m−1)p
p
≤ 2
3
el ee∗l ) ⊗ In
0
p (e
l=1


"
#


2
n−1
+ m−1
Im
0
1
n−1
p
≤ 2
≤6
+
.
3
p m−1
Inp
0
p

29

n
o
p
Conditioned on the event max1≤l≤p kΛl k ≥ 2γ log(mp) , applying (5.32) with t = γ log(np+m)
gives
p
nr 1
X
n−1p
(γ + 1) log(np + m),
+
≤ C0 max
Zl
p m−1
l=1
s
o 1
2γ log(mp)
log(mp) log(np + m) ≤
(γ + 1)
p
2
−γ
−γ+1 and it suffices to require the condition
with
 at least 1 − (np + m) − 2(mp)
 probability
n−1
γ 3 log4 (np + m) ≤ 41 .
C0 1p + m−1

5.3

Blind Calibration from multiple snapshots

Recall that Aw in (3.3) and the only difference from (3.1) is that here all Al are
εl = 0, Aw (excluding the last row) can be factorized into
√

pD √ 0
  diag(Av

A
1)
√
− √m · · ·
0
mIn
kx1 kIm · · ·
0
 0
p


kx1 k


.
.
.
.
.
.
.



..
..
..
..
..
..
..
A0 := 

..
  ...
.

diag(Av
)
p
A
0
· · · kxp kIm
√
0
· · · − √m
p
|
{z
}|
0
0
{z
}
Q
{z
|
mp×(np+m)
Z∈C

equal to A. If
···
···
..
.
···

0
0
..
.
√

mIn
kxp k

P








}

(5.22)
xi
n is the normalized x .
where vi = kx
∈
C
i
ik
Before we proceed to the main result in this section we need to introduce some notation. Let
ai be the i-th column of A∗ , which is a complex Gaussian random matrix; define Zi ∈ C(np+m)×p
to be a matrix whose columns consist of the i-th column of each block of Z ∗ , i.e.,
#
"
√1 hai , v1 iei · · ·
√1 hai , vp iei
p
p
Zi =
− √1m ee1 ⊗ ai · · · − √1m eep ⊗ ai
(np+m)×p

m
where “⊗” denotes Kronecker product and both {ei }m
el }pl=1 ∈ Rp are the standard
i=1 ∈ R and {e
orthonormal basis in Rm and Rp , respectively. In this way, the Zi are independently from one
another. By definition,
" 1 Pp
#
2 e e∗ − √1 e v ∗ (I ⊗ a a∗ )
|
ha
,
v
i
|
i
i
i
p
i
l
i
i
l=1
p
mp
Zi Zi∗ =
1
1
∗
(Ip ⊗ ai a∗i )ve∗i
− √mp
m Ip ⊗ (ai ai )

where
by

Pp

el
l=1 e

⊗ vl = v ∈ Cnp×1




v1
√
 
and v =  ...  with kvk = p. The expectation of Zi Zi∗ is given
vp

#
∗
√1 ei v ∗
e
e
−
i
i
mp
,
E(Zi Zi∗ ) =
1
1
ve∗i
− √mp
m Inp
"

C :=

m
X
i=1

30

"

#
√1 1m v ∗
I
−
m
mp
E(Zi Zi∗ ) =
.
1
v1∗m
Inp
− √mp

Our analysis depends on the mutual coherence of {vl }pl=1 . One cannot expect to recover all
{vl }pl=1 and D if all {vl }pl=1 are parallel to each other. Let G be the Gram matrix of {vl }pl=1 with
p ≤ n, i.e., G ∈ Cp×p and Gk,l = hvk , vl i , 1 ≤ k ≤ l ≤ p and in particular, Gl,l = 1, 1 ≤ l ≤ p. Its
eigenvalues are denoted by {λl }pl=1 with λp ≥ · · · ≥ λ1 ≥ 0. Basic linear algebra tells that
p
X

U GU ∗ = Λ,

λl = p,

(5.23)

l=1




v1∗
 
where U ∈ Cp×p is unitary and Λ = diag(λ1 , · · · , λp ). Let V =  ...  ∈ Cp×n , then there holds
v∗
√
√p
∗
∗
∗
Λ = U GU = U V (U V ) since G = V V . Here 1 ≤ kGk ≤ p and p ≤ kGkF ≤ p. In
√
particular, if hvk , vl i = δkl , then G = Ip ; if | hvk , vl i | = 1 for all 1 ≤ k, l ≤ p, then kGk = p and
kGkF = p.
We are now ready to state and prove the main result in this subsection.
Proposition 5.6. There holds
kZ ∗ Z − Ck =

m
X
i=1

Zi Zi∗ − C ≤

1
2

with probability at least 1 − 2m(np + m)−γ if




kGk kGk2F
n
1
C0 max
,
+
log2 (np + m) ≤
2
p
p
m
16(γ + 1)

(5.24)

and each al is i.i.d. complex Gaussian, i.e., ai ∼ √12 N (0, In ) + √i2 N (0, In ). In particular, if
√
G = Ip and kGkF = p, (5.24) becomes


n
1
1
+
.
(5.25)
C0
log2 (np + m) ≤
p m
16(γ + 1)
Remark 5.7. The proof of Theorem 3.5 follows exactly from that of Theorem 3.3 when Proposition 5.6 holds. Hence we just give a proof of Proposition 5.6.
Proof: [Proof of Proposition 5.6] Let Zi Zi∗ − E(Zi Zi∗ ) =: Zi,1 + Zi,2 , where Zi,1 and Zi,2 are
defined as
 1 Pp

(| hai , vl i |2 − 1)ei e∗i 0
l=1
p
Zi,1 :=
,
0
0
"
#
1
ei v ∗ (Ip ⊗ (ai a∗i − In ))
0
− √mp
Zi,2 :=
.
1
1
∗
− √mp
(Ip ⊗ (ai a∗i − In ))ve∗i
m Ip ⊗ (ai ai − In )
Estimation of k
p
X
l=1

Pm

i=1 Zi,1 k

Following from (5.23), we have
p

2

2

2

1/2

| hai , vl i | = kV ai k = kU V ai k = kΛ
31

−1/2

Λ

1X
2
λl ξi,l
U V ai k =
2
2

l=1

where Λ−1/2 U V is a p×n matrix with orthonormal
rows and hence each ξi,l is Rayleigh distributed.
√
(We say ξ is Rayleigh distributed if ξ = X 2 + Y 2 where both X and Y are standard real Gaussian
variables.)
Due to the simple form of Zi,1 , it is easy to see from Bernstein’s inequality for scalar random
variables (See Proposition 5.16 in [42]) that
m
X
i=1

p

Zi,1

1X
≤ max
| hai , vl i |2 − 1 ≤ C0 max
1≤i≤m p
l=1

with probability 1 − me−t . Here
Therefore,
p
X

Var

1
p

2
λl (ξi,l

l=1

2
max (λl |ξi,l
1≤l≤p

Pp

2
l=1 | hai , vl i |

!

− 1)

− 1|)ψ1

≤

2
Var(ξi,1

1
p

− 1)

tkGk
,
p

Pp

p
X
l=1

2
l=1 λl (ξi,l

√

tkGkF
p

− 1) where

(5.26)
Pp

l=1 λl

= p.

λ2l = c0 kGk2F ,

1≤l≤p

Pm

∗
i=1 E(Zi,2 Zi,2 )k

in order to bound k

Denote zi := (Ip ⊗ (ai a∗i − I))v and there holds,

Estimation of kZi,2 kψ1



≤ c1 max λl = c1 kGk.

Now we only need to find out kZi,2 kψ1 and σ02 = k

kZi,2 k ≤

−1 =



Pm

i=1 Zi,2 k.

1
1
1
1
max{kai k2 , 1} + √ kzi e∗i k =
max{kai k2 , 1} + √ kzi k.
m
mp
m
mp

For kai k2 , its ψ1 -norm is bounded by C0 n due to Lemma 5.14 and we only need to know kzi kψ1 :
v
u p
uX
√
∗
∗
| hai , vl i |2 kai k + p.
kzi k = k(Ip ⊗ (ai ai − In ))vk ≤ kIp ⊗ (ai ai )vk + kvk = t
l=1

Let u =

Pp

2
l=1 | hai , vl i |

=

√
( uv)ψ1

Pp

2
l=1 λl ξi,l

and v = kai k2 . The ψ1 -norm of kIp ⊗ (ai a∗i )vk satifies

√
= sup q −1 (E( uv)q )1/q ≤ sup q −1 (E uq )1/2q (E v q )1/2q
q≥1

≤

r

q≥1

√
sup q −1 (E uq )1/q sup q −1 (E v q )1/q ≤ C1 np
q≥1

q≥1

where the second inequality follows from the Cauchy-Schwarz inequality, kukψ1 ≤ C1 p, and kvkψ1 ≤
√
C1 n. Therefore, kzi kψ1 ≤ C2 np and there holds
kZi,2 kψ1 ≤ C0



r 

√ 
np
n
n
n
+√
+
.
≤ C0
m
mp
m
m

32

(5.27)

Pm
∗
Estimation of σ02 Note that σ02 :=
i=1 E(Zi,2 Zi,2 ) . Let Ai,1,1 and Ai,2,2 be the (1, 1)-th and
∗ respectively, i.e.,
(2, 2)-th block of Zi,2 Zi,2
Ai,1,1 =
Ai,2,2 =
=

1
k(Ip ⊗ (ai a∗i − In ))vk2 ei e∗i ,
mp
1
1
(Ip ⊗ (ai a∗i − In ))vv ∗ (Ip ⊗ (ai a∗i − In )) + 2 Ip ⊗ (ai a∗i − In )2
mp
m
1
1
[(ai a∗i − In )vk vl∗ (ai a∗i − In )]1≤k≤l≤p + 2 Ip ⊗ (ai a∗i − In )2 .
mp
m

(5.28)

∗ is a positive semi-definite matrix, Lemma 5.10 implies
Since Zi,2 Zi,2

σ02



m
X
Ai,1,1
0
E
.
≤2
0
Ai,2,2

(5.29)

i=1

So we only need to compute E(Ai,1,1 ) and E(Ai,2,2 ).
E k(Ip ⊗

(ai a∗i

2

− In ))vk =

p
X
l=1

E k(ai a∗i

where E(ai a∗i − In )2 = nIn . Now we have E Ai,1,1 =

2

− In )vl k = n

n
∗
m ei ei .

p
X
l=1

kvl k2 = np

For E(Ai,2,2 ), note that

E((ai a∗i − In )vk vl∗ (ai a∗i − In )) = Tr(vk vl∗ )In = hvl , vk i In = Gk,l In

which follows from (5.35). By (5.28), (5.29) and Lemma 5.10, there holds,

n



m n
X
ei e∗i
0
0
n
kGk
2
m
m
≤2
+
.
≤2
σ0 :≤ 2
n
n
1
Inp
0
0 1p G ⊗ In + m
p
m
mp G ⊗ In + m2 Inp
i=1

P
n
∗ ) is positive
each E(Zi,2 Zi,2
One lower bound of σ02 is σ02 ≥ k m
i=1 E(Ai,1,1 )k = m because

√ 
√
√
≤ log 2C0 √ n
≤ log(2C0 m) where R :=
semi-definite. Therefore, we have log σmR
0
n/m
pn
n
max1≤i≤m kZi,2 kψ1 ≤ CP
0( m +
m ) and m ≥ n.
Applying (5.34) to m
Z
i=1 i,2 with (5.29) and (5.27), we have
m
n n r n 
X
+
(t + log(np + m)) log(np + m),
Zi,2
≤ C0 max
m
m
i=1
s

o
kGk
n
+
(log(np + m) + t)
(5.30)
p
m
with
1 − e−t . By combining (5.30) with (5.26) and letting t = γ log(np + m), we have
Pm probability
∗
k i=1 Zi Zi − Ck ≤ 21 with probability 1 − 2m(np + m)−γ if
p


γ log(np + m)kGkF
kGk
1
n
1
≤ , C0 (γ + 1)
+
log2 (np + m) ≤
C0
p
4
p
m
16

n
o

kGk2F
n
1
or equivalently, C0 max kGk
+m
.
log2 (np + m) ≤ 16(γ+1)
p , p2
33

5.4

Proof of the spectral method

We provide the proof of the spectral method proposed in Section 2.3. The proof follows two steps: i)
we provide an error bound for the recovery under noise by using singular value/vector perturbation.
The error bound involves the second smallest singular value of S0 and the noise strength kδSk; ii)
we give a lower bound for σ2 (S0 ), which is achieved with the help of Proposition 5.2, 5.4 and 5.6
respectively for three different models.
The first result is a variant of perturbation theory for the singular vector corresponding to the
smallest singular value. A more general version can be found in [45, 36].
Lemma 5.8. Suppose S = S0 + δS where S0 is the noiseless part of S and δS is the noisy part.
Then
kδSk
kα0 ẑ − z0 k
(I − ẑ ẑ ∗ )z0
≤
=
.
min
kz
k
kz
k
[σ
(S
α0 ∈C
0
0
2 0 ) − kδSk]+
where σ2 (S0 ) is the second smallest singular value of S0 , z0 satisfies S0 z0 = 0, and ẑ is the right
singular vector with respect to the smallest singular value of S, i.e., the solution to (2.10).
 
s
Proof: By definition, there holds S0 z0 = 0 where z0 = 0 is the ground truth and also the
x0
right singular vector corresponding to the smallest singular value. Without loss of generality, we
assume kz0 k = kẑk = 1. For S, we denote its singular value decomposition as
S := S − σ1 (S)ûẑ ∗ + σ1 (S)ûẑ ∗
{z
} | {z }
|
B1

B0

where σ1 (S), û and ẑ are the smallest singular value/vectors of S. By definition, the vector ẑ is
also the solution to (2.10).
First note that I − ẑ ẑ ∗ = B1∗ (B1∗ )† where (B1∗ )† is the pseudo-inverse of B1 . Therefore, we
have
(I − ẑ ẑ ∗ )z0
kz0 k

= kz0 z0∗ B1∗ (B1∗ )† k = kz0 z0∗ (S0∗ + (δS)∗ − (B0 )∗ )(B1∗ )† k
= kz0 z0∗ (δS)∗ (B1∗ )† k ≤ kδSkk(B1∗ )† k ≤

kδAk
kδSk
≤
σ2 (S)
[σ2 (S0 ) − kδSk]+

where S0 z0 = 0 and (B0 )∗ (B1∗ )† = 0. And the last inequality follows from |σ2 (S0 ) − σ2 (S)| ≤ kδSk
and σ2 (S) ≥ [σ2 (S0 ) − kδSk]+ .
The second smallest singular value σ2 (S0 ) is estimated by using the proposition 5.2, 5.4 and 5.6
and the following fact:
Lemma 5.9. Suppose P is an invertible matrix, and A is a positive semi-definite matrix with the
dimension of its null space equal to 1. Then the second smallest singular value of P AP ∗ is nonzero
and satisfies
σ2 (P AP ∗ ) ≥ σ2 (A)σ12 (P )
where σ1 (·) and σ2 (·) denotes the smallest and second smallest singular values respectively.

34

Proof: The proof is very straightforward. Note that k(P AP ∗ )† k =
deficient. Also from the property of pseudo inverse, there holds
k(P AP ∗ )† k = k(P −1 )∗ A† P −1 k ≤ kA† kkP −1 k2 =

1
σ2 (P AP ∗ )

since A is rank-1

1
1
.
2
σ2 (A) σ1 (P )

Hence, σ2 (P AP ∗ ) ≥ σ2 (A)σ12 (P ).
Proof: [Proof of Theorem 3.9] Combined with Lemma 5.8, it suffices to estimate the lower
bound of the second smallest singular P
value of S0 for the proposed three models. We start with
∗
(a). From (5.1), we know that S0 S0 = pl=1 P Zl Zl∗ P ∗ where
"
#

 ∗
Λl
D kxk
0
(m+n)×m
√
∈ C(m+n)×(m+n) .
Zl :=
∈
C
,
P
:=
− √1m A∗l
mIn
0
P
From Proposition 5.2, we know that the second smallest eigenvalue of pl=1 Zl Zl∗ is at least
it is also rank-1 deficient. Applying Lemma 5.9 gives
!
p
X
√
p
Zl Zl∗ σ1 (P P ∗ ) ≥ (min{ m, dmin kxk})2
σ2 (S0∗ S0 ) ≥ σ2
2

p
2

and

l=1

and hence σ2 (S0 ) ≥

q

p
2

√
min{ m, dmin kxk}.

Since (b) and (c) are exactly the same, it suffices to show (b). From (5.15), we know that S0
can be factorized into A0 := QZP and Proposition 5.4 implies
"
#
√1 1m v ∗
I
−
m
1
mp
.
kZ ∗ Z − Ck ≤ , C := E(Z ∗ Z) =
1
− √mp
v1∗m
Inp
2
Therefore, S0∗ S0 = QZ ∗ P P ∗ ZQ∗ and σ2 (Z ∗ Z) ≥ 12 . Applying Lemma 5.9 leads to
σ2 (S0∗ S0 ) ≥ σ2 (Z ∗ P P ∗ Z)σ12 (Q) ≥ σ2 (Z ∗ Z)σ12 (P )σ12 (Q)

√ 2
√
m
1 2
pdmin ,
≥ xmin min
2
xmax
where xmin = min{kxl k} and xmax = max{kxl k}.

Appendix

S11 S12
, there
Lemma 5.10. For any Hermitian positive semi-definite matrix S with S :=
∗
S12
S22
holds,




S11 0
S11 S12
2
.
∗
S22
S12
0 S22


In other words, kSk ≤ 2 max{kS11 k, kS22 k}.
35

Lemma 5.11. Corollary 7.21 in [17]. Let a ∈ CM and ε = (ε1 , · · · , εM ) be a Rademacher sequence,
then for u > 0,


M
X
P
εj aj ≥ kaku ≤ 2 exp(−u2 /2).
j=1

For Gaussian and random Hadamard cases, the concentration inequalities are slightly different.
The following theorem is mostly due to Theorem 6.1 in [39] as well as due to [42].

Theorem 5.12. Consider a finite sequence of Zl of independent centered random matrices
PL with
dimension M1 × M2 . We assume that kZl k ≤ R and introduce the random matrix S := l=1 Zl .
Compute the variance parameter
σ02

L
L
n X
o
X
∗
= max k
E(Zl Zl )k, k
E(Zl∗ Zl )k ,
l=1

then for all t ≥ 0
kSk ≤ C0 max{σ0

p

(5.31)

l=1

t + log(M1 + M2 ), R(t + log(M1 + M2 ))}

(5.32)

with probability at least 1 − e−t where C0 is an absolute constant.
The concentration inequality is slightly different from Theorem 5.12 if kZl k is a sub-exponential
random variable. Here we are using the form in [24]. Before presenting the inequality, we introduce
the sub-exponential norm k · kψ1 of a matrix, defined as
kZkψ1 := inf {u : E[exp(kZk/u)] ≤ 2}.
u≥0

(5.33)

One can find more details about this norm and norm on Orlicz spaces in [42] and [40].
Theorem 5.13. For a finite sequence of independent M1 × M2 random matrices Zl with R :=
max1≤l≤L kZl kψ1 and σ02 as defined in (5.31), we have the tail bound on the operator norm of S,
!
√
p
LR
(t + log(M1 + M2 ))}
(5.34)
kSk ≤ C0 max{σ0 t + log(M1 + M2 ), R log
σ0
with probability at least 1 − e−t where C0 is an absolute constant.
The estimation of the ψ1 -norm of a sub-exponential random variable easily follows from the
following lemma.
Lemma 5.14 (Lemma 2.2.1 in [40]). Let z be a random variable which obeys P{|z| > u} ≤ ae−bu ,
then kzkψ1 ≤ (1 + a)/b.
Remark 5.15. A direct implication of Lemma 5.14 gives the ψ1 -norm of kak2 where a ∼ N (0, 12 In )+
iN (0, 21 In ) is a complex Gaussian random vector, i.e., (kak2 )ψ1 ≤ C0 n for the absolute constant
C0 .

36

Lemma 5.16. For a ∼ N (0, 12 In ) + iN (0, 21 In ), there holds
E(kak2 aa∗ ) =
∗

∗

E((a Xa)aa ) =

(n + 1)In ,
X + Tr(X)In ,

(5.35)

for any fixed x ∈ Cn and X ∈ Cn×n . In particular, we have
E | ha, xi |2 aa∗

(5.36)

4

= kxk2 In + xx∗ ,

= 2kxk ,

(5.37)

2

= nIn .

(5.38)

E | ha, xi |
∗

E(aa − In )

4

Lemma 5.17. Suppose that a ∈ Rn is a Rademacher sequence and for any fixed x ∈ Cn and
X ∈ Cn×n , there holds
E(kak2 aa∗ )
E((a∗ Xa)aa∗ )

= nIn ,

(5.39)

= Tr(X)In + X + X T − 2

n
X

Xkk Ekk

(5.40)

k=1

where Ekk is an n × n matrix with only one nonzero entry equal to 1 and at position (k, k). In
particular, setting X = xx∗ gives
E | ha, xi |2 aa∗
E | ha, xi |4

= kxk2 In + xx∗ + xx∗ − 2 diag(x) diag(x̄)  3kxk2 In ,
= 2kxk4 +

n
X

2

x2i

k=1

−2

n
X
k=1

|xk |4 ≤ 3kxk4 .

(5.41)
(5.42)

Proof: Since a is a Rademacher sequence, i.e, each ai takes ±1 independently with equal probability, this implies a2i =P1 and
kak2 = n. Therefore, E(kak2 aa∗ ) = n E(aa∗ ) = nIn . The (k, l)-th
P
entry of (a∗ Xa)aa∗ is ni=1 nj=1 Xij ai aj ak al .
1. If k = l,



n
n
n X
X
X
2

E(Xii |ai |2 |ak |2 ) = Tr(X)
Xij ai aj |ak |
=
E
i=1

i=1 j=1

where E(ai aj |ak |2 ) = 0 if i 6= j.
2. If k 6= l,



n X
n
X
Xij ai aj ak al  = Xkl E(|ak |2 |al |2 ) + Xlk E(|ak |2 |al |2 ) = Xkl + Xlk .
E
i=1 j=1

Hence, we have E((a∗ Xa)aa∗ ) = Tr(X)In + X + X T − 2
Lemma 5.18. There holds
E(ΛA(ΛA)∗ ) =

Pn

k=1 Xkk Ekk .

(n − 1)(mIm − 1m 1∗m )
+ 1m 1∗m
m−1
37

where A = HM , Λ = diag(Av) and v = (v1 , · · · , vn )T ∈ Cn is a deterministic unit vector.
H ∈ Cm×n is a random partial Fourier/Hadamard matrix with H ∗ H = mIn and m ≥ n, i.e., the
columns of H are uniformly sampled without replacement from an m × m DFT/Hadamard matrix;
M is a diagonal matrix with entries random sampled from ±1 with equal probability; moreover, we
assume M and H are independent from each other. In particular, if m = n = 1, E(ΛA(ΛA)∗ ) = 1.
Proof: We only prove the case when A is a random Fourier matrix since the Hadamard case is
essentially the same modulo very minor differences. By definition,
ΛAA∗ Λ∗ = ΛHH ∗ Λ∗ = diag(HM v)HH ∗ diag(HM v).
Let hi be the i-th column of H ∗ and the (i, j)-th entry of diag(HM v)HH ∗ diag(HM v) is
hhi , M vi hhj , M vi hhi , hj i where hu, vi = u∗ v. The randomness of ΛAA∗ Λ∗ comes from both H
and M and we first take the expectation with respect to M .
E(hhi , M vi hhj , M vi |H)

= h∗i E(M vv ∗ M |H)hj = h∗i diag(v) diag(v̄)hj

= (H diag(v) diag(v̄)H ∗ )i,j

where E(M vv ∗ M ) = diag(v) diag(v̄) follows from each entry in M being a Bernoulli random
variable. Hence, E((ΛAA∗ Λ∗ )i,j |H) = (H diag(v) diag(v̄)H ∗ )i,j · (HH ∗ )i,j .
Let uk be the k-th columnPof H and 1 ≤ k ≤ n and “⊙” denotes the Hadamard (pointwise)
product. So we have HH ∗ = nk=1 uk u∗k . There holds,
!
n
X
|vk |2 ūk ū∗k ⊙ HH ∗
E(ΛAA∗ Λ∗ |H) = (H diag(v) diag(v̄)H ∗ ) ⊙ HH ∗ =
k=1

=

n
X
k=1

=

|vk |2 diag(ūk )HH ∗ diag(uk ) =

X

1≤k6=l≤n

2

|vk |

n
n X
X
k=1 l=1

diag(ūk )ul u∗l diag(uk ) +

|vk |2 diag(ūk )ul u∗l diag(uk )

1m 1∗m

(5.43)

where the third equation follows from linearity of the Hadamard product and from
ūk ū∗k ⊙ HH ∗ = diag(ūk )HH ∗ diag(uk ).
The last one uses the fact that diag(ūk )uk = 1m if uk is a vector from the DFT matrix or
Hadamard matrix. By the property of conditional expectation, we know that E(ΛAA∗ Λ∗ ) =
E(E(ΛAA∗ Λ∗ )|H) and due to the linearity of expectation, it suffices to find out for k 6= l,
E(diag(ūk )ul u∗l diag(uk )) where uk and ul , by definition, are the k-th and l-th columns of H
which are sampled uniformly without replacement from an m × m DFT matrix F . Note that
(uk , ul ) is actually an ordered pair of random vectors sampled without replacement from columns
of F . Hence there are in total m(m − 1) different choices of diag(ūk )ul and
P(uk = φi , ul = φj ) =

1
,
m(m − 1)

38

∀1 ≤ i 6= j ≤ m, ∀1 ≤ k 6= l ≤ n

where φi is defined as the i-th column of an m × m DFT matrix F . Now we have, for any k 6= l,
E(diag(ūk )ul u∗l diag(uk )) =

=

=

X
1
diag(φ̄i )φj φ∗j diag(φi )
m(m − 1)
i6=j


X
1

diag(φ̄i )φj φ∗j diag(φi ) − m1m 1∗m 
m(m − 1)
1≤i,j≤m
!
m
X
1
∗
diag(φ̄i ) diag(φi ) − 1m 1m
m−1
i=1

=

mIm − 1m 1∗m
.
m−1

(5.44)

P
∗
∗ ∗
where diag(φ̄i )φi = 1m and m
i=1 φi φi = mIm . Now we return to E(ΛAA Λ ). By substituting (5.44) into (5.43), we get the desired formula:


X
E(ΛAA∗ Λ∗ ) = E(E(ΛAA∗ Λ∗ |H)) = E 
|vk |2 diag(uk )ul u∗l diag(uk ) + 1m 1∗m
1≤k6=l≤n

=

where

P

2
1≤k6=l≤n |vk |

1m 1∗m

mIm −
m−1

=n−

Pn

X

1≤k6=l≤n

2
k=1 |vk |

|vk |2 + 1m 1∗m =

= n − 1 follows from

(n − 1)(mIm − 1m 1∗m )
+ 1m 1∗m ,
m−1

Pm

2
k=1 |vk |

= 1.

References
[1] A. Ahmed, A. Cosse, and L. Demanet, A convex approach to blind deconvolution with
diverse inputs, in 2015 IEEE 6th International Workshop on Computational Advances in MultiSensor Adaptive Processing (CAMSAP), IEEE, 2015, pp. 5–8.
[2] A. Ahmed and L. Demanet, Leveraging diversity and sparsity in blind deconvolution, arXiv
preprint arXiv:1610.06098, (2016).
[3] A. Ahmed, F. Krahmer, and J. Romberg, Empirical chaos processes and blind deconvolution, arXiv preprint arXiv:1608.08370, (2016).
[4] A. Ahmed, B. Recht, and J. Romberg, Blind deconvolution using convex programming,
IEEE Transactions on Information Theory, 60 (2014), pp. 1711–1732.
[5] S. Bahmani and J. Romberg, Lifting for blind deconvolution in random mask imaging:
Identifiability and convex relaxation, SIAM Journal on Imaging Sciences, 8 (2015), pp. 2203–
2238.
[6] L. Balzano and R. Nowak, Blind calibration of sensor networks, in Proceedings of the 6th
International Conference on Information processing in sensor networks, ACM, 2007, pp. 79–88.
[7] L. Balzano and R. Nowak, Blind calibration of networks of sensors: Theory and algorithms,
in Networked Sensing Information and Control, Springer, 2008, pp. 9–37.
39

[8] C. Bilen, G. Puy, R. Gribonval, and L. Daudet, Convex optimization approaches for
blind sensor calibration using sparsity, IEEE Transactions on Signal Processing, 62 (2014),
pp. 4847–4856.
[9] V. Cambareri and L. Jacques, A greedy blind calibration method for compressed sensing
with unknown sensor gains, arXiv preprint arXiv:1610.02851, (2016).
[10] V. Cambareri and L. Jacques, A non-convex blind calibration method for randomised
sensing strategies, in Compressed Sensing Theory and its Applications to Radar, Sonar and
Remote Sensing (CoSeRa), 2016 4th International Workshop on, IEEE, 2016, pp. 16–20.
[11] V. Cambareri and L. Jacques, Through the haze: A non-convex approach to blind calibration for linear random sensing models, arXiv preprint arXiv:1610.09028, (2016).
[12] P. Campisi and K. Egiazarian, Blind Image Deconvolution: Theory and Applications, CRC
press, 2007.
[13] E. J. Candes, Y. C. Eldar, T. Strohmer, and V. Voroninski, Phase retrieval via
matrix completion, SIAM Review, 57 (2015), pp. 225–251.
[14] S. Curtis, J. Lim, and A. Oppenheim, Signal reconstruction from one bit of fourier transform phase, in 1984 IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), vol. 9, IEEE, 1984, pp. 487–490.
[15] M. A. Davenport and J. Romberg, An overview of low-rank matrix recovery from incomplete observations, IEEE Journal of Selected Topics in Signal Processing, 10 (2016), pp. 608–
622.
[16] J. R. Fienup, Reconstruction of an object from the modulus of its fourier transform, Optics
Letters, 3 (1978), pp. 27–29.
[17] S. Foucart and H. Rauhut, A Mathematical Introduction to Compressive Sensing, Springer,
2013.
[18] B. Friedlander and T. Strohmer, Bilinear compressed sensing for array self-calibration,
in 2014 48th Asilomar Conference on Signals, Systems and Computers, Asilomar, 2014.
[19] L. Gan, T. T. Do, and T. D. Tran, Fast compressive imaging using scrambled block
hadamard ensemble, in 2008 16th European Signal Processing Conference, IEEE, 2008, pp. 1–
5.
[20] R. Gribonval, G. Chardon, and L. Daudet, Blind calibration for compressed sensing by
convex optimization, in 2012 IEEE International Conference on Acoustics, Speech and Signal
Processing (ICASSP), IEEE, 2012, pp. 2713–2716.
[21] G. Harikumar and Y. Bresler, Perfect blind restoration of images blurred by multiple
filters: Theory and efficient algorithms, IEEE Transactions on Image Processing, 8 (1999),
pp. 202–219.

40

[22] A. Ito, A. C. Sankaranarayanan, A. Veeraraghavan, and R. G. Baraniuk, Blurburst: Removing blur due to camera shake using multiple images, ACM Trans. Graph., Submitted, 3 (2014).
[23] M. Kech and F. Krahmer, Optimal injectivity conditions for bilinear inverse problems with
applications to identifiability of deconvolution problems, SIAM Journal on Applied Algebra and
Geometry, 1 (2017), pp. 20–37.
[24] V. Koltchinskii et al., Von Neumann entropy penalization and low-rank matrix estimation,
The Annals of Statistics, 39 (2011), pp. 2936–2973.
[25] K. Lee, Y. Li, M. Junge, and Y. Bresler, Blind recovery of sparse signals from subsampled
convolution, IEEE Transactions on Information Theory, 63 (2017), pp. 802–821.
[26] X. Li, S. Ling, T. Strohmer, and K. Wei, Rapid, robust, and reliable blind deconvolution
via nonconvex optimization, arXiv preprint arXiv:1606.04933, (2016).
[27] Y. Li, K. Lee, and Y. Bresler, Identifiability in blind deconvolution with subspace or
sparsity constraints, IEEE Transactions on Information Theory, 62 (2016), pp. 4266–4275.
[28] Y. Li, K. Lee, and Y. Bresler, Optimal sample complexity for blind gain and phase calibration, IEEE Transactions on Signal Processing, 64 (2016), pp. 5549–5556.
[29] S. Ling and T. Strohmer, Blind deconvolution meets blind demixing: Algorithms and performance bounds, arXiv preprint arXiv:1512.07730, (2015).
[30] S. Ling and T. Strohmer, Self-calibration and biconvex compressive sensing, Inverse Problems, 31 (2015), p. 115002.
[31] J. Lipor and L. Balzano, Robust blind calibration via total least squares, in 2014 IEEE
International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, 2014,
pp. 4244–4248.
[32] V. I. Morgenshtern, E. Riegler, W. Yang, G. Durisi, S. Lin, B. Sturmfels, and H. Bölcskei, Capacity pre-log of noncoherent SIMO channels via Hironaka’s Theorem, IEEE Transactions on Information Theory, 59 (2013), pp. 4213–4229,
http://www.nari.ee.ethz.ch/commth/pubs/p/mrydlsb12.
[33] M. Pollefeys, R. Koch, and L. Van Gool, Self-calibration and metric reconstruction inspite of varying and unknown intrinsic camera parameters, International Journal of Computer
Vision, 32 (1999), pp. 7–25.
[34] Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, 2003.
[35] P. J. Shin, P. E. Larson, M. A. Ohliger, M. Elad, J. M. Pauly, D. B. Vigneron,
and M. Lustig, Calibrationless parallel imaging reconstruction based on structured low-rank
matrix completion, Magnetic Resonance in Medicine, 72 (2014), pp. 959–970.
[36] G. W. Stewart, Perturbation theory for the singular value decomposition, Technical Report
CS-TR 2539, University of Maryland, (September 1990).

41

[37] G. Tang and B. Recht, Convex blind deconvolution with random masks, in Computational
Optical Sensing and Imaging, Optical Society of America, 2014, pp. CW4C–1.
[38] L. Tong, G. Xu, B. Hassibi, and T. Kailath, Blind identification and equalization based
on second-order statistics: A frequency domain approach, IEEE Transactions on Information
Theory, 41 (1995), pp. 329–334.
[39] J. A. Tropp, User-friendly tail bounds for sums of random matrices, Foundations of Computational Mathematics, 12 (2012), pp. 389–434.
[40] A. Van der Vaart and J. Wellner, Weak Convergence and Empirical Processes: with
Applications to Statistics, Springer Series in Statistics, Springer-Verlag, New York, 1996.
[41] S. Verdú and S. Shamai, Spectral efficiency of cdma with random spreading, IEEE Transactions on Information theory, 45 (1999), pp. 622–640.
[42] R. Vershynin, Introduction to the non-asymptotic analysis of random matrices, in Compressed Sensing: Theory and Applications, Y. C. Eldar and G. Kutyniok, eds., Cambridge
University Press, 2012, ch. 5.
[43] L. Wang and Y. Chi, Blind deconvolution from multiple sparse inputs, IEEE Signal Processing Letters, 23 (2016), pp. 1384–1388.
[44] L. Wang, A. Singer, and Z. Wen, Orientation determination of cryo-em images using
least unsquared deviations, SIAM journal on imaging sciences, 6 (2013), pp. 2450–2483.
[45] P.-Å. Wedin, Perturbation bounds in connection with singular value decomposition, BIT Numerical Mathematics, 12 (1972), pp. 99–111.
[46] M. Wei, The perturbation of consistent least squares problems, Linear Algebra and its Applications, 112 (1989), pp. 231–245.
[47] A. J. Weiss and B. Friedlander, Eigenstructure methods for direction finding with sensor
gain and phase uncertainties, Circuits, Systems, and Signal Processing, 9 (1990), pp. 271–300.

42

