Additive Poisson Process: Learning Intensity of Higher-Order Interaction in Poisson Processes
Authors: Simon Luo
Feng Zhou
lamiae azizi
Mahito Sugiyama
Venue: ICLR 2022 Submitted (2022)
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 Poisson processes using projections into lower-dimensional space. 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 Poisson 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 2022
ADDITIVE POISSON PROCESS: LEARNING INTENSITY OF HIGHER-ORDER INTERACTION IN POISSON PROCESSES
Anonymous authors Paper under double-blind review
ABSTRACT
We present the Additive Poisson Process (APP), novel framework that can model the higher-order interaction effects of the intensity functions Poisson processes using projections into unwer-dimensional space. Our model combines the tech- niques from information geometry to model higher-order interactions on statis- tical manifold and in generalized additive models to use lower-dimensional pro- jections to overcome the effects from the curse of dimensionality. Our approach solves convex optimization problem by minimizing the KL divergence from a sample distribution in lower-dimensional projections to the distribution modeled by intensity function in the Poisson process. Our empirical results show that our model is able to use samples observed in the lower dimensional space to esti- mate the higher-order intensity function with extremely sparse observations.
INTRODUCTION
The Poisson process is counting process used in wide range of disciplines such as spatial- temporal sequential data in transportation (Zhou et al. 2018 (, finance (Ialan 2016) and ecol- ogy (Thompson 1955) to model the arrival rate by learning an intensity function. For given time interval, the integral of the intensity function represents the average number of events occurring in that interval.
The intensity function can be generalized to multiple dimensions. However, for most practical appli- cations, learning the multi-dimensional intensity function becomes a challenge because of the sparse observations. Despite the recent advances of Poisson processes, most current Poisson process mod- els are unable to learn the intensity function of a multi-dimensional Poisson process. Our research question is, "Are there any good ways of approximating the high dimensional intensity function?" Our proposed model, the Additive Poisson Process (APP), provides a novel solution to this problem.
Throughout this paper, will use running example in spatial-temporal setting. Say want to learn the intensity function for a taxi to pick up customers at a given time and location. For this setting, each event is multi-dimensional! that is, (I,, yi, W.), where pair of and yi represents two spatial coordinates and W; represents the day of the week. In addition, observation time is associated with each component of the event, which is represented as t, (tx., ty:; fw.). In this example, it is reasonable to assume tu, ty = fw, as they are observed simultaneously. For any given location or time, we can expect most few pick-up events, which makes difficult for any model to learn the low-valued intensity function. Figure 26 visualizes this problem. In this problem setup, if we would like to learn the intensity function at given location (æ., y) and day of the week W. the naïve approach would be to learn the intensity at (ar, y, W) directly from observations. This extremely difficult because there could be only few events for given location and day.
However, there is useful information in dower-dimensional space; for example, the marginalized ob- servations at the location (x,, y:) across all days of the week, or on the day W; at all locations. This information can be included into the model to improve the estimation of the joint intensity function. Using the information in diver-dimensionall space provides a structured approach to include prior information based on the location day of the week to improve the estimation of the joint intensity function. For example, a given location could be a shopping center or a hotel, where it is common for taxis to pick up passengers, and therefore expect more passengers at this location. There

---

Under review as a conference paper at ICLR 2022
could also be additional patterns that could be uncovered based on the day of the week. We can then use the observations of events to update our knowledge of the intensity function.
Previous approaches, such as kernel density estimation (KDE) (Rosenblatt 1956), able to learn the joint intensity function by using information in lower dimensions. However, KDE suffers from the curse of dimensionality, which means that it requires large size of samples build accurate model. In addition, the complexity of the model expands exponentially with respect to the number of dimensions, which makes it infeasible compute. Bayesian approaches, such as using mixture of beta distributions with Dirichlet prior (Kottas 2006 Kottas & Sansó (007)) and Reproducing Ker- nel Hilbert Space (RKHS) (Flaxman et | 77. have been proposed quantify the uncertainty for the intensity function. However, these approaches are often non-convex, making it difficult to obtain the globally optimal solution. Besides, if observations are sparse, it is hard for these approaches to learn reasonable intensity function. Additional related work about Bayesian inference for Poisson processes and Poisson factorization can be found Appendix A
In this paper, we propose a novel framework to learn the higher-order interaction effects of intensity functions in Poisson processes. Our model combines the techniques introduced by Luo & Sugiyama (2019) to model higher-order interactions between Poisson processes and by |Fiedmann & Stuetzle (198I) in generalized additive models to learn the joint intensity function using samples in a lower dimensional space. Our proposed approach is to decompose
Poisson process into lower-dimensional representations. For example, we have points (zz) in the z--imension and (w:) in the y-dimension. Such data in the lower-dimensional space can be used to improve the estimation of the joint intensity function. This is different from the traditional approaches where only the joint occurrence of events is used to learn the joint intensity.
We first show the connection between generalized additive models and Poisson processes, and then provide the connection between generalized additive models and the log-linear model (Agresti 2012), which has well-established theoretical background in information geometry (Amari [016). We draw parallels between the formulation of the generalized additive models and the binary log- linear model on partially ordered set (poset) (Sugiyama 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 Kulbback-Leibler (KL) divergence from the sample distri- bution in a power-dimensions space to the distribution modeled by the learned intensity function. This connection provides remarkable properties our model: the ability to learn higher-ordes inten- sity functions using lower-dimensional projections, thanks to the Kolmogorov-Aroolo representation theorem. This property makes advantageous to use our proposed approach for cases where there no observations, missing samples, or low event rates. Our model is flexible because it can cap- ture the interaction effects between events in a Poisson process as partial order structure in the log-ling linllinear model and the parameters of the model are fully customizable meet the requirements of the application. Our empirical results show that our model effectively uses samples projected onto lower dimensional space estimate the higher-order intensity function. More importantly, our model is also robust to various sample sizes.
2 FORMULATION
We start this section by introducing the technical background in the Poisson process and its extension to 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 tight link to the Poisson process modeled by GAMs. We show that the learning of APP can be achieved via convex optimization using natural gradient.
The Poisson process is characterized by an intensity function X Poisson process is an extension of homogeneous Poisson process, where the arrival rate changes with time. The process with time-changing intensity X(t) is defined as a counting process N(t), which has an independent increment property. For any time & > 0 and infinitesimal interval 8 > 0. the probability of events count is 6NN(t+
R - R. An inhomogeneous
= 0) = 1- (ooo00()()))))) ((((t+))-N(( = =)
8X(t)=00 - ((()) and ((N((t)) N(t) 2 2) o(6). where o(-) denotes little-o notation (Daley & Vere-
Jones 2007). This formulation can be generalized into the (inhomogeneous) multi-dimensional Poisson process, which can be used to learn the intensity function À RDD IR from realization
2

---

Under review as a conference paper at ICLR 2022
of timestamps t1. t22 ..... ty with t; [0, Each is the time of occurrence for the i-th event across D dimensions and T is the observation duration. The likelihood for the multi-dimensional Poisson process (Daley & Vere-Jones 2007) is given by
N (-p(-/XUUg A(US)) I
p ((( (A()) exp exp
(I)
where t = = (t
t (D)) E RR. We define the functional prior on A(t) as
X(t) 9 (f(t)) exp (f(t))
(2)
The function g(-) 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 the multi-dimensional Poisson process with the functional prior described as
log logr (ttr(tt) (A(E)) (A()) = =E ((tt)) -f (f (t)) dt.
In the following sections, we introduce the generalized additive models and propose to model by the log-ling log-linear model to learn f(t) and the normalizing term (ʃ exp (f (t)) dt).
2.1 GENERALIZED ADDITIVE MODEL.
In this section, we present the connection between Poisson processes and the Generalized Additive Model (GAM) proposed by Friedman & Stuetzle (198I). The GAM projects higher-dimennsinaal 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 com- pared to directly using smoothing in a nigher-dimensionss
space. For a given set processes J [D] (1....,D) D.. the traditional GAM using one-dimensional projections defined as (((o)) ((((()
log (t) = jje
B, with some smoothing function fj-
In this paper, extend it to include higher-order interactions between features in GAM by introduc- ing terms that represent the multiplicative product between events. The k-th order GAM is defined
log (t) (t)(
Eo (t(6)) ,..., ((0)) B3 fE
jEJ
=
Σ
f((t
B3,
(4)
where (4) |R!41 denotes the subvector jee) of with respect to C [D]. The function fr R!"l + IR is a smoothing function to fit the data, and the normalization constant 3 for the intensity function is obtained as B, = f A)(t))dt = exp( Eigh11 gk. (t() )dt. The definition of the additive model in the same form as Equation (3). particular, if we compare Equation (3) and Equation (4), we can see that the smoothing function f in Equation is realized as the summation over lower-dimensional projections in Equation (4).
Learning of continuous function using lower-dimensional projections is well known because of the Kolmogorov-Arnon representation theorem, which states that:
Theorem 1 (Kolmogorov-Anunn Arnold Representation Theorem (Braun & Griebel 2009 Kol-
mogorov 1957)). Any multivariate continuous function can be represented as superposition of one-dimensional functions; that is, f (f,,...,t,,) = Eant1 f, (E 1_p119m 9q. (t,))
Braun (2009) showed that the GAM is approximation to the general form presented in Kolmogorov- representation theorem by replacing the range
{l...., 2n + 1} with C J and the inner function Iu by the identity if q p and zero otherwise, yielding f(t) S(t)))
3

---

Under review as a conference paper at ICLR 2022
T=M
( s)
r=1
of of
(1.233)
(122
bettwan
== "j
Time
Example for k
X
Figure 1: Partial order structured sam- ple space (22 with D
(a)
(b)
Each
Figure 2: a A visualisation of the truncated parame-
node represents state and and the di- ter space to approximate the joint intensity function with rected edge represents the direction of the partial ordering.
4 and k = 2. (b) A visualization of the input datasets, where the blue points represent events with two spatial di- mensions and one time dimension.
Interestingly, the canonical form for additive models Equation 4 can be rearranged to in the same form as Kolmogorov-Arool representation theorem. By letting f(t) Eigs fi(t!))
g (A(t)) and g(-) = exp(-), we have
(t) exp (B,) exp Ei (ti)) exp (Eis fi (t()).
1
(5)
if [7] >> K for the k-th order model and 1/ exp(33) is the normal- ization term for the intensity function. Based on the Kolmogorov-Arool representation theorem, generalized additive models are able learn the intensity of the higher-order interaction between Poisson processes by using projections into lower dimensional spaces. The log-likelihood function for a kth-order model is obtained by substituting (t) in Equation (4p into A(t) in Equation (3)
where we assume ()((
log (( = Σ Σ (1)( -8',
i=1 1EI.
where 3' is constant given by 8' d(t) By. In the following subsection introduce log-lineat formulation equipped with partially ordered sample space, which aligns with the GAM formulation in Equation (5.
2.2 ADDITIVE POISSON PROCESS
We introduce our key technical contribution in this section. We introduce log-linear formulation called the additive Poisson process to estimate the parameters for the higher-order interactions in equation 6 We begin by discretizing the time window [0, T) into M bins and treat each bin natural number т E [M] (1,2...., M) for each process. The discretization avoids the need to compute the intractable integral in the likelihood function in Equation (6). The discretization approach tends to perform better in nigh-dimension compare to the alternative approaches such as variational inference. We assume that M is predetermined by the user. First we introduce structured space for the Poisson process to incorporate interactions between processes. Let Ω =
(J.T) J E 2|D| \ 0, E [M] C {(1,0)) We define the partial order < Davey & Priestley 2002
on Ω as
- J and r TST.
for each w, w' E S2
(7)
and (1, 3 for all w E Ωt which illustrated in Figure) 1 The relation J C J' is used to model anyerrder interactions between Poisson processes (Luo & Sugiyama, 2019) (Amari 2016 Section 6.8.4) and each t in (J, represents the auto-regressive component ("time") ") in our model. The symbol I. denotes the least element in the partial order structure, which is required normalize probabilities the log-linear model. Each node w in the partially ordered set (poset) represents the
state of the sample space and the arrows in Figure represeet the partial order relationship between two nodes?) that is, if + w', then w 3 w'
In information geometry, this corresponds hypergraphs which have simplicial complex structure Ay et al. 2018 Section 2.9).
"Thi graph structure should not be confused with the graph structure studied in graphical models, where the nodes typically represent random variable and the arrows represent the relationship hetween the two random variables.
4

---

Under review as a conference paper at ICLR 2022
Intuitively, the greatest node for each E
Algorithm 1 Additive Poísson Process (APP)
[M.] which is (11.2,3}, T) in Figure [ rep-
1: Function APP((t.) S, M, ht):
resents the multi-dimensional Poisson process.
2: Initialize Ω with the number M of bins
Other nodes represent projections onto lower-
3: Apply Gaussian Kernel with bandwidth
dimensional space that correspond to the marginal-
on [t:] to compute
ized observations: for example, {{I}, (2). { and {[1,2), {1, 3), 333 represent the first- and second-order processes. Using our example in In- troduction, where we wanted to estimate the in- tensity function of pick-up event of taxi, {1} and (2) correspond to spacial coordinates æ and y. respectively, and {3) to the day of the week W. and r represents the (discretized) observation time. We can then update our belief to model the second-order intensity function using obser- vations of the second order events. For example,
4: Compute n) (ñs))es from ÿ
5: Initialize 0 = (0,)sEs (randomly or 0g
0)
6: repeat
7: Compute P using the current (0s)sES
8:
Compute = (Ms)ss from p
9:
An t 1-1
10:
Compute the Fisher information matrix
G using Equation (II)
{1, 2), {1, 3), {2, 3) represents an event occurring at (I, {ær, W). and (y, W). We can then con- 12:
11:
0 0 a An
until convergence of 0 (f.)ses
tinue this process to an arbitrary order of interac- 13: End Function
tions. Later on in this section we introduce the mathematies to estimate the higher-ordee function using a restricted number of nower-dimensiondl projections.
The domain of can be generalized from [M]] to [M]] 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 demonstrate that we can still accurately estimate the density of processes. Our model can be applied not only to time-series data, but to any sequential data.
On any set equipped with partial order, we can introduce a log-linear model (Sugiyama et al., 2017). Let us assume that a parameter domain S _ St is given. For a partially ordered set (S2, 3), the log- log-linear model with parameters (0*)sES is introduced as
log logp(w;0)
Σ SES 1s - (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 Ewea exp( Eees 1 (s = - (.0)- A special case of this formulation coincides with the density function of the Boltamann machines (Sugiyamal
et al. 2018 Luo & Sugiyama 2019).
Here there is clear correspondence between the log-linear formulation and that in the form of Kolmogorov- representation theorem in Equation (5) we rewrite Equation X as
p(u; 0) cxp (0) exp (Eees (ssww)..) x exp < 1 1s
We call this model with (S2, =) defined in Equation 7 the additive Poisson process, which rep- resents the intensity À as the joint distribution across all possible states. The intensity of the multi- Poisson process given via the GAM in Equation (5) is fully modeled (parameter- ized) by Equation (8) and each intensity fr(-) is obtained as 0g..)- To consider the k-th order model, we consistently use the parameter domain S, given as S (J. e Ω [J] < & }.. where is an input parameter the model that specifies the upper bound of the order of interactions. This means that 0, = 0 for all 8 & S. Note that our model is well-defined for any subset S C Ω and the user can arbitrary domain in applications. A visualization of the truncated parameter space is shown in Figure 2a
For given J C [D] and each bin 7 with w (J,T), the empirical probability p(w) of input observations is given as
N
(() - t
1 Σ (T(())
p(w) = Z
= Σ p(w),
0T(T) = Nhj
K
(10)
hr
ICJ
wERt
for each discretized state w = (J, where T = (T,...,77) T) E R.. The function Ø, performs smoothing time stamps t1,-.., tN, which is the kernel smoother proposed by Buja al. 1989).
5

---

Under review as a conference paper at ICLR 2022
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 SSchbee (993)..
2.3
OPTIMIZATION
Given empirical distribution P defined in Equation 10 the task is to learn the parameter (0g sesE such that the distribution via the log-linear model in Equation (8) is as close to p as much as possible. Let us define Gş = {p 0, 0 if S); which is the of distributions that be represented by the log-linear model using the parameter domain S. Then the objective function is given as minpees DKL (p. p). where (P. p) = Ewen P log(p/p) is the KL divergence from á p. In this optimization, let p" be the learned distribution from the sample with infinitely large sample size and let p be the learned distribution for each sample. Then we can lower bound the uncertainty (variance) (p' (p',ppp))) p)] by |8|/2N (Barron & Hengartner 1998).
Thanks the well-developed theory of information geometry (Amani 2016) for the log-lineas model (Amari] 2001), it is known that this problem can be solved by e-projection, which coincides with the maximum likelihood estimation and is always convex optimization (Amaril 2016] Chapter 2.8.3). The gradient with respect each parameter 0g is obtained by (0/00,) (p.2)p (p.p) = 7/+ fa.
where T/ = Ewen |wzzs)((w). The value 17s is known as the expectation parameter (Sugiyama et al. 2017) and 1is is obtained by replacing with á in the above equation. If ńr 0 for some 8 E S, we remove from S to ensure that the model is well-defined.
Let S {s....., } (|s)) and 0 [0........ 0g 1",
Meis]]" IT. We always 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 tnext A G G- f)), where the
[M . ,....,
Fisher information matrix G is obtained as
Pij = 00,,00,, DKL(F.p) =
Σ - WEB
no. Nog-
(11)
Theoretically the Fisher information matrix is numerically stable perform matrix inversion. However, computationally. floating point errors may cause the matrix to become indefinite. To over- come 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 computationally positive semi-definite NNeal 1999).
The pseudocode for APP is shown in Algorithm The time complexity of computing line 7 is ((9[[S)). This means when implementing the model using gradient descent, the time complexity of the model is (1 to update the parameters in S for each iteration. For natural gradient, the cost of inverting the Fisher information matrix G is ((S))): therefore, the time complexity to update the parameters in S is ( + [S][[]] for each iteration. The time complexity for natural gradient is significantly higher higher because of the requirement to invert the fisher information matrix; if the number of parameters is small, it is more efficient to use natural gradient because it requires significantly fewer iterations. However, if the number of parameters is large, is more efficient to use gradient descent.
EXPERIMENTS
We perform experiments using two-dimensional synthetic data, higher-dimensional synthetic data, and real-world data to evaluate the performance of our proposed approach. Our code is imple- mented in 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 [ In experiments with synthetic data, we simulate random events using Equation ( We generate an intensity function using a mixture of Gaussians, where the mean is drawn from 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 sample from the probability density function. An arbitrary number of samples is drawn from the mixture of Gaussians. We
"The code is available in the supplementary material and will be publicly available online.

---

Under review as a conference paper at ICLR 2022
Presse
Presest Dt
Dense observations N 100, 000.
(b) Sparse observations, N = 1000.
Figure 3: KL Divergence for second-order Poisson process. The order of the model (color of the line) represents the k-th order model, i.e., l: (blue) and k = 2 (orange).
(a) Dense observations, 1000, h 0.4.
(b) Sparse observations, N 100, 000, à 0.2.
Figure 4: Intensity function of two dimensional processes. Dots represent observations. Left: Repre- sents marginalized observation of the first dimension. Middle: Represents marginalized observation of the second dimension. Right: The joint observation of dimension and 2. The order of the model (color of the line) represents the k-th order model, i.e., k = 1 (blue) and k = 2 (orange).
The lowest KL divergence from the ground truth distribution to the obtained distri- bution on two types of single processes ([1] and [2]) and joint process of them ([1,2]). APP- # represents the order of the Additive Poisson Process. Missing values mean that the compu- tation did not finish within two days.
Table 1:
Table 2: Negative test log-likeli 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.
APP-1
RKHS
[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
RKHS
DP-beta
[T.]]
249.60
380.22
259.22 299.29
4.98c3
4.98c55
2.8e.4
[T]
713.43
713.43
755.71
795.61
765.76
Dense
[2]
2836-5
283e-5
1.17c4
Feb
[W]
738.66
738.66
773.65
811.34
792.10
[1.2]
2.98--2
1.27c-3
6,33e-4
4.09e-2
4.51e2
[T.]
328.84
244.21
307.
334.31
326.52
7.26e-4
7,26e-4
8.83e-1
1.96e22
2,62c-3
[T]
716.72
716.72
733.74
755
741.28
Sparse
2.28c-4
2.76e-4
2.35e-3
2.49e-3
Mar
[W]
738.06
738.06
$16:99
853.33
832.43
[1.2]
2.886-2
1.77e-2
3.676-3
1.84e-2
3.68e-2
291. 20
246.19
289.60
328.47
300.36
then run our models and compare with Kernel Density Estimation (KDE) (Rosenblati 1956), an
inhomogeneous Poisson process whose intensity is estimated by reproducing kernel Hilbert space formulation (RKHS) (Flaxman et al.] | (007)) and a Dirichlet process mixture of Beta distributions (DP-beta) (Kottas 2006 Kottas & Sansó 2007). The hyper-parameters M and f in proposed model are selected using grid search and vross-validation. For situations where a validation set is not available, then h could be selected using rule of thumb approach such as Scott's Rule (Scott 2013) and M could be selected empirically from the input data by computing the time interval of the joint observation.
3.1 EXPERIMENTS ON D--NNNINSIONL PROCESSES
For our experiment, use 20 Gaussian components and simulate dense case with 100,000 obser- vations and a sparse case with 1,000 observations within the time frame 10 seconds. We consider that a joint event occurs if the two events oceur 0.1 seconds apart. Figure Ba and Figure 3b com- pare the KL divergence between the first- and second-order models and plots in Figure are the corresponding intensity functions. In the first-order processes, both first- and second-order models
7

---

Under review as a conference paper at ICLR 2022
Prasem
1..3
1-2.311
Fromuy D.2.0.0
Teld Sharpis
(a) Dense observations,
N
10"
D1
(b) Sparse observations, N 10
Figure 5: KL Divergence for fourth-order Poisson process. We selected four representative examples for experimental results, full results available in the supplementary material. The line color signifies the order of the model, i.e., k = (blue), k = 2 (orange), 3 (green) and = 4 (red).
(a) Dense observations, N 10"
Proces EDI
Poces 13
Poome ts 1n
Prices
(b) Sparse observations, 10
Figure 6: Intensity function of higher dimensional processes. Dots represent observations. We have selected four representative examples for our experimental results, full results available in the supplementary material. The order of the model (color of the line) represents the k-th order model, i.e., k = 1 (blue), &- (orange), k (green) and k (red).
have the same performance. This is expected, as both of the models can treat first-order interactions and are 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 wwo-dimensional space. In contrast, the first-order model must approximate the second-order process using the observations from the first-order pro- cesses. In the sparse case, the second-order model performs better when the correct bandwidth is selected.
Table compares our approach APP with other satee-of-thearr 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 two days for the dense experiment. In the second- order process, our approach was outperformed by KDE, while the second-order APP was able to outperform both RKHS and DP-beta process for both sparse and dense experiments. Figures 3a and 3bb 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-DIMENSIONN PROCESSES
We generate a fourth-order process to simulate the behavior of the model in higher dimensions. The model is generalizable to higher dimensions, but 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 sample with the size of 10" for the dense case and that with the size of 10 for the sparse case. We consider the joint event to be the time frame of 0.1 seconds.
We not able when there are no or few joint observations in third- and fourth-order processes. In addition, the time complexity is too high learn from joint observations in first- and second-order processes because
comparison experiments with other models because they unable to learn
8

---

Under review as a conference paper at ICLR 2022
all the other models have their time complexity proportional to the number of observations. The time complexity for KDE is O(N') for the dimensionality with D, while DP-beta is O(N²K), where K is the number of clusters, and RKHS is ( for each iteration with respect the sample size N, where DP-beta and RKHS are applied directly on the joint observation as they cannot use the projections in lower-dimensional space. KDE is able to make an estimation of the intensity function using projections in lower-dimensional space, but it was too computationally expensive to complete running the experiment. By contrast, 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 (( where [2] MP and |S] = Σ -و= (2). Our model scales combinatorially with respect to the number of dimensions. However, this is unavoidable for any model that 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 Sal we observe similar behavior 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 predict better on the fourth-order process. This is because the observation shown in Figure 6aa is largely sparse and learning from the observations directly may overfit. A lower-dimensional approximation is able provide a better result in the third-order model. Similar trends can be seen in the sparse case, as shown in Figure 5b where a second-order model is able to produce better estimation in third- and fourth-order processes. The observations are extremely sparse, as seen in Figure [b] 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 which is a open source dataset with a CCO: Public Domain licences. We are interested in finding the common pick-up patterns across Tuesdays and Wednesdays. 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 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 Tuesdays and Wednesdays in January 2016 as our training and validation sets and the 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- lig-likel 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, 1PP-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.
4 CONCLUSION
We have proposed novel framework, called Additive Poisson Process (APP), to learn the intensity function of the higher-order interaction between Poisson processes using observations projected into lower-dimensional spaces. We formulated our proposed model using the log-linear model and opti- mize it using information geometric structure of the distribution space. We drew parallels between our proposed model and generalized additive model and showed the ability learn from lower dimensional projections via the Kolmogorov-Aruol representation theorem. Our empirical results show the superiority of our method when learning the higher-order interactions between Poisson processes and when there are no or extremely sparse joint observations. 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 is enormous potential to apply APP to real-world applications, where higher order interaction effects need to be modeled such as transportation, fi- nance and ecology.

---

Under review as a conference paper at ICLR 2022
ETHICS STATEMENT
This paper is intended to be written neutral tone. Any suggestion of any application that may be considered unethical have negative social impact is only coincidental and not our intention. The results of this work should only be used to have a positive impact in society.
REPRODUCIBILITY STATEMENT
Our model is implemented in Python 3.7.5 with NumPy version 1.8.2. The code to reproduce experiments available in the supplementary material and will be made publicly available GitHub after the paper is accepted.
REFERENCES
Alan Agresti. Categorical Data Analysis. Wiley, 3 edition, 2012.
Shun-Ichi Amari. Natural gradient works efficiently learning. Neural Computation, 10(2):251- 276, 1998.
Shun-Ichi Amari. Information geometry on hierarchy of probability distributions. IEEE Transac- fions on Information Theory, 47(5) 2001.
Sun-Ichi Amari. Information Geometry and Its Applications. Springer, 2016.
Nihat Ay, Paolo Gibilisco, and F Matus. Information Geometry and Its Applications. Springer, 2018.
A. Barron and N. Hengartner. Information theory and superefficiency. The Annals of Statistics, 26 (0:1800-1825 1998.
Jürgen Braun. An application of Kolmogorov's superposition theorem to function reconstruction in higher dimensions. PhD thesis, Universitäts-und Landesbibliothek Bonn, 2009.
Jürgen Braun and Michael Griebel. On a constructive proof of Kolmogorov's '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.
Eric C Chi and Tamara G Kolda. On tensors, sparsity, and nonnegative factorizations. SIAM Journal on Matrix Analysis and Applications, 33(4): 1221119299999 2012.
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 kernels. Electronic Journal of Statistics, 11(2)5081-51044 2017.
Jerome H Friedman and Werner Stuetzle. Projection pursuit regression. Journal of the American Statistical Association, 166376) 1981.
Deniz Ilalan. A poisson process with random intensity for modeling financial stability. The Spanish Review of Financial Economics, 14(2):43-00, 2016.
Andrei Nikolaevich Kolmogorov. On the representation of continuous functions of many variables by superposition of continuous functions of one variable and addition. Doklady Akademii Nauk, 114(5):953-99 95569 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 Machine Learning (ICML), volume 47, 2006.
10

---

Under review as a conference paper at ICLR 2022
Athanasios Kottas and Bruno Sansó. Bayesian mixture modeling for spatial poisson process inten- sities, with applications to extreme value analysis. Journal of Statistical Planning and Inference, 137(10):3151-3166 2007.
Chris Lloyd, Tom Gunter, Michael Osborne, and Stephen Roberts. Variational inference for gaussian process modulated poisson processes. In International Conference on Machine Learning, pp- 1814-1822. PMLR, 2015.
Simon Luo and Mahito Sugiyama. ing higher-order feature interactions. Proceedings of the 33rd AAAI Conference on Artificial Intelligence, pp. 4488 0188-449,444 2019.
trade-off in hierarchical probabilistic models
Radford M Neal. Regression and classification using gaussian process priors. Bayesian Statistics, 6:475-501, 1999.
Krzysztof Nowicki and Tom A B Snijders. Estimation and prediction for stochastic blockstructures. Journal of the American statistical association, 96(455):1077-1088, 2001.
Yosihiko Ogata. On Lewis' simulation method for point processes. IEEE Transactions Informa- tion Theory, 27(1):23-31, 1981.
Murray Rosenblatt, Remarks on some nonparametric estimates of density function. The Annals of Mathematical Statistics, pp. 832-837, 1956.
Yves-Laurent Kom Samo and Stephen Roberts. Scalable nonparametric bayesian inference point processes with gaussian processes. In International Conference on Machine Learning, pp. 2227- 2236. PMLR, 2015.
H Schäbe. Nonparametric estimation of intensities of nonhomogeneous poisson processes. Statisti- cal Pupers, 34(1):113-131, 1993.
Aaron Schein, Mingyuan Zhou, David Blei, and Hanna Wallach. Bayesian poisson tucker decompo- sition for learning the structure of international relations. In International Conference on Machine Learning, pp. 2810-2819. PMLR, 2016.
David Scott. Multivariate density estimation: theory, practice, and visualization. John Wiley & Sons, 2015.
Mahito Sugiyama, Hiroyuki Nakahara, and Koji Tsuda. Tensor balancing on statistical manifold. In Proceedings of the 34th International Conference on 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-8835, 2018.
HR Thompson. Spatial point processes, with applications ecology. Biometrika, 42(1/2)):02-115,5 1955.
Feng Zhou, Zhidong Li, Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen. A refined MISD algorithm based on Gaussian process regression. In Pacific-Asia Conference on Knowledge Dis- covery and Data Mining, pp. 584 596. Springer, 2018.
Feng Zhou, Simon Luo, Zhidong Li. Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen. Efficient em on-variational inference for nonparametrie hawkes process. Statistics and Computing, 1114):1-11,021
11

---

Under review as a conference paper at ICLR 2022
APPENDIX
A RELATED WORK
A.1 BAYESIAN INFERENCE FOR POISSON PROCESS
Learning an intensity function from sparse high-dimensional datasets is challenging problem. It often comes a trade-off between computational complexity and numerical accuracy. For applications which are sparse with higher dimensionality, numerical accuracy is often sacrificed to approximate the intensity function.
The na've approach using MCMC has cubic complexity O(n)) ) for each dimension and increases exponentially with respect with additional input dimension, that is, (())))), where n refers to the number of observations and d is the number of dimensions. One example of this approach is using a Dirichlet process mixture of Beta distribution as prior for the intensity function of the Poisson process (Kottas) 2006;; Kottas & Sansó 2007). The solution from MCMC is often accurate and asymtopically converges to the true posterior intensity function. However, due to its computational complexity, it is infeasible estimate any high-dimensional intensity function.
More recent techniques have attempted to scale up these approaches by using Gaussian processes as functional prior the intensity function (Samo & Roberts (2015). However, the model complexity is (n²kk k) for each dimension and is exponential with respect to the number of input dimensions, which is still infeasible estimate any high-dimensiioal intensity function.
Variational inference (Lloyd et al. 2015) approaches that can be used to make Bayesian inference of Poisson process much more efficient scaling up the linear complexity ((n) for each dimension. However variational inference is not guaranteed to asymptotically converge to the true posterior distribution.
Our approach uses the discretization approach to scale up the model to higher dimensions. We use a graph (partial order) structure allow the flexibility for domain expertise to specific how each dimension is treated and which interaction effects should be included into the model.
A.2 POISSON FACTORIZATION
Our work is closely related to Poisson Factorization (Chi & Kolda 2012). where random variables in a tensor are represented with a Poisson distributed or Poisson process likelihood. The tensor is usually used represent some high-dimensional dataset such contingency tables or other collec- tion of counting datasets, which are often large and sparse. The objective of Poison Factorization decompose the high-ddmensiioal sparse matrices into lower dimensional space, where we can find some meaningful latent structure.
The effectiveness of Poisson factorization for high dimensional datasets makes it ideal analyze problems consisting of sparse count data. One example of this work is Bayesian Poisson Tucker decomposition (BPTD) (Schein et al. 2016). where dataset interaction events is represented as a set of N events, each of which consists of a pair of token e; that encodes certain features and time, that is, (e;, t.)) BPTD uses an MCMC inference algorithm to learn the latent structure, which is based 2001) with a Poisson likelihood.
extension of stochastic block models (SBM) (Nowicki & Snijders)
Our approach provides generalization of this idea of Poisson Factorization by using Legendre tensor decomposition (Sugiyama et al. 2018) and demonstrating its ability on a teatial-temporal problem. Our optimization is much more efficient it is guided by gradients to minimize the KL- divergence. Our approach also contains a graph structure which allows domain experts to encode certain properties into the model.
12

---

Under review as a conference paper at ICLR 2022
Algorithm 2 Thinning Algorithm for non-homogenous Poisson Process
1: Function Thinning Algorithm (A (t), T):
2: n m 0, lo = 80 0, suposigTA (1)
3: repeat
4: #f ~ uniform (0, 1)
5:
w w=
{w exponential {X)}
6:
8m+1
+
7:
D ~ uniform (0,
D < M(mm1 then
8:
9:
bn+1 5m+1
10:
=n++
else
12: m m +
13: end if
14: if In < T then
15:
return {tx k=,,...
else
17: return {tkJR= 122
18:
end if
19:
until 8m ST
20:
End Function
B ADDITIONAL EXPERIMENTS
B. BANDWIDTH SENSITIVITY ANALYSIS
demonstrate the ability for our proposed model learn an intensity func- tion from samples. We generate Bernoulii process with probably of 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 model. In Figure 7al we observe that apply- ing no kernel, 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 seconds and 80 seconds mark, can be seen when observations have overlapping kernels, the intensity function becomes larger in magnitude.
Our first experiment
B.2 ONE DIMENSIONAL POISSON PROCESS
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 0 intensity and peaks at 2. Figure 7d 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 7b shows that the optimal bandwidth h =
13

---

Under review as a conference paper at ICLR 2022
No Karnel
0.125
0.3
0.8
Divergence
0.4
0.075
Karnel,
0.050
0.025
0.0
0.000
-0.3
20
40
60
8.0
100
Time
(a) Toy Example
(b) KL divergence dense experi- (c) divergence sparse experi- ment
ment
h2
(d) Ogata's thinning algorithm with high intensity
(c) Ogata's 's thinning algorithm with low intensity
Figure 7: One dimensional experiments
14

---

Under review as a conference paper at ICLR 2022
OO An TA LAT
25 OF OO DO LEE TN L. BRNEEE
Prosee I. 2.4)
LI DO ES LEE Baseeit
KImH Samml Berteeth
M
Froomet
Prooees: NI
OO ES LE
ES
Figure 8: KL Divergence for four-order Poisson process.
ZZT SO 3xE
Lo 2n Zie
Prosse 1..)
Ler baroedsa
LET TO Sn basdedd dwidd
Kern Berenits
1.30
PRICEE
Prices
Dense observations.
(b) Sparse observations.
aoaddaeg X
S00dsm00 TX
assndme0
#XUndamAD TM
12 TM
sooddeeo TM
12 OF On TH 19 311 31
OF TY L/ SR
OE AT OF TE 1I
Preae 41
Komel Bunde inddwii
Karal Hindwets
4)
2,344
Prices::
D0 OE Do TH 21 LIT
OT TO DE LTE
50 OI DE DO TOR EET L12 BEAACC
Prea (3)
Pro 20
Kuraa Bantaith
Drssnn 22.4 41
Preeeass (3)
LA
Rarial aentiich 21
Koral Setticc
43
Brich
21
15

---

Under review as a conference paper at ICLR 2022
Ground Trush
Order:
Process:
Time
Process:
Time
Tima
Order:
Order:
Order:
Process:
100000
50000
2000
1008
22
Process:
Time
Process:
Time
Time
Time
Process:
(a) Dense observations.
1000
1
10
Process:
Time
Process:
Time
Process:
Time
Time
Process::
100000
2000
1090
8
Process:
Time
Process:
Time
Process:
Time
Process:
Time
1000
1
16

---

Under review as a conference paper at ICLR 2022
10
10
Onler: Ground
Process: 12 31
Process: 2.
Process: (4)
Onter: Onter: Order: 3
Time
Time
Time
Figure 10: Intensity function of higher dimensional processes. Dots represent observations.
16
18
16
10
Proces: 2. 4,
Process: Dt 4]
Process: 2.
Process: 33
Time
Time
Time
Time
(a) Sparse observations.
8.06
intaniity
10
16
Process: 3.3.4 4)
Proces:: f1 31
Process: (3 4]
Process: [2]
Time
Time
Time
Time
1000
0
10
10
10
Process: [1, 21
Process: f1, 4)
Process: [1]
Process: (2.
Time
Time
Time
Time
1000
17


Reviews:
[{"review_id": "f5WK0PQt-R", "reviewer": "Reviewer_s5rK", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "I found the main idea of the paper novel and interesting. I believe the paper is well written and easy to read even for someone not familiar with information geometry. I particularly appreciated the theoretical guarantees that come with this model and inference scheme. The experimental section is extensive but could be improved (see comments below). Addressing the following comments could help clarify the benefits and contributions of the proposed approach:\n\n1) It seems to me that the point of this paper is not to model multiple correlated poisson processes but rather to learn the intensity of a single poisson process by “decomposing” it into sub processes over lower dimensional input spaces. Can you confirm this is correct?\n2) The authors speak about general high dimensional input spaces in the introduction also making the example of space-time locations. However, the input definition is not clear in Section 2. From Section 2, the authors speak only about timestamps where each time step seems to have a D dimensionality. This is confusing. It would be useful to repeat your initial example when you define $\\mathbf{t}_1, …. \\mathbf{t}_N$ to give examples of what these time stamps are. Ultimately it seems to me that you could just call $\\mathbf{x}_i = [x_i^1, …., x_i^{D-1}, t_i ] $ the input vector in a D dimensional space which also includes the time stamp. \n3)Section 2.3. In terms of computational complexity, what is generally the cardinality of $\\mathcal{S}$? Is this not given by the number of terms considered in Eq (6)? In other words the number of lower dimensional projections in Eq 4? If this is the case it does not seem the computational complexity is significantly lower compared to the methods mentioned in Section A.1 in the appendix. \n4) I believe the authors should comment an important approximation introduced when discretising the input space thus getting rid of the intractable integral in the likelihood? There are indeed approaches in the literature that avoid that and despite resorting to variational inference (thus no guarantees) are scalable and seem to be working fine in experimental settings (see [1] and [2] below).\n5) Figure 4. It seems like there is a typo in the captions. For the sparse case N should be 1000? Also it does not seem like the observations are sparse in the two left plots of the second row.   \n6) Experimental comparison. Given the number of approaches existing in the literature and using GP to model the intensity function (and same link function used in Eq (2)) I think the authors should compare the proposed method to one of them (one example could be the LGCP model) especially cause for some of them the approximated inference scheme is fast and stochastic gradient descent method can be used to reduce the computational complexity. In addition there exist out of the box tools that can be used to train such models making it a fast default choice. They would probably fail in sparse settings but I think the paper would significantly benefit from that comparison. \n\n\n**Minor comments**\n1) Just above Eq (6) I believe the authors refer to Eq (3) and not Eq (1)\n2) The figures are very scatters in the paper and I would prefer to have them reorganised in order to have them closer to the point where they are discussed in the text. In addition it would be useful to see Fig 4 together with fig 3 and fig 5 together with fig 6 as they refer to the same experiment. \n\n[1] Aglietti, Virginia, et al. \"Structured variational inference in continuous cox process models.\" arXiv preprint arXiv:1906.03161 (2019).\n\n[2] Lloyd, Chris, et al. \"Latent point process allocation.\" Artificial Intelligence and Statistics. PMLR, 2016.\n", "summary_of_the_paper": "This paper proposes a new model, called additive Poisson process, for estimating the intensity of a high dimensional Poisson process. By representing the intensity as a sum of lower dimensional projections, the method can deal with sparse data while allowing for an efficient inference scheme based on the information geometric structure of the distribution space. The assumed intensity form is motivated by the Kolmogorov-Arnold theorem while the estimation procedure has some formal guarantees of convergence.", "summary_of_the_review": "I found the main idea of the paper novel and interesting. I believe the paper is well written and I particularly appreciated the theoretical guarantees that come with this model and inference scheme. However there are some parts of the paper that are confusing, even in the basic formulation section. Therefore, I believe the paper is borderline and could be significantly improved by clarifying the points above. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "c6_YoXx6uW-", "reviewer": "Reviewer_cFZz", "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": "The paper is overall well written and I read it with interest. The major issue I have with this paper is that the proposed approach has long existed in the statistics literature for decades, under the context of generalized additive models. In particular, in chapter 5.6 of  [1], the tensor product bases are more flexible than the approach proposed in this paper. In my opinion, the lower-dimensional approximation proposed in this paper is essentially a piecewise constant approximation of the intensity function (in lower dimension), which is a special case of the tensor product basses.  If comparisons should be made, it should be between the proposed method and the GAM with tensor product bases, instead of the RKHS formulation. \n\nIn addition, the partial order in the formulation (7) seems to rely on the assumption that all dimensions have the same common support [0, T], which is not the case for the spatial-temporal processes considered in this paper. How is this issue addressed?\n\nFor the proposed method, the most difficult task should be to determine the value of hyper-parameters M and h. If M is too large, it can easily lead to overfitting. If M is too small, there will be a large bias in the estimated intensity functions. Please provide more details on how cross-validation is performed to choose an optimal combination of h and M.\n\nAnother difficult task with the proposed method is how to justify the use of lower-dimensional approximation. For example, is it justifiable to assume that there is no need to consider the spatial and temporal interaction in a spatial-temporal point process? In GAM, one can do this by performing an analysis of variance (ANOVA), see, e.g., section 5.6.3 of [1]. Could this be done with the proposed approach?\n\nMinor points: \n1. In the motivating example, the time coordinate is the day of the week, which is discrete. Such an example is not coherent with the definition of the Poisson process, which should be with continuous coordinates.\n2.  in equation (6), should it be $t_i^{(I)}$ instead of $t^{(I)}$ inside the function f_{I}?\n\nReference:\n\n[1]. Wood, S.N., 2017. Generalized additive models: an introduction with R. CRC Press. (2nd Edition)", "summary_of_the_paper": "This paper proposes a method to estimate the multi-dimensional intensity function of a Poisson process with a lower-dimensional projection. The proposed method is motivated by perspectives of information geometry and generalized additive models in statistics. The performance of the proposed method is extreme through synthetic examples and a real data analysis.", "summary_of_the_review": "I think the paper is well written, but I do not see the methodology proposed to be better than the existing ones. ", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "CkG-ffERTN", "reviewer": "Reviewer_wbmT", "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": "This approach is interesting but it would deserve more works to further develop the idea. The paper is not motivated and clear enough which leaves the reader perplexed. Here are some remarks: \n\n\n- Several statements in the paper are not well-justified by references, e.g. on KDE in the introduction.\n\n- More details should be added on how GAM approximates the Kolmogorov-Arnold representation and how well it performs. For example, 'based on the K-A representation theorem, GAM are able to learn the intensity of the higher-order interaction between Poisson processes' It would be nice to add references and examples since it is the key of the paper.\n\n- Motivations for introducing APP are not well explained. Background on Poisson Processes are not clear enough and do not give any insight about PP.\n\n- The presentation of APP, in Section 2.2, should be clarified. It is not clear how to match the log-linear model (8) and the intensity defined in Sections 2 and 2.1, e.g. the taxi toy example at the top of page 5 makes things even more confusing. \n\n- I would like to see experiments on the computation time of the approach compared to existing methods as it could be an advantage of this approach.\n\nMinor comments:\n\n- The notation \\mathbf{t}_i in Section 2 is very confusing. \n\n- Figure 2 is cited in the text before Figure 1.\n\n- Section 2.2 The symbol $\\bot$ should be defined before, or just after,  using it. Furthermore,  $\\bot$ denotes the least element in the partial order structure but is used before 'partial order' definition. Also, I think $J$ belongs to $ 2^{|[D]|}$.\n\n- Equation (5) is weird since you basically inverse log-function on equation (4) but the index $k$ disappears and appears again in (6). \n\n- Figures are too small, it is difficult to see anything.\n\n- \"to arrive at\": \"to converge to\"", "summary_of_the_paper": "The paper propose to apply generalized additive models (GAM) to learn the intensity of the multi-dimensional Poisson processes (PP). ", "summary_of_the_review": "The paper is not motivated and clear enough. Furthermore, the proposed approach is empirically poorly studied while no theoretical guarantees are investigated. Several claims are not supported by references or proof.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "3zryyCAccY", "reviewer": "Reviewer_AyVi", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "•\tThe idea of restricting your model to lower-order interactions is well-motivated and flexible. \n\n•\tThe paper is well written (I only found one typo!). \n\n•\tIn synthetic data experiments, why not compare \\lambda with the estimators directly? \n\n•\tIn Figure 2 (a), why does it read \\eta = \\hat{\\eta}? I thought the purpose of the figure was to show that the parameters for order-3 or -4 interactions were set to zero?\n\n•\tCould you please give a real-life example of a four-dimensional (or higher) Poisson process? It seems to me like most real-life examples are at most three-dimensional.\n\n•\tYou said that KDE is too slow for d >= 4. How does it do for d = 3?\n\n•\tHow was the bandwidth chosen in the taxi experiment?\n\n•\tI like how you related the GAM to your estimation problem. It is helpful to include the Kolmogorov-Arnold representation theorem to explain how the restricted functional form is qualitatively approximating the form given in Theorem 1. However, given that this relationship is not quantified, statements like “Based on the Kolmogorov-Arnold representation theorem, generalized additive models …” (pg 4) should be scaled back. There is a similar statement in the conclusion.\n\n•\tFrom Figure 3, one could argue that KDE is a better method than APP for two-dimensional processes, despite its sensitivity to the kernel bandwidth.\n\n•\tIn Figure 4, I can’t tell which method has the best performance. It would be helpful to give a numeric comparison, perhaps some norm of \\lambda - \\hat{\\lambda}.\n\n•\tPlease define t^{j}.\n\n•\tPlease justify the use of the exponential function in the functional prior on \\lambda(t) in Equation (2).\n\n•\tPlease give some indications of computation time.\n\n•\tTypo on pg. 9: “The third-order method is able to period better”\n", "summary_of_the_paper": "This paper considers the problem of estimating the intensity function of a multi-dimensional Poisson process. Inspired by the generalized additive model, the authors propose a certain non-parametric form for the log-likelihood of the intensity function, involving multiple functions of groups of dimensions. Further, the higher-order interactions among the dimensions of the Poisson process are restricted to be within a small order k, so only functions of up to k of the dimensions are considered. Considering lower-order interactions naturally leads to a partial order, and hence the authors propose a log-linear model. For each element in the discretized state space, there is a parameter to be estimated. In turn, searching for the parameters amounts to solving a KL divergence minimization problem, which can be solved using convex optimization. Experiments on synthetic data as well as real data (New York taxis) are included. The experiments investigate the performance of the new approach, compared with three existing methods from the literature.", "summary_of_the_review": "The idea to decompose the function f into lower-order contributions appears to be novel, and the link to GAM’s is interesting. From the empirical evaluation, I am not convinced that the APP approach is better than existing methods for d = 2, 3, while higher-dimensional applications are not motivated.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "2: Several of the paper’s claims are incorrect or not well-supported.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a novel approach for estimating the high-dimensional intensity function of a Poisson process. The proposed approach builds on generalized additive models, using lower-dimensional projections. 

The reviewers noted that, although the paper is well written, the position of this paper compared to earlier related work is unclear, and the empirical evaluation of the method should be strenghtened. The authors clarified some points in their response, but the paper would still require some more modifications to be ready for publication. I therefore recommend this paper to be rejected.
Decision: Reject