Additive Poisson Process: Learning Intensity of Higher-Order Interaction in Stochastic Processes
Authors: Simon Luo
Feng Zhou
Lamiae Azizi
Mahito Sugiyama
Venue:  (2021)
Keywords: Poisson Process
Log-Linear Model
Energy-Based Model
Generalized Additive Models
Information Geometry

Abstract: We present the Additive Poisson Process (APP), a novel framework that can model the higher-order interaction effects of the intensity functions in stochastic processes using lower dimensional projections. Our model combines the techniques in information geometry to model higher-order interactions on a statistical manifold and in generalized additive models to use lower-dimensional projections to overcome the effects from the curse of dimensionality. Our approach solves a convex optimization problem by minimizing the KL divergence from a sample distribution in lower dimensional projections to the distribution modeled by an intensity function in the stochastic process. Our empirical results show that our model is able to use samples observed in the lower dimensional space to estimate the higher-order intensity function with extremely sparse observations.

Under review as a conference paper at ICLR 2021
ADDITIVE POISSON PRocESS: LEARNING INTENSITY OF HIGHER-ORDER INTERACTION IN POINT PROCESSES
Anonymous authors Paper under double-bling review
ABSTRACT
We present the Additive Poisson Process (APP), novel framework that can model the higher-order interaction effects of the intensity functions in point processes using lower dimensional projections. Our model combines the techniques in in- formation geometry to model higher-order interactions on a statistical manifold and in generalized additive models to use lower-dimensional projections to over- come the effects from the curse of dimensionality. Our approach solves a convex optimization problem by minimizing the KL divergence from a sample distribu- tion in lower dimensional projections to the distribution modeled by an intensity function in the point process. Our empirical results show that our model is able to use samples observed in the lower dimensional space to estimate the higher-order intensity function with extremely sparse observations.
INTRODUCTION
Consider two point processes which are correlated with arrival times for an event. For given time interval, what is the probability of observing an event from both processes? Can we learn the joint intensity function by just using the observations from each individual processes? Our proposed model, the Additive Poisson Process (APP), provides novel solution to this problem.
The Poisson process a counting process used in wide range of disciplines such as time-space se- quence data including transportation (Zhou et al.) 2018), finance (Tlalan 2016). ecology (Thompson) 1955), and violent crime (Tadyy 2010) to model the arrival times for a single system by learning an intensity function. For given time interval of the intensity function, it represents the probability of a point being excited at a given time. Despite the recent advances of modeling of the Poisson processes and its wide applicability, majority of the point processes model do not consider the cor- relation between two or more point processes. Our proposed approach learns the joint intensity function of the point process which is defined to be the simultaneous occurrence of events.
For example in a spatial-temporal problem, we want to learn the intensity function for taxi to pick up customers at given time and location. For this problem, each point is multi-dimensional, that is (x, y, t) where a pair of æ and represents two spatial dimensions and represents the time dimension. For any given location time, we can only expect very few pick-up events occurring. therefore making it difficult for any model to learn the low valued intensity function.
Previous approaches such as Kernel density estimation (KDE) (Rosenblati 1956) are able to learn the joint intensity function. However, KDE suffers from the curse of dimensionality, which means that KDE requires large size sample or a high intensity function to build an accurate model. In ad- dition, the complexity of the model expands exponentially with respect to the number of dimensions, which makes it infeasible to compute. Bayesian approaches such as using a mixture of beta distribu- tions with Dirichlet prior (Kottas) 2006) and Reproducing Kernel Hilbert Space (RKHS) (Flaxman et al. 2017) have been proposed to quantify the uncertainty with a prior for the intensity function. However, these approaches are often non-convex, making difficult to obtain the global optimal solution. In addition, if observations are sparse, it is hard for these approaches to learn a reasonable intensity function.

---

Under review as a conference paper at ICLR 2021
All previous models are unable to efficiently and accurately learn the intensity of the interaction between point processes. This is because the intensity of the joint process is often low, leading to sparse samples or, in an extreme case, no direct observations of the simultancous event at all, making it difficult to learn the intensity function from the joint samples. In this paper, we propose a novel framework to learn the higher-order interaction effects of intensity functions in point processes. Our model combines the techniques introduced by Luo & Sugiyama (2019) to model higher-order interactions between point processes and by Friedman & Stuetzle (1981) generalized additive models to learn the intensity function using samples in a lower dimensional space. Our proposed approach is to decompose For example, in the x-dimension we have points (axa) N in the y-dimension, we have points (y+) and in the time dimension we have (tt) The data in these lower dimensional space can be used to improve the estimate of the joint intensity function. This different from the traditional approach where we only use the simultaneous events to learn the joint intensity function.
point process into nower-dimensionn representations.
We first show the connection between generalized additive models and Poisson processes. We then provide the connection between generalized additive models and the log-lineat model (Agresti) 2012), which has a well-established theoretical background information geometry (Amari [Amari|[2110 2016). We draw parallels between the formulation of the generalized additive models and the binary log- linear model on a partially ordered set (poset) (Sugiyama et al. 2017). The learning process in our model is formulated as a convex optimization problem to arrive at unique optimal solution using natural gradient, which minimizes the kullback-Leibie (KL) divergence from the sample distribu- tion in lower dimensional space to the distribution modeled by the learned intensity function. This connection provides remarkable properties to our model: the ability to learn higher-order intensity functions using lower dimensional projections, thanks to the Kelmogorov-Arrole representation the- orem. This property makes it advantageous to use our proposed approach for the cases where there are, no observations. missing samples, or low event rate. Our model flexible because it can capture interaction between processes as a partial order structure in the log-linear model and the parameters of the model are fully customizable to meet the requirements for the application. Our empirical results show that our model effectively uses samples projected onto a lower dimensional space to estimate the higher-order intensity function. Our model is also robust to various sample sizes.
FORMULATION
In this section we first introduce the technical background in the Poisson process and its extension to multi linli-dimensional Poisson process. We then introduce the Generalized Additive Model (GAM) and its connection to the Poisson process. This is followed by presenting our novel framework, called Additive Poisson Process (APP), which is our main technical contribution and has tight link to the Poisson process modelled by GAMs. We show that learning of APP can be achieved via convex optimization using natural gradient.
The Poisson process is characterized by an intensity function A: - IR where we assume multiple D processes. An inhomogeneous Poisson process is general type of processes, where the arrival intensity changes with time. The process with time-changing intensity A(t) is defined as a counting process N(t). which has an independent increment property. For all time 2 0 and changes in time
5 2 O, the probability p for the observations is given as p (N((++)) N(t) = 0) = 1 - SX(L) (M()++ + o(0). p(N(1 + 8) - N(C) = 1) 5X(1) + o()), and p(N(1 5) - N(I) 2 2) o()). where o(-) denotes little-o notation (Daley & Vone-Jones] 2007). Given a realization of timestamps t1, .,,,.., ty with t, [0, T] from an inhomogencous imulti-dimensional) Poisson process with the intensity A. Each t; is the time of occurrence for the i-th event across D processes and T' is the observation duration. The likelihood for the Poisson process (Daley & Vere-Jones) 2007) is given by
p (( A(t)) exp (-fA(O)as di Itt,,,
where t = (((1),...
(D) E RDD We define the functional prior on X(t) as X(t) g (f(t)) exp (f())
(2)
The function 9(-) is a positive function to guarantee the non- nog-atittyv of the intensity which we choose to be the exponential function, and our objective is to learn the function s()). The log-
2

---

Under review as a conference paper at ICLR 2021
likelihood of the multi-dimennional Poisson process with the functional prior is described as
N
log (ft:) N (tt) E - exp (f (t)) dt.
(3)
In the following sections, we introduce generalized additive models and propose to model by the log-linear model to learn f(t) and the normalizing term.
2.1 GENERALIZED ADDITIVE MODEL
In this section we present the connection between Poisson processes with Generalized Additive Model (GAM) proposed by Friedman & Stuetzle (1981). GAM projects higher-dimensions fea- tures into lower-dimensional space to apply smoothing functions to build a restricted class of non-parametric regression models. GAM is less affected by the curse of dimensionality com- pared to directly using smoothing in a higher-dimensional space. For a given set of processes (l....,D), D}. the traditional GAM using one-dimensional projections defined as
S [2] 1 log -je f(((6)) - with some smoothing function fj-
In this paper, extend it include higher-order interactions between features in GAM. The k-th order GAM is defined as
log (t)
Σ (u)(())))
(tlrs) ((22))
Se (166)) ,...,, (166)....,()))) By
jEJ
JUNEEJ
Σ 1((t))
- 3s=
(4)
IEJ, T|SA
where t(() E |R:41 denotes the subvector )EE of with respect to S [D]. The function
fi R|1 + is a smoothing function to fit the data, and the normalization constant By for the
intensity function is obtained as 33 = (t) exp( (EEp(E fi(t() ( )dt. The definition of the additive model is in the same form as Equation (3). In particular, if we compare Equation 3 and 44. we can see that the smoothing function f in (3) corresponds to the right-hand side of
Learning of a continuous function using lower dimensional projections is well known because the Kolmogorov-Arnole representation theorem, which states follows:
Theorem 1 (Kolmogorov-Aninn
Representation Theorem (Braun & Griebel 2009
Kol-
mogorov 1957)). Any multivariate continuous function can be represented as superposition of 2n+1 fq (Ep=199.0 (g. (tp)
one-dimensional functions, i.e., f ((,....,ta)) = E
Braun (2009) showed that the GAM is approximation to the general form presented in Kolmogorov- representation theorem by replacing the range {l,..., , 2n + 13 with C J and the inner function Iq. by the identity if g p and zero otherwise, yielding
f(t) f((D))
Interestingly, the canonical form for additive models in Equation 44 can be rearranged to be in the same form as Kolmogorov-Amole representation theorem. By letting f(t) = Eigif((t()) =
-- (A(t)) and g(-) exp(.), we have
Az(t) = exp (3r) exp (Eic (t x exp (Erg ((0)),
(5)
where we assume () 0 if [/] > for the k-th order model and 1/ (3)) is the normal- ization term for the intensity function. Based the Kolmogorov-Armold representation theorem, generalized additive models are able to learn the intensity of the higher-order interaction between point processes by using projections into lower dimensional space. The higlliilldood function for a kth-order model is obtained by substituting the Equation (4) into Equation (T.
(((( 3r,
log ((t)
(A
Σ exp
5
1=1
ICJ.
where is a constant given by 3' x(t)dt + Ey In the following section we will detail log-linear formulation that efficiently maximizes this log-likelihood equation.
3

---

Under review as a conference paper at ICLR 2021
2.2 ADDITIVE POISSON PROCESS
We introduce our key technical contribution in this section, the log-linear formulation of the additive Poisson process, and draw parallels between higher-order interactions in the log-linear model and the lower dimensional projections in generalized additive models. In the following, we discretize the time window [0, I' into M bins and treat each bin as natural number for each process. We assume that M is predetermined by the user. First we introduce a structured space for the Poisson process to incorporate interactions between processes. Let f2 (J, [J
[M]] 21....... M)
2(P) 0, T E [M]] } U (((,0))}. We define the partial order 3 (Davey & Priestley 2002 on as
(J. ( + J C and r S T', for each a = (J, = (J',') Ω (6)
and (1.0) 3 (J,T) for all (J.r) E Ω which illustrated in Figure [ The relation J JS J' is used to model any-order interactions between point processes (Luo Sugiyama
t= M
(0, =)
T=2
2019) (Amari, 2016) Section 6.8.4) and each т in (J.T) represents "time" in our model with
r=1
Order interaction
between processes
(1,2,3)
1 denoting the least element in the partial or-
der structure. Note that the domain of r can be
(1.2)
generalized from [M] [Mf] to take different
(3)
time stamps into account, while in the follow-
Time
ing we assume that observed time stamps are always the same across processes for simplic- ity. Our experiments in the next section demon- strates that we can still accurately estimate the density of processes. Our model can be applied not only time-series data but any sequential data.
Or=
Figure Partial order structured sample space (52, 3) with D = 3.
On any set equipped with partial order, we can introduce a log- model (Sugiyama et al.) 2016 2017). Given a parameter domain S - 52. For a partially ordered set (S2, 3).. the log-linear model with parameters (0.)xes is introduced as
logp(w; 0) Σ EES (s (0)
for each w E Ω where [ = 1 if the statement in [+] is true and 0 otherwise, and (0) E R is the par- tition function uniquely obtained as (0) = log Ewen exp( E.es )1ssw) - special case of this formulation coincides with the density function of the Boltzmann machines (Sugiyama
(et al. 2018 Luo & Sugiyama 2019).
Here we have clear correspondence between the log-linear formulation and that in the form of Kolmogorov-Aroo representation theorem in Equation (5) if we rewrite Equation 7 as
p(u; 0) = exp (0)
exp (Σ .ES (2 x exp Σ .ES 11
(8)
We call this model with (S2, >) defined in Equation (6) the additive Poisson process, which rep- resents the intensity À as the joint distribution across all possible states. The intensity of the multi-dimennionall Poisson process given via the GAM in Equation (5) is fully modeled (parameter- ized) by Equation 7 and each intensity fr(-) is obtained as (1,,)- To consider the k-th order model, we consistently use the parameter domain S given as = (J. E Ω | |J| < k: }. where is an input parameter to the model that specifies the upper bound of the order of interactions. This means that 0, = O for all s & S. Note that our model is well-defined for any subset S C Ω and the user use arbitrary domain in applications.
For a given J and each bin 7 with w = (J, the empirical probability p(w). which corresponds to the input observation, is given
N
1 Ear(T) Z I(w,, and (T( = Nlg ((w) 7
Σ K
hr
ICJ
WEB
im1
for each discretized state smoothing on time stamps ti...., ty. which is the kernel smoother proposed by Buja et al. (1989). The function K is kernel and hy is the bandwidth for each projection I S [D]. We use the Gaussian kernel as K to ensure that probability is always nonzero. meaning that the definition of the kernel smoother coincides with the kernel estimator of the intensity function proposed by |Schäbe (993)...
(J, where T = (T. E RP, The function Ø, performs
4

---

Under review as a conference paper at ICLR 2021
2.3 OPTIMIZATION
Given an empirical distribution ş defined in Equation ( the task is to learn the parameter (0,) (0)) such that the distribution via the log-linear model in Equation ( is close á as much as possible. Let us define 6 = (p 0g 0 if & S)) which is the set of distributions that can be represented by the log-linear model using the parameter domain S. Then the objective function is given as m DKL (p. p). where DKL (p. = Ewen P plog(p/p) /)) is the KL divergence from to p. In this optimization, let p" be the learned distribution from the sample with infinitely large sample size and p be the learned distribution for each sample. Then we can lower bound the uncertainty (variance) (pD p)] by |S|/2N (Barron & Hengartner, 1998).
Thanks to the well developed theory of information geometry (Amari 2016) for the log-linear model (Amari 200T), is known that this problem can solved by e-prejection, which coincides with the maximum likelihood estimation, and it is always convex optimization (Amari 2016 Chapter 2.8.3). The gradient with respect to each parameter 0g is obtained by (0/00, DKL (p, p) 7/ ne where 77, Euga ((wz))((w) The value 7/s is known as the expectation parameter (Sugiyamal et al. 2017) and 7/s is obtained by replacing with ý in the above equation. If s) 0 for some E S, remove from S to ensure that the model is well-defined.
Let = (s and the natural gradient (Amari 1998) the closed form solution of the Fisher information matrix is always available Sugiyama et al. 2017). The update step is Onext G G rt), where the Fisher information matrix G is obtained
[0......,88. 8. 1", = I "ais 17. We can always use
Pij = 00,,00, DKL(D.P)
( Ta, Tsg
(10)
WEn
Theoretically the Fisher information matrix is nu- merically stable to perform matrix inversion. However, computationally floating point errors
Algorithm 1 Additive Poisson Process (APP)
1: Function APP([t.) S, M, lt):
may cause the matrix become indefinite. To overcome this issue, small positive value is added along the main diagonal of the matrix. This technique is known as jitter and it is used in ar- eas like Gaussian processes to ensure that the co- variance matrix is computationally positive semi- definite (Neal 1999).
2: Initialize Ω with the number M of bins
3: Apply Gaussian Kernel with bandwidth
on [t:),= compute ò
4: Compute rj = (Ma) from ÿ
5: Initialize 0 = (02).e (randomly 0, 0)
6: repeat
7: Compute using the current (0x)ses
The pseudocode for APP is shown in Algo- rithm [ The time complexity of computing line
8:
Compute 77 (n.).s from
7 is ((2)))))) This means when implement- ing the model using gradient descent, the time up-
An + n-
10:
Compute the Fisher information matrix
complexity of the model is (( date the parameters in S for each iteration. For natural gradient the cost of inverting the Fisher information matrix G is O([S)"), therefore the time complexity update the parameters in is O(S)) + (2] for each iteration. The time complexity for natural gradient is significantly higher to invert the fisher information matrix, the number of parameter small, more efficient to use natural gradient because it requires significantly less iterations. However, if the number of parame- ters is large, more efficient to use gradient descent.
G using Equation (TO)
11:
i + 0 G 'nn
12:
until convergence of 0 = (0g)ses
13:
End Function
EXPERIMENTS
We perform experiments using two dimensional synthetic data, higher dimensional synthetic data, and rea- world data to evaluate the performance of our proposed approach. Our code is implemented on Python 3.7.5 with NumPy version 1.8.2 and the experiments are run on Ubuntu 18 04 LTS with an Intel 17-8700 6c/12t with 16GB of memory In experiments of synthetic data, we simulate
'The code is available in the supplementary material and will be publicly available online after the peer review process.
5

---

Under review as a conference paper at ICLR 2021
Procere
Raried beadwitt
L3 3.
Kdreal Beededt
Kome birenitt
Norell Birdwidh
Dense observations.
Proces:T1)
Tagit EL Sivergacc
(b) Sparse observations.
Figure 2: KL Divergence for four-order Poisson process.
Preesss
Presmes
Presmus It,
(a) Dense observations, h 0.4.
Fresess 12!
Fruses B!
Procest FL.
(b) Sparse observations, 0.2.
Figure 3: Intensity function of two dimensional processes. Dots represent observations.
random events using Equation ( We generate an intensity function using mixture of Gaus- sians, where the mean is drawn from a uniform distribution and the covariance is drawn from an inverted Wishart distribution. The intensity function is then the density function multiplied by the sample size. The synthetic data generated by directly drawing sample from the probability density function. An arbitrary number of samples is drawn from the mixture of Gaussians. We then run our models and compare with Kernel Density Estimation (KDE) (Rosenblat) 1956), an inhomogeneous Poisson process whose intensity is estimated by reproducing kernel Hilbert space formulation (RKHS) (Flaxman et al.) 2017), and a Dirichlet process mixture of Beta distributions (DP-beta) (Kottas 2006). The hyper-parameter M and h in our proposed model are selected us- ing grid search and cross-validation. For situations where a validation set is not available, then could be selected using rule of thumb approach such as Scott's Rule (Scott) (2015) and M could be selected empirically from the input data by computing the time interval of the joint observation.
3.1 EXPERIMENTS ON TWO-DIMENSIONAL PROCESSES
20 Gaussian components and simulate dense case with 100,000 obser- vations and a sparse case with 1,000 observations within the time frame of 10 seconds. We consider that a joint event occurs the two events occur 0.1 seconds apart. Figure and Figure2bjcomparee the KL divergence between the first- and second-orded models and Figure are the corresponding intensity functions. In the first-order processes, both first- and second-order models have the same performance. This is expected both of the model can treat first-order interactions and is able to learn the empirical intensity function exactly which is the superposition of the one-dimensional pro- jection of the Gaussian kernels each observation. For the second-order process, the second-order model performs better than the first-order model because it is able to directly learn the intensity function from the projection onto the wo-diimensional space. In contrast, the first-order model must approximate the second-order process using the observations from the first order-processes. In the sparse case, the second-order model performs better when the correct bandwidth is selected.
For our experiment,

---

Under review as a conference paper at ICLR 2021
Preeent
Procers 21
Frecem: D. A.A.M
Proesse 2.
(a) Dense observations.
Prcees: 131
masil_1 31
Proses 2.33
Pracest Ct, 2.4-4
Tase EL Zirepprcs
gg
6- Sww.Pa
L4 23
2
Kemaa lesdwid
Kamell bandwidd
(b) Sparse observations.
Figure 4: KL Divergence for four-order Poisson process.
113010
1400
(a) Dense observations.
(b) Sparse observations.
Figure 5: Intensity function of higher dimensional processes. Dots represent observations.
Table D compares our approach APP with other state-of-the-art approaches. APP performs the best for first-order processes in both the sparse and dense experiments. Experiments for RKHS and DP-beta were unable to complete running within 2 days for the dense experiment. In the second- order process our approach was outperformed by KDE, while both the second-order APP is able to outperform both RKHS and DP-beta process for both sparse and dense experiments. Figure 2a and Figure 2b show that KDE is sensitive to changes in bandwidth, which means that, for any practical implementation of the model, second-order APP with less sensitive bandwidth is more likely to learn more accurate intensity function when the ground truth is unknown.
3.2 EXPERIMENTS ON HIGHER-DIMENSIOON
PROCESSES
We generate fourth-order process to simulate the behaviour of the model in higher dimensions. The model generalizable to higher dimensions, however it is difficult to demonstrate results for processes higher than fourth-order. For our experiment, we generate an intensity function using 50 Gaussian components and draw a sample with the size of 107 for the dense case and that with the size of 1 for the sparse case. We consider the joint event to be the time frame of 0.1 seconds.
We were not able to run comparison experiments with other models because they are unable to learn when there are
few direct observations in third- and fourth-order processes. In addition, the time complexity is too high to learn from direct observations in first- and second-order processes because all the other models have their time complexity proportional the number of observations. The time complexity for KDE is (NN)) for the dimensionality with D, while DP-beta is O K), where K the number of clusters, and RKHS is (N for each iteration with respect to the sample size N. where De-beta and RKHS are applied to a single dimension as they cannot directly treat multiple dimensions. KDE is able to make an estimation of the intensity function when there are no direct observations of the simultaneous event, however, it was too computationally expensive to complete running the experiment. Differently, our model more efficient because the time com- plexity is proportional to the number of bins in our model. The time complexity of APP for each iteration is (( where |ol M and |S| Σ (D)- Our model scales combinatorially

---

Under review as a conference paper at ICLR 2021
The lowest KL divergence from the ground truth distribution to the obtained distri- bution on two types of single processes ([1] and (2I) and joint process of them ([1.2)). APP-# represents the order the Additive Poisson Pro- cess. Missing values mean that the computation did not finish within two days.
Table 1:
Table Negative test log-likelihood for the New York Taxi data. Single processes ([T) and [W]) and joint process of them (IT.WD. APP-# repre- sents the order of the Additive Poisson Process.
Process
APP-1
APP-2
RKHS
DP-bcta
[T]
714.07
714.07
713.77
728.13
731.01
Jan
[W]
745.60
745.60
745.23
853.42
790.04
Process
APP-1
APP-2
KDE
RKIS
bet
[T.W]
249.60
246.05
380.22
259.29
260.30
[1]
4.98
4.%c-5
2.8ee4
[T]
713.43
713.43
755.71
795.61
765.76
Dense
[2]
2.83-5
2.83-5
1.17e4
Feb
738.66
738.66
773.65
811.34
792.10
[1.2]
2.98c-2
1.27e-3
6.33e-4
4.09e-2
4.54e-2
[T.W]
328.84
244.21
307.86
334,31
326.52
[1]
7.26c-4
7.26e-4
8.83e-4
1.96e
2.62e-3
716.72
716.72
733.74
755.48
741.28
Sparse
[2]
2,28c-4
2.28c-4
2.76e-4
2.36e-3
2.49c-3
Mar
[W]
738.06
738.06
816.99
853.33
832,43
[1.2]
2.88822
1.77622
3,67e-3
188c-2
3.680-2
FT.11
291.20
246.19
289.69
328.47
300.36
with respect to the number of dimensions. However, this is unavoidable for any model which directly takes into account the high-order interactions, For practical applications, the number of dimensions D and the order of the model & is often small, making it feasible compute.
In Figure 4a we observe similar behaviour in the model, where the first-order processes fit precisely to the empiriall distribution generated by the Gaussian kernels. The third-order model is able to pe- riod better on the fourth-order process. This is because the observation shown in Figure Sajis largely sparse and learning from the observations directly may overfit. A lower dimensional approximation is able to provide a better result in the third-order model. Similar trends can be seen in the sparse case as shown in Figure 4b where a second- -oder model is able to produce better estimation in third- and fourth-order processes. The observations are extremely sparse seen in Figure 5b where there are only a few observations or no observations at all to learn the intensity function.
3.3 UNCOVERING COMMON PATTERNS IN THE NEW YORK TAXI DATASET
We demonstrate the capability of our model on the 2016 Green Taxi Trip datase We are interested in finding the common pick up patterns across Tuesday and Wednesday. We define a common pick up time to be within minute intervals of each other between the two days. We have chosen to learn an intensity function using the Poisson process for Tuesday and Wednesday and a joint process for both of them. The joint process uncovers the common pick up patterns between the two days. We have selected to use the first two weeks of Tuesday and Wednesday January 2016 as our training and validation sets and Tuesday and Wednesday of the third week of January 2016 as our testing set. We repeat the same experiment for February and March.
We show our results in Table 2 where we use the negative test log- ligllikeli as an evaluation measure. APP-2 has consistently outperformed all the other approaches for the joint process between Tuesday and Wednesday. In addition, for the individual process, APP-1 and -2 also showed the best result for February and March. These results demonstrate the effectiveness of our model in capturing higher-order interactions between processes, which is difficult for the other existing approaches.
CONCLUSION
We have proposed a novel framework, called Additive Poisson Prucess (APP), to learn the intensity function of the higher-order interaction between point processes using samples from lower dimen- sional projections. We formulated our proposed model using the the log-lineat model and optimize it using information geometric structure of the distribution space. We drew parallels between our proposed model and generalized additive model and showed the ability to learn from lower dimen- sional projections via the Kolmogorov-Arnold representation theorem. Our empirical results show the superiority of our method when learning the higher-order interactions between point processes when there are no or extremely sparse direct observations, and our model is also robust to varying sample sizes. Our approach provides novel formulation to learn the joint intensity function which typically has extremely low intensity. There enormous potential to apply APP to real-world appli- cations, where higher order interaction effects need to be model such as in transportation, finance, ecology, and violent crimes.
8

---

Under review as a conference paper at ICLR 2021
REFERENCES
Alan Agresti. Categorical Data Analysis. Wiley, 3 edition, 2012.
S. Amari. Information geometry on hierarchy of probability distributions. IEEE Transactions on Information Theory, 47(5):1701-1711, 2001.
Shun-Ichi Amari. Natural gradient works efficiently learning. Neural Computation, 10(2):251- 276, 1998.
Sun-Ichi Amari. Information Geometry and Its Applications. Springer, 2016.
A. Barron and N. Hengartner. Information theory and superefficiency. The Annals of Statistics, 26 (5): 15:11800-825, 1998.
Jürgen Braun. An application of Kolmogorov's superposition theorem to function reconstruction in higher dimensions. PhD thesis, Universitäts Landesbibliothek Bonn, 2009.
Jürgen Braun and Michael Griebel. On a constructive proof of Kolmogorov's superposition theorem. Constructive Approximation, 30(3):653, 2009.
Andreas Buja, Trevor Hastie, and Robert Tibshirani. Linear smoothers and additive models. The
Annals of Statistics, 17(2):453-510, 1989.
Daryl Daley and David Vere-Jones. An Introduction to the Theory of Point Processes: Volume II:
General Theory and Structure. Springer, 2007.
Brian A Davey and Hilary A Priestley. Introduction to Lattices and Order. Cambridge University Press, 2002.
Seth Flaxman, Yee Whye Teh, and Dino Sejdinovic. Poisson intensity estimation with reproducing kerells. Electronic Journal of Statistics, 11(2)50881-5104 2017.
Jerome H Friedman and Werner Stuetzle. Projection pursuit regression. Journal of the American Statistical Association, 76(3766:817-823, 1981.
Deniz Ilalan. A poisson process with random intensity for modeling financial stability. The Spanish Review of Financial Economics, 14(2):33500 2016.
Andrei Nikolaevich Kolmogorov. On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition. Doklady Akademit Nauk, 114(5):953 9966,,95656111 1957.
Athanasios Kottas. Dirichlet process mixtures of beta distributions, with applications to density and intensity estimation. In Workshop on Learning with Nonparametric Bayesian Methods, 23rd International Conference on Machine Learning (ICML), volume 47, 2006.
Simon Luo and Mahito Sugiyama. Bias-variance trade-off in hierarchical probabilistic models us- ing oigher-order feature interactions. Proceedings of the 33rd AAAI Conference on Artificial Intelligence, pp. 1488-4495, 2019.
Radford M Neal. Regression and classification using gaussian process priors. Bayesian Statistics,
6:475-501, 1999.
Yosihiko Ogata. On Lewis' simulation method for point processes. IEEE Transactions on Informa- fion Theory, 27(1):33-31 1981.
Murray Rosenblatt. Remarks on some nonparametric estimates of density function. The Annals of Mathematical Statistics, pp. 832-837, 1956.
H Schäbe. Nonparametric estimation of intensities of nonhomogeneous poisson processes. Statisti- cal Papers, 34(1):113-131, 1993.
David W Scott. Multivariate density estimation: theory, practice. and visualization. John Wiley & Sons, 2015.

---

Under review as a conference paper at ICLR 2021
Mahito Sugiyama, Hiroyuki Nakahara, and Koji Tsuda. Information decomposition structured space. In 2016 IEEE International Symposium on Information Theory, pp. 575-579, 2016.
Mahito Sugiyama, Hiroyuki Nakahara, and Koji Tsuda. Tensor balancing on statistical manifold. In Proceedings of the 34th International Conference Machine Learning, volume 70 of Proceed- ings of Machine Learning Research, pp. 3270-3279, 2017.
Mahito Sugiyama, Hiroyuki Nakahara, and Koji Tsuda. Legendre decomposition for tensors. In Advances in Neural Information Processing Systems 31, pp. 8825-8828 8835, 2018.
Matthew A Taddy, Autoregressive mixture models for dynamic spatial poisson processes: Applica- tion to tracking intensity of violent crime. Journal of the American Statistical Association, 105 (492):1403-1417, 2010.
HR Thompson. Spatial point processes, with applications to ecology. Biometrika, 42(1/2):102-115, 1955.
Feng Zhou, Zhidong Li, Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen. refined MISD algorithm based on Gaussian process regression. In Pacific-Asia Conference on Knowledge Dis- covery and Data Mining, pp. 584-596. Springer, 2018.
A ADDITIONAL EXPERIMENTS
A.1 BANDWIDTH SENSITIVITY ANALYSIS
Our first experiment to demonstrate the ability for our proposed model to learn an intensity func- tion from samples. We generate Bernoulli process with probably of p 0.1 to generate samples for every 1 seconds for 100 seconds to create toy problem for our model. This experiment is to observe the behaviour of varying the bandwidth in our model. In Figure 6a| we observe that apply- ing no kemel,, we learn the deltas of each individual observation. When we apply Gaussian kernel, the output of the model for the intensity function is much more smooth. Increasing the bandwidth of the kernel will provide a wider and much smoother function. Between the 60 seconds and 80 seconds mark, can be seen when two observations have overlapping kernels, the intensity function becomes larger in magnitude.
A.2 ONE DIMENSIONAL POISSON PROCESS
A one dimensional experiment is simulated using Ogata's thinning algorithm (Ogata 1981). We gen- erate two experiments use the standard sinusoidal benchmark intensity function with frequency of 20T. The dense experiment has troughs with O intensity and peaks at 201 and the sparse experiment has troughs with intensity and peaks case, our model has no problem learning the intensity function. We compare our results using KI. divergence between the underlying intensity function used to generate the samples to the intensity function generated by the model. Figure 6b shows that the optimal bandwidth is h 1.
Figure 6d shows the experimental results of the dense
10

---

Under review as a conference paper at ICLR 2021
No Karnel
0.125
0.3
0.8
C9
Divergance
0.6
0.075
Karnel,
0.050
0.025
0.0
0.000
-0.2
20
40
60
8.0
160
Time
(a) Toy Example
(b) KL divergence dense experi- (c) KL divergence sparse experi-
ment
ment
(d) Ogata's 's thinning algorithm with high intensity
c Ogata's 's thinning algorithm with low intensity
Figure 6: One dimensional experiments
Algorithm Thinning Algorithm for non-homogenous Poisson Process
1: Function Thinning Algorithm (A (t), T):
2: 0, so 0, X suposegr>
3: repeat
4:
uniform (0, 1)
5: w - In u {w exponential()) {)}
6: Sm+1 9m +
7: D ~ uniform (0,1)
if D M then
8:
X
9:
tn+1 = 8m+1
==n+
11: else
m m
13:
end if
if In S T then
15:
return {4x (1
16: else
return {t+ } n-
end if
19:
until 8m <T
End Function
11

---

Under review as a conference paper at ICLR 2021
0200 12 E.5 Oi E3R 33 D5 1/5 2.1
070 12 eaE 075 SAR 124 100 17% TAN
Kamd LO 175 230
125 13 07, 1HE 13. 15 175 2.18
0N 14E LAN AAN 3/2 7
Procese (4)
Harmil Banteett
Precess: 31
Total EL Dierrpncs
MITH Randwatt
Proces 14
Bandeelts
Pr 31
11. 21,2,44
0.0 LS Bern Fandwatth
Tetal EL pireppnns
KarT Bundwatth
0.0 L2 IS
0.0 5:2
0.0
22
25
TL
0.00 LI E3p 075 331 135 13o 175 2
528 aRe 675 125 1.0 377h 2
0.00 0.25 83 075 33x 125 15 175 3.70
175
040 2.R 0A0 6A 1/8 AAN 1.2%
175
Figure 7: KL Divergence for four four-order Poisson process.
Procese EX
Preceme: E1, 41
Proteet: (1 2,
Bernd antwaats
BRIT Aandrana
Procese D3
Precuss: 41
Proce:: D1. 2,
030 12 Berne Sandwatt
Process: 11
BITH Sandwitta
Frocess: I1,
040 D2 DS
Dense observations.
(b) Sparse observations.
040
TM
Onter
0.30 OZ 03g tis Sa1 225 13o 17s 2.00
0l0 0.28 aND 6.5 1i 125 L.0 17% 2
0,30 0.22 030 875 3at 125 13o 175 3,00
040 FUN LA THE 1.4 1A0 1.7 2
Romel mondoiith 175
Orter
Prscosn D23
Precus:: It, 31
Precess (1 41
ProsEs:: 3,4)
Gomel Mandwilth
Precus: [t, 31
Precuss D. 41
femel mindiiittt
Fronss:: 3,41
Kamel Manovitta
Gan
Orden
Ordmn
Cooeara TM
TX
130 LT 2.00
6:3 6:28 8ND 0.5 1 122 1.0 677 2.00
6J3 6.2 830 075 3al 123 11e L75 2.00
6d8 FUR ADD ONW 1.2 144 1.7 2200
LT9
OJ BI Kemell ainduiihi LA
Parcoso E11
6.2 6.22 B30 0.5 591 125
Precuse: (4, 21
Proous: (2. 41
0J 6.22 Nemel Bindaiiih
Procuss (1, 3.44)
Kamee mnndiii
Porrac DI
Precuss: L1, 21
ProoE:: D2, 41
Precazs (t. 3.4)
Kamal Handutith
12

---

Under review as a conference paper at ICLR 2021
10
Onlerr
Proces: (4)
[2,3)
Process: fL, 11.2244
Process
Time
Time
Ordenn
Procemm:
Time
Time
Orden Orden
Figure 8: Intensity function of higher dimensional processes. Dots represent observations
100000
200
22
1000
500
30
1.0
Process: 12,22,3 31
Procersi TL. 2,3,3.4 4,
Precess [3]
Processo C1,
Processs 2.33
Process
Time
Time
Process:
Time
Time
Process:
(a) Dense observations.
(b) Sparse observations
2090
0.00
Precuss (4)
[3,4]
Processi F2. 3. 41
(A. 31
Process [3 41
Process 41
Precess:
Time
Time
Proceur:
Time
Time
Time
Time
4000
2000
1000
3
Agsuaguu
Pro-esso rb, 47
Proceas: [A]
+]
Process: [1.
Process: 12. 41
Time
Process:
Time
Time
Timee
Process:
Time
Time
Proces::
Time
2000
Aesuagn
0.29
0.09
Aspuagua
Agsusqui
13


Reviews:
[{"review_id": "v8pbpRzocBO", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The paper under review proposes a new model for multi-dimensional temporal Point processes, allowing efficient estimation of high order interactions. This new model, called additive Poisson process, relies on a log-linear structure of the intensity function that is motivated thanks to the Kolmogorov-Arnold theorem. Such structure is then linked to generalized additive models, a result that is used to devise an efficient estimation procedure with formal guarantees  of convergence.\n\nThe paper is well written but some part lacks of details making the exposition sometime difficult to follow. Overall the contribution is original and promising, but in its current state suffers from imprecision making difficult to assess its quality. The following points summarizes these remarks:\n\n\n## Major points\n\n* The experimental design is unclear: what is the structure of the intensity function? Is it log-linear? If yes, why APP, which is designed for this case is outperformed by KDD? Why is the interaction between processes 1 and 2, is summarized by a univariate function ? Does this imply that only simultaneous interactions are considered in this work, as suggested by the second sentence of Section 3.1? \n* How would the method compare to others when the log-linearity is not satisfied, e.g., $\\lambda(t_1,t_2) = \\lambda_1(t_1)\\lambda_2(t_2)$? Such cases, if not studied in detail, should at least be mentioned and discussed as a limitation of the current model.\n* The application is used to showcase the usefulness of the new model, but model comparison relies on a negative log-likelihood criteria which does not account for the complexity of the models. If possible, some more general criteria such as AIC should be considered. If such calculations are not possible/does not apply, then differences in model's complexity should be at least discussed.\n\n## Minor points\n\n* p.1: the first sentence is really unclear. Not all stochastic process have an 'arrival time'. Note that correlation does not mean simultaneity of the events.\n* p.1: 'intensity function is multiplied ...'\n* p.1 : twice 'the'.\n* p,1 'Differently' sounds awkward.\n* p.1: 'if observations are sparse'.\n* p.2: 'intensity is often low' would be better to say 'lower than lower-orders'.\n* p.2: of the model 'are'.\n* p.2: 'the model is robust to' this is not clear and never discussed in the paper. No misspecified experiment is presented.\n* p.2 'arrival intensity'?\n* p.2: 'for time t > 0, ...' is this sentence needed? Why not simply giving the classical setting with integration over a non-negligible region?\n* p.3: the notation $g_{p,q}$ is never defined.\n* p.3: the link between Theorem 1 and the log-linear model is unclear: how the two summations become only one? The summation in the Theorem uses the number of sample $n$ while the summation in (5) depends on the dimensionality of the subset. This is particularly unclear.\n* p.3: notation $\\perp$ is not defined.\n* p.4: the variable $Z$ is defined but never used.\n* p5.: 'We can always use the natural gradient ...' this sentence is unclear.\n* p.5: 'uniform distribution' on $[0,1]$? What are the 'expected number of samples' that are used?\n* p.8: Table b: What are $T$ and $W$? How is it possible to have a joint process for APP-1?\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "o3BsNZZ568E", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "## Summary\n\nThis work claims to address the problem of learning the structure of interactions between multiple point processes. This can, e.g., be used to predict co-occurrences between different types of events. The authors develop a method that (a) is theoretically principled and (b) scales to large datasets.\n\n## Reasons for score\n\nThe overall problem is clearly relevant, and some aspects of the contribution appear to be promising (connection to GAMs, approach to inference). However, in its current form, the paper is excessively unclear, to the point where I am not sure that I understand the precise problem the authors are trying to address (despite my familiarity with the area).\n\n## Pros\n\n- if the problem is that of understanding correlations between multiple point processes, then that is an important problem that is clearly relevant to the ML community.\n- The approach builds on generalized additive models, which gives their approach a solid theoretical foundation\n- The inference algorithm based on minimizing KL-divergence appears to be interesting.\n\n## Cons\n\n- the paper lacks clarity throughout, making it difficult to understand the problem, let alone the contribution itself (examples below).\n- the experiments are not explained in sufficient details (NYC taxi dataset for example), and as such they are not supporting the contribution.\n\n## Comments\n\nReferring to \"stochastic processes\" is too general—maybe talking about \"point processes\" would be better. For example, a Gaussian process is a stochastic process, but doesn't have anything to do with the processes considered in this work.\n\nSection 2, \"Formulation\" is confusing. Usually, a multi-dimensional / spatial poisson process is a *single* point process defined on $\\mathbf{R}^D$, and a realization of that process is a collection of points in D-dimensional space. This is what, e.g., Flaxman et al. (2017) consider.\n\nIn this paper, the authors seem to assume that every dimension is a different process, and somehow events in each dimension are grouped together based on their arrival order (\"Each $t_i$ is the time of occurrence for the $i$-th event across $D$ processes and $T$ is the observation duration\"). I am not even sure this is always possible, since it seems to assume that every process / dimension has exactly $N$ events during the interval $[0, T]$. Thus, I fail to understand even the basic foundation on which their work builds.\n\nIt seems to me that the authors are really considering a *multivariate* Poisson process, where the intensity function of each process are dependent. In that case, the likelihood would be different from (1). Can the authors comment on this?\n\nExperiments:\n\n- \"the sample size is randomly chosen by the mixture of Gaussians\" - what does this mean? a GMM does not define a distribution over number of samples.\n- where is the notion of timestamps in the synthetic GMM experiments?\n- Figure 3 does not seem to be referenced in the text. How can the intensity function be 1D for if the experiments are with 2D Gaussians?\n- NYC taxi dataset: how is space taken into account? Defining a \"common pick-up\" as one that happened within at 1 min interval but at opposite ends of the city seems pointless.\n\nOther comments / questions:\n\n- How does your work compare / contrast to multivariate Hawkes processes?\n- In what way is Flaxman et al. 2017 non-convex?\n- \" ...where the different events are marked so that a particular subset of events can be used to generate the intensity function\" what does this mean?\n\n## POST-REBUTTAL\n\nThank you for rebuttal which has helped me, to some extent, to understand your model. I still believe the formulation needs to be clarified before the contributions can be assessed objectively.\n\nFor example, on page 2:\n\n> Given a realization of timestamps $\\mathbf{t}_1, \\mathbf{t}_2, \\ldots ,\\mathbf{t}_N$ with $\\mathbf{t}_i \\in [0, T]^D$ from an inhomogeneous (multi-dimensional) Poisson process with the intensity $\\lambda$. Each $\\mathbf{t}_i$ is the time of occurrence for the $i$-th event across $D$ processes and $T$ is the observation duration.\n\nThis is confusing: only a single one of the dimensions will likely represent \"time\". Other dimensions will be space, etc. Furthermore, why would every dimension need to be within $[0, T]$? Different dimensions might use different units and be in different ranges.\n\nI am upgrading my score slightly, because I sense that your contributions might be interesting once properly explained, but in the present state I believe the paper is not ready for publication.", "extra_scores": {}, "extra_text": {}}, {"review_id": "bGN3uwrEQP", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper proposes a generalized additive model to learn joint intensity functions for multiple Poisson processes. The main contribution is to consider the partial order of high-order interactions when parameterizing the model.\n\nThe proposed method is a reasonable solution, but there are several points in the paper that are particularly weak:\n\nFirst, the paper does not find a proper scenario to apply the additive Poisson process model. The experimental part lacks real data analysis. In Section 3.3, they only show a quantitative result regarding two Poisson processes [T, W]. It is hard to summarize that the performance gain comes from high-order modeling without explaining it in detail. For example:\ni) Why are KDE and RKHS predictions inconsistent across January through March?\nii) How do you define M?\niii) How sparse is the data?\n\nSecond, jumping from section 2.1 to 2.2 is confusing, since section 2.2 seems to have nothing to do with the Poisson process. So my first question is:\ni) Can you define the likelihood in a single Poisson process model in order to consider all the interactions you mentioned in section 2.2?\nAfter carefully reading the setup in section 2.2, I feel that the entire method is based on a crude approximation of Poisson processes. For example, I’m curious about:\nii) Does your method really depend on discretizing the time window? What if you let M goes to infinity? I’m particularly interested in this question because the original Poisson process likelihood model does not require discretization.\n\nTo conclude, this paper proposes a reasonable method to model high-order correlations. However, its experimental result is weak and the model construction is of a question to me.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes a method for modeling higher-order interactions in Poisson processes. Unfortunately, the reviewers do not feel that the paper, in its current state, meets the bar for ICLR. In particular, reviewers found the descriptions unclear and the justifications lacking. While the responses did aid the reviewers understanding, the paper would benefit from rewriting and more careful thought given to the experimental design.
Decision: Reject