arXiv:1709.04836v1 [stat.ML] 14 Sep 2017

Informed Non-convex Robust Principal
Component Analysis with Features
Niannan Xue, Jiankang Deng, Yannis Panagakis, Stefanos Zafeiriou
Department of Computing, Imperial College London, UK, SW7 2AZ
Abstract
We revisit the problem of robust principal component analysis with
features acting as prior side information. To this aim, a novel, elegant,
non-convex optimization approach is proposed to decompose a given observation matrix into a low-rank core and the corresponding sparse residual. Rigorous theoretical analysis of the proposed algorithm results in
exact recovery guarantees with low computational complexity. Aptly designed synthetic experiments demonstrate that our method is the first to
wholly harness the power of non-convexity over convexity in terms of both
recoverability and speed. That is, the proposed non-convex approach is
more accurate and faster compared to the best available algorithms for the
problem under study. Two real-world applications, namely image classification and face denoising further exemplify the practical superiority of
the proposed method.

1

Introduction

Many machine learning and artificial intelligence tasks involve the separation of
a data matrix into a low-rank structure and a sparse part capturing different
information. Robust principal component analysis (RPCA) Candes et al. [2011],
Chandrasekaran et al. [2011] is a popular framework that logically characterizes
this matrix separation problem.
Nevertheless, prior side information, oftentimes in the form of features, may
also be present in practice. For instance, features are available for the following
tasks:
– Collaborative filtering: apart from ratings of an item by other users, the
profile of the user and the description of the item can also be exploited in
making recommendations Chiang et al. [2015];
– Relationship prediction: user behaviours and message exchanges can assist
in finding missing links on social media networks Xu et al. [2013];
– Person-specific facial deformable models: an orthonormal subspace learnt
from manually annotated data captured in-the-wild, when fed into an im1

age congealing procedure, can help produce more correct fittings Sagonas
et al. [2014].
It is thus reasonable to investigate how propitious it is for RPCA to exploit the
available features. Indeed, recent results Liu et al. [2017] indicate that features
are not redundant at all. In the setting of multiple subspaces, RPCA degrades
as the number of subspaces grows because of the increased row-coherence. On
the other hand, the use of feature dictionaries allows accurate low-rank recovery by removing the dependency on row-coherence. Despite the theoretical and
practical merits of convexified RPCA with features, such as LRR Liu et al.
[2010] and PCPF Chiang et al. [2016], convex relaxations of the rank function
and l0 -norm necessarily lead into algorithm weakening Chandrasekarana and
Jordan [2013].
On a separate note, recent advances in non-convex optimization algorithms
continue to undermine their convex counterparts Gong et al. [2013], Ge et al.
[2016], Kohler and Lucchi [2017]. In particular, non-convex RPCA algorithms
such as fast RPCA Yi et al. [2016] and AltProj Netrapalli et al. [2014] exhibit better properties than the convex formulation. Most recently, Niranjan
et al. [2017] embedded features into a non-convex RPCA framework known as
IRPCA-IHT with faster speed. However, it remains unclear as to whether features have been effectively incorporated into non-convex RPCA and the benefits
of accuracy, speed and so on have been exploited as much as possible.
In this work, we give positive answers to the above questions by proposing
a novel, non-convex scheme that fully leverages features, which reveal true row
and column subspaces, to decompose an observation matrix into a core matrix
with given rank and a residual part with informed sparsity. Even though the
proposed algorithm is inspired by the recently proposed fast RPCA Yi et al.
[2016], our contributions are by no means trivial, especially from a theoretical
perspective. First, fast RPCA cannot be easily extended to consistently take account of features. Second, as we show in this paper, incoherence assumptions on
the observation matrix and features play a decisive role in determining the corruption bound and the computational complexity of the non-convex algorithm.
Third, fast RPCA is limited to a corruption rate of 50% due to their choice
of the hard threshold, whereas our algorithm ups this rate to 90%. Fourth,
we prove that the costly projection onto factorized spaces is entirely optional
when features satisfy certain incoherence conditions. Although our algorithm
n
) and complexity of O(rn2 log( 1 ))
maintains the same corruption rate of O( r1.5
as fast RPCA, we show empirically that massive gains in accuracy and speed
can still be obtained. Besides, the transfer of coherence dependency from observation to features means that our algorithm is capable of dealing with highly
incoherent data.
Unavoidably, features adversely affect tolerance to corruption in IRPCAIHT (O( nd )) compared to its predecessor AltProj (O( nr )). This is not always
true with our algorithm in relation to fast RPCA. And when the underlying
rank is low but features are only weakly informative, i.e. r  d, which is often
the case, our tolerance to corruption is arguably better. IRPCA-IHT also has
2

a higher complexity of O((dn2 + d2 r) log( 1 )) than that of our algorithm. Although feature-free convex and non-convex algorithms have higher asymptotic
error bounds than our algorithm, we show in our experiments that this does not
translate as accuracy in reality. Our algorithm still has the best performance in
recovering accurately the low-rank part from highly corrupted matrices. This
may be attributed to the fact that our bounds are not tight. Besides, PCPF and
n3
) and O(r2 n2 log( 1 ))) than ours.
AltProj have much higher complexity (O( √

For PCPF, there does not exist any theoretical analysis under the deterministic
sparsity model. Nonetheless, we show in our experiments that our algorithm
is superior with regard to both recoverability and running time. The overall
contribution of this paper is as follows:
• A novel non-convex algorithm integrating features with informed sparsity
is proposed in order to solve RPCA problem.
• We establish theoretical guarantees of exact recovery under different assumptions regarding the incoherence of features and observation.
• Extensive experimental results on synthetic data indicate that the proposed algorithm is faster and more accurate in low-rank matrix recovery
than the compared state-of-the-art convex and non-convex methods for
RPCA (with and without features).
• Experiments on two real-world datasets, namely MNIST and Yale B database
demonstrate the practical merits of the proposed algorithm.

2

Notations

Lowercase letters denote scalars and uppercase letters denote matrices, unless
otherwise stated. Ai· and A · j represent the ith row and the j th column of A.
Projection onto support set Ω is given by ΠΩ . |A| is the element-wise absolute
value of matrix A. For norms of matrix A, kAkF is the Frobenius norm; kAk∗
is the nuclear norm; kAk2 is the largest singular value; otherwise, kAkp is the
lp -norm of vectorized A; and kAk2,∞ is the maximum of matrix row l2 -norms.
Moreover, hA, Bi represents tr(AT B) for real matrices A, B. Additionally, σi
is the ith largest singular value of a matrix.
The Euclidean metric is not applicable here because of the non-uniqueness
of the bi-factorisation L∗ = A∗ B∗T , which corresponds to a manifold rather
than a point. Hence, we define the following distance between (A, B) and any
of the optimal pair (A∗ , B∗ ) such that L∗ = A∗ B∗T :
q
(1)
d(A, B, A∗ , B∗ ) = min kA − A∗ Rk2F + kB − B∗ Rk2F ,
R

where R is an r × r orthogonal matrix.

3

3

Related Work

RPCA concerns a known observation matrix M which we are seeking to decompose into matrices L∗ , S∗ such that L∗ is low-rank and S∗ is sparse and
of arbitrary magnitude. Conceptually, it is equivalent to solving the following
optimization problem:
min rank(L) + γkSk0
L,S

subject to L + S = M,

(2)

for appropriate γ. This problem, regrettably, is NP-hard.
PCP Wright et al. [2009] replaces (2) with convex heuristics:
min kLk∗ + γkSk1
L,S

subject to L + S = M,

(3)

for some γ. In spite of the simplification, PCP can exactly recover the solution
of RPCA under the random model Candes et al. [2011] and the deterministic
model Chandrasekaran et al. [2011], Hsu et al. [2011].
If feasible feature dictionaries, X and Y, regarding row and column spaces
are available, PCPF Chiang et al. [2016] makes use of these to generalize (3) to
the below objective:
min kHk∗ + γkSk1
H,S

subject to XHYT + S = M,

(4)

for the same γ as in (3). Convergence to the RPCA solution has only been
established for the random sparsity model.
AltProj Netrapalli et al. [2014] addresses RPCA by minimizing an entirely
different objective:
min kM − L − SkF
L,S

subject to L ∈ set of low-rank matrices

(5)

S ∈ set of sparse matrices,
where the search consists of alternating non-convex projections. That is, during
each cycle, hard-thresholding takes place first to remove large entries and projection of appropriate residuals onto the set of low-rank matrices with increasing
ranks is carried out next. Exact recovery has also been established.
Fast RPCA Yi et al. [2016] follows yet another non-convex approach to solve
RPCA. After an initialization stage, fast RPCA updates bilinear factors U, V
such that L = UVT through a series of projected gradient descent and sparse
estimations, where U, V minimize the following loss:
1
1
min kUVT + S − Mk2F + kUT U − VT Vk2F ,
U,V 2
8

(6)

for U, V properly constrained. Recovery guarantee is ensured.
IRPCA-IHT Niranjan et al. [2017] includes features X, Y in an iterative
4

non-convex projection algorithm. Similar to AltProj, at each step, a new sparse
estimate is calculated from hard thresholding via a monotonically decreasing
threshold. After that, spectral hard thresholding takes place to attain the lowrank estimate. IRPCA-IHT provably converges to the solution of RPCA.
We also mention here several works of non-convex objectives Oh et al. [2015],
Shang et al. [2017], though exact recovery guarantees are lacking.

4

Problem Setup

Suppose that there is a known data matrix M ∈ Rn1 ×n2 , which can be decomposed into a low-rank component L∗ and a sparse error matrix S∗ of compatible
dimensions. Our aim is to identify these underlying matrices and hence robustly
recover the low-rank component with the help of available side information in
the form of feature matrices X and Y.
Concretely, let L∗ = U∗ Σ∗ V∗T be the singular value decomposition and
1
1
∗
P = XT U∗ Σ∗ 2 and Q∗ = YT V∗ Σ∗ 2 . S∗ follows the random sparsity model.
That is, the support of S∗ is chosen uniformly at random from the collection
of all support sets of the same size. Furthermore, let us be informed of the
proportion of non-zero entries per row and column, denoted by α. Assume that
there are also available features X ∈ Rn1 ×d1 and Y ∈ Rn2 ×d2 such that they are
feasible, i.e. col(X)⊇col(U∗ ) and col(Y)⊇col(V∗ ) where col(A) is the column
space of A and XT X = YT Y = I1 .
In this paper, we discuss robust low-rank recovery using the above mentioned
q
features and three different incoherence conditions: (i) kU∗ k2,∞ ≤ µn11r and
q
q
q
kV∗ k2,∞ ≤ µn12r ; (ii) kXk2,∞ ≤ µn2 d1 1 and kYk2,∞ ≤ µn2 d2 2 ; (iii) both (i)
and (ii), where r is the given rank of L∗ and µ1 , µ2 are constants.

5

Algorithm

We use a non-convex approach to achieve the above objective. The algorithm
consists of an initialization phase followed by a gradient descent phase. At each
stage, we keep track of the factors P, Q such that L = XPQT YT .

5.1

Hard-thresholding

We first introduce the sparse estimator via hard-thresholding which is used in
both phases. Given a threshold θ, Tθ (A) removes elements of A that are not
among the largest θ-fraction of elements in their respective rows and columns,
breaking ties arbitrarily for equal elements:
(
0
if |Aij | ≤ Aθ i· and |Aij | ≤ Aθ ·j,
Tθ (A)ij =
(7)
Aij
otherwise,
1 This

can always achieved via orthogonalisation.

5

where Aθ i·, Aθ ·j are the (n2 θ)th and (n1 θ)th largest element in absolute value
in row i and column j respectively.

5.2

Initialization

S is first initialized as S0 = Tα (M). Next, we obtain U0 Σ0 V0T as the rtruncated SVD of L0 , which is calculated via L0 = M − S0 . We can then
1
1
construct P0 = XT U0 Σ02 and Q0 = YT V0 Σ02 . Such an initialization scheme
gives P, Q the desirable properties for use in the second phase.

5.3

Gradient Descent

In case (i), we need the following sets:
P = {A ∈ Rd1 ×r |kXAk2,∞ ≤
Q = {A ∈ R

d2 ×r

r

r
|kYAk2,∞ ≤

2µ1 r
kP0 k2 },
n1

(8)

2µ1 r
kQ0 k2 }.
n2

(9)

Otherwise, we can simply take P as Rd1 ×r and Q as Rd2 ×r .
To proceed, we first regularise P0 and Q0 :
P = ΠP (P0 ), Q = ΠQ (Q0 ).

(10)

At each iteratiion, we first update S with the sparse estimator using a threshold of α + min(10α + 0.1):
S = Tα+min(10α+0.1) (M − XPQT YT ).

(11)

For P, Q, we define the following objective function
L(P, Q) =

1
1
kXPQT YT + S − Mk2F + kPT P − QT Qk2F .
2
64

(12)

P and Q are updated by minimizing the above function subject to the constraints imposed by the sets P and Q. That is,
P = ΠP (P − η∇P L),

(13)

Q = ΠQ (Q − η∇Q L),

(14)

where the step size η is determined analytically below. With properly initialized
P and Q, such an optimization design converges to P∗ and Q∗ . The procedure
is summarized in Algorithm 1.

6

Algorithm 1 Non-convex solver for robust principal component analysis with
features
Input: Observation M, features X, Y, rank r, corruption approximation α and
step size η.
Initialization:
1: S = Tα (M)
2: UΣVT = r-SVD(M − S)
1
3: P = XT UΣ 2
1
T
4: Q = Y VΣ 2
Gradient descent:
5: P = ΠP (P)
6: Q = ΠQ (Q)
7: while not converged do
8:
S = Tα+min(10α+0.1) (M − XPQT YT )
9:
P = ΠP (P − η∇P L)
10:
Q = ΠQ (Q − η∇Q L)
11: end while
Return: L = XPQT YT , S

6

Analysis

We first provide theoretical justification of our proposed approach. Then we
evaluate its computational complexity. The proofs can be found in the supplementary material.

6.1

Convergence

The initialization phase provides us with the following guarantees on P and Q.
Theorem 1. In cases (i) and (iii), if α ≤

1
16κrµ1 ,

d(P0 , Q0 , P∗ , Q∗ ) ≤ 18αrµ1
In case (ii), if α ≤

1√
,
16κµ2 d1 d2

we have

p

rκσ1∗ .

(15)

rd1 d2 κσ1∗ ,

(16)

we have

d(P0 , Q0 , P∗ , Q∗ ) ≤ 18αµ2

p

where κ is the condition number of L∗ and d is a distance metric defined in the
appendix.
1
Theorem 2. For η ≤ 192kL
, there exist constants c1 > 0, c2 > 0, c3 > 0,
0 k2
c1
c4 > 0, c5 > 0 and c6 > 0 such that, in case (i), when α ≤
3 , we have the
µ1 (κr) 2

following relationship
d(Pt , Qt , P∗ , Q∗ )2 ≤ (1 − c2 ησr∗ )t d(P0 , Q0 , P∗ , Q∗ )2 ,

7

(17)

in case (ii), when α ≤

c3

1

3

µ2 dr 2 κ 2

, we have

d(Pt , Qt , P∗ , Q∗ )2 ≤ (1 − c4 ησr∗ )t d(P0 , Q0 , P∗ , Q∗ )2 .
and in case (iii), when α ≤ c5 min( µ21dκ ,

1
3

µ1 (κr) 2

), we have

d(Pt , Qt , P∗ , Q∗ )2 ≤ (1 − c6 ησr∗ )t d(P0 , Q0 , P∗ , Q∗ )2 .

6.2

(18)

(19)

Complexity

From Theorem 2, it follows that our algorithm converges at a linear rate under
assumptions (ii) and (iii). To converge below  of the initial error, O(log( 1 ))
iterations are needed. At each iteration, the most costly step is matrix multiplication which takes O(rn2 ) time. Overall, our algorithm has total running time
of O(rn2 log( 1 )).

7

Experimental results

We have found that when the step size is set to 0.5, reasonable results can be
obtained. For all algorithms in comparison, we run a total of 3000 iterations or
until kM − L − SkF /kMkF < 10−7 is met.

7.1

Phase transition

Here, we vary the rank and the error sparsity to investigate the behavior of both
our algorithm and existing state-of-art algorithms in terms of recoverability.
True low-rank matrices are created via L∗ = JKT , where 200 × r matrices J, K
have independent elements drawn randomly from a Gaussian distribution of
mean 0 and variance 5 · 10−3 so r becomes the rank of L∗ . Next, we corrupt
each column of L∗ such that α of the elements are set independently with
r
magnitude U(0, 40
). However, this does not guarantee α row corruption. We
thus select only matrices whose maximum row corruption does not exceed α +
6.5% but we still feed α to the algorithms in order to demonstrate that our
algorithm does not need the exact value of corruption ratio. We consider two
types of signs for error: Bernoulli ±1 and sgn(L∗ ). The resulting M thus
becomes the simulated observation. In addition, let L∗ = UΣVT be the SVD
of L∗ . Feature X is formed by randomly interweaving column vectors of U with
5 arbitrary orthonormal bases for the null space of UT , while permuting the
expanded columns of V with 5 random orthonormal bases for the kernel of VT
forms feature Y. Hence, the feasibility conditions are fulfilled: col(X) ⊇col(L0 ),
col(Y) ⊇col(LT0 ). For each (r, α) pair, three observations are constructed. The
recovery is successful if for all these three problems,
kL − L∗ kF
< 10−3
kL∗ kF
8

(20)

Figure 1: Domains of recovery by various algorithms: (a) for random signs and
(b) for coherent signs.
from the recovered L.
Figures 1(I) plot results from algorithms incorporating features. Besides,
our algorithm contrasts with fast RPCA in Figure 1(II). Other feature-free algorithms are investigated in Figure 1(III). Figures 1(a) illustrate the random
sign model and Figures 1(b) for the coherent sign model. All previous nonconvex attempts fail to outperform their convex equivalents. IRPCA-IHT is

9

unable to deal with even moderate levels of corruption. The frontier of recoverability that has been advanced by our algorithm over PCPF is phenomenal,
massively ameliorating fast RPCA. The anomalous asymmetry in the two sign
models is no longer observed in non-convex algorithms.

7.2

Running Time
700
600

Time

500

(i)
PCP
PCPF
fast RPCA
our algorithm
AltProj

400

800

(ii)
IRPCA-IHT
our algorithm

600

400

300
200

200
100

0
1500 1750 2000 2250 2500 2000 4000 6000 8000 10000

Dimension

Dimension

Figure 2: (i) Running times for observation matrices of increasing dimensions
for (i) PCP, PCPF, fast RPCA, AltProj, our algorithm and (ii) IRPCA-IHT
∗
kF
and our algorithm when kL−L
≤ 1%.
kL∗ kF
Next, we highlight the speed of our algorithm for large-scale matrices, typical of video sequences Xiong et al. [2016]. 1500×1500 to 2500×2500 random
observation matrices are generated, where the rank is chosen to be 20% of the
column number and random sign error corrupts 11% of the entries, with features
X, Y having a dimension of 50% of the column number. The running times of
all algorithms except IRPCA-IHT are plotted in 2 (i) because IRPCA-IHT is
∗
kF
not able to achieve a relative error ( kL−L
kL∗ kF ) less than 1% for larger matrices.
For fair comparison, we have relaxed the rank to 0.3% of the column number
and error rate to 0.1% to compare our algorithm with IRPCA-IHT for matrices
ranging from 2000×2000 to 10000×10000. We have used features X, Y having
a dimension of 80% of the column number to speed up the process. The result
is shown in Figure 2 (ii). All times are averaged over three trials. It is evident
that, for large matrices, our algorithm overtakes all existing algorithms in terms
of speed. Note that features in PCPF even slow down the recovery process.

7.3

Image Classification

Once images are denoised, classification can be performed on them. The classification results directly reflect the image denoising ability. For a set of correlated
10

α
10
15
20
25
30
35

clean

89.65

noisy
30.45
25.1
23.15
18.65
18.6
16.95

PCP
82.75
82.95
83.5
81.35
77.95
71.2

PCPF
83.35
83.4
84
82.65
79
73.4

AltProj
81.4
81.15
79.3
74.05
71.5
67.75

IRPCA-IHT
65.2
49.65
37.8
30.35
24.1
21.05

fast RPCA
81.1
79.65
78.65
75.3
72.9
71.45

our algorithm
86.9
84.8
83.8
83.15
82.05
79.05

Table 1: Classification results obtained by a linear SVM.
α
10
15
20
25
30
35

clean

92.25

noisy
87
75.85
64.35
55.85
47.15
40.55

PCP
87.25
87.15
87.6
87
81.15
74.8

PCPF
87.3
87.4
87.55
86.95
81.55
75.7

AltProj
86.45
86.75
84.65
79.4
76.75
71

IRPCA-IHT
89.3
82.85
71.2
62.35
53.5
47.4

fast RPCA
89.25
87.2
85.55
82.65
78.3
76.75

our algorithm
90.3
89.8
88.55
87.8
85.65
85.15

Table 2: Classification results obtained by an SVM with RBF kernel.
0.95
PCP
PCPF
fast RPCA
our algorithm
AltProj
IRPCA-IHT

0.85

Relative Error

0.75

0.65

0.55

0.45

0.35

0.25
0.1

0.15

0.2

0.25

0.3

0.35

Sparsity

∗

kF
Figure 3: Relative error ( kL−L
kL∗ kF ) for sparsity values: 10%, 15%, 20%,
25%, 30%, 35%.

images, low-rank algorithms are normally used to remove noise that is sparse.
The same classifier is thus able to compare the different low-rank models.
The MNIST dataset is such an example which contains hand-written digits
divided into training and testing sets. Let the observation matrix be composed
of 2000 vectorized random images from the test set stacked column-wise. In this

11

case, the left feature obtained from the training set is also applicable to the test
set because of the Eigendigit nature. This imparts our algorithm to supervised
learning where there are clean related training samples available. The right
feature does not posses such property and is set to the identity matrix. We add
a range of sparse noise to the test set separately where the noise sets the pixel to
255. For PCPF, we take d = 300 as in Chiang et al. [2016] and for IRPCA-IHT
and our algorithm we use d = 150 instead.
The relative error between the recovered matrix by the competing algorithms
and the clean test matrix is plotted in Figure 3. Our algorithm is most accurate
in removing the added artificial noise. To evaluate how classifiers perform on the
recovered matrices, we train the linear and kernel SVM using the training set
and test the corresponding models on the recovered images. Table 1 tabulates
the linear SVM. Table 2 tabulates the kernel SVM. Both classifiers confirm
the recovery result obtained by various models corroborating our algorithm’s
pre-eminent accuracy.

7.4

Face denoising

It is common practice to decompose raw facial images as a low-rank component
for faithful face representation and a sparse component for defects. This is because the face is a convex Lambertian surface which under distant and isotropic
lighting has an underlying model that spans a 9-D linear subspace Basri and
Jacobs [2003], but theoretical lighting conditions cannot be realised and there
are unavoidable occlusion and albedo variations in real images. We demonstrate
that there can be a substantial boost to the performance of facial denoising by
leveraging dictionaries learnt from the images themselves.
The extended Yale B database is used as our observation which consists
images under different illuminations for a fixed pose. We study all 64 images
of a randomly chosen person. A 32556 × 64 observation matrix is formed by
vectorizing each 168 × 192 image. For fast RPCA and our algorithm, a sparsity
of 0.2 is adopted. We learn the feature dictionary as in Xue et al. [2017]. In a
nutshell, the feature learning process can be treated as a sparse encoding problem. More specifically, we simultaneously seek a dictionary D ∈ Rn1 ×c and a
sparse representation B ∈ Rc×n2 such that:
minimize
D,B

kM − DBk2F

(21)

subject to γi ≤ t for i = 1 . . . n2 ,
where c is the number of atoms, γi ’s count the number of non-zero elements in
each sparsity code and t is the sparsity constraint factor. This can be solved
by the K-SVD algorithm. Here, feature X is the dictionary D, feature Y corresponds to a similar solution using the transpose of the observation matrix as
input. We set c to 40, t to 40 and used 10 iterations.
As a visual illustration, recovered images from all algorithms are exhibited in
Figure 4. For this challenging scenario, our algorithm totally removed all shadows. PCPF is smoother than PCP but still suffers from shade. AltProj and
12

(i)

(ii)

(iii)

(iv)

(v)

(vi)

(vii)

Figure 4: (i) original; (ii) PCPF; (iii) our algorithm; (iv) IRPCA-IHT; (v) PCP;
(vi) fast RPCA; (vii) AltProj.
fast RPCA both introduced extra artefacts. Although IRPCA-IHT managed to
remove the shadows but brought back a severely distorted image. To quantitatively verify the improvement made by our proposed method, we examine the
structural information contained within the denoised eigenfaces. Singular values
of the recovered low-rank matrices from all algorithms are plotted in Figure 5.
All non-convex algorithms are competent in incorporating the rank information
to keep only 9 singular values, vastly outperforming convex approaches. Among
them, our algorithm has the most rapid decay that is found naturally Wright
et al. [2011].

13

PCP
PCPF
fast RPCA
AltProj
IRPCA-IHT
our method

5

3

log 10

i

4

2

1

0
1

5

9

13

17

21

25

29

i

Figure 5: Log-scale singular values of the denoised matrices.

8

Conclusion

This work proposes a new non-convex algorithm to solve RPCA with the help
of features when the error sparsity is roughly known. Exact recovery guarantee
has been established for three different assumptions about the incoherence conditions on features and the data observation matrix. Simulation experiments
suggest that our algorithm is able to recover matrices of higher ranks corrupted
by errors of higher sparsity than previous state-of-the-art approaches. Large
synthetic matrices also show that our algorithm scales best with observation
matrix dimension. MNIST and Yale B datasets further justify that our algorithm leads other approaches by a fair margin. Future work may involve finding
a more accurate initialization scheme.

14

Appendices
A

Convex Projection

Given P, the problem of finding ΠP (P) can be seen as projection
T onto
T the
intersection of a series of closed convex
sets
P
,
that
is
P
=
P
·
·
·
Pd1 ,
i
1
q

where Pi = {A ∈ Rd1 ×r ||Xi· A|2 ≤ 2µn11 r kP0 k2 }. We have emperically found
that the Cyclic Dykstra algorithm Reich and Zaslavski [2012] has the fastest
rate of convergence. Let A0 = P, and B−(d1 −1) = B−(d1 −2) = · · · = B−1 =
B0 = 0 ∈ Rd1 ×r , the Cyclic Dykstra algorithm updates, at each iteration,
Ak+1 = ΠPk+1 mod d1 (Ak + Bk+1−d1 ) and Bk+1 = Ak + Bk+1−d1 − Ak+1 .
For ΠPi (P), we formulate the equivalent optimisation problem below
r
2µ1 r
min kA − Pk2F s.t. |Xi· A|2 =
kP0 k2 ,
A
n1
q
for |Xi· P|2 > 2µn11 r kP0 k2 . Its solution is given by
|Xi· P|2
( q 2µ
1r

A = (Id1 ×d1 +

n1

kP0 k2

(22)

− 1)XTi· Xi·

|Xi· |22

)−1 P.

(23)

For Q, ΠQ (Q) follows similarly.
We have also run experiments to see how much improvement can be gained by
convex projection. 200×200 high-incoherence matrices are created with ranks
from 140 to 155 and corrupted by 10% random sign errors. Our algorithm is
applied with projection several times. Each uses a different number of iterative
steps ranging from 0 to 2000. Recoverability is plotted against the number of
iterative projections in Figure 6. There is hardly any noticeable improvement
so we do not use convex projection in our comparison experiments. Further
analysis is demanded to justify the redundency of convex projection.

15

151
149

rank

147
145
143
141
0

400

800

1200

1600

2000

iterations

None

our algorithm

Figure 6: Effectiveness of convex projection.

B

Proofs

For simplicity, we assume that n1 = n2 = n, d1 = d2 = d.

B.1

Proof of Theorem 1

We first declare some lemmas that will be essential to our result.
Lemma B.1. Let S0 be obtained from the initialisation phase, we have
kM − S0 − L∗ k∞ ≤ 2kL∗ k∞ .

(24)

Proof. See Yi et al. [2016] theorem 1.
Lemma B.2. For any matrix A ∈ Rn×n for which the proportion of non-zero
entries per row and column is β, we have
kAk2 ≤ βnkAk∞ .

(25)

Proof. See Netrapalli et al. [2014] lemma 4.
Lemma B.3. For two rank r matrices L1 and L2 of the same dimension whose
compact SVDs are L1 = U1 Σ1 V1T and L2 = U2 Σ2 V2T , we have
1

1

1

1

d(U1 Σ12 , V1 Σ12 , U2 Σ22 , V2 Σ22 )2 ≤ √
16

2
kL1 − L2 k2F
,
2 − 1 σr (L2 )

(26)

provided kL1 − L2 k2 ≤ 12 σr (L2 ).
Proof. See Tu et al. [2016] lemma 5.14.
Lemma B.4. For any matrices A and B of consistent sizes, we have
kABk2,∞ ≤ kAk2,∞ kBk2 .

(27)

Proof. See Liu et al. [2017] lemma 4.2.
Lemma B.5. For any matrix A with compact SVD A = UΣVT ,
kAk∞ ≤ kΣk2 kUk2,∞ kVk2,∞ .

(28)

Proof. See Yi et al. [2016] theorem 1.
Lemma B.6. Let U0 , V0 , Σ0 , S0 be obtained from the initialisation phase, we
have
kU0 Σ0 V0T − M + S0 k2 ≤ kM − S0 − L∗ k2
(29)
Proof. Weyl’s theorem tells us that, for 1 ≤ i ≤ n, |σi (L∗ )−σi (M−S0 )| ≤ kM−
S0 −L∗ k2 . When i = r +1, σi (L∗ ) = 0 and σi (M−S0 ) = kU0 Σ0 V0T −M+S0 k2
because L∗ has rank r and U0 Σ0 V0T = r-SVD(M − S0 ).
Lemma B.7. For A, B, C, D ∈ Rd×r
d(XT A, YT B, XT C, YT D) ≤ d(A, B, C, D).

(30)

Proof.
d(XT A, YT B, XT C, YT D)
q
= min kXT (A − CR)k2F + kYT (B − DR)k2F
R
q
≤ min kXT k22 k(A − CR)k2F + kYT k22 k(B − DR)k2F
R
q
= min k(A − CR)k2F + k(B − DR)k2F

(31)

R

= d(A, B, C, D).

We begin by deriving a bound on kM − S0 − L∗ k2 ,
kM − S0 − L∗ k2 ≤ 2αnkM − S0 − L∗ k∞ ≤ 4αnkL∗ k∞
≤ 4αnkΣ∗ k2 kU∗ k2,∞ kV∗ k2,∞ ,

17

(32)

where the first inequality follows from Lemma B.2 with β = 2α, the second from Lemma B.1 and the third from Lemma B.5. Next, we look at
kU0 Σ0 V0T − L∗ k2 :
kU0 Σ0 V0T − L∗ k2
≤ kU0 Σ0 V0T − M + S0 k2 + kM − S0 − L∗ k2

(33)

≤ 2kM − S0 − L∗ k2
≤ 8αnkΣ∗ k2 kU∗ k2,∞ kV∗ k2,∞ ,

where we have used Lemma B.6 and (32).
1
In cases (i) and (iii), the condition α ≤ 16κµ
gives kU0 Σ0 V0T −L∗ k2 ≤ 12 σr (L∗ )
1r
and we have
d(P0 , Q0 , P∗ , Q∗ )2
1

1

1

1

= d(XT U0 Σ02 , YT V0 Σ02 , XT U∗ Σ∗ 2 , YT V∗ Σ∗ 2 )2
1

1

1

1

≤ d(U0 Σ02 , V0 Σ02 , U∗ Σ∗ 2 , V∗ Σ∗ 2 )2

(34)

kU0 Σ0 V0T − L∗ k2F
2
σr (L∗ )
2−1
128r3 α2 κσ1∗ µ21
2r kU0 Σ0 V0T − L∗ k22
√
≤√
≤
,
∗
σr (L )
2−1
2−1
≤√

using Lemma B.7, Lemma B.3 and (33). So, we have
p
d(P0 , Q0 , P∗ , Q∗ ) ≤ 18µ1 αr rκσ1∗ .

(35)

In case (ii), we have
r
∗

T

∗

T

∗

∗

T

∗

T

∗

kU k2,∞ = kXX U k2,∞ ≤ kX U k2 kXk2,∞ ≤
r
kV k2,∞ = kYY V k2,∞ ≤ kY V k2 kYk2,∞ ≤
The condition α ≤
similar to (35)

1
16κµ2 d

gives kU0 Σ0 V0T − L∗ k2 ≤

µ2 d
,
n

(36)

µ2 d
.
n

(37)

1
∗
2 σr (L )

and we have

p
d(P0 , Q0 , P∗ , Q∗ ) ≤ 18µ2 αd rκσ1∗ .

B.2

(38)

Proof of Theorem 2

To ease our exposition, we define the following auxiliary quantities.
Let the solution set be
E = {(A, B) ∈ Rd×r × Rd×r |d(A, B, P∗ , Q∗ ) = 0}.

18

(39)

For any (P, Q) ∈ Rd×r × Rd×r , the corresponding solution is given by
(P† , Q† ) ∈ arg min kP − Ak2F + kQ − Bk2F .

(40)

(A,B)∈E

Let ∆P = P − P† , ∆Q = Q − Q† and δ = k∆Pk2F + k∆Qk2F , from which we
have
2k∆PkF k∆QkF ≤ k∆Pk2F + k∆Qk2F ,
√
k∆PkF + k∆QkF ≤ 2δ,
√
k∆Pk2F + k∆QkF k∆PkF ≤ 2δk∆PkF ,
(41)
√
k∆Qk2F + k∆QkF k∆PkF ≤ 2δk∆QkF ,
4k∆QkF k∆PkF ≤ k∆Pk2F + k∆Qk2F + 2k∆QkF k∆PkF
√
≤ 2δ(k∆QkF + k∆PkF ).
Let H = 21 kXPQT YT + S − Mk2F and ∆M = ∇L H(P, Q), we have
∇L H(P, Q) = XPQT YT + S − M = L + S − L∗ − S∗ .

(42)

We also have
∇P H(P, Q) = XT ∇L HYQ,
T

T

∇Q H(P, Q) = (X ∇L HY) P.
Let G(P, Q) =

1
T
64 kP P

(43)
(44)

− QT Qk2F , we have

1
P(PT P − QT Q),
16
1
∇Q G(P, Q) =
Q(QT Q − PT P).
16
 
 †
P
P
Let F =
, F† =
and ∆F = F − F† , then we have δ = k∆Fk2F .
Q
Q†
∇P G(P, Q) =

(45)
(46)

We now state several lemmas that will help us construct the proof.
Lemma B.8. For any P ∈ Rd×r and Q ∈ Rd×r , we have
√
p
2δ 2
∗ 2
∗
kL − L kF ≤ 2δ( σ1 +
) .
4
Proof.

(47)

kL − L∗ kF = kXPQT YT − XP† Q†T YT kF
= kX(P† ∆QT + ∆PQ†T + ∆P∆QT )YT kF
≤ kP† ∆QT kF + k∆PQ†T kF + k∆P∆QT kF
≤ kP† k2 k∆QkF + k∆PkF kQ† k2 + k∆PkF k∆QkF
√
p
p
2δ
∗
∗
≤ σ1 k∆QkF + σ1 k∆PkF +
(k∆QkF + k∆PkF )
4
√
√
p
p
2δ
2δ √
∗
∗
≤ ( σ1 +
)(k∆QkF + k∆PkF ) ≤ ( σ1 +
) 2δ.
4
4
19

(48)

q

3µ rσ ∗

1
1
Lemma B.9. For 1 ≤ i, j ≤ n, in case (i), if kXPk2,∞ ≤
and
2n
q
∗
3µ1 rσ1
kYQk2,∞ ≤
2n , then
r
r
3
1 µ1 rσ1∗
∗
|(L − L )ij | ≤
(3 +
)(k(X∆P)i· k2 + k(Y∆Q)j· k2 )
(49)
2
n
2
q
q
3µ2 dσ1∗
3µ2 dσ1∗
and
kYQk
≤
and in cases (ii) and (iii), if kXPk2,∞ ≤
2,∞
2n
2n ,
then
r
r
1 µ2 dσ1∗
3
∗
|(L − L )ij | ≤
(3 +
)(k(X∆P)i· k2 + k(Y∆Q)j· k2 ).
(50)
2
n
2

Proof.
|(L − L∗ )ij |
= |(XPQT YT − XP† Q†T YT )ij |
≤ |(XP† )i· (Y∆Q)T·j | + |(X∆P)i· (YQ† )T·j | + |(X∆P)i· (Y∆Q)T·j |
≤ k(XP† )i· k2 k(Y∆Q)j· k2 + k(X∆P)i· k2 k(YQ† )j· k2 + k(X∆P)i· k2 k(Y∆Q)j· k2
≤ kXP† k2,∞ k(Y∆Q)j· k2 + kYQ† k2,∞ k(X∆P)i· k2
1
1
+ kX∆Pk2,∞ k(Y∆Q)j· k2 + kY∆Qk2,∞ k(X∆P)i· k2
2
2
1
†
≤ ((3kXP k2,∞ + kXPk2,∞ )k(Y∆Q)j· k2
2
+ (3kYQ† k2,∞ + kYQk2,∞ )k(X∆P)i· k2 ),
(51)
where we have used kX∆Pk2,∞ ≤ kXP† k2,∞ + kXPk2,∞ and kY∆Qk2,∞ ≤
kYQ† k2,∞ + kYQk2,∞ .
In case (i),
r
r
1
µ1 rσ1∗
3µ1 rσ1∗
|(L − L )ij | ≤ (3
+
)(k(X∆P)i· k2 + k(Y∆Q)j· k2 )
2
n
2n
r
r
1 µ1 rσ1∗
3
=
(3 +
)(k(X∆P)i· k2 + k(Y∆Q)j· k2 ).
2
n
2
∗

In cases (ii) and (iii),
r
r
1
µ2 dσ1∗
3µ2 dσ1∗
∗
+
)(k(X∆P)i· k2 + k(Y∆Q)j· k2 )
|(L − L )ij | ≤ (3
2
n
2n
r
r
1 µ2 dσ1∗
3
=
(3 +
)(k(X∆P)i· k2 + k(Y∆Q)j· k2 ).
2
n
2

20

(52)

(53)

Lemma B.10. For any α ∈ (0, 1), suppose the support index set Ω ⊆ [n] × [n]
satisfies |Ωi· | ≤ αn for all i ∈ [n] and |Ω·j | ≤ αn for all j ∈ [n] where Ωi· =
{(i, j) ∈ Ω|j ∈ [n]} and Ω·j = {(i, j) ∈ Ω|i ∈ [n]}. In case (i), we have

∗

kΠΩ (L − L

)k2F

αµ1 rσ1∗
(3 +
≤
2

r

3 2
) (k∆Pk2F + k∆Qk2F ),
2

(54)

3 2
) (k∆Pk2F + k∆Qk2F ).
2

(55)

|(L − L∗ )ij |2 .

(56)

and in cases (ii) and (iii), we have
∗

kΠΩ (L − L

)k2F

αµ2 dσ1∗
≤
(3 +
2

r

Proof.
kΠΩ (L − L∗ )k2F =

X
i,j∈Ω

Using Lemma B.9, in case (i),
X
|(L − L∗ )ij |2
i,j∈Ω

X µ1 rσ ∗
1
≤
(3 +
4n
i,j∈Ω

r

3 2
) (k(X∆P)i· k2 + k(Y∆Q)j· k2 )2
2

r
X µ1 rσ ∗
3 2
1
(3 +
) (k(X∆P)i· k22 + k(Y∆Q)j· k22 )
≤
2n
2
i,j∈Ω
r
X
µ1 rσ1∗
3 2 X
≤
(3 +
) (
k(X∆P)i· k22 +
k(Y∆Q)j· k22 )
2n
2
i,j∈Ω
i,j∈Ω
r
∗
X
X
X X
µ1 rσ1
3 2
≤
(3 +
) (
k(X∆P)i· k22 +
k(Y∆Q)j· k22 )
2n
2
i j∈Ωi·
j i∈Ω·j
r
∗
3 2
αµ1 rσ1
≤
(3 +
) (kX∆Pk2F + kY∆Qk2F )
2
2
r
αµ1 rσ1∗
3 2
≤
(3 +
) (k∆Pk2F + k∆Qk2F ).
2
2

21

(57)

and in cases (ii) and (iii),
X
|(L − L∗ )ij |2
i,j∈Ω

X µ2 dσ ∗
1
≤
(3 +
4n
i,j∈Ω

r

3 2
) (k(X∆P)i· k2 + k(Y∆Q)j· k2 )2
2

r
X µ2 dσ ∗
3 2
1
(3 +
) (k(X∆P)i· k22 + k(Y∆Q)j· k22 )
2n
2
i,j∈Ω
r
X
µ2 dσ1∗
3 2 X
≤
(3 +
) (
k(X∆P)i· k22 +
k(Y∆Q)j· k22 )
2n
2
i,j∈Ω
i,j∈Ω
r
∗
X
X
X X
µ2 dσ1
3 2
(3 +
) (
k(X∆P)i· k22 +
k(Y∆Q)j· k22 )
≤
2n
2
i j∈Ωi·
j i∈Ω·j
r
∗
αµ2 dσ1
3 2
≤
(3 +
) (kX∆Pk2F + kY∆Qk2F )
2
2
r
αµ2 dσ1∗
3 2
≤
(3 +
) (k∆Pk2F + k∆Qk2F ).
2
2
≤

(58)

Lemma B.11. Given that S = Tα+min(10α,0.1) (M − XPQT YT ), we have in
case (i)
hXT ∇L H(P, Q)Y, PQT − P† Q†T + ∆P∆QT i ≥ kL − L∗ k2F
r
√
p
µ1 rσ1∗ δ
3 2
2αδ
2δ 2
∗
−
((4+β)α+2 min(10α, 0.1))(3+
) −
( σ1 +
)
4
2
β min(10α, 0.1)
4
q
√
α
√
2 + 2 min(10α,0.1)
√ p
2δ
∗
3
−
δ ( σ1 +
), (59)
2
4
and in cases (ii) and (iii)
hXT ∇L H(P, Q)Y, PQT − P† Q†T + ∆P∆QT i ≥ kL − L∗ k2F
r
√
p
3 2
2αδ
2δ 2
µ2 dσ1∗ δ
∗
((4+β)α+2 min(10α, 0.1))(3+
) −
( σ1 +
)
−
4
2
β min(10α, 0.1)
4
q
√
α
√
2 + 2 min(10α,0.1)
√ p
2δ
∗
3
−
δ ( σ1 +
). (60)
2
4
Proof.
hXT ∇L H(P, Q)Y, PQT − P† Q†T + ∆P∆QT i
= hL + S − L∗ − S∗ , L − L∗ + X∆P∆QT YT i
≥ kL − L∗ k2F − |hS − S∗ , L − L∗ i| − |hL + S − L∗ − S∗ , X∆P∆QT YT i|.
(61)
22

Following Yi et al. [2016] lemma 2, we have
β
α
|hS−S∗ , L−L∗ i| ≤ kΠΩ (L−L∗ )k2F +(1+ )kΠΩ∗ \Ω (L−L∗ )k2F +
kL−L∗ k2F ,
2
β min(10α, 0.1)
(62)
where β > 0, Ω and Ω∗ are supports of S and S∗ respectively.
On the other hand,
|hL + S − L∗ − S∗ , X∆P∆QT YT i| ≤ |hΠΩ∗c ∩Ωc (L − L∗ ), X∆P∆QT YT i|
+ |hΠΩ∗ ∩Ωc (∆M), X∆P∆QT YT i|, (63)
because ∆M has support Ωc . From Cauchy-Swartz inequality, we have
|hΠΩ∗c ∩Ωc (L − L∗ ), X∆P∆QT YT i|
≤ kΠΩ∗c ∩Ωc (L − L∗ )kF kX∆P∆QT YT kF
≤ kL − L∗ kF kX∆P∆QT YT kF

(64)

≤ kL − L∗ kF k∆PkF k∆QkF
δ
≤ kL − L∗ kF .
2
From Yi et al. [2016] lemma 2, we have
s

2α
kL − L∗ kF kX∆P∆QT YT kF
min(10α, 0.1)
r
α
kL − L∗ kF .
≤δ
2 min(10α, 0.1)
(65)

|hΠΩ∗ ∩Ωc (∆M), X∆P∆QT YT i| ≤

So,
r
δ
α
kL − L∗ kF + δ
kL − L∗ kF
2
2 min(10α, 0.1)
s
δ
2α
)kL − L∗ kF .
≤ (1 +
2
min(10α, 0.1)
(66)

|hL + S − L∗ − S∗ , X∆P∆QT YT i| ≤

Together, we have
hXT ∇L H(P, Q)Y, PQT −P† Q†T +∆P∆QT i ≥ kL−L∗ k2F −kΠΩ (L−L∗ )k2F
s
β
α
δ
2α
∗ 2
∗ 2
−(1+ )kΠΩ∗ \Ω (L−L )kF −
kL−L kF − (1+
)kL−L∗ kF .
2
β min(10α, 0.1)
2
min(10α, 0.1)
(67)
23

From Lemma B.8, we have
hXT ∇L H(P, Q)Y, PQT −P† Q†T +∆P∆QT i ≥ kL−L∗ k2F −kΠΩ (L−L∗ )k2F
q
√
α
√
√
2 + 2 min(10α,0.1)
√ p
p
2δ
2δ
β
2αδ
∗ 2
2
∗
∗
3
−(1+ )kΠΩ∗ \Ω (L−L )kF −
( σ1 +
) −
δ ( σ1 +
).
2
β min(10α, 0.1)
4
2
4
(68)
Since ΠΩ (L − L∗ ) and ΠΩ∗ \Ω (L − L∗ ) have at most α + min(10α, 0.1)-fraction
and α-fraction non-zero entries per row and column respectively, from Lemma
B.10, we have in case (i)
hXT ∇L H(P, Q)Y, PQT − P† Q†T + ∆P∆QT i
r
3 2
αµ1 rσ1∗
∗ 2
(3 +
) (k∆Pk2F + k∆Qk2F )
≥ kL − L kF −
2
2
r
min(10α, 0.1)µ1 rσ1∗
3 2
−
(3 +
) (k∆Pk2F + k∆Qk2F )
2
2
r
αµ1 rσ1∗
β
3 2
−
(1 + )(3 +
) (k∆Pk2F + k∆Qk2F )
2
2
2
q
√
α
√
√
2
+
2
p
p
min(10α,0.1) √
2αδ
2δ 2
2δ
∗
∗
−
( σ1 +
) −
δ 3 ( σ1 +
)
β min(10α, 0.1)
4
2
4
r
µ1 rσ1∗ δ
3 2
∗ 2
≥ kL − L kF −
((4 + β)α + 2 min(10α, 0.1))(3 +
)
4
2
q
√
α
√
√
2 + 2 min(10α,0.1)
√ p
p
2αδ
2δ 2
2δ
∗
∗
( σ1 +
) −
),
−
δ 3 ( σ1 +
β min(10α, 0.1)
4
2
4
(69)

24

and in cases (ii) and (iii)
hXT ∇L H(P, Q)Y, PQT − P† Q†T + ∆P∆QT i
r
αµ2 dσ1∗
3 2
∗ 2
≥ kL − L kF −
(3 +
) (k∆Pk2F + k∆Qk2F )
2
2
r
min(10α, 0.1)µ2 dσ1∗
3 2
(3 +
) (k∆Pk2F + k∆Qk2F )
−
2
2
r
αµ2 dσ1∗
β
3 2
−
(1 + )(3 +
) (k∆Pk2F + k∆Qk2F )
2
2
2
q
√
α
√
√
2
+
2
p
p
min(10α,0.1) √
2δ 2
2δ
2 min(10α, 0.1)αδ
δ 3 ( σ1∗ +
−
( σ1∗ +
) −
)
β
4
2
4
r
3 2
µ2 dσ1∗ δ
∗ 2
((4 + β)α + 2 min(10α, 0.1))(3 +
)
≥ kL − L kF −
4
2
q
√
α
√
√
2 + 2 min(10α,0.1)
√ p
2 min(10α, 0.1)αδ p ∗
2δ 2
2δ
∗
−
( σ1 +
) −
δ 3 ( σ1 +
).
β
4
2
4
(70)

Lemma B.12. When kF−F† k2 ≤
q ∗
3σ1
2 we have

√

2σr∗ , given that kPk2 ≤

q

3σ1∗
2

and kQk2 ≤

1
kPT P − QT Qk2F
64 √
√
1
2 + 3p ∗ 3
1 p ∗
2
∗ 2
σ1 δ . (71)
+ (2 σr δ − δ) − kL − L kF −
64
16
32

h∇P G(P, Q), P − P† i + h∇Q G(P, Q), Q − Q† i ≥

Proof.
1

1

P†T P† = (XT U∗ Σ∗ 2 R)T (XT U∗ Σ∗ 2 R)
1

1

= RT Σ∗ 2 T U∗T XXT U∗ Σ∗ 2 R
1

1

= RT Σ∗ 2 T U∗T U∗ Σ∗ 2 R
1

1

= RT Σ∗ 2 T Σ∗ 2 R
1

(72)

1

= RT Σ∗ 2 T V∗T V∗ Σ∗ 2 R
1

1

= RT Σ∗ 2 T V∗T YYT V∗ Σ∗ 2 R
1

1

= (YT V∗ Σ∗ 2 R)T (YT V∗ Σ∗ 2 R)
= Q†T Q† .

25

Then, following Yi et al. [2016] lemma 3, we have
h∇P G(P, Q), P − P† i + h∇Q G(P, Q), Q − Q† i
1
1
=
kPT P − QT Qk2F + hPT P − QT Q, ∆PT ∆P − ∆QT ∆Qi. (73)
32
32

≤
≤
≤
≤
≤
≤
≤
≤
≤
≤

1
hPT P − QT Q, ∆PT ∆P − ∆QT ∆Qi
32
1
|hPT P − QT Q, ∆PT ∆P − ∆QT ∆Qi|
32
1
kPT P − QT QkF k∆PT ∆P − ∆QT ∆QkF
32
1
kPT P − QT QkF (k∆Pk2F + k∆Qk2F )
32
1
kPT P − P†T P† + Q†T Q† − QT QkF δ
32
1
(kPT P − P†T P† kF + kQT Q − Q†T Q† kF )δ
32
1
(kPT P − PT P† + PT P† − P†T P† kF + kQT Q − QT Q† + QT Q† − Q†T Q† kF )δ
32
1
(kPT ∆P + ∆PT P† kF + kQT ∆Q + ∆QT Q† kF )δ
32
1
((kPk2 + kP† k2 )k∆PkF + (kQk2 + kQ† k2 )k∆QkF )δ
32
r
3σ1∗
1 p ∗
( σ1 +
)(k∆P kF + k∆QkF )δ
32
2
√
√
2 + 3p ∗ 3
σ1 δ .
32
(74)

Following Yi et al. [2016] lemma 3, we have
1
1
1 √
1
kPT P−QT Qk2F ≥
kPT P−QT Qk2F + ( 2k∆FF†T kF −δ)2 − kL−L∗ k2F ,
32
64
64
16
(75)
where we have used the fact that −kPQT − P† Q†T k2F ≥ −kXPQT YT −
XP† Q†T YT k2F .
 T ∗
 †
 T ∗
X U
X U
P
∗ 12
†
=
Σ R. If we let E =
, then
We know that F =
Y T V∗
Y T V∗
Q†
 T ∗
X U
ET E = [U∗T X V∗T Y ]
= 2I ∈ Rr×r . So
Y T V∗
√ 

2 X T U∗ √ ∗ 1
)( 2Σ 2 )R,
F =(
2 Y T V∗
†

26

(76)

is the SVD of F† . Therefore,
1
1
1
1 p
kPT P − QT Qk2F ≥
kPT P − QT Qk2F + (2 σr∗ δ − δ)2 − kL − L∗ k2F .
32
64
64
16
(77)
Thus, altogether we have
1
kPT P − QT Qk2F
64 √
√
1 p
1
2 + 3p ∗ 3
σ1 δ . (78)
+ (2 σr∗ δ − δ)2 − kL − L∗ k2F −
64
16
32

h∇P G(P, Q), P − P† i + h∇Q G(P, Q), Q − Q† i ≥

Lemma B.13. When S = Tα+min(10α,0.1) (M − XPQT YT ), given that kPk2 ≤
q ∗
q ∗
3σ1
3σ1
and
kQk
≤
2
2
2 , we have
s
k∇L H(P, Q)k2F ≤ (1 +

2α
)2 kL − L∗ k2F ,
min(10α, 0.1)

k∇P G(P, Q)k2F + k∇Q G(P, Q)k2F ≤

3σ1∗ T
kP P − QT Qk2F .
256

(79)

(80)

Proof.
k∇P G(P, Q)k2F + k∇Q G(P, Q)k2F
1
1
= k P(PT P − QT Q)k2F + k Q(QT Q − PT P)k2F
16
16
1
2
2
T
≤
(kPk2 + kQk2 )kP P − QT Qk2F
256
1 3σ1∗
3σ ∗
≤
(
+ 1 )kPT P − QT Qk2F
256 2
2
3σ1∗ T
kP P − QT Qk2F .
≤
256

(81)

From Lemma B.11, we have
s
k∇L H(P, Q)kF ≤ (1 +
so

s
k∇L H(P, Q)k2F ≤ (1 +

2α
)kL − L∗ kF ,
min(10α, 0.1)

(82)

2α
)2 kL − L∗ k2F .
min(10α, 0.1)

(83)

27

We conduct the proof of Theorem 2 by induction.
If α is small, then from Theorem 1 we have kU0 Σ0 V0T − L∗ k2 ≤
Weyl’s theorem, we have
r
1
3σ1∗
2
,
kU0 Σ0 k2 ≤
2
r
1
3σ1∗
2
kV0 Σ0 k2 ≤
,
2
and
r
1
1
3σ1∗
T
2
2
kP0 k2 ≤ kX U0 Σ0 k2 ≤ kXk2 kU0 Σ0 k2 ≤
,
2
r
1
1
3σ1∗
.
kQ0 k2 ≤ kYT V0 Σ02 k2 ≤ kYk2 kV0 Σ02 k2 ≤
2

1 ∗
2 σ1 .

By
(84)
(85)

(86)
(87)

In case (i), we thus have
r
kXΠP (P0 )k2,∞ ≤
r
kYΠQ (Q0 )k2,∞ ≤

2µ1 r
kP0 k2 ≤
n

r

2µ1 r
kQ0 k2 ≤
n

r

3σ1∗ µ1 r
,
n

(88)

3σ1∗ µ1 r
.
n

(89)

And it also follows that d(ΠP (Pt ), ΠQ (Qt ), P∗ , Q∗ ) ≤ d(Pt , Qt , P∗ , Q∗ ).
By definition,
kP − P † k2F ≤ δ,

(90)

1

kP − P † k2 ≤ kP − P † kF ≤ δ 2 .
q
p
1
And from Weyl’s theorem, if δ 2 ≤ ( 32 − 1) σ1∗ , we have
r
kP k2 ≤

3σ1∗
.
2

(91)

3σ1∗
.
2

(92)

Similarly, we also have
r
kQk2 ≤
In cases (ii) and (iii), we have
r
kXPk2,∞ ≤ kPk2 kXk2,∞ ≤
r
kYQk2,∞ ≤ kQk2 kYk2,∞ ≤

28

3σ1∗
×
2

r

3σ1∗
×
2

r

µ2 d
≤
n

r

µ2 d
≤
n

r

3dµσ1∗
,
2n

(93)

3dµσ1∗
.
2n

(94)

Now, we verify that δ diminishes:
δt+1 = kPt+1 − P†t+1 k2F + kQt+1 − Q†t+1 k2F
≤ kPt+1 − P†t k2F + kQt+1 − Q†t k2F
= kPt − η∇P Ht − η∇P Gt − P†t k2F + kQt − η∇Q Ht − η∇Q Gt − Q†t k2F
= δt − 2ηh∇P Ht + ∇P Gt , Pt − P†t i − 2ηh∇Q Ht + ∇Q Gt , Qt − Q†t i
+ η 2 k∇P Ht + ∇P Gt k2F + η 2 k∇Q Ht + ∇Q Gt k2F
= δt + η 2 k∇P Ht + ∇P Gt k2F + η 2 k∇Q Ht + ∇Q Gt k2F − 2ηh∇P Gt , Pt − P†t i
T
T
− 2ηh∇Q Gt , Qt − Q†t i − 2ηh∇L Ht , X(Pt QTt − P†t Q†T
t + ∆Pt ∆Qt )Y i

≤ δt − 2ηh∇P Gt , Pt − P†t i − 2ηh∇Q Gt , Qt − Q†t i
T
T
− 2ηh∇L Ht , X(Pt QTt − P†t Q†T
t + ∆Pt ∆Qt )Y i

+ η 2 (k∇P Ht kF + k∇P Gt kF )2 + η 2 (k∇Q Ht kF + k∇Q Gt kF )2
≤ δt − 2ηh∇P Gt , Pt − P†t i − 2ηh∇Q Gt , Qt − Q†t i
T
T
− 2ηh∇L Ht , X(Pt QTt − P†t Q†T
t + ∆Pt ∆Qt )Y i

+ 2η 2 (kXT ∇L Ht YQk2F + k∇P Gt k2F + k(XT ∇L Ht Y)T Pk2F + k∇Q Gt k2F )
≤ δt − 2ηh∇P Gt , Pt − P†t i − 2ηh∇Q Gt , Qt − Q†t i
T
T
− 2ηh∇L Ht , X(Pt QTt − P†t Q†T
t + ∆Pt ∆Qt )Y i

+ 2η 2 (k∇Q Gt k2F + k∇P Gt k2F )
+ 2η 2 (kQk22 kXT ∇L Ht Yk2F + kPk22 k(XT ∇L Ht Y)T k2F )
≤ δt − 2ηh∇P Gt , Pt − P†t i − 2ηh∇Q Gt , Qt − Q†t i
T
T
− 2ηh∇L Ht , X(Pt QTt − P†t Q†T
t + ∆Pt ∆Qt )Y i

+ 2η 2 (k∇Q Gt k2F + k∇P Gt k2F )
+ 2η 2 (kQk22 kXk22 kYk22 k∇L Ht k2F + kPk22 kXk22 kYk22 k∇L Ht k2F )
≤ δt − 2ηh∇P Gt , Pt − P†t i − 2ηh∇Q Gt , Qt − Q†t i
T
T
− 2ηh∇L Ht , X(Pt QTt − P†t Q†T
t + ∆Pt ∆Qt )Y i

+ 2η 2 (k∇Q Gt k2F + k∇P Gt k2F )
3σ ∗
3σ ∗
+ 2η 2 ( 1 k∇L Ht k2F + 1 k∇L Ht k2F )
2
2
≤ δt − 2ηh∇P Gt , Pt − P†t i − 2ηh∇Q Gt , Qt − Q†t i
T
T
− 2ηh∇L Ht , X(Pt QTt − P†t Q†T
t + ∆Pt ∆Qt )Y i

+ 2η 2 (k∇Q Gt k2F + k∇P Gt k2F + 3σ1∗ k∇L Ht k2F ).
(95)

29

Applying Lemma B.13, we get
δt+1 ≤ δt − 2ηh∇P Gt , Pt − P†t i − 2ηh∇Q Gt , Qt − Q†t i
T
T
− 2ηh∇L Ht , X(Pt QTt − P†t Q†T
t + ∆Pt ∆Qt )Y i
s
3σ ∗
2α
+ η 2 ( 1 kPTt Pt − QTt Qt k2F + 6(1 +
)2 σ1∗ kLt − L∗t k2F ).
128
min(10α, 0.1)
(96)

Applying Lemma B.12, we get
1
1
δt+1 ≤ δt + η( kLt − L∗t k2F − kPTt Pt − QTt Qt k2F )
8√
32
√ q
1 p
2+ 3
σ1∗ δt3 − (2 σr∗ δt − δt )2 )
+ η(
16
32
T
T
T
− 2ηh∇L Ht , X(Pt Qt − P†t Q†T
t + ∆Pt ∆Qt )Y i
s
3σ ∗
2α
+ η 2 ( 1 kPTt Pt − QTt Qt k2F + 6(1 +
)2 σ1∗ kLt − L∗t k2F ).
128
min(10α, 0.1)
(97)
Applying Lemma B.11, we have in case (i)
δt+1
1
15
≤ δt − η( kPTt Pt − QTt Qt k2F + kLt − L∗t k2F )
32
8
√
√ q
√
r
q p
p
√
2+ 3
1
α
2δt
+ η(
σ1∗ δt3 − (2 σr∗ δt − δt )2 + ( 2 + 2
) δt3 ( σ1∗ +
))
16
32
min(10α, 0.1)
4
r
√
p
3 2
4αδt
2δt 2
µ1 rσ1∗ δt
∗
((4 + β)α + 2 min(10α, 0.1))(3 +
) +
( σ1 +
) )
+ η(
2
2
β min(10α, 0.1)
4
s
3σ ∗
2α
)2 σ1∗ kLt − L∗t k2F ),
+ η 2 ( 1 kPTt Pt − QTt Qt k2F + 6(1 +
128
min(10α, 0.1)
(98)

30

and in cases (ii) and (iii)
δt+1
1
15
≤ δt − η( kPTt Pt − QTt Qt k2F + kLt − L∗t k2F )
32 √
8
√
√
r
q p
q
p
√
1
2+ 3
α
2δt
σ1∗ δt3 − (2 σr∗ δt − δt )2 + ( 2 + 2
) δt3 ( σ1∗ +
))
+ η(
16
32
min(10α, 0.1)
4
r
√
p
3 2
2δt 2
µ2 dσ1∗ δt
4αδt
+ η(
((4 + β)α + 2 min(10α, 0.1))(3 +
) +
( σ1∗ +
) )
2
2
β min(10α, 0.1)
4
s
∗
2α
2 3σ1
T
T
2
+η (
kPt Pt − Qt Qt kF + 6(1 +
)2 σ1∗ kLt − L∗t k2F ),
128
min(10α, 0.1)
(99)
If 10α < 0.1, then min(10α, 0.1) = 10α.
Therefore, we have in case (i)
δt+1
15
1
≤ δt − η( kPTt Pt − QTt Qt k2F + kLt − L∗t k2F )
32
8
r q
√
√ q
√
p
√
p
2+ 3
1
2
2δt
2
∗
3
+ η(
σ1 δt − (2 σr∗ δt − δt ) + ( 2 +
) δt3 ( σ1∗ +
))
16
32
5
4
r
√
µ1 rσ1∗ δt
3 2 2δt p ∗
2δt 2
+ η(
((4 + β)α + 20α)(3 +
) +
( σ1 +
) )
2
2
5β
4
r
∗
1 2 ∗
T
T
2
2 3σ1
kP Pt − Qt Qt kF + 6(1 +
) σ1 kLt − L∗t k2F )
+η (
128 t
5
15
1
≤ δt − η( kPTt Pt − QTt Qt k2F + kLt − L∗t k2F )
32 √
8
√
√
q
p
√
√ q p
2+ 3
2δt
1
2
∗
3
∗
3
∗
+ η(
σ1 δt − (2 σr δt − δt ) + ( 2 + 2 10) δt ( σ1 +
))
16
32 r
4
√
µ1 rσ1∗ δt
3 2 40δt p ∗
2δt 2
+ η(
((24 + β)α)(3 +
) +
( σ1 +
) )
2
2
β
4
√
3σ ∗
+ η 2 ( 1 kPTt Pt − QTt Qt k2F + 6(1 + 20)2 σ1∗ kLt − L∗t k2F ),
128
(100)

31

and in cases (ii) and (iii)
δt+1
1
15
≤ δt − η( kPTt Pt − QTt Qt k2F + kLt − L∗t k2F )
32 √
8
√
√
q
p
√
√ q p
1
2+ 3
2δt
σ1∗ δt3 − (2 σr∗ δt − δt )2 + ( 2 + 2 10) δt3 ( σ1∗ +
))
+ η(
16
32 r
4
√
µ2 dσ1∗ δt
3 2 40δt p ∗
2δt 2
+ η(
((24 + β)α)(3 +
) +
( σ1 +
) )
2
2
β
4
√
3σ ∗
+ η 2 ( 1 kPTt Pt − QTt Qt k2F + 6(1 + 20)2 σ1∗ kLt − L∗t k2F ).
128
(101)
On the other hand, we have min(10α, 0.1) = 0.1 if 10α ≥ 0.1.
Then, we have in case (i)
δt+1
15
1
≤ δt − η( kPTt Pt − QTt Qt k2F + kLt − L∗t k2F )
32 √
8
√
√
q p
q
p
√
√
1
2+ 3
2δt
2
∗
3
3
∗
∗
σ1 δt − (2 σr δt − δt ) + ( 2 + 2 10α) δt ( σ1 +
))
+ η(
16
32
4
r
√
µ1 rσ1∗ δt
3 2 40αδt p ∗
2δt 2
+ η(
((4 + β)α + 0.2)(3 +
) +
( σ1 +
) )
2
2
β
4
√
3σ ∗
+ η 2 ( 1 kPTt Pt − QTt Qt k2F + 6(1 + 20α)2 σ1∗ kLt − L∗t k2F ).
128
(102)
But α ≤ 1, so
δt+1
15
1
≤ δt − η( kPTt Pt − QTt Qt k2F + kLt − L∗t k2F )
32 √
8
√
√
q
p
√
√ q p
2+ 3
1
2δt
2
∗
3
3
∗
∗
+ η(
σ1 δt − (2 σr δt − δt ) + ( 2 + 2 10) δt ( σ1 +
))
16
32
4
r
√
µ1 rσ1∗ δt
3 2 40δt p ∗
2δt 2
+ η(
((4 + β)α + 20α)(3 +
) +
( σ1 +
) )
2
2
β
4
√
3σ ∗
+ η 2 ( 1 kPTt Pt − QTt Qt k2F + 6(1 + 20)2 σ1∗ kLt − L∗t k2F ).
128
(103)

32

And, similarly, we have in cases (ii) and (iii)
δt+1
1
15
≤ δt − η( kPTt Pt − QTt Qt k2F + kLt − L∗t k2F )
32 √
8
√
√
q
p
√
√ q p
1
2+ 3
2δt
2
∗
3
3
∗
∗
))
+ η(
σ1 δt − (2 σr δt − δt ) + ( 2 + 2 10) δt ( σ1 +
16
32
4
r
√
µ2 dσ1∗ δt
3 2 40δt p ∗
2δt 2
+ η(
((4 + β)α + 20α)(3 +
) +
( σ1 +
) )
2
2
β
4
√
3σ ∗
+ η 2 ( 1 kPTt Pt − QTt Qt k2F + 6(1 + 20)2 σ1∗ kLt − L∗t k2F ).
128
(104)
If η ≤

5
√
,
16(1+ 20)2 σ1∗

we have in case (i)

r
√
3 2 40δt p ∗
2δt 2
µ1 rασ1∗ δt
δt+1 ≤ δt + η(
(24 + β)(3 +
) +
( σ1 +
) )
2
2
β
4
√ q
√
√
√
√ q p
2+ 3
2δt
1 p ∗
2
∗
3
∗
3
+ η(
σ1 δt − (2 σr δt − δt ) + ( 2 + 2 10) δt ( σ1 +
)).
16
32
4
(105)
and in cases (ii) and (iii)
r
√
µ2 dασ1∗ δt
3 2 40δt p ∗
2δt 2
(24 + β)(3 +
) +
( σ1 +
) )
δt+1 ≤ δt + η(
2
2
β
4
√
√ q
√
√
√ q p
2+ 3
1 p
2δt
+ η(
σ1∗ δt3 − (2 σr∗ δt − δt )2 + ( 2 + 2 10) δt3 ( σ1∗ +
)).
16
32
4
(106)
If δt ≤ 2σr∗ , we have in case (i)
√
√
√
√
25 2 + 3 + 48 10 p ∗
3− 2 ∗
δt+1 ≤ δt + η(
δ t σ1 δ t −
σr δ t )
16
16
r
µ1 rσ1∗ δt α
3 2 90δt σ1∗
+ η(
(24 + β)(3 +
) +
)
2
2
β
√
√
√
√
25 2 + 3 + 48 10 p ∗
3− 2 ∗
≤ δt (1 + η(
σ1 δ t −
σr )
16
16
r
µ1 rσ1∗ α
3 2 90σ1∗
(24 + β)(3 +
) +
)),
+ η(
2
2
β

(107)

and in cases (ii) and (iii)
δt+1

√
√
√
√
25 2 + 3 + 48 10 p ∗
3− 2 ∗
≤ δt (1 + η(
σ1 δ t −
σr )
16
16
r
3 2 90σ1∗
µ2 dσ1∗ α
+ η(
(24 + β)(3 +
) +
)),
2
2
β
33

(108)

In case (i), if α ≤

1
16κrµ1 ,

we have

1

δ 2 = d(P0 , Q0 , P∗ , Q∗ ) ≤ 18αrµ1

p

rκσ1∗ ,

(109)

which leads to
δt+1

r
√
µ1 rσ1∗ α
3 2 3− 2 ∗
(24 + β)(3 +
) −
σr )
≤ δt (1 + η(
2
2
16
√
√
√
√
90σ1∗
225 2 + 9 3 + 432 10
+ η(
+
αrµ1 σ1∗ rκ)),
β
8
1
16κdµ2 ,

In case (ii), if α ≤

(110)

we have

1

δ 2 = d(P0 , Q0 , P∗ , Q∗ ) ≤ 18αdµ2

p

rκσ1∗ ,

(111)

which leads to
δt+1

r
√
µ2 dσ1∗ α
3 2 3− 2 ∗
(24 + β)(3 +
) −
σr )
≤ δt (1 + η(
2
2
16
√
√
√
√
90σ1∗
225 2 + 9 3 + 432 10
+ η(
+
αdµ2 σ1∗ rκ)),
β
8

In case (iii), if α ≤

1
16κrµ1 ,

(112)

we have

1

δ 2 = d(P0 , Q0 , P∗ , Q∗ ) ≤ 18αrµ1

p

rκσ1∗ ,

(113)

which leads to
δt+1

r
√
µ2 dσ1∗ α
3 2 3− 2 ∗
≤ δt (1 + η(
(24 + β)(3 +
) −
σr )
2
2
16
√
√
√
√
225 2 + 9 3 + 432 10
90σ1∗
+
αrµ1 σ1∗ rκ)),
+ η(
β
8

(114)

In case (i), we require that
r
√
µ1 rσ1∗ α
3 2 3 − 2 ∗ 90σ1∗
(24 + β)(3 +
) −
σr +
2
2
16
β
√
√
√
√
225 2 + 9 3 + 432 10
+
αrµ1 σ1∗ rκ ≤ 0, (115)
8
which leads to
α≤

√
3− 2
16
µ1 rκ
2 (24

+ β)(3 +

q

3 2
2)

+

90κ
β
√
√
√
√ .
225 2+9 3+432 10
κrµ
rκ
1
8

+

(116)

Since other constraints on α are milder, for β large enough, there exist c1 and
c1
c2 such that if α ≤
3 ,
µ1 (κr) 2

δt ≤ (1 − c2 ησr∗ )t δ0 .
34

(117)

In case (ii), we require that
µ2 dσ1∗ α
(24 + β)(3 +
2

r

√
3 2 3 − 2 ∗ 90σ1∗
) −
σr +
2
16
β
√
√
√
√
225 2 + 9 3 + 432 10
+
αdµ2 σ1∗ rκ ≤ 0, (118)
8

which leads to
α≤

√
3− 2
16
µ2 dκ
2 (24

+ β)(3 +

q

3 2
2)

+

90κ
β
√
√
√
√ .
225 2+9 3+432 10
rκ
κdµ
2
8

+

(119)

Since other constraints on α are milder, for β large enough, there exist c1 and
c3
c2 such that if α ≤
1
3 ,
µ2 dr 2 κ 2

δt ≤ (1 − c4 ησr∗ )t δ0 .

(120)

In case (iii), we require that
µ2 dσ1∗ α
(24 + β)(3 +
2

r

√
3 2 3 − 2 ∗ 90σ1∗
) −
σr +
2
16
β
√
√
√
√
225 2 + 9 3 + 432 10
+
αrµ1 σ1∗ rκ ≤ 0, (121)
8

which leads to
α≤

√
3− 2
16
µ2 dκ
2 (24

+ β)(3 +

q

3 2
2)

+

90κ
β
√
√
√
√ .
225 2+9 3+432 10
κrµ
rκ
1
8

+

(122)

Since other constraints on α are milder, for β large enough, there exist c5 and
1
c6 such that if α ≤ c5 min( µ21dκ ,
3 ),
µ1 (κr) 2

δt ≤ (1 − c6 ησr∗ )t δ0 .

(123)

References
R. Basri and D. W. Jacobs.
TPAMI, 25:218–233, 2003.

Lambertian reflectance and linear subspaces.

E. J. Candes, X. Li, Y. Ma, and J. Wright. Robust principal component analysis?
Journal of the ACM, 58:11:1–11:37, 2011.
V. Chandrasekaran, S. Sanghavi, P. A. Parrilo, and A. S. Willsky. Rank-sparsity
incoherence for matrix decomposition. SIAM J. Optim., 21:572–596, 2011.
35

V. Chandrasekarana and M. I. Jordan. Computational and statistical tradeoffs
via convex relaxation. PNAS, 110:1181–1190, 2013.
K. Chiang, C. Hseih, and I. Dhillon. Matrix completion with noisy side information. In NIPS, 2015.
K. Chiang, C. Hsieh, and I. Dhillon. Robust principal component analysis with
side information. In ICML, 2016.
R. Ge, J. Lee, and T. Ma. Matrix completion has no spurious local minimum.
In NIPS, 2016.
P. Gong, C. Zhang, Z. Lu, J. Huang, and J. Ye. A general iterative shrinkage and
thresholding algorithm for non-convex regularized optimization problems. In
ICML, 2013.
D. Hsu, S. M. Kakade, and T. Zhang. Robust matrix decomposition with sparse
corruptions. TIT, 2011.
J. M. Kohler and A. Lucchi. Sub-sampled cubic regularization for non-convex
optimization. In ICML, 2017.
G. Liu, Z. Lin, and Y. Yu. Robust subspace segmentation by low-rank representation. In ICML, 2010.
G. Liu, Q. Liu, and P. Li. Blessing of dimensionality: Recovering mixture data
via dictionary pursuit. TPAMI, 39:47–60, 2017.
P. Netrapalli, N. U. N, S. Sanghavi, A. Anandkumar, and P. Jain. Non-convex
robust pca. In NIPS, 2014.
U.N. Niranjan, A. Rajkumar, and T. Tulabandhula. Provable inductive robust
pca via iterative hard thresholding. In UAI, 2017.
T. Oh, Y. Tai, J. Bazin, H. Kim, and I. S. Kweon. Partial sum minimization
of singular values in robust pca: Algorithm and applications. TPAMI, 38:
744–758, 2015.
S. Reich and A. Zaslavski. Infinite Products of Operators and Their Applications.
the AMS and Bar-Ilan University, 2012.
C. Sagonas, Y. Panagakis, S. Zafeiriou, and M. Pantic. Raps: Robust and efficient automatic construction of person-specific deformable models. In CVPR,
2014.
F. Shang, J. Cheng, Y. Liu, Z. Luo, and Z. Lin. Bilinear factor matrix norm
minimization for robust pca: Algorithms and applications. TPAMI, PP, 2017.
S. Tu, R. Boczar, M. Simchowitz, M. Soltanolkotabi, and B. Recht. Low-rank
solutions of linear matrix equations via procrustes flow. In ICML, 2016.

36

J. Wright, A. Ganesh, S. Rao, and Y. Ma. Robust principal component analysis:
Exact recovery of corrupted low-rank matrices via convex optimization. In
NIPS, 2009.
J. Wright, A. Ganesh, A. Yang, Z. Zhou, and Y. Ma. Sparsity and robustness
in face recognition. arXiv:1111.1014, 2011.
H. Xiong, T. Liu, and D. Tao. Diversified dynamical gaussian process latent
variable model for video repair. In AAAI, 2016.
M. Xu, R. Jin, and Z. Zhou. Speedup matrix completion with side information:
application to multi-label learning. In NIPS, 2013.
N. Xue, Y. Panagakis, and S. Zafeiriou. Side information in robust principal
component analysis: Algorithms and applications. In ICCV, 2017.
X. Yi, D. Park, Y. Chen, and C. Caramanis. Fast algorithms for robust pca via
gradient descent. In NIPS, 2016.

37

