Regularized Gradient Descent: A Nonconvex Recipe
for Fast Joint Blind Deconvolution and Demixing ∗
Shuyang Ling† Thomas Strohmer‡
November 29, 2017

arXiv:1703.08642v2 [cs.IT] 28 Nov 2017

Abstract
We study the question of extracting a sequencePof functions {fi , gi }si=1 from observing
s
only the sum of their convolutions, i.e., from y = i=1 fi ∗ gi . While convex optimization
techniques are able to solve this joint blind deconvolution-demixing problem provably and
robustly under certain conditions, for medium-size or large-size problems we need computationally faster methods without sacrificing the benefits of mathematical rigor that come with
convex methods. In this paper we present a non-convex algorithm which guarantees exact
recovery under conditions that are competitive with convex optimization methods, with the
additional advantage of being computationally much more efficient. Our two-step algorithm
converges to the global minimum linearly and is also robust in the presence of additive
noise. While the derived performance bounds are suboptimal in terms of the informationtheoretic limit, numerical simulations show remarkable performance even if the number of
measurements is close to the number of degrees of freedom. We discuss an application of the
proposed framework in wireless communications in connection with the Internet-of-Things.

1

Introduction

The goal of blind deconvolution is the task of estimating two unknown functions from their
convolution. While it is a highly ill-posed bilinear inverse problem, blind deconvolution is
also an extremely important problem in signal processing [1], communications engineering [39],
imaging processing [5], audio processing [24], etc. In this paper, we deal with an even more
difficult and more general variation of the blind deconvolution problem, in which we have
to extract multiple convolved signals mixed together in one observation signal. This joint
blind deconvolution-demixing problem arises in a range of applications such as acoustics [24],
dictionary learning [2], and wireless communications [39].
We briefly discuss one such application in more detail. Blind deconvolution/demixing problems are expected to play a vital role in the future Internet-of-Things. The Internet-of-Things
will connect billions of wireless devices, which is far more than the current wireless systems can
technically and economically accommodate. One of the many challenges in the design of the
Internet-of-Things will be its ability to manage the massive number of sporadic traffic generating devices which are most of the time inactive, but regularly access the network for minor
updates with no human interaction [43]. This means among others that the overhead caused by
the exchange of certain types of information between transmitter and receiver, such as channel
estimation, assignment of data slots, etc, has to be avoided as much as possible [36, 27].
Focusing on the underlying mathematical challenges, we consider a multi-user communication scenario where many different users/devices communicate with a common base station, as
illustrated in Figure 1. Suppose we have s users and each of them sends a signal gi through an
unknown channel (which differs from user to user) to a common base station,. We assume that
∗

The authors acknowledge support from the NSF via grants dtra-dms 1322393 and DMS 1620455.
Courant Institute of Mathematical Sciences, New York University (Email: sling@cims.nyu.edu).
‡
Department of Mathematics, University of California at Davis (Email: strohmer@math.ucdavis.edu).
†

1

the i-th channel, represented by its impulse response fi , does not change during the transmission of the signal gi . Therefore fi acts as convolution operator, i.e., the signal transmitted by
the i-th user arriving at the base station becomes fi ∗ gi , where “∗” denotes convolution. The

User
1

⋮
User
𝑖

⋮
User
𝑠

𝑓$ : channel

𝐸𝑠𝑡𝑖𝑚𝑎𝑡𝑒 (𝑓$ , 𝑔$ )

𝑔$ : signal
𝑓$ ∗ 𝑔$

𝑓1 : channel 𝑓1 ∗ 𝑔1

𝑦 = ∑516$ 𝑓1 ∗ 𝑔1 + 𝑛

𝐸𝑠𝑡𝑖𝑚𝑎𝑡𝑒 (𝑓1 , 𝑔1 )

𝑔1 : signal
𝑓5 ∗ 𝑔5

𝑓5 : channel

Decoder

𝑔5 : signal

𝐸𝑠𝑡𝑖𝑚𝑎𝑡𝑒 (𝑓5 , 𝑔5 )

Figure 1: Single-antenna multi-user communication scenario without explicit channel estimation: Each of the s users sends a signal gi through an unknown channel fi to a common
P base station. The base station measures the superposition of all those signals, namely, y = si=1 fi ∗ gi
(plus noise). The goal is to extract all pairs of {(fi , gi )}si=1 simultaneously from y.
antenna at the base station, instead of receiving each individual component fi ∗ gi , is only able
to record the superposition of all those signals, namely,
y=

s
X

fi ∗ gi + n,

(1.1)

i=1

where n represents noise. We aim to develop a fast algorithm to simultaneously extract all pairs
{(fi , gi )}si=1 from y (i.e., estimating the channel/impulse responses fi and the signals gi jointly)
in a numerically efficient and robust way, while keeping the number of required measurements
as small as possible.

1.1

State of the art and contributions of this paper

A thorough theoretical analysis concerning the solvability of demixing problems via convex
optimization can be found in [26]. There, the authors derive explicit sharp bounds and phase
transitions regarding the number of measurements required to successfully demix structured
signals (such as sparse signals or low-rank matrices) from a single measurement vector. In
principle we could recast the blind deconvolution/demixing problem as the demixing of a sum
of rank-one matrices, see (2.3). As such, it seems to fit into the framework analyzed by McCoy
and Tropp. However, the setup in [26] differs from ours in a crucial manner. McCoy and Tropp
consider as measurement matrices (see the matrices Ai in (2.3)) full-rank random matrices, while
in our setting the measurement matrices are rank-one. This difference fundamentally changes
the theoretical analysis. The findings in [26] are therefore not applicable to the problem of
joint blind deconvolution/demixing. The compressive principal component analysis in [42] is
also a form of demixing problem, but its setting is only vaguely related to ours. There is a
large amount of literature on demixing problems, but the vast majority does not have a “blind
deconvolution component”, therefore this body of work is only marginally related to the topic
of our paper.

2

Blind deconvolution/demixing problems also appear in convolutional dictionary learning,
see e.g. [2]. There, the aim is to factorize an ensemble of input vectors into a linear combination
of overcomplete basis elements which are modeled as shift-invariant—the latter property is why
the factorization turns into a convolution. The setup is similar to (1.1), but with an additional
penalty term to enforce sparsity of the convolving filters. The existing literature on convolutional
dictionary learning is mainly focused on empirical results, therefore there is little overlap with
our work. But it is an interesting challenge for future research to see whether the approach
in this paper can be modified to provide a fast and theoretically sound solver for the sparse
convolutional coding problem.
There are numerous papers concerned with blind deconvolution/demixing problems in the
area of wireless communications [31, 37, 20]. But the majority of these papers assumes the
availability of multiple measurement vectors, which makes the problem significantly easier.
Those methods however cannot be applied to the case of a single measurement vector, which is
the focus of this paper. Thus there is essentially no overlap of those papers with our work.
Our previous paper [23] solves (1.1) under subspace conditions, i.e., assuming that both fi
and gi belong to known linear subspaces. This contributes to generalizing the pioneering work
by Ahmed, Recht, and Romberg [1] from the “single-user” scenario to the “multi-user” scenario.
Both [1] and [23] employ a two-step convex approach: first “lifting” [9] is used and then the
lifted version of the original bilinear inverse problems is relaxed into a semi-definite program.
An improvement of the theoretical bounds in [23] was announced in [29].
While the convex approach is certainly effective and elegant, it can hardly handle large-scale
problems. This motivates us to apply a nonconvex optimization approach [8, 21] to this blinddeconvolution-blind-demixing problem. The mathematical challenge, when using non-convex
methods, is to derive a rigorous convergence framework with conditions that are competitive
with those in a convex framework.
In the last few years several excellent articles have appeared on provably convergent nonconvex optimization applied to various problems in signal processing and machine learning, e.g.,
matrix completion [17, 16, 34], phase retrieval [8, 11, 33, 3], blind deconvolution [19, 4, 21],
dictionary learning [32], super-resolution [12] and low-rank matrix recovery [35, 41]. In this
paper we derive the first nonconvex optimization algorithm to solve (1.1) fast and with rigorous
theoretical guarantees concerning exact recovery, convergence rates, as well as robustness for
noisy data. Our work can be viewed as a generalization of blind deconvolution [21] (s = 1) to
the multi-user scenario (s > 1).
The idea behind our approach is strongly motivated by the nonconvex optimization algorithm for phase retrieval proposed in [8]. In this foundational paper, the authors use a two-step
approach: (i) Construct a good initial guess with a numerically efficient algorithm; (ii) Starting
with this initial guess, prove that simple gradient descent will converge to the true solution. Our
paper follows a similar two-step scheme. However, the techniques used here are quite different
from [8]. Like the matrix completion problem [7], the performance of the algorithm relies heavily
and inherently on how much the ground truth signals are aligned with the design matrix. Due
to this so-called “incoherence” issue, we need to impose extra constraints, which results in a different construction of the so-called basin of attraction. Therefore, influenced by [17, 34, 21], we
add penalty terms to control the incoherence and this leads to the regularized gradient descent
method, which forms the core of our proposed algorithm.
To the best of our knowledge, our algorithm is the first algorithm for the blind deconvolution/blind demixing problem that is numerically efficient, robust against noise, and comes with
rigorous recovery guarantees.

1.2

Notation

For a matrix Z, kZk denotes its operator norm and kZkF is its the Frobenius norm. For a
vector z, kzk is its Euclidean norm and kzk∞ is the `∞ -norm. For both matrices and vectors,
Z ∗ and z ∗ denote their complex conjugate transpose. z̄ is the complex conjugate of z. We
equip the matrix space CK×N with the inner product defined by hU , V i := Tr(U ∗ V ). For a

3

given vector z, diag(z) represents the diagonal matrix whose diagonal entries are z. For any
z ∈ R, let z+ = z+|z|
2 .

2

Preliminaries

Obviously, without any further assumption, it is impossible to solve (1.1). Therefore, we impose
the following subspace assumptions throughout our discussion [1, 23].
• Channel subspace assumption: Each finite impulse response fi ∈ CL is assumed to
have maximum delay spread K, i.e.,
 
h
fi = i .
0
Here hi ∈ CK is the nonzero part of fi and fi (n) = 0 for n > K.
• Signal subspace assumption: Let gi := Ci x̄i be the outcome of the signal x̄i ∈ CN
encoded by a matrix Ci ∈ CL×N with L > N , where the encoding matrix Ci is known
and assumed to have full rank1 .
Remark 2.1. Both subspace assumptions are common in various applications. For instance
in wireless communications, the channel impulse response can always be modeled to have finite
support (or maximum delay spread, as it is called in engineering jargon) due to the physical
properties of wave propagation [14]; and the signal subspace assumption is a standard feature
found in many current communication systems [14], including CDMA where Ci is known as
spreading matrix and OFDM where Ci is known as precoding matrix.
The specific choice of the encoding matrices Ci depends on a variety of conditions. In this
paper, we derive our theory by assuming that Ci is a complex Gaussian random matrix, i.e.,
each entry in Ci is i.i.d. CN (0, 1). This assumption, while sometimes imposed in the wireless
communications literature, is somewhat unrealistic in practice, due to the lack of a fast algorithm
to apply Ci and due to storage requirements. In practice one would rather choose Ci to be
something like the product of a Hadamard matrix and a diagonal matrix with random binary
entries. We hope to address such more structured encoding matrices in our future research. Our
numerical simulations (see Section 4) show no difference in the performance of our algorithm
for either choice.
Under the two assumptions above, the model actually has a simpler form in the frequency
domain. We assume throughout the paper that the convolution of finite sequences is circular
convolution2 . By applying the Discrete Fourier Transform DFT) to (1.1) along with the two
assumptions, we have
s

X
1
1
√ Fy =
diag(F hi )(F Ci x̄i ) + √ F n
L
L
i=1
where F is the L × L normalized unitary DFT matrix with F ∗ F = F F ∗ = IL . The noise
is
assumed to be additive white complex Gaussian noise with n ∼ CN (0, σ 2 d20 IL ) where d0 =
pP
s
∗
s
2
2
i=1 khi0 k kxi0 k , and {(hi0 , xi0 )}i=1 is the ground truth. We define di0 = khi0 xi0 kF and
assume without
loss of generality that khi0 k and kxi0 k are of the same norm, i.e., khi0 k =
√
kxi0 k = di0 , which is due to the scaling ambiguity3 . In that way, σ12 actually is a measure of
SNR (signal to noise ratio).
1

Here we use the conjugate x̄i instead of xi because it will simplify our notation in later derivations.
This circular convolution assumption can often be reinforced directly (for example in wireless communications
the use of a cyclic prefix in OFDM renders the convolution circular) or indirectly (e.g. via zero-padding). In
the first case replacing regular convolution by circular convolution does not introduce any errors at all. In the
latter case one introduces an additional approximation error in the inversion which is negligible, since it decays
exponentially for impulse responses of finite length [30].
3
Namely, if the pair (hi , xi ) is a solution, then so is (αhi , α−1 xi ) for any α 6= 0.
2

4

Let hi ∈ CK be the first K nonzero entries of fi and B ∈ CL×K be a low-frequency DFT
matrix (the first K columns of an L × L unitary DFT matrix). Then a simple relation holds,
B ∗ B = IK .

F fi = Bhi ,

We also denote Ai := F Ci and e := √1L F n. Due to the Gaussianity, Ai also possesses complex
Gaussian distribution and so does e. From now on, instead of focusing on the original model,
we consider (with a slight abuse of notation) the following equivalent formulation throughout
our discussion:
s
X
y=
diag(Bhi )Ai xi + e,
(2.1)
i=1
σ 2 d20

where e ∼ CN (0, L IL ). Our goal here is to estimate all {hi , xi }si=1 from y, B and {Ai }si=1 .
Obviously, this is a bilinear inverse problem, i.e., if all {hi }si=1 are given, it is a linear inverse
problem (the ordinary demixing problem) to recover all {xi }si=1 , and vice versa. We note that
there is a scaling ambiguity in all blind deconvolution problems that cannot be resolved by
any reconstruction method without further information. Therefore, when we talk about exact
recovery in the following, then this is understood modulo such a trivial scaling ambiguity.
Before proceeding to our proposed algorithm we introduce some notation to facilitate a more
convenient presentation of our approach. Let bl be the l-th column of B ∗ and ail be the l-th
column of A∗i . Based on our assumptions the following properties hold:
L
X

bl b∗l = IK ,

kbl k2 =

l=1

K
,
L

ail ∼ CN (0, IN ).

Moreover, inspired by the well-known lifting idea [9, 1, 6, 22], we define the useful matrix-valued
linear operator Ai : CK×N → CL and its adjoint A∗i : CL → CK×N by
Ai (Z) := {b∗l Zail }L
l=1 ,

A∗i (z) :=

L
X

zl bl a∗il = B ∗ diag(z)Ai

(2.2)

l=1

for each 1 ≤ i ≤ s under canonical inner product over CK×N . Therefore, (2.1) can be written
in the following equivalent form
s
X
y=
Ai (hi x∗i ) + e.
(2.3)
i=1

Hence, we can think of y as the observation vector obtained from taking linear measurements
with respect to a set of rank-1 matrices {hi x∗i }si=1 . In fact, with a bit of linear algebra (and
ignoring the noise term for the moment), the l-th entry of y in (2.3) equals the inner product
of two block-diagonal matrices:

  ∗

∗
0
···
0
bl a1l
0
···
0 +
* h1,0 x1,0
 0

h2,0 x∗2,0 · · ·
0 
bl a∗2l · · ·
0 

  0

yl = 
,
(2.4)


.
.
..  + el ,
..
..
..
..
..
.
.




.
.
.
.
.
.
.
.
0
0
· · · bl a∗sl
0
0
· · · hs0 x∗s0
|
{z
}
defined as X0

Ps

where yl = i=1 b∗l hi0 x∗i0 ail + el , 1 ≤ l ≤ L and X0 is defined as the ground truth matrix. In
other words, we aim to recover such a block-diagonal matrix X0 from L linear measurements
with block structure if e = 0.
By stacking all {hi }si=1 (and {xi }si=1 , {hi0 }si=1 , {xi0 }si=1 ) into a long column, we let
 




 
h1
x1
x1,0
h1,0
 


 


h :=  ...  , h0 :=  ...  ∈ CKs , x :=  ...  , x0 :=  ...  ∈ CN s .
(2.5)
hs

xs

hs0
5

xs0

We define H as a bilinear operator which maps a
matrix in CKs×N s , i.e.,

h1 x∗1
0
 0
h2 x∗2

H(h, x) :=  .
..
 ..
.
0
0

pair (h, x) ∈ CKs × CN s into a block diagonal
···
···
..
.

0
0
..
.

···

hs x∗s




 ∈ CKs×N s .


(2.6)

Let X := H(h, x) and X0 := H(h0 , x0 ) where X0 is the ground truth as illustrated in (2.4).
Define A(Z) : CKs×N s → CL as
s
X
A(Z) :=
Ai (Zi ),
(2.7)
i=1

where Z = blkdiag(Z1 , · · · , Zs ) and blkdiag is thePstandard MATLAB function to construct
block diagonal matrix. Therefore, A(H(h, x)) = si=1 Ai (hi x∗i ) and y = A(H(h0 , x0 )) + e.
The adjoint operator A∗ is defined naturally as

 ∗
A1 (z)
0
···
0
 0
A∗2 (z) · · ·
0 


Ks×N s
∗
,
(2.8)
A (z) :=  .
.
..  ∈ C
.
.
.
.
 .
.
.
. 
0

···

0

A∗s (z)

which is a linear map from CL to CKs×N s . To measure the approximation error of X0 given by
X, we define δ(h, x) as the global relative error:
qP
sP
s
∗
∗ 2
s
2 2
i=1 khi xi − hi0 xi0 kF
kX − X0 kF
i=1 δi di0
P
δ(h, x) :=
(2.9)
=
=
s
2 ,
kX0 kF
d0
i=1 di0
where δi := δi (hi , xi ) is the relative error within each component:
δi (hi , xi ) :=

khi x∗i − hi0 x∗i0 kF
.
di0

Note that δ and δi are functions of (h, x) and (hi , xi ) respectively and in most cases, we just
simply use δ and δi if no possibility of confusion exists.

2.1

Convex versus nonconvex approaches

As indicated in (2.4), joint blind deconvolution-demixing can be recast as the task to recover
a rank-s block-diagonal matrix from linear measurements. In general, such a low-rank matrix
recovery problem is NP-hard. In order to take advantage of the low-rank property of the
ground truth, it is natural to adopt convex relaxation by solving a convenient nuclear norm
minimization program, i.e.,
min

s
X

kZi k∗ ,

s.t.

i=1

s
X

Ai (Zi ) = y.

(2.10)

i=1

The question of when the solution of (2.10) yields exact recovery is first answered in our
previous work [23]. Late, [29, 15] have improved this result to the near-optimal bound L ≥
C0 s(K + N ) up to some log-factors where the main theoretical result is informally summarized
in the following theorem.
Theorem 2.2 (Theorem 1.1 in [15]). Suppose that Ai are L × N i.i.d. complex Gaussian
matrices and B is an L × K partial DFT matrix with B ∗ B = IK . Then solving (2.10) gives
exact recovery if the number of measurements L yields
L ≥ Cγ s(K + N ) log3 L
with probability at least 1 − L−γ where Cγ is an absolute scalar only depending on γ linearly.
6

While the SDP relaxation is definitely effective and has theoretic performance guarantees,
the computational costs for solving an SDP already become too expensive for moderate size
problems, let alone for large scale problems. Therefore, we try to look for a more efficient
nonconvex approach such as gradient descent, which hopefully is also reinforced by theory. It
seems quite natural to achieve the goal by minimizing the following nonlinear least squares
objective function with respect to (h, x)
2

F (h, x) := kA(H(h, x)) − yk =

s
X

2

Ai (hi x∗i )

−y

.

(2.11)

i=1

In particular, if e = 0, we write
F0 (h, x) :=

s
X

2

Ai (hi x∗i

−

hi0 x∗i0 )

.

(2.12)

i=1

As also pointed out in [21], this is a highly nonconvex optimization problem. Many of the
commonly used algorithms, such as gradient descent or alternating minimization, may not
necessarily yield convergence to the global minimum, so that we cannot always hope to obtain
the desired solution. Often, those simple algorithms might get stuck in local minima.

2.2

The basin of attraction

Motivated by several excellent recent papers of nonconvex optimization on various signal processing and machine learning problem, we propose our two-step algorithm: (i) Compute an
initial guess carefully; (ii) Apply gradient descent to the objective function, starting with the
carefully chosen initial guess. One difficulty of understanding nonconvex optimization consists
in how to construct the so-called basin of attraction, i.e., if the starting point is inside this
basin of attraction, the iterates will always stay inside the region and converge to the global
minimum. The construction of the basin of attraction varies for different problems [8, 3, 34].
For this problem, similar to [21], the construction follows from the following three observations.
Each of these observations suggests the definition of a certain neighborhood and the basin of
attraction is then defined as the intersection of these three neighborhood sets Nd ∩ Nµ ∩ N .
1. Ambiguity of solution: in fact, we can only recover (hi , xi ) up to a scalar since (αhi , α−1 xi )
and (hi , xi ) are both solutions for α 6= 0. From a numerical perspective, we want to avoid
the scenario when khi k → 0 and kxi k → ∞ while khi kkxi k is fixed, which potentially leads
to numerical instability. To balance both the norm of khi k and kxi k for all 1 ≤ i ≤ s, we
define
p
p
Nd := {{(hi , xi )}si=1 : khi k ≤ 2 di0 , kxi k ≤ 2 di0 , 1 ≤ i ≤ s},
which is a convex set.
2. Incoherence: the performance depends on how large/small the incoherence µ2h is, where µ2h
is defined by
LkBhi0 k2∞
.
µ2h := max
1≤i≤s
khi0 k2
The idea is that: the smaller the µ2h is, the better the performance is. Let us consider an
extreme case: if Bhi0 is highly sparse or spiky, we lose much information on those zero/small
entries and cannot hope to get satisfactory recovered signals. In other words, we need the
ground truth hi0 has “spectral flatness” and hi0 is not highly localized on the Fourier domain.
A similar quantity is also introduced in the matrix completion problem [7, 34]. The larger µ2h
is, the more hi0 is aligned with one particular row of B. To control the incoherence between
bl and hi , we define the second neighborhood,
p
√
Nµ := {{hi }si=1 : LkBhi k∞ ≤ 4 di0 µ, 1 ≤ i ≤ s},
(2.13)
where µ is a parameter and µ ≥ µh . Note that Nµ is also a convex set.
7

3. Close to the ground truth: we also want to construct an initial guess such that it is close
to the ground truth, i.e.,


khi x∗i − hi0 x∗i0 kF
s
≤ ε, 1 ≤ i ≤ s
(2.14)
N := {(hi , xi )}i=1 : δi =
di0
1
where ε is a predetermined parameter in (0, 15
].

Remark 2.3. To ensure δi ≤ ε, it suffices to ensure δ ≤ √εsκ where κ :=
because
s
1 X 2
ε2
2
δ
≤
δ
≤
i
sκ2
sκ2

max di0
min di0

≥ 1. This is

i=1

which implies max1≤i≤s δi ≤ ε.
Remark 2.4. When we say (h, x) ∈ Nd , Nµ or N , it means for all i = 1, . . . , s we have
(hi , xi ) ∈ Nd , Nµ or N respectively. In particular, (h0 , x0 ) ∈ Nd ∩ Nµ ∩ N where h0 and x0
are defined in (2.5).

2.3

Objective function and Wirtinger derivative

To implement the first two observations, we introduce the regularizer G(h, x), defined as the
sum of s components
G(h, x) :=

s
X

Gi (hi , xi ).

(2.15)

i=1

For each component Gi (hi , xi ), we let ρ ≥ d2 + 2kek2 , 0.9d0 ≤ d ≤ 1.1d0 , 0.9di0 ≤ di ≤ 1.1di0
for all 1 ≤ i ≤ s and




 X

L
h
L|b∗l hi |2 i
khi k2
kxi k2
,
Gi := ρ G0
G0
+ G0
+
2di
2di
8di µ2
l=1
|
{z
} |
{z
}
Nd

(2.16)

Nµ

where G0 (z) = max{z − 1, 0}2 . Here both d and {di }si=1 are data-driven and well approximated
by our spectral initialization procedure; and µ2 is a tuning parameter which could be estimated
if we assume a specific statistical model for the channel (for example, in the widely used Rayleigh
fading model, the channel coefficients are assumed to be complex Gaussian). The idea behind
Gi is quite straightforward though the formulation is complicated. For each Gi in (2.16), the
first two terms try to force the iterates to lie in Nd and the third term tries to encourage the
iterates to lie in Nµ . What about the neighborhood N ? A proper choice of the initialization
followed by gradient descent which keeps the objective function decreasing will ensure that the
iterates stay in N .
Finally, we consider the objective function as the sum of nonlinear least squares objective
function F (h, x) in (2.11) and the regularizer G(h, x),
Fe(h, x) := F (h, x) + G(h, x).

(2.17)

Note that the input of the function Fe(h, x) consists of complex variables but the output is
real-valued. As a result, the following simple relations hold
∂ Fe
∂ Fe
=
,
∂hi
∂ h̄i

∂ Fe
∂ Fe
=
.
∂ x̄i
∂xi

Similar properties also apply to both F (h, x) and G(h, x).

8

Therefore, to minimize this function, it suffices to consider only the gradient of Fe with
respect to h̄i and x̄i , which is also called Wirtinger derivative [8]. The Wirtinger derivatives of
F (h, x) and G(h, x) w.r.t. h̄i and x̄i can be easily computed as follows
∇Fhi = A∗i (A(X) − y) xi = A∗i (A(X − X0 ) − e) xi ,
∇Fxi =
∇Ghi =
∇Gxi =
where A(X) =

∗

− y)) hi = (A∗i (A(X − X0 ) − e))∗ hi ,




L
i
L X 0 L|b∗l hi |2
ρ h 0 khi k2
∗
G0
hi + 2
b
b
h
,
G0
i
l
l
2di
2di
4µ
8di µ2
l=1


ρ 0 kxi k2
G0
xi ,
2di
2di

(A∗i (A(X)

(2.18)
(2.19)
(2.20)
(2.21)

Ps

∗
i=1 Ai (hi xi )

and A∗ is defined in (2.8). In short, we denote




∇Fh1
∇Gh1




∇Feh := ∇Fh + ∇Gh , ∇Fh :=  ...  , ∇Gh :=  ...  .
∇Fhs
∇Ghs

(2.22)

Similar definitions hold for ∇Fex , ∇Fx and Gx . It is easy to see that ∇Fh = A∗ (A(X) − y)x
and ∇Fx = (A∗ (A(X) − y))∗ h.

3

Algorithm and Theory

3.1

Two-step algorithm

As mentioned before, the first step is to find a good initial guess (u(0) , v (0) ) ∈ CKs × CN s such
that it is inside the basin of attraction. The initialization follows from this key fact:
 

s
X
E(A∗i (y)) = E A∗i 
Aj (hj0 x∗j0 ) + e = hi0 x∗i0 ,
j=1

where we use B ∗ B =

PL

∗
l=1 bl bl

= IK , E(ail a∗il ) = IN and

E(A∗i Ai (hi0 x∗i0 )) =
E(A∗j Ai (hi0 x∗i0 )) =

L
X
l=1
L
X

bl b∗l hi0 x∗i0 E(ail a∗il ) = hi0 x∗i0 ,
bl b∗l hi0 x∗i0 E(ail a∗jl ) = 0,

∀j 6= i.

l=1

Therefore, it is natural to extract the leading singular value and associated left and right singular
vectors from each A∗i (y) and use them as (a hopefully good) approximation to (di0 , hi0 , xi0 ).
This idea leads to Algorithm 1, the theoretic guarantees of which are given in Section 6.5.
The second step of the algorithm is just to apply gradient descent to Fe with the initial guess
(0) (0)
(0)
{(ui , vi , di )}si=1 or (u(0) , v (0) , {di }si=1 ), where u(0) stems from stacking all ui into one long
vector4 .
Remark 3.1. For Algorithm 2, we can rewrite each iteration into
u(t) = u(t−1) − η∇Feh (u(t−1) , v (t−1) ),

v (t) = v (t−1) − η∇Fex (u(t−1) , v (t−1) ),

where ∇Feh and ∇Fex are in (2.22), and
 (t) 
u1
 .. 
(t)
u :=  .  ,

v (t)

(t)

 (t) 
v1
 .. 
:=  .  .
(t)

us

vs

4

It is clear that instead of gradient descent one could also use a second-order method to achieve faster
convergence at the tradeoff of increased computational cost per iteration. The theoretical convergence analysis
for a second-order method will require a very different approach from the one developed in this paper.

9

Algorithm 1 Initialization via spectral method and projection
1: for i = 1, 2, . . . , s do
2:
Compute A∗i (y).
3:
Find the leading singular value, left and right singular vectors of A∗i (y), denoted by
(di , ĥi0 , x̂i0 ).
4:
Solve the following optimization problem for 1 ≤ i ≤ s:
p
p
√
(0)
ui := argminz∈CK kz − di ĥi0 k2 s.t. LkBzk∞ ≤ 2 di µ.
√
(0)
Set vi = di x̂i0 .
6: end for
(0) (0)
7: Output: {(ui , vi , di )}si=1 or (u(0) , v (0) , {di }si=1 ).

5:

Algorithm 2 Wirtinger gradient descent with constant stepsize η
1:
2:
3:
4:
5:
6:
7:

3.2

Initialization: obtain (u(0) , v (0) , {di }si=1 ) via Algorithm 1.
for t = 1, 2, . . . , do
for i = 1, 2, . . . , s do
(t)
(t−1)
(t−1) (t−1)
ui = ui
− η∇Fehi (ui
, vi
),
(t)
(t−1)
(t−1) (t−1)
e
vi = vi
− η∇Fxi (ui
, vi
),
end for
end for

Main results

Our main findings are summarized as follows: Theorem 3.2 shows that the initial guess given
by Algorithm 1 indeed belongs to the basin of attraction. Moreover, di also serves as a good
approximation of di0 for each i. Theorem 3.3 demonstrates that the regularized Wirtinger
gradient descent will guarantee the linear convergence of the iterates and the recovery is exact
in the noisefree case and stable in the presence of noise.
Theorem 3.2. The initialization obtained via Algorithm 1 satisfies
\ 1
\
1
√ Nµ N √2ε
(u(0) , v (0) ) ∈ √ Nd
5 sκ
3
3

(3.1)

and
0.9di0 ≤ di ≤ 1.1di0 ,

0.9d0 ≤ d ≤ 1.1d0 ,

(3.2)

holds with probability at least 1 − L−γ+1 if the number of measurements satisfies
L ≥ Cγ+log(s) (µ2h + σ 2 )s2 κ4 max{K, N } log2 L/ε2 .

(3.3)

1
Here ε is any predetermined constant in (0, 15
], and Cγ is a constant only linearly depending
on γ with γ ≥ 1.

Theorem 3.3. Starting with the initial value z (0) := (u(0) , v (0) ) satisfying (3.1),the Algorithm 2
creates a sequence of iterates (u(t) , v (t) ) which converges to the global minimum linearly,
√
εd0
(1 − ηω)t/2 + 60 skA∗ (e)k
kH(u(t) , v (t) ) − H(h0 , x0 )kF ≤ √
2sκ2

(3.4)

with probability at least 1 − L−γ+1 where ηω = O((sκd0 (K + N ) log2 L)−1 ) and
s
γs(K + N )(log2 L)
kA∗ (e)k ≤ C0 σd0
L
if the number of measurements L satisfies
L ≥ Cγ+log(s) (µ2 + σ 2 )s2 κ4 max{K, N } log2 L/ε2 .
10

(3.5)

Remark 3.4. Our previous work [23] shows that the convex approach via semidefinite programming (see (2.10)) requires L ≥ C0 s2 (K + µ2h N ) log3 (L) to ensure exact recovery. Later, [15]
improves this result to the near-optimal bound L ≥ C0 s(K + µ2h N ) up to some log-factors. The
difference between nonconvex and convex methods lies in the appearance of the condition number
κ in (3.5). This is not just an artifact of the proof—empirically we also observe that the value
of κ affects the convergence rate of our nonconvex algorithm, see Figure 5.
Remark 3.5. Our theory suggests s2 -dependence for the number of measurements L, although
numerically L in fact depends on s linearly, as shown in Section 4. The reason for s2 -dependence
will be addressed in details in Section 5.2.
Remark 3.6. In the theoretical analysis, we assume that Ai (or equivalently Ci ) is a Gaussian
random matrix. Numerical simulations suggest that this assumption is clearly not necessary.
For example, Ci may be chosen to be a Hadamard-type matrix which is more appropriate and
favorable for communications.
Remark 3.7. If e = 0, (3.4) shows that (u(t) , v (t) ) converges to the ground truth at a linear
rate. On the other hand, if noise exists, (u(t) , v (t) ) is guaranteed to converge to a point within a
small neighborhood of (h0 , x0 ). More importantly, if the number of measurements L gets larger,
kA∗ (e)k decays at the rate of O(L−1/2 ).

4

Numerical simulations

In this section we present a range of numerical simulations to illustrate and complement different
aspects of our theoretical framework. We will empirically analyze the number of measurements
needed for perfect joint deconvolution/demixing to see how this compares to our theoretical
bounds. We will also study the robustness for noisy data. In our simulations we use Gaussian
encoding matrices, as in our theorems. But we also try more realistic structured encoding
matrices, that are more reminiscent of what one might come across in wireless communications.
While Theorem 3.3 says that the number of measurements L depends quadratically on the
number of sources s, numerical simulations suggest near-optimal performance. Figure 2 demonstrates that L actually depends linearly on s, i.e., the boundary between success (white) and
failure (black) is approximately a linear function of s. In the experiment, K = N = 50 are
fixed, all Ai are complex Gaussians and all (hi , xi ) are standard complex Gaussian vectors.
For each pair of (L, s), 25 experiments are performed and we treat the recovery as a success
if

kX̂−X0 kF
kX0 kF

≤ 10−3 . For our algorithm, we use backtracking to determine the stepsize and the

iteration stops either if kA(H(h(t+1) , x(t+1) ) − H(h(t) , x(t) ))k < 10−6 kyk or if the number of
iterations reaches 500. The backtracking is based on the Armijo-Goldstein condition [25]. The
1
initial stepsize is chosen to be η = K+N
. If Fe(z (t) − η∇Fe(z (t) )) > Fe(z (t) ), we just divide η by
two and use a smaller stepsize.
We see from Figure 2 that the number of measurements for the proposed algorithm to
succeed not only seems to depend linearly on the number of sensors, but it is actually rather
close to the information-theoretic limit s(K + N ). Indeed, the green dashed line in Figure 2,
which represents the empirical boundary for the phase transition between success and failure
corresponds to a line with slope about 32 s(K + N ). It is interesting to compare this empirical
performance to the sharp theoretical phase transition bounds one would obtain via convex
optimization [10, 26]. Considering the convex approach based on lifting in [23], we can adapt
the theoretical framework in [10] to the blind deconvolution/demixing setting, but with one
modification. The bounds in [10] rely on Gaussian widths of tangent cones related to the
measurement matrices Ai . Since simply analytic formulas for these expressions seem to be
out of reach for the structured rank-one measurement matrices used in our paper, we instead
compute the bounds for full-rank Gaussian random matrices, which yields a sharp bound of
about 3s(K + N ) (the corresponding bounds for rank-one sensing matrices will likely have a
constant larger than 3). Note that these sharp theoretical bounds predict quite accurately

11

the empirical behavior of convex methods. Thus our empirical bound for using a non-convex
methods compares rather favorably with that of the convex approach.

Number of measurements: L, from 100 to 1250

L vs. s, K = N = 50, Regularized GD, Gaussian
1

1250

0.9
0.8

1000

0.7
0.6

750

0.5
0.4

500

0.3
0.2

250

0.1

1

2

3

4

5

6

7

0

8

s: 1 to 8

Figure 2: Phase transition plot for empirical recovery performance under different choices of
(L, s) where K = N = 50 are fixed. Black region: failure; white region: success. The red solid
line depicts the number of degrees of freedom and the green dashed line shows the empirical
phase transition bound for Algorithm 2.
Similar conclusions can be drawn from Figure 3; there all Ai are in the form of Ai = F Di H
where F is the unitary L × L DFT matrix, all Di are independent diagonal binary ±1 matrices
and H is an L × N fixed partial deterministic Hadamard matrix. The purpose of using Di is to
enhance the incoherence between each channel so that our algorithm is able to tell apart each
individual signal and channel. As before we assume Gaussian channels, i.e., hi ∼ CN (0, IK )
Therefore, our approach does not only work for Gaussian encoding matrices Ai but also for
the matrices that are interesting to real-world applications, although no satisfactory theory has
been derived yet for that case. Moreover, due to the structure of Ai and B, fast transform
algorithms are available, potentially allowing for real-time deployment.
Plot: K = N = 50; Hadamard matrix
L = 512
L = 786
L = 1024
L = 1280

1

Empirical probability of success

0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

1

2

3

4

5

6

7

8

9

10

11

12

s: number of sources

Figure 3: Empirical probability of successful recovery for different pairs of (L, s) when K =
N = 50 are fixed.
Figure 4 shows the robustness of our algorithm under different levels of noise. We also run
25 samples for each level of SNR and different L and then compute the average relative error.
It is easily seen that the relative error scales linearly with the SNR and one unit of increase in
SNR (in dB) results in one unit of decrease in the relative error.
Theorem 3.3 suggests that the performance and convergence rate actually depend on the
di0
condition number of X0 = H(h0 , x0 ), i.e., on κ = max
min di0 where di0 = khi0 kkxi0 k. Next we
12

K = N = 64, s = 6, Gaussian

K = N = 64, s = 6, Hadamard

10

10
L = 2s(K+N)
L = 4s(K+N)

Average Relative Error of 10 Samples(dB)

Average Relative Error of 10 Samples(dB)

L = 2s(K+N)
L = 4s(K+N)
0

−10

−20

−30

−40

−50

−60

0

5

10

15

20

25

30

35

40

45

0

−10

−20

−30

−40

−50

−60

50

0

5

10

15

20

SNR(dB)

25

30

35

40

45

50

SNR(dB)

Figure 4: Relative error vs. SNR (dB): SNR = 20 log10



kyk
kek



.

demonstrate that this dependence on the condition number is not an artifact of the proof, but is
indeed also observed empirically. In this experiment, we let s = 2 and set for the first component
d1,0 = 1 and for the second one d2,0 = κ for κ ∈ {1, 2, 5}. Here, κ = 1 means that the received
signals of both sensors have equal power, whereas κ = 5 means that the signal received from the
second sensor is considerably stronger. The initial stepsize is chosen as η = 1, followed by the
backtracking scheme. Figure 5 shows how the relative error decays with respect to the number
of iterations t under different condition number κ and L.
The larger κ is, the slower the convergence rate is, as we see from Figure 5. This may result
from two reasons: our spectral initialization may not be able to give a good initial guess for those
weak components; moreover, during the gradient descent procedure, the gradient directions for
the weak components could be totally dominated/polluted by the strong components. Currently,
we still have no effective way of how to deal with this issue of slow convergence when κ is not
small. We have to leave this topic for future investigations.
K = N = 64, s = 2, L = 1024, Gaussian

K = N = 64, s = 2, L = 512, Gaussian
0

0

−1
−0.5

−1

Relative Error in log10

Relative Error in log10

−2

−1.5

−2

−2.5

−3

−4

−5

−6

−7

−8
−3

−3.5

kappa = 1
kappa = 2
kappa = 5
0

10

20

kappa = 1
kappa = 2
kappa = 5

−9

30

40

50

60

70

80

90

−10

100

Number of iterations

0

10

20

30

40

50

60

70

80

90

100

Number of iterations

Figure 5: Relative error vs. number of iterations t.

5

Convergence analysis

Our convergence analysis relies on the following four conditions where the first three of them
are local properties. We will also briefly discuss how they contribute to the proof of our main
theorem. Note that our previous work [21] on blind deconvolution is actually a special case
(s = 1) of (2.1). The proof of Theorem 3.3 follows in part the main ideas in [21]. The readers
may find the technical parts of [21] and this manuscript share many similarities. However,
13

there are also important differences. After all, we are now dealing with a more complicated
problem where the ground truth matrix X0 and measurement matrices are both rank-s blockdiagonal matrices, as shown in (2.4), instead of rank-1 matrices in [21]. The key is to understand
the properties of the linear operator A applying to different types of block-diagonal matrices.
Therefore, many technical details are much more involved while on the other hand, some of
results in [21] can be used directly. During the presentation, we will clearly point out both the
similarities to and differences from [21].

5.1

Four key conditions

s(K+N ) and ∇F
e(z) :=
Condition
"
# 5.1. Local regularity condition: Let z := (h, x) ∈ C
∇Feh (z)
∈ Cs(K+N ) , then
∇Fex (z)
k∇Fe(z)k2 ≥ ω[Fe(z) − c]+
(5.1)

for z ∈ Nd ∩ Nµ ∩ N where ω =

d0
7000

and c = kek2 + 2000skA∗ (e)k2 .

We will prove Condition 5.1 in Section 6.3. Condition 5.1 states that Fe(z) = 0 if k∇Fe(z)k =
0 and e = 0, i.e., all the stationary points inside the basin of attraction are global minima.
Condition 5.2. Local smoothness condition: Let z = (h, x) and w = (u, v) and there
holds
kFe(z + w) − Fe(z)k ≤ CL kwk
(5.2)
for z + w and z inside Nd ∩ Nµ ∩ N where CL ≈ O(d0 sκ(1 + σ 2 )(K + N ) log2 L) is the Lipschitz
constant of Fe over Nd ∩ Nµ ∩ N . The convergence rate is governed by CL .
The proof of Condition 5.2 can be found in Section 6.4.
Condition 5.3. Local restricted isometry property: Denote X = H(h, x) and X0 =
H(h0 , x0 ). There holds
3
2
kX − X0 k2F ≤ kA(X − X0 )k2 ≤ kX − X0 k2F
3
2
uniformly all for (h, x) ∈ Nd ∩ Nµ ∩ N .

(5.3)

Condition 5.3 will be proven in Section 6.2. It says that the convergence of the objective
function implies the convergence of the iterates.
Remark 5.4 (Necessity of inter-user incoherence). Although Condition 5.3 is seemingly
the same as the one in our previous work [21], it is indeed very different. Recall that A is
a linear operator acting on block-diagonal matrices and its output is the sum of s different
components involving Ai . Therefore, the proof of Condition 5.3 heavily depends on the inter-user
incoherence whereas this notion of incoherence is not needed at all for the single-user scenario.
At the beginning of Section 2, we discuss the choice of Ci (or Ai ). In order to distinguish
one user from another, it is essential to use sufficiently different5 encoding matrices Ci (or
Ai ). Here the independence and Gaussianity of all Ci (or Ai ) guarantee that kPTi A∗i Aj PTj k
is sufficiently small for all i 6= j where Ti is defined in (6.1). It is a key element to ensure the
validity of Condition 5.3 which is also an important component to prove Condition 5.1. On the
other hand, due to the recent progress on this joint deconvolution and demixing problem, one is
also able to prove a local restricted isometry property with tools such as bounding the suprema
of chaos processes [15] by assuming {Ai }si=1 as Gaussian matrices.
Condition 5.5. Robustness condition: Let ε ≤

1
15

be a predetermined constant. We have

kA∗ (e)k = max kA∗i (e)k ≤
1≤i≤s

where e ∼ CN (0,
5

σ 2 d20
L )

εd0
√
,
10 2sκ

if L ≥ Cγ κ2 s2 (K + N )/ε2 .

Suppose all Ci are the same, there is no hope to recover all pairs of {(hi , xi )}si=1 simultaneously.

14

(5.4)

We will prove Condition 5.5 in Section 6.5. We now extract one useful result based on
Conditions 5.3 and 5.5. From these two conditions, we are able to produce a good approximation
of F (h, x) for all (h, x) ∈ Nd ∩ Nµ ∩ N in terms of δ in (2.9). For (h, x) ∈ Nd ∩ Nµ ∩ N , the
following inequality holds
2 2 2
εδd2
3
εδd2
δ d0 − √ 0 + kek2 ≤ F (h, x) ≤ δ 2 d20 + √ 0 + kek2 .
3
2
5 sκ
5 sκ

(5.5)

Note that (5.5) simply follows from
F (h, x) = kA(X − X0 )k2F − 2 Re(hX − X0 , A∗ (e)i) + kek2 .
Note that (5.3) implies
cross-term,

2 2 2
3 δ d0

≤ kA(X − X0 )k2F ≤

3 2 2
2 δ d0 .

Thus it suffices to estimate the

| Re(hX − X0 , A∗ (e)i)| ≤ kA∗ (e)kkX − X0 k∗ = kA∗ (e)k

s
X

khi x∗i − hi0 x∗i0 k∗

i=1

≤

√

2kA∗ (e)k

s
X

khi x∗i − hi0 x∗i0 kF

i=1

≤

√

2skA∗ (e)kkX − X0 kF ≤

εδd20
√
10 sκ

(5.6)

where k · k∗ and k · k are a pair of dual norms and kA∗ (e)k comes from (5.4).

5.2

Outline of the convergence analysis

For the ease of proof, we introduce another neighborhood:


2 d2
ε
2
0
+ kek .
NFe = (h, x) : Fe(h, x) ≤
3sκ2
Moreover, another reason to consider NFe is based on the fact that gradient descent only allows
one to make the objective function decrease if the step size is chosen appropriately. In other
words, all the iterates z (t) generated by gradient descent are inside NFe as long as z (0) ∈ NFe .
On the other hand, it is crucial to note that the decrease of the objective function does
not necessarily imply the decrease of the relative error of the iterates. Therefore, we want to
construct an initial guess in N ∩ NFe so that z (0) is sufficiently close to the ground truth and
then analyze the behavior of z (t) .
In the rest of this section, we basically try to prove the following relation:
1
1
√ Nd ∩ √ Nµ ∩ N √2ε ⊂
5 sκ
3
3
|
{z
}
Initial guess

N ∩ NFe
| {z }

{z (t) }t≥0 in N ∩NFe

⊂

N d ∩ Nµ ∩ N 
|
{z
}

.

Key conditions hold over Nd ∩Nµ ∩N

Now we give a more detailed explanation of the relation above, which constitutes the main
structure of the proof:
1. We will show

√1 Nd ∩ √1 Nµ ∩ N √
2ε
3
3
5 sκ

⊂ N ∩ NFe in the proof of Theorem 3.3 in Section 5.3,

which is quite straightforward.
2. Lemma 5.6 explains why it holds that N ∩NFe ⊂ Nd ∩Nµ ∩N and where the s2 -bottleneck
comes from.
3. Lemma 5.8 implicitly shows that the iterates z (t) will remain in N ∩ NFe if the initial
guess z (0) is inside N ∩NFe and Fe(z (t) ) is monotonically decreasing (simply by induction).
Lemma 5.9 makes this observation explicit by showing that z (t) ∈ N ∩NFe implies z (t+1) :=
z (t) − η∇Fe(z (t) ) ∈ N ∩ NFe if the stepsize η obeys η ≤ C1L . Moreover, Lemma 5.9
guarantees sufficient decrease of Fe(z (t) ) in each iteration, which paves the road towards
the proof of linear convergence of Fe(z (t) ) and thus z (t) .
15

Remember that Nd and Nµ are both convex sets, and the purpose of introducing regularizers
Gi (hi , xi ) is to approximately project the iterates onto Nd ∩ Nµ . Moreover, we hope that once
the iterates are inside N and inside a sublevel subset NFe , they will never escape from NFe ∩ N .
Those ideas are fully reflected in the following lemma.
Lemma 5.6. Assume 0.9di0 ≤ di ≤ 1.1di0 and 0.9d0 ≤ d ≤ 1.1d0 . There holds NFe ⊂ Nd ∩ Nµ ;
moreover, under Conditions 5.3 and 5.5, we have NFe ∩ N ⊂ Nd ∩ Nµ ∩ N 9  .
10

If (h,
/ Nd ∩ Nµ , by the definition of G in (2.15), at least one component in G
 x)∈
2di0
exceeds ρG0 di . We have
Proof:


Fe(h, x) ≥ ρG0

2di0
di



2

2

≥ (d + 2kek )



2
2di0
−1
di

≥ (2/1.1 − 1)2 (d2 + 2kek2 )
1 2
ε2 d20
+ kek2 ,
≥
d0 + kek2 >
2
3sκ2
where ρ ≥ d2 + 2kek2 , 0.9d0 ≤ d ≤ 1.1d0 and 0.9di0 ≤ di ≤ 1.1di0 . This implies (h, x) ∈
/ NFe
and hence NFe ⊂ Nd ∩ Nµ .
Note that (h, x) ∈ Nd ∩ Nµ ∩ N if (h, x) ∈ NFe ∩ N . Applying (5.5) gives
2 2 2
εδd2
ε2 d20
δ d0 − √ 0 + kek2 ≤ F (h, x) ≤ Fe(h, x) ≤
+ kek2
3
3sκ2
5 sκ
which implies that δ ≤

9 √ε
10 sκ .

By definition of δ in (2.9), there holds

Ps
Ps
2
δi2 d2i0
81ε2
1
2
i=1
i=1 δi
P
≥
δ
=
≥
≥ 2 max δi2 ,
s
2
2
2
100sκ
sκ
sκ 1≤i≤s
i=1 di0
which gives δi ≤

9
10 ε

(5.7)

and (h, x) ∈ N 9 ε .
10

Remark 5.7. The s2 -bottleneck comes from (5.7). If δ ≤ ε is small, we cannot guarantee that
each δP
i is also smaller than ε. Just consider the simplest case when all di0 are the same: then
d20 = si=1 d2i0 = sd2i0 and there holds
s

ε2 ≥ δ 2 =

1X 2
δi .
s
i=1

Obviously, we cannot conclude that max δi ≤ ε but only say that δi ≤
require δ = O( √εs ) to ensure δi ≤ ε, which gives s2 -dependence in L.

√

sε. This is why we

Lemma 5.8. Denote z1 = (h1 , x1 ) and z2 = (h2 , x2 ). Let z(λ) := (1 − λ)z1 + λz2 . If z1 ∈ N
and z(λ) ∈ NFe for all λ ∈ [0, 1], we have z2 ∈ N .
Proof: Note that for z1 ∈ N ∩ NFe , we have z1 ∈ Nd ∩ Nµ ∩ N 9 ε which follows from
10
the second part of Lemma 5.6. Now we prove z2 ∈ N by contradiction. Let us suppose
that z2 ∈
/ N and z1 ∈ N . There exists z(λ0 ) := (h(λ0 ), x(λ0 )) ∈ N for some λ0 ∈ [0, 1]
kh x∗ −h x∗ k
such that max1≤i≤s i i di0i0 i0 F = . Therefore, z(λ0 ) ∈ NFe ∩ N and Lemma 5.6 implies
max1≤i≤s

khi x∗i −hi0 x∗i0 kF
di0

≤

9
10 ,

which contradicts max1≤i≤s

khi x∗i −hi0 x∗i0 kF
di0

= .

Lemma 5.9. Let the stepsize η ≤ C1L , z (t) := (u(t) , v (t) ) ∈ Cs(K+N ) and CL be the Lipschitz
constant of ∇Fe(z) over Nd ∩ Nµ ∩ N in (5.2). If z (t) ∈ N ∩ NFe , we have z (t+1) ∈ N ∩ NFe
and
Fe(z (t+1) ) ≤ Fe(z (t) ) − ηk∇Fe(z (t) )k2
(5.8)
(t+1)
(t)
(t)
where z
= z − η∇Fe(z ).
16

Remark 5.10. This lemma tells us that once z (t) ∈ N ∩ NFe , the next iterate z (t+1) = z (t) −
η∇Fe(z (t) ) is also inside N ∩ NFe as long as the stepsize η ≤ C1L . In other words, N ∩ NFe
is in fact a stronger version of the basin of attraction. Moreover, the objective function will
decay sufficiently in each step as long as we can control the lower bound of the ∇Fe, which is
guaranteed by the Local Regularity Condition 5.3.
Proof:

Let φ(τ ) := Fe(z (t) − τ ∇Fe(z (t) )), φ(0) = Fe(z (t) ) and consider the following quantity:
τmax := max{µ : φ(τ ) ≤ Fe(z (t) ), 0 ≤ τ ≤ µ},

where τmax is the largest stepsize such that the objective function Fe(z) evaluated at any point
over the whole line segment {z (t) − τ Fe(z (t) ), 0 ≤ τ ≤ τmax } is not greater than Fe(z (t) ). Now
we will show τmax ≥ C1L . Obviously, if k∇Fe(z (t) )k = 0, it holds automatically.
Consider k∇Fe(z (t) )k =
6 0 and assume τmax < 1 . First note that,
CL

d
φ(τ ) < 0 =⇒ τmax > 0.
dτ
By the definition of τmax , there holds φ(τmax ) = φ(0) since φ(τ ) is a continuous function w.r.t.
τ . Lemma 5.8 implies
{z (t) − τ ∇Fe(z (t) ), 0 ≤ τ ≤ τmax } ⊆ N ∩ NFe .
Now we apply Lemma 6.20, the modified descent lemma, and obtain
2
Fe(z (t) − τmax ∇Fe(z (t) )) ≤ Fe(z (t) ) − (2τmax − CL τmax
)kFe(z (t) )k2 ≤ Fe(z (t) ) − τmax kFe(z (t) )k2

where CL τmax ≤ 1. In other words, φ(τmax ) ≤ Fe(z (t) − τmax ∇Fe(z (t) )) < Fe(z (t) ) = φ(0) contradicts φ(τmax ) = φ(0).
Therefore, we conclude that τmax ≥ C1L . For any η ≤ C1L , Lemma 5.8 implies
{z (t) − τ ∇Fe(z (t) ), 0 ≤ τ ≤ η} ⊆ N ∩ NFe
and applying Lemma 6.20 gives
Fe(z (t) − η∇Fe(z (t) )) ≤ Fe(z (t) ) − (2η − CL η 2 )kFe(z (t) )k2 ≤ Fe(z (t) ) − ηkFe(z (t) )k2 .

5.3

Proof of Theorem 3.3

Combining all the considerations above, we now prove Theorem 3.3 to conclude this section.
Proof:

The proof consists of three parts:

Part I: Proof of z (0) := (u(0) , v (0) ) ∈ N ∩ NFe . From the assumption of Theorem 3.3,
\ 1
1
√ Nµ ∩ N √2ε .
z (0) ∈ √ Nd
5 sκ
3
3
First we show G(u(0) , v (0) ) = 0: for 0 ≤ i ≤ s and the definition of Nd and Nµ ,
(0)

(0)

L|b∗l ui |2
kui k2
2di0
L
16di0 µ2
2di0
≤
< 1,
≤
·
≤
< 1,
2
2
2di
3di
8di µ
8di µ
3L
3di
√
√
√
(0)
(0)
9
where kui k ≤ 2 √d3i0 , LkBui k∞ ≤ 4 √d3i0 µ and 10
di0 ≤ di ≤ 11
10 di0 . Therefore
(0)

G0

kui k2
2di

!

(0)

= G0

kvi k2
2di
17

!

(0)

= G0

L|b∗l ui |2
8di µ2

!
=0

for all 1 ≤ l ≤ L and G(u(0) , v (0) ) = 0.
For z (0) = (u(0) , v (0) ) ∈ N
δ(z (0) ) ≤

2ε
√
5 sκ

2ε
√
5 sκ

, we have δ(z (0) ) :=

√Ps

2 2
i=1 δi di0
d0

≤

2ε
√
.
5 sκ

By (5.5), there holds

and G(u(0) , v (0) ) = 0,

3
εδ(z (0) )d20
ε2 d20
√
Fe(u(0) , v (0) ) = F (u(0) , v (0) ) ≤ kek2 + δ 2 (z (0) )d20 +
≤ kek2 +
2
3sκ2
5 sκ
T
and hence z (0) = (u(0) , v (0) ) ∈ N NFe .
Part II: The linear convergence of the objective function Fe(z (t) ). Denote z (t) :=
(u(t) , v (t) ). Note that z (0) ∈ N ∩ NFe , Lemma 5.9 implies z (t) ∈ N ∩ NFe for all t ≥ 0 by
induction if η ≤ C1L . Moreover, combining Condition 5.1 with Lemma 5.9 leads to
h
i
Fe(z (t) ) ≤ Fe(z (t−1) ) − ηω Fe(z (t−1) ) − c , t ≥ 1
+

with c = kek2 + akA∗ (e)k2 and a = 2000s. Therefore, by induction, we have
i
h
i
h
i
ε2 d20
Fe(z (t) ) − c ≤ (1 − ηω) Fe(z (t−1) ) − c ≤ (1 − ηω)t Fe(z (0) ) − c ≤
(1 − ηω)t
3sκ2
+
+
+
h
i
 1 2 2

ε2 d2
ε2 d20
∗
2
where Fe(z (0) ) ≤ 3sκ02 + kek2 and Fe(z (0) ) − c ≤ 3sκ
≤
2 ε d0 − akA (e)k
2 . Now we
3sκ
+
+
h
i
conclude that Fe(z (t) ) − c converges to 0 linearly.
h

+

Part III: The linear convergence of the iterates (u(t) , v (t) ).
δ(z (t) ) :=

Denote

kH(u(t) , v (t) ) − H(h0 , x0 )kF
.
d0

Note that z (t) ∈ N ∩NFe ⊆ Nd ∩Nµ ∩N and over Nd ∩Nµ ∩N , there holds F0 (z (t) ) ≥ 23 δ 2 (z (t) )d20
which follows from Local RIP Condition in (5.3) and F0 (z (t) ) defined in (2.12). Moreover


Fe(z (t) ) − kek2 ≥ F0 (z (t) ) − 2 Re hA∗ (e), H(u(0) , v (0) ) − H(h0 , x0 )i
√
2 2 (t) 2
≥
δ (z )d0 − 2 2skA∗ (e)kδ(z (t) )d0
3
where G(z (t) ) ≥ 0 and the second inequality follows from (5.6). There holds
h
i
√
2 2 (t) 2
ε2 d20
δ (z )d0 − 2 2skA∗ (e)kδ(z (t) )d0 − akA∗ (e)k2 ≤ Fe(z (t) ) − c ≤
(1 − ηω)t
3
3sκ2
+
and equivalently,
√
3 2 ∗
δ(z )d0 −
kA (e)k
2

2

(t)

ε2 d20
≤
(1 − ηω)t +
2sκ2



3
9
a+
2
2



kA∗ (e)k2 .

Solving the inequality above for δ(z (t) ), we have
!
r
√
εd0
3 2
3
9
t/2
δ(z )d0 ≤ √
(1 − ηω) +
+
a+
kA∗ (e)k
2
2
2
2
2sκ
√
εd0
≤ √
(1 − ηω)t/2 + 60 skA∗ (e)k
(5.9)
2sκ2
q
Ps
(t) 2
(t) 2
where a = 2000s. Let d(t) :=
i=1 kui k kvi k for t ∈ Z≥0 . By (5.9) and triangle inequality,
√
we immediately obtain |d(t) − d0 | ≤ √εd0 2 (1 − ηω)t/2 + 60 skA∗ (e)k.
(t)

2sκ

18

6

Proof of the four conditions

This section is devoted to proving the four key conditions introduced in Section 5. The local
smoothness condition and the robustness condition are relatively less challenging to deal with.
The more difficult part is to show the local regularity condition and the local isometry property.
The key to solve those problems is to understand how the vector-valued linear operator A
in (2.7) behaves on block-diagonal matrices, such as H(h, x), H(h0 , x0 ) and H(h, x)−H(h0 , x0 ).
In particular, when s = 1, all those matrices become rank-1 matrices, which have been well
discussed in our previous work [21].
First of all, we define the linear subspace Ti ⊂ CK×N along with its orthogonal complement
for 1 ≤ i ≤ s as
Ti := {Zi ∈ CK×N : Zi = hi0 vi∗ + ui x∗i0 , ui ∈ CK , vi ∈ CN },

 


hi0 h∗i0
xi0 x∗i0
⊥
K×N
Ti :=
IK −
Zi IN −
: Zi ∈ C
di0
di0
√
where khi0 k = kxi0 k = di0 . In particular, hi0 x∗i0 ∈ Ti for all 1 ≤ i ≤ s.

(6.1)

The proof also requires us to consider block-diagonal matrices whose i-th block belongs to
Ti (or Ti⊥ ). Let Z = blkdiag(Z1 , · · · , Zs ) ∈ CKs×N s be a block-diagonal matrix and say Z ∈ T
if
T := {blkdiag({Zi }si=1 )|Zi ∈ Ti }
and Z ∈ T ⊥ if
T ⊥ := {blkdiag({Zi }si=1 )|Zi ∈ Ti⊥ }
where both T and T ⊥ are subsets in CKs×N s and H(h0 , x0 ) ∈ T.
Now we take a closer look at a special case of block-diagonal matrices, i.e., H(h, x) and
calculate its projection onto T and T ⊥ respectively and it suffices to consider PTi (hi x∗i ) and
PT ⊥ (hi x∗i ). For each block hi x∗i and 1 ≤ i ≤ s, there are unique orthogonal decompositions
i

hi := αi1 hi0 + h̃i ,

x := αi2 xi0 + x̃i ,

(6.2)

where hi0 ⊥ h̃i and xi0 ⊥ x̃i . It is important to note that αi1 = αi1 (hi ) =
hxi0 ,xi i
di0

hhi0 ,hi i
di0

and αi2 =

αi2 (xi ) =
and thus αi1 and αi2 are functions of hi and xi respectively. Immediately, we
have the following matrix orthogonal decomposition for hi x∗i onto Ti and Ti⊥ ,
hi x∗i − hi0 x∗i0 = (αi1 αi2 − 1)hi0 x∗i0 + αi2 h̃i x∗i0 + αi1 hi0 x̃∗i +
|
{z
}
belong to Ti

h̃i x̃∗i
| {z }

(6.3)

belongs to Ti⊥

where the first three components are in Ti while h̃i x̃∗i ∈ Ti⊥ .

6.1

Key lemmata

From the decomposition in (6.2) and (6.3), we want to analyze how kh̃i k, kx̃i k, αi1 and αi2
kh x∗ −h x∗ k
depend on δi = i i di0i0 i0 F if δi < 1. The following lemma answers this question, which can
be viewed as an application of singular value/vector perturbation theory [40] applied to rank-1
matrices. From the lemma below, we can see that if hi x∗i is close to hi0 x∗i0 , then PT ⊥ (hi x∗i ) is
i
in fact very small (of order O(δi2 di0 )).
Lemma 6.1. (Lemma 5.9 in [21]) Recall that khi0 k = kxi0 k =
1, we have the following useful bounds
|αi1 | ≤

khi k
,
khi0 k

|αi1 αi2 − 1| ≤ δi ,

19

√

di0 . If δi :=

khi x∗i −hi0 x∗i0 kF
di0

<

and

δi
δi2
δi
khi k, kx̃i k ≤
kxi k, kh̃i kkx̃i k ≤
di0 .
1 − δi
1 − δi
2(1 − δi )
√
√
√
√
T
Moreover,
if khi k ≤ 2 di0 and LkBhi k∞ ≤ 4µ di0 , i.e., hi ∈ Nd Nµ , we have LkBi h̃i k∞ ≤
√
6µ di0 .
kh̃i k ≤

Now we start to focus on several results related to the linear operator A.
Lemma 6.2. (Operator norm of A). For A defined in (2.7), there holds
p
kAk ≤ s(N log(N L/2) + (γ + log s) log L)

(6.4)

with probability at least 1 − L−γ .
Proof:

Note that Ai (Zi ) := {b∗l Zi ail }L
l=1 in (2.2). Lemma 1 in [1] implies
p
kAi k ≤ N log(N L/2) + γ 0 log L
0

with probability at least 1 − L−γ . By taking the union bound over 1 ≤ i ≤ s,
p
max kAi k ≤ N log(N L/2) + (γ + log s) log L
with probability at least 1 − sL−γ−log s ≥ 1 − L−γ .
For A defined in (2.7), applying the triangle inequality gives
v
u s
s
s
X
X
u X
√
kA(Z)k =
Ai (Zi ) ≤
kAi kkZi kF ≤ max kAi kts
kZi k2F = s max kAi kkZkF
i=1

i=1

1≤i≤s

1≤i≤s

i=1

where Z = blkdiag(Z1 , · · · , Zs ) ∈ CKs×N s . Therefore,
p
√
kAk ≤ s max kAi k ≤ s(N log(N L/2) + (γ + log s) log L)
1≤i≤s

with probability at least 1 − L−γ .
Lemma 6.3. (Restricted isometry property for A on T ). The linear operator A restricted
on T is well-conditioned, i.e.,
1
kPT A∗ APT − PT k ≤
(6.5)
10
where PT is the projection operator from CKs×N s onto T , given L ≥ Cγ s2 max{K, µ2h N } log2 L
with probability at least 1 − L−γ .
Remark 6.4. Here APT and PT A∗ are defined as
APT (Z) =

s
X

Ai (PTi (Zi )),

PT A∗ (z) = blkdiag(PT1 (A∗1 (z)), · · · , PTs (A∗s (z)))

i=1

respectively where Z is a block-diagonal matrix and z ∈ CL .
As shown in the remark above, the proof of Lemma 6.3 depends on the properties of both
PTi A∗i Ai PTi and PTi A∗i Aj PTj for i 6= j. Fortunately, we have already proven related results
in [23] which are written as follows:
Lemma 6.5 (Inter-user incoherence, Corollary 5.3 and 5.8 in [23]). There hold
kPTi A∗i Aj PTj k ≤

1
,
10s

∀i 6= j;

kPTi A∗i Ai PTi − PTi k ≤

1
,
10s

∀1 ≤ i ≤ s

with probability at least 1 − L−γ+1 if L ≥ Cγ s2 max{K, µ2h N } log2 L log(s + 1).
20

(6.6)

1
Note that kPTi A∗i Aj PTj k ≤ 10s
holds because of independence between each individual
random Gaussian matrix Ai . In particular, if s = 1, the inter-user incoherence kPTi A∗i Aj PTj k ≤
1
10s is not needed at all. With (6.6), it is easy to prove Lemma 6.3.

Proof of Lemma 6.3. For any block diagonal matrix Z = blkdiag(Z1 , · · · , Zs ) ∈ CKs×N s
and Zi ∈ CK×N ,
X
hAi PTi (Zi ), Aj PTj (Zj )i − kPT (Z)k2F
hZ, PT A∗ APT (Z) − PT (Z)i =
=

1≤i,j≤s
s
X

hZi , PTi A∗i Ai PTi (Zi ) − PTi (Zi )i +

i=1

X
hAi PTi (Zi ), Aj PTj (Zj )i.
i6=j

(6.7)
Using (6.6), the following two inequalities hold,
kZi k2F
,
10s
kZi kF kZj kF
|hAi PTi (Zi ), Aj PTj (Zj )i| ≤ kPTi A∗i Aj PTj kkZi kF kZj kF ≤
.
10s

|hZi , PTi A∗i Ai PTi (Zi ) − PTi (Zi )i| ≤ kPTi A∗i Ai PTi − PTi kkZi k2F ≤

After substituting both estimates into (6.7), we have
X kZi kF kZj kF
1
|hZ, PT A APT (Z) − PT (Z)i| ≤
≤
10s
10s
∗

1≤i,j≤s

s
X

!2
kZi kF

i=1

≤

kZk2F
.
10

Finally, we show how A behaves when applied to block-diagonal matrices X = H(h, x). In
particular, the calculations will be much simplified for the case s = 1.
Lemma 6.6. (A restricted on block-diagonal matrices with rank-1 blocks).
Consider X = H(h, x) and
2
σmax
(h, x) := max

1≤l≤L

s
X

|b∗l hi |2 kxi k2 .

(6.8)

i=1

Conditioned on (6.4), we have
q
4
2 (h, x)(K + N ) log L+8sσ 2 (h, x)(K +N ) log L, (6.9)
kA(X)k2 ≤ kXk2F +2 2skXk2F σmax
max
3
uniformly for any h ∈ CKs and x ∈ CN s with probability at least 1 −
P
L ≥ Cγ s(K + N ) log L. Here kXk2F = kH(h, x)k2F = si=1 khi k2 kxi k2 .

1
γ

exp(−s(K + N )) if

2 (h, x). Note that kA(X)k2
Remark 6.7. Here are a few more explanations and facts about σmax
is the sum of L sub-exponential 6 random variables, i.e.,

2

kA(X)k =

s
L X
X

2

b∗l hi x∗i ail

.

(6.10)

l=1 i=1
2 (h, x) corresponds to the largest expectation of all those components in kA(X)k2 .
Here σmax
2 (h, x), without loss of generality, we assume kx k = 1 for 1 ≤ i ≤ s and let h ∈ CKs
For σmax
i
P
be a unit vector, i.e., khk2 = si=1 khi k2 = 1. The bound

1
K
2
≤ σmax
(h, x) ≤
L
L

(6.11)

6
For the definition and properties of sub-exponential random variables, the readers can find all relevant
information in [38].

21

PL Ps
2 (h, x) ≥
∗
2
2
follows from Lσmax
l=1
i=1 |bl hi | = khk = 1.
2
Moreover, σmax (h, x) and σmax (h, x) are both Lipschitz functions w.r.t. h. Now we want to
determine their Lipschitz constants. First note that for kxi k = 1, σmax (h, x) equals
σmax (h, x) = max k(Is ⊗ b∗l )hk
1≤l≤L

where ⊗ denotes Kronecker product. Let u ∈ CKs be another unit vector and we have
|σmax (h, x) − σmax (u, x)| = max k(Is ⊗ b∗l )h − max k(Is ⊗ b∗l )uk
1≤l≤L

1≤l≤L

= max |k(Is ⊗
1≤l≤L

b∗l )hk

− k(Is ⊗ b∗l )uk|

≤ max k(Is ⊗ b∗l )(h − u)k ≤ kh − uk

(6.12)

1≤l≤L

where kIs ⊗ b∗l k = kbl k

q

K
L

2 (h, x),
< 1. For σmax

2
2
|σmax
(h, x) − σmax
(u, x)| ≤ (σmax (h, x) + σmax (u, x)) · |σmax (h, x) − σmax (u, x)|
2K
kh − uk ≤ 2kh − uk.
(6.13)
≤
L
P
Proof of Lemma 6.6. Without loss of generality, let kxi k = 1 and si=1 khi k2 = 1. It suffices
to prove f (h, x) ≤ 43 for all (h, x) ∈ CKs × CN s in (2.5) where f (h, x) is defined as
p
2 (h, x)(K + N ) log L − 8sσ 2 (h, x)(K + N ) log L.
f (h, x) := kA(X)k2 − 2 2sσmax
max

Part I: Bounds
of kA(X)k2 for any fixed (h, x). From (6.10), we already know that Y =
P
2L
kA(X)k2F = i=1 ci ξi2 where {ξi } are i.i.d. χ21 random variables and c = (c1 , · · · , c2L )T ∈ R2L .
More precisely, we can determine {ci }2L
i=1 as
s
X

2

b∗l hi x∗ ail

s

2
2
,
+ c2l ξ2l
= c2l−1 ξ2l−1

c2l−1 = c2l =

i=1

1X ∗ 2
|bl hi |
2
i=1

Ps

Ps

because i=1 b∗l hi x∗i ail ∼ CN 0, i=1 |b∗l hi | .
By the Bernstein inequality, there holds




t2
t
P(Y − E(Y ) ≥ t) ≤ exp −
∨ exp −
8kck2
8kck∞

2

(6.14)

where E(Y ) = kXk2F = 1. In order to apply the Bernstein inequality, we need to estimate kck2
and kck∞ as follows,
s

kck∞ =

X
1 2
1
max
|b∗l hi |2 = σmax
(h, x),
2 1≤l≤L
2
i=1

s
L
1X X ∗ 2
|bl hi |
kck22 =
2
l=1 i=1

2

1
≤
2

s X
L
X

!
|b∗l hi |2

i=1 l=1

max

1≤l≤L

s
X
i=1

1 2
|b∗l hi |2 ≤ σmax
(h, x).
2

Applying (6.14) gives

P(kA(X)k ≥ 1 + t) ≤ exp −
2

t2
2 (h, x)
4σmax




∨ exp −

t
2
4σmax (h, x)


.

In particular, by setting
t = g(h, x) := 2

p
2 (h, x)(K + N ) log L + 8sσ 2 (h, x)(K + N ) log L,
2sσmax
max

we have

P kA(X)k2 ≥ 1 + g(h, x) ≤ e−2s(K+N )(log L) .
So far, we have shown that f (h, x) ≤ 1 with probability at least 1 − e−2s(K+N )(log L) for a fixed
pair of (h, x).
22

Part II: Covering argument. Now we will use a covering argument to extend this result
for all (h, x) and thus prove that f (h, x) ≤ 43 uniformly for all (h, x).
We start with defining K and Ni as 0 -nets of S Ks−1 and S N −1 for h and xi , 1 ≤ i ≤ s,
respectively. The bounds |K| ≤ (1 + 20 )2sK and |Ni | ≤ (1 + 20 )2N follow from the covering
numbers of the sphere (Lemma 5.2 in [38]). Here we let N := N1 × · · · × Ns . By taking the
union bound over K × N , we have that f (h, x) ≤ 1 holds uniformly for all (h, x) ∈ K × N with
probability at least
1 − (1 + 2/0 )2s(K+N ) e−2s(K+N ) log L = 1 − e−2s(K+N )(log L−log(1+2/ε0 )) .
For any (h, x) ∈ S Ks−1 × |S N −1 × ·{z
· · × S N −1}, we can find a point (u, v) ∈ K × N satisfying
s times

kh − uk ≤ ε0 and kxi − vi k ≤ ε0 for all 1 ≤ i ≤ s. Conditioned on (6.4), we know that
kAk2 ≤ s(N log(N L/2) + (γ + log s) log L) ≤ s(N + γ + log s) log L.
Now we aim to evaluate |f (h, x) − f (u, v)|. First we consider |f (u, x) − f (u, v)|. Since
2 (u, v) if kx k = kv k = kuk = 1 for 1 ≤ i ≤ s, we have
= σmax
i
i

2 (u, x)
σmax

|f (u, x) − f (u, v)| =

kA(H(u, x))k2F − kA(H(u, v))k2F

≤ kA(H(u, x − v))k · kA(H(u, x + v))k
v
v
u s
u s
X
u
uX
2t
2
2
≤ kAk
kui k kxi − vi k t
kui k2 kxi + vi k2
i=1

i=1

2

≤ 2kAk ε0 ≤ 2s(N + γ + log s)(log L)ε0
where the first inequality is due to ||z1 |2 − |z2 |2 | ≤ |z1 − z2 ||z1 + z2 | for any z1 , z2 ∈ C.
We proceed to estimate |f (h, x) − f (u, x)| by using (6.13) and (6.12),
|f (h, x) − f (u, x)| ≤ J1 + J2 + J3
p
≤ (2kAk2 + 2 2s(K + N ) log L + 16s(K + N ) log L)ε0
≤ 25s(K + N + γ + log s)(log L)ε0
where (6.13) and (6.12) give
kA(H(h, x))k2F − kA(H(u, x))k2F ≤ kA(H(h − u, x))k kA(H(h + u, x))k ≤ 2kAk2 ε0 ,
p
p
= 2 2s(K + N ) log L · |σmax (h, x) − σmax (u, x)| ≤ 2 2s(K + N ) log Lε0 ,

J1 =
J2

2
2
J3 = 8s(K + N )(log L) · |σmax
(h, x) − σmax
(u, x)| ≤ 16s(K + N )(log L)ε0 .

Therefore, if 0 =

1
81s(N +K+γ+log s) log L ,

there holds

4
f (h, x) ≤ f (u, v) + |f (u, x) − f (u, v)| + |f (h, x) − f (u, x)| ≤
{z
} 3
|
≤27s(K+N +γ+log s)(log L)ε0 ≤ 13

for all (h, x) uniformly with probability at least 1 − e−2s(K+N )(log L−log(1+2/ε0 )) . By letting
L ≥ Cγ s(K + N ) log L with Cγ reasonably large and γ ≥ 1, we have log L − log (1 + 2/ε0 ) ≥
1
1
2 (1 + log(γ)) and with probability at least 1 − γ exp(−s(K + N )).

6.2

Proof of the local restricted isometry property

Lemma 6.8. Conditioned on (6.5) and (6.9), the following RIP type of property holds:
2
3
kX − X0 k2F ≤ kA(X − X0 )k2 ≤ kX − X0 k2F
3
2
uniformly for all (h, x) ∈ Nd ∩ Nµ ∩ N with µ ≥ µh and  ≤
for some numerical constant Cγ .
23

1
15

if L ≥ Cγ µ2 s(K + N ) log2 L

Proof: The main idea of the proof follows two steps: decompose X − X0 onto T and T ⊥ ,
then apply (6.5) and (6.9) to PT (X − X0 ) and PT ⊥ (X − X0 ) respectively.
1
For any X = H(h, x) ∈ N with δi ≤ ε ≤ 15
, we can decompose X − X0 as the sum of two
block diagonal matrices U = blkdiag(Ui , 1 ≤ i ≤ s) and V = blkdiag(Vi , 1 ≤ i ≤ s) where each
pair of (Ui , Vi ) corresponds to the orthogonal decomposition of hi x∗i − hi0 x∗i0 ,

hi x∗i − hi0 x∗i0 := (αi1 αi2 − 1)hi0 x∗i0 + αi2 h̃i x∗i0 + αi1 hi0 x̃∗i + h̃i x̃∗i
|
{z
} | {z }
Ui ∈Ti

(6.15)

Vi ∈Ti⊥

which has been briefly discussed in (6.2) and (6.3). Note that A(X − X0 ) = A(U + V ) and
kA(U )k − kA(V )k ≤ kA(U + V )k ≤ kA(U )k + kA(V )k.
Therefore, it suffices to have a two-side bound for kA(U )k and an upper bound for kA(V )k
where U ∈ T and V ∈ T ⊥ in order to establish the local isometry property.
Estimation of kA(U )k:

For kA(U )k, we know from Lemma 6.3 that
r
r
9
11
kU kF ≤ kA(U )k ≤
kU kF
10
10

(6.16)
δ2

i
and hence we only need to compute kU kF . By Lemma 6.1, there also hold kVi kF ≤ 2(1−δ
di0
i)
and δi − kVi kF ≤ kUi kF ≤ δi + kVi kF , i.e.,




δi2
δi2
di0 ≤ kUi kF ≤ δi +
di0 , 1 ≤ i ≤ s.
δi −
2(1 − δi )
2(1 − δi )




Ps
ε
ε
2 , it is easy to get δd
1
−
kU
k
With kU k2F =
≤
kU
k
≤
δd
1
+
0
i F
0
F
i=1
2(1−ε)
2(1−ε) .
Combined with (6.16), we get
r 
r 


9
11
ε
ε
1−
δd0 ≤ kA(U )k ≤
1+
δd0 .
(6.17)
10
2(1 − ε)
10
2(1 − ε)

Estimation of kA(V )k: Note that V is a block-diagonal matrix with rank-1 block. So
applying Lemma 6.6 gives us
q
4
2
2
2 (h̃, x̃)(K + N ) log L + 8sσ 2 (h̃, x̃)(K + N ) log L
kA(V )k ≤ kV kF + 2 2skV k2F σmax
max
3
(6.18)
 
h̃1
 .. 
2 (h̃, x̃) to
where V = H(h̃, x̃) and h̃ =  .  . It suffices to get an estimation of kV kF and σmax
h̃s
bound kA(V )k in (6.18).
Lemma 6.1 says that kh̃i kkx̃i k ≤
kx̃i k ≤

δi2
2(1−δi ) di0

δi
2δi p
kxi k ≤
di0 ,
1 − δi
1 − δi

≤
√

ε
2(1−ε) δi di0

if ε < 1. Moreover,

p
LkB h̃i k∞ ≤ 6µ di0 ,

if (h, x) belongs to Nd ∩ Nµ ∩ N . For kV kF ,
v
v
u s
u s
uX
uX
2
t
kVi kF = t
kh̃i k2 kx̃i k2 ≤
kV kF =
i=1

i=1

24

1≤i≤s

εδd0
.
2(1 − ε)

(6.19)

2 (h̃, x̃) by using (6.19),
Now we aim to get an upper bound for σmax
2
σmax
(h̃, x̃)

= max

1≤l≤L

s
X

|b∗l h̃i |2 kx̃i k2

≤ C0

i=1

µ2

Ps

2 2
i=1 δi di0

L

= C0

µ2 δ 2 d20
.
L

2 (h̃, x̃) into (6.18)
By substituting the estimations of kV kF and σmax
√ 2 2r
ε2 δ 2 d20
2εδ d0 C0 µ2 s(K + N ) log L 8C0 µ2 δ 2 d20 s(K + N ) log L
2
kA(V )k ≤
+
+
. (6.20)
2
3(1 − ε)
1−ε
L
L

By letting L ≥ Cγ µ2 s(K + N ) log2 L with Cγ sufficiently large and combining (6.20) and (6.17),
we have
r
r
2
3
δd0 ≤ kA(U )k − kA(V )k ≤ kA(U + V )k ≤ kA(U )k + kA(V )k ≤
δd0 ,
3
2
which gives 32 kX − X0 k2F ≤ kA(X − X0 )k2 ≤ 23 kX − X0 k2F .

6.3

Proof of the local regularity condition

We first introduce a few notations: for all (h, x) ∈ Nd ∩ N , consider αi1 , αi2 , h̃i and x̃i defined
in (6.2) and define
∆hi = hi − αi hi0 , ∆xi = xi − α−1
i xi0
where

(
(1 − δ0 )αi1 ,
αi (hi , xi ) =
1
(1−δ0 )αi2 ,

if khi k2 ≥ kxi k2
if khi k2 < kxi k2

with

δ
.
(6.21)
10
The function αi (hi , xi ) is defined for each block of X = H(h, x). The particular form of
αi (h, x) serves primarily for proving the Lemma 6.11, i.e., local regularity condition of G(h, x).
We also define




h1 − α1 h1,0
x1 − α1 x1,0




..
..
Ks
Ns
∆h := 
 ∈ C , ∆x := 
∈C .
.
.
δ0 :=

hs − αs hs0

xs − αs xs0

The following lemma gives bounds of ∆xi and ∆hi .
Lemma 6.9. For all (h, x) ∈ Nd ∩ N with  ≤

1
15 ,

there hold

max{k∆hi k22 , k∆xi k22 } ≤ (7.5δi2 + 2.88δ02 )di0 ,
1
k∆hi k22 k∆xi k22 ≤ (δi2 + δ02 )d2i0 .
26
√
√
Moreover, if we assume (hi , xi ) ∈ Nµ additionally, we have LkB(∆hi )k∞ ≤ 6µ di0 .
Proof: We only consider khi k2 ≥ kxi k2 and αi = (1 − δ0 )α1i , and the other case is exactly
the same due to the symmetry. For both ∆hi and ∆xi , by definition,
∆hi = hi − αi hi0 = δ0 αi1 hi0 + h̃i ,


1
1
∆xi = xi −
xi0 = αi2 −
xi0 + x̃i ,
(1 − δ0 )αi1
(1 − δ0 )αi1

(6.22)
(6.23)

where hi = αi1 hi0 + h̃i and xi = αi2 xi0 + x̃i come from the orthogonal decomposition in (6.2).

25

We start with estimating k∆hi k2 . Note that khi k22 ≤ 4di0 and kαi1 hi0 k22 ≤ khi k22 since
(h, x) ∈ Nd ∩ Nµ . By Lemma 6.1, we have
!
2

δi
2
2
2
2
2
k∆hi k2 = kh̃i k2 + δ0 kαi1 hi0 k2 ≤
+ δ0 khi k22 ≤ (4.6δi2 + 4δ02 )di0 .
(6.24)
1 − δi
Then we calculate k∆xi k: from (6.23), we have
k∆xi k2 = αi2 −

2

1
(1 − δ0 )αi1

di0 + kx̃i k2 ≤ αi2 −

1
(1 − δ0 )αi1

δi
2δi √
1−δi kxi k2 ≤ 1−δi di0 for (h, x)
αi2 − (1−δ10 )αi1 , which satisfies

where Lemma 6.1 gives kx̃i k2 ≤
So it suffices to estimate
αi2 −

1
1
1
δ0
=
αi1 αi2 − 1 −
≤
(1 − δ0 )αi1
|αi1 |
1 − δ0
|αi1 |

2

di0 +

4δi2 di0
,
(1 − δi )2

∈ Nd ∩ N .


|(αi1 αi2 − 1)| +

δ0
1 − δ0


.

Lemma 6.1 implies that |αi1 αi2 − 1| ≤ δi , and (6.2) gives




1
δi2
δi2
1
2
(khi k2 − kh̃i k2 ) ≥
1−
kh
k
≥
1
−
(1 − ε)
|αi1 |2 =
i
di0
di0
(1 − δi )2
(1 − δi )2

(6.25)

(6.26)

δi
where kh̃i k ≤ 1−δ
khi k and khi k2 ≥ khi kkxi k ≥ (1 − ε)di0 if khi k ≥ kxi k. Substituting (6.26)
i
into (6.25) gives

1
αi2 −
(1 − δ0 )αi1

≤

1
√
1−ε



δi2
1−
(1 − δi )2

−1/2 
δi +

δ0
1 − δ0


≤ 1.2(δi + δ0 ).

Then we have
k∆xi k22


≤

4δi2
1.44(δi + δ0 ) +
(1 − δi )2
2



di0 ≤ (7.5δi2 + 2.88δ02 )di0 .

(6.27)

δi2 di0
2(1−δi )

and |αi1 | ≤ 2.

Finally, we try to bound k∆hi k2 k∆xi k2 . Lemma 6.1 gives kh̃i k2 kx̃i k2 ≤
Combining them along with (6.22), (6.23), (6.24) and (6.27), we have

2
1
k∆hi k22 k∆xi k22 ≤ kh̃i k22 kx̃i k22 + δ02 |αi1 |2 khi0 k22 k∆xi k22 + αi2 −
kxi0 k22 k∆hi k22
(1 − δ0 )αi1


δi4
2
2
2
2
2
2
≤
+ 4δ0 (7.5δi + 2.88δ0 ) + 1.44(δi + δ0 ) (4.6δi + 4δ0 ) d2i0
4(1 − δi )2
(δ 2 + δ02 )d2i0
.
≤ i
26

By symmetry, similar results hold for the case khi k2 < kxi k2 and max{k∆hi k, k∆xi k} ≤
(7.5δi2 + 2.88δ02 )di0 .
√
Next,
under
the
additional
assumption
(h,
x)
∈
N
,
we
now
prove
LkB(∆hi )k∞ ≤
µ
√
6µ di0 :
Case 1: khi k2 ≥ kxi k2 and αi = (1 − δ0 )αi1 . By Lemma 6.1 gives |αi1 | ≤ 2, which implies
√
√
√
LkB(∆hi )k∞ ≤ LkBhi k∞ + (1 − δ0 )|αi1 | LkBhi0 k∞
p
p
p
≤ 4µ di0 + 2(1 − δ0 )µh di0 ≤ 6µ di0 .
Case 2: khi k2 < kxi k2 and αi =

1
(1−δ0 )αi2 .

Using the same argument as (6.26) gives


|αi2 | ≥ 1 −
2

δi2
(1 − δi )2
26


(1 − ε).

Therefore,
√
√
LkB(∆hi )k∞ ≤ LkBhi k∞ +

√
1
LkBh0 k∞
(1 − δ0 )|αi2 |
√
−1/2

p
p
δi2
µh d0
√
≤ 6µ d0 .
≤ 4µ d0 + 1 −
2
(1 − δi )
(1 − δ0 ) 1 − ε

Lemma 6.10. (Local Regularity for F (h, x)) Conditioned on (5.3) and (6.9), the following
inequality holds
Re (h∇Fh , ∆hi + h∇Fx , ∆xi) ≥

√
δ 2 d20
− 2 sδd0 kA∗ (e)k,
8

uniformly for any (h, x) ∈ Nd ∩ Nµ ∩ N with  ≤
numerical constant C.
Proof:

1
15

if L ≥ Cµ2 s(K + N ) log2 L for some

First note that for
I0 = h∇Fh , ∆hi + h∇Fx , ∆xi =

s
X
h∇Fhi , ∆hi i + h∇Fxi , ∆xi i.
i=1

For each component, recall that (2.18) and (2.19), we have
h∇Fhi , ∆hi i + h∇Fxi , ∆xi i = hA∗i (A(X − X0 ) − e)xi , ∆hi i + h(A∗i (A(X − X0 ) − e))∗ hi , ∆xi i
= hA(X − X0 ) − e, Ai ((∆hi )x∗i + hi (∆xi )∗ )i .
Define Ui and Vi as
Ui := αi hi0 (∆xi )∗ + αi −1 (∆hi )x∗i0 ∈ Ti ,

Vi := ∆hi (∆xi )∗ .

(6.28)

Here Vi does not necessarily belong to Ti⊥ . From the way of how ∆hi , ∆xi , Ui and Vi are
constructed, two simple relations hold:
hi x∗i − hi0 x∗i0 = Ui + Vi ,
(∆hi )x∗i + hi (∆xi )∗ = Ui + 2Vi .
Define U := blkdiag(U1 , · · · , Us ) and V := blkdiag(V1 , · · · , Vs ). I0 can be simplified to
I0 =

s
s
X
X
h∇Fhi , ∆hi i + h∇Fxi , ∆xi i =
hA(U + V ) − e, Ai (Ui + 2Vi )i
i=1

i=1

= hA(U + V ), A(U + 2V )i − he, A(U + 2V )i .
|
{z
} |
{z
}
I01

I02

Now we will give a lower bound for Re(I01 ) and an upper bound for Re(I02 ) so that the lower
bound of Re(I0 ) is obtained. By the Cauchy-Schwarz inequality, Re(I01 ) has the lower bound
Re(I01 ) ≥ (kA(U )k − kA(V )k)(kA(U )k − 2kA(V )k).

(6.29)

In the following, we will give an upper bound for kA(V )k and a lower bound for kA(U )k.

27

Upper bound for kA(V )k: Note that V is a block-diagonal matrix with rank-1 blocks, and
applying Lemma 6.6 results in
s

p
4X
2
kV k2F +2σmax (∆h, ∆x)kV kF 2s(K + N ) log L+8sσmax
(∆h, ∆x)(K+N ) log L.
3
i=1
√
√
By using Lemma 6.9, we have k∆hi k2 ≤ (7.5δi2 + 2.88δ02 )di0 and LkB(∆hi )k∞ ≤ 6µ di0 .
2 (∆h, ∆x) gives
Substituting them into σmax
!
s
s
X
36µ2 X
272µ2 δ 2 d20
2
∗
2
2
σmax (∆h, ∆x) = max
|bl ∆hi | k∆xi k
≤
(7.5δi2 + 2.88δ02 )d2i0 ≤
.
1≤l≤L
L
L
kA(V )k2 ≤

i=1

i=1

For kV kF , note that k∆hi k2 k∆xi k2 ≤
kV k2F =

s
X

1
2
26 (δi

+ δ02 )d2i0 and thus
s

k∆hi k2 k∆xi k2 ≤

i=1

δ 2 d20
1 X 2
1
· 1.01δ 2 d20 =
.
(δi + δ02 )d2i0 ≤
26
26
25
i=1

1
Then by δ ≤ ε ≤ 15
and letting L ≥ Cµ2 s(K + N ) log2 L for a sufficiently large numerical
constant C, there holds

kA(V )k2 ≤

δ 2 d20
δd0
=⇒ kA(V )k ≤
.
16
4

(6.30)

1
Lower bound for kA(U )k: By the triangle inequality, kU kF ≥ δd0 − 15 δd0 ≥ 45 δd0 if  ≤ 15
since kV kF ≤ 0.2δd0 . Since U ∈ T , by Lemma 6.3, there holds
r
9
3
kA(U )k ≥
kU kF ≥ δd0 .
(6.31)
10
4

With the upper bound of A(V ) in (6.30), the lower bound of A(U ) in (6.31), and (6.29), we
δ 2 d2
get Re(I01 ) ≥ 8 0 .
Now let us give an upper bound for Re(I02 ),
kI02 k ≤ kA∗ (e)kkU + 2V k∗ = kA∗ (e)k

s
X
i=1

≤
≤

√
√

2kA∗ (e)k

s
X

k Ui + 2Vi k∗
| {z }
rank-2

kUi + 2Vi kF

i=1

√
2skA∗ (e)kkU + 2V kF ≤ 2 sδd0 kA∗ (e)k

where k · k and k · k∗ are a pair of dual norms and
kU + 2V kF ≤ kU + V kF + kV kF ≤ δd0 + 0.2δd0 ≤ 1.2δd0 .
Combining the estimation of Re(I01 ) and Re(I02 ) above leads to
Re(h∇Fh , ∆hi + h∇Fx , ∆xi) ≥

√
δ 2 d20
− 2 sδd0 kA∗ (e)k.
8

T
Lemma 6.11. (Local Regularity for G(h, x) For any (h, x) ∈ Nd N with  ≤
9
11
9
11
10 d0 ≤ d ≤ 10 d0 , 10 di0 ≤ di ≤ 10 di0 , the following inequality holds uniformly
Re (h∇Ghi , ∆hi i + h∇Gxi , ∆xi i) ≥ 2δ0

p
δp
ρGi (hi , xi ) =
ρGi (hi , xi ),
5

1
15

and

(6.32)

where ρ ≥ d2 + 2kek2 . Immediately, we have
Re (h∇Gh , ∆hi + h∇Gx , ∆xi) =

s
X

Re (h∇Ghi , ∆hi i + h∇Gxi , ∆xi i) ≥

i=1

28

δp
ρG(h, x). (6.33)
5

Remark 6.12. For the local regularity condition for G(h, x), we use the results from [21] when
s = 1. This is because each component Gi (h, x) only depends on (hi , xi ) by definition and thus
the lower bound of Re (h∇Ghi , ∆hi i + h∇Gxi , ∆xi i) is completely determined by (hi , xi ) and
δ0 , and is independent of s.
Proof:

For each i : 1 ≤ i ≤ s, ∇Ghi (or ∇Gxi ) only depends on hi (or xi ) and there holds
Re (h∇Ghi , ∆hi i + h∇Gxi , ∆xi i) ≥ 2δ0

p
δp
ρGi (hi , xi ) =
ρGi (hi , xi ),
5

which follows exactly from Lemma 5.17 in [21]. For (6.33), by definition of ∇Gh and ∇Gx
in (2.22),
Re (h∇Gh , ∆hi + h∇Gx , ∆xi) =

s
X

Re (h∇Ghi , ∆hi i + h∇Gxi , ∆xi i)

i=1
s

≥

δp
δ Xp
ρGi (hi , xi ) ≥
ρG(h, x)
5
5
i=1

where G(h, x) =

Ps

i=1 Gi (hi , xi ).

Lemma 6.13. (Proof of the Local Regularity Condition) Conditioned on (5.3), for the
objective function Fe(h, x) in (2.17), there exists a positive constant ω such that
i
h
(6.34)
k∇Fe(h, x)k2 ≥ ω Fe(h, x) − c
+

with c = kek2 + 2000skA∗ (e)k2 and ω =
ρ ≥ d2 + 2kek2 .
Proof:

d0
7000

for all (h, x) ∈ Nd ∩ Nµ ∩ N . Here we set

Following from Lemma 6.10 and Lemma 6.11, we have
Re(h∇Fh , ∆hi + h∇Fx , ∆xi) ≥
Re(h∇Gh , ∆hi + h∇Gx , ∆xi) ≥

√
δ 2 d20
− 2 sδd0 kA∗ (e)k
8
δd p
9δd0 p
G(h, x) ≥
G(h, x)
5
50

9
for all (h, x) ∈ Nd ∩ Nµ ∩ N where ρ ≥ d2 + 
2kek2 ≥ d2 and 10
d0 ≤ d ≤ 11
10 d0 . Adding
e
e
them together gives Re h∇Fh , ∆hi + h∇Fx , ∆xi on the left side. Moreover, Cauchy-Schwarz
inequality implies


p
Re h∇Feh , ∆hi + h∇Fex , ∆xi ≤ 4δ d0 k∇Fe(h, x)k

where both k∆hk2 and k∆xk2 are bounded by 8δ 2 d0 in Lemma 6.9 since
k∆hk2 =

s
X
i=1

k∆hi k2 ≤

s
X
(7.5δi2 + 2.88δ02 )di0 ≤ 8δ 2 d0 .
i=1

Therefore,
p
p
√
δ 2 d20 9δd0 G(h, x)
+
− 2 sδd0 kA∗ (e)k ≤ 4δ d0 k∇Fe(h, x)k.
8
50
Dividing both sides of (6.35) by δd0 , we obtain
4
√ k∇Fe(h, x)k ≥
d0
≥

√
δd0
9p
δd0
+
− 2 skA∗ (e)k
G(h, x) +
12
50
24
p
√
δd
1 p
√ [ F0 (h, x) + G(h, x)] + 0 − 2 skA∗ (e)k
24
6 6
29

(6.35)

1
0
√
where the Local RIP condition (5.3) implies F0 (h, x) ≤ 23 δ 2 d20 and hence δd
12 ≥ 6 6
where F0 (h, x) is defined in (2.12).
Note that (5.6) gives
√
√
q √
q
6δd0 4 s ∗
∗
∗
2 [Re(hA (e), X − X0 i)]+ ≤ 2 2skA (e)kδd0 ≤
+ √ kA (e)k.
4
6

p
F0 (h, x),

(6.36)

By (6.36) and Fe(h, x) − kek2 ≤ F0 (h, x) + 2[Re(hA∗ (e), X − X0 i)]+ + G(h, x), there holds
4
√ k∇Fe(h, x)k ≥
d0

q

p
1 h p
√
F0 (h, x) + 2 [Re(hA∗ (e), X − X0 i)]+ + G(h, x)
6 6
!
√
√
√
δd0
1
6δd0 4 s ∗
+
− √
+ √ kA (e)k − 2 skA∗ (e)k
24
4
6 6
6
rh

i
√
1
∗
2
e
√
− 1000skA (e)k .
≥
F (h, x) − kek
+
6 6
p
p
For any nonnegative real numbers a and b, we have [ (x − a)+ − b]+ + b ≥ (x − a)+ and
it implies
p
p
(x − a)+
(x − a)+ − b]2+ + b2 ) =⇒ [ (x − a)+ − b]2+ ≥
− b2 .
2
√
Therefore, by setting a = kek2 and b = 1000skA∗ (e)k, there holds
"
#
e(h, x) − kek2
d
F
0
k∇Fe(h, x)k2 ≥
− 1000skA∗ (e)k2
3500
2
+
i
d0 h e
≥
F (h, x) − (kek2 + 2000skA∗ (e)k2 ) .
7000
+
(x − a)+ ≤ 2([

6.4

Local smoothness

Lemma 6.14. Conditioned on (5.3), (5.4) and (6.4), for any z := (h, x) ∈ C(K+N )s and
w := (u, v) ∈ C(K+N )s such that z and z + w ∈ N ∩ NFe , there holds
k∇Fe(z + w) − ∇Fe(z)k ≤ CL kwk,
with



2ρ
2L
CL ≤ 10kAk d0 +
5+ 2
min di0
µ
p
where ρ ≥ d2 + 2kek2 and kAk ≤ s(N log(N L/2) + (γ + log s) log L) holds with probability at
least 1 − L−γ from Lemma 6.2.
In particular, L = O((µ2 + σ 2 )s(K + N ) log2 L) and kek2 = O(σ 2 d20 ) follows from kek2 ∼
σ 2 d20 2
2L χ2L



2

and (6.14). Therefore, CL can be simplified to
CL = O(d0 sκ(1 + σ 2 )(K + N ) log2 L)

by choosing ρ ≈ d2 + 2kek2 .
Proof: By Lemma 5.6, we know that both z = (h, x) and z+w = (h+u, x+v) ∈ Nd ∩Nµ ∩N .
Note that
∇Fe = (∇Feh , ∇Fex ) = (∇Fh + ∇Gh , ∇Fx + ∇Gx ),
where (2.18), (2.19), (2.20) and (2.21) give ∇Fh , ∇Fx , ∇Gh and ∇Gx . It suffices to find out
the Lipschitz constants for all of those four functions.
30

Step 1: We first estimate the Lipschitz constant for ∇Fh and the result can be applied to
∇Fx due to symmetry.
∇Fh (z + w) − ∇Fh (z) = A∗ A(H(h + u, x + v))(x + v) − [A∗ A(H(h, x))x + A∗ (y)v]
= A∗ (A(H(h + u, x + v) − H(h, x)))(x + v)
+ A∗ A(H(h, x) − H(h0 , x0 ))v − A∗ (e)v
= A∗ (A(H(h + u, v) + H(u, x)))(x + v)
+ A∗ A(H(h, x) − H(h0 , x0 ))v − A∗ (e)v.
pPs

2
2
i=1 khi k kxi k

≤ khkkxk and z, z + w ∈ Nd directly implies
p
kH(u, x) + H(h + u, v)kF ≤ kukkxk + kh + ukkvk ≤ 2 d0 (kuk + kvk)
√
where kh + uk ≤ 2 d0 . Moreover, (5.3) implies
Note that kH(h, x)kF ≤

kH(h, x) − H(h0 , x0 )kF ≤ d0
√
since z ∈ Nd ∩ Nµ ∩ N . Combined with kA∗ (e)k ≤ εd0 in (5.4) and kx + vk ≤ 2 d0 , we have
k∇Fh (z + w) − ∇Fh (z)k ≤ 4d0 kAk2 (kuk + kvk) + εd0 kAk2 kvk + εd0 kvk
≤ 5d0 kAk2 (kuk + kvk).

(6.37)

Due to the symmetry between ∇Fh and ∇Fx , we have,
k∇Fx (z + w) − ∇Fx (z)k ≤ 5d0 kAk2 (kuk + kvk).

(6.38)

In other words,
√
k∇F (z + w) − ∇F (z)k ≤ 5 2d0 kAk2 (kuk + kvk) ≤ 10d0 kAk2 kwk
√
where kuk + kvk ≤ 2kwk.
Step 2: We estimate the upper bound of k∇Gxi (zi +wi )−∇Gxi (zi )k. Implied by Lemma 5.19
in [21], we have
k∇Gxi (zi + wi ) − ∇Gxi (zi )k ≤

5di0 ρ
kvi k.
d2i

(6.39)

We estimate the upper bound of k∇Ghi (z + w) − ∇Ghi (z)k. Denote
 


 
khi + ui k2
khi k2
ρ
0
0
∇Ghi (z + w) − ∇Ghi (z) =
G0
(hi + ui ) − G0
hi
2di
2di
2di
|
{z
}

Step 3:

j1

ρL
+
8di µ2
|

L 
X
l=1

G00



L|b∗l (hi + ui )|2
8di µ2



b∗l (hi

+ ui ) −

G00



L|b∗l hi |2
8di µ2



b∗l hi

{z

5di0 ρ
kui k,
d2i

kj2 k ≤

31

3ρLdi0
kui k.
2d2i µ2

bl .
}

j2

Following the same estimation of j1 and j2 in Lemma 5.19 of [21], we have
kj1 k ≤



(6.40)

Therefore, combining (6.39) and (6.40) gives
v
u s
uX
k∇G(z + w) − ∇G(z)k = t
(k∇Ghi (z + w) − ∇Ghi (z)k2 + k∇Gxi (z + w) − ∇Gxi (z)k2 )
i=1

v
v
 uX
 uX

u s
u s
5d
ρ
i0
t
t
kui k2 + max
kvi k2
≤ max
2
d
i
i=1
i=1




5di0 ρ
5di0 ρ 3ρLdi0
+
kvk
kuk + max
≤ max
d2i
2d2i µ2
d2i


2ρ
2L
≤
5 + 2 kwk.
min di0
µ


5di0 ρ 3ρLdi0
+
d2i
2d2i µ2

In summary, the Lipschitz constant CL of Fe(z) has an upper bound as follows:
k∇Fe(z + w) − ∇Fe(z)k ≤ k∇F (z + w) − ∇F (z)k + k∇G(z + w) − ∇G(z)k



2L
2ρ
2
5+ 2
kwk.
≤ 10kAk d0 +
min di0
µ

6.5

Robustness condition and spectral initialization

In this section, we will prove the robustness condition (5.4) and also Theorem 3.2. To prove (5.4),
it suffices to show the following lemma, which is a more general version of (5.4).
Lemma 6.15. Consider a sequence of Gaussian independent random variable c = (c1 , · · · , cL ) ∈
λ2
CL where cl ∼ CN (0, Li ) with λi ≤ λ. Moreover, we assume Ai in (2.2) is independent of c.
Then there holds
kA∗ (c) = k max kA∗i (c)k ≤ ξ
1≤i≤s

with probability at least 1 − L−γ if L ≥ Cγ+log(s) ( λξ +
Proof:

λ2
) max{K, N } log L/ξ 2 .
ξ2

It suffices to show that max1≤i≤s kA∗i (c)k ≤ ξ. For each fixed i : 1 ≤ i ≤ s,
A∗i (c)

=

L
X

cl bl a∗il

l=1

The key is to apply the matrix Bernstein inequality (6.53)
and we need to estimate kZl kψ1 , and
√
PL
λ KN
∗
the variance of l=1 Zl . For each l, kcl bl ail kψ1 ≤
follows from (6.58). Moreover, the
L
variance of A∗i (c) is bounded by

λ2 max{K,N }
L

E[A∗i (c)(A∗i (c))∗ ] =
E[(A∗i (c))∗ (A∗i (c))] =

L
X
l=1
L
X

since

E(|cl |2 kail k2 )bl b∗l =

L
N X 2 ∗ λ2 N
λl bl bl 
,
L
L

kbl k2 E(|cl |2 ail a∗il ) =

l=1

l=1
L
X

K
L2

l=1

λ2i IN 

λ2 K
.
L

Letting t = γ log L and applying (6.53) leads to
( √
)
r
2
2 max{K, N } log L
C
λ
λ
KN
log
L
γ
kA∗i (c)k ≤ C0 max
,
≤ ξ.
L
L
Therefore, by taking the union bound over 1 ≤ i ≤ s,
kA∗i (c)k ≤ ξ
with probability at least 1 − L−γ if L ≥ Cγ+log(s) ( λξ +
32

λ2
) max{K, N } log2 L.
ξ2

The robustness condition is an immediate result of Lemma 6.15 by setting ξ =
λ = σd0 .
Corollary 6.16. [Robustness Condition] For e ∼ CN (0,
kA∗i (e)k ≤

εd0
√
,
10 2sκ

with probability at least 1 − L−γ if L ≥ Cγ ( s
Lemma 6.17. For e ∼ CN (0,

σ 2 d20
L IL ),

2 κ2 σ 2
ε2

εd
√0
10 2sκ

and

σ 2 d20
L IL )

∀1 ≤ i ≤ s

+

sκσ
ε ) max{K, N } log L.

there holds

kA∗i (y) − hi0 x∗i0 k ≤ ξdi0 ,

∀1 ≤ i ≤ s

(6.41)

with probability at least 1 − L−γ if L ≥ Cγ+log(s) sκ2 (µ2h + σ 2 ) max{K, N } log L/ξ 2 .
Remark 6.18. The success of the initialization algorithm completely relies on the lemma above.
As mentioned in Section 3, E(A∗i (y)) = hi0 x∗i0 and Lemma 6.41 confirms that A∗i (y) is close to
hi0 x∗i0 in operator norm and hence the spectral method is able to give us a reliable initialization.
Proof:

Note that
A∗i (y) = A∗i Ai (hi0 x∗i0 ) + A∗i (wi )

where
wi = y − Ai (hi0 x∗i0 ) =

X

Aj (hj0 x∗j0 ) + e

(6.42)

j6=i

is independent of Ai . The proof consists of two parts: 1. show that kA∗i Ai (hi0 x∗i0 ) − hi0 x∗i0 k ≤
ξdi0
ξdi0
∗
2 ; 2. prove that kAi (wi )k ≤ 2 .
Part I:

Following from the definition of Ai and A∗i in (2.2),
A∗i Ai (hi0 x∗i0 ) − hi0 x∗i0 =

L
X
l=1

bl b∗l hi0 x∗i0 (ail a∗il − IN ) .
|
{z
}
defined as Zl

where B ∗ B = IK . The sub-exponential norm of Zl is bounded by
kZl kψ1 ≤ max

1≤l≤L

kbl k|b∗l hi0 |k(ail a∗il

− IN )xi0 kψ1

√
µ KN di0
≤
L

√
and k(ail a∗il − IN )xi0 kψ1 ≤ N di0 follows from (6.56).
P
We proceed to estimate the variance of L
l=1 Zl by using (6.55) and (6.57):

where kbl k =

L
X

q

K
L,

maxl |b∗l hi0 |2 ≤

E(Zl Zl∗ )

=

E(Zl∗ Zl )

=

X

µ2 di0
L

|b∗l hi0 |2 x∗i0 E(ail a∗il − IN )2 xi0 bl b∗l ≤

l=1
L
X
l=1

µ2 N d2i0
,
L

L
K X ∗
Kd2i0
|bl hi0 |2 E [(ail a∗il − IN )xi0 x∗i0 (ail a∗il − IN )] ≤
.
L
L
l=1

P
max{K,µ2h N }d2i0
Therefore, the variance of L
. By applying matrix Bernstein
l=1 Zl is bounded by
L
inequality (6.53) and taking the union bound over all i, we prove that
kA∗i Ai (hi0 x∗i0 ) − hi0 x∗i0 k ≤

ξdi0
,
2

∀1 ≤ i ≤ s

holds with probability at least 1 − L−γ+1 if L ≥ Cγ+log(s) max{K, µ2h N } log L/ξ 2 .
33

Part II:

For each 1 ≤ l ≤ L, the l-th entry of wi in (6.42), i.e., (wi )l =

is independent of b∗l hi0 x∗i0 ail and obeys CN (0,
σil2 = L E |(wi )l |2 = L

X

2
σil

L

∗
∗
j6=i bl hj0 xj0 ajl + el ,

P

). Here

|b∗l hj0 |2 kxj0 k2 + σ 2 kX0 k2F

j6=i

µ2h

≤

X

2

khj0 k kxj0 k2 + σ 2 kX0 k2F ≤ (µ2h + σ 2 )kX0 k2F .

j6=i

This gives maxi,l σil2 ≤ (µ2h + σ 2 )kX0 k2F . Thanks to the independence between wi and Ai ,
applying Lemma 6.15 results in
ξdi0
kA∗i (wi )k ≤
(6.43)
2


√ 2 2
µh +σ kX0 kF
(µ2h +σ 2 )kX0 k2F
−γ+1
,
max{K, N } log L.
with probability 1 − L
if L ≥ C max
ξdi0
ξ 2 d2
i0

Therefore, combining (6.42) with (6.43), we get
kA∗i (y) − hi0 x∗i0 k ≤ kA∗i Ai (hi0 x∗i0 ) − hi0 x∗i0 k + kA∗i (wi )k ≤ ξdi0
for all 1 ≤ i ≤ s with probability at least 1 − L−γ+1 if
L ≥ Cγ+log(s) (µ2h + σ 2 )sκ2 max{K, N } log L/ξ 2
where kX0 kF /di0 ≤

√

sκ.

Before moving to the proof of Theorem 3.2, we introduce a property about the projection
onto a closed convex set.
√
√
Lemma 6.19 (Theorem 2.8 in [13]). Let Q := {w ∈ CK | LkBwk∞ ≤ 2 dµ} be a closed
nonempty convex set. There holds
Re(hz − PQ (z), w − PQ (z)i) ≤ 0,

∀ w ∈ Q, z ∈ CK

where PQ (z) is the projection of z onto Q.
With this lemma, we can easily see
kz−wk2 = kz−PQ (z)k2 +kPQ (z)−wk2 +2 Re(hz−PQ (z), PQ (z)−wi) ≥ kPQ (z)−wk2 (6.44)
for all z ∈ CK and w ∈ Q. It means that projection onto nonempty closed convex set is
non-expansive. Now we present the proof of Theorem 3.2.
Proof of Theorem 3.2. By choosing L ≥ Cγ+log(s) (µ2h +σ 2 )s2 κ4 max{K, N } log L/ε2 , we have
kA∗i (y) − hi0 x∗i0 k ≤ ξdi0 ,

∀1 ≤ i ≤ s

(6.45)

where ξ = 10√ε2sκ .
By applying the triangle inequality to (6.45), it is easy to see that
(1 − ξ)di0 ≤ di ≤ (1 + ξ)di0 ,
which gives

9
10 di0

≤ di ≤

11
10 di0

|di − di0 | ≤ ξdi0 ≤

di0
εdi0
√
<
,
10
10 2sκ

where di = kA∗i (y)k according to Algorithm 1.

34

(6.46)

√
√
(0)
Part I: Proof of (u(0) , v (0) ) ∈ √13 Nd ∩ √13 Nµ Note that vi = di kx̂i0 k = di where x̂i0 is
the leading right singular vector of A∗i (y). Therefore,
p
p
p
2 p
(0)
kvi k = di kx̂i0 k = di ≤ (1 + ξ)di0 ≤ √
di0 , ∀1 ≤ i ≤ s
3
(0)

√1 Nd .
3
(0)
that ui

which implies {vi }si=1 ∈

(0)
√1 Nd ∩ √1 Nµ by Lemma 6.19. By Algorithm 1, u
is
i
3
3
√
√
√
the minimizer to the function f (z) = 21 kz − di ĥi0 k2 over Qi := {z| LkBzk∞ ≤ 2 di µ}.
√
(0)
(0)
Obviously, by definition, ui is the projection
of di ĥi0 onto Qi . Note that ui ∈ Qi implies
√
p
√
√
(0)
(0)
LkBui k∞ ≤ 2 di µ ≤ 2 (1 + ξ)di0 µ ≤ 4 √d3i0 µ and hence ui ∈ √13 Nµ .

Now we will prove

∈

Moreover, due to (6.44), there holds
p
(0)
k di ĥi0 − wk2 ≥ kui − wk2 ,

∀w ∈ Qi

(6.47)

In particular, let w = 0 ∈ Qi and immediately we have
(0)

kui k2 ≤ di ≤
(0)

(0)

In other words, {(ui , vi )}si=1 ∈

√1 Nd
3

5
(0)

≤

(0)

We will show kui (vi )∗ − hi0 x∗i0 kF ≤ 4ξdi0 for

2ε
√
sκ

(0)

kui (vi )∗ −hi0 x∗i0 kF
di0

√1 Nµ .
3
(0)

Part II: Proof of (u(0) , v (0) ) ∈ N
1 ≤ i ≤ s so that

∩

4
1
(0)
=⇒ ui ∈ √ Nµ .
3
3

2ε
√
.
5 sκ

First note that σj (A∗i (y)) ≤ ξdi0 for all j ≥ 2, which follows from Weyl’s inequality [28] for
singular values where σj (A∗i (y)) denotes the j-th largest singular value of A∗i (y). Hence there
holds
kdi ĥi0 x̂∗i0 − hi0 x∗i0 k ≤ kA∗i (y) − di ĥi0 x̂∗i0 k + kA∗i (y) − hi0 x∗i0 k ≤ 2ξdi0 .
(6.48)
On the other hand, for any i,




hi0 h∗i0
hi0 h∗i0
IK −
ĥi0 =
ĥi0 x̂∗i0 x̂i0 ĥ∗i0
IK −
di0
di0



1
hi0 h∗i0
hi0 x∗i0
∗
∗
∗
=
((Ai (y) − di ĥi0 x̂i0 ) + ĥi0 x̂i0 −
x̂i0 ĥ∗i0
IK −
di0
di0
di0
1
di
kA∗i (y) − hi0 x∗i0 k +
− 1 ≤ 2ξ
=
di0
di0
where (IK −

hi0 h∗i0
∗
di0 )hi0 xi0

= 0 and (A∗i (y) − di ĥi0 x̂∗i0 )x̂i0 ĥ∗i0 = 0. Therefore, we have

ĥi0 −
√

h∗i0 ĥi0
hi0 ≤ 2ξ,
di0

k

p
p
di ĥi0 − ti0 hi0 k ≤ 2 di ξ,

(6.49)

p
√
and |ti0 | ≤ di /di0 < 2. If we substitute w by ti0 hi0 ∈ Qi into (6.47),
p
(0)
k di ĥi0 − ti0 hi0 k ≥ kui − ti0 hi0 k.
(6.50)
√
√
√
where ti0 hi0 ∈ Qi follows from L|ti0 |kBhi0 k∞ ≤ |ti0 | di0 µh ≤ 2di0 µ.

where ti0 =

di h∗i0 ĥi0
di0

(0)

(0)

Now we are ready to estimate kui (vi )∗ − hi0 x∗i0 kF as follows,
(0)

(0)

(0)

(0)

(0)

(0)

kui (vi )∗ − hi0 x∗i0 kF ≤ kui (vi )∗ − ti0 hi0 (vi )∗ kF + kti0 hi0 (vi )∗ − hi0 x∗i0 kF
di
(0)
(0)
≤ kui − ti0 hi0 kkvi k +
hi0 h∗i0 ĥi0 x̂∗i0 − hi0 x∗i0 .
di0
|
{z
}
F
|
{z
}
I1
I2

35

(0)

Here I1 ≤ 2ξdi because kvi k =
For I2 , there holds

√

√
(0)
di and kui −ti0 hi0 k ≤ 2 di ξ follows from (6.49) and (6.50).


hi0 h∗i0 
di ĥi0 x̂∗i0 − hi0 x∗i0
di0

I2 =

√
≤ kdi ĥi0 x̂∗i0 − hi0 x∗i0 kF ≤ 2 2ξdi0 ,
F

which follows from (6.48). Therefore,
√
2εdi0
(0) (0)
kui (vi )∗ − hi0 x∗i0 kF ≤ 2ξdi + 2 2ξdi0 ≤ 5ξdi0 ≤ √ .
5 sκ

Appendix
Descent Lemma
Lemma 6.20 (Lemma 6.1 in [21]). If f (z, z̄) is a continuously differentiable real-valued function
with two complex variables z and z̄, (for simplicity, we just denote f (z, z̄) by f (z) and keep in
the mind that f (z) only assumes real values) for z := (h, x) ∈ N ∩ NFe . Suppose that there
exists a constant CL such that
k∇f (z + t∆z) − ∇f (z)k ≤ CL tk∆zk,

∀0 ≤ t ≤ 1,

for all z ∈ N ∩ NFe and ∆z such that z + t∆z ∈ N ∩ NFe and 0 ≤ t ≤ 1. Then
f (z + ∆z) ≤ f (z) + 2 Re((∆z)T ∇f (z)) + CL k∆zk2
where ∇f (z) :=

∂f (z,z̄)
∂z

is the complex conjugate of ∇f (z) =

∂f (z,z̄)
∂ z̄ .

Concentration inequality
We define the matrix ψ1 -norm via
kZkψ1 := inf {E[exp(kZk/u)] ≤ 2}.
u≥0

Theorem 6.21. [18] Consider a finite sequence of Zl of independent centered random matrices
with dimension M1 × M2 . Assume that R := max1≤l≤L kZl kψ1 and introduce the random matrix
S=

L
X

Zl .

(6.51)

l=1

Compute the variance parameter
σ02 = max

L
n X

E(Zl Zl∗ ) ,

l=1

L
X

)
E(Zl∗ Zl )

,

(6.52)

l=1

then for all t ≥ 0
√
p
kSk ≤ C0 max{σ0 t + log(M1 + M2 ), R log

LR
σ0

!
(t + log(M1 + M2 ))}

with probability at least 1 − e−t where C0 is an absolute constant.

36

(6.53)

Lemma 6.22 (Lemma 10-13 in [1], Lemma 12.4 in [23]). Let u ∈ Cn ∼ CN (0, In ), then
kuk2 ∼ 21 χ22n and
kkuk2 kψ1 = khu, uikψ1 ≤ Cn
(6.54)
and
E(uu∗ − In )2 = nIn .

(6.55)

Let q ∈ Cn be any deterministic vector, then the following properties hold
√
k(uu∗ − I)qkψ1 ≤ C nkqk,

(6.56)

E[(uu∗ − I)qq ∗ (uu∗ − I)] = kqk2 In .

(6.57)

Let v ∼ CN (0, Im ) be a complex Gaussian random vector in Cm , independent of u, then
√
kkuk · kvkkψ1 ≤ C mn.
(6.58)

Acknowledgement
S.Ling would like to thank Felix Krahmer and Dominik Stöger for the discussion about [29], and
also thank Ju Sun for pointing out the connection between convolutional dictionary learning
and this work.

References
[1] A. Ahmed, B. Recht, and J. Romberg. Blind deconvolution using convex programming.
IEEE Transactions on Information Theory, 60(3):1711–1732, 2014.
[2] H. Bristow, A. Eriksson, and S. Lucey. Fast convolutional sparse coding. In Proceedings of
the IEEE Conference on Computer Vision and Pattern Recognition, pages 391–398, 2013.
[3] T. T. Cai, X. Li, Z. Ma, et al. Optimal rates of convergence for noisy sparse phase retrieval
via thresholded wirtinger flow. The Annals of Statistics, 44(5):2221–2251, 2016.
[4] 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.
[5] P. Campisi and K. Egiazarian. Blind Image Deconvolution: Theory and Applications. CRC
press, 2007.
[6] E. Candès, Y. Eldar, T. Strohmer, and V. Voroninski. Phase retrieval via matrix completion. SIAM Journal on Imaging Sciences, 6(1):199–225, 2013.
[7] E. Candès and B. Recht. Exact matrix completion via convex optimization. Foundations
of Computational Mathematics, 9(6):717–772, 2009.
[8] E. J. Candes, X. Li, and M. Soltanolkotabi. Phase retrieval via Wirtinger flow: Theory
and algorithms. IEEE Transactions on Information Theory, 61(4):1985–2007, 2015.
[9] E. J. Candès, T. Strohmer, and V. Voroninski. Phaselift: Exact and stable signal recovery
from magnitude measurements via convex programming. Communications on Pure and
Applied Mathematics, 66(8):1241–1274, 2013.
[10] V. Chandrasekaran, B. Recht, P. Parrilo, and A. Willsky. The convex geometry of linear
inverse problems. Foundations of Computational Mathematics, 12(6):805–849, 2012.
[11] Y. Chen and E. Candes. Solving random quadratic systems of equations is nearly as easy
as solving linear systems. In Advances in Neural Information Processing Systems, pages
739–747, 2015.
37

[12] A. Eftekhari and M. B. Wakin. Greed is super: A fast algorithm for super-resolution. arXiv
preprint arXiv:1511.03385, 2015.
[13] R. Escalante and M. Raydan. Alternating Projection Methods, volume 8. SIAM, 2011.
[14] A. Goldsmith. Wireless Communications. Cambridge University Press, 2005.
[15] P. Jung, F. Krahmer, and D. Stöger. Blind demixing and deconvolution at near-optimal
rate. arXiv preprint arXiv:1704.04178, 2017.
[16] R. Keshavan, A. Montanari, and S. Oh. Matrix completion from noisy entries. In Advances
in Neural Information Processing Systems, pages 952–960, 2009.
[17] R. H. Keshavan, A. Montanari, and S. Oh. Matrix completion from a few entries. IEEE
Transactions on Information Theory, 56(6):2980–2998, 2010.
[18] V. Koltchinskii et al. Von Neumann entropy penalization and low-rank matrix estimation.
The Annals of Statistics, 39(6):2936–2973, 2011.
[19] K. Lee, Y. Li, M. Junge, and Y. Bresler. Blind recovery of sparse signals from subsampled
convolution. IEEE Transactions on Information Theory, 63(2):802–821, 2017.
[20] X. Li and H. H. Fan. Direct blind multiuser detection for CDMA in multipath without
channel estimation. IEEE Transactions on Signal Processing, 49(1):63–73, 2001.
[21] X. Li, S. Ling, T. Strohmer, and K. Wei. Rapid, robust, and reliable blind deconvolution
via nonconvex optimization. arXiv preprint arXiv:1606.04933, 2016.
[22] S. Ling and T. Strohmer. Self-calibration and biconvex compressive sensing. Inverse
Problems, 31(11):115002, 2015.
[23] S. Ling and T. Strohmer. Blind deconvolution meets blind demixing: Algorithms and
performance bounds. IEEE Transactions on Information Theory, 63(7):4497–4520, July
2017.
[24] J. Liu, J. Xin, Y. Qi, F.-G. Zheng, et al. A time domain algorithm for blind separation of convolutive sound mixtures and L1 constrained minimization of cross correlations.
Communications in Mathematical Sciences, 7(1):109–128, 2009.
[25] D. G. Luenberger and Y. Ye. Linear and Nonlinear Programming, volume 228. Springer,
2015.
[26] M. B. McCoy and J. A. Tropp. Achievable performance of convex demixing. Technical
report, Caltech, 2017, Paper dated Feb. 2013. ACM Technical Report 2017-02.
[27] M. Shafi, A. F. Molisch, P. J. Smith, T. Haustein, P. Zhu, P. De Silva, F. Tufvesson,
A. Benjebbour, and G. Wunder. 5G: A tutorial overview of standards, trials, challenges,
deployment, and practice. IEEE Journal on Selected Areas in Communications, 35(6):1201–
1221, 2017.
[28] G. W. Stewart. Perturbation theory for the singular value decomposition. Technical Report
CS-TR 2539, University of Maryland, September 1990.
[29] D. Stöger, P. Jung, and F. Krahmer. Blind deconvolution and compressed sensing. In
Compressed Sensing Theory and its Applications to Radar, Sonar and Remote Sensing
(CoSeRa), 2016 4th International Workshop on, pages 24–27. IEEE, 2016.
[30] T. Strohmer. Four short stories about Toeplitz matrix calculations. Linear Algebra Appl.,
343/344:321–344, 2002. Special issue on structured and infinite systems of linear equations.

38

[31] P. Sudhakar, S. Arberet, and R. Gribonval. Double sparsity: Towards blind estimation
of multiple channels. In Latent Variable Analysis and Signal Separation, pages 571–578.
Springer, 2010.
[32] J. Sun, Q. Qu, and J. Wright. Complete dictionary recovery over the sphere I: Overview and
the geometric picture. IEEE Transactions on Information Theory, 63(2):853–884, 2017.
[33] J. Sun, Q. Qu, and J. Wright. A geometric analysis of phase retrieval. Foundations of
Computational Mathematics, Aug 2017.
[34] R. Sun and Z.-Q. Luo. Guaranteed matrix completion via non-convex factorization. IEEE
Transactions on Information Theory, 62(11):6535–6579, 2016.
[35] S. Tu, R. Boczar, M. Simchowitz, M. Soltanolkotabi, and B. Recht. Low-rank solutions
of linear matrix equations via procrustes flow. In Proceedings of The 33rd International
Conference on Machine Learning, pages 964–973, 2016.
[36] R. Vannithamby and S. Talwar. Towards 5G: Applications, Requirements and Candidate
Technologies. John Wiley & Sons, 2017.
[37] S. Verdu. Multiuser Detection. Cambridge University Press, 1998.
[38] R. Vershynin. Introduction to the non-asymptotic analysis of random matrices. In Y. C.
Eldar and G. Kutyniok, editors, Compressed Sensing: Theory and Applications, chapter 5.
Cambridge University Press, 2012.
[39] X. Wang and H. V. Poor. Blind equalization and multiuser detection in dispersive CDMA
channels. IEEE Transactions on Communications, 46(1):91–103, 1998.
[40] P.-Å. Wedin. Perturbation bounds in connection with singular value decomposition. BIT
Numerical Mathematics, 12(1):99–111, 1972.
[41] K. Wei, J.-F. Cai, T. F. Chan, and S. Leung. Guarantees of riemannian optimization for
low rank matrix recovery. SIAM Journal on Matrix Analysis and Applications, 37(3):1198–
1222, 2016.
[42] J. Wright, A. Ganesh, K. Min, and Y. Ma. Compressive principal component pursuit.
Information and Inference, 2(1):32–68, 2013.
[43] G. Wunder, H. Boche, T. Strohmer, and P. Jung. Sparse signal processing concepts for
efficient 5G system design. IEEE Access, 3:195—208, 2015.

39

