Empirical Bayes Matrix Completion

TL;DR

This paper introduces an empirical Bayes algorithm for matrix completion that leverages singular value shrinkage, offering a parameter-tuning free, efficient, and accurate solution especially effective for matrices with large row-column disparities.

Contribution

The paper presents a novel EB algorithm for matrix completion based on singular value shrinkage, eliminating heuristic tuning and improving performance over existing methods.

Findings

Achieves a good balance between accuracy and efficiency.

Performs well when the number of rows and columns differ greatly.

Demonstrates practical utility on real datasets.

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 essentially 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.

Tables

Table 1. Table 1: Performance of the matrix completion algorithms ( p = 1000 𝑝 1000 p=1000 , q = 100 𝑞 100 q=100 , r = 10 𝑟 10 r=10 , | Ω | / ( p ​ q ) = 0.5 Ω 𝑝 𝑞 0.5 |\Omega|/(pq)=0.5 , and σ 2 = 1 superscript 𝜎 2 1 \sigma^{2}=1 )

	error1	error2	time
EB	0.21	0.18	4.63
SVT	0.28	0.31	5.12
SOFT-IMPUTE	0.28	0.31	2.75
OPTSPACE	0.16	0.17	8.06

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

	error	time
EB	0.85	52.21
SOFT-IMPUTE	0.82	100.11
OPTSPACE	0.83	68.71

Equations

$$L(M,\hat{M}(Y))=\Vert \hat{M}(Y)-M{\Vert }_{\mathrm{F}}^2=\sum _{i=1}^p\sum _{j=1}^q({\hat{M}}_{ij}(Y)-M_{ij}{)}^2.$$

$${\hat{M}}_{\mathrm{EM}}=Y\left(I_q-(p-q-1)(Y^{\top }Y{)}^{-1}\right)$$

$${\sigma }_i({\hat{M}}_{\mathrm{EM}})=\left(1-\frac{p-q-1}{{\sigma }_i(Y{)}^2}\right){\sigma }_i(Y)(i=1,\dots ,q),$$

$$\hat{U}=U,\hat{V}=V.$$

$$M\sim {\mathrm{N}}_{p,q}(0,I_p,\Sigma ),$$

$$Y\mid M\sim {\mathrm{N}}_{p,q}(M,I_p,I_q).$$

$${\hat{M}}^{\Sigma }(Y)=Y(I_q-(I_q+\Sigma {)}^{-1}).$$

$${\pi }_{\mathrm{SVS}}(M)=\det (M^{\top }M{)}^{-(p-q-1)/2}$$

$${\pi }_{\mathrm{SVS}}(M)\propto \int p(M\mid I_p,\Sigma )\mathrm{d}\Sigma ,$$

$$A=\left(\begin{array}{ccc}1& 2& 3\\ 4& 5& 6\\ 7& 8& 9\end{array}\right),R=\{1,2\},C=\{1,3\},$$

$$A[R;C]=\left(\begin{array}{cc}1& 3\\ 4& 6\end{array}\right).$$

$$M\sim {\mathrm{N}}_{p,q}(0,I_p,\Sigma ),$$

$$Y\mid M\sim {\mathrm{N}}_{p,q}(M,{\sigma }^2{I}_p,I_q).$$

$$\log p(Y_{\Omega }\mid \Sigma ,{\sigma }^2)=$$

$$-\frac{1}{2}\sum _{i=1}^{p}Y[\{i\};{\Omega }_i]({\sigma }^2{I}_{|{\Omega }_i|+\Sigma [{\Omega }_i,{\Omega }_i]}{)}^{-1}Y[\{i\};{\Omega }_i{]}^{\top }.$$

$$\underset{\hat{M}}{\mathrm{minimize}}$$

$$\mathrm{subject}\mathrm{to}$$

$$\underset{\hat{M}}{\mathrm{minimize}}$$

$$\mathrm{subject}\mathrm{to}$$

$$\underset{\hat{M}}{\mathrm{minimize}}\frac{1}{2}\sum _{(i,j)\in \Omega }(Y_{ij}-{\hat{M}}_{ij}{)}^2+\lambda \Vert \hat{M}{\Vert }_{\ast },$$

$$error1:=\frac{\Vert \hat{M}-M{\Vert }_{\mathrm{F}}}{\Vert M{\Vert }_{\mathrm{F}}},$$

$$error2:=\frac{({\sum }_{(i,j)\notin \Omega }({\hat{M}}_{ij}-M_{ij}{)}^2{)}^{1/2}}{({\sum }_{(i,j)\notin \Omega }{M}_{ij}^2{)}^{1/2}}.$$

$$\mathrm{error}:=\frac{({\sum }_{(i,j)\in {\Omega }_0\setminus \Omega }({\hat{M}}_{ij}-M_{ij}{)}^2{)}^{1/2}}{({\sum }_{(i,j)\in {\Omega }_0\setminus \Omega }{M}_{ij}^2{)}^{1/2}}.$$

$$C=({\sigma }^2{I}_k+\Sigma [{\gamma }_k,{\gamma }_k]{)}^{-1},{\gamma }_k=\{1,\cdots ,k\}.$$

$$p(m\mid Z=z)$$

$$\propto \exp \left(-\frac{1}{2}{m}^{\top }{\Sigma }^{-1}m-\frac{{\sum }_{i=1}^k(z_i-m_i{)}^2}{2{\sigma }^2}\right)$$

$$\propto \exp \left(-\frac{1}{2}(m-\hat{m}{)}^{\top }{R}^{-1}(m-\hat{m})\right),$$

$$R={\left({\Sigma }^{-1}+{\sigma }^{-2}\left(\begin{array}{cc}{I}_k& O\\ O& O\end{array}\right)\right)}^{-1}=\Sigma -\Sigma \left(\begin{array}{cc}C& O\\ O& O\end{array}\right)\Sigma ,$$

$$\hat{m}={\sigma }^{-2}R\left(\begin{array}{c}z\\ 0\end{array}\right).$$

$$({\Sigma }^{-1}+U^{\top }U{)}^{-1}=\Sigma -\Sigma U^{\top }(I+U\Sigma U^{\top }{)}^{-1}U\Sigma .$$

$$P[\omega ,\omega ]=({\sigma }^2{I}_k+\Sigma [\omega ,\omega ]{)}^{-1},P_{j_1{j}_2}=0(j_1\notin \omega \mathrm{or}{j}_2\notin \omega ),$$

$$b_j=\{\begin{array}{ll}{Z}_l& (j={\omega }_l)\\ 0& (j\notin \omega )\end{array}.$$

$$p(m\mid Z=z)$$

Full text

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 (Candès and Recht, 2008; Recht, 2011) and developed efficient algorithms (Srebro, 2005; Cai, Candè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.

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=(y_{1},\cdots,y_{p})^{\top}\in\mathbb{R}^{p\times q}$, whose row vectors $y_{i}^{\top}\in\mathbb{R}^{q}$ have the distribution $y_{i}\sim{\rm N}_{q}(m_{i},I_{q})$ independently, where $I_{q}$ is the $q$-dimensional identity matrix. Here, $m_{i}^{\top}\in\mathbb{R}^{q}$ is the $i$-th row vector of $M=(m_{1},\cdots,m_{p})^{\top}\in\mathbb{R}^{p\times q}$. In the notation of Dawid (1981), this situation is denoted as $Y\sim{\rm N}_{p,q}(M,I_{p},I_{q})$. We assume $p\geq q$. We consider estimation of $M$ under the Frobenius loss:

[TABLE]

Let $Y=U\Lambda V^{\top}$ be the singular value decomposition of a matrix $Y$, where $U\in O(p)$, $V\in O(q)$, and $\Lambda=\left[{\rm diag}(\sigma_{1}(Y),\ldots,\sigma_{q}(Y))\;O_{q,p-q}\right]^{\top}$. Here, $O(k)$ is the $k$-dimensional orthogonal group, $O_{lm}$ is the zero matrix of size $l\times m$, and $\sigma_{1}(Y)\geq\cdots\geq\sigma_{q}(Y)\geq 0$ are the singular values of $Y$. Similarly, let $\hat{M}=\hat{U}\hat{\Lambda}\hat{V}^{\top}$ be the singular value decomposition of an estimator $\hat{M}$ of $M$, where $\hat{U}\in O(p)$, $\hat{V}\in O(q)$, $\hat{\Lambda}=\left[{\rm diag}(\sigma_{1}(\hat{M}),\cdots,\sigma_{q}(\hat{M}))\;O_{q,p-q}\right]^{\top}$, and $\sigma_{1}(\hat{M})\geq\cdots\geq\sigma_{q}(\hat{M})\geq 0$ are the singular values of $\hat{M}$.

Efron and Morris (1972) proposed the estimator

[TABLE]

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

[TABLE]

[TABLE]

Therefore, $\hat{M}_{{\rm EM}}$ shrinks the singular values of $Y$ for each and preserves the singular vectors of $Y$.

The Efron–Morris estimator $\hat{M}_{{\rm EM}}$ in (1) was derived as an empirical Bayes estimator based on the following hierarchical model:

[TABLE]

[TABLE]

The above model assumes that each row vector $m_{i}^{\top}\in\mathbb{R}^{q}$ of $M=(m_{1},\cdots,m_{p})^{\top}\in\mathbb{R}^{p\times q}$ has the distribution $m_{i}\sim{\rm N}_{q}(0,\Sigma)$ independently. If $\Sigma$ is given, then the Bayes estimator of $M$ is written as

[TABLE]

To obtain an empirical Bayes estimator, we estimate $\Sigma$ from $Y$. Since the marginal distribution of $Y$ is $Y\sim{\rm N}_{p,q}(0,I_{p},I_{q}+\Sigma)$, the marginal distribution of $Y^{\top}Y$ is $Y^{\top}Y\sim W_{q}(p,I_{q}+\Sigma)$. From the property of the Wishart distribution, we have ${\rm E}(Y^{\top}Y)^{-1}=(p-q-1)^{-1}(I_{q}+\Sigma)^{-1}$. Thus, we can estimate $\Sigma$ by $(I_{q}+\hat{\Sigma})^{-1}=(p-q-1)(Y^{\top}Y)^{-1}$. By substituting this estimate into (2), the Efron–Morris estimator $\hat{M}_{{\rm EM}}$ in (1) is obtained. We note that $\hat{M}_{{\rm EM}}$ is not a generalized Bayes estimator.

Recently, Matsuda and Komaki (2015) developed the singular value shrinkage prior

[TABLE]

and proved its superharmonicity. When $q=1$, the singular value shrinkage prior $\pi_{{\rm SVS}}$ coincides with the Stein prior $\pi_{{\rm S}}(\mu)=\|\mu\|^{2-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 $\hat{M}_{{\rm EM}}$ in (1). This is an extension of the relationship between the James–Stein estimator and the Stein prior. We note that the singular value shrinkage prior has the integral representation

[TABLE]

where $p(M\mid I_{p},\Sigma)$ is the probability density function of $M\sim{\rm N}_{p,q}(0,I_{p},\Sigma)$ and ${\rm d}\Sigma$ is the Lebesgue measure on the space of $q\times 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 $\hat{M}_{{\rm EM}}$ in (1).

We assume that the data matrix $Y\in\mathbb{R}^{p\times q}$ has the distribution $Y\sim{\rm N}_{p,q}(M,\sigma^{2}I_{p},I_{q})$, where $p\geq q$ and $\sigma^{2}$ is an unknown variance. Namely, each row vector $y_{i}^{\top}\in\mathbb{R}^{q}$ of $Y=(y_{1},\cdots,y_{p})^{\top}\in\mathbb{R}^{p\times q}$ has the distribution $y_{i}\sim{\rm N}_{q}(m_{i},I_{q})$ independently, where $m_{i}^{\top}\in\mathbb{R}^{q}$ is the $i$-th row vector of $M=(m_{1},\cdots,m_{p})^{\top}\in\mathbb{R}^{p\times q}$. We observe only part of the entries of $Y$. Let $\Omega\subset\{1,\cdots,p\}\times\{1,\cdots,q\}$ be the set of indices of the observed entries and $\Omega_{i}=\{j\mid(i,j)\in\Omega\}\ (i=1,\cdots,p)$ be the set of indices of the observed entries in the $i$-th row. We denote the observed entries of $Y$ by $Y_{\Omega}$. We denote the submatrix of a matrix $A$ with row indices $R$ and column indices $C$ as $A[R;C]$. For example, if

[TABLE]

then we have

[TABLE]

Our goal is to estimate $M$ from the observed entries $Y_{\Omega}$ of $Y$. We tackle this problem using an empirical Bayes approach based on the following hierarchical model:

[TABLE]

[TABLE]

Namely, we assume that each row vector $m_{i}^{\top}\in\mathbb{R}^{q}$ of $M=(m_{1},\cdots,m_{p})^{\top}\in\mathbb{R}^{p\times q}$ has the distribution $m_{i}\sim{\rm N}_{q}(0,\Sigma)$ independently. Note that the Efron–Morris estimator $\hat{M}_{{\rm EM}}$ in (1) was also derived from this model with $\sigma^{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 $\Sigma$ and $\sigma^{2}$. As a result, the EB algorithm is described as Algorithm 1. Derivation of the EB algorithm is given in Appendix B. Empirically, this algorithm converges in less than 20 iterations for most cases. We note that the log-likelihood $\log p(Y_{\Omega}\mid\Sigma,\sigma^{2})$ is obtained as

[TABLE]

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, Candè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, Candès and Shen, 2010) solves the following optimization problem:

[TABLE]

where $\|\cdot\|_{*}$ denotes the nuclear norm and $E_{ij}$ is a tolerance parameter. In the same manner as the original paper, we set $E_{ij}$ equal to the standard deviation of the observation noise $\sigma$. We adopt the default settings of the algorithm parameters: $\tau=5\sqrt{pq}$, $\delta=1.2pq/|\Omega|$, $k_{\max}=1000$, $\epsilon=10^{-4}$.

The SOFT-IMPUTE algorithm (Mazumder, Hastie and Tibshirani, 2010) solves the following optimization problem:

[TABLE]

which is rewritten as

[TABLE]

where $\lambda\geq 0$ is a regularization parameter. In the same manner as the original paper, we select $\lambda$ from $K$ candidate values by cross-validation using $\eta$ % of the observed entries as the validation set. We adopt the default settings of the algorithm parameters: $K=20$, $\eta=20$, $\epsilon=10^{-4}$.

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\in\mathbb{R}^{p\times r}$ and $V\in\mathbb{R}^{r\times q}$ whose entries are sampled from the standard normal distribution ${\rm N}(0,1)$ independently and put $M=UV$. Here, $r$ denotes the rank of $M$. Then, we generate $Y=M+E$, where the entries of $E$ are sampled from ${\rm N}(0,\sigma^{2})$ independently. The indices of the observed entries $\Omega$ 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

[TABLE]

and the normalized error for the unobserved entries

[TABLE]

We also compare the efficiency of the matrix completion algorithms by the computation time in seconds.

For the EB algorithm, we set $\varepsilon_{1}=10^{-3}$ and $\varepsilon_{2}=10^{-4}$. Also, we set $\sigma_{0}^{2}=\sigma^{2}$ except for Figure 6. In Figure 6, we investigate how the selection of $\sigma_{0}^{2}$ affects the performance of the EB algorithm.

Table 1 shows the results when $p=1000$, $q=100$, $r=10$, $\sigma^{2}=1$, and $|\Omega|/(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.

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 $|\Omega|/(pq)$, and the observation noise variance $\sigma^{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\in[10^{2},10^{5}]$, where $q=100$, $r=10$, $\sigma^{2}=1$, and $|\Omega|/(pq)=0.5$. Here, error2 decreases with $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\leq 10^{3}$ but becomes worse when $p>10^{3}$. 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\in[10^{2},10^{3}]$, where $q=p$, $r=10$, $\sigma^{2}=1$, and $|\Omega|/(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\in[5,50]$, where $p=1000$, $q=100$, $\sigma^{2}=1$, and $|\Omega|/(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 time111In SVT, we set $\epsilon=10^{-3}$, since the default setting of $\epsilon=10^{-4}$ often caused divergence. When $|\Omega|/(pq)\leq 0.3$, SVT did not converge in $k_{\max}=1000$ iterations and so the computation time is extremely large. In OPTSPACE, we omit the results for $|\Omega|/(pq)=0.1$, since the step of guessing the rank sometimes did not finish. as a function of $|\Omega|/(pq)\in[0.1,0.9]$, where $p=1000$, $q=100$, $r=10$, and $\sigma^{2}=1$. Here, error2 decreases with $|\Omega|/(pq)$ for all algorithms. EB has the best accuracy for all values of $|\Omega|/(pq)$ and OPTSPACE also attains the best accuracy when $|\Omega|/(pq)\geq 0.4$. The computation time of SOFT-IMPUTE is almost constant with $|\Omega|/(pq)$, whereas those of EB and OPTSPACE increase with $|\Omega|/(pq)$ and almost converges at $|\Omega|/(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 $\sigma^{2}\in[10^{-1},10^{2}]$, where $p=1000$, $q=100$, $r=10$, and $|\Omega|/(pq)=0.5$. For all algorithms, error2 increases with $\sigma^{2}$. Whereas OPTSPACE has the best accuracy when $\sigma^{2}\leq 1$, error2 is almost the same for all algorithms when $\sigma^{2}\geq 10$. The computation time of EB and SOFT-IMPUTE is almost constant with $\sigma^{2}$, whereas that of SVT decreases with $\sigma^{2}$ and that of OPTSPACE decreases rapidly when $\sigma^{2}\geq 10$.

Finally, we investigate the effect of $\sigma_{0}^{2}$ on the performance of EB. Figure 6 plots error2 and the computation time of EB as a function of $\sigma_{0}^{2}\in[10^{-2},10^{2}]$, where $p=1000$, $q=100$, $r=10$, $\sigma^{2}=1$, and $|\Omega|/(pq)=0.5$. From Figure 6, we see that error2 is almost constant with $\sigma_{0}^{2}$ except for the case $\sigma_{0}^{2}=100$. Therefore, the accuracy of EB does not depend on $\sigma_{0}^{2}$ provided that $\sigma_{0}^{2}$ is not significantly larger than $\sigma^{2}$. On the other hand, the computation time increases as $\sigma_{0}^{2}$ becomes farther from $\sigma^{2}$. Interestingly, the selection $\sigma_{0}^{2}=\sqrt{10}$ provides better efficiency than the correct selection of $\sigma_{0}^{2}=\sigma^{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 $\varepsilon_{1}=10^{-3}$, $\varepsilon_{2}=10^{-4}$, and $\sigma_{0}^{2}=10$. In SVT, SOFT-IMPUTE, and OPTSPACE, we adopt the default settings of the algorithm parameters.

In the original data set, $|\Omega_{0}|=1810455$ among $pq=2498300$ entries were observed. We randomly sampled $|\Omega|=5\times 10^{5}$ 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 $|\Omega_{0}|-|\Omega|=1310455$ entries:

[TABLE]

Table 2 shows the results. We only present results for EB, SOFT-IMPUTE, and OPTSPACE since SVT diverged222We also tried $E_{ij}=\hat{\sigma}$ instead of $E_{ij}=\sigma_{0}$, where $\hat{\sigma}^{2}=18.21$ is the estimate of $\sigma^{2}$ by EB, but it diverged again. Changing $\epsilon$ from $10^{-4}$ to $10^{-3}$ was not useful either.. 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.

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.

Appendix A Posterior of Normal Mean from Missing Observation

Suppose that $m\sim{\rm N}_{q}(0,\Sigma)$ and $X\mid m\sim{\rm N}_{q}(m,\sigma^{2}I_{q})$. 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=(X_{1},\cdots,X_{k})$. Let $C\in\mathbb{R}^{k\times k}$ be defined by

[TABLE]

Then, the posterior distribution of $m$ given $Z=z$ is obtained as

[TABLE]

where

[TABLE]

[TABLE]

Here, we used the Sherman-Morrison-Woodbury formula (Golub and van Loan, 1996)

[TABLE]

Therefore, the posterior distribution of $m$ given $Z=z$ is the normal distribution ${\rm N}(\hat{m},R)$.

The above result is straightforwardly extended to the general case where we observe $k$ entries of $X$ with indices $\omega=\{\omega_{1},\cdots,\omega_{k}\}\subset\{1,\cdots,q\}$. We denote the observed subvector of $X$ by $Z=(X_{\omega_{1}},\cdots,X_{\omega_{k}})$. Let ${P}\in\mathbb{R}^{q\times q}$ and $b\in\mathbb{R}^{q}$ be defined by

[TABLE]

[TABLE]

Then, the posterior distribution of $m$ given $Z=z$ is obtained as

[TABLE]

where

[TABLE]

[TABLE]

Therefore, the posterior distribution of $m$ given $Z=z$ is the normal distribution ${\rm N}(\hat{m},R)$.

Appendix 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\mid\theta)$ (Dempster, Laird and Rubin, 1972). Here, $X$ denotes the observed variables and $Z$ denotes the latent variables. Each iteration in the EM algorithm is described as

[TABLE]

where

[TABLE]

In the present model (3) and (4), $\theta=(\Sigma,\sigma^{2})$ is the parameter, $Y_{\Omega}$ is the observed variables, and $M$ is the latent variables. The log-likelihood function is calculated as

[TABLE]

From the covariance structure in (3) and (4), row vectors $m_{1},\cdots,m_{p}$ of $M=(m_{1},\cdots,m_{p})^{\top}$ are independent given $Y_{\Omega}$. Using the results in Appendix A, the posterior distribution of each $m_{i}$ given $Y_{\Omega}$ is obtained as

[TABLE]

where

[TABLE]

[TABLE]

Here, ${P}_{i}\in\mathbb{R}^{q\times q}$ and $b_{i}\in\mathbb{R}^{q}$ are defined by

[TABLE]

[TABLE]

Therefore,

[TABLE]

where

[TABLE]

[TABLE]

By maximizing $Q_{1}(\widetilde{\Sigma},\Sigma)$ and $Q_{2}(\widetilde{\sigma}^{2},\sigma^{2})$ with respect to $\widetilde{\Sigma}$ and $\widetilde{\sigma}^{2}$ respectively, we obtain

[TABLE]

[TABLE]

Thus, Algorithm 1 is obtained.