Fast Bayesian Estimation of Point Process Intensity as Function of Covariates
Authors: Hideaki Kim
Taichi Asami
Hiroyuki Toda
Venue: NeurIPS 2022 Accept (2022)
Keywords: Point Process
Gaussian Cox Process
Kernel method
Permanental Process
Gaussian Process

Abstract: In this paper, we tackle the Bayesian estimation of point process intensity as a function of covariates. We propose a novel augmentation of permanental process called augmented permanental process, a doubly-stochastic point process that uses a Gaussian process on covariate space to describe the Bayesian a priori uncertainty present in the square root of intensity, and derive a fast Bayesian estimation algorithm that scales linearly with data size without relying on either domain discretization or Markov Chain Monte Carlo computation. The proposed algorithm is based on a non-trivial finding that the representer theorem, one of the most desirable mathematical property for machine learning problems, holds for the augmented permanental process, which provides us with many significant computational advantages. We evaluate our algorithm on synthetic and real-world data, and show that it outperforms state-of-the-art methods in terms of predictive accuracy while being substantially faster than a conventional Bayesian method.

Fast Bayesian Estimation of Point Process Intensity as Function of Covariates
Hideaki Kim
Taichi Asami
Hiroyuki Toda*
NTT Human Informatics Laboratories
NTT Corporation
(hideaki. kin. en, taichi asami. &hco .nt..c.. co. jp, hirotodadacm. org
Abstract
In this paper, tackle the Bayesian estimation of point process intensity function covariates. We propose novel augmentation of permanental pro- called augmented permanental process, doubly-stochastic point process that uses a Gaussian process on covariate space to describe the Bayesian a pri- ori uncertainty present in the square root of intensity, and derive fast Bayesian estimation algorithm that scales linearly with data size without relying on either domain discretization or Markov Chain Monte Carlo computation. The proposed algorithm is based on non-trivial finding that the representer theorem, one of the most desirable mathematical property for machine learning problems, holds for the augmented permanental process, which provides with many significant computational advantages. We evaluate our algorithm on synthetic and real-world data, and show that it outperforms statc-of-the of-tth-artttttt methods in terms of predictive accuracy while being substantially faster than a conventional Bayesian method.
Introduction
Point processes are widely-used statistical tools to analyze point patterns observed on temporal, spatial, and spatio-temporal domains. They have wide ranging list of applications in neuroscience [9], finance [6), epidemiology [14],, ecology [II]], and others.
The essential task in point processes is to estimate from event data an intensity function, that is, an instantaneous probability of events occurring each point in the observation domain. Generally, there are two kinds of scenarios in the estimation task: one is to estimate intensities functions of every point in the observation domain; the other is to estimate intensitics as functions of covariates with the assumption that covariates are given at every point in the observation domain. For example, in ecological applications where the event points are the locations trees, the former estimates the value of intensity at each geographical point (longitude and latitude), while the latter estimates the intensity as function of spatial covariates such as slope, elevation, and humidity, each of which varies across the geographical domain. The latter scenario has greater practical importance because en- ables us to investigate factors that may be causing events. But unfortunately, it has been overlooked in the recent machine learning community. This paper focuses using the function of covariates to solve the estimation problem of intensity.
The simplest and most classical model used to estimate intensities functions of covariates is the Cox proportional hazards model [8, 43, [0] (or the log-linear model), which forms the logarithm of the intensity function easy to handle, fails to address the nonlinear dependence of point patterns on covariates, which is often the case in real-world data. To rectify this limitation, two nonparametric approaches have been
linear combination of covariates. Although this parametric model is very
36th Conference on Neural Information Processing Systems (NeurIPS 2022).
"Current affiliation is Yokohama City University.

---

proposed: Kernel-based intensity estimator (KIE) [ IR D], which is rescaled version of kernel density estimator, the gold standard in estimating intensity as function of covariates in a non- Bayesian way; it is easy to implement and rapid in computation. A Bayesian alternative is Gaussian Cox process (GCP) [299 37, 22, 21, 28, IIIl in which flexible prior over latent intensity es- tablished by Gaussian process
random field [[5] through a positive link function. GCPs are promising approaches as they allow the estimation of the uncertainty of intensity and the hyper- parameter optimization in statistically formal way. However, the Bayesian estimation for GCPs is highly challenging as the likelihood function depends on the functional form of the latent intensity distribu-
over compact domain, which imposes the intractable computation of tions. Most algorithms published to date have dealt with the difficulty by using domain discretization or/and Markov Chain Monte Carlo (MCMC). But domain discretization suffers from poor scaling in terms of domain dimension, as well as being sensitive to discretization size. Moreover, MCMC has significant computation demands, making GCPs less attractive than KIEs.
In this paper, for the problem of estimating intensity as function of covariates, we propose a novel Bayesian estimation scheme for GCPs with quadratic link functions without recourse to either do- main discretization or MCMC. GCPs with quadratic link functions are called permanental processes [27] when the Gaussian process for intensity is defined in the observation domain. In our proposed model. the Gaussian process for intensity should be defined in the covariate domain, and thus we call it the augmented permanental process (APP). We show that the representer theorem [_ one of the most desirable mathematical property for machine learning problems, holds for the APP, which is one of the main contributions of this paper. The finding provides with many significant com- putational advantages, resulting in a fast Bayesian estimation scheme that scales linearly with the number of data points. Furthermore, based on Laplace approximation in the functional space, we derive the predictive distribution and the marginal likelihood in a feasible form, which enables us to implement the uncertainty evaluation and the hyper-parameter optimization
Bayesian manner.
In Section E we introduce APP and construct a scalable Bayesian estimation scheme for it. In Section B we outline related work. In Section H. we compare APP against KIE- and GCP-based algorithms on synthetic and real-world data, and confirm that our scheme achieves higher predictive accuracy than the state-of-the-at algorithms well as substantial speed improvements over the conventional GCP-based algorithm. Finally, Section 5 states our conclusions.
2 Methods
2.1 Augmented Permanental Processes
Suppose that a dataset of N point events and a map of Dy-dimensional covariate, denoted by D = (tn, E TT) and (t) - c CRD, respectively, are observed for a compact Dy-dimensional 1=
domain, T RD, We consider variant of Gaussian Cox processes which assumes that latent function, x0y) y + IR. is generated from a Gaussian process (GP) on V, and point events are generated on T from a point process with intensity function, X(t), so that
(V)) (P()()) p(a(y))D) J Drlly) P (D)((y)) (P()())k))
N
log (Dpp(D(x(8) (z)) = Σ log A(tn) - / A(t)dt, X(t) = (00=x(x(y(())) n(z(v(t))).
(2)
n=1
where (P(a(y)) kb) represents GP with kernel function k(y. y²), p(D) xr (r))) is the likelihood func- tion of the point process, k(z)= R + E4 non-negative function called the link function, while f Tr(y) in the denominator represents the integral over the function or the infinite-dimension variable 2gy). Among popular link functions such as exponential and sigmoidal, this paper adopts the quadratic form of the link function, x(x) x² for its computing benefits. We call the model defined by (I-2) with the quadratic link function the augmented permanental process (APP), be- cause it reduces to the ordinary permanental process when the covariate map is the identity mapping (t) = t. Intensity X(t) represents an instantaneous probability of events occurring at each point on T, and our goal is to estimate, based on APP, the functional form of the intensity over covariate domain, x(x(y))= = z² 22(),, from the observed data, {T, (t,, S3 in the supplementary material.
1N 1 y()). For illustration, see Figure
2

---

In general, Bayesian estimation on normal Gaussian Cox processes, equivalent to (I-22 with y(t) often referred as Kloubly-intractaile function æ(y) as well as computing the intractable posterior probability of an infniteedimensionall variable. Furthermore, APP is challenged by one more difficulty: the domain on which the GP defined ()) differs from that over which the intensity integral is performed (T), which degrades mathematical tractability. Below, we address this "triply-intractall" problem with the help of path integral representation, mathematical tool that has been proposed most recently [23]].
because requires solving the integral of the stochastic
2.2 Maximum A Posteriori Estimator
We consider the problem of obtaining the maximum posteriori (MAP) estimator of x(y) that maximizes the posterior probability (I)) First, we solve the discrepancy between the domains for the GP in (II) and the intensity integral in the likelihood function (LD, as it makes functional analyses problematic. If the covariate map (t) is bijective, it is trivial that the intensity integral over 'T in (L) can be transformed into the integral over the covariate domain as, f dt Jy - Idy. where Jy is the Jacobian matrix of the map yt). We can generalize the transformation for a non-bijective map, which often happens in practice (e.g. the value of covariate varies periodically over time), as follows:
12==--200000 ((((00)) |
[e y))y = ((())=(()dd)
where 5(-) is the Dirac delta function, and p(v) = f ((y(t) - y d). We call p(y) the weight distribution of the covariate. Actually, our proposed algorithm does not need to evaluate p(y) directly, and thus do not go into detail about how compute p(y)- Rather, emphasize that the following transformation holds for arbitrary functions of covariates, f(y):
J. plu))(w)ds | f(y(())dt.
(4)
Next, we derive the MAP estimator of sty), denoted by #(y), by using the path integral representa- tion of GP [24]]
(P(x(y)Ik) 2x(y)= V [C] exp - 1 k²" k" y ('))(()))())))))/o)))) "dydy' Dx(y),
(5)
where K. and Jy (yy, 3/) -dy are the integral operator with kernel function k(y, y') and its inverse operator, respectively: x()) = Jy k y (ay')dy' KCk (y, y') 6(y - y'). Here [K] represents the function determinant [|.S of K., defined by the product of its eigenvalues. Using the
relation (B) and the representation (5), we write the posterior of APP (MI2) in the following functional form,
p(r(y)))) 1
((D) exp (((()))) Φr,
(6)
where 5(x(y), z()) is the action integral, defined by
((x(v),z(u)) = / x ((v))(() ply)a (v)
Σ
N
log x bnd) dy + log [K], (7)
yn y(t)), and r(y) = fy (y. y )a(() dyy Then we apply calculus of variations the action integral, where the functional derivative of (vx()))())) on the MAP estimator f(y) should be 6s equal to zero: o dz( form,
si(y) O. which results in the exact MAP estimator in feasible kā(w)
69
((v) = 2 Σ h(y yn),,, Un (1
(8)
where hyy y²) is a transformed kernel function defined by Fredholm integral equation of the second kind [32].
My. v') + 2 S kly, s))(s))(s g'dds k(y, v").
3

---

See the supplementary material ($1) for the detailed derivations. Equation (S) shows that the MAP estimator augmented permanental process involves the representer theorem under the transformed kernel function hy, y²), and thus the Bayesian estimation reduces to miiie-dimensionll optimiza- tion problem. Here, we call h(y, 3/' the equivalent kernel following studies by Flaxman et al. [[3]]
and Walder & Bishop [H Given the equivalent kernel, the unknown coefficients Un in (S) solve the simultaneous quadratic equations derived from (S).
N
A
- 2 Cn
Σ n==11
hyyn: 3rd Jn O, 1.2,...,
(10)
In this paper, we estimate set of coefficients, {En by solving minimization problem of the mean of the squared residuals, Σ AA.l//,, with popular gradient descent algorithm, Adam [25]] We discuss how to obtain hy,, y in Section 2
Computational Complexity
The objective function to be minimized MAP estimation is Σ AnI?/NN which, in gradi- descent algorithms, naively demands the computation of ((N2) for each evaluation. How- ever, when the equivalent kernel is given in degenerate form with rank M (< N) such that h(y.v') Σm=1 Pm Pm (my'), the computational complexity reduces O(N linear com- putation against data size N. See Section 233 for the construction the equivalent kernel h y²).
2.3
Construction of Equivalent Kernels
2.3.1 Naive Approach
The equivalent kernel hyy, y²) solves the Fredholm integral equation (H), which includes the weight distribution of covariate, p(y). Although it is possible evaluate p(y) numerically, we take an approach that avoids to treat p(y) directly. Through the relation (H)), we rewrite the integral term in (9) as fy k(y, (hy)) y )ds = S k(y, ((y))) y d'), and apply to it the Nyström method [30], a popular method solve integral equations that approximates the integral operator by J-point numerical integration:
hly, v') + wy wyk(y y: (y s') k(y.v'), y - ytt;)
(11)
where T and wy E IR represent the integration points and weights, respectively, It is easily shown that the equation (LI) can be solved an analytical form as follows,
h(y. y²) = k(y.))) k(y) T (w K) k(y'),
(12)
where k(y) (k(y. 81) .1)....... k(y T. Wjr = 2 and K3 = k(y By for S J.j <
J. We call (122) the naive approach to constructing equivalent kernel, and the complexity of eval- uating ky., y²) is ((J3). ). this paper, we use the quasi-Monte Carlo method [[]], where {4,)/=1 low-discrepancy sequence (Sobol sequence was used here) and w, = [T/J for the volume of the continuous domain T. However, other techniques like Monte Carlo methods and random quadrature approximation [S0] are also possible.
2.3.2 Degenerate Approach
When the kernel function of GP has degenerate form with rank M(< oo),
k(y.v') Σ (m (b(('))
(13)
is well known that the Fredholm integral equation (g) can be solved analytically [B] as,
hly y') x (IM + 2A) Φlv'),
(14)

---

where )(y) = (0) (x())) and A is the MxM matrix defined by A = Jy (())((y)(())) dy. As in the naive approach, we rewrite and approximate the integral oper- ator A by using relation (+) and a J-point numerical integration, resulting in
A Σ w w, y/ y(t,),
(15)
where i, T and Wj E I represent the integration points and weights, respectively. We call (IE (.5) the degenerate approach for constructing an equivalent kernel. The computational complexity of evaluating hly, y') is (M3 JM²), where ((MM) is given by numerical integration (I.))
should be emphasized here that when k(y, y²) has degenerate form with rank M, the relation (TE)) shows that the equivalent kernel hyy, y²) also has degenerate form obtained through Cholesky decomposition:
hyy,y') (Lo(y)) (LO()')) (Lo()))(Lo('')))
L'L (IM + 2A)-1
(16)
The degenerate equivalent kernel (IG) offers fast Bayesian estimation that scales linearly with N, the number of data points (see Section 2 and Section ME!)
There are several ways of preparing degenerate kernel functions 94, 45]] approximates popular kernel functions, such as shift invariant kernels and additive kernels, by the inner products of randomized feature map vectors. Also, Nyström approximation for kernel methods [M 3 finds degenerate form of the function that best approximates an arbitrary kernel function at a finite subset of data points. We apply the two approximation methods to a Gaussian kernel in our experiments.
Random feature map [53]
2.3.3 Discussion about Approximation Error of Integral Operator
We approximate the integral equation (g) through J-point numerical integration, which generally requires a large J when the dimension of the covariate space, V, is high. However, it is worth noting that in our approaches, the required number of integration points, J, does not grow with the dimension of y (covariate domain), but does with the dimension of T (observation domain). It is easily verified by the fact that the integral operator over y in (H) can be transformed into the operator over 'I through the relation (H)) and that the integration points are selected not on ) but on T. Point processes are usually utilized to analyze at most dimensiaaan data (T C ), and thus the required J practically stays within feasible size regardless of the dimension of V An intuitive explanation for the requirement of / is that observed covariate yt) exists on a lower-dimensional manifold in high-dimensional domain V, and the effective volume of integral region is rather small.
The discussion far is about the approximation error of the integral operator. Next, we discuss about the approximation eтror of the equivalent kernel under the approximate integral. Let the integral operator in (E) and its approximation through J-point numerical integration be denoted by K' and K'',, respectively. Then the equivalent kernel to solve (9) has the relation (for detailed derivation, see [BI],
h. < (I 21C') (K'h K'rh), h h where h and h, are the solutions of (H) under KC' and K',, respectively. The relation shows that the speed of convergence of h, is at least as rapid as the speed of convergence of the integral operator, K h to K'h. Therefore, the required J to obtain a good approximation of the equivalent kernel also stays within a feasible size, as long as the dimension of T low.
(17)
2.4 Predictive Distribution and Marginal Likelihood
One of the advantages of GP models over non-Baycsian approaches is that they can provide pre- dictive distributions and marginal likelihoods, which enable to perform uncertainty evaluations, hyper-parameter optimization, and model selection in a Bayesian manner. We apply to APP (I-2) Laplace approximation in the functional space, and find the approximate form of the predictive distribution and the marginal likelihood.
We now know the mode of the posterior, k(y), and consider a Taylor expansion of functional action potential S (z(()) I(y)) centered on the mode such that
S (((()) = S((()) (()) + 2 f/fl yxx
o (w -a(y -a('))dyyy''' (18)
5

---

is the second derivative of S. The first term in the Taylor pansion vanishes due to the stationary condition. The quadratic approximation of the action integral corresponds to the approximation of the posterior process by GP, and the predictive covariance or the kernel function for the posterior GP, denoted by (y, y²), can be obtained by the functional inversion of o o"(,, y'). which results in
where o" (y.y") )n(y)az(yy) l==8
S(a,)))
ovy. v)) = h(,,)') - -(y) (Z + H) 'n(y',
(19)
where
naw = 8ymn', Haw = h(yn, yr,) h = (h(y.y))
My).
(20)
The full derivation of (199) is given in the supplementary material ($2). When the latent function (ty) follows posterior GP with mean of ity) and the kernel (y y"), it is easily verified that the value of the squared function, À= x² (z), on each point of the covariate domain y E y follows a
Gamma distribution,
24 to
( + o)² (2 + +a)
(p. T(V) X- exp(-2/p), 1 == 2a
(21)
2 to
where è and are abbreviations of (yy) and o (y y). respectively. distribution of the intensity as function of covariates.
represents the predictive
Furthermore, under Laplace approximation (IES, we can obtain the marginal likelihood, p(D), in (G) by performing the path integral as,
log (DD) = -log / exp|-S(>()))
1 2x -5(2(g). + log |1),
(22)
where |X]] the functional determinant of integral operator I Sy- ay y dyy' We can rewrite the result in a more tractable form by substituting (77 (9] into (222):
log = log| Z - log IN + Z-' + #2 (log | - log [K]) + (log 2 - 1)M,
(23)
where IN is the identity matrix with size N, and [H] is the functional determinant of integral oper- ator H = Jy- hy, y d'y'' The full derivation of (23) is given in the supplementary material ($3), and the method of estimating the functional determinant term is discussed in the end of this section.
Computational Complexity
Given the equivalent kernel h(y,y'), the evaluation of the predictive covariance (LL)) for covariate pair, (y. y'), naively needs the computation of O(N3 + N2) due to the matrix inversion in the second term. Fortunately, the complexity reduces O(M + N2) given a kernel with degenerate
M5 form, h(y. y') Em=1 Pm Py) (3'), through the Woodbury matrix identity: The N x X matrix
H can be decomposed into product of N M matrix R. defined by Rum = Pm 4m (Yn), and its transpose as H RR the matrix inversion is transformed as (Z + RR )-1 Z-' "R(IM + R Z R) 'R Z-1. Note that Z RR costs O(NM) because Z is a diagonal matrix. Therefore, the computational complexity of the predictive distribution (2I) is O(N3 + MN2) naively and ((M3 + NM2 with degenerate equivalent kernel.
Z-1
In computing the marginal likelihood (23), the complexity of the first and the second terms is O(N3 MN2) with a naive computation, but reduces to ((M + NM2) with degenerate equiv- alent kernel. The complexity of the third term, which is related to the functional determinants, depends on how the equivalent kernel hy yy) is constructed. Under the naive approach (I2), the functional determinant of operator, H Sy h(y,y')dy y''dy', can be obtained as product of the functional determinant of IC and J J matrix determinant: |
[ + K (see the supplementary material ($4)). Thus the functional determinant term in PX is given by the tractable determinant: (log [7[| - log [K])/2 11/2 log I,, + W K), of which computational complexity is OJ + .M 22). Under the degenerate approach (LE)) the functional determinant of operator H. can be obtained as a product of the functional determinant of IC and M M matrix determinant: |H| |K] ar 2A| (see the supplementary material ($4)). Thus the functional determinant term (2K) is given by the tractable determinant: (log ///- - log [Kl)/2 = 1 log Tas + 2AA,, of which computational complexity is -MM3-JM² where ( JMM 2) is given by numerical integration (I)).

---

APP INLA
INLA
True
True
1.0
Distance to R
Distance to R
Figure 1: (A) Simulated data on 9ı (da) with a 0.5 and contours of distance to letter "R". (B) The estimated intensity functions. Solid lines and shaded areas represent predictive medians and [0.025, 0.975] prediction intervals, respectively.
Related Work
Kernel Intensity Estimators (KIEs) are the de facto standard nonparametric approaches for esti- mating point process intensities as functions of covariates in a non-Bayesian manner [177, I& []]- They construct intensity estimators by rescaling well-studied kernel (probability) density estimators for biased sample data. KIES offer the lowest computational complexity as well as sophisticated methods for bandwidth selection. Analogue estimators, denoted by LLEs, have been also proposed based on local likelihood density estimators [5]] KIES and LLEs have a variety of variants, some of which are readily available from library spatstat in R [[]. It is notable that the estimators also provide approximate point-wise wise confidence intervals for estimated intensities, which allows for uncertainty analysis.
Gaussian Cox Processes (GCPs) are the Bayesian alternatives KIEs, where latent intensities are modeled by Gaussian processes or Gaussian random fields through positive link functions [29]. They provide us methodology to perform the interval estimation of intensities as well as the hyper- parammerr optimization in a Bayesian manner. GCPs have been intensively devel- estimating intensities as functions each point on the observation space [M, 2 1X I22 III I1. AR 2TC 36, g2]. being simpler case (y(t) L) of what we consider in this paper. The estimation of intensities as functions general covariates, which varies across observation do- main, has recently been tackled by using domain discretization [2T, 27, 28. GY) 51]] it partitions the observation domain into finite disjoint sub-areas and assumes the intensity to be constant over each sub-area. Although domain discretization is applicable to general GCPs regardless of link/kernel functions, it suffers from poor scaling in terms of domain dimension, as well as being sensitive to discretization size. Fitting GCPs data are often based time-consuming Markov chain Monte Carlo methods, but the integrated nested Laplace approximation (INLA) [211, 22), 37] can alleviate the fitting costs substantially. Fortunately, INLA also accessible through the library INLA in R. INLA and KIE seem to be old methods published in 2008-2012, but to the best of our knowledge, there has been feasible methods proposed to address the problem "estimating intensity function of covariates" other than the methods. Thus, they are SOTA baselines for the problem we tackled.
oped for the problem
Permanental Processes (PPs) are special GCPs that adopt quadratic forms of link functions [[M]]. of which computational advantages have been highlighted by recent researches [TER] 23, 26, [7]] The fact that PPs involve representer theorem when yt) t was exploited through the RKHS theory [13]], the Mercer's theorem [17]] and the path integral formulation [24]] Therefore, our proposed model is generalization of their works for yt)) +
Neural Network- based Point Processes to explorer the dependence of point patterns on covariates have rarely been proposed. An exception is the deep mixture point process [3I], which locates the finite representative points on the observation domain and models the intensity value on the point by neural network model that eats as input the covariates the point. This model provides a feasible estimation algorithm, but cannot perform the interval estimation of intensities.
Experiments
We examined the validity and the potential efficiency of our proposed model by evaluating it against the state-of-the-art Bayesian and non-Bayesian approaches on synthetic and open real-world data. As non-Bayese Bayesian approaches, we adopted KIEs (KIE) and LLEs (LLE) with the ratio (ra), the
7

---

Table 1: Results on two types of synthetic data across 100 trials with standard errors brackets. cpu is the CPU time second, and N is the average data size. The underlines represent the best predictive performances on each metric, and the performances not significantly (p < 10-2) different from the best one under the Mann- Whhnney U test with Holm method [20) are shown in bold.
91 (dn)
a
0.2
(N
232)
0.5
(N = 586)
a
1.0 (N
1180)
Lo2
197s
cpu
1,25
1.5
L975
cpu
Lo25
1.5
1975
cpu
0.15
1.16
0.21
0.05
0.29
2.26
0.36
0.06
0.43
3.45
0,51
0.06
KlEtrs
(0.05)
(0.35)
(0.10)
(0.00)
(0.20)
(0.52)
(0.10)
(0.00)
(0.26)
(0.84)
(0.12)
(0.00)
0.17
1.41
0.23
0.07
0.43
2.91
0.57
0.08
1.11
5.05
0.92
0.08
LLEra
(0.09)
(0.39)
(0.05)
(0.00)
(0.34)
(0.80)
(0.29)
(0.00)
(0.85)
(1.11)
(0.45)
(0.00)
0.27
1.31
0.62
101
0.54
2.54
1.00
108
0.89
3.71
1.51
110
INLA
(0.03)
(0.30)
(0.13)
(6.01)
(0.06)
(0.51)
(0.14)
(5.19)
(0.07)
(0.70)
(0.20)
(6.10)
0.14
1.08
0.18
9.30
0.24
1.93
0.34
14.1
0.38
2.74
0.45
25.5
APP
(0.06)
(0.32)
(0.09)
(0.52)
(0.10)
(0.63)
(0.22)
(0.80)
(0.26)
(0.78)
(0.33)
(1.05)
0.14
1.08
0.18
0.61
0.24
1.93
0.34
2.62
0.38
2.72
0.44
2.81
APPNYS
(0.05)
(0.31)
(0.09)
(0.09)
(0.10)
(0.62)
(0.22)
(0.19)
(0.27)
(0.78)
(0.33)
(0.35)
0.13
1.06
0.18
2.42
0.24
1.90
0.34
2.48
0.37
2.65
0.43
2.56
APPREM
(0.07)
(0.33)
(0.13)
(0.17)
(0.10)
(0.60)
(0.24)
(0.07)
(0.27)
(0.75)
(0.34)
(0.05)
APP
0.17
1.38
0.24
0.30
0.31
2.48
0.45
0.32
0.49
3.66
0.60
0.35
(0.04)
(0.32)
(0.07)
(0.01)
(0.12)
(0.64)
(0.23)
(0.01)
(0.24)
(0.82)
(0.29)
(0.01)
g2(dR)
0.2
(N
322)
0.5
(N 801)
1.0 (N 1610)
Lo25
L.
1,75
cpu
1 Lo2
1975
cpu
1,025
1.s
1975
cpu
0.18
1.47
0.24
0.06
0.38
2.76
0.40
0.06
0.83
4.77
0.64
0.06
KlEtrss
(0.05)
(0.35)
(0.10)
(0.00)
(0.20)
(0.52)
(0.10)
(0.00)
(0.26)
(0.84)
(0.12)
(0.00)
0.24
1.67
0.28
0.07
0.74
3.47
0.54
0.08
1.95
6.17
1.02
0.08
LLEra
(0.09)
(0.39)
(0.05)
(0.00)
(0.34)
(0.80)
(0.29)
(0.00)
(0.85)
(1.11)
(0.45)
(0.00)
0.45
1.73
1.34
112
0.96
3.31
2.22
115
1.68
5.04
3,45
116
INLA
(0.03)
(0.30)
(0.13)
(6.01)
(0.06)
(0.51)
(0.14)
(5.19)
(0.07)
(0.70)
(0.20)
(6.10)
0.15
1.06
0.17
10.1
0.25
1.87
0.27
17.7
0.40
2.75
0.44
40.1
APP
(0.06)
(0.32)
(0.09)
(0.52)
(0.10)
(0.63)
(0.22)
(0.80)
(0.26)
(0.78)
(0.33)
(1.05)
0.15
1.08
0.18
1.02
0.25
1.88
0.27
2.72
0.39
2.79
0,44
3.15
APPNYS
(0.05)
(0.31)
(0.09)
(0.09)
(0.10)
(0.62)
(0.22)
(0.19)
(0.27)
(0.78)
(0.33)
(0.35)
0.14
0.99
0.17
2.40
0.22
1.76
0.26
2.57
0.36
2.53
0.40
2.78
APPREM
(0.07)
(0.33)
(0.13)
(0.17)
(0.10)
(0.60)
(0.24)
(0.07)
(0.27)
(0.75)
(0.34) (0.05)
APP
0.20
1.57
0.28
0.30
0.37
2.76
0.46
0.34
0.57
4.06
0.65
0.40
(0.04)
(0.32)
(0.07)
(0.01)
(0.12)
(0.64)
(0.23)
(0.01)
(0.24)
(0.82)
(0.29)
(0.01)
reweighted (re), and the transformation (tr) forms [55], and selected the kernel bandwidths of KIEs by Silverman's 's rule of thumb [+0]] (SI) and a method by Sheather & Jones [B9] (SJ), resulting in the nine non-Ban -on-Baysiann variants (e.g. K]Etrgs). Some of them are not implemented for multivariate covariates in spatstat, and they were applied only to synthetic data, where Y C IR. As Bayesian ap- proach, we employed the log-Gaussian Cox process with INLA [22]] (INLA) with 100 x 100 domain discretization, where the random walk model of order was used as the prior process, and the hyper- parameters were optimized by maximizing the marginal likelihood through the 16-points grid search. For our proposal, we adopted APPs with the naive approach (APPNai) and with the degenerate ap- proaches of the random Fourier map [XS]] (APPREM) and the Nyström approximation [+9]] (APPNYs)
a We applied to the APPs multiplicative Gaussian kernel, k(y, y') = O I (Oalua =()))²
where the byper-parameter 9 .o,.....6 (D,)) was optimized for each data by maximizing the marginal likelihood (23) through the 25-points grid search. A MacBook Pro with 4-core CPU (2.8
GHz Intel Core 17) was used. For details of the model configurations, see the supplementary material (85).
4.1 Synthetic Data
Following [S], we created 2D data sets (T C R.2 c R) generated from two types of intensity functions: A(t) 91 (dR(t)) = a exp(5 (d)(t)) and A(t) g2(dR(t))= a exp(5 4d (t))), each of which has 100 trial sequences. Here dt((t) denotes the shortest distance from given location t to the set of lines arranged in the shape of the letter "R" (see Figure IA). The coefficient
"Code and data provided at https::/ h/ppst/ttttttbb www/Hiiimiiii....oocc Http/H/iKKm/APP....
8

---

Table 2: Results on real-world data across 10 trials. opu is the CPU time in second, and the under- lines represent the best predictive performances on each metric. Notations follow Table MI
copper Ches
bei Chust
clmf ires
Uesa
cpu
Uust
cpu
Uust
Chese
cpu
5.65
1.64
0.32
2.92
1.31
0.26
2.39
1.30
0.24
KlFragi
(0.37)
(0.09)
(0.01)
(0.00)
(0.07)
(0.00)
(0.00)
(0.08)
(0.01)
5.65
1.64
0.33
2.92
1.31
0.25
2.40
1.38
0.24
KlErest
(0.37)
(0.09)
(0.01)
(0.00)
(0.07)
(0.00)
(0.00)
(0.08)
(0.01)
5.64
1.66
82.9
2.90
1.13
131
2.39
1.18
129
INLA
(0.34)
(0.10)
(3.00)
(0.01)
(0.06)
(3.37)
(0.01)
(0.06)
(4.61)
5.65
1.64
0.60
2.87
1.01
15.8
2.36
1.15
17.9
APPNYs
(0.37)
(0.09)
(0.06)
(0.01)
(0.05)
(0.41)
(0.01)
(0.07)
(1.06)
5.65
1.64
2.12
2.88
1.02
4.64
2.36
1.16
4.28
APPREM
(0.37)
(0.09)
(0.10)
(0.01)
(0.05)
(0.21)
(0.01)
(0.07)
(0.22)
was set 0.2, 0.5, and 1.0. The predictive performance was evaluated based on the integrated p-quantile loss [[8], defined as 1, - Sn ( ( (ola())>0(y) (g-p) dy, where I. ((y), and ( denote the indicator, the predicted p-quantile of the intensity function covariate domain, and the true one, respectively. Here, we adopted Lo25. ( (integrated absolute error) and 1975.
Table displays the predictive performances on various synthetic data, where only the methods with relatively good accuracy are listed due to space limitations (see the supplementary material ($5) for the full results). It shows that our approaches (APPs) outperformed the reference methods across all predictive metrics, where the performance gaps were not marginal (e.g., the improve- ments of L,5 were 9%~455)). Furthermore, APPs with the degenerate approaches (APPNYs and APPREM) were performed several tens of times faster than INLA, feasible Bayesian method based on domain discretization. Table I empirically suggests that our Laplace approximation on the func- tional space (see Section Mell can provide reasonable approximation of the predictive distribution: APPs' good scores regarding l (o22 and (975 demonstrate that the approximation can recover the pos- terior distribution accurately; APP optimized with the approximate marginal likelihood (APPREM) outperformed APP with an initial hyper-parameee approximate marginal likelihood for determining (FPP))M). showing the practical utility of the Figure II implies that INLA tends to overestimate the predictive errors. Also, we ran an additional experiment based on larger data sets and synthetic data sets with 2D covariate (see the supplementary material ($6-$7)).
4.2 Real-world Data
We examined the validity of our approach based on real-world spatial (T R²2) data provided by spatstat. data in R (GPL-3) [M]] copper consists of 67 location points of copper ore deposits and 146 line segments representing faults. The covariate of interest the shortest distance from given location t to the set of faults () C IR); bei consists of locations of 3605 trees of the species Beilschniedia pendula and geo-information in tropical rain forest. The covariates of interest are the terrain elevation and the terrain slope () C RR);; clmf ires consists of locations of forest fires in the Castilla- Mancha region of Spain and the geographical information. We restricted the analysis to events happening within rectangular region (180, 60)- (330, 360), resulting in 4241 location points. The covariates of interest are the terrain elevation and the terrain slope ()² C R²). For each data set, we randomly split the data points into 10 subsets, assigned one to test and the others to training data, and conducted 10-fold cross validations of the predictive performances based on the negative test log likelihood of point patterns (Mess) and the negative test likelihood of counts (cles)). Details of the metrics are given in the supplementary material ($5). Note that to eliminate the effect of data assignment bias, the evaluation multiplied estimated intensity by a factor, NNwsa/Nraain
Table displays the predictive performances on the three real-world data sets, which shows that our approach (APP.) outperformed the reference methods in most cases: On bei, APP achieved the significant 10%~20% improvements of chest against the reference methods; On clmfires, the APP's improvements of predictive performance against INLA were small, but being statistically significant for Uess: On copper, the predictive performances of all the methods were found not to be significantly different from each other, which might be due to the too small data size (N 67).

---

x10
B
15
0.3
0.3-
1.5
0.2
14
1.3
0.1
0.1-
Elevation
Slope
KlEras
INLA
APP
120
140
160
Elevalion
Figure 2: Estimation results on bei. (A) Covariate maps on 2D spatial domain with observed points (black dots). (B) The estimated intensities functions of 2D covariates.
It should be noted that the CPU times of INLA are less affected by the size of data due to domain discretization. Thus INLA remains a valuable estimation approach when the size of data is large and the dimension of observation domain is low. Figure E displays the estimation results on bei, which implies that KIES tend to underestimate the intensity modulation than Bayesian methods. Here, the estimated intensity by INLA seems less smooth than the others, which is due to the additive assumption in the R package INLA: X(t) x explgi (elev((t) g2(slppe(t))]
5 Conclusions
We have proposed novel Bayesian estimation method of point process intensity as function of co- variates. We augmented permanental process such that latent intensities can be defined on covariate space, and showed that representer theorem holds for the augmented process. By exploiting the representer theorem, we have derived fast estimation scheme which scales linearly with data size without using either domain discretization MCMC. Based on synthetic and open real-world data, we confirmed that our method outperforms state-of-the-art methods in terms of predictive accuracy while being substantially faster than conventional Bayesian method. Although the key derivation of the scheme are based on the methodology proposed in [244]. we believe that the result makes a non -on-incremental impact on the Bayesian estimation community for point processes. Also, our find- ing makes non-trivial contribution to the reproducing kernel Hilbert space community because it shows a novel and practically important likelihood function for which the representer theorem holds under the regularization term of the squared Hilbert space norm.
Future work & limitations: This time, optimized the arper-parameter of model by max- imizing the marginal likelihood through grid search, but it is clearly beneficial if gradient descent algorithms could be utilized. Actually, the maximization of the marginal likelihood, which has an argmin (see Section (.2), belongs to bi-level optimization, and we can compute gradients of the marginal likelihood regarding the found by implementing the gradient descent methods.
exactly [[6]. Better hyper-parameters could be
We achieved scalable algorithm under degenerate form of equivalent kernel function, but an equivalent kernel function with a small rank might lead to overconfident predictive covariance [S4]]. A promising approach the problem is employing variational Bayesian (VB) approximations with inducing points, which usually provide a full-rank predictive covariance. VB-based methods for the scenario "estimating intensity function of covariates" have not been proposed so far, and thus it is an important next step of our study.
As in ordinary permanental processes, the nodal line problem could arise in the augmented perma- nental process (APP): the posterior distribution of the latent variable z(-) has many local modes since =x(+)) can Iead to similar intensity functions À z²(.). and artificial zero crossings of zx-) could happen especially on a location where the intensity is row, To alleviate the problem, John and Hensman [23]]] have proposed to extend the quadratic link function to include offset parameter 3, so that k(x)= = (x + 3) which is also valid for our APP.
One promising extension of our model is for the survival analysis methods, which has been actively developed in medical field. In it, observation domain time (one-dimensional) and the observation is finished when an event occurs for the first time, which our model cannot be applied directly due to the dependence of the observation domain on the event data. However, it might be possible to construct feasible algorithm through the similar methodology proposed in our paper.
10

---

References
[1] Ryan Prescott Adams, lain Murray, and David J.C. MacKay. Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities. In International Conference on Machine Learning, pages 9-16, 2009.
[2] Virginia Aglietti, Edwin V. Bonilla, Theodoros Damoulas, and Sally Cripps. Structured vari- ational inference in continuous Cox process models. In Advances Neural Information Pro- cessing Systems 32, 2019.
[3] Kendall Atkinson. A personal perspective on the history of the numerical analysis of Fredholm integral equations of the second kind. In The Birth of Numerical Analysis, pages 53-72. World Scientific, 2010.
[4] Adrian Baddeley and Rolf Turner. Spatstat: an package for analyzing spatial point patterns. Journal of Statistical Software, 12:1-42, 2005,
[5] Adrian Baddeley, Ya-Mei Chang, Yong Song, and Rolf Turner. Nonparametric estimation of the dependence of a spatial point process on spatial covariates. Statistics and Its Interface, (2):221-236, 2012.
[6] Sankarshan Basu and Angelos Dassios. A Cox process with log-normal intensity. Insurance: Mathematics and Economics, 31(2):297-322, 2002.
[7] Russel E Caflisch. Monte Carlo and quasi-Monte Carlo methods. Acta Numerica, 7:1-49, 1998.
[8] David Cox. Regression models and life-tables. Journal of the Royal Statistical Society: Series B (Methodological), 34(2):187-202, 1972.
[9] John Cunningham, M. Yu Byron, Krishna Shenoy, and Maneesh Sahani. Inferring neural firing rates from spike trains using Gaussian processes. In Advances in Neural Information Processing Systems 20, 2007.
[10] Steve Cygu, Jonathan Dushoff, and Benjamin M Bolker. pcoxtime: Penalized Cox proportional hazard model for time-dependent covariates. arXiv preprint arXiv:2 102iv2210..02297, 2021.
[11] Peter J. Diggle, Paula Moraga, Barry Rowlingson, and Benjamin M. Taylor. Spatial and spatio- temporal log- Gag-Gaassinn Cox processes: extending the geostatistical paradigm. Stafistical Sci- ence, 28(4):542-563 2013.
[12] Christian Donner and Manfred Opper. Efficient Bayesian inference of sigmoidal Gaussian Cox processes. Journal of Machine Learning Research, 19:1-34, 2018.
[13] Seth Flaxman, Yee Whye Teh, and Dino Sejdinovic. Poisson intensity estimation with repro- ducing kernels. In Artificial Intelligence and Statistics, pages 270-279, PMLR, 2017.
[14] Anthony C. Gatrell, Trevor C. Bailey, Peter J. Diggle, and Barry S. Rowlingson. Spatial point pattern analysis and its application in geographical epidemiology. Transactions of the Institute of British Geographers, pages 256-274, 1996.
[15] Israel Gohberg, Seymour Goldberg, and Nahum Krupnik. Traces and Determinants of Linear Operators, volume 116. Birkhäuser, 2012.
[16] Stephen Gould, Basura Fernando, Anoop Cherian, Peter Anderson, Rodrigo Santa Cruz, and Edison Guo. On differentiating parameterized argmin and argmax problems with application to bi-level optimization. arXiv preprint arXiv:1607.0544/ 2016.
[17] Yongtao Guan. On consistent nonparametric intensity estimation for inhomogeneous spatial point processes. Journal of the American Statistical Association, 103(483): 2008.
[18] Yongtao Guan and Hansheng Wang- Sufficient dimension reduction for spatial point processes directed by Gaussian random fields. Journal of the Royal Statistical Society: Series (Statis- tical Methodology), 72(3):367-387, 2010.

---

[19] Tom Gunter, Chris Lloyd, Michael A, Osborne, and Stephen J. Roberts. Efficient Bayesian nonparametric modelling of structured point processes. In Uncertainty Artificial Intelli- gence, 2014.
[20] Sture Holm. simple sequentially rejective multiple test procedure. Scandinavian Journal of Statistics, pages 65-70, 1979.
[21] Janine Illian, Sigrunn Sørbye, and Hávard Rue. toolbox for fitting complex spatial point process models using integrated nested Laplace approximation (INLA). The Annals of Applied Statistics, ((4):1499-155) 1530, 2012.
[22] Janine Barbel Illian, Sigrumn Holbek Sørbye, Håvard Rue, and Ditte Katrine Hendrichsen. Using INLA to fit complex point process model with temporally varying effects-a case study. Journal of Environmental Statistics, 3(7), 2012.
[23] S. T. John and James Hensman. Large-scale Cox process inference using variational Fourier features. In International Conference Machine Learning, volume 80, pages 2362-2370. PMLR, 2018.
[24] Hideaki Kim. Fast Bayesian inference for Gaussian Cox processes via path integral formula- tion. In Advances in Neural Information Processing Systems 34, 2021.
[25] Diederik P. Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint arXiv: 142.6990, 2014.
[26] Chris Lloyd, Tom Gunter, Michael Osborne, and Stephen Roberts. Variational inference for Gaussian process modulated Poisson processes. In International Conference on Machine Learning, volume 37, pages 1812-1822. PMLR, 2015.
[27] Peter McCullagh and Jesper Møller. The permanental process. Advances in Applied Probabil- ity, 387 888, 2006.
[28] Jesper Møller and Rasmus Plenge Waagepetersen. Statistical Inference and Simulation for Spatial Point Processes. CRC Press, 2003.
[29] Jesper Møller, Anne Randi Syversveen, and Rasmus Plenge Waagepetersen. Log Gaussian Cox processes. Scandinavian Journal of Statistics, 25(3):451-485 1998.
[30] E. J. Nyström. Über die praktische Auflösung von Integralgleichungen mit Anwendungen auf Randwertaufgaben. Acta Mathematica, 54:185 204, 1930.
[31] Maya Okawa, Tomoharu Iwata, Takeshi Kurashima, Yusuke Tanaka, Hiroyuki Toda, and Naonori Ueda. Deep mixture point processes: Spatio-temporal event prediction with rich contextual information. In ACM SIGKDD International Conference Knowledge Discovery & Data Mining, pages 373-383, 2019.
[32] Andrei D. Polyanin and Alexander Manzhirov. Handbook Integral Equations. CRC press, 1998.
[33] Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. In Ad- vances in Neural Information Processing Systems 20, 2007.
[34] Carl Edward Rasmussen and Joaquin Quinoñero Quineñero-Candela. Healing the relevance vector ma- chine through augmentation. In International Conference on Machine Learning, pages 689- 696. PMLR, 2005.
[35] Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. MIT press, 2006.
[36] Stephen L. Rathbun and Noel Cressie. Asymptotic properties of estimators for the parameters of spatial inhomogeneous Poisson point processes. Advances Applied Probability, 26(1), 122-154, 1994.
12

---

[37] Hävard Rue, Sara Martino, and Nicolas Chopin. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 7122)319-3392, 2009.
[38] Matthias Seeger, David Salinas, and Valentin Flunkert. Bayesian intermittent demand fore- casting for large inventories. In Advances Neurul Information Prucessing Systems 29, 2016.
[39] Simon Sheather and Michael C. Jones. A reliable data-based bandwidth selection method for kernel density estimation. Journal of the Royal Statistical Society: Series B (Methodological), 53(3):683 690, 1991.
[40] Bernard W. Silverman. Density Estimation for Statistics and Data Analysis. Routledge, 2018.
[41] Dietrich Stoyan and Antti Penttinen. Recent applications of point process methods in forestry statisties. Statistical Science, pages 61-78, 2000.
[42] Ming Teng, Farouk Nathoo, and Timothy D. Johnson. Bayesian computation for log-Gaussian Cox processes: A comparative analysis of methods. Journal of Statistical Computation and Simulation, $7:2227-2252. 2017.
[43] Terry Therneau, Cindy Crowson, and Elizabeth Atkinson. Using time dependent covariates and time dependent coefficients in the Cox model. Survival Vignettes, 2:3, 2017.
[44] Andrea Vedaldi and Andrew Zisserman. Efficient additive kernels via explicit feature maps.
IEEE Transactions on Pattern Analysis and Machine Intelligence, 34(3):480-492, 2012.
[45] Sreekanth Vempati, Andrea Vedaldi, Andrew Zisserman, and C. V. Jawahar. Generalized RBF feature maps for efficient detection. In BMVC, pages 1-11, 2010.
[46] Grace Wahba. Spline Models for Observational Data, volume 59. SIAM, 1990.
[47] Christian J. Walder and Adrian N. Bishop. Fast Bayesian intensity estimation for the permanen- tal process. International Conference on Machine Learning, volume 70, pages 3579-3588.
PMLR, 2017.
[48] Christopher Williams and Matthias Seeger. The effect of the input density distribution on kernel-based classifiers. In International Conference on Machine Learning, pages 1159-1166. Citeseer, 2000.
[49] Christopher Williams and Matthias Seeger. Using the Nyström method to speed up kernel machines. In Advances in Neural Information Processing Systems 13, 2000.
[50] Ganggang Xu, Rasmus Waagepetersen, and Yongtao Guan. Stochastic quasi-likelihood for case-control point pattern data. Journal of the American Statistical Association, 114(526): 631-644, 2019.
[51] Yu Ryan Yue and Ji Meng Loh. Bayesian semiparametric intensity estimation for inhomoge- neous spatial point processes. Biometrics, 67(3):937-946, 2011.
Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes]
(c) Did you discuss any potential negative societal impacts of your work? [Yes] See supplementary materials
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes] See Section D
13

---

(b) Did you include complete proofs of all theoretical results? [Yes] Also see supplemen- tary materials
3. If you experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material or as a URL)? [Yes] Full instruction of experiments
in supplementary materials, and code is also submitted.
(b) Did you specify all the training details (e.g., data splits, hyper-parameters, how they were chosen)? [Yes] See Section E and supplementary materials
(c) Did you report enror bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes] See Section
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [Yes] See Section !
4. If you using existing assets (e.g., code, data, models) curating/releasing new assets...
(a) If your work uses existing assets, did you cite the creators? [Yes] See Section
(b) Did you mention the license of the assets? [Yes] See Section !!
(c) Did you include any assets either in the supplemental material or URL? [No]
(d) Did you discuss whether and how consent was obtained from people whose data you're 'r using/curating? [No]
(e) Did you discuss whether the data you are using/curating contains personally identifi- able information or offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, applicable? [N/A]
(b) Did you describe any potential participant risks, with links Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent on participant compensation? [N/A]
14


Reviews:
[{"review_id": "kfhD2c1POAO", "reviewer": "Reviewer_V1pm", "summary": "This paper proposes an approach for making Maximum A Posteriori (MAP) inference on the intensity function of a Poisson process, regarded as a function of covariates. The representer theorem, which the authors show applies when a Gaussian Process prior is placed on the square root of the intensity function, is used to achieve a complexity linear in the sample size.", "questions": "\nWhy do we need Augmented Permanental Processes? Why not use Permanental Processes in the covariates space?", "limitations": "N/A", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "\n**It is unclear why this work is needed. Do we really need Augmented Permanental Processes?**\n\nThere is nothing in the theory of Poisson processes that requires the input space to be time, space or space-time. \n\nIf the authors want to infer the intensity function as a function of covariates, they could simply define the space of points to be the space of covariates, namely $\\mathcal{Y}$. \n\nThe Poisson likelihood would simply read $p(\\mathcal{D} \\vert  x(\\textbf{y}))= \\sum_{i=1}^n \\log \\lambda(\\textbf{y}_i) - \\int_\\mathcal{Y} \\lambda(\\textbf{y}) d\\textbf{y}, ~~~ \\lambda(\\textbf{y}) := \\kappa(f(\\textbf{y})),$ and the prior would read $f \\sim \\mathcal{GP}(0, k).$\n\nFor all practical purposes, this would work even if the mapping function $g: \\textbf{t} \\to \\textbf{y}$ is not injective.\n\nWhich begs the question: what contribution does this work make to the literature on Permanental Processes [14, 24, 46]?\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "cfmUwfHeEz3", "reviewer": "Reviewer_wjWf", "summary": "This paper proposed an augmented permanental process whose intensity is modeled by a Gaussian process on covariate space, and derive a efficient Bayesian inference method without the need of discretizing the domain. The proposed algorithm is based on path integration and Laplace approximation, which facilitates the computation of predictive distribution (uncertainty) and marginal distribution (hyper-parameter). Also the paper designed different methods, including naive methods, Nyström approximation and random features, to construct the equivalent kernel. ", "questions": "I suggest to provide more experimental comparison results with SOTA baseline models, e.g., sigmoidal Cox process [2], permanental process [3]. Provide more analysis expeirmentally or theoretically on the ' Nodal Lines' problem which is a common issue in permanental processes. \n\n\n[2] Efficient Bayesian Inference of Sigmoidal Gaussian Cox Processes. JMLR 2018. \n[3] Large-Scale Cox Process Inference using Variational Fourier Features. ICML 2018. ", "limitations": "The work does not have potential negative societal impact. \nThe authors did not mention the limitations of their work if I did not miss any important points. ", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "strengths:\nThe paper is writting very well and easy to follow. The methodology is solid and seems correct to me but I did not check all the derivations carefully. \n\nweaknesses:\nHowever, I have two concerns about the current version, which makes me hasitate on the rating of this submission. \nThe first is the experimental comparision with baseline models is rather weak. The proposed methods (NAI/NYS/RFM) are only compared with very old works. For example, KIE and INLA are published in 2012. I did not find any SOTA baselines in this work, which cannot convince me of the merits of this work. \nThe second is the contribution of this work is rather incremental. The key components of the proposed inference method all come from previous works, e.g., the MAP estimate using path integral is from [1], let alone the Laplace approximation of posterior and Nyström or random features approximation of kernel. \n\n[1] Fast Bayesian inference for Gaussian Cox processes via path integral formulation. Neurips 2021.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "NtcbmM6Qql", "reviewer": "Reviewer_oHGh", "summary": "This paper studies a fast Bayesian inference of point process intensity as a function of covariate.\n- Using the path integral formulation of Gaussian processes, the maximum a posteriori (MAP) estimate can be obtained by minimizing the action integral, which results in an optimization problem with O(N^2) time complexity for each evaluation. \n- In particular, as the quadratic link function is used in this context, the Nystroem method, and degenerate kernel method are employed to solve integral equations.", "questions": "Can you provide the running time comparisons of your method against the conventional Bayesian method, such as INLA, for the final experiments? ", "limitations": "The authors claim that the proposed method may harm people's privacy by predicting their fine temporal patterns.\nI expect the authors to present more in-depth analysis about the limitations of the method.", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Originality:\n-The novelty of this paper lies in the nonlinear dependence of the latent intensity on covariates, \\lambda( K( y(t) ) ), compared to the previous work [Hideaki Kim 2021], in which the latent intensity is given by \\lambda( K( t ) ).\n\n- In particular, this paper employs quadratic link function, and thus enables to obtain a linear integral equation, which is first solved using Nystroem method, and degenerate kernel by this paper.\n\n- I consider that most results in deriving MAP estimate, marginal likelihood and predictive distribution in Sec 2.2 and 2.3 are a bit incremental, given the results by [Hideaki Kim 2021]. \n\nQuality:\nThe work appears to be sound.\n\nClarity:\nOverall, the methodology is well documented, and the main derivations are easy to follow.\n\nSignificance:\nThe work comes up a faster Bayesian inference method for covariates-dependent point processes, compared to previous MCMC and VB methods using domain discretizations. Thus, I would think it will make some impacts on the Bayesian inference community for point processes although the novelty is not that strong.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ["I don’t know"]}, "extra_text": {}}, {"review_id": "J0vEthBQxIU", "reviewer": "Reviewer_bJU1", "summary": "This paper proposes an efficient Beyian algorithm to estimate the intensity function as a function of Covariates in a Reproducing Hilbert Kernel Space (RHKS). The purposed algorithm is based on an augmented permanent process and enjoys significant computational advantages. The effectiveness of the proposed algorithm is demonstrated through numerical experiments.", "questions": "The proposed approach is essentially an extension of Flaxman et al (2017), which focuses on estimating intensity function as a function of time/spatial locations. But I think the extension is meaningful and solves an important problem. My only concern is about the computation of the Equivalent Kernel in Section 2.4. Because of the change of support from t to y, the integral equation (9) indeed can be computationally intensive to solve when the dimension of covariates is high. This somehow offsets some benefits of kernel ridge regression in RHKS, which can be used to fit high-dimensional functions using kernels such as the multiplicative Gaussian kernel. The algorithm in 2.4 essentially proposes to approximate the integral equation (9) through J-point numerical integration. However, when the dimension of y is large, the required J would grow exponentially, and the proposed approach cannot be applied to the multiplicative Gaussian kernel. One possible solution is perhaps to use a stochastic approximation of (9) in a similar fashion to those in  [1]. It will be helpful to give some discussion on this issue.\n\nIn addition, an important class of nonparametric models in Statistics is the additive model. So if the reproducing kernel K(.,.) is additive, would the Equivalent Kernel also be additive? If this is the case, the computation of the equivalent kernel can be greatly reduced for the additive model, which would be a great contribution as well.\n\nThe simulation study only covers 1-D functions. It would be better if at least 2-D functions can be studied, just to be consistent with the real data analysis.\n\nFinally, the title seems to emphasize \"inference\" using the proposed method. But I do not see any inferential results in numerical studies. What kind of inferences can be carried? Can you give some examples?\n\n[1] Xu, G., Waagepetersen, R. and Guan, Y., 2019. Stochastic quasi-likelihood for case-control point pattern data. Journal of the American Statistical Association, 114(526), pp.631-644.", "limitations": "none", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: the paper is well written, the presentation is clear,  the motivations are valid and the proposed approach is novel and effective.\nWeaknesses: there might still be computational bottlenecks when the number of covariates is large; the synthetic data experiment can be improved.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "8Vz-T75jb1", "reviewer": "Reviewer_Gqug", "summary": "The paper proposes a new inference method for a variant of a Gaussian Cox process, which the authors call augmented permanental process (APP). The model is described as a point process that uses a Gaussian process on the covariate space with an intensity function based on a square link function. The main focus of the paper is to develop a tractable Bayesian inference to estimate the posterior distribution of the latent intensity function that is based on a novel application of the represent theorem. In the next step the predictive distribution is approximated using a Laplace approximation of the posterior. The method is evaluated on a synthetic and three real-world datasets.\n\nUpdate after rebuttal: increased score to 5.", "questions": "Could you provide a discussion of the accuracy of your approximation methods? How much does the Laplace approximation affect the results? What is the effect of approximating the kernel?\n\nFor me the discussion of the \"degenerate approach\" came a bit late. In the section before it was always given that we have degenerate kernel. It would be helpful to emphasize in the beginning that we can obtain it by different approximation techniques.\n\nHere you could also discuss the connection to inducing point GPs (which are related to the Nystrom approximation)?", "limitations": "The authors discussed the limitations.", "rating": 5, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "*Strengths:*\nThe paper is well written and the math seems to be sound.\n\nThe paper develops interesting theory deriving a representer theorem for the latent intensity function. Using this connection the authors obtain a MAP estimator that scales linearly in the number of datapoint. The theoretical insights could be used beyond the scope of the paper.\n\n*Weaknesses:*\nThe empirical evaluation of the method is not overwhelming. \n\nFirst, on all four considered datasets the benefits of using the new method seem to be rather marginal. The performance gaps to the reference methods are rather small. For instance, when comparing KIE_ra_SI and APP_NYS on the real-world datasets we get the log-likelihood values  5.56 vs 5.56, 2.92 vs. 2.87, and 2.39 vs. 2.36. It is not clear to me if these improvement are of any practical value. If the authors belief so, could you find another metric or experiment demonstrating the benefit. If the benefit is speed (at least in comparison to INLA), it would be helpful demonstrating the benefit clearer. For instance, you could show a plot with dataset size on the x-axis and inference speed/prediction performance on the y-axis.\n\nSecond, the datasets seem rather small. It would be interesting to see how those algorithms scale with increased dataset size (see above). Is this an area where you method can shine? For instance, due to the linear scaling can we use APP on huge datasets where competing methods fail?\n\nAs stated above the paper makes interesting theoretical contributions which is focused on the proposed model. If not demonstrating clear empirical benefits, another way where the paper could have impact would be to extend the theory and to show how similar ideas could be applied beyond the scope of the specific model.\n\nThe presentation is quite dense and might be improved by some visual guiding, e.g., a model diagram.\n\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept