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):
The work describes an application of a spatial point process for solving problems with missing data. The authors introduce a novel method based on a non-parametric definition of point intensities for the multivariate case. The method incorporates VAE framework to effectively handle missing points via smooth intensity estimation and enjoys amortized inference for efficient computations and quick prediction generation. Using a sequence of mild assumptions, the authors show connection to a popular VAE-based collaborative filtering model, which turns out to be a special case of their approach. This is a rigorous study providing theoretically justified evidence on the effectiveness of the proposed approach. Apart from the issues in the last part of experiments with classical collaborative filtering task (which will be detailed below), the work presents a solid research. I would therefore vote for accepting it. The text is well structured, and all key points are clearly explained. The problem solved by the authors is well described, and the motivation for this work is convincing. The way point process theory is applied constitutes a rigorous probabilistic approach. The authors convincingly justify the need for all approximations and simplifications made in the model. One of key results making the entire model feasible is supported by the corresponding theorem proved by the authors. I haven’t carefully verified all the derivations, though. My major concern is related to the last part with experiments on the Movielens data. As the authors state, “applications without explicit spatial information, we embed each event into a latent space as a vector.” “No spatial information” is exactly the case with the standard collaborative filtering task, which the authors attempt to solve. This leads to an introduction of an additional model like GNN, which is unrelated to the main approach. As GNN is involved it’s not immediately obvious that the improvement over standard VAE architecture, observed in the experiments on ML-100K and ML-1M, is due to a better point process modelling. No evidence is provided to argue that this is not simply due to a good compression or a good data preprocessing achieved by a GNN architecture itself. Therefore, the results on a pure recommender systems part are not convincing. What would happen if GNN was trained and fed into another (simpler) algorithm? Maybe a simple KNN based algorithm would produce comparable or even better results? As indicated by the work of [Dacrema, Cremonesi, Jannach 2019] on “A Worrying Analysis of Recent Neural Recommendation Approaches”, VAE-CF (along with several other recently proposed neural network-based methods) is inferior to even properly tuned kNN-based models. I would not be surprised, if a kNN model trained on GNN output would produce even better results than the proposed VAE-SPP. Another related question is how incorporating GNN affects the training time? Is it comparable to that of VAE-CF or is it much worse? Computational performance is an important part in making practical decisions and should be also considered. Furthermore, both ML datasets used for tests are too small and not very representative to make any generalized conclusions. Even on a larger ML-20M dataset an optimal SVD-based model can be trained within several minutes on a standard CPU on a laptop (according to my experiments, VAE-CF would take at least twice longer on Tesla K80). Therefore, it can hardly be considered a realistic example. In practice, there could be millions and hundreds of millions of items. The authors even mention it in the in the introduction, using it as a vehicle to motivate their approach. However, computing similarities between that many entities can be a laborious task on its own, which adds an extra layer of complexity and again is not directly related to the main approach. It can easily become a bottleneck or make further computations inefficient. More efficient similarity computations may in contrast reduce the resulting accuracy. The issue can get even worse, because, unlike classical MF methods, there’s still no proper support for sparse operations in NN frameworks. In the VAE framework it means that, during the training, user batches will be converted into dense arrays and may become inefficient to work with in terms of memory and CPU utilization (a few non-zero entries vs. hundreds of millions of explicitly stored zeroes). In spite of all this, I’d also suggest rephrasing “We validate these bene?ts through extensive experiments” as it sounds a bit exaggerated (if we are considering real recommender systems applications). I agree that the proposed approach is potentially applicable in real cases for recommender systems, however, there’s still not enough evidence for this. In fact, I don’t even think that completely removing the part with ML-100K and ML-1M datasets would make the whole work any worse. Clearly stating the region of applicability of the proposed approach would be enough. Right now some statements in this section in contrast are raising concerns rather than convincing the reader. The wording should be at least changed, so that readers do not get an impression that the case with classical CF task is solved purely by the proposed VAE-SPP approach. Other remarks to help improve the text: 1) “… points are more likely to … form clusters than the simple Poisson process …” the sentence seems to be inconsistent. 2) “The generative process of our model can be described as follow:” -> … as follows: 3) Page 4, last paragraph, line 6 – shouldn’t the upper bound for summation be N_u instead of just N? References: Dacrema, Maurizio Ferrari, Paolo Cremonesi, and Dietmar Jannach. *Are we really making much progress? A worrying analysis of recent neural recommendation approaches.* In Proceedings of the 13th ACM Conference on Recommender Systems, pp. 101-109. ACM, 2019.