A Tutorial Review of Bayesian Optimization with Gaussian Processes to Accelerate Stationary Point Searches

Rohit Goswami

stat.MLMar 11, 2026

A Tutorial Review of Bayesian Optimization with Gaussian Processes to Accelerate Stationary Point Searches

Authors: Rohit Goswami

Organizations: Institute IMX and Lab-COSMO, École polytechnique fédérale de Lausanne · Institute IMX and Lab-COSMO, École polytechnique fédérale de Lausanne (EPFL), Station 12, CH-1015 Lausanne, Switzerland

Abstract

Building local surrogates to accelerate stationary point searches on potential energy surfaces spans decades of effort. Done correctly, surrogates can reduce the number of expensive electronic structure evaluations by roughly an order of magnitude while preserving the accuracy of the underlying theory, with the gain depending on oracle cost, search distance, and the availability of analytical forces. We present a unified Bayesian optimization view of minimization, single-point saddle searches, and double-ended path searches: all three share one six-step surrogate loop and differ only in the inner optimization target and the acquisition criterion. The framework uses Gaussian process regression with derivative observations, inverse-distance kernels, and active learning, and we develop optional extensions for production use, including farthest-point sampling with the Earth Mover's Distance, MAP regularization, an adaptive trust radius, and random Fourier features for scaling. Accompanying pedagogical Rust code demonstrates that all three applications use the same Bayesian optimization loop, bridging the gap between theoretical formulation and practical execution.

Figures & tables

Figure 1: Geometry of the dimer and Householder reflection. The dimer pair (

\mathbf{R}_{1},\mathbf{R}_{2}

) straddles the midpoint

\mathbf{R}_{0}

with axis

\hat{\mathbf{N}}

. The true force

\mathbf{F}

(blue) is reflected about the hyperplane perpendicular to

\hat{\mathbf{N}}

, producing the modified force

\mathbf{F}^{\dagger}

(coral) that climbs along the minimum mode while relaxing perpendicular to it.

Figure 2: GP conditioning in three panels. ( Left ) Before any data, the prior (Eq. 3.1 ) admits a wide family of smooth functions around a chosen mean function. ( Center ) Oracle evaluations supply energies and forces at selected configurations. ( Right ) Conditioning on the data collapses the posterior near training points while preserving wide uncertainty elsewhere; the posterior mean serves as the surrogate surface

V_{\text{GP}}

. The optional nonzero prior-mean branch shown in the schematic is one later extension point in the design space, not the main implementation thread of this tutorial.

Figure 3: GP surrogate fidelity as a function of training set size on the Muller-Brown surface. Each panel shows the GP posterior mean contours after training on

N=3,8,15,30

Latin hypercube-sampled configurations (white markers). With three points the surrogate captures only crude basin structure; by 30 points the contours closely match the true PES (Figure 12 ) in the sampled region.

Figure 4: GP predictive variance on the Muller-Brown surface after 20 training evaluations clustered near minimum A and saddle S1 (black dots). The variance is near zero close to training data and grows with distance, reaching a maximum (coral diamond) in the unexplored region. This landscape is a map of sampling density in the kernel geometry, not of accuracy against the true surface: it tells the active-learning loop where the surrogate has seen least data, which is the quantity an acquisition rule needs, but a low-variance prediction does not by itself imply a small error.

Panel of Figure 5 — Figure 5: Block structure of the full covariance matrix $K_{\text{full}}$ . The base kernel in feature space generates four Cartesian-space blocks through differentiation via the feature Jacobian J. Darker shading indicates higher computational cost.

Figure 6: The inverse-distance feature map. Cartesian coordinates (

\mathbb{R}^{3N}

, not invariant) are mapped to pairwise inverse distances (

N(N{-}1)/2

features, invariant to rotation and translation). The SE kernel operates in this feature space. The Jacobian

\mathbf{J}=\partial\boldsymbol{\phi}/\partial\mathbf{x}

propagates through the kernel via the chain rule to produce the derivative blocks needed for force predictions.

Figure 7: Hyperparameter sensitivity on the Muller-Brown surface. Each panel shows a 1D slice at

y=0.5

with the true PES (black dashed), the GP posterior mean (teal), and the

\pm 2\sigma

confidence band (light blue), for nine combinations of length scale

\ell\in\{0.05,0.3,2.0\}

(columns) and signal variance

\sigma_{f}\in\{0.1,1.0,100.0\}

(rows). Small

\ell

produces noisy interpolation; large

\ell

over-smooths and misses barrier structure. The center cell (

\ell=0.3

\sigma_{f}=1.0

) shows well-calibrated behavior where the confidence band tightly encloses the true surface near training data and widens appropriately in data-sparse regions.

Figure 8: Visual overview of the Bayesian surrogate loop (Algorithm 4 ). Numbered steps proceed clockwise: (1) train the GP, (2) optimize on the surrogate, (3) check trust constraints, (4) evaluate the oracle, (5) select the next query point, (6) update the training set. The oracle (coral) is the only expensive step; all others operate on the cheap surrogate. Method-specific annotations indicate how each algorithm instantiates the inner optimization and acquisition steps. The prior-mean box marks one broader design-space extension that can be studied within the same outer loop without becoming the focus of the present tutorial.

Step	Minimization	GP-dimer	GP-NEB
1. Select subset	Global FPS	Global FPS	Per-bead FPS
2. Train $\boldsymbol{\theta}$	SCG on MAP	SCG on MAP	SCG on MAP
3. Build model	Exact/RFF	Exact/RFF	Exact/RFF
4. Inner optimization	L-BFGS	CG rotate + L-BFGS	NEB relaxation
5. Trust clip	EMD/Euclidean	EMD	EMD
6. Acquisition	Implicit	Implicit	MaxVariance/UCB (OIE) or exhaustive (AIE)

Table 1: Instantiation of Algorithm 4 across methods

Method	Mode	Selection Criterion	Calls/iter
GP-minimization	Implicit	Trust-clipped step	1
GP-dimer	Implicit	Trust-clipped step	1
GP-NEB OIE	Explicit	MaxVariance or UCB	1
GP-NEB AIE	Exhaustive	All images	$P$

Table 2: Acquisition criteria across GP methods. GP-minimization and GP-dimer use implicit acquisition (trust-clipped step from inner loop), GP-NEB OIE uses explicit one-image selection from unevaluated images (MaxVariance or UCB; UCB is the default in chemgp-core ), and GP-NEB AIE uses exhaustive evaluation of all images.

Figure 9: Convergence comparison of the standard dimer, GP-dimer, and OT-GP dimer (OTGPD) on a molecular system (C 3H 5 allyl radical, 8 atoms, 24 DOF) via eOn serve mode. The maximum per-atom force (eV/Å) is plotted against oracle evaluations on a logarithmic scale. The OTGPD variant reaches the convergence threshold (gray dashed, 0.1 eV/Å) with the fewest oracle calls; the GP-dimer shows oscillations from surrogate retraining instabilities that the OT-GP extensions suppress.

Figure 10: Trust region mechanism. ( Left ) A GP-proposed step (coral) that exceeds the trust boundary

d_{\text{EMD}}\leq\Theta

is clipped to the boundary; the oracle evaluates at the clipped location. ( Right ) The trust radius grows with accumulated data via an exponential saturation curve (

\Theta_{\text{earned}}

), capped by a system-size-dependent physical ceiling (

\Theta_{\text{phys}}

Figure 11: Trust region mechanism on a 1D slice (

y=0.5

) of the Muller-Brown surface. The GP posterior mean (teal) and

\pm 2\sigma

confidence band (light blue) are accurate near the training data (black dots) but diverge from the true surface (black dashed) outside the trust boundaries (magenta dotted verticals). A hypothetical GP-proposed step at

x=1.0

(coral cross, labeled "GP step") falls outside the trust region, where the surrogate is unreliable; the algorithm instead evaluates the true PES at the trust boundary (teal star, "Oracle fallback").

Figure 12: Muller-Brown potential energy surface with NEB path overlay. Filled contours show the energy landscape with three local minima (A, B, C) and two saddle points (S1, S2). Eleven NEB images (coral circles, numbered) trace the minimum energy path from A to B through S2. The climbing image (highest-energy interior image) approximates the saddle point. Energy values are reported in the conventional Muller-Brown reduced units and clipped to the range

[-200,50]

for visualization.

Figure 13: LEPS potential energy surface with NEB path overlay. The collinear atom transfer reaction

A+BC\to AB+C

is plotted as a function of bond distances

r_{AB}

and

r_{BC}

. Seven interior NEB images and two fixed endpoints (nine path points, coral circles numbered; endpoints also marked by yellow stars) are optimized in the full 9-dimensional coordinate space and projected onto the

(r_{AB},r_{BC})

plane. The climbing image converges to the saddle region near

r_{AB}\approx r_{BC}\approx 1.0

Å. Contour spacing is 0.5 eV; energies clipped to

[-5,5]

eV.

Figure 14: Convergence of NEB variants on the LEPS surface. Maximum per-atom force versus oracle evaluations on a logarithmic scale. Standard NEB (156 calls), AIE (100 calls), and OIE (42 calls) all reach the convergence threshold (dashed, 0.1 eV/Å). The OIE variant evaluates only the highest-variance image per cycle (Eq. 4.3 ) and converges fastest.

Figure 15: Convergence comparison of the GP-minimizer and classical L-BFGS on the LEPS surface. With a force convergence threshold of

10^{-2}

eV/Å, the GP surrogate reaches the threshold in 9 oracle calls, compared with 57 for direct L-BFGS on the same starting configuration. Force values plotted on a logarithmic scale.

Panel of Figure 16 — Figure 16: Decision flow of the OT-GP framework. The training pipeline (FPS, MAP regularization) and adaptive trust region (EMD-based) address the failure modes of the basic GP-dimer. The optional prior-mean branch shown here is included to place later extensions in the same design space, not to redefine the main algorithmic thread discussed in the text.

Figure 17: ( Left ) Farthest point sampling selects a geometrically spread-out subset (coral) from the full training set (gray) by greedily maximizing the minimum distance to the existing selection. ( Right ) The Earth mover’s distance measures structural dissimilarity between two molecular configurations as the optimal transport cost of matching their atom-pair distance distributions; it is invariant to rotation and atom indexing.

Figure 18: Farthest point sampling (FPS) on the LEPS surface. Approximately 50 candidate configurations from a GP-NEB run are projected onto their first two principal components in inverse-distance feature space. FPS selects 20 points (teal diamonds) that maximally cover the feature space; pruned points (gray circles) lie near already-selected configurations and would add redundancy to the training set without improving kernel matrix conditioning.

Figure 19: RFF approximation quality on the LEPS surface (H 3 , 3 inverse-distance features). Energy MAE (top) and gradient MAE (bottom) between RFF and exact GP predictions on held-out test points, plotted against

D_{\text{rff}}

. The kernel is well approximated by

D_{\text{rff}}\sim 100

Figure 20: Convergence of GP-accelerated minimization on a real molecular system (PET-MAD potential). The maximum per-atom force is plotted against the number of oracle evaluations on a logarithmic scale.

Figure 21: Illustrative convergence comparison of GP-NEB variants on a 9-atom cycloaddition (PET-MAD surface, 27 DOF). Climbing-image force is plotted against oracle evaluations on a logarithmic scale for the classical NEB, AIE, and OIE update patterns. The figure is used here to compare the qualitative behavior of the three outer-loop choices on the same molecular path; the tighter CI-targeted literature-style benchmark is summarized separately in the Supporting Information.

Figure 22: Energy profiles along the converged MEP for the cycloaddition system on the PET-MAD surface. All three NEB variants recover the same barrier and exothermic product basin. The AIE profile is slightly shifted near the saddle region but converges to the same endpoints.

Figure 23: Reaction valley projection [ 94 ] of the NEB paths on the PET-MAD surface. Reaction progress $s$ and orthogonal deviation $d$ are projected RMSD coordinates. Standard NEB (blue stars) and GP-NEB AIE (orange stars) saddle points are compared; both lie near the true saddle (black square). Dashed contours show GP variance ( $\sigma^{2}$ ); the paths remain within the well-sampled region where the surrogate has accumulated data, which in combination with the per-step trust region is what keeps them close to the true saddle. Bottom: molecular snapshots at key points along the reaction coordinate.

Appendix figures & tables5 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 26: MAP-regularized negative log marginal likelihood landscape on the LEPS surface. The filled contours show the NLL as a function of the log-hyperparameters

\ln\sigma^{2}

and

\ln\theta

, computed on a 40

\times

40 grid from 5 training points near the reactant. The dashed black contours show the gradient norm. The coral star marks the MAP optimum that SCG converges to. Regions where the Cholesky factorization fails (hyperparameters that produce a non-positive-definite covariance matrix) are masked.

Parameter	Default	Meaning
Kernel
$\sigma_{c}^{2}$	1.0	Constant offset (eV-scale energies); set to 0.0 for centered model surfaces
$\sigma_{f}^{2}$	Free (MLL)	Signal variance (optimized via MAP-NLL)
$l_{\phi(i,j)}$	Free (MLL)	Inverse-distance length scales (per atom-type pair, optimized)
Noise and Regularization
$\sigma_{E}^{2}$	$10^{-8}$	Energy noise (Tikhonov regularizer, not physical noise)

Appendix

Table 4: Default hyperparameters for GP-accelerated saddle point searches. Parameters marked with † are documented in [ 26 ] .

Concept	chemgp-core (Rust)	gpr_optim (C++)
Inverse-distance kernel	kernel.rs ( MolInvDistSE )	InvDistSE
Derivative kernel blocks	kernel.rs ( molinvdist_kernel_blocks )	GPKernelBlocks
GP training (SCG)	train.rs + scg.rs + nll.rs	GPModel
Dimer method	dimer.rs ( gp_dimer )	DimerSearch
NEB method	neb.rs + neb_oie.rs	NEBSearch
OT-GP Dimer	otgpd.rs ( otgpd )	OTGPDimer

Appendix

Table 5: Conceptual mapping between mathematical ideas, chemgp-core Rust modules, and gpr_optim C++ classes. Both implementations share the same algorithmic structure; the Rust code runs the benchmarks reported in this review.

Equation	chemgp-core Function	Purpose
Eq. 2.6	curvature()	Dimer curvature estimate
Eq. 2.7	rotational_force()	CG rotation force
Eq. 2.8	translational_force()	Householder reflection
Eq. 2.10	neb_force()	NEB total force
Eq. 3.6 – 3.7	molinvdist_kernel_blocks()	Derivative covariance blocks
Eq. 3.8	build_full_covariance()	Assemble $\mathbf{K}_{\text{full}}$

Appendix

Table 6: Equation-to-function mapping for chemgp-core. Each equation in the main text maps to a specific function in the Rust implementation. The annotations in the code listings connect variable names to the mathematical symbols in the equations.

Task	System	Variants	Baseline	Variant A	Variant B
Minimization	LEPS (9D)	gp_minimize vs L-BFGS	57	9	-
Minimization	Muller-Brown	gp_minimize vs L-BFGS	34	50	-
Minimization	PET-MAD molecule	gp_minimize vs L-BFGS	45	13	-
Single saddle	Muller-Brown	gp_dimer vs std. dimer	6	7	-
Single saddle	LEPS (9D)	gp_dimer / OTGPD vs dimer	10	8	8
Single saddle	d000 metatomic molecular run	standard dimer vs OTGPD	50	-	13

Appendix

Table 7: Oracle-call counts to convergence for the illustrative benchmarks used in this review. Each row corresponds to a chemgp-core example whose execution trace drives a figure in the main text or a literature-aligned molecular benchmark in the accompanying ChemGP harness. "Baseline" denotes the classical reference algorithm for that task. "Variant A" and "Variant B" are method-specific comparison columns so the table can mix minimization, dimer, and CI-targeted path-search examples without forcing them into one GP/OT-GP taxonomy. Smaller is better.

Explore similar work

Adaptive Candidate Point Thompson Sampling for High-Dimensional Bayesian Optimization
Apr 10, 2026Donney Fan, Geoff PleissBayesian Optimization Thompson Sampling
University of British Columbia Vector Institute
Transfer Learning Architectures for Scalable Multi-Fidelity Bayesian Optimization
Jul 26, 2026Jaewook Lee, Ethan Errington, Christian D. Lorenz +1Multi-Fidelity Optimization Surrogate Modeling
Department of Engineering, King’s College London, London, WC2R 2LS, United Kingdom.
Frugal Bayesian Optimization: Scalable Surrogates for Data- and Resource-Limited Discovery
Jul 31, 2026Panagiotis Krokidas, Christoforos Rekatsinas, Vassilis Sioros +3Surrogate Modeling Bayesian Optimization
Institute of Informatics and Telecommunications, National Centre for Scientific Research "Demokritos", Agia Paraskevi, Greece. · Department of Mechanical Engineering and Aeronautics, University of Patras, Patras Greece. · Department of Informatics and Telecommunications, National and Kapodistrian University of Athens, Athens, Greece. +4

Apr 10, 2026J/K move · Enter open · S save

Adaptive Candidate Point Thompson Sampling for High-Dimensional Bayesian Optimization

Donney Fan, Geoff Pleiss

University of British Columbia Vector Institute

Bayesian Optimization Thompson Sampling Adaptive Sampling

Open paper AskPDF arXiv

In Bayesian optimization, Thompson sampling selects the evaluation point by sampling from the posterior distribution over the objective function maximizer. Because this sampling problem is intractable for Gaussian process (GP) surrogates, the posterior distribution is typically restricted to fixed discretizations (i.e., candidate points) that become exponentially sparse as dimensionality increases. While previous works aim to increase candidate point density through scalable GP approximations, our orthogonal approach increases density by adaptively reducing the search space during sampling. Specifically, we introduce Adaptive Candidate Thompson Sampling (ACTS), which generates candidate points in subspaces guided by the gradient of a surrogate model sample. ACTS is a simple drop-in replacement for existing TS methods -- including those that use trust regions or other local approximations -- producing better samples of maxima and improved optimization across synthetic and real-world benchmarks.

Abstract

Explore similar work

Adaptive Candidate Point Thompson Sampling for High-Dimensional Bayesian Optimization

Transfer Learning Architectures for Scalable Multi-Fidelity Bayesian Optimization

Frugal Bayesian Optimization: Scalable Surrogates for Data- and Resource-Limited Discovery