Title: VARIATIONAL AUTOENCODERS FOR HIGHLY MULTI-VARIATE SPATIAL POINT PROCESSES INTENSITIES

Abstract: Multivariate spatial point process models can describe heterotopic data over space. However, highly multivariate intensities are computationally challenging due to the curse of dimensionality. To bridge this gap, we introduce a declustering based hidden variable model that leads to an efficient inference procedure via a variational autoencoder (VAE). We also prove that this model is a generalization of the VAE-based model for collaborative filtering. This leads to an interesting application of spatial point process models to recommender systems. Experimental results show the method's utility on both synthetic data and real-world data sets.

1. INTRODUCTION:

Multivariate point processes are widely used to model events of multiple types occurring in an n dimensional continuum. This paper focuses on multivariate spatial point processes (SPP), which can uncover hidden connections between subprocesses based on the correlations of their spatial point patterns. Often we encounter missing data problems, where some subprocesses are not fully observable. The underlying connections could further contribute to the prediction of these subprocesses over the unobserved areas. Moreno-Muñoz et al. (2018) has shown the effectiveness of this joint model for Gaussian processes with heterotopic data. Multi-output models in Lian et al. (2015) such as coregionalization and cokriging can outperform independent predictions. However, there is limited literature on the statistical methodology of the highly multivariate spatial point processes, according to the very recent paper (Choiruddin et al., 2019) .

Inference for multivariate spatial point processes intensities is still a challenging problem (Taylor et al., 2015) , especially with a large number of subprocesses. For popular Gaussian processes-based approaches (Williams & Rasmussen, 2006) , the multivariate intensity often consists of independent and multi-output Gaussian processes. The complexity of the models and the curse of dimensionality hinder this approach for highly multivariate data, such as friendship networks and recommender systems with millions of users. In these problems, we only partially observe the events (e.g. users interact with items, locations) for each subprocess (user). It is necessary to jointly infer the preference of each user based on their hidden correlations. For example, a common approach in recommender systems, collaborative filtering (He et al., 2017) , predicts the item interests of each user with the help of the collection of item preferences for a large number of users.

To address these problems, we propose a multivariate spatial point process model with a nonparametric intensity. We extend the well-known kernel estimator in Diggle (1985) to the multivariate case. This generalization is achieved through the introduction of hidden variables inspired by stochastic declustering (Zhuang et al., 2002) . The latent variables naturally lead to a variational Bayesian inference approach, which is different from the frequentist point estimation in the kernel estimator. To reduce the complexity in the highly multivariate case, we consider an alternative set of hidden variables that are designed to work well as latent variables for a variational autoencoder (VAE) (Kingma & Welling, 2014) . This amortized inference (Gershman & Goodman, 2014) approach leads to fast inference once the model is fully trained. Further, we show the equivalence for these two different settings of hidden variables using the properties of the spatial point process. This efficient approach makes it possible to apply multivariate spatial point processes in many areas, including locationbased social networks and recommender systems with many users. Moreover, the nonparametric method for analyzing spatial point data patterns is not related to specific parametric families of models, which only requires the intensity to be well-defined.

Our approach is not a direct replacement for current inference methods on few-variate spatial point processes (Jalilian et al., 2015) . In contrast to the classical methodology, VAE requires a large number of training data. The highly multivariate data that are widely available in social networks and recommender systems can be ideal applications for our approach. In fact, it can be shown that our model is a generalization of a state-of-the-art VAE-based collaborative filtering model (Liang et al., 2018) . Our model nonparametrically fits the underlying intensity function. Compared with the multinomial distribution used in Liang et al. (2018) , this leads to not only a smoother intensity over space but also better predictions in terms of ranking-based losses. Compared to a univariate model, such as the trans-Gaussian Cox processes (Williams & Rasmussen, 2006) , our multivariate model enhances the predictive ability on missing or unobserved areas, which is consistent with the results of heterogeneous multi-output Gaussian processes (Moreno-Muñoz et al., 2018) .

The contributions of this paper are three-fold. We first build a novel multivariate spatial point process model and find a direct connection with the VAE-based collaborative filtering through detailed theoretical analysis. Secondly, this connection introduces amortized inference for an efficient multivariate point process estimation. Finally, point processes generalize the discrete distribution used in (Liang et al., 2018) and lead to a better modeling of spatial heterogeneity. We validate these benefits through experiments for multiple multivariate data sets, showing improvement over classic SPP methods and potentials on collaborative filtering applications.

2. PRELIMINARIES:

Spatial point process A point process (PP) is a random counting measure N (x) on a complete separable metric space R (here we always assume that R ⊂ R n ) that takes values on {0, 1, 2, ...} {∞}. While the major theory of point processes centers around the temporal dynamics, spatial point process models (Diggle et al., 1983) are established in forestry and seismology, focusing on the stationary and isotropic case. We focus on the (first-order) intensity function λ(x), which is the expected rate of the accumulation of points around a particular spatial location x. We write

EQUATION

where ∆x is a small ball in the metric space, e.g. the Euclidean space R n , with the centre x and measure |∆x| . The second-order intensity function is naturally defined as

EQUATION

measuring the chance of points co-occurring in both ∆x and ∆y. Normalizing this leads to the pair-correlation function g(x, y) = λ (2) (x, y)/λ(x)λ(y). g(x, y) > 1 indicates that points are more likely to form clusters than the simple Poisson process where g(x, y) = 1.

Common models in SPPs include the Poisson process with a non-stationary rate λ(x), and the Cox process with a nonnegative-valued intensity process Λ(x), which is also a stochastic process. Cox processes conditional on a realization of the intensity process Λ(x) = λ(x) are Poisson processes with intensity λ(x). To model the aggregated points patterns, Poisson cluster (Neyman-Scott) processes generate parent events from a Poisson process. Then each parent independently generates a random number of offsprings. The relative positions of these offsprings to the parent are distributed according to some p.d.f K σ (x) in space (Diggle et al., 1983) . Many point process models, including most Cox processes, are in fact Poisson cluster processes. The duality between Cox processes and cluster processes is widely used to construct Cox process models. For example, the kernel-based intensity process Λ(x) = ∞ i=1 K σ (x -x i ) with x i from a Poisson process, is essentially a Poisson cluster process. The number of offsprings is from a Poisson distribution with λ = 1 and the relative position distribution is K σ (x). Repulsive SPPs, on the other hand, model that nearby points of the process tend to repel each other. Higher order intensities are often considered in this case, such as determinantal PPs.

Alternatively, if we are more interested in the realization intensity λ(x) than the mechanical interpretation, the trans-Gaussian Cox process provides a tractable way to construct the Cox process using a nonlinear transformation on a Gaussian process S(x). Popular choices for Λ(x) include the log-Gaussian Cox process (LGCP) with exp(S(x)) and the permanental process with S(x) 2 . Recent works on Cox processes have been extensively focused on the cases that are modulated via the Gaussian random field, due to its capability in modeling the intensity and pair-correlations between subprocesses. In the next section, we develop a more explicit approach to model interactions for fast inference and the generalization ability for new subprocesses.

Inference for Point Processes Inference methods for point processes are mainly based on the order statistics or likelihood function. The order statistics are often estimated nonparametrically, such as the kernel estimator (Diggle, 1985) of the intensity function. For the likelihood-based inference, we assume that one observes events X = {x i } N i=1 of the underlying spatial point process over the area R. The log-likelihood for the inhomogeneous Poisson process over space R is

log p(X|Θ) = N i=1 log(λ(x i )) - R λ(x) dx .

(3)

The integration term is the log void probability and can be viewed as a normalization term for the likelihood. For Cox processes, the likelihood is the expectation over the Poisson likelihood. It is difficult to directly integrate over the distribution of Λ. Monte Carlo methods (Adams et al., 2009) are commonly used to approximate the expectation. To improve the scalability of the expensive sampling, many methods such as variational inference (Lloyd et al., 2015) , Laplace approximation (Williams & Rasmussen, 2006) and reproducing kernel Hilbert spaces (Flaxman et al., 2017) are proposed.

Variational Autoencoder As a stochastic variational inference algorithm, VAE (Kingma & Welling, 2014) is maximizing the evidence lower bound (ELBO) of the log-likelihood function

EQUATION

The hidden variables z have a simple multivariate Gaussian prior p(z) = N (z; 0, I). The true posterior, which is often intractable as in the Cox process, is approximated via a multivariate Gaussian q φ (z|X) = N (z; µ φ (X), σ φ (X)). The KL divergence term in the ELBO can be calculated analytically. VAE uses a multilayer perceptron (MLP) to learn the mean and variance of the approximated posterior directly from the data. The most related work here is a recent VAE-based model for collaborative filtering (VAE-CF) (Liang et al., 2018) . They assume that each user is a multinomial distribution over items with the log-likelihood log p θ (X u |z u ) = N i=1 X iu log π i (z u ) for each user u. Here X u is the observed data of user clicking items, π i (z u ) is the probability that user u clicks the item i and X iu is an indicator function on whether the user u clicked the item i.

3. MULTIVARIATE SPATIAL POINT PROCESSES:

Here we consider a multivariate case of the SPP, with U interdependent univariate point processes on the sample space R. The intensity function is measured in a similar way as the univariate case via λ u (x) = lim |∆x|↓0 (E [N u (∆x)] /|∆x|) , where N u (∆x) is the number of events within a set ∆x for the subprocess u.

3.1. A NONPARAMETRIC MODEL:

The observed data of multivariate SPP include the location of N u events X u = {x u i } Nu i=1 associated with each subprocess u. For each u, the observed event locations follow a Poisson process with spatial intensity λ u (x), which is a realization of the random intensity Λ u (x). Using the nonparameteric kernel estimator, the intensity of the subprocess u is estimated by

EQUATION

Here K σ (x) is a kernel function and we usually adopt the radial basis function kernel (RBF) where K σ (x) = exp(-x 2 /2σ 2 ). We ignore the end-correction (Diggle, 1985) for now.

In real-world applications, however, one often encounters the missing data problem, where we cannot directly observe points in certain areas for some subprocesses. Instead, we seek to infer the hidden data from other fully observed subprocesses. In our model, we assume that each subprocess reflects the stochastic and heterogeneous patterns. For example, users in an e-commerce platform usually prefer different categories. As in Poisson cluster processes, this naturally leads to events clustering in specific areas. Another real-world example is the aggregation of check-in activities around the home and workplaces for social network users (Cho et al., 2011) . Note that N = U u=1 N u is the total number of events. We introduce hidden variables Y u i for each event x i = 1, ..., N and subprocess u = 1, ..., U , where Y u i = 1 if the subprocess u includes event x i and Y u i = 0 otherwise. EY u i = p u i is the probability that event x i is from the subprocess u. Then the intensity process for our multivariate SPP model is

EQUATION

for each subprocess u. This model generalizes the kernel density-based intensity to the missing data case. Similarly to the original method, it can be applied to estimate the intensity for both cluster processes such as Cox processes and repulsive ones like determinantal PPs. In order to incorporate prior information and model the data uncertainty, we adopt a variational inference approach for the hidden variables.

3.2. VARIATIONAL INFERENCE:

A major drawback of current inference methods for SPP is the introduction of a large number of parameters in the highly multivariate case. For our model, we use an amortized inference approach -VAE (Kingma & Welling, 2014) to avoid the computational complexity of directly estimating the posterior for each subprocess u.

The generative process of our model can be described as follows: For each subprocess u, it has a K-dimensional hidden variable z u with a multivariate normal prior z u ∼ N (0, I K ). Here we use a low-dimensional representation and then a nonlinear mapping f θ (z u ) = {p u i } N i=1 transforms z u so that it has the same dimension as the number of events N . Finally the spatial points of the subprocess u are sampled according to the intensity λ u (x) = N i=1 p u i K σ (x -x i ). We approximate the intractable posterior distribution of z, q(z|X) with a multivariate Gaussian N (µ φ (X), σ φ (X)). As in Liang et al. (2018) , we use MLPs to learn the nonlinear function f θ (z) with parameters θ and the mean and variance with parameters φ. The variational bound of our multivariate Cox process model is then

EQUATION

The first term in L is essentially a complete likelihood function. For each subprocess u, it has the following (expected) intensity function

EQUATION

and a Poisson process log-likelihood function from ( 8) and ( 3)

E q φ (zu|Xu) log p θ (X u |z u ) = Nu i=1 log( Nu i=1 p u i K σ (x -x u i )) - R Nu i=1 p u i K σ (x -x u i )dx. (9)

For applications without explicit spatial information, we embed each event into a latent space as a vector. First, we obtain a similarity graph for all events. Then the embedding x i of i th event in this graph is obtained via graph neural networks (GNNs) such as GraphSAGE (Hamilton et al., 2017) .

See Figure 1 for an illustration of our framework. Both x i and z u are learned jointly combining the information of item embedding and user hidden variable. 

3.3. ALTERNATIVE MODEL:

Recall that the hidden variables Y u i describe whether the event x i is from the subprocess u. By definition, we have U u=1 Y u i = 1 and U u=1 p u i = 1 for any i. During the training process, it is difficult to normalize the probability p u i over all subprocesses (have to use the full data). Moreover, this constraint leads to g uv < 1 for u = v, implying mutual-inhibition behaviors between subprocesses. Instead, we consider an alternative model where p u i is the probability that the subprocess u generates an event x i i.e. Nu i=1 p u i = 1 for each u. During the training, the total number of events N u is not viewed as a hidden variable for each subprocess. Thus the alternative model essentially normalizes λ u by a constant. With the reparameterization trick in Kingma & Welling (2014), we sample the log-likelihood function using all events within a mini user batch and compute the gradient. This approach incorporates all information about the user so that negative sampling is not needed. See Algorithm 1 for our training procedure. For the model prediction, the normalized intensity of a new subprocess can be efficiently calculated in O(N ) using the approximated posterior q φ (z|X new ) and nonlinear function f θ (z) with parameters θ, φ inferred from data. We can further reduce the computational challenge (Liang et al., 2018) for large N due to f θ (z) by discretizing the space. Now we show the equivalence of our multivariate model and the alternative one. There are two probabilities to consider. The first one is the conditional probability of observed events X u in the subprocess u with an intensity function λ u (x), given that there are N u events within the metric space R. The second one is the probability of sample X u of size N u from the normalized density h u (x) = λ u (x)/ R λ u (s)ds. For general SPPs data, we have Theorem 1. A spatial point process on a measurable set R ⊂ R n with an intensity function λ u (x) is equivalent to N u i.i.d samples within R with p.d.f h u (x) = λ u (x)/ R λ u (s)ds, given we know N u = λ u (s)ds, which is the number of points within R for the point process model.

Proof. See Section B in Appendix.:

According to this theorem (see supplementary material), we can replace the log-likelihood function (9) in the ELBO with

EQUATION

)

Here C is related to the log-likelihood on the number of events N u , which is a constant because R λ u (s)ds = N u is observed. One drawback of this approach is that, for the prediction of actual missing data, we cannot infer the number of missing points. Instead, our VAE-based model generates the normalized intensity predicting the possible locations for the missing events. We use the alternative definition of p u i from now on.

This result shows that VAE-CF is a special case of this multivariate SPP model over a discrete space X of events. In fact, VAE-CF is the alternative model with a delta function as the kernel (h u (x i ) = λ u (x i )/ X λ u (s)ds = p u i ), which is equivalent to the SPP model according to the theorem. To better model the spatial heterogeneity of events, one can replace the delta function with other kernels or use more advanced SPP intensities. We simply use a RBF kernel here, resulting in 

h u (x) = N i=1 p u i exp( x -x i 2 /2σ 2 ).

b = u∈U b X u ; forall u ∈ U b do Sample z u ∼ N (µ φ (X u ), σ φ (X u )) with reparameterization trick; Compute f θ (z u ) = {p u i } xi∈X b ; forall x ∈ X u do Compute sampled normalized intensity h u (x) ≈ xi∈X b p u i K σ (x -x i ); end

Compute noisy gradients of the ELBO L w.r.t θ and φ end Average noisy gradients over batch; Update θ and φ with the Adam optimizer (Kingma & Ba, 2015) ; end One benefit of this alternative model is its resulted consistency. The nonparametric kernel estimation for the point process intensity is unbiased. To see this, for any measurable set R, we take the expectation of the estimated intensity λ(x) over the Poisson point process distribution

EQUATION

where ρ(y) is the true intensity function. Then Eλ(x) = ρ(x) under mild conditions, e.g., a spatially continuous assumption on ρ. But it is inconsistent due to the non-vanishing variance without normalization. For our alternative model, the normalized intensity function h u (x) is still unbiased. And according to the standard theory of the multivariate kernel density estimation (KDE), the consistency of h u (x) is also guaranteed. Another benefit of using this alternative form can be seen from the cross pair-correlation function. For the alternative model, we remove the undesirable restriction of negative correlations between all users (g uv < 1 for u = v) and can incorporate more diverse relationships between users. To see this, we first consider the auto and cross pair-correlation function g uv = EΛ u Λ v /EΛ u EΛ v . For our original model, it is straightforward to prove that g uu > 1 and g uv < 1, u = v (see supplementary material). The auto pair-correlation functions show that our model is more aggregate than the simple Poisson process.

4. EXPERIMENTS:

We compare our model (with the RBF kernel, VAE-SPP) with both VAE-CF (Liang et al., 2018) and univariate spatial point process models using a standard KDE (Diggle, 1985) or TGCP (Williams & Rasmussen, 2006) as intensity functions. We adopt the experiment setting in VAE-CF. We split the data into training, validation and testing sets. For the multivariate model, the training data is used to learn the parameters θ, φ. For KDE and TGCP models, we omit the training data because different subprocesses are assumed to be independent and also because of the computational complexity of fitting a highly multivariate TGCP. We assume that only 80% of the events in the validation and test sets are observed. The remaining 20% are viewed as missing data to be inferred by different models. Hyperparameters are selected on the validation data as in Liang et al. (2018) . Finally, we compare the prediction performance of different models on the missing data given the partially-observed events. We use standard ranking losses such as NDCG@K and Recall@K defined in Appendix D.1.

4.1. MULTIVARIATE SPP ON SPATIAL DATA:

Synthetic data sets We simulate two different data sets using multiexponential and multisine models. For the multiexponential data set, we simulate 5,000 Poisson processes with Location-based Social Network. We consider the Gowalla data set (Cho et al., 2011) in New York City (NYC) and California (CA). We use a bounding box of -124.4096, 32.5343, -114.1308, 42.0095 for CA and -74.0479, 40.6829, -73.9067, 40.8820 for NYC (both from flickr1 ). Each user with at least 20 events (check-ins) is viewed as a subprocess. There are 673,183 events and 6,728 users for Gowalla-CA. We randomly select 500 users as the validation set and 500 users as the testing set. We use the remaining users for training. For Gowalla NYC, there are 86,703 events from 1,171 users. We set the size of both validation and testing sets to 100. For the spatial tessellation, we use uniform grids (32 × 32 for NYC and 64 × 64 for CA). Both our model and VAE-CF can work without grids. We further compare the performance of our model with VAE-CF by viewing each location as an item. In Table 1 , we summarize the performance of both multivariate and univariate models on the simulation data sets. It is clear that the multivariate models outperform the univariate ones. Moreover, testing on multivariate models takes less time because it only evaluates the posterior probability and intensity function. This illustrates the power of multivariate models using amortized inference. Within the multivariate models, our continuous model further improves upon the discrete VAE-CF. This is due to the fact that these simulation intensities are continuous over R. For real-world applications, the results on the location-based social network prediction and recommendation with and without grids are presented in Table 2 and 3 . We observe the same pattern in both NYC and CA. We stop using univariate models from now on due to their inferior performances, especially for collaborative filtering applications. Moreover, our model improves discrete VAE-CF regardless of the choice of spatial grids. For visualization purposes, in Figure 2 , we plot a user's check-in locations in Gowalla-NYC and intensities estimated via different methods. Comparing with VAE-CF, our model generates a continuous intensity. The univariate models overfit the training data and lead to inferior predictions of the missing data.

λ k (x) = a k e -b k x , k = 1, ...,

4.2. MULTIVARIATE SPP WITH A LATENT SPACE:

MovieLens data sets (ML-100K and ML-1M) include the movie (item) rating by users and we binarize the rating with a threshold of 4. In the spatial point process setting, we view each user as a subprocess over the latent space of item embeddings. Here the item embedding is generated via a GNN. This framework is a natural generalization of the multimodal distribution over items. The item-item graph is constructed based on item-item similarities. We use the Jaccard distance to measure the similarities between items, which are further viewed as the sampling probabilities for GNN. Currently, we only consider 1-hop connections. Both GNN and VAE are trained jointly, which is more expensive than VAE-CF but leads to better performance compared to separate training (see Appendix D). For movie recommendation tasks, we compare the discrete VAE-CF to our joint model with GNN. The results in Table 4 show again the improvement of our model over the baseline. 

5. CONCLUSION:

In this paper, we introduce a novel spatial point process model for efficient inference on the highly multivariate case. Through amortized inference, our model makes it possible to investigate correlations between a myriad of point patterns based on a large number of training data, and the theoretical analysis on the density and intensity for SPPs builds the connection between our model and VAE-CF. There are many promising directions of future works including the extension for multivariate spatiotemporal PPs (Mohler et al., 2011; Yuan et al., 2019) and using features as covariances. There are multiple ways to estimate the mean rate (h u (x)) of a spatial point process overall events, including Gaussian mixture models, Gaussian processes and flow-based models. For future work, we can investigate the connections between our model and other density-based estimations for point processes. Another interesting application is to handle real-world recommender systems via improving the joint training efficiency and comparing thoroughly with simpler algorithms as in Dacrema et al. (2019) .

Using 2), we can evaluate this conditional probability The theorem follows from G c (v) = G(v) as the probability generating functional completely determines the probability structure of the point process.

We show that (10) holds in the main paper.

Corollary 1.1.

EQUATION

)

Proof. Define λ u (x) = E q φ (zu|Xu) Λ u (x).

E q φ (zu|Xu) log p(X u |z u ) = log p(N (R) = N u )p(x u 1 , ..., x u Nu |N u ) 

C AUTO AND CROSS PAIR-CORRELATION FUNCTIONS:

We show that g u,v (x, y) > 1 for u = v and g u,v (x, y) < 1 for u = v for our orginal model.

EQUATION

We have

EQUATION

Similarly,

EQUATION

Note that U u=1 Y u i = 1 and U u=1 p u i = 1. When i = j, we have EY u i Y v j = p u i p v j for any u, v. When i = j, EY u i Y u i = p u i > (p u i ) 2 for u = v and EY u i Y v i = 0 < p u i p v i . Then it is easy to see g u,v (x, y) > 1 for u = v and g u,v (x, y) < 1 for u = v.

:

https://www.flickr.com/places/info/

https://github.com/dawenl/vae_cf

Figures and Tables:

Figure 1: Visual illustration of our spatial point process model via VAE during the training.

Training VAE SPP with stochastic gradient descent. Input: Traning subprocesses u ∈ U T with their point locations X u Result: Parameters θ and φ Initialize θ and φ randomly; while not converged do Sample a subprocesses batch U b from U T and their points X

5000, x ∈ [0, 30] as training data. Here a k and b k are uniformly sampled between [5, 10] and [0.1, 0.2] separately. 500 validation and 500 test subprocesses are generated in the same way with parameters sampled from a k and b k . The multisine data set is generated via replacing the intensity function with λ k (x) = max(a k * sin(b k x), 5) and sampling a k and b k uniformly between[5, 10] and [1, 2]  separately. Each realization of the spatial point process is discretized using a uniform grid over x with grid spacing 0.01.

Figure 2: Estimated density functions for a Gowalla user in NYC (log scale). The first row from left to right: observed check-in locations (in red), held-out check-in locations (in blue, as missing data) and the estimated intensity from VAE-SPP. The second row from left to right: the estimated intensity (or density) from VAE-CF, KDE and TGCP.

s))N (ds))|N (R) = n] = ( R λ u (s)v(s)ds λ u (s)ds ) n(15)Using 2) again and because the point process observation probability isp(ω) = p(N (R) = N u )p(x 1 , ..., x Nu |N u ) v(s)ds) n ) = exp( R λ(s)(v(s) -1)ds). (17)

function of N u .

Testing results on the simulation data sets. Both the mean and variance are percentages (same below).

Table content: <table><tr><td/><td/><td>Multiexp</td><td/><td/><td>Multisine</td><td/></tr><tr><td>Name</td><td colspan="3">NDCG@100 Recall@50 Recall@100</td><td colspan="3">NDCG@100 Recall@50 Recall@100</td></tr><tr><td>VAE-CF</td><td>6.78(0.28)</td><td>7.25(0.40)</td><td>14.5(0.52)</td><td>3.30(0.15)</td><td>2.49(0.13)</td><td>4.64(0.18)</td></tr><tr><td colspan="2">VAE-SPP 7.11(0.31)</td><td>7.34(0.40)</td><td>14.9(0.54)</td><td>3.53(0.15)</td><td>2.58(0.13)</td><td>4.90(0.18)</td></tr><tr><td>KDE</td><td>5.27(0.15)</td><td>5.85(0.12)</td><td>11.8(0.17)</td><td>3.23(0.15)</td><td>2.29(0.12)</td><td>4.55(0.27)</td></tr><tr><td>TCGP</td><td>3.11(0.14)</td><td>3.32(0.11)</td><td>6.44(0.11)</td><td>3.77(0.14)</td><td>1.88(0.11)</td><td>3.92(0.17)</td></tr></table>

Testing results on the Gowalla data sets with uniform grids.

Table content: <table><tr><td/><td/><td>CA</td><td/><td/><td>NYC</td><td/></tr><tr><td>Name</td><td colspan="3">NDCG@100 Recall@50 Recall@100</td><td colspan="3">NDCG@100 Recall@50 Recall@100</td></tr><tr><td>VAE-CF</td><td>41.8(1.5)</td><td>64.8(2.0)</td><td>70.0(2.0)</td><td>43.6(2.3)</td><td>73.9(2.9)</td><td>86.2(2.2)</td></tr><tr><td colspan="2">VAE-SPP 42.3(1.5)</td><td>65.2(2.0)</td><td>70.2(1.9)</td><td>44.8(2.4)</td><td>74.5(2.9)</td><td>86.2(2.2)</td></tr><tr><td>KDE</td><td>34.5(1.5)</td><td>59.2(2.0)</td><td>64.0(2.0)</td><td>41.2(1.5)</td><td>69.9(2.0)</td><td>83.6(2.0)</td></tr><tr><td>TCGP</td><td>31.8(1.3)</td><td>56.5(2.0)</td><td>60.9(2.0)</td><td>37.3(2.3)</td><td>59.9(3.3)</td><td>75.9(2.8)</td></tr></table>

Testing results on the Gowalla data sets without discretization.

Table content: <table><tr><td/><td/><td>CA</td><td/><td/><td>NYC</td><td/></tr><tr><td>Name</td><td colspan="3">NDCG@100 Recall@20 Recall@100</td><td colspan="3">NDCG@100 Recall@20 Recall@100</td></tr><tr><td>VAE-CF</td><td>21.3(0.77)</td><td>16.6(0.74)</td><td>32.8(0.97)</td><td>16.0(1.7)</td><td>13.2(1.7)</td><td>26.3(2.4)</td></tr><tr><td colspan="2">VAE-SPP 21.6(0.77)</td><td>17.0(0.80)</td><td>33.5(0.76)</td><td>16.1(1.7)</td><td>13.7(1.8)</td><td>27.1(2.5)</td></tr></table>

Testing results on the MovieLens data sets.

Table content: <table><tr><td/><td/><td>ML100K</td><td/><td/><td>ML1M</td><td/></tr><tr><td>Name</td><td colspan="3">NDCG@100 Recall@20 Recall@100</td><td colspan="3">NDCG@100 Recall@20 Recall@100</td></tr><tr><td>VAE-CF</td><td>40.8(2.8)</td><td>32.3(2.8)</td><td>57.6(3.3)</td><td>41.6(0.76)</td><td>33.1(0.81)</td><td>56.8(0.88)</td></tr><tr><td colspan="2">VAE-SPP 41.5(2.9)</td><td>31.3(2.7)</td><td>59.0(3.5)</td><td>42.3(0.77)</td><td>33.9(0.82)</td><td>57.6(0.88)</td></tr></table>

Peer review (iclr2020):
In this paper, the authors propose a VAE model for spatial point processes. The model generalizes the kernel density-based intensity and applies variational inference. The model is applied to synthetic datasets, a location-based social network dataset, and a recommender system dataset. The paper is well motivated and clearly written. I found the probabilistic modeling interesting. My major concern is that, with added complexity, the experimental results suggest that VAE-SPP is not significantly better than the existing VAE-CF on most tasks. For the MovieLens task, it seems that the GNN is an important component in the pipeline, but no further detail is provided about it, including the network architecture. The authors might also want to provide the definitions of NDCG@k and Recall@k, at least in the appendix. --- I appreciate the author*s detailed response and updated paper. I have changed my rating.