This paper proposes a hierarchical probabilistic model for ordinal matrix factorization by actively modelling the ordinal nature of ranking data, which is typical of large-scale collaborative filtering tasks. Two algorithms are presented for inference, one based on Gibbs sampling and one based on variational Bayes. The model is evaluated on a collaborative filtering task, where users have rated a collection of movies and the system is asked to predict their ratings for other movies. The Netflix data set is used for evaluation, which consists of around 100 million ratings. Using root mean-squared error (RMSE) as an evaluation metric, results show that the suggested model improves similar factorization techniques. Results also show how Gibbs sampling outperforms variational Bayes on this task, despite the large number of ratings and model parameters.
Matrix factorization is a highly effective technique for both predictive and explanatory data modeling. An observed data set is represented as an M ×N matrix of results, R, where rows represent the observation number and columns represent the variables of interest. This matrix is then decomposed as R = U⊤V + ǫ, where U is a K × M matrix, V is a K × N matrix, and ǫ is a noise term. U and V represent the values of explanatory variables which, when multiplied and added, give a predictor of the values in R.
This paper discusses a model for matrix factorization when the observed data
is a collection of ranks, also known as ordinal data. It contributes:
– An extended probabilistic model for ordinal matrix factorization. The
additional parameters express a rich model, through which the prior distributions
over factors are flexible, but remain coupled through shared hyper-priors.
Salakhutdinov and Mnih’s hierarchical framework of matrix factorization [
– An efficient variational approach for inference in the model.
– An efficient Gibbs sampling algorithm for inference in the model.
The model is tested on the Netflix data set, which is a collection of around
100 million movie ratings. The matrix is sparse in that many of the movie-user
pairs have no rating. Performance is computed by estimating missing values
and computing the root mean-squared error (RMSE) on a held-out collection of
ratings. The model is shown to improve or equal alternative matrix factorization
approaches on this metric. There exists a rich ensemble of similar probabilistic
factorization models for large-scale collaborative filtering [
Ordinal regression arises when the possible observed values in R are ranked. For example, in a collaborative filtering task an item with a five-star rating (r = 5) is regarded as superior to one with a four-star rating (r = 4), which in turn is better than one with a three-star rating. Instead of directly modelling the ranks, r = 1, . . . , R, we will model a collection of hidden variables, h. The probability distribution of the ranks will depend on the position of h compared to a collection of fixed boundaries br,
−∞ = b1 < b2 < . . . < bR+1 = +∞ .
Let rmn denote the rank for row m = 1, . . . , M and columns n = 1, . . . , N , or user n’s rank for item m. The probability of rmn in terms of a hidden variable hmn is
p(rmn|hmn) = Φ(hmn − br) − Φ(hmn − br+1) , (1) where Φ(x) = R−x∞ N (z; 0, 1) dz is the cumulative Gaussian density or probit function. The hidden variables, hmn, are modelled by a Gaussian distribution, with mean equal to the dot product of a row factor um and a column factor vn, p(hmn|um, vn, γ) = N (hmn ; u⊤mvn, γ−1) . (2) Factors um and vn can be viewed as item m and user n’s latent traits, each of which has a Normal prior distribution. The um’s are conditionally independent given a shared mean and covariance,
p(um|μu, Ψ u) = N (um; μu, Ψ u−1) .
A completely analogous model is used for the user factors vn. The mean and precision matrix (inverse covariance) is modelled with a conjugate Normal-Wishart prior p(μu, Ψ u) = N (μu; μ0, (β0Ψ u)−1) W(Ψ u; W0, ν0) . (3) ν0 W0
The Wishart distribution W(Ψ ; W, ν) ∝ |Ψ |(ν−K+1)/2 exp(− 21 tr[W−1Ψ ]) over symmetric positive definite matrices is parameterized with a scale matrix W and ν degrees of freedom. The inverse variance parameter γ is non-negative and is modelled with its conjugate Gamma distribution p(γ; a0, b0) = Γ (γ; a0, b0), where Γ (γ; a, b) ∝ γa−1 exp(−γ/b). The hyperparameters {β0, W0, ν0} and {a0, b0} have to be specified by the user. Figure 1 summarizes the joint distribution of the data and model parameters θ = {H, U, V, γ, μu, Ψ u, μv, Ψ v} in a Bayesian network. In Section 3 we will sample from, and in Section 4 approximate the posterior distribution p(θ|D). 3
Gibbs sampling is a MCMC method that sequentially samples from the conditional distributions of the model. We briefly describe two sampling steps here. Factors. With Ω(m) being the set of users that rated item m, the conditional distribution for each item factor um is Gaussian
um ∼ N um; Σm Ψ uμu + γ
X n∈Ω(m)
hmnvn , Σm with covariance Σm = (Ψ u + γ Pn∈Ω(m) vnvn⊤)−1. Notice that the distribution for factor m only requires knowledge of γ and the variables directly connected with m: {(vn, hmn|n ∈ Ω(m)}. Latent variables. Sampling the conditional for the latent variable p(hmn|rmn, um, vn, γ) ∝ p(rmn|hmn) p(hmn|um, vn, γ) requires one evaluation of a unit interval random number, one random Normal number, two evaluations of Φ and one of Φ−1. To derive the sampler we introduce the “noise-free” latent variable Φ(h − b) = R N (f ; h, 1) Θ(f − b) df . For any m and n, which is omitted here for brevity, the joint marginal distribution of r, f , and h, given μ = u⊤v and γ, is
p(r|f ) p(f |h) p(h|μ, γ) = hΘ(br+1 − f ) − Θ(br − f )i N (f ; h, 1) N (h; μ, γ−1) . (5) The density f, h|r, μ, γ in (5) can be sampled from in two steps, f |r, μ, γ and h|f, μ, γ. The distribution f |r, μ, γ is a truncated Normal p(f |r, μ, γ) =
N (f ; μ, 1 + γ−1) Θ(br+1 − f ) − Θ(br − f ) Φmax − Φmin with Φmax = Φ((br+1 − μ)/p1 + γ−1) and Φmin = Φ((br − μ)/p1 + γ−1)). A sample can be drawn from p(f |r, μ, γ) using the cumulative distribution f = μ + p1 + γ−1 Φ−1(Φmin + rand (Φmax − Φmin)), where “rand” gives a uniform random number between zero and its argument. The desired sample is obtained from p(h|f, μ, γ) = N h ; (f + γμ)/(1 + γ), (1 + γ)−1 . 4
Variational Bayes (VB) is one popular algorithm for obtaining an approximation q(θ) to p(θ|D). At each stage in the algorithm, one of the approximating factors is chosen, all other factors are fixed and the algorithm minimizes the Kullback-Leibler divergence KL(q(θ)kp(θ|D)) by varying the given factor. An approximating family q(θ) =
Y q(hmn) Y q(um) Y q(vn) q(μu, Ψ u) q(μv, Ψ v) q(γ) (m,n) m n is chosen, with the item and user factors having Gaussian approximations q(um) = N (um; humi, Σm) and q(vn) = N (vn; hvni, Ξn). The factorized approximations for the hierarchical parameters follow the prior’s Normal-Wishart form, while q(γ) is chosen to be a Gamma density. We present two example VB updates: Factors. The full update for each of the item factors is
" q(um) = N um; Σm hΨ uμui + hγi
X n∈Ω(m)
# hhnmihvni , Σm , with Σm = (hΨ ui + hγi Pn∈Ω(m)(Ξn + hvnihvn⊤i))−1. Latent variables. The mean and variance of hmn are determined when needed in other updates from q(hmn) ∝ p(rmn|hmn) exp
D log p(hmn|um, vn, γ)E q(um) q(vn) q(γ) .
Define μ = hu⊤mihvni, and let γ be short for hγi, and N (z) short for N (z; 0, 1). If br and br+1 are the boundaries associated with rmn, and zr = (μ−br)/p1 + γ−1, the explicit expression for hhnmi is (hh2nmi can be similarly evaluated) hhmni = μ + γ−1
N (zr) − N (zr+1) .
p1 + γ−1 Φ(zr) − Φ(zr+1) 5
The proposed model is evaluated by testing its predictive performance on the Netflix data set. The data set contains around 100 million ratings from N = 480, 189 users on M = 17, 770 movie titles. Each rating is a number of stars 1 to 5, and is used as the ranked observation at the relevant point in R. A test or “qualifying” set of almost three million user–movie pairs for which the ratings are withheld. Algorithms are benchmarked by their RMSE computed over an unknown half of the test set.
Performance results We compare inference with Gibbs sampling and VB for a
number latent dimensions, K = 50, 100, 200, and γ settings γ = 0.8, 0.9, 0.10,
and one setting in which γ is also inferred. The results are summarized in Figure
2. When the latent factor dimensionality K is increased, a higher precision γ
gives better performance. The Gibbs-sampled results provide further evidence
that proper regularization in models with far more parameters than the data set
size |D| is possible with Bayesian averaging [
This paper has proposed a hierarchical model for ordinal matrix factorization.
Comparing our results to the regular factor model of [
Two standard methods of inference were compared: Gibbs sampling and variational Bayes. The comparison of the two approaches on such a large task gives rise to several conclusions that are applicable to other similar settings (factor models with relatively high noise levels). Variational inference for the smaller factor sizes (K = 50) showed promising results, but unfortunately overfitted for
Gibbs sampling 100
150 γ = 0.09 γ = 0.1 200 larger models. On this problem, Gibbs sampling performed better and would generally be recommended.
In this paper a minimal model was used to keep the message as clean as
possible. The results can clearly be improved by blending with other models,
as was shown by [