Global Analytic Solution
for Variational Bayesian Matrix Factorization

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

Ryota Tomioka
The University of Tokyo
Tokyo 113-8685, Japan
tomioka@mist.i.u-tokyo.ac.jp

Abstract
Bayesian methods of matrix factorization (MF) have been actively explored recently as promising alternatives to classical singular value decomposition. In thispaper, we show that, despite the fact that the optimization problem is non-convex,
the global optimal solution of variational Bayesian (VB) MF can be computed
analytically by solving a quartic equation. This is highly advantageous over a
popular VBMF algorithm based on iterated conditional modes since it can only
ﬁnd a local optimal solution after iterations. We further show that the global optimal solution of empirical VBMF (hyperparameters are also learned from data) can
also be analytically computed. We illustrate the usefulness of our results throughexperiments.

1 Introduction
The problem of ﬁnding a low-rank approximation of a target matrix through matrix factorization(MF) attracted considerable attention recently since it can be used for various purposes such as
reduced rank regression [19], canonical correlation analysis [8], partial least-squares [27, 21],
multi-class classiﬁcation [1], and multi-task learning [7, 29].
Singular value decomposition (SVD) is a classical method for MF, which gives the optimal lowrank approximation to the target matrix in terms of the squared error. Regularized variants of SVD
have been studied for the Frobenius-norm penalty (i.e., singular values are regularized by the ℓ2 penalty) [17] or the trace-norm penalty (i.e., singular values are regularized by the ℓ1 -penalty) [23].
Since the Frobenius-norm penalty does not automatically produce a low-rank solution, it should be
combined with an explicit low-rank constraint, which is non-convex. In contrast, the trace-norm
penalty tends to produce sparse solutions, so a low-rank solution can be obtained without explicit
rank constraints. This implies that the optimization problem of trace-norm MF is still convex, andthus the global optimal solution can be obtained. Recently, optimization techniques for trace-norm
MF have been extensively studied [20, 6, 12, 25].
Bayesian approaches to MF have also been actively explored. A maximum a posteriori (MAP)
estimation, which computes the mode of the posterior distributions, was shown [23] to correspond to
the ℓ1 -MF when Gaussian priors are imposed on factorized matrices [22]. The variational Bayesian(VB) method [3, 5], which approximates the posterior distributions by factorized distributions, hasalso been applied to MF [13, 18]. The VB-based MF method (VBMF) was shown to perform wellin experiments, and its theoretical properties have been investigated [15].
1

M
L

U

H
=L

B

M
A⊤

H

Figure 1: Matrix factorization model. H ≤ L ≤ M . A = (a1 , . . . , aH ) and B = (b1 , . . . , bH ).
However, the optimization problem of VBMF is non-convex. In practice, the VBMF solution iscomputed by the iterated conditional modes (ICM) [4, 5], where the mean and the covariance of theposterior distributions are iteratively updated until convergence [13, 18]. One may obtain a local
optimal solution by the ICM algorithm, but many restarts would be necessary to ﬁnd a good local
optimum.
In this paper, we ﬁrst show that, although the optimization problem is non-convex, the global optimal solution of VBMF can be computed analytically by solving a quartic equation. This is highly
advantageous over the standard ICM algorithm since the global optimum can be found without any
iterations and restarts. We next consider an empirical VB (EVB) scenario where the hyperparameters (prior variances) are also learned from data. Again, the optimization problem of EVBMF is
non-convex, but we still show that the global optimal solution of EVBMF can be computed analytically. The usefulness of our results is demonstrated through experiments.
Recently, the global optimal solution of VBMF when the target matrix is square has been obtainedin [15]. Thus, our contribution to VBMF can be regarded as an extension of the previous result togeneral rectangular matrices. On the other hand, for EVBMF, this is the ﬁrst paper that gives theanalytic global solution, to the best of our knowledge. The global analytic solution for EVBMF isshown to be highly useful in experiments.

2

Bayesian Matrix Factorization

In this section, we formulate the MF problem and review a variational Bayesian MF algorithm.
2.1

Formulation

The goal of MF is to approximate an unknown target matrix U (∈ RL×M ) from its n observations
V n = {V (i) ∈ RL×M }ni=1 .
We assume that L ≤ M . If L > M , we may simply re-deﬁne the transpose U ⊤ as U so that L ≤ M
holds. Thus this does not impose any restriction.
A key assumption of MF is that U is a low-rank matrix. Let H (≤ L) be the rank of U . Thenthe matrix U can be decomposed into the product of A ∈ RM ×H and B ∈ RL×H as follows (seeFigure 1):
U = BA⊤ .
Assume that the observed matrix V is subject to the following additive-noise model:
V = U + E,
where E (∈ RL×M ) is a noise matrix. Each entry of E is assumed to independently follow theGaussian distribution with mean zero and variance σ 2 . Then, the likelihood p(V n |A, B) is given by
Ã
!
n
X
1
p(V n |A, B) ∝ exp − 2∥V (i) − BA⊤ ∥2Fro ,2σ i=1
where ∥ · ∥Fro denotes the Frobenius norm of a matrix.
2

2.2

Variational Bayesian Matrix Factorization

We use the Gaussian priors on the parameters A = (a1 , . . . , aH ) and B = (b1 , . . . , bH ):
Ã H
!
Ã H
!
X ∥ah ∥2
X ∥bh ∥2
φ(U ) = φA (A)φB (B), where φA (A) ∝ exp −and φB (B) ∝ exp −
.
2c2ah
2c2bh
h=1

h=1

c2ah and c2bh are hyperparameters corresponding to the prior variance. Without loss of generality, we
assume that the product cah cbh is non-increasing with respect to h.
Let r(A, B|V n ) be a trial distribution for A and B, and let FVB be the variational Bayes (VB) freeenergy with respect to r(A, B|V n ):À¿r(A, B|V n ),FVB (r|V n ) = logp(V n , A, B) r(A,B|V n )where 〈·〉p denotes the expectation over p.
The VB approach minimizes the VB free energy FVB (r|V n ) with respect to the trial distributionr(A, B|V n ), by restricting the search space of r(A, B|V n ) so that the minimization is computationally tractable. Typically, dissolution of probabilistic dependency between entangled parameters (Aand B in the case of MF) makes the calculation feasible:1
H
Yr(A, B|V n ) =rah (ah |V n )rbh (bh |V n ).
(1)
h=1

b VB is given by the VB
The resulting distribution is called the VB posterior. The VB solution U
posterior mean:b VB = 〈BA⊤ 〉r(A,B|V n ) .
U
By applying the variational method to the VB free energy, we see that the VB posterior can be
expressed as follows:H
Yr(A, B|V n ) =NM (ah ; µah , Σah )NL (bh ; µbh , Σbh ),h=1where Nd (·; µ, Σ) denotes the d-dimensional Gaussian density with mean µ and covariance matrix
Σ. µah , µbh , Σah , and Σbh satisfy
³ nβ
´−1
³ nα
´−1
h
h
µah = Σah Ξ⊤+c−2
IM , Σbh =+c−2
IL , (2)h µbh , µbh = Σbh Ξh µah , Σah =
ah
bh
2
2
σ
σ
where Id denotes the d-dimensional identity matrix, andαh = ∥µah ∥2 + tr(Σah ), βh = ∥µbh ∥2 + tr(Σbh ),n´
Xn³
1 X (i)Ξh = 2 V −,µbh′ µ⊤
V
=
V .
ah′
σ
n i=1
′
h ̸=h

The iterated conditional modes (ICM) algorithm [4, 5] for VBMF (VB-ICM) iteratively updatesµah , µbh , Σah , and Σbh by Eq.(2) from some initial values until convergence [13, 18], allowing one
to obtain a local optimal solution. Finally, an estimator of U is computed as
HX
VB−ICMbU=µbh µ⊤ah .
h=1
2

When the noise variance σ is unknown, it may be estimated by the following re-estimation formula:


°2
°
n °
H ³
H
°
´
X
X
X
1
1
°
°

σ2 = 2µbh µ⊤+
αh βh − ∥µah ∥2 ∥µbh ∥2  ,°V (i) −ah °°°σ LM n
i=1

h=1

Fro

h=1

which corresponds to the derivative of the VB free energy with respect to σ 2 set to zero (see Eq.(4)in Section 3). This can be incorporated in the ICM algorithm by updating σ 2 from some initial value
by the above formula in every iteration of the ICM algorithm.
1
Although a weaker constraint, r(A, B|V n ) = rA (A|V n )rB (B|V n ), is sufﬁcient to derive a tractable iterative algorithm [13], we assume the stronger one (1) used in [18], which makes our theoretical analysis tractable.

3

2.3

Empirical Variational Bayesian Matrix Factorization

In the VB framework, hyperparameters (c2ah and c2bh in the current setup) can also be learned fromdata by minimizing the VB free energy, which is called the empirical VB (EVB) method [5].
By setting the derivatives of the VB free energy with respect to c2ah and c2bh to zero, the following
optimality condition can be obtained (see also Eq.(4) in Section 3):
c2ah = αh /M and c2bh = βh /L.

(3)

The ICM algorithm for EVBMF (EVB-ICM) is to iteratively update c2ah and c2bh by Eq.(3), in addition to µah , µbh , Σah , and Σbh by Eq.(2). Again, one may obtain a local optimal solution by thisalgorithm.

3

Analytic-form Expression of Global Optimal Solution of VBMF

In this section, we derive an analytic-form expression of the VBMF global solution.
The VB free energy can be explicitly expressed as follows.
!
Ã
H
X
1
1
αh
L
βh
nLM
M
2
2
2
n
log σ +
log cah− log |Σah |+ 2 + log cbh− log |Σbh |+ 2
FVB (r|V ) =2222cah 2
2
2cbh
h=1
°
°
2
n
H
H
°
´
1 X°
n X³
° (i) X
°
2
2
α
β
−
∥µ
∥
∥µ
∥
,
(4)
+ 2
µbh µ⊤
+
°V −
°
h h
ah
bh
ah
°
2σ i=1 °
2σ 2
h=1

Fro

h=1

where | · | denotes the determinant of a matrix. We solve the following problem:Given(c2ah , c2bh ) ∈ R2++ (∀h = 1, . . . , H), σ 2 ∈ R++ ,min FVB ({µah , µbh , Σah , Σbh ; h = 1, . . . , H})s.t.

Lµah ∈ RM , µbh ∈ RL , Σah ∈ SM++ , Σbh ∈ S++ (∀h = 1, . . . , H),

where Sd++ denotes the set of d × d symmetric positive-deﬁnite matrices. This is a non-convex
optimization problem, but still we show that the global optimal solution can be analytically obtained.
Let γh (≥ 0) be the h-th largest singular value of V , and let ω ah and ω bh be the associated rightand left singular vectors:2
LX
V =γ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,(5)where the coefﬁcients are deﬁned byÃ!
(L − M )2 γh(L2 + M 2 )b
ηh2
2σ 4
ξ3 =
, ξ2 = − ξ3 γh ++ 2 2 2,LMLMn cah cbh
!2
Ã
µ
¶µ
¶
4
σ2 L
σ2 M
σ
2
2
, ηbh = 1 −1−γh2 .
ξ0 = ηbh − 2 2 2n cah cbh
nγh2
nγh2

Let

ξ1 = ξ3

p
ξ0 ,

v
vÃ
u
!2
u
u
u (L + M )σ 2
44u (L + M )σ 2σ
σ
LM σ 4tγeh = t+ 2 2 2 ++ 2 2 2.
−2n
2n cah cbh
2n2n cah cbhn2(6)b VB as in the following theorem.
Then we can analytically express the VBMF solution U
2In our analysis, we assume that V has no missing entry, and its singular value decomposition (SVD) iseasily obtained. Therefore, our results cannot be directly applied to missing entry prediction.

4

Theorem 1 The global VB solution can be expressed asb VB =UH
X

γbhVB ω bh ω ⊤bhVB =
ah , where γ

h=1

½
γ
bh
0

if γh > γeh ,otherwise.

Sketch of proof: We ﬁrst show that minimizing (4) amounts to a reweighed SVD and any minimizer
is a stationary point. Then, by analyzing the stationary condition (2), we obtain an equation withrespect to γbh as a necessary and sufﬁcient condition to be a stationary point (note that its quadratic
approximation gives bounds of the solution [15]). Its rigorous evaluation results in the quartic equation (5). Finally, we show that only the second largest solution of the quartic equation (5) lies within
the bounds, which completes the proof.
The coefﬁcients of the quartic equation (5) are analytic, so γ
bh can also be obtained analytically3 ,e.g., by Ferrari’s method [9] (we omit the details due to lack of space). Therefore, the global VB
solution can be analytically computed. This is a strong advantage over the standard ICM algorithm
since many iterations and restarts would be necessary to ﬁnd a good solution by ICM.
Based on the above result, the complete VB posterior can also be obtained analytically as follows.
Corollary 2 The VB posteriors are given byrA (A|V n ) =HY

NM (ah ; µah , Σah ),rB (B|V n ) =

h=1

H
Y

NM (bh ; µbh , Σbh ),h=1where, for γ
bhVB being the solution given by Theorem 1,qqVBbbh δh · ω ah , µbh = ± γ
bhVB δbh−1 · ω bh ,µah = ± γ
Ã ¡
!
p
¢
− nb
ηh2 − σ 2 (M − L) + (nbηh2 − σ 2 (M − L))2 + 4M nσ 2 ηbh2
Σah =
IM ,
2nM (bγhVB δbh−1 + n−1 σ 2 c−2
ah )
!
Ã ¡
¢ p
ηh2 + σ 2 (M − L))2 + 4Lnσ 2 ηbh2
− nb
ηh2 + σ 2 (M − L) + (nb
IL ,
Σbh =
2nL(b
γhVB δbh + n−1 σ 2 c−2
bh )
q
4 LM
n(M − L)(γh − γ
bhVB ) + n2 (M − L)2 (γh − γ
bhVB )2 + 4σ
c2a c2bhhδbh =,2σ 2 M c−2ah(
ηh2
if γh > γ
eh ,
2
ηbh =
σ2
otherwise.
nca cbhh

When the noise variance σ 2 is unknown, one may use the minimizer of the VB free energy withrespect to σ 2 as its estimate. In practice, this single-parameter minimization may be carried out
numerically based on Eq.(4) and Corollary 2.

4Analytic-form Expression of Global Optimal Solution of Empirical VBMF

In this section, we solve the following problem to obtain the EVBMF global solution:Given

σ 2 ∈ R++ ,min FVB ({µah , µbh , Σah , Σbh , c2ah , c2bh ; h = 1, . . . , H})
s.t.

L
2
2
2
µah ∈ RM , µbh ∈ RL , Σah ∈ SM++ , Σbh ∈ S++ , (cah , cbh ) ∈ R++ (∀h = 1, . . . , H),

where Rd++ denotes the set of the d-dimensional vectors with positive elements. We show that, although this is again a non-convex optimization problem, the global optimal solution can be obtained
analytically. We can observe the invariance of the VB free energy (4) under the transform
©
ª
©
ª
−2
−2 2
2
2 2
(µah , µbh , Σah , Σbh , c2ah , c2bh ) → (sh µah , s−1
h µbh , sh Σah , sh Σbh , sh cah , sh cbh )3R.
In practice, one may solve the quartic equation numerically, e.g., by the ‘roots’ function in MATLAB⃝

5

3

5

3.5

Global solution

2.5

4.5

3.25
2

Global solution

0

1

2

4

Global solution
3

0

1

2

3

ch

ch

(a) V = 1.5

0

1

2

3

ch

(b) V = 2.1(c) V = 2.7

Figure 2: Proﬁles of the VB free energy (4) when L = M = H = 1, n = 1, and σ 2 = 1 forobservations V = 1.5, 2.1, and 2.7. (a) When V = 1.5 < 2 = γ h , the VB free energy is monotone
increasing and thus the global solution is given by ch → 0. (b) When V = 2.1 > 2 = γ h , a local
minimum exists at ch = c̆h ≈ 1.37, but ∆h ≈ 0.12 > 0 so ch → 0 is still the global solution. (c)When V = 2.7 > 2 = γ h , ∆h ≈ −0.74 ≤ 0 and thus the minimizer at ch = c̆h ≈ 2.26 is the global
solution.
for any {sh ̸= 0; h = 1, . . . , H}. Accordingly, we ﬁx the ratios to cah /cbh = S > 0, and refer toch := cah cbh also as a hyperparameter.
Let



sµ
¶
2
2 2
4
1
(L
+
M
)σ
(L
+
M
)σ
4LM
σ
γh2 −
,c̆2h =
γh2 −
+
−
2LM
n
n
n2
√
√
√
γ h = ( L + M )σ/ n.

(7)Then, we have the following lemma:
Lemma 3 If γh ≥ γ h , the VB free energy function (4) can have two local minima, namely, ch → 0
and ch = c̆h . Otherwise, ch → 0 is the only local minimum of the VB free energy.
Sketch of proof: Analyzing the region where ch is so small that the VB solution given ch is γbh = 0,
we ﬁnd a local minimum ch → 0. Combining the stationary conditions (2) and (3), we derive aquadratic equation with respect to c2h whose larger solution is given by Eq.(7). Showing that thesmaller solution corresponds to saddle points completes the proof.
Figure 2 shows the proﬁles of the VB free energy (4) when L = M = H = 1, n = 1, and σ 2 = 1
for observations V = 1.5, 2.1, and 2.7. As illustrated, depending on the value of V , either ch → 0
or ch = c̆h is the global solution.
Let

³ nγ

´

³ nγ

´

¢
n ¡
−2γh γ̆hVB + LM c̆2h , (8)2σVBwhere γ̆h is the VB solution for ch = c̆h . We can show that the sign of ∆h corresponds to that ofthe difference of the VB free energy at ch = c̆h and ch → 0. Then, we have the following theorem
and corollary.
∆h := M log

h VBγ̆ + 1
M σ2 h+ L logh VBγ̆ + 1
Lσ 2 h+Theorem 4 The hyperparameter bch that globally minimizes the VB free energy function (4) is givenby b
ch = c̆h if γh > γ h and ∆h ≤ 0. Otherwise b
ch → 0.
Corollary 5 The global EVB solution can be expressed as
(
H
X
γ̆hVB
b EVB =Uγ
bhEVB ω bh ω ⊤bhEVB :=
ah , where γ
0h=1

if γh > γ h and ∆h ≤ 0,otherwise.

Since the optimal hyperparameter value b
ch can be expressed in a closed-form, the global EVB
solution can also be computed analytically using the result given in Section 3. This is again a strong
advantage over the standard ICM algorithm since ICM would require many iterations and restarts toﬁnd a good solution.
6

5

Experiments

In this section, we experimentally evaluate the usefulness of our analytic-form solutions using artiRﬁcial and benchmark datasets. The MATLAB⃝
code will be available at [14].
5.1

Artiﬁcial Dataset

PH ∗∗We randomly created a true matrix V ∗ = h=1 b∗h a∗⊤h with L = 30, M = 100, and H = 10,
where every element of {ah , bh } was drawn independently from the standard Gaussian distribution.
We set n = 1, and an observation matrix V was created by adding independent Gaussian noise withvariance σ 2 = 1 to each element. We used the full-rank model, i.e., H = L = 30. The noise
variance σ 2 was assumed to be unknown, and estimated from data (see Section 2.2 and Section 3).
We ﬁrst investigate the learning curve of the VB free energy over EVB-ICM iterations. We createdthe initial values of the EVB-ICM algorithm as follows: µah and µbh were set to randomly createdorthonormal vectors, Σah and Σbh were set to identity matrices multiplied by scalars σa2h and σb2h ,respectively. σa2h and σb2h as well as the noise variance σ 2 were drawn from the χ2 -distribution with
degree-of-freedom one. 10 learning curves of the VB free energy were plotted in Figures 3(a). The
value of the VB free energy of the global solution computed by our analytic-form solution was alsoplotted in the graph by the dashed line. The graph shows that the EVB-ICM algorithm reduces the
VB free energy reasonably well over iterations. However, for this artiﬁcial dataset, the convergence
speed was quite slow once in 10 runs, which was actually trapped in a local minimum.
Next, we compare the computation time. Figure 3(b) shows the computation time of EVB-ICM overiterations and our analytic form-solution. The computation time of EVB-ICM grows almost linearlywith respect to the number of iterations, and it took 86.6 [sec] for 100 iterations on average. On theother hand, the computation of our analytic-form solution took only 0.055 [sec] on average, including the single-parameter search for σ 2 . Thus, our method provides the reduction of computation
time in 4 orders of magnitude, with better accuracy as a minimizer of the VB free energy.
Next, we investigate the generalization error of the global analytic solutions of VB and EVB, meab − V ∗ ∥2 /(LM ). Figure 3(c) shows the mean and error bars (min and max)
sured by G = ∥U
Froover 10 runs for VB with various hyperparameter values and EVB. A single hyperparameter value
was commonly used (i.e., c1 = · · · = cH ) in VB, while each hyperparameter ch was separately
optimized in EVB. The result shows that EVB gives slightly lower generalization errors than VB
with the best common hyperparameter. Thus, automatic hyperparameter selection of EVB works
quite well.
Figure 3(d) shows the hyperparameter values chosen in EVB sorted in the decreasing order. Thisshows that, for all 10 runs, ch is positive for h ≤ H ∗ (= 10) and zero for h > H ∗ . This implies thatthe effect of automatic relevance determination [16, 5] works excellently for this artiﬁcial dataset.
5.2

Benchmark Dataset

MF can be used for canonical correlation analysis (CCA) [8] and reduced rank regression (RRR)[19] with appropriately pre-whitened data. Here, we solve these tasks by VBMF and evaluate theperformance using the concrete slump test dataset [28] available from the UCI repository [2].
The experimental results are depicted in Figure 4, which is in the same format as Figure 3. The
results showed that similar trends to the artiﬁcial dataset can still be observed for the CCA task withthe benchmark dataset (the RRR results are similar and thus omitted from the ﬁgure). Overall, theproposed global analytic solution is shown to be a useful alternative to the popular ICM algorithm.

6

Discussion and Conclusion

Overcoming the non-convexity of VB methods has been one of the important challenges in theBayesian machine learning community, since it sometimes prevented us from applying the VB methods to highly complex real-world problems. In this paper, we focused on the MF problem with no
missing entry, and showed that this weakness could be overcome by computing the global optimal
solution analytically. We further derived the global optimal solution analytically for the EVBMF
7

0.3

120

EVB-Analytic
EVB-ICM

1.95

EVB-Analytic
EVB-ICM
100

0.28

80

0.26

VB-Analytic
EVB-Analytic

1.6
1.4
1.2

0.24

1.91

40

0.22

1.9

20

0.2

1.89
0

0
0

0.18

ch

60

1

^

1.92

G

Time(sec)

F VB /(LM)1.93

EVB-Analytic

1.8

1.94

0.8
0.6
0.4
0.2

50
Iteration

100(a) VB free energy

50
Iteration0(b) Computation time

0
0

1

10
√ ch

100

10

10

20

30

h

(c) Generalization error(d) Hyperparameter value

Figure 3: Experimental results for artiﬁcial dataset.
1.4

EVB-Analytic
EVB-ICM−63.52200

EVB-Analytic
EVB-ICM

1.2

0.25

VB-Analytic
EVB-Analytic

EVB-Analytic
0.2

1−63.55

ch

0.6

0.1

100

0.4
0.05

0.2

−63.5600.15^−63.54150

0.8
G

Time(sec)−63.5350

100 150
Iteration

200(a) VB free energy

250

0
0

50

50

100
150
Iteration

200

250

(b) Computation time

−3

10

√ ch

10

−1

1

10

(c) Generalization error

0
0

1

2
h

3

(d) Hyperparameter value

Figure 4: Experimental results of CCA for the concrete slump test dataset.
method, where hyperparameters are also optimized based on data samples. Since no hand-tuning
parameter remains in EVBMF, our analytic-form solution is practically useful and computationallyhighly efﬁcient. Numerical experiments showed that the proposed approach is promising.
When cah cbh → ∞, the priors get (almost) ﬂat and the quartic equation (5) is factorized aslim

cah cb →∞“““” ”“ “” ”“ “
” ”“
” ”
2
2
2
2
fh (t) = t + M
1− σnγL2 γh t + 1− σnγM
1− σnγL2 γh = 0.
t − 1− σnγM
t−M
2 γh
2 γh
L
L

h

h

h

h

h

Theorem 1 states that its second largest solution gives the VB estimator for γh > limcah cbh →∞ γeh =p
2
M σ /n. Thus we have
µ µ
¶¶
M σ2
VB
lim γ
bh = max 0, 1 −γh .
cah cbh →∞nγh2This is the positive-part James-Stein (PJS) shrinkage estimator [10], operated on each singular component separately, and this coincides with the upper-bound derived in [15] for arbitrary cah cbh > 0.
The counter-intuitive fact—a shrinkage is observed even in the limit of ﬂat priors—can be explainedby strong non-uniformity of the volume element of the Fisher metric, i.e., the Jeffreys prior [11], inthe parameter space. We call this effect model-induced regularization (MIR), because it is induced
not by priors but by structure of model likelihood functions. MIR was shown to generally appear in
Bayesian estimation when the model is non-identiﬁable (i.e., the mapping between parameters anddistribution functions is not one-to-one) and the parameters are integrated out at least partially [26].
Thus, it never appears in MAP estimation [15]. The probabilistic PCA can be seen as an exampleof MF, where A and B correspond to latent variables and principal axes, respectively [24]. The
MIR effect is observed in its analytic solution when A is integrated out and B is estimated to be themaximizer of the marginal likelihood.
Our results fully made use of the assumptions that the likelihood and priors are both spherical Gaussian, the VB posterior is column-wise independent, and there exists no missing entry. They werenecessary to solve the free energy minimization problem as a reweighted SVD. An important future work is to obtain the analytic global solution under milder assumptions. This will enable us tohandle more challenging problems such as missing entry prediction [23, 20, 6, 13, 18, 22, 12, 25].
Acknowledgments
The authors appreciate comments by anonymous reviewers, which helped improve our earlier
manuscript and suggested promising directions for future work. MS thanks the support from the
FIRST program. RT was partially supported by MEXT Kakenhi 22700138.
8

References[1] Y. Amit, M. Fink, N. Srebro, and S. Ullman. Uncovering shared structures in multiclass classiﬁcation. In
Proceedings of International Conference on Machine Learning, pages 17–24, 2007.
[2] A. Asuncion and D.J. Newman. UCI machine learning repository, 2007.
[3] H. Attias. Inferring parameters and structure of latent variable models by variational Bayes. In Proceedings of the Fifteenth Conference Annual Conference on Uncertainty in Artiﬁcial Intelligence (UAI-99),
pages 21–30, San Francisco, CA, 1999. Morgan Kaufmann.
[4] J. Besag. On the Statistical Analysis of Dirty Pictures. J. Royal Stat. Soc. B, 48:259–302, 1986.
[5] C. M. Bishop. Pattern Recognition and Machine Learning. Springer, New York, NY, USA, 2006.
[6] J.-F. Cai, E. J. Candes, and Z. Shen. A singular value thresholding algorithm for matrix completion. SIAM
Journal on Optimization, 20(4):1956–1982, 2008.
[7] O. Chapelle and Z. Harchaoui. A Machine Learning Approach to Conjoint Analysis. In Advances inneural information processing systems, volume 17, pages 257–264, 2005.
[8] D. R. Hardoon, S. R. Szedmak, and J. R. Shawe-Taylor. Canonical correlation analysis: An overview
with application to learning methods. Neural Computation, 16(12):2639–2664, 2004.
[9] M. Hazewinkel, editor. Encyclopaedia of Mathematics. Springer, 2002.
[10] W. James and C. Stein. Estimation with quadratic loss. In Proceedings of the 4th Berkeley Symposium on
Mathematical Statistics and Probability, volume 1, pages 361–379. University of California Press, 1961.
[11] H. Jeffreys. An Invariant Form for the Prior Probability in Estimation Problems. In Proceedings of theRoyal Society of London. Series A, volume 186, pages 453–461, 1946.
[12] S. Ji and J. Ye. An accelerated gradient method for trace norm minimization. In Proceedings of International Conference on Machine Learning, pages 457–464, 2009.
[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. Matlab Code for VBMF, http://sites.google.com/site/shinnkj23/, 2010.
[15] S. Nakajima and M. Sugiyama. Implicit regularization in variational Bayesian matrix factorization. In
Proceedings of 27th International Conference on Machine Learning (ICML2010), 2010.
[16] R. M. Neal. Bayesian Learning for Neural Networks. Springer, 1996.
[17] A. Paterek. Improving Regularized Singular Value Decomposition for Collaborative Filtering. In Proceedings of KDD Cup and Workshop, 2007.
[18] T. Raiko, A. Ilin, and J. Karhunen. Principal Component Analysis for Large Sale Problems with Lots ofMissing Values. In Proc. of ECML, volume 4701, pages 691–698, 2007.
[19] G. R. Reinsel and R. P. Velu. Multivariate reduced-rank Regression: Theory and Applications. Springer,
New York, 1998.
[20] J. D. M. Rennie and N. Srebro. Fast maximum margin matrix factorization for collaborative prediction.
In Proceedings of the 22nd International Conference on Machine learning, pages 713–719, 2005.
[21] R. Rosipal and N. Krämer. Overview and recent advances in partial least squares. In Subspace, Latent
Structure and Feature Selection Techniques, volume 3940, pages 34–51. Springer, 2006.
[22] 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, 2008.
[23] N. Srebro, J. Rennie, and T. Jaakkola. Maximum Margin Matrix Factorization. In Advances in NIPS,
volume 17, 2005.
[24] M.E. Tipping and C.M. Bishop. Probabilistic Principal Component Analysis. Journal of the Royal Statistical Society: Series B, 61(3):611–622, 1999.
[25] R. Tomioka, T. Suzuki, M. Sugiyama, and H. Kashima. An efﬁcient and general augmented Lagrangianalgorithm for learning low-rank matrices. In Proceedings of International Conference on Machine Learning, 2010.
[26] S. Watanabe. Algebraic Geometry and Statistical Learning. Cambridge University Press, Cambridge,
UK, 2009.
[27] K. J. Worsley, J-B. Poline, K. J. Friston, and A. C. Evanss. Characterizing the Response of PET and fMRI
Data Using Multivariate Linear Models. NeuroImage, 6(4):305–319, 1997.
[28] I-Cheng Yeh. Modeling slump ﬂow of concrete using second-order regressions and artiﬁcial neural networks. Cement and Concrete Composites, 29(6):474–480, 2007.
[29] K. Yu, V. Tresp, and A. Schwaighofer. Learning Gaussian Processes from Multiple Tasks. In Proc. of
ICML, page 1019, 2005.

9