PaperlayerView exact revision on arXiv

Source-linked AI summary

ISBO: Scalable Spatio-Temporal Bayesian Optimization with Log Gaussian Cox Process Models via the INLA-SPDE Approach

Kaichuang Yang, Håvard Rue, Jakob Zeitler

arXiv:2610.12213v1stat.MLcs.LG

TL;DR

Standard Bayesian optimization does not fit spatio-temporal point events generated by a latent intensity, and repeated Cox-process inference can be costly. ISBO combines LGCP models with INLA-SPDE inference and region-aware acquisition; experiments report improved optimization performance and computational efficiency over existing baselines.

  • Problem

    Standard Gaussian-process BO is ill-suited to spatio-temporal point events from a random latent intensity, and repeated inference for sequential optimization can be computationally expensive.

    Source: §1 · p. 1

  • Method

    ISBO combines an LGCP surrogate with INLA-SPDE inference and sequential region-aware optimization using time-varying UCB and masking.

    Source: §1 · p. 2, Our contributions · p. 2

  • Results

    Across synthetic and real-world experiments, ISBO consistently improves over existing baselines in optimization performance and computational efficiency.

    Source: §5 · p. 8

  • Takeaways & Limitations

    ISBO supports efficient optimization of latent intensity across temporal, spatial, and spatio-temporal settings while maintaining uncertainty quantification and computational tractability as events are revealed.

    Source: §5 · p. 8

Abstract

from arXiv · show

Bayesian Optimization (BO) is a popular method for efficiently optimizing expensive black-box objectives. However, BO utilizing standard Gaussian Processes is ill-suited for doubly stochastic Cox Processes that are often used in spatio-temporal problem spaces. We introduce INLA-SPDE Spatio-Temporal Bayesian Optimization (ISBO): the first scalable BO framework for spatio-temporal data, that models the log-intensity with a Log-Gaussian Cox Process(LGCP) and performs inference via Integrated Nested Laplace Approximation and Stochastic Partial Differential Equations (INLA-SPDE) approach. Using a Matern field on meshes yields a sparse Gaussian Markov Random Field, where INLA provides fast and accurate posterior inference throughout sequential optimization. ISBO stably locates high-intensity regions and the peak of the latent intensity with minimal evaluations. A time-varying Upper Confidence Bound acquisition with masking avoids revisits, while penalized-complexity priors regularize early rounds. Experiments on synthetic and real-world spatio-temporal datasets show accurate peak discovery, intensity recovery, and substantial speedups over an RKHS-based baseline, positioning ISBO as a practical choice for BO with point-process data.

1 INTRODUCTION

Bayesian optimization is designed for costly black-box evaluations, but its standard Gaussian-process setup does not fit spatio-temporal point events generated by a random latent intensity. ISBO combines LGCP inference with Bayesian optimization for these data and reports improvements over existing baselines across synthetic and real-world experiments.

1 INTRODUCTION

  • Standard BO models real-valued, possibly noisy query outcomes, unlike spatio-temporal point events generated by a latent intensity.The intensity is random, and events are conditionally generated by an inhomogeneous Poisson process.Source: §1 · p. 1
  • Repeatedly refitting a Cox-process model makes sequential BO computationally demanding, while acquisition depends on reliable posterior uncertainty in the latent field.Bayesian inference methods such as MCMC and variational inference can provide uncertainty but may be expensive.Source: §1 · p. 1
  • ISBO combines a log-Gaussian Cox Process with INLA and SPDE for scalable BO over temporal, spatial, and spatio-temporal point-process data.The paper adopts Mei et al. (2024a) as a baseline and reports improving over it.Source: §1 · p. 2

Our contributions

  • The framework optimizes latent intensity from sequentially revealed event information using time-varying UCB, masking, and region-based regret measures.The acquisition strategy is tailored to point-process observations.Source: Our contributions · p. 2
  • Across synthetic and real-world datasets, ISBO improves optimization performance and computational scalability over existing baselines.Ablation studies assess the contributions of its key model components.Source: Our contributions · p. 2

2 RELATED WORK

Related work distinguishes standard BO from optimization over point-process data and reviews LGCP inference approaches suited to latent Gaussian models. It positions ISBO against an RKHS-based BO baseline and a fast approximate LGCP inference method.

Log-Gaussian Cox Process Models

  • LGCPs model the logarithm of intensity as a Gaussian process and provide a flexible class for spatial prediction, though intensity estimation is harder because Cox processes are doubly stochastic.The passage describes LGCPs as the point-process analogue of linear Gaussian models in geostatistics.Source: §2 · p. 2

Integrated Nested Laplace Approximation

  • INLA provides deterministic posterior approximations for LGCP inference and can save computation relative to MCMC while retaining accurate uncertainty estimates.This speed–accuracy trade-off is relevant when the surrogate must be repeatedly refitted after queries.Source: §2 · p. 2

Baseline 1: Laplace Approximation with Kernel Transformation in Reproducing Kernel Hilbert Space (RKHS)

  • Mei et al. (2024a) reformulate a Laplace-approximated MAP objective in an RKHS and use Nyström approximation for posterior mean and covariance estimation.The authors report variability, initialization and hyperparameter sensitivity, limited scalability, and slower computation in their reproductions; a large-n assumption may not hold early in BO.Source: §2 · p. 2

Baseline 2: Fast Approximate LGCP Inference

  • Dovers et al. (2023) use reduced-rank interpolation, likelihood approximation, and automatic differentiation for fast approximate LGCP inference.Their reported speed-ups over the INLA implementation available at the time came with a small loss in model fit; the paper says this method is slower than modern INLA and its runtime increases with event count.Source: §2 · pp. 2–3

Bayesian Optimization through Acquisition Functions

  • Most BO work uses GP surrogates and continuous domains; Mei et al. (2024a) is identified as a baseline, while the paper says existing BO research otherwise does not address GCP models.The focus here is discrete temporal or spatial point events.Source: §2 · p. 3

3 METHODOLOGY

ISBO couples an off-grid LGCP surrogate inferred with modern INLA–SPDE to sequential Bayesian optimization over latent intensity. Its acquisition varies exploration over time and masks already revealed regions.

3.1 LGCP Surrogate Model with Modern INLA–SPDE

  • Off-grid SPDE represents the field with piecewise-linear mesh basis functions, retaining exact event locations and yielding a sparse Gaussian Markov random field.Integration points approximate the likelihood integral, while observed events remain at their original locations.Source: §3.1 · pp. 3–4
  • Modern INLA preserves sparse conditional posterior precision; PC priors regularize the latent field with few revealed events, and a low-rank variational Bayes correction improves its posterior mean.The PC priors are particularly relevant during early BO iterations.Source: §3.1 · p. 4

3.2 Bayesian Optimization over Estimated Intensity

  • ISBO sequentially fits an LGCP surrogate as queries reveal point-process events within local regions, then uses the posterior to select subsequent queries.The objective is to identify a maximizer of latent intensity or a high-intensity region with as few revealed regions as possible.Source: §3.2 · p. 5
  • Because the exponential link preserves maximizers, ISBO acquires on the latent scale with a decreasing-exploration UCB and masks regions whose events were already revealed.Under a fixed SPDE mesh and model specification, the paper states that refitting time remains constant as revealed events increase.Source: §3.2 · p. 5

4 EXPERIMENTS

Experiments assess intensity estimation and Bayesian optimization on synthetic and real-world temporal, spatial, and spatio-temporal point-process data. They report lower regret for ISBO than its baseline and near-constant INLA–SPDE fitting time as event counts grow.

4.1 Evaluation on One Dimensional Data

  • ISBO’s end-to-end temporal fit and projection took 2.27s, compared with 35.74s for the baseline.Both methods showed peaks around days 30 and 50; mesh padding and PC priors avoided boundary collapse in ISBO’s estimate.Source: Figure 1 · p. 6, §4.1.1 · p. 6

4.1.2 Bayesian Optimization over 1D Synthetic Intensity

  • By iterations 14–15, ISBO mean simple regret reached zero; its PC-prior and time-varying-UCB ablations were less robust, and all ISBO variants beat Mei et al. (2024a).The fixed-UCB variant reduced regret more slowly with wider intervals, while the Gaussian-prior variant retained positive regret.Source: §4.1.2 · pp. 6–7, Figure 2 · p. 7

4.2 Evaluation on Two Dimensional Data

  • INLA–SPDE yielded more detailed Porto intensity surfaces and near-constant fitting time through 55,000 events, while Nyström and SCAMPR runtimes grew.The constant-runtime pattern is reported under a fixed SPDE mesh and model specification; INLA–SPDE was faster than both SCAMPR approximations in the tested large-scale regime.Source: §4.2.1 · p. 7, Figure 3 · p. 8

4.2.2 Bayesian Optimization on Spatio-Temporal Data (Washington, D.C. Crime)

  • On Washington, D.C. crime data, ISBO reduces simple regret faster and approaches zero, with lower cumulative regret and narrower later intervals, while enabling direct two-dimensional spatial BO.The comparison with Mei et al. (2024a) uses the same temporal BO setting.Source: §4.2.2 · pp. 7–8, Figure 4 · p. 8

5 DISCUSSION

ISBO is presented as a scalable framework for point-process optimization, with experiments showing improved optimization performance and computational efficiency. The discussion highlights its practical and methodological implications, while noting possible extensions beyond the experiments.

5 DISCUSSION

  • Across synthetic and real-world experiments, ISBO improves over existing baselines in optimization performance and computational efficiency.The ablations quantify contributions from its main design components.Source: §5 · p. 8
  • ISBO supports efficient optimization of latent intensity across temporal, spatial, and spatio-temporal settings while maintaining uncertainty quantification as events are revealed.The discussion identifies computational tractability as part of this practical scope.Source: §5 · p. 8
  • Coupling posterior intensity estimation with region-aware acquisition and masking provides a framework for adaptive acquisition under structured event observations.The discussion presents this integration as a methodological contribution and suggests further work on richer acquisition strategies.Source: §5 · p. 8
  • The authors suggest that scalable INLA inference may extend to machine-learning problems involving uncertainty quantification and sequential decision making, including multi-agent reinforcement learning.The proposed extensions include hyperparameter tuning and adaptive exploration.Source: §5 · p. 8

AI use statement

  • Generative AI was used for limited writing support and routine coding assistance, not for the paper’s core research ideas or scientific conclusions.The authors report reviewing and revising AI-assisted text and code where necessary.Source: AI use statement · p. 9

CHECKLIST

This section documents the INLA–SPDE and latent Gaussian machinery underlying ISBO, including low-rank corrections to posterior inference and the role of PC priors in early optimization. It also provides reproducibility materials and checklist items.

CHECKLIST

  • The checklist marks model descriptions, theoretical assumptions and proofs, and empirical reporting as complete, while human-subject requirements are not applicable; the anonymous repository contains code, benchmarks, data, maps, and result files.The checklist also reports source-code availability and marks asset consent and sensible-content items as addressed.Source: CHECKLIST · p. 11, A CODE AND REPRODUCIBILITY · p. 12

B.1 Latent Gaussian Models

  • The appendix defines LGMs by Gaussian latent fields, conditionally independent likelihoods, and hyperpriors; LGCPs use a Gaussian linear predictor for stochastic intensity, and Laplace approximation is normal around the mode.Sparse precision in the LGM supports INLA inference.Source: §B.1 · p. 12, §B.3 · p. 13, §B.4 · p. 14

B.2 Gaussian Markov Random Fields

  • GMRF conditional independence is encoded by zeros in the precision matrix, enabling sparse Cholesky factorization with typical lattice costs O(n) in one dimension, O(n^3/2) in two, and O(n^2) in three.Fill-reducing orderings can preserve sparsity in the factorization.Source: §B.2 · p. 13

B.7 The Stochastic Partial Differential Equation Approach

  • For point-referenced spatial data, SPDE representation replaces dense GRF inference (O(n^3) time, O(n^2) memory) with a mesh-based Matérn GMRF having sparse precision; temporal data use a one-dimensional mesh.Finite-element basis functions have local support, yielding sparse precision; meshes can be locally refined for irregular boundaries.Source: §B.5 · p. 14, §B.7 · p. 15

Meaning of the parameters

  • For the Matérn field, σ² sets marginal variance, ν controls smoothness, and range ρ describes correlation decay; PC priors are placed on range and marginal standard deviation.The range prior can use a fraction of domain size, while the standard-deviation prior is chosen on the linear-predictor scale.Source: §B.6 · p. 14, Meaning of the parameters · pp. 14–15

C.1 Reformulation of the Latent Gaussian Model

  • The modern INLA formulation defines linear predictors as a sparse design-matrix transformation of the latent field, removing n predictors and improving scalability and numerical stability in data-rich settings.Conditional on the predictors, observations are independent.Source: §B.1 · p. 12, §C.1 · p. 16

C.2 Approximating ˜p(θ | y)

  • INLA approximates the hyperparameter posterior using a Gaussian conditional approximation at its mode, integrates over hyperparameters for marginals, and recovers linear-predictor marginals using a sparse selected inverse.In the modern formulation, η = Ax, so predictor marginals are recovered rather than computed directly as latent-field marginals.Source: §C.2 · p. 16, §C.3 · p. 16, §C.4 · pp. 16–17

C.5 Low-rank Variational Bayes Correction

  • For non-Gaussian likelihoods, a low-rank variational Bayes correction addresses biased Gaussian-approximation posterior means by retaining QX(θ) and modifying only the mean (Van Niekerk and Rue, 2024).The correction uses a restricted Gaussian variational family.Source: §C.5 · p. 17
  • The low-rank correction optimizes a variational objective using a Taylor approximation and Gauss–Hermite quadrature, achieving accuracy comparable to full nested Laplace at negligible added cost.The objective combines expected negative log-likelihood under the Gaussian approximation with a quadratic penalty from the prior precision; the resulting mean improves marginal posteriors while retaining INLA efficiency.Source: §C.5 · p. 17

D KEY PROPERTIES AND USAGE OF PC PRIORS

  • PC priors shrink toward a base model, remain robust with weak or noisy data, and matter in early BO when data are scarce; INLA hyperparameter posteriors support practical tuning and checks.They are also invariant under one-to-one reparameterization and are presented as default hyperpriors for the latent Gaussian model block.Source: D KEY PROPERTIES AND USAGE OF PC PRIORS · pp. 17–18

E PROOF OF THEOREM 1

The proof establishes that maximizing the latent log-intensity and maximizing the intensity select the same locations. It also characterizes the transformed posterior uncertainty and motivates defining acquisition on the latent Gaussian scale.

E PROOF OF THEOREM 1

  • The exponential transformation preserves the complete ordering of candidate locations, so the full set of maximizers is the same for latent log-intensity and intensity.The argument also allows multiple maximizers.Source: E PROOF OF THEOREM 1 · p. 18
  • A Gaussian approximation for latent log-intensity induces an approximately log-normal intensity posterior, whose scale increases nonlinearly with the latent variance.Source: E PROOF OF THEOREM 1 · p. 19
  • The preserved maximizer ordering and transformed uncertainty motivate defining the UCB-style acquisition directly on the latent Gaussian scale.Source: E PROOF OF THEOREM 1 · p. 19

F PROOF OF THEOREM 2

With a fixed SPDE mesh and model specification, modern INLA keeps its principal sparse system at fixed dimension as events accumulate. The mesh-based representation explains why the dominant refitting operations do not grow with the revealed event count.

F PROOF OF THEOREM 2

  • Additional events change the numerical entries of the system but do not enlarge the latent field, because observations project onto the same mesh-based field.Local support of the finite-element basis functions limits each observation to basis functions in its mesh element.Source: F PROOF OF THEOREM 2 · p. 20
  • The dominant sparse Cholesky factorization and linear solves have complexity independent of the increasing number of revealed events under the fixed-mesh formulation.The sparse system retains fixed dimension and mesh-controlled sparsity at each BO iteration.Source: F PROOF OF THEOREM 2 · p. 20
  • Likelihood evaluation and assembling AᵀDₜAₜ may still depend on the event count, even though they do not increase the principal sparse system’s dimension.Source: F PROOF OF THEOREM 2 · p. 20

G.1 Data Preparations

The experiments use reproducible synthetic data and prepared Washington, D.C. crime and Porto taxi datasets, with choices about integration domains and query windows tailored to the task. For temporal crime data, a padded mesh and PC priors support a smooth latent intensity whose posterior indicates sustained rather than isolated changes.

G.1 Data Preparations

  • The synthetic process has 596 events and a known multimodal reference intensity; real-data experiments use May–June firearm incidents and sampled Porto taxi trips.The Washington, D.C. subset contains 342 incidents, and Porto trips are sampled with a fixed seed because the baseline subset was unavailable.Source: §G.1 · pp. 20–21

G.2 Experimental Note on IPS usage.

  • The temporal integration point set alone receives a 5% padding because the experiments report larger edge effects in time than in space.A spatial integration point set could also be added for tasks with strong spatial edge effects.Source: §G.2 · p. 21

G.3 Experimental Note on Query Radius.

  • The query radius is kept moderate so each query reveals only a limited domain portion, preserving sequential decision-making rather than making optimization trivial.Its practical scale should match the application's feasible observation window and resolution.Source: §G.3 · p. 21

G.4 Experimental Note on Computing Resources

  • All methods run locally on a MacBook Air M2 with 16 GB RAM, without HPC resources, to evaluate scalability on commodity hardware.The paper frames this setting as relevant to resource-constrained applications such as UAV sensing.Source: §G.4 · p. 21

H 1D TEMPORAL INTENSITY ESTIMATION (WASHINGTON, D.C. CRIME)

  • The temporal crime model uses a 60-day axis and a 300-node mesh padded by 10%, with a separately padded integration mesh for quadrature.The mesh discretization uses piecewise-linear basis functions to yield a sparse GMRF for INLA inference.Source: H 1D TEMPORAL INTENSITY ESTIMATION… · p. 22
  • PC priors favor longer Matérn ranges and smaller marginal variability on the linear-predictor scale, discouraging short-scale oscillations while retaining sensitivity to sustained changes.The prior parameters are the range ρ and marginal standard deviation σ.Source: H 1D TEMPORAL INTENSITY ESTIMATION… · p. 22
  • The posterior range peaks around 14–16 days and the variance mode is near 0.1, supporting a smooth field with moderate amplitude.The Matérn correlation remains noticeable to about 8–10 days and is close to zero by 25 days; this smoothness reduces overfitting to isolated events while retaining sensitivity to sustained changes.Source: H 1D TEMPORAL INTENSITY ESTIMATION… · p. 23

I BAYESIAN OPTIMIZATION OVER 1D SYNTHETIC INTENSITY

The one-dimensional synthetic experiment evaluates sequential optimization against a known intensity curve, using shared initial regions and repeated seeds. ISBO first explores broadly, then identifies the main peaks and stabilizes its posterior estimate as uncertainty contracts.

I BAYESIAN OPTIMIZATION OVER 1D SYNTHETIC INTENSITY

  • The experiment uses a reproducible synthetic point process with known ground truth, allowing region-based simple regret to be computed directly.All methods share three initial centres within each of ten seeds.Source: I BAYESIAN OPTIMIZATION OVER 1D… · pp. 23–24
  • The synthetic LGCP uses an intercept and a one-dimensional Matérn SPDE field, with empirical Bayes inference and predictions on a common candidate grid.Previously revealed regions and query centres are masked; the fixed-UCB ablation sets κ_t = 6.Source: I BAYESIAN OPTIMIZATION OVER 1D… · pp. 23–24
  • The ablations compare ISBO with a standard Matérn prior, fixed rather than time-varying UCB, and random querying.Each configuration changes a specified component while retaining the others; random query centres are sampled without replacement.Source: I BAYESIAN OPTIMIZATION OVER 1D… · p. 24
  • By iteration 5, ISBO has identified the main peaks, and by iteration 13 its posterior mean closely tracks the ground truth with substantially contracted uncertainty.The sequence illustrates broad early exploration followed by increasingly concentrated sampling of high-intensity regions.Source: I BAYESIAN OPTIMIZATION OVER 1D… · p. 24, Figure 9 · p. 25

J PSEUDOCODE FOR 1D TEMPORAL BO VIA INLA–SPDE–LGCP

  • Table 1 presents the detailed one-dimensional implementation of ISBO.The supplied table caption identifies its scope but does not specify the algorithm steps.Source: Table 1 · p. 25

K 2D SPATIAL INTENSITY ESTIMATION (PORTO TAXI TRAJECTORIES)

The Porto spatial analysis fits a Matérn-SPDE LGCP to taxi starting locations and predicts both latent and response-scale intensity. It also benchmarks inference procedures across event-set sizes, including subsets up to 55,000 events.

K 2D SPATIAL INTENSITY ESTIMATION (PORTO TAXI TRAJECTORIES)

  • The Porto data contain 106,732 trips, with spatial modeling and prediction restricted to a rescaled unit-square domain.The visualization fits use subsets of 200 and 66,666 events.Source: K 2D SPATIAL INTENSITY ESTIMATION… · p. 25

0.01. The latent spatial field follows a Mat´ern–SPDE model with PC priors

  • The spatial LGCP uses an intercept and Matérn field, with empirical Bayes inference and predictions on a pixel grid for both η(s) and λ(s) = exp{η(s)}.The 200-event Nyström baseline fit is shown for comparison; Figure 12 presents posterior mean surfaces at 200 and 66,666 events.Source: §0.01 · pp. 25–26, Figure 11 · p. 26, Figure 10 · p. 26
  • The runtime benchmark uses common event subsets from 200 to 2,500 for four methods and from 5,000 to 55,000 for three methods.All methods receive the same sampled events at each size, and runtimes are mean fitting times across repeated runs on a local MacBook Air M2 with 16 GB RAM.Source: §0.01 · p. 26

L BAYESIAN OPTIMIZATION OVER 2D SPATIAL INTENSITY

Across ten seeds, spatial simple regret decreases over 25 iterations as ISBO approaches high-intensity regions of a full-data reference surface. The visual sequence likewise shows the hotspot structure emerging and the estimated surface sharpening over successive queries.

L BAYESIAN OPTIMIZATION OVER 2D SPATIAL INTENSITY

  • At each iteration, the spatial LGCP is refitted using currently revealed events, then previously queried disks are masked before selecting the next location.Candidate locations are spatial pixels generated from the mesh and restricted to the domain.Source: L BAYESIAN OPTIMIZATION OVER 2D… · p. 27
  • The reference surface is the LGCP intensity fitted to the full dataset, and each queried disk is valued by its maximum reference intensity.Simple and cumulative regret substitute this reference for unavailable true spatial intensity.Source: L BAYESIAN OPTIMIZATION OVER 2D… · pp. 27–28
  • Average simple regret decreases steadily over 25 iterations, bringing queried regions closer to the maximum of the full-data reference intensity.The curves average ten seeds and show 95% Student-t confidence intervals; cumulative regret rises because it aggregates regret across iterations.Source: L BAYESIAN OPTIMIZATION OVER 2D… · p. 28, Figure 13 · p. 28
  • By around iteration 11, the main hotspot structure has emerged, and by iteration 25 the dominant high-intensity areas are clearly identified.The posterior surface sharpens as more regions are queried, consistent with the regret curves.Source: L BAYESIAN OPTIMIZATION OVER 2D… · p. 28, Figure 14 · p. 28

M PSEUDOCODE FOR 2D SPATIAL BO VIA INLA–SPDE–LGCP

  • Table 2 presents the detailed two-dimensional spatial implementation of ISBO.The supplied table caption identifies its scope but does not specify the algorithm steps.Source: Table 2 · p. 29
Loading 2610.12213v1…