1

High Dimensional Low Rank plus Sparse Matrix
Decomposition

arXiv:1502.00182v3 [cs.NA] 16 Mar 2017

Mostafa Rahmani, Student Member, IEEE and George K. Atia, Member, IEEE

Abstract—This paper is concerned with the problem of lowrank plus sparse matrix decomposition for big data. Conventional
algorithms for matrix decomposition use the entire data to
extract the low-rank and sparse components, and are based
on optimization problems with complexity that scales with the
dimension of the data, which limits their scalability. Furthermore,
existing randomized approaches mostly rely on uniform random
sampling, which is quite inefficient for many real world data
matrices that exhibit additional structures (e.g. clustering). In
this paper, a scalable subspace-pursuit approach that transforms
the decomposition problem to a subspace learning problem
is proposed. The decomposition is carried out using a small
data sketch formed from sampled columns/rows. Even when
the data is sampled uniformly at random, it is shown that the
sufficient number of sampled columns/rows is roughly O(rµ),
where µ is the coherency parameter and r the rank of the lowrank component. In addition, adaptive sampling algorithms are
proposed to address the problem of column/row sampling from
structured data. We provide an analysis of the proposed method
with adaptive sampling and show that adaptive sampling makes
the required number of sampled columns/rows invariant to the
distribution of the data. The proposed approach is amenable to
online implementation and an online scheme is proposed.
Index Terms—Low-rank Matrix, Subspace Learning, Big Data,
Matrix Decomposition, Column Sampling, Sketching

I. I NTRODUCTION

S

UPPOSE we are given a data matrix D ∈ RN1 ×N2 , which
can be expressed as
D = L + S,

(1)

where L is a low rank (LR) matrix and S a sparse matrix with
arbitrary unknown support, whose entries can have arbitrarily
large magnitude. Many important applications in which the
data under study can be naturally modeled using (1) were discussed in [1]. The cutting-edge Principal Component Pursuit
approach developed in [1], [2], directly decomposes D into
its LR and sparse components by solving the convex program
min λkṠk1 + kL̇k∗
L̇,Ṡ

subject to L̇ + Ṡ = D

(2)

where k.k1 is the `1 -norm, k.k∗ is the nuclear norm and λ determines the trade-off between the sparse and LR components
[2]. The convex program (2) can precisely recover both the
LR and sparse components if the columns and rows subspace
of L are sufficiently incoherent with the standard basis and the
This material is based upon work supported by NSF CAREER Award CCF1552497 and NSF grant No. CCF-1320547.
The authors are with the Department of Electrical and Computer Engineering, University of Central Florida, Orlando, FL 32816 USA (e-mail:
mostafa@knights.ucf.edu, george.atia@ucf.edu).

non-zero elements of S are sufficiently diffused [2]. Although
the problem in (2) is convex, its computational complexity
is intolerable with large volumes of high-dimensional data.
Even the efficient iterative algorithms proposed in [3], [4] have
prohibitive computational and memory requirements in highdimensional settings.
Contributions: This paper proposes a new randomized decomposition approach, which extracts the LR component in
two consecutive steps. First, the column-space (CS) of L
is learned from a small subset of the columns of the data
matrix. Second, the representation of the columns of L with
respect to the learned CS is obtained from a small subset
of the rows. Unlike conventional decomposition that uses the
entire data, we only utilize a small data sketch, and solve two
low-dimensional optimization problems in lieu of one highdimensional matrix decomposition problem (2) resulting in
significant running time speed-ups.
To the best of our knowledge, it is shown here for the
first time that the sufficient number of randomly sampled
columns/rows scales linearly with the rank r and the coherency
parameter of L even with uniform random sampling. Also,
in contrast to the existing randomized approaches [5], [6],
which use blind uniform random sampling, we propose a
new methodology for efficient column/row sampling. When
the columns/rows of L are not distributed uniformly in the
CS/row space (RS) of L, which prevails much of the real world
data, the proposed sampling approach is shown to achieve
significant savings in data usage compared to uniform random
sampling-based methods that require remarkable portions of
the data. The analysis presented shows that the proposed
adaptive sampling procedure can make the required number
of sampled columns/rows invariant to the data distribution. In
addition, the proposed sampling algorithms can be independently used for feature selection from high-dimensional data.
In the presented approach, once the CS is learned, each
column is decomposed efficiently and independently using the
proposed randomized vector decomposition method. Unlike
most existing approaches, which are batch-based, this unique
feature enables applicability to online settings. The presented
vector decomposition method can be independently used in
many applications as an efficient vector decomposition algorithm or for efficient linear decoding [7].
A. Notation and definitions
We use bold-face upper-case letters to denote matrices and
bold-face lower-case letters to denote vectors. Given a matrix
L, kLk denotes its spectral norm, kLkF its Frobenius norm,

2

and kLk∞ the infinity norm, which is equal to the maximum
absolute value of its elements. In an N -dimensional space,
ei is the ith vector of the standard basis (i.e., the ith element
of ei is equal to one and all the other elements are equal to
zero). The notation A = [ ] denotes an empty matrix and the
matrix A = [A1 A2 ... An ] is the column-wise concatenation
of the matrices {Ai }ni=1 . Random sampling refers to sampling
without replacement.
II. BACKGROUND AND R ELATED W ORK
A. Exact LR plus sparse matrix decomposition
The incoherence of the CS and RS of L is an important requirement for the identifiability of the decompostion problem
in (1) [1], [2]. For the LR matrix L with rank r and compact
SVD L = UΣVT (where U ∈ RN1 ×r , Σ ∈ Rr×r and
V ∈ RN2 ×r ), the incoherence condition is typically defined
through the requirements [1], [2]
µr
µr
, maxkVT ei k22 ≤
maxkUT ei k22 ≤
i
i
N1
N2
r
(3)
µr
T
and kUV k∞ ≤
N2 N1
for some parameter µ that bounds the projection of the
standard basis {ei } onto the CS and RS. Other useful measures
for the coherency of subspaces are given in [8] as,
p
p
γ(U) = N1 max |U(i, j)|, γ(V) = N2 max |V(i, j)|, (4)
i,j

i,j

where γ(U) and γ(V) bound the coherency of the CS and
the RS, respectively. When some of the elements of the
orthonormal basis of a subspace are too large, the subspace is
coherent with the standard vectors.
The decomposition of a data matrix into its LR and sparse
components was analyzed in [1], [2], and sufficient conditions
for exact recovery using the convex minimization (2) were
derived. In [1], the sparsity pattern of the sparse matrix is
selected uniformly at random following the so-called Bernoulli
model to ensure that the sparse matrix is not LR with overwhelming probability. In this model, which is also used in
this paper, each element of the sparse matrix can be nonzero independently with a constant probability. Without loss
of generality (w.l.o.g.), suppose that N2 ≤ N1 . The following
lemma states the main result of [1].
Lemma 1 (Adapted from [1]). Suppose that the support set of
S follows the Bernoulli model with parameter ρ. The convex
program (2) with λ = √1N yields the exact decomposition
1
with probability at least 1 − c1 N1 −10 provided that
r ≤ ρr N2 µ−1 (log(N1 ))

−2

,

ρ ≤ ρs

(5)

where ρs , c1 and ρr are numerical constants.
The optimization problem in (2) is convex and can be solved
using standard techniques such as interior point methods [2].
Although these methods have fast convergence rates, their
usage is limited to small-size problems due to the high
complexity of computing a step direction. Similar to the
iterative shrinking algorithms for `1 -norm and nuclear norm
minimization, a family of iterative algorithms for solving the

optimization problem (2) were proposed in [3], [4]. However,
they also require working with the entire data. For example,
the algorithm in [4] requires computing the Singular Value
Decomposition (SVD) of an N1 ×N2 matrix in every iteration.
B. Randomized approaches
Owing to their inherent low-dimensional structures, robust
principal component analysis (PCA) and low rank plus sparse
matrix decomposition can be conceivably solved using small
data sketches, i.e., a small set of random observations of the
data [6], [9]–[14]. In [9], it was shown based on simple degreeof-freedom analysis that the LR and sparse components can
be precisely recovered using a small set of random linear
measurements of D. A convex program was proposed in [9]
to recover these components using random matrix embedding
with a polylogarithmic penalty factor in sample complexity,
albeit the formulation also requires solving a high-dimensional
optimization problem.
The iterative algorithms which solve (2) have complexity
O(N1 N2 r) per iteration since they compute the partial SVD
of N1 × N2 matrices [4]. To reduce complexity, GoDec [15]
uses a randomized method to efficiently compute the SVD,
and the decomposition algorithm in [16] minimizes the rank
of ΦL instead of L, where Φ is a random projection matrix.
However, these approaches do not have provable performance
guarantees and their memory requirements scale with the full
data dimensions. Another limitation of the algorithm in [16]
is its instability since different random projections may yield
different results.
The divide-and-conquer approach in [5] (and a similar
algorithm in [17]), can achieve super-linear speedups over fullscale matrix decomposition. This approach forms an estimate
of L by combining two low-rank approximations obtained
from submatrices formed from sampled rows and columns of
D using the generalized Nyström method [18]. Our approach
also achieves super-linear speedups in decomposition, yet is
fundamentally different from [5] and offers several advantages. First, our approach is a subspace-pursuit approach that
focuses on subspace learning in a structure-preserving data
sketch. Once the CS is learned, each column of the data
is decomposed independently using a proposed randomized
vector decomposition algorithm. Second, unlike [5], which
is a batch approach that requires to store the entire data,
the structure of the proposed approach naturally lends itself
to online implementation (c.f. Section VI), which could be
very beneficial when the data comes in on the fly [19]–[22].
Third, while the analysis provided in [5] requires roughly
O(r2 µ2 max(N1 , N2 )) random observations to ensure exact
decomposition with high probability (whp), we show that the
order of sufficient number of random observations depends
linearly on the rank and the coherency parameter even if
uniform random sampling is used. Fourth, the structure of the
proposed approach enables us to leverage adaptive sampling
strategies for challenging and realistic scenarios in which
the columns and rows of L are not uniformly distributed
in their respective subspaces, or when the data exhibits additional structures (e.g. clustering) (c.f. Sections IV-B,V). It

3

is shown that the proposed adaptive sampling scheme can
make the number of randomly sampled columns/rows required
for exact decomposition invariant to the data distribution.
In such settings, the uniform random sampling used in [5]
requires significantly larger amounts of data to carry out the
decomposition.
III. S TRUCTURE OF THE P ROPOSED A PPROACH AND
T HEORETICAL R ESULT
In this section, the structure of the proposed randomized
decomposition method is presented. A step-by-step analysis
of the proposed approach is provided and sufficient conditions
for exact decomposition are derived. Theorem 5 stating the
main theoretical result of the paper is presented at the end of
this section. The proofs of the lemmas and the theorem are
deferred to the appendix.
Let us rewrite (1) as D = UQ + S, where Q = ΣV. The
representation matrix Q ∈ Rr×N2 is a full row rank matrix that
contains the expansion of the columns of L in the orthonormal
basis U. The first step of the proposed approach aims to learn
the CS of L using a subset of the columns of D, and in the
second step the representation matrix is obtained using a subset
of the rows of D.
Let U denote the CS of L. Fundamentally, U can be obtained
from a small subset of the columns of L. However, since we
do not have direct access to the LR matrix, a random subset
of the columns of D is first selected. Hence, the matrix of
sampled columns Ds1 can be written as Ds1 = DS1 , where
S1 ∈ RN2 ×m1 is the column sampling matrix and m1 is the
number of selected columns. The matrix of selected columns
can be written as
Ds1 = Ls1 + Ss1 ,

(6)

where Ls1 and Ss1 are its LR and sparse components, respectively. The idea is to decompose the sketch Ds1 into its LR and
sparse components to learn the CS of L from the CS of Ls1 .
Note that the columns of Ls1 are a subset of the columns of L
since Ls1 = LS1 . Should we be able to decompose Ds1 into
its exact LR and sparse components (c.f. Lemma 3), we also
need to ensure that the columns of Ls1 span U. The following
lemma establishes that a small subset of the columns of D
sampled uniformly at random contains sufficient information
(i.e., the columns of the LR component of the sampled data
span U) if the RS is incoherent.
Lemma 2. Suppose m1 columns are sampled uniformly at
random from the matrix L with rank r. If

 
3
,
(7)
m1 ≥ rγ 2 (V) max c2 log r, c3 log
δ
the selected columns of the matrix L span the CS of L with
probability at least 1 − δ for constants c2 and c3 .
Thus, if γ(V) is small (i.e., the RS is not coherent), a
small set of randomly sampled columns can span U. Based
on Lemma 2, if m1 satisfies (7), then L and Ls1 will have
the same CS whp. The following optimization problem (of

dimensionality N1 m1 ) is solved to decompose Ds1 into its
LR and sparse components.
1
min √ kṠs1 k1 + kL̇s1 k∗
N1

L̇s1 ,Ṡs1

(8)

subject to L̇s1 + Ṡs1 = Ds1 .
Thus, the columns subspace of the LR matrix can be recovered
by finding the columns subspace of Ls1 . Our next lemma
establishes that (8) yields the exact decomposition using
roughly m1 = O(µr) randomly sampled columns. To simplify
the analysis, in the following lemma it is assumed that the
CS of the LR matrix is sampled from the random orthogonal
model [23], i.e., the columns of U are selected uniformly at
random among all families of r-orthonormal vectors.
Lemma 3. Suppose the columns subspace of L is sampled
from the random orthogonal model, Ls1 has the same column
subspace of L and the support set of S follows the Bernoulli
model with parameter ρ. In addition, assume that the columns
of Ds1 were sampled uniformly at random. If
r 0
(9)
m1 ≥ µ (log N1 )2 and ρ ≤ ρs ,
ρr
then (8) yields the exact decomposition with probability at
least 1 − c8 N1−3 , where


0
c7 max(r, log N1 )
2
2
, 6γ (V), (c9 γ(V) log N1 ) (10)
µ = max
r
and c7 , c8 and c9 are constant numbers provided that N1 is
greater than the RHS of the first inequality of (9).
Therefore, according to Lemma 3 and Lemma 2, the CS of
L can be obtained using roughly O(rµ) uniformly sampled
data columns. Note that m1  N1 for high-dimensional data
as m1 scales linearly with r. Hence, the requirement that N1
is also greater than the RHS of the first inequality of (9) is by
no means restrictive and is naturally satisfied.
Define Û as an orthonormal basis for the learned CS. An
arbitrary column di of D can be written as di = Uqi + si ,
where qi and si are the corresponding columns of Q and S,
respectively. Thus, di − Uqi is a sparse vector. This suggests
that qi can be learned using the minimization
minkdi − Ûq̂i k1 ,
qˆi

(11)

where the `1 -norm is used as a surrogate for the `0 -norm to
promote a sparse solution [7], [8]. The optimization problem
(11) is similar to a system of linear equations with r unknown
variables and N1 equations. Since r  N1 , the idea is to
learn qi using only a small subset of the equations. Thus, we
propose the following vector decomposition program
minkST2 di − ST2 Ûq̂i k1 ,
qˆi

(12)

where S2 ∈ RN1 ×m2 selects m2 rows of Û (and the corresponding m2 elements of di ).
First, we have to ensure that the rank of ST2 U is equal to
the rank of U, for if q∗ is the optimal point of (12), then
Uq∗ will be the LR component of di . According to Lemma
2, m2 = O(rγ 2 (U)), is sufficient to preserve the rank of U

4

when the rows are sampled uniformly at random. In addition,
the following lemma establishes that if the rank of U is equal
to the rank of ST2 U, then the sufficient value of m2 for (12)
to yield the correct columns of Q whp is linear in r.
Lemma 4. Suppose that the rank of ST2 U is equal to the
rank of L and assume that the CS of L is sampled from the
random orthogonal model. The optimal point of (12) is equal
0
to qi with probability at least (1 − 3δ ) provided that
0.5

 ,
ρ≤
rβ c6 κ log Nδ01 + 1
 
 2rβ(β − 2) log 10 

N1
δ
(13)
c6 κ log 0 + 1 ,
m2 ≥ max
3(β − 1)2
δ
r 
N1 2 6 3
c5 (log 0 ) ,
δ
δ0
where κ = logrN1 , c5 and c6 are constant numbers and β can
be any real number greater than one.
Therefore, we can obtain the LR component of each column
using a random subset of its elements. Since (11) is an `1 -norm
minimization, we can write the representation matrix learning
problem as
minkDs2 − Ûs2 Q̇k1 ,

(14)

Q̇

where Ûs2 = ST2 Û. Thus, (14) learns Q using a subset of the
rows of D as ST2 D is the matrix formed from m2 sampled
rows of D.
As such, we solve two low-dimensional subspace pursuit
problems (8) and (14) of dimensions N1 m1 and N2 m2 ,
respectively, instead of an N1 N2 -dimensional decomposition
problem (2), and use a small random subset of the data to learn
U and Q. The table of Algorithm 1 explains the structure of
the proposed approach.
We can readily state the following theorem which establishes sufficient conditions for Algorithm 1 to yield exact
decomposition.
Theorem 5. Suppose the CS of the LR matrix is sampled from
the random orthogonal model and the support set of S follows
the Bernoulli model with parameter ρ. Also, it is assumed
that Algorithm 1 samples the columns and rows uniformly at
random. If for any small δ > 0, m1 satisfies the inequalities
0
(7) and (9), ρ satisfies inequality (13) with δ = δ/N2 , m2
0
satisfies inequality (13) with δ = δ/N2 and also
0
0
3
(15)
m2 ≥ r log N1 max c2 log r, c3 log
, ρ ≤ ρs
δ
0

0

0

where {ci }9i=1 , c2 and c3 are constant numbers, µ is equal to
(10), κ = logrN1 , and β can be any real number greater than
one, then the proposed approach (Algorithm 1) yields exact
decomposition with probability at least (1 − 5δ − 3rN1−7 −
c8 N1−3 ) provided that N1 is greater than the RHS of the first
inequality of (9).
Theorem 5 guarantees that the LR component can be
obtained using a small subset of the data. The randomized approach has two main advantages. First, it significantly reduces

the memory/storage requirements since it only uses a small
data sketch and solves two low-dimensional optimization problems versus one large problem. Second, the proposed approach
has O(max(N1 , N2 )×max(m1 , m2 )×r) per-iteration running
time complexity, which is significantly lower than O(N1 N2 r)
per iteration for full scale decomposition (2) [3], [4] implying
remarkable speedups for big data. For instance, consider U
and Q sampled from N (0, 1), r = 5, and S following the
Bernoulli model with ρ = 0.02. For values of N1 = N2 equal
to 500, 1000, 5000, 104 and 2 × 104 , if m1 = m2 = 10r,
the proposed approach yields the correct decomposition with
90, 300, 680, 1520 and 4800 - fold speedup, respectively, over
directly solving (2).
Algorithm 1 Structure of Proposed Approach
Input: Data matrix D ∈ RN1 ×N2
1. Initialization: Form column sampling matrix S1 ∈ RN2 ×m1 and row
sampling matrix S2 ∈ RN1 ×m2 .
2. CS Learning
2.1 Column sampling: Matrix S1 samples m1 columns of the given data
matrix, Ds1 = DS1 .
2.2 CS learning: Matrix L̂s1 is obtained as the LR component of Ds1 (8).
2.3 CS calculation: Matrix Û is formed as an orthonormal basis for the CS
of L̂s1 .
3. Representation Matrix Learning
3.1 Row sampling: Matrix S2 samples m2 rows of the given data matrix,
Ds2 = ST
2 D.
3.2 Matrix Q̂ is obtained as the optimal point of (14).
Output: If Û is an orthonormal basis for the learned CS and Q̂ is the obtained
representation matrix, then L̂ = ÛQ̂ is the obtained LR component.

IV. E FFICIENT C OLUMN /ROW S AMPLING
In sharp contrast to randomized algorithms for matrix
approximations rooted in numerical linear algebra (NLA) [24],
[25], which seek to compute matrix approximations from sampled data using importance sampling, in matrix decomposition
and robust PCA we do not have direct access to the LR matrix
to measure how informative particular columns/rows are. As
such, the existing randomized algorithms for matrix decomposition and robust PCA [5], [6], [10] have predominantly relied
upon uniform random sampling of columns/rows.
In Section IV-A, we briefly describe the implications of nonuniform data distribution and show that uniform random sampling may not be favorable for data matrices exhibiting some
structures that prevail much of the real datasets. In Section
IV-B, we demonstrate an efficient column sampling strategy
which will be integrated with the proposed decomposition
method. The decomposition method with efficient column/row
sampling is presented in Section V.
A. Non-uniform data distribution
When data points lie in a low-dimensional subspace, a small
subset of the points can span the subspace. However, uniform
random sampling is only effective when the data points are
distributed uniformly in the subspace. To clarify, Fig. 1 shows
two scenarios for a set of data points in a two-dimensional

5

7 confirms that the upper bound is tight by establishing a
converse. In these lemmas, it is assumed that the columns lie
in a union of linear subspaces as per the following assumption.

Rank of the Sampled Columns

Fig. 1. Data distributions in a two-dimensional subspace. The red points are
the normalized data points.
60
n=1
n=2
n=4
n = 10
n = 20
n = 30
n = 60

50
40
30
20
0

1000

2000

3000

4000

5000

6000

Number of Randomly Sampled Columns

Assumption 1. The matrix L can be represented as L =
[U1 Q1 ... Un Qn ]. The CS of {Ui ∈ RN1 ×r/n }ni=1 are
random r/n-dimensional subspaces in RN1 . The RS of
r/n-dimensional subspaces
{Qi ∈ Rr/n×ni }ni=1 are random
Pn
in {Rni }ni=1 , respectively, i=1 ni = N2 , and min ni  r/n.
i

Lemma 6. Suppose L follows Assumption 1. If the rank of L
is equal to r, then

#
"
n
X
1
N
c
rϕ
2
7
1

 ≤2
P max keTi Vk22 >
n−3
i
i
N2
n min nk
i=1
k


(16)
n
X
3r
N
2
−7
2
≤
n
P γ (V) > 20 log(max ni )
i
min ni
n i=1 i
i

Fig. 2. The rank of a set of uniformly sampled columns for different number
of clusters.

where ϕ1 =
subspace. In the left plot, the data points are distributed
uniformly at random. In this case, two randomly sampled data
points can span the subspace whp. In the right plot, 95 percent
of the data lie on a one-dimensional subspace, thus we may
not be able to capture the two-dimensional subspace from a
small random subset of the data points.
In practice, the data points in a low-dimensional subspace
may not be uniformly distributed, but rather exhibit some
additional structures. A prevailing structure in many modern
applications is clustered data [26]. For example, user ratings
for certain products (e.g. movies) in recommender systems
are not only LR due to their inherent correlations, but also
exhibit additional clustering structures owing to the similarity
of the preferences of individuals from similar backgrounds
(e.g. education, culture, or gender) [26].
To further show that uniform random sampling falls short
when the data points are not distributed uniformly in the
subspace, consider a matrix G ∈ R2000×6150 generated as
G = [G1 G2 ... Gn ]. For 1 ≤ i ≤ n2 , Gi = Ui Qi , where
r
r
200r
Ui ∈ R2000× n , Qi ∈ R n × n . For n/2 + 1 ≤ i ≤ n,
r
r
5r
Gi = Ui Qi , where Ui ∈ R2000× n , Qi ∈ R n × n . The
elements of Ui and Qi are sampled independently from a
normal N (0, 1) distribution. The parameter r is set equal to
60, thus the rank of G is equal to 60 whp. Fig. 2 illustrates
the rank of the randomly sampled columns versus the number
of sampled columns for different number of clusters n. When
n = 60, it turns out that we need to sample more than half
of the columns to span the CS, so we cannot evade highdimensionality with uniform random column/row sampling.
Thus, when the distribution of the clustered data is less
uniform, more randomly sampled columns are needed to
capture the CS. This is indeed confirmed by Lemma 2, which
established that the sufficient number of randomly sampled
columns is proportional to the RS coherency. In the following
lemmas, it is shown that the RS coherency increases if the
distribution of the columns is less uniform. Lemma 6 provides
an upper bound on the coherency of the RS and Lemma

max(r/n,log max nk )
k

r/n

.

Lemma 7. If L follows Assumption 1, the rank of L is equal
to r, r/n ≥ 18 log max ni and ni ≥ 96 nr log ni , 1 ≤ i ≤ n,
i
then



n
X
0.5
r
1
N
2 

P max kVT ei k22 <
≤2
n−5
i . (17)
i
N2
n min ni
i=1
i

Based on Lemma 6, the RS coherency of L is linear with
confirming that the RS coherency increases if the

N2
min ni ,
i

distribution of the columns within the CS is less uniform.
B. Efficient column sampling method
Column sampling is widely used for dimensionality reduction and feature selection [12], [27]. In the column sampling
problem, the LR matrix (or the matrix whose span is to be
approximated with a small set of its columns) is available.
Thus, the columns are sampled based on their importance,
measured by the so-called leverage scores [24], as opposed to
blind uniform sampling. We refer the reader to [12], [27] and
references therein for more information about efficient column
sampling methods.
Next, we present a sampling approach to be used in Section
V where the proposed decomposition algorithm with efficient
sampling is presented. The proposed sampling strategy is
inspired by the approach in [27] in the context of volume sampling. Algorithm 2 details the presented sampling procedure.
Given a matrix A with rank rA , the algorithm aims to sample
a small subset of columns that span its CS. The first column
is sampled uniformly at random or based on a judiciously
chosen probability distribution [24]. The next columns are
selected sequentially so as to maximize the novelty to the span
of the selected columns. As shown in step 2.2 of Algorithm
2, a design threshold τ is used to decide whether a given
column brings sufficient novelty to the sampled columns by

6

Algorithm 3 Decomposition with Adaptive Sampling

Fig. 3. Visualization of the matrices defined in Section V. Matrix Dw is
selected randomly or using Algorithm 3 described in Section V-B.

thresholding the `2 -norm of its projection on the complement
of the span of the sampled columns. The threshold τ is set
to zero in a noise-free setting. Once the selected columns are
believed to span the CS of A, they are removed from A. This
procedure is repeated C times (using the remaining columns).
In each time, the algorithm finds rA columns spanning the CS
of A. After every iteration, the rank of the matrix of remaining
columns is bounded above by rA . As such, the algorithm
samples approximately m1 ≈ CrA columns in total. In the
proposed decomposition method with efficient column/row
sampling (c.f. Sec. V), we set C large enough to ensure that
the selected columns form a low rank matrix.
Algorithm 2 Efficient Sampling from LR Matrices
Input: Matrix A.
1. Initialize
1.1 The parameter C is chosen as an integer greater than or equal to one.
The algorithm finds C sets of linearly dependent columns.
1.2 Set I = ∅ as the index set of the sampled columns and set v = τ ,
B = A and C = [ ].
2. Repeat C Times
2.1 Let b be a non-zero randomly sampled column from B with index ib .
Update C and I as C = [C b], I = {I , ib }.
2.2 While v ≥ τ
2.2.1 Set E = Pc B , where Pc is the projection matrix onto the complement
space of span(C).
2.2.2 Define f as the column of E with the maximum `2 -norm with index
if . Update C, I and v as C = [C f ] , I = {I , if } and v = kf k2 .
2.2 End While
2.3 Set C = [ ] and set B equal to A with the columns indexed by I set to
zero.
2. End Repeat
Output: The set I contains the indices of the selected columns.

1. Informative Column Sampling
1.1 Sample Cr r̂ rows of D uniformly at random to form Dw where r̂ is a
known upper bound on r.
1.2 Obtain L̂w via (2) as the LR component of Dw .
1.3 Apply Algorithm 2 to L̂w .
1.4 Form the matrix Ds1 from the columns of D corresponding to the
sampled columns of L̂w .
2. CS Learning
2.1 The convex program (8) is applied to Ds1 .
2.2 Obtain Û as an orthonormal basis for the CS of the calculated LR
component of Ds1 .
3. Representation Learning
3.1 Apply Algorithm 2 to ÛT and define Ûs2 as the matrix of sampled
rows.
3.2 Form the matrix Ds2 as the rows of D corresponding to the sampled
rows of Û.
3.3 Define Q̂ as the optimal point of (14).
Output: The matrix Û is the obtained basis for the column space of L and
L̂ = ÛQ̂ is the obtained LR matrix.

sense that they do not align along any specific directions,
such that Cr r rows of L sampled uniformly at random
span its RS whp, for some constant Cr . In Section V-B, we
dispense with this assumption. The proposed decomposition
algorithm consists of three steps detailed next. We support the
idea of each step with some examples and theoretical analysis.
1) Informative column sampling: The first important idea
underlying the proposed sampling approach is to start sampling along the dimension that has the better distribution.
For instance, consider an extreme scenario where only two
columns of G ∈ R1000×1000 are non-zero. In this case, with
random sampling we need to sample almost all the columns
to ensure that the sampled columns span the CS of G. But,
if the non-zero columns are non-sparse, a small subset of
randomly chosen rows of G will span its row space. As
another example, suppose that the distribution of the columns
of L follows Assumption 1, i.e., the distribution of the columns
admits a clustering structure. The following lemmas compare
the sufficient numbers of randomly sampled columns and rows
to capture the CS and RS, respectively.
Lemma 8. Suppose L follows Assumption 1. If m1 columns
of L are sampled uniformly at random with replacement, the
rank of L is equal to r, and


3
2n ξmax N2
m1 ≥ 2 +
log
(18)
ξmin
δ
min ni
i

V. P ROPOSED DECOMPOSITION ALGORITHM WITH
EFFICIENT SAMPLING

In this section, we present a modified decomposition
algorithm (Algorithm 3) that replaces uniform random
sampling with the efficient column/row sampling method
(Algorithm 2). We subsequently provide an analysis showing
that the required number of sampled columns using the
proposed approach can be invariant to the distribution of
the columns – even if their distribution is highly nonuniform. The table of Algorithm 3 details the proposed
method with efficient sampling. In this algorithm, it is
assumed that the rows of L are well distributed, in the

where
2r
δ
(19)
2r
ξmax = 10 c7 max(r/n, log max ni ) log
,
i
δ
then the sampled columns spanPthe column space of L with
n
probability at least 1 − 2δ − 2 i=1 n−3
i .
ξmin = 10 c7 max(r/n, log min ni ) log
i

Lemma 9. Suppose L follows Assumption 1 and m2 rows of
L are sampled uniformly at random with replacement. If the
rank of L is equal to r and
2r
m2 ≥ 10 c7 rϕ2 log ,
(20)
δ

7

then the sampled rows span the row space of L with
N1 )
probability at least 1 − δ − 2N1−3 , where ϕ2 = max(r,log
.
r

consisting of Cri columns sampled uniformly at random from
each submatrix {Gi }ni=1 .

The sufficient number of sampled rows indicated in Lemma
9 is roughly O(r) while the sufficient
 number of sampled

Remark 1. The obtained low rank component of Dw , L̂w , is
only used to locate the indices of the informative columns.
Hence, it is not imperative to obtain exact decomposition
of Dw and we can still locate the informative columns if
L̂w approximates Lw . However, for analysis purposes in the
following lemma it is assumed that L̂w = Lw .
If the rank of Lw is equal to r and L̂w approximates Lw
well enough, the only difference between the actual and the
imaginary sampler is that the former samples the informative
columns of Li , and the latter samples the columns corresponding to the i-th cluster uniformly at random from the columns of
Li . However, sampling random columns from a given cluster
is no different that sampling informative columns since r/n
randomly sampled columns of Li span its CS whp.

N2
columns in Lemma 8 is of order O r min
ni . Thus, the
i

sufficient number of randomly sampled rows is invariant to
the distribution of the columns, while the number of columns
grows linearly in N2 / mini ni , and in turn increases if the
distribution of the columns is less uniform.
In Algorithm 3, the columns can admit a clustering
structure. Thus, we start the sampling with row sampling. Let
r̂ denote a known upper bound on r. Such knowledge is often
available as side information depending on the particular
application. For instance, facial images under varying
illumination and facial expressions are known to lie on a
special low-dimensional subspace [28]. For visualization, Fig.
3 provides a simplified illustration of the matrices defined in
this section. We sample Cr r̂ rows of D uniformly at random.
Let Dw ∈ R(Cr r̂)×N2 denote the matrix of sampled rows.
We choose Cr sufficiently large to ensure that the non-sparse
component of Dw is a LR matrix. Define Lw , assumably with
rank r, as the LR component of Dw . If we locate a subset
of the columns of Lw that span its CS, the corresponding
columns of L would span its CS. To this end, the convex
program (2) is applied to Dw to extract its LR component
denoted L̂w . Then, Algorithm 2 is applied to L̂w to find a set
of informative columns by sampling m1 ≈ Cr r̂ columns. The
matrix Ds1 is formed using the columns of D corresponding
to the sampled columns of L̂w .
2) CS learning: Similar to the CS learning step of Algorithm 1, we obtain the CS of L by decomposing Ds1 .
Adaptive sampling makes the number of sampled columns
required to capture the CS of L invariant to the distribution
of the columns since only informative columns are sampled.
Yet, another important advantage of adaptive sampling is that
it also makes the sufficient number of sampled columns to
ensure correct decomposition of Ds1 almost invariant to the
distribution of the columns.
To show this fact, we present the following example along
with its theoretical footing. Suppose L follows Assumption
1. Thus, per Lemma 9 a small number of randomly sampled
rows can capture the RS. Accordingly, the RS of Lw (the
LR component of Dw ) is equal to the RS of L. Thus, if
Dw is decomposed correctly and Algorithm 2 is applied to
L̂w , it samples Cr nr columns from each matrix Li = Ui Qi
whp. Since Algorithm 2 is deterministic, it is hard to analyze
the decomposition algorithm with deterministic sampling and
upper bound the coherency parameters. Instead, we consider
an imaginary randomized sampling algorithm whose performance emulates that of the proposed sampling approach when
L follows Assumption 1.
Imaginary sampler: Suppose matrix G can be represented
as G = [G1 ... Gn ]. The CS of matrices {Gi }ni=1 are
independent subspaces with dimensions {ri }ni=1 , respectively.
Pn
The imaginary sampler applied to G samples C i=1 ri

The following two lemmas confirm that with adaptive column
sampling, the sufficient number of sampled columns to ensure
correct decomposition of Ds1 is invariant to the distribution
of the columns.
Lemma 10. Suppose L follows Assumption 1, the imaginary
sampler is applied to Lw to form Ds1 , Ls1 has the same CS
of L, and the support set of S follows the Bernoulli model
with parameter ρ. If the rank of L is equal to r,
r 00
m1 ≥ µ (log N1 )2 ,
(21)
ρr
and ρ ≤ ρs , then (8) yields the exact decomposition with
probability at least 1 − 2n(m1 /n)−3 − c8 N1−3 , where ϕ3 =
max(r/n,log m1 /n)
and
r/n


00
c7 max(r, log N1 )
c9
2
, c7 ϕ3 , c7 ϕ3 ( √ log N1 ) . (22)
µ = max
r
6
Lemma 11. Suppose L follows Assumption 1, the columns of
D are sampled uniformly at random to form Ds1 , Ls1 has the
same CS of L, and the support set of S follows the Bernoulli
model with parameter ρ. If the rank of L is equal to r,
r
m1 ≥ µ̇(log N1 )2
(23)
ρr
where

c7 max(r, log N1 )
2
µ̇ = max
, 6µv , µv (c9 log N1 )
,
r


max(r/n , log min ni ) (24)
c7  N2 
i
µv = ϕ1
, ϕ1 =
,
n
min nk
r/n


k

then (8) decomposes
Pn Ds1 correctly with probability at least
1 − c8 N1−3 − 2 i=1 n−3
i .
Remark 2. The sufficient value for m1 in Lemma 10 is O(r),
N2
versus O(r min
ni ) in Lemma 11. Hence, through adaptive
i

column sampling not only is the number of sampled columns
required to capture the CS of L invariant to the distribution of
the columns, but also the number of sampled columns required

8

for exact decomposition. If we sample from a dataset uniformly
at random, the resulting data sketches will have the same
structure of the data whp. For example, if most of the columns
of L are aligned along a given direction, then whp most of
the randomly sampled columns of L will be aligned along
that direction as well. By contrast, adaptive sampling samples
the informative columns regardless of the population of the
clusters, thereby balances the distribution of the sampled data.
Thus, the distribution of adaptively sampled columns is closer
to a uniform distribution. For instance, the adaptive sampler
applied to the data on the right plot of Fig. 1 samples an
equal number of data points from each cluster even though
the number of data points in one cluster is notably greater
than the other. This can also be observed by comparing the RS
coherency of matrix Ls1 formed with adaptive versus uniform
random sampling. The analysis provided in the proof of Lemma
10 and Lemma 11 shows that the RS coherency of Ls1
following adaptive and uniform sampling is roughly c7 and
N2
c7 min
ni , respectively. Thus, adaptive sampling significantly
i

improves the coherency of the matrix of sampled columns.
3) Representation matrix learning: In this step, unlike
Algorithm 1, the rows are not sampled randomly. Instead,
we leverage the information embedded in U to select the
informative rows. Algorithm 2 is applied to UT to locate
m2 ≈ Cr rows of U. Thus, we form the matrix Ds2 from
the rows of D corresponding to the selected rows of U. Then,
the representation matrix is learned as the optimal point of
(14). Subsequently, the LR matrix can be obtained from the
learned CS and the representation matrix. Since in Algorithm
3 it is assumed that the rows do not follow a clustering
structure, the analysis of this step is similar to the analysis
of the corresponding step in Algorithm 1.
We can readily state the following theorem which supports
the performance of Algorithm 3. In this theorem, it is assumed
that the imaginary sampler is used to sample the columns of
D. As mentioned earlier, the performance of the imaginary
sampler is approximately equivalent to our column sampling
procedure if L̂w = Lw , or L̂w approximates Lw well enough.

A. An alternative approach to the CS Learning Step
In this section, we present an alternative approach to the
CS learning step of Algorithm 3. We utilize the information
embedded in matrix L̂sw (the sampled columns of L̂w ) to obtain
U. In particular, if Dw is decomposed correctly, the RS of
L̂sw will be the same as that of Ls1 given that the rank of Lw
is equal to r. Let Vs1 be an orthonormal basis for the RS
of L̂sw . Thus, to learn the CS of Ds1 we only need to solve
T
minkDs1 − ÛVs1
k1 .
Û

Remark 3. The convex algorithm (2) may not always yield
accurate decomposition of Dw since structured data may not
be sufficiently incoherent, suggesting that the decomposition
step can be further improved. Let Dsw be the matrix consisting
of the columns of Dw corresponding to the columns selected
from L̂w to form L̂sw . According to our investigations, an
improved Vs1 can be obtained by applying the decomposition
algorithm presented in [29] to Dsw and use the RS of L̂sw as an
initial guess for the RS of the non-sparse component of Dsw .
Since Dsw is low-dimensional (roughly O(r) × O(r) matrix),
this extra step is a low complexity operation.
B. Column/Row sampling from sparsely corrupted data

In Algorithm 3, we assumed that the LR component of Dw
has rank r. However, if the rows are not well-distributed, a
reasonably sized random subset of the rows may not span
the RS of L. Here, we present a sampling approach which
can find the informative columns/rows even when both the
columns and the rows exhibit clustering structures such that a
small random subset of the columns/rows of L cannot span its
CS/RS. The algorithm presented in this section (Algorithm 4)
can be independently used as an efficient sampling approach
from big data. In this paper, we use Algorithm 4 to form Dw
if both the columns and rows exhibit clustering structures.
The table of Algorithm 4, Fig. 4 and its caption provide the
details of the proposed sampling approach and the definitions
of the used matrices. We start the cycle from the position
marked “I” in Fig. 4 with Dw formed according to the
initialization step of Algorithm 4. For ease of exposition,
Theorem 12. Suppose L follows Assumption 1, the support assume that L̂ = L and L̂ = L , i.e., D and D
w
w
c
c
w
c
set of S follows the Bernoulli model with parameter ρ, the are decomposed correctly. The matrix L̂s is the informative
w
rank of Lw is equal to r, the imaginary sampler is applied to columns of L̂ . Thus, the rank of L̂s is equal to the rank
w
w
Lw to locate the columns forming Ds1 , and applied to Û to of L̂ . Since L̂ = L , L̂s is a subset
of the rows of Lc .
w
w
w
w
locate the rows forming Ds2 . If m2 and ρ satisfy (13) with If the rows of L exhibit a clustering structure, it is likely
0
δ = δ/n,
that rank(L̂sw ) < rank(Lc ). Thus, rank(Lw ) < rank(Lc ).
!2
!
We continue one cycle of the algorithm by going through
r
r
r
2n
r 00
2
m1 ≥ max n
+ 2 log
, µ (log N1 ) , (25) steps 1, 2 and 3 of Fig. 4 to update Dw . Using a similar
n
δ
ρr
argument, we see that the rank of an updated Lw will be
greater than the rank of Lc . Thus, if we run more cycles
where {ci }9i=1 are constant numbers, κ = logrN1 , β any of the algorithm – each time updating Dw and Dc – the
00
real number greater than one, µ is equal to (22), and rank of Lw and Lc will increase. As detailed in the table
m1 /n)
ϕ3 = max(r/n,log
, then Algorithm 3 yields exact decom- of Algorithm 4, we stop if the dimension of the span of the
r/n
position with probability at least 1−4δ−n(m1 /n)−3 −c8 N1−3 , obtained LR component does not change in T consecutive
provided that N1 is greater than the RHS of the first inequality iterations. While there is no guarantee that the rank of Lw
will converge to r (it can converge to a value smaller than
of (25).
r), our investigations have shown that Algorithm 4 performs
quite well and the RS of Lw converges to the RS of L in few

9

Fig. 4. Visualization of Algorithm 4. We run few cycles of the algorithm
and stop when the rank of the LR component does not change over T
consecutive steps. One cycle of the algorithm starts from the point marked
“I” and proceeds as follows. I: Matrix Dw is decomposed and L̂w is the
obtained LR component of Dw . II: Algorithm 2 is applied to L̂w to select
the informative columns of L̂w . L̂sw is the matrix of columns selected from
L̂w . III: Matrix Dc is formed from the columns of D that correspond to
the columns of L̂sw . 1: Matrix Dc is decomposed and L̂c is the obtained LR
component of Dc . 2: Algorithm 2 is applied to L̂T
c to select the informative
rows of L̂c . L̂sc is the matrix of rows selected from L̂c . 3: Matrix Dw is
formed as the rows of D corresponding to the rows used to form L̂sc .

steps. We have also found that adding some randomly sampled
columns (rows) to Dc (Dw ) can effectively avert converging
to a lower dimensional subspace. For instance, some randomly
sampled columns can be added to Dc , which was obtained by
applying Algorithm 2 to L̂w .
Algorithm 4 Efficient Column/Row Sampling from Sparsely
Corrupted LR Matrices
1. Initialization
Form Dw ∈ RCr r̂×N2 by randomly choosing Cr r̂ rows of D. Initialize
k = 1 and set T equal to an integer greater than 1.
2. While k > 0
2.1 Sample the most informative columns
2.1.1 Obtain L̂w via (2) as the LR component of Dw .
2.1.2 Apply Algorithm 2 to L̂w with C = Cr .
2.1.3 Form the matrix Dc from the columns of D corresponding to the
sampled columns of L̂w .
2.2 Sample the most informative rows
2.2.1 Obtain L̂c via (2) as the LR component of Dc .
2.2.2 Apply Algorithm 2 to L̂T
c with C = Cr .
2.2.3 Form the matrix Dw from the rows of D corresponding to the sampled
rows of L̂c .
2.3 If the dimension of the RS of L̂w does not increase in T consecutive
iterations, set k = 0 to stop the algorithm.
2. End While
Output: The matrices Dw and L̂w can be used for column sampling in the
first step of the Algorithm presented in Section V.

Algorithm 4 was found to converge in a very small number
of iterations (typically less than 4). Thus, even when Algorithm
4 is used to form the matrix Dw , the order of complexity of
the proposed decomposition method with efficient column/row
sampling is O(max(N1 , N2 )r2 ).
VI. O NLINE I MPLEMENTATION
The proposed decomposition approach consists of two main
steps, namely, learning the CS of the LR component then
decomposing the columns independently. This structure lends
itself to online implementation, which could be very beneficial

in settings where the data arrives on the fly. The idea is to
first learn the CS of the LR component from a small batch
of the data and keep tracking the CS. Since the CS is being
tracked, any new data column can be decomposed based on
the updated subspace. The table of Algorithm 5 details the
proposed online matrix decomposition algorithm, where dt
denotes the tth received data column.
Algorithm 5 uses a parameter nu which determines the rate
at which the algorithm updates the CS of the LR component.
For instance, if nu = 20, then the CS is updated every 20
new data columns (step 2.2 of Algorithm 5). The parameter
nu has to be set in accordance with the rate of the change
of the subspace of the LR component; a small value for nu
is used if the subspace is changing rapidly. The parameter ns
determines the number of columns last received that are used
to update the CS. If the subspace changes rapidly, the older
columns may be less relevant to the current subspace, hence a
small value for ns is used. On the other hand, when the data
is noisy and the subspace changes at a slower rate, a larger
value for ns can lead to more accurate estimation of the CS.
Algorithm 5 Online Implementation
1. Initialization
1.1 Set the parameters nu and ns equal to integers greater than or equal to
one.
1.2 Form D0 ∈ RN1 ×(Cr r̂) as D0 = [d1 d2 ... dCr r̂ ]. Decompose D0
using (2) and obtain the CS of its LR component. Define Uo as the learned
CS, Qo the appropriate representation matrix and Ŝ the obtained sparse
component of D0 .
1.3 Apply Algorithm 2 to UT
o to construct the row sampling matrix S2 .
2. For any new data column dt do
2.1 Decompose dt as
T
minkST
2 dt − S2 Uo q̂t k1 ,

(26)

q̂t

and update
Qo ← [Qo q∗t ], Ŝ ← [Ŝ (dt − Uo q∗t )], where q∗t is the optimal point of
(26).
2.2 If the remainder of nt is equal to zero, update Uo as
u

minkDt − Ûo Qto k1 ,

(27)

Ûo

where Qto is the last ns r̂ columns of Qo and Dt is the matrix formed from
the last ns r̂ received data columns. Apply Algorithm 2 to the new UT
o to
update the row sampling matrix S2 .
2. End For
Output The matrix Ŝ as the obtained sparse matrix, L̂ = D − Ŝ as the
obtained LR matrix and Uo as the current basis for the CS of the LR
component.

A. Noisy data
In practice, noisy data can be modeled as D = L + S + N ,
where N is an additive noise component. In [30], it was shown
that the program
min λkŜk1 + kL̂k∗
L̂,Ŝ

s. t.

L̂ + Ŝ − D

F

≤ n ,

(28)

can recover the LR and sparse components with an error bound
that is proportional to the noise level. The parameter n has
to be chosen based on the noise level. This modified version
can be used in the proposed algorithms to account for the

10

noise. Similarly, to account for the noise in the representation
learning problem (14), the `1 -norm minimization problem can
be modified as follows:
minkST2 D − ST2 UQ̂ − Êk1

subject to kÊkF ≤ δn . (29)

Q̂,Ê

Ê ∈ Rm2 ×N2 is used to cancel out the effect of the noise and
the parameter δn is chosen based on the noise level [31].
VII. N UMERICAL S IMULATIONS
In this section, we present some numerical simulations to
study the performance of the proposed randomized decomposition method. First, we present a set of simulations confirming
our analysis which established that the sufficient number of
sampled columns/rows is linear in r. Then, we compare the
proposed approach to the state-of-the-art randomized algorithm [5] and demonstrate that the proposed sampling strategy
can lead to notable improvement in performance. We then
provide an illustrative example to showcase the effectiveness
of our approach on real video frames for background subtraction and activity detection. Given the structure of the proposed
approach, it is shown that side information can be leveraged
to further simplify the decomposition task. In addition, a
numerical example is provided to examine the performance
of Algorithm 4. Finally, we investigate the performance of the
online algorithm and show that the proposed online method
can successfully track the underlying subspace.
In all simulations, the Augmented Lagrange multiplier
(ALM) algorithm [1], [4] is used to solve the optimization
problem (2). In addition, the `1 -magic routine [32] is used
to solve the `1 -norm minimization problems. It is important
to note that in all the provided simulations (except in Section
VII-D), the convex program (2) that operates on the entire data
can yield correct decomposition with respect to the considered
criteria. Thus, if the randomized methods cannot yield correct
decomposition, it is because they fall short of acquiring the
essential information through sampling.

Fig. 5. Phase transition plots for various rank and sparsity levels. White designates successful decomposition and black designates incorrect decomposition.

Fig. 6. Phase transition plots for various data matrix dimensions (r = 15, ρ =
0.05).

The top row of Fig. 5 confirms that the sufficient values for
m1 and m2 are roughly linear in r. For instance, when the rank
is increased from 5 to 25, the required value for m1 increases
from 30 to 140. In this experiment, the column and RS of
L are sampled from the random orthogonal model. Thus, the
CS and RS have small coherency whp [23]. Therefore, the
important factor governing the sample complexity is the rank
of L. Indeed, Fig. 6 shows the phase transition for different
sizes of the data matrix when the rank of L is fixed. One
can see that the required values for m1 and m2 are almost
independent of the size of the data confirming our analysis.

A. Phase transition plots
In this section, we investigate the required number of randomly sampled columns/rows. The LR matrix is generated as
a product L = Ur Qr , where Ur ∈ RN1 ×r and Qr ∈ Rr×N2 .
The elements of Ur and Qr are sampled independently from
a standard normal N (0, 1) distribution. The sparse matrix S
follows the Bernoulli model with ρ = 0.02. In this experiment, Algorithm 1 is used and the column/rows are sampled
uniformly at random.
Fig. 5 shows the phase transition plots for different numbers
of randomly sampled rows/columns. In this simulation, the
data is a 1000 × 1000 matrix. For each (m1 , m2 ), we generate
10 random realizations. A trial is considered successful if the
L̂kF
≤ 5×10−3 . It is clear
recovered LR matrix L̂ satisfies kL−
kLkF
that the required number of sampled columns/rows increases
as the rank or the sparsity parameter ρ are increased. When the
sparsity parameter is increased to 0.3, the proposed algorithm
can hardly yield correct decomposition. Actually, in this case
the matrix S is no longer a sparse matrix.

B. Efficient column/row sampling
In this experiment, Algorithm 3 is compared to the randomized decomposition algorithm in [5]. It is shown that
the proposed sampling strategy can effectively reduce the
required number of sampled columns/rows, and makes the
proposed method remarkably robust to structured data. In this
experiment, D is a 2000 × 4200 matrix. The LR component is
generated as L = [G1 G2 ...Gn ]. For 1 ≤ i ≤ n2 , Gi = Ui Qi ,
r
r
130r
where Ui ∈ R2000× n , Qi ∈ R n × n and the elements of Ui
and Qi are sampled independently from a normal distribution
N (0, 1). For n/2 + 1 ≤ i ≤ n, Gi = 13Ui Qi , where
r
r
10r
Ui ∈ R2000× n , Qi ∈ R n × n , and the elements of Ui and Qi
are sampled independently from an N (0, 1) distribution. We
set r equal to 60; thus, the rank of L is equal to 60 whp. The
sparse matrix S follows the Bernoulli model and each element
of S is non-zero with probability 0.02. In this simulation, we
do not use Algorithm 4 to form Dw . The matrix Dw is formed
from 300 uniformly sampled rows of D.

11

Fig. 7. Performance of the proposed approach and the randomized algorithm
in [5]. A value 1 indicates correct decomposition and a value 0 indicates
incorrect decomposition.

Fig. 8. Stationary background.

We evaluate the performance of the algorithms for different
values of n, i.e., different number of clusters. Fig. 7 shows
the performance of the proposed approach and the approach
in [5] for different values of m1 and m2 . For each value
of m1 = m2 , we compute the error in LR matrix recovery
kL−L̂kF
averaged over 10 independent runs, and conclude
kLkF
that the algorithm can yield correct decomposition if the
average error is less than 0.01. In Fig. 7, the values 0, 1
designate incorrect and correct decomposition, respectively. It
can be seen that the presented approach requires a significantly
smaller number of samples to yield the correct decomposition.
This is due to the fact that the randomized algorithm [5]
samples both the columns and rows uniformly at random
and independently. In sharp contrast, we use L̂w to find the
most informative columns to form Ds1 , and also leverage
the information embedded in the CS to find the informative
rows to form Ds2 . When n = 60, [5] cannot yield correct
decomposition even when m1 = m2 = 1800.
C. Vector decomposition for background subtraction
The LR plus sparse matrix decomposition can be effectively
used to detect a moving object in a stationary background
[1], [33]. The background is modeled as a LR matrix and the
moving object as a sparse matrix. Since videos are typically
high dimensional objects, standard algorithms can be quite
slow for such applications. Our algorithm is a good candidate
for such a problem as it reduces dimensionality significantly.
The decomposition problem can be further simplified by
leveraging prior information about the stationary background.
In particular, we know that the background does not change
or we can construct it with some pre-known dictionary.
For example, consider the video from [34], which was also
used in [1]. Few frames of the stationary background are
illustrated in Fig. 8. Thus, we can simply form the CS of
the LR matrix using these frames which can describe the

Fig. 9. Two frames of a video taken in a lobby. The first column displays
the original frames. The second and third columns display the LR and sparse
components recovered using the proposed approach.

stationary background in different states. Accordingly, we just
need to learn the representation matrix. As such, background
subtraction is simplified to a vector decomposition problem.
Fig. 9 shows that the proposed method successfully separates
the background and the moving objects. In this experiment,
500 randomly sampled rows are used (i.e., 500 randomly
sampled pixels) for representation matrix learning (14). While
the running time of our approach is just few milliseconds, it
takes almost half an hour if we use (2) to decompose the video
[1].
D. Alternating algorithm for column sampling
In this section, we investigate the performance of Algorithm
4 for column sampling. The rank of the selected columns is
shown to converge to the rank of L even when both the rows
and columns of L exhibit a highly structured distribution. To
generate the LR matrix L we first generate a matrix G as
in Section IV-A but setting r = 100. Then, we construct the
matrix Ug from the first r right singular vectors of G. We
then generate G in a similar way and set Vg equal to the first
r right singular vectors of G. Let the matrix L = Ug VgT .
For example, for n = 100, L ∈ R10250×10250 . Note that the
resulting LR matrix is nearly sparse since in this simulation
we consider a very challenging scenario in which both the
columns and rows of L are highly structured and coherent.
Thus, in this simulation we set the sparse matrix equal to zero
and use Algorithm 4 as follows. The matrix Dc is formed
using 300 columns sampled uniformly at random and the
following steps are performed iteratively:
1. Apply Algorithm 2 to DTc with C = 3 to sample approximately 3r columns of DTc and form Dw from the rows of D
corresponding to the selected rows of Dc .
2. Apply Algorithm 2 to Dw with C = 3 to sample approximately 3r columns of Dw and form Dc from the columns
of D corresponding to the selected columns of Dc . Fig.
10 shows the rank of Dc after each iteration. It is evident
that the algorithm converges to the rank of L in less than 3
iterations even for n = 100 clusters. For all values of n, i.e.,
n ∈ {2, 50, 60}, the data is a 10250 × 10250 matrix.
E. Online Implementation
In this section, the proposed online method is examined.
It is shown that the proposed scalable online algorithm tracks
the underlying subspace successfully. The matrix S follows the

12

Rank of the Sampled
Columns/Rows

100
90

Number of Clusters = 2
Number of Clusters = 50
Number of Clusters = 100

80
70
60
1

2

3

4

5

6

Iteration Number

Fig. 10. The rank of the matrix of sampled columns.

Bernoulli model with ρ = 0.01. Assume that the orthonormal
matrix U ∈ RN1 ×r spans a random r-dimensional subspace.
The matrix L is generated as follows.
For k from 1 to N2
1. Generate E ∈ RN1 ×r and q ∈ Rr×1 randomly.
2. L = [L Uq] .
3. If (mod(k, n) = 0)
U = approx-r(U + αE).
End If
End For
The elements of qi and E are sampled from standard normal
distributions. The output of the function approx-r is the matrix
of the first r left singular vectors of the input matrix and
mod(k, n) is the remainder of k/n. The parameters α and
n control the rate of change of the underlying subspace. The
subspace changes at a higher rate if α is increased or n is
decreased. In this simulation, n = 10, i.e., the CS is randomly
rotated every 10 new data columns. In this simulation, the
parameter r = 5 and N1 = 400. We compare the performance
of the proposed online approach to the online algorithm
in [35]. For our proposed method, we set C = 20 when
Algorithm 2 is applied to U, i.e., 20r rows of U are sampled.
The method presented in [35] is initialized with√the exact
CS and its tuning parameter is set equal to 1/ N1 . The
algorithm [35] updates the CS with every new data column.
The parameter nu of the proposed online method is set equal
to 4 (i.e., the CS is updated every 4 new data columns) and
the parameter ns is set equal to 5r. Define L̂ as the recovered
LR matrix. Fig. 11 shows the `2 -norm of the columns of
L − L̂ normalized by the average `2 -norm of the columns
of L for different values of α. One can see that the proposed
method can successfully track the CS while it is continuously
changing. The online method [35] performs well when the
subspace is not changing (α = 0), however, it fails to track
the subspace when it is changing.
A PPENDIX
Proof of Lemma 2
The selected columns of L can be written as Ls1 = LS1 .
Using the compact SVD of L, Ls1 can be rewritten as Ls1 =
UΣVT S1 . Therefore, to show that the CS of Ls1 is equal to
that of L, it suffices to show that the matrix VT S1 is a full
rank matrix. The matrix S1 selects m1 rows of V uniformly
at random. Therefore, using Theorem 2 in [8], if


3
m1 ≥ rγ 2 (V) max c2 log r, c3 log
,
(30)
δ

Fig. 11. Performance of the proposed online approach and the online
algorithm in [35].

then the matrix VT S1 satisfies the inequality
kI −

N2 T
1
V S1 ST1 Vk ≤
m1
2

(31)

with probability at least (1 − δ), where c2 , c3 are numerical
constants [8]. Accordingly, if σ1 and σr denote the largest and
smallest singular values of ST1 V, respectively, then
3m1
m1
≤ σ12 ≤ σr2 ≤
2N2
2N2

(32)

T
Therefore,
q the singular values of the matrix V S1 are greater
m1
T
than
2N2 . Accordingly, the matrix V S1 is a full rank
matrix.

Remark 4. A direct application of Theorem 2 in [8] would
in fact lead to the sufficient condition


3
m1 ≥ rγ 2 (R) max c2 log r, c3 log
,
(33)
δ
where R ∈ RN2 ×N2 denotes the matrix of right singular vectors of L. The bound in (30) is slightly tighter
since
it uses the incoherence parameter γ(V) ≤ γ(R) ,
√
N2 maxi,j |R(i, j)| in (33), where V consists of the first r
columns of R. This follows easily by replacing the incoherence
parameter in the step that bounds the `2 -norm of the row
vectors of the submatrix in the proof of ( [8], Theorem 2).
Proof of lemma 3
The sampled columns are written as Ds1 = DS1 = Ls1 +Ss1 .
First, we investigate the coherency of the new LR matrix Ls1 .
Define PST1 V as the projection matrix onto the CS of ST1 V
which is equal to the rows subspace of Ls1 . Therefore, the
projection of the standard basis onto the rows subspace of
Ls1 can be written as
maxkPST1 V ei k22 = maxkST1 V(VT S1 ST1 V)−1 VT S1 ei k22
i

≤
≤
≤

i
T
T
maxkS1 V(V S1 ST1 V)−1 VT ej k22
j
kST1 V(VT S1 ST1 V)−1 k2 kVT ej k22
γ 2 (V)r σ12
γ 2 (V)r 6N2
(6γ 2 (V))r
( 4) =
=
N2
σr
N2 m1
m1

(34)

where (ST1 V(VT S1 ST1 V)−1 VT S1 ) is the projection matrix
onto the CS of ST1 V. The first inequality follows from the
N2
1
fact that {S1 ei }m
i=1 is a subset of {ej }j=1 . The second

13

inequality follows from Cauchy-Schwarz inequality and the
third inequality follows from (4) and (32).
Using lemma 2.2 of [23], there exists numerical constant
µ r
c7 such that maxkUT ei k22 ≤ Np1 , with probability at least
i

N1 )
1 − 2N1−3 and µp = c7 max(r,log
.
r
In addition, we need to find a bound similar to the third
T
condition of (3) for the LR matrix Ls1 . Let Ls1 = Us1 Σs1 Vs1
be the SVD decomposition of Ls1 . Define
T
H = Us1 Vs1
=

r
X

i T
Uis1 (Vs1
)

(35)

i=1

Uis1

i
where
is the i column of Us1 and Vs1
is the ith column
of Vs1 . Given the random orthogonal model of the CS, H has
the same distribution as
r
X
0
i T
H =
i Uis1 (Vs1
)
(36)
th

i=1

where {i } is an independent Rademacher sequence. Using
Hoeffding’s inequality [36], conditioned on Us1 and Vs1 we
have
−t2
 0

2
P |H (i, j)| > t ≤ 2e 2hij ,
r
(37)
X
h2ij =
(Us1 (i, k))2 (Vs1 (j, k))2 .
k=1

Consider the following lemma adapted from Lemma 2.2 of
[23].
Lemma 13 (Adapted from lemma 2.2 of [23]). If the orthonormal
 matrix U follows the random orthogonal model,
1
then P |U(i, j)|2 ≥ 20 logNN
≤ 3N1−8 .
1
Therefore,
|Us1 (i, k)|2 ≤ 20

log N1
N1

(38)

with probability at least 1 − 3N1−8 . Thus, we can bound h2ij
as
log N1
kVs1 ei k22 .
(39)
h2ij ≤ 20
N1
Using (34), (39) can be rewritten as
h2ij ≤ 120

log N1 γ 2 (V)r
.
N1 m1

(40)

√
r
√
ω γ(V)
N1 m 1

Choose t =
for some constant ω. Thus, the unconditional form of (37) can be written as

√ 
−ζω 2
0
γ(V) r
P |H (i, j)| > ω √
≤ 2e log N1 +
N1 m1


(41)
log
N
γ 2 (V)r
1
2
P hij ≥ 120
N1 m1
0

for some numerical constant ζ. Setting ω = ζ log N1 where
0
ζ is a sufficiently large numerical constant gives

√ 
0
γ(V) r
P kH k∞ ≥ c9 log N1 √
≤ 3rN1−7
(42)
N1 m1

for some constant number c9 since (38) should be satisfied for
rN1 random variables.
Therefore, according to Lemma 1, if (9) is satisfied, the
convex algorithm (8) yields the exact decomposition with
probability at least 1 − c8 N1−3 .
Proof of lemma 4
Based on (1), the matrix of sampled rows can be written as
Ds2 = ST2 D = ST2 L + ST2 S = Ls2 + Ss2

(43)

T
Let Ls2 = Us2 Σs2 Vs2
be the compact SVD decomposition
c
c T
of Ls2 and Ls2 = Us2 Σcs2 (Vs2
) its complete SVD. It can
be shown [7] that (12) is equivalent to

min kẑi k1
ẑi

s.t.

⊥ T i
T
(U⊥
s2 ) ẑi = (Us2 ) Ss2 .

(44)

where Sis2 is the ith column of Ss2 and U⊥
s2 is the last (m2 −r)
columns of Ucs2 which are orthogonal to Us2 . In other words,
if q∗i is the optimal point of (12) and z∗i ∈ Rm2 is the optimal
point of (44), then z∗i = ST2 (di − Uq∗i ). Thus, it is enough to
show that the optimal point of (44) is equal to Sis2 .
The columns subspace of Us2 obeys the random orthogonal
model. Thus, U⊥
s2 can be modeled as a random subset of
Ucs2 . Based on the result in [8], if we assume that the sign of
the non-zero elements of Sis2 are uniformly random, then the
optimal point of (44) is Sis2 with probability at least (1 − δ)
provided that



m2 2
m2
i
2
c
, c5 log
(45)
m2 − r ≥ max c4 kSs2 k0 γ (Us2 ) log
δ
δ
for some fixed numerical constants c4 and c5 . The parameter
√
γ(Ucs2 ) = m2 max|Ucs2 (i, j)| and kSis2 k0 is the l0 -norm
i,j

of Sis2 . In this paper, we do not assume that the sign of the
non-zero elements of the sparse matrix S is random. However,
according to Theorem 2.3 of [1] (de-randomization technique)
if the locations of the nonzero entries of S follow the Bernoulli
model with parameter 2ρ, and the signs of S are uniformly
random and if (44) yields the exact solution who, then it is
also exact with at least the same probability for the model in
which the signs are fixed and the locations follow the Bernoulli
model with parameter ρ [1]. Therefore, it suffices to provide
the sufficient condition for the exact recovery of a random sign
sparse vector with Bernoulli parameter 2ρ.
First, we provide sufficient conditions to guarantee that
m2 − r ≥ c4 kSis2 k0 γ 2 (Ucs2 ) log

m2
δ

(46)

with high probability. Using Lemma 13 and the union
2
bound, max|Ucs2 (i, j)|2 ≤ 20 logmm
with probability at least
2
i,j

1 − 3m−6
2 .
Now, we find the sufficient number of randomly sampled
rows, m2 , to guarantee that (46) is satisfied with high probability. It is obvious that m2 < N1 . Define κ = logrN1 . Therefore,
it is sufficient to show that


m2
N1
≥
r
c
κ
log
+
1
(47)
6
δ
kSis2 k0

14

whp, where c6 = 20c4 . Suppose that
ρ≤

1
βr c6 κ log

N1
δ

+1

(48)



where β is a real
 number greater than one. Define α =
r c6 κ log Nδ1 + 1 . According to (48) and the Chernoff Bound
for Binomial random variables [37], we have
!


2
m
−a
2
P kSis2 k0 −
> a ≤ exp
.
(49)
a
2
βα
2( m
αβ + 3 )


If we set a = mα2 1 − β1 , then the inequality (47) is satisfied.
Therefore, (49) can be rewritten as



m2
m2
1
P kSis2 k0 −
>
1−
βα
α
β


(50)
−m22 (β − 1)2
3αβ
≤ 2 exp
.
α2 β 2
2m2 (β + 2)

Proof of Lemma 6
Since the rank of L is equal to r, the column spaces of
{Ui }ni=1 are independent r/n-dimensional subspaces and the
rank of the matrices {Qi }ni=1 is equal to r/n. The matrix L
can be expressed as L = [U1 ... Un ]O, where O is a block
matrix with blocks equal to the matrices {Qi }ni=1 . Suppose
Qok T ∈ Rnk ×r/n is an orthonormal matrix such that the row
space of Qok is equal to the row space of Qk . Thus, according
to Lemma 13,


log nk
≤ 3n−8
P |Qok (i, j)|2 > 20
(54)
k .
nk
Thus, maxi,j |Qok (i, j)|2 ≤ 20 lognknk with probability at least
−7
1 − 3r
n nk . Given the block diagonal structure of O,


γ 2 (V) = N2 max max |Qok (i, j)|2
≤ 20 log(max ni )

Therefore, if
2rβ(β − 2) log
m2 ≥
3(β − 1)2

i,j

k

i

1
δ

N2
min ni

(55)

i





N1
+1 ,
c6 κ log
δ

(51)

then the inequality (47) is satisfied with probability at least
(1 − δ). Accordingly, if (13) is satisfied, then (44) returns the
exact sparse vector with probability at least 1 − 3δ. The factor
0.5 in the numerator of the RHS of the first inequality of (13)
is due to the de-randomization technique [1] to provide the
guarantee for the fixed sign case.
Proof of Theorem 5
The proposed decomposition algorithm yields the exact decomposition if:
1. The sampled columns of the LR matrix span the the columns
subspace of L. Lemma 2 provides the sufficient conditions on
m1 to guarantee that the columns of Ls1 span U with high
probability.
2. The program (8) yields the correct LR and sparse components of Ds1 . Lemma 3 provides the sufficient conditions on
m1 and ρ to guarantee that Ds1 is decomposed correctly with
high probability.
3. The sampled rows of the LR matrix span the rows subspace
of L. Since it is assumed that the CS of L is sampled from
the random orthogonal model, according to Lemma 13


log N1
2
P max |U(i, j)| ≥ 20
≤ 3rN1−7 .
(52)
i,j
N1
Therefore, according to Lemma 2 if

 
0
0
3
m2 ≥ r log N1 max c2 log r, c3 log
,
δ

(53)

then the selected rows of the matrix L span the rows subspace
0
of L with probability at least (1 − δ − 3rN1−7 ) where c2 and
0
c3 are numerical constants.
4. The minimization (14) yields the correct RS. Lemma 4
provides the sufficient conditions to ensure that (12) yields
the correct representation vector. In order to guarantee the
performance of (14), we substitute δ with δ/N2 since (14) has
to return exact representation for the columns of D. Therefore,
P (Incorrect Decomposition) ≤ δ + c8 N1−3 + δ + 3rN1−7 + 3δ.

Pn

−7
with probability at least 1 −
i=1 ni . Now we compute
T
2
maxkV ei k2 . According to Lemma 2.2 in [23],
i


c7 max(r/n, log nk )
P max kQok ei k22 >
< 2n−3
(56)
k ,
i
nk
3r
n

where c7 is a numerical constant. Thus, based on (56) and the
block structure of O,
"


P max keTi Vk22 = max max kQok ei k22 >
i

i

k



#

c7 rϕ1  1 N2 
N2
n min nk

≤2

k

n
X

(57)
n−3
i

i=1

max(r/n,log max nk )

k
where ϕ1 =
.
r/n
Proof of Lemma 7
The proof of this lemma is similar to the proof of Lemma
6. We make use of the following lemma to establish a lower
bound on the coherency of the block diagonal matrix O. The
proof of Lemma 14 is provided at the end of the Appendix.

Lemma 14. Suppose U ∈ RN1 ×r is an orthonormal matrix
which spans a random r-dimensional
h matrix. If r ≥ 18 log
i N1
0.5 r
T
2
and N1 ≥ 96r log N1 , then P max kU ei k2 < N1
≤
2N1−5 .

i

Proof of Lemma 8
To prove Lemma 8, it suffices to ensure that the number of
sampled columns corresponding to each submatrix {Li =
Ui Qi }ni=1 is sufficiently large to span the column space of
{Li }ni=1 . Suppose Qok T ∈ Rnk ×r/n is an orthonormal matrix
and the RS of Qok is equal to the RS of Qk . According to
Lemma 2.2 in [23],


c7 max(r/n, log nk )
o
2
P max kQk ei k2 >
< 2n−3
(58)
k .
i
nk
Define Lsi as the columns sampled from submatrix Li .
According to Lemma 5 in [10] and (58), if the number of columns of Lsi is greater than or equal to ξi =

15

s
10c7 max(r/n, log ni ) log 2r
δ , then the CS of Li is equal to the
CS of Li with probability at least1 − δ/n − 2n−3
k . According
ξi N2
to Lemma 6 in [10], if m1 ≥ 2 + ξ3i log 2n
δ
ni , then
the number of columns sampled from Li is greater than ξi
with probability at least 1 − δ/n. Thus, if (18) is satisfied, the
sampled columns
Pn span the CS of L with probability at least
1 − 2δ − 2 i=1 ni−3 .

Lemma 15. Let A be an N × n matrix whose entries are
independent
standard normal variables. Then for every t ≥
p
2 log 2/δ,
√
√
√
√
N − n − t ≤ σmin (A) ≤ σmax (A) ≤ N + n + t (63)

Proof of Lemma 9
From Lemma 2.2 in [23],


c7 max(r, log N1 )
T
2
< 2N1−3 ,
P max kei Uk2 >
i
N1

Lemma 10 establishes a sufficient condition to guarantee the
performance of the CS learning step. In addition, since the
rows lie in one subspace and the imaginary sampler samples
the rows within each subspace uniformly at random, the
analysis of the representation learning step is similar to the
analysis of the corresponding step of Algorithm 1.

(59)

where U is an orthonormal basis for the CS of L. Thus,
according to Lemma 5 in [10], if (20) is true, the sampled
rows span the RS of L with probability at least 1 − 2N1−3 − δ.

with probability at least 1 − δ, where σmin (A) and σmax (A)
are the minimum and maximum singular values of A.

Proof of Lemma 14
First we review the result below from [40].

Proof of Lemma 10
Define Qsi as the columns of Qi corresponding to the sampled Lemma 16. Let Yr be a chi-squared random variable with r
columns from Li = Ui Qi . First, we establish an upper bound degrees of freedom. Then for each t > 0 and for each  ∈
on the RS coherency of the matrix of sampled columns. Define (0, 1)
h
√ i
2
PQsi T as the projection matrix onto the row space of Qsi . From
P Yr − r ≤ −t 2r ≤ e−t /2
Lemma 2.2 in [23],
(64)


−1
−2 r/4




P
Y
≥
r(1
−
)
≤
e
.
−3
r
r
c7 max(r/n, log Cr/n)
.
<2 C
P maxkPQsi T ei k22 >
i
Cr/n
n
(60)
The distribution of kUT ei k22 is equivalent to the distribution
of YYNr [23]. For each λ < 0, it follows from Lemma 16 that
Since the rank of L is equal to r, the matrices {Ui }ni=1 span n
1
"
√ #


independent r/n-dimensional linear subspaces. Thus, similar
√ YN
r
λ 2r
T
2
to the analysis used in the proof of Lemma 8,
= P Yr ≤ (r + λ 2r) 1
P kU ei k2 −
≤
N1
N1
N1 (65)



−3
h
i
√
c7 rϕ3
Cr
2
2
(61)
P maxkPST1 V ei k2 >
< 2n
,
≤ P Yr ≤ (r + λ 2r)(1 − )−1 + e− N1 /4 ,
i
m1
n
where PST1 V is the projection matrix onto the RS of Ls1 , m1
is the number of columns sampled to form Ds1 (which in this
Cr/n)
.
lemma is equal to Cr), and ϕ3 = max(r/n,log
r/n
Finally, similar to the analysis provided in the proof of
Lemma 3, if (21) is satisfied, then the convex algorithm
(8) yields exact decomposition with probability at least 1 −
00
−3
2N1−3 − 2n mn1
where µ is equal to (22).
Proof of Lemma 11
According to Lemma 6,

#
"
n
X
c
rϕ
1
N
7
1 
2 
P max keTi Vk22 >
≤2
n−3
(62)
i
i
N2
n min nk
i=1
k

Thus, based on (62) and the analysis provided in the proof
of Lemma 3, if (23) is satisfied, then (8) decomposes
Ds1
Pn
correctly with probability at least 1 − 2N1−3 − 2 i=1 n−3
i .
Proof of Theorem 12
The imaginary sampler samples m1 = Cr columns of D.
First, we ensure that the sampled columns corresponding to
each submatrix Li = Ui Qi span the column space of Li .
Define Qsi as the columns of Qi , corresponding to the columns
sampled from Li . Based on the following lemma from [38],
q
p
2
r
2n
[39], if Cr
≥
+
2
log
, then the rank of Qsi is
n
n
δ
r
equal to n with probability at least 1 − δ/n.

where the inequality follows from the
√ second inequality of
(64). Set  = 1+1√r and λ = −3/2 log N1 . If r > 1, then
(1 − )−1 ≤ 2. Thus, from Lemma 16, the first term of the
RHS of (65) can be expanded as



p
√
1
P Yr − r ≤ 2r √ − 3 log N1
≤ e−6 log N1 = N1−6 .
2
h
√ i
Accordingly, P kUT ei k22 − Nr1 ≤ λN12r
≤
N1−6 +
−N
√ 1

e 4( r+1)2 ≤ N1−6 + e
and r ≥ 18 log N1 , then
max kU
i

T

ei k22

−N1
16r

. Thus, if N1 ≥ 96r log N1

√
r
3/2 2r log N1
0.5 r
≥
−
≥
N1
N1
N1

(67)

with probability at least 1 − 2N1−5 .
R EFERENCES
[1] E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal component
analysis?” Journal of the ACM (JACM), vol. 58, no. 3, p. 11, 2011.
[2] V. Chandrasekaran, S. Sanghavi, P. A. Parrilo, and A. S. Willsky,
“Rank-sparsity incoherence for matrix decomposition,” SIAM Journal
on Optimization, vol. 21, no. 2, pp. 572–596, 2011.
[3] X. Yuan and J. Yang, “Sparse and low-rank matrix decomposition via
alternating direction methods,” preprint, vol. 12, 2009.
[4] Z. Lin, M. Chen, and Y. Ma, “The augmented Lagrange multiplier
method for exact recovery of corrupted low-rank matrices,” arXiv
preprint arXiv:1009.5055, 2010.

(66)

16

[5] L. W. Mackey, M. I. Jordan, and A. Talwalkar, “Divide-and-conquer
matrix factorization,” in Advances in Neural Information Processing
Systems (NIPS), 2011, pp. 1134–1142.
[6] X. Li and J. Haupt, “Identifying outliers in large matrices via randomized
adaptive compressive sampling,” IEEE Transactions on Signal Processing, vol. 63, no. 7, pp. 1792–1807, 2015.
[7] E. J. Candès and T. Tao, “Decoding by linear programming,” IEEE
Transactions on Information Theory, vol. 51, no. 12, pp. 4203–4215,
2005.
[8] E. Candès and J. Romberg, “Sparsity and incoherence in compressive
sampling,” Inverse problems, vol. 23, no. 3, p. 969, 2007.
[9] J. Wright, A. Ganesh, K. Min, and Y. Ma, “Compressive principal
component pursuit,” Information and Inference, vol. 2, no. 1, pp. 32–68,
2013.
[10] M. Rahmani and G. Atia, “Randomized robust subspace recovery for
high dimensional data matrices,” Preprint arXiv:1505.05901, 2015.
[11] D. Feldman, M. Monemizadeh, C. Sohler, and D. P. Woodruff, “Coresets
and sketches for high dimensional subspace approximation problems,”
in Proc. of the twenty-first annual ACM-SIAM symposium on Discrete
Algorithms, 2010, pp. 630–649.
[12] N. Halko, P.-G. Martinsson, and J. A. Tropp, “Finding structure with
randomness: Probabilistic algorithms for constructing approximate matrix decompositions,” SIAM review, vol. 53, no. 2, pp. 217–288, 2011.
[13] M. Rahmani and G. K. Atia, “Randomized subspace learning approach
for high dimensional low rank plus sparse matrix decomposition,” in
49th Asilomar Conference on Signals, Systems and Computers, 2015,
pp. 1796–1800.
[14] M. Rahmani and G. Atia, “A subspace learning approach for high
dimensional matrix decomposition with efficient column/row sampling,”
in Proceedings of The 33rd International Conference on Machine
Learning (ICML), 2016, pp. 1206–1214.
[15] T. Zhou and D. Tao, “Godec: Randomized low-rank & sparse matrix
decomposition in noisy case,” in International conference on machine
learning (ICML). Omnipress, 2011.
[16] Y. Mu, J. Dong, X. Yuan, and S. Yan, “Accelerated low-rank visual
recovery by random projection,” in IEEE Conference on Computer
Vision and Pattern Recognition (CVPR), 2011, pp. 2609–2616.
[17] R. Liu, Z. Lin, S. Wei, and Z. Su, “Solving principal component pursuit
in linear time via l 1 filtering,” arXiv preprint arXiv:1108.5359, 2011.
[18] S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, “A theory
of pseudoskeleton approximations,” Linear Algebra and its Applications,
vol. 261, no. 1, pp. 1–21, 1997.
[19] M. Mardani, G. Mateos, and G. B. Giannakis, “Subspace learning
and imputation for streaming big data matrices and tensors,” IEEE
Transactions on Signal Processing, vol. 63, no. 10, pp. 2663–2677, 2015.
[20] M. Mardani and G. B. Giannakis, “Online sketching for big data
subspace learning,” in Signal Processing Conference (EUSIPCO), 2015
23rd European. IEEE, 2015, pp. 2511–2515.
[21] J. Zhan, N. Vaswani, and C. Qiu, “Performance guarantees for reprocscorrelated low-rank matrix entries case,” in IEEE International Symposium on Information Theory, 2014, pp. 2187–2191.
[22] H. Guo, C. Qiu, and N. Vaswani, “An online algorithm for separating
sparse and low-dimensional signal sequences from their sum,” IEEE
Transactions on Signal Processing, vol. 62, no. 16, pp. 4284–4297, 2014.
[23] E. J. Candès and B. Recht, “Exact matrix completion via convex
optimization,” Foundations of Computational Mathematics, vol. 9, no. 6,
pp. 717–772, 2009.
[24] C. Boutsidis, M. W. Mahoney, and P. Drineas, “An improved approximation algorithm for the column subset selection problem,” in Proc.
of the twentieth Annual ACM-SIAM Symposium on Discrete Algorithms,
2009, pp. 968–977.
[25] M. W. Mahoney, “Randomized algorithms for matrices and data,”
Foundations and Trends® in Machine Learning, vol. 3, no. 2, pp. 123–
224, 2011.
[26] G. Liu and P. Li, “Recovery of coherent data via low-rank dictionary
pursuit,” in Advances in Neural Information Processing Systems (NIPS),
2014, pp. 1206–1214.
[27] A. Deshpande and L. Rademacher, “Efficient volume sampling for
row/column subset selection,” in 51st Annual IEEE Symposium on
Foundations of Computer Science (FOCS), 2010, pp. 329–338.
[28] J. Wright, A. Y. Yang, A. Ganesh, S. S. Sastry, and Y. Ma, “Robust face
recognition via sparse representation,” IEEE Transactions on Pattern
Analysis and Machine Intelligence, vol. 31, no. 2, pp. 210–227, 2009.
[29] Q. Ke and T. Kanade, “Robust l 1 norm factorization in the presence
of outliers and missing data by alternative convex programming,” in
IEEE Computer Society Conference on Computer Vision and Pattern
Recognition (CVPR), vol. 1, 2005, pp. 739–746.

[30] Z. Zhou, X. Li, J. Wright, E. Candès, and Y. Ma, “Stable principal
component pursuit,” in IEEE International Symposium on Information
Theory Proceedings (ISIT). IEEE, 2010, pp. 1518–1522.
[31] E. J. Candès, J. K. Romberg, and T. Tao, “Stable signal recovery from
incomplete and inaccurate measurements,” Communications on pure and
applied mathematics, vol. 59, no. 8, pp. 1207–1223, 2006.
[32] E. Candès and J. Romberg, “l1-magic: Recovery of sparse
signals via convex programming,” URL: www. acm. caltech.
edu/l1magic/downloads/l1magic. pdf, vol. 4, p. 14, 2005.
[33] T. Bouwmans, N. S. Aybat, and E.-h. Zahzah, “Handbook of robust lowrank and sparse matrix decomposition: Applications in image and video
processing,” 2016.
[34] L. Li, W. Huang, I. Y.-H. Gu, and Q. Tian, “Statistical modeling of complex backgrounds for foreground object detection,” IEEE Transactions
on Image Processing, vol. 13, no. 11, pp. 1459–1472, 2004.
[35] J. Feng, H. Xu, and S. Yan, “Online robust PCA via stochastic
optimization,” in Advances in Neural Information Processing Systems
(NIPS), 2013, pp. 404–412.
[36] W. Hoeffding, “Probability inequalities for sums of bounded random
variables,” Journal of the American statistical association, vol. 58, no.
301, pp. 13–30, 1963.
[37] C. McDiarmid, “Concentration,” in Probabilistic methods for algorithmic discrete mathematics. Springer, 1998, pp. 195–248.
[38] R. Vershynin, “Introduction to the non-asymptotic analysis of random
matrices,” arXiv preprint arXiv:1011.3027, 2010.
[39] K. R. Davidson and S. J. Szarek, “Local operator theory, random
matrices and Banach spaces,” Handbook of the geometry of Banach
spaces, vol. 1, pp. 317–366, 2001.
[40] B. Laurent and P. Massart, “Adaptive estimation of a quadratic functional
by model selection,” Annals of Statistics, pp. 1302–1338, 2000.

