Year: 2020

Abstract:
We present the Additive Poisson Process (APP), a 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 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 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.

text: We present the Additive Poisson Process (APP), a 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 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 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.

heading: 1 INTRODUCTION
text: Consider two point processes which are correlated with arrival times for an event. For a 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 a novel solution to this problem.
The Poisson process is a counting process used in a wide range of disciplines such as time-space sequence data including transportation (Zhou et al., 2018), finance (Ilalan, 2016), ecology (Thompson, 1955), and violent crime (Taddy, 2010) to model the arrival times for a single system by learning an intensity function. For a 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 correlation 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 two events.
For example in a spatial-temporal problem, we want to learn the intensity function for a taxi to pick up customers at a given time and location. For this problem, each point is multi-dimensional, that is (x, y, t)Ni=1, where a pair of x and y represents two spatial dimensions and t represents the time dimension. For any given location or 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) (Rosenblatt, 1956) are able to learn the joint intensity function. However, KDE suffers from the curse of dimensionality, which means that KDE requires a large size sample or a high intensity function to build an accurate model. In addition, 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 distributions with a 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 it 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.
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 simultaneous 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) in generalized additive models to learn the intensity function using samples in a lower dimensional space. Our proposed approach is to decompose a multi-dimensional point process into lower-dimensional representations. For example, in the x-dimension we have points (xi)Ni=1, in the y-dimension, we have points (yi) N i=1 and in the time dimension we have (ti)Ni=1. The data in these lower dimensional space can be used to improve the estimate of the joint intensity function. This is different from the traditional approach where we only use the simultaneous events to learn the joint intensity function.
We first show the connection between generalized additive models and Poisson processes. We then provide the connection between generalized additive models and the log-linear model (Agresti, 2012), which has a well-established theoretical background in information geometry (Amari, 2016). We draw parallels between the formulation of the generalized additive models and the binary loglinear 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 a unique optimal solution using natural gradient, which minimizes the Kullback-Leibler (KL) divergence from the sample distribution in a 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 Kolmogorov-Arnold representation theorem. 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 is 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.

heading: 2 FORMULATION
text: In this section we first introduce the technical background in the Poisson process and its extension to a multi-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 a 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 λ:RD → R, where we assume multiple D processes. An inhomogeneous Poisson process is a general type of processes, where the arrival intensity changes with time. The process with time-changing intensity λ(t) is defined as a counting process N(t), which has an independent increment property. For all time t ≥ 0 and changes in time δ ≥ 0, the probability p for the observations is given as p(N(t+ δ)−N(t) = 0) = 1− δλ(t) + o(δ), p(N(t + δ)− N(t) = 1) = δλ(t) + o(δ), and p(N(t + δ)− N(t) ≥ 2) = o(δ), where o(·) denotes little-o notation (Daley & Vere-Jones, 2007). Given a realization of timestamps t1, t2, . . . , tN with ti ∈ [0, T ]D from an inhomogeneous (multi-dimensional) Poisson process with the intensity λ. Each ti 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 ( {ti}Ni=1 | λ (t) ) = exp ( − ∫ λ (t) dt ) N∏ i=1 λ (ti) , (1)
where t = [t(1), . . . , t(D)] ∈ RD. We define the functional prior on λ(t) as
λ(t) := g (f(t)) = exp (f(t)) . (2)
The function g(·) is a positive function to guarantee the non-negativity of the intensity which we choose to be the exponential function, and our objective is to learn the function f(·). The log-
likelihood of the multi-dimensional Poisson process with the functional prior is described as log p ( {ti}Ni=1 | λ (t) ) = N∑ i=1 f (ti)− ∫ exp (f (t)) dt. (3)
In the following sections, we introduce generalized additive models and propose to model it by the log-linear model to learn f(t) and the normalizing term.

heading: 2.1 GENERALIZED ADDITIVE MODEL
text: In this section we present the connection between Poisson processes with Generalized Additive Model (GAM) proposed by Friedman & Stuetzle (1981). GAM projects higher-dimensional features 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 compared to directly using smoothing in a higher-dimensional space. For a given set of processes J ⊆ [D] = {1, . . . , D}, the traditional GAM using one-dimensional projections is defined as log λJ(t) = ∑ j∈J fj(t (j))− βJ with some smoothing function fj .
In this paper, we extend it to include higher-order interactions between features in GAM. The k-th order GAM is defined as log λJ(t) = ∑ j∈J f{j}(t (j)) + ∑ j1,j2∈J f{j1,j2}(t (j1), t(j2)) + · · ·+ ∑ j1,...,jk∈J f{j1,...,jk}(t (j1), . . . , t(jk))− βJ
= ∑
I⊆J, |I|≤k
fI(t (I))− βJ , (4)
where t(I) ∈ R|I| denotes the subvector (t(j))j∈I of t with respect to I ⊆ [D]. The function fI : R|I| → R is a smoothing function to fit the data, and the normalization constant βJ for the intensity function is obtained as βJ = ∫ λJ(t)dt = ∫ exp( ∑ I⊆J, |I|≤k fI(t
(I)) )dt. The definition of the additive model is in the same form as Equation (3). In particular, if we compare Equation (3) and (4), we can see that the smoothing function f in (3) corresponds to the right-hand side of (4).
Learning of a continuous function using lower dimensional projections is well known because of the Kolmogorov-Arnold representation theorem, which states as follows: Theorem 1 (Kolmogorov–Arnold Representation Theorem (Braun & Griebel, 2009; Kolmogorov, 1957)). Any multivariate continuous function can be represented as a superposition of one–dimensional functions, i.e., f (t1, . . . , tn) = ∑2n+1 q=1 fq (∑n p=1 gq,p (tp) ) .
Braun (2009) showed that the GAM is an approximation to the general form presented in Kolmogorov-Arnold representation theorem by replacing the range q ∈ {1, . . . , 2n + 1} with I ⊆ J and the inner function gq,p by the identity if q = p and zero otherwise, yielding f(t) = ∑ I⊆J fI(t (I)).
Interestingly, the canonical form for additive models in Equation (4) can be rearranged to be in the same form as Kolmogorov-Arnold representation theorem. By letting f(t) = ∑ I⊆J fI(t (I)) = g−1(λ(t)) and g(·) = exp(·), we have
λJ(t) = 1
exp (βJ) exp (∑ I⊆J fI ( t(I) )) ∝ exp (∑ I⊆J fI ( t(I) )) , (5)
where we assume fI(t(I)) = 0 if |I| > k for the k-th order model and 1/ exp(βJ) is the normalization term for the intensity function. Based on the Kolmogorov-Arnold 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 log-likelihood function for a kth-order model is obtained by substituting the Equation (4) into Equation (1),
log p ( {t}Ni=1|λ (t) ) = N∑ i=1 exp  ∑ I⊆J, |I|≤k fI ( t(I) )− β′,
where is a constant given by β′ = ∫ λ(t)dt + ∑ I⊆J βJ . In the following section we will detail a log-linear formulation that efficiently maximizes this log-likelihood equation.

heading: 2.2 ADDITIVE POISSON PROCESS
text: 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, T ] into M bins and treat each bin as a natural number τ ∈ [M ] = {1, 2, . . . ,M} 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 Ω = { (J, τ)|J ∈ 2[D] \ ∅, τ ∈ [M ] } ∪ {(⊥, 0)}. We define the partial order (Davey & Priestley, 2002) on Ω as
(J, τ) (J ′, τ ′) ⇐⇒ J ⊆ J ′ and τ ≤ τ ′, for each ω = (J, τ), ω′ = (J ′, τ ′) ∈ Ω, (6)
and (⊥, 0) (J, τ) for all (J, τ) ∈ Ω, which is illustrated in Figure 1. The relation J ⊆ J ′ is used to model any-order interactions between point processes (Luo & Sugiyama, 2019) (Amari, 2016, Section 6.8.4) and each τ in (J, τ) represents “time” in our model with ⊥ denoting the least element in the partial order structure. Note that the domain of τ can be generalized from [M ] to [M ]D to take different time stamps into account, while in the following we assume that observed time stamps are always the same across processes for simplicity. Our experiments in the next section demonstrates that we can still accurately estimate the
density of processes. Our model can be applied to not only time-series data but any sequential data.
On any set equipped with a partial order, we can introduce a log-linear model (Sugiyama et al., 2016; 2017). Given a parameter domain S ⊆ Ω. For a partially ordered set (Ω, ), the log-linear model with parameters (θs)s∈S is introduced as
log p(ω; θ) = ∑
s∈S 1[s ω]θs − ψ(θ) (7)
for each ω ∈ Ω, where 1[·] = 1 if the statement in [·] is true and 0 otherwise, and ψ(θ) ∈ R is the partition function uniquely obtained as ψ(θ) = log ∑ ω∈Ω exp( ∑ s∈S 1[s ω]θs ) = −θ(⊥,0). A special case of this formulation coincides with the density function of the Boltzmann machines (Sugiyama et al., 2018; Luo & Sugiyama, 2019).
Here we have a clear correspondence between the log-linear formulation and that in the form of Kolmogorov-Arnold representation theorem in Equation (5) if we rewrite Equation (7) as
p(ω; θ) = 1
expψ(θ) exp (∑ s∈S 1[s ω]θs ) ∝ exp (∑ s∈S 1[s ω]θs ) . (8)
We call this model with (Ω, ) defined in Equation (6) the additive Poisson process, which represents the intensity λ as the joint distribution across all possible states. The intensity λ of the multi-dimensional Poisson process given via the GAM in Equation (5) is fully modeled (parameterized) by Equation (7) and each intensity fI(·) is obtained as θ(I,·). To consider the k-th order model, we consistently use the parameter domain S given as S = { (J, τ) ∈ Ω | |J | ≤ k }, where k is an input parameter to the model that specifies the upper bound of the order of interactions. This means that θs = 0 for all s /∈ S. Note that our model is well-defined for any subset S ⊆ Ω and the user can use arbitrary domain in applications.
For a given J and each bin τ with ω = (J, τ), the empirical probability p̂(ω), which corresponds to the input observation, is given as
p̂(ω) = 1
Z ∑ I⊆J σI(τ ), Z = ∑ ω∈Ω p̂(ω), and σI(τ ) := 1 NhI N∑ i=1 K
( τ (I) − t(I)i
hI
) (9)
for each discretized state ω = (J, τ), where τ = (τ, . . . , τ) ∈ RD. The function σI performs smoothing on time stamps t1, . . . , tN , which is the kernel smoother proposed by Buja et al. (1989). The functionK is a kernel and hI is the bandwidth for each projection I ⊆ [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 Schäbe (1993).

heading: 2.3 OPTIMIZATION
text: Given an empirical distribution p̂ defined in Equation (9), the task is to learn the parameter (θs)s∈S such that the distribution via the log-linear model in Equation (7) is close to p̂ as much as possible. Let us define SS = {p | θs = 0 if s 6∈ 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 minp∈SS DKL(p̂, p), where DKL(p̂, p) = ∑ ω∈Ω p̂ log(p̂/p) is the KL divergence from p̂ 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) E[DKL(p∗, p)] by |S|/2N (Barron & Hengartner, 1998). Thanks to the well developed theory of information geometry (Amari, 2016) for the log-linear model (Amari, 2001), it is known that this problem can be solved by e-projection, 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 θs is obtained by (∂/∂θs)DKL(p̂, p) = ηs − η̂s, where ηs = ∑ ω∈Ω 1[ω s]p(ω). The value ηs is known as the expectation parameter (Sugiyama et al., 2017) and η̂s is obtained by replacing p with p̂ in the above equation. If η̂s = 0 for some s ∈ S, we remove s from S to ensure that the model is well-defined. Let S = {s1, . . . , s|S|} and θ = [θs1 , . . . , θs|S| ]T , η = [ηs1 , . . . , ηs|S| ]T . We can always use the natural gradient (Amari, 1998) as the closed form solution of the Fisher information matrix is always available (Sugiyama et al., 2017). The update step is θnext = θ −G−1(η − η̂), where the Fisher information matrix G is obtained as
gij = ∂
∂θsi∂θsj DKL(p̂, p) = ∑ ω∈Ω 1[ω si]1[ω sj ]p(ω)− ηsiηsj . (10)
Algorithm 1 Additive Poisson Process (APP) 1: Function APP({ti}Ni=1, S, M , h): 2: Initialize Ω with the number M of bins 3: Apply Gaussian Kernel with bandwidth h
on {ti}Ni=1 to compute p̂ 4: Compute η̂ = (η̂s)s∈S from p̂ 5: Initialize θ = (θs)s∈S (randomly or θs = 0)
6: repeat 7: Compute p using the current θ =
(θs)s∈S 8: Compute η = (ηs)s∈S from p 9: ∆η ← η − η̂
10: Compute the Fisher information matrix G using Equation (10) 11: θ ← θ −G−1∆η 12: until convergence of θ = (θs)s∈S 13: End Function
Theoretically the Fisher information matrix is numerically stable to perform a matrix inversion. However, computationally floating point errors may cause the matrix to become indefinite. To overcome this issue, a small positive value is added along the main diagonal of the matrix. This technique is known as jitter and it is used in areas like Gaussian processes to ensure that the covariance matrix is computationally positive semidefinite (Neal, 1999).
The pseudocode for APP is shown in Algorithm 1. The time complexity of computing line 7 is O(|Ω||S|). This means when implementing the model using gradient descent, the time complexity of the model is O(|Ω||S|2) to update the parameters in S for each iteration. For natural gradient the cost of inverting the Fisher information matrix G is O(|S|3), therefore the time complexity to update the parameters in S is O(|S|3 + |Ω||S|) for each iteration. The time complexity for natural gradient is significantly higher to invert the fisher information matrix, if the number of parameter is small, it is more efficient to use natural gradient because it requires significantly less iterations. However, if the number of parameters is large, it is more efficient to use gradient descent.

heading: 3 EXPERIMENTS
text: 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 i7-8700 6c/12t with 16GB of memory 1. In experiments of synthetic data, we simulate
1The code is available in the supplementary material and will be publicly available online after the peer review process.
random events using Equation (1). We generate an intensity function using a mixture of Gaussians, 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 is generated by directly drawing a 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) (Rosenblatt, 1956), an inhomogeneous Poisson process whose intensity is estimated by a reproducing kernel Hilbert space formulation (RKHS) (Flaxman et al., 2017), and a Dirichlet process mixture of Beta distributions (DP-beta) (Kottas, 2006). The hyper-parameters M and h in our proposed model are selected using grid search and cross-validation. For situations where a validation set is not available, then h could be selected using a 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.

heading: 3.1 EXPERIMENTS ON TWO-DIMENSIONAL PROCESSES
text: For our experiment, we use 20 Gaussian components and simulate a dense case with 100,000 observations and a sparse case with 1,000 observations within the time frame of 10 seconds. We consider that a joint event occurs if the two events occur 0.1 seconds apart. Figure 2a and Figure 2b compares the KL divergence between the first- and second-order models and Figure 3 are the corresponding intensity functions. In the first-order processes, both first- and second-order models have the same performance. This is expected as 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 projection of the Gaussian kernels on 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 two-dimensional 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.
Table 1 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 secondorder 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 a less sensitive bandwidth is more likely to learn a more accurate intensity function when the ground truth is unknown.

heading: 3.2 EXPERIMENTS ON HIGHER-DIMENSIONAL PROCESSES
text: We generate a fourth-order process to simulate the behaviour of the model in higher dimensions. The model is 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 105 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 no or 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 to the number of observations. The time complexity for KDE isO(ND) for the dimensionality withD, while DP-beta isO(N2K), whereK is the number of clusters, and RKHS isO(N2) for each iteration with respect to the sample size N , where DP-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 is more efficient because the time complexity is proportional to the number of bins in our model. The time complexity of APP for each iteration is O(|Ω||S|), where |Ω| = MD and |S| = ∑k c=1 ( D c ) . Our model scales combinatorially
Table 2: Negative test log-likelihood for the New York Taxi data. Single processes ([T] and [W]) and joint process of them ([T,W]). APP-# represents the order of the Additive Poisson Process.
Process APP-1 APP-2 KDE RKHS DP-beta
Jan [T] 714.07 714.07 713.77 728.13 731.01 [W] 745.60 745.60 745.23 853.42 790.04
[T,W] 249.60 246.05 380.22 259.29 260.30
Feb [T] 713.43 713.43 755.71 795.61 765.76 [W] 738.66 738.66 773.65 811.34 792.10
[T,W] 328.84 244.21 307.86 334.31 326.52
Mar [T] 716.72 716.72 733.74 755.48 741.28 [W] 738.06 738.06 816.99 853.33 832.43
[T,W] 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 k is often small, making it feasible to compute.
In Figure 4a we observe similar behaviour in the model, where the first-order processes fit precisely to the empirical distribution generated by the Gaussian kernels. The third-order model is able to period better on the fourth-order process. This is because the observation shown in Figure 5a is 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-order model is able to produce better estimation in thirdand fourth-order processes. The observations are extremely sparse as seen in Figure 5b, where there are only a few observations or no observations at all to learn the intensity function.

heading: 3.3 UNCOVERING COMMON PATTERNS IN THE NEW YORK TAXI DATASET
text: We demonstrate the capability of our model on the 2016 Green Taxi Trip dataset2. We are interested in finding the common pick up patterns across Tuesday and Wednesday. We define a common pick up time to be within 1 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 in 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-likelihood 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.

heading: 4 CONCLUSION
text: We have proposed a novel framework, called Additive Poisson Process (APP), to learn the intensity function of the higher-order interaction between point processes using samples from lower dimensional projections. We formulated our proposed model using the the log-linear 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 dimensional 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 a novel formulation to learn the joint intensity function which typically has extremely low intensity. There is enormous potential to apply APP to real-world applications, where higher order interaction effects need to be model such as in transportation, finance, ecology, and violent crimes.
2https://data.cityofnewyork.us/Transportation/2016-Green-Taxi-Trip-Data/hvrh-b6nb

heading: A ADDITIONAL EXPERIMENTS
text: A.1 BANDWIDTH SENSITIVITY ANALYSIS
Our first experiment is to demonstrate the ability for our proposed model to learn an intensity function from samples. We generate a Bernoulli process with probably of p = 0.1 to generate samples for every 1 seconds for 100 seconds to create a 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 applying no kernel, we learn the deltas of each individual observation. When we apply a 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, it 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 generate two experiments use the standard sinusoidal benchmark intensity function with a frequency of 20π. The dense experiment has troughs with 0 intensity and peaks at 201 and the sparse experiment has troughs with 0 intensity and peaks at 2. Figure 6d shows the experimental results of the dense case, our model has no problem learning the intensity function. We compare our results using KL 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.
Algorithm 2 Thinning Algorithm for non-homogenous Poisson Process 1: Function Thinning Algorithm (λ (t), T ): 2: n = m = 0, t0 = s0 = 0, λ̄ = sup0≤t≤Tλ (t) 3: repeat 4: u ∼ uniform (0, 1) 5: w = − 1
λ̄ lnu {w ∼ exponential(λ̄)}
6: sm+1 = sm + w 7: D ∼ uniform (0, 1) 8: if D ≤ λ(sm+1)
λ̄ then
9: tn+1 = sm+1 10: n = n+ 1 11: else 12: m = m+ 1 13: end if 14: if tn ≤ T then 15: return {tk}k=1,2,...,n 16: else 17: return {tk}k=1,2,...,n−1 18: end if 19: until sm ≤ T 20: End Function
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 6 10 − 5 10 − 4 10 − 3
KL Divergence
Pr oc
e : [ 1]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 6 10 − 5 10 − 4 10 − 3
KL Divergence
Pr oc
e : [ 2]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 6 10 − 5 10 − 4 10 − 3
KL Divergence
Pr oc
e : [ 3]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 6 10 − 5 10 − 4 10 − 3
KL Divergence
Pr oc
e : [ 4]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 4 10 − 3 10 − 2 KL Divergence
Pr oc
e : [ 1, 2 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 4 10 − 3 10 − 2 KL Divergence
Pr oc
e : [ 1, 3 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 4 10 − 3 10 − 2
KL Divergence
Pr oc
e : [ 1, 4 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 4 10 − 3 10 − 2
KL Divergence
Pr oc
e : [ 2, 3 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 4 10 − 3 10 − 2
KL Divergence
Pr oc
e : [ 2, 4 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 4 10 − 3 10 − 2 KL Divergence
Pr oc
e : [ 3, 4 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 3 10 − 2
KL Divergence
Pr oc
e : [ 1, 2 , 3
]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 3 10 − 2 10 − 1
KL Divergence
Pr oc
e : [ 1, 2 , 4
]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 3 10 − 2 10 − 1
KL Divergence
Pr oc
e : [ 1, 3 , 4
]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 3 10 − 2 10 − 1
KL Divergence
Pr oc
e : [ 2, 3 , 4
]
O rd
er : 1
O rd
er : 2
O rd
er : 3
O rd
er : 4
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 2 10 − 1
KL Divergence
Pr oc
e : [ 1, 2 , 3
, 4 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B an
dw id th
10 − 2 10 − 1 10 0
KL Divergence
To ta l K
L D iv er ge
nc e
(a )D
en se
ob se
rv at
io ns
.
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 4 10 − 3
KL Divergence
Pr oc
es s:
[ 1]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 4 10 − 3 KL Divergence
Pr oc
es s:
[ 2]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 4 10 − 3 KL Divergence
Pr oc
es s:
[ 3]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 4 10 − 3 KL Divergence
Pr oc
es s:
[ 4]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 3 10 − 2
KL Divergence
Pr oc
es s:
[ 1,
2 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 3 10 − 2
KL Divergence
Pr oc
es s:
[ 1,
3 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 3 10 − 2 KL Divergence
Pr oc
es s:
[ 1,
4 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 3 10 − 2 KL Divergence
Pr oc
es s:
[ 2,
3 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 3 10 − 2 KL Divergence
Pr oc
es s:
[ 2,
4 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 2
KL Divergence
Pr oc
es s:
[ 3,
4 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 2 10 − 1
KL Divergence
Pr oc
es s:
[ 1,
2 , 3
]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 1
KL Divergence
Pr oc
es s:
[ 1,
2 , 4
]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 2 10 − 1
KL Divergence
Pr oc
es s:
[ 1,
3 , 4
]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 1 KL Divergence
Pr oc
es s:
[ 2,
3 , 4
]
O rd
er : 1
O rd
er : 2
O rd
er : 3
O rd
er : 4
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 1 10 0 KL Divergence
Pr oc
es s:
[ 1,
2 , 3
, 4 ]
0. 00
0. 25
0. 50
0. 75
1. 00
1. 25
1. 50
1. 75
2. 00
Ke rn
el B
an dw
id h
10 − 1 10 0
KL Divergence
To a
l K L
D iv
er ge
nc e
(b )S
pa rs
e ob
se rv
at io
ns .
Fi gu
re 7:
K L
D iv
er ge
nc e
fo rf
ou r-
or de
rP oi
ss on
pr oc
es s.
0 2
4 6
8 10
Ti m e
0
10 00
00 Intensity Pr oc
es s: [ 1]
0 2
4 6
8 10
Ti m e
0
10 00
00
Intensity
Pr oc
es s: [ 2]
0 2
4 6
8 10
Ti m e
0
10 00
00
Intensity
Pr oc
es s: [ 3]
0 2
4 6
8 10
Ti m e
0
50 00
0
10 00
00
Intensity
Pr oc
es s: [ 4]
0 2
4 6
8 10
Ti m e
0 20 00 Intensity
Pr oc
es s: [ 1, 2 ]
0 2
4 6
8 10
Ti m e
0 20 00 40 00 Intensity
Pr oc
es s: [ 1, 3 ]
0 2
4 6
8 10
Ti m e
0 10 00 20 00 Intensity
Pr oc
es s: [ 1, 4 ]
0 2
4 6
8 10
Ti m e
0 10 00 20 00
Intensity
Pr oc
es s: [ 2, 3 ]
0 2
4 6
8 10
Ti m e
0 10 00 20 00 Intensity
Pr oc
es s: [ 2, 4 ]
0 2
4 6
8 10
Ti m e
0 10 00 20 00
Intensity
Pr oc
es s: [ 3, 4 ]
0 2
4 6
8 10
Ti m e
05010 0
Intensity
Pr oc
es s: [ 1, 2 , 3
]
0 2
4 6
8 10
Ti m e
02550 Intensity
Pr oc
es s: [ 1, 2 , 4
]
0 2
4 6
8 10
Ti m e
050 Intensity
Pr oc
es s: [ 1, 3 , 4
]
0 2
4 6
8 10
Ti m e
050 Intensity
Pr oc
es s: [ 2, 3 , 4
]
0 2
4 6
8 10
Ti m e
024 Intensity
Pr oc
es s: [ 1, 2 , 3
, 4 ]
O rd er : 1 O rd er : 2 O rd er : 3
O rd er : 4 G ro un
d Tr ut h
(a )D
en se
ob se
rv at
io ns
.
0 2
4 6
8 10
Ti m e
0 10 00 Intensity
Pr oc
es s: [ 1]
0 2
4 6
8 10
Ti m e
0 10 00 Intensity
Pr oc
es s: [ 2]
0 2
4 6
8 10
Ti m e
0 10 00 Intensity
Pr oc
es s: [ 3]
0 2
4 6
8 10
Ti m e
0 50 0 10 00 Intensity
Pr oc
es s: [ 4]
0 2
4 6
8 10
Ti m e
020 Intensity
Pr oc
es s: [ 1, 2 ]
0 2
4 6
8 10
Ti m e
02040 Intensity
Pr oc
es s: [ 1, 3 ]
0 2
4 6
8 10
Ti m e
01020 Intensity
Pr oc
es s: [ 1, 4 ]
0 2
4 6
8 10
Ti m e
01020 Intensity
Pr oc
es s: [ 2, 3 ]
0 2
4 6
8 10
Ti m e
01020 Intensity
Pr oc
es s: [ 2, 4 ]
0 2
4 6
8 10
Ti m e
01020 Intensity
Pr oc
es s: [ 3, 4 ]
0 2
4 6
8 10
Ti m e
0. 0 0. 5 Intensity Pr oc
es s: [ 1, 2 , 3
]
0 2
4 6
8 10
Ti m e
0. 0 0. 2 0. 4 Intensity
Pr oc
es s: [ 1, 2 , 4
]
0 2
4 6
8 10
Ti m e
0. 00 0. 25 0. 50 Intensity
Pr oc
es s: [ 1, 3 , 4
]
0 2
4 6
8 10
Ti m e
0. 0 0. 5 1. 0
Intensity
Pr oc
es s: [ 2, 3 , 4
]
0 2
4 6
8 10
Ti m e
0. 00 0. 05 Intensity
Pr oc
es s: [ 1, 2 , 3
, 4 ]
O rd er : 1 O rd er : 2 O rd er : 3
O rd er : 4 G ro un
d Tr ut h
(b )S
pa rs
e ob
se rv
at io
ns .
Fi gu
re 8:
In te
ns ity
fu nc
tio n
of hi
gh er
di m
en si
on al
pr oc
es se
s. D
ot s
re pr
es en
to bs
er va
tio ns
.

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