Nonconvex Matrix Factorization from Rank-One Measurements
Yuanxin Li∗

Cong Ma†

Yuxin Chen‡

Yuejie Chi§

arXiv:1802.06286v1 [] 17 Feb 2018

February 2018

Abstract
We consider the problem of recovering low-rank matrices from random rank-one measurements,
which spans numerous applications including covariance sketching, phase retrieval, quantum state tomography, and learning shallow polynomial neural networks, among others. Our approach is to directly
estimate the low-rank factor by minimizing a nonconvex quadratic loss function via vanilla gradient
descent, following a tailored spectral initialization. When the true rank is small, this algorithm is guaranteed to converge to the ground truth (up to global ambiguity) with near-optimal sample complexity
and computational complexity. To the best of our knowledge, this is the first guarantee that achieves
near-optimality in both metrics. In particular, the key enabler of near-optimal computational guarantees
is an implicit regularization phenomenon: without explicit regularization, both spectral initialization and
the gradient descent iterates automatically stay within a region incoherent with the measurement vectors.
This feature allows one to employ much more aggressive step sizes compared with the ones suggested
in prior literature, without the need of sample splitting.

Keywords: matrix factorization, rank-one measurements, gradient descent, nonconvex optimization

1 Introduction
This paper is concerned with estimating a low-rank positive semidefinite matrix M ♮ ∈ Rn×n from a few
rank-one measurements. Specifically, suppose that the matrix of interest can be factorized as
M ♮ = X ♮ X ♮⊤ ∈ Rn×n ,
♮
where X ♮ ∈ Rn×r (r ≪ n) denotes the low-rank factor. We collect m measurements {yi }m
i=1 about M
taking the form
♮
⊤ ♮ 2
i = 1, · · · , m,
yi = a⊤
i M ai = ai X 2 ,

⊤ m
where {ai ∈ Rn }m
i=1 represent the measurement vectors known a priori. One can think of {ai ai }i=1 as a
♮
1
♮
set of linear sensing matrices (so that yi = hai a⊤
i , M i), which are all rank-one . The goal is to recover M ,
♮
or equivalently, the low-rank factor X , from a limited number of rank-one measurements. This problem
spans a variety of important practical applications, with a few examples listed below.
∗

Department of Electrical and Computer Engineering, The Ohio State University, Columbus, OH 43210, USA; Email:
li.3822@osu.edu
†
Department of Operations Research and Financial Engineering, Princeton University, Princeton, NJ 08544, USA; Email:
congm@princeton.edu
‡
Department of Electrical Engineering, Princeton University, Princeton, NJ 08544, USA; Email: yuxin.chen@princeton.edu
§
Department of Electrical and Computer Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, USA; Email:
yuejiechi@cmu.edu
1
Given that yi is a quadratic function with respect to both X ♮ and ai , the measurement scheme is also referred to as quadratic
sampling.

1

• Covariance sketching. Consider a zero-mean data stream {xt }t∈T , whose covariance matrix M ♮ :=
E[xt x⊤
t ] is (approximately) low-rank. To estimate the covariance matrix, one can collect m aggregated quadratic sketches of the form
yi =

1 X ⊤ 2
(ai xt ) ,
|T |
t∈T

2
⊤
♮
which converges to E[(a⊤
i xt ) ] = ai M ai as the number of data instances grows. This quadratic
covariance sketching scheme can be performed under minimal storage requirement and low sketching
cost. See [1] for detailed descriptions.

• Phase retrieval and mixed linear regression. This problem subsumes as a special case the phase retrieval problem [2], which aims to estimate an unknown signal x♮ ∈ Rn from intensity measurements
♮ 2
(which can often be modeled or approximated by quadratic measurements of the form yi = (a⊤
i x ) ).
This problem has found numerous applications in X-ray crystallography, optical imaging, astronomy,
etc. Another related problem in machine learning is mixed linear regression with two components,
where the data one collects are generated from one of two unknown regressors; see [3] for precise
formulation.
• Quantum state tomography. Estimating the density operator of a quantum system can be formulated
as a low-rank positive semidefinite matrix recovery problem using rank-one measurements, when the
density operator is almost pure [4]. A problem of similar mathematical formulation occurs in phase
space tomography [5], where the goal is to reconstruct the correlation function of a wave field.
• Learning shallow polynomial neural networks. Taking {ai , yi }m
i=1 as training data, our problem is
equivalent to learning a one-hidden-layer, fully-connected neural networkP
with a quadratic activation
function [6, 7, 8], where the output of the network is expressed as y = ri=1 σ(a⊤ x♮i ) with X ♮ =
[x♮1 , x♮2 , · · · , x♮r ] ∈ Rn×r and the activation function σ(z) = z 2 .

1.1 Main Contributions
Due to the quadratic nature of the measurements, the natural least-squares empirical risk formulation is
highly nonconvex and in general challenging to solve. To be more specific, consider the following optimization problem:
m

1 X
2 2
yi − a⊤
X
,
(1)
minimizeX∈Rn×r f (X) :=
i
2
4m
i=1

which aims to optimize a degree-4 polynomial in X and is NP hard in general. The problem, however,
may become tractable under certain random designs, and may even be solvable using simple methods like
gradient descent. Our main finding is the following: under i.i.d. Gaussian design (i.e. ai ∼ N (0, In )),
vanilla gradient descent combined with spectral initialization achieves appealing performance guarantees
both statistically and computationally.

• Statistically, we show that gradient descent converges exactly to the true factor X ♮ (modulo unrecoverable global ambiguity), as soon as the number of measurements exceeds the order of O(nr 4 log n).
When r is fixed independent of n, this sample complexity is near-optimal up to some logarithmic
factor.
• Computationally, to achieve ǫ-accuracy, gradient descent requires an iteration complexity of O(r 2 log(1/ǫ))
(up to logarithmic factors), with a per-iteration cost of O(mnr). When r is fixed independent of m
2

and n, the computational complexity scales linearly with mn, which is proportional to the time taken
to read all data.
These findings significantly improve upon existing results that require either resampling (which is not
sample-efficient and is not the algorithm one actually runs in practice [9, 10, 8]), or high iteration complexity
(which results in high computation cost [11]). In particular, our work is most related to [11] that also
studied the effectiveness of gradient descent. The results in [11] require a sample complexity on the order
of nr 6 log2 n, as well as an iteration complexity of O(n4 r 2 log(1/ǫ)) (up to logarithmic factors) to attain
ǫ-accuracy. In comparison, our theory improves the sample complexity to O(nr 4 log n) and, perhaps more
importantly, establishes a much lower iteration complexity of O(r 2 log(1/ǫ)) (up to logarithmic factor). To
the best of our knowledge, this work is the first nonconvex algorithm (without resampling) that achieves
both near-optimal statistical and computational guarantees with respect to n.

1.2 Surprising Effectiveness of Gradient Descent
Recently, gradient descent has been widely employed to address various nonconvex optimization problems
due to its appealing efficiency from both statistical and computational perspectives. Despite the nonconvexity of (1), [11] showed that within a local neighborhood of X ♮ , where X satisfies

1 σr2 X ♮
♮
,
(2)
X −X F ≤
24 kX ♮ kF
f (X) behaves like a strongly convex function, at least along certain descending directions. However, this
region itself is not enough to guarantee computational efficiency, and consequently, the smoothness parameter derived in [11] is as large as n2 (even ignoring additional polynomial factors in r), leading to a step size
as small as O(1/n4 ) and an iteration complexity of O(n4 log(1/ǫ)). These are fairly pessimistic.
In order to improve computational guarantees, it might be tempting to employ appropriately designed
regularization operations — such as truncation [12] and projection [13]. These explicit regularization operations are capable of stabilizing the search direction, and make sure the whole trajectory is in a basin
of attraction with benign curvatures surrounding the ground truth. However, such explicit regularizations
complicate algorithm implementations, as they introduce more tuning parameters.
Our work is inspired by [14], which uncovers the “implicit regularization” phenomenon of vanilla gradient descent for nonconvex estimation problems such as phase retrieval and low-rank matrix completion. In
words, even without extra regularization operations, vanilla gradient descent always follows a path within
some region around the global optimum with nice geometric structure, at least along certain directions.
The current paper demonstrates that a similar phenomenon persists in low-rank matrix factorization from
rank-one measurements.
To describe this phenomenon in a precise manner, we need to specify which region enjoys the desired
geometric properties. To this end, consider a local region around X ♮ where X is “incoherent”2 with all
sensing vectors in the following sense:


σr2 X ♮
1p
⊤
♮
.
(3)
log n ·
max al X − X 2 ≤
1≤l≤m
24
kX ♮ kF

We term the intersection of (2) and (3) the Region of Incoherence and Contraction (RIC). The nice feature
of the RIC is this: within this region, the loss function f (X) enjoys a smoothness parameter that scales
as O(max{r, log n}) (namely, k∇2 f (x)k . max{r, log n}, which is much smaller than O(n2 ) provided
2
⊤
♮
√This is called incoherent because if X is aligned (and hence coherent) with the sensing vectors, al X − X
O( n) times larger than the right-hand side of (3).

3



2

can be

in [11]). As is well known, a region enjoying a smaller smoothness parameter enables more aggressive
progression of gradient descent.
A key question remains as to how to prove that the trajectory of gradient descent never leaves the RIC.
This is, unfortunately, not guaranteed by standard optimization theory, which only ensures contraction of the
Euclidean error. To address this issue, we resort to the leave-one-out trick [14, 15, 16] that produces auxiliary
trajectories of gradient descent that use all but one sample. This allows us to establish the incoherence
condition by leveraging the statistical independence of the leave-one-out trajectory w.r.t. the corresponding
sensing vector that has been left out. Our theory refines the leave-one-out argument and further establishes
linear contraction in terms of the entry-wise prediction error.

1.3 Notations
We use boldface lowercase (resp. uppercase) letters to represent vectors (resp. matrices). We denote by kxk2
the ℓ2 norm of a vector x, and X ⊤ , kXk and kXkF the transpose, the spectral norm and the Frobenius norm
of a matrix X, respectively. The kth largest singular value of a matrix X is denoted by σk (X). Moreover,
the inner product between two matrices X and Y is defined as hX, Y i = Tr Y ⊤ X , where Tr (·) is
the trace. We also use vec(V ) to denote vectorization of a matrix V . The notation f (n) . g(n) or
f (n) = O(g(n)) means that there exists a universal constant c > 0 such that |f (n)| ≤ c|g(n)|. In addition,
we use c and C with different subscripts to represent positive numerical constants, whose values may change
from line to line.

2 Algorithms and Main Results
To begin with, we present the formal problem setup. Suppose we are given a set of m rank-one measurements as follows
♮ 2
i = 1, · · · , m,
(4)
yi = a⊤
i X 2,

where ai ∈ Rn is the ith sensing vector composed of i.i.d. standard Gaussian entries, i.e. ai ∼ N (0, In ),
for i = 1, · · · , m. The underlying ground truth X ♮ ∈ Rn×r is assumed to have full column rank but not
necessarily having orthogonal columns. Define the condition number of M ♮ = X ♮ X ♮⊤ as

σ12 X ♮
.
(5)
κ= 2
σr (X ♮ )
Throughout this paper, we assume the condition number is bounded by some constant independent of n and
r, i.e. κ = O(1). Our goal is to recover X ♮ , up to (unrecoverable) orthonormal transformation, from the
measurements y = {yi }m
i=1 in a statistically and computationally efficient manner.

2.1 Vanilla Gradient Descent
The algorithm studied herein is a combination of vanilla gradient descent and a judiciously designed spectral
initialization. Specifically, consider minimizing the squared loss:
m

1 X
yi − a⊤
f (X) :=
i X
4m
i=1


2 2
,
2

(6)

which is a nonconvex function. We attempt to optimize this function iteratively via gradient descent
Xt+1 = Xt − µt ∇f (Xt ) ,
4

t = 0, 1, · · · ,

(7)

Algorithm 1: Gradient Descent with Spectral Initialization
m
Input: measurements y = {yi }m
i=1 , and sensing vectors {ai }i=1 .
Parameters: step size µt , rank r, and number of iterations T .
1/2
Initialization: set X0 = Z0 Λ0 , where the columns of Z0 ∈ Rn×r contain the normalized
eigenvectors corresponding to the r largest eigenvalues of the matrix
m

Y =

1 X
yi ai a⊤
i ,
2m

(9)

i=1

and Λ0 is an r × r diagonal matrix, with the entries on the diagonal given as

[Λ0 ]i = λi Y − λ, i = 1, · · · , r,
1 Pm
where λ = 2m
i=1 yi and λi (Y ) is the ith largest eigenvalue of Y .
Gradient loop: for t = 0 : 1 : T − 1, do
m

Xt+1

1 X ⊤
ai Xt
= Xt − µt ·
m
i=1

2
2



− yi ai a⊤
i Xt .

(10)

(11)

Output: XT .
where Xt denotes the estimate in the tth iteration, µt is the step size/learning rate, and the gradient ∇f (X)
is given by
m

1 X ⊤
2
(8)
ai X 2 − yi ai a⊤
∇f (X) =
i X.
m
i=1

[11],3

For initialization, similar to
we apply the spectral method, which sets the columns of X0 as the
top-r eigenvectors — properly scaled — of a matrix Y as defined in (9). The rationale is this: the mean of
Y is given by
1
2
E [Y ] =
X ♮ F In + X ♮ X ♮⊤ ,
2
and hence the principal components of Y form a reasonable estimate of X ♮ , provided that there are sufficiently many samples. The full algorithm is described in Algorithm 1.

2.2 Performance Guarantees
Before proceeding to our main results, we specify the metric used to assess the estimation error of the run
⊤
ning iterates. Since X ♮ P X ♮ P
= X ♮ X ♮⊤ for any orthonormal matrix P ∈ Rr×r , X ♮ is recoverable
up to orthonormal transforms. Hence, we define the error of the tth iterate Xt as

dist Xt , X ♮ = Xt Qt − X ♮ F ,
(12)

where Qt is given by

Qt := argminP ∈Or×r Xt P − X ♮

F

(13)

with Or×r denoting the set of all r ×r orthonormal matrices. Accordingly, we have the following theoretical
performance guarantees of Algorithm 1.
3

Compared with [11], when setting the eigenvalues in (10), we use the sample mean λ rather than λr+1 (Y ) to estimate

1
kX ♮ k2F .
2

5

Theorem 1. Suppose that m ≥ cnr 4 κ3 log n with some large enough constant c > 0, and that the step size
c4
−7 ), the iterates satisfy
obeys 0 < µt := µ = (rκ+log n)
2 σ 2 (X ♮ ) . Then with probability at least 1 − O(mn
r

dist Xt , X

Xt Qt − X

♮

for all t ≥ 0. In addition,
max

1≤l≤m

a⊤
l



♮




t σ 2 X ♮ 
r
2
♮
,
≤ c1 1 − 0.5µσr (X )
kX ♮ kF

2



t p
σr2 X ♮
2
♮
log n ·
≤ c2 1 − 0.5µσr (X )
,
kX ♮ kF

(14)

(15)

for all 0 ≤ t ≤ c3 n5 . Here, c1 , · · · , c4 are some universal positive constants.

Remark 1. The precise expression of required sample complexity in Theorem 1 can be written as m ≥
kX ♮ k6

c σ6 (X ♮F) nr log (κn) with some large enough constant c > 0. By adjusting constants, with probability at
r

least 1 − O(mn−7 ), (15) holds for 0 ≤ t ≤ O(nc5 ) in any power c5 ≥ 1.
Theorem 1 has the following implications.
• Near-optimal sample complexity when r is fixed: Theorem 1 suggests that spectrally-initialized vanilla
gradient descent succeeds as soon as m = O(nr 4 log n). When r = O(1), this leads to near-optimal
sample complexity up to logarithmic factor. In fact, once the spectral initialization is finished, a
sample complexity at m = O(nr 3 log n) can guarantee the linear convergence to the global optima.
To the best of our knowledge, this outperforms all performance guarantees in the literature obtained
for any nonconvex method without requiring resampling.

♮ ≤ ǫkXk ,
• Near-optimal computational complexity: In order to achieve ǫ-accuracy,
i.e.
dist
X
,
X
t
F

it suffices to run gradient descent for T = O r 2 poly log(n) log(1/ǫ)
iterations.
This
results
in
a
total

computational complexity of O mnr 3 poly log(n) log(1/ǫ) .

• Implicit regularization: Theorem 1 demonstrates that both the spectral initialization and the gradient
♮
descent updates provably control the entry-wise error max1≤l≤m a⊤
, and the iterl Xt Qt − X
2
ates remain incoherent with respect to all the sensing vectors. In fact, the entry-wise error decreases
linearly as well, which is not characterized in [14].
Theorem 1 is established using a fixed step size. According to our theoretical analysis, the incoherence
condition (15) has a significant impact on the convergence rate. After a few iterations, the incoherence
condition can be bounded independent of log n, which leads to a larger step size and faster convergence.
Specifically, we have the following corollary.
Corollary 1. Under the same setting of Theorem 1, after Ta = c6 max{κ2 r 2 log n, log3 n} iterations, the
step size can be relaxed as 0 < µt := µ = r2 κ2 σc27 (X ♮ ) , with some universal constant c6 , c7 > 0, then the
r
iterates satisfy

t σ 2 X ♮ 

r
2
♮
♮
dist Xt , X ≤ c1 1 − 0.5µσr (X )
,
(16)
kX ♮ kF
for all t ≥ Ta , with probability at least 1 − O(mn−7 ).

6

3 Related Work
Instead of directly estimating X ♮ , the problem of interest can be also solved by estimating M ♮ = X ♮ X ♮⊤
in higher dimension via nuclear norm minimization, which requires O(nr) measurements for exact recovery
[1, 17, 4, 18]. See also [19, 20, 21, 22] for the phase retrieval problem. However, nuclear norm minimization,
often cast as the semidefinite programming, is in general computationally expensive to deal with large-scale
data.
On the other hand, nonconvex approaches have drawn intense attention in the past decade due to their
ability to achieve computational and statistical efficiency all at once. Specifically, for the phase retrieval
problem, Wirtinger Flow (WF) and its variants [2, 12, 23, 14, 24, 25, 26] have been proposed. As a twostage algorithm, it consists of spectral initialization and iterative gradient updates. This strategy has found
enormous success in solving other problems such as low-rank matrix recovery and completion [13, 27],
blind deconvolution [28], and spectral compressed sensing [29]. We follow a similar route but analyze a
more general problem that includes phase retrieval as a special case.
The paper [11] is most close to our work, which studied the local convexity of the same loss function and
developed performance guarantees for gradient descent using a similar, but different spectral initialization
scheme. As discussed earlier, due to the pessimistic estimate of the smoothness parameter, they only allow
a diminishing learning rate (or step size) of O(1/n4 ), leading to a high iteration complexity. We not only
provide stronger computational guarantees, but also improve the sample complexity, compared with [11].
Algorithms with resampling
AltMin-LRROM [9]
gFM [10]
EP-ROM [8]
AP-ROM [8]

Sample complexity
O(nr 4 log2 n log ( 1ǫ ))
O(nr 3 log ( 1ǫ ))
O(nr 2 log4 n log ( 1ǫ ))
O(nr 3 log4 n log ( 1ǫ ))

Computational complexity
O(mnr log ( 1ǫ ))
O(mnr log ( 1ǫ ))
O(mn2 log ( 1ǫ ))
O(mnr log n log ( 1ǫ ))

Algorithms without resampling
Convex [1]
GD [11]
GD (Algorithm 1, Ours)

Sample complexity
O(nr)
O(nr 6 log2 n)
O(nr 4 log n)

Computational complexity
O(mn2 √1ǫ )
O(mn5 r 3 log4 n log ( 1ǫ ))
O(mnr max{log2 n, r 2 } log ( 1ǫ ))

Table 1: Comparisons with existing results in terms of sample complexity and computational complexity to
reach ǫ-accuracy. The top half of the table is concerned with algorithms that require resampling, while the
bottom half of the table covers algorithms without resampling.
Several other existing works have suggested different approaches for low-rank matrix factorization from
rank-one measurements, of which the statistical and computational guarantees to reach ǫ-accuracy are summarized in Table 1. We note our guarantee is the only one that achieves simultaneous near-optimal sample
complexity and computational complexity. Iterative algorithms based on alternating minimization or noisy
power iterations [9, 10, 8] require a fresh set of samples at every iteration, which is never executed in practice,
and the sample complexity grows unbounded for exact recovery.
Many nonconvex methods have been proposed and analyzed recently to solve the phase retrieval problem, including the Kaczmarz method [30, 31, 32] and approximate message passing [33]. In [34], the
Kaczmarz method is generalized to solve the problem studied in this paper, but no theoretical performance
guarantees are provided.
The local geometry studied in our paper is in contrast to [35], which studied the global landscape of
phase retrieval, and showed that there are no spurious local minima as soon as the sample complexity is
above O(n log3 n). It will be interesting to study the landscape property of the generalized model in our
7

paper.
Our model is also related to learning shallow neural networks. [36] studied the performance of gradient
descent with resampling and an initialization provided by the tensor method for various activation functions,
however their analysis did not cover quadratic activations. For quadratic activations, [6] adopts a greedy
learning strategy, and can only guarantee sublinear convergence rate. Moreover, [7] studied the optimization
landscape for an over-parameterized shallow neural network with quadratic activation, where r is larger than
n.

4 Outline of Theoretical Analysis
This section provides the proof sketch of the main results, with the details deferred to the appendix. Our
theoretical analysis is inspired by the work of [14] for phase retrieval and follows the general recipe outlined
in [14], while significant changes and elaborate derivations are needed. We refine the analysis to show that
both the signal reconstruction error and the entry-wise error contract linearly, where the latter is not revealed
by [14]. In below, we first characterize a region of incoherence and contraction that enjoys both strong
convexity and smoothness along certain directions. We then demonstrate — via an induction argument —
that the iterates always stay within this nice region. Finally, the proof is complete by validating the desired
properties of spectral initialization.

4.1 Local Geometry and Error Contraction
We start with characterizing a local region around X ♮ , within which the loss function enjoys desired restricted strong convexity and smoothness properties. This requires exploring the property of the Hessian of
f (X), which is given by
m

i

1 X h ⊤
2
X
⊗ ai a⊤
(17)
ai X 2 − yi Ir + 2X ⊤ ai a⊤
∇2 f (X) =
i
i .
m
i=1

Here, we use ⊗ to denote the Kronecker product and hence ∇2 f (X) ∈ Rnr×nr . Now we are ready to state
the following lemma regarding this local region, which will be referred to as the region of incoherence and
contraction (RIC) throughout this paper. The proof is given in Appendix B.
kX ♮ k4F
nr log (nκ) for some
♮
r (X )
me−1.5n − mn−12 , we have

Lemma 1. Suppose the sample size obeys m ≥ c σ4
c > 0. Then with probability at least 1 − c1 n−12 −

sufficiently large constant

vec (V )⊤ ∇2 f (X)vec (V ) ≥ 1.026σr2 (X ♮ ) kV k2F ,

and
∇2 f (X) ≤ 1.5σr2 (X ♮ ) log n + 6 X ♮

2
F

hold simultaneously for all matrices X and V satisfying the following constraints:

1 σr2 X ♮
♮
X −X F ≤
,
24 kX ♮ kF


σr2 X ♮
1p
⊤
♮
log n ·
,
max al X − X
≤
1≤l≤m
24
kX ♮ kF
2

(18)
(19)

(20a)
(20b)

and V = T1 QT − T2 satisfying


2 X♮
σ
1
r
,
T2 − X ♮ ≤
24 kX ♮ k
where QT := argminP ∈Or×r kT1 P − T2 kF . Here, c1 is some absolute positive constant.
8

(21)

The condition (20) on X formally characterizes the RIC, which enjoys the claimed restricted strong
convexity (see (18)) and smoothness (see (19)). With Lemma 1 in mind, it is easy to see that if Xt lies
within the RIC, the estimation error shrinks in the presence of a properly chosen step size. This is given in
the lemma below whose proof can be found in Appendix D.
kX ♮ k4F
nr log (nκ) for some sufficiently large constant
♮
r (X )
−12
c1 n
− me−1.5n − mn−12 , if Xt falls within the RIC as

Lemma 2. Suppose the sample size obeys m ≥ c σ4

c > 0. Then with probability at least 1 −
described in (20), we have

 

dist Xt+1 , X ♮ ≤ 1 − 0.513µσr2 (X ♮ ) dist Xt , X ♮ ,
provided that the step size obeys 0 < µt ≡ µ ≤

constant.

1.026σr2 (X ♮ )

1.5σr2 (X ♮ ) log n+6kX ♮ k2F

2 . Here, c1 > 0 is some universal

Assuming that the iterates {Xt }, stay within the RIC (see (20)) for the first Tc iterations, according to
Lemma 2, we have, by induction, that

√
Tc +1
 

log n σr2 X ♮
1
♮
2
♮
♮
·
dist XTc +1 , X ≤ 1 − 0.513µσr (X )
dist X0 , X ≤ √ · √
n
kX ♮ kF
24 6
as soon as
)
(
♮ 4
X
F
log n,
(22)
Tc ≥ c max log2 n, 4
σr (X ♮ )
for some large enough constant c. The iterates when t ≥ Tc are easier to deal with; in fact, it is easily seen
that Xt+1 stays in the RIC since

♮
max a⊤
≤ max al 2 Xt+1 Qt+1 − X ♮
l Xt+1 Qt+1 − X
1≤l≤m
1≤l≤m
2

√
√
log n σr2 X ♮
1
(23)
·
≤ 6n · √ · √
kX ♮ kF
n
24 6

σr2 X ♮
1p
=
log n ·
,
24
kX ♮ kF

where (23) follows
from Lemma 10 for all t ≥ Tc . Consequently, contraction of the estimation error

♮
dist Xt , X can be guaranteed by Lemma 1 for all t ≥ Tc with probability at least 1−c1 n−12 −me−1.5n −
mn−12 .

4.2 Introducing Leave-One-Out Sequences
It has now become clear that the key remaining step is to verify that the iterates {Xt } satisfy (20) for the
first Tc iterations, where Tc is on the order of (22). Verifying (20b) is conceptually hard since the iterates
{Xt } are statistically dependent with all the sensing vectors {ai }m
i=1 . To tackle this problem, for each
(l)
1 ≤ l ≤ m, we introduce an auxiliary leave-one-out sequence {Xt }, which discards a single measurement
(l)
from consideration. Specifically, the sequence {Xt } is the gradient iterates operating on the following
leave-one-out function

1 X
2 2
.
(24)
yi − a⊤
f (l) (X) :=
i X 2
4m
i:i6=l

See Algorithm 2 for a formal definition of the leave-one-out sequences. Again, we want to emphasize that
Algorithm 2 is just an auxiliary procedure useful for the theoretical analysis, and it does not need to be
implemented in practice.
9

Algorithm 2: Leave-One-Out Versions
Input: Measurements {yi }i:i6=l , and sensing vectors {ai }i:i6=l .
Parameters: Step size µt , rank r, and number of iterations T .
(l)
(l) (l)1/2
(l)
Initialization: X0 = Z0 Λ0
, where the columns of Z0 ∈ Rn×r contain the normalized
eigenvectors corresponding to the r largest eigenvalues of the matrix
Y (l) =

1 X
yi ai a⊤
i ,
2m

(25)

i:i6=l

(l)

and Λ0 is an r × r diagonal matrix, with the entries on the diagonal given as
h
i

(l)
Λ0
= λi Y (l) − λ(l) , i = 1, · · · , r,
i

(26)


1 P
(l) is the ith largest eigenvalue of Y (l) .
where λ(l) = 2m
i:i6=l yi and λi Y
Gradient loop: For t = 0 : 1 : T − 1, do
(l)

(l)

Xt+1 = Xt − µt ·

1 X  ⊤ (l)
ai Xt
m

2
2

i:i6=l


(l)
− yi ai a⊤
i Xt .

(27)

(l)

Output: XT .

4.3 Establishing Incoherence via Induction
Our proof is inductive in nature with the following induction hypotheses:
Xt Qt − X
(l)

♮
F

(l)

max

Xt Qt − Xt Rt

max



Xt Qt − X ♮
a⊤
l

1≤l≤m

1≤l≤m

F

2

  t σ 2 X ♮ 

r
2
,
≤ C1 1 − 0.5σr X ♮ µ
kX ♮ kF
  t r log n σ 2 X ♮ 

2
≤ C3 1 − 0.5σr X ♮ µ
· r ♮ ,
n
kX kF

  t p

2 X♮
σ
r
≤ C2 1 − 0.5σr2 X ♮ µ
log n ·
,
kX ♮ kF

(l)

(l)

where Rt = argminP ∈Or×r Xt Qt − Xt P
C1 + C3 ≤

1
,
24

C2 +

F

(28a)
(28b)
(28c)

, and the positive constants C1 , C2 and C3 satisfy

√
1
6C3 ≤ ,
24

5.86C1 + 5.86C3 +

√

6C3 ≤ C2 .

(29)

Furthermore, the step size µ is chosen as
µ=

c0 σr2 X ♮



σr2 (X ♮ ) log n + kX ♮ k2F

2

(30)

with appropriate universal constant c0 > 0.
Our goal is to show that if the tth iteration Xt satisfies the induction hypotheses (28), then the (t + 1)th
iteration Xt+1 also satisfies (28). It is straightforward to see that the hypothesis (28a) has already been
established by Lemma 2, and we are left with (28b) and (28c). We first establish (28b) in the following
(l)
lemma, which measures the proximity between Xt and the leave-one-out versions Xt , whose proof is
provided in Appendix E.
10

kX ♮ k4F
nr log (nκ)
♮
r (X )

Lemma 3. Suppose the sample size obeys m ≥ c σ4

for some sufficiently large constant

c > 0. If the induction hypotheses (28) hold for the tth iteration, with probability at least 1 − c1 n−12 −
me−1.5n − mn−12 , we have
  t+1 r log n σ 2 X ♮ 

(l)
(l)
· r ♮ ,
max Xt+1 Qt+1 − Xt+1 Rt+1 ≤ C3 1 − 0.5σr2 X ♮ µ
1≤l≤m
n
kX kF
F
as long as the step size obeys (30). Here, c1 > 0 is some absolute constant.
(l)

In addition, the incoherence property of Xt+1 with respect to the lth sensing vector al is relatively easier
to establish, due to their statistical independence. Combined with the proximity bound from Lemma 3, this
allows us to justify the incoherence property of the original iterates Xt+1 , as summarized in the lemma
below, whose proof is given in Appendix F.
kX ♮ k4F
nr log (nκ)
♮
r (X )

Lemma 4. Suppose the sample size obeys m ≥ c σ4

for some sufficiently large constant

c > 0. If the induction hypotheses (28) hold for the tth iteration, with probability exceeding 1 − c1 n−12 −
me−1.5n − 2mn−12 ,

  t+1 p



σr2 X ♮
♮
2
♮
⊤
log n ·
max al Xt+1 Qt+1 − X
≤ C2 1 − 0.5σr X µ
1≤l≤m
kX ♮ kF
2
holds as long as the step size satisfies (30). Here, c1 > 0 is some universal constant.

4.4 Spectral Initialization
Finally, it remains to verify that the induction hypotheses hold for the initialization, i.e. the base case when
t = 0. This is supplied by the following lemma, whose proof is given in Appendix G.
6

kX ♮ k
Lemma 5. Suppose that the sample size exceeds m ≥ c σ6 X ♮F nr log n for some sufficiently large constant
)
r(
c > 0. Then X0 satisfies (28) with probability at least 1 − c1 n−12 − me−1.5n − 3mn−12 , where c1 is some
absolute positive constant.

5 Conclusions
In this paper, we have shown that low-rank positive semidefinite matrices can be recovered from a nearminimal number of random rank-one measurements, via the vanilla gradient descent algorithm following
spectral initialization. Our results significantly improve upon existing results in several ways, both computationally and statistically. In particular, our algorithm does not require resampling at every iteration (and
hence requires fewer samples). The gradient iteration can provably employ a much more aggressive step
size than what was suggested in prior literature (e.g. [11]), thus resulting in much smaller iteration complexity and hence lower computational cost. All of this is enabled by establishing the implicit regularization
feature of gradient descent for nonconvex statistical estimation, where the iterates remain incoherent with
the sensing vectors throughout the execution of the whole algorithm.
There are several problems that are worth exploring in future investigation. For example, our theory
♮
reveals the typical size of the fitting error of Xt (i.e. yi − ka⊤
i X k2 ) in the presence of noiseless data, which
would serve as a helpful benchmark when separating sparse outliers in the more realistic scenario. Another
direction is to explore whether implicit regularization remains valid for learning shallow neural networks
[36]. Since the current work can be viewed as learning a one-hidden-layer fully-connected network with
a quadratic activation function σ(z) = z 2 , it would be of great interest to study if the techniques utilized
herein can be used to develop strong guarantees when the activation function takes other forms.
11

Acknowledgements
The work of Y. Li and Y. Chi is supported in part by AFOSR under the grant FA9550-15-1-0205, by
ONR under the grant N00014-18-1-2142, and by NSF under the grants CAREER ECCS-1818571 and CCF1704245.

Appendices
A

Technical Lemmas

In this section, we document a few useful lemmas that are used throughout the proof.
Lemma 6. [27, Lemma 5.4] For any matrices X, U ∈ Rn×r , we have
q √
XX ⊤ − U U ⊤ F ≥ 2( 2 − 1)σr (X) dist(X, U ).

Lemma 7 (Covering number for low-rank matrices). [37, Lemma 3.1] Let Sr = {X ∈ Rn1 ×n2 , rank(X) ≤
r, kXkF = 1}. Then there exists an ǫ-net S̄r ⊂ Sr with respect to the Frobenius norm obeying S̄r ≤
(9/ǫ)(n1 +n2 +1)r .
Lemma 8. [38, 2] Suppose x1 , · · · , xm are i.i.d.
 real-valued random variables obeying xi ≤ b for some
deterministic number b > 0, E [xi ] = 0, and E x2i = d2 . Setting σ 2 = m · max{b2 , d2 }, we have
!




 
m
X
t2
t
xi ≥ t ≤ min exp − 2 , 25 1 − Φ
P
,
2σ
σ
i=1

where Φ(·) is the cumulative distribution function of a standard Gaussian variable.
Lemma 9. [39, Theorem 5.39] Suppose the ai ’s are i.i.d. random vectors following ai ∼ N (0, In ), i =
1, · · · , m. Then for every t ≥ 0 and 0 < δ ≤ 1,
m

1 X
ai a⊤
≤δ
In −
i
m
i=1

pn
2
holds with probability at least 1 − 2e−ct , where δ = C m
+ √tm . On this event, for all W ∈ Rn×r , there
exists
m
1 X ⊤
2
ai W 2 − kW k2F ≤ δ kW k2F .
m
i=1

Lemma 10. [2] Suppose the ai ’s are i.i.d. random vectors following ai ∼ N (0, In ), i = 1, · · · , m. Then
with probability at least 1 − me−1.5n , we have
√
max kai k2 ≤ 6n.
1≤i≤m

Lemma 11. Fix W ∈ Rn×r . Suppose the ai ’s are i.i.d. random vectors following ai ∼ N (0, In ), i =
1, · · · , m. Then with probability at least 1 − mrn−13 , we have
p
log n kW kF .
max a⊤
W
≤
5.86
i
2
1≤i≤m

12

2
Pr
2
⊤
Proof. Define W = [w1 , w2 , · · · , wr ], then we can write a⊤
i W 2 =
k=1 ai wk . Recognize that
2

wk
follows the χ2 distribution with 1 degree of freedom. It then follows from [40, Lemma 1] that
a⊤
i kwk k
2

P



wk
a⊤
i
kwk k2

2

√
≥ 1 + 2 t + 2t

!

≤ exp (−t),

for any t > 0. Taking t = 13 log n yields


2
2
⊤
ai wk ≤ 34.3 kwk k2 log n ≥ 1 − n−13 .
P
Finally, taking the union bound, we obtain
2
max a⊤
i W 2
1≤i≤m

≤

r
X
k=1

34.3 kwk k22 log n = 34.3 kW k2F log n

with probability at least 1 − mrn−13 .
Lemma 12. Suppose a ∼ N (0, In ). Then for any fixed matrices X, H ∈ Rn×r , we have
i
h
2
2
2
2
2
E a⊤ H 2 a⊤ X 2 = H F X F + 2 H ⊤ X F ;
h
2 i 
2

2
E a⊤ HX ⊤ a
= Tr H ⊤ X
+ Tr H ⊤ XH ⊤ X + HX ⊤ F .



≤ ck kHk2k
Moreover, for any order k ≥ 1, we have E ka⊤ Hk2k
2
F , where ck > 0 is a numerical constant
that depends only on k.
Proof. Let X = [x1 , x2 , · · · , xr ] and H = [h1 , h2 , · · · , hr ]. Based on the simple facts
i
h
E (x⊤ a)2 aa⊤ = kxk22 In + 2xx⊤ ,
i
h
⊤
⊤
E (a⊤ xi )(a⊤ xj )aa⊤ = xi x⊤
j + xj xi + xi xj In ,

we can derive

E



⊤

a H

2
2

⊤

a X

2
2




r
r X
2 
2 
X
⊤
⊤
a xj
E a hi
=
=

i=1 j=1
r 
r X
X
i=1 j=1

khi k22 kxj k22

+2

= kHk2F kXk2F + 2 H ⊤ X



h⊤
i xj

2
,
F

2 

and
E



a⊤ HX ⊤ a

2 



r 



2 X 

2 
X
a ⊤ xi +
a ⊤ h i a ⊤ xi a ⊤ h j a ⊤ xj 
a ⊤ hi
= E
i=1

=

r 
X
i=1

i6=j



khi k22 kxi k22 + 2 h⊤
i xi
13

2 

+

i

 

 

X h
⊤
⊤
x⊤
x⊤
h⊤
h⊤
i hj
i xj + h i xj
j xj + h i h j
i xi
i6=j

Finally, to bound E

h


2
= Tr H ⊤ X
+ HX ⊤

a⊤ H

2k
2

i

2
F


+ Tr H ⊤ XH ⊤ X .

for an arbitrary H ∈ Rn×r , we write the singular value decomposition

of H as H = U ΣV ⊤ , where U = [u1 , u2 , · · · , ur ] ∈ Rn×r , Σ = diag {σ1 , σ2 , · · · , σr }, and V ∈ Rr×r .
This gives
r
X
2
σi2 (a⊤ ui )2 .
a⊤ H 2 =
i=1


2
Let bi = σi
i for i = 1, · · · , r, which are independent random variables obeying bi ∼ N 0, σi due to
= σi2t (2t − 1)!! ≤ ck σi2t for any i = 1, · · · , r and t = 1, · · · , k, where
the fact U ⊤ U = Ir . Since E b2t
i
ck is some large enough constant depending only on k, we arrive at

!k
!k 
r
r
X
X
σi2
= ck kHk2k
b2i  ≤ ck
E
F
a⊤ u

i=1

i=1

as claimed.

Lemma 13. Fix X ♮ ∈ Rn×r . Suppose the ai ’s are i.i.d. random vectors following ai ∼ N (0, In ), i =
1, · · · , m. For any 0 < δ ≤ 1, suppose m ≥ cδ−2 n log n for some sufficiently large constant c > 0. Then
we have
m
1 X ⊤ ♮ 2
2
♮ ♮⊤
♮ 2
≤ δ X ♮ F,
ai X 2 ai a⊤
i − X F In − 2X X
m
i=1

with probability at least 1 − c1 rn−13 , where c1 > 0 is some absolute constant.

Proof. This proof adapts the results of [2, Lemma 7.4] with refining the probabilities. Let a(1) be the first
element of a vector a ∼ N (0, In ). Based on [41, Theorem 1.9], we have
!
m
1/2
2
1 X
2
(ai (1)) − 1 ≥ δ ≤ e2 · e−(c1 δ m) ;
P
m
i=1
!
m
X
1/4
2
1
P
(ai (1))4 − 3 ≥ δ ≤ e2 · e−(c2 δ m) ;
m
i=1
!
m
1/6
2
1 X
6
P
(ai (1)) − 15 ≥ δ ≤ e2 · e−(c3 δ m) .
m
i=1

So, by setting m ≫ δ−2 n, we have
m

m

m

i=1

i=1

i=1

1 X
1 X
1 X
(ai (1))2 − 1 ≤ δ,
(ai (1))4 − 3 ≤ δ, and
(ai (1))6 − 15 ≤ δ,
m
m
m

(31)

with probability at least 1 − c4 n−13 for some constant c4 > 0. Moreover, following [40, Lemma 1], we
know


√
P (ai (1))2 ≥ 1 + 2 t + 2t ≤ exp (−t),
14

which gives



P (ai (1))2 ≥ 36.5 log m ≤ exp (−14 log m) = m−14 ,

if setting t = 14 log m. Therefore, as long as m ≥ cn, we have
p
max |ai (1)| ≤ 36.5 log m,

(32)

1≤i≤m

with probability at least 1 − c5 n−13 for some constant c5 > 0.
With (31) and (32), the results in [2, Lemma 7.4] imply that for any 0 < δ ≤ 1, as soon as m ≥
cδ−2 n log n for some sufficiently large constant c, with probability at least 1 − c1 n−13 ,
m

1 X ⊤ 2
2
⊤
ai x ai a⊤
≤ δ kxk22
i − kxk2 I − 2xx
m
i=1

holds for any fixed vector x ∈ Rn . Let X ♮ = [x♮1 , x♮2 , · · · , x♮r ]. Instantiating the above bound for the set of
vectors x♮k , k = 1, . . . , r and taking the union bound, we have
m

1 X ⊤ ♮
ai X
m
i=1

2
a a⊤
2 i i

− X♮

2
I
F

r
X

− 2X ♮ X ♮⊤ ≤

m

1 X ⊤ ♮ 2
♮
a i xk a i a ⊤
i − xk
m
i=1

k=1

≤δ

r
X

x♮k

2
2

= δ X♮

2
I
2

− 2x♮k x♮⊤
k

2
.
F

k=1

B Proof of Lemma 1
The crucial ingredient for proving the lower bound (18) is the following lemma, whose proof is provided in
Appendix C.
X♮

4

Lemma 14. Suppose m ≥ c σ4

F
nr log (nκ) with some large enough positive constant c, then with
(X ♮ )
−12
−1.5n
probability at least 1 − c1 n
− me
, we have


vec (V )⊤ ∇2 f (X)vec (V ) ≥ 2Tr X ♮⊤ V X ♮⊤ V + 1.204σr2 (X ♮ ) kV k2F ,
(33)
r

for all matrices X and V where X satisfies X − X ♮
constant.

F

≤

♮
2
1 σr (X )
24 kX ♮ k .
F

Here, c1 > 0 is some universal

With Lemma 14 in place, we are ready to prove (18). Let V = T1 QT − T2 satisfy the assumptions in
Lemma 1, then we can demonstrate that


Tr X ♮⊤ V X ♮⊤ V

⊤ 
⊤
= Tr X ♮ − T2 + T2 V X ♮ − T2 + T2 V





⊤
⊤ 
⊤
= Tr X ♮ − T2 V X ♮ − T2 V + 2Tr X ♮ − T2 V T2⊤ V + Tr T2⊤ V T2⊤ V


2
≥ Tr T2⊤ V T2⊤ V − X ♮ − T2 kV k2F − 2 X ♮ − T2 kT2 k kV k2F
15

2

2
F

− X ♮ − T2 kV k2F − 2 X ♮ − T2 kT2 k kV k2F
!
 !2


σr2 X ♮
σr2 X ♮
σr2 X ♮
1
1
1
+ X ♮  kV k2F
≥ −
+2·
·
24 kX ♮ k
24 kX ♮ k
24 X ♮

= T2⊤ V


(34)
(35)

≥ −0.0886σr2 (X ♮ ) kV k2F ,

(36)

where (34) follows from the fact that T2⊤ V ∈ Rr×r is a symmetric matrix [42, Theorem 2], (35) arises
2
from the fact T2⊤ V F ≥ 0 as well as the assumptions of Lemma 1, and (36) is based on the fact X ♮ ≥
σr (X ♮ ). Combining (36) with Lemma 14, we establish the lower bound (18).
To prove the upper bound (19) asserted in the lemma, we make the observation that the Hessian in (17)
satisfies
∇2 f (X)

m


i 
1 X h ⊤
⊤
⊤
⊤
♮ 2
I
+
2X
a
a
X
⊗
a
a
kai Xk22 − ka⊤
X
k
=
r
i
i
i
i
i
2
m
i=1
 
m 
⊤



1 X
2
♮
⊤
♮
⊤
≤
ai Ir + 2 ai X 2 Ir ⊗ ai a⊤
X −X
ai X + X
i
m

1
m

≤

i=1
m h
X



⊤
⊤ ♮
♮
ka⊤
Xk
+
ka
X
k
2
2 · ai X − X
i
i

i=1

2

m


1 X h ⊤
⊤
♮
♮
kai Xk2 + ka⊤
=
i X k2 · ai X − X
m

2

i=1

m

+
3
m

≤
|

1 X
♮
2 a⊤
i X
m



2
2

i=1
m
X



−
2
ai a⊤
X♮
i

2
I
F n

{z

+2



2



ai a⊤
i

:=B1

m

+2
|

1 X ⊤ ♮
ai X
m
i=1


− X♮
ai a⊤
i

2
2

2
I
F n

{z

:=B2

2
2

a⊤
i X

i

ai a⊤
i

2
2

− 2X ♮ X ♮⊤ + 2
|
}



}

X♮

(37)

♮
− a⊤
i X



+ 2X ♮ X ♮⊤ + 2 X ♮



⊤
⊤ ♮
♮
ka⊤
Xk
+
ka
X
k
2
2 · ai X − X
i
i

i=1

+ 2 a⊤
i X

2
I
F n

2
2


i 
ai a⊤
i

+ 2X ♮ X ♮⊤



2
I
F n

+ 2X ♮ X ♮⊤ ,
{z
}

:=B3

where (37) follows from the fact kI ⊗ Ak = kAk. It is seen from Lemma 13 that

2
B2 ≤ δ X ♮ F ≤ 0.02σr2 X ♮ ,

σr2 X ♮
when setting δ ≤ 0.02
2 . Moreover, it is straightforward to check that
X♮

F

B3 ≤ 6 X ♮

2
.
F

With regards to the first term B1 , note that by Lemma 11 and (20b), we can bound
a⊤
i X

2

≤

♮
a⊤
i X

2

+

a⊤
i (X

♮

−X )

2

p
≤ 5.86 log n X ♮
16


σr2 X ♮
1p
log n ·
+
F
24
kX ♮ kF

(38)

for 1 ≤ i ≤ m, and therefore,
B1 ≤

1.471σr2

X

♮



m


1 X
log n
ai a⊤
≤ 1.48σr2 X ♮ log n,
i
m

(39)

i=1

where the last line follows from Lemma 9. The proof is then finished by combining (38) with the preceding
bounds on B1 , B2 and B3 .

C

Proof of Lemma 14

Without loss of generality, we assume kV kF = 1. Write
vec (V )⊤ ∇2 f (X)vec (V )

 

m

2
1 X
⊤
⊤
⊤
⊤
⊤
vec (V )
ai X − yi Ir + 2X ai ai X ⊗ ai ai
=
vec (V )
m
2
i=1

m 
m




2
1 X
1 X
⊤
⊤
⊤
⊤
⊤
vec (V )⊤ vec 2ai a⊤
V
X
a
a
X
ai X − yi vec (V ) vec ai ai V +
=
i
i
i
m
m
2
i=1
i=1

m 
2 

2
2
2
1 X
⊤
⊤
⊤ ♮
⊤
⊤
=
ai V + 2 ai XV ai
ai X − ai X
.
(40)
m
2
2
2
i=1

In what follows, we let X = X ♮ + t
X − X♮

F

≤

1 σr (X )
24 kX ♮ k ,
♮

2

σr2 (X ♮ )

kX ♮ kF

H with t ≤ 1/24 and kHkF = 1 which immediately obeys

and express the right-hand side of (40) as

F

p (V , H, t)
m 
1 X
:=
a⊤
i X
m
i=1
|

2

a⊤
i V

2

2
2

+2

{z



⊤
a⊤
i XV ai

2 

m

−

}

:=q(V ,H,t)

1 X ⊤ ♮
ai X
m
i=1

2
2

a⊤
i V

2
2

.

(41)

The aim is thus to control p (V , H, t) for all matrices satisfying kHkF = 1 and kV kF = 1, and for all t
obeying t ≤ 1/24.
We first bound the second term in (41). Let V = [v1 , v2 , · · · , vr ], then by Lemma 13,
m

1 X ⊤ ♮
ai X
m

2
2

a⊤
i V

2
2

i=1

m

=

1 X ⊤ ♮
ai X
m

2
2

i=1

≤
=
≤

r
X
k=1

r
X

k=1
r
X
k=1

1
m

m
X
i=1

vk⊤
kvk k22

− X♮

2
F

V

r 
2
X
a⊤
v
− X♮
i k

1 X ⊤ ♮
ai X
m
i=1
m
X

1
m

i=1

− 2 X ♮⊤ V
2
F

2
a a⊤
2 i i

♮
a⊤
i X

−

2
a a⊤
2 i i

2
X♮ F

2
F

kvk k22 − 2

2
kvk k22
F

− X♮

17

r
X
k=1

k=1


2
⊤
⊤ ♮ 2
ai X 2 ai vk − X ♮
m

2
F

!

− 2X ♮ X ♮⊤

X ♮⊤ vk

k=1

− 2 X ♮⊤ vk

− 2X ♮ X ♮⊤
2
F

r
X

vk

2
2

2
2

≤ δ X♮
By setting δ ≤

2
F

♮
2
1 σr (X )
24 kX ♮ k2 ,
F

m

r
X
k=1

kvk k22 = δ X ♮

2
kV
F

k2F .

we see that with probability at least 1 − c1 rn−13 ,

1 X ⊤ ♮
ai X
m

2
2

a⊤
i V

i=1

2
2

≤ X♮

2
F

V

2
F

+ 2 X ♮⊤ V

2
F

+


1 2
σr X ♮ kV k2F ,
24

(42)

4

kX ♮ kF
n log n.
σr4 (X ♮ )
Next, we turn to the first term q (V , H, t) in (41), and we need to accommodate all matrices satisfying
kHkF = 1 and kV kF = 1, and all scalers obeying t ≤ 1/24. The strategy is that we first establish the
bound of q (V , H, t) for any fixed H, V and t, and then extend the result to a uniform bound for all H, V
and t by covering arguments.
holds simultaneously for all matrices V , as long as m &

C.1 Bound with Fixed Matrices and Scaler
Recall that
m

q (V , H, t) =

1 Xh ⊤
ai X
m
|

2
2

a⊤
i V

i=1

2 i
⊤
+ 2 a⊤
XV
a
.
i
i
{z
}

2
2

:=Gi

We will start by assuming that X and V are both fixed and statistically independent of {ai }m
i=1 . In view of
Lemma 12,
i
h
h
2 i
2
2
⊤
⊤
⊤
+
2E
a
V
a
XV
a
E [Gi ] = E a⊤
X
i
i
i
i
2
2


2
2
2
+ 2 XV ⊤ F + 2Tr X ⊤ V X ⊤ V
= kXk2F kV k2F + 2 X ⊤ V F + 2 Tr X ⊤ V
≤ kXk2F kV k2F + 2 kXk2 kV k2F + 2 kXk2F kV k2F + 2 kXk2 kV k2F + 2 kXk2 kV k2F

2
σr2 X ♮
2
2
♮
H
≤ 9 kXkF kV kF = 9 X + t
kX ♮ kF
F
!

4
♮
σ X
2
2
2
≤ 18
X ♮ F + t2 r
≤ 18.002 X ♮ F ,
2 kHkF
♮
kX kF

where (43) follows kV kF = 1 and X = X ♮ + t

σr2 (X ♮ )

(43)

(44)

H, and (44) arises from the calculations with
kX ♮ kF
kHkF = 1 and t ≤ 1/24. Therefore, if we define Ti = E [Gi ] − Gi , we have E [Ti ] = 0 and
Ti ≤ E [Gi ] ≤ 18.002 X ♮

due to Gi ≥ 0. In addition,
 
 
 
E Ti2 = E G2i − (E [Gi ])2 ≤ E G2i


2  2
2
2
⊤
⊤
⊤
+
2
a
XV
a
a
V
a⊤
X
=E
i
i
i
i
2
2
18

2
,
F

=E

h

a⊤
i X

4
2

h

a⊤
i V
4

4
2

i
4

+ 4E
i

⊤
≤ 9E a⊤
i X 2 ai V 2
r h
i h
8
≤ 9 E a⊤
a⊤
i X 2 E
i V

4

h
8
2

⊤
a⊤
i XV ai

4 i

i

+ 4E

h

⊤
a⊤
i XV ai

2

a⊤
i X

2
2

a⊤
i V

2
2

i

(45)
(46)

≤ 9c4 kXk4F kV k4F = 9c4 X F

4
σr2 X ♮
♮
= 9c4 X + t
H . X♮
kX ♮ kF

(47)
4
,
F

F

where (45) follows from the Cauchy-Schwarz inequality, (46) comes from the Hölder’s inequality, and (47)
is a consequence of Lemma 12. Apply Lemma 8 to arrive at
!
!
m
mσr4 X ♮
1 2  ♮
1 X
≤ exp −c
Ti ≥ σr X
P
,
(48)
4
m
24
kX ♮ kF
i=1

which further leads to
m

m

1 X
1 X
Gi = E [Gi ] −
Ti
m
m
i=1
i=1

1
≥ E [Gi ] − σr2 X ♮
24
 



2
1
2
2
= kXk2F kV k2F + 2 X ⊤ V F + 2 Tr X ⊤ V
+ 2 XV ⊤ F + 2Tr X ⊤ V X ⊤ V − σr2 X ♮
24


1 2  ♮
2
2
2
⊤
⊤
⊤ 2
⊤
≥ kXkF kV kF + 2 X V F + 2 XV F + 2Tr X V X V − σr X .
(49)
24

q (V , H, t) =

Substituting X = X ♮ + t

σr2 (X ♮ )

H for X, and using the facts kHkF = 1, kV kF = 1 and t ≤ 1/24,
kX ♮ kF
we can calculate the following bounds:

 

4
♮
σr2 X ♮
2
2
2 σr X
♮ 2
♮⊤
kXkF = X F + t
Tr
X
H
kHk
+
2t
F
2
kX ♮ kF
kX ♮ kF

 
σ2 X ♮
1
2
2
≥ X ♮ F − 2t r ♮
X ♮ F kHkF ≥ X ♮ F − σr2 X ♮ ;
kX kF
12

 

4
♮
2
♮
2
2
σ X
σr X
2
⊤
♮⊤
⊤
+ t2 r
X ⊤ V F = X ♮⊤ V
Tr
V
HX
V
+
2t
H
V
2
kX ♮ kF
F
F
kX ♮ kF

 
2
2
σ2 X ♮
1
− 2t r ♮
≥ X ♮⊤ V
− σr2 X ♮ ;
X ♮ kHk kV k2F ≥ X ♮⊤ V
kX kF
12
F
F

 
4
♮
2
♮

2
2
σ X
σr X
2
⊤
⊤ ♮ ⊤
XV ⊤ F = X ♮ V ⊤ + t2 r
+
2t
HV
Tr
V
H
X
V
2
kX ♮ kF
F
F
kX ♮ kF

2
2
σr2 X ♮
1 2  ♮
2
♮ ⊤
♮
♮ ⊤
− 2t
≥ X V
σ X ;
−
X
kHk
kV
k
≥
X
V
F
kX ♮ kF
12 r
F
F
 
 






4
♮
σr2 X ♮
⊤
♮⊤
2 σr X
⊤
⊤
⊤
⊤
♮⊤
♮⊤
Tr
H
V
X
V
+
t
Tr
H
V
H
V
Tr X V X V = Tr X V X V + 2t
2
kX ♮ kF
kX ♮ kF
19





4
♮
σr2 X ♮
2
2
2
♮
2 σr X
♮⊤
♮⊤
X kHk kV kF − t
≥ Tr X V X V − 2t
2 kHk kV kF
kX ♮ kF
kX ♮ kF
  

 1
1
♮⊤
♮⊤
+
σr2 X ♮ ,
≥ Tr X V X V −
12 242
which, combining with (49), yields
q (V , H, t)
  

  15
1
♮⊤
♮⊤
+ 2Tr X V X V −
+2 X V
+2 X V
≥
+
σr2 X ♮
24 12 · 24
F
F
  
   15


2
1
♮
♮⊤
♮⊤
2
♮⊤
♮ 2
+
σr2 X ♮
+ 2Tr X V X V + 2σr X −
≥ X F+2 X V
24 12 · 24
F
 


2
2
+ 2Tr X ♮⊤ V X ♮⊤ V + 1.371σr2 X ♮ .
≥ X ♮ F + 2 X ♮⊤ V
2
X♮ F

♮⊤

2

♮

⊤

2

F

C.2 Covering Arguments
Since we have obtained a lower bound on q (V , H, t) for fixed V , H and t, we now move on to extending
it to a uniform bound that covers all V , H and t simultaneously. Towards this, we will invoke the ǫ-net
√
covering arguments for all V , H and t, respectively, and will rely on the fact max1≤i≤m kai k2 ≤ 6n
asserted in Lemma 10. For notational convenience, we define
g (V , H, t) = q (V , H, t)
− X♮

2
F

− 2 X ♮⊤ V

2
F

 


− 2Tr X ♮⊤ V X ♮⊤ V − 1.371σr2 X ♮ .

First, consider the ǫ-net covering argument for V . Suppose V1 and V2 are such that kV1 kF = 1,
kV2 kF = 1, and kV1 − V2 kF ≤ ǫ. Then, since
X ♮⊤ V1

2
F

− X ♮⊤ V2

2
F

≤



X ♮⊤ V1

F

+ X ♮⊤ V2

F



X ♮⊤ (V1 − V2 )

F

≤ 2 X♮

2

ǫ,

and




Tr X ♮⊤ V1 X ♮⊤ V1 − Tr X ♮⊤ V2 X ♮⊤ V2








≤ Tr X ♮⊤ V1 X ♮⊤ V1 − Tr X ♮⊤ V1 X ♮⊤ V2 + Tr X ♮⊤ V1 X ♮⊤ V2 − Tr X ♮⊤ V2 X ♮⊤ V2
≤ X♮

2

kV1 kF kV1 − V2 kF + X ♮

2

kV2 kF kV1 − V2 kF ≤ 2 X ♮

2

ǫ,

we have
|g (V1 , H, t) − g (V2 , H, t)|

2

2

≤ |q (V1 , H, t) − q (V2 , H, t)| + 2 X ♮⊤ V1 − X ♮⊤ V2
F
F




♮⊤
♮⊤
♮⊤
♮⊤
+ 2 Tr X V1 X V1 − Tr X V2 X V2
m 
m 
2 

2
2
1 X
1 X
⊤
⊤
⊤
⊤
ai V1 + 2 ai XV1 ai
−
ai X
a⊤
≤
i X
m
m
2
2
i=1

+ 8 X♮

i=1

2

ǫ

20

2
2

a⊤
i V2

2
2

+2



⊤
a⊤
i XV2 ai

2 

m

≤
≤

1 X
m
1
m
+

a⊤
i X

i=1
m
X

2
m

i=1

a⊤
i V1

2
2

a⊤
i X

i=1
m
X

2

2

·



2
2

a⊤
i V1

− a⊤
i X

2

2

a⊤
i V2

2

+ a⊤
i V2

2



2
2

+

m
2 
2
2 X  ⊤
⊤
ai XV1⊤ ai − a⊤
XV
a
+ 8 X♮
i
i
2
m
i=1

· a⊤
i (V1 − V2 )

2

⊤
⊤
♮
⊤
a⊤
i X (V1 + V2 ) ai · ai X (V1 − V2 ) ai + 8 X

√
√
≤ 6n · kXk2 · 2 6n · 6n · ǫ + 2 · 12n · kXk · 6n · kXk ǫ + 8 X ♮

2
σr2 X ♮
2
2
♮
H + 8 X♮ ǫ
= 216ǫn X + t
♮
kX kF
!

4 X♮
σ
2
2
kHk2 + 8 X ♮ ǫ
X ♮ + t2 r
≤ 432ǫn2
2
kX ♮ kF
 

1
2
≤ 432.75n2 + 8 ǫ X ♮ ≤ σr2 X ♮ ,
24
as long as ǫ =

σr2 (X ♮ )

2

10584n2 kX ♮ k

2

ǫ
2

ǫ

. Based on Lemma 7, the cardinality of this ǫ-net will be

 (n+r+1)r
9
=
ǫ

9 · 10584n2 X ♮
σr2 (X ♮ )

2

!(n+r+1)r

≤ exp (cnr log (nκ)).

Secondly, consider the ǫ-net covering argument for H. Suppose H1 and H2 obey kH1 kF = 1,
kH2 kF = 1, and kH1 − H2 kF ≤ ǫ. Then one has
|g (V , H1 , t) − g (V , H2 , t)|

= |q (V , H1 , t) − q (V , H2 , t)|

!2 
! 2
!


m
2 X♮
2 X♮
2
σ
σ
1 X ⊤
X ♮ + t r ♮ H1 V ⊤ ai 
+ 2 a⊤
=
a⊤
ai X ♮ + t r ♮ H1
i
i V
m
kX kF
kX kF
2
i=1
2

!2 
! 2
!


m
2 X♮
2 X♮
2
σ
σ
1 X ⊤
X ♮ + t r ♮ H2 V ⊤ ai 
+ 2 a⊤
−
a⊤
ai X ♮ + t r ♮ H2
i
i V
m
kX kF
kX kF
2
i=1

m

1 X ⊤
≤
ai V
m
i=1

2

2
2

·

a⊤
i


♮

σ2 X
X ♮ + t r ♮ H1
kX kF
!

!

!2

2

2

−

a⊤
i

!

2 X♮
σ
X ♮ + t r ♮ H2
kX kF



σr2 X ♮
σr2 X ♮
♮
⊤
⊤
♮
⊤
H1 V ai − ai X + t
H2
ai X + t
kX ♮ kF
kX ♮ kF
i=1


√
√
σr2 X ♮
σr2 X ♮
25
25
♮
ǫ · 2 6n ·
ǫ · 12n ·
X + 2 · 6n · t
X♮
≤ 6n · 6n · t
♮
♮
kX kF
24
kX kF
24

75 2 σr2 X ♮
1 2  ♮
♮
≤ ǫn
σ X ,
X
≤
8
kX ♮ kF
24 r
m

2 X
+
m

21

!

2

2

V ⊤ ai

!2

2

ǫ

as long as ǫ =

1
225n2

·

kX ♮ kF
. Based on Lemma 7, the cardinality of this ǫ-net will be
kX ♮ k

 (n+r+1)r
9
=
ǫ

X♮
9 · 225n2 ·
kX ♮ kF

!(n+r+1)r

≤ exp (cnr log n).

Finally, consider the ǫ-net covering argument for all t, such that t ≤ 1/24. Suppose t1 and t2 satisfy
t1 ≤ 1/24, t2 ≤ 1/24 and |t1 − t2 | ≤ ǫ. Then we get
|g (V , H, t1 ) − g (V , H, t2 )|

= |q (V , H, t1 ) − q (V , H, t2 )|

!2 
 !
 ! 2
m
2
♮
2
♮
X
2
σ X
σ X
1
 a⊤
+ 2 a⊤
a⊤
X ♮ + t1 r ♮ H V ⊤ a i 
=
X ♮ + t1 r ♮ H
i
i V
i
m
kX kF
kX kF
2
i=1
2

!2 
 !
 ! 2
m
2 X♮
2 X♮
X
2
σ
σ
1
 a⊤
−
X ♮ + t2 r ♮ H V ⊤ a i 
+ 2 a⊤
X ♮ + t2 r ♮ H
a⊤
i
i V
i
m
kX kF
kX kF
2
i=1

m

2

1 X ⊤
ai V
≤
m
i=1

2
2

·

 !
σr2 X ♮
X + t1
H
kX ♮ kF
!


♮

m

 !
σr2 X ♮
X + t2
H
kX ♮ kF

2

♮

a⊤
i

2

−

!2


♮

σr2 X
σr2 X
♮
⊤
⊤
♮
X
+
t
−
a
H
V
a
H
X
+
t
a⊤
2
i
1
i
i
kX ♮ kF
kX ♮ kF
i=1


√
√
σr2 X ♮
σr2 X ♮
25
25
♮
ǫ · 2 6n ·
ǫ · 12n ·
X + 2 · 6n ·
X♮
≤ 6n · 6n ·
kX ♮ kF
24
kX ♮ kF
24

2 X♮
σ
1 2  ♮
r
♮
≤ 225ǫn2
σ X ,
X
≤
kX ♮ kF
24 r
2 X
+
m

as long as ǫ =

1
5400n2

·

kX ♮ kF
. The cardinality of this ǫ-net will be
kX ♮ k
X♮

2

♮

a⊤
i

4

1/24
ǫ

≤ cn2 ·

!

2

V ⊤ ai

!2

kX ♮ k
.
kX ♮ kF

Therefore, when m ≥ c σ4

F
nr log (nκ) with some large enough constant c, for all matrices V and
(X ♮ )
♮
2
1 σr (X )
X such that X − X ♮ F ≤ 24
, we have
♮
kX kF
 


2
2
(50)
q (V , H, t) ≥ X ♮ F + 2 X ♮⊤ V F + 2Tr X ♮⊤ V X ♮⊤ V + 1.246σr2 X ♮ ,
r

with probability at least 1 − e−c1 nr log (nκ) − me−1.5n .

C.3 Finishing the Proof
Combining (42) and (50), we can prove
vec (V )⊤ ∇2 f (X)vec (V ) ≥ X ♮

2

+ 2 X ♮⊤ V

F
m

1 X ⊤ ♮
−
ai X
m
i=1

2
2

22

2
F

 


+ 2Tr X ♮⊤ V X ♮⊤ V + 1.246σr2 X ♮

a⊤
i V

2
2

 


+ 2Tr X ♮⊤ V X ♮⊤ V + 1.246σr2 X ♮
F
F
 
2
2
1
− σr2 X ♮
− X ♮ − 2 X ♮⊤ V
 F
 F 24  
♮⊤
♮⊤
≥ 2Tr X V X V + 1.204σr2 X ♮
2

≥ X♮

as claimed.

D

+ 2 X ♮⊤ V

2

Proof of Lemma 2

We first note that
Xt+1 Qt+1 − X ♮

2
F

≤ Xt+1 Qt − X ♮

2

(51)

F

= (Xt − µ∇f (Xt )) Qt − X ♮

2
F
2

= Xt Qt − µ∇f (Xt Qt ) − X ♮
F

 
= xt − x♮ − µ · vec ∇f (Xt Qt ) − ∇f X ♮

where we write

xt := vec (Xt Qt )

(52)
2
2

,

(53)


x♮ := vec X ♮ .

and

to the identity ∇f (Xt ) Qt =
Here, (51) follows from the definition of Qt+1 (see (13)), (52) holds owing

r×r
♮
∇f (Xt Qt ) for Qt ∈ O , and (53) arises from the fact that ∇f X = 0. Let

Xt (τ ) = X ♮ + τ Xt Qt − X ♮ ,

where τ ∈ [0, 1]. Then, by the fundamental theorem of calculus for vector-valued functions [43],
Z 1
2


RHS of (53) = xt − x♮ − µ ·
∇2 f (Xt (τ )) xt − x♮ dτ
0


Z
=
I −µ·

1

0

(54)

2



xt − x♮
∇ f (Xt (τ )) dτ

2

2

2

2 

2
xt − x♮
∇ f (Xt (τ )) dτ
0



⊤ Z 1
2
2
♮
♮
xt − x♮
∇ f (Xt (τ )) dτ
− 2µ · xt − x
= xt − x
Z
⊤ 

♮
= xt − x
I −µ·

1

2

0

2 

xt − x♮
∇ f (Xt (τ )) dτ
0


⊤ Z 1

2
2
♮
♮
xt − x♮
∇ f (Xt (τ )) dτ
− 2µ · xt − x
≤ xt − x
⊤ Z

2
♮
+ µ · xt − x

1

2

2

+ µ2 ·

Z

0

1

0

2

∇2 f (Xt (τ )) dτ

xt − x♮

It is easy to verify that Xt (τ ) satisfies (20) for any τ ∈ [0, 1], since
Xt (τ ) − X

♮
F

= τ Xt Qt − X
23

♮
F

2
2

.


1 σr2 X ♮
≤
,
24 kX ♮ kF

(55)

and
max

1≤l≤m



Xt (τ ) − X ♮
a⊤
l

Lemma 1 then implies that
⊤ Z

♮
xt − x
and

Z

1
0

1
0



= τ · max

1≤l≤m

2

2

∇ f (Xt (τ )) dτ



Xt Qt − X ♮
a⊤
l



2


σr2 X ♮
1p
log n ·
≤
.
24
kX ♮ kF

 

xt − x♮ ≥ 1.026σr2 X ♮ xt − x♮

 
∇2 f (Xt (τ )) dτ ≤ 1.5σr2 X ♮ log n + 6 X ♮

2
2

,

2
.
F

Substituting the above two inequalities into (53) and (55) gives
Xt+1 Qt+1 − X ♮

2
F


 
+ µ · 1.5σr2 X ♮ log n + 6 X ♮
xt − x
X
− 2µ ·
≤ xt − x
2
2
#
"


 
 
2
2 2
µ2 Xt Qt − X ♮
= 1 − 2.052σr2 X ♮ µ + 1.5σr2 X ♮ log n + 6 X ♮
♮

2

1.026σr2



♮



♮

2

2

with the proviso that µ ≤

2

F

F

xt − x♮

,

1.026σr2 (X ♮ )

 .
2 2
1.5σr2 (X ♮ ) log n+6kX ♮ kF

Xt+1 Qt+1 − X ♮

F

2

F

F

  

≤ 1 − 1.026σr2 X ♮ µ Xt Qt − X ♮

2

This allows us to conclude that


 
≤ 1 − 0.513σr2 X ♮ µ Xt Qt − X ♮

F

.

E Proof of Lemma 3
Recognizing that
(l)

(l)

Xt+1 Qt+1 − Xt+1 Rt+1

(l)

F

(l)

≤ Xt+1 Qt+1 − Xt+1 Rt Q⊤
t Qt+1
= Xt+1 −
(l)

(l)
(l)
Xt+1 Rt Q⊤
t

(l)

we will focus on bounding Xt+1 Qt − Xt+1 Rt

F

F
(l)

F

(l)

= Xt+1 Qt − Xt+1 Rt

. Since




(l)
(l)
(l)
(l)
(l)
Rt
Xt+1 Qt − Xt+1 Rt = (Xt − µ∇f (Xt )) Qt − Xt − µ∇f (l) Xt


(l) (l)
(l)
(l)
= Xt Qt − Xt Rt − µ∇f (Xt ) Qt + µ∇f (l) Xt Rt

m 
2
1 X
(l) (l)
⊤
ai Xt − yi ai a⊤
= Xt Qt − Xt Rt − µ
i Xt Qt
m
2
i=1



m 
2
2
1
1 X
(l) (l)
⊤ (l)
⊤ (l)
⊤ (l) (l)
− yl al a⊤
− yi ai ai Xt Rt − µ
al Xt
ai Xt
+µ
l Xt Rt
m
m
2
2
i=1

24

F

,

2
2




1
(l)
(l) (l)
(l) (l)
a⊤
= Xt Qt − Xt Rt − µ∇f (Xt Qt ) + µ∇f Xt Rt − µ
l Xt
m
|
{z
} |
(l)

:=St,1

(l)

(l)

we aim to control St,1

F

and St,2

We first bound the term

(l)
St,2 F ,

F


(l) (l)
− yl al a⊤
l Xt Rt ,
2
{z
}
2

(l)

:=St,2

separately.

which is easier to handle. Observe that by Cauchy-Schwarz,

⊤


(l) (l)
(l) (l)
− yl = a⊤
Xt Rt − X ♮ Xt Rt + X ♮ al
l
2




(l) (l)
(l) (l)
♮
♮
⊤
−
X
+
X
R
R
X
X
a
≤ a⊤
l
l
t
t
t
t

(l) 2

a⊤
l Xt

2

The first term in (56) can be bounded by


(l) (l)
♮
−
X
R
X
a⊤
l
t
t
2




(l) (l)
♮
⊤
X
Q
−
X
+ a⊤
≤ al Xt Rt − Xt Qt
t t
l
2

2

.

(56)

2


  t p

2 X♮
σ
(l) (l)
r
≤ 6n Xt Rt − Xt Qt + C2 1 − 0.5σr2 X ♮ µ
log n ·
kX ♮ kF
r


  t p

  t log n σ 2 X ♮

√
σr2 X ♮
r
♮
2
♮
2
≤ 6nC3 1 − 0.5σr X µ
+ C2 1 − 0.5σr X µ
·
log n ·
n
kX ♮ kF
kX ♮ kF

  t p

√
σ2 X ♮
(57)
= ( 6C3 + C2 ) 1 − 0.5σr2 X ♮ µ
log n r ♮ ,
kX kF
√

where we have used the triangle inequality, Lemma 10, as well as the induction hypotheses (28c) and (28b).
Similarly, the second term in (56) can be bounded as




(l) (l)
(l) (l)
♮
♮
⊤
♮
X
R
−
X
+ 2 a⊤
a
X
R
+
X
≤
a⊤
l X
l
l
t
t
t
t
2
2
2

√
p
p
σr2 X ♮
≤
+
11.72
6C3 + C2
log n
log n X ♮
kX ♮ kF
F
p
√
≤
6C3 + C2 + 11.72
log n X ♮ F ,
(58)


2
where we have used (57), Lemma 11, and σr2 X ♮ ≤ X ♮ F . Similarly, we can also obtain
p
√
(l)
6C3 + C2 + 5.86
log n X ♮ F .
≤
a⊤
l Xt
2

Substituting (57) and (58) into (56), and using the above inequality, we get
(l)

St,2

F

1
(l)
(l) 2
− yl · al a⊤
· a⊤
l Xt
l Xt
m
2
F
  t



1
(l)
≤ C42 1 − 0.5σr2 X ♮ µ · µ · σr2 X ♮ log n · kal k2 a⊤
l Xt
m
2
  t
 
p
√ 3
√
1
♮
2
♮
♮
2
≤ 6C4 1 − 0.5σr X µ · µ · σr X log n · n X F log n
m
  t √n · (log n)3/2  
√ 3
σr2 X ♮ X ♮ F ,
= 6C4 1 − 0.5σr2 X ♮ µ · µ
m
=µ

25

(59)

√
where C4 := 6C3 + C2 + 11.72.
(l)
Next, we turn to St,1 . By defining
F

we can write

(l)
(l) 
st,1 = vec St,1 ,

xt = vec (Xt Qt ) ,

and

(l) 

(l)

(l)

xt = vec Xt Rt

,



(l) (l)
(l)
(l)
st,1 = xt − xt − µ · vec ∇f (Xt Qt ) − ∇f (Xt Rt )
Z 1



(l)
(l)
(l)
∇2 f Xt (τ ) xt − xt dτ
= xt − xt − µ ·
0

Z 1

 

(l)
(l)
2
xt − xt .
∇ f Xt (τ ) dτ
= I −µ·
0

Here, the second line follows from the fundamental theorem of calculus for vector-valued functions [43],
where


(l) (l)
(l) (l)
(l)
(60)
Xt (τ ) = Xt Rt + τ Xt Qt − Xt Rt ,
for τ ∈ [0, 1]. Using very similar algebra as in Appendix D, we obtain
(l) 2
St,1
F

≤ xt −

(l) 2
xt
2

+µ

2

Z

1

0

 Z

(l) ⊤
− 2µ · xt − xt



(l)
∇2 f Xt (τ ) dτ
1

0

2

(l) 2

xt − xt

2


 

(l)
(l)
2
xt − xt .
∇ f Xt (τ ) dτ

(61)

It is easy to verify that for all τ ∈ [0, 1],
(l)

Xt (τ ) − X ♮

F



(l) (l)
= (1 − τ ) Xt Rt − Xt Qt + Xt Qt − X ♮
(l) (l)
Xt Rt

F
♮

− Xt Qt + Xt Qt − X
≤ (1 − τ )
F
r


2
♮
σr2 X ♮
log n σr X
≤ C3
+ C1
·
n
kX ♮ kF
kX ♮ kF
!
r


σr2 X ♮
1 σr2 X ♮
log n
= C3
≤
,
+ C1
n
kX ♮ kF
24 kX ♮ kF

F

(62)
(63)

where (62) follows from the induction hypotheses (28a) and (28b), and (63) follows as long as C1 +C3 ≤
Further, for all 1 ≤ l ≤ m, by the induction hypothesis (28b) and (28c),






(l) (l)
(l)
♮
⊤
⊤
♮
X
Q
−
X
a
−
X
Q
+
R
X
a
(τ
)
−
X
≤
(1
−
τ
)
X
a⊤
t
t
t
t
l
l
l
t
t
t
2
2
2

2 X♮
p
σ
(l) (l)
≤ kal k2 Xt Rt − Xt Qt + C2 log n · r ♮
kX kF
r


p
√
σr2 X ♮
log n σr2 X ♮
≤ 6nC3
+ C2 log n ·
·
n
kX ♮ kF
kX ♮ kF


p
√
σr2 X ♮
σr2 X ♮
1p
log n ·
≤
,
6C3 + C2
log n ·
≤
kX ♮ kF
24
kX ♮ kF
26

1
24 .

√
as long as 6C3 + C2 ≤
be further bounded by

1
24 .

(l)

Therefore, Lemma 1 holds for Xt (τ ), and similar to Appendix D, (61) can

(l)

St,1
as long as µ ≤

F

  

(l) (l)
≤ 1 − 0.513σr2 X ♮ µ Xt Qt − Xt Rt

1.026σr2 (X ♮ )

 .
2 2
1.5σr2 (X ♮ ) log n+6kX ♮ kF
(l)

(l)

Xt+1 Qt+1 − Xt+1 Rt+1

Consequently, combining (59) and (64), we can get

(l)

F

(64)

F

(l)

≤ St,1 + St,2
F
 F 

(l) (l)
2
≤ 1 − 0.513σr X ♮ µ Xt Qt − Xt Rt
F
  t √n · (log n)3/2  
√ 3
+ 6C4 1 − 0.5σr2 X ♮ µ · µ
σr2 X ♮ X ♮ F
m
  t+1 r log n σ 2 X ♮ 

≤ C3 1 − 0.5σr2 X ♮ µ
· r ♮ ,
(65)
n
kX kF
2

kX ♮ k
where (65) follows from the induction hypothesis (28b), as long as m ≥ c σ2 X ♮F n log n for some large
)
r(
enough constant c > 0.

F

Proof of Lemma 4
(l)

For any 1 ≤ l ≤ m, by the statistical independence of al and Xt+1 and by Lemma 11, we have


p
(l)
(l)
(l)
(l)
♮
≤ 5.86 log n Xt+1 Rt+1 − X ♮ .
X
R
−
X
a⊤
l
t+1 t+1
2

F

Further, by the triangle inequality, Lemma 10, Lemma 3 and Lemma 2, we can deduce that






(l)
(l)
(l)
(l)
♮
⊤
⊤
♮
X
R
−
X
X
Q
−
X
R
a
a
+
X
Q
−
X
≤
a⊤
t+1 t+1
t+1 t+1
l
l
l
t+1 t+1
t+1 t+1
2
2
2
p
(l)
(l)
(l)
(l)
♮
≤ kal k2 Xt+1 Qt+1 − Xt+1 Rt+1 + 5.86 log n Xt+1 Rt+1 − X
F
p
√
(l)
(l)
(l)
(l)
≤ 6n Xt+1 Qt+1 − Xt+1 Rt+1 + 5.86 log n · Xt+1 Qt+1 − Xt+1 Rt+1
F
p
♮
+ 5.86 log n · Xt+1 Qt+1 − X
F
  t+1 r log n σ 2 X ♮ 
 
√
p
2
6n + 5.86 log n C3 1 − 0.5σr X ♮ µ
≤
· r ♮
n
kX kF





2
♮
p
t+1 σr X
·
+ 5.86 log nC1 1 − 0.513σr2 X ♮ µ
kX ♮ kF
!
r

  t+1 p
√
σr2 X ♮
log n 
♮
2
6C3 + 5.86C1 + 5.86C3
log n ·
≤
1 − 0.5σr X µ
n
kX ♮ kF

  t+1 p

σ2 X ♮
log n · r ♮ ,
≤ C2 1 − 0.5σr2 X ♮ µ
kX kF
√
where the last line follows as long as 6C3 +5.86C1 +5.86C3 ≤ C2 . The proof is then finished by applying
the union bound for all 1 ≤ l ≤ m.
27

G

Proof of Lemma 5

Define
Σ0 = diag {λ1 (Y ) , λ2 (Y ) , · · · , λr (Y )} = Λ0 + λI
o




n 
(l)
(l)
= Λ0 + λ(l) I,
Σ0 = diag λ1 Y (l) , λ2 Y (l) , · · · , λr Y (l)
(l)

(l)

1 ≤ l ≤ m,

(l)

then by definition we have Y Z0 = Z0 Σ0 , Y (l) Z0 = Z0 Σ0 , and
(l)

(l)

(l)

Σ0 Z0⊤ Z0 − Z0⊤ Z0 Σ0 =

1
(l)
yl Z0⊤ al a⊤
l Z0 .
2m

(l)

(66)

(l)

Moreover, let Z0,c and Z0,c be the complement matrices of Z0 and Z0 , respectively, such that both
h
i
(l)
(l)
[Z0 , Z0,c ] and Z0 , Z0,c are orthonormal matrices. Below we will prove the induction hypotheses (28) in
the base case when t = 0 one by one.

G.1 Proof of (28a)
From Lemma 6, we have
X0 Q0 − X ♮

F

≤q

2

√

1

X0 X0⊤ − X ♮ X ♮⊤

♮
2 − 1 σr (X )

F

1

Z0 Λ0 Z0⊤ − X ♮ X ♮⊤

F
2 2 − 1 σr (X ♮ )
√
r
≤q √
Z0 Σ0 Z0⊤ − X ♮ X ♮⊤ − λZ0 Z0⊤ .

♮
2 2 − 1 σr (X )
=q

√

(67)

The last term in (67) can be further bounded as
Z0 Σ0 Z0⊤ − X ♮ X ♮⊤ − λZ0 Z0⊤

1
1
1
⊤
kX ♮ k2F Z0 Z0⊤ − λZ0 Z0⊤
+
≤ Y − kX ♮ k2F I − X ♮ X ♮⊤ + Z0 Σ0 Z0⊤ − Y + kX ♮ k2F Z0,c Z0,c
2
2
2
≤ δ X♮

2
F

+ δ X♮

2
F

+ δ X♮

2
F

= 3δ X ♮

2
F

,

(68)

where (68) follows from
kY − E[Y ]k = Y −

1
X♮
2

2
F

I − X ♮ X ♮⊤ ≤ δ X ♮

via Lemma 13, the Weyl’s inequality, and
|λ − E [λ]| = λ −

1
X♮
2

2
F

≤ δ X♮

2
F

via Lemma 9. Plugging (68) into (67), we have
X0 Q0 − X

♮
F

√
2
δ r X♮ F
3
≤q √
 · σr (X ♮ ) ,
2 2−1
28

2
F

6

kX ♮ k
for a sufficiently small constant c, i.e. m & σ6 X ♮F nr log n, we get X0 Q0 − X ♮
Setting δ =
)
r(
σr2 (X ♮ )
σ2 (X ♮ )
(l) (l)
≤ C1 X ♮ . Following similar procedures, we can also show X0 Q0 − X ♮ ≤ C1 rX ♮ .
k kF
k kF
F
σ 3 (X ♮ )
c√ r ♮ 3
r kX kF

F

G.2 Proof of (28b)
Following Weyl’s inequality, by (28a), we have

 
σr2 X ♮
♮
≤ C1
σi (X0 ) − σi X
,
kX ♮ kF



σ2 (X ♮ )
(l)
and similarly, σi X0 − σi X ♮ ≤ C1 rX ♮ , for i = 1, · · · , r. Combined with Lemma 6, there
k kF
exists some constant c such that
(l)

(l)

X0 Q0 − X0 R0
c
σr (X ♮ )
c
=
σr (X ♮ )
c
=
σr (X ♮ )
c
≤
σr (X ♮ )
≤

F

≤q

2

(l)

√

1

(l)

F

(l)⊤

X0 X0⊤ − X0 X0

F

Z0 Λ0 Z0⊤ −

(l) (l) (l)⊤
Z0 Λ0 Z0

Z0 Σ0 Z0⊤

(l) (l) (l)⊤
Z0 Σ0 Z0

−

(l)⊤

X0 X0⊤ − X0 X0

2 − 1 σr (X0 )

(l)

(l)

(l)⊤

Z0 Σ0 Z0⊤ − Z0 Σ0 Z0

F
(l)

(l)⊤

− λZ0 Z0⊤ + λ(l) Z0 Z0

F

c
(l) (l)⊤
+
λZ0 Z0⊤ − λ(l) Z0 Z0
σr (X ♮ )
F

We will bound each term in (69), respectively. For the first term, we have
h
i
(l)
(l) (l)
(l)
(l) (l) (l)⊤
Z0 Σ0 Z0⊤ − Z0 Σ0 Z0
= Z0 Σ0 Z0⊤ Z0 − Z0 Σ0 , Z0 Σ0 Z0⊤ Z0,c
F

≤

≤

(l)
Z0 Σ0 Z0⊤ Z0

(l) (l)
− Z0 Σ0

+

F
⊤ (l)
⊤ (l) (l)
Z0 Σ0 Z0 Z0 − Z0 Z0 Z0 Σ0

1
(l)
≤ Z0 ·
yl Z0⊤ al a⊤
l Z0
2m

.

(69)

F

(l)
Z0 Σ0 Z0⊤ Z0,c

F

+

F
⊤ (l) (l)
Z0 Z0 Z0 Σ0
(l)

F

F

(l)⊤

+ Z0 Z0⊤ − Z0 Z0

(l)

(l)

(l)

− Z0 Σ0

F

+ kY k Z0⊤ Z0,c
(l)

F

Y (l) + kY k Z0⊤ Z0,c

F

,

F

(70)

where the last line follows from (66). Note that the first term in (70) can be bounded as
Z0 ·

1
(l)
yl Z0⊤ al a⊤
l Z0
2m

F

2
1
(l)
♮
a⊤
a⊤
a⊤
l X
l Z0
l Z0
2m
2
2
√
√
n · (log n)3/2 · r
2
X ♮ F,
.
m

≤

2

(71)

which follows Lemma 10 and Lemma 11. The second term in (70) can be bounded as
(l)

(l)⊤

Z0 Z0⊤ − Z0 Z0

F






(l) (l) ⊤
(l) (l)
(l) (l) ⊤
= Z0 Z0 − Z0 T0
+ Z0 − Z0 T0
Z0 T0

≤ 2 Z0 −

(l) (l)
Z0 T0

29

F

√
(l)
≤ 2 2 Z0⊤ Z0,c

F

F

(l)

(l)

(l)

(l)

where Tt = argminP ∈Or×r Zt − Zt P , and the last line follows from the fact Z0 − Z0 T0
F
√
(l)
2 Z0⊤ Z0,c [44]. Putting this together with the third term in (70), we have

F

≤

F

(l)

(l)⊤

Z0 Z0⊤ − Z0 Z0


 √
(l)
≤ 2 2 Y (l) + kY k Z0⊤ Z0,c
F
 (l)
1
⊤
m yl al al Z0
(l)

Y (l) + kY k Z0⊤ Z0,c

F

. X♮

.
.

2

F

σr2 (X ♮ )

F
2
2

♮
a⊤
l X

(l)

a⊤
l Z0

2

kal k2

m
√

n · (log n)3/2 ·
m

√
r

F

(72)

2

X♮ F
σr2 (X ♮ )

4

X♮ F
,
σr2 (X ♮ )

(73)

where (72) follows from Lemma 13 and the Davis-Kahan sin Θ theorem [45], and (73) follows from
Lemma 10 and Lemma 11.
For the second term in (69), we have
(l)

(l)⊤

λZ0 Z0⊤ − λ(l) Z0 Z0

(l)

F

(l)⊤

= λZ0 Z0⊤ − λZ0 Z0

(l)

(l)⊤

− λ(l) Z0 Z0

(l)

(l) (l)⊤
Z0 Z0

+ λZ0 Z0

(l) (l)⊤
Z0 Z0

Z0 Z0⊤

(l)

≤λ·
−
+ λ−λ ·
F
√
√
4
yl √
n · (log n)3/2 · r X ♮ F
+
r
.
2
♮
m
σr (X ) 2m
√
√
√
4
n · (log n)3/2 · r X ♮ F
r · log n
+
X♮
.
m
σr2 (X ♮ )
m

(l)⊤
F
F

(74)
2
F

,

(75)

where the first term of (74) is bounded similarly as (73), and (75) follows from Lemma 11. Combining (71),
(73), and (75), we obtain
(l) (l)
X0 Q0 − X0 R0

F

.

√

n · (log n)3/2 ·
m

√

4

r X♮ F
.
σr3 (X ♮ )

r


log n σr2 X ♮
,
·
n
kX ♮ kF

5

where the last inequality holds as long as m &

kX ♮ kF √
n r log n = nr 3 log n.
σr5 (X ♮ )

G.3 Proof of (28c)
For every 1 ≤ l ≤ m, from (28a) and (28b), we have
(l) (l)
X0 R0

−X

♮
F

≤

(l) (l)
X0 R0

− X0 Q0

F

+ X0 Q0 − X

♮
F


σr2 X ♮
.
.
kX ♮ kF

This further gives


♮
X
Q
−
X
a⊤
0 0
l
1≤l≤m
2


(l)
(l)
X
Q
−
X
R
≤ max a⊤
0
0
l
0
0
max

1≤l≤m

2

30

+ max

1≤l≤m



(l) (l)
♮
X
R
−
X
a⊤
l
0
0

2

(l)

(l)

≤ max kal k2 X0 Q0 − X0 R0
1≤l≤m

√

+ max

1≤l≤m

p

(l) (l)
X0 R0



(l) (l)
X0 R0 − X ♮
a⊤
l

n · max X0 Q0 −
+ log n · max
1≤l≤m
1≤l≤m
r


2
♮
2
♮
p
√
σ X
log n σr X
+ log n · r ♮
·
. n·
♮
n
kX kF
kX kF

2
♮
p
σ X
. log n · r ♮ ,
kX kF

.

(l) (l)
X0 R0

−X

2
♮
2

(76)
(77)

where (76) follows from Lemma 10 and Lemma 11, and (77) follows from (28b).

G.4 Finishing the Proof
The proof of Lemma 5 is now complete by appropriately adjusting the constants .

References
[1] Y. Chen, Y. Chi, and A. J. Goldsmith, “Exact and stable covariance estimation from quadratic sampling
via convex programming,” IEEE Transactions on Information Theory, vol. 61, no. 7, pp. 4034–4059,
2015.
[2] E. J. Candès, X. Li, and M. Soltanolkotabi, “Phase retrieval via Wirtinger flow: Theory and algorithms,”
IEEE Transactions on Information Theory, vol. 61, no. 4, pp. 1985–2007, 2015.
[3] Y. Chen, X. Yi, and C. Caramanis, “A convex formulation for mixed regression with two components:
Minimax optimal rates,” in Conference on Learning Theory, 2014, pp. 560–604.
[4] R. Kueng, H. Rauhut, and U. Terstiege, “Low rank matrix recovery from rank one measurements,”
Applied and Computational Harmonic Analysis, vol. 42, no. 1, pp. 88–116, 2017.
[5] L. Tian, J. Lee, S. B. Oh, and G. Barbastathis, “Experimental compressive phase space tomography,”
Optics express, vol. 20, no. 8, pp. 8296–8308, 2012.
[6] R. Livni, S. Shalev-Shwartz, and O. Shamir, “On the computational efficiency of training neural networks,” in Advances in Neural Information Processing Systems, 2014, pp. 855–863.
[7] M. Soltanolkotabi, A. Javanmard, and J. D. Lee, “Theoretical insights into the optimization landscape
of over-parameterized shallow neural networks,” arXiv preprint arXiv:1707.04926, 2017.
[8] M. Soltani and C. Hegde, “Towards provable learning of polynomial neural networks using low-rank
matrix estimation,” in Proceedings of the 21st International Conference on Artificial Intelligence and
Statistics, 2018.
[9] K. Zhong, P. Jain, and I. S. Dhillon, “Efficient matrix sensing using rank-1 Gaussian measurements,”
in International Conference on Algorithmic Learning Theory. Springer, 2015, pp. 3–18.
[10] M. Lin and J. Ye, “A non-convex one-pass framework for generalized factorization machine and rankone matrix sensing,” in Advances in Neural Information Processing Systems, 2016, pp. 1633–1641.
[11] S. Sanghavi, R. Ward, and C. D. White, “The local convexity of solving systems of quadratic equations,” Results in Mathematics, vol. 71, no. 3-4, pp. 569–608, 2017.
31

[12] Y. Chen and E. J. Candès, “Solving random quadratic systems of equations is nearly as easy as solving
linear systems,” Comm. Pure Appl. Math., vol. 70, no. 5, pp. 822–883, 2017. [Online]. Available:
http://dx.doi.org/10.1002/cpa.21638
[13] Y. Chen and M. J. Wainwright, “Fast low-rank estimation by projected gradient descent: General
statistical and algorithmic guarantees,” arXiv preprint arXiv:1509.03025, 2015.
[14] C. Ma, K. Wang, Y. Chi, and Y. Chen, “Implicit regularization in nonconvex statistical estimation:
Gradient descent converges linearly for phase retrieval, matrix completion and blind deconvolution,”
arXiv preprint arXiv:1711.10467, 2017.
[15] Y. Zhong and N. Boumal, “Near-optimal bounds for phase synchronization,” arXiv preprint
arXiv:1703.06605, 2017.
[16] Y. Chen, J. Fan, C. Ma, and K. Wang, “Spectral method and regularized MLE are both optimal for
top-k ranking,” arXiv preprint arXiv:1707.09971, 2017.
[17] T. Cai and A. Zhang, “ROP: Matrix recovery via rank-one projections,” The Annals of Statistics, vol. 43,
no. 1, pp. 102–138, 2015.
[18] Y. Li, Y. Sun, and Y. Chi, “Low-rank positive semidefinite matrix recovery from corrupted rank-one
measurements,” IEEE Transactions on Signal Processing, vol. 65, no. 2, pp. 397–408, 2017.
[19] E. J. Candes, T. Strohmer, and V. Voroninski, “Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming,” Communications on Pure and Applied Mathematics,
vol. 66, no. 8, pp. 1241–1274, 2013.
[20] E. J. Candès and X. Li, “Solving quadratic equations via phaselift when there are about as many
equations as unknowns,” Foundations of Computational Mathematics, vol. 14, no. 5, pp. 1017–1026,
2014.
[21] L. Demanet and P. Hand, “Stable optimizationless recovery from phaseless linear measurements,” Journal of Fourier Analysis and Applications, vol. 20, no. 1, pp. 199–221, 2014.
[22] I. Waldspurger, A. d?Aspremont, and S. Mallat, “Phase recovery, maxcut and complex semidefinite
programming,” Mathematical Programming, vol. 149, no. 1-2, pp. 47–81, 2015.
[23] 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, vol. 44, no. 5, pp. 2221–2251, 2016.
[24] H. Zhang, Y. Zhou, Y. Liang, and Y. Chi, “A nonconvex approach for phase retrieval: Reshaped
wirtinger flow and incremental algorithms,” Journal of Machine Learning Research, 2017.
[25] M. Soltanolkotabi, “Structured signal recovery from quadratic measurements: Breaking sample complexity barriers via nonconvex optimization,” arXiv preprint arXiv:1702.06175, 2017.
[26] G. Wang, G. B. Giannakis, and Y. C. Eldar, “Solving systems of random quadratic equations via
truncated amplitude flow,” IEEE Transactions on Information Theory, 2017.
[27] 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 International Conference on Machine Learning-Volume 48. JMLR. org, 2016, pp. 964–973.

32

[28] X. Li, S. Ling, T. Strohmer, and K. Wei, “Rapid, robust, and reliable blind deconvolution via nonconvex
optimization,” CoRR, vol. abs/1606.04933, 2016. [Online]. Available: http://arxiv.org/abs/1606.04933
[29] J.-F. Cai, T. Wang, and K. Wei, “Spectral compressed sensing via projected gradient descent,” arXiv
preprint arXiv:1707.09726, 2017.
[30] K. Wei, “Solving systems of phaseless equations via Kaczmarz methods: A proof of concept study,”
Inverse Problems, vol. 31, no. 12, p. 125008, 2015.
[31] Y. S. Tan and R. Vershynin, “Phase retrieval via randomized kaczmarz: Theoretical guarantees,” arXiv
preprint arXiv:1706.09993, 2017.
[32] H. Jeong and C. S. Güntürk, “Convergence of the randomized kaczmarz method for phase retrieval,”
arXiv preprint arXiv:1706.10291, 2017.
[33] J. Ma, J. Xu, and A. Maleki, “Optimization-based AMP for phase retrieval: The impact of initialization
and ℓ2 -regularization,” arXiv preprint arXiv:1801.01170, 2018.
[34] Y. Chi and Y. M. Lu, “Kaczmarz method for solving quadratic equations,” IEEE Signal Processing
Letters, vol. 23, no. 9, pp. 1183–1187, 2016.
[35] J. Sun, Q. Qu, and J. Wright, “A geometric analysis of phase retrieval,” in Information Theory (ISIT),
2016 IEEE International Symposium on. IEEE, 2016, pp. 2379–2383.
[36] K. Zhong, Z. Song, P. Jain, P. L. Bartlett, and I. S. Dhillon, “Recovery guarantees for one-hidden-layer
neural networks,” in Proceedings of the 34th International Conference on Machine Learning, vol. 70,
2017, pp. 4140–4149.
[37] E. J. Candes and Y. Plan, “Tight oracle inequalities for low-rank matrix recovery from a minimal
number of noisy random measurements,” IEEE Transactions on Information Theory, vol. 57, no. 4, pp.
2342–2359, 2011.
[38] V. Bentkus, “An inequality for tail probabilities of martingales with differences bounded from one
side,” Journal of Theoretical Probability, vol. 16, no. 1, pp. 161–173, 2003.
[39] R. Vershynin, “Introduction to the non-asymptotic analysis of random matrices,” Compressed Sensing,
Theory and Applications, pp. 210 – 268, 2012.
[40] B. Laurent and P. Massart, “Adaptive estimation of a quadratic functional by model selection,” Annals
of Statistics, pp. 1302–1338, 2000.
[41] W. Schudy and M. Sviridenko, “Concentration and moment inequalities for polynomials of independent random variables,” arXiv preprint arXiv:1104.4997, 2011.
[42] J. M. Ten Berge, “Orthogonal procrustes rotation for two or more matrices,” Psychometrika, vol. 42,
no. 2, pp. 267–276, 1977.
[43] S. Lang, “Real and functional analysis, volume 142 of graduate texts in mathematics,” Springer-Verlag,
New York,, vol. 10, pp. 11–13, 1993.
[44] Y. Yu, T. Wang, and R. J. Samworth, “A useful variant of the davis–kahan theorem for statisticians,”
Biometrika, vol. 102, no. 2, pp. 315–323, 2014.
[45] C. Davis and W. M. Kahan, “The rotation of eigenvectors by a perturbation. iii,” SIAM Journal on
Numerical Analysis, vol. 7, no. 1, pp. 1–46, 1970.

33

