Bayesian Optimization through Gaussian Cox Process Models for Spatio-temporal Data
Authors: Yongsheng Mei
Mahdi Imani
Tian Lan
Venue: ICLR 2024 poster (2024)
Keywords: Bayesian optimization
Gaussian Cox process
Primary area: probabilistic methods (Bayesian methods, variational inference, sampling, UQ, etc.)
TL;DR: Bayesian optimization method on Gaussian Cox process stochastic models for spatial temporal data analysis.

Abstract: Bayesian optimization (BO) has established itself as a leading strategy for efficiently optimizing expensive-to-evaluate functions. Existing BO methods mostly rely on Gaussian process (GP) surrogate models and are not applicable to (doubly-stochastic) Gaussian Cox processes, where the observation process is modulated by a latent intensity function modeled as a GP. In this paper, we propose a novel maximum *a posteriori* inference of Gaussian Cox processes. It leverages the Laplace approximation and change of kernel technique to transform the problem into a new reproducing kernel Hilbert space, where it becomes more tractable computationally. It enables us to obtain both a functional posterior of the latent intensity function and the covariance of the posterior, thus extending existing works that often focus on specific link functions or estimating the posterior mean. Using the result, we propose a BO framework based on the Gaussian Cox process model and further develop a Nyström approximation for efficient computation. Extensive evaluations on various synthetic and real-world datasets demonstrate significant improvement over state-of-the-art inference solutions for Gaussian Cox processes, as well as effective BO with a wide range of acquisition functions designed through the underlying Gaussian Cox process model.

Published as a conference paper at ICLR 2024
BAYESIAN OPTIMIZATION THROUGH GAUSSIAN COX PROCESS MODELS FOR SPATIO TEMTORA DATA
Yongsheng Mei', Mahdi Imani?, Tian Lan'
'The George Washington University
Northeastern Univeristy
{ysmei, clan)egwu... m. m..iaanl@nrrthassteen.ccm
ABSTRACT
Bayesian optimization (BO) has established itself as a leading strategy for effi- ciently optimizing expensive-to-evalle -evallute functions. Existing BO methods mostly rely on Gaussian process (GP) surrogate models and are not applicable to (doubly- stochastic) Gaussian Cox processes, where the observation process is modulated by a latent intensity function modeled as a GP. In this paper, we propose a novel maximum a posteriori inference of Gaussian Cox processes. It leverages the Laplace approximation and change of kernel technique to transform the prob- lem into a new reproducing kernel Hilbert space, where it becomes more tractable computationally. It enables us to obtain both a functional posterior of the latent in- tensity function and the covariance of the posterior, thus extending existing works that often focus specific link functions estimating the posterior mean. Using the result, we propose BO framework based on the Gaussian Cox process model and further develop Nyström approximation for efficient computation. Exten- sive evaluations on various synthetic and real- worldworld datasets demonstrate signifi- cant improvement over state-of-the-a inference solutions for Gaussian Cox pro- cesses, as well as effective BO with a wide range of acquisition functions designed through the underlying Gaussian Cox process model.
INTRODUCTION
Bayesian optimization (BO) has emerged a prevalent sample-efficient scheme for global opti- mization of expensive multimodal functions. It sequentially samples the space by maximizing an acquisition function defined according to past samples and evaluations. The BO has shown suc- cess in various domains, including device tuning (Dalibard et al., 2017). drug design (Griffiths) & Hernández-Lobato 2020), and simulation optimization (Acerbi & Jil 2017). Existing BO ap- proaches are often built on the Gaussian process (GP) regression model account for correlation across the continuous samples/search space and provide the prediction the unknown functions in terms of the mean and covariance. Such models become insufficient when the observation process ee.g., consisting of point events in a Poisson process - is generated from a latent intensity function, which itself established by a GP through a non-negative link function.
To this end, doubly-stochastic point process models such as Gaussian Cox processes are commonly adopted in analyzing spatio-temporal data in Bayesian manner and have been successfully ap- plied to many problems in engineering, neuroscience, and finance (Cunningham et al. 2008 Basu
& Dassios 2002p. The problem is typically formulated as a maximum posteriori (MAP) infer-
ence of the latent intensity function over a compact domain. Its dependence on functional form of the latent intensity function makes an exact solution intractable due to integrals over infinite- dimensional distributions. Existing works often focus on the inference problem for specific link functions, e.g., sigmoidal link function (Adams et al. 2009) Gunter et al. (2014) and quadratic link function (Lloyd et al. 2013) Walder & Bishop (017). or leverage approximation techniques, e.g., variational Bayesian approximation (Aglietti et al. 2019), mean-field approximation (Donner Op- per 2018) and path integral approximations of GP (Kim 2021). However, MAP estimation of the latent intensity function alone is not enough to support BO since optimization strategies in BO must be constructed through acquisition functions that quantify the underlying uncertainty of inference

---

Published as a conference paper at ICLR 2024
through both posterior mean and covariance. This requires building a stochastic model of the latent intensity function to enable BO with such doubly-stochastee models.
In this paper, we consider Gaussian Cox process models with general smooth link functions and formulate a MAP inference of the functional posterior of the latent intensity function and the co- variance of the posterior. We show that by leveraging Laplace approximation and utilizing a change of kernel technique, the problem can be transformed into a new reproducing kernel Hilbert space (RKHS) (Flaxman et al., 2017), where a unique solution exists following the representer theorem and becomes more tractable computationally regarding a new set of eigenfunctions in the RKHS. The proposed approach does not rely on variational Bayesian or path integral approximations in previous work (Aglietti et al. 2019 Kim (021).
We apply the inference model (of both posterior mean and covariance of the latent intensity function) to propose novel BO framework and further develop Nyström approximation when closed-form kernel expansion unavailable. The computation highly efficient since it leverages approxima- tions in previous BO steps to obtain new ones incrementally. To the best of our knowledge, this is the first work on BO using Gaussian Cox process models with kernel transformation in RKHS. enables to consider a wide range of acquisition functions, such as detecting peak intensity, idle time, change point, and cumulative arrivals, through the underlying Gaussian Cox process model. Extensive evaluations conducted using both synthetic functions in the literature (Adams et al., 2009) and real-world spatio-temporal datasets, including DC crime incidents (DC.gov 2022). 2I) neuronal data (Sargolini al.) (006), and taxi data of Porto (O (O'Connell et al., 2015). The proposed method shows significant improvement over existing baselines.
The rest of the paper is organized as follows: Section 2 discusses related work and their differences with our approach. JectionBJdevelog our solution for both posterior mean and covariance, followed by the proposed BO framework with different acquisition function designs using the Gaussian Cox process model. Section presents extensive evaluations on multiple datasets. Finally. Section states our conclusions and underlines potential future work.
RELATED WORK
Gaussian Cox Process Models Doubly stochastic Gaussian Cox processes connecting the Poisson
process with an underlying GP through link function have been the golden standards in analyzing
and modeling spatio-temporal data in many domains (Cunningham et al. 2008 Basu & Dassios
2002]. The MAP intensity estimation problem is much more challenging than estimating a deter-
ministic intensity function (Flaxman et al., 2017). and discretization (Møller et al.
1998 Diggle
al. 2013) is often needed to ensure tractability. Recent works have exploited specific link func-
tion structures and considered Markov chain Monte Carlo methods (Adams et al. 2009: Gunter
cet al. 2014| Nava et al. (022). mean-field approximations (Donner & Opper 2018). Langevin dy- namics (Mutny& & Krause (2022), sparse variational Bayesian algorithms with link functions like sigmoidal functions (Aglictti et al. 2019) and quadratic functions (Lloyd et al., 2015), as well as
transformations into permanental processes (Flaxman et al. 2017 Walder & Bishop 2017). Re- cently, path integral formulation (Kim) 2021) has been proposed as method for effective posterior mean predicative covariance estimations without being constrained by particular link functions. In this work, we proposed a new method for quantifying both posterior mean and covariance of the latent intensity function, which is required to support BO and various acquisition function designs.
Bayesian Optimization through Acquisition Functions Bayesian optimization (BO) can op- timize objective functions that are expensive to evaluate, e.g., neural network hyperparameter tuning (Mei et al.) 2023; Masum et al.] 202T) and drug discovery (Stanton et al. 2022). BO a sequential approach to finding a global optimum of an unknown objective function f()): x= = arg max f (xx, where X' C Rd is compact set. It at each iteration i obtains new sample #x+11 to evaluate f(*++1), update the model of f (often using a posterior GP), and se- lects the next sample using acquisition function
X R from the new model, until reach- ing the optimum x". Various BO strategies have been proposed, including expected improvement
(EI) (Močkus 1975). knowledge gradient (KG) (Frazier et al.] 2200), probability of improvement
(PI) (Kushner 1964,, upper- coper-confidence bounds (UCB) (Lai & Robbins 1985), and entropy search (ES) Heenández-LLobato al. 2014). Despite the success of BO in many practical problems, most
2

---

Published as a conference paper at ICLR 2024
Table I: Common link functions x(x) and their derivatives and inverses.
x(x)
k(x)
K(x))
x (x))
Exponential
exp(x)
exp(x)
exp(x)
log
Quadratic
x²
2r
2
v
Sigmoidal
(1 + exp(-x))
exp(-x)
(L-oxp(a)) exp(a)
log(2 -1
(I+exp(-2)]
(expl()) +
Sofiplus
log(1 exp(x))
(1 exp(-xx)
expi
log(exp(xx) 1)
[1+exp[-2) =])
existing work have not considered BO with Gaussian Cox process models, where point events in Poisson processes are modulated by latent GP intensity functions.
METHODOLOGY
3.1 GAUSSIAN Cox PROCESS MODEL
We consider a Lebesgue measurable compact observation space S C Rd and latent random func- tion g(-) + following Gaussian Process (GP), denoted by (g( μr, Σ) with mean and covariance 4 - (n. Σ). respectively. Observations are generated from a point process modu- lated by the latent random function (t) (through deterministic, nen-negative link function that connects g(t) with latent intensity A(t) = x(g(t))) Given set of n observed point events {t,) in the observation space S, we consider the likelihood function (in terms of GP parameters 8, that is formulated as an expectation over the space of latent random functions, i.e.,
p( p((t:)=-11
p p(tt)) (l)())),))(
g(t)
where p(g(t) p,, is the Gaussian prior. Within equation specific random function g(t) takes the following form:
I the pog-brbbayilityy conditioned
Σ log A(t,) -s X(t)
log p( ((tt))7-199(
where latent intensity function X(t) = x(g(t)) is obtained using the deterministic non-negative link function s<() R + Rt For simplicity of notations, we will use short form g to represent (t)) in the rest of the paper. Table shows some common link functions. Existing works on Gaussian Cox process estimation often focus on specific link functions.
To develop a BO framework using Gaussian Cox process models, we consider the problem of es- timating both the mean and covariance of the Gaussian Process (GP) posterior. Equation ( and equation (2) outline the probabilistic surrogate model aligning with this purpose. In the following, first propose a solution to the posterior mean and covariance estimation problem using Laplace approximation and change of kernel into new RKHS. Following existing work (Møller et al. 1998) Diggle al., 2013). we consider discretization of the observation space to ensure tractability of numerical computation and further propose a Nyström approximation. Then, we will present our BO framework built upon these results, along with new acquisition function designs based on the Gaussian Cox process model.
3.2 ESTIMATING POSTERIOR MEAN AND COVARIANCE
Let (g) A p( ([ g)((g) (,,x) denote the multiplication term inside the integral of equation
We apply the Laplace approximation Illiaa et al.) 2012) estimate the mean and variance of the Gaussian posterior p(g (t, )2)) which requires maximizing log f(g) The result using Laplace approximation is provided in the following lemma (proof in Appendix A)
Lemma (Laplace approximation). For the d d-dimensional multivariate distribution regarding e
Rd given mode such that ğ arg max, log f fgg), the likelihood (7) can be approximated as:
(2TT)
f(g) dg f(g)
|
3

---

Published as a conference paper at ICLR 2024
where A - -
V--g log is the Hessian matrix.
Lemma ! provides the Laplace approximation to the model posterior p(g (t, centered on the
MAP estimate problem:
|A| p(gl(tilia)) J (g) dg (2+)
f(g)
2 2 - - Agg 6) ~ N.gggg,AA (4)
exp
where 9 and A are the mean i and covariance Σ of the GP posterior, respectively.
Before we start solving the posterior mean and covariance, we utilize Lemma to rewrite the Gaus- sian Cox process log-likelihood in equation to obtain:
log p( (t+) 1 ,),x) = -log ( | (g) + log (ggr())), + d log 2x
1 log A
2
(5)
log [t:) |= g) const.,
where approximation (a) follows from the Bayesian information criterion (BIC) (Schwarz 1978). BIC further simplifies the Laplace approximation by assuming that when the event number n is large, the prior p(g px)) is independent of n and the term p({+ (=1)==19) will dominate the rest. allows us treat other independent terms as constant in this inference problem.
We can compute the estimate 9 via a maximization problem max, log ((g) in Lemma 1 which equivalent to minimizing the negative log-likelihood derived in equation (5). Based on the likelihood of the Poisson point process, we have:
min Σ log (gg(((()))++ m(g(t) mlo(t)) 3
(6)
To tackle the optimization problem in equation (6). we utilize the concept of reproducing kernel Hilbert space (RKHS), which a Hilbert space of functions : S + R where point evaluation is continuous linear functional. Given non-empty domain S and symmetric positive definite kernel l S x + R, a unique RKHS 74 can be constructed. If we can formulate a regularized empirical risk minimization (ERM) problem
(/((()/(tt) ml,,, where R(-) denotes the empirical risk of h, 7 is the penalty factor, and (2()) is a non non-decreasing error in the RKHS norm, a unique optimal solution exists given by representer theorem (Schölkopt et al.) 2001) as (+) = Ei= o;k((z,)) :), and the optimization can be cast regarding dual coefficients a E RD.
minkeR. R((h(t)));
Since s(() is non-negative smooth link function, we define l(tt) A (g((t) so that 1:2 k(g(t)) : S + RR. This definition will provide us with access to the property of L2 L2-norm that simplifies the problem-solving in the next. Then, we can formulate minimization problem of the penalized negative log-likelihood as regularized ERM, with a penalty factor y, given by:
Σ log (t,) (²( dt 7l1((t) min k(tt)
/
However, equation 77) does not allow direct application of the representer theorem due to the existence of an extra integral term Ss (t) dt in the optimization. To this end, we show that the term can be merged into the square norm term 7))/((t)]? by a change of kernel technique (resulting in a new RKHS) and using the Mercer's theorem (Rasmussen & Williams) (0066). Specifically, note that the integral term actually defines a L2 -orm of h(t). Taking Mercer's ' expansion of the integral term and the square norm term in Hg. we can add these two expansions together and view the result as Mercer's 's expansion of the RKHS square norm regarding new kernel. It transforms the optimization in equation 7 into a new RKHS (concerning the new kernels), thus enabling the application of the representer theorem. The result is stated in the following lemma (proof in Appendix B).
Lemma (Kernel transformation). The minimization objective J (h) in equation can be written
using a new kernel (t(t, t') as:
J(h) =
Σ log (t,) + |
(8)
4

---

Published as a conference paper at ICLR 2024
where the new kernel function defined by Mercer's theorem is k k= = Ei= m7(n: given that {m. {t) represent the eigenvalues and orthogonal eigenfunctions of the ker- nel function k, respectively.
The new objective is now solvable through the application of the representer theorem. Hence, we can derive the solution h((t) for the original minimization problem in equation 7 and further obtain the g to maximize the likelihood [I. which leads to the following theorem (proof in Appendix C Theorem 1 (Posterior mean). Given observations {ti in S, the posterior mean μ of GP is the solution 9 of the minimization problem 8 laking the form و (h²(()))) where h(.) Ei= ogk(ta,)).
Theorem
I gives the solution posterior mean estimate using the kernel transformation technique in Lemma 2 Similarly, posterior covariance in equation can be obtained in the next theorem (proof Appendix D
Theorem (Posterior covariance). Given observations
[4:]==1 in S, the posterior covariance ma-
trix A-1 is given by
A- E diag {
{1
B(os)alg:) ((() (9 i= IT (9(
1-x2(6 - (O))
where 9s.9. represent g(t.). f(t;), and j and At are from the m-partition Riemann sum of the second term in equation 16 as Is (i(g) dt Ej=1 At.
We note that in Theorem 2 the diagonal values of the covariance matrix exhibit two distinct patterns, depending on whether the observation point &; overlaps with the partition of the Riemann sum 4j. Results in this section establish posterior mean and covariance estimates for the Gaussian Cox process model, which lays the foundation for BO framework.
3.3 NUMERICAL KERNEL APPROXIMATION
The computation of posterior mean and covariance requires solving the new kernel function i in Lemmaβ] which may not yield closed-fome solution for kernels that cannot be expanded explicitly by Mercer's 's theorem. To tackle this, we discretize space S with uniform grid (æ; ) and propose an approximation of the kernel matrix Ktt via a m x m symmetric positive definite Gram matrix K _m + = using Nyström approximation (Baker & Taylor 1979). This method is highly efficient since it leverages approximations in previous BO steps and obtains the next approximation from new samples in incremental fashion. We integrate the grid into kernel matrix
Ket
K
K 1arK Kat.
mat u.uu K = UAUT Σ
(10)
where Ammat, U, are the eigenvalue and eigenvector of K Ken, and A diag(A .A, Am) is the diagonal matrix of eigenvalues. The Nyström method provide us with the estimates of the eigenvalue and eigenfunction of the Mercer's 's expansion, i.e.,
X Amat
de(t)
Kt Uei
(11)
mat
where ktz = (k(t, x.)) This approximation enables efficient computation of posterior mean and covariance using the proposed method while achieving superior estimates than existing baselines as shown later in the evaluation, given by the following lemma (proof in Appendix E Also, since BO requires step-wise intensity estimation, Nyström method allows iterative approximation based on Ktt obtained in previous steps to facilitate the estimation on a larger observation space in the next.
Lemma 3 (Nyström approximation). Based on eigenvalue and eigenfunction in equation 77 the new kernel matrix can be approximated as:
Kut = KiuU ( (3AA+taa)) UT m
(12)
where KT- (kr.)) Kut X. n. matrices.

---

Published as a conference paper at ICLR 2024
3.4 BAYESIAN OPTIMIZATION OVER ESTIMATED INTENSITY
In this section, we introduce our proposed BO framework built upon the Gaussian Cox process model and introduce new acquisition function designs guided by the estimated posterior mean and covariance f that provides a stochastic surrogate model over T C R.. Specifically, in each step i, given the current observed region S; C 'T, our BO obtains a Gaussian Cox process model using posterior mean and covariance estimates and then samples the next observation (e.g., events in a small interval т from unknown parts of 7) according to an acquisition function to expand the observations to S events to E Sy and continues to update the Gaussian Cox process model (as new samples are iteratively obtained) until the BO terminates. The acquisition function is designed capture the desired optimization objective (e.g., finding a peak through expected improvement) relying on the proposed Gaussian Cox process model.
{-} U, This process starts from initial So C T with known
The Gaussian Cox process model introduced in this paper enables a wide range of BO acquisition functions. We note that having only the posterior mean estimate is not enough support BO acqui- sition functions, which require estimating the uncertainty of the underlying surrogate model, c.g., through the notion of probability of improvement (PI) (Kushner) 1964). upper-confidence bounds (UCB) (Lai & Robbins, I985), and expected improvement (El) (Močkus 1975). We demonstrate the design of four acquisition functions enabled by our model in this section, while the BO frame- work can be readily generalized to wide range of BO problems.
Peak intensity prediction: This can be achieved by many forms of acquisition functions (such as UCB, Pl, and El) based on our proposed model. We take UCB as an example, i.e.,
(tcn(t;wi) ((t) + wio(t),
(13)
where denotes the posterior mean, W1 is the scaling coefficient, and a represents the standard deviation obtained from posterior covariance. This formulation contains explicit exploitation and exploration terms, allowing effective identification of the latent peak in Speak-
Maximum idle time: We can calculate the distribution of the number of arrivals of a Gaussian Cox process model:
Proposition 1. Let the point process (t) be modulated by latent intensity function k(g). then the probability of the number of arrivals in the region [a, b] is given by (r,N(a, = n)
k(g) (f" k(g) dt)".
We can plug in the posterior mean and variance derived in Theorems I [and] into the above Propo- sition, i.e., 9 + w2a (t), where We is the scaling coefficient. Therefore, the acquisition function for maximum idle time detection takes the form:
Ogall (t;e) Pr(N(a, s)),
(14)
where is a pre-defined small threshold.
Maximum cumulative arrivals: Finding an interval with maximum cumulative arrivals can be useful for many data analytic tasks. This acquisition function also differs from peak intensity detection since it considers arrivals within an interval. According to Proposition [I] by replacing g + wyyo (t) with the scaling coefficient w3, the cumulative arrival detection acquisition function can be defined as:
Geam (t: = Pr(N b) 2 E).
(15)
where E is a threshold and can be updated based on samples observed as in existing BO.
Change point detection: The occurrence of a sudden change in intensity often indicates latent incident in the real-world scenario. We can integrate the Bayesian online change point detection algorithm into our framework using Gaussian Cox process models. The following lemma extends the results in Adams & MacKay (007).
Lemma 4 (Bayesian online change point probability). Considering the run length Tk with time step k since the last change point given the arrivals so far observed, when Tk reduces to zero, the change point probability is Pr(rk- = 0, tx) = =Erz, Pr(rk-1. -t-1ttt--11111 H (Tk-1),, where is the underlying predictive probability and H(-) is the hazard function.

---

Published as a conference paper at ICLR 2024
The estimated posterior mean and variance can be directly plugged into Lemma 4 to obtain change point acquisition function, i.c.,
GCPD(t;;) = Pr(r 0,
(16)
important to note that the posterior mean and variance are re-estimated each time a new sam- ple is obtained. Thus, the change point acquisition function continuously improved using the knowledge of observed ground-truth arrivals in our BO framework.
The proposed BO framework for intensity estimation with four mentioned acquisition functions implemented on the real-world dataset. The representative results are provided in Section 4.3
EXPERIMENTS
We evaluate our proposed solutions - i.... Gaussian Cox process inference well as BO based on the model - using both synthetic functions in the literature (Adams et al., 2009) and real-world spatio-temporal datasets, including DC crime incidents (DC.gov 2. 2D neuronal data (Sargolini al.) 2006), and taxi trajectory data in city of Porto (O (O'Connell et a 2015)). More evaluation results about acquisition function, additional dataset, and sensitivity study, along with further details about evaluation setups provided in the appendix.
4.1
EVALUATION USING SYNTIETIC DATA
4.1.1
ESTIMATION OF SYNTHETIC INTENSITY FUNCTIONS
10.0
8XbS
KFVS
7.5
STV8
0.5
25
0.0
I11
0.0
||
10
20
40
50
20
40
60
80
100
Input:
Inputtt
Irput:
(a) Au
(b) (t
(c) (t)
Figure Mean estimations comparison on three types of synthetic data.
Following commonly used examples in the literature (Adams et al. (2009): Kim (201)) for Gaus- sian Cox process inference, we consider synthetic intensity functions: X
2exp(-1/15)) + exp(-1(t 25) 2) for E [0, 50). (2 = sin + for E [0. 5]. and Xg that is piecewise linear over (0, 20), (25, 3), (50, 1), (75, 2.5), (100, 3) for E [0, 100]. For three given functions, we generate 44, 27, and 207 synthetic events, respectively. We compare our intensity estimation result with those of existing baselines, including RHKS (Flaxman et al., 2017), MFVB (Donner Opper 2018), STVB (Aglietti et al. 2019) and PIF [K21). The radial basis function (RBF) kernel and quadratic link function are used in this experiment.
The intensity estimation by different methods is visualized in Fig. I Table 2 provides the quan- titative results three metrics, consisting of la-norm to the ground-truth intensity function and integrated p-quantile loss (IQL,,) (Seeger et al. 2016) where we adopt TQT..0 (mean absolute er- ror) and IQL s5 (0. (0.85--uanttle). The smaller numbers of these metrics indicate better performances. Since the PIF code has not made available, we borrow the reported numbers in the paper recorded in IQL only. In the table, our method outperforms the baselines in 7 out of all 9 settings. In partic- ular, when considering function λ2 our result regarding the 12-norm surpasses the STVB by 1.43. For function Ag. the difference between our result and MFVB concerning the IQL .5 is 2.75. This demonstrates the superior performance of the proposed method for Gaussian Cox process inference.
4.1.2 BAYESIAN OPTIMIZATION OVER SYNTHETIC INTENSITY
As the central focus of our work, the proposed BO framework with UCB is performed on the Gaus- sian Cox process model for effective identification of the intensity peaks. The proposed BO frame-
IQL, - S& 2/A(4) - >(W))Kel + (1- (1-p)4xa))< di, where I, A, and Á denote the indicator function, ground-truth intensity, and MAP estimation, respectively.

---

Published as a conference paper at ICLR 2024
Table 2: Average performance on synthetic data regarding // norm and IQL.
(t)
(t)
Aa(t)
Baselines
l2
IQL.50 .50
IQL 85
l2
IQL,5 .5
IQL
l/
IQL 50
IQL,ss
Ours
2.98
11.44
7.38
28.58
12.56
8,59
2.94
30.02
24.32
RKHS
4.37
16.19
12.29
44.78
18.84
11.27
6.63
54.31
53.39
MFVB
3.15
14.36
10.88
32.60
14.41
9.40
3.82
32.77
17.89
STVB
2.96
13.64
10.26
30.01
12.86
8.66
4.19
36.68
21.86
12.50
9.00
13.05
8.65
30.81
20.03
rit Our
Irput:1
Irput:1
(a) Initial step.
(b) Step 8,
(c) Step 14.
(d) Step 20.
Figure 2: Step-wise visualization of BO on synthetic intensity function.
work facilitates effective intensity estimation and computation of the posterior covariance, enabling the step-wise optimization of next samplinglobeervation region through UCB acquisition functions, as shown Fig. 2 The BO process in this test has a budget of 25 steps. As shown in Fig. 2a the algorithm keeps sampling by maximizing UCB acquisition function and then improving the estima- tion based on new samples observed. Three peaks are detected at Step 8 (Fig. 2b Step 14 (Fig. 2c
and Step 20 (Fig. 2db. The completed figure showcasing the consecutive BO procedure is provided
in Appendix L We will present BO with other acquisition functions in Section
4.2 EVALUATION USING REAL-WORLD DATA
4.2.1
INTENSITY ESTIMATION ON 2D SPATIAL DATA
2015
STVB
MPVB
Curs
RKHS
STVB
MFV8
Ours
(a) neuronal data.
(b) Porto data.
Figure 3: Mean estimations comparison on two types of 2D real-world data.
We have compared our method with selected baselines on two 2D real-world spatial datasets and provided the qualitative results in Fig. 3 The first dataset is 2D neuronal data in which event loca- tions relate to the position of mouse moving in 2006). The other dataset consists of the trajectories of taxi travels in 2013 111--21144 in Porto, where we consider their starting locations within the coordinates (41.15, -8.63) and (41.18, -8.60).
with recorded cell firing (Sargolini et al.)
As depicted in Fig. 3. our method successfully recovered the intensity functions. In Fig. 3a our approach exhibits an enhanced ability capture structural patterns reflecting the latent gathering locations of events in comparison to other baselines. Meanwhile, in Fig. 3bb the predicted intensity by our approach is smoother than other baselines, while important structures identified. The performance demonstrated in these 2D comparative experiments underscores the validity of our approach for supporting effective BO using the Gaussian Cox models.
the map are all
4.2.2 BAYESIAN OPTIMIZATION ON TEMPORAE DATA
We utilize public spatio-temporal dataset of crime incidents (i.e., locations and time) in Washing- ton, DC, USA, for 2022. We specifically filtered out 343 crime events involving firearm violence and applied the proposed BO and Gaussian Cox process model those events. Our initial obser- vations consist of the events that took place on May 20th and June 10th. We employ UCB as the acquisition function in this experiment.
8

---

Published as a conference paper at ICLR 2024
1.0
0.58
6.
0.8
0.56
0.7
0.6
034
15
(a) Initial step.
(b) Step 3.
(c) Step 7.
(d) Step 10.
(c) Final step.
Figure 4: Step-wise visualization of BO on 2022 DC crime incidents data.
The results are visualized in Fig. where our proposed method successfully identified the general temporal-spatiee intensity pattern after 10 steps, observing a small fraction of incidents. Since the BO continues sample regions with the highest intensity, it keeps exploring and pinpointing pri- mary latent central locations on the map, which are south-east DC area at step (Fig. (4b), downtown area at step 7 (Fig. 4c), and east DC area at step 10 (Fig. 4d). After Step 10, the resulingg intensity pattern does not undergo major changes, even as more observations were incorporated. ultimately stabilizing by step 25 (Fig. 4e). The entire BO process of 25 steps took only 91.56 seconds. This demonstrates the effectiveness of our BO framework in efficiently identifying spatial patterns and optimizing unknown functions.
4.3 EXPERIMENTS WITH DIFFERENT ACQUISITION FUNCTIONS
(a) BO with O,dle. From left: step 3, 5, and 12.
1.0
65
100
(b) BO with aCPD. From left: step 2, 4, and 12.
Figure 5: BO with Oidle and aCPD coal mining disaster data.
As outlined in Section 3.4 our Gaussian Cox process model and inference solution enable novel acquisition functions beyond the standard ones (c.g., dUCB, dpI, and agI), to address tasks like idle time and change point detection. We evaluate these novel acquisition functions aster data (Jarrett 1979), comprising 190 mining disaster events in the UK with at least 10 casualties recorded between March 15, 1851 and March 22, 1962.
coal mining dis-
Fig. 5a shows the results for Oidee which prioritizes the exploration of regions with low mean and high variance from the current prediction, implying higher probability of fewer arrivals (idle time). When #CPD is adopted, the BO leans toward sampling regions where significant changes in inten- sity are likely to occur based on the posterior mean and variance, as demonstrated Fig. (bb These experiments showcase the versatility of our method in addressing different BO objectives by intro- ducing corresponding acquisition functions.
CONCLUSION
In this paper, we provide novel framework for estimating the posterior mean and covariance the Gaussian Cox process model by employing Laplace approximation and transforming the problem in a new reproducing kernel Hilbert space. The results enable novel BO framework based on Gaus- sian Cox process models, allowing the design of various acquisition functions. Our experimental results on synthetic and real-world datasets demonstrate significant improvement over baselines in estimating apatio-tempore data and supporting BO of unknown functions. The work paves the path to considering more complex scenarios, such as models with the time-variam latent intensity, sparse input spaces, and high dimensions.

---

Published as a conference paper at ICLR 2024
ACKNOWLEDGMENTS
This research based on work supported by the Army Research Office (ARO) under the grant 19F221100299
REFERENCES
Luigi Acerbi and Wei Ji. Practical Bayesian optimization for model fitting with Bayesian adaptive direct search. Advances in Neural Information Processing Systems, 30: 18:1866-184,, 2017.
Ryan Prescott Adams and David JC MacKay, Bayesian online changepoint detection. arXiv preprint arXiv:0710.374 2007.
Ryan Prescott Adams, lain Murray, and David JC MacKay. Tractable nonparametric Bayesian inference in Poisson processes with Gaussian process intensities. In Proceedings of the 26th International Conference on Machine Learning, pp. 9-16. PMLR, 2009.
Virginia Aglietti, Edwin V Bonilla, Theodoros Damoulas, and Sally Cripps. Structured variational inference in continuous Cox process models. Advances in Neural Information Processing Sys- tems, 32, 2019.
Christopher TH Baker and RL Taylor. The numerical treatment of integral equations. Journal of Applied Mechanics, 46(4):969, 1979.
Sankarshan Basu and Angelos Dassios. A Cox process with log-normal intensity. Insurance: Math- ematics and Economics, 31(2):997-302, 2002.
John Cunningham, Krishna V Shenoy, and Maneesh Sahani. Fast Gaussian process methods for point process intensity estimation. In Proceedings of the 25th International Conference on Machine Learning, pp. 192-199, 2008
Valentin Dalibard, Michael Schaarschmidt, and Eiko Yoneki. BOAT: Building auto-tuners with structured Bayesian optimization. In Proceedings of the 26th International Conference World Wide Web, pp. 479-488, 2017.
DC..ov. 2022 Crime incidents in Washington, DC, 2022. URL https//openndat..d...wddaaaa [Accessed: August 02, 2023].
Peter J Diggle, Paula Moraga, Barry Rowlingson, and Benjamin M Taylor. Spatial and spatio- temporal log- Gaussian Cox processes: Extending the geostatistical paradigm. Statistical Science, 28(4):542-563 2013.
Christian Donner and Manfred Opper. Efficient Bayesian inference of sigmoidal Gaussian Cox processes. The Journal of Machine Learning Research, 19(1):2110-7433 2018.
Seth Flaxman, Yee Whye Tch, and Dino Sejdinovic. Poisson intensity estimation with reproducing kernels. In International Conference Artificial Intelligence and Statistics, pp. 270-279. PMLR, 2017.
Peter Frazier, Warren Powell, and Savas Dayanik. The knowledge-gradiig policy for correlated normal beliefs. INFORMS Journal on Computing, 21(():599-611 613, 2009.
Whan-Rhhs Griffiths and José Miguel Hernández- Constrained Bayesian optimization for automatic chemical design using variational autoencoders. Chemical Science, 11(2):577-586, 2020.
Tom Gunter, Chris Lloyd, Michael A Osborne, and Stephen J Roberts. Efficient Bayesian nonpara- metric modelling of structured point processes. In Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, pp. 310-319, 2014.
José Miguel Hernández- Matthew W Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of black-box functions. Advances in Neural Information Processing Systems, 27, 2014.
10

---

Published as a conference paper at ICLR 2024
Janine Illian, Sigrumn H Sørbye, and Håvard Rue. A toolbox for fitting complex spatial point pro- cess models using integrated nested Laplace approximation (INLA). Annals of Applied Statistics, 6(4):1499-1500, 2012.
Richard G Jarrett. note on the intervals between coal-mining disasters. Biometrika, 66(1):191- 193, 1979.
Hideaki Kim. Fast Bayesian inference for Gaussian Cox processes via path integral formulation. Advances in Neural Information Processing Systems, 14:26130-26144 2021.
Harold J. Kushner. A new method of locating the maximum point of arbitrary multipeak curve in the presence of noise. Journal of Basic Engineering, 86:97-106, 1964.
Tze Leung Lai and Herbert E. Robbins. Asymptotically efficient adaptive allocation rules. Advances in Applied Mathematics, 6:4-22, 1985.
PA W Lewis and Gerald Shedler. Simulation of nonhomogeneous Poisson processes by thinning. Naval Research Logistics Quarterly, 26(3):403-413 1979.
Chris Lloyd, Tom Gunter, Michael Osborne, and Stephen Roberts. Variational inference for Gaus- sian process modulated Poisson processes. In Proceedings of the 32nd International Conference on Machine Learning, pp. 1814-1822. PMLR, 2015.
Mohammad Masum, Hossain Shahriar, Hisham Haddad, Md Jobair Hossain Faruk, Maria Valero, Md Abdullah Khan, Mohammad A Rahman, Muhaiminul I Adnan, Alfredo Cuzzocrea, and Fan Wu. Bayesian hyperparameter optimization for deep neural network-based network intrusion detection. 202 IEEE International Conference on Big Data (Big Data), pp. 5413-5419 IEEE, 2021.
Yongsheng Mei, Tian Lan, Mahdi Imani, and Suresh Subramaniam. Bayesian optimization frame- work for finding local optima in expensive multi-modal functions. In ECAI 2023 26th European Conference on Artificial Intelligence, volume 372 of Frontiers in Artificial Intelligence and Ap- plications, pp. 1711-1711. IOS Press, 2023.
Jonas Močkus. On Bayesian methods for seeking the extremum. In Proceedings of the 7th IFIP conference, pp. 400-404440 404. Springer, 1975.
Jesper Møller, Anne Randi Syversveen, and Rasmus Plenge Waagepetersen. Log Gaussian Cox processes. Scandinavian Journal of Statistics, 25(3):451-482, 1998.
Mojmir Mutny and Andreas Krause. Sensing cox processes via posterior sampling and positive bases. In International Conference on Artificial Intelligence and Statistics, pp. 6968-6989. PMLR, 2022.
Elvis Nava, Mojmir Mutny, and Andreas Krause. Diversified sampling for batched bayesian op- timization with determinantal point processes. In International Conference on Artificial Intelli- gence and Statistics, pp. 7031-7054. PMLR, 2022.
Meghan O'Connell, Luis Morcira-Matias, and Wendy Kan. ECML/KKDD 15: Taxi Trajectory Pre- diction, 2015, URL http://www......oommmmom//ooonne hirpeeittt------rrrrtttttttiiiie trajector-- [Accessed: August 10, 2023].
Carl Edward Rasmussen and Christopher KI Williams. Gaussian processes for machine learning. volume 1. Springer, 2006.
Francesca Sargolini, Marianne Fyhn, Torkel Hafting, Bruce I. McNaughton, Menno P Witter, May- Britt Moser, and Edvard Moser. Conjunctive representation of position, direction, and velocity in entorhinal cortex. Science, 312 2006.
Bernhard Schölkopf, Ralf Herbrich, and Alex Smola. A generalized representer theorem. In International Conference Computational Learning Theory, pp. 416-426. Springer, 2001.
Gideon Schwarz. Estimating the dimension of a model. The Annals of Statistics, pp. 461 661-464, 1978.

---

Published as a conference paper at ICLR 2024
Matthias W Seeger, David Salinas, and Valentin Flunkert. Bayesian intermittent demand forecasting for large inventories. Advances in Neural Information Processing Systems, 29, 2016.
Samuel Stanton, Wesley Maddox, Nate Gruver, Phillip Maffettone, Emily Delaney, Peyton Green- side, and Andrew Gordon Wilson. Accelerating Bayesian optimization for biological sequence design with denoising autoencoders. In Proceedings of the 39th International Conference on Machine Learning. pp. 20459-20478. PMLR, 2022.
Christian J Walder and Adrian N Bishop. Fast Bayesian intensity estimation for the permanental process. In Proceedings of the 34th International Conference on Machine Learning, pp. 3579- 3588. PMLR, 2017.
12

---

Published as a conference paper at ICLR 2024
A PROOF OF LEMMA I
Proof. Assuming a mode such that 9 arg maxg log (g), the condition V log f() = 0 must
be satisfied. Therefore, when expanding log fg) using Taylor's 's formula at 9. we have:
9) + 1
log f(g) = log [g) + V log (
5
v) log f()) g)
(17)
و=و log g g). log [g) [g 9) /)"
Let A é - be approximated as:
9=9 log f()). Based on the derived expansion the likelihood in equation
can
1 9
[(o)dg= f(0) / exp 2 9 9) A(g 9) dg f((g)
(2m) 2 |
(18)
This concludes the proof.
B PROOF OF LEMMA 2
Proof.
With the uniform convergence on S S, the Mercer's representation of the kernel function
is:
k(t,t') Σ (t)((( t,t' E S,
(19)
where {m)) and {t) are the corresponding eigenvalues and orthogonal eigenfunctions of kernel function k, respectively.
The second integral term in the original objective 7 is the L2-norm, which is:
[ (t dt ( [|((t)ll
(20)
Since h Hgg we write it as h Ei= 6,0% using Mercer's theorem, leading to | ( EE b3 and ||A||3 Σ 6,n, Here, {6,] signifies the eigenvalues, {w), represents the orthogonal eigenfunctions of h, and {m corresponds to the eigenvalues of the original kernel function k. Therefore, the original objective in equation 7 can be rewritten as:
J(h) Σ
log 122 (4.) + (1(()) + 7 (l1
log h2 (t,) Σ
b3 +7 Σ 83
j=1
(21)
log h2 (t,)) +
Σ
63 + mm,
80
Σ log ((.) +
75 (1g + 7
Considering a eigenvalue of a new kernel ie satisfying its eigenvalue î, ny (ny +7)-1, the second and third terms equation (7) can be merged into a single square norm in H.. In this case, compared to the original expansion in equation (19) the Mercer's expansion of the new RHKS kernel is:
i(t,t') = Σ 7
(t))), t,t' Et,
(22)
ns +
13

---

Published as a conference paper at ICLR 2024
and thus the simplified objective is:
J(h) = Σ log (t.)) |
(23)
This concludes the proof.
C PROOF OF THEOREM 1
if and only if Ej= B3 + +7)n' that is h E Hr is equivalent to h Hi, two RHKS spaces then correspond to exactly the same set of functions. Thus, optimization over H4 equals the optimization over Hg. In this case, we can apply the representer theorem to the simplified objective (8) after kernel transformation into
Proof. Since Σ j=1 b3
Given {t. J any h E Hg can be decomposed into a part that lives in the span of the k and part orthogonal it, i.e.:
h(-) = Σ a,t(tt, v,
(24)
where, for IR" and Hi. we have:
(v, (v,b(t,,)) 0, Vt, E.
(25)
The first term in objective 8 is independent of and only related h(-). Using equation 25 the
reproducing property, application of h(-) to other points tú yields:
h(t.) = (h, k(tt, .)
or a, -) u. k(ts.:)
34
(26)
a t,)) + (v. -),
Σ oak(ti, tz).
As for the second term, since the orthogonal property in equation 25 we get:
[n(t)]
Σ a,t
+
H2
Σ a;i((ti)) :)
+ |1111
(27)
Σ (tt.) .)
Pl.g
where, in the last inequality, the equality case occurs if and only if v = 0. Setting the = 0 thus does have any effect on the first term equation (8). while strictly reducing the second term. Hence, the minimizer must have 0 and the solution takes the form:
h(() = Σ a.k((t..)
(28)
14

---

Published as a conference paper at ICLR 2024
Since define h(t) x² (g(t)), the mean μ of GP should be:
27
g = K-
Σ (rk((i;tt))
(29)
The theorem proven.
D PROOF OF THEOREM
Proof. Since p ~ (H, Σ), we can expand and rewrite log- log-liiel of equation as:
log f(ğ) log p - ( log ((0(0
(30)
log (t. _= 9
2 -
log [2
2
and its second- -odde gradient with respect 9
2 log f(ğ) = vy log
-ー
(31)
31 denote the inhomogeneous Poisson process log-likelihood in equa- tion 2 We approximate its integral term with m-partition Riemann for all dimensions and compute the second-order gradient, resulting:
The first term in equation
V2 logg(((t))
Σ x(g(t())
Σ rlo(t)))
0 log (4.) 10)
log 17-
Ogn 0gi
dym 0gi
(32)
8 log p(t. (=))) 10)
log log (t: Yi 15)
0gi Gg...
8gm Gg...
R(G. R(..)x (m(..) (
N22 (g.)
i=
diag
(#))))
( At
i
Given A = 2 Vz log ((), conclude the proof by providing the variance A in equation (4) with respect to Gaussian mean by applying subtraction:
A- (Σ-1- V3 log p( (4:)
(9)
(33)
where the result of V3 log ((t,)) V (g) is derived in equation 32 To align the dimension between the covariance and gradient results, we discretize each dimension of the covariance matrix via a size m uniform grids to ensure a lossless computation, where we have denoting Kronecker product.
Σ x ... X Eg with X
is worth noting that & and fú may not always be exactly identical in practical implementations. To address this, we can consider them the same if their absolute difference falls within small error range.
PROOF OF LEMMA 3
Proof. We use the ith estimated eigenfunction di(t) and eigenvalue î compute the approximation of the new kernel function in Lemma 2
as:
mat
it t') = Σ
(,(t)),(tt)
A.m mat
im
m
Amat
m ken ue
kir
(34)
mat
+ my (Amatyz
ktt
Σ
(X + Amy)a uis
m
kTT
15

---

Published as a conference paper at ICLR 2024
Since the A in equation 100 is the diagonal matrix of every ith eigenvalue Ammi, by applying the to all input observations, the estimation of the new kernel
derived kernel function estimation matrix Ket becomes:
34)
m
Ket Ktv
m
(A.mary2 + 72. -u;
Kat
1=1
(35)
U
A² + +2a)
U Kzt
m
which concludes the proof.
F PROPOSED ALGORITHMS
We summarize the workflow of the proposed BO framework on Gaussian Cox Process Models in Algorithm [ The numerical approximation approach for posterior mean and variance estimation provided in Algorithm
Algorithm BO on the Gaussian Cox Process
1: initialize Starting observation region Sg, acquisition function a(t), total step T
2: for 0 (T 1) do
3: Estimate posterior mean and covariance of the Gaussian Cox process 04,Σ using Theorem
and2basee on observations in region S;
4: Find observations in next sampling region {r} through the acquisition function a(t: p, regarding current mean and covariance
5: S,+1 +- (7) U.
6: end for
Algorithm Estimations using Nyström Approximation
1: initialize Kernel matrix K, observed events {t:]i-n, discretized grid {x,). random dual coefficient ag. learning rate 8
2: Perform eigendecompossiton on grid K = UAUT E Amat u.u,
3: Derive the eigenvalue and eigenfunction of k(tt,, t2) using equation
4: Update dual coefficient a for new kernel function using gradient descent:
5: Compute the gradient of optimization objective J in equation J
6: for = 0: (n 1)
7: Of++1 o'k 8V.
8: end for
9: Compute new kernel matrix Ke2 regarding & based on Lemma 3
Compute the posterior mean 9 and covariance respectively
G DETAILED EVALUATION SETUPS
We conducted our experiments on the Ubuntu 20.04 system, with Intel(R) Core(TM) i7-6700 4- core CPU (3.4 GHz) and 16.0 GB RAM. The algorithm is implemented in Python 3.8, using main Python libraries NumPy 1.22.3 and Pandas 2.0.3. The code has been made available on GitHub via https://gthub..cm//yse1177/////
To evaluate the BO over synthetic intensity Section 4.1.2 a random ground-truth intensity func- tion is generated (black curve in Fig. (2). where the samples are shown by thinned events (Lewis & Shedler 1979) in the bottom of plots, marked by black vertical bars. In the experiment, we initialize the region centers t (25, 60) in the time domain of [0, 100) with region radius of 2, i.e.. region size of 4. Observations in selected regions are highlighted in red vertical bars. We use the UCB ac- quisition function for identifying intensity peaks. To expedite the BO process, set the acquisition
16

---

Published as a conference paper at ICLR 2024
values of explored regions prevent the algorithm from becoming trapped within the same region.
For different acquisitions experiments in Section 4 utilize coal mining disaster data (Jarrtti 1979), comprising 190 mining disaster events the UK with at least casualties recorded between March 15, 1851 and March 22, 1962. We conducted the experiments with total budget of 12 steps. The initial region centers are (20, 60), and the region radius For visualization purposes, we let the year 1851 be the starting point (i.e., 0) of the time axis, shown in Fig. 5
H ADDITIONAL EXPERIMENTS ON ACQUISITION FUNCTIONS
H.1 INTENSITY PEAK PREDICTION
we provide measures of the l2 distance between our predictions and the ground-truth intensity the steps where intensity maxima are identified. All three intensity peaks in Fig. (from left to right) correspond to the respective maxima in the table. We include information about the step at which each maximum identified, along with current l2 norm, considering different ac- quisition functions for each maximum. We adopt common acquisition functions for deciding the optimum solution, i.c., UCB, PI, and EI. These results illustrate that all three acquisition functions can effectively identify the peaks in the synthetic intensity. Furthermore, we present the l2 norm every 8 steps, demonstrating a decreasing trend in the distances between the ground truth and pre- dicted intensities over time. This experiment underscores the capability of our method efficiently locate peaks in a given point process when equipped with an appropriate acquisition function.
In Table 3
Table 3: Step-wise evaluation to intensity peak regarding different acquisition functions (AFs).
I: maximum Step
2nd maximum
maximum Step
Ig-norm every steps
AFs
lz
Step
f2
12
Step 8
Step 16
Step 24
UCB
10.79
8.91
20
9.16
10.79
10.21
9.27
EI
12
13.68
23
8.07
8
9.92
9.92
10.54
8.51
Pl
9
17.01
12
14.50
25
11.17
17.34
11.79
11.15
H.2 CUMULATIVE ARRIVAL DETECTION
15
1.0
00
100
53
100
00
200
30
100
(a) Initial step.
(b) Step
(c) Step
(d) Step 10.
Figure 6: BO with Geum coal mining disaster data.
Fig. shows the results for aeum, which indicates the region that contains the most cumulative arrivals. This acquisition function differs from peak intensity detection since the selected region by the former might not necessarily have the global maxima but most point events. As showeased in the figure, the designed acquisition function aeum can successfully pinpoint the region of interest. where more arrivals are possibly gathered.
I BAYESIAN OPTIMIZATION ON 2022 USA TORNADO DATA
We utilize a public spatio-temporal dataset of tornado disasters (i.e., locations and time) in the USA in 2022. We specifically filtered out 1126 tornado events with actual loss and applied the proposed BO and Gaussian Cox process model on those events with a budget of 20 steps. We employ the UCB acquisition function in this experiment. The results visualized in Fig. where proposed method successfully identified the general cemporal-spatio intensity pattern after 6 steps. After Step 12, the resulting intensity pattern does not change significantly as more observations fall in
17

---

Published as a conference paper at ICLR 2024
(a) Initial step.
(b) Step 4.
(c) Step 6.
(d) Step 14.
(e) Final step.
Figure 7: Step-wise visualization of BO 2022 USA tornado data.
Table 4: Sensitivity results regarding different link functions in l2- normm IQL, and errors.
(1)
(t)
(a
Baselines
l2
IQL, .50
IQL ss
l
IQL
IQL
l
IQL
IQL .85
2,98
11.44
7.38
28.58
12.56
8.59
2.94
30.02
24.32
Ours (q)
(1.32)
(2.87)
(3.01)
(7.64)
(1.01)
(2.34)
(0.98)
(5.73)
(3.04)
3.11
12.46
8.56
31,22
15.63
12.84
4.00
31.98
26.41
Ours (c)
(1.01)
(3.44)
(3.55)
(8.14)
(1.73)
(2.11)
(1.25)
(6.54)
(6.03)
2.91
11.98
7.21
29.73
15.32
11.37
4.17
33.77
24.15
Ours (8)
(1.27)
(3.21)
(3.18)
(7.87)
(1.69)
(2.85)
(1.88)
(7.91)
(5.34)
4.37
16.19
12.29
44.78
18.84
11.27
6.63
54.31
53.39
RKHS
(1.64)
(3.21)
(3.04)
(11.56)
(2.81)
(2.34)
(2.78)
(8.20)
(7.89)
3.15
14.36
10.88
32.60
14.41
9.40
3.82
32.77
17.89
MFVB
(1.23)
(2.84)
(2.98)
(9.01)
(2.45)
(1.89)
(2.31)
(5.94)
(3.01)
2,96
13.64
10.26
30.01
12.86
8.66
4.19
36.68
21.86
STVB
(1.43)
(2.91)
(2.94)
(8.57)
(2.30)
(2.41)
(2.01)
(6.22)
(3.16)
12.50
9.00
13.05
8.65
30.81
20.03
PIF (q)
(4.40)
(4.29)
(1.88)
(1.70)
(5.89)
(4.12)
9.30
15.68
12.35
32.24
21.05
PIF (e)
(4.33)
(3.74)
(1.67)
(1.99)
(7.14)
(4.79)
11.58
7.62
14.46
10.96
32.59
20.21
(8)
(3.26)
(3.22)
(0.83)
(3.20)
(7.24)
(6.94)
predicted high-intensity region, ultimately stabilizing by step (Fig. 7el. The evolution intensity of this experiment shows that the southern USA will experience more tornado disasters, demonstrating the effectiveness of our BO framework.
J SENSITIVITY EXPERIMENTS REGARDING DIFFERENT LINK FUNCTIONS
We change the link functions used in our framework and show the sensitivity results estimating the synthetic intensity function Table In this table, we use the abbreviations q. and to represent the quadratic, exponential, and softplus link functions, respectively. Detailed results of baseline methods, including statistical errors, are also given in the table. Notably, our results consistently outperform the other baseline methods in 8 out of 9 cases. In particular, for the synthetic function Å1. our approach with the softplus link function yields superior results, while the quadratic link function proves effective in handling other scenarios.
K RUNTIME EXPERIMENTS
We conducted runtime tests of our BO framework on three distinct datasets, each with trial budget of 20 steps. The step-wise and total runtime, along with the error in parenthesis for each trial, are presented in Table Notably, the test involving 1126 events from the USA tornado dataset requires the most time to complete, whereas the tests on DC crime data with 343 events and coal mining disaster data with 190 events need comparatively less time. Additionally, the step-wise runtime exhibits noticeable increasing trend as the number of observed events grows time. This is due to the corresponding increase in estimation complexity as more point events are incorporated into the analysis.
18

---

Published as a conference paper at ICLR 2024
Table BO runtime evaluation (in second) for 20 steps regarding different datasets.
Datasets
Step 5 time
Step 10 time
Step 15 time
Step 20 time
Total time
1.42
2.55
4.22
5.95
68.61
crimes
(0.56)
(0.44)
(0.84)
(1.28)
(1.45)
14.71
58.80
77.89
194.23
820.56
USA tornadoes
(2.56)
(8.73)
(9.75)
(14.67)
(21.40)
0.57
0.96
1,45
2.71
25.96
Coal mining disasters
(0.21)
(0.26)
(0.32)
(0.51)
(0.33)
L VISUALIZATION THE BAYESIAN OPTIMIZATION PROCEDURE
Dur Vantori
PANTUR- 30
29
50
90
362
20
150
100
20
50
382
20
60
100
23
100
50
60
382
20
C0
100
23
100
m0
-0
100
013
80
392
20
100
27
100
Figure 8: Complete BO procedure from step (top left) to 20 (bottom right) on a synthetic intensity.
In this section, we provide complete visualization of Fig. in Section
In the figure, the acquisition function depicted by the green curve anticipates the region of interest for the next step based on our posterior estimates of the current step. For instance, step 1, the acquisition function suggests the region centered at 100. Thus, the algorithm will observe the indicated region at the step though no events exist in that region.
M SENSITIVITY EXPERIMENTS REGARDING WEIGHTS IN ACQUISITION FUNCTIONS
Introduced in Section W1, 42, Wl are weights for mean and covariance terms to balance their contributions the acquisition functions aUCB, aidles and acum- respectively. When these hyper- parameter values are small, BO will place higher emphasis on the mean, favoring exploitation in the process. In contrast, when they are large, BO favors exploration by placing higher emphasis on large variance. In previous experiments, these hyperparameters were fixed at 0.8. In this section, we conducted the sensitivity experiments regarding W3 as the illustrative example. Comparing the results in Fig. and Fig. 6 (where wa 0.8), when reduce the value of w3 0.6, the unexplored
19

---

Published as a conference paper at ICLR 2024
100
-30
wa 0.6. From left: step 1, 5, and
800
1986
(b) wa 1.0. From left: step 1, 2, 5, and 10.
Figure 9: Sensitivity study about tuning w3 for deum-
areas with high mean and low uncertainty will be prioritized, such as region centered at = 30 in step 2, Fig 9a Conversely, increasing w/ to 1.0, shown in Fig. 9b shifts the focus to areas with high uncertainty.
N BAYESIAN OPTIMIZATION ON SPARSE SYNTHETIC DATA
Tyth
Varisnca
LCB
Irpot:1
1111
11
20
Figure 10: Complete BO procedure from step 1 (top left) to 18 (bottom right) sparse data.
In Section 3.2 we utilized an approximation assuming the number of events being large. Therefore, we experiment over a sparse synthetic data with 48 events to demonstrate the capability of the proposed framework. The results of 18 steps are shown in Fig. I0 Initially, the estimation deviates significantly from the ground truth due to the limited observations. However, as the number of observations accumulates during the BO process, the approximation becomes more accurate after each step. the figure, the method demonstrates ability to predict latent intensity peaks with relatively sparse data.
20

---

Published as a conference paper at ICLR 2024
(a) A1-
(b) A2-
(c) Às-
Figure 11: Visualization of tabular results regarding no-m.
O VISUALIZATION OF TABULAR SYNTHETIC INTENSITY RESULTS
In this section, we visualize the tabular results using box plots for (2-norm metric with three synthetic functions Section 4.1.1 In the experiments, we adopt ten replicates for each available baseline. The results are provided in Fig. [] where the red horizontal lines show the medians.
21


Reviews:
[{"review_id": "4qqIMcTbW3", "reviewer": "Reviewer_kMEd", "summary": "In this work, the authors propose a novel method to estimate the posterior mean and covariance of the gaussian cox process model.\n\nThey do this by first approximating the posterior $p(g|{t_i})$ via Laplace approximation, and then using BIC to further simplify the computation. This is in terms of $\\hat g$, which must be solved for by minimizing Eq. 6. To do this, they use RKHS along with a transformation of kernel to make the problem computationally cheap to solve. Once this is done, the posterior mean and covariance can be estimated by $\\hat g$ and the expression in Eq. (9). For kernels that cannot be expanded explicitly, they also discretize and use a Nystrom approximation.\n\nWith a way of estimating posterior mean and covariance, one now is free to choose an acquisition function for the specific problem being solved. The authors discuss various settings in which different acquisition can be applied within this framework.\n\nExperiments are carried out showing both the modelling of the latent intensity, as well as the full framework applied in various spatiotemporal settings.", "questions": "-I'm slightly confused about Section 3.4 in that it seems like one can choose any acquisition function that would solve their problem. What about using Gaussian Cox enables us to do this in contrast to standard BO?\n\n-Were there any experiments done which could find the posterior mean and covariance in closed form (without Nystrom approximation)? I don't have good intuition for how much expressivity is lost in doing this approximation.\n\n-The paper analyzes this method on spatial-temporal data, but couldn't I use this method with any temporal data?", "limitations": "", "rating": 6, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 3, "strengths": "This paper claims to be the first work on BO using Gaussian Cox Process models. I could not disprove this claim through a short search, and if true, I think shows a clear strength in its originality. Every claim seemed technically sound and I could not find any glaring problems, and there were a myriad of experiments demonstrating the method in various synthetic and real world settings. The results present are qualitatively and quantitatively compelling, and the whole paper is relatively clear to understand and well written.", "weaknesses": "Because many other people have not used Gaussian Cox Process models for BO before, I wonder how much modelling the latent intensity actually helps. I did not see any results or discussion on this, but it feels like a useful comparison to make to show that using GCP is actually more performant than standard BO.", "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": "", "extra_scores": {}}, {"review_id": "Z9BhsCMzZD", "reviewer": "Reviewer_bJz1", "summary": "The paper presents a novel framework for conducting BO by leveraging Cox processes. This approach hinges on a Laplace approximation of the likelihood and uses kernel techniques to transform the optimization problem into a RKHS. The framework is empirically evaluated across a range of scenarios, encompassing well-known synthetic functions and real-world databases. The results of numerical experiments indicate that this approach exhibits competitive performance in comparison to other state-of-the-art methods. Unlike the other frameworks, it stands out by enabling BO within the context of Cox process-based models.", "questions": "**Questions**\n- Are the results in Table 2 consistent, i.e. similar results are obtained for a different seed? If no, the authors must consider several random replicates and provide the mean +- std of the results\n- In Figure 2, at the initial step, the UCB acquisition function suggests adding new events at $t > 90$ (since we seek to maximize such criterion) but they are added somewhere else. Similarly, in step 14, the UCB targets the instants around $t = 40$ but events are again added somewhere else. Besides the authors argue that \"the algorithm keeps sampling by maximizing UCB acquisition function and then improving the estimation based on new samples observed\", the plots do not validate their point. Can the authors further explain the results while clarifying my concern? Is it possible to add extra plots at consecutive steps (e.g. steps 1 and 2) for a better understanding of the BO's choice?\n- In the experiments, the choice of the hyperparameters $w_1, w_2, w_3$ is not discussed. Can the authors precise their values in each experiment and explain how they were tuned? \n- The authors approximate the integral $\\int_{\\mathcal{S}} \\kappa(g(t)) dt$ using an $m$-partition Riemann sum to obtain a closed-form of the posterior covariance. Since such approximation depends on $m$, can the authors discuss the quality of the approximation in terms of $m$ and precise how they tune that value in the experimental setup? Can they also discuss the scalability of the approximation when $d$ increases?\n- The limitations of the proposed framework are not discussed in the paper. Can the authors add a remark on this subject?\n\n**Other minor remarks**\n- Page 3, Table 1: the derivatives of the link functions need to be checked. For instance, $\\dot{\\kappa}(x) = 2x$ (quadratic case), $\\dot{\\kappa}(x) = \\frac{e^{-x}}{(1+e^{-x})^2}$ (sigmoidal), $\\ddot{\\kappa}(x) = \\frac{e^{-x}}{(1+e^{-x})^2}$ (softplus), ...\n- Page 3, Section 3.1: $\\Sigma$ is a **CO**variance\n- Page 3, Section 3.1, after Eq. (2): $\\lambda(t) = \\kappa(g(t)) \\to \\lambda(t)$ (it has been already defined before Eq.(1) ) \n- Page 4, after Eq. (7): However, Equation equation (7)\n- Page 4, after Eq. (8): $\\eta_i$ and $\\phi_i(\\cdot)$ need to be defined in the main part of the paper (they were defined in the supplementary material)\n- Page 5, Eq. (10): $\\Lambda = \\operatorname{diag}(\\lambda_1, \\ldots, \\lambda_m)$ needs to be defined\n- Page 7, Section 4: To precise that further details on the \"evaluation setup\" are given in Appendix G\n- Page 7, Section 4.1.1: to indicate the number of events considered in each toy example\n- Page 9, Figure 5: to indicate the iteration step in each panel \n- In Appendix C, Eq. (26): $h(t_j) = \\langle h, \\tilde{k}(t_j, \\cdot) \\rangle_{\\mathcal{H}_{\\tilde{k}}}$ ($j$ rather than $i$)\n- In Appendix C, Eq. (32): the first line must be $\\sum_{i=1}^{n} \\log(\\kappa(g(t_i))) - \\sum_{j=1}^m \\kappa(g(t)) \\Delta t$. Then, the sign of $\\ddot{\\kappa}^2(\\hat{g}_i) \\Delta$ must be inverted.\n- In Appendix C, Eq. (32): given the proposed notation, it is not clear that the dimension of $\\nabla_{\\hat{g}}^{2} \\Psi(\\hat{g})$ matches the dimension of the $d \\times d$ matrix $\\Sigma$. Can the authors clarify this and/or propose a more readable notation?\n- In the References: laplace $\\to$ Laplace (Illian et al., 2012), bayesian $\\to$ Bayesian (Kim, 2021), to add all the authors in (Lai et al., 1985), to complete the reference (Stanton et al., 2022), to be consistent with the names of the journals and conferences and the style of displaying them.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "Theoretical contributions bring together techniques from the machine learning and functional analysis communities. Lemmas and other theoretical developments can be easily verified thanks to the clarity of the discussions. The diversity of examples, which take into account well-studied synthetic functions and real databases, makes it possible to assess the competitiveness of the framework in relation to the literature. The paper is generally well-written and well-organized.", "weaknesses": "Although the strengths lie in both the theoretical and numerical aspects, the motivation for performing BO in point processes lacks practical utility. For example, in the spatio-temporal application describing tornadoes in the USA, I fail to see how new events (tornadoes involving damage) can be sampled sequentially to promote active learning of the intensity function $\\lambda$. I assume that for illustrative purposes, the authors considered adding the \"closest event\" available in the database that matches the BO's suggestion. Is this correct? If this is the case, and if the size of the database allows tractable implementations, we can consider all events for inference of $\\lambda$. If the model cannot handle the whole database, the BO schema is an interesting idea that promotes a threshold between inference quality and the number of observed point events. However, what can be done if no similar events are recorded in the database?  Can the authors give further details on the practical utility of their framework?\n\nThe authors have suggested publishing the Python codes in a Github repository, but there is no evidence of their existence. I suggest sharing an anonymous repository (e.g. via https://anonymous.4open.science/) for further examination.", "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": "", "extra_scores": {}}, {"review_id": "mO6GAzjO9b", "reviewer": "Reviewer_dgpT", "summary": "The paper seeks to build a novel Bayesian Optimization approach that is based on Gaussian Cox Processes for spatio-temporal data. The authors emphasize that Gaussian Cox Processes have never been used within BO settings. The posterior distribution is computed using Laplace approximation and a change of kernel that enables to transform the inference problem into a kernel regression problem. The latter kernel is computed using Nyström approximation. Then the authors present several acquisition functions that are built using posterior mean and posterior variance of the Gaussian Cox Process. Finally the authors illustrate the efficiency of their approach on two types of experiments. First, they show the quality of the mean estimation on synthetic data (4.1.1) and real word data (4.2.1). Second, they show how the BO with UCB function performs for some synthetic dataset (4.1.2) and real world dataset (4.2.2).", "questions": "See questions above.\nIs the code publicly available?\nDepending on author responses, I would change my score.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is well written and the overall approach is scientifically sound and compelling. \nThe motivation is clear, which is to provide a BO framework with Gaussian Cox processes.\nAlthough Laplace approximation is a standard tool, the authors provide an elegant way to derive the maximum of the log likelihood through a trick that transform the problem into a standard kernel regression problem so that they can use the representer theorem. The main results are clearly detailed. The use of Nystrom approximation is standard, but Lemma 3 helps the reader to figure out what has been implemented. \nThe numerical experiments section show that the authors made great effort to compare all the components of their methodology with some of the state-of-the-art approaches. In addition, this has been done for simulated data and real-world data.", "weaknesses": "Although the paper is well written and has many interesting components, there are a couple of points that need to be detailed.\n\nMajor comments:\n\n- My first question is general. In a standard BO setting, we aim at minimizing a function which is costly to evaluate. It seems that this is not the objective of the presented BO problem. It would have been good to make a clear distinction between the two problems.  \n\n - In the literature review, the authors write that \"existing works mostly concentrate on the mean estimation\". I am not an expert in Cox models, but it seems that this has been investigated in [1]. This means that this approach could have also been tested in the numerical experiments with Bayesian Optimization. It also seems that this cited paper uses some tools (Laplace approximation, eigenfunctions decomposition) similar to the ones in the paper. Would it be possible to highlight the main modeling differences with this paper? \nWhat is the theoretical/computational benefit of the paper's approach compared to the cited paper? \n\n-  In Equation (5), the authors claim that the rest of the terms of the likelihood are dominated by the first term when $n$ is large. However, this assumption is not always true in a BO setting (it is not in the numerical experiments). Could the authors comment that point and provide more details?\n\n- The authors claim that they can use the Nystrom approximation to compute the next approximation from new samples in an incremental fashion. Is this step used in the experiment? How does it work in practice? \n\nMinor comments:\n\n- After reading the proof of Lemma 2, there is a point I did not understand. Why does the function $h$ belong to $\\mathcal{H}$? This fact does not look straightforward to me. Perhaps I missed something in the development.\n\n- In the numerical experiments, the authors do not report the results of PIF in Table 2. Is there a reason for that? I've not found it in the paper.\n\n- In figure 2c, it looks like the maximum of the acquisition function is around t=40. This means that we have to sample around $t=40$. However I don't see any sample around $t=40$ in figure 2d. \n\n\n[1] Hideaki Kim. Fast Bayesian Inference for Gaussian Cox Processes via Path Integral Formulation. Advances in Neural Information Processing Systems  2021.", "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": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)