Fast Bayesian Inference for Gaussian Cox Processes via Path Integral Formulation
Authors: Hideaki Kim
Venue: NeurIPS 2021 Spotlight (2021)
Keywords: Gaussian Cox processes
point processes
time series analysis
path integral
Gaussian processes

Abstract: Gaussian Cox processes are widely-used point process models that use a Gaussian process to describe the Bayesian a priori uncertainty present in latent intensity functions. In this paper, we propose a novel Bayesian inference scheme for Gaussian Cox processes by exploiting a conceptually-intuitive {¥it path integral} formulation. The proposed scheme does not rely on domain discretization, scales linearly with the number of observed events, has a lower complexity than the state-of-the-art variational Bayesian schemes with respect to the number of inducing points, and is applicable to a wide range of Gaussian Cox processes with various types of link functions. Our scheme is especially beneficial under the multi-dimensional input setting, where the number of inducing points tends to be large. We evaluate our scheme on synthetic and real-world data, and show that it achieves comparable predictive accuracy while being tens of times faster than reference methods.

Fast Bayesian Inference for Gaussian Cox Processes via Path Integral Formulation
Hideaki Kim
NTT Human Informatics Laboratories
NTT Corporation
hideaki .iin. cn@hco. catiicotttt.ooo ct..co...
Abstract
Gaussian Cox processes are widely-used point process models that use Gaussian process to describe the Bayesian priori uncertainty present in latent intensity functions. In this paper, we propose novel Bayesian inference scheme for Gaus- sian Cox processes by exploiting a conceptually-intuiie tion. The proposed scheme does not rely on domain discretization, scales linearly with the number of observed events, has lower complexity than the state-of-the- art variational Bayesian schemes with respect to the number of inducing points, and applicable
path integral formula-
wide range of Gaussian Cox processes with various types of link functions. Our scheme is especially beneficial under the multi-dimensional input setting, where the number of inducing points tends to large. We evaluate our scheme on synthetic and real-world data, and show that it achieves comparable predictive accuracy while being tens of times faster than reference methods.
1 Introduction
Gaussian Cox processes constitute a class of doubly-stochastic point process models which flexible prior over a latent intensity function is established by Gaussian process through positive link function under which point events are generated. Gaussian Cox processes are the gold stan- dard in analyzing event data in Bayesian manner, and have wide ranging list of applications in neuroscience [ТКС], finance [5]], and spatio-temporal analysis [6U]]
The observation process, i.e., the likelihood function of the Gaussian Cox process, depends on the functional form of the latent intensity function over a compact domain, which makes Bayesian inference challenging as it imposes the intractable computation of infiniteddimensional distribu- tions. Many algorithms have been proposed to deal with this difficulty. They include Markov Chain Monte Carlo (MCMC) sampling, domain discretization, and variational Bayesian (VB) approxima- tion. MCMC approaches [2. 22 provide exact inferencing for a specific link function, but they have excessive computation costs. Although domain discretization || [K] achieves sub-linear computa- tional scaling with discretization size, it suffers from poor scaling with data size and the dimension of the domain. To overcome these limitations, state-of-the-un art algorithms employ sparse VB ap- proximation with inducing points [G. BII;; they achieve improved scaling with the dimension of the domain as well as computation cost that is linear to data size. However, the computation costs of VB approaches scale at the rate of O(N + L3) where L and are the number of inducing points and data points, respectively, making them problematic in the multi-dimensional domain setting, where tends to be large. Also, VB approaches achieve good approximations by aggressively focusing on the unique structures of individual kernel functions link functions, and thus are not applicable to a wide range of Gaussian Cox processes.
In this paper, we introduce a cnneeptuall--intitive path integral formulation of Gaussian processes, by which we derive a novel inference scheme for tackling intractable Gaussian Cox processes. The
35th Conference on Neural Information Processing Systems (NcurIPS 2021).

---

formulation predicts that the maximum posteriori (MAP) estimation of a latent intensity function coincides with the problem of finding the most likely path that maximizes the functional posterior, whose solution can be obtained by solving an explicit integral equation. The approximate predictive covariance is also acquired by solving the corresponding integral equation. We propose an efficient algorithm to solve the integral equations; it offers O( + T.²) computation scaling, that is, lower complexity than the state-of-the-arr variational Bayesian schemes with respect to the number of inducing points, as well as linear computation scaling with data size. Furthermore, our inference scheme via path integral holds for any kernel function and a wide range of link functions, and thus provides a practical methodology against Gaussian Cox processes that conventional approaches have yet to tackle successfully.
In Section 2, we introduce the path integral formulation of Gaussian processes, and construct a new Bayesian inference scheme for Gaussian Cox processes. In Section EL we provide practical methodology to perform inferencing. In Section we outline related work. In Section ES, we perform comparative evaluations against reference methods on synthetic and real-world data, and confirm that our scheme achieves comparable predictive accuracy with substantial speed improvements over state-of-the-art algorithms. Finally, Section B states our conclusions.
2 Path Integral Formulation
2.1 Gaussian Cox Process
We assume that latent function on compact space, T C RD, Gaussian Process (GP), denoted by (P(x(t)) /u, kit,, f')), and observations are generated from a pro- cess under latent function xtt), where j and are the mean value and the kernel/covariance function for (t), respectively. Given data set D observed within T, we consider the problem of obtaining the maximum a posteriori (MAP) estimator of (t)) that maximizes the posterior probability
RR is generated from
p (t)) (P(x(t))p, k) p(x(t))D) f (x(() p(D) (t)) (P(x(t)u,kk)
(1)
where p(D)xxt)) is the likelihood function, and f 9r(() in the denominator represents the integral over the function or the diinite-dimensionll variable (t).
The GP prior, (P())tt) (u, k)), be regarded as normal distribution for an variable (see Section (.2)., and thus if the likelihood function depends only on a set of latent functions
on finite N training points, p(D]x(tt) = p(D ((t.)))))=)).. then MAP estimation for an infinite- dimensional variable reduces to that for finite N-dimensional variable by marginalizing out the variables on the points other than {=x(tn)) } (z=(t.))_\==i this inference is tractable and has been examined extensively (e.g. GP regression, GP classification). Here, we focus on a more challenging situation where the likelihood function includes a factor that depends on an integral of latent function over T as follows,
N
- / A(t)dt, A(t) = <(((t))
log pggp(D)aa (()()(t)) = Σ log A(tn) ==1
(2)
where D = [tt is the set of point events occurring in the observation region T, k(x): : R + R4 is a non-negative function called the link function, and >(t) is the intensity function or the instantaneous probability of events occurring at each point in T. The probabilistic model defined by (I-22 is called the Gaussian Cox procexs.
Bayesian inference (I) on situation (2) is often referred as "douublyinntrcttbbl" because it requires solving the integral of the stochastic function ætt) well as computing the intractable posterior probability of an infinitedimmensional variable. Below, we address the problem directly by introduc- ing path integral, a mathematical tool familiar in the field of quantum physics, treat probability distributions of stochastic functions properly and intuitively.
2.2 Gaussian Process Prior via Path Integral
In order to derive an explicit description of GP prior, we first review well-known formulation of GP [Z], and reconstruct its explicit form of probability density distribution by using path integral representation.

---

Consider model defined in terms of a linear combination of M fixed basis functions, fóm (t) pm such that
M
x(t) = + Σ
am (t) p J (tt), tET,
(3)
m=1
where a A (a1,... a'm) T is the M-dimensional weight vector, ((t) & (t), (t)),
and μ is a mean parameter. If a is generated from an M-dimensionnl Gaussian with diagonal covariance matrix of the form
p(a) = -
(4)
the joint probability distribution for any finite (even infinite) set of variables, a - (a(t1),
axt.))) for J S oc, follows the J-dimensional Gaussian given by
p(æ)dæ= = -N (11,,1, K)dæ, K3 = k(),, =c (()) (())))
(5)
where k(t,t) t') is semi-definite kernel function, is the precision parameter, S j,J S J, and 1 (1.....1) T)² Because (5)) holds for arbitrary numbers of dimensions, J, it provides an implicit representation of GP prior.
Next, based on the expression (D), we derive a more explicit representation of GP prior. For sim- plicity of explanation. we assume that the domain is one-dimensional, = and T [0, 7]. and mean yu is zero. Let IC be the integral operator corresponding to k(t, E'), and KC be the inverse operator of K.,
(())(()) /, (1,)
K*K((,8) = Sf k (t, (1,1))((t' (t))) = 5(1 8).
(6)
where 8(.) represents the Dirac delta function. Note that 1C" and k* (t, should be described by a differential operator because relation (①) indicates that k(t, s) is the Green function [TS]] for the operator KC, but we do not go into detail here. Under the finite scenario (α) with [.o [0, 2] and J > 1, the second equation in (8) can be approximated by A' A-1, where A jj' (t, - 153r) Kig kn (ty.ty)) tyy, and 658' is the Kronecker delta. Then we can write down the J-dimensional Gaussian (5) by using K as
AK* A
p(a))dx
exp - (Ax) K" (*(Ax) dar
(2T)
By taking the limit of the division number + oa (A + 0) for the dimee-dimensional Gaus- sian measure (Z), achieve the following path integral representation of GP prior over the latent function zt)):
/ TXT
GP(z(t)O,k)2 & lim p(z) 1
II dx(t,) =
exp
-
(t, s))(()(()))) 2x, (8)
KC|
j=1
where ||| A lim path integral holds,
KA), 2x A lim I1/ V(t, (tt-ty-1)/10 (2T (tx and the following
[xp[[2 1-11-/111 (t, (t,ss)) ()a(()=(s)dddd) gr = VKk. TXT
exp
[A]] represents the Fredholm determinant or functional determinant [T>]]]]]]] of the integral operator KC, which is defined by the product of its cigenvalues. See the supplementary material ($1) for the derivation. The expression for GP(((tt)) [1c, k) is easily recovered by sx(() + (() - r and / + in (
The path integral representation (8-M) has a clear advantage over the standard one (5) in that the distribution of the latent function z(t) is written in terms of the explicit integral of z(t) over the domain T, as with the likelihood function (2), which makes it possible to apply functional analyses, calculus of variation, to the Bayesian estimation (I)) To the best of our knowledge, this is the first proposal of path integral representation for GP prior in the machine learning community.
It should be noted that our derivation of the path integral representation is based on an intuitive view reminiscent Feynman's 's path integral formulation of the quantum field theory [[Z)), and thus not mathematically rigorous. In this paper, we confirm that our path integral expression is valid by using experiments rather than taking more mathematically rigorous approach; we discuss this briefly in Section M!
3

---

2.3 Maximum A Posteriori Estimation
Using the representation (8) and the likelihood function (2), functional.
we can rewrite the posterior (II as the
p(a(t))))z
P(D) exp x((() Ør.
(10)
where ((x(())z(()) is an action integral, defined by
((z(t),z(t)) / ( - - (z((t) + +x(x(t))
1 log (z()))))))( d log [C], (11)
and st)) - (L)) p) K p)).
(IT indicates that MAP estimation of sxt)
equivalent to the problem of finding the most likely function path z(t) that minimizes the
action integral. Thus we now apply calculus of variations to the action integral, where the functional
derivative of (ax (t), .t)) on the MAP estimator & (t) equal to zero: 82(2) dá (t)+ 68 5t(() - 0.
82(0)
This leads us to realize the following equation for deriving the MAP estimator 2(t),
a(t) + / kett, k(t,,t)) ())(((t = +
N
((((tT)) <(((tx)))
Σ ket, tn
TET,
(12)
where k(x)) A /dr())/dr. See the supplementary material ($2) for the detailed derivations (12))
2.4 Predictive Covariance
One of the advantages of GP models over non-Bayesian approaches is that they can provide predic- tive distributions. We apply a Laplace approximation to GP models (II), and find the approximate form of the predictive covariance.
We now know the mode of the posterior, it)), and consider Taylor expansion of functional action potential S (()))..z(() centered the mode such that
1 s(x(t),x(t)) = S(f(t), 2(t)) + 2 / o'
2.
((s))))dss
(13)
TXT
where o² (tt, s) (22(2)))
bz(t)om(() l== is the second derivative of S. The first term in the Taylor expansion vanishes due to the stationary condition. The quadratic approximation of the action integral corre- sponds the approximation of the posterior process by GP, and the predictive covariance or the kernel function for the posterior GP, denoted by at, s), can be obtained by the functional inversion of a² (t, 8) (see Equation (8),, which results in
(t, h(t, 8) h(t) (Z where Zyn'' 8an (x)) (()))x(()(CC) {(=)((())))))))))))) >=)(()))"' k(x)
H) h(ss,
(14)
68 (a)
dx((x) doa
H & hitn, tw))
lxtt)
-
(h(t, t1).... ,,,,tttt (x)) and ltt,, 8) is defined by the following Fredholm integral equation of the second kind [36]].
s(,,s) + [
s)dt' k(t,s) s).
(15)
Note that Zun' has infinite value for k(x) e", which leads to the relation, at,, 8) hit, The full derivation (14-LS) is given in the supplementary material ($3).
2.5 Marginal Likelihood
Let Σ and H be the integral operators for (t, s) and h(t, s), respectively: Σ A fr (o)t,'))dt' H A +h(,tt'd Under Laplace approximation (I3) we can obtain the marginal likelihood evidence of Gaussian Cox processes, p(D), by performing the following path integral (
log p(D) = log
expl -apl-S(ann S (x(t) Zr = -S(((t) (())) + log Σ],
(16)

---

Table 1: Link functions k(z) and their derivation.
k(x))
k(x))
i(x) x
7(6) = 1s/16
quadratic exponential
x²
2x
2
4/h
exp(x)
exp(x)
exp(x)
softplus
log(1+expp())
(1+expl- (x)xx x))
exp(-x))(1 expl (x))2
- -log(1 A)
where log
| log [ - log IIN + 1--zz1110
HI, and and
N
(((((t)) ((((tt))))
((a(t),((t) log [K| + Σ
(b(t))
log logx(((t.))
f ((((t)
1
(3(t)-)(((((()
(17)
The full derivation is given in the supplementary material ($4).
3
How to Solve Integral Equations
To obtain the MAP estimator at) and the predictive covariance at, s). we need to solve the corre- sponding integral equations. We provide efficient algorithms to solve them (also see Suppl. ($8)).
3.1
MAP Estimator
Equation (I22) for the MAP estimator is nonlinear integral equation of the second kind, sometimes called the Hammerstein integral equation [23]. We can obtain approximation solution of the equation through the most widely used class of methods, projection methods [H]], which approximate
the solution its derivation by choosing an approximation from given finite L-dimensional linear
subspace of functions, denoted by Z. We employ the collocation method [29]] variant of the projection method:
I 1i11 arg min Σ (((p)). (BiHE)
A(s(t)) Σ Bipi(t),
(18)
(1)) A
Σ ket,, kn)7
Σ Bispi (tn))
Σ
Br f kit, (') +w) Σ Scp((!).
( is basis for 2, {p ET) are the collocation points, and (K) represents k(z)/n({ being re-written as function of k(x) (see Table M.. As the basis {5p: H= we adopt the cigenfunctions of the kernel operator K. which reduces the problematic integral term in the residual into a tractable one as fy k(t, (')'''t = (t)) where {Ar) is the set of the eigenvalues of I. The procedure for finding the eigenfunctions provided at the end of this section. For fair comparison with VB-based approaches of Gaussian Cox processes [GL BT], we solve the minimization problem of the sum of the squared residuals by using popular gradient descent algorithm, Adam [2.7]. We can obtain the MAP estimator i(t) by substituting (LE) into (I2): att) = Σ,, k(t, tn 1 (E, (tpr(tx)) E, (1))
where r(t) is the residual derived from (12),
The sparsely located collocation points, (mr) reduce the computational complexity of MAP estimation substantially, and play a similar role in the Bayesian inference scheme to that played by the inducing points in the variational framework for sparse GP regression [В4]] In this paper, adopt a regular grid over T as the collocation points in accordance with the VB-based references BL 3, but note that other location choices are possible. For simplicity of explanation, we hereafter call (pi Hp inducing points.
It should be noted here that our approach (18) with the collocation method can be applied to cases where >() = {())/<(a) defined properly as function of A(x). Table MI shows that 7(k) can be defined for popular link functions which include quadratic [[T [T], exponential [02, 30 33) and
softplus [28] functions. But our approach is not applicable to link functions whose derivatives, A(x)),
not monotonic with respect to æ, such as sigmoidal link functions [2) R. 22).
5

---

3.2 Predictive Covariance
Given the MAP estimator ((),, the predictive covariance (L) can be obtained by solving the linear integral equation (LS) about fit s). We apply the Gelerkin method [2]] variant of the projection method, to solve the equation, which results in
h(t, 8) Σ (t))
N WI + XE E
- f (z
(19)
See the supplementary material ($5) for the detailed derivation. We estimate the integral in {E via the Monte Carlo method, but the quadratic link function (x<(z) x²)) an exceptional case in
which the integration can be performed analytically. E, 2, due to the orthogonal relation of the
eigenfunctions, J7 (pr(())ppr (t)dt ōu. It should be noted that (I99 is the Mercer expansion [6]]
of function (tt, s), where {wi) is set of its eigenvalues.
3.3 Marginal Likelihood
We can evaluate the marginal likelihood by substituting (E) into (IHTT),
Σ log K (x(tn)) - - -f ((((t))
p(g) =
Σ Sn
(20)
1 Σ xx((t)) x(((t)))
((((t)))
Σ ((Brpi(th) (((tt))-)) (ft(t)t))))()))))))))))))))))))) 2
1 log II Z-
1
[C] 7H
log 2
where the second term, fy dttt))d is estimated by the Monte Carlo method. The last term associated with the functional determinant (see the supplementary material ($ 1)), which can be [C] calculated by the product of its eigenvalues as, log [ = E'= log dL = Σ 1og(1+AAE)) A/))), where wi we approximate the kernel function in terms of the Mercer expansion up to the top L eigenvalues, ktt, 8) Σ (t)r
3.4 Computational Complexity
The objective function to be minimized in MAP estimation is [t=]r())) which needs the com- putation of O(N L + L2) for each evaluation (see Eq. (L the computational complexities of fmm B [E (t, H and (En, k(Pi, ) 1 are both ((L)), and {E, AuBroi ) or EE (rspr(p:))/ } needs the computation of O(L²). When gradient descent algorithm with automatic diffcrentiation [6i] is employed, the computational complexity of the MAP estimation per gradient descent step is equal to O(N L + L²), which is one-Lth the complexity of the VB-based approaches [3 G] We adopt the algorithm in Section H
Furthermore, in view of the shape of the objective function, we can substantially improve the compu- tational efficiency the MAP estimation by using a mini-batch gradient descent (MGD) algorithm: MGD reduces the computational complexity to O(NL), which is due to [7m)W We verify the effectiveness of MGD in the supplementary material ($10). Also, when {pr) are located regu- larly over T and k(t, t') is separable across the dimensions, matrix & defined as Drr = (pr) has the Kronecker structure (see the the supplementary material ($7)), which can further reduce the complexity to O(N L) under (batch) gradient descent. Note that the cost of the VB-based methods [LI [ become ON when the Kronecker structure the gram matrix is exploited.
Once the MAP estimator or the set of coefficient {B)) is obtained, estimating the predictive co- variance (LE) needs the computation of O((N + L) min(N2 I.2) + Nm. where Nme is the number of samples for the Monte Carlo method. The cost O(Nme comes from the computa- tion of h(t, s) (II). and the rest is due to the matrix operation. Due to the fact that h(t, s) is degenerate kernel of rank L, the N N matrix H can be decomposed into a product of N L (tx)) and its transpose as H RR The matrix opera- tion h(t) (Z + RR (h(s) costs O(LN2 + N3) in a direct manner, but O(L2N + L3) if the Woodbury matrix identity is used. Because the predictive covariance computation is not part of the iterative optimization in MAP estimation, its computation cost is negligible compared to that of MAP estimation. Also, the selection of syper-parameter can be performed based on the marginal likelihood (20). which needs the computation of O(N + L) min(N22 22 + Nme).
matrix R, defined by R2

---

3.5 Eigenfunctions of Kernel Operator
The pairs of eigenvalues and eigenfunctions, {^r, 9:) for the kernel operator I. are obtained by solving a homogeneous Fredholm integral equation, J k(t, (s) = (tp)(t), where the set of eigenvalues is at most countable [[6]]. Here we assume a multiplicative kernel, k(t, s) = s(d)), and hyper-rectangulas domain with interval T(4) = [0, 7(4) in each dimen- sion d. The multidimennional integral equation can then be reduced to a set of unidimensional ones, (d) (
I
UTr,, (t((), s) (s)ds = Ai Pis (t()))a. via the separation of variables, which results in multiplicatiie form of solution, (Ae; (t()t) (Ia Aig
(d) II Pia (d) )). where (X .4%% is the la-th solution of the integral equation in dimension d, and 1 (41...,(p)) Ip). Thus the computa- tion complexity of solving the integral equation scales linearly with the number of dimensions, D.
We adopt the cigenfunctions with the top La cigenvalues in each dimension d as the basis for the projection method, resulting in a set of L = II La basis functions.
The integral equation for each dimension d can rarely be solved analytically, but we can approximate the solution by the Nyström method [35]]: The integral term approximated by a J J-point numerical integration such that Σj=1 k(t, 85 ()((s)))))))))))))))))))))) y = (t))((t), where (Az, (i) is the approximation solution, ( and w - f ldt/ denote the nodes and the weight, respectively, and index dl ís omitted for simplicity of explanation (e.g. (t((), + (tt,, s)); the approximation solution is then obtained by using the eigenvalues et and eigenvectors of the JxJ matrix K, defined by Kjy - k(s,, sy as
X = www,
(i(t) = k(t) ( (t/(er/w),
(21)
where k() A (I:(, ... k(,, 8z)) Under the Nyström method, the computational complexity of finding the eigenfunctions is ((D.3)). We set = 1000 in this paper, which is negligible compared to the cost of the MAP estimation.
4 Related Work
Various algorithms have been proposed to deal with the "doubly- intractable" difficulty of inferencing in Gaussian Cox processes (GCPs). The classical approach is domain discretization [T 11 53, (3X, |3];; it approximates the integral in the likelihood function using a Riemann sum by assuming that the intensity function is constant over each grid. Although domain discretization is applicable to general GCPs regardless of link/kernel function, it suffers from poor scaling in terms of the data size and domain dimension, as well as being sensitive to dis- cretization size. The MCMC approaches [2, 22) provide exact inferencing by exploiting the unique structure of the sigmoidal link function, but they are not applicable to large-scale datasets because demands (N computations. State-of-the-ar algorithms employ the sparse VB approximation with inducing points [E. [i] and achieve computation cost linear to the data size, O(N + L*), where the number of inducing points. However, the VB approaches remain problematic in multi-dimensional domain setting, where the number of inducing points tends to be large. Also, they exploit the unique structures of individual link functions (i.e. quadratic [[I]] and sigmoidal [[I) to obtain feasible algorithms. which limits their application. An approach closely related to ours is the sparse MAP-Laplace method developed by Donner and Opper [[3]]; restricts the link function to being sigmoidal and demands Monte Carlo integration to perform MAP estimation. By focusing on the unique structure of quadratic link function and Matern kernel function, John and Hensman proposed a MCMC- method [24]] that has the same complexity as our model. The merit of [2a] is that it allows for a non non-Gaussian posterior over latent function to improve the approxima- tion quality, while our method is beneficial in that it's deterministic and applicable to various link and kernel functions. Our focus is on a deterministic approach with the most standard Gaussian kernel function, but it is an important next step to make comparative analysis with the scalable MCMC-based method. A summary of the related works about the complexity is provided in the supplementary material ($9).
Gaussian Cox Processes:
Under the quadratic link function, our MAP estimation (12) involves a representer theorem, and the exact MAP estimator with the equivalent kernel [#2] can be derived (see the supplementary material ($6)). The same MAP estimator was derived by Flaxman et al. for the regularized maximum like- lihood problem [[S], while Walder and Bishop derived its predictive covariance through Mercer's
7

---

theorem [95]. Our scheme, with its path integral formulation, is generalization of their works in that our scheme holds (i) for general link functions and (ii) for non-zere mean parameter μ *
(I()) is almost the same as that of quantum mechanics (QM) as posited by Richard Feynman [[Z7], except for the imaginary factor in QM's action potential. In QM, the MAP estimator equation equals the Euler-Lagrange
Path Integral Formulation: The path integral formulation of GCPs
equation, and the Laplace approximation (T3) for functional distribution corresponds to a kind of semi-classical calculation of QM solutions. The connection between GCPs (or generally GP models)
and QM suggests that physics-inspired techniques developed in QM could be applied Bayesian inferencing in machine learning. In statistical physics, for example, Chang et al. have proposed to apply a path integral approach the Bayesian inference problem in membrane biophysics [9]]
Our GP formulation ( agrees with the probability function expressed by the Onsager-Machlup functional (OMF) I&, [66] for linear stochastic differential equation (SDE), which indicates that the inverse kernel function k. (t, t') could be described by using the linear differential operator, L = En (x(ddt)"" of the corresponding SDE. This is consistent with the fact that the kernels in GPs are strongly related to differential operators [[Z]]_ However, our approach differs from OMF in terms of how they represent the MAP estimator: our estimator ((22 is described by means of kernel function k(t, t'), while the OMF approach obtains the estimator by maximizing the posterior described by L. The advantage of the former is, as in the case of kernel trick in the linear regression problem, that it does not need to know the explicit form of L, which would otherwise require an infinite sum of differentials for the closed form of k(t, t). The benefit of the latter is that it is applicable to nonlinear SDEs, which is out of scope for GP.
A possible approach to make our derivation of path integral formulation (8) clearer in measure- theoretic terms is introducing the abstract Wiener measure [2I]] which is closely related to the reproducing kernel Hilbert space [25] and can be used define Gaussian measure in an infinite- dimensional space. Alternatively, the Radon-Nikodym derivative is a promising tool for justifying the path integral expression with the aid of the fertile stream of works by Andrew Stuart and co- authors [[
5 Experiments
A
16
Helt
20
40
30
PFF
VEPP
30
10
15
30
10
|||
10
huy
20
10
30
30
40
30
20
20
10
10
10
20
10
20
10
20
Number of Inducing Points
Number Inducing Points
Number Inducing Points
STVB
True
0
110
20
Figure 1: Results on three types of synthetic data. (A) The predictive performances and the CPU times. Lower values are better. The error bars represent the standard deviations. (B) The estimated intensity functions (L 20). Solid lines and shaded areas represent MAP estimators and 10.15, 0.85] prediction intervals, respectively.
We examined the validity and the potential efficiency of our path integral formulation (PIF) by evaluating it against deterministic and scalable VB-based inference approaches on synthetic and open real-world data. We adopted the state-of-the-art structured variational inference of sig- moidal GCP [3] (STVB), and the variational inference of GCP with quadratic link function [[]]] (VBPP) as the references. For our proposal, we employed GCPs with exponential (PIF.), quadratic (PIF,). and softplus (PIF,) link functions. This time we used a multiplicative Gaussian ker-
nel, ktt,s) = II. (0(22=a))) where the hyper-parameter 0 was optimized by maximizing the
8

---

marginal likelihood (IT)) or the evidence lower bound [3. BI] through grid search. For fair compar- ison, we performed estimations using Adam [27]] with 5,000 iterations across all methods. Each of the CPU times reported is the amount of time required to calculate the MAP, the predictive covari- ance, and the marginal likelihood given hyper-paramest
where computing the eigenfunctions included. A MacBook Pro with 4-core CPU (2.8 GHz) was used. For details, see the supplementary material ($10).
5.1 Synthetic Data
In accordance with [R we created ID data sets generated from three types of intensity func-
tions: (1) = 2 exp( + exp(-[( 25)/10]3) for E [0, 50]; (2 = sin((?) 6 for 1. [0:5]: (t) is piecewise linear function in Figure 1 for E [0, 100), each of which has 20 trial sequences. For each trial, the intensity function was estimated by each of the meth- ods, and the performance was evaluated based on the integrated p-quantile loss [HI defined as IQL, f indicator, the predicted p-quantile of the intensity function, and the true one, respectively. Here, we adopted , (integrated absolute error) and IQL,
X(t)) (plx()> - [I- dt, where I, X(t), and X(t) denote the
Figure MA displays the predictive performances as functionn of the number of inducing points L; it shows that our approaches (PIF) matched the performances of the VB-based methods across the three data sets. In particular, the comparison between PIF, and VBPP is informative because both used the same kernel and link functions but relied on different approximations. VB-based approx- imations are usually better at estimating posterior distributions than Laplace approximations, but Figure IIA shows that our PIF, with Laplace approximation comparable to VBPP under L ž 5 with regard to both IQL,s, 85 and IQL..o- which demonstrates the practical utility for recovering the posterior distribution well point estimation.
Figure MA also plots the CPU times needed for estimation, showing that PIF can be performed several times faster than the scalable VB-based methods. Figure IA-- shows that when L is large, the quadratic link function (PIF.. VBPP) is preferable for the highly-mouulated scenario of (2()). while the sigmoidal and softplus ones (STVB, PIF,) work better on X1 (t) and (t). Link functions should be carefully selected depending on the data, and thus is an advantage of our PIF that supports scalable estimation algorithms for various link functions within the same scheme.
5.2 Real-world Data
Test Log Likelhood (x10
B
Test Log Likelhhood [x10"
PIF.
PIFa
PIFs
13
8TV
42 -6
[+]
VBP
STVB
VBPR
10"
10
16
10
10
10
10
20
10"
15'
Number of Indusing Points
Numbar of Inducing Points
Figure 2: The estimated intensity functions, the predictive performances, and the CPU times on open real-world data. Higher values of log likelihoods and lower values of CPU times are better. (A) 2D neuronal data. (B) 3D spatio-temporal taxi data. The CPU times of STVB and VBPP exceeded 10 hours with L = 153, and the estimations were given up.
real-world data sets. One is a 2D neuronal data, where event locations correspond to the position of a mouse moving in an arena with recorded cell firing [3 BQJ (CC-BY); We adopted the training (Nirain = 583) and the test (Nist 3D taxi dataset in the city of Porto [34]], where the pick-up locations and times are considered as observations (no personally identifiable information is included, CC-BY); We adopted the training (Numin = 1000) and the test (Nies = 3401) data used by [3]. We ran the models on the training data, and evaluated their predictive
We examined the validity of our approach based on
29127) data as randomly assigned by [G]. The other is
performances based on the test log likelihood. Note that to climinate the effect of data assignment

---

bias, the evaluation multiplied estimated intensity by a factor, Nua//Nain We used regular grid over the multi-dimeniional domain T as the inducing points.
Figure Q displays the predictive and time performances as functions of I,, and shows that our PIF- based approaches achieved comparable predictive accuracy while being tens of times faster than the VB-based methods. All the methods typically required a relatively large number of inducing points, L, to achieve good performance, which suggests that our approaches are especially beneficial in the multi-dimensional domain settings because they have lower complexity with L than the VB- based alternatives. Stochastic optimization algorithms could further reduce the complexity, which is demonstrated in the supplementary material ($ 10). Note here that the memory complexity of our proposed method, O(NL++ is the same as that of the VB-based approaches. In the 3D taxi dataset, PIF, output the intensity function with the smallest volatility among the compared methods (see Figure [B), which led to the poor performance. This implies that the softplus link function might tend to underestimate the intensity function modulation given sparsely located data points, which highlights the importance of appropriately selecting the link function. Also, to clarify that our method is scalable to data size, we ran additional experiments a larger (N = 108) taxi dataset, see the supplementary material ($10).
Discussions
We have proposed a novel path integral formulation of Gaussian Cox processes (GCPs), by which we have derived a scalable inference scheme that holds for wide range of GCPs with various types of link functions. Based on synthetic and open real-world data, we con- firmed that our scheme achieves comparable predictive accuracy while being substantially faster than conventional alternatives.
Conclusions:
Future work & limitations: We focused on investigating the practicality of our scheme, and our derivation of the path integral representation is intuitive but not mathematically rigorous. Making our derivation of path integral formulation clearer, especially in measure- measure-theoretic terms, is an important task for future work. Also, we adopted the collocation method to solve the MAP (integral) equation, but exploiting other algorithms for solving integral equations, such as the Nyström method [G)] and the degenerate kernel method [26]), might yield further improvements in computational efficiency.
We did not consider optimizing typer-parameters with gradient descent algorithms in this paper, but VB-based approaches allow for optimization with gradient methods, which are particularly benefi- cial when the kernel function has several kinds of hyper-parameters. It is worth noting that our ap- proach can also employ gradient methods: because the objective function, L(0, (()) log loggppppp (), has an argmin, (0) = ( (0)) (see Eqs. (I)), the hyper-parameter optimization in our ap- proach belongs to the bi-level optimization, and the exact computation of the gradient of L(O, (0)) can be executed [20]] Here we show only the result, d(.(0.ā(0)) aL (V (V3s à Vaf).
where /0, = Er(r(p)))'' When automatic differentiation is employed, the complexity of com- puting the gradient is equal to the sum of the complexity of the marginal likelihood and O(L"), where ( comes from matrix inversion. Here, we note that the gradient of the eigenfunctions and eigenvalues obtained by the Nyström method with respect to the byper-paramed analytically. When gradient descent algorithm is employed for hyper- parem--ereetter optimization, our approach needs to alternate between the optimization in the L-dimensional space (MAP estimation) and that in the (iim)mdniiesio space, while VB-based approaches perform the optimization in higher (1,2 + d0)) (ii((0)-dimena -dimensional space. Although it is out-of-scope of this paper, a comparative analysis between our approach and VB-based alternatives to investigate the practical utility is an important next step in future work.
can be evaluated
Our proposed scheme via the path integral representation can be applied to not only GCPs, but also wide range of Gaussian process models whose likelihood functions include problematic integrals of latent functions. Gaussian process density models (GPDMs) [I]. [2]] for example, come under such models, and our scheme would provide new tangible inference algo- rithm for GPDMs different from MCMC or VB-based algorithms. Although our scheme itself does not contain any ethical problems nor negative societal impacts, applying GCPs for predicting fine spatio-temporal patterns of people's 's behaviors might harm their privacy in some cases, and thus great care should be taken to protect personal information.
Broader impact:
10

---

References
[1] Ryan Prescott Adams, lain Murray, and David J. C. MacKay. The Gaussian process density sampler. In Advances in Neural Information Processing Systems 21, 2008.
[2] Ryan Prescott Adams, lain Murray, and David J.C. MacKay. Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities. In Proceedings of the 26th Annual International Conference on Machine Learning, pages 9-16, 2009.
[3] 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. http://gituub..oo//iii/////////
[4] Kendall E. Atkinson. A survey numerical methods for solving nonlinear integral equations. The Journal of Integral Equations and Applications, 4(1):15 46, 1992.
[5] Sankarshan Basu and Angelos Dassios. A Cox process with log-normal intensity. Insurance: Mathematics and Economics, 31(2):297-302 2002.
[6] Atulum Güneş Baydin, Barak A. Pearlmutter, Alexey Andreyevich Radul, and Jeffrey Mark Siskind. Automatic differentiation in machine learning: survey. Journal of Machine Learning Research, 18:1-43, 2018
[7] Christopher M. Bishop. Pattern Recognition and Machine Learning. Springer, 2006.
[8] Vladimir Igorevich Bogachev. Gaussian Measures, volume 62. American Mathematical Soci- ety, 1998.
[9] Joshua C. Chang, Van M. Savage, and Tom Chou. A path-integral approach to Bayesian inference for inverse problems using the semiclassical approximation. Journal of Statistical Physics, 157(3):582-666 602, 2014.
[10] John P. Cunningham, M. Yu Byron, Krishma V. Shenoy, and Maneesh Sahani. Inferring neural firing rates from spike trains using Gaussian processes. Advances Neural Information Processing Systems 20, 2007.
[11] Masoumeh Dashti and Andrew M. Stuart. The Bayesian approach to inverse problems. In Handbook of Uncertainty Quantification, pages 311-428, Springer, Cham, 2017.
[12] Peter Diggle, Paula Moraga. Barry Rowlingson, and Benjamin M. Taylor. Spatial and spatio- temporal log- aag-Gausiian Cox processes: extending the geostatistical paradigm. Statistical Sci- ence, 284455542- 2013.
[13] Christian Donner and Manfred Opper. processes, Journal of Machine Learning Research, 19:1-34, 2018.
Efficient Bayesian inference of sigmoidal Gaussian Cox
[14] Christian Donner and Manfred Opper. Efficient Bayesian inference for a Gaussian process density model. In Uncertainty in Artificcaa Intelligence, 2018.
[15] Dean G. Duffy. Green's 's Functions with Applications. CRC Press, 2015.
[16] Dimas Abreu Dutra, Bruno Otávio Soares Teixcira, and Luis Antonio Aguirre. Maximum posteriori state path estimation: Discretization limits and their interpretation. Automatica, 50 15):136011368, 2014.
[17] Richard Phillips Feynman. Space-time approach to non-relativistic quantum mechanics. Re- views of Modern Physics, 20:367-387, 1948.
[18] Seth Flaxman, Yee Whye Teh, and Dino Sejdinovic. Poisson intensity estimation with repro- ducing kernels. In Artificial Intelligence and Statistics, pages 279-279. PMLR. 2017.
[19] Israel Gohberg, Seymour Goldberg, and Nahum Krupnik, Traces and determinants of linear operators, volume 116. Birkhäuser, 2012.

---

[20] 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.05444 2016.
[21] Leonard Gross. Abstract Wiener spaces. In Proceedings of the Fifth Berkeley Symposium Mathematical Statistics and Probability, volume 2, pages 31-42, 1965.
[22] 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.
[23] A. Hammerstein. Nichtlineare Integralgleichungen nebst Anwendungen. Acta Mathematica, 54:117-176, 1930.
[24] 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.
[25] G. Kallianpur. Abstract Wiener processes and their reproducing kernel Hilbert spaces. Z. Wahrschhinnicckketstteerr und verw. Gebiete, 17:113-123, 1971.
[26] Hideaki Kaneko and Yuesheng Xu. Degenerate kernel method for Hammerstein equations. Mathematics of Computation, 56(193):141-148, 1991.
[27] Diederik P. Kingma and Jimmy Ba. Adam: method for stochastic optimization. arXiv preprint 1 2014.
[28] Young-Jun Ko and Matthias W. Seeger. Expectation propagation for rectified linear Poisson
regression. In Asian Conference on Machine Learning, pages 253-268. PMLR, 2016.
(29) Sunil Kumar and Ian H. Sloan. A new collocation-type method for Hammerstein integral equations. Mathematics of Computation, 48(118):585-553 1987.
[30] Ye Li, Patrick Brown, Dionne C. Gesink, and Håvard Rue. Log Gaussian Cox processes and spatially aggregated disease incidence data. Statistical Methods in Medical Research, 21(5): 479-507, 2012.
[31] Chris Lloyd, Tom Gunter, Michael Osborne, and Stephen Roberts. Variational inference for Gaussian process modulated Poisson processes. In International Conference Machine Learning, volume 37, pages 1814-1822. PMLR, 2015.
[32] James Mercer. Functions of positive and negative type, and their connection with the theory of integral equations. Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character, 209:415-4466 1909.
[33] Jesper Møller, Anne Randi Syversveen, and Rasmus Plenge Waagepetersen. Log Gaussian Cox processes. Scandinavian Journal of Statistics, 25(3):451-422, 1998.
[34] Luis Moreira-Matias, João Gama, Michel Ferreira, João Mendes-Moreira, and Luis Damas. Predicting taxi-passenger demand using streaming data. IEEE Transactions on Intelligent Transportation Systems, 393-1122,, 2013.
[35] E. J. Nyström. Über die praktische Auflösung von Integralgleichungeen mit Anwendungen auf Randwertaufgaben. Acta Mathematica, 14:185-204, 1930.
[36] Andrei D. Polyanin and Alexander Manzhirov. Handbook Integral Equations. CRC press, 1998.
[37] Carl Edward Rasmussen and Christopher I. Williams. Gaussian Processes for Machine Learning. MIT press, 2006.
[38] 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

---

[39] Francesca Sargolini, Marianne Fyhn, Torkel Hafting, Bruce L. McNaughton, Menno P. Witter, May-Britt Britt Moser, and Edvard Moser. Conjunctive representation of position, direction, and velocity in entorhinal cortex. Science, 312 2006.
[40] Fransesca Sargolini. Grid cell data Sargolini et al 2006. Centre for the Biology of Memory,
2014. https://doi.org/10.. 1188222014.00003
[41] Matthias Seeger, David Salinas, and Valentin Flunkert. Bayesian intermittent demand fore- casting for large inventories. In Advances in Neural Information Processing Systems 29, 2016.
[42] Peter Sollich and Christopher K. Williams. Understanding Gaussian process regression using the equivalent kernel. In International Workshop on Deterministic and Statistical Methods in Machine Learning, pages 211-228. Springer, 2004.
[43] Ming Teng, Farouk Nathoo, and Timothy D. Johnson. Bayesian computation for Sog-Gaussian Cox processes: A comparative analysis of methods. Jourmal of Statistical Computation and Simulation, 87:2227-2233 2252, 2017.
[44] Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics, pages 567-574. PMLR, 2009.
[45] Christian J. Walder and Adrian N. Bishop. Fast Bayesian intensity estimation for the permanen- tal process. In International Conference on Machine Learning, volume 70, pages 3579-3588. PMLR, 2017.
13


Reviews:
[{"review_id": "r5COpvEdiTE", "reviewer": "Reviewer_SdDJ", "summary": "The paper addresses the inference of Gaussian Cox processes (GCP), and provides the path integral formulation to derive a MAP estimate and a Laplace approximation for Bayesian posterior inference. The strongest contribution is, that the formulation holds for general GP prior kernels, and for a variety of link functions. Comparison of the new methods is done on simulated and real datasets.", "questions": "", "limitations": "", "rating": 7, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I found the paper nicely written, and the derivations elegant, and consider them as a good contribution to the active research problem of GCP inference. The formulation allows usage of various link functions and kernels. However, when it comes to the algorithm, I feel that the authors are not discussing all practical challenges.\n\n+ A major point is the optimization of hyperparameters. I know, that is always an ugly task, but it is important for practical applications, and not discussed at all in the manuscript (at least I didn’t find it). I have the feeling, that it is difficult with the proposed task, because you would need to recompute the eigenfunctions of the kernel, every time hyperparameters are updated. That is potentially the reason, why the authors resorted to grid search. However, other approaches allow for optimization with gradient methods (e.g. the VBPP). Also, it is not clear to me, whether the grid search is included in the CPU times, that are reported.\n+ I like the complexity analysis, but I miss the discussion of the complexity of finding the eigenfunctions of the kernel operator. Though, this has to be done only the beginning (as long as kernel is fixed), I think, it should be discussed here. Also, that for finding the basis, D one dimensional integrals have to solved numerically, should be part of the story. In general, I’d appreciate, if the authors could include a bit more details about finding the eigenfunctions in the main text.\n\n## Update\n\nI thank the authors for their response, and I raised my score given the promised changes of the manuscript. I think the final manuscript will be a strong publication.\nI believe the paper is interesting, and I don’t consider any of the points above as a deathblow to the paper, but I would like to see them discussed.\n\nSome minors:\n+ Would be interesting to see comparison also with [24] and [13] in the comparisons, since the first has same complexity as the presented method and the latter is a very related method, as the authors state themselves. If code is available, that is.\n+ I didn’t like that t is used for the space variable and transpose operation. Maybe use \\dagger or \\top for the latter.\n+ L.125: an action integral\n", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "For limitations see above. Societal impact is briefly discussed at the end of the manuscript.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5h"}, "extra_text": {}}, {"review_id": "YJatyGWELMV", "reviewer": "Reviewer_SxYf", "summary": "   The manuscript tackles the problem of scalable and tractable approximate inference for Gaussian Cox models using a path integral approach. Experiments on synthetic 1D data and on 2D neural data illustrate certain properties of the approach.\n\n   Strengths\n   + Gaussian Cox processes are a relevant class of Bayesian models that are computationally intractable and hence a meaningful research objective.\n\n\n   Weaknesses\n   - The paper does not come with code; hence it will be hard to reproduce the results.\n   - The data sets considered in the main paper are rather small compared to [24]. A direct comparison to [24] in the main paper would have been instructive.\n   - The empirical evaluation does not support a solid claim on the runtime and memory footprint of the method.\n   - An empirical comparison to previous methods e.g. [24] is only partly available.", "questions": "", "limitations": "", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "   1) Clarity\n\n   The technical content of the paper on page 3 in particular which is concluded by the strange sentence in lines 117-118 is hard to grasp as the authors chose to use an \"intuitive\" way of deriving the results. A concise summary of the key expressions and the relevant computations e.g. in an algorithm box would have been beneficial.\n\n   2) Originality\n\n   The use of the path integral methodology seems not explored beforehand.\n\n   3) Significance\n\n   The method is specific to the Gaussian Cox process and the benefits w.r.t. to previous methods in terms of theory and empirical performance are not fully established. Hence, the significance of the work is limited.\n\n   4) Reproducibility\n\n   As the code to run the experiments is not included in the submission, it should be rather difficult to reproduce the results.\n\n   5) Empirical analysis\n\n   The manuscript contains experiments on synthetic and real-world data. But the evidence backing up the \"fast\" from the title is limited. It is hard to assess how accurate the approximation is and how it compares to other competing methods e.g. [24]. An experiment backing up the claim in the abstract on the benefits in the \"multidimensional input setting\" where \"the number of inducing inputs tends to be large\" should be added.\n\n   6) Minor Issues and Typos\n\n   line 3: In this paper, we propose\n   line 49: inferencing, lines 80 and 92: of GP prior, line 133: is that it can\n   Figure 2: over10\n   Capitalisation in References: Integralgleichungen, Anwendungen", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"ethical_concerns": "No concerns.", "needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "NKCmXzfe1p0", "reviewer": "Reviewer_jgS8", "summary": "This paper gives a path integral formulation of a Gaussian Cox process (GCP) with an accompanying inference framework for efficiently resolving the integrals. The proposed method is empirically validated on a series of synthetic and real-world data examples where the results seem promising.\n", "questions": "", "limitations": "", "rating": 7, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I found this article to be an enjoyable read. The authors do a very good job of motivating the topic before presenting the proposed model in an intuitive step-by-step framework. There is a lengthy related work discussion that effectively places the proposed model within the GCP literature. Finally, the experimental results are compelling with the proposed path-integral GCP obtaining comparable performance to existing methods, but with substantial computational accelerations. It is of note that some of the concepts in this paper were new to me, however, the accompanying appendix does a good job of providing background information for readers like myself. I have some minor points that I will outline below.\n\nMain comments:\n\n1. Are you able to provide error bars in Figure 2? It would be interesting to see how stable each of the methods is on real-world data.\n2. What is the memory footprint of your method in comparison to STVB and VBPP methods? Given the method's computational runtime performance, this result would be of interest. \n3. Could you clarify what is meant in lines 94-96 of the Appendix? You mention that the kernel's hyperparameters are learned through optimisation of the marginal likelihood, but then write that the same hyperparameter is learned through grid search which confuses me.\n4. I am interested to see how this approach scales to large data sets. I can see that you have considered an example with up to 110705 points in Section 8.3 of the appendix but you only report timings. Do you have any metrics to show the quality of your model here, such as test log-likelihood, as I would be interested to see if the model's capacity to recover the underlying process is maintained with such large datasets?\n\n\nMinor comments:\n\n1. Can you explain the behaviour of PIFq in Figure 2? The model seems unstable \n2. Typo in the caption of Figure 2 - \"over10 hours\".\n3. Given the similarity of all five models' behaviour in Figure 1, I found the plots quite hard to review. My opinion is that these results may be better suited in a table.\n4. I agree with the authors that providing a measure-theoretic formulation of the path integral would be a useful next step for future work in this area.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors do highlight two limitations of this work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "5"}, "extra_text": {}}, {"review_id": "BAcrUjnIVK-", "reviewer": "Reviewer_rjT5", "summary": "This outstanding paper presents an extremely novel approach (for the machine learning community) to the Gaussian process formulation. This formulation is interesting in its own right. But, the paper goes beyond the pure GP interpretation aspect and introduces this formulation for a very specific and precise purpose: to handle the challenging Poisson process likelihood function. This likelihood function is non-i.i.d. as it involves and integral term (equation 2 of the submission). The integral formulation leads to sneaky analytical tricks on the Poisson likelihood, and provides a new way to handle the interesting \"doubly stochastic\" inference problem.\n\nThe authors import some highly technical and novel concepts from physics, to attack the latent GP intensity function of the Poisson process. This leads to a highly efficient and principled algorithm. Moreover, the high degree of novelty is such that the technique may be applied to other problems with latent GP intensity. For example, their method would also allow to handle the likelihood function of survival analysis with latent GP + link function for the hazard function, which also includes an integral term.\n\nThe writing is crisp and compelling. This was a pleasure to behold. The experiments are more than sufficient. It is impressive that the authors made this work in practice, and I hope they see fit to release their code, but if they don't it is okay and the paper more than stands on its own.\n", "questions": "", "limitations": "", "rating": 9, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "My very positive appraisal is in the overview section. Here I give some very minor constructive feedback. \n\nMy main comment regards the discussion of computational complexity. The abstract is slightly misleading, since the stated big O pertains to one iteration, but the number of iterations required clearly depends on N. It would be great to see a table of computational complexities per iteration, for the various alternatives algorithms. This could be broken down by \"predictive mean, predictive variance, marginal likelihood\", and in some cases with an alternative for primal and dual formulations. Also, the O(N^2) on line 256 is a little unclear as the methods under question can be implemented in primal or dual form (scaling in # basis functions of # data), and in the primal form do obtain a similar scaling to the proposed method. \n\nThese are minor constructive criticisms and would be trivial to handle in the camera ready version. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "2"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper presents a path integral formulation for inference in Gaussian Cox process models, which is general for any prior kernel and a variety of link functions. Three out of four reviewers recommend acceptance, with one of the reviewers assessing the paper as a top 15% of accepted papers at NeurIPS. The main criticism of the paper (raised by Reviewer SxYf and shared by other reviewers and myself)  is that the lack of commitment from the authors to make the code publicly available will impede reproducibility. I agree with this point but I understand there are cases where this is not possible. I believe the technical contribution of the paper is strong enough to be beneficial to the NeurIPS community but encourage the authors to improve the reproducibility of their approach. For example, they can include an algorithm (in pseudocode) in the main paper detailing the exact computations required in their method. Besides addressing the reviewers’ concerns, I also encourage the authors to discuss the limitations of their approach a bit further, for example by comparing hyper-parameter estimation in their approach versus other methods such as those based on variational inference.
Decision: Accept (Spotlight)