VARIATE SPATIAL POINT PROCESSES INTENSITIES
Year: 2020

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.

heading: 1 INTRODUCTION
text: 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-Muñ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-Muñ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.

heading: 2 PRELIMINARIES
text: 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 ⊂ Rn) 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
λ(x) = lim |∆x|↓0 E [N(∆x)] |∆x| , (1)
where ∆x is a small ball in the metric space, e.g. the Euclidean space Rn, with the centre x and measure |∆x| . The second-order intensity function is naturally defined as
λ(2)(x, y) = lim |∆x|,|∆y|↓0 E [N(∆x)N(∆y)] |∆x||∆y| , (2)
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=1Kσ(x−xi) with xi 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 = {xi}Ni=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(λ(xi))− ∫ 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
log p(X|Θ) ≥ Eqφ(z|X)[log(pθ(X|z)]−KL(qφ(z|X)|p(z)). (4) 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θ(Xu|zu) = ∑N i=1Xiu log πi(zu) for each user u. Here Xu is the observed data of user clicking items, πi(zu) is the probability that user u clicks the item i and Xiu is an indicator function on whether the user u clicked the item i.

heading: 3 MULTIVARIATE SPATIAL POINT PROCESSES
text: 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 [Nu(∆x)] /|∆x|) , where Nu(∆x) is the number of events within a set ∆x for the subprocess u.

heading: 3.1 A NONPARAMETRIC MODEL
text: The observed data of multivariate SPP include the location of Nu events Xu = {xui } 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
λu(x) = Nu∑ i=1 Kσ(x− xui ). (5)
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=1Nu is the total number of events. We introduce hidden variables Y u i for each event xi = 1, ..., N and subprocess u = 1, ..., U , where Y ui = 1 if the subprocess u includes event xi and Y ui = 0 otherwise. EY ui = pui is the probability that event xi is from the subprocess u. Then the intensity process for our multivariate SPP model is
Λu(x) = N∑ i=1 Y ui Kσ(x− xi), (6)
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.

heading: 3.2 VARIATIONAL INFERENCE
text: 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 zu with a multivariate normal prior zu ∼ N (0, IK). Here we use a low-dimensional representation and then a nonlinear mapping fθ(zu) = {pui }Ni=1 transforms zu 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 iKσ(x−xi). 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
log p(Xu|Θ) ≥ Eqφ(zu|Xu)[log(pθ(Xu|zu)]−KL(qφ(zu|Xu)|p(zu)) = L. (7)
The first term in L is essentially a complete likelihood function. For each subprocess u, it has the following (expected) intensity function
Eqφ(zu|Xu)Λu(x) = Nu∑ i=1 puiKσ(x− xui ) (8)
and a Poisson process log-likelihood function from (8) and (3)
Eqφ(zu|Xu) log pθ(Xu|zu) = Nu∑ i=1 log( Nu∑ i=1 puiKσ(x− xui ))− ∫ R Nu∑ i=1 puiKσ(x− xui )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 xi of ith 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 xi and zu are learned jointly combining the information of item embedding and user hidden variable.

heading: 3.3 ALTERNATIVE MODEL
text: Recall that the hidden variables Y ui describe whether the event xi 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 pui over all subprocesses (have to use the full data). Moreover, this constraint leads to guv < 1 for u 6= v, implying mutual-inhibition behaviors between subprocesses. Instead, we consider an alternative model where pui is the probability that the subprocess u generates an event xi i.e. ∑Nu i=1 p u i = 1 for each u. During the training, the total number of events Nu 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|Xnew) 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 Xu in the subprocess u with an intensity function λu(x), given that there are Nu events within the metric space R. The second one is the probability of sample Xu of size Nu from the normalized density hu(x) = λu(x)/ ∫ R λu(s)ds. For general SPPs data, we have
Theorem 1. A spatial point process on a measurable set R ⊂ Rn with an intensity function λu(x) is equivalent to Nu i.i.d samples within R with p.d.f hu(x) = λu(x)/ ∫ R λu(s)ds, given we know
Nu = ∫ λ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
Eqφ(zu|Xu) log pθ(Xu|zu) = Nu∑ i=1 log(hu(x u i )) + C. (10)
Here C is related to the log-likelihood on the number of events Nu, which is a constant because∫ R λu(s)ds = Nu 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 pui 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 (hu(xi) = λu(xi)/ ∫ 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 hu(x) = ∑N i=1 p u i exp(‖x− xi‖2/2σ2).
Algorithm 1: Training VAE SPP with stochastic gradient descent. Input: Traning subprocesses u ∈ UT with their point locations Xu Result: Parameters θ and φ Initialize θ and φ randomly; while not converged do
Sample a subprocesses batch Ub from UT and their points Xb = ⋃ u∈Ub Xu ; forall u ∈ Ub do Sample zu ∼ N (µφ(Xu), σφ(Xu)) with reparameterization trick; Compute fθ(zu) = {pui }xi∈Xb ; forall x ∈ Xu do
Compute sampled normalized intensity hu(x) ≈ ∑ xi∈Xb p u iKσ(x− xi);
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
E ∫ R λ(x)dx = ∫ R E N∑ i=1 Kσ(x− xi)dx = ∫ R ∫ R Kσ(x− y)ρ(y)dydx = ∫ R ρ(y)dy, (11)
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 hu(x) is still unbiased. And according to the standard theory of the multivariate kernel density estimation (KDE), the consistency of hu(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 (guv < 1 for u 6= v) and can incorporate more diverse relationships between users. To see this, we first consider the auto and cross pair-correlation function guv = EΛuΛv/EΛuEΛv . For our original model, it is straightforward to prove that guu > 1 and guv < 1, u 6= v (see supplementary material). The auto pair-correlation functions show that our model is more aggregate than the simple Poisson process.

heading: 4 EXPERIMENTS
text: 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.

heading: 4.1 MULTIVARIATE SPP ON SPATIAL DATA
text: 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 λk(x) = ake −bkx, k = 1, ..., 5000, x ∈ [0, 30] as training data. Here ak and bk 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 ak and bk. The multisine data set is generated via replacing the intensity function with λk(x) = max(ak ∗ sin(bkx), 5) and sampling ak and bk 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.
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.
1https://www.flickr.com/places/info/
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.

heading: 4.2 MULTIVARIATE SPP WITH A LATENT SPACE
text: 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.

heading: 5 CONCLUSION
text: 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 VAECF. 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 (hu(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).

heading: ACKNOWLEDGMENTS
text: We would like to thank the comments from Frederic P. Schoenberg and George Mohler. Andrea L. Bertozzi and Baichuan Yuan want to thank the support of NIJ fellowship 2018-R2-CX-0013 and NSF DMS-1737770.

heading: A TABLE OF NOTATIONS

heading: B PROOF OF THEOREM 1
text: Proof. We define our model as a point process on R with the intensity function λu(x). The alternative model is Nu i.i.d samples within R with p.d.f hu(x), given that we know Nu =∫ λu(s)ds is the number of points within the point process model.
1) Our model has the following probability generating functional G(v) = exp(− ∫
Rd [1− v(x)] Λ(dx)) (12)
2) Given Nu,
p(x1, ..., xNu |Nu) = Nu∏ i=1 hu(xi) (13)
3) Our alternative model (a counting r.v. N(x) with locations according to hu(x)) has the following characteristic functional
Gc(v) = ∞∑ n=0 p(N(R) = n)E[exp( ∫ R log(v(s))N(ds))|N(R) = n] (14)
Using 2), we can evaluate this conditional probability E[exp( ∫ R log(v(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 is
p(ω) = p(N(R) = Nu)p(x1, ..., xNu |Nu) = 1
n! [ n∏ i=1 λu(xi)] exp(− ∫ R λ(x)dx), (16)
we have Gc(v) = exp(− ∫ R λ(x)dx)(1 + ∞∑ n=1 1 n! ( ∫ R λ(s)v(s)ds)n) = exp( ∫ R λ(s)(v(s)− 1)ds). (17)
The theorem follows from Gc(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.
Eqφ(zu|Xu) log p(Xu|zu) = Nu∑ i=1 log(hu(x u i )) + C. (18)
Proof. Define λu(x) = Eqφ(zu|Xu)Λu(x). Eqφ(zu|Xu) log p(Xu|zu) = log ( p(N(R) = Nu)p(x u 1 , ..., x u Nu |Nu) ) (19)
= Nu∑ i=1 log(hu(x u i )) + log(p(N(R) = Nu)) (20)
log(p(N(R) = Nu)) = n log( ∫ R λ(x)dx)− log(n!)− ∫ R λ(x)dx (21)
is only a function of Nu.

heading: C AUTO AND CROSS PAIR-CORRELATION FUNCTIONS
text: We show that gu,v(x, y) > 1 for u = v and gu,v(x, y) < 1 for u 6= v for our orginal model.
gu,v(x, y) = EΛu(x)Λv(x) EΛu(x)EΛv(x)
(22)
We have
EΛu(x)Λv(x) = E( N∑ i=1 Y ui Kh(x− xi))( N∑ j=1 Y vj Kh(y − xj)) (23)
= N∑ i=1 N∑ j=1 EY ui Y vj Kh(x− xi)Kh(y − xj). (24)
Similarly,
EΛu(x)EΛv(x) = (E N∑ i=1 Y ui Kh(x− xi))(E N∑ j=1 Y vj Kh(y − xj)) (25)
= N∑ i=1 N∑ j=1 pui p v jKh(x− xi)Kh(y − xj). (26)
Note that ∑U u=1 Y u i = 1 and ∑U u=1 p u i = 1. When i 6= j, we have EY ui Y vj = pui pvj for any u, v. When i = j, EY ui Y ui = pui > (pui )2 for u = v and EY ui Y vi = 0 < pui pvi . Then it is easy to see gu,v(x, y) > 1 for u = v and gu,v(x, y) < 1 for u 6= v.

heading: D MORE ON EXPERIMENTS
text: D.1 METRICS DEFINITION
The ranking performance is evaluated through recall at K (Recall@K) and normalized discounted cumulative gain at K (NDCG@K). In our VAE-SPP model, the predicted rank of the held-out items Iu for each user u are obtained from sorting the intensity function λu(x).
Here we keep the definition in Liang et al., (2018). Recall@K is defined as Recall@K = ∑K i=1 I(ri ∈ Iu) |Iu| . (27)
NDCG@K is calculated by normalizing discounted cumulative gain (DCG@K) with ideal DCG@K (IDCG@K). The definition are as follows
DCG@K = K∑ i=1 2I(ri∈Iu) − 1 log2(i+ 1) , NDCG@K = DCG@K IDCG@K , (28)
where I is the indicator function and ri is the ith item among held-out items; IDCG@K is the ideal DCG@K when the ranked list is perfectly ranked.
D.2 HYPERPARAMETERS
We implement our models in Tensorflow based on VAE-CF2. We keep the same MLP architectures and hyperparameters for both of them. We use β-VAE with as suggested. The only additional hyperparameters for our model is the σ2 in the kernel function, which is determined using grid search on the validation set. We conducted the experiments on a single GTX 1080 TI 11GB GPU.
For simulation data, we train both models for 200 epochs using Adam optimizer with β = 0.2, lr = 5 × 10−5. We use mini-batches of size 20. Our architectures consist of a one layer MLP with K = 50. For VAE-SPP, σ2 = 0.001. For Gowalla-NYC data, we train both models for 200 epochs using Adam optimizer with β = 0.2, lr = 5 × 10−4. We use mini-batches of size 20. Our architectures consist of a one layer MLP with K = 50. For VAE-SPP, σ2 = 1 × 10−5. For Gowalla-LA data, we train both models for 200 epochs using Adam optimizer with β = 0.2, lr = 1× 10−3. We use mini-batches of size 20. Our architectures consist of a one layer MLP with K = 50. For VAE-SPP, σ2 = 0.001. For ML-1M data, we train both models for 100 epochs using Adam optimizer with β = 0.2, lr = 1 × 10−3. We use mini-batches of size 5. Our architectures consist of a one layer MLP with K = 200. For VAE-SPP, σ2 = 1 × 10−5. For ML-100K data, we train both models for 100 epochs using Adam optimizer with β = 0.2, lr = 1 × 10−3. We use mini-batches of size 5. Our architectures consist of a one layer MLP with K = 200. For VAESPP, σ2 = 1 × 10−5. The one-layer GNN in ML data is trained using GraphSAGE, for which the embedding dimension is 32 and the number of neighborhood is 10 for item and 5 for user. The graph is consist of the edges between users and items as well as the edges between items based on their Jaccard similarity.
We use python statsmodel for the KDE and GPy for TGCP. Bandwidth for KDE is selected automatically. The hyperparameters for TGCP are determined with a grid search on the validation set. For simulation data sets, we set an RBF kernel with variance=1, lengthscale=0.1 for TGCP. For Gowalla-CA data sets, we set an Matern32 kernel with variance=1e-3, lengthscale=0.1 for TGCP. For Gowalla-NYC data sets, we set an Matern32 kernel with variance=1e-4, lengthscale=0.01 for TGCP.
D.3 ADDITIONAL EXPERIMENTS
On the training of VAE and GNN, we tried different training settings (separately or jointly) and choose to train them jointly. We also tested the point estimate version of the VAE-CF called DAECF (Mult-DAE in Liang et al., (2018), with the same setting), which can improve the result under
2https://github.com/dawenl/vae_cf
certain metrics. One can easily extend our work to a DAE-SPP to obtain a point estimation for the SPP intensity.
We also did experiments on the MLPs for VAE. For Movie Lens 1M, the larger network in VAE-CF leads to a 40.3 NDCG@100 for VAE-CF and 41.9 for VAE-SPP. As a result, we use the smaller one instead.

Review summaries:
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.