1Sparse Bayesian Methods for Low-Rank Matrix
Estimation
S. Derin Babacan, Member, IEEE, Martin Luessi, Member, IEEE, Rafael Molina, Member, IEEE,
Aggelos K. Katsaggelos, Fellow, IEEE

Abstract
Recovery of low-rank matrices has recently seen significant activity in many areas of science and engineering,
motivated by recent theoretical results for exact reconstruction guarantees and interesting practical applications.
In this paper, we present novel recovery algorithms for estimating low-rank matrices in matrix completion androbust principal component analysis based on sparse Bayesian learning (SBL) principles. Starting from a matrix
factorization formulation and enforcing the low-rank constraint in the estimates as a sparsity constraint, we develop
an approach that is very effective in determining the correct rank while providing high recovery performance. We
provide connections with existing methods in other similar problems and empirical results and comparisons withcurrent state-of-the-art methods that illustrate the effectiveness of this approach.

I. I NTRODUCTION
Recently, there has been a significant interest in problems involving the estimation of low-rank matrices. Thisis motivated by recent theoretical advances [1]–[4], as well as interesting practical problems where the underlying
data resides in a low-dimensional linear subspace. Incorporating a low-rank constraint on the data to be processedleads to new and powerful modeling options for many applications in science and engineering.
A typical example is the matrix completion problem, where an unknown (approximately) low-rank matrix isestimated from its limited set of observed entries. Although this problem is not new [5], interesting and challenging
problems (e.g., the Netflix prize) along with recently developed theoretical recovery guarantees [1], [2] created a
rapidly growing interest in this area. Matrix completion finds application in many areas of engineering, including
system identification [6], sensor networks [7], machine learning [8], computer vision [9], [10], and medical imaging[11].
A second important problem is robust principal component analysis (RPCA), where the high dimensional data is
assumed to lie in a lower dimensional subspace with a small number of the data points corrupted with (arbitrarily)large errors. Widely used classical methods, such as principal component analysis (PCA), often fail to provide
meaningful results in these cases. Some earlier methods attempt to overcome these issues using robust statistics[12]–[17]. Recently, theoretical performance guarantees for RPCA have been developed in [3], where it is shownthat a data matrix can be decomposed into its low-rank and sparse components via convex optimization. Robust
PCA has many important applications, such as video surveillance (background/foreground separation in video), facerecognition [18], latent semantic indexing [19], image alignment [20], among many others.
Mathematically, problems involving the estimation of low-rank matrices can be formulated in a common framework as follows. Let X ∈ Rm×n be an unknown matrix with rank r ≪ min(m, n). Suppose that one is given an
observation matrix Y which is a function f (X) of X. In matrix completion, the observation is a subset Ω of its
entries, that is, {Yij = Xij : (i, j) ∈ Ω}. In other words, the observation Y is a projection PΩ of X on a subset Ωof its entries, such that the (i, j)th component of Y is equal to Xij if (i, j) ∈ Ω and zero otherwise. In RPCA, theS. Derin Babacan is with the Beckman Institute for Advanced Science and Technology, University of Illinois at Urbana-Champaign,
Urbana, IL, USA. e-mail: dbabacan@illinois.edu
Martin Luessi is with the Department of Radiology, Martinos Center for Biomedical Imaging, Massachusetts General Hospital, Harvard
Medical School, Boston, MA, USA e-mail: mluessi@nmr.mgh.harvard.edu
Rafael Molina is with the Departamento de Ciencias de la Computación e I.A. Universidad de Granada, Spain. e-mail: rms@decsai.ugr.es
Aggelos K. Katsaggelos is with the Department of Electrical Engineering and Computer Science, Northwestern University, IL, USA.
e-mail: aggk@eecs.northwestern.edu
This work has been supported in part by the Beckman Institute postdoctoral fellowship, ”Ministerio de Ciencia e Innovacion” under contract
TIN2010-15137, the Spanish research program Consolider Ingenio 2010: MIPRCV (CSD2007-00018), and a grant from the Department ofEnergy (DE-NA0000457).

2

observation can be expressed as Y = X + E, where E is a sparse error matrix where only a very small numberof coefficients are non-zero with (arbitrarily) large magnitudes.
In both cases, most matrices can be recovered by solving the affine rank minimization problem1 [1]–[4]minimize rank(X)subject to Y = f (X).

(1)

Although this optimization guarantees exact recovery of X under a set of conditions [1], [3], it is NP-hard andno known polynomial-time algorithms exist (analogous to the l0 -norm based recovery approaches in compressive
sensing). A popular approach is to utilize convex relaxation based on the nuclear norm, given byminimize kXk∗
subject to Y = f (X),(2)

where kXk∗ is equal to the sum of the singular values of X. Under some conditions the solutions of these two
problems coincide and recovery guarantees exist (see, for example, [1], [3], [21]). Subsequent works [2], [22], [23]improved on the theoretical recovery guarantees for the matrix completion problem.
If the observed entries are corrupted by dense (non-sparse) noise, the problem in (2) becomesminimize kXk∗
subject to k Y − f (X) k2F < ǫ,(3)

where k · kF denotes the Frobenius norm. Both nuclear norm based optimization problems in (2) and (3) can berecast as a semidefinite program, and can be solved with interior-point solvers [6], [24]. Although they provide
good empirical results, these methods can be inefficient when the matrix size is large.
A number of methods have been developed consequently for different problems involving low-rank estimation.
For matrix completion, singular value thresholding [25] and projection methods [26] are attractive in terms ofcomputation, while they nearly optimize (2). FPCA [27] introduced an efficient nuclear norm-based regularized
least-squares method, whereas OPTSPACE [22] developed a method based on optimization over the Grasmann
manifold with a theoretical performance guarantee for the noiseless case. Similarly to the approaches for compressive
sensing recovery, greedy approaches have been proposed for matrix completion [28]. Finally, Bayesian methods
have also been developed [29]–[34]: a nonparametric approach for symmetric positive definite matrices is proposedin [29], and a variational Bayes method is developed for collaborative filtering in [31]. The method in [32] is basedon beta-Bernoulli processes for modeling and Gibbs sampling for inference.
For robust PCA, the original work in [3] proposed iterative thresholding methods with low complexity, but
their convergence is generally very slow. Lin et al. [35] proposed accelerated proximal gradient (APG) methods
which are faster and generally more accurate. The augmented Lagrange Multiplier Method (ALM) [36] is, to thebest of our knowledge, the state-of-the-art method for robust PCA in terms of both speed and accuracy. However,algorithm parameters need to be tuned carefully to obtain the best performance. The Bayesian method proposed in[37] addresses this issue by simultaneously estimating the necessary parameters along with the unknowns, but theresulting algorithm uses sampling for inference and has high computational complexity.
In this paper, we present a novel Bayesian formulation for low-rank matrix recovery based on the sparse
Bayesian learning principles. We specifically consider the matrix completion and robust principal component analysis
problems, but the proposed framework can be translated to other problems involving low-rank structures. Based on
the low-rank factorization of the unknown matrix, we employ independent sparsity priors on the individual factors
with a common sparsity profile which favors low-rank solutions. Other elements in the problems are also modeledusing a hierarchical Bayesian framework for simultaneous and automated estimation.
The proposed Bayesian formulation offers several advantages over deterministic approaches. Firstly, prior knowledge on the rank of the matrix is not required; the proposed formulation implicitly estimates the rank of the unknown
matrix similarly to the automatic relevance determination principle in machine learning [38]. This property is notpresent in most of the existing deterministic approaches. Second, algorithmic parameters are treated as stochastic
quantities in the proposed approach, and are handled with the combination of prior distributions and fully-Bayesianinference procedures. In this regard, this type of formulation frees the user from extensive parameter-tuning and data1

A sparsity term is also incorporated in the objective function in the robust PCA case, which is omitted here for generality.

3

and application-dependent supervision. Finally, empirical results demonstrate that the proposed methods provide
very good reconstruction performance compared to existing methods while accurately estimating the unknown
effective rank.
This work is closely related to some probabilistic formulations in collaborative filtering [31], [39] and nonnegative
matrix factorization [40], regarding the modeling of the unknown low-rank components. The works [31], [39] arevariational Bayesian approaches to matrix completion, where the low-rank matrix is modeled via two factors andtwo sets of independent hyperparameters. On the other hand, [40] proposed to relate these factors by a single setof hyperparameters, and a maximum a posteriori approach is used for inference. One of the main contributionsof this work is to combine the modeling approaches in these works to obtain an intuitive modeling structure, andto provide a variational algorithm using this modeling which renders heuristic measures unnecessary and enables
fully-automated estimation without free parameters. Another contribution of this work is to extend the applicationof this low-rank modeling to the robust PCA problem by including additional hierarchical modeling for the sparse
errors with arbitrarily large coefficients. The relation with prior art is discussed in more detail later in this paper(Section IV-A).
The rest of this paper is organized as follows. We present the proposed Bayesian modeling in Section II. Section III
develops the estimation algorithms based on variational Bayesian inference. We present an analysis of the proposed
approach in Section IV and empirical results with synthetic and real data in Section V, and finally conclude in
Section VI.
II. BAYESIAN M ODELING
In order to simultaneously estimate all latent variables, we make use of a hierarchical Bayesian framework where
all observed and unknown quantities are treated as stochastic quantities and their joint probability distribution is
specified. For tractable mathematical modeling, this distribution is given in a factorized form using a generativemodel where each factor is a prior or a conditional distribution used to model a specific quantity. We provide thedescription of each distribution used in this work in the following sections.
A. Proposed Low-Rank Modeling
Our modeling is based on the low-rank parametrization of the unknown matrix X, given by
X = ABT ,(4)

where A is an m × r matrix, and B an n × r matrix, such that rank(X) = r ≤ min(m, n). Any matrix of rank r
can be decomposed in this form, as can be seen by considering the singular value decomposition



(5)
X = USVT = US1/2 S1/2 VT ,where U and V are respectively m × r and n × r matrices with orthogonal columns, and S is a r × r diagonal
matrix of the non-zero singular values. Algorithms based on this factorization are commonly used for nonnegative
matrix factorization [41] and matrix completion [42], which generally aim to find solutions tominimize k A k2F + k B k2F
subject to k Y − f (X) k2F < ǫ.

(6)The equivalence of this optimization problem to (3) can be shown (see [21]). We formulate the problem in (6)using the Bayesian methodology as follows. It is clear from X = ABT that X is the sum of outer-products of thecolumns of A and B, that is,X=kX

a·i b·i T ,

(7)

i=1

to denote the ith column and row of A, respectively. Notice that each outer-

where k ≥ r and we use a·i and ai·
product contributes at most one to the rank to X. Since a low-rank estimate of X is sought, our goal is to achieve

4column sparsity in A and B, such that most columns in A and in B are set equal to zero. To enforce this constraint,
we associate the columns of A and B with Gaussian priors of precisions (inverse variances) γi , that is,p(A|γ) =p(B|γ) =

k
Y

i=1
k
Y
i=1


N a·i |0, γi−1 Im ,(8)
N b·i |0, γi−1 In ,(9)where Im denotes the m × m identity matrix. Thus, the columns of A and B have the same sparsity profile
enforced by the common precisions γi . As shown later, many of the precisions γi will assume very large values
during inference, which effectively removes the corresponding outer-products from X, and hence reduces the rankof the estimate. This formulation is therefore the analog of sparse Bayesian learning formulation (or automatic
relevance determination) [38], [43] successfully utilized for compressive sensing reconstruction, where sparsityinducing Gaussian priors are employed on each of the coefficients of the unknown vector.
In addition to (8) and (9), we incorporate the conjugate Gamma hyperprior on the precisions γi
1p(γi ) = Gamma(a, ) ∝ γia−1 exp (−b γi ) .
b(10)

The parameters a and b are treated as deterministic whose values are set to small values (e.g., 10−5 ) to obtainbroad hyperpriors.
B. Observation and Noise Models
In this work, the prior structure in (8), (9) and (10) is used as a common low-rank matrix model for X in thematrix completion and robust PCA problems. The descriptions of the distributions used to model other latent andobserved variables are provided in the following sections.
1) Matrix Completion: In matrix completion, the observations are generated according to
Yij = Xij + Nij ,(i, j) ∈ Ω,

(11)or in a more compact form as
Y = PΩ (X + N) ,(12)

where N is the dense error matrix with coefficients Nij . The cardinality of the set Ω is pmn, with p the fraction ofobserved coefficients. Using this model, we follow the standard assumption and incorporate white Gaussian noisein the observations, such thatYp(Y|A, B, β) =N Yij |Xij , β −1 ,(13)(i,j)∈Ω

with β = 1/ǫ the noise precision. The noise precision β is assigned the noninformative Jeffrey’s priorp(β) = β −1 .

(14)p(Y, A, B, γ, β) = p(Y|A, B, β) p(A|γ) p(B|γ)p(γ) p(β) .

(15)

The joint distribution, therefore, is expressed as

52) Robust PCA: In this case, the generative model can be expressed as Y = X + E + N, where E is the sparse
error matrix with arbitrarily large coefficients, and N is the dense error matrix with relatively smaller coefficients.
Using white Gaussian noise modeling on N, we obtain the following conditional distribution for the observations
p(Y|A, B, E, β) = N Y|ABT + E, β −1 Imn


β
∝ exp
k Y − ABT − E k2F .
(16)
2

As in the matrix completion case, we assign the Jeffrey’s prior in (14) to β . The modeling of the sparse component
E is done by employing independent Gaussian priors on each of the coefficients Eij of the matrix E, that is,p(E|α) =nm YY

i=1 j=1


,
N Eij |0, α−1
ij

(17)

where α = {αij } and αij is the precision of the Gaussian on the (i, j)th coefficient. As with the noise precision,
we use Jeffrey’s priors on αij
p(αij ) = α−1
ij , ∀i, j .

(18)

Notice that when an individual precision goes to infinity, i.e., α−1
ij → 0, the corresponding coefficient Eij goes tozero. Hence, the sparsity in E is achieved when a large number of precision variables are set to high values. As
in the original formulation of sparse Bayesian learning, this is achieved in this work by simultaneously estimating
the coefficients Eij and the precision variables αij , as shown later.
Finally, the joint distribution is expressed asp(Y, A, B, E, γ, α, β) = p(Y|A, B, E, β) p(A|γ) p(B|γ)× p(E|α)p(γ) p(α) p(β).

(19)

III. A PPROXIMATE BAYESIAN I NFERENCE
As is widely known, exact full-Bayesian inference using joint distributions such as (15) and (19) is intractable,
since p(y) cannot be computed by marginalizing all latent variables. Therefore, approximation methods must be
utilized. Common approximations include maximum a posteriori (MAP) estimation, evidence-based analysis andvariational Bayes. Although all of these methods can only provide local minima, Bayesian inference (where atleast one variable is integrated out) is generally more effective in avoiding undesired local minima compared todeterministic methods such as MAP. This is mainly due to the fact that Bayesian methods approximate the fullposterior distributions instead of providing point-estimates of its modes. Although in theory sampling methods
can provide the optimal approximation to the posteriors, the computational complexity is significantly higher for
high-dimensional data than that of other methods, and convergence is generally hard to assess.
In this work, we present an inference procedure based on mean field variational Bayes [44], [45]. Our goal
is to compute posterior distribution approximations by minimizing the Kullback-Leibler (KL) divergence in analternating fashion for each latent variable. Let z be the vector of all latent variables such that z = (A, B, γ, β)for the matrix completion case, and z = (A, B, E, γ, α, β) for robust PCA. The posterior approximation q(zk ) of
each latent variable zk ∈ z is found using
log q(zk ) = h log p(Y, z)iz\zk + const,(20)

where z\zk denotes the set z with zk removed. The distribution p(Y, z) is the joint probability distribution givenin (15) for the matrix completion problem, and in (19) for robust PCA.
Q
Using mean field approximation, we employ the posterior factorization q(z) = q(zk ) such that the posterior
distribution of each unknown is estimated by holding the others fixed using their most recent distributions. Thus,
for each latent variable, the expectations of all parameters (excluding the current one) in the joint distribution aretaken with respect to their most recent distributions, and the result is normalized to find the approximate posterior
distribution. Since all distributions in the hierarchical model presented in the previous section are in the conjugate
exponential family, the form of each posterior approximation can be found without major difficulties. We present
the update rules resulting from this inference scheme for each problem in the following subsections.

6

A. Inference for Matrix Completion
1) Estimation of factors A and B: With some algebra, it follows from (20) that the approximation to theposterior distributions of A and B decompose as independent distributions of their rows. By combining the priorin (8) and the observation model in (13), the posterior density of the ith row ai· of A is found asq(ai· ) = N (ai· |hai· i, Σai ) ,(21)

hai· iT = hβi Σai hBi iT yi· T ,−1,Σai = hβi hBTi Bi i + Γ(22)with mean and covariance(23)

where the matrix Bi contains only the j th rows of B for which (i, j) ∈ Ω, such that,XX hbj· T bj· i =hbj· T ihbj· i + Σbj ,hBTi Bi i =j:(i,j)∈Ω(24)

j:(i,j)∈Ωwith Σbj the posterior covariance of the j th row of B. Additionally, the row vector yi· contains the observed entries
in the ith row of Y . Similarly, by combining the prior in (9) and the observation model in (13), the posterior
density of the j th row bj· of B is found as a normal distribution(25)
q(bj· ) = N bj· |hbj· i, Σbjwith parameters

hbj· iT = hβi Σbj hAj iT y·j ,Σbj = hβi hATj Aj i + Γ

(26)
−1

,

(27)

where Aj contains the ith rows of A for which (i, j) ∈ Ω, and the vector y·j is constructed from the observed
entries in the j th column of Y . It can be observed that the covariances Σbi of the estimate of B are incorporated
in the estimation of A (and vice versa).
2) Estimation of hyperparameters γ : By combining p(A|γ), p(B|γ) and p(γi ), the posterior density of γi
becomes a Gamma distribution
2b + ha·i T a·i i + hb·i T b·i i(a−1+ m+n2 )q(γi ) ∝ γi
exp −γi(28)2with meanhγi i =

2a + m + n.
2b + ha·i T a·i i + hb·i T b·i i(29)

The required expectations are given by
ha·i T a·i i = ha·i iT ha·i i +

X
j

T

T

hb·i b·i i = hb·i i hb·i i +Σaj ii ,

X
j

Σbj



ii

(30)
.

(31)

3) Estimation of noise precision β : The Bayesian methodology allows for the estimation of the noise precision
as well. From (20), the posterior approximation assumes a Gamma distribution with mean
pmn
hβi =
.
(32)
h k Y − PΩ (ABT ) k2F i
In summary, the algorithm proceeds by first estimating the rows of A and B using (22) and (26), respectively,
followed by the estimation of the precisions γi using (29), and (if desired) the noise precision β using (32). By
the properties of the variational Bayes methods, the algorithm is guaranteed to converge to a local minimum of thevariational bound [45].

7

B. Inference for Robust PCA
1) Estimation of factors A and B: The approximations to the posterior distributions of A and B take forms
similar to (21) and (25) with the same factorization over the rows of A and B, respectively. However, as opposedto the matrix completion case, the covariances Σai of the rows of A are equal since there are no missing values(the same applies to B). The posterior approximation of the ith row of A is given byq(ai· ) = N ai· |hai· i, ΣA ,(33)

with mean and covariancehai· iT = hβi ΣA hBiT (yi· − ei· )T ,
−1.
ΣA = hβi hBT Bi + Γ

(34)
(35)

Similarly, the posterior approximation of bj· is another multivariate normal distribution given byq(bj· ) = N bj· |hbj· i, ΣB(36)with parameters

hbj· iT = hβi ΣB hAiT (y·j − e·j ) ,−1.
ΣB = hβi hAT Ai + Γ(38)

hAT Ai = hAiT hAi + mΣA ,

(39)

T

(40)

(37)

The required expectations can be found asT

B

hB Bi = hBi hBi + nΣ .
T

Using these updates, the estimate of X is then found by X = hAihBi .
2) Estimation of E: Using (20), the posterior distribution approximation of E is found to be factorized on eachcoefficient Eij with distributions
q(Eij ) = N Eij |hEij i, ΣE
(41)
ij ,
with parameters

hEij i = hβi ΣE
ij
ΣE
ij =



Yij − hai· ihbj· iT1.
hβi + hαij i



,

(42)
(43)

Notice that (42) can be rewritten as
hEij i =hβi
Yij − hai· ihbj· iT ,hβi + hαij i(44)where the first term is at most 1, and hence this estimation of Eij corresponds to a shrinkage of the difference
between the observations and the low-rank estimates controlled by the noise precision β and the hyperparameters αij .
When a specific hyperparameter goes to infinity, that is, hαij i−1 → 0, the mean and variance of the corresponding
coefficient Eij become zero, resulting in sparse estimates of E.
3) Estimation of hyperparameters γ : Similarly to the above, the posterior density of γ is found as a Gamma
distribution with mean given in (29). The only difference is in the calculation of the expectations, which are givenby
ha·i T a·i i = ha·i iT ha·i i + m ΣA ii ,
(45)

T
T
B
hb·i b·i i = hb·i i hb·i i + n Σ ii .
(46)4) Estimation of hyperparameters α: The posterior density of hyperarameters αij is found as a Gamma distribution with mean1.
(47)hαij i =2
hEij i + ΣEij85) Estimation of noise precision β : Finally, the posterior approximation of the noise precision assumes a Gamma
distribution with mean
mn
hβi =
,
(48)
h k Y − ABT − E k2F i
where
h k Y − ABT − E k2F i =k Y − hAihBiT − hEi k2F+ Tr n hAiT hAiΣB + Tr m hBiT hBiΣA
n
m X
 X
ΣE
+ Tr mn ΣA ΣB +
ij .

(49)

i=1 j=1

In summary, the proposed algorithm estimates the low rank component X by estimating its factors using (34) and(37), followed by the estimation of the sparse matrix E using (42), and finally the estimation of all hyperparameters
using (29), (47) and (48), until convergence.
IV. D ISCUSSION
A. Related Prior Art
The methodology presented in this work is closely related to some methods developed for collaborative filtering,probabilistic principal component analysis (PCA) and (nonnegative) matrix factorization. In collaborative filtering
methods proposed in [31], [39], independent Gaussian priors are placed on the columns of A and B with separate
sets of variances, and a variational Bayesian analysis is employed for inference. Although these models are similar
to our approach, the columns of A and B are not coupled through the use of common precisions as in our work.
Employing common parameters is of crucial importance in removing redundant components from the estimatedmatrix and determining the effective rank. In theory, the modeling in (8) and (9) with common precisions is used
to represent the correlation between the columns of A and B, and it also removes possible scale problems due tothe use of separate sets of precisions. To cope with scalability issues, [31] uses fixed, heuristically selected values
for one set, and estimates the other hyperparameter set. Finally, in contrast to our work, [31] has reported that nosparsity in the precisions occurs during the application of their algorithm.
The idea of coupling the columns of A and B is also used in [40], which aims at solving the nonnegative matrix
factorization problem. This work, however, employs nonnegative priors on A and B, and resorts to a multiplicative
MAP based estimation procedure for the sake of maintaining nonnegativity. Note also that this method has notbeen developed to handle the missing values as in the matrix completion problem, or the large sparse errors as
in the robust PCA problem. Some statistical approaches [16], [17], [46] use heavy-tailed distributions for robust
estimation against outliers, but these do not include explicit modeling of sparse errors and hence cannot separate
these from dense errors.
The Bayesian PCA methods [47]–[49] also have some similarity with our approach (with a different priorstructure); these methods can be seen as marginalizing the matrix B out from the joint distribution and estimating
A only (or vice versa). Although a similar approach can be developed in our formulation, i.e., marginalize one
matrix factor to estimate the other, estimation of the common precisions γi becomes problematic since A and B
cannot be integrated out together from the joint distribution.
Finally, another Bayesian modeling and inference strategy is proposed in [37] for the robust PCA. The work
uses four distinct factors for the low-rank component, two of which are modeled using Gaussian priors similar tothis work, and the remaining two is used to model the sparse eigenvalues of the low-rank matrix. Sparseness is
explicitly imposed using a beta-Bernoulli hierarchical prior such that irrelevant components can be removed. This
is in contrast to our work and the approaches presented above, where the model does not lead to exact pruning, butrather to a “soft” pruning (by driving components to values numerically indistinguishable from machine precision).
The sparse component is modeled in a similar fashion in [37] with a combination of a beta-Bernoulli and normalGamma prior hierarchies. Due to the complex modeling, the posterior distributions can only be inferred using
sampling strategies.

9

B. Estimating the effective rank
The proposed algorithm enforces low-rank solutions by enforcing column sparsity in A and B. During inference,
most of the hyperparameters γi are driven to very large values, which will force the posterior means of the columns
to go to zero, effectively removing them from the model and reducing the rank. In our implementation, columnsof A and B were declared irrelevant at convergence if the corresponding γi−1 assumes a very small value (e.g.,
10−16 ).
C. Sparsity of the estimate of E
As discussed in Section III-B2, the update procedure (42) of the coefficients Eij is in fact a shrinkage procedure,
where the amount of shrinkage is controlled by the estimates of both the noise precision β and the hyperparameters
αij . This resembles closely the automatic relevance determination in the original work of relevance vector machines[38]. During the iterative procedure, many of the estimated precisions αij will approach very high values, whichmakes the corresponding posteriors in (41) very sharply peaked at zero. In the limit of αij → ∞, the posterior isinfinitely peaked at zero, leading to a zero estimate of hEij i in (42). In our implementation, we prune the coefficients
Eij with large corresponding αij values (e.g., 1016 ) via thresholding, leading to a sparse estimate of E. In addition,
we can find another update rule from (47) as(50)
hαij i hEij i2 + ΣE
ij = 1hαij inew hEij i2 + hαij iold ΣE
ij = 1oldhαij inew =1 − hαij ihEij i2

ΣE
ij

,

(51)

(52)

which is a fixed-point update for hαij i. We have also observed empirically that using these updates instead of (47)
leads to much faster convergence and enhanced sparsity, although no theoretical convergence guarantees exist. Notethat this update is also used in the original formulation of sparse Bayesian learning in [38].
D. Computational Complexity
While the proposed algorithms have demonstrated good empirical performance for a variety of matrix completion
and robust PCA problems, care must be taken when applied to large scale problems. In matrix completion, thecomputation of the inverse matrices in (23) and (27) can be quite expensive; their computation is O(k3 ), where kis the number of columns in each Aj matrix (or the number of columns in each Bi matrix). k is also equal to theestimated rank at each iteration. However, by construction, many rows of A (B) are removed to obtain Aj (Bi ),such that Aj (Bi ) might possibly have fewer rows than columns. Each Aj has on the average pm rows and kcolumns (recall p is the fraction of observed entries to the matrix size with p < 1). If pm < k, we can utilize the
Woodbury identity [50] to obtain a different form for Σbj , given by−1hAj iΓ−1 ,(53)
Σbj = Γ−1 − Γ−1 hAj iT hβi−1 Ik + hAj iΓ−1 hAj iT

which has the average-case complexity O(p3 m3 ). In practice, we compare the number of columns and rows in
Aj and Bi at each iteration to automatically choose the least complexity update. Overall, the complexity of thealgorithm is O(m·min(p3 n3 , k3 )+n·min(p3 m3 , k3 )). Empirically, however, we observed that convergence is rapid;
most of the precisions assume very large values in the very first iterations, and the norms of the corresponding
columns become numerically equal to zero, so that they can be removed from the model (similarly to [38]). Other
optimizations can also be implemented such as using the conjugate gradient method to solve for posterior meansin (22) and (26), and avoiding the computation of the off-diagonal terms of Σai and Σbj . These optimizations
will lead to decreased computational complexity at the expense of recovery performance. In the robust PCA case,an analysis similar to the above (using similar identities as (53)) gives an overall computational complexity ofO(min(n3 , k3 ) + min(m3 , k3 )) per iteration. However, as in the matrix completion case, the effective rank is
generally reduced rapidly in the first few iterations, therefore resulting in a very efficient inference scheme.

10

3

10

0.6
0.4

150

2

10

1

10

10

−1

10

15

20

25 30 35
True Rank

40

45

10

50

5

10

15

20

(a)

45

0
5

50

10

15

20

(b)

VSBL
OPTSPACESVTFPCA
ADMIRA

25 30 35
True Rank

40

45

50

40

45

50

(c)

4

200

3

10

0.6
0.4

150

2

10

1

10

100

50

0

0.2
0
5

40

Estimated Rank

Relative Error

0.8

25 30 35
True Rank

10

Running Times (s)

1

100

50

0

0.2
0
5

200

Estimated Rank

Relative Error

0.8

4

10

VSBL
OPTSPACESVTFPCA
ADMIRA

Running Times (s)

1

10

−1

10

15

20

25 30 35
True Rank

(d)

40

45

50

10

5

10

15

20

25 30 35
True Rank

40

(e)

45

50

0
5

10

15

20

25 30 35
True Rank(f)Fig. 1. Estimation results with matrices of size 500 × 500 with varying ranks when 20% of the entries are observed. Top: No observationnoise, Bottom: with observation noise. (a,d) Relative recovery error, (b,e) running times, and (c,f) estimated ranks. The legend is common
to all figures.

E. Initialization
Although randomly initializing the matrices A and B generally provided satisfactory results, faster convergenceand better reconstruction performance can be achieved by more carefully selecting the initial values. In our
1
1
implementations, we calculate the SVD of the matrix Y = USVT and set A = US 2 and BT = S 2 VT . With
this choice, the algorithm is initialized with a (near) full-rank matrix Y . On the other hand, one can initialize thealgorithm with a lower rank estimate by removing columns of A and B which correspond to small eigenvaluesof Y . Empirical results show negligible difference in performance if a reasonable initial rank (larger than the truerank) is chosen, whereas the computational complexity can be significantly reduced. Moreover, independently of theinitial rank, the algorithm successfully removes irrelevant components from the estimate and estimates the effective
rank accurately.
V. E MPIRICAL R ESULTS

In this section, we provide experimental results for the matrix completion and robust PCA problems with both synthetically generated and real data sets. To examine the empirical performance of the proposed method, we performedexperiments commonly used in the literature and compared the proposed methods to some existing algorithms. The
source code developed to obtain the results shown in this section can be found at https://netfiles.uiuc.edu/dbabac
A. Matrix Completion
Our first example illustrates the effectiveness of the proposed approach on determining the correct rank. We
generated test matrices X of size 500 × 500 of ranks r = 5, . . . , 50 by randomly sampling 500 × r matrices A
and B from a standard normal distribution N (0, 1) and setting X = ABT . The fraction of observed entries p is0.2, and they are sampled uniformly at random. For each experiment, the relative recovery error is measured as
k X̂ − X kF / k X kF , where X̂ is the estimate of X.

11

We present comparisons with the following algorithms: OPTSPACE [22], SVT [25], FPCA [27] and ADMIRA[28]. All of these are deterministic methods with different optimization strategies: OPTSPACE is based on optimization over the Grasmann manifold, SVT uses nuclear-norm minimization with singular-value thresholding,
FPCA is a fixed-point Bregman iterative method, and ADMIRA is an efficient greedy method which iteratively
adds components during reconstruction.
Our method, developed in Section III-A, is denoted by VSBL. We used the procedure proposed in [22] to estimate
the initial target rank required by ADMIRA and OPTSPACE. On the other hand, other methods automaticallyestimate the rank of the unknown matrix. Notice also that in the proposed method VSBL, all required parameters(including the noise variance) are estimated in an automated fashion. In VSBL, we use kX̂i − X̂i−1 kF /kX̂i−1 kF <10−5 as the convergence criteria, where X̂i and X̂i−1 are estimates of X in the ith and (i − 1)th iterations,
respectively.
We consider two test cases, one with noiseless observations, and one where observed entries are corrupted byzero-mean white Gaussian noise with standard deviation 0.05. Each simulation result is obtained by averaging 20
random instances. Figure 1 shows the relative reconstruction error, running times (on a 3GHz Core2 Duo CPU)and estimated ranks for each algorithm for both test cases. Among all algorithms, VSBL provides the highest
recovery performance for all ranks, and also estimates the correct rank in all cases where the rank r ≤ 35. As
expected, errors in both the recovery and the estimated rank increase as the original rank increases. OPTSPACE and
ADMIRA generally underestimate the rank, whereas FPCA and SVT consistently overestimate it. A similar behavior
is observed in the presence of observation noise: although the recovery performance of all algorithms decreases,
VSBL still exhibits a better ability to recover the original matrix and the correct rank than other methods.
We next consider another set of experimental conditions where 500 × 500 matrices of fixed rank of 10 are
generated, and the number of observed entries is varied according to different oversampling degrees of freedom.
Note that a matrix of size m×n of rank r depends upon r(m+n−r) degrees of freedom (df ), and the oversampling
degrees of freedom (osdf ) is defined as pmn/df with the number of measurements pmn [23]. Experimental resultsfor osdf = 2, 3, . . . , 7 are depicted in Figure 2 for the same noise conditions as above. The corresponding sampling
ratios are p ≈ 0.08, 0.12, 0.16, 0.20, 0.24, 0.28. It is evident that VSBL provides very accurate reconstructions andestimates the correct rank even with very low number of observations. In terms of computation time, ADMIRA
provided the best performance in most of the simulations, whereas execution times for VSBL were stable throughout
the testing conditions and were comparable to those of the other methods.
We next illustrate a real-world application of low-rank matrix completion methods on rating prediction fromexisting ratings. We use the Jester joke2 and MovieLens3 datasets, which are commonly used for testing recommendation systems. The Jester joke dataset contains user ratings on jokes where the ratings range from −10 to 10
with 200 quantization levels. The MovieLens dataset consists of user ratings on movies with integer ratings ranging
from 1 to 5. Most of the entries are not available in these datasets, and the goal is predicting the missing entries
by modeling the dataset as low-rank.
In the Jester joke data set, we generated a full rating matrix by removing all users containing missing entries,
and applied the algorithms to randomly generated subsets of this matrix with different number of users and fractionof observed ratings p. The number of jokes is fixed to 100. As the P
performance measure we use the normalized|Xij −X̂ij |mean absolute error (NMAE), which, for this dataset, is defined as (i,j)∈T20 |T |[27], with X̂ij the estimatedmissing components, T the set of missing entries, and |T | = p. It is known that as with most real data sets, Jester
data set is not low rank or even approximately low rank. To account for this in the proposed algorithm, we used
a fixed, high value for the noise variance (β −1 = 20) to encourage low-rank estimates (other values > 5 provided
very similar results). Numerical results (average of 10 realizations) are shown in Table I for two p values and threedifferent number of users. It can be observed that VSBL achieves a better prediction error than other algorithms
in all test cases.
In the MovieLens data set, we experimented with the 100k dataset with 100,000 ratings from 1000 users on
1700 movies, and the 1M dataset with 1,000,209 ratings from 6040 users on 3900 movies. In both datasets, we
randomly generated subsets of the rating matrices by sampling p = 0.1 and p = 0.5 of the available ratings foreach user. Note that the rating matrices are very sparse: The 100k dataset contains only about 5% of the entries,
2
3

Available at http://eigentaste.berkeley.edu/jester-data/
Available at http://www.grouplens.org/node/73/

12
4

50

−2

10

3

−4

10

−6

10

−8

10

Estimated Rank

Running Times (s)Relative Error

60

10

0

10

2

10

30
20

1

10

10

10
0

−10

10

2

3

4

5

6

10

7

2

3

4

(a)

6

0
2

7

(b)

Running Times (s)VSBL
OPTSPACESVTFPCA
ADMIRA
−1

10

10

100

3

80

10

2

10

1

10

4

5

6

7

10

5

6

7

5

6

7

60
40
20

0

−2

3

4

(c)

Estimated Rank

10

2

3

osdf

4

0

10

5
osdf

osdf

Relative Error

40

2

3

4

5

6

7

0
2

osdf

osdf

(d)

3

4
osdf

(e)

(f)

Fig. 2. Estimation results with matrices of size 500 × 500 of rank 10 with varying oversampling degrees of freedom. Top: No observationnoise, Bottom: with observation noise. (a,d) Relative recovery error, (b,e) running times, and (c,f) estimated ranks. The legend is common
to all figures.
TABLE INMAE VALUES ON THE J ESTER JOKE DATA SET# of usersVSBL
ADMIRA
OPT
SVT
FPCA

100
0.16250.16980.16850.18040.2026p = 0.1300
0.15940.17050.17000.16820.20461000
0.15840.16180.16100.16210.2052100
0.17200.17750.17440.19430.2096p = 0.5300
0.16690.17370.17150.18240.20601000
0.16260.17100.16940.17430.2051and the 1M dataset about 4%. Therefore, these datasets are extremely challenging for matrix completion methods
which do not take any other information into account (such as user and genre information). The NMAE results(average of 10 realizations) are shown in Table II. It can be observed that VSBL provides lower prediction errors
than other methods.
B. Robust PCA
1) Comparison with state-of-the-art: In our first experiment, we demonstrate the performance of the proposed
method using synthetic data in comparison with existing approaches. The low-rank component X is generated as in
Section V-A. The non-zero entries of the sparse matrix E ∈ Rm×n are located uniformly at random and are drawn
from a uniform distribution in the range [−10, 10]. The number of non-zero entries is set equal to 0.05mn. We
consider both a noise-free and a noisy case where white Gaussian noise with variance 10−3 is added to the originaldata. As before, the relative recovery error is measured as k X̂i − X kF / k X kF and k Êi − E kF / k E kF andthe convergence criterion is kX̂i − X̂i−1 kF /kX̂i−1 kF < 10−5 , where X̂i and Êi represent the estimates in the ith
iteration.

13

TABLE IINMAE VALUES ON THE M OVIELENS DATA SETS
Dataset
VSBL
OPT
FPCA

100k
p = 0.1
p=0.5
0.20450.13500.21100.17800.31750.2633

1M
p = 0.1
p=0.5
0.18400.14600.18750.17550.19900.1862

We present comparisons with the Bayesian method proposed in [37] (denoted by BRPCA) and the optimizationbased method in [36] (denoted by ALM). As mentioned before, BRPCA is based on factorizations of both the
low-rank and the sparse components, and low-rank and sparsity constraints are imposed using a combination ofbeta-Bernoulli priors on each component. The inference is performed using a Markov chain Monte Carlo (MCMC)
sampling scheme. On the other hand, ALM is based on soft-thresholding the singular values of the low-rankcomponent and elements of the sparse component. The inference is deterministic and is based on the augmented
Lagrange multiplier method. We use the exact inference method in [36] and manually tuned its parameters to reportits best performance in terms of recovery error. The proposed method, developed in Section III-B, is denoted as
VBRPCA.
Table III shows the relative reconstruction error, running times (on a 3GHz Core2 Duo CPU) and estimatedranks/sparsity levels for each algorithm for both noiseless and noisy cases. The average of 10 random instancesis reported in each experiment. It is clear that all methods provide very good reconstructions with both noiselessand noisy observations; both the low-rank and sparse components are recovered with high accuracy in all test
cases. While the running times of ALM and VBRPCA are very similar, the proposed method generally showedfaster convergence rate, especially in large matrix sizes. BRPCA, on the other hand, has a very high computational
complexity and therefore has longer running times in all test cases.
Although ALM is a very attractive method due to its recovery performance and fast convergence, it does
not provide means to estimate the dense noise level. Therefore, its convergence threshold should be adapted to
the noise variance to achieve the optimal performance, which requires user supervision. We empirically foundout that ALM is very sensitive to this parameter, and generally requires careful tuning (see [37] for a related
discussion). A comparison of ALM and VBRPCA is shown in Fig. 3, where matrices of size m × n withm = n = 500, 1000, 2000, 3000, 4000 are generated with the rank of the low-rank component equal to 0.05m
and the number of non-zeros in the sparse component equal to 0.05mn. Results are shown both with noiselessand noisy observations, where in the latter case noise variance is set equal to 10−3 . It is evident that while ALM
provides very low reconstruction errors in the noiseless case, its performance is significantly decreased and therank is consistently overestimated when dense noise is present. On the other hand, the performance of VBRPCA
is comparable to ALM in the noiseless case and better in the noisy case. In addition, VBRPCA estimates the rankcorrectly in both cases and requires lower computation times than ALM. It should be emphasized that ALM required
careful manual tuning of its convergence parameter in all experiments, while VBRPCA automatically estimates allalgorithmic parameters including the dense noise level. BRPCA has a similar mechanism for automatic noise
estimation through a Bayesian formulation, but its results are generally inferior to the proposed method and its
computational complexity is significantly higher.
Our second example illustrates a real-world application of robust PCA methods. We consider the foreground/backgroundseparation problem in video as in [37]. Each column of the data matrix Y is generated by concatenating pixels
of one video frame into a vector. In this application, the low-rank component corresponds to the background of
the scene, and the sparse component is used to model the moving objects in the foreground. It is clear that for acompletely static background, the ideal estimate of the rank of the background is 1, but in the case of dynamic
backgrounds (e.g., due to illumination changes), the rank can be higher.
All algorithms are applied to the video data4 consisting of 158 frames of size 192 × 144. Example results
obtained by the algorithms in one video frame are shown in Fig. 4. Due to the slow motion of the people, they
can be incorporated by mistake into the low-rank component (i.e., the background), which is the case with the4

The data can be found in http://homepages.inf.ed.ac.uk/rbf/CAVIARDATA1/.

14

TABLE IIIR ELATIVE R ECONSTRUCTION E RRORS , E STIMATED R ANKS AND C OMPUTATION T IMES FOR ROBUST PCA

Method
ALM
ALM
ALM
ALM
BRPCA
BRPCA
BRPCA
BRPCA
VBRPCAVBRPCAVBRPCAVBRPCA

σ
0
0
0
0
0
0
0
0
0
0
0
0

m=n
200
200
500
1000
200
200
500
1000
200
200
500
1000

rank(X)
5
10
25
50
5
10
25
50
5
10
25
50

k E k0
2000
2000
12500
50000
2000
2000
12500
50000
2000
2000
12500
50000

rank(X̂)
5
10
25
50
4
10
25
50
5
10
25
50

ALM
ALM
ALM
ALM
BRPCA
BRPCA
BRPCA
BRPCA
VBRPCAVBRPCAVBRPCAVBRPCA

10−3
10−3
10−3
10−3
10−3
10−3
10−3
10−3
10−3
10−3
10−3
10−3

200
200
500
1000
200
200
500
1000
200
200
500
1000

5
10
25
50
5
10
25
50
5
10
25
50

2000
2000
12500
50000
2000
2000
12500
50000
2000
2000
12500
50000

140
140
349
663
5
10
25
50
5
10
25
50

−3

kX̂−XkF
kXkF−8

3.5 × 10
1.5 × 10−8
0.6 × 10−8
4.4 × 10−8
3.5 × 10−5
5.2 × 10−4
1.1 × 10−4
9.5 × 10−5
0.5 × 10−6
1.8 × 10−6
2.1 × 10−6
3.5 × 10−6

kÊ−EkF
kEkF−7

1.3 × 10
0.7 × 10−7
0.5 × 10−7
6.0 × 10−7
1.3 × 10−5
1.2 × 10−4
7.2 × 10−4
1.9 × 10−4
4.0 × 10−7
2.0 × 10−7
4.8 × 10−7
1.7 × 10−7

time (s)
1.71
2.53
15.3
83.94
41.40
358.203116.33
27160.10
0.81
1.02
7.78
38.47

2.6 × 10−4
2.0 × 10−4
1.3 × 10−4
0.9 × 10−4
4.1 × 10−3
4.2 × 10−3
7.2 × 10−3
8.2 × 10−3
1.0 × 10−5
9.0 × 10−5
6.4 × 10−5
3.7 × 10−5

5.5 × 10−4
6.0 × 10−4
6.0 × 10−4
5.0 × 10−4
1.2 × 10−3
1.7 × 10−3
6.7 × 10−3
5.4 × 10−3
1.7 × 10−4
1.9 × 10−4
1.8 × 10−4
1.8 × 10−4

4.39
4.63
40.12
229.9845.39
370.213360.1027412.21
0.90
1.19
7.83
39.98

4

10

10

2500

−4

−5

10

−6

10

2000
3

10

ALM (no noise)
VBRPCA (no noise)
ALM (noisy)
VBRPCA (noisy)

2

10

−7

10

Estimated Rank

Running Times (s)Relative Error

10

1500
1000
500

−8

10

500

1000 1500 2000 2500 3000 3500 4000
m=n

(a)

500

1000 1500 2000 2500 3000 3500 4000
m=n

(b)

500

1000 1500 2000 2500 3000 3500 4000
m=n

(c)

Fig. 3. Comparison of ALM and VBRPCA with varying matrix sizes for robust PCA. Matrices are of size m × n with m = n, the rank of
the low-rank component is set equal to 0.05m, and the number of non-zeros in the sparse component is set equal to 0.05mn. (a) Relative
reconstruction errors, (b) running times, and (c) estimated ranks. The legend is common to all figures.

ALM algorithm. This is due to overfitting in the low-rank component, which was also observed in the synthetic
experiments with the ALM method in the presence of dense noise. The BRPCA algorithm provides a better result,but parts of the foreground are mistakingly classified as background. The proposed algorithm results in a muchcleaner separation, mainly due to the fact that a lower-rank estimate for the background is enforced compared tothe other methods (the estimated rank in this case is 1). This helps to avoid misclassification of foreground andbackground pixels. In this dataset, the running times of the algorithms were around 10 mins for ALM, 60 mins for
BRPCA, and 11 mins for the proposed method.
2) Comparison of inference methods: Although in this work we developed the algorithms based on variational
Bayesian inference, other inference methods can be employed as well based on the same Bayesian modeling shownin Sec. II. Here we compare VBRPCA with two other inference schemes, namely, maximum a posteriori (MAP)15(a) Original Frame(b) ALM(c) BRPCA(d) VBRPCA
Fig. 4. Video background/foreground separation. (a) Original video frame, the reconstructions by (b) ALM, (c) BRPCA, and (d) VBRPCA.
Left: background reconstruction, right: foreground reconstruction.

estimation and Gibbs sampling. This comparison will provide some insight both on the effectiveness and accuracy ofthe variational Bayesian approach. Moreover, it can be used to assess the accuracy of the variational approximation.
For both methods, using the observation model (16) and the priors given in (8), (9), (10), (14), (17), (18), we16form the conditional posterior distributions as
p(ai· |Y, z\ai· ) = N ai· |ãi· , Σ̃A ,
p(bi· |Y, z\bi· ) = N bi· |b̃i· , Σ̃B ,p(Eij |Y, z\Eij ) = N Eij |Ẽij , Σ̃E
ij ,

(54)
(55)
(56)

a+m+n
2
,
T
2
2b + ã·i ã·i + b̃T·i b̃·i
!
1 2
, 2 ,
2 Ẽijp(γi |Y, z\γi ) = Gammap(αij |Y, z\αij ) = Gammap(β|Y, z\β) = Gamma

2mn,
2 k Y − ÃB̃T − Ẽ k2F

!

,

!

,

(57)
(58)

(59)

where z is the set of all latent variables as before, Gamma(k, θ) is the Gamma distribution with shape parameterk and scale parameter θ (see (10)). The parameters of the distributions above are given byãi· = β̃ Σ̃A B̃T (yi· − ẽi· )T ,
−1
Σ̃A = β̃ B̃T B̃ + Γ̃,b̃Tj· = β̃ Σ̃B ÃT (y·j − ẽ·j ) ,−1Σ̃B = β̃ ÃT Ã + Γ̃,Γ̃ = diag (γ̃i ) ,


T
Y
−
ã
b̃
Ẽij = β̃ Σ̃Eij
i· j· ,
ij
Σ̃E
ij =

1

β̃ + α̃ij

.

(60)
(61)
(62)
(63)
(64)
(65)
(66)

The MAP estimates are found as the mode of these distributions, whereas in Gibbs sampling we sample fromthese distributions in an alternating fashion and collect the sampled values. In both cases, the estimated values aredenoted with a tilde (˜).
For empirical comparison, we create synthetic datasets similar to Section V-B1 where the low-rank component
X is 400 × 400 with coefficients drawn from a N (0, 1) distribution, the sparse matrix E has 8000 non-zero entries(sparsity level 5%) drawn from a uniform distribution [−10, 10] located uniformly at random. We consider bothnoiseless and noisy settings where in the latter case white Gaussian noise with variance 10−3 is added to theobservations. The rank of X is varied from 10 to 60 in steps of 10.
The relative reconstruction errors and rank estimates of X, along with average running times are shown in Fig. 5.
Reconstruction errors in the estimates of the sparse component E are similar to those of X and are not shown.
It can be seen that both variational Bayesian inference and Gibbs sampling provide more accurate estimates than
MAP. The MAP approach is very sensitive to the values of the hyperparameters a and b, and is prone to overand under-fitting depending on their selection. Similar results are obtained even when the correct noise level isprovided to MAP (data not shown), indicating that MAP is unable to avoid undesirable local minima. An importantresult is that the variational Bayesian inference provides results comparable to those of Gibbs sampling in termsof reconstruction error, and for rank ≤ 40, it correctly estimates the unknown rank and the sparsity level. In
addition, its running times are 2-4 orders of magnitude lower than those of Gibbs sampling, although a relatively
low number of iterations are used for the sampling method (10000 for burn-in and 2000 for collection, compared to25000 burn-in and 5000 collection in BRPCA [37]). Its running times are also comparable to the MAP approach,which has lower complexity per iteration but generally requires many more iterations for convergence.
In summary, the inference procedure developed in this work based on variational Bayesian analysis provides very
accurate results compared to sampling while being computationally considerably more efficient.

17
40

0

10

4

10

35
3

VB
Sampling
MAP

−4

10

−6

Running Times (s)Estimated Rank

10
Relative Error

10

30

−2

25
20
15
10

10

2

10

1

10

0

10

5
−8

10

10

20

30
40
True Rank

50

0
10

60

−1

20

30
40
True Rank

50

10

60

20

30
40
True Rank

50

60

20

30
40
True Rank

50

60

4

40

0

10

10

10

35
−1

3

−2

10

−3

10

Running Times (s)30
Estimated Rank

Relative Error

10

25
20
15
10

−4

10

2

10

1

10

10

5
−5

10

10

20

30
40
True Rank

50

60

0
10

0

20

30
40
True Rank

50

60

10
10

Fig. 5. Comparison of inference methods. Top row: noiseless case, bottom row: noisy case. From left to right: relative reconstruction errors,
estimated ranks, and running times.

VI. C ONCLUSIONSIn this paper, we have applied sparse Bayesian learning principles to the low-rank matrix estimation in matrix
completion and robust principal component analysis. We introduced a formulation where the low-rank constraint
is imposed on the estimate by using its sparse representation; starting from the factorized form of the unknown
matrix, we enforce a common sparsity profile on its underlying components using a probabilistic formulation. The
sparse error component in the robust PCA problem is also modeled and effectively inferred by sparse Bayesianlearning principles. We modeled the remaining unknown variables and observations within the hierarchical Bayesianframework and developed inference methods based on mean-field variational Bayes approximating the posteriors
of interest. This inference scheme is shown to be advantageous both in terms of computational requirements andestimation performance compared to other inference schemes. Empirical results suggest that the proposed algorithms
are very effective in pruning irrelevant dimensions and recover the correct number of effective components in thematrix estimate, and they provide competitive, and even higher, performance than current state-of-the-art approaches
in terms of reconstruction performance.
ACKNOWLEDGMENTS
The authors would like to thank Shinichi Nakajima, the associate editor and anonymous reviewers for theircomments and suggestions that helped improve the paper.
R EFERENCES[1] E. J. Candès and B. Recht, “Exact matrix completion via convex optimization,” Found. of Comput. Math., vol. 9, pp. 717–772, 2008.
[2] E. J. Candès and T. Tao, “The power of convex relaxation: Near-optimal matrix completion,” IEEE Trans. Inform. Theory, vol. 56,
no. 5, pp. 2053–2080, 2009.
[3] E. J. Candès, X. Li, Y. Ma, and J. Wright, “Robust principal component analysis?” CoRR, vol. abs/0912.3599, 2009.
[4] Z. Zhou, X. Li, J. Wright, E. Candès, and Y. Ma, “Stable principal component pursuit,” in Proceedings of IEEE International Symposiumon Information Theory (ISIT), 2010.
[5] C. R. Johnson, “Matrix completion problems: a survey,” in Proceed. of Symposia in Applied Mathematics, vol. 40, 1990, pp. 171–198.
[6] Z. Liu and L. Vandenberghe, “Interior-point method for nuclear norm approximation with application to system identification,” SIAM
Journal on Matrix Analysis and Applications, vol. 31, no. 3, pp. 1235–1256, 2009.
[7] S. Oh, A. Karbasi, and A. Montanari, “Sensor network localization from local connectivity: Performance analysis for the MDS-MAP
algorithm,” in IEEE Information Theory Workshop (ITW 2010), 2010.

18[8] N. Srebro, “Learning with matrix factorizations,” Ph.D. dissertation, Massachusetts Institute of Technology, Cambridge, MA, USA,2004.
[9] P. Chen and D. Suter, “Recovering the missing components in a large noisy low-rank matrix: Application to SFM,” IEEE Trans. PatternAnal. Mach. Intell., vol. 26, pp. 1051–1063, 2004.
[10] H. Ji, C. Liu, Z. Shen, and Y. Xu, “Robust video denoising using low rank matrix completion,” in IEEE Conference on Computer
Vision and Pattern Recognition (CVPR 2010), 2010.
[11] Z.-P. Liang, “Spatiotemporal imaging with partially separable functions,” in NFSI-ICFBI 2007, 2007, pp. 181 –182.
[12] P. J. Huber and E. M. Ronchetti, Robust Statistics. Wiley, 2009.
[13] F. D. la Torre and M. J. Black, “A framework for robust subspace learning,” International Journal of Computer Vision, vol. 54, pp.
117–142, 2003.
[14] Q. Ke and T. Kanade, “Robust l1 norm factorization in the presence of outliers and missing data by alternative convex programming,”
in IEEE Conf. Computer Vision and Pattern Recognition, 2005.
[15] M. A. Fischler and R. C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis andautomated cartography,” Communications of the ACM, vol. 24, pp. 381–385, 1981.
[16] J. Gao, “Robust l1 principal component analysis and its Bayesian variational inference,” Neural Computation, vol. 20, pp. 555–578,
2008.
[17] J. Luttinen, A. Ilin, and J. Karhunen, “Bayesian robust PCA for incomplete data,” in International Conference on Independent Component
Analysis and Signal Separation, 2009.
[18] A. Wagner, J. Wright, A. Ganesh, Z. Zhou, H. Mobahi, and Y. Ma, “Towards a practical face recognition system: Robust alignment
and illumination by sparse representation,” submitted to IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), July
2010.
[19] C. Papadimitriou, P. Rghavan, H. Tamaki, and S. Vempala, “Latent semantic indexing, a probabilistic analysis,” Journal of Computerand System Sciences, vol. 61, pp. 217–235, 2000.
[20] Y. Peng, A. Ganesh, J. Wright, W. Xu, and Y. Ma, “RASL: Robust alignment by sparse and low-rank decomposition for linearlycorrelated images,” submitted to IEEE Transactions on Pattern Analysis and Machine Intelligence (PAMI), July 2010.
[21] B. Recht, M. Fazel, and P. A. Parrilo, “Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization,”
SIAM Review, vol. 52, no. 3, pp. 471–501, 2010.
[22] R. H. Keshavan, S. Oh, and A. Montanari, “Matrix completion from a few entries,” submitted to IEEE Trans. Inf. Theory,
arXiv:0901.3150v2, 2009.
[23] E. J. Candès and Y. Plan, “Matrix completion with noise,” Proceedings of the IEEE, vol. 98, no. 6, pp. 925 – 936, 2009.
[24] V. Chandrasekharan, S. Sanghavi, P. Parillo, and A. Wilsky, “Rank-sparsity incoherence for matrix decomposition.” preprint, 2009.
[25] J.-F. Cai, E. J. Candès, and Z. Shen, “A singular value thresholding algorithm for matrix completion,” SIAM J. on Optimization, vol. 20,
no. 4, pp. 1956–1982, 2008.
[26] R. Meka, P. Jain, and I. S. Dhillon, “Guaranteed rank minimization via singular value projection,” arXiv:0909.5457, 2009.
[27] S. Ma, D. Goldfarb, and L. Chen, “Fixed point and Bregman iterative methods for matrix rank minimization,” arXiv:0905.1643v2,
2009.
[28] K. Lee and Y. Bresler, “ADMiRA: Atomic decomposition for minimum rank approximation,” arXiv:0905.0044, 2009.
[29] J. Paisley and L. Carin, “A nonparametric Bayesian model for kernel matrix completion,” in ICASSP 2010, Dallas, USA, April 2010.
[30] N. Ding, Y. A. Qi, R. Xiang, I. Molloy, and N. Li, “Nonparametric Bayesian matrix factorization by Power-EP,” in AISTATS, vol. 9,
2010, pp. 169–176.
[31] Y. J. Lim and Y. W. Teh, “Variational Bayesian approach to movie rating prediction,” in Proceedings of KDD Cup and Workshop,
2007.
[32] M. Zhou, C. Wang, M. Chen, J. Paisley, D. Dunson, and L. Carin, “Nonparametric Bayesian matrix completion,” in Proc. IEEE Sensor
Array and Multichannel Signal Processing Workshop (SAM2010), Oct 2010.
[33] F. Wood and T. L. Griffiths, “Particle filtering for nonparametric Bayesian matrix factorization,” in Advances in Neural Information
Processing Systems 19. Cambridge, MA: MIT Press, 2006, pp. 1513–1520.
[34] R. Salakhutdinov and A. Mnih, “Bayesian probabilistic matrix factorization using Markov chain Monte Carlo,” in Proceedings of the25th international conference on Machine learning, ser. ICML ’08. New York, NY, USA: ACM, 2008, pp. 880–887.
[35] Z. Lin, A. Ganesh, J. Wright, L. Wu, M. Chen, and Y. Ma, “Fast convex optimization algorithms for exact recovery of a corrupted
low-rank matrix,” SIAM J. Optimization (submitted), 2009.
[36] Z. Lin, M. Chen, L. Wu, and Y. Ma, “The augmented Lagrange multiplier method for exact recovery of corrupted low-rank matrices,”
University of Illinois at Urbana-Champaign, Tech. Rep., 2010.
[37] X. Ding, L. He, and L. Carin, “Bayesian robust principal component analysis,” submitted to IEEE Trans. Image Processing, 2010.
[38] M. Tipping, “Sparse Bayesian learning and the relevance vector machine,” Journal of Machine Learning Research, pp. 211–244, 2001.
[39] T. Raiko, A. Ilin, and J. Karhunen, “Principal component analysis for large scale problems with lots of missing values,” in Proceedingsof the 18th European Conference on Machine Learning (ECML 2007), J. Kok, J. Koronacki, R. Mantaras, S. Matwin, D. Mladenic,
and A. Skowron, Eds. Springer Berlin / Heidelberg, 2007, vol. 4701, pp. 691–698.
[40] V. Y. F. Tan and C. Févotte, “Automatic relevance determination in nonnegative matrix factorization,” in Proc. Workshop on Signal
Processing with Adaptive Sparse Structured Representations (SPARS’09), 2009.
[41] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix factorization,” in NIPS 13, 2001, pp. 556–562.
[42] J. Haldar and D. Hernando, “Rank-constrained solutions to linear matrix equations using power factorization,” IEEE Signal Process.
Lett., vol. 16, no. 7, pp. 584–587, 2009.
[43] D. Wipf, J. Palmer, and B. D. Rao, “Perspectives on sparse Bayesian learning,” NIPS 16, 2004.
[44] M. Beal, “Variational algorithms for approximate Bayesian inference,” Ph.D. dissertation, The Gatsby Computational Neuroscience
Unit, University College London, 2003.
[45] C. M. Bishop, Pattern Recognition and Machine Learning. Springer-Verlag, 2006.

19[46] B. Lakshminarayanan, G. Bouchard, and C. Archambeau, “Robust Bayesian matrix factorization,” in Proc. Intl. Conf. on Artificial
Intelligence and Statistics (AISTATS), 2011.
[47] M. E. Tipping, “Sparse kernel principal component analysis,” in NIPS, 2000, pp. 633–639.
[48] C. M. Bishop, “Bayesian PCA,” in NIPS, 1999, pp. 382–388.
[49] ——, “Variational principal components,” in Proceedings Ninth International Conference on Artificial Neural Networks, ICANN’99,
1999, pp. 509–514.
[50] K. Mardia, J. Kent, and J. Bibby, Multivariate analysis. Academic Press; New York, 1979.