Additive Poisson Process: Learning Intensity of Higher-Order Interaction in Poisson Processes
Authors: Simon Luo
Feng Zhou
lamiae azizi
Mahito Sugiyama
Venue: Submitted to ICLR 2023 (2023)
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 2023
ADDITIVE POISSON PROCESS: LEARNING INTENSITY OF HIGHER-ORDER INTERACTION IN POISSON PRO- CESSES
Anonymous authors Paper under double-bling blind review
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 tech- niques from information geometry to model higher-order interactions tical manifold and in generalized additive models to use lower-dimensional pro- jections 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 effectively uses samples observed in lower dimensional space to esti- mate a higher-ordee intensity function with sparse observations.
statis-
INTRODUCTION
The Poisson process is counting process used in wide range of disciplines such spatial- temporal sequential data in transportation (Zhou et al.] 2021). finance [Ialaan 2016) and ecology [Thompson] 1955) model the arrival rate by learning an intensity function. For given time inter- val, 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 prac- tical applications, learning the multi-dimeesionn intensity function is challenge due to the sparsity of observations. Despite the recent advances of Poisson processes, current Poisson process models are unable to learn the intensity function of multi-dimensional Poisson process. Our research ques- tion is. "Are there any good ways of approximating the high dimensional intensity function?" Our proposal, the Additive Poisson Process (APP), provides novel solution to this problem.
Throughout this paper, we use running example a spatial-temporal setting. Say we want to learn the intensity function for taxi to pick up eustomers at given time and location. For this setting, each event is multi-dimensional:: that is, (x, y, W). where pair of x and y represents two spatial coordinates and W represents the day of the week. In addition, observation time & is associated with this event. For any given location or time, we can expect at most a few pick-up events, which makes difficult for any model to learn the low-valued intensity function. Figure visualizes this problem. In this problem setup, if we would like to learn the intensity function at given location (ar, y) and day of the week W. the naïve approach would be to learn the intensity (I, y, W) directly from observations. This is extremely difficult because there could be only few events for given location and day.
However, there is useful information in lower-dimensional space; for example, the marginalized observations the location (xr, y) across all days of the week, or on the day W all locations. This information can be included into the model to improve the estimation of the joint intensity function. Using the information lower-dimensional space provides a structured approach to include prior information based on the location or day of the week to improve the estimation of the joint intensity function. For example, a given location could be shopping center or a hotel, where it is common for taxis to pick up passengers, and therefore could also be additional patterns that could be uncovered based on the day of the week. We can then the observations of events to update knowledge of the intensity function.
expect more passengers at this location. There

---

Under review as a conference paper at ICLR 2023
Previous approaches, such kernel density estimation (KDE) (Rosenblat) 1956), learn the joint intensity function using information in lower dimensions. However, KDE suffers from the curse of dimensionality, which means that it requires large size of samples to build accurate model. In addition, the complexity of the model expands exponentially with respect to the number of dimen- sions, which makes it infeasible to compute. Bayesian approaches, such as using a mixture of beta distributions with a Dirichlet prior (Kottas] 2006; Kottas & Sansó (007) and Reproducing Kernel Hilbert Space (RKHS) (Flaxman al., 2017), have been proposed to 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 Poisson processes. Our model combines the techniques introduced by [Luo & Sugiyama 2019) to model higher-order interactions between Poisson processes and by Friedman & Stuetzle (198T) in generalized additive models to learn the joint intensity function using samples in lower dimensional space. Our proposed approach is to decompose
Poisson process into lower-dimes dinenssona representations. For example, we have points (zz), in the z-dimension
and (v.) 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 differrnt 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 (0166) 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 Kullback-Leibler (KL) divergence from the sample distri- bution in a snwer-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 inten- sity functions using ower-dimensions. projections, thanks to the DolmogorrvvArnolo representation theorem. This property makes advantageous to use our proposed approach for cases where there are 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-lineas model and the parameters of the model are fully customizable to 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.
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 a tight link 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 À : E + R. An inhomogeneous 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 counting process N(t),
which has an independent increment property. For any & 2 O and infinitesimal interval 5 2 0. probability of events count is (t) (N)) = = (X(t) + (6)), ((N((t)) N(t) = 1)
8A(t) + o()), and p(N(t- + 8) N()) > 2) o()). where o(-) denotes little-o notation (Daley & Vere-Jones 2007). To take (inhomogeneous) multi cuulti-dimensional attributes into account, we consider multiple intensity functions, each of which is given as X R - IR associated with a subset J of the domain of possible states [D] = {1...., D). Each J C [D] determines the condition of the occurrence of the event. To flexibly consider any combination of possible states, [D] is composed
2

---

Under review as a conference paper at ICLR 2023
of possible states across all dimensions. For example, in the taxi pick up example in Introduction, domains z, y. W are discretized to disjoint sets as {l...., Dz}. (D₂ +1.... D2 + Dy). {Dx + Dy + +l..... D2 D, + Dw respectively, without loss of generality and gathered as single set (D) with D = D2 + Dy + Dw. For each J - [D), the likelihood of this model is given by
p ((W) As (1)) exp (-/ As(()dtt
II A (t).
We define the functional prior on (t) as
ty = (f)(t)) exp (f)(t)).
(2)
The function g(-) is a positive function to guarantec 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 described as
log p As (A)) E ( (fa (f) dt.
(3)
In the following sections, introduce the generalized additive models and propose model it by the log-linear model to learn ft) and the normalizing term J exp (fa (t)) dt.
2.1 GENERALIZED ADDITIVE MODEL.
We present the connection between Poisson processes and the Generalized Additive Model (GAM) proposed by |Fieedman & Stuetzle (198I). The GAM projects lower-dimensional space apply smoothing functions build restricted class of non paramarmetr regression models. GAM is less affected by the curse of dimensionality compared to directly using smoothing in a higher-dimensional space. For given set of processes C [D] = [ D), the traditional GAM using one-dimensional projections is defined as log A₂(1) = - fz(()- B, with some smoothing function fz-
features into
In this paper, we 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 as
log As(1)= = ESOO()) + E400..00 + .... 50 B, =
E ICA.IS
f7(1) 35,
JEJ
J1J2EJ
(4)
The function f1 : R R smoothing function to the data, and the normalization constant By for the intensity function is obtained as B, (t)dtt fexp( Err,, |T| (t) )dt. The definition of the additive model is in the Equation 3 and Equation we can see that the smoothing function J in Equation 3 realized as the summation over dimer-dimeniionall projections in Equation
form Equation In particular, if compare
Learning of continuous function using lower-dimensional projections is well known because of the Kolmogorov-Arnolu representation theorem, which states that:
Theorem 2.1 (Kolmogorov-Ariin Representation Theorem (Braun & Griebel, 2009 Kol-
mogorov 1957)). Any multivariate continuous function can be represented as a superposition of one dineennnsio functions; that is, (t.....,ty) =Zant fy (Ep=194> (tp)
9q.
Braun (2099) showed that the GAM an approximation to the general form presented in Kolmogorov-Amolk representation theorem by replacing the range {1...., 2n + 1} with C and the inner function Da. by the identity if p and zero otherwise, yielding 15(1) 1-()=En Erg,
Interestingly, the canonical form for additive models in Equation 4 can be rearranged to be in the same form as Kolmogorov-Arool representation theorem. By letting f(t) = EEgg fi(t) = (A,((t) and g(-) = exp(:).
Xt(t) = exp (3,) exp EE (((() exp (i)
(5)
3

---

Under review as a conference paper at ICLR 2023
M
(0,z)
==2
t=1
Order of inter interaction
beoveen processes
(1,2,3)
1p
Example for
(1,2)
Maaysss (1)
(a)
(b)
Time
Figure 2: (a A visualisation of the truncated pa-
rameter space approximate the joint intensity
Figure 1: Partial order structured sample space (N, 3) with D 3. Each node represents state and and the directed edge represents the direction of the partial ordering.
function with = and k: = 2. (b) A visualiza- tion of the imput datasets, where the blue points represent events with two spatial dimensions and one time dimension.
where we assume ft) 0 if if [1] > k for the k-th order model and 1/ exp(3)) is the normal- ization term for the intensity function. Based on the Kolmogorov-Arool representation theorem, generalized additive models are able to 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
log ((t)))
E
(t) SC
where 8' is a constant given by 3' (tt 31. In the following subsection introduce log-linear formulation equipped with partially ordered sample space, which aligns with the GAM formulation Equation
2.2 ADDITIVE POISSON PROCESS
We introduce key technical contribution in this section, a log-linear formulation called the additive Poisson process to estimate the parameters for the higher-order interactions in Equa- tion We begin by discretizing the time window [0, into M bins and treat each bin as natural number E [M] {1, ........ ,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 high-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,r) J E 2(D) \ V, E [M] U {(L,0) where the symbol 1. denotes the least el- ement in our partial order structure, which required normalize probabilities the resulting log-linear model. We define the partial order 3 (Davey & Priestley 2002) on Ω as
w = (J,T) - = JC J and r S -', for each w, w' E S2
and (T, 0) < w for all w E Q which illustrated in Figure The relation J S J' is used to model any-order interactions (Amari 2016, Section 6.8.4) between Poisson processes and each 7 in (J.T) represents the auto- auto-regressive component ("time") in our model. Each node w in the partially ordered set (poset) represents the state of the sample space and the arrows in Figure represent the partial order relationship between two node that is, if w - w', then w 3 w'.
Intuitively, the greatest node for each sents the multi-dimmensionl Poisson process. Other nodes represent projections onto lower- dimensional space that correspond the marginalized observations; for example, {{1), (2), {3}}
[M]. which is ((1,2,3),T) in Figure repre-
and {{1, 2}, {1, 3), (}} represent the first- and second-order processes. Using our example in
In information geometry, this corresponds hypergraphs which have simplicial complex structure Ay et al.
2018 g Section 2.9).
"This 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 2023
Introduction, where we wanted estimate the intensity function of pick-up event of taxi, {1} and (2) correspond to spacial coordinates z and y, respectively, and {3} to the day of the week W. and represents the (discretized) observation time. We can then update our belief to model the second-order intensity function using observations of the second order events. For example, (1,2), (1,3), (2 represents event occurring at {ar, y), {ar, W). and {y. W). We then continue this process to an arbitrary order of interactions. Later on in this section we introduce the mathematics to estimate the higher-order function using restricted number of lower-dimensional projections.
On any set equipped with a partial order, we can introduce log-linear model (Sugiyama et al. 2017). Let us assume that a parameter domain SC is given. For partially ordered set (Ωt 3 ) the log-linear model with parameters (fs)sEs introduced
Algorithm Additive Poisson Process (APP)
Function APP([ S, M, h):
2: Initialize Z with the number M of bins
Apply Gaussian Kernel with bandwidth on {t,),=1 to compute p
Compute r = (fss)Es from p
logp(w; 0) --E.. EE 16 - (0)
(8)
5: Initialize (0s)ses (randomly or 82 = =0) repeat
for each w E 2, where [ = 1 if the state- ment in [-] is true and 0 otherwise, and (0) E R is the partition function uniquely obtained
7: Compute using the current 0 = =(0))es
8:
Compute 17 = =m7)ses from
as (0) = log Ewen exp() Eves 0s) - -0 A special case of this formulation co- incides with the density function of the Boltz-
An f 17 - -)
10:
Compute the Fisher information matrix G using Equation [T
mann machines (Sugiyama et al. 2018 Luo & Sugiyama 2019).
11:
0 <- 0 G An
12:
until convergence of (0.)sEs
Here there is clear correspondence between the log-linear formulation and that in the form of Kolmogorov-Amot representation theorem in Equation if we rewrite Equation 8 as
13: End Function
p(w; 0) exp (0) exp Σ SES O.)
x exp (E 1 1l=
(9)
We call this model with (S, 3) defined in Equation the additive Poisson process, which represents the intensity X the joint distribution across all possible states. The intensity X of the multi- dimensional Poisson process given via the GAM in Equation 5 is fully modeled (parameterized) by Equation & and each intensity fr(.) is obtained
0(,,)- To consider the k-th order model, consistently use the parameter domain S, given as S (J, Ω |J| < k 3. where k is input parameter to the model that specifies the upper bound of the order of interactions. This means that 8, 0 for all S. Note that our model is well-defined for any subset S Ω and the user use an arbitrary domain in applications. A visualization of the truncated parameter space is shown in Figure Ia
For a given J C [D] and each bin 7 with w = (J, T), the empirical probability f(w) of input observations given
N
p(w) 2 Eo(t). Z ICJ
Z Σ p(w), WEN
(T(T)
K
7 1. h1
(10)
Nh
for each discretized state w = (J, The function o, performs smoothing on time stamps fi....,to IN. which is the kernel smoother proposed by Buja et al. (1989). The function K is a kernel and h2 is the bandwidth for each projection / I [D]. We use the Gaussian kernel as K to ensure that probability always nonzero, meaning that the definition of the kernel smoother coincides with the kernel estimator of the intensity function proposed by (Schabe (193). The normalization constant
for the intensity function can be computed by B, = Σ15 h1
2.3 OPTIMIZATION
Given an empirical distribution P defined in Equation 10 the task is to learn the parameter (0,)sEs such that the distribution via the log-linear model in Equation is as close to much as possible.
5

---

Under review as a conference paper at ICLR 2023
Let us define 6g {pl 0, 0f 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 minpess DKI. (p. p). where (p. p) Ewen P log(p/p) /p) is the KL divergence from to p. In this optimization, let p" be the learned distribution from the sample with infinitely large sample size and let p be the learned distribution for each sample. Then can lower bound the uncertainty (variance) [[DKL (p",p) ,p)] by |S|/2N (Barron & Hengartner, 1998).
Thanks the well-developed theory of information geometry (Amari 2016) for the log-linear model (Amari 200T), it is known that this problem can be solved by e-projection, which coincides with the maximum likelihood estimation and is always conver optimization (Amani 2016 Chapter 2.8.3). The gradient with respect to each parameter 0g is obtained by (0/00.) DKL (P. p) 7/ ne.
where 7i = Ewen [[w=))((w). The value 7)s is known as the expectation parameter (Sugiyama et al. 2017) and 1/s 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 = {sI...., } and 0 10................. Os s) I" 7 s......., - 1T. 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 0 0 G G-1(n n). where the
Fisher information matrix G is obtained as
Ô (p.(P = Σ [wes)) - Ts Nsy. Dij a0,,00,, 00.,
(11)
WEN
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 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 is computationally positive semi-definite (Neal 1999).
The pseudocode for APP is shown in Algorithm The time complexity of computing line 7 is (([[[))). This means when implementing the model using gradient descent, the time complexity of the model is (( to update the parameters in S for each iteration. For natural gradient, the cost of inverting the Fisher information matrix G is (( therefore, the time complexity to update the parameters in S is ([SS] + [ 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, it is more efficient to use gradient descent.
3 EXPERIMENTS
We perform experiments using two-dimensional synthetic data, diche--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 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 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) (Rosenblati 1956), an inhomogeneous Poisson process whose intensity is estimated by reproducing kernel Hilbert space formulation (RKHS) (Flaxman et al.) (017). and a Dirichlet process mixture of Beta distributions (Dt- (Kottas 2006 Kottas & Sansó] 2007). The hyper-parameters M and h in our proposed model are selected using grid search and cross vross-valiattion. 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) 2013) and M could be selected empirically from the input data by computing the time interval of the joint observation.
'The code is available in the supplementary material and will be publicly available online.

---

Under review as a conference paper at ICLR 2023
Presse
Presest
(a) Dense observations N
100, 000.
Raroe
(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., (blue) and h: 2 (orange).
Presee
(a) Dense observations,
= 100, 000, 0.4.
Pap mr 11
Paom
21.77
(b) Sparse observations, N 100, h 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 (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 ([I) 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- likellhoo for the New York Taxi data. Single processes ([T) and [W]) and joint process of them ([T,W]). APP-# repre- sents the order of the Additive Poisson Process.
Process
APP-1
APP-2
KDE
RKHS
DP-beta
[T]
713.77
731.01
Jan
[W1
745.60
745.60
745.23
853.42
790.04
[TW]
249.60
246.05
380.22
259.29
260.30
Process
APP-1
APP-2
KDB
RKHS
DP-beta
4.98e-5
4.9c-5
2.81e41
Feb
[W]
[T]
738.66
738.66
773.65
755.71
$11.34
792.10
765.76
Dense
[2]
2.83e-5
2.83-55
L.764
[TW]
328.84
244.21
307.86
334,31
326.52
[122]
2.98--2
1,27e-3
6.33e-4
4.09e-2
4.54e-2
[T]
716.72
716.72
733.74
755.
741.28
Sparse
[21
[1]
2.28e-4
7.26e-4
2.88e-2
2.28c-4
7.26c-4
1.77e-2
3.67633
2.76c-1
8.83e-4
2.35e-3
1.96-22
1.84e-2
2.19c-3
2.62e-3
3.68e-2
[TW]
[W]
738.06
291.20
246.19
738.06
$16.99
289.69
853,33
328.47
832.43
300.36
3.1 EXPERIMENTS ON TWO-DIMENSIONAL PROCESSES
For our experiment, we use 20 Gaussian components and simulate a dense case with 100,000 obser- vations and 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 3a 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 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 is able directly learn
7

---

Under review as a conference paper at ICLR 2023
Prasom
1..3
1-2.311
Fromuy D.2.A.4
Told Sharpis
(a) Dense observations, N 10"
Proset 131
Prome 2
Roae 2
Trall ELDwayen
(b) Sparse observations, N 10
Figure 5: KL Divergence for fourth-order Poisson process. We selected four representative examples for our experimental results, full results available in the supplementary material. The line color signifies the order of the model, i..., k (bluc), As 2 (orange), k 3 (green) and k (red).
(a) Dense observations, N 10?
112.144
Sparse observations, N 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.c., k 1 (blue), k- (orange), k- (green) and k (red).
the intensity function from the projection onto the two- -vo-dimensional space. In contrast, the first-order model must approximate the second-orde 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 state-of-the-ar approaches. APP performs the best for first-order processes both the sparse and dense experiments. Experiments for RKHS and DP-beta were unable 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 Ba 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 leamn a more accurate intensity function when the ground truth is unknown.
3.2 EXPERIMENTS ON HIGHER- DIMEERSIIDEN 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 2023
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 2023
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 10-1800-1225, 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 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 Chi and Tamara Kolda. On tensors, sparsity, and nonnegative factorizations. SIAM Journal on Matrix Analysis and Applications, 33(4):1272-1299) 1299, 2012.
Daryl J 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 Pricstley. 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, 1102):081151044 2017.
Jerome H Friedman and Werner Stuetzle. Projection pursuit regression. Journal of the American
Statistical Association, 76(376))817-823 1981.
Deniz Ilalan. A poisson process with random intensity for modeling financial stability. The Spanish Review of Financial Economics, 14(2):43-00, 2016.
ST John and James Hensman. Large-scale cox process inference using variational fourier features. In International Conference on Machine Learning, pp. 2362-2370, PMLR, 2018.
Andrei Nikolaevich Kolmogorov. On the representation of continuous functions of many variables by superposition continuous functions of one variable and addition. Doklady Akademii Nauk, 114(5):953-99 99665 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.
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):315-331 31633313,, 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. Bias-variance trade-off in hierarchical probabilistic models ing higher-order feature interactions. Proceedings of the 33rd AAAI Conference on Artificial Intelligence, pp. 4488 4495, 2019.
10

---

Under review as a conference paper at ICLR 2023
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-1087 2001.
Yosihiko Ogata. On Lewis' simulation method for point processes. IEEE Transactions on Informa- tion Theory, 27(1):33-31, 1981.
Murray Rosenblatt. Remarks on some nonparametric estimates of density function. The Annals of Mathematical Statistics, pp. 832-837, 1956.
Yves-Laurent Kom Samo and Stephen Roberts. Scalable nonparametric bayesian inference on point processes with gaussian processes. In International Conference on Machine Learning, pp. 2227- 2236. PMLR, 2015.
Schäbe. Nonparametric estimation intensities of nonhomogeneous poisson processes. Statisti- cal Papers, 34(1):113331)) 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 Learming, pp. 2810-2819. PMLR, 2016.
David W 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 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 to ecology. Biometrika, 42(1/2):102-115,, 1955.
Feng Zhou, Simon Luo, Zhidong Li, Xuhui Fan, Yang Wang, Arcot Sowmya, and Fang Chen. Efficient em-variational inference for nonparametric hawkes process. Statistics and Computing. 31(4):1-11, 2021.
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 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 ) for each dimension and increases exponentially with respect with additional input dimension, that is, ((n))d)), where n refers the number of observations and df 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 to estimate any nigh-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

---

Under review as a conference paper at ICLR 2023
is (n?)k for each dimension and is exponential with respect to the number of input dimensions, which is still infeasible to estimate any high- eigh-dimensional intensity function.
Variational inference (Lloyd al. 2015) approaches that can be used to make Bayesian inference of Poisson process much more efficient scaling it up to the linear complexity ((n) for each dimension. One of the most notable work in this is using Fourier Features in the Cox Process (John & Hensman 2018). However variational inference is not guaranteed to asymptotically converge to the true posterior distribution. Because variational inference does not have any theoretical guarantees, it is likely to fail in extremely high dimensions and sparse observations.
Our approach uses the discretization approach scale up the model to higher dimensions. We use a graph (partial order) structure to allow the flexibility for domain expertise to specific how each dimension is treated and which interaction effects should be included into the model. Unlike varia- tional inference based approaches, our estimation procedure of the intensity function has theoretical convergence which is based by the Kolmogorov-Aroold theorem.
A.2 POISSON FACTORIZATION
Our work is closely related to Poisson Factorization (Chi & Kolda 2012), where random variables in tensor are represented with a Poisson distributed or Poisson process likelihood. The tensor is usually used represent some high-dimensional dataset such as contingency tables or other collec- tion of counting datasets, which are often large and sparse. The objective of Poison Factorization decompose the high- nigh-dimensional 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 spatial-tempont problems consisting of sparse count data. One example of this work is Bayesian Poisson Tucker decomposition (BPTD) (Schein et al. 2016). where dataset of interaction events 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 spatial-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.
B ADDITIONAL EXPERIMENTS
B.1 BANDWIDTH SENSITIVITY ANALYSIS
Our first experiment is to demonstrate the ability for our proposed model to learn an intensity func- tion 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 toy problem for our model. This experiment is to observe the behaviour of varying the bandwidth in our model. In Figure 7al we observe that apply- ing no kernel, we learn the deltas of each individual observation. When apply Gaussian kernel, the output of the model for the intensity function is much more smooth. Increasing the bandwidth of the kernel will provide a wider and much smoother function. Between the 60 seconds and 80 seconds mark, it can be seen when two observations have overlapping kernels, the intensity function becomes larger in magnitude.
B.2 ONE DIMENSIONAL POISSON PROCESS
A one dimensional experiment is simulated using Ogata's thinning algorithm (Ogata 1981). We gen- erate two experiments use the standard sinusoidal benchmark intensity function with frequency of 20T. The dense experiment has troughs with 0 intensity and peaks 201 and the sparse experiment has troughs with 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 KI.
12

---

Under review as a conference paper at ICLR 2023
Algorithm 2 Thinning Algorithm for non-homogenous Poisson Process
1: Function Thinning Algorithm (A (t), T):
2: n m 0, lo = 80 0, suPos.sT (1)
3: repeat
4: #f ~ uniform (0, 1)
5:
w w= R {w exponential {)]
6:
8m+1
+
7:
D ~ uniform (0, 1)
D < Ammm1 then
8:
9:
bn+1 5m+1
10:
n =n++
else
12:
m m +
13: end if
14: if In < T then
15:
return {tk]k=1,2....
else
17:
return {tR m-
18:
end if
until 8mm ST
20:
End Function
1.0
No Kerne
Gaussi an Karnel, 2n
Gaussi an Karnel, hn
Divergence evergence
KX Divergence
0.100
KL Divergence
0.2
0.6
Gaussi an Karnel, hw
0.4-
Gaussi an Karnel, he 10
0.050
0.025
-0.1
10
Kernal Bandwith
Kernel Bandwith
(a) Toy Example
(b) KL divergence of dense experi- (c) KL divergence of sparse experi-
ment
ment
Yew 1407
(d) Ogata's thinning algorithm with high intensity
(c) Ogata's 's thinning algorithm with low intensity
Figure 7: One dimensional experiments
divergence between the underlying intensity function used generate the samples the intensity function generated by the model. Figure shows that the optimal bandwidth is h = 1.
13

---

Under review as a conference paper at ICLR 2023
OO An TA LAT
25 OF OO DO LEE TN L. BRNEEE
Prosee I. 2.4)
LI DO ES LEE Baseeit
KImH Samel Berteeth
M
Froomet
Prooees: NI
OO EA LE
ES
Figure 8: KL Divergence for four-order Poisson process.
ZZT SO 3xE
Lo 2n Zie
Prosse 1..)
Ler baroedsa
LET TO Sn basdee dwidd
Kern Berewits
1.30
PRICEE
Prices
Dense observations.
(b) Sparse observations.
#00dSma0 TX
#XUndamAD TM
12 TM
sooddeeo TM
12 OF On TH 19 311 31
OF TY L/ SR Nomel Bunde inddwib
OE AT OF TE 1I
Preae 41
Karal Hindwets
4)
2,34
Prices::
S0 OE Do Tee 21 LIT
OT TO DE LTE
50 OI DE DO TOR EET L12 BEALCC
DA DE DH
Pro 20
Kuraa Bantaith
Drssnn 22.4 41
Preeeass (3)
LA
Rarial aentiich 21
Koral Satticc
43
Brich
21
14

---

Under review as a conference paper at ICLR 2023
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
15

---

Under review as a conference paper at ICLR 2023
10
10
Orler: Ground
Process: 12 31
Process: 2.
Process: (4)
Onten: Onder: Order: 3
Time
Time
Time
imnuiitty
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
M Aapsve
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
16


Reviews:
[{"review_id": "iKbDL1V2MQ", "reviewer": "Reviewer_x4an", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "This paper introduces an inference model for learning the intensity function of a multidimensional Poisson process. The approach is presented as an application of generalized additive models to Poisson processes. One key point comes from information geometry for the log-linear model which shows that the likelihood function of the model at stake is convex.", "summary_of_the_review": "The approach followed in this paper is in my opinion interesting but requires further developments, both theoretical and numerical, to be published in this conference.", "strength_and_weaknesses": "This paper makes an interesting connection between multidimensional Poisson processes, generalized additive models, and information geometry. The approach seems to me original and deserves to be studied.\n\nHowever, the two key points of the paper (the model itself, obtained as an application of generalized additive models, and the convexity of the optimization problem) do not require any deep theoretical or empirical development. This is not necessarily an obstacle but, in this case, I would expect at least some theoretical guarantees on the quality of the estimate obtained by the method, which are lacking.\n", "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": "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": "EhAMs-6CNw", "reviewer": "Reviewer_Cwen", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "I reviewed this same paper last year for ICLR 2022. By comparing the current manuscript and the manuscript from last year, I did not find significant differences. Therefore, I will reiterate my questions from last year.\n\nThis 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 real data analysis.\n\n ", "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": "Strength: The paper is overall well written and the presentation is clear. \nWeakness:  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, they 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 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 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:\n\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.\n\nReference:\n[1]. Wood, S.N., 2017. Generalized additive models: an introduction with R. CRC Press. (2nd Edition)", "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": "BIG7iEY-gN", "reviewer": "Reviewer_2omy", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper proposed Additive Poisson Process (APP), which models higher-order interaction effects of the intensity functions in Poisson processes using projections into lower-dimensional spce. The intensity function is modelled using a log-linear model, whose parameters are learned by minimizing the KL divergence between a sample distribution in lower-dimensional projections and the distribution with an intensity function. The empirical experiments showed promising results on synthetic datasets and a real world dataset. \n", "summary_of_the_review": "The paper presents a framework for modelling the intensity function in high-dimensional Poisson processes. There are potential improvement area in terms of the presentation, more thoroughly analysis/experiments and discussion, but overall I think it is an interesting idea. \n", "strength_and_weaknesses": "Strength:\n* The idea is interesting. The problem of modelling the intensity function in multi-dimensional Poisson processes is challenging due to sparse observations. Approximating the high-dimensional intensity function using low-dimensional projections and framing it in an GAM form is a nice idea.\n\nWeakness:\n* a log-linear model is proposed in APP, which is a bit restrictive in its functional form. How easy it is to extend the method to richer class of (non-linear) models? \n* Is it not obvious to me that such a low-dimensional representation always exists. Related to this, How big of k (i.e., the truncation level of the order) is sufficient? Have you thought about ways to improve/constrain APP such that k needed can be low? On the contrary, as the authors also pointed out, large k may result in overfitting, how should one pick k in practice? \n* The hyperparameters M and h are learned with grid search, which may be expensive in high-dimensional problems. \n* The experimental results are all on synthetic data expect one on the New York taxi dataset. More extensive results will make the paper stronger. \n* the number of parameters grows exponentially with the dimension and there is no structural dependency between the theta's. Will it help to enforce some dependency between them (for examples, theta_{1,2} depends on theta_1 or theta_2)? \n\n", "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": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "6eP6JdfjJqs", "reviewer": "Reviewer_8qDw", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The authors consider the problem of learning the intensity function of an inhomogeneous Poisson process (PP) over higher dimensional spaces:\n 1. The authors propose a certain family of generalized additive models (GAM) to model the log-intensity function of the PP and call their method the additive Poisson process (APP).\n 2. They propose to learn the parameters of the GAM by minimizing the Kullback-Leibler (KL) divergence between a tractable kernel density estimate (KDE) approximation to the empirical intensity and the intensity given by the GAM.\n 3. The authors test their method on two synthetic and one real-world dataset.", "summary_of_the_review": "The authors present a fine idea, but the presentation quality is low, and the empirical results seem weak.", "strength_and_weaknesses": "## Strengths\nThe authors consider a relevant problem, and their proposed solution is sound and holds some promise for estimating the intensity functions of low- to medium-dimensional PPs.\n\n## Weaknesses\nWhile the paper has some good ideas, some of the presented theoretical connections are dubious, and I am not convinced that the empirical results bear out all of the authors' claims. Concretely:\n - The authors choose to model the log-intensity function with a GAM, but I do not see how this modelling choice is supposed to constitute any sort of deeper connection. I am hung up on this point because the authors seem to dedicate a significant part of Section 2 to this claimed connection, which appears to only consist of trivial algebraic manipulations.\n - The authors' claim that their method scales better to higher-dimensional problems is dubious. While strictly speaking, their learning process does scale better than other methods, it still scales exponentially in the dimensionality of the space. Hence, their method is limited to moderate-dimensional problems.\n - Some of the empirical results on problems that are supposed to be favourable for APP do not seem convincing. A well-known issue of GAMs is that due to their flexibility, they can heavily overfit, which appears to be precisely the case in Figure 4. Especially in the right plot, the APP intensity is wholly concentrated on the data points, while the KDE fit seems far more reasonable.\n - Similarly, as shown in Figure 3, KDE outperforms APP on the process with higher-order interactions. The authors note this but claim that \"Figures 3a and 3b 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.\" - However, I don't think Figure 3 demonstrates that APP would be much less sensitive to the kernel bandwidth than KDE. \n - I also think the name for the authors' method is quite misleading. \"Additive Poisson process\" would suggest that the family of processes described is a subclass of PPs for which some nice additive property holds (e.g. stability under addition). However, this is not the case, and I would like to ask the authors to pick a more fitting name.\n - The authors claim, \"Based on the Kolmogorov-Arnold representation theorem, generalized additive models are able to learn the intensity of the higher-order interaction between Poisson processes by using projections into lower dimensional spaces.\" - This cannot be true, as GAMs use a more restricted family of functions, so the K-A theorem does not necessarily apply. It would be necessary to discuss the approximation properties of the particular family of GAMs they use to make this claim precise.\n - The experimental methodology is not properly described in the case of synthetic data experiments. It is unclear what data the methods were tested from which the plots were generated.\n - The authors invoke a lot of mathematical machinery, e.g. the poset formulation for the model and information geometry, without seemingly relying on it too much. I think the method and training process could be described in much simpler terms which would heavily improve the paper's readability.\n - The authors claim to use natural gradients to solve their proposed optimization problem. However, they never demonstrate the utility of this over standard gradient descent.\n\n## Questions\n - \"Our model combines the techniques introduced by Luo & Sugiyama (2019) to model higher-order interactions between Poisson processes and by Friedman & Stuetzle (1981) in generalized additive models to learn the joint intensity function using samples in a lower dimensional space\" - a lower dimensional space compared to what? Also, this sounds very similar to the approach in the paper. How do they actually differ?\n - \"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...\" - the previous sentence is not explaining a connection to anything. What are the authors referring to here?\n - What do the authors mean by \"higher-order intensity function\"? Do they mean higher-dimensional? Perhaps they meant to use it analogously to the phrase \"higher-order interaction\", but I don't think \"higher-order intensity function\" is standard terminology.", "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": {}}]

Author rebuttal:

Decision: Reject