On Bayesian PCA: Automatic Dimensionality Selection and Analytic Solution

NAKAJIMA . S @ NIKON . CO . JP

Shinichi Nakajima
Nikon Corporation, Tokyo, 140-8601 Japan

SUGI @ CS . TITECH . AC . JP

Masashi Sugiyama
Tokyo Institute of Technology, Tokyo 152-8552, Japan
Derin Babacan
Beckman Institute, University of Illinois, Urbana-Champaign, USA

Abstract

DBABACAN @ ILLINOIS . EDU

1. Introduction

In probabilistic PCA, the fully Bayesian estimation is computationally intractable. To copewith this problem, two types of approximation
schemes were introduced: the partially Bayesian
PCA (PB-PCA) where only the latent variables
are integrated out, and the variational Bayesian
PCA (VB-PCA) where the loading vectors arealso integrated out. The VB-PCA was proposed
as an improved variant of PB-PCA for enabling
automatic dimensionality selection (ADS). In
this paper, we investigate whether VB-PCA is really the best choice from the viewpoints of computational efﬁciency and ADS. We ﬁrst show thatADS is not the unique feature of VB-PCA—PBPCA is also actually equipped with ADS. We further show that PB-PCA is more advantageous incomputational efﬁciency than VB-PCA because
the global solution of PB-PCA can be computed
analytically. However, we also show the negative fact that PB-PCA results in a trivial solution in the empirical Bayesian framework. We
next consider a simpliﬁed variant of VB-PCA,
where the latent variables and loading vectors are
assumed to be mutually independent (while theordinary VB-PCA only requires matrix-wise independence). We show that this simpliﬁed VBPCA is the most advantageous in practice because its empirical Bayes solution experimentally works as well as the original VB-PCA, andits global optimal solution can be computed efﬁciently in a closed form.

Principal component analysis (PCA) is a well-establishedtechnique for unsupervised dimensionality reduction(Hotelling, 1933). A decade ago, PCA was given its
probabilistic interpretation as a latent variable model
called probabilistic PCA (PPCA) (Tipping & Bishop,
1999; Roweis & Ghahramani, 1999). In the ﬁrst formulation of PPCA, only the latent variables are integrated out
and the loading vectors are estimated by maximizing themarginal likelihood (Tipping & Bishop, 1999). Since theloading vectors are point-estimated, we refer to this approach as partially Bayesian PCA (PB-PCA) in the following.

Appearing in Proceedings of the 28 th International Conference
on Machine Learning, Bellevue, WA, USA, 2011. Copyright 2011by the author(s)/owner(s).

The purpose of this paper is to revisit PPCA, and investigate whether VB-PCA is really the best choice from theviewpoints of computational efﬁciency and ADS, within a

Subsequently, a variant of PPCA was proposed (Bishop,
1999b), which has the following two features:(A) The fully Bayesian treatment: both the latent variables
and the loading vectors are integrated out. We refer tothis as fully Bayesian PCA (FB-PCA).
(B) The empirical Bayesian procedure: the prior variancesof the loading vectors can also be estimated from observation.
A notable advantage of FB-PCA is that it offers automatic dimensionality selection (ADS). However, since exact FB-PCA is computationally intractable, the Laplace approximation (Bishop, 1999b; Hoyle, 2008), Markov chain
Monte Carlo (Bishop, 1999b), and the variational approximation (Bishop, 1999a) were used for approximate inference in practice. Among them, the variational approximation (which we refer to as VB-PCA) seems to be the mostpopular choice.

On Bayesian PCA: Automatic Dimensionality Selection and Analytic Solution

uniﬁed framework of free energy minimization under different constraints on the posterior distribution.
First, in defense of PB-PCA, we theoretically show the following two positive facts:• Neither (A) nor (B) is essential for PPCA to induce
ADS. PB-PCA is actually equipped with ADS. Thus,
the original unique advantage of VB-PCA is lost.
• The global solution of PB-PCA can be computed very
efﬁciently in a closed form (which is an extension of
the result given in Tipping & Bishop (1999)). On theother hand, VB-PCA requires a number of iterations
to ﬁnd a local optimal solution (Bishop, 1999a).

analytic-form solution for PB-PCA, and discuss its behavior in comparison with VB-PCA. In Section 4, we introduce a simpliﬁed variant of VB-PCA, and compare its behavior with the original VB-PCA. We further discuss various issues of PPCA in Section 5, and conclude in Section 6.

2. Formulation
In this section, we formulate PPCA, and review approximation methods to the fully Bayesian inference. Then, we
describe the empirical Bayesian procedure.
2.1. Probabilistic PCA
PPCA assumes that the observation y ∈ RL is driven by ae ∈ RH as follows:
latent vector a

These facts encourage the use of PB-PCA. However, wenext show that the following negative fact:• PB-PCA has a critical drawback when hyperparameters are learned from data in the empirical Bayesianframework. More precisely, the empirical Bayes
method in PB-PCA results in a trivial solution (i.e.,zero), and thus cannot be used in practice. On theother hand, the empirical VB-PCA works well.
Based on this fact, we conclude that PB-PCA cannot be apractical alternative to VB-PCA unfortunately.
We next consider a simpliﬁed variant of VB-PCA, where
the latent variables and loading vectors are mutually independent (cf. the latent variables and loading matrix are
matrix-wise independent in the original VB-PCA). We refer to this method as simple-VB-PCA. For simple-VB-PCA,
we show the following two positive facts:• The global optimal solution of simple-VB-PCA canbe computed analytically (which can be immediately obtained from the result given in Nakajima et al.
(2010)). On the other hand, the ordinary VB-PCA requires a number of iterations to ﬁnd a local optimal
solution (Bishop, 1999a).
• The mutual independence assumption of the latentvariables and loading vectors is not restrictive in practice. More speciﬁcally, we experimentally show thatthe performance of simple-VB-PCA is comparable tothat of the original VB-PCA.
Based on these observations, we conclude that simple-VBPCA is the most attractive as a Bayesian PCA method.
This paper is organized as follows. In Section 2, we formulate PPCA, and introduce approximation methods tothe fully Bayesian inference. In Section 3, we derive they = Bea + ε,
where B ∈ RL×H speciﬁes the linear relationship bee and y, and ε ∈ RL is a Gaussian noise subtween aject to NL (0, σ 2 IL ). Here, we denote by Nd (µ, Σ) thed-dimensional Gaussian distribution with mean µ and covariance Σ, and by Id the d-dimensional identity matrix.
e is subject to NH (0, IH ).
We assume that the latent vector a
Suppose that we are given M observed samples{y 1 , . . . , y M }, which are generated from the latent vectorse M }. Deﬁne the following matrices:{e
a1 , . . . , ae M ) ∈ RH×M.
Y=(y 1 , . . . , y M ) ∈ RL×M, A⊤=(ea1 , . . . , aThen the PPCA model is written as follows1 :()
1p(Y |A, B) ∝ exp − 2 ∥Y − BA⊤ ∥2Fro ,2σ()
1
φA (A) ∝ exp − ∥A∥2Fro ,
2
(
)
)
1 (
−1 ⊤
φB (B) ∝ exp − tr BCB
B
,
2

(1)
(2)
(3)

where ∥ · ∥Fro and tr(·) denote the Frobenius norm and thetrace of a matrix, respectively, and CB ∈ RH×H is a diagonal matrix with positive entries. The column vectors of
B = (b1 , . . . , bH ) correspond to the loading vectors, andthe diagonal elements of CB = diag(c2b1 , . . . , c2bH ) correspond to the prior variances of the loading vectors. Without loss of generality, we assume that {cbh } are arranged in
non-increasing order.
Note that the above PPCA model can be seen as Bayesianmatrix factorization (Salakhutdinov & Mnih, 2008), if U =BA⊤ is regarded as a low rank matrix approximating Y .
1

Note that, for controlling the regularization effect and introducing an empirical Bayesian variant, we added the prior (3) on
B to the original PPCA formulation given in Tipping & Bishop(1999). When the diagonal elements of CB tend to inﬁnity, themodel (1)–(3) is reduced to the original formulation.

On Bayesian PCA: Automatic Dimensionality Selection and Analytic Solution

Throughout the paper, we denote a column vector of a matrix by a bold smaller letter, and a row vector by a boldsmaller letter with a tilde, namely,

be the singular value decomposition of Y . In many methods, the PCA solution can be written in the following form:⊤

S = spane L ) ∈ RL×M ,Y = (y 1 , . . . , y M ) = (ey1 , . . . , y
⊤

,

(8)

where θ(·) is the function taking one if the event is true andzero otherwise, and γ h is a threshold.
Below, we review two types of approximation methods.

2.2. Approximate Bayesian Inference

2.2.2. PARTIALLY BAYESIAN PCA (PB-PCA)

The Bayes posterior is given byp(Y |A, B)φA (A)φB (B)
,
Z(Y )

)
θ(γh > γ h )ω bh ω ⊤ahh=1e M ) ∈ RM ×H ,A = (a1 , . . . , aH ) = (ea1 , . . . , a)⊤(∈ RL×H .
b1 , . . . , ebL
B = (b1 , . . . , bH ) = ep(A, B|Y ) =

(H
∑

(4)

where Z(Y ) = 〈p(Y |A, B)〉φA (A)φB (B) . Here, 〈·〉p denotes
the expectation over the distribution p. Since this expectation is intractable, it needs to be approximated. Here we
review methods of approximate Bayesian inference.
2.2.1. F REE E NERGY M INIMIZATIONFirst, we describe an approximation framework.
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)F (r|Y ) = − log(5)p(Y |A, B)φA (A)φB (B) r(A,B)〈〉
r(A, B)= − log− log Z(Y ).
p(A, B|Y ) r(A,B)
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 (5) amounts to ﬁnding a distribution
closest to the Bayes posterior in the sense of the KL distance.

In the ﬁrst formulation of PPCA by Tipping & Bishop(1999), only the latent matrix A is integrated out, and theloading matrix B is point-estimated. This amounts to restricting the posterior to the following form:rPB-A (A, B) = rAPB-A (A)δ(B; B ∗ ),(9)where δ(B; B ∗ ) denotes a (pseudo-)Dirac delta function of
B located at B = B ∗ .2
We refer to this method as PB-A-PCA since A is integrated out. An analytic-form solution for PB-A-PCA when{cbh → ∞} is given as follows (this proposition is obtainedby combing Eqs.(6) and (7) in Tipping & Bishop (1999).):Proposition 1 (Tipping & Bishop, 1999) The PB-A-PCA
solution when {cbh → ∞} is given by Eq.(8) with the following threshold:√
γ PB-A= σ M.
h(10)

In the above formulation, A is integrated out and B is pointestimated. On the other hand, we can naturally think of the
opposite variant, i.e., B is integrated out and A is pointestimated (which we refer to as PB-B-PCA):

A general approach to Bayesian approximate inference is
to ﬁnd the minimizer of the free energy (5) with respect tor in some restricted function space.

rPB-B (A, B) = δ(A; A∗ )rBPB-B (B).

Let rb be such a minimizer. In the context of PCA, the solution is the subspace spanned by the estimated loading vectors:()(6)
S = span 〈B〉rb(A,B) ,In this paper, we slightly extend the formulation by
Tipping & Bishop (1999) and allow the posterior to be chosen from (9) or (11) that gives a smaller free energy. We callthis adaptive method partially Bayesian PCA (PB-PCA).

where span(·) denotes the subspace spanned by the column
vectors of a matrix. Let
Y =

H
∑
h=1

γh ω bh ω ⊤
ah

(7)

(11)

In Section 3, we derive an analytic-form solution of PBPCA for arbitrary {cbh }.
2

By a pseudo-Dirac delta function, we mean”
“ an extremely

localized density function, e.g., δ(B; B ∗ ) ∝ exp −∥B−B ∗ ∥2
Fro2c2

with a very small but strictly positive variance c2 , such that its taileffect can be neglected, while χB = −〈log δ(B; B ∗ )〉δ(B;B ∗ )remains ﬁnite.

On Bayesian PCA: Automatic Dimensionality Selection and Analytic Solution

2.2.3. VARIATIONAL BAYESIAN PCA (VB-PCA)

3.1. Analytic-Form Solution for PB-PCA

In the VB approximation, the independence between theentangled parameter matrices A and B is assumed:

Here, we derive an analytic-form solution for PB-PCA anddiscuss its properties.

rVB (A, B) = rAVB (A)rBVB (B).

We ﬁrst consider PB-A-PCA. Substituting Eq.(9) into
Eq.(5), we have〈〉rAPB-A (A)PB-A ∗F (rA , B |Y ) = − log+χBp(Y |A,B ∗ )φA (A)φB (B ∗ ) rPB-A (A)A〈〉
rAPB-A (A)= − log− log Z(Y |B ∗ )φB (B ∗ ) + χB ,p(A|Y, B ∗ ) rPB-A (A)
A
(17)

(12)

With this constraint, an iterative algorithm for minimizing the free energy (5) was derived (Bishop, 1999a;
Lim & Teh, 2007).
The VB posterior can be written asrVB (A, B) =

M
∏

m=1

NH (e
am ; µeam , ΣAe)

L
∏
NH (e
bl ; µebl , ΣBe ),
l=1

where
p(Y |A, B)φA (A),Z(Y |B)
Z(Y |B) = 〈p(Y |A, B)〉φA (A) ,p(A|Y, B) =where the means and the covariances necessarily satisfy
)
ΣAe (
µ
,
.
.
.
,
µ
e
e
b1
bL y m ,
σ2
)
Σe (
el ,
µea1 , . . . , µeaM yµebl = B
2
σ
)−1
( L
)
∑(
+ ΣBe + σ 2 IH,
ΣAe = σ 2
µebl µe⊤bµeam =

l

(13)
(14)
(15)

l=1

(M
)−1
∑(
)
⊤
2 −1
µeam µeam + ΣAe + σ CB
ΣBe = σ. (16)
2

l=1

We may obtain a local minimizer by iterating (13)–(16).
After convergence, the VB-PCA solution is given as()
S VB = span (µeb1 , . . . , µebL )⊤ .

∗

χB = −〈log δ(B; B )〉δ(B;B ∗ ) .

(18)
(19)
(20)

Note that Eq.(17) is a functional of rAPB-A and B ∗ , and χB
is a constant with respect to them. Since only the ﬁrst
term depends on rAPB-A on which we impose no restriction,
Eq.(17) is minimized whenrbAPB-A (A) = p(A|Y, B ∗ )(21)

for any B ∗ . With Eq.(21), the ﬁrst term in Eq.(17) vanishes,
and thus an estimator for B ∗ is given byb PB-A = argmin F PB-A (B ∗ |Y ),
B
B∗

where
F PB-A (B|Y ) = − log Z(Y |B)φB (B) + χB .

2.3. Empirical Bayesian Procedure

We call Eq.(22) the PB-A free energy.

PPCA has a hyperparameter CB in the prior (3), which controls the sparsity of the result (i.e., the PCA dimensions).
A popular way to set the hyperparameter in the Bayesianframework is again based on the minimization of the freeenergy (5).
()
bB = argmin min F (r; CB |Y ) .
C

Minimizing Eq.(22), we obtain the following solution.

(22)

Theorem 1 The PB-A-PCA solution is given by Eq.(8)with the following threshold:√
γ PB-A=σ
M + σ 2 /c2bh .
(23)hWe refer to this method as an empirical Bayes method.

When cbh → ∞, Eq.(23) converges to Eq.(10). Therefore,
Theorem 1 is an extension of Proposition 1 for arbitrary{cbh }.

3. Behavior of PB-PCA

The PB-B-PCA solution with the constraint (11) is similarly obtained by minimizing the PB-B free energy:

CBrAmong the PPCA methods reviewed in the previous section, VB-PCA seems to be a popular choice. The purposeof this paper is to revisit PPCA, and investigate whether
VB-PCA is really the best choice from the viewpoints ofcomputational efﬁciency and automatic dimensionality selection (ADS).
In this section, we investigate properties of PB-PCA.

F PB-B (A|Y ) = − log Z(Y |A)φA (A) + χA .

(24)

The threshold for PB-B-PCA is similarly given by√
L + σ 2 /c2bh .
γ PB-B=σh(25)

By comparing Eqs.(22) and (24), we obtain the following
lemma:

On Bayesian PCA: Automatic Dimensionality Selection and Analytic Solution

Lemma 1 In PB-PCA, the constraint (9) is chosen when
M > L, and the constraint (11) is chosen when M < L.

4. Behavior of Simpliﬁed VB-PCA(Simple-VB-PCA)

Based on this, we can express the PB-PCA subspace analytically as follows (we can also derive the PB-PCA posterior analytically, but we omit the details due to lack ofspace):In this section, we introduce a simpliﬁed variant of VBPCA, and investigate its properties.

Theorem 2 The PB-PCA solution is given by Eq.(8) withthe following threshold:
√
γ PB
=
σ
max(L, M ) + σ 2 /c2bh .
(26)hIn the context of collaborative ﬁltering, Raiko et al. (2007)
proposed a simple VB iterative algorithm under the following stronger decomposability constraint:

Thanks to Theorem 2, the solution of PB-PCA can be computed analytically in a very efﬁcient way. On the other
hand, VB-PCA requires a number of iterations to ﬁnd a
local optimal solution (see Section 2.2.3). Thus, PB-PCA
is computationally more attractive than VB-PCA.
Theorem 2 also shows that PB-PCA has a thresholding effect, i.e., the components with singular values smaller thanγ PB
are eliminated. Moreover, this thresholding effect reh
mains even when the prior is ﬂat (cbh → ∞).3 Actually,such a thresholding effect is also observed in the existing
work (see Proposition 1) for PB-A-PCA with cbh → ∞.
However, it was regarded as an artifact in the original paper by Tipping & Bishop (1999).
Bishop (1999a) proposed VB-PCA (see Section 2.2.3) forthe purpose of enabling ADS. However, our analysis aboveshows that PB-PCA also has the thresholding effect. We
further show in Section 5.2 that PB-PCA behaves similarly
to a simpliﬁed variant of VB-PCA.
3.2. Empirical PB-PCA
A critical drawback of PB-PCA appears when the hyperparameter CB is estimated from data. Speciﬁcally, the empirical Bayesian procedure for PB-PCA always results inthe trivial solution, i.e., cbh → 0. This is because the ﬁrst
term of the free energy (22), as well as (24), is not lowerbounded (i.e., it tends to minus inﬁnity). More generally, inorder to obtain a meaningful solution when a hyperparameter is estimated through the empirical Bayes procedure, thecorresponding parameter must be integrated out.
Consequently, despite the existence of the ADS effect, PBPCA cannot be a practical alternative to VB-PCA unfortunately.
3This applies only to the case when the noise variance σ 2is strictly positive. When σ 2 is unkown, the free energy minimization fails to estimate it in PB-PCA with cbh → ∞, because σ 2 → 0 is the global minimizer (see Appendix A.2 inTipping & Bishop (1999)). We would say that it is essential for
inducing ADS to correctly estimate σ 2 when it is unknown. Whencbh is ﬁnite, σ 2 is estimated to be positive, although we experimentally found that it tends to be underestimated.

4.1. Analytic-Form Solution for simple-VB-PCA

rsimple-VB (A, B) =H∏

rasimple-VB(ah )rbsimple-VB(bh ). (27)
h
h

h=1

That is, all column-vectors of A and B are assumed to
be mutually independent. Under this stronger constraint,
Nakajima et al. (2010) derived the VB global solution in aclosed-form.
From their result, the following theorem can be immediately obtained.
Theorem 3 The simple-VB-PCA solution is given by
Eq.(8) with the following threshold:√√simple-VBγh= σ κh + κ2h − LM ,(28)κh = (L + M )/2 + σ 2 /(2c2bh ).
Theorem 3 shows that the solution of simple-VB-PCA canbe computed analytically in a very efﬁcient way. On theother hand, the ordinary VB-PCA (with matrix-wise independence) requires a number of iterations to ﬁnd a local optimal solution (see Section 2.2.3). Thus, simple-VB-PCA
is computationally more attractive than the ordinary VBPCA.
4.2. Simple Empirical VB-PCA (Simple-EVB-PCA)
As opposed to PB-PCA, the empirical Bayesian variant ofsimple-VB-PCA is well-deﬁned, and its global solution canbe obtained in a closed form as follows:
Theorem 4 The simple-EVB-PCA solution is given by
Eq.(8) with the following threshold:{ √√σ( L + M ) if ∆h < 0,simple-EVB
γh
=
(29)
∞
otherwise,
where

( γ
)
( γ
)
h
h
simple-VB
simple-VB∆h = M log
γ̆
+
1
+
L
log
γ̆
+
1
h
h
M σ2
Lσ 2
)
1 (
simple-VB+ 2 LM c̆2bh − 2γh γ̆h,σ √τh + τh2 − 4LM σ 4c̆2bh =, τh = γh2 − (L + M )σ 2 .
2LM

On Bayesian PCA: Automatic Dimensionality Selection and Analytic Solution
M = 300, L = 100

60

1.2

EVB-PCA
Simple-EVB-PCA
Running times (s)Estimated dimensions Ĥ

50

40

30

20

10

15

M = 300, L = 100

1.4

Table 1. Estimated effective data dimensions of real datasets.
b EVB and Hb simple-EVB denote the PCA dimensions estimated by
H
EVB-PCA and simple-EVB-PCA, respectively.

EVB-PCA
Simple-EVB-PCA

1

0.8

0.6

25

30

35

40

45

50

M

L

b EVB
H

b simple-EVB
H

Chart
Glass
Wine
Optical Digits
Satellite
Segmentation
Letter

600
214
178
5620
6435
2310
20000

60
9
13
64
36
19
16

11
7
7
56
32
12
15

10
7
8
56
31
13
15

0.4

0.2

20

Data set0
15

True dimensions H ∗

20

25

30

35

40

45

50

True dimensions H ∗Figure 1. EVB-PCA vs. simple-EVB-PCA. Performance of the
PCA dimensionality selection is highly comparable to each other(left), while simple-EVB-PCA is computationally more efﬁcientthan EVB-PCA (right).

γ̆hsimple-VB = 〈∥ah ∥∥bh ∥〉rbsimple-VB (A,B) is the simple-VB solution for cbh = c̆bh .

1 in the ﬁrst H ∗ directions, and variances 0.1 and no covariance in the remaining L − H ∗ directions.

Below, we experimentally show that the above analyticform solution of simple-VB-PCA works as well as the ordinary VB-PCA.

Figure 1 shows the estimated PCA dimensions with varying
true dimension H ∗ . It can be observed that EVB-PCA andsimple-EVB-PCA have very similar dimension estimation
performance, while simple-EVB-PCA is computationallymore efﬁcient than EVB-PCA4 . We also performed similar
experiments with various settings (e.g., different levels ofcorrelation, different number of samples M , and different
dimensionality L), and empirically observed similar trends
to Figure 1.5

4.3. Experimental Comparison
In this section, we experimentally compare the performance of EVB-PCA and simple-EVB-PCA on artiﬁcial
and real-world datasets. Through all the experiments, H
is set to be the full rank (H = min(L, M )), and thenoise variance σ 2 is estimated based on the free energy
minimization, similarly to the hyperparameter estimation
through the empirical Bayes method (see Section 2.3).
In EVB-PCA, hyperparameters are updated in each iteration as( )c2bh = ∥µbh ∥2 /L + ΣBe hh ,(30)which is obtained as a stationary point of the free energy
with respect to {c2bh } (Bishop, 1999a). The noise varianceσ 2 is estimated by using the following update rule in each
iteration
∑L
∑H
⊤
⊤ 2
Σ eµeblµµ∥∥Y−
l=1 µebl A
h=1 bh ah Fro
2
σ =
+
LM
L
∑M
⊤
(
)
µ Σ e µea
+ m=1 eam B m + tr ΣAeΣBe ,(31)
Mwhich is obtained again as a stationary point of the freeenergy with respect to σ 2 (Bishop, 1999a).
In simple-EVB-PCA, we use Theorem 4 combined with anaive 1-dimensional search strategy over the free energy
for the estimation of σ 2 (Nakajima et al., 2010).
The strong independence assumption of simple-VB-PCA
given in Eq.(27) naturally leads to the question, whether
simple-VB-PCA can provide accurate estimation when thedata has a correlated structure. To investigate this, we created an artiﬁcial dataset having variances 5 and covariances

We we further tested both methods on seven datasets takenfrom the UCI repository (Asuncion & Newman, 2007).
The speciﬁcation of the dataset as well as the number of
PCA dimensions retained by EVB-PCA and simple-EVBPCA is shown in Table 1. As with the synthetic data, thebehavior of the two methods is very similar to each other.
In summary, we observed that the stronger independence
assumption (27) does not signiﬁcantly change the performance. Given that simple-EVB-PCA and EVB-PCA perform similarly, we conclude that simple-EVB-PCA is more
attractive than EVB-PCA because it has an analytic-form
solution that can be computed very efﬁciently.

5. Discussion
In this section, we further discuss various issues in several
PPCA methods, and give insights.
5.1. Maximum A Posteriori PCA (MAP-PCA)
For comparison purposes, we deﬁne the maximum a posteriori PCA (MAP-PCA), where both A and B are point4

Note that slightly different hyperpriors are used in the original EVB-PCA papers (Bishop, 1999a;b). We also tested themand conﬁrmed that they also behave similarly to the simple-EVBPCA.
5
We also observed that the iterative algorithm (Eqs.(2) and(3) in Nakajima et al. (2010)) for simple-EVB-PCA gives similardimension estimation performance.

On Bayesian PCA: Automatic Dimensionality Selection and Analytic Solution
L=20,M =50

L=20,M =1015

simple-VB
PB
PB-A
MAP

10

ThresholdThreshold

15

5

0
0

simple-VB
PB
PB-A
MAP

10

5

0.5

1
cb h

1.5

0
0

2

0.5

1
cb h

1.5

2

Figure 2. Behavior of the thresholds, γ simple-VB, γ PB, γ PB-A
, and
h
h
h
MAP
2
γ h . The noise variance is set to σ = 1.
L=20,c b h = ∞
10

6
4
2
0
0

simple-VB
PB
PB-A
MAP

8
Threshold

8
Threshold

L=20,c b h=1.010

simple-VB
PB
PB-A
MAP

6
4
2

10

20
M

30

40

0
0

10

20
M

30

40

Figure 3. Behavior of the thresholds when the number M of samples is changed.

estimated. This corresponds to the posterior restricted tothe product of delta functions:rMAP (A, B) = δ(A; A∗ )δ(B; B ∗ ).

(32)

Its solution is given as follows.
Proposition 2 (Srebro et al., 2005) The MAP-PCA solution is given by Eq.(8) with the following threshold:γ MAP= σ 2 /cbh .
h(33)

When cbh → ∞, MAP-PCA is reduced to the classical
PCA, i.e., no automatic dimensionality selection (ADS).
5.2. Comparison between PB-A-PCA, PB-PCA, andsimple-VB-PCA
The threshold (26) becomes simpler when the prior is ﬂat:√
lim γ PB=σmax(L, M ).
(34)hcbh →∞

Interestingly, the simple-VB-PCA solution agrees with the
PB-PCA solution under the ﬂat prior:√
lim γ simple-VB= σ max(L, M ).
(35)hcbh →∞

Below, we numerically investigate the behavior of the
PPCA thresholds in more general settings.
Figure 2 compares the thresholds of simple-VB-PCA (28),
PB-PCA (26), PB-A-PCA (23), and MAP-PCA (33) in
two situations. In the left graph (when the dimensionality of the original space is L = 20 and the number ofsamples is M = 50), PB-PCA and PB-A-PCA coincide,and simple-VB-PCA behaves similarly to them. On theother hand, MAP-PCA behaves differently. This is because of the effect of model-induced regularization (MIR)(Nakajima et al., 2010), which indicates the fact that, whenthe model is non-identiﬁable, even if the prior is almost
ﬂat (cbh → ∞), Bayesian methods have a regularization
effect as long as (a part of) parameters are integrated out.
PB-PCA and VB-PCA belong to this class, but MAP-PCA
does not.
In the right graph (when L = 20 and M = 10), PB-PCA
and simple-VB-PCA behave almost identically. Although
PB-A-PCA also shows its MIR effect (i.e., the threshold
do not converge to zero even in the limit cbh → ∞), asigniﬁcant difference from simple-VB-PCA and PB-PCA
is observed.
Figure 3 shows how the number M of samples affects thethresholds. In the limit cbh → ∞ (left graph), Eqs.(34)
and (35) together state that the solutions of PB-PCA andsimple-VB-PCA agrees with each other. PB-A-PCA also
gives the identical solution to them when M ≥ L. However, its behavior changes at M = L; the threshold of PBA-PCA smoothly goes down as M decreases, while thoseof PB-PCA and simple-VB-PCA make a sudden turn andbecomes constant. The right graph in Figure 3 shows thecase with a non-ﬂat prior (cbh = 1), which also shows asimilar phenomenon.
A question is which is more desirable, simple-VBPCA/PB-PCA (which are with a sudden turn in the threshold curve), or PB-A-PCA (which is with smooth behavior).
When M < L, PB-PCA and VB-PCA more strongly regularize the solutions than PB-A-PCA. We argue that the
behavior of PB-PCA and VB-PCA is more reasonable because of the following reason. Let us consider the casewhere no driving latent variable exists, i.e., the true dimension is H ∗ = 0. In this case, we merely observe pure noise,and the average of the squared singular values of Y over all
the components is given by〈tr(Y Y ⊤ )〉N (0,σ2 )= σ 2 max(L, M ).
min(L, M )Comparing this with the thresholds (34) and (35) of PBPCA and simple-VB-PCA, we ﬁnd that they eliminate thecomponents with singular values equal to the average noise
contribution. The sudden turn in the threshold curve actually follows the behavior of the noise contribution, which
would be reasonable. In this sense, PB-A-PCA is suboptimal since it strongly overﬁts the noise when L ≫ M .
5.3. Changing Behavior of Simple-VB-PCA
In the left graph of Figure 3, the behavior of PB-PCA andsimple-VB-PCA suddenly changes at L = M . Here, wetheoretically investigate the reason for this phenomenon.

On Bayesian PCA: Automatic Dimensionality Selection and Analytic Solution
Table 2. Properties of approximate PPCA methods. Upper methods have weaker constraints on the posterior distribution. ‘⃝’means the method positively possesses the corresponding property, ‘△’ means weakly, and ‘×’ means not.
MethodAutomatic
dimensionality
selection

Empirical
Bayes

Analytic-form
solution

⃝
⃝
⃝
△
×

⃝
⃝
×
×
×

×
⃝
⃝
⃝
⃝

VB
simple-VB
PB
PB-A (PB-B)
MAP

solution (Section 4.1), and its empirical Bayesian version,
simple-EVB-PCA, also has an analytic-form solution (Section 4.2). Experimentally, simple-EVB-PCA was shown toperform similarly to the computationally more demanding
counterpart, EVB-PCA (Section 4.3). Based on these ﬁndings, we concluded that simple-VB-PCA is the most attractive PPCA method.
Our future work is to extend the current discussion to more
complex models such as a mixture of PPCAs and missing
value estimation. Also, we will further investigate the roleof column-wise independence in simple-VB-PCA.

Acknowledgments
As Lemma 1 states, the PB-PCA posterior drasticallychanges depending on L > M or L < M . This leads to thesudden turn in the thresholding curve. We can show that asimilar effect also occurs in simple-VB-PCA, which is explained below. Let σasimple-VB
and σbsimple-VB
be the standardhhdeviations of the simple-VB posteriors of ∥ah ∥ and ∥bh ∥,respectively. Then the following lemma holds:
Lemma 2 When cbh → ∞ and γh = γ simple-VB, it holdshthat (M −L)√+ o (1) = O (1)if M > L, σ M
simple-VB
σah
−1
−1
1
o(cbh ) = O(cbh )if M = L,= cbh +√
σbsimple-VB σ L 2 + o(c−2 ) = O(c−2 ) if M < L.
h
(L−M )cb

h

bh

bh

This lemma implies that with the (almost) ﬂat prior, theshape of the simple-VB posterior suddenly changes. Morespeciﬁcally, the shape is spherical when M > L, while itis strongly elliptical in the bh -space when M ≤ L.

6. Conclusion
In this paper, we revisited approximate Bayesian methods
in PPCA, which are summarized in Table 2.
Although VB-PCA is the most general approximate PPCA
method, it is computationally less efﬁcient because it requires a number of iterations to ﬁnd a local optimal solution. On the other hand, we showed by providing
an analytic-form solution that PB-PCA is computationally more efﬁcient (Section 3.1). The analytic-form solution also showed that, despite the fact that VB-PCA
was proposed to induce automatic dimensionality selection (ADS), PB-PCA is already equipped with ADS. We
also numerically showed that PB-PCA behaves similarly tosimple-VB-PCA (Section 5.2). However, as shown in Section 3.2, PB-PCA has a critical disadvantage that its empirical Bayesian variant is practically useless.
On the other hand, simple-VB-PCA (whose ‘complexity’ isbetween PB-PCA and VB-PCA) still has an analytic-form

The authors thank anonymous reviewers for helpful comments to improve the paper. They also thank Ryota
Tomioka of The University of Tokyo for discussion. MS
was supported by the FIRST program. DB was supported
by a Beckman Postdoctoral Fellowship.

References
Asuncion, A. and Newman, D.J. UCI machine learning repository, 2007.
Bishop, C. M. Variational principal components. ICANN1999,
1999a.
Bishop, C. M. Bayesian principal components. NIPS1998, 1999b.
Hotelling, H. Analysis of a complex of statistical variables into
principal components. Journal of Educational Psychology, 24:
417–441, 1933.
Hoyle, D. C. Autotmaric PCA dimension selection for high dimensional data and small sample sizes. Journal of Machine
Learning Research, 9:2733–2759, 2008.
Lim, Y. J. and Teh, T. W. Variational Bayesian approach to movie
rating prediction. KDD Cup and Workshop 2007.
Nakajima, S., Sugiyama, M., and Tomioka, R. Global analytic solution for variational Bayesian matrix factorization. NIPS2010.
Neal, R. M. Bayesian Learning for Neural Networks. Springer,
1996.
Raiko, T., Ilin, A., and Karhunen, J. Principal component analysis for large scale problems with lots of missing values.
ECML2007.
Roweis, S. and Ghahramani, Z. A unifying review of linear Gaussian models. Neural Computation, 11:305–345, 1999.
Salakhutdinov, R. and Mnih, A. Probabilistic matrix factorization.
NIPS2007.
Srebro, N., Rennie, J., and Jaakkola, T. Maximum margin matrix
factorization. NIPS2004.
Tipping, M. E. and Bishop, C. M. Probabilistic principal component analysis. Journal of the Royal Statistical Society, 61:
611–622, 1999.