Fast Algorithms for Robust PCA via Gradient Descent
Xinyang Yi∗

Dohyung Park∗

Yudong Chen†

The University of Texas at Austin∗

arXiv:1605.07784v2 [] 19 Sep 2016

{yixy,dhpark,constantine}@utexas.edu∗

Constantine Caramanis∗

Cornell University†
yudong.chen@cornell.edu†

Abstract
We consider the problem of Robust PCA in the fully and partially observed settings. Without
corruptions, this is the well-known matrix completion problem. From a statistical standpoint
this problem has been recently well-studied, and conditions on when recovery is possible (how
many observations do we need, how many corruptions can we tolerate) via polynomial-time
algorithms is by now understood. This paper presents and analyzes a non-convex optimization
approach that greatly reduces the computational complexity of the above problems, compared
to the best available algorithms. In particular, in the fully observed case, with r denoting rank
and d dimension, we reduce the complexity from O(r2 d2 log(1/ε)) to O(rd2 log(1/ε)) – a big
savings when the rank is big. For the partially observed case, we show the complexity of our
algorithm is no more than O(r4 d log d log(1/ε)). Not only is this the best-known run-time for a
provable algorithm under partial observation, but in the setting where r is small compared to
d, it also allows for near-linear-in-d run-time that can be exploited in the fully-observed case as
well, by simply running our algorithm on a subset of the observations.

1

Introduction

Principal component analysis (PCA) aims to find a low rank subspace that best-approximates a
data matrix Y ∈ Rd1 ×d2 . The simple and standard method of PCA by singular value decomposition
(SVD) fails in many modern data problems due to missing and corrupted entries, as well as sheer
scale of the problem. Indeed, SVD is highly sensitive to outliers by virtue of the squared-error
criterion it minimizes. Moreover, its running time scales as O(rd2 ) to recover a rank r approximation
of a d-by-d matrix.
While there have been recent results developing provably robust algorithms for PCA (e.g., [8,
28]), the running times range from O(r2 d2 ) to O(d3 )1 and hence are significantly worse than SVD.
Meanwhile, the literature developing sub-quadratic algorithms for PCA (e.g., [16, 15, 3]) seems
unable to guarantee robustness to outliers or missing data.
Our contribution lies precisely in this area: provably robust algorithms for PCA with improved
run-time. Specifically, we provide an efficient algorithm with running time that matches SVD while
nearly matching the best-known robustness guarantees. In the case where rank is small compared
to dimension, we develop an algorithm with running time that is nearly linear in the dimension.
This last algorithm works by subsampling the data, and therefore we also show that our algorithm
solves the Robust PCA problem with partial observations (a generalization of matrix completion
and Robust PCA).
1

For precise dependence on error and other factors, please see details below.

1

1.1

The Model and Related Work

We consider the following setting for robust PCA. Suppose we are given a matrix Y ∈ Rd1 ×d2 that
has decomposition Y = M ∗ + S ∗ , where M ∗ is a rank r matrix and S ∗ is a sparse corruption matrix
containing entries with arbitrary magnitude. The goal is to recover M ∗ and S ∗ from Y . To ease
notation, we let d1 = d2 = d in the remainder of this section.
Provable solutions for this model are first provided in the works of [9] and [8]. They propose to
solve this problem by convex relaxation:
min |||M |||nuc + λkSk1 , s.t. Y = M + S,
M,S

(1)

where |||M |||nuc denotes the nuclear norm of M . Despite analyzing the same method, the corruption
models in [8] and [9] differ. In [8], the authors consider the setting where the entries of M ∗ are
corrupted at random with probability α. They show their method succeeds in exact recovery with
α as large as 0.1, which indicates they can tolerate a constant fraction of corruptions. Work in [9]
considers a deterministic corruption model, where nonzero entries of S ∗ can have arbitrary position,
but the sparsity of each row
√ and column does not exceed αd. They prove that for exact recovery,
it can allow α = O(1/(µr d)). This was subsequently further improved to α = O(1/(µr)), which
is in fact optimal [12, 19]. Here, µ represents the incoherence of M ∗ (see Section 2 for details). In
this paper, we follow this latter line and focus on the deterministic corruption model.
The state-of-the-art solver [21] for (1) has time complexity O(d3 /ε) to achieve error ε, and is
thus much slower than SVD, and prohibitive for even modest values of d. Work in [22] considers the
deterministic corruption model, and improves this running time without sacrificing the robustness
guarantee on α. They propose an alternating projection (AltProj) method to estimate the low
rank and sparse structures iteratively and simultaneously, and show their algorithm has complexity
O(r2 d2 log(1/ε)), which is faster than the convex approach but still slower than SVD.
Non-convex approaches have recently seen numerous developments for applications in low-rank
estimation, including alternating minimization (see e.g. [20, 18, 17]) and gradient descent (see e.g. [4,
11, 25, 26, 31, 32]). These works have fast running times, yet do not provide robustness guarantees.
One exception is [11], where the authors analyze a row-wise `1 projection method for recovering
S ∗ . Their analysis hinges on positive semidefinite M ∗ , and the algorithm requires prior knowledge
of the `1 norm of every row of S ∗ and is thus prohibitive in practice. Another exception is work
[17], which analyzes alternating minimization plus an overall sparse projection. Their algorithm is
shown to tolerate at most a fraction of α = O(1/(µ2/3 r2/3 d)) corruptions. As we discuss in Section
1.2, we can allow S ∗ to have much higher sparsity α = O(1/(µr1.5 )), which is close to optimal.
After the initial post of our paper on arxiv, Cherapanamjeri et al. [14] posted their paper on
solving robust PCA with partial observations by a method modified from AltProj. Their approach is
shown to have optimal robustness. The sample complexity they established is O(µ2 r2 d log2 d log2 (1/ε)),
which depends on the estimation error ε, since the analysis of AltProj under partial observations
requires sampling splitting. In contrast, the sample complexity O(µ2 r2 d log d) of our approach (see
Corollary 2) does not have such dependence. Moreover, their method requires computing rank-r
SVD of the observed matrices in every round of iterations. Our algorithm only needs to compute
rank-r SVD (approximately) once in the initialization step.
It is worth mentioning other works that obtain provable guarantees of non-convex algorithms or
problems including phase retrieval [5, 13, 30], EM algorithms [2, 27, 29], tensor decompositions [1]
and second order method [24]. It might be interesting to bring robust considerations to these works.
2

1.2

Our Contributions

In this paper, we develop efficient non-convex algorithms for robust PCA. We propose a novel
algorithm based on the projected gradient method on the factorized space. We also extend it to
solve robust PCA in the setting with partial observations, i.e., in addition to gross corruptions,
the data matrix has a large number of missing values. Our main contributions are summarized as
follows.2
1. We propose a novel sparse estimator for the setting of deterministic corruptions. For the lowrank structure to be identifiable, it is natural to assume that deterministic corruptions are
“spread out” (no more than some number in each row/column). We leverage this information
in a simple but critical algorithmic idea, that is tied to the ultimate complexity advantages
our algorithm delivers.
2. Based on the proposed sparse estimator, we propose a projected gradient method on the matrix
factorized space. While non-convex, the algorithm is shown to enjoy linear convergence under
proper initialization. Along with a new initialization method, we show that robust PCA can
be solved within complexity O(rd2 log(1/ε)) while ensuring robustness α = O(1/(µr1.5 )). Our
algorithm is thus faster than the best previous known algorithm by a factor of r, and enjoys
superior empirical performance as well.
3. Algorithms for Robust PCA with partial observations still rely on a computationally expensive convex approach, as apparently this problem has evaded treatment by non-convex methods. We consider precisely this problem. In a nutshell, we show that our gradient method
succeeds (it is guaranteed to produce the subspace of M ∗ ) even when run on no more than
O(µ2 r2 d log d) random entries of Y . The computational cost is O(µ3 r4 d log d log(1/ε)). When
rank r is small compared to the dimension d, in fact this dramatically improves on our bound
above, as our cost becomes nearly linear in d. We show, moreover, that this savings and robustness to erasures comes at no cost in the robustness guarantee for the deterministic (gross)
corruptions. While this demonstrates our algorithm is robust to both outliers and erasures, it
also provides a way to reduce computational costs even in the fully observed setting, when r
is small.
4. An immediate corollary of the above result provides a guarantee for exact matrix completion,
with general rectangular matrices, using O(µ2 r2 d log d) observed entries and O(µ3 r4 d log d log(1/ε))
time, thereby improving on existing results in [11, 25].

1.3

Organization and Notation

The remainder of this paper is organized as follows. In Section 2, we formally describe our problem
and assumptions. In Section 3, we present and describe our algorithms for fully (Algorithm 1) and
partially (Algorithm 2) observed settings. In Section 4.1, we establish theoretical guarantees of
Algorithm 1. The theory for partially observed setting are presented in Section 4.2. The numerical
results are collected in Section 5. Sections 6, 7 and Appendix A contain all the proofs and technical
lemmas.
To ease presentation, the discussion here assumes M ∗ has constant condition number, whereas our results below
show the dependence on condition number explicitly.
2

3



For any index set Ω ⊆ [d1 ]×[d2 ], we let Ω(i,·) := (i, j) ∈ Ω j ∈ [d2 ] , Ω(·,j) := (i, j) ∈ Ω i ∈ [d1 ] .
For any matrix A ∈ Rd1 ×d2 , we denote its projector onto support Ω by ΠΩ (A), i.e., the (i, j)-th
entry of ΠΩ (A) is equal to A if (i, j) ∈ Ω and zero otherwise. The i-th row and j-th column of
A are denoted by A(i,·) and A(·,j) . The (i, j)-th entry is denoted as A(i,j) . Operator norm of A is
|||A|||op . Frobenius norm of A is |||A|||F . The `a /`b norm of A is denoted by |||A|||b,a , i.e., the `a norm of
the vector formed by the `b norm of every row. For instance, |||A|||2,∞ stands for maxi∈[d1 ] kA(i,·) k2 .

2

Problem Setup

We consider the problem where we observe a matrix Y ∈ Rd1 ×d2 that satisfies Y = M ∗ + S ∗ , where
M ∗ has rank r, and S ∗ is corruption matrix with sparse support. Our goal is to recover M ∗ and S ∗ .
In the partially observed setting, in addition to sparse corruptions, we have erasures. We assume
that each entry of M ∗ + S ∗ is revealed independently with probability p ∈ (0, 1). In particular, for
any (i, j) ∈ [d1 ] × [d2 ], we consider the Bernoulli model where
(
(M ∗ + S ∗ )(i,j) , with probability p;
Y(i,j) =
(2)
∗,
otherwise.
We denote the support of Y by Φ = {(i, j) | Y(i,j) 6= ∗}. Note that we assume S ∗ is not adaptive to
Φ. As is well-understood thanks to work in matrix completion, this task is impossible in general –
we need to guarantee that M ∗ is not both low-rank and sparse. To avoid such identifiability issues,
we make the following standard assumptions on M ∗ and S ∗ : (i) M ∗ is not near-sparse or “spiky.”
We impose this by requiring M ∗ to be µ-incoherent – given a singular value decomposition (SVD)3
M ∗ = L∗ Σ∗ R∗> , we assume that
r
r
µr
µr
∗
∗
|||L |||2,∞ ≤
, |||R |||2,∞ ≤
.
d1
d2
(ii) The entries of S ∗ are “spread out” – for α ∈ [0, 1), we assume S ∗ ∈ Sα , where
n
o
Sα := A ∈ Rd1 ×d2 kA(i,·) k0 ≤ αd2 for all i ∈ [d1 ] ; kA(·,j) k0 ≤ αd1 for all j ∈ [d2 ] .

(3)

In other words, S ∗ contains at most α-fraction nonzero entries per row and column.

3

Algorithms

For both the full and partial observation settings, our method proceeds in two phases. In the first
phase, we use a new sorting-based sparse estimator to produce a rough estimate Sinit for S ∗ based on
the observed matrix Y , and then find a rank r matrix factorized as U0 V0> that is a rough estimate
of M ∗ by performing SVD on (Y − Sinit ). In the second phase, given (U0 , V0 ), we perform an
iterative method to produce series {(Ut , Vt )}∞
t=0 . In each step t, we first apply our sparse estimator
to produce a sparse matrix St based on (Ut , Vt ), and then perform a projected gradient descent
step on the low-rank factorized space to produce (Ut+1 , Vt+1 ). This flow is the same for full and
partial observations, though a few details differ. Algorithm 1 gives the full observation algorithm,
and Algorithm 2 gives the partial observation algorithm. We now describe the key details of each
algorithm.
3

Throughout this paper, we refer to SVD of rank r matrix by form LΣR> where Σ ∈ Rr×r .

4

Sparse Estimation. A natural idea is to keep those entries of residual matrix Y − M that have
large magnitude. At the same time, we need to make use of the dispersed property of Sα that every
column and row contain at most α-fraction of nonzero entries. Motivated by these two principles, we
introduce the following sparsification operator: For any matrix A ∈ Rd1 ×d2 : for all (i, j) ∈ [d1 ]×[d2 ],
we let
(
(αd )
(αd )
A(i,j) , if |A(i,j) | ≥ |A(i,·)2 | and |A(i,j) | ≥ |A(·,j)1 |
Tα [A] :=
,
(4)
0,
otherwise
(k)

(k)

where A(i,·) and A(·,j) denote the elements of A(i,·) and A(·,j) that have the k-th largest magnitude
respectively. In other words, we choose to keep those elements that are simultaneously among
the largest α-fraction entries in the corresponding row and column. In the case of entries having
identical magnitude, we break ties arbitrarily. It is thus guaranteed that Tα [A] ∈ Sα .
Algorithm 1 Fast RPCA
INPUT: Observed matrix Y with rank r and corruption fraction α; parameters γ, η; number of
iterations T .
// Phase I: Initialization.
1: Sinit ← Tα [Y ]
// see (4) for the definition of Tα [·].
2: [L, Σ, R] ← SVDr [Y − Sinit ] 4
3: U0 ← LΣ1/2 , V0 ← RΣ1/2 . Let U, V be defined according to (7).
// Phase II: Gradient based iterations.
4: U0 ← ΠU (U0 ), V0 ← ΠV (V0 )
5: for t = 0, 1, . . . , T − 1 do


6:
St ← Tγα Y − Ut Vt>

7:
Ut+1 ← ΠU Ut − η∇U L(Ut , Vt ; St ) − 12 ηUt (Ut> Ut − Vt> Vt )
8:
Vt+1 ← ΠV Vt − η∇V L(Ut , Vt ; St ) − 21 ηVt (Vt> Vt − Ut> Ut )
9: end for
OUTPUT: (UT , VT )

Initialization. In the fully observed setting, we compute Sinit based on Y as Sinit = Tα [Y ]. In the
partially observed setting with sampling rate p, we let Sinit = T2pα [Y ]. In both cases, we then set
U0 = LΣ1/2 and V0 = RΣ1/2 , where LΣR> is an SVD of the best rank r approximation of Y − Sinit .
Gradient Method on Factorized Space. After initialization, we proceed by projected gradient
descent. To do this, we define loss functions explicitly in the factored space, i.e., in terms of U, V
and S:
L(U, V ; S) :=
e V ; S) :=
L(U,

1
|||U V > + S − Y |||2F ,
(fully observed)
2


1
|||ΠΦ U V > + S − Y |||2F .
(partially observed)
2p

(5)
(6)

4
SVDr [A] stands for computing a rank-r SVD of matrix A, i.e., finding the top r singular values and vectors of A.
Note that we only need to compute rank-r SVD approximately (see the initialization error requirement in Theorem
1) so that we can leverage fast iterative approaches such as block power method and Krylov subspace methods [23].

5

Recall that our goal is to recover M ∗ that satisfies the µ-incoherent condition. Given an SVD
M ∗ = L∗ ΣR∗> , we expect that the solution (U, V ) is close to (L∗ Σ1/2 , R∗ Σ1/2 ) up to some rotation.
In order to serve such µ-incoherent structure, it is natural to put constraints on the row norms of
U, V based on |||M ∗ |||op . As |||M ∗ |||op is unavailable, given U0 , V0 computed in the first phase, we rely
on the sets U, V defined as
r
r




2µr
2µr
d2 ×r
d1 ×r
|||A|||2,∞ ≤
U := A ∈ R
|||A|||2,∞ ≤
|||U0 |||op , V := A ∈ R
|||V0 |||op . (7)
d1
d2
Now we consider the following optimization problems with constraints:
min

U ∈U ,V ∈V,S∈Sα

min

U ∈U ,V ∈V,S∈Spα

1
L(U, V ; S) + |||U > U − V > V |||2F ,
8
e V ; S) + 1 |||U > U − V > V |||2 .
L(U,
F
64

(fully observed)
(partially observed)

(8)
(9)

The regularization term in the objectives above is used to encourage that U and V have the same
scale. Given (U0 , V0 ), we propose the following iterative method to produce series {(Ut , Vt )}∞
t=0 and
{St }∞
.
We
give
the
details
for
the
fully
observed
case
–
the
partially
observed
case
is
similar.
For
t=0


>
t = 0, 1, . . ., we update St using the sparse estimator St = Tγα Y − Ut Vt , followed by a projected
gradient update on Ut and Vt


1
>
>
Ut+1 = ΠU Ut − η∇U L(Ut , Vt ; St ) − ηUt (Ut Ut − Vt Vt ) ,
2


1
>
>
Vt+1 = ΠV Vt − η∇V L(Ut , Vt ; St ) − ηVt (Vt Vt − Ut Ut ) .
2
Here α is the model parameter that characterizes the corruption fraction, γ and η are algorithmic
tunning parameters, which we specify in our analysis. Essentially, the above algorithm corresponds
to applying projected gradient method to optimize (8), where S is replaced by the aforementioned
sparse estimator in each step.
Algorithm 2 Fast RPCA with partial observations
INPUT: Observed matrix Y with support Φ; parameters τ, γ, η; number of iterations T .
// Phase I: Initialization.
1: Sinit ← T2pα [ΠΦ (Y )]
2: [L, Σ, R] ← SVDr [ p1 (Y − Sinit )]
3: U0 ← LΣ1/2 , V0 ← RΣ1/2 . Let U, V be defined according to (7).
// Phase II: Gradient based iterations.
4: U0 ← ΠU (U0 ), V0 ← ΠV (V0 )
5: for t = 0, 1, . . . , T − 1 do


>
6:
St ← Tγpα Π
Y
−
U
V
t
Φ
t


e t , Vt ; St ) − 1 ηUt (U > Ut − V > Vt )
7:
Ut+1 ← ΠU Ut − η∇U L(U
t
t
16


1
e
8:
Vt+1 ← ΠV Vt − η∇V L(Ut , Vt ; St ) − 16 ηVt (Vt> Vt − Ut> Ut )
9: end for
OUTPUT: (UT , VT )
6

4

Main Results

In this section, we establish theoretical guarantees for Algorithm 1 in the fully observed setting and
for Algorithm 2 in the partially observed setting.

4.1

Analysis of Algorithm 1

We begin with some definitions and notation. It is important to define a proper error metric because
the optimal solution corresponds to a manifold and there are many distinguished pairs (U, V ) that
minimize (8). Given the SVD of the true low-rank matrix M ∗ = L∗ Σ∗ R∗> , we let U ∗ := L∗ Σ∗1/2 and
V ∗ := R∗ Σ∗1/2 . We also let σ1∗ ≥ σ2∗ ≥ . . . ≥ σr∗ be sorted nonzero singular values of M ∗ , and denote
the condition number of M ∗ by κ, i.e., κ := σ1∗ /σr∗ . We define estimation error d(U, V ; U ∗ , V ∗ ) as
the minimal Frobenius norm between (U, V ) and (U ∗ , V ∗ ) with respect to the optimal rotation,
namely
p
d(U, V ; U ∗ , V ∗ ) := min |||U − U ∗ Q|||2F + |||V − V ∗ Q|||2F ,
(10)
Q∈Qr

for Qr the set of r-by-r orthonormal matrices. This metric controls reconstruction error, as
p
|||U V > − M ∗ |||F . σ1∗ d(U, V ; U ∗ , V ∗ ),
(11)
p
when d(U, V ; U ∗ , V ∗ ) ≤ σ1∗ . We denote the local region around the optimum (U ∗ , V ∗ ) with radius
ω as
o
n
B2 (ω) := (U, V ) ∈ Rd1 ×r × Rd2 ×r d(U, V ; U ∗ , V ∗ ) ≤ ω .
The next two theorems provide guarantees for the initialization phase and gradient iterations,
respectively, of Algorithm 1. The proofs are given in Sections 6.1 and 6.2.
Theorem 1 (Initialization). Consider the paired (U0 , V0 ) produced in the first phase of Algorithm
1. If α ≤ 1/(16κµr), we have
√
√ p
d(U0 , V0 ; U ∗ , V ∗ ) ≤ 28 καµr r σ1∗ .
Theorem 2 (Convergence). Consider the second phase of Algorithm 1. Suppose we choose γ = 2
and η = c/σ1∗ for any c ≤ 1/36.
There exist constants c1 , c2 such that when α ≤ c1 /(κ2 µr), given

p
any (U0 , V0 ) ∈ B2 c2 σr∗ /κ , the iterates {(Ut , Vt )}∞
t=0 satisfy

c t 2
d2 (Ut , Vt ; U ∗ , V ∗ ) ≤ 1 −
d (U0 , V0 ; U ∗ , V ∗ ).
8κ
Therefore, using proper initialization and step size, the gradient iteration converges at a linear
rate with a constant contraction factor 1 − O(1/κ). To obtain relative precision ε compared to
the initial error, it suffices to perform
√ O(κ log(1/ε)) iterations. Note that the step size is chosen
according to 1/σ1∗ . When α . 1/(µ κr3 ), Theorem 1 and the inequality (11) together imply that
|||U0 V0> − M ∗ |||op ≤ 12 σ1∗ . Hence we can set the step size as η = O(1/σ1 (U0 V0> )) using being the top
singular value σ1 (U0 V0> ) of the matrix U0 V0>
Combining Theorems 1 and 2 implies the following result, proved in Section 6.3, that provides
an overall guarantee for Algorithm 1.

7

Corollary 1. Suppose that
(
α ≤ c min

1
1
√ 3, 2
µ κr µκ r

)

for some constant c. Then for any ε ∈ (0, 1), Algorithm 1 with T = O(κ log(1/ε)) outputs a pair
(UT , VT ) that satisfies
|||UT VT> − M ∗ |||F ≤ ε · σr∗ .
(12)
Remark 1 (Time Complexity). For simplicity we assume d1 = d2 = d. Our sparse estimator (4)
can be implemented by finding the top αd elements of each row and column via partial quick sort,
which has running time O(d2 log(αd)). Performing rank-r SVD in the first phase and computing the
gradient in each iteration both have complexity O(rd2 ).5 Algorithm 1 thus has total running time
O(κrd2 log(1/ε)) for achieving an  accuracy as in (12). We note that when κ = O(1), our algorithm
is orderwise faster than the AltProj algorithm in [22], which has running time O(r2 d2 log(1/ε)).
Moreover, our algorithm only requires computing one singular value decomposition.
Remark 2 (Robustness). Assuming κ = O(1), our algorithm can tolerate corruption at a sparsity
√
√
level up to α = O(1/(µr r)). This is worse by a factor r compared to the optimal statistical
guarantee 1/(µr) obtained in [12, 19, 22]. This looseness is a consequence of the condition for (U0 , V0 )
in Theorem 2. Nevertheless, when µr = O(1), our algorithm can tolerate a constant α fraction of
corruptions.
Notably, we show that gradient descent works in the case of α = O(1/(µr)) if initialization is
sufficiently close. Accordingly, to provide an algorithm with optimal robustness, it is straightforward
to use a more complicated initial method such as AltProj that can tolerate 1/(µr) fraction of
corruptions while satisfying our initial condition. As gradient descent provides a simple and efficient
way to successively refine the estimation, such combination still gives a better running time than
prior arts.

4.2

Analysis of Algorithm 2

We now move to the guarantees of Algorithm 2. We show here that not only can we handle partial
observations, but in fact subsampling the data in the fully observed case can significantly reduce the
time complexity from the guarantees given in the previous section without sacrificing robustness.
In particular, for smaller values of r, the complexity of Algorithm 2 has near linear dependence on
the dimension d, instead of quadratic.
In the following discussion, we let d := max{d1 , d2 }. The next two results, proved in Sections
6.4 and 6.5, control the quality of the initialization step, and then the gradient iterations.
Theorem 3 (Initialization, partial observations). Suppose the observed indices Φ follow the Bernoulli
model given in (2). Consider the pair (U0 , V0 ) produced in the first phase of Algorithm 2. There
√
exist constants {ci }3i=1 such that for any  ∈ (0, r/(8c1 κ)), if
 2

1
µr
1
log d
α≤
, p ≥ c2
+
,
(13)
2
64κµr

α d1 ∧ d2
then we have
5

p
√
√ p
d(U0 , V0 ; U ∗ , V ∗ ) ≤ 51 καµr r σ1∗ + 7c1  κσ1∗ ,

In fact, it suffices to compute the best rank-r approximation with running time independent of the eigen gap.

8

with probability at least 1 − c3 d−1 .
Theorem 4 (Convergence, partial observations). Suppose the observed indices Φ follow the Bernoulli
model given in (2). Consider the second phase of Algorithm 2. Suppose we choose γ = 3, and
η = c/(µrσ1∗ ) for a sufficiently small constant c. There exist constants {ci }4i=1 such that if
α≤

κ4 µ2 r2 log d
c1
,
and
p
≥
c
2
κ2 µr
d1 ∧ d2

(14)

then with probability at least 1 − c3 d−1 , the iterates {(Ut , Vt )}∞
t=0 satisfy
2

∗

∗

d (Ut , Vt ; U , V ) ≤



c
1−
64µrκ

t

d2 (U0 , V0 ; U ∗ , V ∗ )


 p
for all (U0 , V0 ) ∈ B2 c4 σr∗ /κ .
The above result ensures linear convergence to (U ∗ , V ∗ ) (up to rotation) even when the gradient
iterations are computed using partial observations. Note that setting p = 1 recovers Theorem 2 up
to an additional factor µr in the contraction factor. For achieving ε relative accuracy, now we need
O(µrκ log(1/ε)) iterations.
Putting Theorems 3 and 4 together, we have the following overall guarantee, proved in Section
6.6, for Algorithm 2.
Corollary 2. Suppose that
(
α ≤ c min

1
1
√ 3, 2
µ κr µκ r

)
, p ≥ c0

κ4 µ2 r2 log d
,
d1 ∧ d2

for some constants c, c0 . With probability at least 1 − O(d−1 ), for any ε ∈ (0, 1), Algorithm 2 with
T = O(µrκ log(1/ε)) outputs a pair (UT , VT ) that satisfies
|||UT VT> − M ∗ |||F ≤ ε · σr∗ .

(15)

This result shows that partial observations do not compromise robustness to sparse corruptions:
as long as the observation probability p satisfies the condition in Corollary 2, Algorithm 2 enjoys
the same robustness guarantees as the method using all entries. Below we provide two remarks on
the sample and time complexity. For simplicity, we assume d1 = d2 = d, κ = O(1).
Remark 3 (Sample complexity and matrix completion). Using the lower bound on p, it is sufficient
to have O(µ2 r2 d log d) observed entries. In the special case S ∗ = 0, our partial observation model
is equivalent to the model of exact matrix completion (see, e.g., [7]). We note that our sample
complexity (i.e., observations needed) matches that of completing a positive semidefinite (PSD)
matrix by gradient descent as shown in [11], and is better than the non-convex matrix completion
algorithms in [20] and [25]. Accordingly, our result reveals the important fact that we can obtain
robustness in matrix completion without deterioration of our statistical guarantees. It is known that
that any algorithm for solving exact matrix completion must have sample size Ω(µrd log d) [7], and a
nearly tight upper bound O(µrd log2 d) is obtained in [10] by convex relaxation. While sub-optimal
by a factor µr, our algorithm is much faster than convex relaxation as shown below.
9

d(U, V ; U ∗ , V ∗ )

GD p = 1
GD p = 0.5
GD p = 0.2

10 3

Time(secs)

d(U, V ; U ∗ , V ∗ )

10 1
10 0

10 -1

10 2

10 1

10 -2

10

GD p = 1
GD p = 0.5
GD p = 0.2
GD p = 0.1
AltProj
IALM

0

10 -1

10 -2

10 -3
0

2

4

6

8

10

Iteration count

(a)

10 3

10 4

10 5

Dimension d

(b)

0

20

40

60

80

100

Time(secs)

(c)

Figure 1. Results on synthetic data. (a) Plot of log estimation error versus number of iterations
when using gradient descent (GD) with varying sub-sampling rate p. It is conducted using d =
5000, r = 10, α = 0.1. (b) Plot of running time of GD versus dimension d with r = 10, α = 0.1, p =
0.15r2 log d/d. The low-rank matrix is recovered in all instances, and the line has slope approximately
one. (c) Plot of log estimation error versus running time for different algorithms in problem with
d = 5000, r = 10, α = 0.1.

Remark 4 (Time complexity). Our sparse estimator on the sparse matrix with support Φ can be
implemented via partial quick sort with running time O(pd2 log(αpd)). Computing the gradient
in each step involves the two terms in the objective function (9). Computing the gradient of the
first term Le takes time O(r|Φ|), whereas the second term takes time O(r2 d). In the initialization
phase, performing rank-r SVD on a sparse matrix with support Φ can be done in time O(r|Φ|). We
conclude that when |Φ| = O(µ2 r2 d log d), Algorithm 2 achieves the error bound (15) with running
time O(µ3 r4 d log d log(1/ε)). Therefore, in the small rank setting with r  d1/3 , even when full
observations are given, it is better to use Algorithm 2 by subsampling the entries of Y .

5

Numerical Results

In this section, we provide numerical results and compare the proposed algorithms with existing
methods, including the inexact augmented lagrange multiplier (IALM) approach [21] for solving
the convex relaxation (1) and the alternating projection (AltProj) algorithm proposed in [22]. All
algorithms are implemented in MATLAB 6 , and the codes for existing algorithms are obtained
from their authors. SVD computation in all algorithms uses the PROPACK library.7 We ran all
simulations on a machine with Intel 32-core Xeon (E5-2699) 2.3GHz with 240GB RAM.

5.1

Synthetic Datasets

We generate a squared data matrix Y = M ∗ + S ∗ ∈ Rd×d as follows. The low-rank part M ∗ is
given by M ∗ = AB > , where A, B ∈ Rd×r have entries drawn independently from a zero mean
Gaussian distribution with variance 1/d. For a given sparsity parameter α, each entry of S ∗ is set
to be nonzero with probability α, and the values of the nonzero entries are sampled uniformly from
[−5r/d, 5r/d].
The results are summarized in Figure 1. Figure 1a shows the convergence of our algorithms
for different random instances with different sub-sampling rate p (note that p = 1 corresponds to
6
7

Our code is available at https://www.yixinyang.org/code/RPCA_GD.zip.
http://sun.stanford.edu/~rmunk/PROPACK/

10

Original

GD (49.8s)

GD, 20% sample (18.1s)

AltProj (101.5s)

IALM (434.6s)

Original

GD (87.3s)

GD, 20% sample (43.4s)

AltProj (283.0s)

IALM (801.4s)

Figure 2. Foreground-background separation in Restaurant and ShoppingMall videos. In each line,
the leftmost image is an original frame, and the other four are the separated background obtained
from our algorithms with p = 1, p = 0.2, AltProj, and IALM. The running time required by each
algorithm is shown in the title.

the fully observed setting). As predicted by Theorems 2 and 4, our gradient method converges
geometrically with a contraction factor nearly independent of p. Figure 1b shows the running time
of our algorithm with partially observed data. We see that the running time scales linearly with d,
again consistent with the theory. We note that our algorithm is memory-efficient: in the large scale
setting with d = 2×105 , using approximately 0.1% entries is sufficient for the successful recovery. In
contrast, AltProj and IALM are designed to manipulate the entire matrix with d2 = 4×1010 entries,
which is prohibitive on a single machine. Figure 1c compares our algorithms with AltProj and IALM
by showing reconstruction error versus real running time. Our algorithm requires significantly less
computation to achieve the same accuracy level, and using only a subset of the entries provides
additional speed-up.

5.2

Foreground-background Separation

We apply our method to the task of foreground-background (FB) separation in a video. We use
two public benchmarks, the Restaurant and ShoppingMall datasets.8 Each dataset contains a video
with static background. By vectorizing and stacking the frames as columns of a matrix Y , the FB
separation problem can be cast as RPCA, where the static background corresponds to a low rank
matrix M ∗ with identical columns, and the moving objects in the video can be modeled as sparse
corruptions S ∗ . Figure 2 shows the output of different algorithms on two frames from the dataset.
Our algorithms require significantly less running time than both AltProj and IALM. Moreover, even
with 20% sub-sampling, our methods still appear to achieve better separation quality (note that in
each of the frames our algorithms remove a person that is not identified by the other algorithms).
Figure 3 shows recovery results for several more frames. Again, our algorithms enjoy better
running time and outperform AltProj and IALM in separating persons from the background images.
In Appendix B, we describe the detailed parameter settings for our algorithm.
8

http://perception.i2r.a-star.edu.sg/bk_model/bk_index.html

11

Original

GD (49.8s)

GD, 20% sample (18.1s)

AltProj (101.5s)

IALM (434.6s)

Original

GD (49.8s)

GD, 20% sample (18.1s)

AltProj (101.5s)

IALM (434.6s)

Original

GD (87.3s)

GD, 20% sample (43.4s)

AltProj (283.0s)

IALM (801.4s)

Original

GD (87.3s)

GD, 20% sample (43.4s)

AltProj (283.0s)

IALM (801.4s)

Figure 3. More results of FB separation in Restaurant and ShoppingMall videos. The leftmost
images are original frames. The right four columns show the results from our algorithms with p =
1, p = 0.2, AltProj [22], and IALM [21]. The runtime of each algorithm is written in the title.

12

6

Proofs

In this section we provide the proofs for our main theoretical results in Theorems 1–4 and Corollaries 1–2.

6.1

Proof of Theorem 1

Let Y := Y − Sinit . As Y = M ∗ + S ∗ , we have Y − M ∗ = S ∗ − Sinit . We obtain Y − M ∗ ∈ S2α
because S ∗ , Sinit ∈ Sα .
We claim that kY −M ∗ k∞ ≤ 2kM ∗ k∞ . Denote the support of S ∗ , Sinit by Ω∗ and Ω respectively.
Since Y − M ∗ is supported on Ω ∪ Ω∗ , to prove the claim it suffices to consider the following three
cases.
• For (i, j) ∈ Ω∗ ∩ Ω, due to rule of sparse estimation, we have (S ∗ − Sinit )(i,j) = 0.
∗ | ≤ 2kM ∗ k . Otherwise, we have |Y
∗
• For (i, j) ∈ Ω∗ \ Ω, we must have |S(i,j)
∞
(i,j) | = |(S +
M ∗ )(i,j) | > kM ∗ k∞ . So |Y(i,j) | is larger than any uncorrupted entries in its row and column.
Since there are at most α fraction corruptions per row and column, we have Y(i,j) ∈ Ω, which
violates the prior condition (i, j) ∈ Ω∗ \ Ω.
∗ , trivially we have |(S
• For the last case (i, j) ∈ Ω \ Ω∗ , since (Sinit )(i,j) = M(i,j)
init )(i,j) | ≤
∗
kM k∞ .

The following result, proved in Section 7.1, relates the operator norm of Y − M ∗ to its infinite
norm.
Lemma 1. For any matrix A ∈ Rd1 ×d2 that belongs to Sα given in (3), we have
p
|||A|||op ≤ α d1 d2 kAk∞ .
We thus obtain
|||Y − M ∗ |||op ≤ 2α

p
p
d1 d2 kY − M ∗ k∞ ≤ 4α d1 d2 kM ∗ k∞ = 4αµrσ1∗ .

In the last step, we use the fact that M ∗ satisfies the µ-incoherent condition, which leads to
µr
kM ∗ k∞ ≤ |||M ∗ |||op |||L∗ |||2,∞ |||R∗ |||2,∞ ≤ √
|||M ∗ |||op .
d1 d2

(16)

(17)

We denote the i-th largest singular value of Y by σi . By Weyl’s theorem, we have |σi∗ − σi | ≤
∗
|||Y − M ∗ |||op for all i ∈ [d1 ∧ d2 ]. Since σr+1
= 0, we have σr+1 ≤ |||Y − M ∗ |||op . Recall that U0 V0> is
the best rank r approximation of Y . Accordingly, we have
|||U0 V0> − M ∗ |||op ≤ |||U0 V0> − Y |||op + |||Y − M ∗ |||op
= σr+1 + |||Y − M ∗ |||op ≤ 2|||Y − M ∗ |||op ≤ 8αµrσ1∗ .
1
Under condition αµr ≤ 16κ
, we obtain |||U0 V0> − M ∗ |||op ≤ 21 σr∗ . Applying Lemma 5.14 in [26]
(we provide it as Lemma 15 for the sake of completeness), we obtain

d2 (U0 , V0 ; U ∗ , V ∗ ) ≤ √

10r|||U0 V0> − M ∗ |||2op
2 |||U0 V0> − M ∗ |||2F
≤
.
σr∗
σr∗
2−1

Plugging the upper bound of |||U0 V0> − M ∗ |||op into the above inequality completes the proof.
13

6.2

Proof of Theorem 2

We essentially follow the general framework developed in [11] for analyzing the behaviors of gradient
descent in factorized low-rank optimization. But it is worth to note that [11] only studies the
symmetric and positive semidefinite setting, while we avoid such constraint on M ∗ . The techniques
for analyzing general asymmetric matrix in factorized space is inspired by the recent work [26] on
solving low-rank matrix equations. In our setting, the technical challenge is to verify the local
descent condition of the loss function (8), which not only has a bilinear dependence on U and V ,
but also involves our sparse estimator (4).
We begin with some notations. Define the equivalent set of optimal solution as
o
n
(18)
E(M ∗ ) := (A, B) ∈ Rd1 ×r × Rd2 ×r A = L∗ Σ∗1/2 Q, B = R∗ Σ∗1/2 Q, where Q ∈ Qr .
 p

Given (U0 , V0 ) ∈ B2 c2 σr∗ /κ , by (11), we have |||U0 V0> − M ∗ |||op ≤ 21 σr∗ when c2 is sufficiently
small. By Weyl’s theorem We thus have
q
q
q
q
∗
∗
∗
σ1 /2 ≤ |||U0 |||op ≤ 3σ1 /2, and σ1 /2 ≤ |||V0 |||op ≤ 3σ1∗ /2.
As a result, for U, V constructed according to (7), we have
E(M ∗ ) ⊆ U × V, and U ⊆ U, V ⊆ V,

(19)

where
(
U :=

A ∈ Rd1 ×r

r
|||A|||2,∞ ≤

3µrσ1∗
d1

(

)
, V :=

A ∈ Rd2 ×r

r
|||A|||2,∞ ≤

3µrσ1∗
d2

)
.

We let

1
(20)
G(U, V ) := |||U > U − V > V |||2F .
8
For L(U, V ; S), we denote the gradient with respect to M by ∇M L(U, V ; S), i.e. ∇M L(U, V ; S) =
UV > + S − Y .
The local descent property is implied by combining the following two results, which are proved
in Section 6.7 and 6.8 respectively.
Lemma
2 (Local descent property of
 L). Suppose U, V satisfy (19). For any (U, V ) ∈ (U × V) ∩
p
B2 ( σ1∗ ), we let S = Tγα Y − U V > , where we choose γ = 2. Then we have that for (Uπ∗ , Vπ∗ ) ∈
argmin(A,B)∈E(M ∗ ) |||U − A|||2F + |||V − B|||2F and β > 0,
>
> 2
∗
hh∇M L(U, V ; S), U V > − Uπ∗ Vπ>∗ + ∆U ∆>
V ii ≥ |||U V − Uπ ∗ Vπ ∗ |||F − νσ1 δ − 3

q
σ1∗ δ 3 .

Here ∆U := U − Uπ∗ , ∆V := V − Vπ∗ , δ := |||∆U |||2F + |||∆V |||2F , and ν := 9(β + 6)αµr + 5β −1 .
√
Lemma 3 (Local descent property of G). For any (U, V ) ∈ B2 ( σr∗ ) and
(Uπ∗ , Vπ∗ ) ∈ arg

min

(A,B)∈E(M ∗ )

14

|||U − A|||2F + |||V − B|||2F ,

(21)

we have
hh∇U G(U, V ), U − Uπ∗ ii + hh∇V G(U, V ), V − Vπ∗ ii
r
σ1∗ δ 3 1
1 >
1 ∗
>
2
≥ |||U U − V V |||F + σr δ −
− |||U V > − Uπ∗ Vπ>∗ |||2F ,
8
8
2
2
where δ is defined according to Lemma 2.
As another key ingredient, we establish the following smoothness condition, proved in Section
6.9, which indicates that the Frobenius norm of gradient decreases as (U, V ) approaches the optimal
manifold.


p
Lemma 4 (Smoothness). For any (U, V ) ∈ B2 ( σ1∗ ), we let S = Tγα Y − U V > , where we choose
γ = 2. We have that
|||∇M L(U, V ; S)|||2F ≤ 6|||U V > − M ∗ |||2F ,
(22)
and
|||∇U G(U, V )|||2F + |||∇V G(U, V )|||2F ≤ 2σ1∗ |||U > U − V > V |||2F .

(23)

With the above results in hand, we are ready to prove Theorem 2.
Proof of Theorem 2. We use shorthands
δt := d2 (Ut , Vt ; U ∗ , V ∗ ), Lt := L(Ut , Vt ; St ), and Gt := G(Ut , Vt ).
For (Ut , Vt ), let (Uπt ∗ , Vπt∗ ) := argmin(A,B)∈E(M ∗ ) |||Ut − A|||2F + |||Vt − B|||2F . Define ∆tU := Ut − Uπt ∗ ,
∆tV := Vt − Vπt∗ .
We prove Theorem 2 by induction. It is sufficient
 topconsider
 one step of the iteration. For any
t ≥ 0, under the induction hypothesis (Ut , Vt ) ∈ B2 c2 σr∗ /κ . We find that
δt+1 ≤ |||Ut+1 − Uπt ∗ |||2F + |||Vt+1 − Vπt∗ |||2F
≤ |||Ut − η∇U Lt − η∇U Gt − Uπt ∗ |||2F + |||Vt − η∇V Lt − η∇V Gt − Vπt∗ |||2F
≤ δt − 2η hh∇U Lt + ∇U Gt , Ut − Uπt ∗ ii −2η hh∇V Lt + ∇V Gt , Vt − Vπt∗ ii
|
{z
}
|
{z
}
W1

2

W2

2

2

2

+ η |||∇U Lt + ∇U Gt |||F +η |||∇V Lt + ∇V Gt |||F ,
|
|
{z
}
{z
}
W3

(24)

W4

where the second step follows from the non-expansion property of projection onto U, V, which is
implied by E(M ∗ ) ⊆ U × V shown in (19). Since ∇U Lt = [∇M Lt ] V and ∇V Lt = [∇M Lt ]> U , we
have
t
t>
hh∇U Lt , Ut − Uπt ∗ ii + hh∇V Lt , Vt − Vπt∗ ii = hh∇M Lt , Ut Vt> − Uπt ∗ Vπt>
∗ + ∆U ∆V ii.

Combining Lemma 2 and 3, under condition δt < σr∗ , we have that
1
1
1
W1 + W2 ≥ |||Ut Vt> − M ∗ |||2F + |||Ut> Ut − Vt> Vt |||2F + σr∗ δt − νσ1∗ δt − 4
2
8
8

15

q
σ1∗ δt3 .

On the other hand, we have
W3 + W4 ≤ 2|||∇U Lt |||2F + 2|||∇U Gt |||2F + 2|||∇V Lt |||2F + 2|||∇V Gt |||2F
≤ 2(|||Ut |||2op + |||Vt |||2op )|||∇M Lt |||2F + 2|||∇U Gt |||2F + 2|||∇V Gt |||2F
≤ 36σ1∗ |||Ut Vt> − M ∗ |||2F + 4σ1∗ |||Ut> Ut − Vt> Vt |||2F ,

 p
where the last step is implied by Lemma 4 and the assumption (Ut , Vt ) ∈ B2 c2 σr∗ /κ that leads
p
p
to |||Ut |||op ≤ 3σ1∗ /2, |||Vt |||op ≤ 3σ1∗ /2.
By the assumption η = c/σ1∗ for any constant c ≤ 1/36, we thus have
q
1 ∗
2
∗
−2η(W1 + W2 ) + η (W3 + W4 ) ≤ − ησr δt + 2ηνσ1 δt + 8η σ1∗ δt3 .
4
In Lemma 2, choosingpβ = 320κ and assuming α . 1/(κ2 µr), we can have ν ≤ 1/(32κ). Assuming
1 ∗
δt . σr∗ /κ leads to 14 σ1∗ δt3 ≤ 16
σr δt . We thus obtain


ησr∗
δt+1 ≤ 1 −
δt .
(25)
8
Under initial condition δ0 . σr∗ /κ, we obtain that such condition holds for all t since estimation
error decays geometrically after each iteration. Then applying (25) for all iterations, we conclude
that for all t = 0, 1, . . .,


ησ ∗ t
δt ≤ 1 − r δ0 .
8

6.3

Proof of Corollary 1

We need α . κ21µr due to the condition of Theorem 2. In order to ensure the linear convergence
happens, it suffices to let the initial error shown in Theorem 1 be less than the corresponding
condition in Theorem 2. Accordingly, we need
p
√
√ p
28 καµr r σ1∗ . σr∗ /κ,
which leads to α .

1
√ 3.
µ rκ

Using the conclusion that gradient descent has linear convergence, choosing T = O(κ log(1/ε)),
we have
σ∗
d2 (UT , VT ; U ∗ , V ∗ ) ≤ ε2 d2 (U0 , V0 ; U ∗ , V ∗ ) . ε2 r .
κ
∗
∗
>
Finally, applying the relationship between d(UT , VT ; U , V ) and |||UT VT − M ∗ |||F shown in (11), we
complete the proof.

6.4

Proof of Theorem 3

Let Y := p1 (Y − Sinit ). Similar to the proof of Theorem 1, we first establish an upper bound on
|||Y − M ∗ |||op . We have that
1
1
|||Y − M ∗ |||op ≤ |||Y − ΠΦ M ∗ |||op + ||| ΠΦ (M ∗ ) − M ∗ |||op .
p
p
16

(26)

For the first term, we have Y − p1 ΠΦ M ∗ = p1 (ΠΦ (S ∗ ) − Sinit ) because Y = ΠΦ (M ∗ + S ∗ ). Lemma
d
10 shows that under condition p & α(dlog
, there are at most 32 pα-fraction nonzero entries in each
1 ∧d2 )
row and column of ΠΦ (S ∗ ) with high probability. Since Sinit ∈ S2pα , we have

ΠΦ (S ∗ ) − Sinit ∈ S4pα .

(27)

kΠΦ (S ∗ ) − Sinit k∞ ≤ 2kM ∗ k∞ .

(28)

In addition, we prove below that

Denote the support of ΠΦ (S ∗ ) and Sinit by Ω∗o and Ω. For (i, j) ∈ Ω∗o ∩ Ω and (i, j) ∈ Ω \ Ω∗o , we
∗ , respectively. To prove the claim,
have (ΠΦ (S ∗ ) − Sinit )(i,j) = 0 and (ΠΦ (S ∗ ) − Sinit )(i,j) = −M(i,j)
∗ | < 2kM ∗ k . If this is not true, then we must have
it remains to show that for (i, j) ∈ Ω∗o \ Ω, |S(i,j)
∞
∗
|Y(i,j) | > kM k∞ . Accordingly, |Y(i,j) | is larger than the magnitude of any uncorrupted entries in
its row and column. Note that on the support Φ, there are at most 32 pα corruptions per row and
column, we have (i, j) ∈ Ω, which violates our prior condition (i, j) ∈ Ω∗o \ Ω.
Using these two properties (27), (28) and applying Lemma 1, we have
p
p
1
|||Y − ΠΦ M ∗ |||op ≤ 4α d1 d2 kΠΦ (S ∗ ) − Sinit k∞ ≤ 8α d1 d2 kM ∗ k∞ ≤ 8αµrσ1∗ ,
(29)
p
where the last step follow from (17).
For the second term in (26), we use the following lemma proved in [10].
Lemma 5 (Lemma 2 in [10]). Suppose A ∈ Rd1 ×d2 is a fixed matrix. We let d := max{d1 , d2 }.
There exists a constant c such that with probability at least 1 − O(d−1 ),
s
!
n
o
1
log d
log d
>
||| ΠΦ (A) − A|||op ≤ c
kAk∞ +
max |||A|||2,∞ , |||A |||2,∞
.
p
p
p
Given the SVD M ∗ = L∗ ΣR∗> , for any i ∈ [d1 ], we have
∗
kM(i,·)
k2

=

kL∗(i,·) ΣR∗> k2

≤

σ1∗ kL∗(i,·) k2

≤

σ1∗

r

µr
.
d1

We can bound |||M ∗> |||2,∞ similarly. Lemma 5 leads to
1
log d µr
√
||| ΠΦ (M ∗ ) − M ∗ |||op ≤ c0
σ ∗ + c0
p
p
d1 d2 1

s

log d
p

r

√
µr
σ1∗ ≤ c0 σ1∗ / r
d1 ∧ d2

(30)

2

log d
under condition p ≥ 4µr
2 (d ∧d ) .
1
2
Putting (29) and (30) together, we obtain

√
|||Y − M ∗ |||op ≤ 8αµrσ1∗ + c0 σ1∗ / r.
Then using the fact that U0 V0> is the best rank r approximation of Y and applying Wely’s theorem
(see the proof of Theorem 1 for a detailed argument), we have
|||U0 V0> − M ∗ |||op ≤ |||U0 V0> − Y |||op + |||Y − M ∗ |||op

√
≤ 2|||Y − M ∗ |||op ≤ 16αµrσ1∗ + 2c0 σ1∗ / r
17

√
Under our assumptions, we have 16αµrσ1∗ + 2c0 σ1∗ / r ≤ 12 σr∗ . Accordingly, Lemma 15 gives
d2 (U0 , V0 ; U ∗ , V ∗ ) ≤ √

10r|||U0 V0> − M ∗ |||2op
2 |||U0 V0> − M ∗ |||2F
≤
.
σr∗
σr∗
2−1

We complete the proof by combining the above two inequalities.

6.5

Proof of Theorem 4

In this section, we turn to prove Theorem 4. Similar to the proof of Theorem 2, we rely on
establishing the local descent and smoothness conditions. Compared to the full observation setting,
e V ) := 1 |||U > U − V > V |||2 merely
we replace L by Le given in (6), while the regularization term G(U,
F
64
differs from G(U, V ) given in (20) by a constant factor. It is thus sufficient to analyze the properties
e
of L.
Define E(M ∗ ) according to (18). Under the initial condition, we still have
E(M ∗ ) ⊆ U × V, and U ⊆ U, V ⊆ V.

(31)

We prove the next two lemmas in Section 6.10 and 6.11 respectively. In both lemmas, for any
(U, V ) ∈ U × V, we use shorthands
(Uπ∗ , Vπ∗ ) = arg

min

(A,B)∈E(M ∗ )

|||U − A|||2F + |||V − B|||2F ,

∆U := U − Uπ∗ , ∆V := V − Vπ∗ , and δ := |||∆U |||2F + |||∆V |||2F . Recall that d := max{d1 , d2 }.
e Suppose U, V satisfy (31). Suppose we let
Lemma 6 (Local descent property of L).
i
h

S = Tγpα ΠΦ Y − U V > ,
√
where we choose γ = 3. For any β > 0 and  ∈ (0, 41 ), we define ν := (14β +81)αµr +26 +18β −1 .
There exist constants {ci }2i=1 such that if
 2 2

µ r
1
log d
p ≥ c1
+ +1
,
(32)
2

α
d1 ∧ d2
then with probability at least 1 − c2 d−1 ,
q
3
>
> 2
∗
e V ; S), U V > − Uπ∗ Vπ>∗ + ∆U ∆>
∗ V ∗ ||| − νσ δ − 10
hh∇M L(U,
ii
≥
|||U
V
−
U
σ1∗ δ 3 − 2δ 2 (33)
π
V
π F
1
16
p ∗
for all (U, V ) ∈ (U × V) ∩ B2
σ1 .


e Suppose U, V satisfy (31). Suppose we let S = Tγαp ΠΦ Y − U V >
Lemma 7 (Smoothness of L).
for γ = 3. There exist constants {ci }3i=1 such that for any  ∈ (0, 41 ), when p satisfies
condition
p
(32), with probability at least 1 − c2 d−1 , we have that for all (U, V ) ∈ (U × V) ∩ B2 ( σ1∗ ),
i
h
e V ; S)|||2 + |||∇V L(U,
e V ; S)|||2 ≤ c3 µrσ ∗ |||U V > − Uπ∗ V >∗ |||2 + µrσ ∗ δ(δ + σ ∗ ) .
(34)
|||∇U L(U,
F
1
π F
1
1
F

18

In the remainder of this section, we condition on the events in Lemma 6 and 7. Now we are
ready to prove Theorem 4.
Proof of Theorem 4. We essentially follow the process for proving Theorem 2. Let the following
shorthands be defined in the same fashion: δt , (Uπt ∗ , Vπt∗ ), (∆tU , ∆tV ), Let , Get .
Here we show error decays in one step of iteration. The induction process is the same as the
proof of Theorem 2, and is thus omitted. For any t ≥ 0, similar to (24) we have that
δt+1 ≤ δt − 2η hh∇U Let + ∇U Get , Ut − Uπt ∗ ii −2η hh∇V Let + ∇V Get , Vt − Vπt∗ ii
{z
}
|
{z
}
|
W1

W2

+ η |||∇U Let + ∇U Get |||2F +η 2 |||∇V Let + ∇V Get |||2F .
|
|
{z
}
{z
}
2

W3

W4

We also have
t
t>
hh∇U Let , Ut − Uπt ∗ ii + hh∇V Let , Vt − Vπt∗ ii = hh∇M Let , Ut Vt> − Uπt ∗ Vπt>
∗ + ∆U ∆V ii,

which can be lower bounded by Lemma 6. Note that Ge differs from G by a constant, we can still
leverage Lemma 3. Hence, we obtain that
q
1
1
1
W1 + W2 ≥ |||Ut Vt> − M ∗ |||2F + |||Ut> Ut − Vt> Vt |||2F + σr∗ δt − νσ1∗ δt − 11 σ1∗ δt3 − 2δt2 .
8
64
64
On the other hand, we have
W3 + W4 ≤ 2|||∇U Let |||2F + 2|||∇U Get |||2F + 2|||∇V Let |||2F + 2|||∇V Get |||2F
h
i
≤ c µrσ1∗ |||Ut Vt> − M ∗ |||2F + µrσ1∗ δt (δt + σ1∗ ) + σ1∗ |||Ut> Ut − Vt> Vt |||2F ,
where c is a constant, and the last step is implied by Lemma 4 and Lemma 7.
By the assumption η = c0 /[µrσ1∗ ] for sufficiently small constant c0 , we thus have
q
1
−2η(W1 + W2 ) + η 2 (W3 + W4 ) ≤ − ησr∗ δt + 2ηνσ1∗ δt + 22η σ1∗ δt3 + 4ηδt2 .
32
√
Recall that ν := (14β + 81)αµr + 26  + 18β −1 . By letting β = c1 κ,  = c2 /κ2 and assuming
α ≤ c3 /(µrκ2 ) and δt ≤ c4 σr∗ /κ for some sufficiently small constants {ci }4i=1 , we can have −2η(W1 +
1
W2 ) + η 2 (W3 + W4 ) ≤ − 64
ησr∗ δt , which implies that


ησr∗
δt+1 ≤ 1 −
δt ,
64
and thus completes the proof.

6.6

Proof of Corollary 2

We need α . µκ12 r due to the condition of Theorem 4. Letting the initial error provided in Theorem
3 be less than the corresponding condition in Theorem 4, we have
p
p
√
√ p
51 καµr r σ1∗ + 7c1  κσ1∗ . σr∗ /κ,
19

which leads to

1
1
α. √
, . √ .
µ r3 κ3
κ3

Plugging the above two upper bounds into the second term in (13), it suffices to have
p&

κ3 µr2 log d
.
d1 ∧ d2

Comparing the above bound with the second term in (14) completes the proof.

6.7

Proof of Lemma 2

Let M := U V > . We observe that
∇M L(U, V ; S) = M + S − M ∗ − S ∗ .
Plugging it back into the left hand side of (21), we obtain
∗
∗
∗
>
hh∇M L(U, V ; S), U V > − Uπ∗ Vπ>∗ + ∆U ∆>
V ii = hhM + S − M − S , M − M + ∆U ∆V ii

≥ |||M − M ∗ |||2F − |hhS − S ∗ , M − M ∗ ii| − |hhM + S − M ∗ − S ∗ , ∆U ∆>
V ii| .
{z
} |
{z
}
|
T1

(35)

T2

Next we derive upper bounds of T1 and T2 respectively.
Upper bound of T1 . We denote the support of S, S ∗ by Ω and Ω∗ respectively. Since S − S ∗ is
supported on Ω∗ ∪ Ω, we have
T1 ≤ |hhΠΩ (S − S ∗ ), M − M ∗ ii| + |hhΠΩ∗ \Ω (S − S ∗ ), M − M ∗ ii| .
|
{z
} |
{z
}
W1

W2

Recall that for any (i, j) ∈ Ω, we have S(i,j) = (M ∗ + S ∗ − M )(i,j) . Accordingly, we have
W1 = |||ΠΩ (M − M ∗ )|||2F .

(36)

Now we turn to bound W2 . Since S(i,j) = 0 for any (i, j) ∈ Ω∗ \ Ω, we have
W2 = |hhΠΩ∗ \Ω S ∗ , M − M ∗ ii|.
Let ui be the i-th row of M − M ∗ , and vj be the j-th column of M − M ∗ . For any k ∈ [d2 ], we let
(k)
ui denote the element of ui that has the k-th largest magnitude. Similarly, for any k ∈ [d1 ], we
(k)
let vj denote the element of vj that has the k-th largest magnitude.
From the design of sparse estimator (4), we have that for any (i, j) ∈ Ω∗ \ Ω, |(M ∗ + S ∗ − M )(i,j) |
is either smaller than the γαd2 -th largest entry of the i-th row of M ∗ + S ∗ − M or smaller than
the γαd1 -th largest entry of the j-th column of M ∗ + S ∗ − M . Note that S ∗ only contains at most
α-fraction nonzero entries per row and column. As a result, |(M ∗ + S ∗ − M )(i,j) | has to be less than
(γαd2 −αd2 )

the magnitude of ui

(γαd1 −αd1 )

. Formally, we have for (i, j) ∈ Ω∗ \ Ω,
n
o
(γαd −αd2 )
(γαd −αd1 )
|, |vj 1
| .
|(M ∗ + S ∗ − M )(i,j) | ≤ max |ui 2
|
{z
}
or vj

bij

20

(37)

Furthermore, we obtain
(γαd2 −αd2 ) 2

b2ij ≤ |ui

(γαd1 −αd1 ) 2

| + |vj

| ≤

kvj k22
kui k22
+
.
(γ − 1)αd2 (γ − 1)αd1

(38)

Meanwhile, for any (i, j) ∈ Ω∗ \ Ω, we have
∗
· (M − M ∗ )(i,j) | = |(M ∗ + S ∗ − M − M ∗ + M )(i,j) · (M − M ∗ )(i,j) |
|S(i,j)

≤ |(M − M ∗ )(i,j) |2 + |(M ∗ + S ∗ − M )(i,j) | · |(M − M ∗ )(i,j) |
≤ |(M − M ∗ )(i,j) |2 + bij · |(M − M ∗ )(i,j) |


b2ij
β
|(M − M ∗ )(i,j) |2 +
,
≤ 1+
2
2β

(39)

where β in the last step can be any positive number. Combining (38) and (39) leads to
X
∗
· (M − M ∗ )(i,j) |
W2 ≤
|S(i,j)
(i,j)∈Ω∗ \Ω



X b2ij
β
≤ 1+
|||ΠΩ∗ \Ω (M − M ∗ )|||2F +
2
2β
(i,j)∈Ω∗ \Ω



X  kui k2
kvj k22
β
1
∗ 2
2
≤ 1+
|||ΠΩ∗ \Ω (M − M )|||F +
+
2
2β
(γ − 1)αd2 (γ − 1)αd1
∗
(i,j)∈Ω \Ω


β
1
≤ 1+
|||ΠΩ∗ \Ω (M − M ∗ )|||2F +
|||M − M ∗ |||2F .
2
β(γ − 1)

(40)

In the last step, we use


X
(i,j)∈Ω∗ \Ω

≤

1
1
kui k22 + kvj k22
d2
d1

X X
i∈[d] j∈Ω∗(i,·)

≤α

X
i∈[d]

1
kui k22 +
d2

kui k22 + α

X


≤


X 1
1
2
2
kui k2 + kvj k2
d2
d1
∗

(i,j)∈Ω

X X
j∈[d] i∈Ω∗(·,j)

1
kvj k22
d1

kvj k22 = 2α|||M − M ∗ |||2F .

(41)

j∈[d]

We introduce shorthand δ := |||∆U |||2F + |||∆V |||2F . We prove the following inequality in the end of this
section.
p
(42)
|||M − M ∗ |||F ≤ 5σ1∗ δ.
Combining (36), (40) and (42) leads to


β
5σ1∗ δ
T1 ≤ |||ΠΩ (M − M ∗ )|||2F + 1 +
|||ΠΩ∗ \Ω (M − M ∗ )|||2F +
2
β(γ − 1)
∗
5σ1 δ
≤ 9(2γ + β + 2)αµrσ1∗ δ +
,
β(γ − 1)

(43)

where the last step follows from Lemma 14 by noticing that ΠΩ (M − M ∗ ) has at most γα-fraction
nonzero entries per row and column.
21

Upper bound of T2 . To ease notation, we let C := M + S − M ∗ − S ∗ . We observe that C is
supported on Ωc , we have
>
T2 ≤ |hhΠΩ∗c ∩Ωc (M − M ∗ ), ∆U ∆>
V ii| + |hhΠΩ∗ ∩Ωc C, ∆U ∆V ii| .
|
{z
} |
{z
}
W3

W4

By Cauchy-Swartz inequality, we have
∗
W3 ≤ |||ΠΩ∗c ∩Ωc (M − M ∗ )|||F |||∆U ∆>
V |||F ≤ |||M − M |||F |||∆U |||F |||∆V |||F ≤

q
5σ1∗ δ 3 /2,

where the last step follows from (42) and |||∆U |||F |||∆V |||F ≤ δ/2.
It remains to bound W4 . By Cauchy-Swartz inequality, we have
∗
∗
>
W4 ≤ |||ΠΩ∗ ∩Ωc C|||F |||∆U ∆>
V |||F ≤ |||ΠΩ∗ ∩Ωc (M + S − M )|||F |||∆U ∆V |||F

1/2
s X
2
2
X
(a)
(b)
kvj k2 
kui k2
≤
b2ij |||∆U |||F |||∆V |||F ≤ 
+
|||∆U |||F |||∆V |||F .
(γ − 1)αd2 (γ − 1)αd1
(i,j)∈Ω∗ \Ω
(i,j)∈Ω∗ \Ω
s
r
(c)
2
5σ1∗ δ 3
∗
|||M − M |||F |||∆U |||F |||∆V |||F ≤
,
≤
γ−1
2(γ − 1)

where step (a) is from (37), step (b) follows from (38), and step (c) follows from (41). Combining
the upper bounds of W3 and W4 , we obtain
s
q
5σ1∗ δ 3
T2 ≤ 5σ1∗ δ 3 /2 +
.
(44)
2(γ − 1)
Combining pieces. Now we choose γ = 2. Then inequality (43) implies that
T1 ≤ [9(β + 6)αµr + 5β −1 ]σ1∗ δ.
Inequality (44) then implies that
q
T2 ≤ 3 σ1∗ δ 3 .
Plugging the above two inequalities into (35) completes the proof.
Proof of inequality (42). We find that
hp
i2
σ1∗ (|||∆V |||F + |||∆U |||F ) + |||∆U |||F |||∆V |||F

2
p
1p ∗
1p ∗
∗
≤
σ1 (|||∆V |||F + |||∆U |||F ) +
σ1 |||∆U |||F +
σ1 |||∆V |||F
2
2

|||M − M ∗ |||2F ≤

≤ 5σ1∗ (|||∆U |||2F + |||∆V |||2F ),
where the first step follows from the upper bound of |||M p
− M ∗ |||F shown in Lemma 12, and the
second step follows from the assumption |||∆U |||F , |||∆V |||F ≤ σ1∗ .

22

6.8

Proof of Lemma 3

We first observe that
1
1
∇U G(U, V ) = U (U > U − V > V ), ∇V G(U, V ) = V (V > V − U > U ),
2
2
Therefore, we obtain
hh∇U G(U, V ), U − Uπ∗ ii + hh∇V G(U, V ), V − Vπ∗ ii
1
= hhU > U − V > V , U > U − V > V − U > Uπ∗ + V > Vπ∗ ii
2
1
1 >
= |||U U − V > V |||2F + hhU > U − V > V , U > U − V > V − 2U > Uπ∗ + 2V > Vπ∗ ii
4
4
1 >
1
>
2
>
= |||U U − V V |||F + hhU > U − V > V , U > U − V > V − 2∆>
U Uπ ∗ + 2∆V Vπ ∗ ii,
4
4

(45)

>
>
>
>
>
where the last step follows from ∆>
U Uπ ∗ − ∆V Vπ ∗ = U Uπ ∗ − V Vπ ∗ since Uπ ∗ Uπ ∗ = Vπ ∗ Vπ ∗ . Note
that

U > U − V > V = (Uπ∗ + ∆U )> (Uπ∗ + ∆U ) − (Vπ∗ + ∆V )> (Vπ∗ + ∆V )
>
>
>
>
= Uπ>∗ ∆U + ∆>
U Uπ ∗ + ∆U ∆U − Vπ ∗ ∆V − ∆V Vπ ∗ − ∆V ∆V ,

where we use Uπ>∗ Uπ∗ = Vπ>∗ Vπ∗ again in the last step. Furthermore, since U > U −V > V is symmetric,
we have
>
>
hhU > U − V > V , Uπ>∗ ∆U + ∆>
U Uπ ∗ − Vπ ∗ ∆V − ∆V Vπ ∗ ii
>
= hhU > U − V > V , 2∆>
U Uπ ∗ − 2∆V Vπ ∗ ii.

Using these arguments, for the second term in (45), denoted by T2 , we have
1
>
T2 = hhU > U − V > V , ∆>
U ∆U − ∆V ∆V ii.
4
Furthermore, we have
>
>
>
2
2
4T2 ≤ |hhU > U − V > V , ∆>
U ∆U − ∆V ∆V ii| ≤ |||U U − V V |||F |||∆U |||F + |||∆V |||F


≤ |||U > U − Uπ>∗ Uπ∗ |||F + |||V > V − Vπ>∗ Vπ∗ |||F δ
q
≤ 2 (|||Uπ∗ |||op |||∆U |||F + |||Vπ∗ |||op |||∆V |||F ) δ ≤ 2 2σ1∗ δ 3 .



(46)

It remains to find a lower bound of |||U > U − V > V |||F . The following inequality, which we turn
to prove later, is true:
|||U > U − V > V |||2F ≥ |||U U > − Uπ∗ Uπ>∗ |||2F + |||V V > − Vπ∗ Vπ>∗ |||2F − 2|||U V > − Uπ∗ Vπ>∗ |||2F .

(47)

Proceeding with the first term in (45) by using (47), we get
1 >
1
1
|||U U − V > V |||2F = |||U > U − V > V |||2F + |||U > U − V > V |||2F
4
8
8
1 >
1
1
1
≥ |||U U − V > V |||2F + |||U U > − Uπ∗ Uπ>∗ |||2F + |||V V > − Vπ∗ Vπ>∗ |||2F − |||U V > − Uπ∗ Vπ>∗ |||2F
8
8
8
4
1 >
1
1
>
2
>
>
2
>
>
2
= |||U U − V V |||F + |||F F − Fπ∗ Fπ∗ |||F − |||U V − Uπ∗ Vπ∗ |||F ,
(48)
8
8
2
23

where we let

 


U
Uπ∗
∗
F :=
, Fπ :=
.
V
Vπ ∗

Introduce ∆F := F − Fπ∗ . Recall that δ := |||∆U |||2F + |||∆V |||2F . Equivalently δ = |||∆F |||2F . We have
>
>
|||F F > − Fπ∗ Fπ>∗ |||F = |||∆F Fπ∗
+ Fπ∗ ∆>
F + ∆F ∆F |||F
>
2
>
>
≥ |||∆F Fπ∗
+ Fπ∗ ∆>
F |||F − |||∆F |||F = |||∆F Fπ∗ + Fπ ∗ ∆F |||F − δ.

For the first term, we have
>
2
> 2
>
>
|||∆F Fπ∗
+ Fπ∗ ∆>
F |||F = 2|||∆F Fπ∗ |||F + hh∆F Fπ∗ , Fπ ∗ ∆F ii
>
∗
2
>
>
≥ 2σr (Fπ∗ )2 |||∆F |||2F + hh∆F Fπ∗
, F π ∗ ∆>
F ii = 4σr |||∆F |||F + hh∆F Fπ∗ , Fπ ∗ ∆F ii.

For the cross term, by the following result, proved in [11] (we also provide a proof in Section 7.5 for
> , F ∗ ∆> ii ≥ 0.
the sake of completeness), we have hh∆F Fπ∗
π
F
√ ∗
Lemma 8. When |||F − Fπ∗ |||op < 2σr , we have that ∆>
F Fπ ∗ is symmetric.
p
p
Accordingly, we have |||F F > −Fπ∗ Fπ>∗ |||F ≥ 2 σr∗ δ −δ ≥ σr∗ δ under condition δ ≤ σr∗ . Plugging
this lower bound into (48), we obtain
1 >
1
1
1
|||U U − V > V |||2F ≥ |||U > U − V > V |||2F + σr∗ δ − |||U V > − Uπ∗ Vπ>∗ |||2F .
4
8
8
2
Putting (45), (46) and the above inequality together completes the proof.
Proof of inequality (47). For the term on the left hand side of (47), it is easy to check that
|||U > U − V > V |||2F = |||U U > |||2F + |||V V > |||2F − 2|||U V > |||2F .

(49)

The property Uπ>∗ Uπ∗ = Vπ>∗ Vπ∗ implies that |||Uπ∗ Uπ>∗ |||F = |||Vπ∗ Vπ>∗ |||F = |||Uπ∗ Vπ>∗ |||F . Therefore,
expanding those quadratic terms on the right hand side of (47), one can show that it is equal to
|||U U > |||2F + |||V V > |||2F − 2|||Uπ>∗ U |||2F − 2|||Vπ>∗ V |||2F + 4hhUπ>∗ U , Vπ>∗ V ii − 2|||U V > |||2F .
Comparing inequalities (49) and (50), it thus remains to show that
−2|||Uπ>∗ U |||2F − 2|||Vπ>∗ V |||2F + 4hhUπ>∗ U , Vπ>∗ V ii ≤ 0.
Equivalently, we always have |||Uπ>∗ U − Vπ>∗ V |||2F ≥ 0, and thus prove (47).

6.9

Proof of Lemma 4

First, we turn to prove (23). As
1
1
∇U G(U, V ) = U (U > U − V > V ), ∇V G(U, V ) = V (V > V − U > U ),
2
2
we have
|||∇U G(U, V )|||2F + |||∇V G(U, V )|||2F ≤


1
|||U |||2op + |||V |||2op |||U > U − V > V |||2F .
4
24

(50)

p
p
As (U, V ) p
∈ B2 ( σ1∗ ), we thus have |||U |||op ≤ |||Uπ∗ |||op + |||Uπ∗ − U |||op ≤ 2 σ1∗ , and similarly
|||V |||op ≤ 2 σ1∗ . We obtain
|||∇U G(U, V )|||2F + |||∇V G(U, V )|||2F ≤ 2σ1∗ |||U > U − V > V |||2F .
Now we turn to prove (22). We observe that
∇M L(U, V ; S) = M + S − M ∗ − S ∗ ,
where we let M := U V > . We denote the support of S, S ∗ by Ω and Ω∗ respectively. Based on the
sparse estimator (4) for computing S, ∇M L(U, V ; S) is only supported on Ωc . We thus have
|||∇M L(U, V ; S)|||F ≤ |||ΠΩc \Ω∗ (M − M ∗ )|||F + |||ΠΩc ∩Ω∗ (M − M ∗ − S ∗ )|||F
≤ |||M − M ∗ |||F + |||ΠΩc ∩Ω∗ (M − M ∗ − S ∗ )|||F .
It remains to upper bound the second term on the right hand side. Following (37) and (38), we
have
|||ΠΩc ∩Ω∗ (M − M ∗ − S ∗ )|||2F ≤

X
(i,j)∈Ωc ∩Ω∗

kvj k22
2
kui k22
+
≤
|||M − M ∗ |||2F ,
(γ − 1)αd2 (γ − 1)αd1
γ−1

where the last step is proved in (41). By choosing γ = 2, we thus conclude that
√
|||∇M L(U, V ; S)|||F ≤ (1 + 2)|||M − M ∗ |||F .

6.10

Proof of Lemma 6

We denote the support of ΠΦ (S ∗ ), S by Ω∗o and Ω. We always have Ω∗o ⊆ Φ and Ω ⊆ Φ.
In the sequel, we establish several results that characterize the properties of Φ. The first result,
proved in Section 7.2, shows that the Frobenius norm of any incoherent matrix whose row (or
column) space are equal to L∗ (or R∗ ) is well preserved under partial observations supported on Φ.
Lemma 9. Suppose M ∗ ∈ Rd1 ×d2 is a rank r and µ-incoherent matrix that has SVD M ∗ =
log d
L∗ Σ∗ R∗> . Then there exists an absolute constant c such that for any  ∈ (0, 1), if p ≥ c 2µr
,
(d1 ∧d2 )
then with probability at least 1 − 2d−3 , we have that for all A ∈ Rd2 ×r , B ∈ Rd1 ×r ,


(1 − )|||L∗ A> + BR∗> |||2F ≤ p−1 |||ΠΦ L∗ A> + BR∗> |||2F ≤ (1 + )|||L∗ A> + BR∗> |||2F .
We need the next result, proved in Section 7.3, to control the number of nonzero entries per row
and column in Ω∗o and Φ.
Lemma 10. If p ≥

56 log d
3 α(d1 ∧d2 ) ,

1
|Φ(i,·) | − pd2 ≤ pd2 ,
2

then with probability at least 1 − 6d−1 , we have
1
3
3
|Φ(·,j) | − pd1 ≤ pd1 , |Ω∗o(i,·) | ≤ αpd2 , |Ω∗o(·,j) | ≤ αpd1 ,
2
2
2

for all i ∈ [d1 ] and j ∈ [d2 ].

25

The next lemma, proved in Section 7.4, can be used to control the projection of small matrices
to Φ.
2 2

log d
Lemma 11. There exists constant c such that for any  ∈ (0, 1), if p ≥ c µ2 (dr 1 ∧d
, then with
2)

probability p
at least 1 − O(d−1 ),pfor all matrices Z ∈ Rd1 ×d2 , U ∈ Rd1 ×r and V ∈ Rd2 ×r that satisfy
|||U |||2,∞ ≤ µr/d1 ,|||V |||2,∞ ≤ µr/d2 , we have
p−1 |||ΠΦ (U V > )|||2F ≤ |||U |||2F |||V |||2F + |||U |||F |||V |||F ;

(51)

p−1 |||ΠΦ (Z)V |||2F ≤ 2µr|||ΠΦ (Z)|||2F ;

(52)

p−1 |||U > ΠΦ (Z)|||2F ≤ 2µr|||ΠΦ (Z)|||2F .

(53)

In the remainder of this section, we condition on the events in Lemmas 9, 10 and 11. Now we
are ready to prove Lemma 6.
Proof of Lemma 6. Using shorthand M := U V > , we have
e V ; S) = p−1 ΠΦ (M + S − M ∗ − S ∗ ) .
∇M L(U,
Plugging it back into the left hand side of (21), we obtain
e V ; S), U V > − Uπ∗ Vπ>∗ + ∆U ∆>
hh∇M L(U,
V ii
1
= hhΠΦ (M + S − M ∗ − S ∗ ), M − M ∗ + ∆U ∆>
V ii
p
1
1
1
≥ |||ΠΦ (M − M ∗ ) |||2F − |hhΠΦ (S − S ∗ ), M − M ∗ ii| − |hhΠΦ (M + S − M ∗ − S ∗ ), ∆U ∆>
V ii| .
p
p
p
|
{z
} |
{z
} |
{z
}
T1

T2

T3

(54)
Next we derive lower bounds of T1 , upper bounds of T2 and T3 respectively.
>
>
Lower bound of T1 . We observe that M −M ∗ = Uπ∗∗ ∆>
V +∆U Vπ ∗ +∆U ∆V . By triangle inequality,
we have
>
>
|||ΠΦ (M − M ∗ )|||F ≥ |||ΠΦ (Uπ∗ ∆>
V + ∆U Vπ ∗ )|||F − |||ΠΦ (∆U ∆V )|||F .

Note that when c ≥ a − b for a, b ≥ 0, we always have c2 ≥ 21 a2 − b2 . We thus have
T1 ≥
≥
≥
≥
≥

1
1
> 2
> 2
|||ΠΦ (Uπ∗ ∆>
V + ∆U Vπ ∗ )|||F − |||ΠΦ (∆U ∆V )|||F
2p
p
1
1
> 2
> 2
(1 − )|||Uπ∗ ∆>
V + ∆U Vπ ∗ |||F − |||ΠΦ (∆U ∆V )|||F
2
p
1
2
2
2
∗
(1 − )|||M − M ∗ − ∆U ∆>
V |||F − |||∆U |||F |||∆V |||F − 9σ1 |||∆U |||F |||∆V |||F
2
1
1
(1 − )|||M − M ∗ |||2F − (1 − )|||∆U ∆V |||2F − |||∆U |||2F |||∆V |||2F − 9σ1∗ |||∆U |||F |||∆V |||F
4
2
1
∗ 2
(1 − )|||M − M |||F − 2δ 2 − 5σ1∗ δ.
4

where the second step is implied
by Lemma 9, the third p
step follows from (51) in Lemma 11 by
p
∗
noticing that |||∆U |||2,∞ ≤ 3 µrσ1 /d1 and |||∆V |||2,∞ ≤ 3 µrσ1∗ /d1 , which is further implied by
(31).
26

Upper bound of T2 . Since S − S ∗ is supported on Ω∗0 ∪ Ω, we have
pT2 ≤ |hhΠΩ∗o \Ω (S ∗ ), ΠΩ∗o \Ω (M − M ∗ )ii| + |hhΠΩ (S − S ∗ ), ΠΩ (M − M ∗ )ii|.

(55)

For any (i, j) ∈ Ω, we have (S − S ∗ )(i,j) = (M ∗ − M )(i,j) . Therefore, for the second term on the
right hand side, we have
|hhΠΩ (S − S ∗ ), ΠΩ (M − M ∗ )ii| ≤ |||ΠΩ (M − M ∗ )|||2F ≤ 18γpαµrσ1∗ δ,

(56)

where the last inequality follows from Lemma 14 and the fact that |Ω(i,·) | ≤ γpαd2 , |Ω(·,j) | ≤ γpαd1
for all i ∈ [d1 ], j ∈ [d2 ].
We denote the i-th row of ΠΦ (M − M ∗ ) by ui , and we denote the j-th column of ΠΦ (M − M ∗ )
(k)
(k)
by vj . We let ui denote the element of ui that has the k-th largest magnitude. We let vj denote
the element of vj that has the k-th largest magnitude.
For the first term on the right hand side of (55), we first observe that for (i, j) ∈ Ω∗o \ Ω, |(M ∗ +
∗
S − M )(i,j) | is either less than the γpαd2 -th largest element in the i-th row of ΠΦ (M ∗ + S ∗ − M ),
or less than γpαd1 -th largest element in the j-th row of ΠΦ (M ∗ + S ∗ − M ). Based on Lemma
10, ΠΦ (S ∗ ) has at most 3pαd2 /2 nonzero entries per row and at most 3pαd1 /2 nonzero entries per
column. Therefore, we have
n
o
((γ−1.5)pαd2 )
((γ−1.5)pαd1 )
|(M ∗ + S ∗ − M )(i,j) | ≤ max |ui
|, |vj
| .
(57)
In addition, we observe that
|hhΠΩ∗o \Ω (S ∗ ), ΠΩ∗o \Ω (M − M ∗ )ii|
X
≤
|(M ∗ + S ∗ − M )(i,j) ||(M ∗ − M )(i,j) | + |(M ∗ − M )(i,j) |2
(i,j)∈Ω∗o \Ω


≤

β
1+
2



|||ΠΩ∗o (M ∗ − M )|||2F +

≤ (27 + 14β)pαµrσ1∗ δ +

1
2β

1
2β

X

|(M ∗ + S ∗ − M )(i,j) |2 ,

(i,j)∈Ω∗o \Ω

X

|(M ∗ + S ∗ − M )(i,j) |2 ,

(58)

(i,j)∈Ω∗o \Ω

where the second step holds for any β > 0 and the last step follows from Lemma 14 under the size
constraints of Ω∗o shown in Lemma 10. For the second term in (58), using (57), we have
X
X
((γ−1.5)pαd2 ) 2
((γ−1.5)pαd1 ) 2
|(M ∗ + S ∗ − M )(i,j) |2 ≤
|ui
| + |vj
|
(i,j)∈Ω∗o \Ω

=

X

(i,j)∈Ω∗o

X

((γ−1.5)pαd2 ) 2
|ui
|

i∈[d1 ] j∈Ω∗o(i,·)

≤

X
i∈[d1 ]

+

X

X

((γ−1.5)pαd1 ) 2

|vj

|

j∈[d2 ] i∈Ω∗o(·,j)

X 1.5
1.5
3
kui k22 +
kvj k22 ≤
|||ΠΦ (M − M ∗ )|||2F .
γ − 1.5
γ − 1.5
γ − 1.5

(59)

j∈[d2 ]

Moreover, we have
> 2
> 2
|||ΠΦ (M − M ∗ )|||2F ≤ 2|||ΠΦ (Uπ∗ ∆>
V + ∆U Vπ ∗ )|||F + 2|||ΠΦ (∆U ∆V )|||F
> 2
2
2
∗
≤ 2(1 + )p|||Uπ∗ ∆>
V + ∆U Vπ ∗ |||F + 2p|||∆U |||F |||∆V |||F + 18pσ1 |||∆U |||F |||∆V |||F

≤ 4(1 + )p |||Uπ∗ |||2op |||∆V |||2F + |||Vπ∗ |||2op |||∆U |||2F + 2p|||∆U |||2F |||∆V |||2F + 18pσ1∗ |||∆U |||F |||∆V |||F

≤ (4 + 13)pσ1∗ δ + 2pδ 2 ,

(60)
27

where the second step follows from Lemma 9 and inequality (51) in Lemma 11. Putting (55)-(60)
together, we obtain
T2 ≤ (18γ + 14β + 27)αµrσ1∗ δ +

3[(2 + 7)σ1∗ δ + δ 2 ]
.
β(γ − 1.5)

Upper bound of T3 . By Cauchy-Schwarz inequality, we have
pT3 ≤ |||ΠΦ (M − M ∗ + S − S ∗ )|||F |||ΠΦ (∆U ∆>
V )|||F
q
≤ |||ΠΦ (M − M ∗ + S − S ∗ )|||F p|||∆U |||2F |||∆V |||2F + 9pσ1∗ |||∆U |||F |||∆V |||F
q
≤ |||ΠΦ (M − M ∗ + S − S ∗ )|||F pδ 2 + 5pσ1∗ δ.
where we use (51) in Lemma 11 in the second step.
We observe that ΠΦ (M − M ∗ + S − S ∗ ) is supported on Φ \ Ω. Therefore, we have
|||ΠΦ (M − M ∗ + S − S ∗ )|||F ≤ |||ΠΦ∩Ωc ∩Φ∗c (M − M ∗ )|||F + |||ΠΦ∩Ωc ∩Φ∗ (M − M ∗ − S ∗ )|||F
≤ |||ΠΦ (M − M ∗ )|||F + |||ΠΩc ∩Φ∗ (M − M ∗ − S ∗ )|||F
r
3
∗
≤ |||ΠΦ (M − M )|||F +
|||ΠΦ (M − M ∗ )|||F
γ − 1.5
r

q
3
(4 + 13)pσ1∗ δ + 2pδ 2 ,
≤ 1+
γ − 1.5
where the third step follows from (59), and the last step is from (60). Under assumptions γ = 3,
 ≤ 1/4 and δ ≤ σ1∗ , we have
q
q
q
√
T3 ≤ 3 9σ1∗ δ + 2δ 2 δ 2 + 5σ1∗ δ ≤ 10 σ1∗ δ 3 + 23 σ1∗ δ.
Combining pieces. Under the aforementioned assumptions, putting all pieces together leads to
e V ; S), U V > − Uπ∗ V >∗ + ∆U ∆> ii
hh∇M L(U,
π
V


q
√
3
18
∗ 2
∗
≥ |||M − M |||F − (14β + 81)αµrσ1 δ − 26  +
σ1∗ δ − 10 σ1∗ δ 3 − 2δ 2 .
16
β

6.11

Proof of Lemma 7

Let M := U V > . We find that
e V ; S) = p−1 ΠΦ (M + S − M ∗ − S ∗ ) V,
∇U L(U,
e V ; S) = p−1 ΠΦ (M + S − M ∗ − S ∗ )> U.
∇V L(U,
Conditioning on the event in Lemma 11, since (U, V ) ∈ U × V, inequalities (52) and (53) imply that
e V ; S)|||2F + |||∇V L(U,
e V ; S)|||2F ≤ 12 µrσ1∗ |||ΠΦ (M + S − M ∗ − S ∗ ) |||2F .
|||∇U L(U,
p
28

It remains to bound the term |||ΠΦ (M + S − M ∗ − S ∗ ) |||2F . Let Ω∗o and Ω be the support of ΠΦ (S ∗ )
and S respectively. We observe that
|||ΠΦ (M + S − M ∗ − S ∗ ) |||2F = |||ΠΩ∗o \Ω (M − M ∗ − S ∗ ) |||2F + |||ΠΦ∗c ∩Ωc ∩Φ (M − M ∗ ) |||2F
≤ |||ΠΩ∗o \Ω (M − M ∗ − S ∗ ) |||2F + |||ΠΦ (M − M ∗ ) |||2F .
In the proof of Lemma 6, it is shown in (59) that
|||ΠΩ∗o \Ω (M − M ∗ − S ∗ ) |||2F ≤

3
|||ΠΦ (M − M ∗ )|||2F .
γ − 1.5

Moreover, following (60), we have that
> 2
2
2
∗
|||ΠΦ (M − M ∗ )|||2F ≤ 2(1 + )p|||Uπ∗ ∆>
V + ∆U Vπ ∗ |||F + 2p|||∆U |||F |||∆V |||F + 18pσ1 |||∆U |||F |||∆V |||F

≤ 4(1 + )p|||M − M ∗ |||2F + (6 + 4)p|||∆U |||2F |||∆V |||2F + 18pσ1∗ |||∆U |||F |||∆V |||F
≤ 4(1 + )p|||M − M ∗ |||2F + (6 + 4)pδ 2 + 9pσ1∗ δ.
We thus finish proving our conclusion by combining all pieces and noticing that γ = 3 and  ≤ 1/4.

7

Proofs for Technical Lemmas

In this section, we prove several technical lemmas that are used in the proofs of our main theorems.

7.1

Proof of Lemma 1

We observe that
|||A|||op =

sup

sup x> Ay.

x∈Sd1 −1 y∈Sd2 −1

We denote the support of A by Ω. For any x ∈ Rd1 , y ∈ Rd2 and β > 0, we have
X 1
kAk∞ (β −1 x2i + βyj2 )
2
(i,j)∈Ω
(i,j)∈Ω


X X
X X
1
= kAk∞ 
β −1 x2i +
βyj2 
2

x> Ay =

X

xi A(i,j) yj ≤

i

j∈Ω(i,·)

j

i∈Ω(·,j)


1
≤ kAk∞ αd2 β −1 kxk22 + αd1 βkyk22 .
2
It is thus implied that |||A|||op ≤ 21 α(β −1 d2 + βd1 )kAk∞ . Choosing β =

7.2

p
d2 /d1 completes the proof.

Proof of Lemma 9

We define a subspace K ⊆ Rd1 ×d2 as
n
o
K := X X = L∗ A> + BR∗> for some A ∈ Rd2 ×r , B ∈ Rd1 ×r .

29

Let ΠK be Euclidean projection onto K. Then according to Theorem 4.1 in [6], under our assumptions, for all matrices X ∈ Rd1 ×d2 , inequality
p−1 ||| (ΠK ΠΦ ΠK − pΠK ) X|||F ≤ |||X|||F

(61)

holds with probability at least 1 − 2d−3 .
In our setting, by restricting X = L∗ A> + BR∗> , we have ΠK X = X. Therefore, (61) implies
that
|||ΠK ΠΦ X − pX|||F ≤ p|||X|||F .
For |||ΠΦ X|||2F , we have
|||ΠΦ X|||2F = hhΠΦ X, ΠΦ Xii = hhΠΦ X, Xii
= hhΠK ΠΦ X, Xii ≤ |||ΠK ΠΦ X|||F |||X|||F ≤ p(1 + )|||X|||2F .
On the other hand, we have
|||ΠΦ X|||2F = hhΠK ΠΦ X, Xii = hhΠK ΠΦ X − pX + pX, Xii
= p|||X|||2F − hhX, −ΠK ΠΦ X + pXii
≥ p|||X|||2F − |||X|||F |||ΠK ΠΦ X − pX|||F ≥ p(1 − )|||X|||2F .
Combining the above two inequalities, we complete the proof.

7.3

Proof of Lemma 10

We observe that |Φ(i,·) | is a summation of d2 i.i.d. binary random variables with mean p and variance
p(1 − p). By Bernstein’s inequality, for any i ∈ [d1 ],
!




− 21 (pd2 /2)2
1
3
Pr |Φ(i,·) | − pd2 ≥ pd2 ≤ 2 exp −
≤
2
exp
−
pd
2 .
2
28
d2 p(1 − p) + 13 (pd2 /2)
By probabilistic union bound, we have
"
Pr sup |Φ(i,·) | − pd2
i∈[d1 ]

#


1
3
≥ pd2 ≤ 2d1 exp − pd2 ≤ 2d−1 ,
2
28

log d
where the last inequality holds by assuming p ≥ 56
3 d2 .
The term |Ω∗o(i,·) | is a summation of at most αd2 i.i.d. binary random variables with mean p and
variance p(1 − p). Again, applying Bernstein’s inequality leads to




h
i 1
3
∗
∗
Pr |Ωo(i,·) | − E |Ωo(i,·) | ≥ pαd2 ≤ exp − pαd2 .
2
28

Accordingly, by the assumption p ≥
"

56 log d
3 αd2 ,

we obtain
#


3
1
∗
Pr sup |Ωo(i,·) | − pk ≥ pk ≤ d1 exp − pαd2 ≤ d−1 .
2
28
i∈[d1 ]

The proofs for |Φ(·,j) | and |Ω∗o(·,j) | follow the same idea.
30

7.4

Proof of Lemma 11

d
d1 ×d2 , we
According to Lemma 3.2 in [8], under condition p ≥ c1 µd1log
∧d2 , for any fixed matrix A ∈ R
have
s
d log d
|||A − p−1 ΠΦ A|||op ≤ c2
kAk∞ ,
p

holds with probability at least 1 − O(d−3 ). Letting A be all-ones matrix, then we have that for all
u ∈ Rd1 , v ∈ Rd2 ,
X
p
ui vj ≤ pkuk1 kvk1 + c2 pd log dkuk2 kvk2 .
(i,j)∈Φ

We find that
|||ΠΦ (U V > )|||2F ≤

X

kU(i,·) k22 kV(j,·) k22

(i,j)∈Φ

sX
sX
p
4
≤ p|||U |||F |||V |||F + c2 pd log d
kU(i,·) k2
kV(j,·) k42
2

2

i∈[d1 ]

j∈[d2 ]

p
≤ p|||U |||2F |||V |||2F + c2 pd log d|||U |||F |||V |||F |||U |||2,∞ |||V |||2,∞
s
pµ2 r2 d log d
≤ p|||U |||2F |||V |||2F + c2
|||U |||F |||V |||F .
d1 d2
By the assumption p &

µ2 r2 log d
,
2 (d1 ∧d2 )

we finish proving (51).

d
According to the proof of Lemma 10, if p ≥ c dlog
, with probability at least 1 − O(d−1 ), we
1 ∧d2
3
3
have |Φ(i,·) | ≤ 2 pd2 and |Φ(·,j) | ≤ 2 pd1 for all i ∈ [d1 ] and j ∈ [d2 ]. Conditioning on this event, we
have
X X
|||ΠΦ (Z)V |||2F =
h(ΠΦ (Z))(i,·) , H(·,k) i2
i∈[d1 ] k∈[r]

≤

X X

k(ΠΦ (Z))(i,·) k22

2
V(j,k)

j∈Ω(i,·)

i∈[d1 ] k∈[r]

= |||ΠΦ Z|||2F

X

X

kV(i,·) k22

j∈Ω(i,·)

3
≤ |||ΠΦ Z|||2F pd2 · |||V |||22,∞ ≤ 2µrp|||ΠΦ Z|||2F .
2
We thus finish proving (52). Inequality (53) can be proved in the same way.

7.5

Proof of Lemma 8

Recall that we let F := [U ; V ] and Fπ∗ := [U ∗ ; V ∗ ]Q for some matrix Q ∈ Qr , which minimizes the
following function
|||F − [U ∗ ; V ∗ ]Q|||2F .
(62)
Let F ∗ := [U ∗ ; V ∗ ]. Expanding the above term, we find that Q is the maximizer of hhF, F ∗ Qii =
Tr(F > F ∗ Q). Suppose F > F ∗ has SVD with form Q1 ΛQ>
2 for Q1 , Q2 ∈ Qr . When the minimum
31

diagonal term of Λ is positive, we conclude that the minimizer of (62) is unique and Q = Q2 Q>
1.
To prove this argument, we note that
X
Λ(i,i) hpi , qi i,
Tr(F > F ∗ Q) =
i∈[r]

where
pi is the i-th column of Q1 and qi is the i-th column of Q> Q2 . Hence, Tr(F > F ∗ Q) ≤
P
i∈[r] Λ(i,i) and the equality holds if and only if pi = qi for all i ∈ [r] since every Λ(i,i) > 0. We
have Q1 = Q> Q2 and thus finish proving the argument.
√
Under our assumption |||F − Fπ∗ |||op < 2σr∗ , for any nonzero vector u ∈ Rr , we have
p
kF > Fπ∗ uk2 ≥ kFπ>∗ Fπ∗ uk2 − k(Fπ∗ − F )> Fπ∗ uk2 ≥ ( 2σr∗ − |||Fπ∗ − F |||op )|||Fπ∗ u|||F > 0.
In √
the second step, we use the fact that the singular values of Fπ∗ are equal to the diagonal terms
of 2Σ∗1/2 . Hence, F > Fπ∗ has full rank. Furthermore, it implies that F > F ∗ has full rank and only
contains positive singular values.
Proceeding with the proved argument, we have
>
>
F > Fπ∗ = Q1 ΛQ>
2 Q2 Q1 = Q1 ΛQ1 ,

which implies that F > Fπ∗ is symmetric. Accordingly, we have (F − Fπ∗ )> Fπ∗ is also symmetric.

Acknowledgment
Y. Chen acknowledges support from the School of Operations Research and Information Engineering,
Cornell University.

References
[1]

Animashree Anandkumar, Rong Ge, Daniel Hsu, Sham M. Kakade, and Matus Telgarsky.
“Tensor decompositions for learning latent variable models”. In The Journal of Machine Learning Research 15.1 (2014), pp. 2773–2832.

[2]

Sivaraman Balakrishnan, Martin J. Wainwright, and Bin Yu. “Statistical guarantees for the
EM algorithm: From population to sample-based analysis”. In arXiv preprint arXiv:1408.2156
(2014).

[3]

Srinadh Bhojanapalli, Prateek Jain, and Sujay Sanghavi. “Tighter low-rank approximation
via sampling the leveraged element”. In Proceedings of the Twenty-Sixth Annual ACM-SIAM
Symposium on Discrete Algorithms. SIAM. 2015, pp. 902–920.

[4]

Srinadh Bhojanapalli, Anastasios Kyrillidis, and Sujay Sanghavi. “Dropping convexity for
faster semi-definite optimization”. In arXiv preprint arXiv:1509.03917 (2015).

[5]

Emmanuel J. Candès, Xiaodong Li, and Mahdi Soltanolkotabi. “Phase retrieval via Wirtinger
flow: Theory and algorithms”. In IEEE Transactions on Information Theory 61.4 (2015),
pp. 1985–2007.

[6]

Emmanuel J. Candès and Benjamin Recht. “Exact matrix completion via convex optimization”. In Foundations of Computational mathematics 9.6 (2009), pp. 717–772.
32

[7]

Emmanuel J. Candès and Terence Tao. “The power of convex relaxation: Near-optimal matrix
completion”. In IEEE Transactions on Information Theory 56.5 (2010), pp. 2053–2080.

[8]

Emmanuel J. Candès, Xiaodong Li, Yi Ma, and John Wright. “Robust principal component
analysis?” In Journal of the ACM (JACM) 58.3 (2011), p. 11.

[9]

Venkat Chandrasekaran, Sujay Sanghavi, Pablo A. Parrilo, and Alan S. Willsky. “Ranksparsity incoherence for matrix decomposition”. In SIAM Journal on Optimization 21.2 (2011),
pp. 572–596.

[10]

Yudong Chen. “Incoherence-Optimal Matrix Completion”. In IEEE Transactions on Information Theory 61.5 (2015), pp. 2909–2923.

[11]

Yudong Chen and Martin J. Wainwright. “Fast low-rank estimation by projected gradient
descent: General statistical and algorithmic guarantees”. In arXiv preprint arXiv:1509.03025
(2015).

[12]

Yudong Chen, Ali Jalali, Sujay Sanghavi, and Constantine Caramanis. “Low-rank Matrix Recovery from Errors and Erasures”. In IEEE Transactions on Information Theory 59.7 (2013),
pp. 4324–4337.

[13]

Yuxin Chen and Emmanuel J. Candès. “Solving random quadratic systems of equations is
nearly as easy as solving linear systems”. In Advances in Neural Information Processing Systems. 2015, pp. 739–747.

[14]

Yeshwanth Cherapanamjeri, Kartik Gupta, and Prateek Jain. “Nearly-optimal Robust Matrix
Completion”. In arXiv preprint arXiv:1606.07315 (2016).

[15]

Kenneth L. Clarkson and David P. Woodruff. “Low rank approximation and regression in
input sparsity time”. In Proceedings of the forty-fifth annual ACM symposium on Theory of
computing. ACM. 2013, pp. 81–90.

[16]

Alan Frieze, Ravi Kannan, and Santosh Vempala. “Fast Monte-Carlo algorithms for finding
low-rank approximations”. In Journal of the ACM (JACM) 51.6 (2004), pp. 1025–1041.

[17]

Quanquan Gu, Zhaoran Wang, and Han Liu. “Low-Rank and Sparse Structure Pursuit via
Alternating Minimization”. In Proceedings of the 19th International Conference on Artificial
Intelligence and Statistics. 2016, pp. 600–609.

[18]

Moritz Hardt. “Understanding alternating minimization for matrix completion”. In 2014 IEEE
55th Annual Symposium on Foundations of Computer Science (FOCS). IEEE. 2014, pp. 651–
660.

[19]

Daniel Hsu, Sham M. Kakade, and Tong Zhang. “Robust matrix decomposition with sparse
corruptions”. In IEEE Transactions on Information Theory 57.11 (2011), pp. 7221–7234.

[20]

Prateek Jain, Praneeth Netrapalli, and Sujay Sanghavi. “Low-rank matrix completion using
alternating minimization”. In Proceedings of the forty-fifth annual ACM symposium on Theory
of computing. ACM. 2013, pp. 665–674.

[21]

Zhouchen Lin, Minming Chen, and Yi Ma. “The Augmented Lagrange Multiplier Method
for Exact Recovery of Corrupted Low-Rank Matrices”. In Arxiv preprint arxiv:1009.5055v3
(2013).

33

[22]

Praneeth Netrapalli, UN Niranjan, Sujay Sanghavi, Animashree Anandkumar, and Prateek
Jain. “Non-convex robust PCA”. In Advances in Neural Information Processing Systems. 2014,
pp. 1107–1115.

[23]

Yousef Saad. Numerical Methods for Large Eigenvalue Problems: Revised Edition. Vol. 66.
Siam, 2011.

[24]

Ju Sun, Qing Qu, and John Wright. “When Are Nonconvex Problems Not Scary?” In arXiv
preprint arXiv:1510.06096 (2015).

[25]

Ruoyu Sun and Zhi-Quan Luo. “Guaranteed matrix completion via nonconvex factorization”.
In 2015 IEEE 56th Annual Symposium on Foundations of Computer Science (FOCS). IEEE.
2015, pp. 270–289.

[26]

Stephen Tu, Ross Boczar, Mahdi Soltanolkotabi, and Benjamin Recht. “Low-rank solutions
of linear matrix equations via procrustes flow”. In arXiv preprint arXiv:1507.03566 (2015).

[27]

Zhaoran Wang, Quanquan Gu, Yang Ning, and Han Liu. “High Dimensional EM Algorithm:
Statistical Optimization and Asymptotic Normality”. In Advances in Neural Information Processing Systems. 2015, pp. 2512–2520.

[28]

Huan Xu, Constantine Caramanis, and Sujay Sanghavi. “Robust PCA via Outlier Pursuit”.
In IEEE Transactions on Information Theory 58.5 (2012), pp. 3047–3064. issn: 0018-9448.

[29]

Xinyang Yi and Constantine Caramanis. “Regularized EM Algorithms: A Unified Framework
and Statistical Guarantees”. In Advances in Neural Information Processing Systems. 2015,
pp. 1567–1575.

[30]

Huishuai Zhang, Yuejie Chi, and Yingbin Liang. “Provable Non-convex Phase Retrieval with
Outliers: Median Truncated Wirtinger Flow”. In arXiv preprint arXiv:1603.03805 (2016).

[31]

Tuo Zhao, Zhaoran Wang, and Han Liu. “A Nonconvex Optimization Framework for Low Rank
Matrix Estimation”. In Advances in Neural Information Processing Systems. 2015, pp. 559–
567.

[32]

Qinqing Zheng and John Lafferty. “A convergent gradient descent algorithm for rank minimization and semidefinite programming from random linear measurements”. In Advances in
Neural Information Processing Systems. 2015, pp. 109–117.

34

Appendices
A

Supporting Lemmas

In this section, we provide several technical lemmas used for proving our main results.
Lemma 12. For any (U ∗ , V ∗ ) ∈ E(M ∗ ), U ∈ Rd1 ×r and V ∈ Rd2 ×r , we have
p
|||U V > − U ∗ V ∗> |||F ≤ σ1∗ (|||∆V |||F + |||∆U |||F ) + |||∆U |||F |||∆V |||F ,
where ∆U := U − U ∗ , ∆V := V − V ∗ .
∗> + ∆ ∆> . Hence,
Proof. We observe that U V > − U ∗ V ∗> = U ∗ ∆>
U V
V + ∆U V
∗>
|||U V > − U ∗ V ∗> |||F ≤ |||U ∗ ∆>
|||F + |||∆U ∆>
V |||F + |||∆U V
V |||F

≤ |||U ∗ |||op |||∆V |||F + |||V ∗ |||op |||∆U |||F + |||∆U |||F |||∆V |||F .

Furthermore, assuming (U, V ) ∈ U × V, where U and V satisfy the conditions in (19), we have
the next result.
Lemma 13. For any (i, j) ∈ [d1 ] × [d2 ], we have
r
r
µrσ1∗
µrσ1∗
>
∗ ∗>
|(U V − U V )(i,j) | ≤ 3
k∆V (j,·) k2 + 3
k∆U (i,·) k2
d1
d2

(63)

Proof. We observe that
∗
∗
|(U V > − M ∗ )(i,j) | ≤ |hU(i,·)
, ∆V (j,·) i| + |hV(j,·)
, ∆U (i,·) i| + |h∆U (i,·) , ∆V (j,·) i|
r
r
µrσ1∗
µrσ1∗
1
1
≤
k∆V (j,·) k2 +
k∆U (i,·) k2 + |||∆U |||2,∞ k∆V (j,·) k2 + |||∆V |||2,∞ k∆U (i,·) k2 .
d1
d2
2
2

By noticing that
r

∗

|||∆U |||2,∞ ≤ |||U |||2,∞ + |||U |||2,∞ ≤ 3

r
µrσ1∗
µrσ1∗
∗
, |||∆V |||2,∞ ≤ |||V |||2,∞ + |||V |||2,∞ ≤ 3
,
d1
d2

we complete the proof.
Lemma 13 can be used to prove the following result.
Lemma 14. For any α ∈ [0, 1], suppose Ω ⊆ [d1 ] × [d2 ] satisfies |Ω(i,·) | ≤ αd2 for all i ∈ [d1 ] and
|Ω(·,j) | ≤ αd1 for all j ∈ [d2 ]. Then we have
|||ΠΩ (U V > − U ∗ V ∗> )|||2F ≤ 18αµrσ1∗ (|||∆V |||2F + |||∆U |||2F ).

35

Proof. Using Lemma 13 for bounding each entry of U V > − U ∗ V ∗> , we have that
X
|(U V > − U ∗ V ∗> )(i,j) |2
|||ΠΩ (U V > − U ∗ V ∗> )|||2F ≤
(i,j)∈Ω

X 18µrσ ∗
18µrσ1∗
1
≤
k∆V (j,·) k22 +
k∆U (i,·) k22
d1
d2
(i,j)∈Ω

X X 18µrσ ∗
X X 18µrσ ∗
1
1
≤
k∆V (j,·) k22 +
k∆U (i,·) k22
d1
d2
j

≤

i

i∈Ω(·,j)

18αµrσ1∗ (|||∆V |||2F

j∈Ω(i,·)

2

+ |||∆U |||F ).

Denote the i-th largest singular value of matrix M by σi (M ).
Lemma 15 (Lemma 5.14 in [26]). Let M1 , M2 ∈ Rd1 ×d2 be two rank r matrices. Suppose they have
SVDs M1 = L1 Σ1 R1> and M2 = L2 Σ2 R2> . Suppose |||M1 − M2 |||op ≤ 21 σr (M1 ). Then we have
1/2
1/2
1/2
1/2
d2 (L2 Σ2 , R2 Σ2 ; L1 Σ1 , R1 Σ1 ) ≤ √

B

2 |||M2 − M1 |||2F
.
2 − 1 σr (M1 )

Parameter Settings for FB Separation Experiments

We approximate the FB separation problem by the RPCA framework with r = 10, α = 0.2, µ = 10.
Our algorithmic parameters are set as γ = 1, η = 1/(2σ̂1∗ ), where σ̂1∗ is an estimate of σ1∗ obtained
from the initial SVD. The parameters of AltProj
are kept as provided in the default setting. For
√
IALM, we use the tradeoff paramter λ = 1/ d1 , where d1 is the number of pixels in each frame
(the number of rows in Y ).
Note that both IALM and AltProj use the stopping criterion
|||Y − Mt − St |||F /|||Y |||F ≤ 10−3 .
Our algorithm for the partial observation setting never explicitly forms the d1 -by-d2 matrix Mt =
Ut Vt> , which is favored in large scale problems, but also renders the above criterion inapplicable.
Instead, we use the following stopping criterion
|||Ut+1 − Ut |||2F + |||Vt+1 − Vt |||2F
≤ 4 × 10−4 .
|||Ut |||2F + |||Vt |||2F
This rule checks whether the iterates corresponding to low-rank factors become stable. In fact, our
stopping criterion seems more natural and practical because in most real applications, matrix Y
cannot be strictly decomposed into low-rank M and sparse S that satisfy Y = M + S. Instead of
forcing M + S to be close to Y , our rule relies on seeking a robust subspace that captures the most
variance of Y .

36

