K$^2$IE: Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
Authors: Hideaki Kim
Tomoharu Iwata
Akinori Fujino
Venue: ICML 2025 poster (2025)
Keywords: point processes
kernel methods
kernel intensity estimators
representer theorem
least squares loss
Primary area: general_machine_learning->kernel_methods

Abstract: Kernel method-based intensity estimators, formulated within reproducing kernel Hilbert spaces (RKHSs), and classical kernel intensity estimators (KIEs) have been among the most easy-to-implement and feasible methods for estimating the intensity functions of inhomogeneous Poisson processes. While both approaches share the term "kernel", they are founded on distinct theoretical principles, each with its own strengths and limitations. In this paper, we propose a novel regularized kernel method for Poisson processes based on the least squares loss and show that the resulting intensity estimator involves a specialized variant of the representer theorem: it has the dual coefficient of unity and coincides with classical KIEs. This result provides new theoretical insights into the connection between classical KIEs and kernel method-based intensity estimators, while enabling us to develop an efficient KIE by leveraging advanced techniques from RKHS theory. We refer to the proposed model as the *kernel method-based kernel intensity estimator* (K$^2$IE). Through experiments on synthetic datasets, we show that K$^2$IE achieves comparable predictive performance while significantly surpassing the state-of-the-art kernel method-based estimator in computational efficiency.

K'IE: Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
Hideaki Kim
Tomoharu Iwata 1 Akinori Fujino
Abstract
assess the risk experiencing events at specified domains and forecast the timings/locations of future events. Poisson processes have variety of applications in reliability engi- neering (Lai & Xie, 2006), clinical research (Cox, 1972; Clark et al., 2003; Lánczky & Györffy, 2021), seismology (Ogata, 1988), epidemiology (Gatrell et al., 1996), ecology (Heikkinen & Arjas, 1999), and more.
Kernel method-based intensity estimators, for- mulated within reproducing kernel Hilbert spaces (RKHSs), and classical kernel intensity estimators (KIEs) have been among the most easy-to-implement and feasible methods for esti- mating the intensity functions of inhomogeneous Poisson processes. While both approaches share the term "kernel", they are founded on distinet theoretical principles, each with its own strengths and limitations. In this paper, we propose a novel regularized kernel method for Poisson processes based on the least squares loss and show that the resulting intensity estimator involves spe- cialized variant of the representer theorem: it has the dual coefficient of unity and coincides with classical KIEs. This result provides new theoretical insights into the connection between classical KIES and kernel method- based intensity estimators, while enabling us to develop an ef- ficient KIE by leveraging advanced techniques from RKHS theory. We refer to the proposed model as the kernel method hased kernel inten- sity estimator IE) Through experiments on synthetic datasets, we show that K²IE achieves comparable predictive performance while sig- nificantly surpassing the state-of-the-art kernel method-based estimator in computational effi- ciency.
Kernel intensity estimators (KIEs) are the simplest non- parametric approaches to estimating intensity functions (Ramlau-Hansen, 1983; Diggle, 1985), with advantages that include superior computational efficiency and theoret- ical tractability. They represent the underlying intensity function as a sum of smoothing kernels' evaluated at data points, where rescaled versions of density functions are usually adopted as smoothing kernels to correct the edge effects, that is, the estimation biases around the edges of observation domains.
Recently, Flaxman et al. (2017) developed a feasible Re- producing Kernel Hilbert Space (RKHS) formulation for inhomogeneous Poisson processes. They showed that the representer theorem (Wahba, 1990; Schölkopf et al., 2001) holds for penalized maximum likelihood estimation un- der the constraint that the square root of the intensity func- tion lies in an RKHS: the obtained square root of the in- tensity estimator is given by linear combination of trans- formed RKHS kernels evaluated at data points. The trans- formed RKHS kernels, often referred to as the equivalent RKHS kernels, naturally account for edge effects through likelihood functions, and the intensity estimator has been shown to outperform KIES scenarios where edge effects are prominent, such as in high-dimensional domains. Al- though the kernel method-basee baathdd-basdd intensity estimator has form similar to KIEs, it requires fitting the dual coefficient using gradient descent method, making it less favorable than KIES in terms of computational efficiency.
1. Introduction
Poisson processes have been the gold standard for modeling point patterns that occur randomly in multi- dimensional domains. They are characterized by an inten- sity function, that is, the instantaneous probability of events occurring at any point in the domain, which allows us
In this paper, we propose a penalized least squares loss for- mulation for estimating intensity functions under the con-
'TTT Corporation, Japan. Correspondence to: Hideaki Kim <hideaki. kin kinttaiimminnnt......... ntt.com>..
'Traditionall, the term "kernel" is used to refer both to the weights assigned to data points in kernel density estimation and to the desitiv--defiite kemelss that define RKHS. To avoid am- biguity, this paper employs two distinet terms "smoothing" and "RKHS" kernels, following Flaxman et al. (2017).
Proceedings of the 1nnd International Conference Machine Learning, Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
straint that the intensity function resides in an RKHS. The least squares loss is motivated by the empirical risk min- imization principle (van de Geer, 2000) and has demon- strated notable computational advantages in recent studies on Poisson processes (Hansen et al., 2015; Bacry et al., 2020; Cai et al., 2024). Utilizing advanced variational anal- ysis via path integral representation (Kim, 2021), we show that a specialized variant of the representer theorem holds for the functional optimization problem: the resulting in- tensity estimator, which we call the kernel method-based kernel intensity estimator (K2 has the unit dual coef- ficient and requires no optimization of dual coefficients given a kernel hyper-parane syper-parameter, which is consistent with KIEs; furthermore, it employs the equivalent RKHS ker- nels appeared in Flaxman's model and can effectively ad- dress edge effects based on the RKHS theory. This result not only establishes significant theoretical connection be- tween classical KIEs and kernel sethod-based intensity es- timators but also enables the development of a more scal- able intensity estimator based on kernel methods.
form binary classification to determine the occurrence of future events within S.
2.1. Kernel Intensity Estimator
Kernel smoothing is classical approach to nonparametric intensity estimation (Diggle, 1985), expressed as:
NV g(z, g x m//())) {xx
v(x) / g(x, s)ds, (3)
n=1
where g: g:::K x R+ represents a non non-negative smooth- ing kernel', and v(x) is an edge-correction term². This method, commonly referred to as the kernel intensity es- timator (KIE), is closely related to the well-known kernel density estimator (Parzen, 1962; Davis et al., 2011), espe- cially with bounded support (Jones, 1993). KIEs offer sev- eral advantages, including theoretical tractability, superior computational efficiency, and ease of implementation.
The smoothing kernel g(2, x') involves bandwidth hy- perparameter, which can be optimized using standard tech- niques such as cross-validation (Cronic et al., 2024) or Siverman's rule-of-thumb (Silverman, 2018). It is worth noting that widely-used cross- voossvvaliittion methods involv- ing test log-likelihood functions require the integration of intensity functions over a test region S C X, where KIES need rely on time-consuming Monte Carlo integration because -)v(() in (3) usually cannot be integrated in closed form (e.g., Gaussian smoothing kernels).
In Section 2, we outline related works of intensity estima- tion. In Section 3, we derive IE via functional analysis with path integral representation of RKHS norm. In Sec- tion 4, compare IE with conventional nonparametric intensity estimators on synthetic datasets, and confirm the effectiveness of the proposed method?, Finally, Section 5 states our conclusions.
2. Background
Let set of N point events, D == (I,,) } 1= being observed in dd-iimeniionll compact space, X' C Rd. We consider learning problem of intensity function in the framework of inhomogeneous Poisson processes, where intensity func- tion, A(z)= X + R4, represents an instantaneous proba- bility of events occurring at any point in X:
2.2. Kernel based Intensity Estimator
Let k X' x X + R denote continuous positive semi- definite kernel. Then there exists unique reproducing ker- nel Hilbert space (RKHS) Hr (Scholkopf & Smola, 2018; Shawe- Shawe-Taylor & Cristianini, 2004) associated with kernel k(.,)) Flaxman et al. (2017) modeled the intensity func- tion as the square of latent function in the RKHS,
X(x)= = lim (dx)) / |dz),
(1)
Idw -D
where (dw) is the number of events occurring in dæ C X. and| represents the measure of domain.
1(x) = (xx), ʃ() E Ha:
(4)
and proposed a regularized minimization problem with log- likelihood loss functional as follows:
One of the most significant applications of Poisson pro- cesses lies in evaluating the risk of experiencing events (e.g., traffic accidents and disaster events) within specified regions. Given intensity function X(x), the probability distribution of event counts (i.e., the Poisson distribution) over an arbitrary compact region S C X', denoted by (s).), is calculated as follows:
Σ log (og()(()))) + J Mlxcddo 7 11 }.
min JERL
(5)
where || 11 represents the squared Hilbert space norm, and y represents the regularization hyper-parameter. Through Mercer's theorem (Mercer, 1909). Flaxman et al. (2017) showed that the representer theorem (Wahba, 1990; Schölkopf et al., 2001) does hold in an appropriately trans-
(n( = A
[ )(z))xx,
A
(2)
n!
where f0, 1,,,
Using Equation (2), we can per-
Inttnsiity estimates near the edge of the domain X' biased downwards since no points are observed outside X.
*Codes are available at: https://gitub.com//HikkkKKK
2

---

Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
formed RKHS, resulting in the following solution:
Other Related Works
Gaussian Cox Processes (GCPs) provide Bayesian alter- native to kemnel intensity estimators and kernel method- based models, where Gaussian processes are used model latent intensity functions via positive-valued link functions (Møller et al., 1998). This framework enables principled interval estimation of intensity functions, along with hyper- parameter inference in a fully Bayesian manner (Rathbun & Cressie, 1994; Cunningham et al., 2007; Adams et al., 2009; Diggle et al., 2013; Gunter et al., 2014; Lloyd al., 2015; Teng et al., 2017: Donner & Opper, 2018; John & Hensman, 2018; Aglietti et al., 2019). Although GCPs are typically more computationally expensive than their non-Bayesian counterparts, extending kernel-based mod- els into the Gaussian process framework can yield effi- cient Bayesian alternatives. For instance, Walder & Bishop (2017) proposed Bayesian variant of Flaxman's model, known the permanental process, and Sellier & Dellapor- tas (2023) further extended this approach within the gener- alized stationary kernel framework.
N
f(x)=
anh(a, æn);
(6)
where a = (a1...... oN)" is the dual coefficient, and h XxX - RR is a transformed RKHS kernel defined in terms of the Mercer expansion of k( x') as
M(x
7m
Cm (2)))
1/7 + Am
m=1
(7)
Jx x
k(rr, s)mm(s)ds (8)ds Cm (2));
where {em (e)mm)mm is the eigenfunctions of the integral op- erator 1xx, (1x, s)ds. Recently, Kim et al. (2022) rewrote the definition (7) terms of a Fredholm integral equation of the second kind (Polyanin & Manzhirov, 1998),
1 ((x, a²) + J k(x,, s)h(s, a')ds k( x²),
(8)
While neural network- based methods often trade compu- tational efficiency for expressive power, Tsuchida et al. (2024) recently introduced the squared neural family, a model that simultaneously achieves expressiveness and an- alytical tractability. Like Flaxman's 's model, it ensures non- negativity and closed-form expressiveness of the intensity function through the use of a squared link function-an ele- gant property that merits further attention.
which enables us to utilize the established approxima- tion techniques for solving Fredholm integral equations (Polyanin & Manzhirov, 1998; Atkinson, 2010). We will discuss how to solve Equation (8) in Section 3.2. The trans- formed RKHS kernel is referred to as the equivalent RKHS kernel (Flaxman et al., 2017; Walder & Bishop, 2017; Kim et al., 2022).
The dual coefficient in the intensity estimator (6) solves the following dual optimization problem:
3. Method
3.1. Kernel Method-basse Kernel Intensity Estimator
N
N
min aERN
log Σ
an xwn) + La T Ha) (9)
In this paper, we introduce the least squares loss functional for Poisson processes (Hansen et al., 2015) given by
((xn) + / (ax))dar. dar.
where H [k( Ven Fwn) Mow' The computational com- plexity of solving (9) is naively O(gN2) for q iterations of gradient descent methods, but reduces to OggMN) when the equivalent RKHS kernel is given in degener- ate form with rank M (< N) such that h(xxx') EM (mm(()) (x))mm())))) (x').
(10)
This loss functional comes from the empirical risk mini- mization principle (van de Geer, 2000) and has demon- strated notable computational advantages in recent Pois- son process literature (Hansen et al., 2015; Bacry et al., 2020; Cai et al., 2024). For readers unfamiliar with the loss defined in (10), we briefly explain the origin of the term squares loss in Appendix A. We model the intensity function latent function in RKHS H,, and consider the problem of minimization of the penalized least squares loss as follows:
Flaxman's intensity estimator (6) differs from the classi- cal KIE (3) in that it does not require explicit edge correc- tion. Instead, the equivalent RKHS kernel I)(.-)) naturally accounts for the effects of the finite observation domain X d through the second term on the left-hand side of Equation (8). Flaxman et al. (2017) demonstrated that the intensity estimator (6) outperforms the KIE in terms of predictive performance, especially in high-dimensional settings. However, it demands the model fitting (9) unlike KIEs, which makes KIEs more favorable regarding compu- tational efficiency.
/ )xx) + 1IMIM }. (11)
min XER
2
A ((xx)
+
1=1
where
represents the squared Hilbert space norm, and 7 represents the regularization hyper-parameter.
||
11
3

---

Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
Through variational analysis, Theorem I below shows that the resulting kernel method-based intensity estimator is consistent with the classical KIE (3). While the main proof relies on the path integral representation (Kim, 2021), for completeness, we also provide an alternative derivation based on Mercer's theorem in Appendix B. To the best of our knowledge, this paper is the first to prove that the repre- senter theorem holds for the penalized minimization of the least squares loss for the intensity estimation in RKHS.
equal to the equivalent RKHS kernel h(::-) defined by (8): applying operator 2 to Equation (13) leads the relation,
1 5(x - =x)==')= q(xx + 7 J k- (xx s)g(s, x x')ds,
and applying operator K. to both sides of the relation yields:
20) S. k(x, s)q(s, x' dds + y (((2, ((222 x'),
The solution of the functional optimization problem (11), denoted as X-), involves the representer the- orem under a transformed RKHS kernel h(:--) defined by Equation (8), and its dual coefficient is equal to unity:
Theorem 1.
which is identical to Equation (8).
Theorem 1 demonstrates, under the least squares loss func- tional, strong connection between classical KIEs and modern kernel methods. From the perspective of KIE the- ory, Theorem 1 implies that the equivalent RKHS kernels h(---) are smoothing kernels constructed based on RKHS kemels. Hence, we call the proposed model (12) the kernel method-base based kernel intensity estimator (K² As Flax- man et al. (2017) discussed, the equivalent RKHS kernels implicitly incorporate edge effects in an effective manner. Therefore, our K'IE is expected to combine the compu- tational efficiency of KIEs with the effectiveness of Flax- man's 's kernel method- based estimator.
{(x) = Σ h(z, X.
(12)
Proof. Let (et) k(x, s) be the integral operator with RKHS kernel -), and K (x) = Sx k" (z,s)dd ssds be its inverse operator. Then, through the path integral rep- resentation of Gaussian processes (Kim, 2021), the squared Hilbert space norm can be represented in a functional form,
|||| 1AI
1|xc (xx (x,)))((x) s))( (s)dads.
xxxx
Similar to the conventional kernel method estimator (6), the support of K'IE in Theorem lies within the obser- vation domain X, i.e., concerns interpolation. However, by broadening the support of the RKHS kernel k(;-), the support in Theorem can be naturally extended: In other words, IE defined by Equation (12) can be applied in its current form to extrapolation as well. A proof of this claim is provided in Appendix C.
Using the representation, the objective functional in Equa- tion (11) can be rewritten as follows:
N
S(A) = 2 Σ ((x.) +
xxx (xx s))(()
n=1
where (-.-) is the weighted sum of k (s,)) and the Dirac delta function 8()),
Unlike conventional methods, KI has the limitation of not guarantceing the non-negativity of intensity functions. The equivalent RKHS kernels may generally take negative values, and since K- is constructed as linear combina- tion of the equivalent RKHS kernels, it can yield negative values in certain regions, particularly in areas with no ob- served events. This issue is caused by the fact that KPIE models intensity function by an RKHS function ʃ()) E Ha, while conventional methods by aff(-) for non-negative link function (-).. In practice, KºIE does not have large negative values because the second term of the objective function (11) penalizes them. Thus we can deal with the issue by applying max(u, 0) for intensity-related values u, such as u = X(x) and u = Ss A(xr) da over a domain S.
qº (x. 8) 8(x0 3) k" (x,
(13)
The solution of Equation (11), A(x)), is obtained by solv- ing the equation where the functional derivative of S()) is equal zero:
85
Σ
N
5(x -x,===)++) + +2 / q= s){(s)ss= =0. O.
2
==1
Let Q² (x) = Jx g² (x, s)ds be the integral operator associated with q" (, and (x(x) = Jx q(2, s)ds be its inverse operator. Then applying operator Q to the equation, 65/8İ 0, leads to a representation of the form,
N \(x) = Σ qe, xx,,
3.2. Construction of Equivalent RKHS Kernel
The primary task in K2IE is to derive the equivalent RKHS kernel that satisfies the integral equation (8). The method- ology varies depending on whether the observation do- main X' is infinite or finite, as elaborated in the subsequent
1=1
where the relation, (x)*)(x= = fxx -(xx=()))))) s)ds, was used. Furthermore, the following derivation shows that q(:))) is
4

---

Kernel Method-basee Kernel Intensity Estimators for Inhomogeneous Poisson Processes
sections. Here, we assume that RKHS kernels are shift- invariant, i.e., k(ar, x²) k(ar= 2²), which includes pop- ular RKHS kernels such as Gaussian, Matérn, and Laplace kernels.
using Nyström approximation (Williams & Seeger, 2000), which approximates the integral term through numerical integration. While this approach has the advantage of be- ing applicable to any RKHS kernel, it could potentially de- grade the accuracy of the edge correction because the in- tegral term is critical for the edge correction. To address the issue, we adopt the degenerate approach (Kim et al., 2022), which approximates RKHS kernels using 2M ran- dom Fourier features (Rahimi & Recht, 2007).
3.2. I. INFINITE OBSERVATION DOMAIN
If the observation domain is infinite, i.e., X = R. the inte- gral equation (8) can be solved by using the Fourier trans- form as follows:
2M
k(w) h(x x')== = F- x-1 k((w)
k( x') ~
Σ m=1
pm (2)) = x(x)) "T(()),) (16)
(x: x²),
(14)
Om (z)) M- cos (W x. + mm)
where (1)x)) denotes the inverse Fourier transform, and k(w E ( represents the Fourier transform of the shift- invariant RKHS kernel k(x -()x. Notably, the equivalent RKHS kernel :): is also shift-invariant due the sym- metry of the integral equation (8). Approximation meth- ods are required because the inverse Fourier transform in (14) generally cannot be expressed in closed form. One promising approach is the random feature map (Rahimi & Recht, 2007), where the equivalent RKHS kernel ap- proximated via Monte Carlo sampling from a probability distribution, p(-) x kl.)//6-11 k()),, such that f(x x') (eco:ppp/epp(ieeeeee /exp( T x)) cxp(iw" x')] Another feasible ap- EEw-p()
WmSM ~ k(w), WmnM = Wmn
0mgM = 0, 0mhi = 10r0==-1// -T/2,
and allows the integral term to be handled without any error as follows:
1(2,, x0) = +-(x) ( 12M + A) "d((',,
J Jx, ((x))(x)) dx,
(17)
A = Σ AJ, A
where I2a represents the identity matrix of size 2M. Notably, 2M x 2M matrix A, which involves the edge- correction, can be computed in a closed form:
proach is to apply the fast Fourier transform to (14).
When X' = Idd the edge-correction term in KIE (3) van- ishes, suggesting that the choice of the smoothing kernels g(ar, a²) KIE (3) effectively equivalent to the selec- tion of the equivalent RKHS kernels fa, x²) in IC (12). Through h(æ, æ²), however, we could find smoothing ker- nels more robust the squared error than popular ones such as Gaussian smoothing kernels. It is an interesting topic, but this paper focused on the case of a finite observa- tion domain, where edge correction plays a crucial role.
(A)) )mm' S/g @m (x)())))m (m)))
1
[S](wm + Wm, On + +mm)
2M
+ Cy(um - wmi, On - 0mn))
11 d
(i(w, 0) = cos
w (X. + x +
2 i=1
3.2.2. FINITE OBSERVATION DOMAIN
d
1 2w0 (X.g
x X
sinc
X.min)
Next, we consider a scenario where the observation do- main is expressed as union of finite number of hyper- rectangular regions:
i=1
(18)
where sine(x) = sin(x)/ //x is the unnormalized sinc func- tion, and w = (w') w"" gests that the equivalent kernel h(x, x²) has degenerate form of rank 2.M, which is obtained through Cholesky de- composition as h(x,, x' = (LΦ(x)) (Lo(x')), ([Lo()))((L') where L L = L2 + A) -1 To enhance the approxima- tion accuracy of the random Fourier features, we employed the quasi-Monte Carlo feature maps (Yang et al., 2014) in this paper.
x= U X X,
II[ [ min X.ma
(15)
The relation (17) sug-
where J denotes the number of hyper-rectangular regions, and X, n2y+z = Ø. While prior studies typically assume single hyper-rectangular region, the assumption (15) en- ables us to deal with more complicated observation do- mains, such as disjoint or irregularly shaped regions, often encountered in practical applications.
The degenerate form of equivalent kernel (17) offers an additional advantage. For cross-validation with the least squares loss, KºIE needs to evaluate the integral of the
The Fredholm integral equation (8) generally cannot be solved in closed form, and Flaxman et al. (2017) proposed
5

---

Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
Gaussian kernel
Laplace kernel
Cauchy kemel
((x%))) g(xx%, (v(0)) g(4
1
-2
Figure 1: Examples of smoothing kernels (2,,2.) smoothing kemels with edge-correction (x),z))/v(x) and equivalent kernels 112 Inn) for data points 2n E {-1.8, 0, 1.8). Gaussian, Laplace, and Cauchy RKHS/smoothing kernels are e-1a=aal" st- and 1+|r- respectively. The regularization hyper-parameter, the number of random features, and the observation domain were set as (r, 2M, X) = (2, 500, 22,
squared intensity function, (x, h(x, xx)) )x)) which requires (N2) computation naively or (NN)) computa- tion with ((> 1) points Monte Carlo integration. But the integral of the squared intensity function be obtained analytically with O(M2 + MN) computation under (17)
alent RKHS kernels in K2IE, showing that both the edge- corrected smoothing kernels and the equivalent RKHS ker- nels assign greater weight data points near the boundary of the observation domain (x) ~ 2) compared to those at the center. Interestingly. K applies edge correction more conservatively through the equivalent RKHS kemnlls h(r, Fn), that is, differentiates the weights between the center and the boundary less significantly compared to KIE with g(æ, 2x, 2x)/v(z))
as follows:
1 da: Σ h(D, ain))
= AE
(19)
< = -7 I2 + A)
4. Experiments
Therefore, regarding hyperparameter tuning, KIE and FIE require MC integration and solving a dual optimization problem for each cross-validation, respectively, whereas K²IE requires neither, which beneficial especially in multi-dimensional settings.
We evaluated the validity and the potential efficiency of our proposed KCIE by comparing it with prior nonparametric approaches, including the kernel intensity estimator with edge correction (KIE) (Diggle, 1985) and Flaxman's ker- nel method-based intensity estimator (FIE) (Flaxman al., 2017), using synthetic datasets. For K$IE and FIE, the number of random features 2M was fixed at 500 (see Ap- pendix D for an ablation study on the feature number 2.M).
The comparison of II with KIE suggests that from the viewpoint of KIE, the primary distinction between them lies in how smoothing kernels with edge-correction are constructed: In KIE, the smoothing kernels with edge- correction are constructed by rescaling density functions with their integrals over the observation domain; In con- trast, Ka constructs smoothing kernels the solution to the integral equation (8). During model training. KIE ben- efits from more computational efficiency than K-IE, which requires solving the integral equation. However, K of- fers computational advantages during inference as can perform the intensity function integration needed in pre- dictive tasks (e.g., see Equation (2)) analytically, while KIE relies on Monte Carlo integration. Furthermore, as demon- strated by (Flaxman et al., 2017), the smoothing kernel in II i.e., the equivalent kernel is expected to achieve more effective edge correction, particularly in high-dimensional domain settings.
For both the smoothing and RKHS kernels, we em- ployed a multiplicative Gaussian function, 2(x, x') = e [00(z==)) where B (31,..., Ba)T is the inverse scale hyper-parameter, and o denotes the Hadamard prod- uct. KIE optimized the cross-validation based on the negative log-likeldhood func- tion: FIE optimized the syper-parameters. (3,7), 7), using the same cross-validation procedure KIE; For K2TE, the ppper-parameters, (B, 7), were optimized via 5-fold cross- validation with the least squares loss function (10). For all models, the Monte Carlo cross-validation with p- p-thinning (Cronie et al., 2024) was adopted, where p was fixed at 0.6. A 10 x 10 logarithmic grid search was conducted for 7 [0.1, 100] and B3 [0.1, 1000-00- -B, where = -B,...., ...). for 3, 1/[max, (xnxx,(X.j0)) (X (Xxa min, (Xiii)) For FIE, the gradient descent algorithm Adam (Kingma & Ba, 2014) was employed to solve the dual optimization problem (9).
/ through 5-fold
Figure/ 1 illustrates examples of smoothing kernels with and without edge correction KIE, as well the equiv-
6

---

Kernel Method-basee based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
3.5
FIE
KIE
AiD (X)
AIo (X) I
A
(X)
0.0
111
0.0
11I
1000
10
20
30
20
40
60
80
100
Figure 2: Examples of the estimated intensity functions ID synthetic data. The vertical lines represent the locations of
observed events.
Results on ID synthetic data across 100 trials with standard errors in brackets. N denotes the average data size trial. The performances not significantly (p > 10-2) different from the best one under the Mann-V Mann-Whitney Whitney U test (Holm, 1979) are shown in bold for 2 and LZ].
Table 1:
(x)(x)
N =6
(x)
Ñ = 33
(z) N
= 226
L
pt
cpu
L
[1|
et
cpu
L +
et
cpu
0.09
0.23
12.6
2.97
0.15
0.30
KIB
(0.07)
(0.08)
(2.82)
(0.28)
(0.07)
(0.08)
0.11
0.24
0.34
1.06
13.2
3.04
0.46
0.29
0.17
0.33
0.33
0.86
FIE
(0.11)
(0.09)
(0.17)
(3.41)
(0.28)
(0.30)
(0.09)
(0.09)
(0.39)
0.12
0.26
0.26
0.01
13.9
3.09
0.48
0.01
0.18
0.34
0.31
0.01
KCIE
(0.08)
(0.08)
(0.00)
(5.03)
(0.45)
(0.00)
(0.08)
(0.09)
(0.00)
10x (Z)
N
466
10x (2):
N
328
10 x (x))
Ñ
2250
I.2
|
pt
cpu.4
I.2 +
[L]4
p
cpu
I!? 4
et
cpu
1.43
0.87
289
13.5
2.84
1.29
KIF
(1.03)
(0.29)
(71.3)
(1.92)
(1.68)
(0.34)
1.74
0.93
0.49
1.77
277
13.0
0.64
0.55
2.70
1.25
0.63
0.61
FIE
(1.53)
(0.39)
(0.13)
(80.6)
(2.09)
(0.33)
(1.79)
(0.37)
(0.13)
1.67
0.92
0.49
0.01
266
12.7
0.77
0.01
3.24
1.34
0.47
0.01
K-IE
(0.71)
(0.36)
(0.00)
(74.6)
(1.98)
(0.00)
(2.08)
(0.41)
(0.00)
Predictive performance assessed using the integrated squared error (L and the integrated absolute error (L) (Kowalczuk & Kozlowski, 1998), defined as follows:
we generated 1D datasets based three types of intensity functions:
\D (x) 2e x/11 + e (00--25511
X' = [0, 50),
2 = 1 / (^*(x) (((x)-((ax)) (A'(x)-((z)) dx,
( = 5 (ai) + 6,
X' = [o, 5],
X x
(20)
x (x) (x) = piecewise linear function,
X' = [0, 100),
X / (x) - {xx))
14
(21)
where (x: and \(x) denote the true and estimated in- tensity functions, respectively. Following (Flaxman et al., 2017), the fraction of times that T,2 is smaller than KIE across the trials, denoted by p. was also reported, where P not defined for KIE Efficiency on the CPU time (in seconds), epu, required execute the model fitting given the optimized
where (1) (z) passes through the points: (0,2), (25, 3), (50, 1), (75, 2.5). and (100,3). Furthermore, evaluate the scalability of K and FIE with respect data size, generated ID datasets using intensity functions scaled by factor of ten, denoted by 10 X Aip (z) for 21,2,3). For each intensity function, we simulated 100 trial sequences and performed intensity estimation 100 times using the compared methods.
evaluated based
All models were implemented using 10 and exccuted on a MacBook Pro equipped with 12-core CPU (Apple M2 Max), with the GPU disabled.
Table displays the predictive performance on the ID syn- thetic datasets. It shows that proposal Kº matched the predictive performances of FIE across all three datasets, while achieving significantly faster model fitting in terms of CPU time. K² achieved comparable predictive perfor- mances with KIE on (x).. ) but was outperformed by KIE on (x) and Aip (z)). This result consistent with (Flax-
4.1. 1D Synthetic Data
In accordance with previous studies (Adams et al., 2009; John & Hensman, 2018: Aglietti et al., 2019; Kim, 2021),
7

---

Kernel Method-bascy Kernel Intensity Estimators for Inhomogeneous Poisson Processes
True
KIE
FIE
K"EE
3
3
2
2
2 3
5
2
4
5
5
4-
3-
3-
1
1
0
0
2
3
4
Figure 3:
Examples the estimated intensity functions on 2D synthetic data AnD and A0.8 locations of observed events, and the unobserved regions are indicated by hatched lines.
The black dots represent the
Results on 2D synthetic data across 100 trials with standard errors in brackets. N denotes the average data size per trial. The notation follows Table 1.
Table 2:
A (x)
554
443 A (ar)
Ñ
U.S ADD
(x)
N
128
2
1
pt
cpu
L
11|
et
cpu
L?
pt
cpu
63.3
6.36
63.5
6.35
64.5
6.34
KIF
(8.96)
(0.40)
(8.92)
(0.45)
(10.9)
(0.52)
56.47
5.38
0.80
1.54
59.8
5,53
0.71
1.45
62.3
5.64
0.64
1.50
FIE
(12.2)
(0.60)
(0.34)
(13.4)
(0.62)
(0.33)
(13.5)
(0.65)
(0.34)
53.0
5.54
0.97
0.16
55.1
5,63
0.90
0.14
57.9
5.77
0.85
0.13
K-IE
(10.2)
(0.49)
(0.03)
(11.1)
(0.51)
(0.02)
(12.2)
(0.55)
(0.03)
man et al., 2017), demonstrating that KIE performed very well fow-dimensional settings. is worth noting that the discrepancy of predictive performances between KIE and II became negligible 10 x (x)), where the dataset size increases. Also, Table demonstrates that the CPU time of FIE increases with the data size, while that of I2 remains nearly constant. Figure 2 displays some estimation results.
using the compared methods. The resulting dataset con- tained approximately 540 data points per trial.
In this study, considered scenario where the observa- tion domain was divided into 5x = 25 sub-domains,
X U
Xy, X2 evenly partitioned 2D domain, (22)
4.2. 2D Synthetic Data
with some of the sub domains being missing. For each trial of the dataset, we randomly selected each sub-domain with probability p E {1.0, 0.9, 0.8}. thereby generating three datasets, denoted as x x (xr). and X (2:), respec- tively.
Following the procedure in (Lloyd et al., 2015), we gen- crated a 2D dataset from sigmoidal Gaussian Cox pro- cess. Specifically,, we first sampled 2D function from Gaussian process with an RKHS kernel, k(x,,²) x') 11==-=11//2 over the domain X' = [0, 5] x [0, 5]. The inten- sity function was then obtained by applying sigmoid link function, (2) 50/(1 + 2002), the sampled fune- tion. Using the intensity function, simulated 100 trials of event data and conducted intensity estimation 100 times
Table displays the predictive performance the 2D syn- thetic datasets, which shows that k and FIE consistently outperformed KIE in all datasets. This result suggests that KºIE and FIE could more effectively handle edge effects than KIE in multi-dimmesionall settings. Notably, regard- ing the integrated squared error I,2, KaIE achieved superior
8

---

Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
predictive performance, on average, than FIE, despite both methods employing the same equivalent kernels (with hy- perparameters optimized individually for each model). It might be due to the fact that KIE is based on the minimiza- tion of the least squares loss (see Section 3.1). Another possible explanation for this result is that the optimization of the dual coefficient required in FIE may become unsta- ble. Indeed, John & Hensman (2018) reported that FIE can yield unreasonable solutions for highly modulating inten- sity functions. In contrast, KC is expected to work more robustly, as it does not require the optimization of dual co- efficients. Figure 2 displays some estimation results on AD (z) and X (z).
ing condition that the optimal function (ax) must satisfy:
1 f(x) + ((s)))))(((s)))((((()
= Σ k(a, (n)d'(((xn)) (f(zn)),
where a' o'(v) do (v). When (v) = y. the above equation reduces to a Fredholm integral equation, for which Theo- rem provides a tractable solution. However, if a is non- linear, even for simple cases like (ty) = y², the resulting nonlinear integral equation becomes analytically and nu- merically challenging to solve.
An alternative approach is to impose non-negativity con- straints at finite set of virtual points, which leads to dual optimization problem. Although this approach may reduce the risk of negative estimates at the virtual points, it neither guarantees global non- non-negatvviyy nor preserves the compu- tational simplicity of K2IE, due to the added complexity introduced by the dual optimization.
Additional experiments with a scalable Bayesian approach and on a real-world dataset are provided in Appendix E.
5. Diseussions
We have proposed a novel penalized least squares loss for- mulation for estimating intensity functions that resides in an RKHS. Through the path integral representation of the squared Hilbert space norm, we showed that the optimiza- tion problem encompasses a representer theorem, and de- rived feasible intensity estimator based on kernel meth- ods. We evaluated the proposed estimator on synthetic data, confirming that it achieved comparable predictive ac- curacy while being substantially faster than the state-of- the-art kernel method-based estimator.
Does 1I truly fail to guarantee non-negativity in its orig- inal form? Interestingly. a sufficient condition to ensure the non-negativyty of the equivalent kernels arising in Flax- man's 's estimator (and, of course, in K"IE) has been estab- lished by Kim (2024). Specifically, when RKHS kernels belong to the class of inverse M-kernels (IMKs), the as- sociated equivalent kernels h(x, x²) are guaranteed to be non-negative. This suggests that K2IF,, as defined by a sum of equivalent kernels in Equation (12), may be a non- negative intensity estimator whenever the RKHS kernel is an IMK. In one-dimensional cases, the Laplace kernel is known to be an IMK, but no general construction of IMKs is currently known in higher dimensions-posing an intrigu- ing open problem.
LIMITATIONS AND FUTURE WORK
As noted at the end of Section 3.1, key limitation of our KI lies in its lack of a general guarantee for the non-negativity of the resulting intensity function. To in- vestigate the effect, we conducted an analysis of how fre- quently IE produces negative values using the 2D syn- thetic dataset ^1o Specifically, we evaluated the estimated intensity values at 500 x 500 grid points within the obser- vation domain and computed the ratio of negative values. The mean + standard deviation of this ratio across 100 tri- als was 0.059 + 0.016, indicating that K2 can indeed pro- duce negative estimates in practice particularly in regions with sparse data-highlighting the necessity of a post-hoc correction such as clipping via max(A (Axx), 0) in applica- tions where negative intensity values are not permitted. As a direction for future work, we explore the technical chal- lenges involved in incorporating non-negativity constraints directly into the functional optimization problem defined in Equation (1)).
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential societal consequences of our work, none of which we feel must be specifically highlighted here.
References
Adams, R. P., Murray, I., and MacKay, D. J. Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities. In International Con- ference on Machine Learning, pp. 9-16. 2009.
Aglietti, V., Bonilla, E. V., Damoulas, T., and Cripps, S. Structured variational inference in continuous Cox pro- cess models. In Advances in Neural Information Pro- cessing Systems 32, 2019.
One natural approach is to model the intensity function as a non-negative transformation of (a=) of a latent function f(x) residing in an RKHS. In this setting, the functional analysis of the objective in Equation yields the follow-
Atkinson, K. personal perspective on the history of the

---

Kernel Method-basel Kernel Intensity Estimators for Inhomogeneous Poisson Processes
numerical analysis of Fredholm integral equations of the second kind. In The Birth of Numerical Analysis, pp. 53-72. World Scientific, 2010.
Gatrell, A. C., Bailey, T. C., Diggle, P. J., and Rowlingson, B. S. Spatial point pattern analysis and its application in geographical epidemiology. Transactions of the Institute of British Geographers, pp. 256-274, 1996.
Bacry, E., Bompaire, M., Gaiffas, S., and Muzy, J.-F. Sparse and low-rank multivariate Hawkes processes. Journal of Machine Learning Research, 21(50):1-32, 2020.
Gunter, T., Lloyd, C., Osborne, M. A., and Roberts, S. J. Efficient Bayesian nonparametric modelling of struc- tured point processes. In Uncertainty in Artificial In- telligence, 2014.
Cai, B., Zhang, J., and Guan, Y. Latent network structure learning from nigh-dimensional multivariate point pro- cesses. Journal of the American Statistical Association, 119(545):55108,, 2024.
Hansen, N. R., Reynaud- Leynaud-Bouret, P., and Rivoirard, V. Lasso and probabilistic inequalities for multivariate point pro- cesses. Bernoulli, 21(1):83-143 2015.
Clark, T. G., Bradburn, M. J., Love, S. B., and Altman, D. G. Survival analysis part i: basic concepts and first analyses. British Journal of Cancer, (92)2232-2238 2003.
Heikkinen, J. and Arjas, E. Modeling a Poisson forest in variable elevations: A nonparametric Bayesian ap- proach. Biometricx, 55(3):738-745, 1999.
Holm, S. A simple sequentially rejective multiple test pro- cedure. Scandinavian Journal of Statistics, pp. 65-70, 1979.
Cox, D. R. Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methalological); 34(2):187202, 1972.
Hubbel, P. and Foster, R. B. Diversity of canopy trees in a neotropical forest and implications for conservation. Tropical Rain Forest: Ecology and Management, pp. 25- 41, 1983.
Cronie, O., Moradi, M., and Biscio, C. A. cross- validation-baet based statistical theory for point processes. Biometrika, 111 (2):625-641, 2024.
Cunningham, J. P., Byron, M. Y., Shenoy, V., and Sa- hani, M. Inferring neural firing rates from spike trains using Gaussian processes. In Advances in Neural Infor- mation Processing Systems 20, 2007.
John, S. T. and Hensman, J. Large-scale Cox process infer- ence using variational Fourier features. In International Conference on Machine Learning, volume 80, pp. 2362- 2370. PMLR, 2018.
Daley, D. and Vere-Jones, D. An Introduction the Theory of Point Processes. Springer- Verlag, New York, 1988.
Jones, M. C. Simple boundary correction for kernel density estimation. Statistics and Computing, 3:135-146, 1993.
Kim, H. Fast Bayesian inference for Gaussian Cox pro- cesses via path integral formulation. In Advances in Neu- rul Information Processing Systems 34, 2021.
Davis, R. A., Lii, K.-S., and Politis, D. N. Remarks some nonparametric estimates of density function. Se- lected Works of Murray Rosenblatt, pp. 95-100, 2011.
Kim, H. Inverse M-kernels for linear universal approxima- tors of non-negative functions. In Advances in Neural Information Processing Systems 37, 2024.
Diggle, P. A kernel method for smoothing point process data. Journal of the Royal Statistical Society: Series C (Applied Statistics), 34(2):138-1477 1985.
Kim, H., Asami, T., and Toda, H. Fast Bayesian estimation of point process intensity as function of covariates. In Advances in Neural Information Processing Systems 35, 2022.
Diggle, P. J., Moraga, P., Rowlingson, B., and Taylor, log- Gaussian Cox processes: extending the geostatistical paradigm. Sta- tistical Science, 28(4):542-563 2013.
B. M. Spatial and
Kingma, D. P. and Ba, J. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6990 2014.
Donner, C. and Opper, M. Efficient Bayesian inference of sigmoidal Gaussian Cox processes. Journal of Machine Learning Research, 19:1-34, 2018.
Kowalezuk, Z. and Kozlowski, J. Integrated squared er- ror and integrated absolute error in recursive identifica- tion of continuous-tine plants. In UKACC International Conference on Control '98 (Conf. Publ. No. 455), vol- ume 1, pp. 693 698. IET, 1998.
Flaxman, S., Teh, W., and Sejdinovic, D. Poisson in- tensity estimation with reproducing kernels. Artificial Intelligence and Statisties, pp. 270-279. PMLR, 2017.
10

---

Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
Lai, C. D. and Xie, M. Stochastic Ageing and Depen- dence for Reliability. Springer Science & Business Me- dia, 2006.
Sellier, and Dellaportas, Sparse spectral Bayesian per- manental process with generalized kernel. In Interna- tional Conference Artificial Intelligence and Statis- tics, pp. 2769-2791 PMLR, 2023.
Lánezky, A. and Györffy, B. Web-based survival analy- sis tool tailored for medical research (kmplot): develop- ment and implementation. Journal of Medical Internet Research, 23(7)=e27333, 2021.
Shawe- Shawe-Taylor, and Cristianini, N. Kernel Methods for Pattern Analysis. Cambridge University Press, 2004.
Silverman, B. W. Density Estimation for Statistics and Data Analysis. Routledge, 2018.
Lloyd, C., Gunter, T., Osborne, M., and Roberts, S. Vari- ational inference for Gaussian process modulated Pois- son processes. In International Conference on Machine Learning. volume 37. pp. 1814-1822. PMLR, 2015.
Teng, M., Nathoo, F., and Johnson, T. D. Bayesian com- putation for log- oog-Gausiian Cox processes: A comparative analysis of methods, Journal of Statistical Computation and Simulation, 87:2227-222 2017.
Mercer, J. Xvi. functions of positive and negative type, and their connection the theory of integral equations. Philo- sophicul Transuctions of the Royal Society of London. Series A, containing papers of mathematical or physi- cal character, 209(44)-458):4154444 1909.
Tsuchida, R., Ong, C. S., and Sejdinovic, D. Exact, fast and expressive Poisson point processes via squared neu- ral families. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, pp. 90559-20566 2024.
Møller, J., Syversveen, A. R., and Waagepetersen, R. P. Log Gaussian Cox processes. Scandinavian Journal of Statistics, 25(3):451-435 1998.
van de Geer, S. Empirical Processes in M-estimation, vol- ume 6. Cambridge University Press, 2000.
Ogata, Y. Statistical models for earthquake occurrences and residual analysis for point processes. Journal of the American Statistical Association, 83(401):9-27, 1988.
Wahba, G. Spline Models for Observational Data, vol- ume 59. SIAM, 1990.
Walder, C. J. and Bishop. A. N. Fast Bayesian intensity estimation for the permanental process. In International Conference on Machine Learning, volume 70, pp. 3579- 3588. PMLR, 2017.
Parzen, E. On estimation of probability density function and mode. The Annals of Mathematical Statistics, 33(3): 1065-1076, 1962.
Williams, C. and Seeger, M. Using the Nyström method to speed up kernel machines. In Advances in Neural Infor- mation Processing Systems 13, 2000.
Polyanin, A. D. and Manzhirov, A. Handbbook of Integral Equations. CRC press, 1998.
Rahimi, A. and Recht, B. Random features for large-scale kernel machines. In Advances in Neural Information Processing Systems 20, 2007.
Yang, J., Sindhwani, V., Avron, H., and Mahoney, M. Quasi-Monte Carlo feature maps for shift-invariant ker- nels. International Conference on Machine Learning. pp. 485 493, PMLR, 2014.
Ramlau-Hansen, H. Smoothing counting process intensi- ties by means of kernel functions. The Annals of Statis- tics, pp. 453-466, 466, 1983.
A. Explanation of the Least Squares Loss
Rathbun, S. L. and Cressic, N. Asymptotic properties of estimators for the parameters of spatial inhomogencous Poisson point processes. Advances in Applied Probabil- ity, (6611):222-155 1994.
Let E denote the expectation with respect to data points generated from the true intensity function X X*()). We con- sider the expected integrated squared loss between the es- timator X(x) and the true intensity function X \*(x), defined as:
Scholkopf, B. and Smola, A. J. Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond. MIT press, 2018
E ll/
EE 11/ i? {² ddæ
[}(x)) X2 (x)))
x
2E [/ {(x) (x)) + E 1 X (x))
Schölkopf, B., Herbrich, R., and Smola, A. J. A gen- eralized representer theorem. In International Confer- ence on Computational Learning Theory, pp. 416-426. Springer, 2001.
The third term on the right-hand side independent the estimator and can, therefore, be omitted. The second term

---

Kernel Method- bbsed Kernel Intensity Estimators for Inhomogeneous Poisson Processes
can be decomposed as follows:
where the orthogonality condition, Sx em (ar) (m))m (amm(a)))) (a))da 8mm's is used. The above equation shows that if we define new RKHS kernel q(:,)) as
NV.
21 / {(x))X(a))da)) = 2E
1/x )(x)
Σ 8(x
2n dae
==1
N
g x')
7m
Cm (x) e emn);
+ 2F
/ ((x) (1
Σ
8(x an) doy)
7m
+ 1/77
the optimization problem in Equation (11) reduces to:
where the second term on the right-hand side vanishes due to Campbell's theorem (Daley & Vere-Jones, 1988):
N
-2 HEHz
Σ Max)) + -IAII
min
N
E (xx) X (x))x
Σ
E [ ==1
where 11-112, represents the squared norm of an RKHS Hi associated with g(:s-). By construction, ((:s)) coincides with the equivalent kernel defined in Equation (7). Accord- ing to the classical representer theorem (Schölkopf et al., 2001), the optimal solution to this problem lies in the span of kernel evaluations at the data points:
= =fx e {[[{<)) (())d
/ E [((xa))((0)))) =0 0.
Combining the above expressions yields the following identity:
E x ((x)) - xxxx))
N
{(x) = Σ anq(a, Xnn,
= E [ (x "x(xx))) +C,
where the dual coefficients a = (o1....,aN' minimize the objective. Taking the gradient of the objective with re- spect to a yields:
where C' is a constant. This shows that the least squares loss defined in (10) corresponds the empirical integrated squared loss.
0
2
X(x,.. + 11A1
B. Proof of Theorem 1 via Mercer's Theorem
N
N
We present proof of Theorem I based on Mercer's ' Theo- rem, following an approach similar to that of Flaxman et al. (2017).
2
Σ
(I xnn) + +20p
( = 0,
n'=1
T2
On =1. 1.
This completes the proof.
Proof. Using the Mercer expansion of the RKHS kernel given in Equation (7), any function À E Hg can be ex- pressed as 2(.) = Σm bmem (c), where {bm) are the expansion coefficients and the RKHS norm is given by | -m ///./ /'/m < co. Substituting this into the
C. Extension of Theorem 1
Proposition 2. Let k IRd x I2 + lR be continuous pos- itive defiidefinite kernel. Then the solution X(-) to the op- timization problem (11) admits representer theorem with respect to a transformed RKHS kernel h((,-), which is de- fined via the following Fredholm integral equation:
objective in Equation (11), we obtain:
1 A(xx) +
11
+ f dxx
1 h(2,2') J k(z, x x')ds k(xx, x'),
N
-2)
A(xn) + - Σ
ben/ 'm
7
(CI)
(a, x E Rd x Bd
n=1
m
bm bime [
em (a)()))m (x) dx
Moreover, its dual coefficient is equal to unity
N
N
= 2
A(x..)) + Σ
12n Um + Σ 12
/m
i(x) = Σ h(x In),
x E Rd
n=1
m
m
N
= 2 Σ A(x,,) + [(
Tmm
62
Proof. Let the integral operator associated with the RKHS kernel .(.,.) be defined as (C-(x) = Sge k(x,, 8 8)ds, and its
1m + 1/7
m=
12

---

Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
inverse operator be denoted by KC* -(xr) = SRu k" (ar, s)ds. Using the path integral formulation of Gaussian processes (Kim, 2021), the squared norm in the RKHS can be ex- pressed in the functional form:
Table DI: Predictive performance of I on the 2D syn- thetic data X2 as a function of the number of feature maps. Brackets represent standard errors over 100 trials.
2M
20
100
300
500
11 - (x, (s)
|
147
75.4
53.2
53.0
L
(8.28)
(10.4)
(10.7)
(10.2)
Based on this representation, the objective functional in Equation (11) becomes:
[1]
9.807
6.681
5.56
5.54
(0.26)
(0.45)
(0.52)
(0.49)
N S(A) 2 Σ Maxn) JL (a (s)
shown in Table DI, both the integrated squared error and the integrated absolute error consistently decrease as M increases. These results indicate that IC benefits from more random features, and that the setting 2M = 500, used in Section 4, provides sufficiently accurate and stable esti- mates.
where ht" (:--) is defined in terms of k- (:--,, the Dirac delta function (-)) and the indicator function (c))
h" (ax 8) = ō(œ s)loEx + 1 k² kxxx,, 8).
(C2)
7
E. Additional Experiments
The minimizer X(x) of (A) satisfies the equation obtained by setting the functional derivative to zero:
E.1. Comparison with Variational Bayesian model
We conducted an additional experiment the 2D synthetic dataset AnD to compare against scalable Bayesian model. Here, we adopted variational Bayesian approach based on a Gaussian Cox process with a quadratic link function (Lloyd et al., 2015), where a Gaussian RKHS kernel and 10 10 inducing points were employed. We employed a gradient descent algorithm, Adam (Kingma & Ba, 2014), to perform the model fitting, where the number of itera- tions and the learning parameter were set 5000 and 0.01, respectively. L2, [L]. and cpu achieved by the Bayesian model were 63.9 (12.2), 5.55 (0.46), and 51.8 (32.2), re- spectively, where standard deviations are in brackets. The result highlights the high efficiency of K2
S5
N
5(x xn)) + yd (a s){(s) O. Σ
Define the integral operator corresponding to (;)) by H' (x) Sea h" (", s)ds, and its inverse operator by H (ar) = fga h(x, s)ds. Applying HH to both sides of the functional equation yields:
\(x)
Σ M(x, mn),
E Rd.
where we have used the identity, (HH** SR -((x- s)ds. Furthermore, applying the operator H to Equa- tion (C2) leads to the relation,
E.2. Comparison on a Real-world Dataset
1 / k" 8)h(s o' d's,
We conducted an additional experiment using an open 2D real-world dataset, bei, in the R package spatsta (GPL- 3). It consists of locations of 3605 trees of the species Beilschmiedia pendula in tropical rain forest (Hubbel & Foster, 1983).
8(x - x') h(r, 1') +
and subsequent application of the operator IC results in
1 M( a'), k(0, = / k(æ, 8)h(s, x' )ds +
Following (Cronie et al., 2024), we randomly labeled the data points with independent and identically distributed marks (,,2,3) from a multinomial distribution with pa- rameters (p1,P2, P2, P2,P)) = (0.3, 0.3, 0.7), and assigned the points with label and to training data and test data, respectively;
7
(x, x') E Rd x Rd,
which is identical to Equation (C1).
repeated it 100 times for evaluation. A 10 x 10 logarithmic grid search was conducted for [0.001, 1) and B E [0.1, 100) B3, where 3 = (B 3a)
D. Ablation Study on the Number of Random Features
for 3, = 1/ [max, (XIIX)) min, min,(X(rim)) (X:rji))
We conducted an ablation study to investigate the effect of the number of random features (2M) on the predictive per- formance of K2IF using the 2D synthetic dataset 21 As
Let the observation domain X be regularly devided into 10 x 10 sub-domains as X' - Xz. We evaluated the pre-
13

---

Kernel Method-based Kernel Intensity Estimators for Inhomogeneous Poisson Processes
Table El: Results on the real-world data bei across 100 trials with standard errors in brackets. The notation follows Table
L, +
Let
quu
-5.80
267
KIE
(0.32)
(11.5)
-5.16
287
5.15
FIE
(0.26)
(15.1)
(1.57)
KºIE
-6.16
279
0.17
(0.44)
(13.2)
(0.04)
dictive performance of the estimator X(x) based the test least squares loss (In)) and the test negative logliikl of counts (Lo):
La / (x) 2 Σ {(xa),
DEDi
L. Σ - Nogg do A,= / {(x) dx, X.
where Diss denotes the data, and N, represents the number of test data points observed within Xy. Table El displays the results, showing that K'IE achieved the best performance on L. but was outperformed by KIE on Le. which could be because the hyperparameters were opti- mized based on the least squares loss and the logllikelihood for II and KIE, respectively.
14


Reviews:
[{"review_id": "MQhZVkmEK9", "reviewer": "Reviewer_upvZ", "summary": "The paper develops a new kernel based estimator for intensity in inhomogeneous Poisson processes. The estimator is shown to be a associated with a unique reproducing kernel Hilbert space and is compared to some previous estimation methods in a simulation study. The simulation study shows that the new methods achieves better reconstructions at a lower computational cost compared to existing methods.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "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": "We thank the reviewer for the deep understanding of our model and the constructive comments. We provide a detailed response to each comment.\n\n**Why not include a lambda>=0 constraint?**, **My feeling is that the paper presents a theoretically sound solution to the wrong problem.**\n\nAs pointed out, one can enforce the non-negativity of the estimator by introducing lambda>=0 constraint or by applying a non-negative transformation. However, as far as we know, such constraints prevent us from obtaining efficient estimators like K$^2$IE. The main contribution of our work lies in showing that, by sacrificing strict non-negativity, one can obtain a feasible kernel-based estimator comparable to classical KIEs. \n\nAs discussed in Section 3.1, non-negativity can be enforced in a post-hoc clipping like max($\\hat{\\lambda}(x)$,0). However, we totally agree with the reviewer's point that the functional optimization problem (11) does not explicitly take non-negativity into account, hence not yielding an optimal \"non-negative\" estimator. For future work, we would like to discuss technical issues arising when non-negativity constraints are imposed. \n\nConsider modeling the intensity as a non-negative transformation $\\sigma(x)$ of a latent function $f(x)$ lying in an RKHS. Then, the functional derivative of the objective functional in Eq. (11) leads to the following equation that the optimal $\\hat{f}(x)$ solves (full derivation is omitted): $$\\frac{1}{\\gamma} f(x) + \\int_{X} k(x,s)\\sigma(f(s)) \\sigma'(f(s)) ds = \\sum_n k(x,x_n) \\sigma'(f(x_n)),$$ where $\\sigma'(x)$ is the derivative of $\\sigma(x)$. When $\\sigma(x) = x$, the equation reduces to a Fredholm integral equation for which Theorem 1 provides a feasible solution. However, when $\\sigma(x)$ is nonlinear, even as simple as $\\sigma(x) = x^2$, deriving a feasible solution becomes non-trivial.\nAlternatively, one may consider enforcing non-negativity of intensity at finite virtual points {$q_1,\\cdots, q_R$}, which leads to a dual optimization problem. While this approach may reduce the risk of negative estimates at $q_i$, it does not guarantee non-negativity of intensity elsewhere and undermines the computational advantages inherent in K$^2$IE due to the added complexity of dual optimization.\n\nFinally, we discuss whether the problem in Eq. (11) is truly improper/wrong. (Kim, NeurIPS2024) showed that when RKHS kernels belong to the class of inverse M-kernels (IMKs), the corresponding equivalent kernels $h(x,x')$ are non-negative. This suggests that Eq. (11) may not be inherently improper because K$^2$IE, a sum of equivalent kernels, is also non-negative under IMKs. In one-dimensional cases, the Laplace kernel is known to be an IMK, but no construction is known for IMKs in higher dimensions, posing an interesting challenge. We will include the discussion in the text.\n\n**Can the method be used to reconstruction of intensity at unobserved regions**\n\nYes, K$^2$IE can reconstruct intensity at unobserved regions. At submission, we had assumed the intensity estimation at observed regions. However, in light of your insightful comment, we revisited the model and confirmed that K$^2$IE in Eq. (12) could accept inputs from unobserved regions as well. Accordingly, a minor revision to Theorem 1 is warranted to reflect this more general setting. Specifically, all integral operators and the squared RKHS norm should be defined over the full domain ($X$ -> $\\mathbb{R}^d$), and Eq. (13) should be updated as $q^*(x,s) = \\delta(x-s) 1(s \\in X) +\\frac{1}{\\gamma}k^*(x,s)$, where $1(\\cdot)$ is the indicator. With this modification, Theorem 1, initially defined in $x \\in X$, should be revised to hold for $x \\in \\mathbb{R}^d$. We sincerely appreciate the valuable comment. Due to time constraints, we have not yet conducted experiments to evaluate accuracy in unobserved regions. But, we are willing to include the results in the final version of the paper if you consider it essential.\n\n**How does the method compare to log-gaussian-cox process methods?**\n\nIn response to several reviewers' suggestions, we conducted an additional experiment using the 2D synthetic dataset $\\lambda_{\\text{2D}}^{1.0}$ to include the result of a scalable Bayesian method​. This time, we adopted a variational Bayesian model with a quadratic link function (Lloyd, ICML2015), where 10 x 10 inducing points were employed. We appreciate your kind suggestion of references (we will cite them), but we have not yet become proficient with the R package. $L^2$,  $|L|$, and cpu of the Bayesian method were 58.0(8.13), 5.61(0.47), and 28.4(0.96), respectively, where standard deviations are in brackets. The result highlights the high efficiency of K$^2$IE. Note that the reported results are based on 10 trials. We will include the results based on 100 trials in the text.  \n\n**Does the methods provide uncertainty estimates for the reconstruction**\n\nK$^2$IE does not provide uncertainty estimates, which limits our model compared to Bayesian models.", "claims_and_evidence": "The theoretical derivation seems sound, I have some issues with the simulations study that does not include comparison to any state of the art statistical methods (e.g. https://doi.org/10.1111/2041-210X.13168) and which seems to solve the wrong problem for partially observed processes.\n\nFor a partially observed process we would like reconstructions for the entire domain even when observations are only from part of the domain. For this paper that would entail in equation (11) that the integral over lambda(x) (second term) be over the observed area while the last integral/norm (i.e. the penalty term) be over the entire domain, thus having a penalty on difference between observations and lambda(x) only for observed parts of the domain but a smoothness penalty for the entire domain, including unobserved parts. For the lower row in figure 3 this would imply that the reconstruction of lambda is computed (and evaluated) also for the hash-marked regions.", "methods_and_evaluation_criteria": "The simulation study is reasonable but evaluation is only based on L1 and L2 error between reconstruction and actual intensity in the observed area. No attempt at checking prediction uncertainties or the models ability to predict intensities at unobserved locations is made (see comment under Claims And Evidence). No comparison to methods for log-gaussian-cox processes is made (see Essential References Not Discussed).", "theoretical_claims": "The derivation of theorem 1 seems sound. However, as the authors themself note the solution can have unreasonable results, i.e. negative intensities, this is due to an improper formulation of the minimisation problem in (11), either an additional constraint lambda>0 or a transform of lambda to ensure non-negative intensities should be included. My feeling is that the paper presents a theoretically sound solution to the wrong problem.", "experimental_designs_or_analyses": "See comments under Methods And Evaluation Criteria.", "relation_to_broader_scientific_literature": "The paper extends on existing RKHS theory for estimation of inhomogeneous Poisson processes. It provides some references to spatial statistical development but does not include any methods from spatial statistics in the comparison (see Essential References Not Discussed).", "essential_references_not_discussed": "The paper lacks references to some recent log-gaussian cox process literature and fast numerical methods for these. E.g. https://doi.org/10.1093/biomet/asv064 and https://doi.org/10.1111/2041-210X.13168. Especially the later could be included in the simulation study.", "other_strengths_and_weaknesses": "No additional comments", "other_comments_or_suggestions": "No additional comments", "questions_for_authors": "Key points to consider:\n1) Why not include a lambda>=0 constraint?\n2) Can the method be used to reconstruction of intensity at unobserved regions (e.g. the hash marked areas in Fig 3)\n3) How does the method compare to log-gaussian-cox process methods?\n4) Does the methods provide uncertainty estimates for the reconstruction, e.g. V(lambda(x)|observations), or is only the best reconstruction provided?", "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": "", "extra_scores": {}}, {"review_id": "b6reeDd6wH", "reviewer": "Reviewer_XSJm", "summary": "This paper introduces K2IE, a kernel method-based kernel intensity estimator for inhomogeneous Poisson processes, which formulates the intensity estimation as a penalized least squares loss minimization in RKHS. A key theoretical contribution is the establishment of a specialized representer theorem leading to a computationally efficient estimator with unit dual coefficients, drawing a formal connection between classical kernel intensity estimators (KIEs) and RKHS-based estimators. The method is validated on 1D and 2D synthetic datasets, demonstrating comparable predictive performance to prior methods while offering improved computational efficiency.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "We would like to thank the reviewer for the highly positive and constructive comments, by which we are strongly encouraged. Below, we provide a detailed response to each of the comments.\n\n**Why are real-world datasets not included in the evaluation? This would be critical to support claims of practical relevance.**\n\nWe focused our evaluations on synthetic datasets, which allow for precise error estimation between the true and estimated intensity functions—an appropriate setting to verify the theoretical soundness of our model. However, from the perspective of practical relevance, we fully agree with the reviewer on the importance of validating our approach on real-world datasets. In response to the reviewer’s suggestion, we have conducted an additional experiment using an open 2D real-world dataset, *bei*, in R package *spatsta* (GPL-3). It consists of locations of 3605 trees of the species *Beilschmiedia pendula* in a tropical rain forest (Hubbell & Foster, 1983). Following (Cronie et al., 2024), we randomly labeled the data points with independent and identically distributed marks {1, 2, 3} from a multinomial distribution with parameters $(p_1 = p_2 = 0.3, p_3 = 0.4)$, and assigned the points with label 1 and 2 to training data and test data, respectively; we repeated it ten times for evaluation. We evaluated the predictive performance of the estimators $\\hat{\\lambda}(x)$ based on the test least squares loss ($L_{s}$) and the negative test likelihood of counts ($L_{c}$): $L_{s}$ was computed as $L_{s} = \\int_{X} \\hat{\\lambda}^2(x) dx - 2\\sum_{n \\in D_{\\text{test}}} \\hat{\\lambda}(x_n)$, where $D_{\\text{test}}$ was the test data; the observation domain $X$ was discretized into10 x 10 sub-domain {$X_1, \\dots, X_{100}$}, and $L_{c}$ was computed as $L_{c} = \\sum_{i=1}^{100} (\\Lambda_{X_i}) - N_{X_i} \\log \\Lambda_{X_i} - \\log(N_{X_i}!)$, where $N_{X_i}$ is the number of test data points observed in $X_i$, and $\\Lambda_{X_i} = \\int_{X_i} \\hat{\\lambda}(x) dx$. We obtained the results across 10 trials with standard errors in brackets (the lower, the better): $L_s$ = -5.74(0.39), -6.17(0.52), -5.09(0.31) for KIE, K$^2$IE, FIE; $L_s$ = 265(14.1), 278(10.8), 287(19.7) for KIE, K$^2$IE, FIE. The results show that our K$^2$IE achieved the best performance on $L_s$, but was outperformed by KIE on $L_c$, which could be because the hyperparameters were optimized based on the least squares loss and the log-likelihood for K$^2$IE and KIE, respectively. We will include the results based on 100 trials in the final manuscript.\n\n**How sensitive is the model to the number of random Fourier features (M)? An ablation study would help demonstrate robustness.**\n\nWe thank the reviewer for the constructive feedback. In response to the suggestion, we conducted an ablation study on the number of random features ($2M$) using the 2D synthetic dataset $\\lambda_{\\text{2D}}^{1.0}$​. The integrated squared errors ($L^2$) of our K$^2$IE were 149(6.72), 74.8(12.1), 50.7(7.71), and 49.0(8.80) for $2M$ = $20$, $100$, $300$, and $500$, respectively, where standard deviations are in brackets. Similarly, the integrated absolute errors ($|L|$) were 9.86(0.25), 6.65(0.60), 5.43(0.52), and 5.31(0.57) for the same values of $M$. These results indicate that the predictive performance of K$^2$IE improves with larger $M$, and our chosen setting ($2M=500$) is sufficiently large to yield accurate estimates.\nPlease note that, due to the limited time available during the rebuttal period, the reported results are based on 10 trials. We will include the results based on 100 trials in the final manuscript.  \n\n**Can the non-negativity of the estimator be more rigorously enforced, e.g., through a post-processing projection or transformation?**\n\nAs discussed in the last paragraph of Section 3.1, applying $\\text{max}(\\hat{\\lambda}(x),0)$ or $\\text{ReLU}(\\hat{\\lambda}(x))$ is the simplest way to enforce the non-negativity of the estimator $\\hat{\\lambda}(x)$. Fortunately, this operation improves the estimator's accuracy in a pointwise sense, as it satisfies $|\\text{max}(\\hat{\\lambda}(x),0) - \\lambda(x)| \\leq |\\hat{\\lambda}(x) - \\lambda(x)|$ for any input $x$ due to the fact that the true intensity function $\\lambda(x)$ is always non-negative. However, when evaluating the integral $\\int_{S} \\lambda(x) dx$ over a compact region $S$, applying the max operation hinders closed-form integration and necessitates the use of computationally expensive Monte Carlo methods. Therefore, whether or not to employ the max⁡ operation depends on the specific application. We will incorporate the above discussion into the text. For clarity, we note that the max⁡ operation was not used in our experiments.\n\n**overlooks some recent work in Bayesian nonparametric methods...**\n\nWe added a study with a scalable Bayesian model (see our 3rd response to Reviewer upvZ). We promise to cite deep kernel/GP approaches, but it would be helpful if you could point us to a few references worth citing.", "claims_and_evidence": "The paper claims that:\n1. K2IE is theoretically consistent with classical KIEs under least squares loss.\n2. It provides comparable predictive performance to state-of-the-art methods while being computationally more efficient.\n3. The proposed estimator handles edge effects effectively via RKHS-derived equivalent kernels.\n\nThese claims are largely supported by the theoretical derivations and empirical experiments. However, the experiments are limited to synthetic data, and evidence on real-world applicability or robustness to noise and irregular event distributions is missing.", "methods_and_evaluation_criteria": "The least squares loss within RKHS is a reasonable and novel formulation for this problem, especially given its computational advantages over log-likelihood loss. The comparison with KIE and Flaxman’s estimator (FIE) is appropriate, and the use of metrics like integrated squared and absolute error (L2, |L|), along with CPU time, provides a fair evaluation.", "theoretical_claims": "The theoretical results are sound and well-supported through rigorous derivation. The connection to Fredholm integral equations and path integral representations is novel and mathematically grounded.", "experimental_designs_or_analyses": "The experiments are designed well for demonstrating performance on a range of synthetic intensities, with appropriate variations in data sparsity and observation domains. Use of both low-dimensional and moderate-dimensional settings is appreciated. Still, the omission of real-world datasets or more complex 3D/temporal domains limits broader validation.\n\nAlso, while hyperparameters are tuned via cross-validation, there is no ablation to show sensitivity to the number of random features (2M), which could impact approximation quality.", "relation_to_broader_scientific_literature": "This work is well-positioned in the literature on nonparametric Poisson intensity estimation, kernel methods, and RKHS theory. It builds directly on foundational work by Flaxman et al. (2017) and distinguishes itself by moving from maximum likelihood to least squares loss, bridging classical and modern techniques.", "essential_references_not_discussed": "The paper overlooks some recent work in Bayesian nonparametric methods for Poisson processes beyond what is cited, including approximate inference in deep Gaussian processes or deep kernel learning which could serve as competitive baselines.", "other_strengths_and_weaknesses": "Strengths:\n1. Theoretical originality in bridging classical KIE and RKHS estimators.\n2. Analytical tractability due to representer theorem and Fourier-based approximation.\n3. High computational efficiency and clear reproducibility through open-source code.\n\nWeaknesses:\n1. Limited to synthetic data.", "other_comments_or_suggestions": "na", "questions_for_authors": "1. Why are real-world datasets not included in the evaluation? This would be critical to support claims of practical relevance.\n2. How sensitive is the model to the number of random Fourier features (M)? An ablation study would help demonstrate robustness.\n3. Can the non-negativity of the estimator be more rigorously enforced, e.g., through a post-processing projection or transformation?", "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": "", "extra_scores": {}}, {"review_id": "3ixrenFxyH", "reviewer": "Reviewer_MUu1", "summary": "**Summary:**\nThe authors consider modelling the intensity function of a Poisson point process as belonging to an RKHS, and then fitting this based on a regularised squared error objective. This yields a method which is similar to other kernel intensity estimators (which were not motivated via RKHS), in that no \"model fitting\" or \"parameters\" need to be found once the data and hyperparameters are fixed. The technique achieves comparable predictive performance while being more computationally efficient than existing methods.\n\n## update after rebuttal\nI'm happy with the authors response --- the comments about the clipping and squared loss were particularly helpful. It looks like all reviewers are leaning accept or accept. I will maintain my current score of accept.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "We would like to thank the reviewer for the highly positive comments. Below, we provide a detailed response to each comment.\n\n**I list some possible related works for the authors' consideration, however these are definitely not essential to cite.....**\n\nWe appreciate the suggestion of including these important references. \"Sparse Spectral Bayesian Permanental Process with Generalized Kernel\" and \"Exact, Fast and Expressive Poisson Point Processes via Squared Neural Families\" both propose intensity estimation methods that employ the quadratic link function, $\\sigma(\\cdot) = (\\cdot)^2$, to ensure non-negativity of the estimators. We will appropriately cite these works in the Other Related Works section.\n\n**Is your squared loss related to equation (12) of \"PSD Representations for Effective Probability Models\"?**\n\nThis reference proposes a quadratic kernel model of the form $f_{\\text{PSD}}(x) = \\sum_{n,n'=1}^N A_{nn'}k(x_n,x)k(x_{n'},x)$, which leverages the non-negativity property of the positive semi-definite matrix $A$. The model is applicable to both density estimation and intensity estimation tasks. While the reference also considers learning the model within the framework of penalized least squares loss minimization, it focuses on fitting the parameter $A$ using iterative optimization methods and does not address whether the functional form of $f_{\\text{PSD}}(x)$ is optimal for that loss--indeed, it likely is not. In fact, (Marteau-Ferey et al., NeurIPS2020) showed that $f_{\\text{PSD}}(x)$ minimizes a certain functional loss, $L(f(x_1),\\dots, f(x_N))$, with appropriate regularization. However, this functional loss does not include the least squares loss, which involves the integral of the latent function. Although identifying the functionally optimal estimator under the least squares loss with similar regularization is an intriguing and important question, the work \"PSD Representations for Effective Probability Models\" does not appear essential in this specific context.\n\n**Weaknesses: As pointed out by authors in section 3.1, the method leads to intensity functions which can be negative, .... The authors use a post-hoc clipping of the intensity function below zero, after the fitting procedure, however this then clearly breaks the optimality of the solution.**\n\nIn fact, applying post-hoc clipping via $\\text{max}(\\hat{\\lambda}(x),0)$ consistently improves the accuracy of the estimator, since for any input $x$, the inequality $|\\text{max}(\\hat{\\lambda}(x),0) - \\lambda(x)| \\leq |\\hat{\\lambda}(x) - \\lambda(x)|$ holds due to the non-negativity of the true intensity function $\\lambda(x)$. The reason why the clipping, despite the fact that it breaks the optimality of the solution, improves the accuracy of K$^2$IE is clear: the non-negativity condition of intensity function is not taken into consideration in the problem of minimization of the penalized least squares loss (11). For more details, see our 3rd response to Reviewer XSJm.\n\n**Question: Could you add some more discussion for equation (10), beyond providing the citations?**\n\nIn response to the reviewer's suggestion, we will include the following explanation of the least squares loss (10), which we believe is satisfactory. Let $\\mathbb{E}$ be the expectation regarding data points generated from the true intensity function $\\lambda(x)$. Then we consider the expectation of the integrated squared loss between the estimator, $\\hat{\\lambda}(x)$, and the true intensity $\\lambda(x)$, defined by \n$$\\mathbb{E} \\Bigl[ \\int_{\\mathcal{X}} \\bigl|\\hat{\\lambda}(x) - \\lambda(x) \\bigr|^2 dx \\Bigr] = \\mathbb{E} \\Bigl[ \\int_{\\mathcal{X}} \\hat{\\lambda}^2(x) dx \\Bigr] - 2\\mathbb{E} \\Bigl[  \\int_{\\mathcal{X}} \\hat{\\lambda}(x) \\lambda(x) dx \\Bigr] + \\mathbb{E} \\Bigl[ \\int_{\\mathcal{X}} \\lambda^2(x) dx \\Bigr].$$ \nThe third term can be omitted, as it is independent of the estimator. The second term can be decomposed into two parts as\n$$2\\mathbb{E} \\Bigl[ \\int_{\\mathcal{X}} \\hat{\\lambda}(x) \\lambda(x) dx \\Bigr] = 2\\mathbb{E} \\Bigl[ \\int_{\\mathcal{X}} \\hat{\\lambda}(x) \\sum_{n=1}^N \\delta(x-x_n) dx \\Bigr] + 2\\mathbb{E} \\Bigl[ \\int_{\\mathcal{X}} \\hat{\\lambda}(x) \\bigl( \\lambda(x) - \\sum_{n=1}^N \\delta(x-x_n) \\bigr) dx \\Bigr],$$ \nwhere the second term vanishes due to Campbell’s theorem:\n\\begin{equation}\n2\\int_{\\mathcal{X}} \\mathbb{E} [\\hat{\\lambda}(x)] \\lambda(x) dx - 2\\sum_{n=1}^N \\mathbb{E} [\\hat{\\lambda}(x_n)] = 2\\int_{\\mathcal{X}} \\mathbb{E} [\\hat{\\lambda}(x)] \\lambda(x) dx - 2\\int_{\\mathcal{X}} \\mathbb{E} [\\hat{\\lambda}(x)] \\lambda(x) dx = 0.\n\\end{equation}\nPutting everything together, we obtain the following (correct) interpretation of Eq. (10): \n$$\\mathbb{E} \\Bigl[ \\int_{\\mathcal{X}} \\bigl|\\hat{\\lambda}(x) - \\lambda(x) \\bigr|^2 dx \\Bigr] = \\mathbb{E} \\Bigl[ \\int_{\\mathcal{X}} \\hat{\\lambda}^2(x) dx - 2\\sum_{n=1}^N \\hat{\\lambda}(x_n) \\Bigr] + C,$$ where $C$ is the constant term.", "claims_and_evidence": "The claims are clear and supported by evidence.", "methods_and_evaluation_criteria": "Yes, the evaluation criteria make sense and are standard for this area. The benchmark with the rectangular domains is particularly nice.", "theoretical_claims": "- This paper bridges an important theoretical and conceptual gap between two \"kernel\" approaches to intensity estimation --- RKHS (e.g. Flaxman et al.) and kernel intensity estimators (like kernel density estimators). In particular, with a regularised least squared loss, RKHS methods actually give a variant of kernel intensity estimation methods! This is actually quite an inspiring result - I wonder what other loss functions with RKHS models yield?", "experimental_designs_or_analyses": "Yes, they are sound.", "relation_to_broader_scientific_literature": "This is well-placed within broader scientific literature -> stats/ml -> point processes and machine learning. I list some possible related works for the authors' consideration, however these are definitely not essential to cite.\n\t- In  Bayesian context, the intensity is modelled as a squared Gaussian process and Random Fourier features are used in \"Sparse Spectral Bayesian Permanental Process with Generalized Kernel\", as well as a Laplace approximation to the posterior.\n\t- Yet another type of \"kernel method\" for estimating intensity functions of inhomogeneous Poisson point processes is using \"\"Exact, Fast and Expressive Poisson Point Processes via Squared Neural Families. These are squared neural networks with closed-form integrated intensity functions, trained using maximum likelihood.\n\t- Is your squared loss related to equation (12) of \"PSD Representations for Effective Probability Models\"? They are looking at density estimation rather than intensity estimation, but apart from that it looks a little bit similar.", "essential_references_not_discussed": "All essential works are cited, to the best of my knowledge.", "other_strengths_and_weaknesses": "**Strengths:**\n- This paper bridges an important theoretical and conceptual gap between two \"kernel\" approaches to intensity estimation --- RKHS (e.g. Flaxman et al.) and kernel intensity estimators (like kernel density estimators). In particular, with a regularised least squared loss, RKHS methods actually give a variant of kernel intensity estimation methods! This is actually quite an inspiring result - I wonder what other loss functions with RKHS models yield? \n- Text, equations and figures are easy to follow, and seem to be without any major errors.\n- I particularly liked the composite domains built from multiple rectangles in Figure 3.\n\n**Weaknesses:**\n- As pointed out by authors in section 3.1, the method leads to intensity functions which can be negative, due to the fact that they do not utilise a nonnegative \"link function\". This can lead to undesirable effects from a modelling perspective (e.g. predict a negative number of events). The authors use a post-hoc clipping of the intensity function below zero, after the fitting procedure, however this then clearly breaks the optimality of the solution. \n- The equivalent kernel has to be approximated (e.g. Fourier features, MC), leading to approximate estimation procedure.", "other_comments_or_suggestions": "See above.", "questions_for_authors": "- Question: Could you add some more discussion for equation (10), beyond providing the citations? It looks like the squared L2 distance between the intensity function and a constant 1 function, which is then approximated by replacing the integrated intensity by a sum of intensities evaluated at the data. Is that the \"correct\" interpretation? Why is this a helpful loss function?", "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": "", "extra_scores": {}}, {"review_id": "yzqktee8hC", "reviewer": "Reviewer_eKVi", "summary": "The paper proposes K²IE, a kernel method-based intensity estimator for inhomogeneous Poisson processes, combining the computational efficiency of classical kernel intensity estimators (KIEs) with edge-correction capabilities from reproducing kernel Hilbert spaces (RKHS). By reformulating the problem using a penalized least squares loss in an RKHS, the authors derive an estimator that matches the form of classical KIEs but uses equivalent RKHS kernels h(⋅,⋅) to implicitly handle edge effects. Theoretical analysis shows the solution satisfies a specialized representer theorem with unit dual coefficients. Experiments on synthetic 1D/2D datasets demonstrate comparable accuracy to Flaxman’s kernel method-based estimator (FIE) but with significantly improved computational efficiency.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "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": "We thank the reviewer for giving valuable comments. Below, we provide a detailed response to each question. We will include all discussions in the revised manuscript.  \n\n**The equivalence between ( q(\\cdot,\\cdot) ) and ( h(\\cdot,\\cdot) ) in Equation 8 is asserted but lacks a proof ...  Can this be rigorously shown without invoking path integral ...?**\n\nYes, the equivalence between $q(\\cdot,\\cdot)$ and $h(\\cdot,\\cdot)$ can alternatively be established via Mercer's theorem, which aligns with the \"rigorous\" approach proposed by (Flaxman, 2017). Due to space limitations, we provide only a brief overview of the derivation here. Let $\\lambda(x)$ be a function in the RKHS associated with the kernel $k(x,x')$, and let $||\\lambda||^2_{H_k}$​ denote its squared RKHS norm. Then, following the notation in Eq. (7) in Section 2.2, we can express $\\lambda(x) = \\sum_{m=1}^{\\infty} b_m \\eta_m$ and $||\\lambda||^2_{H_k} = \\sum_m b_m^2 /\\eta_m$, where $b_m$ are coefficients. Substituting this into the penalized least squares loss in Eq. (11), we obtain: $-2 \\sum_{n} \\lambda(x_n) + \\sum_m \\frac{\\eta_m+1/\\gamma}{\\eta_m} b_m^2$. We can see that the 2nd term corresponds to the squared RKHS norm under a rescaled kernel defined by $q(x,x') = \\sum_m \\frac{\\eta_m}{\\eta_m + 1/\\gamma}e_m(x) e_m(x')$. It is evident that $q(\\cdot,\\cdot)$ is consistent with $h(\\cdot,\\cdot)$ in Eq. (8). A full derivation of Theorem 1 will be included in the text.\n\n**How prevalent are negative values in practice?**\n\nFollowing the reviewer’s question, we conducted an analysis of how frequently K$^2$IE produces negative values using the 2D synthetic dataset $\\lambda_{\\text{2D}}^{1.0}$. Specifically, we evaluated the estimated intensity values at 500 x 500 grid points within the observation domain and computed the ratio of negative values. The mean ± standard deviation of this ratio across 100 trials was $0.059 \\pm 0.016$. This result indicates that K$^2$IE can indeed produce negative estimates in practice—particularly in regions with sparse data—highlighting the necessity of post-hoc clipping like $\\max(\\hat{\\lambda}(x), 0)$ in applications where negative intensity values are not permitted. It is also worth noting that when Laplace RKHS kernels are used, the equivalent kernel $h(x, x')$ is functionally non-negative in one-dimensional input settings (see Section 3.2 in (Kim, NeurIPS2024)). \n\n**What fundamentally new theoretical insight does K$^2$IE offer?**\n\n1. To the best of our knowledge, our paper is the first to prove that minimizing the least squares loss with a squared RKHS norm regularizer yields kernel intensity estimators (KIEs).\n\n2. K$^2$IE demonstrates that the equivalent kernel used in Flaxman (2017) can serve as an edge-corrected smoothing kernel for KIE. This insight enhances computational efficiency, as it eliminates the need to solve the dual optimization problem (9) in Flaxman's model.\n\n**The claim \"no model fitting\" is misleading ... How does hyperparameter tuning affect computational efficiency claims?**\n\nWe will rephrase \"no model fitting\" as \"no optimization of dual coefficients\". Regarding hyperparameter tuning, K$^2$IE offers a significant advantage over the reference models. Specifically, KIE and FIE require MC integration and solving a dual optimization problem for each cross-validation, respectively, whereas K$^2$IE requires neither, which is beneficial especially in multi-dimensional settings.\n\n**How does RFF rank ( M ) trade off edge-correction accuracy?**\n\nPlease refer to our 2nd response to Reviewer XSJm regarding the ablation study on $M$. We want to emphasize that the edge correction (i.e., the integral over the domain in Eq. (8)) is performed exactly under the RFF approach (see the 2nd paragraph of Sec. 3.2.2). Rather, $M$ governs the approximation accuracy of shift-invariant kernels.\n\n**How does (\\gamma) balance edge effects and overfitting?**\n\nWe appreciate the important question. In classical KIEs, tanking $\\beta \\to \\infty$ ($\\beta^{-1}$ is the scale parameter) leads to $\\hat{\\lambda}(x) = \\sum_n \\delta(x-x_n)$. K$^2$IE exhibits the same behavior regardless of $\\gamma$, but taking $\\gamma \\to \\infty$ also yields the same solution regardless of $\\beta$ in K$^2$IE. Although $\\gamma$ as well as $\\beta$ controls the degree of overfitting, it seems that $\\gamma$ acts globally and $\\beta$ locally. A thorough investigation into the distinct roles of $\\gamma$ and $\\beta$ remains an important direction for future work. We believe the explanation is satisfactory to the reviewer.\n\n**Experiments are limited to synthetic 1D/2D data.**\n\nPlease see our 1st response to Reviewer XSJm.\n\n**Standardize notation**\n\nWe will standardize the notation according to the suggestion.\n\n**How does K$^2$IE compare to modern Bayesian nonparametrics?**\n\nPlease see our 3rd response to Reviewer upvZ.\n\n**does not clarify whether KIE’s edge correction was optimized**\n\nAs in Sec 4, KIE's hyperparameter was optimized similarly to K2IE but based on test likelihood.", "claims_and_evidence": "Supported Claims:1. K²IE achieves computational efficiency comparable to classical KIEs (evidenced by CPU time in Tables 1-2).\n2.Edge correction via equivalent RKHS kernels works effectively in multi-dimensional settings (supported by 2D results in Table 2).\n\nProblematic Claims:1.The assertion that K²IE \"combines the computational efficiency of KIEs with the effectiveness of Flaxman’s estimator\" is partially unsubstantiated. While K²IE is faster than FIE, its edge-correction superiority over KIE is only qualitatively shown (Figure 1), lacking quantitative comparison.\n2.The claim that K²IE \"does not require model fitting\" is misleading, as hyperparameter tuning (e.g., γ, β) is still necessary.", "methods_and_evaluation_criteria": "Strengths:\n1.The degenerate approximation using random Fourier features (Eq. 17) is a practical approach for solving the Fredholm equation.\n2.Evaluation metrics (L^2, ∣L∣) align with standard practices for intensity estimation.\n\nWeaknesses:\n1.Experiments are limited to synthetic data with high SNR. Real-world datasets are absent, raising concerns about generalizability.\n2.The cross-validation setup uses p-thinning, which may introduce bias if events are temporally/spatially correlated.", "theoretical_claims": "Theorem 1’s proof relies on path integral representations and operator inversions (Eq. 13–17). While the derivation is logically consistent, the critical step of connecting the least squares loss to the Gaussian process representation lacks rigor. Specifically: 1.The path integral representation of the RKHS norm (Kim, 2021) is cited but not explicitly justified in the context of Poisson processes. 2.The equivalence between q(⋅,⋅) and h(⋅,⋅) (Eq. 8) is asserted without proving uniqueness or convergence under finite-dimensional approximations.", "experimental_designs_or_analyses": "1.The comparison with KIE uses edge-corrected kernels for KIE but does not clarify whether KIE’s edge correction was optimized similarly to K²IE’s hyperparameters.\n2.The regularization parameter γ’s impact on edge correction is not analyzed. For example, does γ→0 degrade h(⋅,⋅) to an uncorrected kernel?\n3.Negative intensity values (Section 3.1) are dismissed without quantitative analysis (e.g., frequency/severity of negative estimates in experiments).", "relation_to_broader_scientific_literature": "The work bridges classical kernel smoothing (Diggle, 1985) and modern RKHS-based methods (Flaxman et al., 2017). By connecting the least squares loss to the representer theorem, it extends the theoretical framework of Walder & Bishop (2017) for Cox processes. However, it does not engage with recent advances in neural point processes or scalable Bayesian methods.", "essential_references_not_discussed": "Scalable Bayesian Methods: Lloyd et al., Variational Inference for Gaussian Process Modulated Poisson Processes (ICML 2015) – omitted despite being a key prior work.", "other_strengths_and_weaknesses": "1.The connection between least squares loss and unit dual coefficients is novel, though the core idea (equivalent kernels for edge correction) builds directly on Flaxman et al. (2017).\n2.Provides a computationally efficient alternative to FIE but does not surpass classical KIE in low dimensions.\n3.The writing is dense, with inconsistent notation (e.g., λ used for both true and estimated intensity functions).", "other_comments_or_suggestions": "The authors need to carefully check all the details making sure no typos included.", "questions_for_authors": "1. The equivalence between \\( q(\\cdot,\\cdot) \\) and \\( h(\\cdot,\\cdot) \\) in Equation 8 is asserted but lacks a proof of uniqueness or convergence under finite-dimensional approximations. Can this be rigorously shown without invoking path integral heuristics (Kim, 2021)?  \n2. Section 3.1 acknowledges potential negative estimates but dismisses them without quantitative analysis (e.g., frequency/severity in experiments). How prevalent are negative values in practice?  \n3. The edge-correction mechanism via solving Equation 8 directly extends prior work. What fundamentally new theoretical insight does K²IE offer?  \n4. The claim \"no model fitting\" is misleading since \\(\\gamma\\) and \\(\\beta\\) require cross-validation. How does hyperparameter tuning affect computational efficiency claims?  \n5.  Equation 16 uses random Fourier features (RFF) to approximate \\( k(\\cdot,\\cdot) \\). How does RFF rank \\( M \\) trade off edge-correction accuracy?  \n6.  Does \\(\\gamma \\to 0\\) reduce \\( h(\\cdot,\\cdot) \\) to an uncorrected kernel? How does \\(\\gamma\\) balance edge effects and overfitting?  \n7. Experiments are limited to synthetic 1D/2D data. How does K²IE perform on real-world data with irregular domains?  \n8. Figure 1 qualitatively compares kernels but lacks quantitative metrics (e.g., \\( L^2 \\) error near boundaries). Why?  \n9.\\( \\lambda(\\cdot) \\) denotes both true and estimated intensity functions (Section 3.1). Standardize notation (e.g., \\( \\lambda^*(\\cdot) \\) vs. \\( \\hat{\\lambda}(\\cdot) \\)).  \n10.Scalable Bayesian methods (e.g., Lloyd et al., ICML 2015) are ignored. How does K²IE compare to modern Bayesian nonparametrics?\n\nI will be willing to improve the scores if the authors give me good answer of the above questions.", "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": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This paper demonstrates a correspondence between classical kernel intesnity estimators in non-homogeneous Poisson proceesses and those based on repreoducing kernel Hilbert space theory.

The reviewers broadly agree that the methodology is theoretically sound. Though one reviewer requests a more detailed derivation without path-integral heuristics.
Additionally, all reviewers remark on the issue that the estimated intensity might become negative. It was also remarked that the error in the random Fourier feature approximation was not investigated sufficiently.

Furthermore, most reviewers agree that there are gaps in the literature review and that more comparisons with recent statistical/ML approaches would be beneficial.
Additionally some reviewers would like an empirical investigation on the possible negativity of the estimated intensity and on the performance in extrapolation.

Some of the key-points that ought to be addressed in a camera-ready version are:

* Discuss related work more thoroughly in accordance with reviewers' suggestions.
* Discuss "clipping" the intensity estimator and inverse M-kernels in more detail.
* Empirical investigation of the sign of the intensity estimator.
Decision: Accept (poster)