Paper deep dive
Out-Of-The-Loop Multi-Fidelity Bayesian Optimization
Gustavo Sutter, Hao Wang, Luis Ricardez-Sandoval, Pascal Poupart, Agustinus Kristiadi
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 92%
Last extracted: 8/9/2026, 2:56:05 AM
Summary
The paper addresses the 'Out-Of-The-Loop' Multi-Fidelity Bayesian Optimization (OOL-MF-BO) problem, where the highest-fidelity objective function is too expensive to query during the optimization loop. Standard multi-fidelity BO algorithms suffer from irreducible regret in this setting because they cannot capture latent processes unique to the high-fidelity target. The authors propose a solution that incorporates historical high-fidelity data from previous tasks, using a multi-task multi-fidelity kernel with deep embeddings of task descriptors to transfer knowledge and mitigate suboptimality.
Entities (7)
Relation Signals (5)
Out-Of-The-Loop Multi-Fidelity Bayesian Optimization → isavariantof → Multi-Fidelity Bayesian Optimization
confidence 95% · The paper formalizes OOL-MF-BO as a specific challenging scenario within the broader MF-BO framework where the highest fidelity is not in the loop.
Historical Data → mitigates → Out-Of-The-Loop Suboptimality
confidence 93% · We mitigate this problem by incorporating historical high-fidelity data accompanied by task descriptors.
Out-Of-The-Loop Multi-Fidelity Bayesian Optimization → suffersfrom → Irreducible Regret
confidence 92% · Theorem 3.4 states that whenever the unobserved component shifts the location of the optimum, the algorithm incurs irreducible regret.
Linear Model of Coregionalization → causes → Suboptimality
confidence 90% · Figure 2 and Section 3 illustrate the pitfalls of out-of-the-loop multi-fidelity optimization under the LMC kernel.
Deep Kernel → isusedtomodel → Task and Fidelity Relationships
confidence 88% · The deep kernel operates on fidelity indices and task features to learn complex, nonlinear relationships between tasks and fidelities.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Black-box optimization is a ubiquitous problem in science and engineering, often dealing with expensive objective functions with cheaper lower-fidelity proxies available. Multi-fidelity Bayesian optimization (MF-BO) is a principled approach to this problem, leveraging correlations across different fidelities when querying the objective. However, for many important MF-BO tasks, the true highest-fidelity function is prohibitively expensive to be part of the optimization loop. Nevertheless, practitioners often have gold standard data (observations of the highest-fidelity function) obtained from previous experiments that might provide information for the current task. For instance, in molecular optimization, chemists often pick the top-$k$ candidate molecules using various computer simulations, and later reveal their true objective function values. In this work, we demonstrate the suboptimality of standard MF-BO algorithms in the real-world scenarios above, even under ideal assumptions. Next, we mitigate this problem by incorporating historical high-fidelity data accompanied by task descriptors---which can be explicitly given or extracted from unstructured metadata. We demonstrate the effectiveness of our methods on synthetic functions, as well as real-world problems in chemistry and hyperparameter optimization.
Tags
Links
- Source: https://arxiv.org/abs/2608.04113v1
- Canonical: https://arxiv.org/abs/2608.04113v1
PDF not stored locally. Use the link above to view on the source site.
Full Text
66,777 characters extracted from source content.
Expand or collapse full text
Out-Of-The-Loop Multi-Fidelity Bayesian Optimization Gustavo Sutter1,2,†, Hao Wang3, Luis Ricardez-Sandoval3,4, Pascal Poupart1,2, Agustinus Kristiadi2,5 Abstract Black-box optimization is a ubiquitous problem in science and engineering, often dealing with expensive objective functions with cheaper lower-fidelity proxies available. Multi-fidelity Bayesian optimization (MF-BO) is a principled approach to this problem, leveraging correlations across different fidelities when querying the objective. However, for many important MF-BO tasks, the true highest-fidelity function is prohibitively expensive to be part of the optimization loop. Nevertheless, practitioners often have gold standard data (observations of the highest-fidelity function) obtained from previous experiments that might provide information for the current task. For instance, in molecular optimization, chemists often pick the top-k candidate molecules using various computer simulations, and later reveal their true objective function values. In this work, we demonstrate the suboptimality of standard MF-BO algorithms in the real-world scenarios above, even under ideal assumptions. Next, we mitigate this problem by incorporating historical high-fidelity data accompanied by task descriptors—which can be explicitly given or extracted from unstructured metadata. We demonstrate the effectiveness of our methods on synthetic functions, as well as real-world problems in chemistry and hyperparameter optimization. $ $$ $footnotetext: Corresponding author: gsutterp@uwaterloo.ca 1 Introduction Figure 1: Visualization of the out-of-the-loop multi-fidelity setting, emphasizing the distinction between the in-the-loop and out-of-the-loop fidelities. (I) The online Bayesian optimization loop actively queries lower-fidelity functions. (I) The offline stage takes the final top-k suggestions out of the loop to be observed on the highest-fidelity objective. A wide array of problems in chemistry (Griffiths and Hernández-Lobato 2020; Greenaway et al. 2023; Muthyala et al. 2025), biology (Romero et al. 2013; Ruberg et al. 2023; Martens et al. 2025), and engineering (Ament et al. 2023; Feurer et al. 2022; Lam et al. 2018) consists of optimizing black-box functions. Such objectives often have prohibitively large input domains and are expensive to compute. This creates the need for data-efficient black-box optimization algorithms that are able to find the optimal input point by querying the objective function as few times as possible. Bayesian optimization (BO; Močkus 1975; Garnett 2023) is a principled way to solve this problem following the principles of Bayesian decision theory. One solution to expensive objective function evaluations is the use of functions with different cost-fidelity trade-offs as proxies for the true, prohibitive objective. This setting is addressed by multi-fidelity Bayesian optimization (MF-BO; Kennedy and O’Hagan 2000; Huang et al. 2006) algorithms, where costs and correlations between proxies and the true objective are considered when selecting the next point to query in the BO loop. For example, one can employ numerical simulators or machine learning models instead of performing wet-lab experiments (Guan et al. 2022). However, in many applications, the highest fidelity level is prohibitively expensive, so practitioners often do not include it in the BO loop. Instead, the highest level available during online exploration is not the exact gold standard that is being optimized. Usually, practitioners run optimization on the best available fidelity and evaluate the obtained maximizers on the real objective function of interest. Notably, the results of previous experiments are saved for future reference and experiments, which can be exploited by future optimization campaigns. For example, the highest fidelity in material design is the physical experiment in the laboratory, but often BO is run considering density functional theory (DFT; Argaman and Makov 2000) as the objective, with lower fidelities corresponding to lighter DFT configurations and/or machine learning interatomic potential models, such as UMA (Wood et al. 2025). Only once the BO algorithm is done, its top suggestions are used in lab experiments, the real objective functions are computed, and the results are saved. In this work, we first characterize the suboptimality of standard MF-BO under this setting, even assuming knowledge of the true correlation between fidelities—which is not available in real-world applications. This highlights the need for extra information when the highest fidelity is not available in the loop. Next, we propose a solution based on transferring information from previous optimization tasks to address the suboptimality problem. We incorporate information from offline data collected at the end of past experiments via a multi-task multi-fidelity kernel. Importantly, the proposed deep kernel can leverage task-specific features that are either explicitly provided or extracted from unstructured metadata using modern foundation models (Bommasani et al. 2021), making use of the rich textual context available in real-world applications. We validate our approach on synthetic benchmarks as well as real-world applications in chemistry and hyperparameter tuning, demonstrating the importance of knowledge transfer in the OOL-MF-BO setting. (a) Latent processes (b) Objective functions (c) Target decomposition Figure 2: Illustration of the pitfalls of out-of-the-loop multi-fidelity optimization under the LMC kernel, with M=2M=2. (a) Three latent functions u1,u2,u3u_1,u_2,u_3 sampled from GPs with different lengthscales. (b) Higher and lower fidelity functions, with f(1)=u1+u2f^(1)=u_1+u_2 and f(2)=2u1+u2+u3f^(2)=2\,u_1+u_2+u_3. That is, the process u3u_3 is present in the target but not in the observable fidelity. (c) The observable component f∥(2) f^(2)_\| does not have the same maximizer as the target f(2) f^(2) due to the shift imposed by f⟂(2) f^(2)_ . The contributions of this work are as follows: • We are, to the best of our knowledge, the first to formalize the OOL-MF-BO problem. • We introduce theoretical results characterizing the suboptimality of OOL-MF-BO even if the true correlations between fidelities are known. • We introduce a method that leverages historical experimental data and task-specific features—extracted from either structured metadata or foundation models—to escape the aforementioned suboptimality, leading to a reduction in cumulative regret in multiple benchmarks. • We show empirical results in the OOL-MF-BO setting across standard synthetic functions and on real-world chemistry and hyperparameter optimization benchmarks. 2 Preliminaries Let f:→ℝf:X→R denote an unknown objective function on a space ⊂ℝdX⊂R^d. The goal of black-box optimization is to find an optimal point ∗∈argmax∈f() x_*∈ *arg\,max_ x∈Xf( x) while assuming that f is expensive to compute and the input domain cannot be explored exhaustively. The objective function is accessed via noisy observations y=f()+ϵy=f( x)+ε with ϵ∼(0,σn2)ε N(0, _n^2), stored in the observation dataset t=(i,yi)i=1tD_t= \( x_i,y_i)\_i=1^t. Bayesian Optimization Bayesian optimization (Močkus 1975; Garnett 2023) provides a principled solution for black-box optimization problems based on two main components: a surrogate model p(f|t)p(f|D_t) and an acquisition function α(;t)α( x;D_t). The surrogate model represents the current posterior belief over the unknown objective function, which is used by the acquisition function to evaluate candidate input locations and select the next query point according to t+1=argmaxα(;t) x_t+1= *arg\,max_Xα( x;D_t). The optimization loop terminates after a predefined number of iterations T is reached. Then, the algorithm returns its final choice ^T=argmaxi∈Tyi x_T= *arg\,max_ x_i _Ty_i, given by the observed point with the highest objective value. Multi-Fidelity Bayesian Optimization The multi-fidelity extension of the Bayesian optimization framework assumes that we have access to several functions in f(m):→ℝm=1M \f^(m):X \_m=1^M with the highest fidelity function f(M)f^(M) being the objective to be maximized, i.e., we aim to find ∗∈argmax∈f(M)() x_*∈ *arg\,max_ x∈Xf^(M)( x) The observations dataset is a collection of triples (i,yi,mi)( x_i,y_i,m_i), where mi∈1,…,Mm_i∈\1,...,M\ indicates the fidelity queried on step i. For each fidelity, there is an associated querying cost c(m)≥0c^(m)≥ 0 (assumed to be independent of x). In the multi-fidelity regime, the acquisition function is given by α(,m;t)α( x,m;D_t), which also incorporates the fidelity as an input. The algorithm runs until the optimization budget Λ is exhausted. Both the surrogate model and the acquisition function are adapted to deal with a set of correlated functions. Analogous to the single-fidelity case, at termination, the algorithm makes a final suggestion. Usually, this is done by suggesting the observed point from the highest fidelity with the highest objective value, i.e. ^T=argmaxi∈T,myi x_T= *arg\,max_ x_i _T,my_i where T,m=(i,yi,mi):mi=MD_T,m=\( x_i,y_i,m_i):m_i=M\. Another option is to suggest ^T=argmax∈[f(M)()|T] x_T= *arg\,max_ x∈XE[f^(M)( x)|D_T], allowing for points that were not observed in the highest fidelity level or not observed at all. Gaussian Processes Gaussian processes (GPs; Rasmussen and Williams 2006) are the standard choice for surrogate modelling in BO. They define a distribution over functions specified by a mean function μ:→ℝμ:X and a covariance kernel k:×→ℝk:X×X . Given a dataset tD_t with inputs and noisy observations, a GP prior f∼(μ,k)f (μ,k) induces a multivariate Gaussian distribution over the unseen function values and observed data tD_t. Conditioning on the data yields a GP posterior whose predictive distribution at a test point x is Gaussian with closed-form mean μt() _t( x) and variance σt() _t( x). Vector-valued Gaussian Process When dealing with vector-valued functions f:→ℝMf:X ^M, it is still possible to use GPs. Under this setting, it is assumed that f(x)=(f1(x),…,fM(x))⊤f(x)=(f_1(x),…,f_M(x)) follows a vector-valued GP (Álvarez et al. 2012), fully characterized by a mean function :→ℝM μ:X ^M and a matrix-valued kernel :×→ℝM×MK:X×X ^M× M. The entries ((,′))m,m′(K( x, x ))_m,m in the matrix (,′)K( x, x ) correspond to the covariance between f(m)()f^(m)( x) and f(m′)(′)f^(m )( x ). A popular choice is the linear model of coregionalization (LMC; Journel and Huijbregts 1978; Goovaerts 1997): [(,′)]m,m′=∑q=1Qm,m′(q)kq(,′),[K( x, x )]_m,m = _q=1^Q B^(q)_m,m k_q( x, x ), (1) where (q)∈ℝM×M B^(q) ^M× M is a positive definite coregionalization matrix and kqk_q is a kernel for q∈1,…,Qq∈\1,…,Q\. This corresponds to modelling each output as a linear combination of latent processes, f(m)()=∑q=1Qaq(m)uq()f^(m)( x)= _q=1^Qa^(m)_qu_q( x), where uq∼(0,kq)u_q (0,k_q) are i.i.d. latent processes. When all latent processes have the same kernel, the resulting model is called an intrinsic coregionalization model (ICM; Goovaerts 1997), for which the covariance is: [(,′)]m,m′=m,m′k(,′),[K( x, x )]_m,m = B_m,m k( x, x ), (2) Another popular kernel, the multi-information source (MISO) model of Poloczek et al. (2017), can be written in the LMC form by taking Q=MQ=M and setting f(M)=u1f^(M)=u_1 and f(m)=u1+um+1f^(m)=u_1+u_m+1 for m=1,…,M−1m=1,…,M-1. Note that this requires distinct kernels for each process, placing it in the LMC family but outside ICM. 3 Pitfalls of Out-of-the-Loop Fidelity Figure 3: Proposed solution for the OOL-MF-BO setting based on data from historical tasks for which the highest fidelity was observed. The surrogate model has access to a dataset of historical tasks alongside the data being collected from the available fidelities of the current task. We start by characterizing what happens when we only have access to the lower fidelities in the BO loop. As noted before, in many complex scientific applications the true objective function is prohibitively expensive—both computationally and financially—to be accessed within the loop. Consider a multi-fidelity optimization problem with functions f(m)m=1M \f^(m)\_m=1^M drawn from a joint Gaussian process with zero mean and an LMC kernel, that is, each fidelity is a linear combination of the same set of latent processes uqq=1Q \u_q\_q=1^Q. Let the costs c(m)m=1M \c^(m)\_m=1^M be such that c(M)=∞c^(M)=∞, and let the budget be finite, Λ>0 >0. Therefore, we have mi<Mm_i<M for all elements in t=(i,yi,mi)i=1ND_t=\( x_i,y_i,m_i)\_i=1^N. The final suggestion of the algorithm is given by ^T=argmax∈[f(M)()∣T], x_T= *arg\,max_ x∈XE[f^(M)( x) _T], (3) since we are unable to query our target fidelity. First, we partition the set of latent functions based on their presence in the lower fidelities. Next, using the two sets in this partition, we decompose the objective function into observable and unobservable components. Definition 3.1. Let obs=q|aq(m)≠0for somem<MS_obs=\q\;|\; a^(m)_q≠ 0\;for some\;m<M\ denote the set of indices of latent processes that are observable in at least one of the lower fidelities. In addition, define its complement obs=q|aq(m)=0for allm<MS_obs c=\q\;|\; a^(m)_q=0\;for all\;m<M\, capturing the latent processes that are not present in any observable fidelity. Definition 3.2. Let the target fidelity be f(M)()=∑q=1Qaq(M)uq()f^(M)( x)= _q=1^Qa^(M)_qu_q( x). Splitting this sum according to the partition obs,obs\S_obs,S_obs c\ of 1,…,Q\1,…,Q\ from Definition˜3.1 yields f(M) f^(M) =f∥(M)+f⟂(M),where = f^(M)_ + f^(M)_ , f∥(M) f^(M)_ :=∑q∈obsaq(M)uq,f⟂(M):=∑q∈obsaq(M)uq. = _q∈S_obsa^(M)_qu_q, f^(M)_ = _q∈S_obs ca^(M)_qu_q. This decomposition of the objective function allows us to reason about what can be reconstructed from out-of-the-loop observations. The key object to study under these conditions is the posterior mean μt(M)=[f(M)∣t]μ^(M)_t=E[f^(M) _t], which is used to make suggestions. Proposition 3.3. Under the LMC kernel defined in (1), the highest-fidelity posterior mean equals the posterior mean of the observable component: μt(M)=[f∥(M)∣t]μ^(M)_t=E[ f^(M)_ _t]. The unobservable component f⟂(M) f^(M)_ contributes nothing. Crucially, this inability to capture information from f⟂(M) f^(M)_ is not resolved as more low-fidelity data is collected. As a result, whenever the unobserved component shifts the location of the optimum, the algorithm incurs irreducible regret. Theorem 3.4. Assume [f∥(M)∣t]→f∥(M)E[ f^(M)_ _t]→ f^(M)_ uniformly as t→∞t→∞ and that the lower-fidelity observations become dense in X. Given the candidate ^t x_t selected via (3), whenever argmaxf(M)≠argmaxf∥(M) *arg\,maxf^(M)≠ *arg\,max f^(M)_ the regret converges to a strictly positive constant: limt→∞(f(M)(∗)−f(M)(^t))>0 _t→∞(f^(M)( x_*)-f^(M)( x_t))>0. Remark 3.5. Under MISO, every latent process appears in some observable fidelity (obs=∅S_obs c= ), therefore f⟂(M)=0 f^(M)_ =0 and Thm.˜3.4 does not hold. However, this is a structural assumption on the fidelity hierarchy, and when it fails, the misspecified prior restores the suboptimality. A concrete instance of the LMC problem is presented in Figure˜2. In this bi-fidelity example, the model has direct access to the true linear coefficients. However, since the lower-fidelity function does not include the higher-frequency component, the posterior collapses to the observed component and the regret converges to a positive value. The suboptimality of the LMC just presented relies on the set of non-observable latents being non-empty. For the ICM setting, however, we show that the regret is likely irreducible even when obs=∅S_obs c= , as the target posterior mean converges to a fixed combination of the lower fidelities that in general does not share the optimum location with the objective. Corollary 3.6. Under the ICM kernel defined in (2), the highest-fidelity posterior mean corresponds to a fixed linear combination of the lower-fidelity posterior means, i.e., μt(M)()=∑m=1M−1wmμt(m)()μ^(M)_t( x)= _m=1^M-1w_mμ^(m)_t( x), where the weights w do not depend on the input location x. Furthermore, if ∗≠argmax∈∑m=1M−1wmf(m)() x_*≠ *arg\,max_ x∈X _m=1^M-1w_mf^(m)( x), the simple regret converges to a strictly positive constant: limt→∞(f(M)(∗)−f(M)(^t))>0 _t→∞(f^(M)( x_*)-f^(M)( x_t))>0. This effect is clear in the bi-fidelity setting, which reduces to the high-fidelity posterior mean being a simple rescaling of the lower-fidelity one. Therefore, whenever the lower-fidelity maximizer or minimizer is not aligned with the true optima, the model will have irreducible regret. 4 Handling Out-of-the-Loop Fidelity The root of the suboptimality stated in the previous section is the absence of the target fidelity information within the loop. Therefore, to escape this regime, it is necessary to incorporate a higher-fidelity signal in the surrogate model. As, by definition, the highest fidelity of the current objective is not available, the solution is to incorporate data from outside the current task. Fortunately, in many science and engineering problems, practitioners do not start their optimization campaign in a vacuum. Instead, they often possess a repository of historical data from previous experiments—datasets where the highest-fidelity outcomes from those experiments have already been observed and recorded. Note, however, that the relationship between those previous experiments and the current experiment is often ambiguous, as past experiments may come from partially relevant configurations or settings. Based on this observation, we propose a solution for the aforementioned problem in OOL-MF-BO. Formally, consider N sets of functions fn(m)m=1M \f^(m)_n\_m=1^M for n=1…Nn=1...N, where the subscript n indicates the task index (i.e., previous experiments), and we are interested in optimizing task N (i.e., current experiment). Assume we have historical observations ℋn=(t,y(mt),mt)t=1TnH_n= \( x_t,y^(m_t),m_t)\_t=1^T_n, for which the highest fidelity is available. We denote by ℋ:=ℋnn=1N−1H:= \H_n\_n=1^N-1 the set of all historical data. In addition, assume each task has a task context vector ene_n that the optimization algorithm can access. We assume that the previous tasks cannot be queried; only their past observations are available. Surrogate modelling To incorporate the information from the previous tasks, we replace the fidelity kernel to also include tasks: Cov(fn(m)(), (f^(m)_n( x), fn′(m′)(′))=k((,m,n),(′,m′,n′)) f^(m )_n ( x ))=k(( x,m,n),( x ,m ,n )) =kdata(,′)⋅kθ((m,n),(m′,n′)). =k_data( x, x )· k_θ((m,n),(m ,n )). The newly introduced term kθ((m,n),(m′,n′))k_θ((m,n),(m ,n )) is a deep kernel (Wilson et al. 2016) that operates on fidelity indices and task features: kθ((m,n),(m′,n′))=kRBF(ℰ(m,en;θ),ℰ(m′,en′;θ)),k_θ((m,n),(m ,n ))=k_RBF(E(m,e_n;θ),E(m ,e_n ;θ)), where ℰE is a neural network parametrized by θ. We use a deep kernel as it can flexibly learn complex, nonlinear relationships between tasks and fidelities. This is particularly relevant in our setting, as the task features—whether manually defined or extracted via foundation models—can interact with fidelity indices in nontrivial ways. This allows the model to reliably identify which information from the previous tasks is relevant to the current objective, as not all previous data is equally relevant. Figure˜3 illustrates how offline data from previous tasks can improve the modelling of the incumbent function and break the linear combination effect previously stated. Data from historical tasks is taken into account in the modelling of the current objective and its lower-fidelity functions. Acquisition function The proposed method does not impose restrictions on the choice of acquisition function. Once the historical dataset ℋH is considered in the surrogate model, any multi-fidelity acquisition function such as MF-MES (Takeno et al. 2020) or MF-EI (Huang et al. 2006) can be used. Importantly, the selected acquisition function should remain formulated to optimize the unobserved highest-fidelity objective. At each step of the optimization loop, the next point and fidelity queried are given by t+1,mt+1=argmax∈,m∈[M−1]α(,m;t,ℋ). x_t+1,m_t+1= *arg\,max_ x∈X,m∈[M-1]α( x,m;D_t,H). The main change from the usual MF-BO is the exclusion of the highest fidelity level from the optimization domain. Task Features Although some optimization problems may have natural task features, this is not the case for all problems. In fact, many real-world problems can be better explained in words than through a hand-crafted feature vector. A natural solution to this is to make use of foundation models trained on a wide variety of modalities and domains. An example of a task description that we use in our experiments is provided in Figure˜4. Notably, we extract only one vector per task, representing a negligible cost. Hydrocarbon task | Solvation energy benchmark Context: Non-polar, hydrophobic solutes. Dominant interactions: weak Van der Waals dispersive forces. Features: Zero dipole moment, lack of hydrogen bond donors or acceptors. Chemical space: Alkanes, alkenes, and aromatics. Solvation driven by cavity formation energy. Figure 4: Task description example from solvation energy benchmark. To generate the task embeddings, the description is fed to an LLM that maps it to a fixed-dimensional vector. 5 Related work Figure 5: Performance comparison of our proposed methods on synthetic benchmark functions. Results are averaged over 20 independent trials with shading indicating the standard error. Multi-task Bayesian Optimization Multi-task BO was originally proposed by Swersky et al. (2013) as a solution to transfer knowledge across tasks to speed up the optimization of a related task. In their work, they use an ICM kernel, learning the correlation between tasks directly by maximum likelihood estimation. Alternatively, Feng et al. (2020) proposed the latent embedding multi-output kernel, in which the external task embeddings are fed to an RBF kernel function. When the task embeddings are extracted by a neural network, this is closely related to deep kernel methods (Wilson et al. 2016; Zhang et al. 2025a); however, in this case, the network weights are frozen. More recently, BOLT (Zeng et al. 2025) was introduced as a framework for using large language models (LLMs) to improve BO transfer across tasks. In their framework, an LLM is fine-tuned to produce a good initialization from task-specific context. Zhang et al. (2025d) is the only work to combine multi-fidelity and multi-task BO. Their work differs from ours in three main ways: (i) the highest fidelity is always observable, (i) the acquisition function used also weighs exploration for future tasks when selecting the next query point, and (i) correlations among tasks arise from the i.i.d. assumption over task parameters rather than an explicit task kernel acting on the task contexts. Bayesian Optimization Without Access to the Objective Function The central characteristic of our problem setting is not having access to the highest fidelity in the BO loop. In contrast, prior work has examined other settings, characterized by different constraints and structural assumptions on the objective function. Zhang et al. (2025b) investigated the setting where, instead of directly selecting a point, one selects an action that is linked to the input space via a conditional distribution. Related to the multi-fidelity setting, Mikkola et al. (2023) explored the problem of unreliable fidelities and proposed a robust algorithm for MF-BO. Furthermore, another related area of work is offline model-based optimization (Kim et al. 2026), which leverages powerful surrogate and generative models to optimize a function in a completely offline setting. In contrast, we study the distinct challenge of OOL-MF-BO, in which the highest fidelity of the objective function is removed from the optimization loop, and lower fidelities are queried directly. 6 Experiments Figure 6: (Top) Performance comparison of all methods on three molecular optimization benchmarks: Xe/Kr selectivity, solvation energy, and polarizability. Results are averaged over 10 trials. (Bottom) Cumulative regret curves for logistic regression and SVM design spaces. Curves are averaged across 80 trials (8 datasets, each with 10 trials). We perform OOL-MF-BO experiments on both synthetic functions, molecular benchmarks, and hyperparameter optimization tasks. To evaluate performance, we visualize the cumulative regret at each step, i.e., Rt=∑τ=1tf(M)(∗)−f(M)(^τ)R_t= _τ=1^tf^(M)( x_*)-f^(M)( x_τ). We consider the following methods: • Random Search: Suggestions are made randomly by picking ^t∼Uniform() x_t (X). • Next-Best: The surrogate models the current task up to the highest observable fidelity. The suggested point is given by ^t=argmax∈[f(M−1)()|t] x_t= *arg\,max_ x∈XE[f^(M-1)( x)|D_t]. • Single-Task: Surrogate models all fidelities of the current task. The suggested point is given by ^t=argmax∈[f(M)()|t] x_t= *arg\,max_ x∈XE[f^(M)( x)|D_t]. • Multi-Task: Our proposed method in which the out-of-the-loop fidelity is maximized. The suggested point is given by ^t=argmax∈[f(M)()|t,ℋ] x_t= *arg\,max_ x∈XE[f^(M)( x)|D_t,H]. • Multi-Task Next-Best: Variation of our proposed method in which the next-best fidelity is maximized. The suggested point is given by ^t=argmax∈[f(M−1)()|t,ℋ] x_t= *arg\,max_ x∈XE[f^(M-1)( x)|D_t,H]. For all methods, except Random Search, we use a GP with a constant mean and an RBF covariance kernel scaled by the dimension (Hvarfner et al. 2024) as the data kernel. For Single-Task and Next-Best, we show results with ICM, LMC, and MISO, covering popular multi-fidelity kernels used in practice. For the Next-Best variants, Appendix C repeats these experiments using the observed incumbent rather than the posterior-mean maximizer, with no qualitative change. The deep kernel network ℰE is parameterized as a two-layer network with 16 hidden units and an 8-dimensional output. The acquisition function used is MF-MES (Takeno et al. 2020). All experiments are implemented using the BoTorch package (Balandat et al. 2020). Synthetic Functions We now evaluate the methods on synthetic functions commonly used in the MF-BO literature. Specifically, given a set of M multi-fidelity objective functions g(m):ℝd→ℝg^(m):R^d for m=1,…,Mm=1,…,M, we generate a family of related tasks by sectioning the original functions along a specific input dimension. Let the input space be decomposed such that the first d′=d−1d =d-1 dimensions represent the design variables ∈ℝd−1 x ^d-1, and the d-th dimension represents a task-defining parameter. By selecting a discrete set of task values =z1,…,zNZ=\z_1,…,z_N\, we define the n-th task at the m-th fidelity level, fn(m):ℝd−1→ℝf^(m)_n:R^d-1 , as fn(m)()=g(m)(x1,…,xd−1,zn) f^(m)_n( x)= g^(m)(x_1,…,x_d-1,z_n). In this framework, znz_n serves as a contextual feature. This construction ensures that tasks are naturally correlated, as they are derived from slices of the same underlying d-dimensional response surface g(m)g^(m). Under this setting, we perform experiments on multi-fidelity versions of eight synthetic functions111https://w.sfu.ca/˜ssurjano/optimization.html: Branin, Michalewicz, Park, Styblinski-Tang, Levy, Ackley, Rosenbrock, and Hartmann. For all functions, we perform 20 trials, in each trial sampling M∼Uniform(2,3,4,5)M (\2,3,4,5\), N∼Uniform(2,3,4,5)N (\2,3,4,5\), |Hn|∼Uniform(20,30,40,50)|H_n| (\20,30,40,50\). Given the number of tasks, the task-defining parameters are sampled uniformly. The costs are set to c(m)=2m−1c^(m)=2^m-1 for m≤M−1m≤ M-1 and c(M)=∞c^(M)=∞, capturing a situation where lowering the fidelity corresponds to halving the data size or tolerance parameter (Eggensperger et al. 2021). Points in the historical observations are sampled uniformly with fidelity ratio 2−(m−1)2^-(m-1). More information is presented in Appendix B. The results are presented in Figure 5. Across all functions, Multi-Task approaches improve the results and are able to perform better than Next-Best and Single-Task methods. Interestingly, the Multi-Task Next-Best method is consistently better than its single-task counterpart. This gap in performance shows the clear advantage of incorporating previous tasks in the OOL-MF-BO setting. As expected, the Single-Task method shows performance comparable to Random Search and, at best, on par with Next-best, highlighting the impossibility of optimizing the unavailable fidelity directly without auxiliary data. For further analysis of the different components of the setting, we refer to the ablation study exploring the effects of number of tasks, task relevance, and historical data size provided in Appendix D. Real-World Applications Molecular Optimization We experiment with three molecular optimization benchmarks, adapting the MF-BO problems used in Sabanza-Gil et al. (2025) to our OOL-MF-BO setting. The problems are (i) Xe/Kr selectivity (Gantzler et al. 2023), (i) solvation energy (Mobley and Guthrie 2014), and (i) polarizability (Ghahremanpour et al. 2018). For all benchmarks M=2M=2 and the number of steps is set to T=30T=30 for all problems. We use the Qwen3 Embedding 8B (Zhang et al. 2025c) model to extract task features from text descriptions of the tasks. Figure˜4 shows one of the textual descriptions used. Refer to Appendix B for information on the benchmarks and tasks. Hyperparameter Optimization We evaluate the method on two hyperparameter optimization (HPO) benchmarks from HPOBench (Eggensperger et al. 2021): (i) logistic regression (LR) and (i) support vector machine (SVM). For each benchmark, the hyperparameter space is provided alongside the validation accuracy for different fidelity levels, with M=5M=5 and M=3M=3 for LR and SVM, respectively. Following the complexity scaling of the benchmark fidelities, costs are c(m)=3c(m−1)c^(m)=3c^(m-1) with c(1)=1c^(1)=1 and the budget is set to Λ=50 =50. In both benchmarks, tasks correspond to different training and validation datasets, for a total of N=8N=8 tasks. We perform 10 trials per dataset in the role of the current task, randomly sampling only 3 of the remaining tasks to be present in the historical data for each trial. Refer to Appendix B for more information. The cumulative regret results for both families of benchmarks are presented in Figure˜6. Our proposed methods Multi-Task and Multi-Task Next Best show strong performance, obtaining the lowest cumulative regrets for three out of five benchmarks. In the polarizability benchmark, the observed gains start from the strong initialization due to the knowledge transfer from historical tasks and are maintained thereafter—resulting in significantly lower cumulative regret throughout the iterations. On the other hand, on both HPOBench tasks (LR and SVM) our methods show consistently lower instantaneous regret during the entire run, leading to better results. The benefits from historical data are clear when comparing Multi-Task and Multi-Task Next Best with the ICM versions of Single Task and Next Best. Even though our proposed deep kernel has a structure very similar to the ICM kernel, namely a separable data kernel, the inclusion of historical information makes the results go from equivalent to Random Search to competitive in all evaluated benchmarks. The strong performance of the MISO-based methods follows from Remark˜3.5, and Next-Best LMC on Xe/Kr selectivity is the same effect in another guise: when the maximizers of the target and next-best fidelities are closely aligned, a method that refuses to extrapolate is difficult to improve upon. Their subpar performance on the HPO and synthetic benchmarks shows that this condition is neither universal nor checkable in advance. Gains from modelling the out-of-the-loop fidelity directly are thus realized precisely in the regime characterized in Section˜3. 7 Conclusion In this work, we explored the out-of-the-loop multi-fidelity Bayesian optimization setting. We started from its formal definition and proceeded by characterizing the suboptimality of tackling problems of this family using standard MF-BO methods with both LMC and ICM—even if the fidelity covariance matrix is known. Next, we demonstrated how historical data with the highest fidelity information can be used to improve the performance in the OOL-MF-BO setting. We showed the benefits of the proposed method across both synthetic and real-world objective functions, demonstrating the effectiveness of our approach. Limitations and Future Work Our analysis of the OOL suboptimality is developed for the LMC family of surrogates, which, although broad, does not cover all possible kernels; extending it to more general covariance structures remains open. Furthermore, performance depends on the relevance of the available historical tasks, which in practice is hard to know in advance. A natural next step is to apply OOL-MF-BO to challenging science and engineering problems using domain-specific foundation models as low-fidelity functions, and to develop more efficient methods for handling historical data that avoid the cubic complexity of GPs. References M. A. Álvarez, L. Rosasco, and N. D. Lawrence (2012) Kernels for Vector-Valued Functions: A Review. Foundations and Trends in Machine Learning 4 (3). Cited by: §2. S. Ament, A. Witte, N. Garg, and J. Kusuma (2023) Sustainable Concrete via Bayesian Optimization. arXiv 2310 (18288). Cited by: §1. N. Argaman and G. Makov (2000) Density Functional Theory: An Introduction. American Journal of Physics 68 (1). Cited by: §1. M. Balandat, B. Karrer, D. Jiang, S. Daulton, B. Letham, A. G. Wilson, and E. Bakshy (2020) BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization. In NeurIPS, Cited by: Appendix B, Appendix B, §6. R. Bommasani, D. A. Hudson, E. Adeli, R. Altman, S. Arora, S. von Arx, M. S. Bernstein, J. Bohg, A. Bosselut, E. Brunskill, E. Brynjolfsson, S. Buch, D. Card, R. Castellon, N. S. Chatterji, A. S. Chen, K. A. Creel, J. Davis, D. Demszky, C. Donahue, M. Doumbouya, E. Durmus, S. Ermon, J. Etchemendy, K. Ethayarajh, L. Fei-Fei, C. Finn, T. Gale, L. E. Gillespie, K. Goel, N. D. Goodman, S. Grossman, N. Guha, T. Hashimoto, P. Henderson, J. Hewitt, D. E. Ho, J. Hong, K. Hsu, J. Huang, T. F. Icard, S. Jain, D. Jurafsky, P. Kalluri, S. Karamcheti, G. Keeling, F. Khani, O. Khattab, P. W. Koh, M. S. Krass, R. Krishna, R. Kuditipudi, A. Kumar, F. Ladhak, M. Lee, T. Lee, J. Leskovec, I. Levent, X. L. Li, X. Li, T. Ma, A. Malik, C. D. Manning, S. P. Mirchandani, E. Mitchell, Z. Munyikwa, S. Nair, A. Narayan, D. Narayanan, B. Newman, A. Nie, J. C. Niebles, H. Nilforoshan, J. F. Nyarko, G. Ogut, L. Orr, I. Papadimitriou, J. S. Park, C. Piech, E. Portelance, C. Potts, A. Raghunathan, R. Reich, H. Ren, F. Rong, Y. H. Roohani, C. Ruiz, J. Ryan, C. R’e, D. Sadigh, S. Sagawa, K. Santhanam, A. Shih, K. P. Srinivasan, A. Tamkin, R. Taori, A. W. Thomas, F. Tramèr, R. E. Wang, W. Wang, B. Wu, J. Wu, Y. Wu, S. M. Xie, M. Yasunaga, J. You, M. A. Zaharia, M. Zhang, T. Zhang, X. Zhang, Y. Zhang, L. Zheng, K. Zhou, and P. Liang (2021) On the Opportunities and Risks of Foundation Models. arXiv 2108 (07258). Cited by: §1. K. Eggensperger, P. Müller, N. Mallik, M. Feurer, R. Sass, A. Klein, N. Awad, M. Lindauer, and F. Hutter (2021) HPOBench: a Collection of Reproducible Multi-Fidelity Benchmark Problems for HPO. In Thirty-fifth Conference on Neural Information Processing Systems Datasets and Benchmarks Track (Round 2), Cited by: Appendix B, §6, §6. Q. Feng, B. Letham, H. Mao, and E. Bakshy (2020) High-Dimensional Contextual Policy Search with Unknown Context Rewards using Bayesian Optimization. In NeurIPS, Cited by: §5. M. Feurer, K. Eggensperger, S. Falkner, M. Lindauer, and F. Hutter (2022) Auto-sklearn 2.0: hands-free AutoML via Meta-Learning. JMLR 23 (261). Cited by: §1. N. Gantzler, A. Deshwal, J. R. Doppa, and C. M. Simon (2023) Multi-Fidelity Bayesian Optimization of Covalent Organic Frameworks For Xenon/Krypton Separations. Digital Discovery 2 (6). Cited by: Appendix B, §6. R. Garnett (2023) Bayesian Optimization. Cambridge University Press. Cited by: §1, §2. M. M. Ghahremanpour, P. J. Van Maaren, and D. Van Der Spoel (2018) The Alexandria Library, A Quantum-chemical Database of Molecular Properties for Force Field Development. Scientific Data 5 (1). Cited by: Appendix B, §6. P. Goovaerts (1997) Geostatistics for Natural Resources Evaluation. Oxford University Press. Cited by: §2, §2. R. L. Greenaway, K. E. Jelfs, A. C. Spivey, and S. N. Yaliraki (2023) From Alchemist to AI Chemist. Nature Reviews Chemistry 7 (8). Cited by: §1. R. Griffiths and J. M. Hernández-Lobato (2020) Constrained Bayesian Optimization for Automatic Chemical Design Using Variational Autoencoders. Chem. Sci. 11 (2). Cited by: §1. Y. Guan, D. Chaffart, G. Liu, Z. Tan, D. Zhang, Y. Wang, J. Li, and L. Ricardez-Sandoval (2022) Machine Learning in Solid Heterogeneous Catalysis: Recent Developments, Challenges and Perspectives. Chemical Engineering Science 248. Cited by: §1. D. Huang, T. T. Allen, W. I. Notz, and R. A. Miller (2006) Sequential Kriging Optimization Using Multiple-Fidelity Evaluations. Structural and Multidisciplinary Optimization 32 (5). Cited by: §1, §4. C. Hvarfner, E. O. Hellsten, and L. Nardi (2024) Vanilla Bayesian Optimization Performs Great in High Dimensions. In ICML, Cited by: §6. A. G. Journel and C. J. Huijbregts (1978) Mining Geostatistics. Academic Press. Cited by: §2. M. C. Kennedy and A. O’Hagan (2000) Predicting the Output From a Complex Computer Code When Fast Approximations Are Available. Biometrika 87 (1). Cited by: §1. M. Kim, J. Gu, Y. Yuan, T. Yun, Z. Liu, Y. Bengio, and C. Chen (2026) Offline Model-Based Optimization: Comprehensive Review. TMLR. External Links: ISSN 2835-8856 Cited by: §5. R. Lam, M. Poloczek, P. Frazier, and K. E. Willcox (2018) Advances in Bayesian Optimization with Applications in Aerospace Engineering. In 2018 AIAA Non-Deterministic Approaches Conference, Cited by: §1. A. Martens, M. Neufang, A. Butté, M. von Stosch, A. del Rio Chanona, and L. M. Helleckes (2025) Holistic Bioprocess Development Across Scales Using Multi-Fidelity Batch Bayesian Optimization. arXiv 2508 (10970). Cited by: §1. P. Mikkola, J. Martinelli, L. Filstroff, and S. Kaski (2023) Multi-Fidelity Bayesian Optimization with Unreliable Information Sources. In AISTATS, Cited by: Appendix B, Appendix B, §5. D. L. Mobley and J. P. Guthrie (2014) FreeSolv: A Database of Experimental and Calculated Hydration Free Energies, with Input Files. Journal of computer-aided molecular design 28 (7). Cited by: Appendix B, §6. J. Močkus (1975) On Bayesian Methods for Seeking the Extremum. In Optimization Techniques IFIP Technical Conference, Cited by: §1, §2. M. R. Muthyala, F. Sorourifar, T. Tan, Y. Peng, and J. A. Paulson (2025) Generative Multiobjective Bayesian Optimization with Scalable Batch Evaluations for Sample-Efficient De Novo Molecular Design. Industrial & Engineering Chemistry Research 65 (1). Cited by: §1. M. Poloczek, J. Wang, and P. Frazier (2017) Multi-Information Source Optimization. In NIPS, Cited by: §2. C. E. Rasmussen and C. K. I. Williams (2006) Gaussian Processes for Machine Learning. The MIT Press. Cited by: §2. P. A. Romero, A. Krause, and F. H. Arnold (2013) Navigating the Protein Fitness Landscape with Gaussian Processes. Proceedings of the National Academy of Sciences 110 (3). Cited by: §1. S. J. Ruberg, F. Beckers, R. Hemmings, P. Honig, T. Irony, L. LaVange, G. Lieberman, J. Mayne, and R. Moscicki (2023) Application of Bayesian Approaches in Drug Development: Starting a Virtuous Cycle. Nature Reviews Drug Discovery 22 (3). Cited by: §1. V. Sabanza-Gil, R. Barbano, D. Pacheco Gutiérrez, J. S. Luterbacher, J. M. Hernández-Lobato, P. Schwaller, and L. Roch (2025) Best Practices for Multi-Fidelity Bayesian Optimization in Materials and Molecular Research. Nature Computational Science 5 (7). Cited by: Appendix B, §6. K. Swersky, J. Snoek, and R. P. Adams (2013) Multi-Task Bayesian Optimization. In NeurIPS, Cited by: §5. S. Takeno, H. Fukuoka, Y. Tsukada, T. Koyama, M. Shiga, I. Takeuchi, and M. Karasuyama (2020) Multi-Fidelity Bayesian Optimization with Max-value Entropy Search and its Parallelization. In ICML, Cited by: §4, §6. J. Vanschoren, J. N. van Rijn, B. Bischl, and L. Torgo (2014) OpenML: Networked Science in Machine Learning. SIGKDD Explor. Newsl. 15 (2). Cited by: Appendix B. A. G. Wilson, Z. Hu, R. Salakhutdinov, and E. P. Xing (2016) Deep Kernel Learning. In AISTATS, Cited by: §4, §5. B. M. Wood, M. Dzamba, X. Fu, M. Gao, M. Shuaibi, L. Barroso-Luque, K. Abdelmaqsoud, V. Gharakhanyan, J. R. Kitchin, D. S. Levine, K. Michel, A. Sriram, T. Cohen, A. Das, A. Rizvi, S. J. Sahoo, Z. W. Ulissi, and C. L. Zitnick (2025) UMA: A Family of Universal Models for Atoms. arXiv 2506 (23971). Cited by: §1. S. Xiong, P. Z. G. Qian, and C. F. J. Wu (2013) Sequential Design and Analysis of High-Accuracy and Low-Accuracy Computer Codes. Technometrics 55 (1). Cited by: Appendix B. Y. Zeng, N. Maus, H. T. Jones, J. Tao, F. Wan, M. Der Torossian Torres, C. de la Fuente-Nunez, R. Marcus, O. Bastani, and J. R. Gardner (2025) Large Scale Multi-Task Bayesian Optimization with Large Language Models. arXiv 2503 (08131). Cited by: §5. F. Zhang, T. Desautels, and Y. Chen (2025a) Robust Multi-Fidelity Bayesian Optimization with Deep Kernel and Partition. In AISTATS, Cited by: §5. M. Zhang, S. Bouabid, C. S. Ong, S. Flaxman, and D. Sejdinovic (2025b) Indirect Query Bayesian Optimization with Integrated Feedback. arXiv 2412 (13559). Cited by: §5. Y. Zhang, M. Li, D. Long, X. Zhang, H. Lin, B. Yang, P. Xie, A. Yang, D. Liu, J. Lin, F. Huang, and J. Zhou (2025c) Qwen3 Embedding: Advancing Text Embedding and Reranking Through Foundation Models. arXiv. Cited by: Appendix B, §6. Y. Zhang, S. Park, and O. Simeone (2025d) Multi-Fidelity Bayesian Optimization with Across-Task Transferable Max-Value Entropy Search. IEEE Transactions on Signal Processing 73 (). Cited by: §5. Appendix A Proofs See 3.3 Proof. Using the decomposition of the objective function and the linearity of expectation, the posterior mean can be written as μt(M)=[f(M)∣t]=[f∥(M)+f⟂(M)∣t]=[f∥(M)∣t]+[f⟂(M)∣t]. splitμ^(M)_t&=E[f^(M) _t]\\ &=E[ f^(M)_ + f^(M)_ _t]\\ &=E[ f^(M)_ _t]+E[ f^(M)_ _t]. split By Definition˜3.2, the second term can be written as [f⟂(M)∣t]=∑q∈obsaq(M)[uq∣t]E[ f^(M)_ _t]= _q∈S_obs ca^(M)_qE[u_q _t]; however, uq⟂tu_q _t for all q∈obsq∈S_obs c. Therefore, [uq∣t]=[uq]=0E[u_q _t]=E[u_q]=0, given that all latent processes have a zero-mean prior. ∎ See 3.4 Proof. By Prop.˜3.3, the posterior mean of the target fidelity evaluates exactly to the posterior mean of the observable component, meaning [f(M)()∣t]=[f∥(M)()∣t]E[f^(M)( x) _t]=E[ f^(M)_ ( x) _t]. Consequently, the candidate point selected at step t via (3) can be equivalently written as ^t=argmax∈[f∥(M)()∣t]. x_t= *arg\,max_ x∈XE[ f^(M)_ ( x) _t]. By the assumption of uniform convergence, [f∥(M)∣t]→f∥(M)E[ f^(M)_ _t]→ f^(M)_ as t→∞t→∞. Since the lower-fidelity observations become dense in X, the maximizer of the posterior mean converges to the maximizer of the limiting function. Let ∥∗=argmax∈f∥(M)() x^*_ = *arg\,max_ x∈X f^(M)_ ( x). It follows that limt→∞^t=∥∗. _t→∞ x_t= x^*_ . Let ∗=argmax∈f(M)() x_*= *arg\,max_ x∈Xf^(M)( x) be the true global optimum. We can now evaluate the limit of the simple regret as t→∞t→∞. Assuming f(M)f^(M) is continuous, we have limt→∞(f(M)(∗)−f(M)(^t))=f(M)(∗)−f(M)(∥∗). _t→∞ (f^(M)( x_*)-f^(M)( x_t) )=f^(M)( x_*)-f^(M)( x^*_ ). By the premise that argmaxf(M)≠argmaxf∥(M) *arg\,maxf^(M)≠ *arg\,max f^(M)_ , the point ∥∗ x^*_ is not a global maximizer of the true objective function f(M)f^(M). Therefore, the value of the target fidelity at the true optimum is strictly greater than its value at the observable optimum: f(M)(∗)>f(M)(∥∗).f^(M)( x_*)>f^(M)( x^*_ ). Subtracting f(M)(∥∗)f^(M)( x^*_ ) from both sides yields f(M)(∗)−f(M)(∥∗)>0f^(M)( x_*)-f^(M)( x^*_ )>0, which concludes the proof that the simple regret converges to a strictly positive constant. ∎ See 3.6 Proof. We begin by writing the general formula for the posterior mean for all m≤Mm≤ M: μt(m)()=m()⊤(+σn2)−1,μ^(m)_t( x)= k_m( x) ( K+ _n^2 I)^-1 y, where m()=[Cov(f(m)(),f(mi)(i))]i=1N k_m( x)=[Cov(f^(m)( x),f^(m_i)( x_i))]_i=1^N is the cross-covariance vector, =[Cov(f(mi)(i),f(mj)(j))]i=1,j=1N,N K=[Cov(f^(m_i)( x_i),f^(m_j)( x_j))]_i=1,j=1^N,N is the covariance matrix of the observation dataset, and =[yi]i=1N y=[y_i]_i=1^N is a vector collecting all observed responses. Recall that, as we are dealing with the ICM kernel, Cov(f(m)(),f(m′)(′))=m,m′k(,′)Cov(f^(m)( x),f^(m )( x ))= B_m,m k( x, x ). Let =<M,<M−1<M,M w= B_<M,<M^-1 B_<M,M, where <M,<M B_<M,<M is the block of the coregionalization matrix corresponding to the observable fidelities and <M,M B_<M,M is the block of cross-terms between the observable and the non-observable fidelity. Since B is positive definite, <M,<M−1 B_<M,<M^-1, and hence w, are well defined. Note that this definition of w implies ∑m<Mwmm,m′=M,m′ _m<Mw_m B_m,m = B_M,m for all observable m′m . Because mi<Mm_i<M holds for all points in the dataset tD_t, we have M()=∑m<Mwmm(). k_M( x)= _m<Mw_m k_m( x). (4) Finally, we write the posterior mean of the highest fidelity and substitute (4) for the cross-covariance term: μt(M)()=M()⊤(+σn2)−1=[∑m<Mwmm()]⊤(+σn2)−1=∑m<Mwm[m()⊤(+σn2)−1]=∑m<Mwmμt(m)(), splitμ^(M)_t( x)&= k_M( x) ( K+ _n^2 I)^-1 y\\ &= [ _m<Mw_m k_m( x) ] ( K+ _n^2 I)^-1 y\\ &= _m<Mw_m [ k_m( x) ( K+ _n^2 I)^-1 y ]\\ &= _m<Mw_mμ^(m)_t( x), split which concludes the first part of the proof. It remains to establish the irreducible regret. Assume that, as t→∞t→∞, each lower-fidelity posterior mean converges uniformly to its function, μt(m)→f(m)μ^(m)_t→ f^(m) for all m<Mm<M, which holds as the lower-fidelity observations become dense in X. Combining this with the first part of the proof, the highest-fidelity posterior mean converges to the fixed combination μt(M)()→∑m<Mwmf(m)()=:g().μ^(M)_t( x)→ _m<Mw_mf^(m)( x)=:g( x). Consequently, the candidate point selected via (3) satisfies limt→∞^t=argmax∈g()=:g∗. _t→∞ x_t= *arg\,max_ x∈Xg( x)=: x^*_g. Let ∗=argmax∈f(M)() x_*= *arg\,max_ x∈Xf^(M)( x) be the true global optimum. Assuming f(M)f^(M) is continuous, the limit of the simple regret is limt→∞(f(M)(∗)−f(M)(^t))=f(M)(∗)−f(M)(g∗). _t→∞ (f^(M)( x_*)-f^(M)( x_t) )=f^(M)( x_*)-f^(M)( x^*_g). By the premise ∗≠argmax∈∑m<Mwmf(m)()=g∗ x_*≠ *arg\,max_ x∈X _m<Mw_mf^(m)( x)= x^*_g, the point g∗ x^*_g is not a global maximizer of f(M)f^(M), so f(M)(∗)>f(M)(g∗)f^(M)( x_*)>f^(M)( x^*_g). Therefore f(M)(∗)−f(M)(g∗)>0f^(M)( x_*)-f^(M)( x^*_g)>0, concluding the proof. ∎ Appendix B Experimental Details Synthetic Functions Here, we provide the formula for all the multi-fidelity functions considered in our experiments. Some multi-fidelity extensions are taken from previous work (indicated accordingly), while others are introduced in this study. All functions have a continuous fidelity parameter m∈[0,1]m∈[0,1], where m=1m=1 recovers the original function. For each benchmark, the fidelity is discretized into M randomly chosen levels. Throughout, we define the fidelity shift δ(m)=0.1(1−m).δ(m)=0.1(1-m). Multi-fidelity Branin (Mikkola et al. 2023): f(,m) f( x,m) =(x2−(b−δ(m))x12+cx1−r)2 = (x_2-(b-δ(m))x_1^2+cx_1-r )^2 +10(1−t)cos(x1)+10, +0(1-t) (x_1)+0, where ∈[−5,10]×[0,15],b=5.14π2,c=5π,r=6,t=18π. gathered x∈[-5,10]×[0,15],\\ b= 5.14π^2, c= 5π,\\ r=6, t= 18π. gathered Multi-fidelity Hartmann (Mikkola et al. 2023): f(,m) f( x,m) =−(α1−δ(m))exp(−∑j=16A1j(xj−P1j)2) =-( _1-δ(m)) \! (- _j=1^6A_1j(x_j-P_1j)^2 ) −∑i=24αiexp(−∑j=16Aij(xj−Pij)2). - _i=2^4 _i \! (- _j=1^6A_ij(x_j-P_ij)^2 ). where ∈[0,1]6,=(1.0,1.2,3.0,3.2), x∈[0,1]^6, α=(1.0,1.2,3.0,3.2), and A and P are the standard Hartmann constants. Multi-fidelity Park (Xiong et al. 2013): f(,m) f( x,m) =x12(1+x4(x2+x32)x12−1) = x_12 ( 1+ x_4(x_2+x_3^2)x_1^2-1 ) +(x1+(3−1.5δ(m))x4)exp(1+sinx3). + (x_1+(3-5δ(m))x_4 ) (1+ x_3). where ∈[0,1]4 x∈[0,1]^4. Multi-fidelity Rosenbrock: (adapted from BoTorch (Balandat et al. 2020)) f(,m)=∑i=1d−2[100(xi+1−xi2+δ(m))2+(xi−1+δ(m))2].f( x,m)= _i=1^d-2 [100(x_i+1-x_i^2+δ(m))^2\\ +(x_i-1+δ(m))^2 ]. where ∈[−5,10]d x∈[-5,10]^d. Multi-fidelity Ackley: f(,m) f( x,m) =−aexp(−b1d−1∑i=1d−1(xi−δ(m))2) =-a \! (-b 1d-1 _i=1^d-1(x_i-δ(m))^2 ) −exp(1d−1∑i=1d−1cos(c(xi−δ(m)))) - \! ( 1d-1 _i=1^d-1 \! (c(x_i-δ(m)) ) ) +a+e, +a+e, where ∈[−32.768,32.768]d,(a,b,c)=(20,0.2,2π). x∈[-32.768,32.768]^d, (a,b,c)=(20,0.2,2π). Multi-fidelity Levy: f(,m) f( x,m) =sin2(πw1) = ^2(π w_1) +∑i=1d−2(wi−1)2(1+10sin2(πwi+1)) + _i=1^d-2(w_i-1)^2 (1+0 ^2(π w_i+1) ) +(wd−1−1)2(1+sin2(2πwd−1)), +(w_d-1-1)^2 (1+ ^2(2π w_d-1) ), where ∈[−10,10]d,wi=1+xi−1+δ(m)4. x∈[-10,10]^d, w_i=1+ x_i-1+δ(m)4. Multi-fidelity Michalewicz: f(,m) f( x,m) =−∑i=1d−1sin(xi−δ(m)) =- _i=1^d-1 (x_i-δ(m)) ×[sin(i(xi−δ(m))2π)]2A, × [ \! ( i(x_i-δ(m))^2π ) ]^2A, where ∈[0,π]d,A=10. x∈[0,π]^d, A=10. Multi-fidelity Styblinski–Tang: f(,m)=12∑i=1d−1[(xi−δ(m))4−16(xi−δ(m))2+5(xi−δ(m))].f( x,m)= 12 _i=1^d-1 [(x_i-δ(m))^4\\ -16(x_i-δ(m))^2+5(x_i-δ(m)) ]. where ∈[−5,5]d x∈[-5,5]^d. Chemistry Applications The three benchmarks are taken from Sabanza-Gil et al. (2025), where the high- and low-fidelity values are provided directly in tabular form. For each dataset, we partition it into separate tasks and generate a textual description of each task using Gemini 3 Pro222https://gemini.google.com/. Such textual descriptions are passed to Qwen3 Embedding 8B (Zhang et al. 2025c), producing task context vectors. For the three benchmarks, we use 30 examples per historical task, sampled uniformly on each trial. We now explain the partitions and provide the task descriptions for each benchmark. Solvation energy (Mobley and Guthrie 2014): For this dataset, we split the tasks based on the functional groups present in the molecules. The four resulting tasks are Nitrogenous, Oxygenated, Halogenated, and Hydrocarbon. We establish a priority hierarchy (Nitrogen > Oxygen > Halogens) to assign each molecule to one of four mutually exclusive tasks. The following task descriptions are used: Hydrocarbon: “Context: Non-polar, hydrophobic solutes. Dominant interactions: weak Van der Waals dispersive forces. Features: Zero dipole moment, lack of hydrogen bond donors or acceptors. Chemical space: Alkanes, alkenes, and aromatics. Solvation driven by cavity formation energy.” Oxygenated: “Context: Polar, hydrophilic solutes. Dominant interactions: Strong hydrogen bonding and permanent dipole-dipole attraction. Features: High electronegativity difference, presence of hydroxyl (OH) or carbonyl (C=O) groups. Chemical space: Alcohols, ethers, ketones.” Nitrogenous: “Context: Basic, polar solutes. Dominant interactions: Proton acceptance and electrostatic contributions. Features: Presence of lone pairs on Nitrogen, variable hybridization (sp2/sp3). Chemical space: Amines, amides, nitriles. Solvation driven by specific H-bond networks.” Halogenated: “Context: Lipophilic, polarizable solutes. Dominant interactions: Electrostatic sigma-hole interactions and halogen bonding. Features: Large atomic radii, high surface area, ’soft’ electron clouds. Chemical space: Alkyl halides. Solvation driven by entropy and polarizability.” Xe/Kr selectivity (Gantzler et al. 2023): We partition this dataset based on the pore diameter of each instance. This results in 3 tasks: Small Pore (Å<15 A<15, Medium Pore (15≤Å≤2515≤ A≤ 25), and Large Pore (Å>25 A>25). The following task descriptions are used: Small Pore: “Context: High confinement regime (< 15 Angstroms). Mechanism: Steric sieving and overlap of potential energy surfaces. Geometry: Dense packing, high framework density, restricted void fraction. Separation driven by repulsive forces and size exclusion.” Medium Pore: “Context: Optimal adsorption regime (15-25 Angstroms). Mechanism: Strong guest-host affinity without steric penalty. Geometry: Mesoporous, balanced void fraction. Separation driven by attractive Van der Waals wells and surface interactions.” Large Pore: “Context: Bulk fluid regime (> 25 Angstroms). Mechanism: Diffusion-dominated transport, weak confinement. Geometry: Open frameworks, extremely high void fraction, low crystal density. Separation driven by pore wall surface area rather than pore size.” Polarizability (Ghahremanpour et al. 2018): For this dataset, we split the tasks based on the elemental composition of the molecules. The four resulting tasks are Hetero-Halogens, CHN Compounds, CHO Compounds, and Hydrocarbons. We established a priority hierarchy (Hetero-Halogens [F, Cl, Br, I, S, P] > Nitrogen > Oxygen) to assign each molecule to one of four mutually exclusive tasks. The following task descriptions are used: Hydrocarbon: “Context: Saturated carbon frameworks. Electronic environment: Hard electron clouds, high HOMO-LUMO gap. Polarizability source: Strictly volumetric scaling (size-dependent). Absence of permanent dipoles or lone pairs.” CHO Compound: “Context: Oxygenated carbon frameworks. Electronic environment: High electronegativity contrast, distorted electron density. Polarizability source: Localized electron density around Oxygen atoms, permanent dipole contributions.” CHN Compound: “Context: Nitrogenous carbon frameworks. Electronic environment: Lone pair donation, resonance effects. Polarizability source: Mobile pi-electrons and orbital hybridization changes (sp/sp2/sp3).” Hetero-Halogens: “Context: Heavy atom substituted frameworks. Electronic environment: Diffuse, soft orbitals (d-orbital participation for S/P). Polarizability source: High deformability of large valence shells (I, Br, S, P). Significant Van der Waals radii overlap.” Figure 7: Synthetic benchmarks with the Next-Best methods reporting the best observed point instead of the posterior-mean maximizer. All other methods are unchanged from Figure˜5. Hyperparameter Optimization Applications The logistic regression (LR) and support vector machine (SVM) from Eggensperger et al. (2021) are used. The LR design space consists of the regularization strength α and the initial learning rate η0 _0, both within [10−5,100][10^-5,10^0]; the SVM design space consists of the penalty parameter C and the kernel coefficient γ, both within [2−10,210][2^-10,2^10]. All four hyperparameters are optimized in log space. Both benchmarks provide a single fidelity dimension, discretized by the original authors. For LR it is the number of SGD iterations, taking the five values 10,37,111,333,1000\10,37,111,333,1000\, with the subsample fidelity fixed to its default value of 11; for SVM it is the fraction of the training data used, taking the three values 1/9,1/3,1\1/9,1/3,1\. We use the following 8 tasks: 10101, 53, 146818, 146821, 9952, 146822, 31, and 3917. Each task is represented by the meta-features of its underlying OpenML (Vanschoren et al. 2014) dataset, retaining only those meta-features that are available for all tasks and non-constant across them. The resulting 20-dimensional feature vector is then standardized to zero mean and unit variance per dimension across the task set. Each of the 8 tasks serves as the target task in turn, and for each run we uniformly sample 3 of the remaining 7 tasks as the historical tasks, so that no method is given access to the full task set at once. Each historical task contributes 30 observations, uniformly distributed across its fidelity levels, and the target task is seeded with 5 observations at its lowest fidelity. The historical subset and all initial observations are resampled independently for each of the 10 runs. Computational Resources Experiments are implemented in Python using the BoTorch (Balandat et al. 2020) library version 0.16.1. The accompanying code contains all the software requirements and installation instructions. All experiments are executed on a single NVIDIA L40S using float64, as suggested by BoTorch’s official documentation. 1 Appendix C Additional Results Throughout the paper, all methods report the maximizer of the posterior mean of their target fidelity. An alternative convention is to report the best point observed so far, ^t=argmaxi∈tyi x_t= _ x_i _ty_i. This rule is only available to methods that observe the fidelity they report on, i.e., the Next-Best variants; the remaining methods target the out-of-the-loop fidelity, which is never observed in the loop. Figures˜7 and 8 repeat the experiments of Section˜6 with the Next-Best methods using the observed incumbent. All other methods and all experimental settings are unchanged. Figure 8: Real-world benchmarks under the observed-incumbent rule. All other methods are unchanged from Figure˜6. Appendix D Ablation Studies Figure 9: Results of the ablation study of our proposed methods, Multi-task (MT) and Multi-task Next-Best (MT-NB) on Branin, Park and Hartmann synthetic functions. We vary the historical dataset size (top), the number of tasks (middle), and historical tasks’ relevance (bottom). Results are averaged over 10 independent trials. In this section, we evaluate the impact of different components of our method on the Branin, Park, and Hartmann synthetic functions described in Section 6. Unless stated otherwise, we use the same setting as before and perform 10 trials per configuration. Historical dataset size We first analyze the effect of the historical dataset size (first row of Figure˜9). We vary |Hn|∈20,30,40|H_n|∈\20,30,40\ for all previous tasks. For lower-dimensional problems (Branin and Park), increasing the amount of historical data consistently reduces cumulative regret. For Hartmann, however, the results are mixed and do not show a clear trend. Number of tasks Next, we vary the number of available historical tasks, N∈2,3,4N∈\2,3,4\. The task-defining parameters are sampled znz_n uniformly. For Branin and Park, having more historical tasks clearly improves performance. For Hartmann, the effect is less conclusive, likely because task diversity increases with dimensionality. Auxiliary Task Relevance Finally, we examine the effect of task similarity. We sample the incumbent task parameter as zN∼Uniform(0,1)z_N (0,1) and draw the remaining parameters z1,…,zN−1z_1,…,z_N-1 under low, medium, or high relevance. Low relevance is defined by plow(z∣zN)∝|z−zN|2p_low(z z_N) |z-z_N|^2, medium by pmid(z∣zN)=Uniform(0,1)p_mid(z z_N)=Uniform(0,1), and high by phigh(z∣zN)=(zN,0.12)p_high(z z_N)=N(z_N,0.1^2). As shown in the last row of Figure˜9, highly relevant historical tasks substantially improve optimization performance. Moreover, as dimensionality increases, the differences between low, medium, and high relevance become more pronounced.