BoGrape: Bayesian optimization over graphs with shortest-path encoded
Authors: Yilin Xie
Shiqiang Zhang
Jixiang Qing
Ruth Misener
Calvin Tsay
Venue: ICLR 2026 Poster (2026)
Keywords: Bayesian optimization
graph optimization
mixed-integer programming
shortest-path
Primary area: optimization
TL;DR: We introduce Bayesian optimization for functions over general graph spaces by encoding shortest-paths in mixed-integer programming.

Abstract: Graph-structured data are central to many scientific and industrial applications where the goal is to optimize expensive black-box objectives defined over graph structures or node configurations---as seen in molecular design, supply chains, and sensor placement. Bayesian optimization offers a principled approach for such settings, but existing methods largely focus on functions defined over nodes of a fixed graph. Moreover, graph optimization is often approached heuristically, and it remains unclear how to systematically incorporate structural constraints into BO. To address these gaps, we build on shortest-path graph kernels to develop a principled framework for acquisition optimization over unseen graph structures and associated node attributes. Through a novel formulation based on mixed-integer programming, we enable global exploration of the combinatorial domain over graph structures and explicit embedding of problem-specific constraints. We demonstrate that our method, BoGrape, is competitive both on general synthetic benchmarks and representative molecular design case studies with application-specific constraints.

Published as a conference paper at ICLR 2026
BOGRAPE: BAYESIAN OPTIMIZATION OVER GRAPHS WITH SHORTEST-PATH ENCODED
Yilin Xie & &hiqiian Zhang*
Department of Computing, Imperial College London, UK
Jyilin. xie22, zhang21 @imperii allaa.aa. .ac.uk
Jixiang Qing
Department of Computing, Imperial College London, UK School Mathematical Sciences, Lancaster University, UK j. g.qinge iiteerippial.. ac uk
Ruth Misener & Calvin Tsay
Department of Computing, Imperial College London, UK {r rr.mssnne,, c c.tsay) simperiall. awmprraalla....m aa..c.ukkkkk
ABSTRACT
Graph-structured data are central to many scientific and industrial applications where the goal is to optimize expensive black-box objectives defined over graph structures node configurations as seen in molecular design, supply chains, and sensor placement. Bayesian optimization offers a principled approach for such settings, but existing methods largely focus on functions defined over nodes of fixed graph. Moreover, graph optimization is often approached heuristically, and it remains unclear how to systenatically incorporate structural constraints into BO. To address these gaps, we build on shortest- -ooot graph kernels to develop a principled framework for acquisition optimization over unseen graph structures and associated node attributes. Through novel formulation based on mixed-intege programming, we enable global exploration of the combinatorial domain over graph structures and explicit embedding of problem-specifie constraints. We demonstrate that our method, BoGrape, is competitive both on general synthetic benchmarks and representative molecular design studies with inpliiatiionspeciine constraints.
INTRODUCTION
Graph-structured data are playing an emerging role across scientific and industrial fields, giving rise to series of decision decision-making problems graph domains, such graph-based molecular design (Korovina et al.] 22200] Mercado et al. 2021] Yang et al.) 2024) and neural architecture search (Elsken) ct al. 2019; W hite et al. (023). Broadly speaking, there two classes of graph optimization problems (Wan et al. 2023): : (i) optimizing over nodes, with given (unknown) graph as the search nodes as the objective, and (ii) optimizing over graphs, with the entire (constrained) graph domain as the search space and function over graphs as the objective. The latter case, which this work studies, usually more challenging since the graph structure itself is optimized, resulting in a complicated combinatorial optimization task.
space and a function
For both aforementioned scenarios, the objective function can black-box. and, when expensive to evaluate, discourages gradient- and population- based methods. These characteristics motivate several works to extend Bayesian optimization (BO) (Frazier 2018
Garneti 2023) to graph domains (Cul & Yangl 2018: Oh et al.) 2019 Wan et al.) 2023 Liang et al. 2024) given its potential sample efficiency. B0 relies on two main components: surrogate model, e.g., Gaussian processes (GPs), trained on available data approximate the underlying function, and an acquisition function used
"Equal contributions.
1 Equal contributions. Corresponding authors.

---

Published as a conference paper at ICLR 2026
Graph Kernel
Graph Gaussian Process (xx()),
Acquisition Optimization min LCB
x = (G),Pl)
x = FG)
~ f(x)~
Kax Kxx
LCB
xxx)
ko G*)
S KEx
Kex Kidk Rxx
X
1ption
Ksp
u.s.
K
MIP encoding
Option 2
Kssp
800 V.S.
K
Next query point
Figure 1: Key components of BoGrape The graph kernel comprises kg and kp on the graph and feature levels, resp. The graph GP is subsequently trained using the chosen kernel, and its posterior is used to build the acquisition function, e.g., LCB. Note that graph GP includes discrete graph domains; the continuous domain is only for illustration purposes. Acquisition optimization is formulated and solved as a MIP using the encoding of shortest paths and graph kernels, giving the next query point.
LCB.
to suggest the next sample. To translate BO to graph domains, one needs a surrogate model over graph inputs with suitable uncertainty quantification, leading existing approaches to adapt GPs with various graph kernels (Ramachandram et al. [017]] Borovitskiy et al., |021 Ru et al. 2021 Zhi et al.,
2023). However, a general graph BO framework IS missing, since existing works either (1) limit the
searchable graph set given fixed graph (Oh et al. 2019 Wan al. 2023 Liang al. 2024).
directed labeled graphs (Ru et al.] [202] Wan et [l.]1021 White et al. (021). unlabeled graphs (Cui
& Yang 2018), etc. or (ii) rely on task-specific similarity metrics (Kandasamy et al.) 2018).
When optimizing over graphs, the search space includes both continuous and discrete variables, thus limiting the choice of optimization techniques. For example, acquisition function optimization in graph BO is mostly performed using evolutionary algorithms (Kandasamy et al., 2018 Wan
al. 2021) or sampling (Ru al., [021] Wan al. 2023), which are incapable of (i) effectively exploring the search domain, (ii) embedding problem-specific constraints, and (iii) guaranteeing optimality in terms of acquisition function, which is essential for optimization convergence. To mitigate these issues, this paper explores mixed-integer programming (MIP) as alternative to represent an analytic expression of the graph function. The challenges this paper addresses are to manage both the black-box setting and MIP encodings of surrogates for graph BO.
Recent advances on applying MIP to optimize trained machine learning (ML) models (Ceccon et al.. 2022 Schweidtmann et al. 2022 Thebelt et al. 2022b) suggest pathways to address these challenges. By equivalently encoding surrogates, c.g., GPs (Schweidtmann et al. 202T). trees (Misie) 2020
Mistry et al. 2021 Ammari al. 2023) neural networks (NNs) (Fischetti & Jo 2018 Anderson et al. 2020 Isa et al. 2021 Wang et al. 2023). as constraints in larger decision-making problems, several MiP-based BO methods are proposed, allowing global optimization over mixed-feature
domains (Thebelt et al. 2021 2022a Papalexopoulos et al. 20222 Xie et al. (2024). Moreover, some works develop MIP-based techniques to handle optimization problems constrained by graph neural networks (GNNs), with applications to molecular design (Zhang et al. 2023) McDonald et al. 2024 Zhang et al. 2024) and robustness certification (Hojny et al.] 2024 Gaines et al.) 2025). However, given the data requirements of GNNs, the computational cost of solving the large resulting MIPs, and the lack of uncertainty quantification, GNNs are impractical surrogates for graph BO.
This paper proposes BoGrape, a MIP-based graph BO method to optimize functions over connected graphs with attributes. GP surrogates are formulated and optimized using global acquisition function optimization techniques introduced in Xie et al. (2024). We develop four variants of the classic shortest-pant patt graph kernel (Borgwardt & Kriegel 2005), as well as MIP encodings of graph search space, for in BoGrape By introducing representation of the shortest paths decision variables, the acquisition function optimization is formulated as a MIP with a mixed-feature search space, graph kernel, and relevant problem-specific constraints. While we choose shortest path kernels for their
2

---

Published as a conference paper at ICLR 2026
ability to model both undirected and directed graphs, BoGrape focuses on undirected graphs. Figure iillustrttes the BoGrape pipeline. Our main contributions include:
We propose graph representations with their corresponding shortest paths, and theoretically prove that the feasible domain of our formulation is equivalent to the graph space consisting of all connected graphs.
We formulate shortest-path graph kernels, node attribute kernels, and GP posterior in- formation based on our graph encoding as MIP constraints, enabling global acquisition optimization.
We provide a principled BO framework over graph spaces from a discrete optimization viewpoint, BoGrape is compatible with problem-specific constraints over graph structures, node attributes, and their interactions.
2
PRELIMINARIES
2.1 BAYESIAN OPTIMIZATION (BO)
BO Frazier 2018 is derivative-free optimization framework iteratively approach the optimum of an sypeessveeft--valuall black-box function. At each iteration, a surrogate model, usually a GP (Schulz al. 2018), trained on the current observed dataset. With the surrogate constructed and trained, an acquisition function is then formulated based on the posterior information, c.g., probability of improvement (PI) (Kushner) [1644). expected improvement (EI) (Jones al. [1998). lower confidence bound (LCB) (Srinivas et 20100,, predictive Entropy search (PES) (Hernández) Lobato al. 2)14), etc.. Optimizing the acquisition function returns the next query, whose function value is evaluated to form the next data point. This process repeats until meeting a stopping criterion.
2.2 GLOBAL OPTIMIZATION OF ACQUISITION FUNCTIONS
Most theoretical results for regret bounds in BO rely on the global optimization over acquisitions (Srinivas al., 2012), i.e., they assume the global minimizer/maximizer of the acquisition function is found at each step, which may not be satisfied using gradient- and sample-based optimizers, Xie et al. (2024) introduce PK-MIQP, a global acquisition optimization framework based on mixed-integer quadratic programming (MIQP). The core of PK- PK-MIQP is the piecewise linearization of stationary or dot- dot-product kernel, e.g., RBF, Matérn, etc., based on which the acquisition optimization is then formulated as an MIQP. PK-MIQP is useful because of its (i) compatibility with various kernels (note the piecewise linearization is unnecessary if the kernel can be expressed linearly), and (ii) theoretical guarantee on regret bounds. We present its formulation for the LCB acquisition function here:
min - -3/12
+ LCB acquisition
(la)
s.t. = K2X XxXy 0² < Kuz Kax K Kx
+ GP posterior mean
(lb)
+ GP posterior variance
(1c)
K,x: k(x,X X'), VI Sict
kernel function
(1d)
a dEA X
search space
(Ic)
2.3 SHORTEST-PATE GRAPH KERNELS
Graph kernels extend the concept of kernels graph domains and are used to measure the similarity between two graphs Mathematically, a graph kernel k(--)) G x G + R is given by k k(G, G') (o(G), ((G") R4 where Φ : 9 -> H a feature map from graph domain G to a reproducing kernel Hilbert space 7 with inner product (+ (Kriege et al.) 2020). Past research develops graph kernels using a variety of graph patterns, e.g., neighborhoods, subgraphs, walks, paths. We refer the reader to /Vishwanathan et al.) 2010 Borgwardt et al. 2020 Kriege et al. 2020 Nikolentzos et al. 2021 for more details on graph kermels. Several works also use graph kernels to optimize over nodes (Oh
et al. 2019; [Borovitskiy et al., 2021 Wan et al. 2023 Liang al.] (024). but the involved kerness measure the similarity of two nodes on one given graph and do not support optimizing over graphs (see Section [ for this distinction). We focus on the shortest-path (SP) kemel (Borgwardt & Kriegell
3

---

Published as a conference paper at ICLR 2026
2005) in this paper, owing its ability (i) handle both directed and undirected graphs, (ii) consider node labels, and (iii) capture the relationship between non non-adjacen graph nodes, making it more general than kernels based subgraph patterns (Shervashidze et al. 2009 Costa & Grave, 2010). We further discuss the choice of kernels Appendix [B.]] For graph G. denote by as the label of node u, Cu,v as the shortest path from distance from node to (which is unique). The SP kernel between graphs G (V), EE) and G = (V2 (V?,E²) is defined as:
(which may not be unique), and du, as the shortest
ksp(G', C
Σ
(o(laysuug) Ki (1 (2) kolbo,
(2)
where k, is kernel comparing node labels and ke is a kernel comparing path lengths.
3
METHODOLOGY
3.1 VARIANTS OF THE SHORTEST-PATE GRAPH KERNELS
We build on Eq. (2) and develop variants of the thortest-path kernel. Both k, and kg in Eq. (2) are usually chosen as Dirac kernels, giving the explicit representation of the SP kernel as:
(G G²)
Σ
(lm, -lay dui,oy = dug,vas Lor Loz).
(SP)
nin2
WI PVE
where ni normalizing coefficient with n1, f2 as the node numbers of G¹, G2, resp.
Each node may additionally have problem-specife features beyond a single label. From here on, we use X (G, F) to denote an attributed graph with as the underlying labeled graph and F as node features. Intuitively, we can compare the features of two nodes instead of labels in kig. However, this could unnecessarily reduce the number of matching paths between two graphs, as requiring identical node features is restrictive and may introduce additional subgraph information into path comparison. Another option is to use a more complicated kernel ky that measures similarity between features of two nodes, which may significantly increase the computational cost of optimization (similarity is computed for all node pairings). Therefore, we borrow from (Cui Yang, 2018) the idea to separate the implicit and explicit information of graphs, i.e., the kernel value between attributed graphs X1, X hecomes:
A(X = Rok B 3 (FFFF,
where kg any graph kernel, e.g., (SP), kF any kernel over features, and a, 3 are learnable parameters controlling the trade-off between graph similarity and feature similarity.
Since node label is usually included as a node feature and considered in kF term, and comparing labels in Eq. (SP) increases the complexity of our upcoming optimization formulations, we further propose a simplified shortest-pate (SSP) kernel corresponding to an unlabeled SP kernel:
kssp(G' G2)
Σ
1dd,, ddg.vz))
(SSP)
n
n
W.. NO VVV2
Lemma 3.1. SP and SSP kernels are positive definite (PD).
Proof. Borgwardt & Kriegel (2005) prove the SP kernel PD. The SSP kernel is a special case of the SP kernel where all nodes have the same label, hence is also PD.
Observe that both the SP and SSP kernels linear kemnels if pre-compute all shortest paths in each graph and count the number of occurrence for each shortest path length. Such linearity simplifies the optimization step (which still requires the non-trivial representation of shortest paths), but reduces the representation ability of the kernels and limits the maximal rank of the Gram matrix. Motivated by the practically strong performance of exponential kernels such as RBF, Matérn, graph diffusion kernel (Oh et al. 2019), etc., we propose two nonlinear graph kernels based on SP and SSP kernels:
KESP (G', G2) exp(ksp(CC 1//) kESSP (G2, G2) exp(kssp G2)) where variance of is added control the magnitude of kernel value.
(ESP)
(ESSP)
4

---

Published as a conference paper at ICLR 2026
Table 1: List of variables introduced represent the shortest path, where is the number of nodes.
variables
type
description
Au, E {0, 1], u, [x]
binary
the existence of edge from node u to v the length of shortest path from node to
dujo E [x], u, N,D. In]
integer
u,v E 10, u, U, u. [n]
binary
the presence of node w on the shortest path from u to v
Lemma 3.2. ESP and ESSP kernels are PD.
Proof. SP and SSP kernels can be rewritten into linear forms, so ESP and ESSP are exponential kernels, which are known to be PD (Fukumizu 2010).
Remark 3.3. The nonlinear kernels introduce additional difficulties for optimization, discussed later Section 3.4 but may demonstrate better empirical performance compared to their linear counterparts, owing to increased representation ability.
3.2 GLOBAL ACQUISITION FUNCTION OPTIMIZATION
To extend the prior formulation in Eq. to optimize LCB acquisition in graph space (see Appendix fo2 for the applicability other acquisition functions), replace Eq. (Td) by our graph kernels in Section 3.1 and define a combinatorial graph search space for Eq. (Te) as x. (G, F) E = 9 x F. (G,
To maintain consistency with the general BO setting, we denote
the next sample and X ((',, F'),)) yf == as the previous samples the 1-th iteration. The difference that now need to optimize over both the graph domain G E G and the feature domain F E F. W.l.o.g., assume that each node has M features F E [= and the first L features denote the one-hot encoding of its label, i.e. EiEry E = 1, where [n] denotes set {0, 1,..., n - 1}. This modified formulation allows: (i) discrete variables, which is key challenge of graph optimization, (ii) problem-speciff constraints over graph domain, and (iii) theoretical guarantees regret bounds.
binary adjacency matrix is sufficient to represent the graph domain; however, encoding correspond- ing shortest- path information (for unknown graph) is not straightforward and comprises main technical contribution of this work. We first introduce the formulation of shortest paths in Section 3.3
and then explicitly derive Eq. (Idh in Section 3.4 the graph kernels in Section 3.1
3.3 ENCODING OF THE SHORTEST PATHS AS OPTIMIZATION CONSTRAINTS
For the sake of exposition, we first consider all connected graphs G with fixed size, i.e., node number is given (Appendix A.. discusses formulations for graphs of unknown size). Table summarizes the optimization variables. Since our formulations involve constant graph information and their variable counterparts, for each variable Var, we use Var(G) denote value on a given graph G. For example, (G,,,()) is the shortest distance from node node in graph G.
If grgrphh G is given, all variables in Table 1 can computed using classic shortest-pale -attest--pahhaae algorithms, such the WarshaWarshall algorithm (Floyd 1962). In graph optimization tasks, however, need to encode the relationships between these variables as constraints. Motivated by the Floyd-Warshe Warshall algorithm, first present the constraints in Eq. (5) of Appendix A.2 and then prove there exists bijective between the feasible domain given by these constraints and all connected graphs with size n. Here directly give the final encoding of the shortest paths in the following linear MIP (details in
5

---

Published as a conference paper at ICLR 2026
Appendix A.2
Au, du, 0, 6e,0 1(w
Vu, [n]
du 1 1 + n (1- (1-Au,,),
Vu, E [n]] #
du 2 2 Au,o
Vu, U E [n]] u u
du S du + duiv - (1 - 8 0,
Vu, v, [n]]
du 2 du + dw.v 2n (1 8 0,,)),
Vu, v, [ra]
(MIP-SP)
8 u =1,
Vu, E [na], u
Σ
0 2 2 + (n - 2) (1 An,))
Vu, E [n]] u+
wE(m) [v]
Σ wE(n)
8 2 (1 Au,o),
Vu, E [n]], a #
Lemma 3.4. (Au,o(G), (., du (G), (G)) is a feasible solution of Eq. MIP-SP with size n = ntG) given any connected graph G.
Proof. Trivial to verify by definition.
Theorem 3.5. Given any for any feasible solution (Au,es due ou, ow,)) of Eq. (MIP-SP with
size n, there exists a unique graph such that:
(Au, (G), du. (G), &g, (Auar duner &w,u),
i.e. there is bijection between the feasible domain of Eq. (MIP-SPI with size and the set consisting of all connected graphs with n nodes.
The formulation becomes more complicated when the graph size is unknown (but bounded). Denote no and n as the minimal and maximal node numbers, resp., and A,, to represent the existence of node Variables du, and 6%, need to be properly assigned when either u or does not exist, Moreover, we extend the domain of du. from [n] to [n + 1] and use n to denote infinity. Eq. (MIP) SPPpluus Appendix A.3 presents the encoding and Theorem 3.6 Theorem 3.5 unknown size.
Theorem 3.6. There a bijection between the feasible domain of Eq. MIP-SP pPpss with size (No. n] and all connected graphs with number of nodes in [Na, n).
See Appendix A.4 for proofs of Theorems 3.5 and which guarantee the equivalence of our encoding for directed and strong connected graphs. Appendix Hypendix[A.6(shoolk how to further simplify our encoding for undirected graphs. Appendix B.3 discusses the effectiveness our encoding.
3.4 ENCODING OF GRAPH KERNELS AS OPTIMIZATION CONSTRAINTS
We now rewrite Eq.
ld using Eq.
as:
Kxx; k(x, X,)
ko - 1.r-B-rrFFFFF (F,
Given that kr is independent of the choice of graph kernel kig and that kernels continuous features are studied in Xie et al. (2024), here we focus on formulating kg. See Appendix A.5 kernel encoding with binary features.
for respective
Formulating KG ko (G,GG) Gi) is straightforward for SP and SSP kernels:
kssp (G,
Σ Σ
day
Awa.ng
Σ D.(G)) di.v>
n2
n2n u.t,SEim)
[i (2,02E(D(G*))
where == n(Gi) is the node number of Gi, di 1(du,v 5) are indicator variables:
Σ 8 dz. du Vu, E [n].
sE(n+1]
and D. is the number of shortest paths with length in G:: D. = ((u, u [nd], du = =s}l-

---

Published as a conference paper at ICLR 2026
Algorithm 1 BoGrape at /-th iteration.
1: Input: dataset X ((G", y H==i, hyperparameter St, graph kernel
2: Model training: kernel parameters a, B, o
graph GP fit to X
3: Acquisition formulation:
4: represent xx and K e in Eqs.
1b
Section
3.4
5: search space in Eq. (Te)
> problem-specine global optimization
6: Optimization: initialize and solve MIP Eq.
7: Output: proposed sample (G', Ft)
dege is not used evaluating the kernel, since it means the shortest path does not exist. Remark 3.7.
Similarly, introducing indicator variables Pu,0 f2 as:
pe 2 = 1Fuu,, = du,e Fvd = 1), vu, [n], 4.4 E [L].
and counting the numbers of each type of paths in G¹:
P.. (G)) = (( v) 22 [mi], (u())) = ,,, (G') =s, (.())) 1211
the SP kernel formulated as:
ksp(G, n'n d.
Σ
Pali,a(G") pulluie
There are several ways to handle the exponential kernels: (i) directly use (local) nonlinear solvers, losing optimality guarantees, (ii) piecewise linearize the exponential function following Xie al. (2024), or (iii) utilize nonlinear MIP functionalities in established solvers such as Gurobi (Gurobi Optimization, LLC. 2024) SCIP (Vigerske & Gleixner 2018). In our experiments, choose to use Gurobi, which by default employs a dynamic outer approximation of the exponential function given an error tolerance.
is noteworthy that Kgz in Eq. ITc is not constant with a non-stationary kernel, making it the most complicated term in the whole formulation. By definition, kssp(G, G) has quadratic form:
kSsp kssp(G, (G,G)
Σ D. sE(n)
n
where D== Em di Vs E [n]. Reusing the indicator trick and introducing D2 = 1(D, = c). the quadratic form is equivalently linearized as:
1 Σ De
KssP(G,G)= nd
where indicator variables D, vs E |n].cE In² + 1] should satisfy:
Σ D =1
D. Da, Vs [n].
+E(n?+1]
+1]
Repeating the procedure for the SP kernel, we have:
Ksp (G,(G,
1
Σ
2 pe Poti.da
n
where indicator variables
P2 1( c), vs In], 4 [L], e [n²
satisfy:
Σ
P Pa
Σ
P Panits e [n], 4 (T.)
cE[n++1]
With the above graph GP model and optimization encodings, we have presented all the pieces needed to implement end-to-end graph BO procedure. Algorithm Tottlnes BoGrape
7

---

Published as a conference paper at ICLR 2026
SM
WLOA
Ideal Fit
Ideal Fit
Ideal Fit
Ideal Fit Test Data
Test Data
Test Data
Test Data
Ideal Fit
Ideal Fit Test Data
Ideal Fit Test Data
Ideal Fit Test Data
Test Data
RW*
SM+
WL
WLOA
Ideal
Ideal Test Data
Ideal
Ideal
Test Data
Test Data
Test Data
Figure 2: Compare predictive performance of GP with different kernels. indicates linear combination of given kernel and feature kernel. 100 samples are randomly chosen from the QM7 dataset with various graph sizes, 30 of which are used for training. The predictive mean with one standard deviation (predicted y) of the remaining 70 graphs are plotted against their real values (true y).
EXPERIMENTS
All experiments are performed on a 4.2 GHz Intel Core i7-7700K CPU with GB memory. We use GPîlow (Matthews et al. 2017) to implement GP models, GraKel (Siglidis et al. (020) to implement the classic graph kernels, PyG (Fey & Lenssen 2019) to implement GNNs, and Gurobi (Gurobi
Optimization. LLC. 2024 to solve MIPs.
There are few synthetic benchmark functions f F + R for general graph domains since most graph BO works focus specific types of graphs. W.o.lgg, consider GNNs graph functions that maps general labeled connected graph to real values. We conduct experiments considering two settings: (i) randomly initialized GNNs which serves as random synthetic functions, and (ii) GNNs trained on molecular datasets as graph property predictor in real-world case studies. is noteworthy that the architectures of GNNs, the training mechanism, and the choice of datasets not major components of our work, since they merely serve as benchmarks of black-box graph functions.
4.1 MODEL PERFORMANCE
Before conducting the optimization tasks, we first compare the performance of graph GPs with various graph kernels on randomly sampled molecules from the QM7 dataset (Blum & Reymond, 2009 Rupp et al. 2012). Figure 2 visualizes the results and more detailed discussion can be found in Appendix C.1 Figure shows that four shortest-path kernels have comparable prediction accuracy against other classic graph kernels, while the two exponential kernels quantify uncertainty more accurately (also supported by Table Appendix C.T.. Additionally, the introduction of the additional feature kernel term generally improves the predictive performance of all kernels. For larger graph sizes, Table
8

---

Published as a conference paper at ICLR 2026
Random
RW-rann
Wi-rand
-avol
SSP
FSSP
best
240
340
iter
iter
(a) GAT, N 10
(h) GCN, N = 10
(c) GraphSAGE, N 10
0150
240
beryy
2400
340
(d) GAT, N 20
(e) GCN, N = 20
(l) GraphSAGE, N 20
Figure 3: Bayesian optimization results on synthetic benchmarks with N E {10,20}. Best objective value is plotted at each iteration. Mean with 0.5 standard deviation over 10 replications is reported.
of Appendix C.I shows that the more complicated kernels, i.e., SP and ESP, are generally better at predicting graph properties, since they impose stronger criteria on comparing shortest-paths between two graphs.
Remark 4.1. See Appendix C. for similar illustrations of other graph kernels and more evaluations, and Appendix B.4 complexity analysis of different graph kernels.
4.2 OPTIMIZATION OF SYNTHETIC BENCHMARKS
We first evaluate the performance BoGrape over synthetic benchmarks, i.e., arbitrary functions over graphs. Specifically, given their ability universal approximators, we cover the range of black-box graph functions using randomly initialized GNNs including GAT (Veličković et al., 2017), GCN (Kipf & Welling 2017). and GraphSAGE (Hamilton et al., 2017). Each GNN consists of two convolutional layers to learn graph embeddings, and two linear layers. Each hidden layer has 64 features. The search space is the set of all connected, undirected graphs with N nodes, and each node has one-hot features with length L, = 5 as its label. We propose these functions as benchmark problems, given: (i) there no existing synthetic benchmarks in graph BO literature, (ii) these benchmarks impose neither cobbeem-ppecifi constraints nor assumptions over the graph space (except for connectivity). making them suitable for comparison of a wide class of methods. These benchmark functions for graph BO are available at: [link to be added after peer review].
BoGrape is compared against the following baselines: (i) Random: random sampling, i.e., randomly sample one connected graph at each iteration, (ii) RW-rand: use graph GP with random walk (RW) kernel as surrogate, and sampling-
acquisition optimization, i.e., choosing the sample with the best LCB value among 20 random graphs, (iii) WI--and: use Weisfeiler-Lehhe RW-rand, (iv) WL-evvll: use evolutionary algorithm for acquisition optimization WL.-rand.
(WL.) kernel in
Remark 4.2. WL-rand and WL-evol are adapted from Ru et al. (2021), which is specifically designed for neural architecture search (NAS). WL-evol could be regarded as the state-of-the- method in graph BO.
For each benchmark with size N = 10, 20, we conduct BO with 10 initial samples and 50 iterations. When solving Eq. ( we observed good solutions to be found early (since more time is spent on proving optimality) and set 600s as the MIP time limit. As shown in Figure 3 BoGrape with all kernel variants outperforms baselines in most cases. When the graph size is small, SP and ESP perform better, since they are more expressive. For larger sizes, using simpler kernels reduces model
9

---

Published as a conference paper at ICLR 2026
Limeade
RW-lime
WL-llme
SSP
SP
ESSP
ESP
besty
Ber
Rer
Ecr
(a) QM7, N = 10
(b) QM7, N = 20
(c) QM7, N = 30
besty
Eor
(d) QM9, N = 10
(e) QM9, N = 20
(i) QM9, N 30
Figure 4: Bayesian optimization results on QM7 and QM9 with N E (10, 20, 30). Best objective value is plotted at each iteration. Mean with 0.5 standard deviation over 10 replications is reported.
complexity and produces better solutions within the given time limit. The trade-off is between the expressiveness of kernel and the complexity of the resulting optimization problem. We further discuss this computational limitation in Appendix B.5 and possible solutions in Appendix Bnd kernel selection in Appendix B.8
4.3 REAL-WORLD CASE STUDY
real- world case study. Following Zhang et al. (2023), we train two GNNs on dataset QM7 Blum & Reymond 2009 Rupp et al. 2012) and QM9 (Ruddigkeit et al. 2012 Ramakrishnan et al. 2014) as oracle predictors, I.e., the functions that we seek to optimize. The additional challenge of this task compared Section 4.2 is that molecules are not arbitrary labeled graphs: the molecular graph should be compatible with atom features. We found these structural constraints to effectively prevent the sampling- and evolutionary-based methods used in Section 4.2 from producing feasible solutions. Therefore, we modify Random to only consider randomly generated feasible molecules from Limeade (Zhang et al., (025). and we remove WL-evol from comparison. We also adapt molecular fcasibility constraints from Limeade and add them to our MIP formulation ensure only valid molecules are considered during the optimization. Further related details are given in Appendix
We next consider optimal molecular design (McDonald et al. 2024 Zhang al. 2024
The computational results in Figure show that BoGrape again generally outperforms baselines regardless of which kernel is used. We discuss this application further in Appendix B.7
CONCLUSION
This work proposes BoGrape to optimize black-box functions over graphs. Four shortest-path graph kernels are presented and tested on both prediction and Bayesian optimization tasks. The underlying mixed-integer formulation provides flexible and general platform including mixed-feature search spaces, graph kernels, acquisition functions, and problem-specific constraints. Our results show promising performance and suggest trade-offs between query-efficiency and computational time when choosing a suitable kernel. Future work may further simplify the formulations BoGrape and relax the requirement graph connectivity.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
The authors gratefully acknowledge support from Department of Computing Scholarship (YX), BASF SE, Ludwigshafen am Rhein (SZ), Engineering and Physical Sciences Research Council Igrant numbers EP/W003317/1 and EP/X025292/1] (RM, CT. JQ), Research England's Expanding Excellence in England (E3) fund awarded to MARS (Mathematics for Al in Real-world Systems) at Lancaster University (JQ), a BASF/RAEng Research Chair Data-Driven Optimisation (RM), BASF/RAEng Senior Research Fellowship (CT).
REPRODUCIBILITY STATEMENT
We take the following measures to facilitate the reproducibility of our work. Theoretical contributions are detailed and explained in both Section the main paper and Appendix Theoretical claims are supported by formal proof' provided in Appendix[A.4 To support the replication of empirical findings, we provide code implementations including our method, synthetic graph functions mentioned in Section 4.2 and models used in experiments in the supplementary materials.
REFERENCES
Bashar L. Ammari, Emma S. Johnson, Georgia Stinchfield, Tachun Kim, Michael Bynum, William E. Hart, Joshua Pulsipher, and Carl D. Laird. Linear model decision trees as surrogates optimization of engineering applications. Computers & Chemical Engineering. 178, 2023.
Ross Anderson, Joey Huchette, Will Ma, Christian Tjandraatmadja, and Juan Pablo Vielma. Strong mixed-integer programming formulations for trained neural networks. Mathematical Programming. 183(1):3-39, 2020.
Lorenz C. Blum and Jean-Louis Reymond. 970 million druglike small molecules for virtual screening in the chemical universe database GDB-13. Journal of the American Chemical Society, 131(25):
8732-8733 2009.
Karsten Borgwardt, Elisabetta Ghisu, Felipe Llinares-Lópezz Leslie O' O'Bray, Bastian Rieck, et al. Graph kernels: State-of-the-an and future challenges. Foundations and Trends0 Machine Learning, 13(5-6):511-712, 2020.
Karsten M Borgwardt and Hans-Peter Kriegel. Shortest-path kernels graphs. In International Conference on Data Mining, 2005.
Viacheslav Borovitskiy, Iskander Azangulov, Alexander Terenin, Peter Mostowsky, Marc Peter Deisenroth, and Nicolas Durrande. Matern Gaussian processes on graphs. International Conference Artificial Intelligence and Statistics, 2021.
Francesco Ceccon, Jordan Jalving, Joshua Haddad, Alexander Thebelt, Calvin Tsay, Carl D. Laird, and Ruth Misener. OMLT: Optimization & machine learning toolkit. Journal of Machine Learning Research, 23())15529-15836, 2022.
Fabrizio Costa and Kurt De Grave. Fast neighborhood subgraph pairwise distance kernel. In ICML, 2010.
Jiaxu Cui and Bo Yang. Graph Bayesian optimization: Algorithms, evaluations and applications. arXiv preprini arXiv: 1805.01 2018.
Thomas Elsken, Jan Hendrik Metzen, and Frank Hutter. Neural architecture search: of Machine Learning Research, 2019.
survey. Journal
Matthias Fey and Jan E. Lenssen. Fast graph representation learning with PyTorch Geometric. In ICLR 2019 Workshop on Representation Learning on Graphs and Manifolds, 2019.
Matteo Fischetti and Jason Jo. Deep neural networks and mixed integer linear optimization. Con- straints, 23(3):296-309, 2018.
Robert W Floyd. Algorithm 97: Shortest path. Communications the ACM, 5(6):345-345 1962.
11

---

Published as a conference paper at ICLR 2026
Peter Frazier. A tutorial Bayesian optimization. arXiv preprint arXiv: 2018.
Kenji Fukumizu. Kernel method: Data analysis with positive definite kernels. Graduate University of Advanced Studies, 2010.
Blake B Gaines, Chunjiang Zhu, and Jinbo Bi. Explaining graph neural networks with mixed-integer programming, Neurocomputing, pp. 130214, 2025.
Roman Garnett. Bayesian Optimization. Cambridge University Press, 2023.
Thomas Gärtner, Peter Flach, and Stefan Wrobel. On graph kernels: Hardness results and efficient alternatives. In Learning Theory and Kernel Machines, 2003.
Gurobi Optimization, LLC. Gurobi optimizer reference manual, 2024. URL https: /www gurobi com
Will Hamilton, Zhitao Ying, and Jure Leskovec. Inductive representation learning on large graphs. In NIPS, 2017.
José Miguel Bernández-Lobato Matthew W Hoffman, and Zoubin Ghahramani. Predictive entropy search for efficient global optimization of black-box functions. Neur/PS, 2014.
Christopher Hojny, Shiqiang Zhang, Juan S Campos, and Ruth Misener. Verifying message-passinee neural networks via topology-based bounds tightening. In ICML, 2024.
Donald R. Jones, Matthias Schonlau, and William J. Welch. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13:455 492, 1998.
Kirthevasan Kandasamy, Willie Neiswanger, Jeff Schneider, Barnabas Poezos, and Eric Xing. Neural architecture search with Bayesian optimisation and optimal transport. NeurlPS, 31, 2018.
Thomas N. Kipf and Max Welling. Sumi-super classification with graph convolutional networks. In ICLR, 2017.
Ksenia Korovina, Sailun Xu, Kirthevasan Kandasamy, Willie Neiswanger, Barnabas Poczos, Jeff Schneider, and Bric Xing. Chembo: Bayesian optimization of small organic molecules with synthesizable recommendations. In AISTATS, 2020.
Nils Kriege and Petra Mutzel. Subgraph matching kernels for attributed graphs. In ICML, 2012.
Nils M Kriege, Picrre-Louis Giscard, and Richard Wilson. On valid optimal assignment kernels and applications to graph classification. NewrIPS, 2016.
Nils M Kriege, Fredrik Johansson, and Christopher Morris. A survey on graph kernels. Applied Network Science, 5:1-42, 2020.
Harold J Kushner. A new method of locating the maximum point an arbitrary multipeak curve in the presence of noise. Journal of Basic Engineering, 86(1):97-106, 1964.
Huidong Liang, Xingchen Wan, and Xiaowen Dong. Bayesian optimization of functions over node subsets in graphs. NeurlPS, 2024.
Alexander G. de G. Matthews, Mark van der Wilk, Tom Nickson, Keisuke. Fujii, Alexis Boukouvalas, Pablo León- Villagil,,,,,, Zoubin Ghahramani, and James Hensman. GPflow: A Gaussian process library using TensorFlow. Journal of Machine Learning Research, 18(40): 2017.
Tom McDonald, Calvin Tsay, Artur M. Schweidtmann, and Neil Yorke-Smith. Mixed-integer optimisation of graph neural networks for computer-aided molecular design. Computers & Chemical Engineering, 185:108660, 2024.
Rocío Mercado, Tobias Rastemo, Edvard Lindelöf, Günter Klambauer, Ola Engkvist, Hongming Chen, and Esben Jannik Bjerrum. Graph networks for molecular design. Machine Learning: Science and Technology, 2021.
Velibor V. MiKić. Optimization of tree ensembles. Operations Research, 68(5):6605-1122 2020.
12

---

Published as a conference paper at ICLR 2026
Miten Mistry, Dimitrios Letsios, Gerhard Krennrich, Robert M. Lee, and Ruth Misener. Mixed- integer convex nonlinear optimization with muddient-boosted trees embedded. INFORMS Journal Computing, 33(3): 103-)111001119,,,, 2011.
Giannis Nikolentzos, Giannis Siglidis, and Michalis Vazirgiannis. Graph kernels: A survey. Journal of Artificial Infelligence Research, 72:943-1027, 2021.
Changyong Oh, Jakub Tomezak, Efstratios Gavves, and Max Welling. Combinatorial Bayesian optimization using the graph cartesian product. NeurlPS, 2019.
Theodore P Papalexopoulos, Christian Tjandraatmadjj, Ross Anderson, Juan Pablo Vielma, and David Belanger. Constrained discrete black-box optimization using mixed-integer programming. In International Conference on Machine Learning, pp. 17295-17322. PMLR, 2022.
Joel Paulson and Calvin Tsay. Bayesian optimization as a flexible and efficient design framework for sustainable process systems. Current Opinion in Green and Sustainable Chemistry, pp. 100983, 2024.
Dhanesh Ramachandram, Michal Lisicki, Timothy J Shields, Mohamed R Amer, and Graham W Taylor. Structure optimization for deep multimodal fusion networks using graph-induced kernels. arXiv preprint arXiv: 1.XXiv117700070000 2017.
Raghunathan Ramakrishnan, Pavlo O. Dral, Matthias Rupp, and Anatole Von Lilienfeld. Quantum chemistry structures and properties of 134 kilo molecules. Scientific Data, 1(1):1-7, 2014.
Binxin Ru, Xingchen Wan, Xiaowen Dong, and Michael Osborn.. Interpretable neural architecture search via Bayesian optimisation with Weisteiler-Lehman kernels. In ICLR, 2021.
Lars Ruddigkeit, Ruud Van Deursen, Lorenz C. Blum, and Jean-Louis Reymond. Enumeration of 166 billion organic small molecules in the chemical universe database gdb- 17. Journal of Chemical Information and Modeling, 52(1 ((111)22644222 2875, 2012.
Matthias Rupp. Alexandre Tkatchenko, Klaus-Robert Müller, and O. Anatole Von Lilienfeld. Fast and accurate modeling of molecular atomization energies with machine learning. Physical Review Letters, 10865)0058301, 2012.
Eric Schulz, Maarten Speekenbrink, and Andreas Krause. A tutorial on Gaussian process regression: Modelling, exploring, and exploiting functions. Journal of mathematical psychology, 85, 2018.
Artur M Schweidtmann, Dominik Bongartz, Daniel Grothe, Tim Kerkenhoff, Xiaopeng Lin, Jaromit Najman, and Alexander Mitsos. Deterministic global optimization with Gaussian processes embedded. Mathematical Programming Computation, 13(3):553-588 5811 2021.
Artur M Schweidtmann, Dominik Bongartz, and Alexander Mitsos. Optimization with trained machine learning models embedded. In Encyclopedia of Optimization, pp. 1-8. Springer, 2022.
Nino Shervashidze, SVN Vishwanathan, Tobias Petri, Kurt Mehlhorn, and Karsten Borgwardt. Efficient graphlet kernels for large graph comparison. In AISTATS, 2009.
Nino Shervashidze, Pascal Schweitzer, Erik Jan Van Leeuwen, Kurt Mehlhorn, and Karsten M Borgwardt. Weisfeiler-lehman graph kernels. Journal of Muchine Learning Research, 12(9), 2011.
Giannis Siglidis, Giannis Nikolentzos, Stratis Limnios, Christos Giatsidis, Konstantinos Skianis, and Michalis Vazirgiannis. GraKel: A graph kernel library in Python. Journal of Machine Learning Research, 21(54):155, 2020.
Niranjan Srinivas, Andreas Krause, Sham Kakade, and Matthias Seeger. Gaussian process optimiza- tion in the bandit setting: No regret and experimental design. ICML, 2010.
Niranjan Srinivas. Andreas Krause, Sham M. Kakade, and Matthias W. Seeger. informition-theoretiic regret bounds for Gaussian process optimization in the bandit setting. JEEE Transactions on Information Theory, 58:3250-3265. 2012.
13

---

Published as a conference paper at ICLR 2026
Alexander Thebelt, Jan Kronqvist, Miten Mistry, Robert M. Lee, Nathan Sudermann- and Ruth Misener. Entmoot: A framework for optimization over ensemble tree models. Computers & Chemical Engineering, 151:107343, 2021.
Alexander Thebelt, Calvin Tsay, Robert M. Lee, Nathan Sudermann-Merx, David Walz, Tom Tranter, and Ruth Misener. Multi-objective constrained optimization for energy applications via tree ensembles. Applied Energy, 306: 18061, 2022a.
Alexander Thebelt, Johannes Wiebe, Jan Kronqvist, Calvin Tsay, and Ruth Misener. Maximizing information from chemical engineering data sets: Applications to machine learning. Chemical Engineering Science, 252:1 17469, 2022b.
Calvin Tsay, Jan Kronqvist, Alexander Thebelt, and Ruth Misener. Partition-based formulations for mixed-integer optimization of trained ReLU neural networks. In NeurlPS, 2021.
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph attention networks. In ICLR, 2017.
Stefan Vigerske and Ambros Gleixner. SCIP: Global optimization of mixed-integer nonlinear programs in a branch-and-cut framework. Optimization Methods and Software, 33(3):563-553 2018.
S Vichy N Vishwanathan, Nicol N Schraudolph, Risi Kondor, and Karsten M Borgwardt. Graph kernels. Journal of Machine Learning Research, 11:1201-1242 2010.
Xingchen Wan, Henry Kenlay, Binxin Ru, Arno Blaas, Michael Osbore,, and Xiaowen Dong. Adversarial attacks on graph classifiers via Bayesian optimisation. In NeurlPS, 2021.
Xingchen Wan, Binxin Ru, Pedro M Esperança, and Zhenguo Li. On redundancy and diversity in cell-based neural architecture search. In ICLR, 2022.
Xingchen Wan, Pierre Osselin, Henry Kenlay, Binxin Ru, Michael Oshorne, and Xiaowen Dong. Bayesian optimisation of functions on graphs. NeurIPS, 2023.
Keliang Wang, Leonardo Lozano, Carlos Cardonha, and David Bergman. Optimizing over an ensemble of trained neural networks. INFORMS Journal on Computing, 2023.
Colin White, Willie Neiswanger, and Yash Savani. Bananas: Bayesian optimization with neural architectures for neural architecture search. In AAAI, 2021.
Colin White, Mahmoud Safari, Rhea Sukthanker, Binxin Ru, Thomas Elsken, Arber Zela, Debadeepta Dey, and Frank Hutter. Neural architecture search: insights from 1000 papers. arXiv preprint arXiv:2301.0777 2023.
Yilin Xie, Shiqiang Zhang, Joel Paulson, and Calvin Tsay. Global optimization of Gaussian process acquisition functions using piccewise-linee kernel approximation. arXiv preprint arXiv:2410.16893, 2024.
Nianzu Yang, Huaijin Wu, Kaipeng Zeng, Yang Li, Siyuan Bao, and Junchi Yan. Molecule generation for drug design: a graph learning perspective. Fundamental Research, 2024.
Shiqiang Zhang, Juan S. Campos, Christian Feldmann, David Walz, Frederik Sandfort, Miriam Mathea, Calvin Tsay, and Ruth Misener. Optimizing over trained GNNs via symmetry breaking. In NeurIPS, 2023.
Shiqiang Zhang, Juan S Campos, Christian Feldmann, Frederik Sandfort, Miriam Mathea, and Ruth Misener. Augmenting optimization- molecular design with graph neural networks. Computers & Chemical Engineering, 186:108684, 2024.
Shiqiang Zhang, Christian W Feldmann, Frederik Sandfort, Miriam Mathea, Juan S Campos, and Ruth Misener. Limeade: Let integer molecular encoding aid. Computers Chemical Engineering, pp. 1091 15, 2025.
Yin-Cong Zhi, Yin Cheng Ng, and Xiaowen Dong. Gaussian processes graphs via spectral kernel learning. IEEE Transactions on Signal and Information Processing over Networks, 2023.
14

---

Published as a conference paper at ICLR 2026
A ENCODING OF GRAPH KERNELS
A.1 NOTATIONS
We provide details for all variables introduced in this paper in Table Recall that the search domain considered here consists of all connected graphs with node number ranging from ng to n, each node has M binary features with the first L node features as the one-hot encoding of node label.
Table 2: All variables introduced in the optimization formulation for graph kemels..
Variables
Domain
Description
Au,n, u, E [n]
[0, 1}
the existence of edge from u to v the length of shortest path from w appears the shortest path from u to
du,s u, V. E [n]
[n 1]
0,,0 u. U. w E [n]
[0, 1]
die Liwn U, E [n], E [n+1]
{0, 1 [n²+1]
indicator: [dd..
DR, s E [n]
shortest paths with length 8
DS. S [n], [re* + 1]
[0.1]
indicator: 1(D, = c)
Pu.o da u.u, E I1 E [4]
[0.1]
indicator: 1(Fa,,, = L, du,. 8, Fu,1a = 1)
Padidar SE 11), E [L]
[n2 [n]+1]
# shortest paths with length and labels 4 la indicator: 1 (P = c)
da E 1,, E [4], E [n² +
{0,1}
N WI [M] No. m E [M], C E [M + 1]
[N +1] [0.
sum of m-thh feature over all nodes indicator: 1NN,.. e)
A.2 SHORTEST PATH ENCODING FOR GRAPHS WITH FIXED SIZE
We first present necessary conditions that Au, du 8t should satisfy in Eq. (5)
Au, = 1,
Ve E [<]
(5a)
du 0,
Ve E [n]
(5b)
1,
Au, =
du,
Vu, E [m]. u +
(5c)
> 1, Au, = 0
dew dww 00, 1
du,
Vu,u E [m], u +
(5d)
duw + dw,
du. = 0'
1,
8
Ve E [n]
(5e)
0,
w + U
8"
=1, 1,
Vu, [n], u +
(5l)
2,
An, =
8
Vu, E |m], +
(5g)
2,
Au,u =
Eq. 5
restricts An,, d, o, in the following rules:
Eq. 5a)
initializes the diagonal elements.
Eq. 5b
initializes the shortest distance from itself.
Eq. 5ch forces the shortest distance from node u and e be if edge u - exists, and larger
than 1 otherwise.
Rewrite Eq. (5c) as:
dy S 1 + n (1 Ana), Vu, v E n), u + v
dy 2 - Au,or
Vu,u E [n], +
where is a big-M coefficient using du, S n 1.
15

---

Published as a conference paper at ICLR 2026
Eq. 5d is the triangle inequality for distance matrix d. Rewrite Eq. (5d) as:
du, I du,o + dw. - (1 - 0 pm du, 2 du,o + dw. 2n (1 8u,
Vu, v, w E [n]
v,, v, E [n]
where 2n is big-M coefficient since du,w + du.o < 2n.
Eq.
Seh
initializes 60 by definition.
Eq.
5i
initializes &h and 6i by definition.
ensures that there is no edge from node u to Otherwise, no node except for u and could appear at the shortest path from u to v.
Eq. 5
least one node at the shortest path from node u if there
Rewrite Eq. g5 as:
Σ WE(n]
8th 2 (n 2). An,u); Vu, [n]. u
Σ
ou 2 2 2 (1 - -A,,)),
Vu [n]. #
WEIn]
where n 2 is big-M coefficient since Ewe(n) 8 S n.
Replacing disjunctive constraints accordingly in Eq. 5 gives the final formulation Eq. MMPP-SP
A.3 SHORTEST PATH ENCODING FOR GRAPH WITH UNKNOWN SIZE
We extend constraints Eq. (5) to handle changeable graph size. Full constraints are as follows:
A Au41,041
Va [n 1]
(6a)
A.p 2 n0,
(6b)
DE[n]
2A S An, + An,os
Vu, E [n], u +
(6c)
dg 0,
E [n]
(6d)
1, Au,
du
Vu, v [n], t
(6c)
> 1, An, = 0)
n1 Au,u Au,
dy
Vu, [n], u
(6f)
n, min An,u: Av,
= du + dm
8w = Ou, 0
du.
Vu, [n], u u
(6g)
duiw + dw
1, w = v
8"
Ve E [n]
(6h)
0,
w +
8M
8
=
Vu, [n],
(6i)
2,
An, =
Σ
> 2,
Au, = 0, Au, = Ag, = 1,
Yu, E [n], u +v
(6)
a.v
W
2,
mint Au,us A 0
Eq. 6
restricts An, d, 8 in the following rules:
Eq. 6a)
forces nodes with smaller indexes exist.
Eq. 6b
gives the lower bound of the number of existed nodes.
Eq. (c
means that there is no edge from node u to any of them does not exist.
Eq.
6d
initializes the shortest distance from one node itself, it does not exist.
16

---

Published as a conference paper at ICLR 2026
Eq. (6e) forces the shortest distance from node and be 1 if there is edge from and larger that 1 otherwise.
Rewrite Eq. (6e) as:
dup S + n (1 Au,), Vu,u E [n], u + Vu, E [n], +
du,v 2 Au,o
where a big-M coefficient using the fact that du,v I
Eq. (6f) sets the shortest distance from node u to v as n, i.e., oc, if any of them does not
exist. Otherwise, the shortest distance is less than n.
Rewrite Eq. (6) as:
duv 2 (1 Au,,) Vu, E [n]. u u
du 2 n (1 Aun), Vu, E [n], u +
Eq. 66g is the triangle inequality for the distance matrix d. Rewrite Eq.
66 as:
du, S du, dw.w (1 - o0,
Vu, u, E [n]
du,v du,w dw. 2n (1 0,,
Vu, E [n]
where 2n is big-M coefficient since du + dur S 2n.
Eq.
6h
initializes d..,. by definition, even node v does not exist.
Eq. 6i
initializes 84 and 6., by definition, even node u or w does not exist.
Eq. [jj makes sure that there is at least on node at the shortest path from node u to v if there is no edge from node and and these two nodes both exist. Otherwise, only 8ú, and 8 equal to 1.
Rewrite Eq. 66 as:
Σ
u, 2 (n 2) - -A,,),
Vu, E [n], +
WE
Σ
< 2 (n -2) :). Au,u
Vu, [m], +
Σ
u,v < 2 + (n -2) A
Vu, E [n], #
Σ
8m > An, + A + (1 An,), Vu, [m], u +
#
m)
where 2 is big-M coefficient since Emejn) 8 < n.
To conclude, the formulation for shortest paths of all connected graphs with least nn nodes and at most nodes is:
17

---

Published as a conference paper at ICLR 2026
Av,o 2 Ast1.+1
Ve E [n 1]
Av,o 2 no:
MEin]
2 I Au,u Av,
Vu, u E [n], +
d2 = 0,
Ve E [n]
du 5 + (1 Au,)):
Vu, U E [n], u +
dup 2 2 - -An,
Vu, V E [n], u # v
du 2 n (1 Au,u)
Vu, E [n], u u
du 2 n (. Ann)
Vu, v E [n], u u
du < dm m0 + duni (1 ou,r);
Vu, U, [n]
du 2 du + duu.
2n (1 0g,
Vu, u, w [n]
(MIP-SP-plus)
60,
11
W = v
Ve E [n]
10
+ +v²
0000 8 1,
Vu, E [n], +
Σ
u.v 2 + (n - 2) (1 Au,)),
V [n], u #
Gup 2 (n (n-2) 2) Augu 2)).
Vu,ue E [n], u #
Ou 2 (n 22)- An,os
Vu,v E [n], +
[
Σ WE(n]
t 2 Au, + Au, + (1 Au,o),
VU E [n], u u
A.4 PROOFS THEOREMS
Proof of Theorem 3.5 If such G exists, it is unique since Au, gives the existence of every edge. Thus it suffices to show that (G),, 88 Sw, (G)) = (du,e, 8 for G defined with A,
We are going to prove by induction on the shortest distance sd from node u to in graph G. Specifically, we want to show that for any 0 < sd < and for any pair of (u,v) such that min(du, (G), du,v) = sd, we have (G, = du,v and (G, = Sw Vu E [n].
For sd = O, (G) du,)) = and only if u = v. For any e E [m], it obvious to have:
de (G) = 0 = d
8 (G) = = 8"
8 (G) = 0 = 8th
Vu +
For sd consider every pair (u, v) such that (G) 1, we have Au,o (G. 1, then is easy to obtain:
du,, (G) 1 du,"
8w, (G) = 1 = sw Vi E {u, t}
8w, (G' = 0 = sw Vw & {u,
where 8"r 0, Ver @ {u, v} since:
Σ [m]
50 0. u,e
a,e
a,
mg (u,))
On the contrary, du,v gives Au, = 1, thus An, (G) 1 and 8w, (G) 5w ou.
Vi by definition.
18

---

Published as a conference paper at ICLR 2026
Now assume that for any pair of (u, v) such that min(d (G), du,u) I sd, we have (G. du, and (G) = Su,e: V.. Since 8w,,, (G) = 6u,ps Vw E {u, 2) always holds by definition, we only consider w & {u.u}.
Part 1: We first consider every pair of (u, v) such that (G) sd + 1. Since sd + 1 > 2, we know that Au, = (G, = and there exists w & {u, 2 on the shortest path from node u to v in graph G.
Case 1.1: For every w & {u, 2t such that (G) = 1, since dG, I sd and ( & sd, we have:
du S dy + dw = d (G) + dw, = ( = sd + The equality has to hold, otherwise, du,u S sd gives du duu, sd by assumption. Therefore, o = 1 = 8 (G).
Case 1.2: For every w & (u, u} such that 83,,, (G) = 0, if o. = 1, then du, + dm, = du, = sd sd+11 1, which means that du S sd and duu,o S sd. By assumption, we have du, (G) = du du (G) d and then:
du, + duu. = du, + duu, = due = dG,
which contradicts to W. (G) O. Thus ow, O.
Part 2: Then consider every pair of (u, such that du, sd + 1. Similarly, we have Au, = (G)
since du,w sd and du,o sd, have Case 2.1: For every & {u, such that 8
d (G) d and (G) dew, then:
d S du, (G) + d = du,w + dw, du,e = sd +
This equality also has to hold, otherwise, Bu,((G) S sd, by assumption du,o = d I sd, which a contradiction.
Case 2.2: For every w & {u, e) such that ou, = O, if 8 (G) =1, then (G) dew (G) du,v (G) = sd + 1, which means that du (G) S sd and dhe (G)) I sd. Therefore, du,w dw (G) + dww (G) du,v(G) du,us
which contradicts to ōu = 0.
Proof of Theorem 3.6 Fix the node number as n1 with nu S n2 S n. Eqs.
6a (6b)
force:
A = {
1,
v E [m]]
o. e [n]][a]]
substituting which to other constraints give us:
Au, Av = 0,
Vu E [M1]. eE [n]]] u u
du, 0,
Ve e [n]]m]]
du du,u =
Vu e [m1]. e [n]]]]] u
[1, {u,
8 Di
Vu [m1], v [m]]m]]
10,
w & (u,
One can easily check that all constraints associated with non-existed nodes are satisfied. Removing those constraints turns Eq. MMP-SP-pllus into Eq. (MIP-SP) with size ni.
A.5 ENCODING FOR KERNEL OVER BINARY FEATURES
Assume that each graph G has binary feature matrix F E (0, (MM)))M we need to formulate (F, F*) and (F(F,P F) properly. kr could be defined in multiple ways, here we propose permutationallinnarrt kernel considering the pair-wise similarity among node features. Given two feature matrices F1, 222 corresponding to graphs G¹, G² resp., define kF as:
kr (F(FFF2))
F2)
1
Σ
FL
F2
Σ Nm N (F²),
M
B1n2 M mE(M]]
19

---

Published as a conference paper at ICLR 2026
where Nm (F')
Σ
Femm Vm
[M], and ningM is the normalized coefficient.
re/n.]
Similar to Section
3.4
we have:
(F F*)
Σ
(Fn
Nm
nn,M m[[M]]
where Nn Σ
Fu
Vm [M]], and:
vEjn
kr(F,F)
1
Σ
N2 1
E
2 N
112
n²M ME(M-]
mE[M]
where indicators No. 1(Nm c). Vm [M], c E [n 1] satisfy:
Σ No 1, Σ
N Nm, Vm E [M].
cE[n+1]
eE[n+1]
A.6 SIMPLIFY PATH ENCODING OVER UNDIRECTED GRAPHS
For undirected graphs, we first add the following constraints to guarantee symmetry:
Au,e Au,u
Vu, v [n], Yu, E [n], <<u Vu, U, [n]. u u<v
du, = de
6
Since the inverse of any shortest path from node u to also shortest path from node to u, for SSP and ESSP kernels, D,, Vs [n] are even and we can fix odd indicators as zero:
D 1, is even
Vs E [n]. o [n² 1].
0, is odd
Similarly, for SP and ESP kernels, we have:
Ps Padadh Vs E [n], fi, [Z].
B DISCUSSION
B. CHOICE GRAPHI KERNELS
Various graph kernels are proposed to better fit graph data. However, none of them could be incorporated as optimization constraints (nor are they designed for this purpose). Thus, current graph BO works mostly use evolutionary algorithms that generate candidates and then evaluate them, where graph kernels are used as graph functions to calculate the posterior mean and variance. The major difference here is that computing k(G', G²) given both G and G2 is quite easy, but representing A(G G2) only given G super challenging since G2 could be any arbitrary graph. BoGrape is built upon theoretical contributions connected graphs. Therefore, it is not that we chose shortest-par -artest-paat&& kernel first for specific reasons then developed necessary formulations, but that the fundamental advances graph optimization led shortest-path kernels.
encoding shortest paths into decision variables for arbitrary
B.2 CHOICE OF ACQUISITION FUNCTIONS
BoGrape formulates acquisition optimization as a MIP, and LCB is chosen as a representative acquisition function given its popularity in BO and relatively simple form. BoGrape could be easily applied to other acquisitions functions in linear forms w.r.t. posterior mean and variance. For nonlinear acquisition functions, one could either nonlinear solvers to optimize the resulting MINLP or linearize the acquisition functions. Since the acquisition only appears in objective, all graph-relevant constraints still work as before.
20

---

Published as a conference paper at ICLR 2026
B.3 EFFECTIVENESS ENCODING
Different shortest-path algorithms might affect the time complexity when computing the graph kernels, but they should be asymptotically similar since cost dominated by computation of the shortest distance between any pair of nodes. For example, if we use Dijkstra's s algorithm, which a single-source shortest path algorithm, then we need to repeat it n times and the complexity is O(n(e + n nt) with the number of edges. Most importantly, the choice of shortest path algorithm is irrelevant to our shortest path encoding. Although our encoding is motivated by Floyd's algorithm, all constraints in our encoding are the necessary conditions that the shortest paths should satisfy no matter what algorithm is used. For the optimality of the encoding, our encoding builds bijection between all connected graphs and all feasible solutions of Eq. (MIP-SP) as shown in Theorem 3.5 meaning it is optimal in terms of representations.
B.4 COMPLEXITY ANALYSIS
The time complexity of computing all shortest paths for graph with nodes is O(n²). When computing the covariance between two graphs (assume they both have n nodes for simplicity), naive implementation of shortest path kernel is Onn*), but implementation is O with L being the number of node labels after storing P Is (-) as defined Section 3.4 For other graph kernels, Random Walk (RW) (Gärtner et 2003) is O(n²), Subgraph Matching (Kriege Mutzel 2012) ((kn*+1) with k as the size of subgraphs considered, and Weisfeiler-Lehmar (WL) (Shervashidze et al. 2011 is (hum) with / as the number of iterations and m as the number of edges [7]. There graph kernels with lower complexity than shortest-par partesst-pathhaee kernels, but the complexity of calculating kernels in graph BO is less important than encoding the graph kernel as optimization constraints.
B.5 LIMITATIONS
The major limitation of BoGrape (and most MIP-based methods) is the computational complexity. The BoGrape complexity stems from solving the MIP rather than computing kernel values. The tradeoff is: (i) BoGrape represents the whole, unavoidably large, search space precisely, and (ii) solving MIP to global optimality is cons-cons since proving optimality of a solution usually takes much more time than finding this solution. To better demonstrate this tradeoff, we perform an ablation study by varying the MIP time limit among {60, 600, 1200] seconds on molecular design case study on QM7 dataset with graph size N 10. As Figure 5 illustrates, extending the computational time does not improve BO performance significantly. Nevertheless, Figure shows that increasing time limits results in smaller MIPgap, i.e. gives the the solver to more time prove a solution's 's optimality. In other words, finding good feasible solutions is easier (and important for practical BO performance), while closing the MIPgap (important for theoretical BO convergence) requires more computational effort.
B.6 SCALIBILITY OF BOGRAPE
Scalability issues always exist for combinatorial optimization, since the search space grows quickly, For BoGrape, there are several ways to improve scalability: (i) reduce the search space, e.g., only consider graphs that are similar previous graphs in trust-region fashion (similar in spirit muta- tion over existing samples in evolutionary algorithms (Ru et al.) (202)). or adversarial perturbations with limited budgets (Wan al., (021)) (ii) limit the solving time as we did in experiments, letting the MIP solver return the current best solution, (iii) develop computational heuristics for specific problems to identify promising candidates earlier, (iv) decompose large graphs into functional groups and optimize the graph structure over groups instead of nodes. Note that (iv) is frequently applied in graph tasks, e.g., cell-based neural architecture search (Wan et al. 2022 fragment- based molecular design (Zhang et al. 2024), etc..
B.7 CHOICE OF APPLICATION
We choose the optimal molecular design task since (i) molecules can be represented attributed, connected graphs, (ii) molecular properties, either measured or predicted, are suitable functions over graphs, and (iii) the biP-based framework for molecular design is well-established. The baselines
21

---

Published as a conference paper at ICLR 2026
600
600
0.3
1200
0.3
1200
best
-0.4
-04
-0.5
0.5
10
20
30
40
50
0
10
20
30
40
50
Iter
tter
(a) SSP kernel
(b) SP kemel
-0.25
00.5
60
60
600
600
1200
1200
.0.35
0.35
best y
besty
-0.45
0.45
0.5
40
50
Iter
tter
ESSP kemel
(d) ESP kernel
Figure 5: Performance of varying the time limit for BoGrape over QM7 datasets with graph size N 10. Best objective plotted at each iteration. Mean with 0.5 standard deviation over 10 replications is reported.
SSP
8
SP
Average MIPgap
6
4
60
1200
time limit (s)
Figure 6: Comparison of the average MIPgap over all iterations when varying the time limit. Experi- ment conducted on QM7 dataset with graph size N = 10.
22

---

Published as a conference paper at ICLR 2026
used in experiments are less tailored to molecular design, and there definitely more advanced methods. But the purpose of this case study is not showing BoGrape is a state-of-the art method in molecular design, but investigating this problem from constrained discrete optimization perspective. Meanwhile, although molecular design is a promising and important application area for BO (Paulson) & Tsay 2024), our proposed BoGrape procedure is general for any setting with functions defined over connected graphs.
B.8 KERNEL SELECTION
Kernel selection is an interesting question explored in BoGrape As discussed in Section there is a trade-off between the kernel's 's expressiveness and the complexity of resulting optimization problems. With sufficient computational resources. more expressive kernels like ESP is preferred. But simpler kernels like SSP yield optimization problems that are easier to solve. As the graph size increases, linear kernels usually achieve better performance due the overhead associated with formulating exponential kernels.
C ADDITIONAL NUMERICAL RESULTS
C.1 KERNEL PERFORMANCE
Besides four shortest-pan paatesst-pathhhe graph kernels, we also test the performance of several classic graph kernels, including Random Walk (RW) (Gärtner et al. 2003), Subgraph Matching (SM) (Kriege & Mutzel 2012), Weisfeiler-Lehhan (WL) (Shervashidze et al., 2011), and Weisfeiler-Lehman Optimal Assignment (WLOA) (Kriege al. 2016) kernels. To justify the effectiveness of the feature component in Eq. 33. we also test the combination of these four kernels with the same feature kernel used in shortest-path kernels. All GPs are trained by maximizing the log marginal likelihood. During GP training, we set bounds for kernel parameters, i.e., o, 3, of to [0.01, 100] with 1 as their initial values, and set noise variance o? as 10-0 3./²2 defined in Eq. lab is set as 1.
Datasets QM7 (Blum & Reymond 2009 Rupp et al. 2012) and QM9 (Ruddigkeit et al., 2012 Ramakrishnan al. 2014g are used test the kernel performance and train GNNs as graph functions. Each dataset consists of molecules with quantum mechanic properties. Following the setting in Zhang al. (2024), we represent each molecule as graph with included, use the same structural constraints, and train GNN predictor for each dataset. The trained GNN on QM7 has train and test errors 0.0356 and 0.0337 respectively. Both the train and test errors of the GNN QM9 are 0.0082. We provide example of the node feature and label in such molecular graph to better distinguish the difference in their definitions:
node features with I, = labels
Example. In the molecular design example QM7 dataset, followed the same setting as in Zhang et al. 2024). Each node (atom) has one label from {C, N, S,S) and a feature vector with length M 15, e.g. (1, 0, 0, 0,0,0, 1, O, 0, 0, 1, 0, 0, 0) where the first four elements indicate the atom has label C, the 5th 8th elements indicate the atom has two neighbors. the 9th 13th elements indicate the atom connected 2 hydrogen atoms, the 14th element indicates the atom included in double bond and the t5th element indicates it is not included in a triple bond. More detailed definitions can be found the Table 2 & of Zhang et al. (2024).
Based on the molecular size N, we consider two settings (a) if the dataset includes molecules of size N, we randomly choose molecules from the dataset and use their real properties, and (b) for larger N, we use Limeade to generate molecules and use the trained GNN to predict their properties. To show the performance of different kernels on representing similarity between graphs, we apply setting (a) and perform property prediction task using GPs equipped with the various kernels, shown in Figure 2 For larger graph sizes, we apply setting (b). The root mean square errors (RMSE) is reported in Figure & and Table 3 the mean negative log likelihood (MNLL) is reported in Table 4
Two observations from these results are: (i) adding feature component largely improves the per- formance of all kernels in terms of predictive accuracy and uncertainty quantification, and (ii) our shortest-pant -ath kernels achieve comparable performance comparing to other graph kennels, which further justifies our choice of these kernels for global acquisition optimization.
23

---

Published as a conference paper at ICLR 2026
C.2 ABLATION STUDIES ON THE CHOICE OF 3t
BoGrape leverages the classical LCB acquisition function, where exploration and exploitation are balanced through its coefficient 3. (Srinivas et al. 2010). We set 3 1 in our experiments for simplicity, as using constant values for 3e is a standard approach in BO literature, e.g., (Thebelt et al.) (2021) Although changing 3, does not affect the complexity of the acquisition optimization, we provide an ablation study on varying the value of 3 as an investigation on the convergence behavior of BoGrape regarding different exploration cntprratioonnndd--ee and-exploitation factor. We consider the same setup of the real-world case study on the QM7 dataset using SSP kernel with graph size N = 10 as in Section
4.3 We include three common choices of 8./2 in BO literature: (1) as this work; (2) 1.96, e.g. Thebelt et (202)); (3) a time-dependent schedule of 3 00.5 + 1)) as suggested by recent graph BO literature (Ru et al. 2021). We present the average minimal objective value (with 0.5 standard deviation brackets) found over iterations with 10 replications in Figure This study confirms that, though there may indeed be value in tuming B, for a particular setting, BoGrape variants with different choices for 3e largely exhibit similar convergence behavior. This observation justifies the choice on Br and further proves the robustness of BoGrape over key hyperparameters.
C.3 DETAILS FOR CASE STUDY
Random sampling is common baseline, but is excluded in Section since it rarely produces even feasible solutions. Randomly sample feasible graphs is untrivial in molecular generation task because the graph structure and features should be reasonable and compatible with each other, e.g.. satisfying structural feasibility, dataset-speciffc constraints, etc.. Here we consider random sampling over QM7 and QM9, to guarantee the feasibility of samples and compare it with Limeade (Zhang ctt al., 2025). Limeade is proposed as a feasible molecule generator, which is further enhanced by incorporating the composition constraints and symmetry -breakknn constraints (Zhang et al. 2023 Figure plots the regret curve over 50 iterations for both sample methods. In all cases, Limeade outperforms random sampling, showing the limitations of random sampling. Therefore, we choose Limeade as our sampling baseline.
For evolutioaary algorithm, we apply the same random sampling and mutation procedure as in Ru et al. 2 First, none of 10 random samples is feasible, which is expected since the sample domain is 2 N (~> 10 for = 10, ~ 10 for N = 20), while the feasible domain is relatively much smaller. Then we give evolutionary algorithms 10 feasible molecules generated by Limeade as an initial population and mutate each one 100 times, but only 0.35% mutations are feasible for N 10, and out of 106 are feasible for N 20. Although tailored evolutionary algorithm could be designed with better performance, it is neither the main focus of this work nor compatible with general settings. Therefore, we exclude WL-evol from baselines.
Given the trained GNNs used in Section C.I as unknown graph functions, we conduct BoGrape to optimize them. In our experiments, 10 random molecules sampled by Limeade are used as the initial dataset, and 50 BO iterations are performed. We set PoolSearchModes feasible solutions using Limeade. For each BO run, we show the mean with 0.5 standard deviation of the best objective value over 10 replications. For the two baselines where we use Limeade as sampling-bast solver for the acquisition functions, we conduct an ablation study by increasing the number of candidates from 20 to 100. Figure shows the results. While multiplying the number of candidates evaluated by five in each acquisition optimization step indeed improves the performance of the sampling- based baselines, BoGrape (which only proposes one sample per iteration) still outperforms Limeade. This ablation study emphasizes the importance of global acquisition function optimization. Besides the results shown in Figure 4 experimental results for N E {15, are reported in Figure II supporting our analysis in Section 4.3
in Gurobi to generate
24

---

Published as a conference paper at ICLR 2026
SSP
ESSP
ESP
RW-lime
RW-llme
WL-line
Wi-lime
2.6
2.6
2.8
2.8
3.0
best
3.4
3.4
-3.6
-3.6
10
10
20
50
Iter
Iter
(a) QM7, N 30
(b) QM7, N 30
-0.4
-0.4
best
-0.6
-0.6
-0.8
-08
-1.0
-1.0
0
10
20
30
40
50
10
20
30
40
50
iter
iter
(c) QM9, = 30
(d) QM9, N = 30
Figure Performance of varying the number of samples used acquisition optimization for baselines over QM7 and QM9 datasets with graph size N = 30. indicates 100 candidates used in each iterations. Best objective is plotted each iteration. Mean with 0.5 standard deviation over 10 replications is reported.
25

---

Published as a conference paper at ICLR 2026
Voln piot
1.2
1560 alss t9b0 1lwab0
Bolii (4-50))
1.35
04
6.10
Figure 8: Violin plots demonstrate GP regression performance on different graph kernels and graph sizes N. * indicates linear combination of given kernel and feature kernel. Each violin plots 25% percentile, median and 75% percentile of the RMSEs over the 100 replications.
26

---

Published as a conference paper at ICLR 2026
Table 3: Model performance of GPs equipped with different graph kernels. indicates linear combination of given kernel and feature kemel. For each graph size N, we use Limeade random generate 20 and 100 molecules for training and testing, respectivelyy root mean square error (RMSE) of predictive error is reported over 100 replications.
1.81((3911))
1.389 (2758)))
1.492(2.750)
1.122(3 005
(62((4907)
(157(19980) 1.574(4.516)
0.995 2.655)
1.278 (3777)
1.006 2.64)
1.049 2..4))
2.816(6.667
N=30
1.073
1.255 .499
1.678(3.041)
0.934(2.085 (2.05)))))
0.682(0.508) 1.00((2.055)
0.638(1 701)
(440())602) .602)
(.766(1.61)
0.523(1 1.523(1.723 723)
0.593(1 01593(1.887) 487)
0.278(0.286) (1.62()
0.608 (.608(1.702) 702)
N=25
0.540(1
1.456() (.56(10074) .074) 0.819(0.79)) (0.792
1.398(4.284)
0.682(0.501)
0.682(0.683)
0.354
0.436(0.550)
0.341(0.304)
0.520(0.613)
0.290(0.301)
0.289(0:293)
N=20
1.693(1 .22)
1.475( .057)
(612(1.084)
1.810(2.812 (00(2.819)
1.785(0.640)
0.928 (3.732)
(790()
0. 0.984(0. 920)
0.765 699)
(0.731)
0.866 (066(1100) .002)
N=15
1.55((09655)
1.301 945)
1.185 894)
1.256(0.930))
SSP 0.710(0.68) 680)
0.709 (1709(1.150) 150)
0.619(0.772)
0.539
0.826(0.922 (2.92)
0.519(0.670)
0.498(0.636)
0.545(0.710)
N=10
Kernel
RW
SM
WL
WLOA
SP
ESSP
ESP
RW*
SM
WL*
WLOA*
27

---

Published as a conference paper at ICLR 2026
Table 4: Model performance of GPs equipped with different graph kerelss. indicates linear combination of given kernel and feature kernel. For each graph size N. we use Limeade random generate 20 and 100 molecules for training and testing, respectively, mean negative log likelihood (MNLL) reported over I00 replications.
36426.89 19126.89111288155 (120815.092) 5.092)
19793. 19793.178 178.17886633... 8..82)2)))
19794.5 528 69333.022 3.027)
12862 222.2266182 84932.3333999
227.393(2336.9)) 968)
248.222(2283 3,,444
1438. 11381175112278.55 7378.339)
102.551(6853.100) 109)
725.048 228)
13835. (8433884912..5
N=30
NA
NA
3381 10T.999599999.. .46)
2944 (1444284(9372.551) (2..525510)
2946. (63.16009937880 893)
1447 017744866778822 234)
1090.3 000.37004677.81777
38 08.819(214.390 4,990)
2.597(7 (.597(7.894) .944)
9.611 1153)
3.782 6433)
3. 3.987(17.66) .698)
N=25
NA
NA
438. 088.756(1797)) (97.4462)
304.2460 246( 1144246((133.77)) 13377.500, (777750)))
304.281 (14.281(1333.6666 (685)
1787 (1..8666333333 1777 177)
56. 0..368(357.666 .663)
40.077 703)
0.658 (1.00)
0.84 0.845(2.916)
0.477(2.357)
0.309(2.049)
N=20
NA
NA
7210.262(1 10225.22220 188)
521 001 (20.08222222777 726)
521 2481 221.248(2233.72) 2243.762 769)
793. 193.739(3287.277
83..86)6))) 9866
24.316(83.98))
530.289(3421.044
1.877(3.650)
2.974(9.078)
2.764(9.
2.737(9.016)
N=15
NA
NA
15688. 598( 70411.0445535 15.488)
15688. (08.665((04411 704155362 1.662)
15920. 29w.222222222222 469)
15690. 000.6677770 36745555..316
276. 276.745(980.00) 0334)
53.093(154.100) 109)
4.962(17.396)
2.572(9.418)
3 3.623(14.137)
2.863(9.559)
N=10
NA
NA
Kernel
RW
SM
WL
WLOAA
SSP
SP
ESSP
ESP
RW*
SM
WL*
WLOA*
28

---

Published as a conference paper at ICLR 2026
dynamic a
A 1196
-0.2
best
-0.3
-0.5
10
20
40
50
Figure 9: Bayesian optimization results on QM7 with N = 10 and different values of 3//22 Best objective value is plotted at each iteration. Mean with 0.5 standard deviation over 10 replications is reported.
Uimende
0.0
Ren
simple rogret
Simple
nter
iter
iter
(a) QM7, N =
(b) QM7, N
(c) QM7, N =
simple regiet
simple regret
Ttar
Itar
Itar
(d) QM7, 7
(c) QM9, N 4
(f) QM9, N
simple tegret
simple reavet
simple
60
(g) QM9,
(h) QM9, =7
QM9,
Figure 10: Performance of random sampling and Limeade over QM7 and QM9 datasets with different graph size N. Simple regret is plotted at each iteration. Mean with 0.5 standard deviation over 10 replications is reported.
29

---

Published as a conference paper at ICLR 2026
Limeado
lime
SSP
SP
ESP
40
best y
-10
-13
30
45
Iter
iter
(a) QM7, N = 15
(b) QM7, N 25
0.2
best y
best y
-01
-02
-03
iter
iter
(c) QM9,
(d) QM9, N
Figure 11: Bayesian optimization results QM7 and QM9 with N E (15, 25). Best objective value
plotted at each iteration. Mean with 0.5 standard deviation over 10 replications is reported.
30


Reviews:
[{"review_id": "Xfx0PRfCAe", "reviewer": "Reviewer_shL1", "summary": "This paper introduces BoGrape, a novel framework for Bayesian optimization (BO) over graph-structured domains. The key idea is to encode shortest-path relations in a mixed-integer programming (MIP) formulation, enabling global optimization of acquisition functions over combinatorial graph spaces while incorporating explicit structural constraints. The method is evaluated on synthetic benchmarks and molecular design tasks, showing superior performance to several baselines including random search, evolutionary algorithms, and existing graph BO methods.", "questions": "1. How would BoGrape perform on larger graphs (e.g., n>50)?\n\n2. In line 52, the paper states that existing methods are “incapable of efficiently exploring the search domain.”However, based on later sections, the proposed method does not appear to demonstrate clear improvements in efficiency. \nCould the authors provide quantitative comparisons of training and inference time costs for all methods across different tasks to substantiate this claim?\n\n3. Why were recent neural graph optimization methods not included for a more comprehensive and up-to-date comparison? For instance, in the molecular experiments (Section 4.3), it is recommended that the authors extend the comparison beyond Bayesian optimization and evolutionary algorithms to include recent domain-specific neural graph optimization methods. Incorporating these approaches would provide a more comprehensive and up-to-date evaluation and better demonstrate the practical potential of the proposed method in real-world applications.", "limitations": "", "rating": 4, "confidence": 4, "soundness": 2, "presentation": 3, "contribution": 3, "strengths": "1.The proposed approach is novel. The paper presents the first unified framework that formulates BO directly over graph structures. The use of MIP for global acquisition optimization is technically innovative and allows explicit incorporation of hard structural constraints, which is critical for applications such as molecular design.\n\n2.The authors provide formal proofs ensuring the feasible domain of the formulation is equivalent to the graph space consisting of all connected graphs. This gives BoGrape a solid theoretical grounding compared to heuristic approaches.\n\n3.The paper is generally well-written, with a clear motivation and logical flow from problem definition to methodology and results.", "weaknesses": "Overall, the main issues of the paper lie in two aspects.\n\nScalability concerns:\nThe proposed MIP formulation scales poorly with the number of nodes O(n3) variables due to shortest-path encoding. The largest experiments involve small graphs (n < 20), limiting its practical applicability. See questions 1 and 2.\n\nLimited baselines:\nComparisons mainly include random, evolutionary, and one BO-based baseline.\nMissing comparisons with recent neural graph optimization methods. see question 3.", "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": "yGJpqxDj6X", "reviewer": "Reviewer_bwgA", "summary": "This paper proposes a graph-based Bayesian optimization that efficiently solves graph-level optimization problems using shortest-path encoding and mixed-integer programming (MIP). It establishes a bijective mapping between shortest-path MIP encoding and the set of connected graphs, thereby addressing the challenge of modeling graph structural constraints. Additionally, it introduces four variants of shortest-path kernels and a mixed kernel formulation to balance structural and feature similarity, and achieves constrained global acquisition optimization over mixed discrete–continuous spaces, with validated effectiveness on both synthetic and molecular design tasks.", "questions": "1. The current MIP encoding only supports undirected connected graphs with no more than 30 nodes, and its scalability to large graphs has not been fully verified. Additionally, the proposed extensions for directed and disconnected graphs lack empirical support. It is recommended that the authors include corresponding validation experiments to enhance the method's generality and practicality.\n2. The current experimental benchmarks do not include domain-specific Bayesian optimization methods, such as those for molecular design. The authors are encouraged to incorporate domain-relevant baselines to more comprehensively evaluate the performance and practical applicability of the proposed approach.", "limitations": "", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. This paper proposes a shortest-path-based MIP encoding strategy that establishes a strict bijection between connected graphs and the optimization constraint space.\n2. It designs four positive-definite graph kernel variants and a mixed kernel formulation to adapt to different graph scales and nonlinear characteristics, integrating structural and attribute similarities.\n3. It achieves global acquisition optimization within the MIP framework, explicitly modeling problem constraints and avoiding the local optima issues of traditional heuristic methods.", "weaknesses": "1. Since the computational complexity of MIP solving increases exponentially with the number of nodes, the current method only supports graphs with ≤30 nodes, making it challenging to handle large-scale graph scenarios.\n2. The experiments mainly focus on undirected connected graphs, with insufficient validation of adaptability to directed and disconnected graphs, and a lack of related experiments, as well as detailed discussion and verification.", "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": "nptLe8cbUu", "reviewer": "Reviewer_Tm62", "summary": "This paper presents BoGrape, a new Bayesian optimization framework designed for black-box functions defined over graphs. The authors introduce four shortest-path-based graph kernels: SP, SSP, ESP, and ESSP; and then use them in BO with a global acquisition optimization scheme formulated as a mixed-integer optimization problem. Experiments are carried out on both synthetic and real-world scenarios, where the authors validate their proposed kernels and use them for BO on QM7 and QM9 datasets.", "questions": "- On Figure 2, why is there no error bar on the first plot for the SSP kernel?\n\n- Instead of comparing among your proposed kernels in Figure 2, I think it will be better to show the performance of other graph kernels in the main text (I notice that they are in Appendix C) by making each subplot a bit smaller. The reason is that the audience will be more interested in how your proposed kernels compare to previous graph kernels, rather than solely comparing among the new ones.\n\n- A minor suggestion: it's better to briefly mention the testing accuracy of the GNNs that are used as ground-truth.\n\n- I am a bit curious about the learned values of $\\alpha$ and $\\beta$ in your kernels, e.g., results in Section 4.1 Model Performance. (I think they should reflect the \"importance\" of graph structure and node feature in the underlying function?)\n\n- Another minor suggestion: please consider using a larger font size for x and y labels, ticks on x and y axes, titles, and legends in Figures 2, 3, and 4.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "- This paper is well-written with clear logic and is very easy to follow.\n\n- Unlike previous works on graph BO, which use evolutionary algorithms to explore the search space, this paper proposes a global optimization framework over the feasible graph space based on mixed Integer programming. This idea sounds very novel to me.\n\n- The proposed graph kernels based on shortest-path look solid to me.\n\n- I am happy to see that the authors compared the performance of their proposed kernels to previous graph kernels (but why are the results in the appendix?)\n\n- I am also glad to see that the authors use baselines based on graph BO rather than simply comparing to random sampling. \n\n- The authors also provide further details on the computational complexity, scalability, and limitations of their proposed framework, which help the audience better understand the pros and cons of this new method.", "weaknesses": "While the authors said their method can be easily applied to other common acquisition functions, they only consider LCB in their methodology and experiments. Is it possible to extend BoGrape to UCB (this one should be straightforward) or EI?\n\nAnother limitation, as mentioned by the authors in the appendix, is the scalability issue when the underlying graph is large. The maximum size of the graphs used in the experiments is 25, but in many real-world applications (e.g., social networks and infrastructure networks), the network size may go up to $10^5$. In these cases, I am unsure whether using a global optimization based on MIP is more efficient than adopting evolutionary algorithms. \n\nBut I don't think this is necessarily a \"weakness\" of the current work, since most of the graph-level optimization algorithms suffer from the scalability problem. Overall, I have a positive evaluation of this work.", "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": "mhdWPx2EzQ", "reviewer": "Reviewer_uSM7", "summary": "- This paper proposes a Bayesian optimization (BO) framework for optimizing black-box functions whose inputs are graphs, which commonly arise in domains such as molecular design, supply chains, and sensor placement.\n- The authors study four variants of the shortest-path graph kernel (SP, SSP, ESP, and ESSP) and integrate them within a Gaussian process surrogate model.\n- The main methodological contribution is the formulation of the acquisition function optimization as a mixed-integer programming (MIP) problem, which enables global optimization using off-the-shelf MIP solvers.\n- The MIP has the acquisition function (in this paper, the authors consider the lower confidence bound, LCB) as its objective and includes constraints that ensure the predictive mean and variance from the GP model are consistent with the graph kernel, as well as additional constraints enforcing valid shortest-path relationships within the graph.\n- The authors theoretically prove that each feasible MIP solution corresponds bijectively to a valid connected graph in the search space.\n- The proposed framework, named BoGrape, is evaluated on both synthetic graph benchmarks and real-world molecular design tasks, demonstrating competitive or superior performance compared to baseline methods.", "questions": "- In the introduction, it is stated that there are two scenarios of graph optimization problems: optimizing over nodes and optimizing over graph structure, and this paper considers only the second case. If the algorithm is not applicable to the first case, this point should be specifically mentioned in the abstract to avoid overclaiming the ability of the framework. If the framework can be generalized to the first scenario, a brief idea on how to do it should be discussed.\n- It is difficult to distinguish between node labels and node features. An explanation through an example would be helpful.\n- The notation $\\beta$ is used twice: in the objective of the MIP in Eq. (1a) and in the graph kernel weights in Eq. (3). I suggest changing the notation for clarity.\n- In the last two constraints of (MIP-SP), what indices are the summations over? $w$?\n- Should the $\\alpha$ and $\\beta$ in Eq. (3) be positive?\n- Line 290: superscript $I$.\n- Line 411: What is NAS? Network Architecture Search?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "- This paper proposes a method to optimize an acquisition function in BO for a graph-input framework. This area of research is under-explored and has only a limited amount of literature.\n- The method proposed here, i.e., formulating the acquisition function optimization problem as a mixed-integer programming (MIP) problem, is novel in the BO literature for graph functions.\n- Formulating the acquisition function as an MIP also generalizes the capability of BO for graph functions, extending it from optimizing functions defined over graph nodes to functions defined over graph structures.\n- The authors also provide a fruitful discussion on how to generalize the framework to other scenarios, such as graph structures with varying but bounded sizes and other nonlinear acquisition functions.", "weaknesses": "- The presentation of the paper needs improvement. Many typos and a wrong reference in the main part cause some confusion. For example, there is no Eq. (5) in the paper, and it seems that the authors were trying to refer to Eq. (MIP-SP).\n- The main limitation of this paper is the high computational cost of solving the mixed-integer programming (MIP) formulation, which may restrict scalability to large graph sizes or complex kernels. More discussion on this limitation is required. This includes considering larger problems and presenting the comparison in a more informative way, for example by showing the Pareto front between acquisition optimization runtime and BO performance.\n- There is a lack of explanation on the choices of parameters used in the competitor methods. For example, for methods that sample random graphs, why did the authors consider only 20 candidate graphs in each iteration? Since random sampling is cheap, it would be nice to see an ablation study with higher numbers.\n- It is unclear how feasibility constraints of valid molecules are encoded in the method; more explanation on this part is needed.\n- The authors considered several kernels but did not discuss the pros and cons of each of them. This discussion would help readers choose kernels that encode the properties relevant to their applications.", "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)