Bayesian Optimization through Gaussian Cox Process Models for Spatio-temporal Data

Thank you for your constructive suggestions. Our answers to your questions are provided below.

Reviewer: 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.

Response: Thank you for the question. We would like to clarify that this is indeed our objective. However, expensive-to-evaluate functions relating to discrete point observations over time and space, such as arrivals, events, and occurrences, must be modeled as a Gaussian Cox process, which is a Poisson process modulated by a latent Gaussian process. Standard BO relies on a Gaussian prior and cannot model such expensive-to-evaluate functions with respect to point observation and spatio-temporal data. It requires using doubly stochastic process models, i.e., Gaussian Cox process models considered in this paper.

We clarify that the paper makes two related contributions: (1) a Gaussian Cox process model for expensive-to-evaluate functions relating to discrete point observations over time and space, and (2) a BO framework that builds on such models.

We will also add two illustrative examples to demonstrate the unique advantages of using Gaussian Cox process models for BO. (1) Consider a problem where we need to find the time and location of maximum truck arrivals from video feeds collected by tens of thousands of geo-distributed cameras in a city. This is an expensive-to-evaluate function since it requires performing analytics on a huge volume of video feeds. It is natural to model the truck arrivals (over time and location) as a Gaussian Cox process. We can sample a small set of video clips and use them to estimate the latent intensity, thus supporting efficient BO. (2) Consider the problem of finding the time and location of maximum bird-related incidents in a city. Again, this is an expensive-to-evaluate function, and we have to estimate a latent intensity function from limited samples. We can use the proposed method to model the incidents through a Gaussian Cox process model and then support a BO to solve this problem. For problems requiring estimating a latent intensity function from point observations, we need to use BO with Gaussian Cox process models as proposed in this paper.

Reviewer: 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? What is the theoretical/computational benefit of the paper's approach compared to the cited paper?

Response: Thanks for the question. In the paper, we have compared our results with those reported in [1]. Specifically, this is included in Table 2 with respect to three commonly used synthetic functions and the IQL metric. Regarding the approaches and modeling differences, we would like to clarify that the previous paper [1] uses path integral to obtain an approximation of the solution. In contrast, our approach is to recast the problem into a new RKHS through a change of kernel technique so it becomes easier to handle. It also allows applying the representer theorem to ensure the uniqueness of the estimated mean and covariance. Additionally, our approach enables the use of Nystrom approximation for efficient numerical computation (in Section 3.3), given the change of kernel and the new RKHS. These fundamental differences allow our proposed solution to achieve more favorable results than [1], as demonstrated in Table 2, and make our framework align more with the BO design.

Thank you for your comment. We will certainly consider discussing the approximation related to a few samples and intelligent acquisition strategy with a low number of data points in future work. For now, we have added an experiment starting with a small number of events to observe each step in Fig. 10, Appendix N. The experimental results show that our method can be applied to cases where input events are low without significant estimation deviation from the true intensity.

Thank you for pointing out and correcting the minor mistakes in our manuscript, which will definitely make the paper much more robust. We have updated our manuscript after checking all the remarks you kindly provided. Besides, we provide our responses regarding your concerns as follows.

Reviewer: 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. 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?

Response: Thank you for the question. We will certainly clarify our motivation in the final version. BO based on Gaussian Cox process models is useful when the target function (or dataset, as the reviewer mentioned) is expensive to evaluate.

We will add two illustrative examples where the target function is hard to evaluate (either because it requires performing expensive analytics on a known dataset or due to the fact that obtaining the complete dataset is expensive). (1) Consider a problem where we need to find the time and location of maximum truck arrivals from video feeds collected by tens of thousands of geo-distributed cameras in a city. This is an expensive-to-evaluate function since it requires performing analytics on a huge volume of video feeds. It is natural to model the truck arrivals (over time and location) as a Gaussian Cox process. We can sample a small set of video clips and use them to estimate the latent intensity, thus supporting efficient BO. (2) Consider the problem of finding the time and location of maximum bird-related incidents in a city. Again, this is an expensive-to-evaluate function, and we have to estimate a latent intensity function from limited samples. We can use the proposed method to model the incidents through a Gaussian Cox process model and then support a BO to solve this problem.

We also need to clarify that our proposed BO leverages Gaussian Cox process models to estimate the target function (or dataset as mentioned by the reviewer) from limited samples available. It is a Bayesian method and does not require similar data. The results support the development of acquisition functions that balance exploration (by collecting new samples) and exploitation (to optimize performance objectives) in BO.

Reviewer: 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.

Response: Thank you for your suggestion. Table 4 was included in the appendix as the completed version of Table 2 with standard deviation. Due to the page limit, we shortened the table in the main text to save space. We will also add a description in the corresponding section for clarification and refer readers to the full version of Table 4 in the appendix.

Reviewer: 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. 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?

Response: Sorry about the confusion. In Fig. 2, the black vertical bars at the bottom of each subplot denote the ground-truth events that are not visible to BO. Only the red vertical bars are regions sampled already and thus visible in BO. For instance, at the initial step, our proposed algorithm tries to build a Gaussian Cox process model based on the observed events shown by red vertical bars ($t=25$ and $t=60$).

The green curve depicting the UCB acquisition function anticipates the region of interest for the next step based on our posterior estimates of the current step. In Fig. 2(c), the maximum UCB acquisition value indicates that the algorithm plans to explore around $t=40$ in the upcoming step 15. However, since no events are observed near $t=40$ when we sample, no new observations (i.e., vertical red bars) are added at $t=40$. The Gaussian Cox process model is adjusted based on this new information, and the UCB acquisition function is updated accordingly to guide future samples/steps. We have added the complete consecutive step-wise figure (Fig. 8 in Appendix L, providing a comprehensive demonstration of the BO procedure.

Reviewer: In the experiments, the choice of the hyperparameters $\omega_1, \omega_2, \omega_3$ is not discussed. Can the authors precise their values in each experiment and explain how they were tuned?

Response: $\omega_1, \omega_2, \omega_3$ are weights for mean and covariance terms to balance their contributions in the acquisition function. When these 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 our experiments, we maintain the fixed values for all these parameters, specifically as 0.8.

Furthermore, additional experiments were conducted to assess the impact. In the experiments, we add two extra levels (0.6 and 1.0) for $\omega_3$ to compare with the original $\omega_3=0.8$ in Fig. 6 as an example. The results are shown in Fig. 9 in Appendix M, where if we increase the $\omega_3$ value, the acquisition function will encourage exploring areas with high uncertainty earlier.

Reviewer: The authors approximate the integral $\int_\mathcal{S} \kappa(g(\boldsymbol{t})) d\boldsymbol{t}$ 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?

Response: Since the exact closed-form solutions are not tractable, we discretize the integral term for obtaining the computable closed-form covariance. In the experiment, we set $m$ as the grid size. Therefore, a higher $m$ leads to a finer granularity but higher computation complexity; otherwise a lower $m$ will increase the computational efficiency but with slightly lower quality. To explore the impact of this value, we conduct an additional experiment and show the results in the following table, where we used the synthetic intensity and ran BO for 10 steps. In the table, the runtime will increase with a higher $m$, while the performance does not deviate much.

Regarding the scalability of the approximation, our solution computes the kernel matrix of each dimension separately and then applies the Kronecker product to construct the final kernel matrix as $K = K_1 \otimes \dots \otimes K_d$. For example, given $d=2$ where we have two kernel matrices with size $m_1$ and $m_2$, the time complexity will be $O(m_1^2 m_2^2)$.

Reviewer: The limitations of the proposed framework are not discussed in the paper. Can the authors add a remark on this subject?

Response: Thanks for the question. For future work, we will analyze more complex scenarios with the high-dimensional Gaussian Cox process model, and we will try to add some discussions about this in the final version.

The code is made available via the link.

Thank you for your appreciation of our work. We provide our answers to your questions as follows.

Reviewer: Because many other people have not used Gaussian Cox Process models for BO before, I wonder how much modeling 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.

Response: Thanks for your question. The Gaussian Cox process is the suitable model for addressing problems involving discrete point observations over time and space, such as arrivals, events, and occurrences. This process, characterized by a Poisson process modulated by a latent Gaussian process, proves essential for scenarios where the standard BO relying on a Gaussian prior falls short. The latter cannot effectively model expensive-to-evaluate functions in the context of point observation and spatio-temporal data. Thus, the estimation of latent intensity from discrete point observations becomes necessary, requiring doubly stochastic process models. Gaussian Cox process models are considered as the golden standard for spatio-temporal data [1,2]. This is what motivates our work.

We will provide two illustrative examples demonstrating the unique advantages of using Gaussian Cox process models for BO. (1) Consider a problem where we need to find the time and location of maximum truck arrivals from video feeds collected by tens of thousands of geo-distributed cameras in a city. This is an expensive-to-evaluate function since it requires analyzing a huge volume of video feeds. It is natural to model the truck arrivals (over time and location) as a Gaussian Cox process, as proposed in this paper. We can sample a small set of video clips and use them to estimate the latent intensity, thus supporting efficient BO. (2) Consider the problem of finding the time and location of maximum bird-related incidents in a city. Again, this is an expensive-to-evaluate function, and we need to estimate a latent intensity function from limited samples. We can use the proposed method to model the incidents through a Gaussian Cox process model and then support a BO to solve this problem. For problems requiring estimating a latent intensity function from point observations, we need to use the Gaussian Cox process models proposed in this paper.

Reviewer: 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?

Response: As mentioned in the previous answer, for problems relating to discrete point observations over time and space, such as arrivals, events, and occurrences, they must be modeled as a (doubly stochastic) Gaussian Cox process. Gaussian process models in Standard BO are insufficient. In this paper, we first develop a Gaussian Cox process model to estimate the latent intensity from point observations and then, in Section 3.4, show that various acquisition functions can be developed in our proposed BO framework, supporting standard acquisition functions like UCB as well as more specialized acquisition functions defined on the latent intensity function, such as maximum idle time detection and change point detection, fully exploiting the temporal properties of the underlying data.

Reviewer: Were there any experiments done which could find the posterior mean and covariance in closed form (without Nystrom approximation)? I don't have a good intuition for how much expressivity is lost in doing this approximation.

Response: This is a great question. We have conducted additional experiments to assess the approximation gap. It uses a synthetic set of point events and an RBF kernel, comparing the mean intensity estimations with and without Nystrom approximation. As shown in the following table, the approximation is quite accurate.

We would also like to note that Nystrom approximation can significantly reduce the computational complexity from $O(n^3)$ to $O(n^2m+m^3)$, where $m$ is the number of grid size and smaller than the number of events $n$ [3]. The theoretical analysis of the approximation gap has been discussed in [4]. Additionally, the Nystrom approximation proves beneficial when obtaining a closed-form solution for a specific kernel is not feasible through explicit Mercer’s expansion. We will cite these results and explain them in the final version.