Dynamic Bayesian Logistic Matrix Factorization for Recommendation with Implicit Feedback

1 Introduction

Personalized recommendation systems have been extensively used in our daily lives. A lot of algorithms and real systems have been developed [Shi et al., 2014]. Practically, matrix factorization models are the most successful approaches that have been widely applied to solve various recommendation problems [Hu et al., 2014; Liu et al., 2016; 2017]. From the data perspective, two types of data, i.e., users' explicit feedback and implicit feedback, are usually used to build recommendation models. Explicit feedback refers to users' explicit ratings to items, e.g., the 1 to 5 scale business ratings in Yelp [Hu et al., 2014]. Recommendation models based on explicit feedback data mainly focus on predicting the exact rating value that a user would like to assign to an item. The implicit feedback is derived from users' interactions with items, e.g., the number of times a user has visited a location [Liu et al., 2014]. Recommendation models based on implicit feedback data usually focus on predicting whether a user would interact with an item. However, explicit feedback data is often very limited, where implicit feedback may be useful for learning users' preferences. Thus, in this work, we exploit implicit feedback data to build recommendation model.

Various matrix factorization methods have been proposed to build recommendation models based on implicit feedback data [Rendle et al., 2009; Pan and Chen, 2013; Johnson, 2014; Liu et al., 2015; He et al., 2016], where the latent embeddings of users and items are assumed to be static. In practice, this assumption is generally not valid, since users' preferences and item attributes are constantly drifting over time. To exploit the temporal dynamics for recommendation, different dynamic matrix factorization models have been proposed [Xiong et al., 2010; Gultekin and Paisley, 2014; Charlin et al., 2015; Lian et al., 2016; Hosseini et al., 2017]. The key assumptions of these methods are that the latent embeddings of users and items are homogeneously evolving with the same diffusion rate.

However, different dimensions of the user and item latent embeddings shall evolve heterogeneously with inconsistent diffusion rates. More concretely, certain dimension of embeddings that reflects users' long-term preferences (e.g., payment method) or items' stable attributes (e.g., item category) may evolve slowly with time. On the other side, the dimensions of embeddings reflecting users' short-term preferences (e.g., usage of social media) or items' dynamic attributes (e.g., monthly popularity) should quickly change over time. In this sense, it would be more desirable to simultaneously estimate the latent embeddings and their corresponding diffusion rates, where more accurate estimations of latent embeddings can be expected. This problem, however, has not yet been adequately investigated in existing literature.

In this paper, we have proposed a novel matrix factorization method, named Dynamic Bayesian Logistic Matrix Factorization (DBLMF), to bridge this gap. More specifically, DBLMF extends logistic matrix factorization [Johnson, 2014] to model the dynamic likelihood that a user would interact with an item at a given timestamp. To characterize the temporal dynamics, a temporal diffusion process is developed for learning dynamic latent embeddings of users and items. In particular, the latent embeddings and the diffusion process are modeled in a probabilistic manner with model conjugacy. Therefore, the latent embeddings and the diffusion rate of each dimension can both be desirably inferred from users' implicit feedback data. To infer all the unknown parameters, we have proposed an efficient variational Bayesian inference algorithm, which adopts a sampling strategy to scale up the proposed method. In addition, we have also conducted extensive experiments on real datasets to demonstrate the effectiveness of the proposed algorithm. The experimental results indicate the proposed DBLMF method usually achieves better results than several state-of-the-art recommendation algorithms, in terms of Precision, Recall, and Normalized Discounted Cumulative Gain (NDCG).

2 RelatedWork

In the literature, different dynamic matrix factorization methods have been developed to explore users' dynamic preferences. For example, a set of time-aware recommendation algorithms based on users' explicit feedback have been proposed, which assigned higher weights to recent observations [Liu et al., 2010], or incorporated the temporal information into latent factor models [Koren, 2009; Xiong et al., 2010]. Moreover, Gultekin and Paisley proposed the collaborative Kalman filter that modeled each low-dimensional latent embedding as a multidimensional Brownian motion [Gultekin and Paisley, 2014]. Yu et al. proposed a temporal regularized matrix factorization framework for temporal learning and forecasting based on high dimensional time series data with missing values [Yu et al., 2016]. For implicit feedback data, the Gaussian process had been proposed to capture users' dynamic preferences [Liu, 2015]. In addition, Charlin et al. developed a dynamic Poisson factorization model that exploited Kalman filter to model evolving latent embeddings and used Poisson distribution to model the user-item interactions [Charlin et al., 2015]. Du et al. developed a convex optimization framework that connected self-exciting point processes and low-rank models to exploit recurrent user activity patterns for recommendation [Du et al., 2015]. Recently, Wang et al. proposed a latent feature process method to study the co-evolving patterns of user and item features [Wang et al., 2016]. Hosseini et al. proposed a recurrent Poisson factorization framework that utilized Poisson process to model users' recurrent activity patterns from their implicit feedback data for recommendation [Hosseini et al., 2017].

3 The Proposed RecommendationModel

In this work, we consider the dynamic recommendation problem based on the implicit feedback data generated by M users {ui}Mi=1 over N items {vj}Nj=1 . Firstly, we use a binary matrix Y t ∈ RM×N to describe the interactions between users and items at time t. If a user ui has interacted with an item vj at time t, we set yt ij to 1; otherwise, we set yt ij to 0. Moreover, we denote the set of observed user-item interactions at time t by Ot . We use E t i = {vj|yt ij = 1}Nj=1 to denote the set of items that have interactions with ui at time t. Similarly, we use F t j = {ui|yt ij = 1}Mi=1 to denote the set of users that have interactions with vj at time t. Following traditional matrix factorization methods, we use a vector uti ∈ R1×d to denote the embedding of ui at time t, where d is the dimensionality of latent space and d ≪ min(M,N). We denote the latent embedding of vj at time t by vt j ∈ R1×d . Following [Johnson, 2014], we define the probability that ui would like to interact with vj at time t as follows:

where σ(x) = 1/(1 + e-x ). The likelihood of the observed user-item interactions at time t can be defined as follows:

where wt ij is a dynamic weight, and I[·] is an indicator function. Similar to [Johnson, 2014], we empirically set wt ij > 1. In other words, we assign larger weights to positive user-item interactions (i.e., yt ij = 1) than "negative" user-item interactions (i.e., yt ij = 0), to balance the positive and "negative" interactions for better results. More sophisticated approaches, e.g., the re-weighted probabilistic models proposed in [Wang et al., 2017], can also be integrated into Eq. (2).

As discussed earlier, the user and item embeddings should evolve with time in a heterogeneous way. Thus, we develop a temporal diffusion process to model the heterogeneous evolution patterns of latent embeddings. Specifically, the following priors are imposed on the latent vectors ut+1 i and vt+1 j ,

where N (µ, σ) is the Gaussian distribution with mean µ and variance σ. For simplicity, we set α ∈ Rd×d and β ∈ Rd×d to be diagonal matrices. The diagonal elements are {αk}dk=1 and {βk}dk=1 . Particularly, the kth diagonal entries in α and β control the drifting rates for the kth dimension of users' and items' embeddings, respectively. At time t = 1, we define p(u1i |u0i ) = N(0, σ0I) and p(v1 j |v0 j ) = N(0, σ0I). In addition, we also place the following Gamma priors on the kth diagonal elements of α and β,

where Γ(a, b) is the Gamma distribution with shape parameters a and inverse scale parameter b.

Let Y ∈ RT×M×N denote the tensor of observed interactions across all times, where T is the number of timestamps. Similarly, U ∈ RT×M×d and V ∈ RT×N×d are used to denote the tensors of user and item latent embeddings across all times, respectively. Then, the joint probability distribution of the observations can be defined as follows:

Based on the above-described probabilistic model, we employ a full Bayesian treatment to infer the unknown parameters from this model.

3.1 Variational Bayesian Inference

For brevity, we denote the parameter set {U, V ,α, β} by Θ. In Bayesian methodology, the objective is to obtain the following full posterior,

However, the normalization term is intractable due to multidimensional integral involved. We resort to variational inference (VI) [Blei et al., 2017] to approximate the posterior distribution p(Θ|Y ) with factorisable assumption on the approximated posterior q(U, V , α, β), which can be expressed as follows:

The objective is to minimize the Kullback-Leibler (KL) divergence between the two distributions q(U, V ,α, β) and p(Θ|Y ). Based on derivation, the approximated posterior minimizing the KL divergence can be expressed as follows:

where E[·]j̸=i represents the expectation taken with respect to (w.r.t.) all Θ except for Θi [Blei et al., 2017].

For probabilistic models that are conjugate, the update rules can be derived based on Eq. (8). However, in the proposed model, the likelihood function is a logistic likelihood function, which is not conjugate with Gaussian priors on the embeddings U and V . To address this issue, the following Gaussian lower bound on the logistic function is employed to approximate the likelihood [Jaakkola and Jordan, 1997],

where λ(ξ) = 1 2ξ (σ(ξ) - 12 ), and ξ is an auxiliary variable that needs to be adjusted to make the bound tight at x = ±ξ. By substituting Eq. (9) into the likelihood in Eq. (2), the lower bound of likelihood function log p(Y t |U t , V t ) is as follows:

where ϕ(ξ t ij ) = logσ(ξ t ij ) - ξt ij 2 + λ(ξ t ij )(ξ t ij )2 . The joint probability distribution can be approximated by replacing the likelihood function with its lower bound.

Based on the updating rule in Eq. (8), the model parameters at time t can be updated in the following iterative fashion: Updating U t and V t : To obtain the optimal approximated posterior, we first assume q(U t ) is fully factorized as q(U t ) = ∏Mi=1 ∏dk=1 q(utik ). Then, we substitute the joint probability into Eq. (8), and obtain the optimal variational distribution for the element utik as follows:

where τ t ijk = ∑ℓ̸=k utiℓ vt jℓ . Due to model conjugacy, it can be known from Eq. (11) that the variational distribution shall follow a Gaussian distribution N(utik |˜utik ,Σtik ), where the mean and variance are as follows:

where ˜ τ t ijk = ∑ℓ̸=k ˜utiℓ ˜ vt jℓ . Similarly, we also assume q(V t ) is fully factorized as q(V t ) = ∏Nj=1 ∏dk=1 q(vt jk ). Then, we can obtain the variational distribution q(vt jk ) in a similar way. It can be seen that q(vt jk ) also follows a Gaussian distribution N (vt jk |˜ vt jk , Ωtjk ) with mean and variance as follows:

Updating α and β: Based on model conjugacy, we can see that the variational distribution for the kth element of α shall follow a Gamma distribution: q(αk) ∼ Γ(˜a0k, ˜b0k), where

The parameter αk can be updated by using the mean as ˜a0k/˜b0k. Similarly, the variational distribution of βk also follows a Gamma distribution as: q(βk) ∼ Γ(˜a1k, ˜b1k), where

Then, βk can be updated using the mean as ˜a1k/˜b1k.

Updating ξ t ij : Since there is no prior assigned to the auxiliary parameter ξ t ij , the optimal value of ξ t ij can be obtained by maximizing the expected log-likelihood function, which can be expressed as follows:

Taking derivative of the expected log-likelihood w.r.t. ξ t ij and setting it to zero, the optimal value of ξ t ij can be obtained as

The details of the proposed optimization algorithm are summarized in Algorithm 1. In this algorithm, the time complexity of each iteration is O(M · N · d + d · (M + N)). Thus, it is not scalable to large datasets. To remedy this issue, we propose the following sampling based inference algorithm.

3.2 Inference by Uniform Sampling

The updating rule in Eq. (11) can be rewritten as follows:

where ρtijk = λ(ξ t ij )[utik (vt jk )2 +2vt jk τ t ijk ]. In practice, a user ui usually only interacts with a small fraction of items, i.e., |E t i | ≪ N, where | · | is the cardinality of a set. Therefore, the main computation portion in Eq. (18) is the expectation ψt ik over all items. To improve computation efficiency, we propose to approximate ψt ik by uniformly sampling over all items as ψt ik ≈ N |S t ik | ∑vj∈S t ik E[ vt jk 2 + ρtijk ], where S t ik denotes the set of sampled items used for optimizing q(utik ).

The approximations of mean and variance of the variational Gaussian distribution for utik in Eq. (12) are as follows:

Algorithm 1: Variational Inference for DBLMF

Similarly, the approximations of mean and variance of the variational Gaussian distribution for vt jk are as follows:

where T t jk denotes the set of sampled users used for optimizing q(vt jk ). Considering model stability, we set S t ik = S t i , T t jk = T t j , 1 ≤ k ≤ d. The final updating rules for the variational distributions of elements in U t and V t are as follows:

where θ is the mean or variance of a variational Gaussian distribution, φ(θ) is the approximation based on uniform sampling, and γ ∈ (0, 1] is a constant which is empirically set at 0.95. Algorithm 2 summarizes the details of the inference algorithm based on uniform sampling. In the experiments, we empirically set |S t i | = |E t i | and |T t j | = |F t j |. Then, the time complexity of each iteration in Algorithm 2 is O(|Ot | · d).

Algorithm 2: Variational Inference for DBLMF by Uniform Sampling

Input : {Ot }Tt=1, d, γ, σ0, a0, b0, a1, b1 Output: {q∗ (U t )}Tt=1, {q∗ (V t )}Tt=1 1 Set αk = 0, βk = 0, ∀1 ≤ k ≤ d; 2 for t = 1, 2, . . . , T do 3 Initialize ˜utik, ˜ vt jk using N(0, 1/d), set Σtik = 1, Ωtjk = 1, ∀1 ≤ k ≤ d, 1 ≤ i ≤ M, 1 ≤ j ≤ N; 4 for iter = 1, 2, · · · ,max iter do 5 for i = 1, 2, . . . ,M do 6 Uniformly sampled items to construct S t i ; 7 Compute ξt ij, ∀vj ∈ E t i ∪S t i using Eq. (17); 8 for k = 1, 2, . . . , d do 9 Compute φ(˜utik), φ(Σtik) using Eq. (19); 10 ˜utik ← (1 - γ)˜utik + γφ(˜utik); 11 Σtik ← (1 - γ)Σtik + γφ(Σtik); 12 for j = 1, 2, . . . , N do 13 Uniformly sampled users to construct T t j ; 14 Compute ξt ij, ∀ui ∈ F t j ∪ T t j using Eq. (17); 15 for k = 1, 2, . . . , d do 16 Compute φ(˜ vt jk), φ(Ωtjk) using Eq. (20); 17 ˜ vt jk ← (1 - γ)˜ vt jk + γφ(˜ vt jk); 18 Ωtjk ← (1 - γ)Ωtjk + γφ(Ωtjk); 19 Update the Gamma distributions for α and β using Eq. (14) and Eq. (15), resepectively;

4 Experiments

4.1 Experimental Setting

Datasets and Setups: We extracted four datasets from Movielens-20M 1 for evaluation. Following [Pan and Chen, 2013; Liu et al., 2015], we kept ratings larger than 3 as implicit feedback. Dataset 2005-3Y was extracted using the following conditions: 1) only ratings between 2005-01-01 and 2010-12-31 were included; 2) ratings before and on 200712-31 were used for training, and ratings after 2007-12-31 were used for testing; and 3) only users with at least 5 ratings in training data and at least one rating in testing data were included. Dataset 2002-6Y was generated in a similar way by using the rating data between 2002-01-01 and 2007-1231 as training data, and the rating data between 2008-01-01 to 2010-12-31 as testing data. As the proposed algorithms focus on modeling users' preferences over time, it is better to have users' activities over a long period. To this end, we further restrict that users must have at least one rating before 2005-12-31 on dataset 2005-3Y, and at least one rating before 2002-12-31 on dataset 2002-6Y. This leads to two more datasets: 2005-3Y-S and 2002-6Y-S. Table 1 summarizes the statistics of these datasets.

The recommendation accuracies of a model are measured by: Precision@K, Recall@K, and NDCG@K, which have been widely applied to evaluate the performances of recommendation algorithms based on implicit feedback [Shi et al., 2014]. In the experiments, we set K=5 and 10. For each metric, we first compute the accuracy for each user, and then report the averaged accuracy over all users.

1 https://grouplens.org/datasets/movielens/

Evaluated Recommendation Models: We compare the following recommendation methods: (1) BPR [Rendle et al., 2009]: This is a pairwise learning-to-rank method developed for recommendation with implicit feedback. (2) LMF [Johnson, 2014]: This is the logistic matrix factorization method that employs MAP (i.e., maximum a posterior) to learn user and item embeddings. (3) eALS [He et al., 2016]: This matrix factorization method weights the missing data based on item popularity and uses element-wise alternating least squares technique to learn latent embeddings. (4) dPF [Charlin et al., 2015]: This method uses Poisson matrix factorization to model users' implicit feedback and utilizes Kalman filter to learn time evolving embeddings. (5) DBLMF: This is the proposed method that uses Algorithm 1 to learn dynamic embeddings. (6) sDBLMF: This is the proposed method that uses Algorithm 2 to learn dynamic embeddings.

For matrix factorization methods, we set the dimensionality of the latent space d at 32. In DBLMF and sDBLMF, we set σ0 = 10, a0 = a1 = 10-4 , and b0 = b1 = 10-4 . Moreover, we define the dynamic weight as wt ij = 1+δrt ij , where rt ij denotes the rating that ui assigned to vj at time t, and we empirically set δ = 1 in the experiments.

4.2 Experimental Results

We have tested two different settings for DBLMF and sDBLMF. The first setting is to have only one time period over the whole dataset (i.e., DBLMF-full and sDBLMF-full). The second setting is to set the time granularity to be 3 months for the proposed methods (i.e., DBLMF-3m and sDBLMF-3m) and dPF. Table 2 summarizes the performances on different datasets. We make the following observations:

• Without considering temporal dynamics, the proposed methods (i.e., DBLMF-full and sDBLMF-full) achieve comparable or even better results than existing recommendation algorithms (i.e., BPR, LMF, and eALS) that are based on implicit feedback data, in terms of all metrics. Moreover, DBLMF-full usually outperforms LMF. This indicates that Bayesian inference can more accurately estimate latent embeddings than MAP by considering higher statistics such as covariance.
• The dynamic matrix factorization methods (i.e., dPF, DBLMF-3m, and sDBLMF-3m) usually achieve better results than BPR, LMF, and eALS, which do not consider temporal information. For instance, on 20026Y-S dataset, DBLMF-3m and sDBLMF-3m outperform eALS by 139.16% and 145.07% respectively, in terms of NDCG@5. This indicates that the recommendation accuracy can be significantly improved by considering temporal dynamics. Additionally, the accuracy improvements obtained by dynamic matrix factorization methods

on 2005-3Y-S and 2002-6Y-S datasets are much more significant than that obtained on other two datasets. The potential reason is that users' activities on 2005-3Y-S and 2002-6Y-S datasets are over a longer period. Thus, these two datasets contain more sufficient temporal information to learn users' dynamic preferences.

• Compared with dPF, DBLMF-3m and sDBLMF-3m achieve better performances on 2005-3Y-S and 2002-6Y-S datasets. For example, in terms of NDCG@5, DBLMF-3m and sDBLMF-3m outperform dPF by 18.56% and 21.50%, on 2002-6Y-S dataset. This is because dPF assumes different dimensions of user/item embeddings are homogeneously evolving with the same diffusion rate. However, the proposed DBLMF methods enable different dimensions of latent embeddings to heterogeneously evolve with inconsistent diffusion rates, so that better performances can be achieved.
• On all datasets, sDBLMF methods (i.e., sDBLMF-full and sDBLMF-3m) achieve comparable results with DBLMF methods (i.e., DBLMF-full and DBLMF-3m). However, sDBLMF methods are tens of times faster than DBLMF methods. This demonstrate the effectiveness of the sampling inference algorithm.

In addition, we have also evaluated the sensitivity to parameters for sDBLMF, because it usually achieves comparable results with DBLMF and is scalable to large datasets. The time granularity is varied in {1, 3, 6, 9, 12, 18, 24, 36, 72} months. The dimensionality of latent space d is adjusted in {8, 16, 32, 64, 128}. The parameter of prior distributions σ0 is varied in {1, 10, 50, 100, 150, 200, 250}. Moreover, we also adjust the number of iteration times from 1 to 200. The results on 2002-6Y-S dataset are summarized in Figure 1. Observed that smaller time granularity generally achieves better results. When the time granularity is set to 3 months, sD- BLMF achieves the best results. We also notice that sDBLMF is relatively stable with respect to different settings of d and σ0. The recommendation accuracy generally increases with the increase of the number of iterations, particularly when the number of iterations is fewer than 10. When the number of iteration times is larger than 10, further increases of iteration times does not improve the performance much. Due to space limitation, we do not report the performance trends regarding with other model parameters in this paper.

5 Conclusions

In this paper, we propose a novel recommendation method, namely Dynamic Bayesian Logistic Matrix Factorization (DBLMF). The proposed model is capable of inferring the dynamic embeddings of users and items, and adapting the diffusion rates of latent embeddings in a data-driven manner. Remarkably, the proposed Bayesian inference framework can obtain the advantages of model flexibility and statistical information incorporation. In addition, a scalable sampling based inference approach has also been developed, which is particularly useful in presence of large datasets. The superiority of the proposed approaches has been validated by comparing with the state-of-the-art baselines on Movielens-20M dataset. As for the future work, we intend to develop different sampling strategies to improve the recommendation accuracy.