Sparse Probabilistic Matrix Factorization by Laplace Distribution for Collaborative Filtering

1 Introduction

With the emerging of big data, recommender systems play a more and more important role to provide personalized recommendation and improve the user experience [Adomavicius and Tuzhilin, 2005]. In the last decade, the collaborative filtering (CF) methods based on matrix factorization (MF) [Koren et al., 2009] have shown their ability to build accurate prediction models and are widely adopted in commercial world such as Amazon, Google and Netflix [Dror et al., 2012]. MF methods predict the user preference by determining the latent features of users and items. More precisely, MF method factorizes the rating matrix into two low-rank matrices, one for latent user factor and the other for latent item factor. By multiplying these two latent factors, we can complete the original incomplete rating matrix and do recommendation.

Even though a variety of MF methods [Berry et al., 1999; Salakhutdinov and Mnih, 2007; 2008; Chen et al., 2011; Shi et al., 2013; Bauer and Nanopoulos, 2014] have been proposed and lead to promising recommendation results, they are limited in modern recommendation systems due to the sparse and large-scale rating data. The sparsity of data refers to the lack of observed rating data that makes it difficult and unreliable to predict the user preference. One accompanying phenomenon is long-tail effect where a large fraction of items have few ratings [Anderson, 2006]. Although the amount of users relating to each individual tail item is small in absolute numbers, collectively they cover a substantial fraction of all users. Additionally, a user's rarer purchases in e-commerce are also more informative of their tastes than their purchases of popular items. Hence, making use of the tail items is important to predict the user preference in modern recommendation systems. In such sparse observed rating data, the number of items that each user are interested in is very small. Recall the principle of MF methods, each user factor measures how much the user likes items that score high on the corresponding item factor [Koren et al., 2009]. Thus, each user can be characterized well by a few latent features, i.e., the user factor vector in latent feature space should be sparse. This is the main motivation of our work.

For the sake of dealing with sparse rating data, probabilistic matrix factorization (PMF) [Salakhutdinov and Mnih, 2007] was proposed based on the assumption that the latent factors follow Gaussian distribution. This assumption makes sense when the observations are continuous, however it is less justified when the data are on an ordinal scale [Lakshminarayananw et al., 2011]. Later, Bayesian matrix factorization (BPMF) [Salakhutdinov and Mnih, 2008] is estimated by imposing Gaussian-Wishart priors over the factor vectors and tries to generate a non-Gaussian distribution via a product of two multivariate Gaussians. Recently, a sparse covariance prior is adopted in [Shi et al., 2013] to enforce the user and item factors and make each latent feature reflect the semantics more properly. To date, PMF and its variants become the arguably representative CF methods, but they can not guarantee that each user is effectively represented by the most informative latent features. In real applications, obtaining proper user/item representation is crucial to complete the sparse rating matrix especially for the cell values corresponding to the tail items.

We therefore present a sparse probabilistic matrix factorization method (SPMF) by utilizing Laplace distribution to model the item/user factor vector. Laplace distribution [Mohamed et al., 2011] makes the most elements of each factor vector close to zero, which on the one hand is beneficial for SPMF to distinguish the relevant and irrelevant latent features for each user/item. Meanwhile, the tails in Laplace distribution are comparatively heavy, which on the other hand is rewarding to identify the tail items and partially solve the longtail problem. Because Laplace distribution is non-smooth, the Bayesian inference of SPMF is not analytically tractable. In this paper, we express it as a scale-mixture of Gaussian distribution and exponential density, and employ the Markov chain Monte Carlo technique to perform the Bayesian inference. The second distinguishing feature of our work is the distributed Gibbs sampling algorithm for training the SPMF model, which can efficiently deal with large-scale, sparse and very imbalanced dataset such as 1M MovieLens and 100M Netflix data.

The rest of the paper is organized as follows. The related work is described in Section 2. Section 3 gives the proposed sparse probabilistic matrix factorization model and the distributed Gibbs sampling algorithm for model inference. A series of experimental results on real world datasets are listed and discussed in Section 4. A brief conclusion and future work are shown in Section 5.

2 RelatedWork

Our proposed model can be taken as a variant of Probabilistic Matrix Factorization. Let us briefly review PMF and Bayesian PMF.

2.1 Probabilistic Matrix Factorization

Probabilistic Matrix Factorization (PMF) is a probabilistic linear model with Gaussian observation noise [Salakhutdinov and Mnih, 2007] . Figure 1(a) shows the graphical model of PMF using the plate convention. The conditional distribution of observed rating data R ∈ ℜM×N on the latent user factor U ∈ ℜD×M and item factor V ∈ ℜD×N is modeled via

where M is the number of users, N is the number of items, D is the number of latent features. N(x|µ, σ2 ) denotes Gaussian distribution with mean µ and variance σ2 . Iij is an indicator function, Iij = 1 if the ith user is associated with the jth item, otherwise Iij = 0. The user factor and item factor are assumed to followNormal distribution, i.e.,

In [Salakhutdinov and Mnih, 2007], maximizing the logposterior of (1) is implemented by minimizing the sum-ofsquares error function with quadratic regularization terms. However, restricting all features to the same regularization level limits the flexibility of the model. Moreover, it is computationally very expensive to search for appropriate values of the regularization parameters.

2.2 Bayesian Probabilistic Matrix Factorization

Bayesian probabilistic matrix factorization (BPMF) [Salakhutdinov and Mnih, 2008] introduces GaussianWishart priors for the hyperparameters, as shown in Figure 1(b). Among them, the prior distributions over U and V are assumed to be Gaussian:

The hyperparameters ΘU = {µU , ΛU }, ΘV = {µV , ΛV } are assumed to be Gaussian-Wishart distribution. By maximizing the log-posterior of the model over both parameters and hyper parameters, BPMF can automatically control the model complexity. Later, Shi et al. [Shi et al., 2013] proposed sparse covariance matrix factorization (SCMF) by imposing Laplace prior on the covariance matrices of U and V to consider the feature correlations and prevent overfitting. Although BPMF and SCMF obtain promising results, they, like PMF, assume that the factors follow Gaussian distribution, which is not justified especially when the data are extremely sparse.

In the next section, we will give a sparse probabilistic matrix factorization model and a distributed inference algorithms based on Gibbs sampling technique.

3 Sparse Probabilistic Matrix Factorization

To achieve sparsity, the sparsity-favoring distributions can be employed which prefer a high excess kurtosis, i.e., a high peak with heavy tails. The set of sparsity-favoring distributions includes spike-and-slab, Student-t, Laplace, and Gamma distribution [Polson and Scott, 2010]. Among them, Laplace distribution is log-concave, leading to a posterior whose log density is a concave function and has a single local maximum, which is essential to design robust and easy-to-use algorithm [Paninski, 2005]. The probability density function (P.D.F.) of Laplace distribution is defined as:

As shown in Figure 2, Laplace distribution has higher excess kurtosis than Gaussian distribution at the mean point, which means that the variable x has greater chance to concentrate close to zero. Meanwhile, Laplace distribution has two comparative heavy tails allowing for occasional large values, which is beneficial for model the tail items in MF-based collaborative filtering methods.

3.1 SPMFModel

In SPMF, the conditional distribution of the observed ratings Rij is still i.i.d. normal distribution with mean UTi Vj and variance σ2 , as shown in (1). To get a full Bayesian approach, we introduce an inverse Gamma distribution (Γ-1 (α, β)) with shape α and rate β to model the Gaussian noise variance σ2 :

The inverse Gamma distribution is already proved to be effective in modeling the unknown variance of a normal distribution [Witkovsky, 2001].

In order to effectively characterize each user and item, we try to select the most informative latent features to represent them. Thus, we employ Laplace distribution to model the latent user factor U with zero mean and ηU scale, and model the latent item factor V with zero mean and ηV scale.

Here, we introduce the Generalized Inverse Gaussian (GIG) distribution to model the scales (ηU and ηV ) of Laplace distribution for obtaining a fully Bayesian treatment of SPMF and enhancing the model robustness. In [Zhang et al., 2012], it has been demonstrated that Laplace mixture with GIG is beneficial to define a regularizer for variable selection, which is useful for our SPMF model to select the most informative latent features for each item or user. Let x = (ηU )ki, then

where Kγ(·) is a modified Bessel function of the second kind with the index γ. The distribution of (ηV )kj has the similar form.

According to the Bayesian theory, the posterior distribution over latent factors U and V can be modeled as

Then, the overall generative process can be summarized as follows.

• Draw scales (ηU )ki and (ηV )kj from GIG(γ, a, b).
• Draw each user factor vector Ui from L(0, (ηU )i) and each item factor vector Vj from L(0, (ηV )j).
• Draw σ2 from Γ-1 (α, β)
• Draw each observed rating Rij from N(UTi Vj, σ2 ).

and the graphical model of SPMF is shown in Figure 1(c).

3.2 Inference

Markov chain Monte Carlo (MCMC)-based methods are widely applied to approximate the predictive distribution (like

(10)) [Neal, 1993]. Its key idea is to construct a Markov chain that will evenly converge to the posterior distribution of the model with the given data. Each state of the Markov chain is used as sample of the desired distribution. When the conditional distributions can be sampled easily, Gibbs sampling is the simplest but efficient algorithm. In SPMF model, however, it is not easy to sample from a non-smooth Laplace distribution. Fortunately, Laplace distribution can be equivalently expressed as a scaled mixture of Gaussians [Andrews and Mallows, 1974], i.e., an infinite Gaussian mixture with an exponential distribution like

Obviously, all priors on the parameters and hyperparameters in SPMF are conjugate, thus, we can develop an efficient Gibbs sampling algorithm to infer SPMF model. Gibbs sampling sample each variable from its distribution conditional on the current values of all other variables, i.e., it samples each variable by fixing all others.

Sample parameters Ui :

For Ui, we can extract all terms related to Ui and use Bayes' rule to obtain

In order to incorporate the alternative expression of Laplace distribution, infinite Gaussian distribution as shown in (11), we introduce a vector (λU )i ∈ RD , where each element (λU )ki is a latent variable with exponential prior, i.e.,

for the corresponding Uki. In this case, we can further express (12) as

where

Sample hyperparameters (λU )i and (ηU )i :

When introducing the vector λU , we have p((λU )i|Ui, (ηU )i ) ∝ p(Ui|(λU )i, (ηU )i)p((λU )i) and each element (λU )ki has an exponential prior. According to the property of exponential distribution, (λU )-1 ki follows an inverse Gaussian distribution (denoted as G-1 ) [Zhang et al., 2012]. Then, we can get the conditional distributions of the hyperparameters (λU )i as follows.

For the scale of Laplace distribution, (ηU )i, it follows GIG distribution as shown in (9), i.e.,

Since it is inefficient to do sampling directly from GIG distribution, we convert the posterior distribution to a special case of GIG distribution by setting γ = -1/2. In this case, GIG distribution becomes an inverse Gaussian distribution, thus we can sampling (ηU )i with

Sample parameters σ2 :

Based on (17), the posterior distribution of σ2 by fixing other variables can be written as

The sampling scheme about parameters Vj and hyperparameters (λV )j and (ηV )j is similar to Ui, (λU )i and (ηU )i respectively.

3.3 Distributed Sampling

Note that the columns of U and V are independent, we can sample them in parallel. More precisely, when sampling U, the i-th column of U, can be obtained by using only the i-th column of rating matrix R, which holds the i-th user's preference, and all the columns of V. In such a setting, we can use a parallel broadcast-join [Blanas et al., 2010] to efficiently infer SPMF model.

As shown in Figure 3,we broadcast the factor V to all participating machines (there are total C machines), which create a hashtable for its contents. The rating data R are stored in the distributed filesystem (DFS) partitioned by its rows (corresponding to users) and forms the input for the map operator. Here R(1) refers to partition 1 of R. The map operator reads the i-th row of R and selects all columns of V from the corresponding hashtable. Next, the map operator can sample Ui via (13), (λU )i via (15), (ηU )i via (16), σ2 via (17), and write back its result. This strategy contains multithreaded mappers that leverage all cores of the worker machines for sampling U.

The computational complexity of inferring the parameter U and its corresponding hyperparameters is O(T( ˆMD3 + ˆ CD2 )), where T is the number of iterations, D is the number of latent features, ˆM and ˆ C are the maximum number of users and the maximum observed ratings in all machines. The sampling process for V works analogously and its computational complexity is O(T( ˆND3 + ˆ CD2 )), here ˆN is the maximum number of items in all machines. Thus, the total cost of inferring process is O(T( ˆMD3 + ˆND3 + ˆ CD2 )). Obviously, the computational complexity of our algorithm is linearly scalable to the size of rating data matrix, thus, it is practical for handling large-scale dataset.

4 Experimental Results and Discussions

In this section, we compare SPMF with other baseline methods PMF and BPMF using the MovieLens and Netflix datasets.

4.1 Methodology

Datasets: In our experiment, we use MovieLens 1 data with 1M ratings and Netflix 2 data with 100M ratings, where the ratings are ordinal values on the scale 1 to 5. These two datasets are extensively used in the literature to test the performance of recommendation systems. The statistical properties of the datasets are summarized in Table 1. Density indicates the percentage of non-zero cells in rating matrix. Netflix is more challenging than MovieLens. In Netflix, there are over 30% of the items having less than 10 ratings.

Parameter setting: In PMF, the regularization parameters are tuned from the candidate set {10-3 , 10-2 , 10-1 , 1, 101 , 102 , 103 }. Following [Salakhutdinov and Mnih, 2008], we set µ0 = 1, v0 = D, W0 = ID×D, and variance σ2 = 1/2 for BPMF. We employ the code 3 given by the authors .

In the proposed SPMF, we adopt an uninformative prior for the noise variance, i.e., initializing the shape and rate of inverse Gamma distribution via α = β = 0. Meanwhile, the prior of inverse Gaussian distribution on hyperparameters ηU is initialized by setting a = 1/D and b = D, here a is a part of scale in (16) which plays an important role to control the sparsity. Smaller a leads to more sparse factor vector. The reason that we use this setting is that the size of latent features (D) in SPMF may affect the sparsity of latent factor vectors. For example, when D is small, i.e., there are few features, then each latent factor vector may be related to most features, which means that the factor vector will be dense. Otherwise, larger D may lead to sparse factor vector. Thus, we can say D is inverse proportional to the sparsity of factor vector.

1 http://www.grouplens.org/node/73

2 www.netflixprize.com

3 http://www.cs.toronto.edu/ rsalakhu/BPMF.html

Evaluation: The Root Mean Square Error (RMSE) is used to measure the performance of our proposed SPMF and the baselines (PMF and BPMF).

where |Ω| is the number of observed ratings in testing data, ̂Rij is the predicated rating value of the i-th user for the j-th item. Smaller RMSE indicates better recommendation.

4.2 Results and Discussions

In experiments, we randomly select 90% data as training data and the remaining ratings as testing data for ten times, and the average results on testing data are recorded. When selecting training data, every item and user appears at least once.

Effect of latent feature size

We are firstly interested in evaluating the effect of the number of latent features (D) on the proposed SPMF and baselines (PMF and BPMF). Table 2 shows the performance of these models with different latent feature sizes on two datasets.

It can be seen that SPMF and BPMF outperforms PMF, which indicates that automatically learning the regularization parameters is beneficial to enhance the model robustness, esp., avoid overfitting problem. With the increase of D, the performances of both SPMF and BPMF steadily improve. To be excited, SPMF consistently performs better than BPMF on all latent feature sizes, and larger D leads to bigger improvement (e.g., SPMF outperforms BPMF by over 1.09% to 1.54% on MovieLens data, 2.73% to 3.46% on Netflix from D = 10 to 150). Furthermore, we show the convergence process along iterations (taking Netflix as an example) in Figure 4, which verifies that SPMF like BPMF does not overfit. The factors of PMF are randomly initialized. The factors in BPMF and SPMF are initialized with the MAP estimations obtained by PMF. Since BPMF and PMF assume that the factors follow Gaussian distribution, so that BPMF starts better than SPMF. The main difference between SPMF and BPMF is the prior over the latent factors, which demonstrates that imposing a sparse-favoring distribution on latent factors does bring performance gain especially on sparse Netflix data.

Even though larger D can lead to better performance for SPMF and BPMF, it also increases the model complexity because running the Gibbs sampling on more parameters is computationally much more expensive. As listed in Table 3, more running time is needed by both SPMF and BPMF when D grows. All experiments are conducted on the PCs with Intel Core i3-2120 3.3GHz processor and 4 GB RAM. Here SPMF-1 indicates non-distributed sampling. SPMF-4 and SPMF-8 use four and eight machines for distributed sampling respectively. The change of RMSE (less than 0.00001) is taken as the stop condition for both methods. Obviously, distributed sampling in SPMF is helpful to speed up the inferring process.

Long tail recommendation

In modern recommendation systems, only a small part of items are popular, while a large part of items have few rating information. However, such unpopular items, also called as long tail items [Anderson, 2006], often play an important role for recommendation, esp., in commercial world. For example, Amazon successfully makes most of their profit from the long tail products rather than the best selling products.

In order to investigate how our proposed SPMF method deal with the long tail items, we firstly ranking items according to their rating frequency, and then almost evenly separated them into ten groups. Figure 5(a) demonstrates the number of ratings in ten item groups for MovieLens dataset, where the x-axis refers to the boundary of rating frequency in each group. Obviously, they follow long tail distribution. Figure 5(b) and (c) give the average RMSE of each item group for MovieLens and Netflix respectively, which are obtained by recommendation results SPMF, BPMF and PMF with setting D = 50. As expected, the proposed SPMF model outperforms two baselines especially for the items in the long tail part.

This result further verifies that SPMF integrating with sparsity-favoring distribution (Laplace distribution) on latent factors has ability to handle sparse rating data. The main reason is that BPMF imposes Gaussian prior on V, while SPMF assumes V is generated from Laplace distribution. Such sparsity prior can effectively partition the elements (in each factor vector) into a large set which is close to zero with high probability and a small set which has significant mass on large values. Thus, SPMF can model the long tail items well and significantly enhance the recommendation performance especially for unpopular items.

5 Conclusions and FutureWork

In this paper we propose a sparse probabilistic matrix factorization model and a distributed Gibbs sampling algorithm to infer the model. The experimental results have demonstrated that the model can be efficiently trained and successfully applied on sparse and large-scale data by comparing with the state-of-the-art MF-based CF methods. One distinguishing characteristic of our model is that it has ability to deal with long tail items. This makes SPMF have potential commercial value, e.g., it can be used to explore the long tail market and boost the one-stop shopping convenience.

The proposed model benefits from the Laplace distribution for modeling the latent factor, however, we have to empirically tune the size of latent features. It will be nice if we can automatically find the proper D. One more thing is that the Gibbs sampling method for log-linear model has to cost large memory, thus, we will transfer it into robust MF problem [Wang et al., 2012] to reduce the space complexity. Another possible direction is to extend the model and make it able to do cold-start recommending which is a challenging issue [Houlsby et al., 2014]. Last but not least, it is interesting to apply the model in other domains such as e-commerce with large-scale and online data.

Acknowledgments

This work was supported in part by the NSFC Grants 61375062 and 61370129, PCSIRT Grant IRT201206 and the Fundamental Research Funds for the Central Universities 2014JBM029 and 2014JBZ005.