# A Recommender System Based on a Double Feature Allocation Model

Qiaohui Lin and Peter Müller

###### Abstract

A collaborative filtering recommender system predicts user preferences by
discovering common features among users and items. We implement such inference
using a Bayesian double feature allocation model, that is, a model for random
pairs of subsets. We use an Indian buffet process (IBP) to link users and
items to features. Here a feature is a subset of users and a matching subset
of items. By training feature-specific rating effects, we predict ratings. We
use MovieLens Data to demonstrate posterior inference in the model and
prediction of user preferences for unseen items compared to items they have
previously rated.

Part of the implementation is a novel semi-consensus Monte Carlo method to
accomodate large numbers of users and items, as is typical for related
applications. The proposed approach implements parallel posterior sampling in
multiple shards of users while sharing item-related global parameters across
shards.

##  1 Introduction

We develop a nonparametric Bayesian model-based approach to collaborative
filtering for random subsets of items and users. The main contributions are
the construction of a suitable prior for pairs of subsets of items and users
(features), the possibility to report coherent inference on such features, and
a consensus Monte Carlo approach to allow practical implementation of
posterior inference.

Collaborative filtering refers to recommender systems that predict
personalized user preferences for products (i.e, ratings, rankings) by
discovering similarity patterns among users and items, and make corresponding
recommendations (Sarwar et al. [2001], Schafer et al. [2007], Koren and Bell
[2015]). It has been widely adopted by e-commerce websites and online
streaming services. The Netflix Prize since 2006 has encouraged more progress
in this field.

The Netflix data (Marlin [2004], Koren and Bell [2015], Sedhain et al. [2015],
Liang et al. [2018]) is a widely used benchmark dataset used in collaborative
filtering research. It is a sparse matrix with 463435 rows for users, 17769
columns for items (movies) and 56.9 million entries based on ratings between
1999 and 2005. On average a movie has 5600 ratings and a user rates 208 movies
(Koren and Bell [2015]). Thus the density of the matrix is as low as 0.69%
(Liang et al. [2018]).

Collaborative filtering methods have progressed from naive nearest neighbor
methods, to well-adopted matrix factorization methods, to probabilistic
Bayesian models with latent factors and, more recently, to generative models.
Nearest neighbour methods (Schafer et al. [2007], Koren and Bell [2015]) are
intuitively appealing. The idea is that the rating $r_{ui}$ of user $u$ for
item $i$ is likely to be close to the ratings $r_{u^{\prime}i}$ of similar
users $u^{\prime}$ for the same item, or the ratings $r_{ui^{\prime}}$ of
similar items $i^{\prime}$ by the same user $u$. The key here is to measure
the similarity of users (or items). The simplest measure uses the correlation
coefficient. Using a similarity measure, we can identify the $k$ items that
are most similar to $i$ rated by $u$. Denote the set of such items by
$S^{k}(i;u)$. We can then predict $r_{ui}$ using a weighted average of the
ratings of items in $S^{k}(i;u)$ that are rated by user $u$, using weights
proportional to the respective similarlities, i.e, the more similar the
neighbour is, the more weight it gets.

Assume now that ratings are metric, facilitating the use of one of the most
popular methods in collaborative filtering based on matrix factorization. The
$m\times n$ rating matrix $R=[r_{ui}]$ with $m$ users and $n$ movies can be
written as the product of a $k\times m$ matrix $P$ and $k\times n$ matrix $Q$,
i.e, $R=P^{T}Q$. Here, $k$ is the number of latent factors $z$, for example,
movie genres. Then the $P$ matrix can be interpreted as representing the
preferences of users for the $k$ genres and the $Q$ matrix can be interpreted
as a classification of the movies relative to these genres. The rating
$r_{ui}$ for a specific movie and user then becomes $r_{ui}=p_{u}^{T}q_{i}$,
using $p_{u}$ and $q_{i}$ to denote the corresponding columns of $P$ and $Q$,
respectively.

The problem then is to find preferences $P$ and $Q$ to best predict observed
ratings $R$ by $\hat{R}=P^{T}Q$, while controlling the number of the latent
factors (the rank of $P$ and $Q$). Thus the objective is

$min_{Q,P}\sum_{u,i}(r_{ui}-p_{u}^{T}q_{i})^{2}+\lambda(\sum_{i}||q_{i}||^{2}+\sum_{u}||p_{u}||^{2}).$
(1)

The parameter $\lambda$ is determined by cross validation. Minimization is
performed by stochastic gradient descent. A popular implementation is reported
in Funk [2006] and Paterek [2007].

Matrix factorization captures the latent pattern in users and items. A
convenient stochastic gradient descent made it a winner of the Netflix prize.
An additional advantage is the ease of incorporating temporal dynamics.
However, the point estimation of the prediction comes without an uncertainty
measure and thus cannot serve a more complicated goal of filtering and
understanding user behavior, one of the reasons why we introduce probabilistic
models in the next section.

##  2 Probabilistic Models

Probabilistic matrix factorization [Mnih and Salakhutdinov, 2008] interprets
(1) from a probabilistic perspective. They show that minimizing the sum of
square errors while penalizing their Frobenius norm is equivalent to
maximizing the log posterior in a probabilistic model with spherical Gaussian
priors.

As before, define $R=[r_{ui}]$ as the $(m\times n)$ rating matrix, and $P$ and
$Q$ are $(k\times m)$ and $(k\times n)$ the low rank latent feature-user and
feature-item matrices. As before let $p_{u}$ and $q_{i}$ denote column $u$ and
$i$ of $P$ and $Q$, respectively, and let $N(x\mid m,V)$ denote a normal p.d.f
for random variable $x$ with moments $m$ and $V$. We assume

$p(R\mid P,Q,\sigma^{2})=\prod^{m}_{u=1}\prod^{n}_{i=1}[N(r_{ui}\mid
p_{u}^{T}q_{i},\sigma^{2})]^{I_{ui}}.$

with $I_{ui}=1$ when user $u$ has rated movie $i$ and $I_{ui}=0$ otherwise.
The model is completed with zero-mean spherical Gaussian priors on $P$ and
$Q$:

$p(P\mid\sigma_{P}^{2})=\prod^{m}_{u=1}N(p_{u}\mid 0,\sigma_{P}^{2}I),\quad
p(Q\mid\sigma_{Q}^{2})=\prod_{i=1}^{n}N(q_{i}\mid 0,\sigma_{Q}^{2}I)$

Maximizing $\log p(P,Q\mid R,\sigma^{2},\sigma_{P}^{2},\sigma_{Q}^{2})$ under
this model is equivalent to minimizing the sum of squares error with quadratic
regularization, as in

$E=\frac{1}{2}\sum_{u=1}^{m}\sum_{i=1}^{n}I_{ui}(r_{ui}-p_{u}^{T}q_{i})^{2}+\frac{\lambda_{P}}{2}\sum_{u=1}^{m}||p_{u}||^{2}+\frac{\lambda_{Q}}{2}\sum_{i=1}^{n}||q_{i}||^{2},$
(2)

where $\lambda_{P}=\frac{\sigma^{2}}{\sigma_{P}^{2}}$,
$\lambda_{Q}=\frac{\sigma^{2}}{\sigma_{Q}^{2}}$.

A minibatch gradient descent is used to find the optimal $P$ and $Q$. The
optimization defines an extension of the classic SVD model, where the modified
SVD is defined as the MAP estimate, and the classic SVD is a special case
where prior variance goes to infinity.

Mnih and Salakhutdinov [2008] also discussed other constraints to be allowed
onto the model. For all proposed models, instead of adopting a full Bayesian
approach and leading to a MCMC posterior simulation, the authors used
stochastic gradient descent. This approach limits meaningful posterior
inference for hyperparameters and uncertainty quantification of the MAP, but
on the other hand has vastly decreased the computational cost.

Some probabilistic models proposed in the recent literature build on the
Bayesian Mallows Model (Liu et al. [2019b], Liu et al. [2019a], Vitelli et al.
[2017]). Different from other models, Mallows model treats the response
variables as ordinal rankings. For example, Liu et al. [2019b] work with the
ranking on $n$ items for a user, $R_{u}=\\{R_{u1},R_{u2},\ldots,R_{un}\\}$,
$R_{ui}\in\\{1,\ldots,n\\}$ and $u=1,\ldots,m$. Mallows model is a probability
model on the space $P_{n}$ of permutations of $n$ items. A basic model uses a
latent consensus ranking $\rho\in P_{n}$ to define,

$p(R_{u}=R\mid\alpha,\rho)\propto\exp\left(-\frac{\alpha}{n}\,d(R,\rho)\right),$
(3)

where $\alpha$ is a scale parameter and
$d(R,\rho)=\sum_{i=1}^{n}|R_{i}-\rho_{i}|$ is a distance between $R$ and
$\rho$. The normalizing constant $Z_{n}(\alpha,\rho)$ in (3) is usually not
analytically tractable. Vitelli et al. [2017], for example, use instead
importance sampling and Metropolis-Hastings posterior simulation schemes.

In a more complicated scenario where users are not homogeneous, $m$ users are
arranged into $C$ clusters; each cluster has its own common consensus
$\rho_{c}$. Latent cluster membership indicator $z_{u}$ assign user $u$ to
cluster $z_{u}$. The model is then

$p(R_{1},\ldots,R_{m}\mid
z_{1},\ldots,z_{m},\alpha_{c},\rho_{c};\;c=1,\ldots,C)=\prod_{u=1}^{m}[Z_{n}(\alpha_{z_{u}})]^{-1}\exp\left\\{-\frac{\alpha_{z_{u}}}{n}d(R_{u},\rho_{z_{u}})\right\\}.$
(4)

with exponential priors on $\alpha_{c}$, uniform prior on $\rho_{c}$. The
prior for the cluster assignments $z_{u}$, $u=1,...,m$ is
$p(z_{1},\ldots,z_{u}\mid\tau_{1},\ldots,\tau_{C})=\prod_{u=1}^{m}\tau_{z_{u}}$
with a Dirichlet prior on $\tau$.

Assume now that a given user $u$ has rated only $m_{u}<m$ (instead of all $m$)
items and the objective is to make $L$ recommendations. This is equivalent to
inferring the unseen items with the $L$ highest rankings. Let $H_{u}$ denote
the top $L$ rankings different from the rankings of the observed items. To
find the top $L$ items for user $u$, we evaluate the posterior probability

$p_{iu}=p(R_{ui}\in H_{u}\mid data)$

for all unrated items $i$ (skipping details of how model (4) is modified to
allow for the observation of $m_{u}<m$ items only). The strength of Mallows
model is the use of a distribution of rankings, allowing inference beyond
point estimation. The main limitation is the need for computation-expensive
posterior MCMC simulation, which is not suitable for large data sets on
sequential updating.

Some approaches to collaborive filtering are based on LDA (latent Dirichlet
allocation) type models. For example, the User Rating Profile model (URP)
proposed in Marlin [2004] represents each user as a mixture as a user
attitudes, and the mixture proportions are distributed according to a
Dirichlet random variable. For any user $u$ we introduce a set of latent item-
specific attitudes $Z_{ui}$. Here $Z_{ui}$ is a user attitude that determines
the rating of $i$. Next, let $\beta_{viz}=P(r_{ui}=v\mid Z_{ui}=z)$ denote
item-specific rating probabilities. Like in the Latent Dirichlet Allocation
(LDA) model, $\theta_{u}$ is a Dirichlet random variable with parameter
$\alpha$, and $p(Z=z)=\theta_{uz}$.

A user profile is thus $r_{u}=[r_{u1},\ldots,r_{un}]$ with a sampling model
$p(r_{u}\mid Z_{u},\beta)$ determined by the described hierarchical model.
However, the posterior distribution for $\theta$ and $Z_{u}$ is intractable.

Another line of research are based on the use neural networks or, more
generally, generative models for collaborative filtering [Sedhain et al.,
2015, Li et al., 2016, He et al., 2017, Liang et al., 2018]. Variational
Autoencoder is one of the popular methods.

The use of variational autoencoder models originated from Kingma and Welling
[2013]. To apply it to the collaborative filtering, following Liang et al.
[2018], we assume the hidden factors are latent variables $z_{u}$, the number
of clicks of a user to all items is $x_{u}=[x_{u1},,,x_{un}]$,
$f_{\theta}(z_{u})$ is a non-linear function to produce a probability
distribution of $\pi(z_{u})$, and $x_{u}$ is multinomial with total number of
clicks $N_{u}$ and probability vector $\pi(z_{u})$.

$\displaystyle z_{u}\sim N($ $\displaystyle
0,I_{k}),\quad\pi(z_{u})\propto\exp[f_{\theta}(z_{u})],$ $\displaystyle
x_{u}\sim\mathrm{Mult}(N_{u},\pi(z_{u})),$

One common approach is to use variational inference to approximate the
intractable posterior $p(z_{u}|x_{u})$ with a variational approximate
distribution $q(z_{u})$. Full algorithm can be seen in Liang et al. [2018] and
Kingma and Welling [2013].

Autorec model (Sedhain et al. [2015]) is another version of this model applied
to rating instead of clicking numbers. Generative Models with variational
inference is an efficient scalable inference well suited when the dataset is
large in size and new users’ data keep flowing in. The problem is the lack of
explanibity, the absence of uncertainty quantification, and the difficulty to
fit the model to a more complicated goal, such as a specific criterion for
recommendation diversity and accuracy tradeoff.

Below we introduce a model that builds on these approaches, aiming to (i)
include learning on underlying structure that determines user preferences, and
(ii) still allows (approximate) full posterior inference. The underlying
structure that can be discovered by the proposed model is an extension of the
clusters that feature in Mallows model by adding matching subsets of items.

##  3 Double Feature Allocation Model

We introduce an alternative generative model for user ratings $r_{ui}$ of
users $u=1,\ldots,m$ for items $i=1,\ldots,n$, using notation as before. We
assume ordinal ratings $r_{ui}\in\\{1,\ldots,5\\}$. The inference goal is to
predict user preferences for un-rated movies by discovering similarity
patterns among users and movies.

The proposed model construction is guided by symmetry assumptions with respect
to items and users. That is, the probability model should be invariant with
respect to arbitrary permutation of user indices, and/or permutations of movie
indices. Models with such structure are also known as separately exchangeable,
and the rating matrix can be characterized as

$\displaystyle r_{1:m,1:n}\overset{d}{=}r_{\pi_{1}(1:m),\pi_{2}(1:n)}$ (5)

for separate permutations $\pi_{1}$ and $\pi_{2}$ of rows and columns,
respectively. Here $X\overset{d}{=}Y$ indicates equality in distribution for
two random variables.

A double feature allocation model is a model for random pairs of subsets,
first proposed in Ni et al. [2019b]. Our model is similar but different from
the Coupled Indian Buffet Process Model proposed by Chatzis [2012]. Chatzis
[2012] used two independent Indian Buffet process (IBP) for two separate
feature allocations of users and items. While in our case, a feature is a
subset of users together with a matching subset of items. We use an Indian
buffet process (IBP) prior to link users and items to features. Suppose there
are $K$ features. We use an $(m\times K)$ binary matrix $A$ to link users to
features, with $A_{uk}=1$ indicating that user $u$ is in feature $k$. Another
$(n\times K)$ binary matrix $B$ links items to features, with $B_{ik}=1$
meaning that movie $i$ is in feature $k$. Figure 1 is a stylized
representation of the double feature allocation model. Same colored block
indicate one feature $k$, a subset of users and a matching subset of items.

| Item 1 | Item 2 | Item 3 | Item 4 | Item 5 | … | Item n-1 | Item n
---|---|---|---|---|---|---|---|---
User 1 |  |  |  |  |  |  |  |
User 2 |  |  |  |  |  |  |  |
User 3 |  |  |  |  |  |  |  |
User 4 |  |  |  |  |  |  |  |
User 5 |  |  |  |  |  |  |  |
… |  |  |  |  |  |  |  |
User m-1 |  |  |  |  |  |  |  |
User m |  |  |  |  |  |  |  |
Figure 1: Illustration of features in Double Feature Allocation Model, with
colored box as features.

The IBP prior on $A$ [Griffiths and Ghahramani, 2011] is defined as follows.
The model includes an unknown number $K$ of features, and can be written as

$p(A)=\frac{\lambda^{K}\text{exp}(-\lambda
H)}{K!}\prod_{k=1}^{K}\frac{\Gamma(m_{k})\Gamma(m-m_{k}+1)}{\Gamma(m+1)},$

where $p(A)$ is without order of columns, $\lambda$ is a fixed hyperparameter
(concentration parameter), $H$ is the harmonic number $H=\sum_{u=1}^{m}1/u$,
$m_{k}$ is the sum of column $k$, $m_{k}=\sum_{u=1}^{m}A_{uk}$. The number of
features $K$ is random and unbounded, and features are exchangeable a priori.
. The IBP is easiest described as a generative model building up $A$ row by
row, starting with $u=1$, and adding columns of $A$, i.e., features, and
indexing features by appearance. Let $K_{u}$ denote the number of features
that are introduced after the first $u$ users, starting with $K_{0}=0$, and
let $m_{u,k}=\sum_{v=1}^{u}A_{vi}$ denote the cardinality of feature $k$ among
the first $u$ users, $k=1,\ldots,K_{u}$. Considering the respective next user
$u$ we then proceed as follows. First we decide inclusion into one of the
existing features, $k=1,\ldots,K_{u-1}$, with probability $p(A_{uk}=1\mid
A_{1\cdots u-1,1\cdots K_{u-1}})=m_{u-1,k}/u$ Then we add a Poisson random
number $K^{+}_{u}\sim\mbox{Poi}(\alpha/u)$ new features with $A_{uk}=1$,
$k=K_{u-1},\ldots,K_{u-1}+K^{+}_{u}$ and increment $K_{u}=K_{u-1}+K^{+}_{u}$.
Implicit in the construction is a constraint of all zeroes in the right upper
corner of $A$, i.e., $A_{vk}=0$, $k>K_{u}$ and $v<u$. We remove the constraint
of indexing items by appearance by using a final step of randomly permuting
the final $K=K_{n}$ columns. In practice, for the purpose of prior sampling,
we only need the conditional probability of $A_{uk}=1$,

$p(A_{uk}=1|A_{-u,k})=m_{-u,k}/m,$

where $A_{-u,k}$ is the $k$th column of $A$ excluding $A_{uk}$ and $m_{-u,k}$
is the sum of column $k$ of $A$ excluding $A_{uk}$.

Given $A$, the item-feature matrix $B$ inherits $K$ features from $A$. For
simplicity, we assume independent Bernoullis:

$p(B_{ik}=1|A)=p,$

with the prior parameter $p$ usually chosen to be a small number to control
the number of features a movie can be in. Such parsimony ensures that features
do not share too many common movies and preserve their differences. In
posterior sampling, $p$ and the dimension of $B$ will be updated each time
after we update $A$.

Note that the IBP $p(A)$ together with $p(B\mid A)$ define a joint model
$p(A,B)$, and therefore also imply a marginal $p(B)$. The current model does
not imply a marginal IBP prior for $p(B)$. If a more symmetric construction
with an IBP marginal prior on $B$ were desired, it could be easily achieved.
Let $p_{\text{IBP}}(B)$ denote an IBP prior on a random binary matrix $B$,
including the number of columns, $K_{B}$, and let $p_{\text{IBP}}(B\mid
K_{B})$ denote the conditional distribution of $B$, conditional on the number
of columns equal to $K_{B}$ under the IBP. Then using $p(B\mid A)=p(B\mid
K_{B}=K_{A})$ would by construction deliver $p(B)=\text{IBP}$, marginally.

Finally, we complete the inference model with a sampling model for the
observed ratings. The model links features to ratings by introducing
probabilities for $r_{ui}$ conditional on the currently imputed features of
which user $u$ and items $i$ are members.

We consider a baseline $b_{0}=2.5$ (between the extremes 1 and 5). Each
feature $k$ of which user $u$ and movie $i$ are a member adds an adjustment
$\theta_{k}$ to this baseline. For example, feature $k$ might be a pair of
subsets of comedy movies and comedy fans. Then $\theta_{k}$ would be a
positive increment from baseline. Similarly, if feature $k$ is a pair of
subsets of comedy movies and comedy haters, then $\theta_{k}$ should be a
negative adjustment. We allow an item and a user to be in multiple subsets
with features having aggregative influences.

We also include global parameters $\rho_{i}$ for each movie regardless of
feature allocation. The parameter $\rho_{i}$ has an interpretation as overall
mean rating for movie $i$ (on the ordinal probit scale). This movie-specific
offset reflects if a movie is generally popular and well-received among
audiences or vice versa. We complete the prior model with independent priors
for $\theta_{k}$ and $\rho_{i}$, assuming $\theta_{k}\sim h_{\theta}$ and
$\rho_{i}\sim h_{\rho}$, i.i.d.

The sampling model is then defined as an ordinal probit model including the
described feature-specific and item-specific parameters. Denote by $Z_{ui}$ a
latent continuous probit score for the rating of user $u$ for movie $i$,
denote by $K^{*}_{ui}=\\{k:\;A_{uk}=B_{ki}=1\\}$ the set of features that
include both, user $u$ and movie $i$. We use an inverse-gamma hyperprior on
$\tau^{2}$, and normal priors for both $\theta$ and $\rho$. We assume

$Z_{ui}\mid A,B,\theta,\tau\sim N(b_{0}+\sum_{k\in
K^{*}}\theta_{k}+\rho_{i},\tau^{2}).$

The probit scores are linked with the data in the usual ordinal probit model
as

$\displaystyle r_{ui}=\left\\{\begin{array}[]{cc}1,&Z_{ui}\leq 1\\\
x,x\in\\{2,3,4\\},&x-1<Z_{ui}\leq x\\\ 5,&Z_{ui}>4\end{array}\right..$ (9)

##  4 Posterior Inference

###  4.1 Posterior Sampling Algorithm

We implement posterior inference using MCMC posterior simulation. Let
$\omega=(A,B,\theta,\tau,\rho)$ denote the currently computed parameters. See
the appendix for transition probabilities to update $B$, $\theta$, $\rho$ and
$\tau$. Only the transition probability for $A$ requires more discussion.

Denote the number of iteration as superscript (t), denote the $u$th row of $A$
and $R$ as $A_{u\cdot}$ and $R_{u\cdot}$, respectively, and denote the $i$th
column of $R$ as $R_{\cdot i}$. The following three steps define a reversible
jump transition probability for $A_{u\cdot}$. Below, let
$(\tilde{A},\tilde{\theta})$ denote proposed new values for $A,\theta$, and
let $\tilde{\omega}=(\tilde{A},B,\tilde{\theta},\tau,\rho)$. The transition
probability to update $A_{u\cdot}$ is defined as follows.

  1. 1.

For all the $k$ with $m_{-u,k}\neq 0$, update
$p(A_{uk}=x\mid\cdot)\propto\frac{m_{-u,k}}{m}p(R_{u\cdot}\mid
A_{uk}=x,A_{-u,k},\theta,\rho,\tau,B)$, $x=\\{0,1\\}$.

  2. 2.

Reversible jump proposal. We refer to all features with $m_{-u,k}=0$ as
singular features. W.l.o.g. assume $k=1,\ldots,K_{0}$ are not singular, and
$k=K_{0}+1,\ldots,K$ are the singular features. We create a proposal by first
dropping all singular features, i.e., retaining in $\tilde{A}$ only the first
$K_{0}$ columns of $A$, proposing $\tilde{A}=A[\;\cdot\;,(1,\ldots,K_{0})]$.
Next we propose $\tilde{K}^{+}_{u}\sim Pois(\lambda/n)$ new (singular)
features, together with (new) feature-specific parameters
$\tilde{\theta}_{k}\sim h_{\theta}$,
$k=K_{0}+1,\ldots,K_{0}+\tilde{K}^{+}_{u}$. We add the new features to
$\tilde{A}$ with $\tilde{A}_{uk}=1$ and $\tilde{A}_{vk}=0$ for
$k=K_{0}+1,\ldots,K_{0}+\tilde{K}^{+}_{u}$ and $v\neq u$.

  3. 3.

Metropolis-Hastings acceptance probability. We denote with
$S=\\{K_{0}+1,\ldots,K\\}$ the indices of the singular features, and find

$\alpha=\frac{p(R_{u\cdot}\mid\tilde{\omega})}{p(R_{u\cdot}\mid\omega)}\cdot\frac{p(\tilde{\theta}_{S})\mbox{Poi}(K^{+}_{u}\mid\frac{\lambda}{n})}{p(\theta_{S})\,\mbox{Poi}(\tilde{K}^{+}_{u}\mid\frac{\lambda}{n})}\cdot\frac{p(\theta_{S})\mbox{Poi}(\tilde{K}^{+}_{u}\mid\frac{\lambda}{n})}{p(\tilde{\theta}_{S})\,\mbox{Poi}(K^{+}_{u}\mid\frac{\lambda}{n})}=\frac{p(R_{u\cdot}\mid\tilde{\omega})}{p(R_{u\cdot}\mid\omega)}$

With probability $\min(1,\alpha)$, we accept the proposal
$\tilde{A},\tilde{\theta}$ and set $A=\tilde{A}$, $\theta=\tilde{\theta}$.
Otherwise we keep $A,\theta$ unchanged.

We evaluate fitted mean ratings as $\bar{r}_{ui}=\int
p(r_{ui}^{mis}=x|\omega)p(\omega|R^{obs})d\omega\approx\frac{1}{T}\sum_{t=1}^{T}p(r_{ui}^{mis}=x|\omega^{t})$
and predict unseen ratings for items $i$ with $I_{ui}=0$ by maximizing

$\hat{r}_{ui}=\operatorname*{\arg\\!\max}_{x\in\\{1,2,3,4,5\\}}\bar{r}_{ui}$
(10)

###  4.2 A Consensus Monte Carlo Method for Large Number of Users

The described posterior simulation can be computationally costly. Large
recommender systems usually have at least thousands, or millions of users,
render the usual MCMC impractical. Thus we use the idea of Consensus Monte
Carlo [Scott et al., 2016], to split data into shards, run MCMC one each shard
in parallel on different machines, and then reconcile posterior inference from
the shards into a reconstruction of posterior inference under the full data.
Ni et al. [2019a] applied Consensus Monte Carlo to Bayesian nonparametric
models on clustering and feature allocation. Their method relies on shared
data points (anchor points) across shards, and they merge random subsets
(clusters or features) on different machines based on the number of common
anchor points in the two sets. This strategy, though appealing, does not work
in our case.

Using a set of users as common anchors to define a criterion for merging
features would ignore the possibility of different sets of movies being paired
with these users in different shards. And similarly for using sets of movies
only. A practicable implementation would need to use a criterion based on
shared users, shared movies and similar imputed feature-specific effects
$\theta_{k}$. A related criterion to merge subsets would require several ad-
hoc choices and tuning parameters. In simulations we found that the involved
approximations left the joint posterior reconstruction of little practical
value. Instead we propose an alternative strategy based on Consensus Monte
Carlo for global parameters, but keeping inference for random subsets local to
each shard. We refer to this strategy as ”semi-local Consensus Monte Carlo”.

We use shards that split the data by subsets of users. Let $s=1,\ldots,S$
index the shards, and let $\bigcup_{s=1}^{S}U_{i}=\\{1,\ldots,m\\}$ denote the
split of users into the $S$ shards. Also, let $R_{s}$ denote the data for the
users in shard $s$. In our model, the random features that are imputed under
posterior inference under shard $s$ naturally include only users from $U_{s}$,
making the features local parameters, while the only global parameters are
$\rho_{i},i=1,\ldots,n$. We use a CMC approximation of the joint posterior for
the global parameter $\rho$ as

$p(\rho\mid R)\approx\prod_{s=1}^{s}p(R_{s}\mid\rho)\,p(\rho)^{1/S}.$

The nature of the approximation is to assume independence of the marginal
distribution of global parameters across shards. Conditional on a posterior
sample $\tilde{\rho}\sim p(\rho\mid R)$ we then use shard-specific posterior
samples of the shard-specific parameters. The latter is implemented by
selecting stored posterior Monte Carlo samples $(A,B,\theta,\rho,\tau)$ with
$|\rho-\tilde{\rho}|<\epsilon$. Here $A$ refers to users in shard $s$ only,
and $B$ and $\theta$ are linked with the subsets that are represented by $A$,
making $(A,B,\theta)$ parameters that are local to each shard only.

There remains the step of creating a Monte Carlo sample for $\rho\sim
p(\rho\mid R)$. For this we use Consensus Monte Carlo for approximate normal
posterior distributions. Let $(\mu_{s},\sigma_{s})$ denote posterior mean and
standard deviation of $\rho_{i}$ in shard $s$ (for a movie $i$). We
approximate $p(\rho\mid R)\approx N(\mu,\sigma)$ with
$1/\sigma^{2}=1/\sigma^{2}_{0}+\sum_{s}1/\sigma^{2}_{s}$ and
$\mu=(\mu_{0}/\sigma_{0}^{2}+\sum_{s}\mu_{s}/\sigma^{2}_{s})/(1/\sigma^{2}_{0}+\sum_{s}1/\sigma^{2}_{s})$,
where $(\mu_{0},\sigma^{2}_{0})$ are the prior moments for $\rho$. We call
this strategy semi-local Consensus Monte Carlo.

The full algorithm is now, run MCMC described in Section 4.1 on each separate
shards, store $A$, $B$, $\theta$, $\tau$, $\rho$, every 5 iterations. Merging
the shard-specific posterior distributions for $\rho$ across shards as
described, we get an approximate global posterior from which we then generate
a posterior Monte Carlo draw $\tilde{\rho}$. In each shard, for each stored
iteration $A$, $B$, $\theta$, $\tau$, resample $\tilde{\rho}$ from the global
distributions, filter the iterations where the stored shard-posterior $\rho$
is close to resampled global $\tilde{\rho}$ and make predictions based on $A$,
$B$, $\theta$, $\tau$ and $\tilde{\rho}$ in these iterations. Note that the
filtering step within the shard is to ensure the closeness of shard posterior
MCMC and global posterior MCMC and to not totally lose the dependence of
$\tilde{\rho}$ and $(A,B,\theta,\tau)$ in the iterations we use for
prediction.

Consensus Monte Carlo Algorithm

  1. 1.

Separate users into S shards, keep entire list of movies.

  2. 2.

In each shard, carry out MCMC simulation for $A$, $B$, $\mathbf{\theta}$,
$\mathbf{\tau}$, $\mathbf{\rho}$ according to the transition probability in
Section 4.1 and appendix. Store after thinning.

  3. 3.

Merge the shard-specific posterior distributions for $\rho$, and derive an
approximate global posterior distribution for $\mathbf{\rho}$ by aggregating
shard precision.

  4. 4.

In each shard, for each stored iteration, resampling $\mathbf{\tilde{\rho}}$
from the (approximate) global posterior, keep those iterations for which
$|\rho-\tilde{\rho}|<\epsilon$, perform prediction using $A$, $B$,
$\mathbf{\theta}$, $\mathbf{\tau}$ and resampled $\tilde{\rho}$ in those
iterations.

##  5 Simulation

We implement the proposed scheme in R. The code is avaiblable in author’s
github. We use this to set up the simulations. In the simulation study, we
generate a $100\times 150$ rating matrix of $m=100$ users and $n=150$ movies
from our model. We first simulate $A$ under the IBP prior with hyperparameter
$\lambda=3$. The (random) number of columns of $A$ determines the number of
features $K$. Next we simulate a $n\times K$ binary matrix $B$ with $p$, using
independent Bernoulli draws with success probability $0.2$. For user $u$ and
movie $i$, identify the subset of features $K_{ui}^{*}={k:A_{uk}=B_{ki}=1}$.
With $b_{0}=2.5$, $\theta_{1:K}\sim N(0,2)$ and $\tau=0.25$, the latent probit
score for rating, $Z_{ui}$ is simulated from

$\displaystyle Z_{ui}|A,B,\theta,\tau\sim N(b_{0}+\sum_{k\in
K^{*}_{ui}}\theta_{k},\tau^{2}).$

Note here we did not introduce global parameter $\rho_{i}$ for movies as the
simulation data size is small and we are not using Consensus Monte Carlo to
split users into shards. Thus a global parameter for movies is optional but
not necessary in this example. The rating is generated from the probit score
as in Eq 9.

We randomly split the simulated data into $80\%$ and $20\%$ for training and
testing. We use the MCMC algorithm proposed in Section 4.1, implementing
$10000$ iterations conditional on the training data. We evaluate an estimated
user-feature relationships $\hat{A}$ following Ni et al. [2019b]. We first
calculate the maximum a posteriori (MAP) estimate $\hat{K}$ from the marginal
posterior distribution of K. Conditional on $\hat{K}$, we follow Dahl [2006]
and compute a point estimate $\hat{A}$ as

$\hat{A}=\operatorname*{\arg\\!\min}_{A^{\prime}}\int
d(A,A^{\prime})dp(A|Z,R,\hat{K}),$

where $d=\min_{\pi}H(A,\pi(A^{\prime}))$, denotes the minimum Hamming distance
between binary matrices $A$ and $A^{\prime}$ over greedy searches of
permutations of $\pi(A^{\prime})$. In Figure 2, we show the estimated user-
feature relationships $\hat{A}$, and conditional on $\hat{A}$, the point
estimates for movie-feature relationships $\hat{B}$ and the rating adjustment
for each feature $\hat{\theta}$, versus the true values $A$, $B$, $\theta$
used in simulation.

Figure 2: Simulation truth (top row) and posterior estimates (bottom row) for
$A$, $B$. Total feature (column) numbers are the same in the truth and
estimates, but features (columns) can be permuted.

We predict rating for training and testing data by maximizing posterior
predictive probability in Eq 10. The predicted rating for the training data
set is $69.12\%$ correct, meaning that in $69.12\%$ of the cases the predicted
rating exactly matches the recorded data. The same for the test data was
$66.89\%$. . We compare to inference under the matrix factorization method as
described in Section 2. We use the R package recosystem. Rank $k$ and sparsity
parameters in Eq 2 are tuned by cross validation. After training, we find a
test RMSE of 0.67 and prediction accuracy for the test data set of $60.80\%$.

Figure 3 shows the boxplots of users’ prediction accuracy for our methods and
Matrix Factorization. The left boxplot shows all users’ accuracy of predicting
ratings (level 1-5) of unseen movies. The right boxplot shows all users’
accuracy of predicting the top 10 unseen movies. Our method has a better
simulation results in both boxplots.

Figure 3: Boxplots of prediction accuracy for our method (Bayesian Double
Feature Allocation, BDFA) to Matrix Factorization (MF). Left shows boxplot of
all users’ accuracy of predicting ratings (level 1-5) of unseen movies. Right
shows boxplots of all users’ accuracy of predicting the top 10 unseen movies.

##  6 Preference Prediction with the MovieLens Data

We use the Movielens dataset 111https://grouplens.org/datasets/ with movie
ratings from 6040 users. We clean the data as in Vitelli et al. [2017], to
keep results comparable. We keep the $200$ most rated movies and users who
rated more than three movies, which yields to a $6040\times 200$ rating matrix
$R$. Our goal include traditional goals for recommender systems, such as
predicting individual user’s ratings to unseen movies, find the top-rated
movies across all users, and a pairwise preference prediction goal similar to
Vitelli et al. [2017]. The latter goal predicts the preference of unseen
movies to rated movies for each user. Here preference is defined as the
following. Denote the predicted rating of user $u$ for movie $i$ as
$\hat{R}_{ui}$, estimated as described in Section 3. We say that user $u$
strictly prefers movie $i$ to $i^{\prime}$, and write $i\succ i^{\prime}$ for
user $u$, if $\hat{R}_{ui}>\hat{R}_{ui^{\prime}}$.

In the $6040\times 200$ rating matrix, only $24.7\%$ entries are observed. We
conduct a pairwise preference prediction test using our model. For each user,
we randomly select one movie that he/she has rated as test, and train the
model based on the remaining data for that user. We compare the user’s
preference of this test movie to other (rated) movies based on the predicted
rating for the test movie to the observed ratings in the training set. This
pairwise preference test is also comparable to Vitelli et al. [2017].

We apply a double IBP prior described in Section 3 initialized with
$\lambda=3$, base line $\theta=2.5$, non-baseline $\theta$s drawn from
$N(0,\sigma_{0}^{2})$, $\sigma_{0}=2$, $\tau$ drawn from an inverse-gamma
distribution with location and scale parameters $(5,1)$. We initialize the
number of features $K$ from matrix factorization algorithm result.

Using the concensus Monte Carlo algorithm, we divided the $6040$ users into
$15$ shards and implement MCMC posterior simulation for each shard separately
in parallel. Each shard with $400$ users and $200$ movies stabilizes at around
$39$ to $43$ features. As described in Section 4, we merge inference across
shards by defining the approximate global posterior for $\rho$. For
simplicity, we set the prior precision for $\rho$ to $0$, i.e.,
$\sigma_{0}=\infty$ and simplify the merge step to evaluating
$\mu=\bar{\mu}_{s}$ as an unweighted average across shards. The after-merge
prediction follows from re-sampling $\rho$ from global posterior and
previously stored MCMC draws of other variables.

The resulting average pairwise preference accuracy in each shard is $79.1\%$
with standard deviation $0.012$. Vitelli et al. [2017] report an accuracy of
$79.6\%$ for the pairwise comparison on the same data. Figure 4 (left panel)
shows the before-merge individual shards MCMC pairwise predication accuracy
and after-merge pairwise predication accuracy. After-merge prediction accuracy
has a higher average and a smaller standard deviation among shards.

Figure 4: Left: Prediction accuracy pf pairwise comparison in each shard
before and after merge. Right: Accuracy of Level (1-5) prediction of unseen
movies within one level distance in each shard after merge and their $95\%$
confidence intervals. Red line indicates the mean over shards.

For exact predicted rating for unseen movies, on the test set within one star
distance, our model has provided $88.5\%$ accuracy. That is saying, we have on
average $88.5\%$ probability that our predicted rating is at or within in one
star difference of the actual rating. Figure 4 (right panel) shows in each
shard the prediction accuracy and its $95\%$ Confidence Interval using
binomial distribution variance and a normal approximation CI,
$(\hat{p}-\mathit{z}_{\alpha}\sqrt{\frac{\hat{p}(1-\hat{p})}{n_{s}}},\hat{p}+\mathit{z}_{\alpha}\sqrt{\frac{\hat{p}(1-\hat{p})}{n_{s}})}$,
where $\hat{p}$ denotes the shard accuracy, $\mathit{z}_{\alpha}$ denotes
$1-\frac{\alpha}{2}$ quantile of a standard normal distribution,
$\alpha=0.05$.

We also present individual predicted rating for users. Figure 5 shows three
randomly selected user’s predicted ratings, $95\%$ credible intervals for all
movies (seen and unseen), sorted by predicted ratings, versus true seen
ratings of train and test movies. Figure 6 shows in one randomly picked shard,
all users’ predicted expectation of probit scores for test movies’ ratings
versus observed test movies’ ratings.

Figure 5: Three randomly selected users’ predicted expectation of probit
scores for ratings of all movies (seen and unseen), sorted from high to low.
The red line indicates the $95\%$ credible intervals. The blue dots are true
ratings in training and green dot is the seen test movie rating. Figure 6: All
users’ predicted expectation of probit scores for ratings of test movies
versus observed ratings of test movies in one randomly picked shard.

For general rating of all movies, we show the top 10 rated movies across all
shards, which is an analysis not available by Mallows Model in Vitelli et al.
[2017]. Table 1 shows top $10$ movies with highest posterior mean of merged
$\rho$.

| Movie Name | Genre
---|---|---
1 | The Godfather (1972) | Action/Crime/Drama
2 | Schindler’s List (1993) | Drama/War
3 | Star Wars Series | Action/Adventure/Drama/Sci-Fi/War
4 | American Beauty (1999) | Comedy/Drama
5 | The Usual Suspects (1995) | Crime/Thriller
6 | Casablanca (1942) | Drama/Romance/War
7 | Raiders of the Lost Ark (1981) | Action/Adventure
8 | The Shawshank Redemption (1994) | Drama
9 | Pulp Fiction (1994) | Crime/Drama
10 | Rear Window (1954) | Mystery/Thriller
Table 1: Top $10$ highest rated movies across shards, selected by highest
posterior mean of merged $\rho$

Finally, we discuss the computational complexity reduction from our Consensus
Monte Carlo method and strategy to choose shard number $S$. Ghahramani and
Griffiths [2006] has established that an Indian Buffet Process with $N$ data
points and linear-Gaussian likelihood model has at least $O(N^{3})$
computation complexity per iteration. To be specific, Doshi-Velez and
Ghahramani [2009] pointed out that in each iteration for $N$ number of
$D$-dimensional data points and a linear-Gaussian likelihood model, given an
$N\times K$ feature assignment matrix in that iteration, the collapsed Gibbs
sampler has complexity $O(N^{3}(K^{2}+KD))$. In our double feature allocation
model with $m\times K$ matrix $A$ and $n\times K$ matrix $B$ ($m>n$) and a
probit likelihood model will at least have complexity $O(m^{3})$ each
iteration. If we split the entire data set into $S$ shards, then for each
shard we will have complexity of $O((m/S)^{3})$ each iteration, (in total
$O(S(m/S)^{3})$ for all shards). With $m=6000$ in Movielens data, we see in
Figure 7 the complexity per iteration decreases and eventually stabilizes with
$S$. We also note that more shards would result in fewer users in a shard and
fewer ratings observed in each feature $k$ in the shard and thus a larger
standard error for $\theta_{k}$. Consider the normal distribution of $\theta$,
heuristically the standard error of $\theta_{k}$ would be proportional to
$\frac{1}{n_{ks}}$ where $n_{ks}$ is the number of ratings observed in feature
$k$ in the shard $s$. Following Korwar et al. [1973], the approximate number
of features in shard $s$ is $\log(m_{s}n)$. Suppose $n_{ks}$ is proportional
to number of users and movies in the shard and inversely proportional to
number of features in shard $s$, i.e., $n_{ks}\propto
m_{s}n(\log(m_{s}n))^{-1}$, then standard of error of $\theta_{k}$ is
proportional to $((m_{s}\times n)/\log(m_{s}\times n))^{-1/2}$. Figure 7 shows
this trade-off of standard error of estimating $\theta_{k}$ increasing with
number of shards $S$ while the computation cost per iteration decreasing with
$S$. Figure 7 can serve as a eyeball guideline to pick suitable $S$.

Figure 7: Complexity per iteration and standard error of $\theta$ with the
number of shards in consensus monte carlo, $\propto$ denotes proportional up
to a constant.

##  7 Conclusion

We introduced a novel model-based approach to collaborative filtering, with
the main features being full posterior inference with interpretable parameters
and structure. Based on a generative model, our inference includes a full
probabilistic description of any desired summary. For example, the same
approach could be used for the performance of learners over problems in a
large on-line course. Inference would allow to identify the subsets of similar
learners and courses.

Limitations in the current implementation is the very approximate nature of
the reconciliation of the shard-specific posterior distribution into a
reconstructed joint posterior distribution. In particular, there is no
borrowing of information about random subsets across shards. Such features
could be added using, for example, common anchors in the Consensus Monte Carlo
method. An important limitation is the lack of using any covariate
information. For example, movies, or generally any items, have known
characteristics like actors, length, origin, year etc. Such information could
be used to include a rudimentary regression in the subset selection of the
generative model.

Interesting applications arise in many other fields beyond marketing. For
example, users could be HIV patients, items could be medications and outcomes
could any ordinally reported health outcomes, for example mental health
outcomes. Including important baseline covariates and classes of medications
could then allow to recommend suitable treatment combinations for future
patients.

## Appendix

##  Appendix A Transition probabilities to update $B$,$\theta$,$\rho$ and
$\tau$.

  * •

Update $B|R,A,\mathbf{\theta},\mathbf{\rho},\tau$,
$p(B_{ik}=x|\cdot)\propto p(B_{ik}=x)p(r_{\cdot
i}|B_{ik}=x,B_{-i,k},\theta,\rho,\tau,A)$, $x=\\{0,1\\}$

  * •

Sample auxiliary variable
$\mathbf{Z}|R,A,B,\mathbf{\theta},\mathbf{\rho},\tau$

$p(Z_{ui}|\cdot)\sim\mathrm{TruncatedNormal}(b_{0}+\sum_{k\in
K^{*}}\theta_{k}+\rho_{i},\tau^{2})$

with lower bound $r_{ui}-1$ ($r_{ui}>0$), upper bound $r_{ui}$ ($r_{ui}\leq
5$), $K^{*}$ is the set of $k$ where $A_{uk}=B_{ik}=1$.

  * •

Update $\mathbf{\tau}|\mathbf{Z},A,B,\mathbf{\theta},\mathbf{\rho}$ from
conjuagate Inverse-Gamma families.

  * •

Update $\mathbf{\theta}|\mathbf{\rho},\mathbf{Z},\tau,A,B$ and
$\mathbf{\rho}|\mathbf{\theta},\mathbf{Z},\tau,A,B$ from conjuagate Normal
families.

## References

  * Chatzis [2012] Sotirios P Chatzis.  A coupled indian buffet process model for collaborative filtering.  In _Asian Conference on Machine Learning_ , pages 65–79. PMLR, 2012\.
  * Dahl [2006] David B Dahl.  Model-based clustering for expression data via a dirichlet process mixture model.  _Bayesian inference for gene expression and proteomics_ , 4:201–218, 2006.
  * Doshi-Velez and Ghahramani [2009] Finale Doshi-Velez and Zoubin Ghahramani.  Accelerated sampling for the indian buffet process.  In _Proceedings of the 26th annual international conference on machine learning_ , pages 273–280, 2009.
  * Funk [2006] S. Funk.  Netflix update: Try this at home.  http://sifter.org/˜simon/journal/20061211.html, 2006.
  * Ghahramani and Griffiths [2006] Zoubin Ghahramani and Thomas L Griffiths.  Infinite latent feature models and the indian buffet process.  In _Advances in neural information processing systems_ , pages 475–482, 2006.
  * Griffiths and Ghahramani [2011] Thomas L Griffiths and Zoubin Ghahramani.  The indian buffet process: An introduction and review.  _Journal of Machine Learning Research_ , 12(Apr):1185–1224, 2011.
  * He et al. [2017] X. He, L. Liao, H. Zhang, L. Nie, Hu., and T. Chua.  Neural collaborative filtering.  In _Proceedings of the 26th international conference on world wide web_ , pages 173–182. International World Wide Web Conferences Steering Committee, 2017.
  * Kingma and Welling [2013] Diederik P Kingma and Max Welling.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Koren and Bell [2015] Y. Koren and R. Bell.  Advances in collaborative filtering.  In _Recommender systems handbook_ , pages 77–118. Springer, 2015\.
  * Korwar et al. [1973] Ramesh M Korwar, Myles Hollander, et al.  Contributions to the theory of dirichlet processes.  _The Annals of Probability_ , 1(4):705–711, 1973\.
  * Li et al. [2016] S. Li, A. Karatzoglou, and C. Gentile.  Collaborative filtering bandits.  In _Proceedings of the 39th International ACM SIGIR conference on Research and Development in Information Retrieval_ , pages 539–548. ACM, 2016\.
  * Liang et al. [2018] D. Liang, R. G. Krishnan, M. D. Hoffman, and T. Jebara.  Variational autoencoders for collaborative filtering.  In _Proceedings of the 2018 World Wide Web Conference_ , pages 689–698. International World Wide Web Conferences Steering Committee, 2018.
  * Liu et al. [2019a] Q. Liu, M. Crispino, I. Scheel, V. Vitelli, and A. Frigessi.  Model-based learning from preference data.  _Annual review of statistics and its application_ , 6,329-354, 2019a.
  * Liu et al. [2019b] Q. Liu, A. H. Reiner, A. Frigessi, and I. Scheel.  Diverse personalized recommendations with uncertainty from implicit preference data with the bayesian mallows model.  _arXiv:1904.03099._ , 2019b.
  * Marlin [2004] B. M Marlin.  Modeling user rating profiles for collaborative filtering.  In _Advances in neural information processing systems_ , pages 627–634, 2004.
  * Mnih and Salakhutdinov [2008] A. Mnih and R.R. Salakhutdinov.  Probabilistic matrix factorization.  In _Advances in neural information processing systems_ , pages 1257–1264, 2008.
  * Ni et al. [2019a] Yang Ni, Yuan Ji, and Peter Mueller.  Consensus monte carlo for random subsets using shared anchors.  _arXiv preprint arXiv:1906.12309_ , 2019a.
  * Ni et al. [2019b] Yang Ni, Peter Müller, and Yuan Ji.  Bayesian double feature allocation for phenotyping with electronic health records.  _Journal of the American Statistical Association_ , pages 1–15, 2019b.
  * Paterek [2007] Arkadiusz Paterek.  Improving regularized singular value decomposition for collaborative filtering.  In _Proceedings of KDD cup and workshop_ , volume 2007, pages 5–8, 2007.
  * Sarwar et al. [2001] Badrul Sarwar, George Karypis, Joseph Konstan, and John Riedl.  Item-based collaborative filtering recommendation algorithms.  In _Proceedings of the 10th international conference on World Wide Web_ , pages 285–295, 2001.
  * Schafer et al. [2007] J. B. Schafer, D. Frankowski, J. Herlocker, and S. Sen.  Collaborative filtering recommender systems.  In _The adaptive web_ , pages 291–324. Springer, 2007.
  * Scott et al. [2016] Steven L Scott, Alexander W Blocker, Fernando V Bonassi, Hugh A Chipman, Edward I George, and Robert E McCulloch.  Bayes and big data: The consensus monte carlo algorithm.  _International Journal of Management Science and Engineering Management_ , 11(2):78–88, 2016.
  * Sedhain et al. [2015] S. Sedhain, A. K. Menon, S. Sanner, and L. Xie.  Autorec: Autoencoders meet collaborative filtering.  In _Proceedings of the 24th International Conference on World Wide Web_ , pages 111–112. ACM, 2015.
  * Vitelli et al. [2017] V. Vitelli, Ø. Sørensen, M. Crispino, A. Frigessi, and E. Arjas.  Probabilistic preference learning with the mallows rank model.  _Journal of Machine Learning Research_ , 18, 158-1, 2017.