Accounting for missing ratings in available training data was recently shown [3, 17] to lead to large improvements in the top-k hit rate of recommender systems, compared to state-of-the-art approaches optimizing the popular rootmean-square-error (RMSE) on the observed ratings. In this paper, we take a Bayesian approach, which lends itself naturally to incorporating background knowledge concerning the missing-data mechanism. The resulting log posterior distribution is very similar to the objective function in [17]. We conduct elaborate experiments with real-world data, testing several variants of our approach under dieffrent hypothetical scenarios concerning the missing ratings. In the second part of this paper, we provide a generalized probabilistic framework for dealing with possibly multiple observed rating values for a user-item pair. Several practical applications are subsumed by this generalization, including aggregate recommendations (e.g., recommending artists based on ratings concerning their songs) as well as collaborative lfitering of sequential data (e.g., recommendations based on TV consumption over time). We present promising preliminary experimental results on IP-TV data.
Algorithms ∗ This work was done while an intern at Bell Labs, AlcatelLucent.
Copyright is held by the author/owner(s). Workshop on the Practical Use of Recommender Systems, Algorithms and Technologies (PRSAT 2010), held in conjunction with RecSys 2010. September 30, 2010, Barcelona, Spain. 1.
The idea of recommender systems is to automatically suggest items to each user that s/he may nfid appealing. The quality of recommender systems can be assessed with respect to various criteria, including accuracy, diversity, surprise / serendipity, and explainability of recommendations.
This paper is concerned with accuracy. The root mean
squared error (RMSE) has become the most popular
accuracy measure in the literature of recommender systems–for
training and testing. Its computational ecffiiency is one of
its main advantages. Impressive progress has been made in
predicting rating values with small RMSE, and it is
impossible to name all approaches, e.g., [
In this paper, we consider the top-k hit rate–based on all (unrated) items–as the natural accuracy measure for recommender systems, as only a few out of all unrated items can be recommended to a user in practice (see Section 3 for exact denfiition of top- k hit rate). While this measure is computationally tractable for testing the predictions of recommender systems, unfortunately it is computationally very costly for training recommender systems. For training, we thus resort to appropriate surrogate objective functions that are computationally ecffiient.
In recent work [
Motivated by the results of [
In our second contribution, we conduct elaborate
experiments on the Netiflx Prize data [
As the third contribution of this paper, we generalize this
probabilistic approach as to account for possibly multiple
observed rating values for a user-item pair in Section 5. This
general framework subsumes several applications in
addition to the one outlined in Section 2. Two of which are
outlined in Section 5: while the training objective function
for a recommender system concerning TV programs seems
motivated in an ad-hoc manner in [
In this section, we outline a probabilistic framework that
allows us to incorporate background knowledge when
training recommender systems on available data. The use of
background knowledge in addition to the available training
data can signicfiantly improve the accuracy of recommender
systems on performance measures like top-k hit rate, recall,
precision, area under the ROC curve, etc. This was
demonstrated for implicit feedback data in [
The Bayesian approach lends itself naturally to this task. We consider the rating matrix R as a matrix of random variables: each element Ri,u concerning item i = 1, ..., i0 and user u = 1, ..., u0 is a random variable with normal distribution, where i0 denotes the number of items and u0 is the number of users. 2.1
We take a collaborative lfitering approach, and use a lowrank matrix-factorization model, which has proven successful in many publications. Like the rating matrix, we consider our model as a matrix of random variables, M . Each random variable Mi,u corresponds to the rating of item i assigned by user u. In matrix notation, it reads
M = roffset + P Q> (1) where roffset ∈ R is an offset value, and P , Q are low-rank matrices of random variables with dimensions i0 × d0 and u0 × d0, respectively, where rank d0 i0, u0. We use upper case symbols to denote random variables (with a Gaussian distribution), and lower case symbols to denote values. 2.2
In our Bayesian approach, we rfist denfie the usual prior
over model parameters, concerning each entry of the low
rank matrices P and Q (see also [
Y Y N (Pi,d|0, σ P2,i) ) ) The vectors of variances σ Q2 = (σ Q2,u)u=1,...,u0 and σ P2 = (σ P2,i)i=1,...,i0 for all users u = 1, ..., u0 and items i = 1, ..., i0 are free parameters of the zero-mean normal prior distribution, denoted by N . There are several ways of denfiing the standard deviations in Eq. 2, eventually resulting in dieffrent kinds of regularization. The obvious choice is to assume that σ P,i = σ Q,u = 1/√2λ 0 ∀ i, u, with λ 0 ∈ R. This results in the regularization term
log p(M |σ P2 , σ Q2) = − λ 0 `||P ||22 + ||Q||22´ + c1, where || · || 2 denotes the Frobenius norm of a matrix, and c1 is an irrelevant constant when training our model.
When optimizing root mean square error on observed data
(like in the Netiflx Prize competition [
X u0(i) X Pi2,d + X i0(u) X Qu,d 2 i d u d
!1
X
X Pi2,d + Qu,d A + c2,
2
d
!
+ c2
(3)
where c2 denotes again an irrelevant constant when
training our model. Note that this choice increasingly regularizes
the model parameters related to the items and users with a
larger number of observed ratings. This may seem
counterintuitive at rfist glance. A theoretical explanation for this
empirical nfiding was recently provided in [
We now incorporate the following background knowledge
into our sequential Bayesian approach: absent rating
values tend to be lower than the observed ratings on average
(see [
Y p(Ri,u = rip,ruior|Mi,u, σ p2rior,i,u), all i all u (4) where rprior denotes the matrix of virtual data points rip,ruior, and σ p2rior the matrix with elements σ p2rior,i,u. We assume that the probabilities in this likelihood are determined by normal distributions with mean Mi,u and variance σ p2rior,i,u. The log likelihood then reads log p(rprior|M, σ p2rior) = − X X wip,ruior “rip,ruior − Mi,u ”2 +c3 all i all u (5) where we denfied the weights of the virtual data points as wip,ruior = 1/(2σ p2rior,i,u); c3 is again an irrelevant constant when training our model.
With Bayes rule, we obtain the posterior distribution of ( · ( i d u d Y Y N (Qu,d|0, σ Q2,u) (2) the model in light of these virtual data points: p(M |rprior, wprior) = p(rprior, wprior|M )p(M ) p(rprior, wprior) This equation combines our prior concerning the elements in the matrices P and Q (for regularization) with our background knowledge on the expected rating values. This serves as our prior when observing the actual rating values in the training data. p(Ri,u = rio,bus|Mi,u, σ o2bs,i,u) Again assuming a normal distribution, the log likelihood reads: log p(robs|M, σ o2bs) = − wio,bus(ri,u − Mi,u)2 + c4
obs
X where we assumed in the denominator that the observed ratings are independent of the chosen prior ratings, i.e., p(robs|rprior) = p(robs). Substituting Eqs. 3, 5, 6 and 7 into Eq. 8, we obtain the following log posterior: log p(M |robs, rprior, wobs, wprior, λ ) = ) ) (9) prior weights to be identical: wprior = wip,ruior for all (i, u). In summary, the three tuning parameters in Eq. 9 are wprior, rprior and λ , which can be chosen as to optimize the performance measure on cross-validation data. 2.6
For computational ecffiiency, our training aims to nfid
the maximum-a-posteriori (MAP) parameter estimate of our
model, i.e., the MAP estimates Pˆ and Qˆ of the matrices P
and Q. We use the alternating least squares approach. The
idea is that one matrix can be optimized exactly while the
other one is assumed xfied. A local maximum of the log
posterior can be found by alternating between the matrices
Pˆ and Qˆ. While local optima exist [
A key challenge in testing recommender systems is that
the observed ratings in the available data typically provide a
skewed picture of the (unknown) true distribution
concerning all ratings [
The top-k hit-rate / recall of relevant items is a
particularly useful performance measure for assessing the accuracy
of recommender systems [
Recall can be calculated for a user by ranking the items according to their scores predicted by the recommender system, and determining the fraction of relevant items that are among the top-k items, i.e., the k items with the highest scores. The value of k ∈ N has to be chosen, e.g., as the number of items that can be recommended to a user. Only small values of k are important in practice. The goal is to maximize recall for the chosen value of k.
Recall has two interesting properties in this context [
X ˆPi2,d + Q2u,d˜ wip,ruior (rip,ruior − Mi,u)2 + λ
X ˆPi2,d + Q2u,d˜ d d − − +c5 We found that using a prior that involves also the regularization term of P and Q in the third line in Eq. 9 leads to a slight improvements in our experimental results. The weights wio,bus and wip,ruior as well as λ 0 are absorbed in λ ; this is a slight but straight-forward generalization of the prior in Section 2.2; c5 is again an irrelevant constant for training.
Eq. 9 serves as our training objective function. For
simplicity, we choose the same value for all virtual rating values
rprior. For computational efficiency, we choose our model
osffet to equal the prior rating values: roffset = rprior. Its
main eeffct is that this retains the sparsity of the observed
rating matrix. Apart from that, it also leads to the
simplicfiation: ( rip,ruior − Mi,u)2 = ((P Q>)i,u)2. For simplicity, we
set wio,bus = 1 for all observed pairs (i, u), and also choose all
recall also has the larger precision. More interestingly,
however, recall can be calculated from the available MNAR data
and provides an unbiased estimate for the recall
concerning the (unknown) complete data (which comprises all
rating values of all users) under the following assumption: the
relevant ratings are missing at random, while an arbitrary
missing-data mechanism may apply to all other rating values
(as long as they are missing with a larger probability than
the relevant ones) [
The expected performance on the (unknown) complete data is important because it is directly related to user experience: the recommender system has to pick a few items from among all items the user has not rated yet (and which may hence be new to the user); one can expect the distribution on all unrated items to be well-approximated by the distribution on all (rated and unrated) items under the (mild) assumption that only a small fraction of the relevant ratings has been observed.
Given that ground truth (i.e., the complete data) is
typically not available (at low cost), the validity of the
assumption in [
We consider the case that the probability of observing a relevant rating depends on the popularity of items. We denfie the popularity of an item by the number Nc+omplete,i of relevant ratings it obtained in the (unknown) complete data. Let No+bs,i be the number of relevant ratings observed in the available data; then the probability of observing a relevant rating regarding item i is pobs(i) =
No+bs,i Nc+omplete,i .
In the rfist scenario, pobs may take only two values: it is
small for unpopular items, and large for popular items. If
the ratio of these two values approaches innfiity, this results
in si → 0 for popular items. Removing the relevant
ratings of the most popular items from the test set is indeed
common practice, e.g., in [
si ∝ 1/(No+bs,i)γ/ (γ +1).
Note that the unknown proportionality factor cancels out in Eq. 12, so that it provides an unbiased estimate of the recall concerning the complete data for the correct polynomial degree γ (and assuming that there are no additional factors underlying the missing data mechanism). If γ = 0, the relevant ratings are missing at random; if γ = 1, the probability of observing relevant ratings increases linearly with item popularity (Nc+omplete,i). The extreme case when γ → ∞ has several interesting properties: rfist, one observes in the available data only relevant ratings of the item with the largest popularity in the complete data, which does not agree with empirical evidence. Second, as γ/ (γ + 1) → 1, the stratification weights si are inversely proportional to the number of observed relevant ratings No+bs,i; this means that every item has the same weight in the recall-measure in Eq. 12, independent of its number of observed relevant ratings. This means that, once an item obtains its rfist relevant rating, it is weighted the same as all other items that may have obtained thousands of relevant ratings. This obviously entails low robustness against statistical noise as well as against manipulation and attacks. As this is the limiting case of γ , we nevertheless provide experimental results for this extreme scenario in Section 4. The (unknown) realistic case can be expected to be less extreme.
If the test set is a random subset of the observed data, then No+bs,i can be determined from the test set. Given that the training set is typically much larger than the test set, it might be more robust to determine No+bs,i based on all the available data, i.e., the training and test set combined. We use the latter in our experiments reported in the next section, but the former choice of No+bs,i leads to very similar results.
Finally, we denfie recall( k) = Pu wurecallu(k) as the
average recall over all users, with normalized weights, P
∝ P u wu =
1, like in [
Obviously, straticfiation like in Eq. 12 carries over
analogously to other measures, like ATOP [
This section summarizes our results on the Netflix Prize
data [
Pi∈Su+,k si recallu(k) = P i∈Su+ si on the available MNAR data; Su+ denotes the set of relevant items of user u; Su+,k is the subset of relevant items ranked in the top k items based on the predictions of the recommender system; the popularity-straticfiation weight for each item i is si =
1 pobs(i) .
If we choose the straticfiation weights si = 1 for all items i, we obtain the usual recall measure. Given that complete data is unavailable (at low cost), pobs in Eq. 11 cannot be calculated in practice. For this reason, we now examine different choices for pobs and their effects on recall in Eq. 12. 0.3 ll ca0.2 e r 0.1 user-pairs. The ratings take integer values from 1 (worst) to 5 (best). The provided data are already split into a training and a probe set. We removed the probe set from the provided data as to obtain our training set.
We consider 5-star ratings as relevant to a user (as denfied
above), and use the popularity-stratiefid recall, as outlined
in Section 3, as our performance measures on the Netiflx
probe set. For all experiments, we chose rank d = 50 of our
low-rank matrix factorization (MF) model (Eq. 1). In the
following, we compared our MF model, trained with dieffrent
prior weights wprior > 0, against the popular MF approach
that ignores the missing-data mechanism (i.e., wprior = 0 in
our notation). The latter achieved a root mean square error
on the observed ratings in the probe set of 0.922.
Additionally, we compared to singular value decomposition (SVD),
where we used the svds function of Matlab, which implicitly
imputes a zero value (with unit weight) for all missing
ratings; and to the bestseller list (BS), which ranks the items
according to the number of ratings in the training set. The
values of tuning parameters in our training objective
function (log posterior) in Eq. 9 are summarized in Figure 1 (b).
Like in [
Figures 1 (c)-(f) shows the performance of these
models under dieffrent test scenarios, concerning our
popularitystratiefid recall measure for the practically important range
of small k values. For computational ecffiiency, we computed
recall by randomly sampling, for a user, 1,000 unrated items
for each relevant rating in the test set, like in [
When comparing the dieffrent graphs, it is obvious that the performance of all the models depends on the (unknown) missing-data mechanism concerning the relevant ratings. In particular, when pobs of relevant ratings is assumed to increase more rapidly with growing popularity (Nc+omplete,i), the expected recall on the (unknown) complete data decreases for all models, cf. Figures 1 (c)→(d)→(e), and (c)→(f).
As pobs of relevant ratings increases more rapidly with
item popularity (compare Figures 1 (c)→(d)→(e)), the
difference in recall among the various MF models decreases.
Training with smaller but positive weights wprior > 0 results
in the best recall on the test set, even in the unrealistic
extreme limit in Figure 1 (e). This suggests that, compared
to the popular approach of ignoring the missing-data
mechanism when training MF models, recall can be improved
by using a small prior weight wprior > 0; its value is
upper bounded by the value that optimizes the (usual) recall
measure on the test set, i.e., under the assumption that the
relevant ratings are missing at random, like in [
The bestseller list (BS) and SVD perform surprisingly well
if relevant ratings are missing at random, see Figure 1 (c),
while the popular MF model with wprior = 0 has low recall
in comparison. This was also found in [
This section outlines a generalization of the Bayesian approach given above. In our probabilistic approach, we consider the rating matrix R as a matrix of random variables. As each entry Ri,u is a random variable (rather than a value), this naturally allows for possibly multiple values concerning each pair (i, u) in the data. This has several advantages over a matrix of values, which has typically been considered in the literature of recommender systems. After developing our generalized probabilistic framework, we outline three special cases / applications in Section 5.2.
Let the given data set be D = {ri,u,j}i,u,j, where i = 1, ..., i0 is the index concerning items, u = 1, ..., u0 is the index regarding users, and j = 1, ... is the index over possibly multiple observed ratings for the same pair (i, u). The likelihood of the model in light of i.i.d. data reads p(D|M ) =
Y p(Riu = ri,u,j|M ). i,u,j (14) Assuming again a normal distribution of the ratings (with standard deviations σ i,u,j, or equivalently weights wi,u,j = 1/(2σ i2,u,j)), the log likelihood of the model is log p(D|M ) = − X X wi,u,j(ri,u,j −
i,u j = − X X wi,u,v(v − i,u v
Mi,u)2 + c5 Mi,u)2 + c5 (15) where the second line is obtained by switching–for each pair (i, u)–from index j (over multiple ratings in the data) to the actual rating values v; the cumulative weight is wi,u,v = Pj wi,u,jIri,u,j=v, where indicator function Iri,u,j=v = 1 if ri,u,j = v and 0 otherwise.
Combining this likelihood with the same kind of prior over the model parameters as in Eq. 2, we obtain the log posterior − − of our model: log p(M |D) = − X X wi,u,v(v − i,u v
Mi,u)2 1 X
X Pi,d −
2 i 2σ P2,i d
1 X
X Qu,d + c6 (16)
2 u 2σ Q2,i d The standard deviations σ Q,i and σ P,i may be chosen as to achieve the desired variant of regularization, as discussed in Section 2.2. 5.1
For computational efficiency, we focus on optimizing the log posterior in Eq. 16. The maximum-a-posteriori (MAP) parameter estimate of our model, i.e., the MAP estimates Pˆ and Qˆ of the matrices P and Q, can be determined by alternating least squares, which alternately optimizes one matrix while the other one is assumed xfied. Using the usual necessary condition for the optimum of Eq. 16, we equate its partial derivative to zero, and obtain the following update equation for each row i of Pˆ (for fixed Qˆ):
" Pi,· = (v¯i,· − roffset )W˜ (i)Qˆ · Qˆ>W˜ (i)Qˆ + ˆ 1 2σ P2,i I #− 1 , (17) where I denotes the identity matrix, and W˜ (i) is a diagonal matrix containing the ith row of the aggregate weight matrix with elements and v¯i,u is the average rating value wi,u =
X wi,u,v, v v¯i,u =
! X vwi,u,v /wi,u.
v Analogously, the update equation for each row u of Qˆ is: " Qˆu,· = (v¯· ,u − roffset )W˜ (u)Pˆ Pˆ>W˜ (u)Pˆ + 1 2σ Q2,u
I #− 1 (18) where W˜u is the diagonal matrix containing the uth column of the aggregate weight matrix.
This derivation shows that optimizing Eq. 16 is equivalent to optimizing log p(M |D) = −
X wi,u(v¯i,u − i,u
Mi,u)2 1 X
X Pi,d −
2 i 2σ P2,i d
1 X
X Qu,d + c6 (19)
2 u 2σ Q2,i d where multiple rating values for an item-user pair are replaced by their mean value v¯i,u and their aggregate weight wi,u. 5.2
This generalized probabilistic approach subsumes several applications as special cases. In addition to the use of virtual ratings in the prior, as outlined in Section 2, we present two additional applications in the following.
The universe of all items available for recommendation may have a structure that goes beyond a aflt list. Items can often be arranged in a hierarchical manner. For instance, songs may be grouped by artist, album, genre, etc. Possibly, there are several layers of hierarchy. Now let us consider the problem that data are available where users have rated individual songs, but the task is to recommend an artist to a user. This problem arises in several situations. For instance, the recommender system may want to suggest a concert to the user, based on the data on individual songs. Another scenario is the release of a new song by an artist: this cold start problem may be overcome by recommending the new song to users who like the artist.
The ratings matrix concerning songs and users can be used to construct a rating matrix regarding artists and users by aggregating the songs of each artist. As a user may have rated several songs of an artist, we now have possibly multiple rating values for an entry in the artist-user matrix. This is exactly the problem solved by our general approach in Section 5. Our framework also shows that, for each user, the rating of each artist can be determined as the weighted average of the ratings of his/her songs, and the weight is the sum of the weights of the songs, as one may have intuitively expected. 5.2.2
IP-TV is much more interactive than traditional TV. Concerning recommender systems, it allows one to collect information on the users’ TV consumption. This can be used to learn preferences of users to TV programs, as to make accurate recommendations of TV shows. Unlike the previous applications described in this paper, we now use implicit feedback data (time spent watching a TV show) in place of the ratings. Our approach carries over immediately. This problem can be cast in our general probabilistic framework as follows: we divide the length of each show into nmax time intervals (of equal length), where nmax is the same large integer for all shows. We consider each time interval associated with a random experiment; a show hence comprises nmax repetition of the experiment. The random variable Ri,u takes the value 1 if user u watched show i for a time-interval, and 0 otherwise. So, if user u watches ni,u ∈ {1, ..., nmax} out of nmax intervals of show i, then we have ni,u observations of value 1 for the random variable Ri,u, and nmax − ni,u observations of value 0; as shown in Section 5, these multiple observations can be substituted equivalently in our least-squares objective function in Eq. 16 by the aggregate weight nmax and the averaged value of the observations: r¯i,u = ni,u/nmax, i.e., the fraction of the show watched. In addition, if a show comprises several (e.g., weekly) episodes, then we assume that the above applies to each episode; the total weight of the show is then ti,unmax, where ti,u ∈ N is the number of episodes of show i watched by user u. Then the averaged observed value is r¯i,u = Ptji=,u1 r¯i,u,j/ti,u, where r¯i,u,j is the fraction of episode j watched (analogous to above). This results in the log likelihood (with an irrelevant constant c7 in our optimization problem): log p(D|M ) = − nmax X ti,u(r¯i,u − Mi,u)2 + c7, (i,u)∈S where the set S contains all pairs (i, u) with shows i that have been watched at least partially by user u. Concerning all shows, we additionally incorporate background knowledge that users tend to not like shows with some small condfience / weight. This weight is small, as it can also be interpreted as the variance of our prior, which is large as there are many reasons for not watching a show. Analogously to Section 2, we use virtual observations of value 0, with weight wprior. This results in the log likelihood in light of the virtual data points Dprior: log p(Dprior|M ) = − nmaxwprior
X (0 − Mi,u)2 + c8 all (i,u) Combining these two likelihoods, together with the prior over the model parameters in Eq. 2 (rfist variant), we obtain the log posterior log p(M |D, Dprior) ∝ −
X ti,u(r¯i,u − Mi,u)2 − wprior X (0 − Mi,u)2 (i,u)∈S
" − λ X d
# X Pi2,d + X Qu,d + c9,
2 i u all (i,u) where we replaced the standard deviations in the prior by λ ∈ R; the proportionality is due to omitting nmax, which is an irrelevant constant when optimizing the posterior; c9 is an irrelevant additive constant in our optimization problem.
In [
# X Pi2,d + X Qu,d ,
2 i u (i,u)6∈S where nit,u = Pj ni,u,j is the total time spent by user u watching show i (including all episodes j); S denotes again the set of pairs (i, u) where show i is at least partially watched by user u; α takes essentially the role of wprior.
Interestingly, their objective function is similar to ours.
The main dieffrence, however, is in the least squares term,
where we fit our model to the fraction ¯ ri,u of the show
watched, while the indicator value 1 is used in [
Preliminary Experiments: We used a (proprietary)
IP-TV data set [
This paper provides three contributions. First, we
outlined a Bayesian framework that naturally allowed us to
insert background knowledge concerning the missing-data
mechanism underlying the observed rating data. The
obtained log posterior probability of the model is very similar
to the training objective function outlined in [
In our second contribution, we conducted experiments where we considered several hypothetical missing-data mechanisms underlying the observed real-world data. Given that the true missing-data mechanism is unknown in the absence of ground truth data, this sensitivity analysis provided valuable information. Our key insight based on these experiments is that the top-k hit-rate or recall can be improved considerably by training recommender systems with an appropriately chosen small positive prior weight concerning background knowledge on the missing-data mechanism. This is in contrast to the popular approach in the literature, which only considers observed ratings.
As third contribution, we provided a generalized probabilistic framework for factorizing a user-item-matrix that possibly has multiple observed rating values associated with each user-item pair. We discussed three important special cases / applications: besides training a top-k recommender system by using virtual data points, we outlined how ratings can be aggregated when items are grouped in a hierarchical manner rather than in a aflt list, and how recommendations can be made using this hierarchical structure. Additionally, this framework enabled us to derive the training objective function for a recommender system on sequential data concerning TV consumption. This derivation not only provides a clear understanding of the assumptions underlying the training objective function, but also led to improvements in the top-k hit rate over state-of-the-art approaches in our preliminary experiments.
We are greatly indebted to Tin Ho for her encouragement and support of this work. We are also very grateful to the anonymous reviewers for their valuable feedback. 7.