Empirical Bayes Matrix Completion

arXiv:1706.01252v2 [stat.ML] 6 Jun 2017

Takeru MATSUDA1 and Fumiyasu KOMAKI1, 2
1 Department of Mathematical Informatics
Graduate School of Information Science and Technology
The University of Tokyo
2 RIKEN Brain Science Institute
{matsuda,komaki}@mist.i.u-tokyo.ac.jp

June 7, 2017

Abstract
We develop an empirical Bayes (EB) algorithm for the matrix completion problems. The EB algorithm is motivated from the singular value shrinkage estimator
for matrix means by Efron and Morris (1972). Since the EB algorithm is derived
as the EM algorithm applied to a simple model, it does not require heuristic parameter tuning other than tolerance. Numerical results demonstrated that the EB
algorithm achieves a good trade-off between accuracy and efficiency compared to
existing algorithms and that it works particularly well when the difference between
the number of rows and columns is large. Application to real data also shows the
practical utility of the EB algorithm.

1

Introduction

In various applications, we encounter problems of estimating the unobserved entries of
a matrix from the observed entries. For example, in the famous Netflix problem, we
have a matrix of movie ratings by users and aim to predict the preference for movies of
each user for recommendation. This problem is called the matrix completion problem
and many studies have investigated its theoretical properties (Candès and Recht, 2008;
Recht, 2011) and developed efficient algorithms (Srebro, 2005; Cai, Candès and Shen,
2010; Keshavan, Montanari, and Oh, 2010; Mazumder, Hastie and Tibshirani, 2010).
In the matrix completion problems, the low-rank property of the underlying matrix
plays a central role. For example, in the Netflix problem, the rank is interpreted as the
number of latent factors in the movie preference and it is believed to be small. Indeed,
existing matrix completion algorithms succeed in estimating the unobserved entries by
assuming the low-rankness. Note that low rank matrices have sparse singular values
since the rank of a matrix is equal to the number of its nonzero singular values. The
sum of singular values of a matrix is called the trace norm or nuclear norm, which is
employed by many existing algorithms for regularization.

1

In practice, the data matrix often contains observation noise and we aim to recover
the true underlying matrix. If the data matrix is fully observed with the Gaussian
observation noise, then the matrix completion problem reduces to the estimation of
the mean matrix parameter of a matrix-variate normal distribution. For this problem,
Efron and Morris (1972) developed an empirical Bayes estimator and proved that it is
minimax and dominates the maximum likelihood estimator under the Frobenius loss.
Later, Stein (1974) pointed out that this estimator shrinks the singular values of the
observed matrix for each. Therefore, this estimator performs well when the true value
of the mean matrix parameter has low rank. Based on this idea, Matsuda and Komaki
(2015) developed singular value shrinkage priors as a natural generalization of the Stein
prior. The singular value shrinkage priors are superharmonic and the Bayes estimators
based on them are minimax estimators with similar properties to the Efron–Morris
estimator.
In this study, we develop an empirical Bayes (EB) algorithm for matrix completion.
The EB algorithm is a natural extension of the Efron–Morris estimator. Since the EB
algorithm is essentially the EM algorithm applied to a simple model, it does not require
heuristic parameter tuning other than tolerance. Numerical experiments demonstrate
the effectiveness of the EB algorithm compared with existing algorithms. Specifically,
the EB algorithm works well when the difference between the number of rows and
columns is large. Application to real data also shows the practical utility of the EB
algorithm.
This paper is organized as follows. Section 2 reviews the previous results on the
empirical Bayes estimation of matrix means. Section 3 provides details of the EB
algorithm. Section 4 presents the results of the numerical experiments while Section 5
applies the EB algorithm to real data.

2

Empirical Bayes Estimation of Matrix Means

In this section, we review the empirical Bayes estimator by Efron and Morris (1972)
for the mean matrix parameter of a matrix-variate normal distribution. We extend
this estimator to matrix completion problems in the next section.
Suppose that we have a matrix observation Y = (y1 , · · · , yp )⊤ ∈ Rp×q , whose
row vectors yi⊤ ∈ Rq have the distribution yi ∼ Nq (mi , Iq ) independently, where Iq
q
is the q-dimensional identity matrix. Here, m⊤
i ∈ R is the i-th row vector of M =
(m1 , · · · , mp )⊤ ∈ Rp×q . In the notation of Dawid (1981), this situation is denoted
as Y ∼ Np,q (M, Ip , Iq ). We assume p ≥ q. We consider estimation of M under the
Frobenius loss:
L(M, M̂ (Y )) = kM̂ (Y ) − M k2F =

p X
q
X
(M̂ij (Y ) − Mij )2 .
i=1 j=1

Let Y = U ΛV ⊤ be the singular value decomposition of a matrix Y , where U ∈
O(p), V ∈ O(q), and Λ = [diag(σ1 (Y ), . . . , σq (Y )) Oq,p−q ]⊤ . Here, O(k) is the kdimensional orthogonal group, Olm is the zero matrix of size l × m, and σ1 (Y ) ≥
· · · ≥ σq (Y ) ≥ 0 are the singular values of Y . Similarly, let M̂ = Û Λ̂V̂ ⊤ be the
2

singular value decomposition of an estimator M̂ of M , where Û ∈ O(p), V̂ ∈ O(q),
i⊤
h
Λ̂ = diag(σ1 (M̂ ), · · · , σq (M̂ )) Oq,p−q , and σ1 (M̂ ) ≥ · · · ≥ σq (M̂ ) ≥ 0 are the
singular values of M̂ .
Efron and Morris (1972) proposed the estimator


M̂EM = Y Iq − (p − q − 1)(Y ⊤ Y )−1

(1)

and proved that it is minimax and dominates the maximum likelihood estimator M̂ =
Y when p − q − 1 > 0. When q = 1, the Efron–Morris estimator M̂EM coincides with
the James–Stein estimator. Stein (1974) pointed out that M̂EM can be represented in
the singular value decomposition form as follows:


p−q−1
σi (M̂EM ) = 1 −
σi (Y ) (i = 1, . . . , q),
σi (Y )2
Û = U,

V̂ = V.

Therefore, M̂EM shrinks the singular values of Y for each and preserves the singular
vectors of Y .
The Efron–Morris estimator M̂EM in (1) was derived as an empirical Bayes estimator based on the following hierarchical model:
M ∼ Np,q (0, Ip , Σ),
Y | M ∼ Np,q (M, Ip , Iq ).

q
⊤
The above model assumes that each row vector m⊤
i ∈ R of M = (m1 , · · · , mp ) ∈
Rp×q has the distribution mi ∼ Nq (0, Σ) independently. If Σ is given, then the Bayes
estimator of M is written as

M̂ Σ (Y ) = Y (Iq − (Iq + Σ)−1 ).

(2)

To obtain an empirical Bayes estimator, we estimate Σ from Y . Since the marginal
distribution of Y is Y ∼ Np,q (0, Ip , Iq +Σ), the marginal distribution of Y ⊤ Y is Y ⊤ Y ∼
Wq (p, Iq + Σ). From the property of the Wishart distribution, we have E(Y ⊤ Y )−1 =
(p−q−1)−1 (Iq +Σ)−1 . Thus, we can estimate Σ by (Iq + Σ̂)−1 = (p−q−1)(Y ⊤ Y )−1 . By
substituting this estimate into (2), the Efron–Morris estimator M̂EM in (1) is obtained.
We note that M̂EM is not a generalized Bayes estimator.
Recently, Matsuda and Komaki (2015) developed the singular value shrinkage prior
πSVS (M ) = det(M ⊤ M )−(p−q−1)/2
and proved its superharmonicity. When q = 1, the singular value shrinkage prior πSVS
coincides with the Stein prior πS (µ) = kµk2−p (Stein, 1974). The generalized Bayes
estimators based on the singular value shrinkage priors are minimax and have similar
properties to the Efron–Morris estimator M̂EM in (1). This is an extension of the

3

relationship between the James–Stein estimator and the Stein prior. We note that the
singular value shrinkage prior has the integral representation
Z
πSVS (M ) ∝ p(M | Ip , Σ)dΣ,
where p(M | Ip , Σ) is the probability density function of M ∼ Np,q (0, Ip , Σ) and dΣ is
the Lebesgue measure on the space of q × q positive-semidefinite matrices. This can
be confirmed by the calculation of the normalization constant in the inverse-Wishart
distribution.

3

The EB Algorithm

In this section, we propose the empirical Bayes (EB) algorithm for the matrix completion problems. This algorithm is motivated from the Efron–Morris estimator M̂EM in
(1).
We assume that the data matrix Y ∈ Rp×q has the distribution Y ∼ Np,q (M, σ 2 Ip , Iq ),
where p ≥ q and σ 2 is an unknown variance. Namely, each row vector yi⊤ ∈ Rq of
Y = (y1 , · · · , yp )⊤ ∈ Rp×q has the distribution yi ∼ Nq (mi , Iq ) independently, where
q
⊤
p×q . We observe only part
m⊤
i ∈ R is the i-th row vector of M = (m1 , · · · , mp ) ∈ R
of the entries of Y . Let Ω ⊂ {1, · · · , p}×{1, · · · , q} be the set of indices of the observed
entries and Ωi = {j | (i, j) ∈ Ω} (i = 1, · · · , p) be the set of indices of the observed
entries in the i-th row. We denote the observed entries of Y by YΩ . We denote the
submatrix of a matrix A with row indices R and column indices C as A[R; C]. For
example, if


1 2 3
A = 4 5 6 , R = {1, 2}, C = {1, 3},
7 8 9
then we have

A[R; C] =




1 3
.
4 6

Our goal is to estimate M from the observed entries YΩ of Y . We tackle this
problem using an empirical Bayes approach based on the following hierarchical model:
M ∼ Np,q (0, Ip , Σ),

(3)

Y | M ∼ Np,q (M, σ 2 Ip , Iq ).

(4)

q
⊤
p×q has
Namely, we assume that each row vector m⊤
i ∈ R of M = (m1 , · · · , mp ) ∈ R
the distribution mi ∼ Nq (0, Σ) independently. Note that the Efron–Morris estimator
M̂EM in (1) was also derived from this model with σ 2 = 1. Here, since only part of
the entries of Y are observed, we use the EM algorithm (Dempster, Laird and Rubin,
1972) to estimate the hyperparameters Σ and σ 2 . As a result, the EB algorithm is
described as Algorithm 1. Derivation of the EB algorithm is given in Appendix B.

4

Empirically, this algorithm converges in less than 20 iterations for most cases. We
note that the log-likelihood log p(YΩ | Σ, σ 2 ) is obtained as
p

|Ω|
1X
log p(YΩ | Σ, σ ) = −
log(2π) −
log det(σ 2 I|Ωi |+Σ[Ωi ,Ωi ] )
2
2
2

i=1

−

p
1X

2

Y [{i}; Ωi ](σ 2 I|Ωi |+Σ[Ωi ,Ωi ] )−1 Y [{i}; Ωi ]⊤ .

i=1

Similarly to the SOFT-IMPUTE algorithm by Mazumder, Hastie and Tibshirani (2010),the
EB algorithm can be viewed as iteratively imputing the missing entries of Y , although
we are updating not Y but M strictly speaking.

4

Numerical Experiments

In this section, we investigate the performance of the EB algorithm by numerical experiments. The EB algorithm is compared with the SVT algorithm by Cai, Candès and Shen
(2010), the SOFT-IMPUTE algorithm by Mazumder, Hastie and Tibshirani (2010),
and the OPTSPACE algorithm by Keshavan, Montanari, and Oh (2010). For these
existing algorithms, we use the MATLAB codes provided by the authors online.
The SVT algorithm (Cai, Candès and Shen, 2010) solves the following optimization
problem:
minimize kM̂ k∗
M̂

subject to

|Yij − M̂ij | ≤ Eij ,

(i, j) ∈ Ω,

where k · k∗ denotes the nuclear norm and Eij is a tolerance parameter. In the same
manner as the original paper, we set Eij equal to the standard deviation of the obser√
vation noise σ. We adopt the default settings of the algorithm parameters: τ = 5 pq,
δ = 1.2pq/|Ω|, kmax = 1000, ǫ = 10−4 .
The SOFT-IMPUTE algorithm (Mazumder, Hastie and Tibshirani, 2010) solves
the following optimization problem:
minimize kM̂ k∗
M̂
X
subject to
(Yij − M̂ij )2 ≤ δ,
(i,j)∈Ω

which is rewritten as
minimize
M̂

1 X
(Yij − M̂ij )2 + λkM̂ k∗ ,
2
(i,j)∈Ω

where λ ≥ 0 is a regularization parameter. In the same manner as the original paper,
we select λ from K candidate values by cross-validation using η % of the observed
entries as the validation set. We adopt the default settings of the algorithm parameters:
K = 20, η = 20, ǫ = 10−4 .
5

Algorithm 1 EB algorithm
Input: set of observation indices Ω, observed entries YΩ , initial value σ02 , and tolerance
ε1 , ε2
Output: M̂
Description: Estimate M from YΩ based on the model (3) and (4)

1:

Initialize Mold ∈ Rp×q by (Mold )ij =

2:

Initialize Σold = (Mold ⊤ Mold )/p

3:

Initialize (σ 2 )old = σ02

4:

while true do

5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:



Yij

0

Set Pi = Oq,q (q × q zero matrix)

−1
2 I
Set Pi [Ωi , Ωi ] = (σold
|Ωi | + Σold [Ωi , Ωi ])

Set Ri = Σold − Σold Pi
Σold

Yij
Set bi ∈ Rq by (bi )j =

0
end for

(j ∈ Ωi )
(j 6∈ Ωi )

2 )−1 (R b , · · · , R b )⊤
Set Mnew = (σold
1 1
p p

Pp
⊤
Set Σnew = Mnew Mnew + i=1 Ri /p

P
2
Set σnew
= (i,j)∈Ω (Yij − (Mnew )ij )2 + (Ri )jj /|Ω|

2 ) − log p(Y | Σ , σ 2 ) < ε then
if log p(YΩ | Σnew , σnew
Ω
1
old old

Set M̂ = Mnew and exit
end if

17:

if kMnew − Mold k2F /kMold k2F < ε2 then
Set M̂ = Mnew and exit

19:

end if

20:

Set Mold = Mnew

21:

Set Σold = Σnew

22:

2 = σ2
Set σold
new

23:

((i, j) 6∈ Ω)

for i = 1 · · · p do

16:

18:

((i, j) ∈ Ω)

end while

6

The OPTSPACE algorithm (Keshavan, Montanari, and Oh, 2010) achieves matrix
completion via spectral techniques and manifold optimization. Here, we use the option
of guessing the rank r from data. We adopt the default setting of algorithm parameters:
maximum number of iterations = 50, tolerance = 10−6 .
We consider the same experimental setting with Mazumder, Hastie and Tibshirani
(2010). We generate U ∈ Rp×r and V ∈ Rr×q whose entries are sampled from the
standard normal distribution N(0, 1) independently and put M = U V . Here, r denotes
the rank of M . Then, we generate Y = M +E, where the entries of E are sampled from
N(0, σ 2 ) independently. The indices of the observed entries Ω are randomly sampled
over pq indices of the matrix. We evaluate the accuracy of the matrix completion
algorithms by the normalized error for the overall matrix
error1 :=

kM̂ − M kF
,
kM kF

and the normalized error for the unobserved entries
P
( (i,j)6∈Ω (M̂ij − Mij )2 )1/2
P
error2 :=
.
( (i,j)6∈Ω Mij2 )1/2

We also compare the efficiency of the matrix completion algorithms by the computation
time in seconds.
For the EB algorithm, we set ε1 = 10−3 and ε2 = 10−4 . Also, we set σ02 = σ 2
except for Figure 6. In Figure 6, we investigate how the selection of σ02 affects the
performance of the EB algorithm.
Table 1 shows the results when p = 1000, q = 100, r = 10, σ 2 = 1, and
|Ω|/(pq) = 0.5. We present error1, error2, and the computation time averaged over
100 simulations. EB and OPTSPACE have the least error2 , whereas EB is faster than
OPTSPACE. SOFT-IMPUTE takes the least computation time, whereas its error2 is
almost twice as large as those of EB and OPTSPACE. Therefore, EB achieves a good
trade-off between accuracy and efficiency under this setting.
EB
SVT
SOFT-IMPUTE
OPTSPACE

error1
0.21
0.28
0.28
0.16

error2
0.18
0.31
0.31
0.17

time
4.63
5.12
2.75
8.06

Table 1: Performance of the matrix completion algorithms (p = 1000, q = 100, r = 10,
|Ω|/(pq) = 0.5, and σ 2 = 1)
Now, we investigate the dependence of the performance of each algorithm on the
size of the matrix (p, q), the true rank r, the proportion of observed entries |Ω|/(pq),
and the observation noise variance σ 2 . We present error2 and the computation time
averaged over 100 simulations for each setting. The qualitative behavior of error1 was
almost the same as that of error2.
Figure 1 plots error2 and the computation time as a function of p ∈ [102 , 105 ],
where q = 100, r = 10, σ 2 = 1, and |Ω|/(pq) = 0.5. Here, error2 decreases with
7

(a)

(b)

EB
SVT
SOFT-IMPUTE
OPTSPACE

-0.3

3
log10 time

log10 error2

-0.2

-0.4
-0.5

2

1
EB
SVT
SOFT-IMPUTE
OPTSPACE

-0.6
0
-0.7
-0.8
2

-1
3

4

5

2

log10 p

3

4

5

log10 p

Figure 1: Plot of (a) error2 and (b) computation time (in seconds) as a function of p,
where q = 100, r = 10, σ 2 = 1, and |Ω|/(pq) = 0.5.
p for EB, SVT, and SOFT-IMPUTE and the computation time increases with p for
all algorithms. EB has the best accuracy, whereas its computation time is almost the
same as those of SVT and SOFT-IMPUTE. The accuracy of OPTSPACE is almost the
same as EB when p ≤ 103 but becomes worse when p > 103 . Also, the computation
time of OPTSPACE grows with p faster than the other algorithms.
Figure 2 plots error2 and the computation time as a function of p ∈ [102 , 103 ],
where q = p, r = 10, σ 2 = 1, and |Ω|/(pq) = 0.5. Here, we are considering the case of
a square matrix. For all algorithms, error2 decreases with p and the computation time
increases with p. OPTSPACE has the best accuracy and efficiency. EB has almost the
same accuracy with SVT and SOFT-IMPUTE, whereas its computation time grows
with p a little faster than the others.
Figure 3 plots error2 and the computation time as a function of r ∈ [5, 50], where
p = 1000, q = 100, σ 2 = 1, and |Ω|/(pq) = 0.5. Here, error2 increases with r for
all algorithms and EB has the best accuracy for all values of r. On the other hand,
the behavior of the computation time varies among the four algorithms. Whereas the
computation time of EB and SOFT-IMPUTE increases with r, that of SVT decreases
with r and that of OPTSPACE is not monotone. As a whole, SVT and SOFT-IMPUTE
has a little better efficiency than EB and OPTSPACE.
Figure 4 plots error2 and the computation time1 as a function of |Ω|/(pq) ∈
[0.1, 0.9], where p = 1000, q = 100, r = 10, and σ 2 = 1. Here, error2 decreases
with |Ω|/(pq) for all algorithms. EB has the best accuracy for all values of |Ω|/(pq)
and OPTSPACE also attains the best accuracy when |Ω|/(pq) ≥ 0.4. The computation time of SOFT-IMPUTE is almost constant with |Ω|/(pq), whereas those of EB
1
In SVT, we set ǫ = 10−3 , since the default setting of ǫ = 10−4 often caused divergence. When
|Ω|/(pq) ≤ 0.3, SVT did not converge in kmax = 1000 iterations and so the computation time is
extremely large. In OPTSPACE, we omit the results for |Ω|/(pq) = 0.1, since the step of guessing the
rank sometimes did not finish.

8

(a)

(b)

3
EB
SVT
SOFT-IMPUTE
OPTSPACE

-0.4

2
log10 time

log10 error2

-0.2

-0.6
-0.8

1
EB
SVT
SOFT-IMPUTE
OPTSPACE

0
-1
-1
2

2.5
log10 p

3

2

2.5
log10 p

3

Figure 2: Plot of (a) error2 and (b) computation time (in seconds) as a function of p,
where q = p, r = 10, σ 2 = 1, and |Ω|/(pq) = 0.5.

(b)

(a)

0

log10 time

-0.2
log10 error2

EB
SVT
SOFT-IMPUTE
OPTSPACE

1.5

-0.4
-0.6
EB
SVT
SOFT-IMPUTE
OPTSPACE

-0.8

1

0.5

-1
10

20

30

40

50

10

r

20

30

40

50

r

Figure 3: Plot of (a) error2 and (b) computation time (in seconds) as a function of r,
where p = 1000, q = 100, σ 2 = 1, and |Ω|/(pq) = 0.5.

9

(a)

(b)

0.2
EB
SVT
SOFT-IMPUTE
OPTSPACE

-0.2

1
0.8
log10 time

log10 error2

0

-0.4
-0.6
-0.8

0.6
0.4
EB
SVT
SOFT-IMPUTE
OPTSPACE

0.2
0

-1

-0.2
0.2

0.4
0.6
|Ω|/(pq)

0.8

0.2

0.4
0.6
|Ω|/(pq)

0.8

Figure 4: Plot of (a) error2 and (b) computation time (in seconds) as a function of
|Ω|/(pq), where p = 1000, q = 100, r = 10, and σ 2 = 1.
and OPTSPACE increase with |Ω|/(pq) and almost converges at |Ω|/(pq) = 0.4. As a
whole, SOFT-IMPUTE has the best efficiency and EB is the second best.
Figure 5 plots error2 and the computation time as a function of σ 2 ∈ [10−1 , 102 ],
where p = 1000, q = 100, r = 10, and |Ω|/(pq) = 0.5. For all algorithms, error2
increases with σ 2 . Whereas OPTSPACE has the best accuracy when σ 2 ≤ 1, error2 is
almost the same for all algorithms when σ 2 ≥ 10. The computation time of EB and
SOFT-IMPUTE is almost constant with σ 2 , whereas that of SVT decreases with σ 2
and that of OPTSPACE decreases rapidly when σ 2 ≥ 10.
Finally, we investigate the effect of σ02 on the performance of EB. Figure 6 plots
error2 and the computation time of EB as a function of σ02 ∈ [10−2 , 102 ], where p =
1000, q = 100, r = 10, σ 2 = 1, and |Ω|/(pq) = 0.5. From Figure 6, we see that error2 is
almost constant with σ02 except for the case σ02 = 100. Therefore, the accuracy of EB
does not depend on σ02 provided that σ02 is not significantly larger than σ 2 . On the other
hand, the computation
time increases as σ02 becomes farther from σ 2 . Interestingly, the
√
2
selection σ0 = 10 provides better efficiency than the correct selection of σ02 = σ 2 = 1.
In summary, EB achieves a good trade-off between accuracy and efficiency and has
comparable performance to SVT, SOFT-IMPUTE and OPTSPACE. In particular, EB
works well when the difference between the number of rows and columns is large.

5

Application to Real Data

In this section, we apply the EB algorithm to the data set from the Jester online
joke recommender system (Goldberg et al., 2001). This data set comprises over 4.1
million continuous ratings (-10.00 to +10.00) of q = 100 jokes from p = 24983 users,
which were collected from April 1999 to May 2003 online. Similarly to the previous
section, we compare EB with SVT, SOFT-IMPUTE, and OPTSPACE. In EB, we set
ε1 = 10−3 , ε2 = 10−4 , and σ02 = 10. In SVT, SOFT-IMPUTE, and OPTSPACE, we
10

(a)

(b)

1.5

EB
SVT
SOFT-IMPUTE
OPTSPACE

log10 time

log10 error2

0

-0.5
EB
SVT
SOFT-IMPUTE
OPTSPACE

-1

-1

0

1

1

0.5

0
-1

2

0

log10 σ

1

2

log10 σ

2

2

Figure 5: Plot of (a) error2 and (b) computation time (in seconds) as a function of σ 2 ,
where p = 1000, q = 100, r = 10, and |Ω|/(pq) = 0.5.

(b)

(a)

1
0.9

-0.3

log10 time

log10 error2

-0.2

-0.4
-0.5
-0.6

0.7
0.6
0.5

-0.7
-0.8
-2

0.8

-1

0

1

0.4
-2

2

log10 σ

-1

0

1

2

log10 σ

2
0

2
0

Figure 6: Plot of (a) error2 and (b) computation time (in seconds) as a function of σ02 ,
where p = 1000, q = 100, r = 10, |Ω|/(pq) = 0.5 and σ 2 = 1.

11

adopt the default settings of the algorithm parameters.
In the original data set, |Ω0 | = 1810455 among pq = 2498300 entries were observed.
We randomly sampled |Ω| = 5 × 105 entries from the observed entries and applied
matrix completion algorithms to the sampled entries. We evaluated the accuracy of
each algorithm by the normalized error for the rest of the |Ω0 | − |Ω| = 1310455 entries:
P
( (i,j)∈Ω0 \Ω (M̂ij − Mij )2 )1/2
P
.
error :=
( (i,j)∈Ω0 \Ω Mij2 )1/2

Table 2 shows the results. We only present results for EB, SOFT-IMPUTE, and
OPTSPACE since SVT diverged2 . Whereas the accuracies of the three algorithms are
almost identical, EB takes the least computation time among the three algorithms.
This result shows the practical utility of EB on real data.
EB
SOFT-IMPUTE
OPTSPACE

error
0.85
0.82
0.83

time
52.21
100.11
68.71

Table 2: Performance of matrix completion algorithms on the Jester data set

6

Conclusion

In this study, we proposed an empirical Bayes (EB) algorithm for matrix completion.
The EB algorithm is motivated from the Efron–Morris estimator for a normal mean
matrix. It is free from heuristic parameter tuning other than tolerance. Numerical
results demonstrated that the EB algorithm achieves a good trade-off between accuracy
and efficiency compared to existing algorithms and that it works particularly well
when the difference between the number of rows and columns is large. In addition,
application to real data showed the practical utility of the EB algorithm.

Acknowledgements
We thank Edward I. George for helpful comments.

A

Posterior of Normal Mean from Missing Observation

Suppose that m ∼ Nq (0, Σ) and X | m ∼ Nq (m, σ 2 Iq ). Here, we derive the posterior
distribution of m given a subvector of X.
First, consider the case where we observe the first k entries of X, which we denote
by Z = (X1 , · · · , Xk ). Let C ∈ Rk×k be defined by
C = (σ 2 Ik + Σ[γk , γk ])−1 ,
2

γk = {1, · · · , k}.

We also tried Eij = σ̂ instead of Eij = σ0 , where σ̂ 2 = 18.21 is the estimate of σ 2 by EB, but it
diverged again. Changing ǫ from 10−4 to 10−3 was not useful either.

12

Then, the posterior distribution of m given Z = z is obtained as
p(m | Z = z) ∝ π(m)p(Z = z | m)

!
Pk
2
1 ⊤ −1
i=1 (zi − mi )
∝ exp − m Σ m −
2
2σ 2


1
⊤ −1
∝ exp − (m − m̂) R (m − m̂) ,
2

where
R=



Σ

−1

+σ

−2






C O
=Σ−Σ
Σ,
O O
 
z
−2
.
m̂ = σ R
0
Ik O
O O

−1

Here, we used the Sherman-Morrison-Woodbury formula (Golub and van Loan, 1996)
(Σ−1 + U ⊤ U )−1 = Σ − ΣU ⊤ (I + U ΣU ⊤ )−1 U Σ.
Therefore, the posterior distribution of m given Z = z is the normal distribution
N(m̂, R).
The above result is straightforwardly extended to the general case where we observe
k entries of X with indices ω = {ω1 , · · · , ωk } ⊂ {1, · · · , q}. We denote the observed
subvector of X by Z = (Xω1 , · · · , Xωk ). Let P ∈ Rq×q and b ∈ Rq be defined by
P [ω, ω] = (σ 2 Ik + Σ[ω, ω])−1 , Pj1 j2 = 0 (j1 6∈ ω or j2 6∈ ω),
(
Zl (j = ωl )
.
bj =
0 (j ∈
6 ω)
Then, the posterior distribution of m given Z = z is obtained as


1
⊤ −1
p(m | Z = z) ∝ exp − (m − m̂) R (m − m̂) ,
2
where
R = Σ − ΣP Σ,
m̂ = σ −2 Rb.

Therefore, the posterior distribution of m given Z = z is the normal distribution
N(m̂, R).

B

Derivation of the EB Algorithm

In general, the EM algorithm is used to iteratively estimate the parameters of statistical
models with latent variables p(X, Z | θ) (Dempster, Laird and Rubin, 1972). Here, X

13

denotes the observed variables and Z denotes the latent variables. Each iteration in
the EM algorithm is described as
θt+1 = arg max Q(θ, θt ),
θ

where
Q(θ, θt ) = EZ|X,θt [log p(X, Z | θ)].

In the present model (3) and (4), θ = (Σ, σ 2 ) is the parameter, YΩ is the observed
variables, and M is the latent variables. The log-likelihood function is calculated as
log p(YΩ , M | θ) = log p(M | Σ) + log p(YΩ | M, σ 2 )
p
1
|Ω|
1 X
= − log det Σ − trΣ−1 M M ⊤ −
(Yij − Mij )2 .
log σ 2 − 2
2
2
2
2σ
(i,j)∈Ω

From the covariance structure in (3) and (4), row vectors m1 , · · · , mp of M = (m1 , · · · , mp )⊤
are independent given YΩ . Using the results in Appendix A, the posterior distribution
of each mi given YΩ is obtained as
mi | YΩ , θ ∼ N(m̂i , Ri ),
where
Ri = Σ − ΣPi Σ,
m̂i = σ −2 Ri bi .

Here, Pi ∈ Rq×q and bi ∈ Rq are defined by

Pi [Ωi , Ωi ] = (σ 2 I|Ωi | + Σ[Ωi , Ωi ])−1 , (Pi )j1 j2 = 0 (j1 6∈ Ωi or j2 6∈ Ωi ),
(
Yij (j ∈ Ωi )
.
(bi )j =
0
(j 6∈ Ωi )

Therefore,

where

e θ) = EM |Y ,θ [log p(YΩ , M | θ)]
e
Q(θ,
Ω
e Σ) + Q2 (e
= Q1 (Σ,
σ 2 , σ 2 ),
p


X 
e Σ) = − p log det Σ
e−1
e −1 m̂i m̂⊤ + Ri ,
Q1 (Σ,
tr Σ
i
2
2
i=1

|Ω|
1 X
Q2 (e
σ2 , σ2 ) = −
log σ
e2 − 2
((Yij − (m̂i )j )2 + (Ri )jj ).
2
2e
σ
(i,j)∈Ω

e Σ) and Q2 (e
e and σ
By maximizing Q1 (Σ,
σ 2 , σ 2 ) with respect to Σ
e2 respectively, we
obtain
p 

X
e=1
m̂i m̂⊤
+
R
Σ
i ,
i
p
i=1
X

1
(Yij − (m̂i )j )2 + (Ri )jj .
σ
e2 =
|Ω|
(i,j)∈Ω

Thus, Algorithm 1 is obtained.

14

References
Cai, J. F., Candès, E. J. & Shen, Z. (2010). A singular value thresholding algorithm
for matrix completion. SIAM Journal on Optimization 20, 1956–1982.
Candès, E. J. & Recht, B. (2008). Exact matrix completion via convex optimization.
Foundations of Computational Mathematics 9, 717–772.
Dawid, A. P. (1981). Some matrix-variate distribution theory: notational considerations and a Bayesian application. Biometrika 68, 265–274.
Dempster, A. P., Laird, N. M. & Rubin, D. B. (1972). Maximum likelihood from
incomplete data via the EM algorithm. Journal of the Royal Statistical Society B
39, 1–38.
Efron, B. & Morris, C. (1972). Empirical Bayes on vector observations: an extension of Stein’s method. Biometrika 59, 335–347.
Goldberg, K., Roeder, T., Gupta, D. & Perkins, C. (2001). Eigentaste: A
constant time collaborative filtering algorithm. Information Retrieval 4, 133–151.
Golub, G. H. & van Loan, C. F. (1996). Matrix Computations. Baltimore, MD:
Johns Hopkins.
Keshavan, R. H., Montanari, A. & Oh, S. (2010). Matrix completion from noisy
entires. Journal of Machine Learning Research 11, 2057–2078.
Matsuda, T. & Komaki, F. (2015). Singular value shrinkage priors for Bayesian
prediction. Biometrika 102, 843–854.
Mazumder, R., Hastie, T. & Tibshirani, R. (2010). Spectral regularization algorithms for learning large incomplete matrices. Journal of Machine Learning Research
11, 2287–2322.
Recht, B. (2011). A simpler approach to matrix completion. Journal of Machine
Learning Research 12, 3413–3430.
Srebro, N., Rennie, J. & Jaakkola, T. (2005). Maximum-margin matrix factorization. In Advances in Neural Information Processing Systems 17, 1329–1336.
Stein, C. (1974). Estimation of the mean of a multivariate normal distribution.
Proceedings of Prague Symposium on Asymptotic Statistics 2, 345–381.

15

