SCMF: Sparse Covariance Matrix Factorization for Collaborative Filtering

1 Introduction

Nowadays, recommender systems complement search engines for information filtering on the Internet. They provide personalized recommendation for items such as products or services by taking into consideration such sources of information as characteristics of the items and preferences of similar users. As two well-known examples, Amazon and Netflix use recommender systems to recommend to their users books and movies, respectively. For this approach to work effectively for personalized advertising which leads to commercial benefit, it is essential that user preference can be predicted as accurately as possible. This explains why research in recommender systems has drawn so much interest from both the academic and commercial worlds.

While recommender systems can use either the contentbased or collaborative filtering (CF) approach (or a hybrid approach incorporating the two), the CF approach is considered particularly promising in this era of the social web (or called Web 2.0) in which social relations are often exploited extensively for a variety of applications. Unlike content-based filtering which defines user or item similarity based on contentbased features extracted from user profiles or item descriptions, the CF approach predicts user preference by analyzing the previous rating data of many users (hence the qualifier 'collaborative').

Memory-based CF methods [Breese et al., 1998; Linden et al., 2003] directly compute user or item similarity using the rating data. These methods have a longer history and have been used by many commercial systems. On the other hand, model-based CF methods make use of machine learning and data mining models to explore patterns in the rating data. They include the clustering model [Xue et al., 2005], aspect model [Hofmann, 2004], latent factor model [Canny, 2002], restricted Boltzmann machine [Salakhutdinov et al., 2007], etc. Over the past few years, extensive empirical studies conducted on large-scale datasets show that model-based CF methods generally have higher accuracy than memorybased methods [T¨oscher et al., 2009].

The most popular and successful model-based CF methods are based on low-rank matrix factorization (MF). The rationale behind the MF approach is that the observed rating data can be explained holistically by a small number of latent factors for both users and items, which, for instance, may be related to user groups and item topics. Given an incomplete and highly sparse rating matrix representing the observed rating data, a typical MF method factorizes the matrix into two low-rank matrices, one for latent user features and the other for latent item features. Multiplying the two matrices thus allows us to complete the original incomplete matrix by filling in the missing entries.

Among these MF methods, probabilistic matrix factorization (PMF) [Salakhutdinov and Mnih, 2008b] is arguably the most representative one due to its high efficiency and accuracy. Many other MF methods are variants of PMF. For example, [Hu et al., 2008] utilized the rated-or-not indicator for implicit feedback to improve the results; [Koren, 2008] unified the latent factor MF model with the neighborhood model under the same framework; and [Rendle and Schmidt-Thieme, 2008] generalized the standard inner product operation for user features and item features to various kernel functions.

Recently, incorporating social network information into MF models has emerged as a hot research direction [Ma et al., 2008; 2009; 2011].

Since the rating matrix contains many missing entries, MF methods usually impose an ℓ2 penalty on both the latent user features and latent item features to avoid overfitting. However, restricting all features to the same regularization level limits the flexibility of the model. One variant [Lakshminarayanan et al., 2011] uses a more general diagonal covariance matrix whose diagonal entries do not have to be the same to give different regularization levels to different features, if we assume that the user and item features are generated from multivariate normal distributions. Nevertheless, it still assumes that the features are uncorrelated. An even more flexible variant uses a full covariance matrix, as in Bayesian probabilistic matrix factorization (BPMF) [Salakhutdinov and Mnih, 2008a]. Although BPMF gives good results, its high computational requirement, due to the use of sampling-based inference techniques, makes it challenging to be applied to very large datasets.

An open question is whether there is indeed the need for a full covariance matrix. To answer this question, let us first take a look at the covariance matrix for the user and item features learned by PMF. Figure 1 shows the covariance matrix and correlation matrix for a movie recommendation dataset (ML-10M with 80% training data). Not only are the diagonal entries in the covariance matrix significantly positive, some non-diagonal entries also differ from 0 significantly although such entries are sparse. Note that the latent features are generally found to be semantically meaningful. For example, the latent features may correspond to such topics as biography, comedy, history, and romance. If two features, e.g., those corresponding to biography and history, are highly correlated, we expect that the corresponding entries in the (symmetric) covariance matrix will also be significantly non-zero. Nevertheless, most other pairs of features are uncorrelated or at most weakly correlated. This explains why the majority of the non-diagonal entries are close to 0. Similar observations have also been noted in some other datasets. This finding inspires us to propose later in this paper a sparse covariance prior which favors a sparse yet non-diagonal covariance matrix for MF models.

Besides its role in reflecting the semantics more properly, we anticipate that imposing sparsity on the covariance matrix also has a side effect of preventing overfitting. When there exist outliers in the very sparse rating matrix, two otherwise uncorrelated features may exhibit high correlation incorrectly. Imposing sparsity on the covariance matrix will likely help to alleviate this problem.

Based on the inspiration above, we propose a probabilistic generative MF model with a novel sparse covariance prior. The same sparse covariance prior is imposed on both the user and item feature vectors. For model inference, we take the maximum a posteriori (MAP) approach for its efficiency and scalability. Our model is compared empirically with two baseline MF methods on some large datasets.

2 Brief Review of PMF

Since our model builds on PMF [Salakhutdinov and Mnih, 2008b], let us first give a quick review of the model formulation and learning algorithm of PMF.

For the sake of presentation, we assume that the recommender system is for movie recommendation. The same model may also be used for recommendation of other items. Let there be M movies and N users. The rating of movie j by user i is denoted by Rij . Like many other MF methods, PMF uses two latent feature matrices U ∈ RD×N and V ∈ RD×M , with column vectors Ui and Vj corresponding to the latent feature vectors of user i and movie j, respectively. It is assumed that the observed ratings are governed by a normal distribution to give the likelihood as follows:

where Zij is an indicator variable introduced to indicate whether Rij is observed, and N(X | μ, σ2 ) means that X satisfies a normal distribution with mean μ and variance σ2 (the same notation is also used for the multivariate normal distribution with mean vector μ and covariance matrix Σ). Zero-mean spherical normal distributions are imposed as priors on the column vectors of U and V independently:

where I is an identity matrix and σ2U and σ2V are variance parameters. Using a MAP estimation approach for model inference, PMF solves a log-posterior maximization problem which is equivalent to the following minimization problem:

where λU = σ2 /σ2U and λV = σ2 /σ2V are regularization parameters. The graphical model of PMF is shown in plate notation in the left panel of Figure 2. To solve the optimization problem in Eq. (3) efficiently, PMF uses stochastic gradient descent to update U and V iteratively.

3 Sparse Covariance Matrix Factorization

We are now ready to present our MF model. We first present it as a probabilistic generative model with a sparse covariance prior. For this reason, we refer to it as sparse covariance matrix factorization (SCMF). As in PMF, model inference in SCMF is formulated as a MAP estimation problem which corresponds to an equivalent minimization problem. We then present an extension of the basic model as well as its complexity analysis.

3.1 Probabilistic Model

The likelihood term of SCMF is exactly the same as that of PMF as given in Eq. (1). The main difference lies in the prior distributions of the latent user and item feature vectors. Specifically, the priors of the feature vectors Up and Vp are zero-mean multivariate normal distributions with a full covariance matrix:

Unlike PMF, using a full covariance matrix allows the features to be correlated. Having the same covariance matrix Σ for both U and V in Eq. (4) provides a bridge to connect the two feature matrices.

As discussed above, it is desirable to impose sparsity on the non-diagonal entries of the covariance matrix both for reflecting the underlying semantics properly and for preventing overfitting. To achieve this, we impose a Laplace prior on every element of Σ:

where λ is a parameter that controls the degree of sparsity, p.d. stands for positive definite, which is a property required for all covariance matrices, 1 and L is the Laplace distribution defined as

1 Although covariance matrices may be positive semi-definite (p.s.d.), we impose the stricter p.d. requirement here because the probability density will be undefined for p.s.d. covariance matrices. Nevertheless, the p.d. requirement is easy to satisfy because the latent dimensionality is typically low.

where |a| denotes the absolute value of the scalar a. The generative model of SCMF is shown in the right panel of Figure 2 to allow us to reveal the similarities and differences between it and PMF.

To understand how model inference is done, we first take a look at the joint posterior distribution of the latent variables and parameters U, V, Σ given the hyperparameters θ = (σ, λ):

subject to the requirement that the p.d. property is satisfied. Maximizing the (log-)posterior in Eq. (7) is equivalent to the following minimization problem:

where Σ should be a p.d. matrix, det(·) denotes the determinant of a matrix, tr(·) denotes the trace of a matrix, P is a matrix with zeros in the diagonal entries and ones in the nondiagonal entries to prevent the shrinking of diagonal entries in Σ, and ⊗ denotes element-wise multiplication.

3.2 Optimization Problem

To solve the optimization problem in Eq. (8), we now present a stochastic gradient descent algorithm which alternates the update of U, V and Σ.

Optimizing with respect to U and V

To optimize with respect to U while keeping V and Σ fixed, the problem in Eq. (8) becomes

For each observed Rij, the gradient with respect to Ui is

Similarly, the gradient with respect to Vj is

With these gradients, we can use stochastic gradient descent to update Ui and Vj iteratively.

Optimizing with respect to Σ

To optimize with respect to Σ while keeping U and V fixed, the problem in Eq. (8) becomes

where

is the sample covariance matrix for the column vectors of U and V . We note that the p.d. property of Σ is preserved since the log-determinant term acts as a barrier function [Boyd and Vandenberghe, 2004].

While it appears challenging to solve the optimization problem in Eq. (12) since the objective function is both nonconvex and non-smooth, it is in fact well structured because tr(Σ-1S) + α‖P ⊗ Σ‖1 is convex with respect to Σ while logdet(Σ) is concave. Problems of this type have been widely studied and can be solved using the majorizationminimization (MM) approach [Hunter and Li, 2005].

Theorem 1. The objective function value of Eq. (12) is nonincreasing under the following update rule:

where α is a step size and T is an element-wise soft-thresholding operator defined by Tij{X, Y } = max(0, |Xij| - Yij). The objective function value remains unchanged under the update rule if and only if Σ is at a stationary point.

Theorem 1 guarantees the convergence of updating via Eq. (14). In what follows, we will give a detailed proof of Theorem 1. We start by defining an auxiliary function in a way similar to that in [Seung and Lee, 2001].

Definition 1. g(x | x′ ) is an auxiliary function for f(x) if the following conditions are satisfied:

Lemma 2. If g(x | xt) is an auxiliary function, then f(x) is non-increasing under the update rule

Therefore, the equality f(xt+1) = f(xt) holds only when xt is a local minimum of g(x | xt). By iteratively updating via Eq. (14), we can reach a local minimum of Eq. (12). We will show in the following that Eq. (14) decreases the value of the auxiliary function of Eq. (12).

Lemma 3. The following is an auxiliary function of Eq. (12):

Proof. It is trivial to prove directly from the definition that g(x | x) = f(x). To prove that g(x | xt) ≥ f(x) with Eq. (17) and Eq. (12), we note that

This is true since logdet(Σt+1) is concave with respect to Σt+1 and log det(Σt)+tr((Σt)-1 (Σt+1-Σt)) is the tangent of logdet(Σt+1) at the point Σt.

Lemma 4. The update rule in Eq. (14) can guarantee a nonincreasing value for the auxiliary function in Eq. (17).

Proof. Eq. (17) consists of a differentiable convex component and a non-differentiable component with an ℓ1 penalty. Thus, we can employ a generalized gradient descent algorithm [Beck and Teboulle, 2009], which is a natural extension of gradient descent to non-differentiable objective functions. The objective function in Eq. (17) is then changed to

Solving Eq. (19) is equivalent to updating with Eq. (14).

Proof of Theorem 1. Since Eq. (17) is an auxiliary function, f(Σ) is non-increasing, according to Lemma 2, when the function value of Eq. (17) decreases. It follows from Lemma 4 that Eq. (14) minimizes the auxiliary function. This completes the proof.

3.3 Incorporating User Bias and ItemBias

Previous study [Koren, 2008] showed that incorporating user bias and item bias terms can lead to further improvement because it can allow for the intrinsic difference between users and the intrinsic difference between items to be represented. Here we present an extension of our basic model to include such terms.

Specifically, the likelihood in Eq. (1) is modified to

where bui , the ith element in a vector bu ∈ RN×1 , represents the bias for user i, and, similarly, bvj is the bias for item j.

The bias vectors bu and bv also follow a zero-mean normal distribution, namely,

Given the probability distribution, the optimization procedure for the extended model can be obtained easily by extending from that in Section 3.2.

3.4 Discussion and Analysis

Positive definite property

The covariance matrix Σ, initialized to the identity matrix I, remains symmetric and p.d. throughout the optimization procedure. For the symmetric property, it is easy to see that every update via Eq. (14) does not change this property. For the p.d. property Σ ≻ 0, it suffices to ensure that Σ ⪰ δI for some δ > 0. Not only can this enhance numerical stability but it can also avoid penalizing some features excessively. Every time when Σ is updated, its eigenvalues are checked. In case the condition Σ ⪰ δI is violated as a result of the update, those (possibly negative) eigenvalues smaller than δ will be rounded up to δ. This corresponds to projecting the solution to the convex set C = {Σ | Σ ⪰ δI} in order to preserve the p.d. property.

Convergence

Our optimization procedure involves three subproblems in each iteration corresponding to optimization with respect to U, V and Σ. Because the subproblems for U and V are convex, convergence to the optimal solutions is guaranteed. From Theorem 1, the subproblem for Σ, though not convex, is guaranteed to decrease the objective function value in each iteration. Given that the objective function in Eq. (8) is bounded below, it always converges to a local minimum.

Complexity analysis

Algorithm 1 summarizes the optimization procedure for SCMF, where β denotes the step size parameter for the gradient descent update of U and V .

The algorithm has a time complexity of O((C + M + N)D2 ) for each epoch, where C is the number of observed ratings. The complexity is linear to the numbers of users, items, and ratings, but quadratic to the latent feature dimensionality. However, in practice D is very small with D ≪ M and D ≪ N and hence the algorithm scales linearly with the data size. This gives the algorithm potential to be applied to very large datasets. How to select D is a commonly encountered problem for MF-based methods. Usually it has positive correlation with the sample size and desired accuracy and has negative correlation with computational efficiency.

The actual performance of our method compared to PMF depends on the batch size for updating Σ. In practice, the time requirement of our method is about 1.2-1.5 times that of the baseline PMF method.

Algorithm 1 Optimization procedure for SCMF

• 1: input: rating matrix R; hyperparameters α, β, θ;
• 2: repeat
• 3: for each observed Rij do
• 4: update Ui via Ui = Ui - β ∂f ∂Ui ;
• 5: update Vj via Vj = Vj - β ∂f ∂Vj ;
• 6: update bui and bvj ;
• 7: end for
• 8: update Σ via Eq. (14) until convergence;
• 9: until validation error stops decreasing;
• 10: return U and V .

4 Experiments

In this section, we will shift our focus to the empirical aspect of the paper. We first conduct an empirical analysis of the covariance sparsity property of SCMF. Then, we compare it with baseline methods PMF [Salakhutdinov and Mnih, 2008b] and biased matrix factorization (biased MF) [Koren, 2008] using the MovieLens 2 and Netflix 3 datasets.

4.1 Experimental Settings

There are three MovieLens datasets: ML-100K with 100K ratings, ML-1M with 1M ratings, and ML-10M with 10M ratings. The Netflix dataset is even larger with about 100M ratings. Since our focus in this paper is on CF using rating data, other sources of information such as implicit feedback and tagging information are discarded.

We use root mean squared error (RMSE) as the evaluation measure. It is defined as

where N is the total number of ratings in the test set, Rij is the ground-truth rating of user i for item j, ˆRij denotes the corresponding predicted rating, and Z T ij is an indicator variable used to select the Rij defined in the test set. For all three methods, we set the regularization parameters via 5-fold cross validation.

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

3 http://www.netflixprize.com/

4.2 Empirical Analysis of Covariance Sparsity

We use the MovieLens ML-1M and ML-10M datasets and train SCMF on 80% of their data. The latent feature dimensionality is set to 10. Figure 3 depicts the correlation matrices derived from the covariance matrices learned. Among the off-diagonal entries, only a very small fraction show relatively high correlation. Compared with the correlation matrix learned using PMF (Figure 1), incorporating a sparse covariance prior increases the degree of sparsity of the covariance matrix learned.

4.3 Quantitative Results of Comparative Study

For the three MovieLens data sets, we try different settings by using different percentages (99%, 80%, 50%) of the data for model training. The remaining data are used for testing to obtain the RMSE results. The training data are randomly selected from each dataset. This process is repeated five times for each percentage setting and the RMSE reported is the average over these runs.

Table 1 summarizes the results. In practice, SCMF ends within 200 iterations. We can see that SCMF consistently outperforms both PMF and biased MF on all settings of all three datasets. As the main difference between SCMF and the other two methods is in the structure of the covariance matrix of the latent features, this shows that imposing a sparse covariance prior does bring about performance gain as anticipated.

For the Netflix dataset, we use the training set data for training and probe data for testing. Since the Netflix dataset is even much larger than the largest MovieLens dataset, we set the latent feature dimensionality to larger values, D = 20, 50, 80, following the settings in [Gantner et al., 2011]. The results are shown as convergence curves in Figure 4. Although SCMF converges a bit more slowly than the other two methods when D is larger, it can always reach a lower RMSE and hence the model learned has higher prediction accuracy.

5 Conclusion and FutureWork

In this paper, we have proposed an MF model with a sparse covariance prior for CF applications. The model does not take either the extreme of assuming that all latent features are uncorrelated or the other extreme that all are correlated. Instead, it allows a small number of correlated feature pairs. Empirical studies show that imposing this particular prior can lead to performance improvement.

It should be noted that incorporating a sparse covariance prior on the latent features is not limited to CF applications that use rating data only. Other sources of information such as content-based features extracted from user profiles or item descriptions may also be incorporated into the MF model. This is one direction in which the current research will be extended. Moreover, the current MAP estimation approach, as a point estimation approach, may not be robust enough for very sparse observed data. A full Bayesian approach to MF [Salakhutdinov and Mnih, 2008a] provides a potential solution to this problem. However, it is very challenging to develop efficient and highly scalable full Bayesian models especially when non-conjugate priors are used. As another possible direction in which this research can be extended, some special structure of the problem will be exploited to devise a full Bayesian model with a sparse covariance prior which is a non-conjugate prior. Last but not least, we will apply SCMF to extremely large datasets by exploiting parallel and distributed computing platforms.

Acknowledgments

This research has been supported by General Research Fund No. 621310, No. 413110 and No. 413212 from the Research Grants Council of Hong Kong.