Global Solution of Fully-Observed
Variational Bayesian Matrix Factorizationis Column-Wise Independent

Shinichi Nakajima
Nikon Corporation
Tokyo, 140-8601, Japan
nakajima.s@nikon.co.jp

Masashi Sugiyama
Tokyo Institute of Technology
Tokyo 152-8552, Japansugi@cs.titech.ac.jp

S. Derin Babacan
University of Illinois at Urbana-Champaign
Urbana, IL 61801, USA
dbabacan@illinois.edu

Abstract
Variational Bayesian matrix factorization (VBMF) efﬁciently approximates theposterior distribution of factorized matrices by assuming matrix-wise independence of the two factors. A recent study on fully-observed VBMF showed that,
under a stronger assumption that the two factorized matrices are column-wise independent, the global optimal solution can be analytically computed. However,it was not clear how restrictive the column-wise independence assumption is. In
this paper, we prove that the global solution under matrix-wise independence isactually column-wise independent, implying that the column-wise independence
assumption is harmless. A practical consequence of our theoretical ﬁnding is thatthe global solution under matrix-wise independence (which is a standard setup)
can be obtained analytically in a computationally very efﬁcient way without any
iterative algorithms. We experimentally illustrate advantages of using our analytic
solution in probabilistic principal component analysis.

1

Introduction

The goal of matrix factorization (MF) is to approximate an observed matrix by a low-rank one. In
this paper, we consider fully-observed MF where the observed matrix has no missing entry1 . This
formulation includes classical multivariate analysis techniques based on singular-value decomposition such as principal component analysis (PCA) [9] and canonical correlation analysis [10].
In the framework of probabilistic MF [20, 17, 19], posterior distributions of factorized matrices are
considered. Since exact inference is computationally intractable, the Laplace approximation [3],the Markov chain Monte Carlo sampling [3, 18], and the variational Bayesian (VB) approximation[4, 13, 16, 15] were used for approximate inference in practice. Among them, the VB approximation
seems to be a popular choice due to its high accuracy and computational efﬁciency.
In the original VBMF [4, 13], factored matrices are assumed to be matrix-wise independent, and alocal optimal solution is computed by an iterative algorithm. A simpliﬁed variant of VBMF (simpleVBMF) was also proposed [16], which assumes a stronger constraint that the factored matrices
1
This excludes the collaborative ﬁltering setup, which is aimed at imputing missing entries of an observed
matrix [12, 7].

1

are column-wise independent. A notable advantage of simpleVBMF is that the global optimal solution can be computed analytically in a computationally very efﬁcient way [15].
Intuitively, it is suspected that simpleVBMF only possesses weaker approximation ability due to its
stronger column-wise independence assumption. However, it was reported that no clear performancedegradation was observed in experiments [14]. Thus, simpleVBMF would be a practically usefulapproach. Nevertheless, the inﬂuence of the stronger column-wise independence assumption was
not elucidated beyond this empirical evaluation.
The main contribution of this paper is to theoretically show that the column-wise independence
assumption does not degrade the performance. More speciﬁcally, we prove that a global optimal
solution of the original VBMF is actually column-wise independent. Thus, a global optimal solution of the original VBMF can be obtained by the analytic-form solution of simpleVBMF—nocomputationally-expensive iterative algorithm is necessary. We show the usefulness of the analyticform solution through experiments on probabilistic PCA.

2

Formulation

In this section, we ﬁrst formulate the problem of probabilistic MF, and then introduce the VB approximation and its simpliﬁed variant.
2.1

Probabilistic Matrix Factorization

The probabilistic MF model is given as follows [19]:¶µ1p(Y |A, B) ∝ exp − 2 ∥Y − BA⊤ ∥2Fro ,
2σ
µ
¶
µ
¶
¢
¢
1 ¡
1 ¡
−1 ⊤
−1 ⊤
p(A) ∝ exp − tr ACA A,
p(B) ∝ exp − tr BCB B
,
2
2

(1)

(2)

where Y ∈ RL×M is an observed matrix, A ∈ RM ×H and B ∈ RL×H are parameter matrices to be
estimated, and σ 2 is the noise variance. Here, we denote by ⊤ the transpose of a matrix or vector, by∥ · ∥Fro the Frobenius norm, and by tr(·) the trace of a matrix. We assume that the prior covariancematrices CA and CB are diagonal and positive deﬁnite, i.e.,
CA = diag(c2a1 , . . . , c2aH ),CB = diag(c2b1 , . . . , c2bH )for cah , cbh > 0, h = 1, . . . , H.

Without loss of generality, we assume that the diagonal entries of the product CA CB are arrangedin the non-increasing order, i.e., cah cbh ≥ cah′ cbh′ for any pair h < h′ .
Throughout the paper, we denote a column vector of a matrix by a bold smaller letter, and a row
vector by a bold smaller letter with a tilde, namely,³´⊤e M )⊤ ∈ RM ×H , B = (b1 , . . . , bH ) = e
A = (a1 , . . . , aH ) = (ea1 , . . . , ab1 , . . . , ebL∈ RL×H .
2.2

Variational Bayesian Approximation

The Bayes posterior is written asp(A, B|Y ) =p(Y |A, B)p(A)p(B)
,
Z(Y )

(3)

where Z(Y ) = 〈p(Y |A, B)〉p(A)p(B) is the marginal likelihood. Here, 〈·〉p denotes the expectation over the distribution p. Since the Bayes posterior (3) is computationally intractable, the VB
approximation was proposed [4, 13, 16, 15].
Let r(A, B), or r for short, be a trial distribution. The following functional with respect to r is called
the free energy:À¿À¿r(A, B)r(A, B)= log− log Z(Y ). (4)F (r|Y ) = log
p(Y |A, B)p(A)p(B) r(A,B)
p(A, B|Y ) r(A,B)
2

In the last equation, the ﬁrst term is the Kullback-Leibler (KL) distance from the trial distribution
to the Bayes posterior, and the second term is a constant. Therefore, minimizing the free energy (4)amounts to ﬁnding the distribution closest to the Bayes posterior in the sense of the KL distance. In
the VB approximation, the free energy (4) is minimized over some restricted function space.
A standard constraint for the MF model is matrix-wise independence [4, 13], i.e.,rVB (A, B) = rAVB (A)rBVB (B).

(5)

This constraint breaks off the entanglement between the parameter matrices A and B, and leads to
a computationally-tractable iterative algorithm. Using the variational method, we can show that,
under the constraint (5), the VB posterior minimizing the free energy (4) is written asrVB (A, B) =MYeb m , ΣA )NH (eam ; a

m=1

L
Y

e
NH (e
bl ; bbl , ΣB ),l=1

where the parameters satisfy
³
´⊤
e
e
b= a
b ΣA ,
b1 , . . . , abMA= Y ⊤B
σ2
µ
¶⊤
e
b
b= e
b ΣB ,
B
b1 , . . . , b
bL
=YA
σ2

³
´−1
b⊤B
b + LΣB + σ 2 C −1
ΣA = σ 2 B
,
A

(6)

³
´−1
b⊤ A
b + M ΣA + σ 2 C −1
ΣB = σ 2 A
.
B

(7)

Here, Nd (·; µ, Σ) denotes the d-dimensional Gaussian distribution with mean µ and covariance mab ΣA , B,b and ΣB by Eqs.(6) and (7) until convergencetrix Σ. Iteratively updating the parameters A,
gives a local minimum of the free energy (4).
When the noise variance σ 2 is unknown, it can also be estimated based on the free energy minimization. The update rule for σ 2 is given by´³
b + LΣB )b + M ΣA )(Bb⊤BbAb⊤ ) + tr (Ab⊤ A∥Y ∥2Fro − tr(2Y ⊤ B
σ2 =
.
(8)
LM
Furthermore, in the empirical Bayesian scenario, the hyperparameters CA and CB are also estimatedfrom data. In this scenario, CA and CB are updated in each iteration by the following formulas:c2bh = ∥bbh ∥2 /L + (ΣB )hh .

c2ah = ∥bah ∥2 /M + (ΣA )hh ,
2.3

(9)

SimpleVB Approximation

A simpliﬁed variant, called the simpleVB approximation, assumes column-wise independence ofeach matrix [16, 15], i.e.,
rsimpleVB (A, B) =H
YrasimpleVB
(ah )
h

h=1

H
Y

rbsimpleVB(bh ).
h

(10)

h=1

This constraint restricts the covariances ΣA and ΣB to be diagonal, and thus necessary memory storage and computational cost are substantially reduced [16]. The simpleVB posterior can be written
asrsimpleVB (A, B) =HYb h , σa2h IM )NL (bh ; b
NM (ah ; abh , σb2h IL ),h=1

where the parameters satisfy

⊤
X
σa2h
b
 b
b h = 2 Y −
b⊤
a
bh′ a
bh ,
h′
σ
h′ ̸=h


X
σb2h
b
b
a
b⊤
bh ,
bh = 2 Y −
bh′ a
h′
σ
′

´−1
³
,
bh ∥2 + Lσb2h + σ 2 c−2
σa2h = σ 2 ∥b
ah

(11)

¡
¢−1
.
σb2h = σ 2 ∥bah ∥2 + M σa2h + σ 2 c−2
bh

(12)

h ̸=h

3

Here, Id denotes the d-dimensional identity matrix. Iterating Eqs.(11) and (12) until convergence,
we can obtain a local minimum of the free energy. Eqs.(8) and (9) are similarly applied if the noise
variance σ 2 is unknown and in the empirical Bayesian scenario, respectively.
A recent study has derived the analytic solution for simpleVB when the observed matrix has no
missing entry [15]. This work made simpleVB more attractive, because it did not only provide substantial reduction of computation costs, but also guaranteed the global optimality of the solution.
However, it was not clear how restrictive the column-wise independence assumption is, beyond its
experimental success [14]. In the next section, we theoretically show that the column-wise independence assumption is actually harmless.

3

Analytic Solution of VBMF under Matrix-wise Independence

Under the matrix-wise independence constraint (5), the free energy (4) can be written as
F = 〈log r(A) + log r(B) − log p(Y |A, B)p(A)p(B)〉r(A)r(B)
2

=

LM
M
|CA | L|CB | ∥Y ∥log σ 2 +
log
+ log
+
+ const.
2
2
|ΣA |
2
|ΣB |
2σ 2
´
³
´
1 n −1 ³ b⊤ b
−1 b ⊤ b+ tr CAA A + M ΣA + CB
B B + LΣB
2
³
³
´³
´´o
b⊤ Y ⊤ Bb+ Ab⊤ A
b + M ΣA Bb⊤B
b + LΣB+σ −2 −2A
.

(13)

b B,
b
Note that Eqs.(6) and (7) together form the stationarity condition of Eq.(13) with respect to A,
ΣA , and ΣB .
Below, we show that a global solution of ΣA and ΣB is diagonal. When the product CA CB is nondegenerate (i.e., cah cbh > cah′ cbh′ for any pair h < h′ ), the global solution is unique and diagonal.
On the other hand, when CA CB is degenerate, the global solutions are not unique because arbitraryrotation in the degenerate subspace is possible without changing the free energy. However, still oneof the equivalent solutions is always diagonal.
Theorem 1 Diagonal ΣA and ΣB minimize the free energy (13).
The basic idea of our proof is that, since minimizing the free energy (13) with respect to A, B, ΣA ,and ΣB is too complicated, we focus on a restricted space written in a particular form that includes
the optimal solution. From necessary conditions for optimality, we can deduce that the solutions ΣAand ΣB are diagonal.
Below, we describe the outline of the proof for non-degenerate CA CB . The complete proof forgeneral cases is omitted because of the page limit.
∗∗(Sketch of proof of Theorem 1) Assume that (A∗ , B ∗ , ΣA
, ΣB
) is a minimizer of the free energy(13), and consider the following set of parameters speciﬁed by an H × H orthogonal matrix Ω:

b = A∗ C −1/2 Ω ⊤ C 1/2 ,AAAΣA = CA ΩCA

b = B ∗ C 1/2 Ω ⊤ C −1/2 ,BAAΣB = CA

−1/2

1/2

−1/2

−1/2

∗
ΣA
CA

Ω ⊤ CA ,
1/2

−1/2

∗
CA Ω ⊤ CAΩCA ΣB
1/2

1/2

.

bA
b⊤ is invariant with respect to Ω, and (A,b B,b ΣA , ΣB ) = (A∗ , B ∗ , Σ ∗ , Σ ∗ ) holds if
Note that BABΩ = IH . Then, as a function of Ω, the free energy (13) can be simpliﬁed as¢ 1/2 ⊤ o
1 n −1 −1
1/2 ¡∗F (Ω) = tr CA
CB ΩCA B ∗⊤ B ∗ + LΣBCA Ω+ const.
2∗∗∗
This is necessarily minimized at Ω = IH , because we assumed that (A , B ∗ , ΣA
, ΣB
) is a mini∗⊤ ∗∗mizer. We can show that F (Ω) is minimized at Ω = IH only if B B + LΣB is diagonal. This∗implies that ΣA(see Eq.(6)) should be diagonal.
Similarly, we consider another set of parameters speciﬁed by an H × H orthogonal matrix Ω ′ :b = A∗ C 1/2 Ω ′⊤ C −1/2 ,A
B
BΣA = CBb=B
B

1/2−1/2 ∗ −1/2 ′⊤ 1/2
ΣB = CB Ω ′ CB ΣB
CB Ω CB .

∗

−1/2
1/2
CB Ω ′⊤ CB ,

−1/2

4

−1/2

∗
Ω ′ CB ΣA
CB Ω ′⊤ CB
1/2

1/2

,

Then, as a function of Ω ′ , the free energy (13) can be expressed as¢ 1/2 ′⊤ o
1 n −1 −1 ′ 1/2 ¡ ∗⊤ ∗∗F (Ω ′ ) = tr CA
CB Ω CB A A + M ΣA
CB Ω+ const.
2∗∗
Similarly, this is minimized at Ω ′ = IH only if A∗⊤ A∗ + M ΣA
is diagonal. Thus, ΣBshould be
diagonal (see Eq.(7)).
¤
The result that ΣA and ΣB become diagonal would be natural because we assumed the independent
Gaussian prior on A and B: the fact that any Y can be decomposed into orthogonal components may
imply that the observation Y cannot convey any preference for singular-component-wise correlation.
Note, however, that Theorem 1 does not necessarily hold when the observed matrix has missing
entries.
Theorem 1 implies that the stronger column-wise independence constraint (10) does not degrade
approximation accuracy, and the VB solution under matrix-wise independence (5) essentially agrees
with the simpleVB solution. Consequently, we can obtain a global analytic solution for VB, by
combining Theorem 1 above with Theorem 1 in [15]:Corollary 1 Let γh (≥ 0) be the h-th largest singular value of Y , and let ω ah and ω bh be theassociated right and left singular vectors:
LXY =γh ω bh ω ⊤
ah .
h=1

Let γ
bh be the second largest real solution of the following quartic equation with respect to t:fh (t) := t4 + ξ3 t3 + ξ2 t2 + ξ1 t + ξ0 = 0,
where the coefﬁcients are deﬁned byÃ!
p(L − M )2 γh(L2 + M 2 )ηh2
2σ 4
ξ3 =
, ξ2 = − ξ3 γh ++ 2 2
, ξ1 = ξ3 ξ0 ,LMLM
cah cbh
Ã
!2
µ
¶µ
¶
σ4
σ2 L
σ2 M
2
2
ξ0 = ηh − 2 2, ηh = 1 − 2
1− 2
γh2 .
cah cbh
γh
γh
Let

v
vÃ
u
!2
u
u
u (L + M )σ 2
44u (L + M )σ 2σσγ
eh = t+ 2 2 +t
+ 2 2
− LM σ 4 .
22cah cbh
22cah cbh

(14)

(15)

Then, the global VB solution under matrix-wise independence (5) can be expressed as½HXγbh if γh > γeh ,VBb VB ≡ 〈BA⊤ 〉rVB (A,B) = B
bA
b⊤ =
=
U
γ
bhVB ω bh ω ⊤
,
where
γ
b
ah
h
0
otherwise.
h=1

Theorem 1 holds also in the empirical Bayesian scenario, where the hyperparameters (CA , CB )are also estimated from observation. Accordingly, the empirical VB solution also agrees with theempirical simpleVB solution, whose analytic-form is given in Corollary 5 in [15]. Thus, we obtain
the global analytic solution for empirical VB:
Corollary 2 The global empirical VB solution under matrix-wise independence (5) is given by(HXγ̆hVB if γh > γ h and ∆h ≤ 0,
EVB
EVB
⊤
EVB
b
U
=
γ
bh ω bh ω ah , where γbh=0otherwise.
h=1

Here,
√
√
γ h = ( L + M )σ,
µ
¶
q
1
2
2
2
2
2
2
4
c̆h =
γh − (L + M )σ + (γh − (L + M )σ ) − 4LM σ ,
2LM
³ γ
´
³ γ
´
¢
1 ¡
h
h
∆h = M log
γ̆hVB + 1 + L log
γ̆hVB + 1 + 2 −2γh γ̆hVB + LM c̆2h ,
2
2
Mσ
Lσ
σ
and γ̆hVB is the VB solution for cah cbh = c̆h .
5

(16)
(17)
(18)

When we calculate the empirical VB solution, we ﬁrst check if γh > γ h holds. If it holds, we
compute γ̆hVB by using Eq.(17) and Corollary 1. Otherwise, γbhEVB = 0. Finally, we check if∆h ≤ 0 holds by using Eq.(18).
When the noise variance σ 2 is unknown, it is optimized by a naive 1-dimensional search to minimize
the free energy [15]. To evaluate the free energy (13), we need the covariances ΣA and ΣB , which
neither Corollary 1 nor Corollary 2 provides. The following corollary, which gives the complete
information on the VB posterior, is obtained by combining Theorem 1 above with Corollary 2 in[15]:
Corollary 3 The VB posteriors under matrix-wise independence (5) are given by
VB
rA
(A) =

H
Y

b h , σa2h IM ),NM (ah ; a

VB
rB
(B) =

h=1

H
Y

NL (bh ; bbh , σb2h IL ),h=1

where, for γ
bhVB being the solution given by Corollary 1,qqbh = ± γabhVB δbh · ω ah , b
bh = ± γ
bhVB δbh−1 · ω bh ,¡¢ p 2
− ηbh2 − σ 2 (M − L) + (bηh − σ 2 (M − L))2 + 4M σ 2 ηbh2
2
σah =
,
2M (b
γhVB δbh−1 + σ 2 c−2
ah )
¡
¢ p 2
− ηbh2 + σ 2 (M − L) + (b
ηh + σ 2 (M − L))2 + 4Lσ 2 ηbh2
2
σbh =
,
2L(b
γhVB δbh + σ 2 c−2
bh )
q
4 LM
bhVB )2 + 4σ
(M − L)(γh − γ
bhVB ) + (M − L)2 (γh − γ
c2a c2bhhδbh =,2σ 2 M c−2
ah
( 2
ηh
if γh > γeh ,ηbh2 =
σ4
otherwise.
c2 c2
ah bh

Note that the ratio cah /cbh is arbitrary in empirical VB, so we can ﬁx it to, e.g., cah /cbh = 1 without
loss of generality [15].

4Experimental Results

In this section, we ﬁrst introduce probabilistic PCA as a probabilistic MF model. Then, we show
experimental results on artiﬁcial and benchmark datasets, which illustrate practical advantages ofusing our analytic solution.
4.1

Probabilistic PCA

e ∈ RHIn probabilistic PCA [20], the observation y ∈ RL is assumed to be driven by a latent vector ain the following form:y = Bea + ε.
e and y, and ε ∈ RL is a Gaussian noiseHere, B ∈ RL×H speciﬁes the linear relationship between asubject to NL (0, σ 2 IL ). Suppose that we are given M observed samples {y 1 , . . . , y M } generatede M }, and each latent vector is subject to ae ∼ NH (0, IH ). Then,from the latent vectors {ea1 , . . . , a
the probabilistic PCA model is written as Eqs.(1) and (2) with CA = IH .
If we apply Bayesian inference, the intrinsic dimension H is automatically selected without predetermination [4, 14]. This useful property is called automatic dimensionality selection (ADS).
4.2

Experiment on Artiﬁcial Data

We compare the iterative algorithm and the analytic solution in the empirical VB scenario withunknown noise variance, i.e., the hyperparameters (CA , CB ) and the noise variance σ 2 are also
6

80

2.8

2.2

40

20

2
1.8
0

60

40

b
H

2.4

50

100
150
Iteration

200

0
0

250

Analytic
Iterative

80

60
Time(sec)

F /( LM )

2.6

100

Analytic
Iterative

Analytic
Iterative(a) Free energy

20

50

100
150
Iteration

200

0
0

250

(b) Computation time

50

100
150
Iteration

200

250

(c) Estimated rank

Figure 1: Experimental results for Artiﬁcial1 dataset, where the data dimension is L = 100, thenumber of samples is M = 300, and the true rank is H ∗ = 20.
80

30

3.2
Analytic
Iterative

25

Analytic
Iterative

Analytic
Iterative

60

2.8

20

b
H

Time(sec)

F /( LM )

3

40

15
10

20

2.6

5
0

50

100
150
Iteration

200

250

0
0

(a) Free energy

50

100
150
Iteration

200(b) Computation time

250

0
0

50

100
150
Iteration

200

250

(c) Estimated rank

Figure 2: Experimental results for Artiﬁcial2 dataset (L = 70, M = 300, and H ∗ = 40).

estimated from observation. We use the full-rank model (i.e., H = min(L, M )), and expect the
ADS effect to automatically ﬁnd the true rank H ∗ .
Figure 1 shows the free energy, the computation time, and the estimated rank over iterations for an
artiﬁcial (Artiﬁcial1) dataset with L = 100, M = 300, and H ∗ = 20. We randomly created true∗∗
matrices A∗ ∈ RM ×H and B ∗ ∈ RL×H so that each entry of A∗ and B ∗ follows N1 (0, 1). An
observed matrix Y was created by adding a noise subject to N1 (0, 1) to each entry of B ∗ A∗⊤ .
The iterative algorithm consists of the update rules (6)–(9). Initial values were set in the following
b and Bb are randomly created so that each entry follows N1 (0, 1). Other variables are set to
way: AΣA = ΣB = CA = CB = IH and σ 2 = 1. Note that we rescale Y so that ∥Y ∥2Fro /(LM ) = 1,
before starting iteration. We ran the iterative algorithm 10 times, starting from different initialpoints, and each trial is plotted by a solid line in Figure 1. The analytic solution consists of applying
Corollary 2 combined with a naive 1-dimensional search for noise variance σ 2 estimation [15]. The
analytic solution is plotted by the dashed line. We see that the analytic solution estimates the trueb = H ∗ = 20 immediately (∼ 0.1 sec on average over 10 trials), while the iterative algorithm
rank Hdoes not converge in 60 sec.
Figure 2 shows experimental results on another artiﬁcial dataset (Artiﬁcial2) where L = 70, M =300, and H ∗ = 40. In this case, all the 10 trials of the iterative algorithm are trapped at local
minima. We empirically observed a tendency that the iterative algorithm suffers from the local
minima problem when H ∗ is large (close to H).

4.3

Experiment on Benchmark Data

Figures 3 and 4 show experimental results on the Satellite and the Spectf datasets available from the
UCI repository [1], showing similar tendencies to Figures 1 and 2. We also conducted experiments
on various benchmark datasets, and found that the iterative algorithm typically converges slowly,
and sometimes suffers from the local minima problem, while our analytic-form gives the global
solution immediately.
7

5

400

Analytic
Iterative

Analytic
Iterative

50
40

4
3.5

300
200

3

100

2.5
0

0
0

50

100
150
Iteration

200

250

30

b
H

Time(sec)

F /( LM )500

Analytic
Iterative

4.5(a) Free energy

20
10
50

100
150
Iteration

200

0
0

250

(b) Computation time

50

100
150
Iteration

200

250

(c) Estimated rank

Figure 3: Experimental results for the Sat dataset (L = 36, M = 6435).
525

Analytic
Iterative

20

4

3.5

3
0

Analytic
Iterative

30
15

20

10

10

5

50

100
150
Iteration

200

250

Analytic
Iterative

40b
H

Time(sec)

F /( LM )

4.5

0
0

(a) Free energy

50

100
150
Iteration

200(b) Computation time

250

0
0

50

100
150
Iteration

200

250

(c) Estimated rank

Figure 4: Experimental results for the Spectf dataset (L = 44, M = 267).

5

Conclusion and Discussion

In this paper, we have analyzed the fully-observed variational Bayesian matrix factorization (VBMF)
under matrix-wise independence. We have shown that the VB solution under matrix-wise independence essentially agrees with the simpliﬁed VB (simpleVB) solution under column-wise independence. As a consequence, we can obtain the global VB solution under matrix-wise independence
analytically in a computationally very efﬁcient way.
Our analysis assumed uncorrelated priors. With correlated priors, the posterior is no longer uncorrelated and thus it is not straightforward to obtain a global solution analytically. Nevertheless, there
exists a situation where an analytic solution can be easily obtained: Suppose there exists an H × H′′
non-singular matrix T such that both of CA= T CA T ⊤ and CB= (T −1 )⊤ CB T −1 are diagonal.
We can show that the free energy (13) is invariant under the following transformation for any T :A → AT ⊤ ,ΣA → T ΣA T ⊤ ,CA → T CA T ⊤ ,B → BT −1 ,

ΣB → (T −1 )T ΣB T −1 ,CB → (T −1 )⊤ CB T −1 .

Accordingly, the following procedure gives the global solution analytically: the analytic solution′′
given the diagonal (CA, CB) is ﬁrst computed, and the above transformation is then applied.
We have demonstrated the usefulness of our analytic solution in probabilistic PCA. On the other
hand, robust PCA has gathered a great deal of attention recently [5], and its Bayesian variant has
been proposed [2]. We expect that our analysis can handle more structured sparsity, in addition to
the current low-rank inducing sparsity. Extension of the current work along this line will allow us togive more theoretical insights into robust PCA and provide computationally efﬁcient algorithms.
Finally, a more challenging direction is to handle priors correlated over rows of A and B. This
allows us to model correlations in the observation space, and capture, e.g., short-term correlationin time-series data and neighboring pixels correlation in image data. Analyzing such a situation, as
well as missing value imputation and tensor factorization [11, 6, 8, 21] is our important future work.
Acknowledgments
The authors thank anonymous reviewers for helpful comments. Masashi Sugiyama was supportedby the FIRST program. S. Derin Babacan was supported by a Beckman Postdoctoral Fellowship.
8

References[1] A. Asuncion and D.J. Newman. UCI machine learning repository, 2007.
[2] S. D. Babacan, M. Luessi, R. Molina, and A. Katsaggelos. Sparse Bayesian methods for low-rank matrix
estimation. arXiv:1102.5288v1 [stat.ML], 2011.
[3] C. M. Bishop. Bayesian principal components. In Advances in NIPS, volume 11, pages 382–388, 1999.
[4] C. M. Bishop. Variational principal components. In Proc. of ICANN, volume 1, pages 514–509, 1999.
[5] E.-J. Candes, X. Li, Y. Ma, and J. Wright. Robust principal component analysis? CoRR, abs/0912.3599,
2009.
[6] J. D. Carroll and J. J. Chang. Analysis of individual differences in multidimensional scaling via an n-way
generalization of ’eckart-young’ decomposition. Psychometrika, 35:283–319, 1970.
[7] S. Funk. Try this at home. http://sifter.org/˜simon/journal/20061211.html, 2006.
[8] R. A. Harshman. Foundations of the parafac procedure: Models and conditions for an ”explanatory”
multimodal factor analysis. UCLA Working Papers in Phonetics, 16:1–84, 1970.
[9] H. Hotelling. Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 24:417–441, 1933.
[10] H. Hotelling. Relations between two sets of variates. Biometrika, 28(3–4):321–377, 1936.
[11] T. G. Kolda and B. W. Bader. Tensor decompositions and applications. SIAM Review, 51(3):455–500,
2009.
[12] J. A. Konstan, B. N. Miller, D. Maltz, J. L. Herlocker, L. R. Gordon, and J. Riedl. Grouplens: Applying
collaborative ﬁltering to Usenet news. Communications of the ACM, 40(3):77–87, 1997.
[13] Y. J. Lim and T. W. Teh. Variational Bayesian approach to movie rating prediction. In Proceedings of
KDD Cup and Workshop, 2007.
[14] S. Nakajima, M. Sugiyama, and S. D. Babacan. On Bayesian PCA: Automatic dimensionality selection and analytic solution. In Proceedings of 28th International Conference on Machine Learning(ICML2011), Bellevue, WA, USA, Jun. 28–Jul.2 2011.
[15] S. Nakajima, M. Sugiyama, and R. Tomioka. Global analytic solution for variational Bayesian matrix factorization. In J. Lafferty, C. K. I. Williams, R. Zemel, J. Shawe-Taylor, and A. Culotta, editors, Advancesin Neural Information Processing Systems 23, pages 1759–1767, 2010.
[16] T. Raiko, A. Ilin, and J. Karhunen. Principal component analysis for large scale problems with lots ofmissing values. In J. Kok, J. Koronacki, R. Lopez de Mantras, S. Matwin, D. Mladenic, and A. Skowron,
editors, Proceedings of the 18th European Conference on Machine Learning, volume 4701 of Lecture
Notes in Computer Science, pages 691–698, Berlin, 2007. Springer-Verlag.
[17] S. Roweis and Z. Ghahramani. A unifying review of linear Gaussian models. Neural Computation,
11:305–345, 1999.
[18] R. Salakhutdinov and A. Mnih. Bayesian probabilistic matrix factorization using Markov chain Monte
Carlo. In International Conference on Machine Learning, 2008.
[19] R. Salakhutdinov and A. Mnih. Probabilistic matrix factorization. In J. C. Platt, D. Koller, Y. Singer,
and S. Roweis, editors, Advances in Neural Information Processing Systems 20, pages 1257–1264, Cambridge, MA, 2008. MIT Press.
[20] M. E. Tipping and C. M. Bishop. Probabilistic principal component analysis. Journal of the Royal
Statistical Society, 61:611–622, 1999.
[21] L. R. Tucker. Some mathematical notes on three-mode factor analysis. Psychometrika, 31:279–311, 1996.

9