Bayesian Factorization Machines

Information Systems & Machine Learning Lab
University of Hildesheim
31141 Hildesheim

Steffen Rendle
Social Network Analysis
University of Konstanz
78457 Konstanz, Germany{freudenthaler, schmidt-thieme}@ismll.de

steffen.rendle@uni-konstanz.de

Christoph Freudenthaler, Lars Schmidt-Thieme

Abstract
This work presents simple and fast structured Bayesian learning for matrix andtensor factorization models. An unblocked Gibbs sampler is proposed for factorization machines (FM) which are a general class of latent variable models subsuming matrix, tensor and many other factorization models. We empirically showon the large Netflix challenge dataset that Bayesian FM are fast, scalable and moreaccurate than state-of-the-art factorization models.

1

Introduction

Many different problems in machine learning such as computer vision, computational biology orrecommender systems deal with complex problems in sparse data. For those problems, e.g., theNetflix prize problem1 , sparse representations such as latent variable models are typically used.
The most basic version are matrix factorization models, extensible to tensor factorization models
for more than two categorical dimensions. In [5], Rendle generalizes these factorization models to
Factorization Machines (FM) and shows how a simple input vector consisting of the appropriate
predictors is enough to mimic most factorization models. The general principle is the same as forstandard linear regression. For instance factorizing a matrix containing the rating of a user u ∈ U
represented in row r = 1, . . . , |U | for a specific item i ∈ I represented in column c = 1, . . . , |I| canbe generated by defining an input vector of length p = |U | + |I| consisting of p indicator variables.
Within the first |U | indicators exactly one indicator is on - the one identifying the current user.
Within the last |I| indicators again exactly one is on - now for the current item.
Typically, the number of latent dimensions k for latent variable models is not known. Too small kunderfit, too high k overfit the data. Structured Bayesian methods using hierarchical priors overcome
brute force grid search for k, e.g., [7, 9] for matrix factorization, [11] for tensor factorization, and[10] multi-relational learning. Each paper derived its own learning algorithm. In this work we comeup with one structured Bayesian learning algorithm for a much wider class of factorization models– among others subsuming matrix, tensor or multi-relational factorization models. Moreover, theruntime complexity of existing algorithms is typically non-linear in k. A first attempt to tackle thisissue is made by Zhu et al. [12]. They adapt the Gibbs sampler of BPMF [7] by unblocking somesampling steps leading to a final runtime complexity of O(nk + k 2 (|U | + |I|)). In our work, wederive a fast and unblocked Gibbs sampler where each iteration is linear in k. The improvement
is based on the work of Rendle et al. [6] where a linear time learning algorithm for FMs using
alternating least-squares is presented. As posterior mean and variance of factors are closely related
to the (conditional) least-squares solution of [6], we can transfer these ideas to derive fast BayesianFMs.
1

http://www.netflixprize.com/

1w , w

 v , vwj

vj

0 , 0 , 

w

v

wj

vj

w
j=1,...,pv
j=1,...,p

x ij
yi

w0
w , w0x ij

w , w



0

yi

w0
i=1,...,n0i=1,...,n0 0 , 0

Figure 1: Comparison of the probabilistic interpretation of standard Factorization Machines (left) toBayesian Factorization Machines (right). BFM extend the standard model by using hyperpriors.

2

Bayesian Factorization Machines

2.1

Factorization Machines (FM)Factorization Machines [5] are a regression model for a target y using p explanatory variables x ∈Rp as:
y(x) := w0 +p
X

w j xj +

p
D X
X

...

d=2 j1 =1

j=1

p
X

wj1 ,...,jd

jd >jd−1

d
Y

xjl

(1)

l=1

where second and higher order parameters wj1 ,...,jd are factorized using PARAFAC[2]:wj1 ,...,jd :=

k Y
d
X
(d)
vjl ,f∀d ≥ 2.

(2)

f =1 l=1

FM subsume (among others) D-way tensor factorization which is representable by engineering the
input features x. Engineering x = (x1 , . . . , xD ) as a sequence of disjoint subsequences xi =xi,1 , . . . , xi,pi where i = 1, . . . , D and ∃!j ∈ xi : xi,j = 1, xi,¬j = 0 ∀i = 1, . . . , D gives a D-way
tensor factorization model. Note that, FM are more complex than tensor factorization models as x
may consist of real valued predictor variables.
2.2

Bayesian Factorization Machines

So far FM lack Bayesian inference. For obvious reasons such as better accuracy, automatic hyperparameter learning (no grid search) or no learning hyperparameter, e.g., learn rate for gradientdescent, we want to enhance FM with structured Bayesian inference. Bayesian theory [1] suggest
Bayesian models with hierarchical hyperpriors as they regularize underlying model parameters. The
intuition behind this is that additional (hierarchical) structure put on model parameters shrinks them
towards each other as long as no evidence in the data clearly indicates the opposite. Without loss ofgenerality but for better readability, we present Bayesian FM (BFM) of order D = 2.
Our hierarchical Bayesian model adds to the standard regularized regression model Gaussian hyperpriors for each mean µθ of all model parameters θ ∈ Θ = {w0 , w1 , . . . , wp , v1,1 , . . . , vp,k } butµw0 as well as Gamma hyperpriors for each prior precision α and λθ but λw0 . The standard L2regularized regression model with hyperparameters ΘH = {λθ , µθ : θ ∈ Θ} is for a sample of nobservations (yi , xi ) ∈ R1+p :p(Θ|y, X, ΘH ) ∝

n
Y
√

2

α

α e− 2 (yi −y(xi ,Θ))

i=1

Yp

λθ

2

λθ e− 2 (θ−µθ )

(3)

θ∈Θ

Figure 1 depicts the graphical representations of the standard and our hierarchically extended
Bayesian model: we define two regularization populations, i.e. one for all linear effects wi andanother one for the latent variable v1,1 , . . . , vp,k .
2

2.3

Inference: Efficient Gibbs Sampling

Since posterior inference for FM is analytically intractable we use Gibbs sampling to draw fromthe posterior. Standard Gibbs sampling is blocked Gibbs sampling, i.e. all inferred variables Θ andΘH b are divided into b disjoint blocks consisting of r1 , . . . , rb variables. However, blocked Gibbs
sampling for our hierarchical Bayes approach leads to higher time complexities in the final Gibbs
sampler as for each block of size ri matrices need to be inverted which is O(ri3 ). Therefore, we
propose single parameter Gibbs sampling. For notational readability, we exploit the multi-linear
nature of FM, i.e. for each θ ∈ Θ eq. 1 can be rewritten as [6]:y(x|Θ) = gθ (x) + θ hθ (x) ∀θ ∈ Θ(4)where gθ , hθ respectively comprise all terms additively, multiplicatively connected to θ (see [6]for details). Using this notation and parameters Θ0 = {γ0 , µ0 , α0 , β0 , αλ , βλ } of the conjugate
hyperpriors the conditional posterior distributions for each θ, µθ , λθ , and α resulting from ourhierarchical Bayesian model are:p(θ|y, X, Θ \ {θ}, ΘH , Θ0 ) = N (µ∗θ , σθ2 ),p(µθ |y, X, Θ, ΘH \ {µθ }, Θ0 ) = N (µµθ , σµ2 θ )p(λθ |y, X, Θ, ΘH \ {λθ }, Θ0 ) = Γ(αθ , βθ ),p(α|y, X, Θ, ΘH \ {α}, Θ0 ) = Γ(αn , βn )with:µ∗θ = σθ2nX!
α(yi − gθ (xi ))hθ (xi ) + µθ λθ

σθ2 =

,

α

!−1
2

hθ (xi ) + λθ

(5)

i=1

i=1

µµθ = σµ2 θ λθ

n
X

p
X

!
θ j + γ 0 µ0,

σµ2 θ = ((p + γ0 )λθ )−1

(6)

j=1

αθ = (αλ + p + 1) /2,

βθ =

p
X

!
2

2

(θj − µθ ) + γ0 (µθ − µ0 ) + βλ

/2

(7)

j=1

αn = (α0 + n) /2,

βn =

n
X
(yi − y(xi , Θ))2 + β0

!
/2

(8)

i=1

The final algorithm iteratively draws samples for hyperparameters ΘH and model parameters Θ.
Sampling of hyperparameters
using eq. 6-8 is simple. With a straight-forward implementation it
Pnis done in O( i=1 k m(Xi. )) = O(k Nz )2 . Efficient computation of posterior mean and posteriorvariance of eq. 5 is more complicated. However the posterior mean is identical to the alternating
least squares solution for which an efficient algorithm has already proposed [6] with an overall
complexity of O(k Nz ) for an update of all model parameters.

3

Evaluation

We evaluate BFM on the Netflix challenge dataset which is well-known in the recommender systems
community. It consists of n = 100, 480, 507 ratings of 480, 189 users for 17, 700 movies. The
train/test split is the same as in the challenge. This allows direct comparison with existing results
reported in [8, 7, 3, 11]. For the parameters of the hyperpriors Θ0 we have chosen trivial values:
α0 = β0 = γ0 = αλ = βλ = 1, µ0 = 0. Please note, due to the high number of explanatory
variables their impact on the Gibbs sampler is negligible (cf. eq. 6-8).
With the derivation of BFM it is possible to learn many different models (e.g. matrix/tensor factorization) Bayesian by appropriate specification of the input data x. Exploiting this flexibility, wetested several models by compiling the feature vector x from the following components: (1) binary
indicators: for each user (u), each item (i), each day of rating (t), each u and t of rating combination(ut), each i and t combination (it)3 , each integer value of the logarithm of the rating user’s rating
frequency on the day of rating (f ), each i and f combination (if ), and for items the current user has
2m(Xi. ) is the number of non-zero elements in the i-th row of the design matrix X. Nz :=is the number of non-zero elements in X.
3
Time is binned in 50 days intervals (following [4]).

3Pni=1 m(Xi. )0.915

●
●
●

RMSE
0.90

●
●
●
●

SGD KNN256
SGD KNN256++
BFM KNN256
BFM KNN256++
SGD PMF
BPMF
BFM (u,i)
BFM (u,i,t,f)
BFM (u,i,ut,it,if)
BFM Ensemble

●

0.88

●

0

50
100
k (#latent Dimensions)

●
●

RMSE
0.895 0.905

●

150

0.885

0.92

Netflix

BPMF k=30
BFM (u,i) k=32
BFM (u,i,t,f) k=200

●
●
●
●

0

●

100

●

200 300 400#Samples●500

Figure 2: Predictive accuracy of various instances of Bayesian Factorization Machines (BFM)
compared to state-of-the-art baselines: SGD PMF [8], BPMF [7], SGD KNN/KNN++ [3]. Rightpanel: Convergence of unblocked Gibbs sampling (BFM) is comparable to blocked Gibbs sampling(BPMF [7]).

co-rated (++). (2) Real-valued variables: the rating values for the co-rated items (KNN)4 . We index
BFM with the features they use, e.g., BFM (u, i) is a Bayesian matrix factorization model using theuser and item indicators u and i. It is similar to BPMF [7]. BFM (u, i, t) uses time information
like the Bayesian PARAFAC model BPTF [11]. BFM KNN and BFM KNN++ are the Bayesianequivalent to SGD KNN(++) [3], etc.
Figure 2 plots the quality with a varying amount of factorization dimensions k. As expected, posterior averages in BPMF and BFM (u, i) outperform posterior modes in standard matrix factorization.
For the same reason BMF KNN(++) outperform SGD KNN(++). Besides, the Bayesian approaches
achieve this without a hyperparameter search. Next, we compare our generic BFM approach to thecorresponding specialized models where Gibbs samplers [7, 11] already exist. Comparing BFM(u, i) to BPMF [7], shows that BFM outperforms BPMF slightly. The reason for this are the additional linear effects in BFM. For the time-aware Bayesian factorization model BPTF [11], the bestvalue the authors report is an RMSE value of 0.9044 for k = 30. For comparison we also run an
BFM (u, i, t) and largely outperform this with 0.8958 (k = 32) and 0.8909 (k = 64). Moreover, incontrast to BPMF and BPTF using blocked Gibbs, BFM achieve these results with a Gibbs samplerof lower computational complexity. Encouraged by the quality of the single factorization models wealso computed a linear ensemble of three BFM: a BFM (u, i, t, f ) with k = 200, a BFM KNN256++
and a BFM (u, i, ut, it, if ) with k = 64. The result of this ensembled BFM is 0.8772.
The right panel of figure 2 depicts the convergence rate of blocked5 and unblocked Gibbs sampling. Surprisingly, the number of iterations for similar models (matrix factorization, k ≈ 30) using
unblocked and blocked Gibbs sampling are comparable.

4

Conclusion

We proposed BFM, a computationally efficient Gibbs sampler for FM. Instead of cubic complexity
BFM are linear in k. Empirically tested on the Netflix challenge, BFM converge as fast as models
using blocked Gibbs sampling. Thus, BFM combine high accuracy and more scalable learning withthe generality of FM which subsume many matrix/tensor factorization models and other (factorized)
polynomial regression models like KNN models.
4
For ++ and KNN we follow the pruning step of [3] and took only the 256 most similar items with respect
to the target item.
5
Parameters are grouped in blocks of size k.

4

References[1] Andrew Gelman, John B. Carlin, Hal S. Stern, and Donald B. Rubin. Bayesian Data Analysis.
Chapman and Hall/CRC, 2nd edition, 2003.
[2] Richard A. Harshman. Foundations of the parafac procedure: models and conditions for an’exploratory’ multimodal factor analysis. UCLA Working Papers in Phonetics, pages 1–84,
1970.
[3] Yehuda Koren. Factorization meets the neighborhood: a multifaceted collaborative filteringmodel. In KDD ’08: Proceeding of the 14th ACM SIGKDD international conference on
Knowledge discovery and data mining, pages 426–434, New York, NY, USA, 2008. ACM.
[4] Yehuda Koren. Collaborative filtering with temporal dynamics. In KDD ’09: Proceedings ofthe 15th ACM SIGKDD international conference on Knowledge discovery and data mining,
pages 447–456, New York, NY, USA, 2009. ACM.
[5] Steffen Rendle. Factorization machines. In Proceedings of the 10th IEEE International Conference on Data Mining. IEEE Computer Society, 2010.
[6] Steffen Rendle, Zeno Gantner, Christoph Freudenthaler, and Lars Schmidt-Thieme. Fast
context-aware recommendations with factorization machines. In Proceedings of the 34th ACM
SIGIR Conference on Research and Development in Information Retrieval. ACM, 2011.
[7] Ruslan Salakhutdinov and Andriy Mnih. Bayesian probabilistic matrix factorization using
Markov chain Monte Carlo. In Proceedings of the 25th International Conference on Machine
Learning, volume 25, 2008.
[8] Ruslan Salakhutdinov and Andriy Mnih. Probabilistic matrix factorization. In Advances in
Neural Information Processing Systems, volume 20, 2008.
[9] Mikkel Schmidt. Linearly constrained bayesian matrix factorization for blind source separation. In Y. Bengio, D. Schuurmans, J. Lafferty, C. K. I. Williams, and A. Culotta, editors,
Advances in Neural Information Processing Systems 22, pages 1624–1632. 2009.
[10] Ilya Sutskever, Ruslan Salakhutdinov, and Joshua Tenenbaum. Modelling relational data using
bayesian clustered tensor factorization. In Y. Bengio, D. Schuurmans, J. Lafferty, C. K. I.
Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems 22,
pages 1821–1828. 2009.
[11] Liang Xiong, Xi Chen, Tzu-Kuo Huang, Jeff Schneider, and Jaime G. Carbonell. Temporalcollaborative filtering with bayesian probabilistic tensor factorization. In Proceedings of the
SIAM International Conference on Data Mining (SDM 2010), pages 211–222. SIAM, 2010.
[12] Shenghuo Zhu, Kai Yu, and Yihong Gong. Stochastic relational models for large-scale dyadic
data using mcmc. In NIPS, pages 1993–2000, 2008.

5