Spatial Processes for Recommender Systems

1 Introduction

Recommender systems (RS) are designed to direct users to personally interesting items in situations where the amount of available information exceeds the users' processing capability [Resnick and Varian, 1997; Burke, 2002]. Typically, such systems (1) use information about a user to predict ratings of items that the user has not yet considered, and (2) recommend suitable items based on these predictions. Collaborative modelling techniques constitute one of the main model classes applied in RS [Albrecht and Zukerman, 2007]. They base their predictions upon the assumption that users who have agreed in their behaviour in the past will agree in the future.

The greatest strength of collaborative approaches is that they are independent of any representation of the items being recommended, and work well for complex objects, for which features are not readily apparent. The two main collaborative approaches are memory-based and model-based. Previous research has mainly focused on memory-based approaches, such as nearest-neighbour models (classic collaborative filtering), e. g., [Herlocker et al., 1999], due to their intuitiveness. The main drawback of memory-based algorithms is that they operate over the entire user database to make predictions. In contrast, model-based approaches use techniques such as Bayesian networks, latent-factor models and artificial neural networks, e. g., [Breese et al., 1998; Bell et al., 2007], to first learn a statistical model in an offline fashion, and then use it to make predictions and generate recommendations. This decomposition can significantly speed up the recommendation generation process.

Spatial processes (random fields) are a subclass of stochastic processes which are applied to domains that have a geospatial interpretation, e. g., [Diggle et al., 1998; Banerjee et al., 2004]. Typical tasks of the field of spatial statistics include modelling spatial associations between a set of observations made at certain locations, and predicting values at locations where no observations have been made.

This paper applies theory from the area of spatial statistics to the prediction of a user's interests or item ratings in RS. We develop a Spatial Process Model (SPM) by adapting a Gaussian spatial process model to the RS scenario, and demonstrate our model's applicability to the task of predicting implicit ratings in the museum domain. The use of spatial processes requires a measure of distance between items in addition to user ratings. This measure, which is non-specific (e. g., it may be a physical or a conceptual distance), can be readily obtained in most cases.

Our application scenario is motivated by the need to automatically recommend exhibits to museum visitors, based on non-intrusive observations of their actions in the physical space. Employing RS in this scenario is challenging due to (1) the physicality of the domain, (2) having exhibit viewing times rather than explicit ratings, and (3) predictions differing from recommendations (we do not want to recommend exhibits that visitors are going to see anyway). We turn the first challenge into an advantage by exploiting the fact that physical distances between exhibits are meaningful, enabling the use of walking distance between exhibits to calculate (content) distance. This supports the direct, interpretable application of spatial processes by using a simple parametric Gaussian spatial process model (with the ensuing low variance in parameter estimates), compared to more complex non-parametric approaches, e. g., [Schwaighofer et al., 2005]. The second challenge, which stems from the variable semantics of viewing times (time t for different exhibits could mean interest or boredom), is naturally addressed by SPM's structure. The third challenge can be addressed by (a) using SPM to build a model of a visitor's interests in exhibits, (b) inferring a predictive model of a visitor's pathway through the museum [Bohnert et al., 2008], and (c) combining these models to recommend exhibits of interest that may be overlooked if the predicted pathway is followed.

Our approach offers advantages over other model-based approaches in that, unlike neural networks (and memorybased techniques), it returns the confidence in a prediction, and its parameters have a clear interpretation; unlike Bayesian networks, our model does not require a domain-specific adaptation, such as designing the network topology. In addition, the distance measure endows our model with capabilities of hybrid RS [Burke, 2002; Albrecht and Zukerman, 2007] by seamlessly supporting the incorporation of other types of models (e. g., content-based). The distance measure also alleviates the cold-start problem. The new-item problem is addressed by utilising the (distance-based) correlation between this item and the other items. The new-user problem is similarly handled through the correlation between items rated by a user and the other items (our model can make useful personalised predictions after only one item has been rated).

SPM was evaluated with a real-world dataset of time spans spent by museum visitors at exhibits (viewed as implicit ratings). We compared our model's performance to that of (1) a baseline model which delivers a non-personalised prediction, and (2) a state-of-the-art nearest-neighbour collaborative filter incorporating performance-enhancing modifications, e. g., [James and Stein, 1961; Herlocker et al., 1999]. Our results show that SPM significantly outperforms both models.

The paper is organised as follows. Our spatial processes approach for modelling and predicting item ratings is described in Section 2. In Section 3, we present the results of our evaluation, followed by our conclusions in Section 4.

2 Using Spatial Processes toModel and Predict Ratings in Recommender Systems

We briefly introduce stationary spatial process models (Section 2.1), before adapting a model based on Gaussian spatial processes to RS (Section 2.2). In Section 2.3, we describe an MCMC-based Bayesian approach for estimating the parameters of the model, and in Section 2.4 we outline the theory for predicting item ratings.

2.1 Stationary Spatial Process Models

Let Y = (Y (s1), . . . , Y (sn)) be a vector of observations at n sites s1, . . . , sn. Assuming stationarity both in the mean and variance, we can define the following basic model to capture spatial associations:

where w(si) are assumed to be realisations from a stationary Gaussian spatial process W (s) with mean 0, variance 1, and isotropic 1 correlation function ρ(‖si - sj‖; φ, ν ) capturing residual spatial association; and ε(si) are realisations from a white-noise process with mean 0 and variance τ 2 , i. e., non-spatial uncorrelated error terms. That is, we assume that Y = (Y (s1), . . . , Y (sn)) are observations from a stationary Gaussian spatial process over s with mean μ, variance σ2 + τ 2 , and correlation function ρ(‖si - sj‖; φ, ν ). Generally, correlation is assumed to approach 0 with increasing distance ‖si - sj‖. Common choices for isotropic correlation functions are the powered exponential

1 A correlation function is isotropic if it depends on the separation vector si - sj only through its length ‖si - sj‖.

where φ > 0 and 0 < ν < 2, or correlation functions from the Mat´ern class [Banerjee et al., 2004]. 2

The model provides a generic framework for modelling spatial associations between locations. Given estimates for μ, σ2 and τ 2 , it can make predictions for new locations. This model can be extended to observations Y that would not naturally be modelled using a Gaussian distribution, e. g., discrete binary variables indicating like/dislike [Diggle et al., 1998].

2.2 Adaptation for Recommender Systems

RS help users find interesting information in a space of many options. Typically, RS identify items that suit the needs of a particular user given some evidence about his/her preferences. A collaborative filter, for example, utilises a set of ratings Y of users U = {u : u = 1, . . . ,m} regarding a set of items I = {i : i = 1, . . . , n} to identify users who are similar to the current user. It then predicts this user's ratings on the basis of the ratings of the most similar users. Ratings of different users are usually considered to be independent, whereas ratings of related items tend to be correlated. Introducing a notion of spatial distance between items in order to functionally specify this correlation structure, we can use spatial process models (Section 2.1) for the prediction task, in a fashion similar to the Gaussian process model for preference prediction described by Schwaighofer et al. [2005]. The assumption made for spatial processes, that correlation between observations increases with decreasing site distance, fits well with RS, where ratings are usually more correlated the closer (i. e., more related) items are. As for the previous section, we use s1, . . . , sn to denote the locations of items i, j = 1, . . . , n in a space providing such a distance measure, i. e., ‖si -sj‖. For example, ‖si - sj‖ could be computed from feature vectors representing the items (similarly to content-based RS), from item-to-item similarities (similarly to item-to-item collaborative filtering [Sarwar et al., 2001]), or from physical distance, as done in this paper.

The following changes are necessary to adapt the model presented in the previous section to RS:

• Multiple ratings. Given a set of users U, we can have multiple (but independent) ratings for a given item i.
• Non-stationarity. Different items can have different rating means and variances. Hence, the underlying process cannot be assumed to be stationary in its mean and variance, as both μ and σ2 depend on an item's location s. We use the notation μ(s) and σ2 (s) to indicate this.
• Item set finiteness. In contrast to traditional geospatial modelling, we require predictions only for a finite set of items, i. e., those at locations s1, . . . , sn. Hence, it is sufficient (and necessary) to know μ(s) and σ2 (s) only at these locations. That is, we do not require a special (functional) structure for μ(s) or σ2 (s), which is usually the case for geospatial models. 3

2 In our experiments (Section 3), we use a powered exponential correlation function, as the Mat´ern class yielded inferior results.

3 In order to compute predictions for a new item i without ratings, μ(si) and σ2 (si) must be externally estimated and supplied to our model, until they can be estimated from observed data.

For a single user with rating vector Y with (Y )i = Y (si), we extend the basic model from Section 2.1 as follows. We set Y (si) to be observations from a non-stationary spatial process, i. e., Y (si) = μ(si)+σ(si)w(si)+ε(si). Exploiting finiteness, μ(si) and σ(si) are respectively components of μ = (μ(s1), . . . , μ(sn)), the vector of mean ratings for items 1, . . . , n, and σ = (σ(s1), . . . , σ(sn)), the vector of standard deviations. Let θ = (μ, σ, τ 2 , φ, ν ) be a vector collecting all model parameters, and W = (w(s1), . . . , w(sn)). Then,

where 1n is the identity matrix of dimension n×n. Note that Y = (Y (s1), . . . , Y (sn)) are mutually independent given θ and W. Marginalising the model over W, we obtain 4

where H(φ, ν ) denotes a correlation matrix with components (H(φ, ν ))ij = ρ(‖si - sj‖; φ, ν ). That is, (H(φ, ν ))ij represents the correlation between the ratings for items i and j.

We are ready to generalise the model to the multi-user case. As above, we denote the set of users with U (cardinality m), and the set of items with I (cardinality n). Typically, for a user u in U, we have ratings for only a subset of I, say for nu items in I. Denoting a rating by user u for item i with Yu(si) and a user's rating vector with Yu, we collect all observed ratings into a vector Y = (Y1, . . . , Ym) of dimension ∑mu=1 nu. Similarly, we structure μ and σ such that μ = (μ1, . . . , μm) and σ = (σ1, . . . , σm), where μu and σu are the vectors of means and standard deviations for those items rated by a user u, respectively. For example, assume that U = {1, 2} and I = {1, 2, 3}. If user 1 rated items 2 and 3, and user 2 rated items 1 and 2, then

Similarly to Equation 1, Y | θ is multivariate normal of dimension ∑mu=1 nu, where H(φ, ν ) is block diagonal with diagonal elements H1(φ, ν ), . . . , Hm(φ, ν ) (due to users u = 1, . . . ,m being independent), and Hu(φ, ν ) denotes a user u's correlation matrix of dimension nu × nu. That is, for all users u = 1, . . . ,m,

Thus, given the model parameters θ= (μ, σ, τ 2 , φ, ν ), our model is fully specified, with μ = (μ(s1), . . . , μ(sn)) and σ=(σ(s1), . . . , σ(sn)).

2.3 Parameter Estimation

This section describes efficient algorithms for estimating the 2n+3 model parameters θ = (μ, σ, τ 2 , φ, ν ). The most popular parameter estimation strategies are maximum-likelihood and Bayesian inference. We opt for a Bayesian solution, as it offers some attractive advantages over the classic frequentist Algorithm 1 Slice Gibbs sampling algorithm

4 Marginalisation over W is possible only in the Gaussian case, not in the more general case described in [Diggle et al., 1998].

• 1: Initialise θ, e. g., by drawing θ from p(θ).
• 2: repeat Updating of auxiliary variable V | θ, Y .
• 3: Draw Z ∼ Exp (1), and set V = l (θ; Y ) + Z.

Component-wise updating of θ | V, Y .

• 4: for k = 1, . . . , |θ| do
• 5: repeat
• 6: Draw the k-th component θk of θ from p(θk), using shrinkage sampling to truncate the domain of p(θk) after each iteration.
• 7: until l (θ; Y ) < V .
• 8: end for
• 9: Keep acquired sample of θ.
• 10: until the number of MCMC samples of θ from p(θ|Y ) is sufficiently large.

approach. For instance, prior knowledge can be formally incorporated into parameter estimation via the prior distribution p(θ), and the uncertainty about the parameters θ is captured by the posterior distribution. Parameter estimates for θ can be obtained from the posterior distribution

where p(Y |θ) is the likelihood of Y given θ. Typically, independent priors are chosen for the different parameters, i. e., in our case, p(θ) = p(μ)p(σ)p(τ 2 )p(φ)p(ν ).

The integrations required to calculate p(θ|Y ) in Equation 3 are generally not tractable in closed form. However, p(θ|Y ) can be approximated numerically using Markov chain Monte Carlo (MCMC) integration methods, such as the Metropolis-Hastings algorithm and the Gibbs sampler. 5 Following Banerjee et al. [2004], we use slice sampling [Neal, 2003], i. e., a slice Gibbs sampler, to sample from the posterior distribution p(θ|Y ). This approach is favoured by Banerjee et al., because it does not require tuning that is tailored to the application, and hence provides an automatic MCMC algorithm for fitting Gaussian spatial process models. Algorithm 1 summarises our sampling procedure. The algorithm consists of two iteratively applied steps: (1) slicing the likelihood, and (2) performing Gibbs updates using draws from the prior along with rejection sampling. When updating a component θk of θ, we use shrinkage sampling to reduce the number of draws required before a point in the slice is found [Neal, 2003].

For our model (Equation 2), the marginal negative loglikelihood l(θ; Y ) = - log p(Y |θ) associated with Y is

where Σu = Σu(σu, τ 2 , φ, ν ) = σu1nu Hu(φ, ν )σu1nu + τ 21nu , and C ≡ const. independent of θ. Given σu, φ and

5 Refer to [Andrieu et al., 2003] for an excellent introduction to MCMC for machine learning.

ν, computing the eigen decomposition of positive-definite σu1nu Hu(φ, ν )σu1nu simplifies the calculation of log |Σu| and Σ-1 u . It also speeds up the sampling procedure when updating τ 2 at a given iteration of the slice Gibbs sampling algorithm. To minimise the number of eigen decompositions, we update φ and ν together. We proceed similarly for the components of σu, and hence for the components of σ.

In the following section, we set σ(si) = √σ2Y (si) - τ 2 , where σ2Y (si) denotes the sample variance of the ratings for item i (at site si), calculated from the ratings Yu(si). This reduces the number of free model parameters from 2n + 3 to n + 3, i. e., θ = (μ, τ 2 , φ, ν ), and significantly speeds up the slice Gibbs sampler.

2.4 Prediction

Given θ, the prediction of a user u's ratings of unseen items, say Yu,1, from a vector of observed ratings Yu,2 is straightforward. That is, we can use standard multivariate normal theory, because Yu = (Yu,1,Yu,2) | θ is normally distributed (Section 2.2, similarly to Equation 1). If we use the following notation

then the conditional distribution p (Yu,1|Yu,2, θ) is normal with mean vector and covariance matrix

The expectation E (Yu,1|Yu,2, θ) represents a personalised prediction of ratings Yu,1, and a measure of confidence can be easily derived from Cov (Yu,1|Yu,2, θ).

3 Evaluation with a Real-World Dataset

This section reports on the results of an evaluation performed with a real-world dataset from the museum domain, including comparison with a state-of-the-art collaborative filter.

In an information-seeking context, viewing time correlates positively with preference and interest. This observation was used in [Bohnert et al., 2008] to propose a formulation of visitors' interests based on viewing times of exhibits. In this paper, we evaluate the predictive accuracy of our model on the basis of its predictions of viewing times.

3.1 Dataset andModel Justification

We obtained the dataset by manually tracking visitors to Melbourne Museum (Melbourne, Australia) from April to June 2008. In general, visitors do not require recommendations to travel between individual, logically related exhibits in close physical proximity. Hence, with the help of museum staff, we grouped the collection exhibited at Melbourne Museum, which comprises a few thousand exhibits, into 126 coherent exhibit areas. We tracked first-time adult visitors travelling on their own, to ensure that neither prior knowledge about the museum nor other visitors' interests influenced a visitor's decisions about what to look at. The resulting dataset comprises 158 complete visitor pathways in the form of time-annotated sequences of visited exhibit areas, with a total visit length of 291:22:37 hours, and a total viewing time of 240:00:28 hours. 6 A total of 8327 exhibit areas were viewed, yielding 52.7 areas per visitor on average (41.8% of the exhibit areas). Hence, 58.2% of the entries are missing from the viewing time (rating) matrix, indicating a potential for pointing a visitor to personally relevant but unvisited exhibit areas. Table 1 summarises further statistics of the dataset.

Clearly, the deployment of an automated RS in a museum requires suitable positioning technologies to non-intrusively track visitors, and models to infer which exhibits are being viewed. Although our dataset was obtained manually, it provides information that is of the same type as information inferable from sensing data. Additionally, the results obtained from experiments with this dataset are essential for model development, as they provide an upper bound for the predictive performance of our model.

The museum space is carefully themed by curatorial staff, such that closely-related exhibits are in physical proximity. Based on this observation, we hypothesise that physical walking distance between exhibits is inversely proportional to their (content) similarity. Thus, in our experiments, we use physical walking distance for measuring (content) distance between exhibits. To calculate walking distances, we employed an SVG file-based representation of Melbourne Museum's site map, mapped onto a graph structure which preserves the physical layout of the museum (i. e., preventing paths from passing through walls or ceilings), and normalised the distances to the interval [0, 1].

We used the Bayesian Information Criterion (BIC) to select the most appropriate family of probability distributions for approximating the distribution of viewing times at each exhibit area. We tested exponential, gamma, normal, lognormal and Weibull distributions. The log-normal family fitted best, with respect to both number of best fits and average BIC score (averaged over all exhibit areas). Hence, we transformed all viewing times to their log-equivalent to obtain normally distributed data.

3.2 Parameter Estimation

We performed slice Gibbs sampling (Algorithm 1) to obtain estimates for the parameters θ of our Spatial Process Model (SPM) (Section 2). For each of the 129 free model parameters, 7 we used (uninformative) independent uniform distance prior distributions. We used every 20-th sample after a burnin phase of 1000 iterations as a sample of θ from p(θ|Y ), and stopped the sampling procedure after 8000 iterations. Thus, in total, this provides 350 samples of θ from p(θ|Y ).

6 For our experiments, we ignore travel time between exhibit areas, and collapse multiple viewing events of one area into one event.

7 We set σ(si) = pσ2Y (si) - τ 2 to speed up the sampling process (Section 2.3), which reduces the number of free parameters from 255 to 129.

Figure 1: ρ(‖si - sj‖; φ = 28.6334, ν = 0.2583)

For visualisation purposes, we ran a slice Gibbs sampler on the complete dataset of 8327 log viewing times, and used the output of this run to obtain posterior mean point estimates of φ (28.6334), ν (0.2583) and τ 2 (0.1578); we omit the estimates for μ and σ due to space limitations. Figure 1 depicts a plot of ρ(‖si -sj‖; φ, ν ) for this parameterisation of φ and μ, showing the shape of the fitted powered exponential correlation function. The dashed lines indicate the correlation functions obtained with φ and ν set to the lower and upper bounds of their 95% credible intervals respectively. The correlation function rapidly drops to values around 0.4, and then more slowly approaches 0.1 as ‖si-sj‖ approaches 1.0. The shape of this function confirms the existence of a relatively high correlation between viewing durations at exhibit areas in close physical proximity.

3.3 Predictive Model Performance

Experimental Setup

To evaluate SPM's predictive performance, we implemented two additional models: Mean Model (MM) and Collaborative Filter (CF). MM, which we use as a baseline, predicts the log viewing time of an exhibit area i to be its (nonpersonalised) mean log viewing time μ(si). For CF, we implemented a nearest-neighbour collaborative filtering algorithm, and added modifications from the literature that improve its performance, such as shrinkage to the mean [James and Stein, 1961] and significance weighting [Herlocker et al., 1999]. Additionally, to ensure that varying exhibit area complexity does not affect the similarity computation for selecting the nearest neighbours (viewing time increases with exhibit complexity), we transformed the log viewing times into z-scores by normalising the values for each of the exhibit areas separately. Visitor-to-visitor differences with respect to their mean viewing durations were neutralised by transforming predictions to the current visitor's viewing-time scale [Herlocker et al., 1999]. We tested several thousand different parameterisations, but in this paper, we report only on the performance of the best one.

Due to the relatively small dataset, we used leave-one-out cross validation to evaluate the performance of the different models. That is, for each visitor, we trained the models with a reduced dataset containing the data of 157 of the 158 visit trajectories, and used the withheld visitor pathway for testing. For SPM, we first obtained the posterior mean of p(θ|Y ) by performing slice Gibbs sampling on the training data (Section 3.2), and then used this posterior mean to compute predictions by conditioning a multivariate normal distribution (Section 2.4). For CF, predictions were computed from the ratings of the nearest neighbours; and for MM, we used μ(si), estimated from the appropriate reduced dataset, as a prediction.

We performed two types of experiments: Individual Exhibit (IE) and Progressive Visit (PV).

• IE evaluates predictive performance for a single exhibit. For each observed visitor-exhibit area pair (u, i), we removed the observation Yu(si) from the vector of visitor u's log viewing durations, and computed a prediction ˆ Yu(si) from the other observations. This experiment is lenient in the sense that all available observations except the observation for exhibit area i are kept in a visitor's viewing duration vector.
• PV evaluates performance as a museum visit progresses, i. e., as the number of viewed exhibit areas increases. For each visitor, we started with an empty visit, and iteratively added each viewed exhibit area to the visit history, together with its log viewing time. We then predicted the log viewing times of all yet unvisited exhibit areas.

For both experiments, we used the mean absolute error (MAE) to measure predictive accuracy as follows:

where Iu denotes a visitor u's set of exhibit areas for which predictions were computed. For IE, we calculated the total MAE for all valid visitor-exhibit area pairs; and for PV, we computed the MAE across the yet unvisited exhibit areas and all visitors for each time fraction of a visit (to account for different visit lengths, we normalised all visits to a length of 1).

Results

Table 2 shows the results for the IE experiment, where SPM achieves an MAE of 0.7548 (stderr 0.0066), outperforming both MM and CF. The performance difference between SPM and the other models is statistically significant with p ≪ 0.01.

The performance of SPM, CF and the baseline MM for the PV experiment is depicted in Figure 2. CF outperforms MM slightly (statistically significantly for visit fractions 0.191 to 0.374 and for several shorter intervals later on, p < 0.05). More importantly, SPM performs significantly better than both MM and CF (statistically significantly for visit fractions 0.019 to 0.922, p < 0.05). Drawing attention to the initial portion of the visits, SPM's MAE decreases rapidly, whereas the MAE for MM and CF remains at a higher level. Generally, the faster a model adapts to a visitor's interests, the more likely it is to quickly deliver (personally) useful recommendations. Such behaviour in the early stages of a museum visit is essential in order to build trust in the RS, and to guide a visitor in a phase of his/her visit where such guidance is most likely needed. As expected, MM performs at a relatively constant MAE level. For CF and SPM, we expected to see an improvement in performance (relative to MM) as the number of visited exhibit areas increases. However, this trend is rather subtle (it can be observed when plotting the models' performance relative to MM). Additionally, for all three models, there is a performance drop towards the end of a visit. We postulate that these phenomena may be explained, at least partially, by the increased influence of outliers on the MAE as the number of exhibit areas remaining to be viewed is reduced with the progression of a visit. This influence in turn offsets potential gains in performance obtained from additional observations. Our hypothesis is supported by a widening in the standard error bands for all models as a visit progresses, in particular towards the end (not shown in Figure 2 for clarity of presentation). However, this behaviour requires further, more rigorous investigation.

4 Conclusions and FutureWork

In this paper, we utilised the theory of spatial processes to develop a model-based approach for predicting users' interests or item ratings in RS. We applied our model to a real-world dataset from the museum domain, where our model attains a higher predictive accuracy than state-of-the-art nearestneighbour collaborative filters. In addition, under the realistic Progressive Visit setting, our model rapidly adapts to a user's ratings (starting from as little as one rating), thus alleviating the new-user problem common to collaborative filtering.

Our dataset is relatively small compared to other real-world RS applications. Although a high number of ratings per user slows down the slice Gibbs sampler due to repeated inversion of matrices of high dimension, employing our model with larger datasets should not represent a problem in practice. This is because the number of ratings per user is usually small compared to the number of users and items, and the computational complexity of evaluating the likelihood function depends only linearly on the number of users in the database.

In the future, we intend to hybridise our model by incorporating content-based item features into our distance measure. We also plan to extend our model to fit non-Gaussian item ratings, e. g., [Diggle et al., 1998; Yu et al., 2006], and consider negative correlations between items.

Acknowledgments

This research was supported in part by grant DP0770931 from the Australian Research Council. The authors thank Carolyn Meehan and her team from Museum Victoria for their assistance; and David Abramson, Jeff Tan and Blair Bethwaite for their help with the computer cluster.