Fast Bayesian Inference for Gaussian Cox Processes via Path Integral Formulation
Year: 2021

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 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-theart 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.

heading: 1 Introduction
text: Gaussian Cox processes constitute a class of doubly-stochastic point process models in which a flexible prior over a latent intensity function is established by a Gaussian process through a positive link function under which point events are generated. Gaussian Cox processes are the gold standard in analyzing event data in a Bayesian manner, and have a wide ranging list of applications in neuroscience [10], finance [5], and spatio-temporal analysis [30].
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 infinite-dimensional distributions. Many algorithms have been proposed to deal with this difficulty. They include Markov Chain Monte Carlo (MCMC) sampling, domain discretization, and variational Bayesian (VB) approximation. MCMC approaches [2, 22] provide exact inferencing for a specific link function, but they have excessive computation costs. Although domain discretization [10, 43] achieves sub-linear computational 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-art algorithms employ sparse VB approximation with inducing points [3, 31]; 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(NL2+L3) where L and N are the number of inducing points and data points, respectively, making them problematic in the multi-dimensional domain setting, where L tends to be large. Also, VB approaches achieve good approximations by aggressively focusing on the unique structures of individual kernel functions or link functions, and thus are not applicable to a wide range of Gaussian Cox processes.
In this paper, we introduce a conceptually-intuitive 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 (NeurIPS 2021).
formulation predicts that the maximum a 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(NL+ L2) computation scaling, that is, lower complexity than the state-of-the-art 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 3, we provide a practical methodology to perform inferencing. In Section 4, we outline related work. In Section 5, 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 6 states our conclusions.

heading: 2 Path Integral Formulation

heading: 2.1 Gaussian Cox Process
text: We assume that a latent function on a compact space, T ⊂ RD, x(t) : T → R, is generated from a Gaussian Process (GP), denoted by GP(x(t)|µ, k(t, t′)), and observations are generated from a process under latent function x(t), where µ and k are the mean value and the kernel/covariance function for x(t), respectively. Given data set D observed within T , we consider the problem of obtaining the maximum a posteriori (MAP) estimator of x(t) that maximizes the posterior probability
p(x(t)|D) = p(D|x(t)) GP(x(t)|µ, k)∫ Dx(t) p(D|x(t)) GP(x(t)|µ, k) , (1)
where p(D|x(t)) is the likelihood function, and ∫
Dx(t) in the denominator represents the integral over the function or the infinite-dimensional variable x(t).
The GP prior, GP(x(t)|µ, k), can be regarded as a normal distribution for an infinite-dimensional variable (see Section 2.2), and thus if the likelihood function depends only on a set of latent functions on finite N training points, p(D|x(t)) = p(D|{x(tn)}Nn=1), then MAP estimation for an infinitedimensional variable reduces to that for a finite N -dimensional variable by marginalizing out the variables on the points other than {x(tn)}Nn=1; 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 a latent function over T as follows,
log p(D|x(t)) = N∑
n=1
log λ(tn)− ∫ T λ(t)dt, λ(t) = κ ( x(t) ) , (2)
where D = {tn}Nn=1 is the set of point events occurring in the observation region T , κ(x) : R → R+ 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 (1-2) is called the Gaussian Cox process.
Bayesian inference (1) on situation (2) is often referred as “doubly-intractable", because it requires solving the integral of the stochastic function x(t) as well as computing the intractable posterior probability of an infinite-dimensional variable. Below, we address the problem directly by introducing path integral, a mathematical tool familiar in the field of quantum physics, to treat probability distributions of stochastic functions properly and intuitively.

heading: 2.2 Gaussian Process Prior via Path Integral
text: In order to derive an explicit description of GP prior, we first review a well-known formulation of GP [7], and reconstruct its explicit form of probability density distribution by using a path integral representation.
Consider a model defined in terms of a linear combination of M fixed basis functions, {ϕm(t)}Mm=1, such that
x(t) = µ+ M∑ m=1 αmϕm(t) = µ+α ⊤ϕ(t), t ∈ T , (3)
where α ≜ (α1, · · · , αM )⊤ is the M -dimensional weight vector, ϕ(t) ≜ (ϕ1(t), · · · , ϕM (t))⊤, and µ is a mean parameter. If α is generated from an M -dimensional Gaussian with a diagonal covariance matrix of the form p(α) = N (α|0, c−1IM ), (4) the joint probability distribution for any finite (even infinite) set of variables, x ≜ (x(t1), . . . , x(tJ))
⊤ for J ≤ ∞, follows the J-dimensional Gaussian given by p(x)dx = N (x|µ1,K)dx, Kjj′ = k(tj , tj′) ≜ c−1ϕ(tj)⊤ϕ(tj′), (5)
where k(t, t′) is a semi-definite kernel function, c is the precision parameter, 1 ≤ j, j′ ≤ J , and 1 ≜ (1, . . . , 1)⊤. Because (5) holds for arbitrary numbers of dimensions, J , it provides an implicit representation of GP prior.
Next, based on the expression (5), we derive a more explicit representation of GP prior. For simplicity of explanation, we assume that the domain T is one-dimensional, t = t and T = [0, T ], and mean µ is zero. Let K be the integral operator corresponding to k(t, t′), and K∗ be the inverse operator of K,
K(∗)x(t) ≜ ∫ T k(∗)(t, t′)x(t′)dt′, K∗k(t, s) = ∫ T k∗(t, t′)k(t′, s)dt′ = δ(t− s), (6)
where δ(·) represents the Dirac delta function. Note that K∗ and k∗(t, t′) should be described by a differential operator because relation (6) indicates that k(t, s) is the Green function [15] for the operator K∗, but we do not go into detail here. Under the finite scenario (5) with [t0, tJ ] = [0, T ] and J ≫ 1, the second equation in (6) can be approximated by K∗∆K = ∆−1, where ∆jj′ = (tj − tj−1)δjj′ , K∗jj′ = k∗(tj , tj′), and δjj′ is the Kronecker delta. Then we can write down the J-dimensional Gaussian (5) by using K∗ as
p(x)dx = √ |∆K∗∆| (2π)J exp [ −1 2 (∆x)⊤K∗(∆x) ] dx. (7)
By taking the limit of the division number J → ∞ (∆ → 0) for the finite-dimensional Gaussian measure (7), we achieve the following path integral representation of GP prior over the latent function x(t):
GP(x(t)|0, k)Dx ≜ lim J→∞ p(x) J∏ j=1 dx(tj)=
√ 1
|K| exp [ −1 2 ∫∫ T ×T k∗(t, s)x(t)x(s)dtds ] Dx, (8)
where |K| ≜ limJ→∞ |K∆|, Dx ≜ limJ→∞ ∏J
j=1
√ (tj−tj−1)/(2π)dx(tj), and the following
path integral holds, ∫ exp [ −1 2 ∫∫ T ×T k∗(t, s)x(t)x(s)dtds ] Dx = √ |K|. (9) |K| represents the Fredholm determinant or functional determinant [19] of the integral operator K, which is defined by the product of its eigenvalues. See the supplementary material (§1) for the derivation. The expression for GP(x(t)|µ, k) is easily recovered by x(t) → x(t)− µ and t → t in (8).
The path integral representation (8-9) has a clear advantage over the standard one (5) in that the distribution of the latent function x(t) is written in terms of the explicit integral of x(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 (1). To the best of our knowledge, this is the first proposal of a 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 of Feynman’s path integral formulation of the quantum field theory [17], and thus is not mathematically rigorous. In this paper, we confirm that our path integral expression is valid by using experiments rather than taking a more mathematically rigorous approach; we discuss this briefly in Section 4.
2.3 Maximum A Posteriori Estimation
Using the representation (8) and the likelihood function (2), we can rewrite the posterior (1) as the functional,
p(x(t)|D)Dx = 1 p(D)
exp [ −S ( x(t), x(t) )] Dx, (10)
where S ( x(t), x(t) ) is an action integral, defined by
S ( x(t), x(t) ) ≜ ∫ T [ 1 2 (x(t)−µ)x(t)+κ ( x(t) ) − N∑ n=1 log κ ( x(t) ) δ(t−tn) ] dt+ 1 2 log |K|, (11)
and x(t) ≜ ∫ T k
∗(t, t′)(x(t′)−µ)dt′ = K∗(x − µ). (10) indicates that MAP estimation of x(t) is equivalent to the problem of finding the most likely function or path x(t) that minimizes the action integral. Thus we now apply calculus of variations to the action integral, where the functional derivative of S(x ( t), x(t) ) on the MAP estimator x̂ ( t) is equal to zero: δSδx̂(t)δx̂(t)+ δS δx̂(t)δx̂(t) = 0. This leads us to realize the following equation for deriving the MAP estimator x̂(t),
x̂(t) + ∫ T k(t, t′)κ̇ ( x̂(t′) ) dt′ = µ+ N∑ n=1 k(t, tn) κ̇ ( x̂(tn) ) κ ( x̂(tn)
) , t ∈ T , (12) where κ̇(x) ≜ dκ(x)/dx. See the supplementary material (§2) for the detailed derivations of (12).

heading: 2.4 Predictive Covariance
text: One of the advantages of GP models over non-Bayesian approaches is that they can provide predictive distributions. We apply a Laplace approximation to GP models (10), and find the approximate form of the predictive covariance.
We now know the mode of the posterior, x̂(t), and consider a Taylor expansion of functional action potential S ( x(t), x(t) ) centered on the mode such that
S ( x(t), x(t) ) ≃ S ( x̂(t), x̂(t) ) + 1
2 ∫∫ T×T σ∗(t, s)(x(t)−x̂(t))(x(s)−x̂(s))dtds, (13)
where σ∗(t, s) ≜ δ 2S(x,x)
δx(t)δx(s) ∣∣ x=x̂
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 corresponds to the approximation of the posterior process by a GP, and the predictive covariance or the kernel function for the posterior GP, denoted by σ(t, s), can be obtained by the functional inversion of σ∗(t, s) (see Equation (8)), which results in
σ(t, s) = h(t, s)− h(t)⊤(Z +H)−1h(s), (14)
where Znn′ ≜ δnn′ κ 2(x) κ̇2(x)−κ(x)κ̈(x) ∣∣ x=x̂(tn) , κ̈(x) ≜ d 2κ(x) dx2 , Hnn′ ≜ h(tn, tn′), h(t) ≜ (h(t, t1),. . . ,h(t, tN ))⊤, and h(t, s) is defined by the following Fredholm integral equation of the second kind [36],
h(t, s) + ∫ T k(t, t′)κ̈ ( x̂(t′) ) h(t′, s)dt′ = k(t, s). (15)
Note that Znn′ has an infinite value for k(x) = ex, which leads to the relation, σ(t, s) = h(t, s). The full derivation of (14-15) is given in the supplementary material (§3).

heading: 2.5 Marginal Likelihood
text: Let Σ and H be the integral operators for σ(t, s) and h(t, s), respectively: Σ ≜∫ T ·σ(t, t ′)dt′, H ≜ ∫ T ·h(t, t
′)dt′. Under Laplace approximation (13), we can obtain the marginal likelihood or evidence of Gaussian Cox processes, p(D), by performing the following path integral (9),
log p(D) = log ∫ exp [ −S ( x(t), x(t) )] Dx ≃ −S ( x̂(t), x̂(t) ) + 1
2 log |Σ |, (16)
where log |Σ | = log |H| − log |IN +Z−1H|, and
S ( x̂(t), x̂(t) ) = 1
2 log |K|+ N∑ n=1 [ 1 2 (x̂(tn)−µ) κ̇ ( x̂(tn) ) κ ( x̂(tn) ) − log κ(x̂(tn))] +
∫ T κ ( x̂(t) ) dt− 1 2 ∫ T (x̂(t)−µ)κ̇ ( x̂(t) ) dt. (17)
The full derivation is given in the supplementary material (§4).

heading: 3 How to Solve Integral Equations
text: To obtain the MAP estimator x̂(t) and the predictive covariance σ(t, s), we need to solve the corresponding integral equations. We provide efficient algorithms to solve them (also see Suppl. (§8)).

heading: 3.1 MAP Estimator
text: Equation (12) for the MAP estimator is a nonlinear integral equation of the second kind, sometimes called the Hammerstein integral equation [23]. We can obtain an approximation solution of the equation through the most widely used class of methods, projection methods [4], which approximate the solution or its derivation by choosing an approximation from a given finite L-dimensional linear subspace of functions, denoted by Z . We employ the collocation method [29], a variant of the projection method:
κ̇ ( x̂(t) ) ≃ L∑ l=1 β̂lφl(t), {β̂l}Ll=1 = arg min {βl}Ll=1 L∑ l=1 [r(pl)] 2,
r(t) ≜ κ̇ (∑
n
k(t, tn)γ [∑
l βlφl(tn) ] − ∑ l βl ∫ T k(t, t′)φl(t ′)dt′ + µ ) − ∑ l βlφl(t), (18)
where r(t) is the residual derived from (12), {φl}Ll=1 is a basis for Z , {pl∈T }Ll=1 are the collocation points, and γ(κ̇) represents κ̇(x)/κ(x) being re-written as a function of κ̇(x) (see Table 1). As the basis {φl}Ll=1, we adopt the eigenfunctions of the kernel operator K, which reduces the problematic integral term in the residual into a tractable one as ∫ T k(t, t ′)φl(t ′)dt′ = λlφl(t), where {λl}Ll=1 is the set of the eigenvalues of K. The procedure for finding the eigenfunctions is provided at the end of this section. For fair comparison with VB-based approaches of Gaussian Cox processes [3, 31], we solve the minimization problem of the sum of the squared residuals by using a popular gradient descent algorithm, Adam [27]. We can obtain the MAP estimator x̂(t) by substituting (18) into (12): x̂(t) = µ+ ∑ n k(t, tn)γ (∑ l β̂lφl(tn) ) − ∑ l λlβ̂lφl(t).
The sparsely located collocation points, {pl}Ll=1, 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 [44]. In this paper, we adopt a regular grid over T as the collocation points in accordance with the VB-based references [31, 3], but note that other location choices are possible. For simplicity of explanation, we hereafter call {pl}Ll=1 inducing points. It should be noted here that our approach (18) with the collocation method can be applied to cases where γ(κ̇) = κ̇(x)/κ(x) is defined properly as a function of κ̇(x). Table 1 shows that γ(κ̇) can be defined for popular link functions which include quadratic [18, 31], exponential [12, 30, 33], and softplus [28] functions. But our approach is not applicable to link functions whose derivatives, κ̇(x), are not monotonic with respect to x, such as sigmoidal link functions [2, 3, 22].

heading: 3.2 Predictive Covariance
text: Given the MAP estimator x̂(t), the predictive covariance (14) can be obtained by solving the linear integral equation (15) about h(t, s). We apply the Gelerkin method [4], a variant of the projection method, to solve the equation, which results in
h(t, s) ≃ ∑ l ωlφl(t)φl(s), ωl = λl 1 + λlΞl , Ξl = ∫ T κ̈ ( x̂(t) ) φ2l (t)dt. (19)
See the supplementary material (§5) for the detailed derivation. We estimate the integral in {Ξl}Ll=1 via the Monte Carlo method, but the quadratic link function (κ(x) = x2) is an exceptional case in which the integration can be performed analytically, Ξl = 2, due to the orthogonal relation of the eigenfunctions, ∫ T φl(t)φl′(t)dt = δll′ . It should be noted that (19) is the Mercer expansion [32] of function h(t, s), where {ωl}Ll=1 is a set of its eigenvalues.

heading: 3.3 Marginal Likelihood
text: We can evaluate the marginal likelihood by substituting (18) into (16-17), log p(D) = ∑ n log κ ( x̂(tn) ) − ∫ T κ ( x̂(t) ) dt− 1 2 ∑ l λlβ 2 l
+ 1
2 ∑ n
κ̇ ( x̂(tn) ) κ ( x̂(tn) )[∑ l λlβlφl(tn)− (x̂(tn)−µ) ] − 1 2 log |IN+Z−1H| − 1 2 log |K| |H| , (20)
where the second term, ∫ T κ ( x̂(t) ) dt, is estimated by the Monte Carlo method. The last term is associated with the functional determinant (see the supplementary material (§1)), which can be calculated by the product of its eigenvalues as, log |K||H| = ∑L l=1 log λl ωl = ∑L
l=1 log(1+λlΞl), where we approximate the kernel function in terms of the Mercer expansion up to the top L eigenvalues, k(t, s) ≃ ∑L l=1 λlφl(t)φl(s).

heading: 3.4 Computational Complexity
text: The objective function to be minimized in MAP estimation is ∑L
l=1[r(pl)] 2, which needs the com-
putation of O(NL + L2) for each evaluation (see Eq. (18)): the computational complexities of {γn≜γ [∑ l βlφl(tn) ] }Nn=1 and { ∑ n k(pl, tn)γn}Ll=1 are both O(NL), and { ∑ l λlβlφl(pl′)}Ll′=1
or { ∑
l βlφl(pl′)}Ll′=1 needs the computation of O(L2). When a gradient descent algorithm with automatic differentiation [6] is employed, the computational complexity of the MAP estimation per gradient descent step is equal to O(NL + L2), which is one-Lth the complexity of the VB-based approaches [31, 3]. We adopt the algorithm in Section 5.
Furthermore, in view of the shape of the objective function, we can substantially improve the computational efficiency of the MAP estimation by using a mini-batch gradient descent (MGD) algorithm: MGD reduces the computational complexity to O(NL), which is due to {γn}Nn=1. We verify the effectiveness of MGD in the supplementary material (§10). Also, when {pl}l are located regularly over T and k(t, t′) is separable across the dimensions, matrix Φ defined as Φll′ = φl(pl′) has the Kronecker structure (see the the supplementary material (§7)), which can further reduce the complexity to O(NL) under (batch) gradient descent. Note that the cost of the VB-based methods [31, 3] become O(NL2) when the Kronecker structure of the gram matrix is exploited. Once the MAP estimator or the set of coefficient {βl}Ll=1 is obtained, estimating the predictive covariance (14) needs the computation of O((N + L) · min(N2, L2) + NmcL), where Nmc is the number of samples for the Monte Carlo method. The cost O(NmcL) comes from the computation of h(t, s) (19), and the rest is due to the matrix operation. Due to the fact that h(t, s) is a degenerate kernel of rank L, the N ×N matrix H can be decomposed into a product of N ×L matrix R, defined by Rnl = √ ωlφl(tn), and its transpose as H = RR⊤; The matrix operation h(t)⊤(Z +RR⊤)−1h(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 hyper-parameters can be performed based on the marginal likelihood (20), which needs the computation of O((N + L) · min(N2, L2) +Nmc).

heading: 3.5 Eigenfunctions of Kernel Operator
text: The pairs of eigenvalues and eigenfunctions, {λl, φl}∞l=1, for the kernel operator K are obtained by solving a homogeneous Fredholm integral equation, ∫ T k(t, s)φl(s)ds = λlφl(t), where the
set of eigenvalues is at most countable [36]. Here we assume a multiplicative kernel, k(t, s) =∏D d=1 k
(d)(t(d), s(d)), and a hyper-rectangular domain with interval T (d) = [0, T (d)] in each dimension d. The multidimensional integral equation can then be reduced to a set of unidimensional ones, { ∫ T (d) k (d)(t(d), s)φ (d) ld (s)ds = λ (d) ld φ (d) ld
(t(d))}d, via the separation of variables, which results in a multiplicative form of solution, (λl, φl(t)) = (∏ d λ (d) ld , ∏ d φ (d) ld (t(d)) ) , where (λ(d)ld , φ (d) ld
) is the ld-th solution of the integral equation in dimension d, and l = (l1, . . . , lD). Thus the computation complexity of solving the integral equation scales linearly with the number of dimensions, D. We adopt the eigenfunctions with the top Ld eigenvalues in each dimension d as the basis for the projection method, resulting in a set of L = ∏ d Ld 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 is approximated by a J-point numerical integration such that ∑J j=1 k(t, sj)φ̃l(sj)w = λ̃lφ̃l(t), where (λ̃l, φ̃l) is the approximation solution,
{sj}Jj=1 and w ≜ ∫ T 1dt/J denote the nodes and the weight, respectively, and index d is omitted for simplicity of explanation (e.g. k(d)(t(d), s(d)) → k(t, s)); the approximation solution is then obtained by using the eigenvalues el and eigenvectors vl of the J×J matrix K̃, defined by K̃jj′ ≜ k(sj , sj′), as
λ̃l = elw, φ̃l(t) = k̃(t) ⊤vl/(el √ w), (21)
where k̃(t) ≜ (k(t, s1), . . . , k(t, sJ))⊤. Under the Nyström method, the computational complexity of finding the eigenfunctions is O(DJ3). We set J = 1000 in this paper, which is negligible compared to the cost of the MAP estimation.

heading: 4 Related Work
text: Gaussian Cox Processes: Various algorithms have been proposed to deal with the “doublyintractable" difficulty of inferencing in Gaussian Cox processes (GCPs). The classical approach is domain discretization [10, 12, 33, 38, 43]; 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 discretization 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 it demands O(N3) computations. State-of-the-art algorithms employ the sparse VB approximation with inducing points [3, 31], and achieve computation cost linear to the data size, O(NL2 + L3), where L is the number of inducing points. However, the VB approaches remain problematic in a 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 [31] and sigmoidal [3]) 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 [13]; it 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 Matèrn kernel function, John and Hensman proposed a MCMC-based method [24] that has the same complexity as our model. The merit of [24] is that it allows for a non-Gaussian posterior over latent function to improve the approximation 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 a comparative analysis with the scalable MCMC-based method. A summary of the related works about the complexity is provided in the supplementary material (§9).
Under the quadratic link function, our MAP estimation (12) involves a representer theorem, and the exact MAP estimator with the equivalent kernel [42] can be derived (see the supplementary material (§6)). The same MAP estimator was derived by Flaxman et al. for the regularized maximum likelihood problem [18], while Walder and Bishop derived its predictive covariance through Mercer’s
theorem [45]. Our scheme, with its path integral formulation, is a generalization of their works in that our scheme holds (i) for general link functions and (ii) for non-zero mean parameter µ ̸= 0. Path Integral Formulation: The path integral formulation of GCPs (10) is almost the same as that of quantum mechanics (QM) as posited by Richard Feynman [17], except for the imaginary factor in QM’s action potential. In QM, the MAP estimator equation equals the Euler-Lagrange equation, and the Laplace approximation (13) 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 to Bayesian inferencing in machine learning. In statistical physics, for example, Chang et al. have proposed to apply a path integral approach to the Bayesian inference problem in membrane biophysics [9].
Our GP formulation (8) agrees with the probability function expressed by the Onsager-Machlup functional (OMF) [8, 16] 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 =∑
n cn(d/dt) n, of the corresponding SDE. This is consistent with the fact that the kernels in GPs are strongly related to differential operators [37]. However, our approach differs from OMF in terms of how they represent the MAP estimator: our estimator (12) 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 measuretheoretic terms is introducing the abstract Wiener measure [21], which is closely related to the reproducing kernel Hilbert space [25] and can be used to define Gaussian measure in an infinitedimensional 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 coauthors [11].

heading: 5 Experiments
text: 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 sigmoidal GCP [3] (STVB), and the variational inference of GCP with quadratic link function [31] (VBPP) as the references. For our proposal, we employed GCPs with exponential (PIFe), quadratic (PIFq), and softplus (PIFs) link functions. This time we used a multiplicative Gaussian kernel, k(t, s) = ∏ d e −(θ(td−sd))2 , where the hyper-parameter θ was optimized by maximizing the
marginal likelihood (16) or the evidence lower bound [3, 31] through grid search. For fair comparison, 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 covariance, and the marginal likelihood given a hyper-parameter, where computing the eigenfunctions is included. A MacBook Pro with 4-core CPU (2.8 GHz) was used. For details, see the supplementary material (§10).

heading: 5.1 Synthetic Data
text: In accordance with [2], we created 1D data sets generated from three types of intensity functions: λ1(t) = 2 exp(−t/15) + exp(−[(t − 25)/10]2) for t ∈ [0, 50]; λ2(t) = 5 sin(t2) + 6 for t ∈ [0, 5]; λ3(t) is a piecewise linear function in Figure 1 for t ∈ [0, 100], each of which has 20 trial sequences. For each trial, the intensity function was estimated by each of the methods, and the performance was evaluated based on the integrated ρ-quantile loss [41], defined as IQLρ ≜ ∫ T 2 ( λ(t)− λ̂(t) )( ρIλ(t)>λ̂(t) − ( 1−ρ ) Iλ(t)≤λ̂(t) ) dt, where I, λ̂(t), and λ(t) denote the indicator, the predicted ρ-quantile of the intensity function, and the true one, respectively. Here, we adopted IQL.5 (integrated absolute error) and IQL.85.
Figure 1A displays the predictive performances as functions 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 PIFq and VBPP is informative because both used the same kernel and link functions but relied on different approximations. VB-based approximations are usually better at estimating posterior distributions than Laplace approximations, but Figure 1A shows that our PIFq with Laplace approximation was comparable to VBPP under L ≳ 5 with regard to both IQL.85 and IQL.50, which demonstrates the practical utility for recovering the posterior distribution as well as point estimation.
Figure 1A also plots the CPU times needed for estimation, showing that PIF can be performed several times faster than the scalable VB-based methods. Figure 1A-B shows that when L is large, the quadratic link function (PIFq, VBPP) is preferable for the highly-modulated scenario of λ2(t), while the sigmoidal and softplus ones (STVB, PIFs) work better on λ1(t) and λ3(t). Link functions should be carefully selected depending on the data, and thus it is an advantage of our PIF that it supports scalable estimation algorithms for various link functions within the same scheme.
5.2 Real-world Data
We examined the validity of our approach based on multi-dimensional 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 [39, 40] (CC-BY); We adopted the training (Ntrain = 583) and the test (Ntest = 29127) data as randomly assigned by [3]. The other is a spatio-temporal 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 (Ntrain = 1000) and the test (Ntest = 3401) data used by [3]. We ran the models on the training data, and evaluated their predictive performances based on the test log likelihood. Note that to eliminate the effect of data assignment
bias, the evaluation multiplied estimated intensity by a factor, Ntest/Ntrain. We used a regular grid over the multi-dimensional domain T as the inducing points. Figure 2 displays the predictive and time performances as functions of L, and shows that our PIFbased 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 VBbased 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+L2), is the same as that of the VB-based approaches. In the 3D taxi dataset, PIFs output the intensity function with the smallest volatility among the compared methods (see Figure 2B), 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 on a larger (N ≃ 105) taxi dataset, see the supplementary material (§10).

heading: 6 Discussions
text: Conclusions: 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 a wide range of GCPs with various types of link functions. Based on synthetic and open real-world data, we confirmed that our scheme achieves comparable predictive accuracy while being substantially faster than conventional alternatives.
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-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 [35] and the degenerate kernel method [26], might yield further improvements in computational efficiency.
We did not consider optimizing hyper-parameters with gradient descent algorithms in this paper, but VB-based approaches allow for optimization with gradient methods, which are particularly beneficial when the kernel function has several kinds of hyper-parameters. It is worth noting that our approach can also employ gradient methods: because the objective function, L(θ, β̂(θ))≡ log p(D|θ), has an argmin, β̂(θ)≡{β̂l(θ)}Ll=1 (see Eqs. (18-20)), the hyper-parameter optimization in our approach belongs to the bi-level optimization, and the exact computation of the gradient of L(θ, β̂(θ)) can be executed [20]; Here we show only the result, dL(θ,β̂(θ))dθ = ∂L ∂θ − ( ∇βL
)⊤(∇2ββf)−1( ∂∂θ∇βf), where f(θ, β̂)≡ ∑ l[r(pl)]
2. When automatic differentiation is employed, the complexity of computing the gradient is equal to the sum of the complexity of the marginal likelihood and O(L3), where O(L3) 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 hyper-parameter can be evaluated analytically. When a gradient descent algorithm is employed for hyper-parameter optimization, our approach needs to alternate between the optimization in the L-dimensional space (MAP estimation) and that in the dim(θ)-dimensional space, while VB-based approaches perform the optimization in higher (L2 + dim(θ))-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.
Broader impact: Our proposed scheme via the path integral representation can be applied to not only GCPs, but also a wide range of Gaussian process models whose likelihood functions include problematic integrals of latent functions. Gaussian process density models (GPDMs) [1, 14], for example, come under such models, and our scheme would provide a new tangible inference algorithm 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 behaviors might harm their privacy in some cases, and thus great care should be taken to protect personal information.

Review summaries:
   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.     Strengths    + Gaussian Cox processes are a relevant class of Bayesian models that are computationally intractable and hence a meaningful research objective.      Weaknesses    - The paper does not come with code; hence it will be hard to reproduce the results.    - 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.    - The empirical evaluation does not support a solid claim on the runtime and memory footprint of the method.    - An empirical comparison to previous methods e.g. [24] is only partly available.

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.

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.  The 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.  The 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.