Automatic Integration for Spatiotemporal Neural Point Processes
Authors: Zihao Zhou
Rose Yu
Venue: NeurIPS 2023 poster (2023)
Keywords: spatiotemporal modeling
neural point processes
integration method
TL;DR: This paper introduces Auto-STPP, a novel paradigm that extends the AutoInt approach to 3D STPP, and showcases its significant advantages in recovering complex intensity functions from irregular spatiotemporal events.

Abstract: Learning continuous-time point processes is essential to many discrete event forecasting tasks. However, integration poses a major challenge, particularly for spatiotemporal point processes (STPPs), as it involves calculating the likelihood through triple integrals over space and time. Existing methods for integrating STPP either assume a parametric form of the intensity function, which lacks flexibility; or approximating the intensity with Monte Carlo sampling, which introduces numerical errors. Recent work by Omi et al. proposes a dual network approach for efficient integration of flexible intensity function. However, their method only focuses on the 1D temporal point process. In this paper, we introduce a novel paradigm: `Auto-STPP` (Automatic Integration for Spatiotemporal Neural Point Processes) that extends the dual network approach to 3D STPP. While previous work provides a foundation, its direct extension overly restricts the intensity function and leads to computational challenges. In response, we introduce a decomposable parametrization for the integral network using ProdNet. This approach, leveraging the product of simplified univariate graphs, effectively sidesteps the computational complexities inherent in multivariate computational graphs. We prove the consistency of `Auto-STPP` and validate it on synthetic data and benchmark real-world datasets. `Auto-STPP` shows a significant advantage in recovering complex intensity functions from irregular spatiotemporal events, particularly when the intensity is sharply localized. Our code is open-source at https://github.com/Rose-STL-Lab/AutoSTPP.

Automatic Integration for Spatiotemporal Neural Point Processes
Zihao Zhou
Rose Yu
Department of Computer Science University California, San Diego
Department of Computer Science University of California, San Diego La Jolla, CA 92092 roseyu@ucsd.
La Jolla, CA 92092 ziz244@uesd. edu
Abstract
Learning continuous-time point processes is essential many discrete event fore- casting tasks. However, integration poses a major challenge, particularly for spatiotemporal point processes (STPPs), as it involves calculating the likelihood through triple integrals over space and time. Existing methods for integrating STPP either assume a parametric form of the intensity function, which lacks flexibility; or approximating the intensity with Monte Carlo sampling, which in- troduces numerical errors. Recent work by Omi al.] [2019] proposes dual network approach for efficient integration of flexible intensity function. However, their method only focuses on the ID temporal point process. In this paper, we introduce a novel paradigm: AutoSTPP (Automatic Integration for Spatiotempo- Neural Point Processes) that extends the dual network approach 3D STPP. While previous work provides foundation, its direct extension overly restricts the intensity function and leads to computational challenges. In response, we introduce decomposable parametrization for the integral network using ProdNet. This approach, leveraging the product of simplified univariate graphs, effectively sidesteps the computational complexities inherent in multivariate computational graphs. We prove the consistency AutoSTPP and validate it on synthetic data and benchmark real-world datasets. AutoSTPP shows significant advantage in recovering complex intensity functions from irregular spatiotemporal events, particularly when the intensity sharply localized. Our code is open-source https::/ /ttps:/://ttttt../// cwwccomm//ooo/ Rose- -TTL-Lab/ Laaassa-tTLLabA
Introduction
Spatiotemporal point process (STPP) DDleyy and Vere-Jones/ 2007 Reinhart 2018 is a continuous
time stochastic process for modeling irregularly sampled events space and time. STPPs are par- ticularly well-suited for modeling epidemic outbreaks, ride-sharing trips, and earthquake occurrences. A central concept in STPP is the intensity function, which captures the expected rates of events occur- rence. Specifically, given the event sequence over space and time He {($1, .)...., (Sn, fn )t, St. the joint fog-likelihood of the observed events is:
-J X" r)dudr
logp(HH) Σ log (s,,
X" is the optimal intensity, S the spatial domain, u and 7 the space and time, and i the time range. Learning STPP requires multivariate integrals of the intensity function, which is numerically challeng- ing. Traditional methods often assume a parametric form of the intensity function, such as the integral
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

can have closed-fome solution Daley and Vere-Jones 2007]. But this also limits the expressive power of the model describing complex spatiotemporal patterns.
Others propose to parameterize the model using neural ODE models Chen et al. (020) and Monte Carlo sampling, but their computation is costly for high-dimensional functions. Recently, work by Zhou et al. [2022] proposes a nonparametric STPP approach. They use kernel density estimation for the intensity and model the parameters of the kernels with deep generative model. Still, their method heavily depends on the number of background data points chosen as a hyper-parameter. Too few background points cause the intensity function to be an inflexible Gaussian mixture, while too many background points may cause overfitting on event arrival.
To reduce computational cost while maintaining expressive power, we propose novel automatic integration scheme based on dual networks to learn STPPs efficiently. Our framework models intensity function the sum of background and influence functions. Instead of relying on closed- form integration or Monte Carlo sampling, we directly approximate the integral of the influence function with a deep neural network (DNN). Taking the partial derivative of DNN results in a new computational graph that shares the same parameters; see Figure [
First, we construct integral network whose derivative is the intensity function. Then, we train the network parameters by maximizing the data likelihood formulated in terms of the linear combinations of the integral networks. Finally, we reassemble the parameters of the integral network to obtain the intensity. This approach leads to the exact intensity function and its an- tiderivative without restricting its parametric forms, thereby making the integration process "automatic".
J.I. Aln d/dnet
W2
W2
W
Derivative Network
Integral Network
A(s,t)
W
Wi
W
Our approach bears a resemblance with the fully NN approach by Omi et al. [019] for ID tempo- ral point processes. There, automatic integration can be easily implemented by imposing mono- tonicity constraints on the integral network (Lin- dell et al. 202T Li et al. 2019]. However, due to triple integration in STPP, imposing mono- tonicity constraints significantly hurdles the ex- pressivity of the integral network, leading to inaccurate intensity. As the experiments show, extending Fully NN to 3D cannot learn complex spatiotemporal intensity functions. Instead, we propose a decomposable parametrization for the integral network that bypasses this restriction.
W
W
X
x
Time
Figure 1: Illustration of AutoSTPP. W denotes the linear layer's 's weight. is the nonlinear activation function. Left shows the intensity network approx- imating A(s, t), and right is the integral network that computes J J, A. The two networks share the same parameters.
experimentt
Our approach can efficiently compute the exact likelihood of any continuous influence function. We validate our approach using synthetic spatiotemporal point processes with complex intensity functions. We also demonstrate the superior performance on several real-world discrete event datasets. Compared FullyNN by Omi et al. (2099) our approach presents a more general higher-order automatic integration scheme and is more effective in learning complex intensity functions. Also, the probability density Omi al. (2019) is ill-defined as its intensity integral does not diverge to infinity [Shchur et al. 2019]. We fix the issue by adding a constant background intensity μ.
To summarize, our contributions include:
We propose the first deep learning framework AutoSTPP to speed up spatiotemporal point process learning with automatic integration. We use dual networks and enforce the non- negativity of the intensity via a monotone integral network.
We show that our automatic integration scheme empirically learns intensity functions more accurately than other integration approaches.
We prove that the derivative network of AutoSTPP a universal approximator of continuous functions and, therefore, is a consistent estimator under mild assumptions.

---

We demonstrate that AutoSTPP can recover complex influence functions from synthetic data, enjoys high training efficiency and model interpretability, and outperforms the state-of- the-art methods on benchmark real-world data.
2 Related Work
Fitting traditional STPP, such as the spatiotemporal Hawkes process, to data points with parametric models can perform poorly if the model is misspecified. To address this issue, statisticians have extensively studied semi- and non-parametric inference for STPP. Early works like Brix and Moller [200]] Brix and Diggle [2001) usually rely on Log- Gaussian Cox processes as a backbone and Epanechnikov kernels as estimators of the pair correlation function. [Adams et al.) [2099] propose a non -on-marmmttticc approach that allows the generation of exact Poisson data from a random intensity drawn from a Gaussian Process, thus avoiding fiite-dimensional approximation. These Cox process models are scalable but assume continuous intensity change over time.
Parametrizing Point Process.
Recently, neural point processes (NPPs) that combine point processes with neural networks have received considerable attention Yan et al. 2018 Upadhhay et al. 2018 Huang et al. 2019 Shang and Sun 2019 Zhang al. 2020 Under this framework, models focus more on approximating discrete set of intensities before and after each event. The continuous intensity comes from interpolating the intensities between discrete events. For example, (Du al.) 2016) uses RNN generate intensities after each event. [Mei and Eisner, 2016) proposes a novel RNN architecture that generates intensities at both ends of each inter-event interval. Other works consider alternative training schema: Xiao iet al. [2017) used Wasserstein distance, Guo et al. [2018) introduced noise-contrastive estimation, and Li et al. [2018) leveraged reinforcement learning. While these NPP models are more expressive than the traditional point process models, they still assume simple (continuous, usually monotonous) inter-event intensity changes and only focus on temporal point processes.
Neural STPP [Chen et al., 2020 Zhou et al. 2022] further generalizes NPP to spatiotemporal data. They use non-negative activation function to map the hidden states to a scalar, i.e., the
temporal intensity immediately after an event, and a conditional spatial distribution. The change of intensity between events is represented by decay function or a Neural ODE. The conditional spatial distribution is represented by a kernel mixture or a normalizing flow. Nevertheless, all models assume continuous transformation of the intensity function and have limited expressivity.
In the context of temporal point processes (TPP), closely related approaches are Omi et al. 2019
and Zhou and Yu [203]]. Both propose using Deep Neural Network (DNN) parameterize the integral of an intensity function. The work by Omi et al. 2019 offers a more flexible formulation yet does not incorporate any specific prior assumptions about the form of the intensity changes. On the other hand, the approach by Zhou and Yu (2023 is capable of capturing more sophisticated influence functions, and is this work that our research builds upon. However, both of these studies focus on the easier problem of learning the derivative with respect to time alone, thereby neglecting the rich amount of other features that may associated with the timestamps.
Integration Methods. Integration methods are largely ignored in NPP literature but are central to a model's 's ability to capture the complex dynamics of a system. Existing works either use an intensity function with an elementary integral [Du et al., [0166] or Monte Carlo integration [Mei and Eisner 2016). However, we will see in the experiment section that the choice of integration method has non-trivial effect on the model performance.
Integration is generally more complicated than differentiation, which can be mechanically solved using the chain rule. Most integration rules, e.g., integration by parts and change of variables, transform antiderivative to another that is not necessarily casier. Elementary antiderivative only exists for a small set of functions, but not for simple composite functions such as exp (x Dunham 2018). The Risch algorithm can determine such elementary antiderivative [Risch] 1969 1970 but has
never been fully implemented due to its complexity. The most commonly used integration methods still numerical: Newton-Cotes Methods, Romberg Integration, Quadrature, and Monte Carlo integration [Davis and Rabinowitz 2007].
Several recent works leverage automatic differentiation to speedup integration, a paradigm known as Automatic Integration (AutoInt). Liu 2020) proposes integrating the Taylor polynomial using the derivatives from Automatic Differentiation (AutoDiff). It requires partitioning of the integral limits
3

---

and choosing the order of Taylor approximation. Though it makes use of AutoDiff, the integration procedure involves a trade-off between runtime and accuracy and is numerical in nature. Li al.) (2199] and Lindell et al. [2021] proposed dual network approach, which we will discuss in detail in Section
This approach guarantees closed-form integral and efficient.
Methodology
We first review the background of Spatiotemporal Point Processes. Then, we introduce the Autolnt technique, which is interpretable and flexible. Lastly, we consider applying the Autolnt technique to the 3D integration in the spatiotemporal point process.
3.1 Spatiotemporal Point Process
spatiotemporal point process (STPP) generalizes TPP to model the number of events ( x (a, b)) that occurred in the Cartesian product of the spatial domain S S Hd (d is the spatial dimensionality) and the time interval (a, b]. It is characterized by non-negative space-time intensity function given the event history He = ((s),.1)..) (Sn, tn)
(B(( As) x (1,1+ A())
X (s = lim
As-00
B(s. As)
which the probability of finding an event in an infinitesimal time interval (t, + At) and an infinitesimal spatial ball S Bss, As) centered at location s. Alternatively, an STPP be seen a TPP with a conditional spatial distribution (s(t),, such that (s, t) = (s)f"
3.2 AutoInt Point Process
Consider the following NPP model that generalizes the spatiotemporal Hawkes process:
ts, -
Σ f - t., H(t))).
(3)
Here μ is the base intensity. fø is positive scalar function that takes space, time, and representations of event history H(s;, f4) as inputs. Each f is approximated by a DNN.
The two main benefits of such design are flexibility and interpretability. fo is neural network that can model complex inter-event change. The additive form allows the decomposition the intensity function for event influence analysis. We extend this model into the spatiotemporal domain,
3.3 Automatic Integration (AutoInt)
One advantage of the neural STPP model in Equation 3 is that we can instantiate automatic integration (AutoInt) and calculate the volume integral fed fo h) := (b h) - Fa h), where h is the latent representation of the event history until generated by a deep sequence model. AutoInt first constructs the integral network Fg and then reorganizes the computational graph of Fo to form the integrant, the derivative network fo. The two networks thus share the same of parameters
Specifically, let x = s @ f ① h, we approximate the integral of the intensity function with DNN of the following form:
Fx(x) = W,, ... (W (W2o(W11))
where n is the number of layers, W Rr + N denotes the weight of the k-th linear layer of the neural network and denotes the elementwise nonlinearity. M₁, and N₂ are the input and output dimension for the k-th layer. Hence, the set of parameters in this neural network is 0 (W M Vk).
The derivative network fe is the partial derivative of the integral network Fe. As long as the activation function is differentiable everywhere, one can compute the intensity recursively.
OF.
fx(x) := at (x) = W Wn''(Wn (Wn-jo w (Wix)) ... W20 (Wix) o WT1.

---

where o indicates the Hadamard product, and W11 is the first column of W1, i.e.
Wi = W11
W12
WI.M.]
Computing (x) involves many repeated operations. For example, the result of W 1i is used for We have developed a program that harnesses the power of dynamical programming compute derivatives efficiently with AutoDiff. See Appendix D for the detailed algorithm.
compute both o (W,x) and o' o'W xx), see Figure [
3.4 AutoInt Point Processes as Consistent Estimators
We show that the universal approximation theorem (UAT) holds for derivative networks. This theorem sigmifies that, given sufficient number of hidden units, derivative networks can theoretically approximate any continuous functions, matter how complex. Therefore, using derivative networks does not limit the range of influence functions that can be approximated.
Proposition 3.1 (Universal Approximation Theorem for Derivative Networks). The set of derivative networks corresponding to two-layer feedforward inlegral networks is dense in C(R) with respect to the uniform norm.
For detailed proof, Appendix With UAT, it is clear that under some mild assumptions, Autolnt Point Processes are consistent estimators of point processes that take the form of Equation
Proposition 3.2 (Consistency of AutoInt Point Process). Under the assumption that the ground truth point process is stationary, ergodic, absolutely continuous and predictable, if the ground truth influence function is truncated (i.e., 3C, f(t) = VI > c), the maximun likelihood estimator fo converges to the true influence function f in probability as T + co.
Our model belongs to the class of linear self-excitation processes, whose maximum likelihood estimator properties were analyzed by Ogata et al. [1978]. Under the assumptions above, two conditions are needed for the proof of consistency:
Assumption 3.3. (Consistency Conditions) For any 0 θ there is neighbourhood U of such that
1. supereu (t, Xz. (1,ww) - 0 in probability as t +
2. supwet |log Ag- (t,w)l has, for some >> 0, finite (2 + a) th moment uniform bounded with respect
The first condition is satisfied by UAT. The second condition depends the rate of decrease of the influence tail and satisfied by truncation. In our experiments, we truncated the history by only including the influences of the previous 20 events. If the ground truth influence function decays over the entire time domain, our estimator may exhibit negligible bias.
3.5
3D Automatic Integration
AutoInt gives an efficient way to calculate line integral over one axis. However, for spatiotemporal point process models, we need calculate the triple integral of the intensity function over space and time. Since we cannot evaluate the integral network with an input of infinity, we assume the spatial domain be rectangle, such that the triple integral is over cuboid. We then convert the triple integral to line integrals using Divergence and Green's 's theorem MMrsden and Tromba 2003).
Define (xr, 3) and == z, we model the spatiotemporal influence fă (zz y,
OP aQ dx dy + Oz
aR
(s,
(4)
The Divergence theorem relates volume and surface integrals. It states that for any P(a,,,;), Qr,, y, =), R(x, y, that is differentiable over the cuboid Ω have
OP aQ
aR du $ Pdydz Qdadr Rdzdy, }ss
n dr By
2z
where E are the six rectangles that enclose Q.
5

---

Ground Truth
Autolnt ProdNet
Constrained Triple Autoint
z=1.14
Figure 2: Comparison of fitting nonnegative function [(x y, 2) sin(x) cos(g) (in(z) + 1. Our proposed Autolnt ProdNet approach can well approximate the ground truth function. In contrast, imposing the constraint through activation function with nonnegative triple derivative "Constrained Triple Autolnt" fails due to overly stringent constraint.
Using R as an example, we begin with two neural network approximators, La and MR. We initialize the two approximators as first-order derivative networks, whose corresponding integral networks are f LRdz and MRdy. We use the zy partial derivatives of the two networks to
0MR ÖLR model R. such that R == dr dy
òR 2²M 0LR
Then the 2 partial derivative
is
dz Ozdz
dydz
OR az can be exactly evaluated as the third derivatives of two integral networks. In the manner, parametrized by six networks: LR, LQ. MR, MQ, and Mp.
8 Mndy 83 Lndd
Ordydz drdydz We can see that
model P and such that the intensity is
We neural network pairs M to parametrize each intensity term in Equation 4 according to Green's 's theorem. states that for any L and M differentiable over the rectangle D,
aL (
0M By
drdy ff. (Ldr Mdy),
dr
where L+ is the counterclockwise path that encloses D. J Rdxdy is then & Lædr MRdy, and can be exactly evaluated using the integral networks.
In practice, observe that the six networks' parametrization of intensity be simplified. LR. LQ- Lp can share the same set weights, and similarly for MR. Mq. Mp. The influence is essentially parametrized by two different neural networks, L and M. That is,
fo(1, h) = 883 My 633 S Ldx
+ (S' Mdz
8th J'L
) Mdx 8' Ldz
ōxōzōz )+
5xdzōz 8zōyōz
ōrōyōz
3.6
Imposing the 3D Non-negativity Constraint
For the NPP model in Equation 3] the derivative network fo needs be non negaativeee.. Imposing the non-negativity constraint for 3D AutoInt is challenging task. It implies that the integral network Fø always has a nonnegative triple derivative.
simple approach is to apply an activation function with a nonnegative triple derivative. An integral network that uses this activation and has nonnegative linear layer weights satisfies the condition. We call this approach "Constrained Triple AutoInt". However, the output of an integral network can grow very quickly with large input, and the gradients are likely explode during training. Moreover, the
FF to be positive. But an activation
non-negative constraint on the influence function only requires
Osot function with a nonnegative triple derivative would also enforce other partial derivatives be positive. F 0F, can both Os0s Otät
Such constraints are overly restrictive for STPPs, whose partial derivatives be negative when the intensity is well-defined.
and
ProdNet. We propose different solution to enforce the 3D non-negative constraint called ProdNet. Specifically, we decompose the influence function fo (81, 82, 1) RR3 -> as (s) (x) (1). the product of three R - IR AutoInt derivative networks. The triple antiderivative of the influence

---

SpatialTemport Temporal Hawkes (STH)
Spatial Temporal Self-Correcting (STSC)
10
10
Auto
20
1/0
Carla
100
STPP
Deep
20
STPP
1.0
30
03
20
02
01
0
140
150
160
170
time
75
80
85
90
time
Figure Comparing the ground truth conditional intensity X \" t) learned intensity on Hawkes Dataset and ST Self-Correcting Dataset 3. Top row: Ground truth. Second row: Our AutoSTPP. Rest of the rows: Baselines. The crosses top represent past events. Larger crosses indicate more recent events.
the product of their respective integral networks. Then we can apply ID non-negative constraint to each of the derivative networks. The other partial derivatives not constrained because F.) 1 F3 can negative.
function is then F2 (1))
One limitation of such decomposition is that it can only learn the joint density of marginally indepen- dent distribution. We circumvent this issue by parameterizing the influence function as the sum of N ProdNet, Σ (S,)), (3.). (t). The formulation is no longer marginally independent since not multiplicative decomposable. Figure shows that the sum of two ProdNets is already sufficient to approximate a function that cannot be written as the sum of products of positive functions. In contrast, the constrained triple AutoInt's intensity is convex everywhere and fails fit ʃ.
Increasing the number of ProdNet improves AutoInt's flexibility at the cost of time and memory. We perform an ablation study of the number of ProdNet in Appendix F
3.7 Model Training
Given the integral network h) := (o 1, h) and the derivative network approximating the influence function fo 8Fa the log-likelihood of an event sequence H., =1( ( (n))
observed in time interval [0, T] with respect to the model
[-1 fols, - s,,t - ty, ti) Σ 1=1
[(Hr) =
Σ log
Σ
Fo - t4 h,) Fo(0,
Obtaining the above from Equation is straightforward by the Fundamental Theorem Calculus. Fo is evaluated using the Divergence theorem. We can learn the parameters 0 in both networks by maximizing the log-likelihood function. experiments, we parametrize so with two AutoInt networks L and M. Each network is a sum of N ProdNets. That is,
83 M,da
83 L1d
fe(s,t) = Σ fe, (s. =
Σ
8x
8z
x=1 (0
We name our method Automatic Spatiotemporal Point Process (AutoSTPP).
4 Experiments
We compare the performances of different neural STPPs using synthetic and real-world benchmark data. For synthetic data, our goal is to validate our AutoSTPP can accurately recover complex

---

Table 1: Test log likelihood (LL) and Hellinger distance of distribution (HD) on synthetic data (LL higher is better, HD lower is better). Comparison between AutoSTPP, NSTPP, Monte Carlo STPP, on synthetic datasets from two types of spatiotemporal point processes.
Spatiosemporal Hawkes process
Spatiotemppral Self Correcting process
L1.
HD
1,1
HD
1.1
HD
1.1.
HD
1.1
HD
1,1.
HD
NSTPP
-5.3110
0.5341
-4.8561
0.5849
-3.7366
0.1498
-2.0759
0.5426
-2.3612
0.3933
-3.0599
0.3097
DSTPP
0.0033
.3.1142
0.1813
-1.3927
Monte Carlo STPP
-4.0066
0.3198
-3.2778
0.3780
-3.7704
0.2587
-100317
0.1224
-1.3681
0.1163
-1.4439
0.2354
AutoSTPP
-3,9588
0.3018
-2.5304
0.1891
33.7700
0.1495
-10026
0.1216
-.3657
0.1119
-1.3979
0.2181
intensity functions. Additionally, we show that the errors resulting from numerical integration lead to higher variance in the learned intensity than closed-form integration. We show that our model performs better on par with the state-of-the-ar methods for real-world data.
4.1 Experimental Setup
Synthetic Datasets. We follow the experiment design Zhou et al. [2022] to validate that our method can accurately recover the true intensity functions of complex STPPs. We use six synthetic point process datasets simulated using Ogata's 's thinning algorithm [Chen] 2016). see Appendix for details. The first three datasets were based spatiotemporal Hawkes processes, while the remaiingg three were based on spatiotemporal self-correcting processes. Each dataset spans a time range of [0, 10000) and is generated using fixed set of parameters. Each dataset divided into training, validation, and testing set in an 8 1: ratio based on the time range.
Spatiotemporal Hawkes process (STH). A spatiotemporal Hawkes process, also known self-exciting process, posits that every past event exerts an additive, positive, and spatially local influence future events. This pattern is commonly observed in social media and earthquakes. The process is characterized by the intensity function [Reinhart 2018):
(²ss == (go(s)
Σ (t, (t, 1>00
9o represents the density of the background event distribution over S. 92 represents the density the event influence distribution centered at s; and over S. 91 describes each event t, t,'s influence decay over time. We implement go and O2 as Gaussian densities and g2 as exponential decay functions.
Spatiotemporal Self-Correcting process (STSC). A spatiotemporal Self-Correcting process assumes that the background intensity always increases between the event arrivals. Each event discretely reduces the intensity in the vicinity. The STSC is often used for modeling events with regular intervals, such as animal feeding times. It is characterized by:
(s, = exp (go(s)
Σ (gg((sss)) a, 3, a
(6)
(8) again represents the density of the background event distribution. g2(s, represents the density of the negative event influence centered 8,.
See the Appendix for the simulation parameters of the six synthetic datasets.
Real-world Datasets. We follow the experiment design of Chen et al. 2020 and use two of the real-world datasets, Earthquake Japan and COVID New Jersey. The first dataset includes information on the times and locations of all earthquakes in Japan between 1990 and 2020, with magnitudes of at least 2.5. This dataset includes 1050 sequences over a [0, 30) time range and is split with a ratio of 950 50 50. The second dataset is published by The New York Times and describes COVID-19 cases in New Jersey at the county level. This dataset includes 1650 sequences over a [0, 7) time range and is split with ratio of 1450 100 100.
Evaluation Metrics. For real-world datasets, we report the average test log-likelihood (L.I.) of events. For synthetic datasets, the ground truth intensities available, report the test likeliioddoood

---

(LL) and the time-av average Hellinger distance between the learned conditional spatial distribution (s(t) and the ground truth distribution. The distributions are estimated as multinomial distributions at k discretized grid points. The Hellinger distance is then
P (p:..... ..... p and Q ((i..... ..., 4s.
calculated H(P, Q)- J V E (vp: v
Baselines. We compare with two state-of-the-art neural STPP models, NSTPP [Chen et al., 2020 and Deep-STPP [Zhou et al. 2022). We also design another baseline, Monte Carlo STPP, which uses the same underlying model as AutoSTPP but applies numerical integration instead of Autolnt to calculate the loss. For fair comparison, Monte Carlo STPP models the influence fğ (8. t) as multi-layer perceptron instead of a derivative network. This numerical baseline aims to demonstrate the benefit of automatic integration.
4.2 Results and Discussion
Comparing baseline and efficient impl of AutoInt, mixed derivative, forward
Exact Likelihood. Not all baseline methods have exact likelihood available. NSTPP [Chen et al. 2020] uses numerical Neural-ODE solver. DeepSTPP [Zhou et al. [022) introduces VAE that optimizes lower bound of likeli- hood. Monte Carlo STPP estimates the triple integral of intensity by Monte Carlo integration. As such, the results presented in Table [ are es- timated for these baseline models, whereas for our model, the results reflect the true likelihood.
baseline
aurs
0.8
relative time
0.6
0.4
2
3
Advantage of Autolnt. We visualize and com- pare two sample sets of intensities from STH Dataset and STSC Dataset 3 in Figure 3 While the influence function approximator in Monte Carlo STPP theoretically approxi- mate more functions than Auto-STPP, the inten- sity it learns is "flatter" than the intensity learned by Auto-STPP. This flatness indicates a lack of information regarding future event locations.
order derivative
Figure 4: Forward mixed (d/da, dæ2 ... dr.) par- tial derivative computation average speed compari- son between our efficient implementation and Py- Torch naive AutoGrad, for a two-layer MLP.
The primary reason behind this flatness can be traced back to numerical errors inherent in the Monte Carlo method. The Monte Carlo integration performs poorly for sharply localized integrants because its samples are homogeneous over the phase space. The intensity of ST-Hawkes is a localized function; it is close to zero in most of the places but high near some event locations. As a result, Monte Carlo STPP can hardly recover the sharpness of the ground truth intensity. NSTPP also uses Monte Carlo integration and suffers the same drawback on STH Dataset 1. In contrast, AutoSTPP evaluates the integral with closed-form integration and alleviates this issue.
Synthetic Datasets Results. Table compares the test LL and the Hellinger distance between AutoSTPP and the baseline models on the six syn- thetic datasets. For the STSC datasets, we can see that AutoSTPP accurately recovers the intensity functions compared to other models. In Figure 3 AutoSTPP is the only model whose peak intensity location is always the same as the ground truth. DeepSTPP does not perform well in learning the dynamics of the STSC dataset; it struggles to align the peak and tends learn the flat intensity overly sharp. The peak intensity location of NSTPP is also biased.
Table 2: Test log likelihood (LLL) comparison for space and time on real-world benchmark data, mean and standard deviation over three runs.
LL
COVID-19 NY
Earthquake JP
NSTPP
2.5566
-4.4949 1172
DSTPP
2.3133
-3.9852
MonteSTPP
2.1070 Lo.000
-3.6085 10.056
AutoSTPP
2.6243+n.500 40.5005
-3.5948 2594482000025
For the STHP datasets, DeepSTPP has clear advantage because it uses Gaussian kemels to ap- proximate the Gaussian ground truth. Table [ AutoSTPP outperforms all other models except DeepSTPP, and its performance is comparable. Figure 3 shows that Monte Carlo STPP and NSTPP

---

can only learn an unimodal function, whereas AutoSTPP can capture multi-moda behavior in the ground truth, especially the small bumps near the mode.
Real-world Datasets Results. Table 2 compares the test LL of AutoSTPP against the baseline models on the earthquakes and COVID datasets. Our model demonstrates superior performance, outperforming all the state-of-the the--e--att--- art methods. One should note that while Monte Carlo STPP shows performance comparable ours on the Earthquake JP dataset, it falls short when applied to the COVID-19 NY dataset. We attribute this discrepancy to the large low-poy www-popuatiiondensiii areas in the COVID-19 NY data, which causes higher numerical error in integration.
4.3
Computational Efficiency
Figure 4 visualizes the benefit of using our implementation of AutoInt instead of the PyTorch naive implementation. More visualizations of the forward and backward computation times can be found in Appendix [ We can see that our implementation can be extended to compute any order of partial derivative. It is significantly faster than the naive autograd. In our AutoSTPP, we calculate the intensity using the product of three first-order derivatives. Our implementation would lead to speedup of up to 68% for computing each first-order derivative.
Conclusion
We propose Automatic Integration for neural spatiotemporal point process models (AutoSTPP) using a dual network approach. AutoSTPP can efficiently compute the exact likelihood of any sophisticated intensity.
We validate the effectiveness of our method using synthetic data and real-world datasets and demon- strate that it significantly outperforms other point process models with numerical integration when the ground truth intensity function localized.
However, like any approach, ours is not without its limitations. While AutoSTPP excels in computing the likelihood, sampling from AutoSTPP is computationally expensive as we only have the expression of the probability density. Closed-fore formm computation of expectations is also not possible; Knowing the form of f x(t). calculating f +t()) is still intractable.
Our work presents a new paradigm for learning continuous- tamin dynamics. Currently, our neural process model takes the form of Hawkes processes (self-exciting) but cannot handle the discrete decreases of intensity after events due to the difficulty of integration. Future work includes relaxing the form of the intensity network with advanced integration techniques. Another interesting direction is to increase the approximation ability of the product network.
Acknowledgement
This work was supported in part by U. S. Army Research Office under Army-ECASH award I 07-R-0003-03, the U.S. Department Of Energy, Office of Science, IARPA HAYSTAC Program, NSF Grants #2205093, #2146343, and #2134274.
References
Ryan Prescott Adams, Iain Murray, and David JC MacKay. Tractable nonparametric bayesian inference in poisson processes with gaussian process intensities. In Proceedings of the 26th annual international conference on machine learning, pages 9-16, 2009.
Anders Brix and Peter Diggle. Spatiotemporal prediction for log-gaussian cox processes. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63(4))823 841, 2001.
Anders Brix and Jesper Moller. Space-time multi type log gaussian cox processes with a view to modelling weeds. Scandinavian Journal of Statistics, 28(3):471-44 1-888, 2001.
Ricky TQ Chen, Brandon Amos, and Maximilian Nickel. Neural spatio-temporal point processes. arXiv preprint arX1v:2011.0458,, 2020.
10

---

Yuanda Chen. Thinning algorithms for simulating point processes. Florida State University, Talla- hassee, FL, 2016.
Daryl Daley and David Vere-Jones. An introduction to the theory of point processes: volume II: general theory and structure. Springer Science & Business Media, 2007.
Hennie Daniels and Marina Velikova. Monotone and partially monotone neural networks. IEEE Transactions Neural Networks, 2166)906-917, 2010
Philip] Davis and Philip Rabinowitz. Methods of numerical integration. Courier Corporation, 2007.
Nan Du, Hanjun Dai, Rakshit Trivedi, Utkarsh Upadhyay, Manuel Gomez-Rodriguez, and Le Song. Recurrent marked temporal point processes: Embedding event history to vector. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pages 1564,1564, 2016.
William Dunham. The calculus gallery. Princeton University Press, 2018.
Ruocheng Guo, Jundong Li, and Huan Liu. Initiator: Noise-contrastive estimation for marked temporal point process. In IJCAI, pages 2191-2197, 2018.
Hengguan Huang, Hao Wang, and Brian Mak. Recurrent poisson process unit for speech recognition. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 6538-6545, 2019.
Haibin Li, Yangtian Li, and Shangjie Li. Dual neural network method for solving multiple definite integrals. Neural computation, 31(1):208 2220008--22211111 2019.
Shuang Li, Shuai Xiao, Shixiang Zhu, Nan Du, Yao Xie, and Song. Learning temporal point
processes via reinforcement learning. arXiv preprint arXiv: 1et:11011.000,, 2018.
David B Lindell, Julien NP Martel, and Gordon Wetzstein. Autoint: Automatic integration for fast neural volume rendering. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 14556-14565, 2021.
Keqin Liu, Automatic integration. arXiv e-prints, pages arXiv-2006, 2020.
Jerrold F. Marsden and Anthony Tromba. Vector Calculus. Macmillan, August 2003. ISBN 978-0-7167-44992-9 Google-Books-ID: LiRLJf2m_dwe
Hongyuan Mei and Jason Eisner. The neural hawkes process: neurally self-modulating multivariate point process. arXiv preprint arXiv:1612.09322 2016.
Yosihiko Ogata et al. The asymptotic behaviour of maximum likelihood estimators for stationary point processes. Annals of the Institute of Statistical Mathematics, 30(1):243-261, 1978.
Takahiro Omi, Naonori Ueda, and Kazuyuki Aihara. Fully neural network based model for general temporal point processes. arXiv preprint arXiv:1905.0963 2019.
Alex Reinhart. A review of self-exciting spatio-temporal point processes and their applications. Statistical Science, 33(3):299-333 318, 2018.
Robert H Risch. The problem of integration in finite terms. Transactions of the American Mathemati- cal Society, 139:167-189, 1969.
Robert H Risch. The solution of the problem of integration in finite terms. Bulletin of the American Mathematical Society, 76(3):605-60) 1970.
Jin Shang and Mingxuan Sun. Geometric hawkes processes with graph convolutional recurrent neural networks. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 33, pages 4878 4885, 2019.
Oleksandr Shchur, Marin Biloš, and Stephan Günneemann. Intensity-free learning temporal point processes, arXiv preprint
2019.

---

Utkarsh Upadhyay. Abir De, and Manuel Gomez-Rodriguez. Deep reinforcement learning of marked temporal point processes. arXiv preprint arXiv:1805.09366, 2018.
Shuai Xiao, Mehrdad Farajtabar, Xiaojing Ye, Junchi Yan, Le Song, and Hongyuan Zha. Wasserstein learning of deep generative point process models. arXiv preprini arXiv: 1705.0805, 2017.
Junchi Yan, Xin Liu, Liangliang Shi, Changsheng Li, and Hongyuan Zha. Improving maximum likelihood estimation of temporal point process via discriminative and adversarial learning. In IJCAI, pages 2948-222 2954,, 2018.
Qiang Zhang, Aldo Lipani, Omer Kimnap, and Emine Yilmaz. Self-attentive hawkes process. In International Conference on Machine Learning, pages 11183-11193 PMLR, 2020.
Zihao Zhou and Rose Yu. Automatic integration for fast and interpretable neural point process. In Learning for Dynamics and Control Conference. PMLR, 2023.
Zihao Zhou, Xingyi Yang, Ryan Rossi, Handong Zhao, and Rose Yu. Neural point process for learning spatiotemporal event dynamics. In Learning for Dynamics and Control Conference, pages 777-789. PMLR, 2022.
12

---

A More Implementation Benchmarks
Comparing baseline and efficient impl AutoInt, mixed derivative, forward
Comparing baseline and efficient impl Autoint, derivative, forward
baseline
1.2
baseinn
reletive time
relative time
0.8
0.2
order of derivative
order of derivative
Figure 5: Forward mixed (d/dz;d (t//dddz...... ddz...drz:) partial derivative computation speed, 3 layers MLP
Figure Forward mixed (d/dr, drr...drs))
partial derivative computation speed, 4 layers MLP
As the number of MLP layers increases, the lower-order derivative computation becomes faster (relative to PyTorch naive implementation). whereas the higher-order derivative computation becomes slower. This performance pattern is because our implementation uses Python for loop. Our approach faster than the baseline in the majority of cases.
Comparing baseline and efficient impl of Autoint, mixed derivative, backward
Comparing baseline and efficient impl of AutoInt, univariate derivative, backward
baseline
baseline
0.8
relative time
0.8
relative time
0.4
order of derivative
order derivative
Figure 7:
Forward Backward mixed
Figure 8: Forward Backward univariate (d/dx ...
(d/da, dxz ... dr&) partial derivative compu- tation speed, layers MLP
... dr 1) partial derivative compu- tation speed, 3 layers MLP
As for both forward and backward computation times, our implementation still consistently surpasses the baseline terms speed. Notably, our method is even more efficient when calculating univariate partial derivatives than mixed partial derivatives. This advantage is primarily due to the reduced number of iterations required by the Python for-loop in the case of univariate derivatives.
B STSC and ST-Hawkes Introduction and Simulation Parameters
We the parameters as Zhou et al. 2022
The STSCP's and the STHP's kernels go (s) and
92 (s, sy) are prespecified to be Gaussian:
gu(s) 1
1 (s [0, (0,0))=00 (s 10. (0.0))²)
Egol 1 exp
2n
2
(2 8,) =
1
g2
1 exp
1 - (s s, 5 (s -s,)²)
2m
2
The STSCP is defined on S 10 x [0. 1), while the STHP is defined on TR.2. The STSCP's
kernel functions normalized according to their cumulative probability S. Table shows the
13

---

simulation parameters. We discretized the STSCP's spatial domain as a 101 X 101 grid during the simulation.
Table 3: Parameter settings for the synthetic dataset
α
B
Ego
Eg²
ST-Hawkes
DS1
.5
1
.2
[.2 0; 0 2]
[0.5 0; 0.5]
DS2
.5
.6
.15
[5 0; 05]
O.1 0.1]
DS3
.3
2
1
[1 0000]
[.1 0; ..]]
ST-Self Correcting
DSI
.2
.2
1
[1
[0.85 0; 0 0.85]
DS2
.3
[,40;
[.3 0; O .3]
DS3
.4
.2
1
[.25 0; 0 .2]]
[.2 0; 0 .2]
Each dataset is single, long sequence that spans over 10,000 time units. We divide each dataset into 50 sequences, each spanning 200 time units. We use 40 sequences for training. 5 for validation, and for testing. Here's summary of the total number events found in each dataset:
Table 4: Number of events in each synthetic dataset
number of events
ST-Hawkes
DSI
3983 9017
DS2
DS3
11693
ST-Self Correcting
DSI
10002
DS2
6668
DS3
5004
The prediction task applies sliding windows each of the datasets. We try use historical events predict the likelihood of the next event in the same sequence.
C Model Setup Details
We detail the specific hyperparameter settings in Table Except for the learning rate, the same set of parameters was applied across all datasets. Despite varying datasets, This consistency in performance demonstrates our model's robustness to hyperparameters.
Name
Value
Description
Optimizer
Adam
Learning rate
Depends on dataset, [0.0002, 0.004]
Momentum
0.9
Adam momentum
Epoch
50 / 100
50 for synthetic dataset and 100 for real-world dataset
Batch size
128
Activation
tanh
Activation function in L and M (intensity parameter networks) in and M
N
1/
Number of product nets for synthetic dataset and 10 for real-world dataset and use bias their linear layers
bias
true
Table 5: Hyperparameter settings for training AutoSTPP on all datasets.
D Forward-pass Algorithm for Automatic Integration
Function: dnforward(f, n, x, dims), partition(n, k) finds all k-subset partitions of
14

---

Data: n, dimension of fz), tensor of shape (batch, dim), dims, list of dimensions to derive, layers, composite functions in f Result: dains f f/dx
Initialize dictionary dnf, mapping from dims to d f //dadine empty list pd if d'mml then
Precompute f())
pd (E(I.))
else
for subdims E combination dims, lení dims)-I) do Precompute dnf anforward(f, 2, subdims) end
pd <- 0
end
for Lager Layers do
layer is linear then
pd append last pd xw W is the linear weight
else if Layer is activation then
if d'ims then
termsum + last pd layer' (f)
else
termsum +
for order 0,..., d'mms do if order 0 then
term t last pd
else
term f
for part partition dims, order /) do
temp 1
for subdims part do temp + temp dnf dnforwarddf... n, æ, subdims) (precomputed)
term +- term + temp
end
termsum t termsum + term
end end
termsum + termsum X layer (f)
end
end pd append termsum
end
return last pd
E Universal Approximation Theorem for Derivative Network
Consider AutoInt integral network with the form
(xx) C (a (A b)),
A A Rk b b C
where a denotes R + continuous non-polynomial function applied elementwise to each input dimension.
The derivative network thus takes the form
g (x) = C (a' oA + Acot),
where Anot C & is column of A that corresponds the deriving dimension.
Recall the universal approximation theorem (Daniels and Velikova, 2010). which says for every compact K C R and C(K, 0, there exist A, b. C such that
sup [If(x) - g(x)) < zek
15

---

Proposition E.1. (Universal Approximation Theorem for Derivative Network) for every compact K C In and f E CR, 0, there exists R E R C R E R. such that
g( = C (ar (Ax b)) 3x
sup f (z) <E
=EK
Proof. Given the mapping J. by UAT, there exists A, b, C that approximate ((x). Construct Ć IR*
and B E IR, such that
Acal.j # 0
C2= 10, 10.
Acol.j 0 and == Σ Cga(by)
Then,
sup =EK
[f(x)) C. (o (
k f(x) C x bb))))
sup zEK
1=1
sup
f(x) Σ Cy(o(A; x -2 Ab)
E
Cjo(by)
=EK
1=1
s1Acolu=0
[1(x)) - Č (o' (A x + b) Acol) -B1,
sup
=EK
Note that Č (a' o (A-x Acot) - d (c (00 (oo b]) Bz). which is the derivative
dx
of two-layer feedforward integral network.
F Relationship between Number of ProdNets and Model Expressivity
"192""
[UL--IL procret
6.01
0.00
0.06
0.00
Figure 9: Training MSE for fitting positive derivative network to sin(z) cos(y) sin(2) + +1
We applied ten summations of positive ProdNets to fit the noo-mutiplicaieeeecompaabe function sin (1) cos(y) sin(2)
Each of these ProdNets consists of three MLP components, each with two hidden layers with 128 dimensions. All models underwent training with consistent learning rate at 0.005.
Our results, shown in Figure indicate that increasing the number of ProdNets generally improves the model's performance in fitting the non-dccomposable function. The model with the best MSE
16

---

uses 10 ProdNets, while the model with ProdNets had the second-wors performance. This result intuitively sensible, more linear terms are typically required express an arbitrary function precisely. However, we observed that employing more ProdNets does not always lead to better performance, as demonstrated by the model with 8 ProdNets.
17


Reviews:
[{"review_id": "QFSV3pFJOv", "reviewer": "Reviewer_HnTF", "summary": "In this paper is proposed an automatic integration scheme for efficient learning for spatiotemporal process models, where the corresponding function is approximated with a deep neural network. A priori, the method can be used to compute the likelihood of any intensity, even for sophisticated ones. The method is validated in both synthetic and real datasets, where the proposed method obtains better results than those provided by state-of-the-art approaches. ", "questions": "As the proposed method is available for learning continuous-time dynamics, I think the authors could consider more challenging scenarios where the applicability is especially relevant in science and technology, such as fluids or continuous mechanics. I would like to see at least a discussion about this type of scenario. \n \nAblation studies by considering the parameters in table 3. The authors should indicate how the best parameters were selected, and if any variation is given in the parameter tuning for both synthetic and real experiments.  ", "limitations": "The authors do not include a proper discussion about the limitations of the method as it is a bit superficial. Moreover, some failure cases could help the reader. To this end, including more complex intensity functions where the method obtains a worse performance could greatly improve the paper. Thanks to that, we can see the real limitations of the method. ", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "The paper is well written and clear enough. In general, the authors provide most of the details, explanations and assumptions to understand their approach and the corresponding contribution. To be honest, I do not have many issues with the current submission. \n \nThe method shows a clear advance in terms of accuracy in comparison with imposing the constraint through activation function with a nonnegative triple derivative. AutoSTPP performance is very competitive with respect to previous neural STPPs. \n \nBoth quantitative and qualitative experiments are provided in the paper. Both of them show the superiority of the method, validating all the claims in the paper. In particular, I like the experiments on synthetic data, as complex intensity functions are considered and the proposed method can capture them correctly. In any case, the real-world scenarios are also needed, and the set of experiments in this line is well executed. Moreover, the method clearly outperforms the rest of competing methods in these experiments. ", "weaknesses": "Some parameters in Eq.(1) are not properly defined. See, u or \\tau, for example. I know many parameters are included in the paper so please, be careful with all of them. \n\nThe authors do not include a proper discussion about the limitations of the method as it is a bit superficial. Moreover, some failure cases could help the reader. To this end, including more complex intensity functions where the method obtains a worse performance could greatly improve the paper. Thanks to that, we can see the real limitations of the method. \n\nI am not sure why the authors only use two real-world datasets of Chen et al. 2020, as more cases were considered there. At least, that could be briefly commented. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for reviewing our paper. We greatly appreciate the reviewer's thorough feedback and the positive evaluation of our paper.\n\n**Q: Some parameters in Eq.(1) are not properly defined.** \n\n**A:** Sorry for not properly defining tau in Equation (1). Thank you for pointing this out! We will make sure to address this in the updated manuscript.\n\n**Q: The authors do not include a proper discussion about the limitations of the method as it is a bit superficial.** \n\n**A:** We agree that a more comprehensive discussion on the limitations and potential failure cases of our method would be valuable. One limitation is the \"same influence\" assumption, which may introduce an inductive bias when the influence of different events significantly differs. However, our empirical findings suggest that removing this assumption can lead to worse performance, presumably due to a lack of estimator consistency. Exploring this limitation further would be an interesting avenue for future research.\n\n**Q: I think the authors could consider more challenging scenarios where the applicability is especially relevant in science and technology, such as fluids or continuous mechanics.** \n\n**A:** Thank you for your suggestion. Fluids or continuous mechanics are not discrete event data, hence are not directly applicable to our paper. The earthquake datasets we used in the paper are challenging datasets for event prediction. \n\nWe stress that the selection of datasets in our experiments do some favors for the baselines (in fact, they are from the baseline papers) models and are relatively challenging for our model. Our synthetic datasets, such as the spatiotemporal Hawkes dataset, feature Gaussian influence functions that are assumed by DSTPP. Incorporating a dataset with more complex influence functions, like a sine wave, while still maintaining the \"same influence\" assumption, could provide a clearer performance gap compared to baseline methods.\n\n**Q: Ablation studies by considering the parameters in Table 3.** \n\n**A:** We follow the simulation parameter settings in the DSTPP paper, and they are only for generating the datasets. One advantage of our proposed method is that it has very few hyperparameters (number of ProdNet and hyperparameters for MLP).\n"}, {"review_id": "dEmftvpvO2", "reviewer": "Reviewer_Hisf", "summary": "The paper introduces a new approach to address the challenge of integrating spatiotemporal point processes in a flexible and efficient manner. The existing methods either lack flexibility or introduce numerical errors. The authors propose AutoSTPP as a novel paradigm that overcomes this limitation and showcases its advantages in recovering complex intensity functions from irregular spatiotemporal events. Its effectiveness is validated through experiments on synthetic data and benchmark real-world datasets.", "questions": "What’s the key contribution of this paper compared to [Zhou et al 2023]?  The authors should carefully differentiate their approach from it.", "limitations": "None", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Clear presentation. The writing is clear and concise.\n2. It provides the proof of consistency for AutoSTPP, which enhances the credibility of the proposed method.\n3. Comprehensive experiments. The authors evaluate the proposed approach on both synthetic datasets and real-world datasets.", "weaknesses": "1. An essential related work is missing: Automatic Integration for Fast and Interpretable Neural Point Processes [1]. \n2. Despite the solution for STPP is very interesting, the proposed technique lacks novelty. Zhou et al.[1] has already proposed a very similar framework for TPPs. The major different part lies in the data dimension. I don’t think just extending 1D to 3D is qualified for an accepted NIPS paper.  \n\n[1] Zhou, Zihao, and Rose Yu. \"Automatic Integration for Fast and Interpretable Neural Point Processes.\" Learning for Dynamics and Control Conference. PMLR, 2023.\n\n## After Rebuttal\n\nThe authors' response has addressed most of my concerns. Therefore, I decide to raise the score to 6 (Weak accept).\n\nHowever, I highly suggest authors carefully discuss links between their work and Zhou et al. [1].\n\nMoreover, recently developed diffusion (or score-based) models seem like a promising workaround for handling likelihood computation in the neural point process. Authors should discuss the following related works that have been published recently.\n- Li, Zichong, et al. \"SMURF-THP: Score Matching-based UnceRtainty quantiFication for Transformer Hawkes Process.\" in ICML 2023.\n- Yuan, Yuan, et al. \"Spatio-temporal Diffusion Point Processes.\" in KDD 2023.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for providing feedback on our paper. We appreciate your positive remarks regarding the soundness and presentation of the paper. Below, we address the concerns you raised regarding the contribution.\n\n**Q: An essential related work is missing: Automatic Integration for Fast and Interpretable Neural Point Processes [1].**\n\n**A:** This paper was scheduled to appear in the 5th L4DC conference (June 14th~16th), but it has not been published or made publicly accessible before the NeurIPS submission deadline (May 11). Hence, we were unable to include it in the related work in our paper prior to submission. According to the NeurIPS [two month policy](https://neurips.cc/Conferences/2023/PaperInformation/NeurIPS-FAQ), Papers appearing less than two months before the submission deadline are generally considered concurrent to NeurIPS submissions. authors are **not** expected to compare to work that appeared only a month or two before the deadline, let alone **after the deadline**. \n\nNevertheless, we recognize the significance of their work and will ensure to include it in the final version of our manuscript. We will also discuss the additional challenges our proposed method addresses compared to their problem settings.\n\n**Q: The major different part lies in the data dimension. I don’t think just extending 1D to 3D is qualified for an accepted NIPS paper.**\n\n**A:** We respectfully disagree. Despite the naming similarity of AutoInt between our approach and the work by [Zhou et al.], the extension from 1D to nD (n>1) is a highly **non-trivial contribution**. Zhou et al. can only learn temporal point processes without markers, whereas our paper focuses on extensible spatiotemporal point processes.  The difference in data dimension makes the problem significantly harder, hence introducing several technical differences. One of the major challenges is the inherent difficulty of constraining the sign of higher-order derivatives without compromising the performance. As detailed in Section 3, if we simply extend the 1D AutoInt technique in [Zhou et al], we would overconstrain the model and fail to recover the ground truth intensity. Another technical challenge is the computational complexity. The computational cost of AutoInt scales exponentially with the number of dimensions. Our proposed method tackles these technical challenges through many novel techniques including ProdNet and Divergence theorem. We circumvent the trivariate computational graph and use simplified univariate graphs to make AutoInt scales with the number of dimensions."}, {"review_id": "OzSO1NQHuQ", "reviewer": "Reviewer_wohe", "summary": "The paper proposes a novel framework to extend the automatic integration of conditional intensity function to spatio-temporal point processes, aiming to overcome the computational challenges of triple integration. Specifically, the proposed model adopts an integral network and converts the triple integral into line integrals using the Divergence theorem and Green’s theorem. The authors use a summation of ProdNet to handle the non-negativity of the conditional intensity function (derivative network) instead of directly imposing a non-negative constraint on the model parameters, which enhances the model expressiveness. The consistency, validity, and computational efficiency of the AutoInt STPP are comprehensively demonstrated.", "questions": "1. Forgive me if I understand wrong, but is there an integral of f over S missed in the integral at line 151?\n2. I would like to know how N (in ProdNet) is chosen in experiments, and if there is a significant performance difference with different N.\n3. It would be great if the authors could explain the advantages of AutoSTPP over DSTPP.\n4. The computational efficiency of the derivative is great, however, the readers or practitioners would like to know the comparison of the computational efficiency of the model with other baselines.", "limitations": "1. Modeling performance and computational efficiency can be further explained or demonstrated.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The proposed framework based on auto-integration is valuable for spatio-temporal point process modeling. It significantly improves the model flexibility and accuracy compared to the close-form or monte-carlo based integration model.\n2. The natural constraint of the spatio-temporal point processes (e.g., non-negativity of the conditional intensity function) is well and not overly satisfied. The triple integral has also been handled in a novel way.\n3. The comprehensive experimental results on both synthetic and real-world data, including the ablation study, well support the claimed advantage of the proposed model.\n4. The consistency of the AutoInt PP estimator is proved under mild constraints.", "weaknesses": "1. Although it does not impact the overall goodness of the paper, I have to admit that the presentation of Sec. 3.5 is not straightforward. The purpose and advantage of using line integral to handle the triple integral are not clear to me. It would be better if the authors could further explain it with several high-level sentences.\n2. More experiments on the spatio-temporal Hawkes process with a non-separable influence function can be included. For example, the model performance on synthetic data generated by the ETAS [1], which is a widely-used STPP model in real-world practice.\n3. I would say the performance advantage of AutoSTPP is not significant or highlighted. For example, it seems like the DSTPP has a similar performance to AutoSTPP on synthetic data, and so does MonteSTPP if the standard deviation is considered on real data.\n\n\n---\n[1] Space-Time Point-Process Models for Earthquake Occurrences. Yosihiko Ogata. Annals of the Institute of Statistical Mathematics, 1998, vol. 50, issue 2, 379-402", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review of our paper. We appreciate your positive feedback on the proposed framework and the comprehensive experimental results supporting our claims. We acknowledge the weaknesses and limitations pointed out and will address them accordingly.\n\n**Q: I have to admit that the presentation of Sec. 3.5 is not straightforward. The purpose and advantage of using line integral to handle the triple integral are not clear to me.** \n\n**A:** Regarding the presentation of Section 3.5, we apologize if it was not clear enough. The triple integral is necessary for computing the likelihood during the training stage. The line integral is not particular; we convert the triple integral over a cuboid to double integrals over 6 squares, then line integrals over 24 lines. It is a consequence of sequential application of Divergence theorem, Green’s theorem and the Fundamental theorem of calculus. AutoInt guarantees that the triple, double and line integral all have closed forms. We will make sure to provide a clearer high-level explanation in the updated draft.\n\n**Q: More experiments on the spatio-temporal Hawkes process with a non-separable influence function can be included. For example, the model performance on synthetic data generated by the ETAS [1].**\n\n**A:** We appreciate the suggestion. However, the intensity of the ETAS model is still additively separable, based on the equation mentioned in the 4th page of [1]. The current synthetic spatio-temporal Hawkes process in our paper is a simplified version of ETAS without the term t^{1+\\beta}. Indeed, adding more point processes with intensity that is not additively separable would be valuable. That is why we included the spatiotemporal extension of the self-correcting process in the synthetic experiments.\n\n**Q: I would say the performance advantage of AutoSTPP is not significant or highlighted.** \n\n**A:** We agree that the advantages of AutoSTPP over DSTPP should be discussed more. One major advantage of AutoSTPP is its ability to learn any intensity without relying on kernel density estimation. This not only reduces the number of hyperparameters but also improves flexibility. When examining the estimated intensity of the self-correcting process, one may notice that DSTPP incorrectly approximates the underlying intensity with a mixture of Gaussian kernels, while AutoSTPP learns a much smoother intensity. We will include a more detailed discussion in the final draft.\n\n**Q: Forgive me if I understand wrong, but is there an integral of f over S missed in the integral at line 151?**\n\n**A:** It's a typo and thank you for pointing it out! \n\n**Q: I would like to know how N (in ProdNet) is chosen in experiments, and if there is a significant performance difference with different N.**\n\n**A:** Please refer to Appendix F: while N=2 could reduce the performance by a margin, further increasing the number of ProdNet has limited effect on performance, and does not always boost the performance. It shows that N=5 can already cover a sufficient number of functions. We selected N by this gridsearch and fixed it to be 10 during experiments.\n\n**Q: The computational efficiency of the derivative is great, however, the readers or practitioners would like to know the comparison of the computational efficiency of the model with other baselines.**\n\n**A:** Our current records show that both DSTPP (\\~1 hour) and AutoSTPP (1\\~3 hour) are much faster than the normalizing-flow-based NSTPP (> 12 hours). AutoSTPP is often slower than DSTPP due to the larger computation graph, which can be influenced by the choice of parameter N. We will include a table containing the multirun averages in the appendix to provide a comprehensive comparison.\n"}, {"review_id": "TIp3zE5JTR", "reviewer": "Reviewer_WN8G", "summary": "In this paper, the authors propose a new method for spatiotemporal neural point processes (STPP). A key innovation of their method is that they introduced AutoInt for efficient integration of the intensity function of STPP. Using the divergence theorem and Green’s theorem, they convert the original triple integral to line integrals which AutoInt can efficiently calculate. They also impose 3D non-negativity constraint on the integral network. Experimental results show efficacy of the proposed method\n", "questions": "1. It’s not clear why we cannot use AutoInt to do triple automatic integration on the intensity function directly. \n2. Compared with traditional methods on integral calculation, e.g. Monte Carlo, what are the advantages of the proposed method from the viewpoint of computational complexity?\n3. To train AutoInt (or corresponding duo neural networks), are there any minimum requirement of the number of events (training samples)?", "limitations": "Yes", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Applying AutoInt to solve the STPP problem looks pretty novel.\n2. Converting the original triple integral to line integrals using the divergence theorem and Green’s theorem and imposing non-negative constraint on the integral network also have strong methodology contribution.\n", "weaknesses": "1. The paper misses an important reference which also uses AutoInt for point processes: Automatic Integration for Fast and Interpretable Neural Point Processes. So the authors are encouraged to clarify contributions of this paper over the published paper above.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your review of our paper. We appreciate your positive feedback on the novelty and methodology contribution of our proposed method.\n\n**Q: So the authors are encouraged to clarify contributions of this paper over the published paper above.**\n\n**A:** This paper was scheduled to appear in the 5th L4DC conference (June 14th~16th), but it has not been published or made publicly accessible before the NeurIPS submission deadline (May 11). Hence, we were unable to include it in our paper prior to submission. According to the NIPS [two month policy](https://neurips.cc/Conferences/2023/PaperInformation/NeurIPS-FAQ), authors are **not** expected to compare to work that appeared only a month or two before the deadline.\n\nNevertheless, we recognize the significance of their work and will ensure to include it in the final version of our manuscript. We will also discuss the additional challenges our proposed method addresses compared to their problem settings.\n\n**Q: It’s not clear why we cannot use AutoInt to do triple automatic integration on the intensity function directly.**\n\n**A:** Figure 2 illustrates why direct extension of AutoInt for triple automatic integration is not feasible due to convexity challenges. If we directly constrain the mixed triple derivative (df / dxdydt) of a neural network to be positive by setting the whole network monotone, all triple derivatives would become positive. This would lead to an overly constrained model. For instance, df / dtdtdt is also positive, leading to an intensity that skyrockets on the time axis. As a result, the intensity function is approximated using a piecewise convex function, as shown in the third plot of Figure 2. To overcome this issue, we propose the ProdNet technique.\n\n**Q: Compared with traditional methods on integral calculation, e.g. Monte Carlo, what are the advantages of the proposed method from the viewpoint of computational complexity?**\n\n**A:** The main advantages include (1) sample complexity and (2) memory efficiency. The computational efficiency comparison with other integration methods heavily depends on the number of Monte Carlo samples. With a small number of samples (10 per dimension), Monte Carlo can be faster than our method but suffers from large numerical errors, as depicted in Figure 3. With a large number of samples, Monte Carlo has high memory usage, and backpropagation needs to be done sequentially for each batch of samples, making it inevitably slower than AutoInt.\n\n**Q: To train AutoInt (or corresponding duo neural networks), are there any minimum requirements for the number of events (training samples)?**\n\n**A:** There is no strict minimum requirement on the number of events (training samples) for training AutoInt networks. The training process can be performed with a small number (e.g. tens of) of event sequences, and one can always add a loss of regularity to avoid the overfitting of intensity function. \n"}]

Author rebuttal:
Once again, we sincerely thank you for your valuable feedback and constructive criticism. We will incorporate the suggested improvements and address the raised concerns in the final version of the manuscript. We will emphasize the differences between our approach and the work by Zhou et al. published after the submission. 

Meta-review: The manuscript proposes an automatic integration scheme for efficient learning for spatiotemporal process models. While the initial reviewers were mixed, the authors quite satisfactorily answered the questions raised by the reviewers. This is also reflected in the final scores of this manuscript. I also stand with the authors on the NeurIPS *two-month policy*. However, it would be greatly appreciated to include a discussion on the differences and links to the work of Zhou et al. in the final version. Considering all the initial reviews, the rebuttal, and the reviewer discussion, I do recommend "Accept".
Decision: Accept (poster)