1 Introduction↩︎

Bayesian optimisation (BO) is a principled framework for optimising expensive black-box functions. Such problems are prevalent in science and engineering applications, including materials discovery, drug design and simulation-based engineering. In these domains, where each evaluation may involve expensive physical experiments or high-resolution simulations, sample efficiency becomes critical. Constraints, often black-box functions themselves, are often equally expensive to evaluate, yet essential to ensure feasible designs. In many applications the feasible region is small, discontinuous, or highly non-linear, making even the discovery of a single feasible point challenging, especially in high-dimensional spaces where data is sparse, as commonly seen in crashworthiness design [1].
 
In many engineering applications, however, practitioners may have access to multiple information sources. Examples include coarse simulations [2], simplified physics-based models [3], [4], or analytical approximations [5], that are cheaper to evaluate but potentially biased or noisy. In drug design, this may include inexpensive simulations alongside costly laboratory experiments. While these information sources are often correlated with the target information source (ground truth), the strength and structure of this correlation can vary significantly. Yet, even when information sources are only weakly correlated with the target, they can still provide valuable information for optimisation and help populate the design space, particularly in data-sparse scenarios.
 
Unlike traditional BO methods that rely on a single information source, multi-source BO leverages multiple models in a cost-aware manner, using cheap sources to guide exploration while reserving expensive evaluations to maintain accuracy. These approaches, often referred to as multi-fidelity in literature, were initially modelled through autoregressive processes [6], [7], with later work extending acquisition strategies to this setting [2], [8]–[10]. While much of the literature uses the term multi-fidelity, our work adopts the more general multi-source perspective, which does not assume a strict fidelity hierarchy.
 
Despite recent progress, the integration of multi-source modelling with constrained BO remains under-explored, as most methods address either constraints or multiple sources in isolation. In particular, leveraging low-cost sources for early-stage exploration, especially when feasible points are unknown, has received little attention. To close this gap, we propose a scalable framework for constrained BO with multiple information sources, targeting problems with expensive target evaluations, hard-to-locate feasible regions, and sparse data. We formalise the constrained optimisation problem with multiple sources as \[\label{eq:optimisation95problem} \begin{align} \min_{\mathbf{x} \in \mathcal{X} \subset \mathbb{R}^d} \quad & f^{(L)}(\mathbf{x}) \\ \text{s.t.} \quad & c_i^{(L)}(\mathbf{x}) \leq 0, \quad i = 1, \ldots, g \end{align}\tag{1}\] where \(f^{(L)} :\mathcal{X} \to \mathbb{R}\) and \(c_i^{(L)}:\mathcal{X} \to \mathbb{R} \;\forall i = 1,...,g\) respectively denote the objective and constraints at the target information source \(L\) (most accurate but costly). We assume that the objective and all constraints are always evaluated jointly at a given source. While this assumption may not hold in all domains, it is natural in computational engineering. The design space \(\mathcal{X} \subset \mathbb{R}^d\) may be high-dimensional. To reduce queries to source \(L\), we additionally leverage cheaper information sources \(f^{(\ell)}\), \(c_i^{(\ell)}\) with \(\ell \in \mathcal{S} = \{0, 1, \ldots, L\}\) which provide biased or noisy estimates at lower cost. These auxiliary sources help explore the design space and identify promising regions, while target source evaluations refine solutions near the feasibility boundary or optima. The goal is to efficiently optimise \(f^{(L)}\) subject to the constraints \(c_i^{(L)} \leq 0\) under a limited budget by balancing information gain and query cost across sources. Our work makes the following key contributions:

• A unified framework for constrained BO with multiple sources, combining cheap auxiliary samples and a trust region heuristic for scalability in to high-dimensional input spaces.

• A systematic comparison of scaling multi-source GP models.

• Extensive experiments on real-world and synthetic benchmarks, showing clear improvements over existing constrained BO methods.

2 Related Work↩︎

Early approaches to multi-source or multi-fidelity modelling used hierarchical Gaussian processes (GPs) with nested fidelities [6], [7], later extended to non-nested settings via autoregressive and discrepancy-based models [11], [12]. In BO, [8] proposed a bandit-based fidelity selection method, while [9] introduced the Multi-Information Source Optimisation (MISO) framework, treating each source as a biased observation of a latent function. Further developments include trace-aware knowledge gradients [2], multi-fidelity MES [10], and safeguards against uninformative sources [13]. Additionally, recent pre-prints [14], [15] explore constrained multi-source BO. The former introduces a simple cost-aware heuristic based on expected improvement with feasibility filtering, while the latter performs sequential optimisation over design and fidelity. However, both approaches remain limited in empirical validation.
 
A parallel line of research has explored information-theoretic acquisition functions which aim to maximise expected information gain about the global optimum. Predictive Entropy Search (PES) [16] and Max-value Entropy Search (MES) [17] estimate information gain about the global optimum. Extensions include Fast Information-theoretic BO (FITBO) [18] and General-purpose Information-based BO (GIBBON) [19]. Constrained variants incorporate feasibility modelling via PES with constraints (PESC) [20] and constrained MES (CMES) [21], [22].
 
For scalability to high dimensions with constraints, proposed strategies include trust-region heuristics [23], its extension of feasibility-aware refinements [24], and scaled lengthscale priors [25], [26], often in combination with log Constrained Expected Improvement (logCEI) [27]. High-dimensional problems with thousands of black-box constraints have been addressed in [28], while [29] employ flow-based ensembles to improve GP scalability.
 
Despite this progress, the joint treatment of constraints and multiple sources to scale BO algorithms remains under-explored. Early discovery of feasible regions is especially critical, yet existing methods rarely exploit auxiliary sources for this purpose. We address this gap by proposing a scalable, unified, information-theoretic framework for constrained BO with multiple sources.

3 Bayesian Optimisation with Black-Box Constraints and Multiple Information Sources↩︎

In this section, we first review BO with unknown constraints and its use of Gaussian processes. We then show how this framework can be extended to incorporate multiple information sources.

3.1 Constrained Bayesian Optimisation↩︎

We consider the constrained optimisation problem in Equation 1 , where the feasible set under the target information source \(L\) is defined as: \[\mathcal{X}_f^{(L)} = \{ \mathbf{x} \in \mathcal{X} \mid c^{(L)}_i(\mathbf{x}) \leq 0, \, i = 1, \ldots, g \}.\] Since constraint predictions may vary across sources, feasibility must be defined with respect to the target source \(L\). BO [30], [31] addresses expensive, black-box problems by learning probabilistic surrogate models, typically GPs, to approximate the objective and constraints, allowing for efficient exploration and exploitation of the design space [32]. An acquisition function \(\alpha(\mathbf{x}; \mathcal{D}_n): \mathcal{X} \to \mathbb{R}\) balances this exploration-exploitation trade-off by selecting the next evaluation point \(\mathbf{x}_+ \in \mathcal{X}\). In the constrained setting, feasibility such that \(\mathbf{x}_+ \in \mathcal{X}_f\), must be accounted for by the acquisition function: \[\mathbf{x}_+ = \mathop{\mathrm{arg\,max}}_{\mathbf{x} \in \mathcal{X}_f} \alpha(\mathbf{x}; \mathcal{D}_n).\] Popular constrained acquisition functions have been proposed in [20], [23], [27], [33], [34]. These methods leverage GP-based uncertainty estimates to guide the search towards promising, feasible regions, enabling efficient optimisation despite limited budgets and unknown constraint landscapes.

3.2 Gaussian Process Regression↩︎

GPs provide a flexible, non-parametric prior over functions. A function \(u : \mathcal{X} \to \mathbb{R}\), where \(\mathcal{X} \subset \mathbb{R}^d\), is said to follow a GP prior if for any finite collection of input points \(\{\mathbf{x}_i\}_{i=1}^n \subset \mathcal{X}\). The corresponding vector of function values follows a multivariate normal distribution \(u \sim \mathcal{GP}(\mu, \Sigma)\) such that \[[u(\mathbf{x}_1), \dots, u(\mathbf{x}_n)]^\top \sim \mathcal{N}(\boldsymbol{\mu}, \mathbf{K}),\] where \(\mu : \mathcal{X} \to \mathbb{R}\) is the mean function with \(\mu_i = \mu(\mathbf{x}_i)\), and \(\Sigma : \mathcal{X} \times \mathcal{X} \to \mathbb{R}\) is the covariance function with \(K_{ij} = \Sigma(\mathbf{x}_i, \mathbf{x}_j)\). In practice, the mean is often taken to be zero, while typical covariance functions include the squared exponential or Matérn kernels, parametrised by hyperparameters that are obtained by optimising the marginal likelihood [35].

3.3 Extending Gaussian Processes to Multiple Data Sources↩︎

To extend GPs to the multi-source setting, we model evaluations indexed by \(\ell \in \mathcal{S} = \{0,1,\dots,L\}\), where \(\ell = L\) denotes the most target source and \(\ell = 0\) the cheapest auxiliary source. The joint domain is then the source–augmented input space \(\mathcal{S} \times \mathcal{X}\), with \(\mathcal{X} \subset \mathbb{R}^d\). Following the MISO framework [9], we consider each source specific function \(u^{(\ell)} : \mathcal{S} \times \mathcal{X} \to \mathbb{R}\): \[\label{eq:miso} u^{(\ell)}(\mathbf{x}) = u^{(L)}(\mathbf{x}) + \Delta^{(\ell)}(\mathbf{x}), \quad \text{with } \Delta^{(L)} \equiv 0.\tag{2}\] where \(u^{(L)}: \mathcal{X} \to \mathbb{R}\) is a latent target function and \(\Delta^{(\ell)} : \mathcal{X} \to \mathbb{R}\) a source-dependent discrepancy. We place independent GP priors over both: \[\begin{align} u^{(L)} &\sim \mathcal{GP}(\mu_L, \Sigma_L), \\ \Delta^{(\ell)} &\sim \mathcal{GP}(\mu_\ell, \Sigma_\ell), \quad \forall\, \ell < L. \end{align}\] We typically assume \(\mu_\ell \equiv 0\), whereas \(\Sigma_L\) and \(\Sigma_\ell\) are chosen from standard kernel families and learned jointly via marginal likelihood maximisation. Since the sum of independent GPs is again a GP, the combined model defines a joint GP over the source-augmented input space: \(u : \mathcal{S} \times \mathcal{X} \to \mathbb{R}\) with \(u(\ell, \mathbf{x}) := u^{(\ell)}(\mathbf{x}) \sim \mathcal{GP}(\mu, \Sigma)\) with mean and covariance functions given by: \[\label{eq:miso95mu95sigma} \begin{align} &\mu(\ell, \mathbf{x}) = \mathbb{E}[u^{(\ell)}(\mathbf{x})] = \mathbb{E}[u^{(L)}(\mathbf{x})] + \mathbb{E}[\Delta^{(\ell)}(\mathbf{x})], \\ &\Sigma\big((\ell, \mathbf{x}), (\ell', \mathbf{x}')\big) = \Sigma_L(\mathbf{x}, \mathbf{x}') + \delta(\ell,\ell') \Sigma_\ell(\mathbf{x}, \mathbf{x}'), \end{align}\tag{3}\] where \(\delta(\ell,\ell')\) is the Kronecker delta. This structure captures inter-source correlations through a shared latent process, while allowing source-specific discrepancies. Unlike naïve multi-output GPs [36] or multi-fidelity models with a strict hierarchy [6], the additive characteristic explicitly encodes a source structure, allowing for joint inference over the latent function \(u^{(L)}\), discrepancies \(\Delta^{(\ell)}\) and hyperparameters from observations \(\mathcal{D}_n = \{(\ell_i, \mathbf{x}_i, y_i)\}_{i=1}^n\), while maintaining flexibility.
 
In practice, we fit a separate MISO model for each black-box function, i.e. one for the objective \(f\) and one for each constraint \(c_i \;\forall i = \{1,...,g\}\), all defined over the same source-augmented domain \(\mathcal{S} \times \mathcal{X}\). Figure 1 illustrates how the model can approximate the target source \(f^{(L)}\) from only two target source observations \(\mathcal{D}^{(L)}\) by leveraging six auxiliary observations \(\mathcal{D}^{(\ell)}\).

4 Constrained Max-Value Entropy Search with Multiple Information Sources↩︎

We propose Multi-Source Constrained Max-value Entropy Search (MS-CMES), extending the Max-value Entropy Search [17] principle to settings where both objective and constraints can be queried at multiple information sources. In each iteration we aim to choose the next pair \((\mathbf{x}_{+}, \ell_{+})\) by balancing three factors: (i) expected information gain about the best feasible solution at the target source, (ii) the correlation between auxiliary and target sources, and (iii) the cost of each source, allowing the method to exploit cheap, approximate sources early on. This is achieved by solving \[\label{eq:acqf95opti95main} (\mathbf{x}_{+}, \ell_{+}) = \mathop{\mathrm{arg\,max}}_{\mathbf{x} \in \mathcal{X}, \ell \in \mathcal{S}} \frac{\alpha_n(\mathbf{x}, \ell)}{\lambda(\mathbf{x},\ell)},\tag{4}\] where \(\lambda(\mathbf{x},\ell)\) denotes the cost of evaluating source \(\ell\). In the following we assume that each function \(f^{(\ell)}\) and constraint \(c_i^{(\ell)} \;\forall i \in \{1,...,g\}\) can be observed through any source \(\ell \in \mathcal{S}\).

4.0.0.1 Mutual Information Gain.

Following the Max-value Entropy Search principle, the function \(\alpha_n\) maximises the expected information gain about the constrained optimum \(f^* = \max_{\mathbf{x} \in \mathcal{X}_f^{(L)}} f^{(L)}(\mathbf{x})\) at the target source \(L\): \[\alpha_n(\mathbf{x}, \ell) = \mathbb{I}\left( \mathbf{u}^{(\ell)}(\mathbf{x}); f^* \right).\] with \(\mathbf{u}^{(\ell)}(\mathbf{x}) = [f^{(\ell)}(\mathbf{x}), c_1^{(\ell)}(\mathbf{x}), \dots, c_g^{(\ell)}(\mathbf{x})] \in \mathbb{R}^{g+1}\). In the following we will omit that \(\mathbf{u}^{(\ell)}\) depends on \(\mathbf{x}\). This expression measures how much querying \((\mathbf{x}, \ell)\) reduces uncertainty over \(f^*\), with all quantities modelled via the joint GP framework described in Section 3.3. As \(f^*\) is unknown, we approximate it from the GP posteriors via: \[\label{eq:f95star95main} f^* := \begin{cases} \max\limits_{\mathbf{x} \in \mathcal{X}^{(L)}_f} f^{(L)}(\mathbf{x}) &\text{if } \mathcal{X}^{(L)}_f \neq \emptyset \\ f^{(L)}\left( \mathop{\mathrm{arg\,min}}\limits_{\mathbf{x} \in \mathcal{X}} \sum_j \bar{c}_j^{(L)}\right) & \text{else} \end{cases}\tag{5}\] with \(\bar{c}_j^{(L)} = \max\big(0, c_j^{(L)}(\mathbf{x}) \big)\). This fallback ensures that, in the absence of any feasible point at the target source (i.e. when \(\mathcal{X}^{(L)}_f = \emptyset\)), the algorithm focuses on reducing uncertainty around the most promising infeasible region. Based on this definition of \(f^*\), mutual information can be written as \[\label{eq:mutual95information95main} \begin{align} \mathbb{I}(\mathbf{u}^{(\ell)}; f^*) &= \mathbb{E}_{f^*} \left[ \mathrm{D}_\mathrm{KL}\left( p(\mathbf{u}^{(\ell)} \mid f^*) \,\|\, p(\mathbf{u}^{(\ell)}) \right) \right] \\ &= \mathbb{E}_{f^*} \left[ \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \left[ \log \frac{p(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)})} \right] \right] \end{align}\tag{6}\] However, \(p(\mathbf{u}^{(\ell)} \mid f^*)\) is intractable to compute directly. Therefore, we adopt the variational lower bound on mutual information as derived by [22], adapted to our multi-source context. By introducing a variational distribution \(q(\mathbf{u}^{(\ell)} \mid f^*)\), we can write the information gain as: \[\label{eq:bounded95MI95main} \begin{align} \mathbb{I}(\mathbf{u}^{(\ell)}; f^*) \geq \underbrace{ \mathbb{E}_{f^*} \left[ \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \log \frac{q(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)})} \right] }_{:= \alpha_n(\mathbf{x}, \ell) }. \end{align}\tag{7}\] This lower bound quantifies the expected information gain from evaluating \((\mathbf{x}, \ell)\), i.e. whether the outcome falls within the feasible region defined by the current value of the constrained optimum \(f^*\). The variational distribution \(q(\mathbf{u} \mid f^*)\) is defined as the normalised posterior over the joint outputs, restricted to the feasible set \(\mathcal{F} := (-\infty, f^*] \times (-\infty, 0]^g \subset \mathbb{R}^{g+1}\): \[\label{eq:q95set95main} q(\mathbf{u}^{(\ell)} \mid f^*) = \begin{cases} \displaystyle\frac{p(\mathbf{u}^{(\ell)})}{\mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F})} & \text{if } \mathbf{u}^{(\ell)} \in \mathcal{F} \\ 0 & \text{otherwise} \end{cases}\tag{8}\] Substituting Equation 8 into the variational lower bound in Equation 7 yields the acquisition function \[\label{eq:alpha95main} \begin{align} \alpha_n(\mathbf{x}, \ell) &:= \mathbb{E}_{f^*} \left[ \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \log \frac{q(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)})} \right]\\ &\approx - \frac{1}{K} \sum_{k=1}^K \log \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F}) \end{align}\tag{9}\] where \[\label{eq:prob95in95F} \begin{align} \mathrm{Pr}&(\mathbf{u}^{(\ell)} \in \mathcal{F}) \\ &= \mathrm{Pr}\left( f^{(\ell)}(\mathbf{x}) \leq f^* \right) \prod_{i=1}^G \mathrm{Pr}\left( c_i^{(\ell)}(\mathbf{x}) \leq 0 \right) \end{align}\tag{10}\] The expectation is estimated by drawing \(K\) samples \(\{f^*_k\}_{k=1}^K\) of the constrained optimum using discrete Thompson sampling from the GP posterior at the target source \(L\). A full derivation is provided in Appendix 7.1.

4.0.0.2 Variance correction.

As shown in Equations 9 and 10 , we are interested in \(\mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F})\). However, this introduces a mismatch: while \(\mathcal{F}\) is defined at the target source \(L\), \(\mathbf{u}^{(\ell)}\) is potentially computed on a different source with \(\ell \neq L\). Thus, a candidate may appear feasible under source \(\ell\), yet violating constraints or optimality at source \(L\), since \(\mathcal{X}^{(L)}_f \neq \mathcal{X}^{(\ell)}_f\). The authors in [37] show that \(\mathbf{u}^{(\ell)} \mid f^{(L)}\) is not a truncated Gaussian but rather an extended skew Gaussian distribution, which they approximate in [19], [37] via a variance-correction by quantifying the correlation \(\rho(\mathbf{x},\ell) \in (0,1]\) between source \(\ell\) and target \(L\). This mechanism ensures that weakly informative sources are automatically penalised, preventing spurious feasibility or infeasibility from dominating early-stage acquisition decisions.
 
Since the MISO model assumes conditional independence between fidelities given the latent, without backward conditioning, we retain the truncated Gaussian approximation in Equation 9 , but correct its variance [19] via: \[\label{eq:sigma95tilde95main} \tilde{\sigma}^{(\ell)}(\mathbf{x}) \approx \sigma^{(L)}(\mathbf{x}) \left(1 - \rho^2(\mathbf{x},\ell) \Psi \left( \gamma^{(L)} \right)\right),\tag{11}\] where \(\gamma^{(L)}=\frac{t-\mu^{(L)}(\mathbf{x})}{\sigma^{(L)}(\mathbf{x})}\) and \(\Psi(\gamma)=\frac{\phi(\gamma)}{\Phi(\gamma)}\Big(\gamma+\frac{\phi(\gamma)}{\Phi(\gamma)}\Big)\). We set \(t=f^*\) for the objective and \(t=0\) for the constraints, respectively. Hence, the adjusted feasibility probabilities from Equation 10 can be computed via: \[\label{eq:var95corr95main} \begin{align} \mathrm{Pr}\left( f^{(\ell)}(\mathbf{x}) \leq f^* \right) &\approx \Phi\!\left(\frac{f^*-\mu_f^{(\ell)}(\mathbf{x})}{\tilde{\sigma}_f^{(\ell)}(\mathbf{x})}\right) \\ \mathrm{Pr}\left( c_i^{(\ell)}(\mathbf{\mathbf{x}}) \leq 0 \right) &\approx \Phi\!\left(-\frac{\mu_{c_i}^{(\ell)}(\mathbf{x})}{\tilde{\sigma}_{c_i}^{(\ell)}(\mathbf{x})}\right). \end{align}\tag{12}\] with \(\Phi(\bullet)\) being the cumulative distribution function of a standard Gaussian. We emphasise that the corrected variance \(\tilde{\sigma}^{(\ell)}\) is always computed with respect to the predictive mean \(\mu^{(L)}\) and variance \(\sigma^{(L)}\) of data source \(L\), since the goal is to reduce uncertainty about the target-source optimum \(f^*\).This follows from the fact that the mutual information gained by evaluating at source \(\ell\) can be interpreted as upper bounded: \[\label{eq:I95inequ} \mathbb{I}(\mathbf{u}^{(\ell)}(\mathbf{x});f^*) \leq \mathbb{I}(\mathbf{u}^{(L)}(\mathbf{x});f^*),\tag{13}\] with reduction scaled by \(\rho^2(\mathbf{x},\ell)\), as illustrated in Figure 2. Importantly, the scaling also prevents overconfidence of auxiliary models which may appear more certain due to being trained on a larger data set.

This correction is crucial in early-stage optimisation, where feasible regions are unknown and auxiliary sources may be misleading. By inflating variances at untrusted sources (when \(\rho^2 \approx 0\)), the model can still assign non-negligible feasibility probability as uncertainty is high. Rather than discarding such points, the acquisition function keeps encouraging exploration near the feasibility boundary. Details on how the correlation is computed in case of the MISO model can be found in Appendix 7.2

4.0.0.3 Trust region acquisition optimisation.

To further increase the efficiency and scalability, we restrict acquisition optimisation to a dynamically updated trust region (TR), centred at the best observed point \(\mathbf{x}^\dagger \in \mathcal{D}_n\), restricted to all points where \(\ell=L\), following the idea of [23]: \[\mathbf{x}^\dagger = \mathop{\mathrm{arg\,max}}_{(\mathbf{x}_i, L)\in\mathcal{D}_n} f^{(L)}(\mathbf{x}_i) \quad \text{s.t.} \quad c_j^{(L)}(\mathbf{x}_i)\le 0,\;\forall j ,\] or, if no feasible point is known, as the one with minimal total constraint violation. The TR is defined as a hypercube around \(\mathbf{x}^\dagger\) with side length \(r\), as \[\mathcal{T}(\mathbf{x}^\dagger, r) = \big\{ \mathbf{x}\in\mathcal{X}\,\big|\, \mathbf{x}^\dagger-\tfrac{r}{2}\le \mathbf{x} \le \mathbf{x}^\dagger+\tfrac{r}{2}\},\] clipped to the domain \(\mathcal{X}=[0,1]^d\). Hence, the \(K\) samples \(\{f^*_{k}\}_{k=1}^K\) are generated within \(\mathcal{T}\), and candidate points are then selected by \[\label{eq:acqf95opti95main95tr} (\mathbf{x}_{+}, \ell_{+}) = \mathop{\mathrm{arg\,max}}_{\mathbf{x}\in \mathcal{T}(\mathbf{x}^\dagger, r),\, \ell\in\mathcal{S}} \;\frac{\alpha_n(\mathbf{x},\ell)}{\lambda(\mathbf{x},\ell)}.\tag{14}\] Depending on the progress of the optimisation, \(r\) shrinks or expands. More details can be found in Appendix 7.3. A summary of the acquisition strategy is provided in Algorithm 3, where the optimisation problem in Line 4 is solved using gradient ascent. The complexity is discussed in Appendix 8.
 
Summarising, the key novelty lies in combining a constrained, entropy-based information gain criterion with explicit multi-source modelling and variance correction. Unlike previous single-source constrained BO methods [22] or multi-source methods without constraints [9], MS-CMES provides the first principled acquisition rule for constrained optimisation that automatically down-weights uninformative sources.

5 Numerical Experiments↩︎

We present numerical experiments that first compare multi-source models and their scalability with increasing dimensionality, then benchmark the proposed MS-CMES acquisition strategy against state-of-the-art methods, followed by an ablation study and parameter sensitivity analysis.

5.1 Scalability of Gaussian Processes with Multiple Data Sources↩︎

We compare the performance of the MISO model (Subsection 3.3) against three representative alternatives: the Kennedy–O’Hagan (KOH) model [6], [7], multi-task Gaussian processes (MTGP) [36], and a standard GP trained only on target-source data. As test function we use the Rosenbrock function [38] in dimensions \(d \in \{10,50,100\}\). The target source contains \(n_L = d\) data points, while the auxiliary source contains \(n_\ell = 4d\). Auxiliary signals with no, weak, and strong correlation are constructed following the procedure in Appendix 10. Figure 4 reports the normalised RMSE averaged over 10 random seeds.

When auxiliary and target sources are strongly correlated, all multi-source models achieve comparable accuracy and consistently outperform the single-source GP. This behaviour is expected in settings such as mesh-refined aerodynamic solvers, where lower-cost simulations preserve much of the structural information of their high-cost counterparts.
 
By contrast, strong correlations cannot always be assumed. In scenarios where simulators represent the same physical system but differ substantially, or where simulation data complements experimental data, the outputs are often weakly correlated or structurally divergent. In such cases, models such as KOH and MTGP, which impose a strict hierarchy or assume shared latent structure, show degraded performance. MISO, by explicitly modelling source-specific discrepancies, remains robust across a wide range of inter-source correlations. Notably, even when auxiliary sources are entirely uninformative or uncorrelated, the performances of MISO and MTGP do not degrade much below that of a single-source GP trained only on high-fidelity data. In contrast, KOH, which assumes a strict hierarchy between sources, performs significantly worse in such settings. This robustness makes MISO particularly well-suited for general multi-source optimisation, especially when the relationships between sources are unknown, weak, or heterogeneous.

5.2 Benchmark Tests↩︎

To evaluate performance and scalability, we select five constrained benchmark problems with dimensionalities ranging from \(d = 4-100\). These include the physics-based Pressure Vessel benchmark [39], the Rosenbrock function with two constraints [38], the (Rotated) Rastrigin and the Different Powers function [40]. More information can be found in Appendix 13. We use the approach presented in Appendix 10 to derive a corrupted, weakly correlated signal for each objective and constraint of the respective function. Moreover, we focus on the practically relevant regime of low to moderate evaluation budgets, specifically up to 200 evaluations of the target information source, as common in scenarios where evaluation costs of the target source are immense [41]. The corresponding code can be found at github.com/released/upon/acceptance.

5.2.0.1 MS-CMES Setup.

In our MS-CMES implementation, we adopt a constant cost model, depending only on the data source, with cost weights \(c_L = 1000\) for the target source and \(c_\ell = 1\) for the auxiliary source (more information on the cost function can be found in Appendix 12). This translates to an evaluation budget of \(c \approx 2 \cdot 10^5\). We use \(K = 32\) samples to approximate \(f^*\) and initialise the optimisation with a ratio of \(n_\ell / n_L = 5\) points during the design of experiments.

5.2.0.2 Baseline methods.

We compare our method against several state-of-the-art single- and multi-source approaches for constrained BO. These include SCBO [23], FuRBO [24], a vanilla BO (VBO) baseline with dimensionality-scaled lengthscale priors [25] using LogCEI [27] as the acquisition function, and a random search baseline. We also compare against CMES-IBO, the information-theoretic constrained BO method proposed by [22], which our work extends. In Appendix 11, we introduce a minor modification to this method and refer to the modified version as CMES-IBO+. Additionally, we include Constrained Multi-Fidelity Bayesian Optimisation (CMFBO) from [14] as a representative fidelity-aware baseline. We do not include PESC [42] in our comparison, as it has previously been shown to be outperformed by SCBO [23].

5.2.0.3 Results.

All GP models use Matérn kernels and are trained jointly across sources. For benchmarks with \(d \geq 40\), we use a batch size of \(q = 5\) and \(n_0 = 50\) initial samples, while for the Pressure Vessel benchmark we use \(q = 1\) and \(n_0 = 10\) initial samples. For multi-source methods we use \(5n_0\) auxiliary source samples. The results are depicted in Figure 5, averaged over ten random seeds with mean and variance reported. We emphasise that in multi-source methods only outputs on the target source \(f^{(L)},c_1^{(L)},...,c_g^{(L)}\) determine the optimality and feasibility of a sample. Moreover, when no feasible point has been found, we assign its value to the highest feasible objective discovered so far.
 

On the Pressure Vessel benchmark, all methods perform similarly, though MS-CMES achieves slightly better results than CMES-IBO+. CMFBO struggles to converge and random search fails to discover a feasible point. On the Rastrigin family (standard, rotated) and the Different Powers benchmark, most methods find it difficult to locate the first feasible point. Here, MS-CMES stands out, significantly outperforming alternatives, whereas CMFBO, CMES-IBO+, VBO, and random never find a single feasible sample. SCBO performs somewhat better but is eventually surpassed by FuRBO. On the Rosenbrock benchmark, CMFBO and random again fail to converge, while VBO and CMES-IBO+ achieve promising final values but are clearly outperformed by MS-CMES.
 
Overall, MS-CMES delivers the strongest performance, particularly on the more challenging high-dimensional problems. On the Pressure Vessel benchmark, performance is broadly comparable across methods, though MS-CMES still surpasses CMES-IBO+, highlighting the benefits of incorporating auxiliary information. The advantage of our approach becomes especially evident in higher dimensions: whereas most baselines fail to identify any feasible point even after many iterations, MS-CMES discovers feasible solutions immediately after the initial design. Although the auxiliary signals are only weakly correlated (see Appendix 10), exploiting them proves highly beneficial, enabling MS-CMES to locate feasible regions early and substantially enhance optimisation performance.

5.3 Ablation and Parameter Study↩︎

We further examine the effect of the trust region (TR) heuristic by comparing MS-CMES with and without TR-constrained acquisition optimisation, see Figure 10 in Appendix 14. In the unconstrained variant, \(f^*\) is sampled and Equation 14 is optimised over the full domain \(\mathcal{X}\) rather than within the adaptive TR (in \(\mathcal{T}\)). Both variants perform similarly on Pressure Vessel and Rosenbrock, but the TR heuristic provides clear efficiency gains on the more challenging high-dimensional problems. This indicates that dynamically restricting the search space helps concentrate exploration and prevents drifting into unpromising regions.
 
We also investigate the impact of the number of Monte Carlo samples \(K\) used to estimate the entropy bound. Figure 9 in Appendix 14 shows that while performance improves slightly for larger \(K\), moderate values as \(K=10\) are already sufficient. This underlines the robustness of the method, consistent with the findings of [22].

6 Conclusion↩︎

We introduced MS-CMES, the first information-theoretic framework for constrained Bayesian optimisation with multiple information sources. By explicitly modelling inter-source correlations and evaluation costs, MS-CMES leverages inexpensive, approximate sources for early exploration while reserving costly target evaluations for refinement. Across a diverse set of benchmarks, we demonstrated that MS-CMES consistently outperforms established baselines in both identifying feasible regions and converging towards optimal solutions. These results highlight the promise of multi-source modelling for constrained optimisation in domains ranging from engineering design to drug discovery.

6.0.0.1 Limitations.

Our current formulation assumes that the objective and all constraints are evaluated jointly at a given source. While this assumption holds naturally in many simulation-based settings, it is restrictive in cases where partial observations are available (e.g. objective-only or constraint-only measurements). Extending MS-CMES to accommodate selective output queries would require acquisition strategies capable of balancing heterogeneous information gains. Moreover, the method currently presumes that auxiliary sources exist for every objective and constraint. In problems where this is not the case, the framework could easily be adapted by defaulting to single-source models for missing signals. Finally, performance depends on the quality of the user-specified cost function \(\lambda(\mathbf{x},\ell)\). If costs are poorly estimated, the optimisation may misallocate evaluations, making careful cost modelling crucial in practice.

6.0.0.2 Future Work.

Several promising directions remain. Extending the framework to support partial-output queries would increase flexibility in settings where objectives and constraints can be decoupled. In addition, scalable generative models trained on archival data offer a particularly promising avenue. We believe that by treating such models as auxiliary information sources, organisations could recycle past experiments and simulations, effectively turning historical data into a reusable asset that accelerates and guides future optimisation.

7 Details on MS-CMES↩︎

This appendix complements Section 4 by providing detailed derivations for the variational bound used in MS-CMES, the associated variance correction mechanism that accounts for discrepancies between information sources and the trust region heuristic.

7.1 Variational Bound and Variance Correction in MS-CMES↩︎

We aim to compute the mutual information between a new observation \(\mathbf{u}^{(\ell)}\) and the unknown constrained optimum \(f^*\). This quantity is defined as: \[\label{eq:mutual95information} \begin{align} \mathbb{I}(\mathbf{u}^{(\ell)}; f^*) &= \mathbb{E}_{f^*} \left[ \mathrm{D}_\mathrm{KL}\left( p(\mathbf{u}^{(\ell)} \mid f^*) \,\|\, p(\mathbf{u}^{(\ell)}) \right) \right] \\ &= \mathbb{E}_{f^*} \underbrace{ \left[ \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \left[ \log \frac{p(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)})} \right] \right] }_{(*)} \end{align}\tag{15}\] This follows from the standard mutual information identity \(\mathbb{I}(X; Y) = \mathbb{E}_{Y} \left[ \mathrm{D}_\mathrm{KL}\left( p(X \mid Y) \,\|\, p(X) \right) \right]\) with \(\mathrm{D}_\mathrm{KL}\) being the Kullback-Leibler divergence. However, \(p(\mathbf{u}^{(\ell)} \mid f^*)\) is intractable to compute directly. We continue to follow the approach of [22] who introduce a variational distribution \(q(\mathbf{u}^{(\ell)} \mid f^*)\): \[\label{eq:bounded95term} \begin{align} \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \left[ \log \frac{p(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)})} \right] &= \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \left[ \log \frac{q(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)})} \right] +\mathrm{D}_\mathrm{KL} \left( q(\mathbf{u}^{(\ell)} \mid f^*) \,\|\, p(\mathbf{u}^{(\ell)} \mid f^*) \right) \\ &\geq \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \left[ \log \frac{q(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)})} \right] \end{align}\tag{16}\] Since \(\mathrm{D}_\mathrm{KL}(\cdot \,\|\, \cdot) \geq 0\) we end up at the information source-dependent lower bound. Inserting Equation 16 into Equation 15 , we can write: \[\label{eq:bounded95MI} \begin{align} \mathbb{I}(\mathbf{u}^{(\ell)}; f^*) &\geq \mathbb{E}_{f^*} \left[ \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \log \frac{q(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)}(\mathbf{x}))} \right] \end{align}\tag{17}\] This lower bound quantifies the expected information gain obtained by observing whether the evaluation \((\mathbf{x}, \ell)\) yields a realisation within the feasible region associated with the current value of the constrained optimum \(f^*\). In the unconstrained case, the variational distribution \(q(\mathbf{u} \mid f^*)\) corresponds to a truncated Gaussian. In the constrained setting, we define \(q(\mathbf{u} \mid f^*)\) as the normalised posterior over the joint outputs, restricted to the feasible set \(\mathcal{F} \subset \mathbb{R}^{g+1}\), i.e. \(\mathcal{F} := (-\infty, f^*] \times (-\infty, 0]^g\). We define the feasibility probability \(\mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F})\) at data source \(\ell\) as: \[\label{eq:pr95feas} \begin{align} \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F}) = \mathrm{Pr}\left( f^{(\ell)}(\mathbf{x}) \leq f^* \right) \prod_{i=1}^G \mathrm{Pr}\left( c_i^{(\ell)}(\mathbf{x}) \leq 0 \right) \\ \end{align}\tag{18}\] where \(\Phi(\bullet)\) is the commulative distribution function of a standard Gaussian. Using this, the variational distribution \(q(\mathbf{u}^{(\ell)} \mid f^*)\) is defined as \[\label{eq:q95set} q(\mathbf{u}^{(\ell)} \mid f^*) = \begin{cases} \displaystyle\frac{p(\mathbf{u}^{(\ell)})}{ \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F})} & \text{if } \mathbf{u}^{(\ell)} \in \mathcal{F} \\ 0 & \text{otherwise} \end{cases}\tag{19}\] Substituting Equation 19 into the lower bound in Equation 17 gives: \[\label{eq:alpha} \begin{align} \alpha_n(\mathbf{x}, \ell) &:= \mathbb{E}_{f^*} \left[ \mathbb{E}_{\mathbf{u}^{(\ell)} \mid f^*} \log \frac{q(\mathbf{u}^{(\ell)} \mid f^*)}{p(\mathbf{u}^{(\ell)})} \right]\\ &= \mathbb{E}_{f^*} \left[ \int p(\mathbf{u}^{(\ell)} | f^*) \log \frac{q(\mathbf{u}^{(\ell)} | f^*)}{p(\mathbf{u}^{(\ell)})} \, d\mathbf{u}^{(\ell)} \right]\\ &= \mathbb{E}_{f^*} \left[ \int p(\mathbf{u}^{(\ell)} | f^*) \log \frac{p(\mathbf{u}^{(\ell)})}{ \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F}) p(\mathbf{u}^{(\ell)})} \, d\mathbf{u}^{(\ell)} \right] \\ &= \mathbb{E}_{f^*} \left[ \int p(\mathbf{u}^{(\ell)} | f^*) \log \frac{1}{ \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F})} \, d\mathbf{u}^{(\ell)} \right] \\ &= \mathbb{E}_{f^*} \left[ - \log \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F}) \int p(\mathbf{u}^{(\ell)} | f^*)\, d\mathbf{u}^{(\ell)} \right] \\ &= \mathbb{E}_{f^*} \left[ - \log \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F}) \right] \end{align}\tag{20}\] Finally, the expected value is approximated with \(K\) samples as \[\mathbb{E}_{f^*} \left[ - \log \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F}) \right] \approx - \frac{1}{K} \sum_{k=1}^K \log \mathrm{Pr}(\mathbf{u}^{(\ell)} \in \mathcal{F})\]

7.2 Variance Correction↩︎

Recalling the setup of the multi-source model in Equations 2 and 3 , we can quantify the correlation [35] with: \[\begin{align} \rho(\mathbf{x}, \ell) &= \frac{\mathrm{Cov}\left(u^{(L)}(\mathbf{x}), u^{(\ell)}(\mathbf{x})\right)}{\sigma_L(\mathbf{x}) \sigma_\ell(\mathbf{x})} \\ &= \frac{\sigma_L(\mathbf{x})}{\sigma_L(\mathbf{x}) + \sigma_\Delta(\mathbf{x})} \in (0,1]. \end{align}\] with \(\mathrm{Cov}\!\big(u^{(\ell)}(\mathbf{x}),u^{(L)}(\mathbf{x})\big)=\mathrm{Var}\!\big(u^{(L)}(\mathbf{x})\big)={\sigma_{L}}^2\). From that, it follows that \(\rho(\mathbf{x},L) \equiv 1\). The discrepancy \(\Delta^{(\ell)}\) increases predictive variance at more inaccurate data sources which changes the shape of the conditional distribution.

7.3 Trust-Region Constrained Optimisation of the Acquisition Function↩︎

We constrain the optimisation of the acquisition function to a dynamically updated trust region (TR) in order to stabilise the search. This idea, inspired by [23], focuses optimisation on promising neighbourhoods around the current best observed solution, which is particularly beneficial in high-dimensional settings.
 
At iteration \(n\), let the dataset of all evaluations be \(\mathcal{D}_n = \{(\mathbf{x}_i, \ell_i, f^{(\ell_i)}(\mathbf{x}_i), c^{(\ell_i)}_1(\mathbf{x}_i), \dots, c^{(\ell_i)}_g(\mathbf{x}_i)) \}_{i=1}^n\). We select \(\mathbf{x}^\dagger\) from the subset of points evaluated at the target source \(L\), i.e.\(\{(\mathbf{x}_i, \ell_i) \in \mathcal{D}_n \mid \ell_i = L\}\). If at least one feasible point has been observed, \(\mathbf{x}^\dagger\) is defined as the best feasible point: \[\mathbf{x}^\dagger = \mathop{\mathrm{arg\,max}}_{(\mathbf{x}_i, L) \in \mathcal{D}_n} \; f^{(L)}(\mathbf{x}_i) \quad \text{s.t.} \quad c_j^{(L)}(\mathbf{x}_i) \leq 0 \;\forall j .\] If no feasible point exists, \(\mathbf{x}^*\) is chosen as the infeasible point with the smallest total constraint violation: \[\mathbf{x}^* = \mathop{\mathrm{arg\,min}}_{(\mathbf{x}_i, \ell_i) \in \mathcal{D}_n,\;\ell_i = L} \; \sum_{j=1}^g \max \!\big(0, c_j^{(L)}(\mathbf{x}_i)\big).\] The trust region is defined as a hypercube centred at \(\mathbf{x}^*\) with side length \(r\): \[\mathcal{T}(\mathbf{x}^\dagger, r) = \Big\{ \mathbf{x} \in \mathcal{X} \,\big|\, x^\dagger_z - \tfrac{r}{2} \leq x_z \leq x^\dagger_z + \tfrac{r}{2}, \; z=1,\dots,d \Big\},\] clipped to the global domain \(\mathcal{X}=[0,1]^d\).
 
The side length \(r\) is adapted dynamically. Let \(s_n\) and \(f_n\) denote counters of consecutive successes and failures, respectively. A success occurs if the new evaluation improves upon the best feasible objective, or reduces total constraint violation relative to \(\mathbf{x}^*\). Otherwise, a failure is recorded. When \(s_n\) reaches a tolerance \(s_{\max}\), the region is expanded as \(r \gets \min(2r, r_{\max})\) and \(s_n\) is reset. Conversely, if \(f_n\) reaches \(f_{\max}\), the region is shrunk as \(r \gets \tfrac{1}{2} r\) and \(f_n\) reset. If \(r\) falls below \(r_{\min}\), the TR is restarted around the current best point.
 
At each iteration, new candidates are chosen by optimising the cost-scaled acquisition function within the TR: \[(\mathbf{x}_{+}, \ell_{+}) = \mathop{\mathrm{arg\,max}}_{\mathbf{x} \in \mathcal{T}(\mathbf{x}^\dagger, r), \, \ell \in \mathcal{S}} \; \frac{\alpha_n(\mathbf{x}, \ell)}{\lambda(\mathbf{x},\ell)} .\] This strategy prevents the optimiser from wasting evaluations in unpromising regions (particularly when auxiliary sources are misleading) and provides a principled balance between local refinement and global exploration through adaptive expansion and contraction of the trust region.

8 Computational Complexity↩︎

We briefly analyse the computational complexity of the proposed method, separating the cost of model training from that of acquisition optimisation.

8.0.0.1 Multi-Source Gaussian Process Model

Let \(n = \sum_{\ell \in \mathcal{S}} n_{\ell}\) denote the total number of observations across sources. The training and inference requires a Cholesky factorisation of the \(n \times n\) covariance matrix, leading to the classical complexities for time and memory, \(\mathcal{O}(n^3)\) and \(\mathcal{O}(n^2)\). Optimisation of the hyperparameters requires repeated gradient evaluations, hence scales identically. Due to the additive structure, additional kernel parameters are introduced due to the discrepancy GP \(\Delta^{(\ell)}\), compared to a standard, single-source GP. For large data sets \(n\), inducing point methods such as Structured Kernel Interpolation (SKI) or Sparse GPs (SGPR) could be applicable.

8.0.0.2 Acquisition function (MS-CMES)

The cost of evaluating the acquisition function decomposes as follows: We draw \(K\) samples from the posterior of the constrained optimum on a candidate set of size \(n_c\) to obtain \(f^*\). This costs \(\mathcal{O}(K n_c)\) plus the cost of GP posterior evaluation on \(n_c\) candidates. Evaluating \(\alpha(\mathbf{x},\ell)\), for each candidate \(\mathbf{x}\), computing feasibility probabilities involves Gaussian CDFs of the posterior mean/variance. With \(g\) constraints, a batch of \(b\) candidates and \(|\mathcal{S}|\) number of sources, the complexity yields \(\mathcal{O}(b K (g+1)|\mathcal{S}|)\). In total this leads to \(\mathcal{O}(K n_c + b K (g+1)|\mathcal{S}|)\).

8.0.0.3 Computational resources in this work.

All experiments were carried out on a compute cluster equipped with Intel Xeon Gold 5218 CPUs (2.3 GHz). All jobs were executed on single CPU nodes without GPU acceleration. Runtime per experiment varied with problem dimensionality and evaluation budget but remained within a few hours for the largest benchmarks, showing that the method does not need extensive compute.

9 Details on Implementation of Models↩︎

9.1 MISO Model↩︎

We implement a variant of the Multi-Information Source Optimisation (MISO) model based on the structure proposed by [9]. This subsection accompanies the mathematical description in Section 3.3. Our model assumes a hierarchical design of experiments, wherein each target source evaluation is accompanied by evaluations from all auxiliary sources at the same input location. This nested design is motivated by the assumption that auxiliary sources are inexpensive to query, allowing evaluations to be collected “for free” whenever a target source point is selected by the acquisition function. The model decomposes the objective into a shared latent function and source-specific discrepancies. Both the base kernel and the discrepancy kernels are instantiated as Matérn-5/2 kernels. Discrepancy kernels are masked to activate only on data from their associated source level. These components are combined additively in the covariance module. To capture the magnitude of inter-fidelity variation, we place a log-normal prior on the outputscale of each discrepancy kernel. The prior’s mean is set to the empirical mean squared discrepancy between target and auxiliary source observations across the training data, extracted from the nested design. This reflects an informed prior belief on the relative strength of the discrepancy signal. Finally, the model exposes a utility function that returns the local squared correlation \(\rho^2(x)\) between each auxiliary source and the target source, based on the kernel structure. We implemented the model with the help of GPyTorch [43] and use it within BoTorch [44] for optimisation.

9.2 KOH Model↩︎

The Kennedy-O’Hagan (KOH) model [6] is a seminal framework for multi-fidelity/ multi-source GP modelling, originally proposed for calibrating computer models using experimental data. It has since become a cornerstone in multi-source BO and surrogate modelling, particularly when relating a cheap but approximate simulator to a more expensive high-fidelity function.
 
In its canonical form, the KOH model assumes a hierarchical structure that expresses the high-fidelity response \(f_{L}(\mathbf{x})\) as a scaled version of the low-fidelity model \(f_{\ell}(\mathbf{x})\), plus a discrepancy term: \[f_{L}(\mathbf{x}) = \rho f_{\ell}(\mathbf{x}) + \Delta(\mathbf{x}),\] where \(\rho \in \mathbb{R}\) is a scalar calibration or scaling parameter that is learned alongside the other hyperparameters defining the model, \(f_{\ell}(\mathbf{x}) \sim \mathcal{GP}(0, k_{\ell}(\mathbf{x}, \mathbf{x}’))\) is a GP representing the low-fidelity simulator, \(\Delta(\mathbf{x}) \sim \mathcal{GP}(0, k_{\Delta}(\mathbf{x}, \mathbf{x}’))\) is a discrepancy GP capturing the systematic mismatch between fidelities. Under this model, the joint prior over \(f_{\ell}(\cdot)\) and \(f_{L}(\cdot)\) is analytically tractable and Gaussian, and the posterior predictions retain closed-form expressions under Gaussian noise. This hierarchical design enforces that high-fidelity predictions are informed by low-fidelity evaluations through the shared term \(f_{\ell}(\cdot)\), while allowing the discrepancy GP to correct for systematic errors. We implemented this model with the help of GPyTorch [43], using the Matérn-5/2 kernel.

9.3 MTGP Model↩︎

The MTGP model proposed by [36] offers a principled approach for jointly modelling multiple related outputs (tasks) using a shared GP framework. In the context of multi-source modelling, tasks correspond to different sources, providing a flexible and non-hierarchical alternative to models like KOH.
 
Let \(f_\ell(\mathbf{x})\) denote the output of task \(\ell \in \{1, \dots, L\}\) at input \(\mathbf{x} \in \mathcal{X}\). The MTGP assumes that the joint function \(f(\mathbf{x}, \ell)\) is drawn from a GP: \[f(\mathbf{x}, \ell) \sim \mathcal{GP}(0, k((\mathbf{x}, \ell), (\mathbf{x}’, \ell’))),\] with a covariance function decomposed as: \[k((\mathbf{x}, \ell), (\mathbf{x}’, \ell’)) = k_{x}(\mathbf{x}, \mathbf{x}’) \cdot k_{\ell}(\ell, \ell’).\] Here \(k_{x}\) is a standard kernel (e.g. RBF or Matérn) defined over the input space, \(k_{\ell}\) is a positive semi-definite covariance matrix \(\mathbf{K}_{\ell} \in \mathbb{R}^{L \times L}\) that encodes task relationships.
 
This structure allows the MTGP to capture correlations across tasks (sources), enabling knowledge transfer from low-fidelity data to improve predictions at high fidelity. Importantly, unlike the KOH model, MTGPs do not impose any explicit hierarchy or discrepancy structure. Instead, task dependencies are entirely captured through the learned task covariance matrix. In this work, we make use of GPyTorch’s [43] Multi-Task Gaussian Process implementation [36], again using the Matérn-5/2 kernel.

10 Extending Benchmarks to Multi-Fidelities↩︎

Let \(u^{(L)}: \mathbb{R}^d \to \mathbb{R}\) denote a target source objective or constraint function. We define a auxiliary source \(u^{(\ell)}\) by applying an input-dependent oscillatory distortion to the function \(u^{(L)}\). The distortion is scaled relative to the estimated magnitude of each function. Let \(S_f > 0\) denote an empirical estimate of the objective’s scale: \[S := \mathbb{E}_{x \sim \mathcal{U}([0,1]^d)} \left[ |u^{(L)}(x)| \right]\] The oscillatory signal is defined as: \[s(x) := \sin\left(\frac{2\pi}{d} \sum_{j=1}^{d} x_j \right)\] For the weak correlation level, set the relative distortion factor to \(\rho := 1\), for a strong correlation \(\rho := 0.1\). We define the auxiliary source \(\tilde{u}^{(\ell)}\) as: \[u^{(\ell)}(x) := u^{(L)}(x) \cdot \left( 1 + \frac{\rho S}{\max\left( | u^{(L)}(x)|, \varepsilon \right)} \cdot s(x) \right) + \xi\] with \(\xi \sim \mathcal{N}(0,\sigma_\xi^2)\).

11 Baseline Methods: Experiment Setup↩︎

11.0.0.1 Multi-Source Constrained Max-value Entropy Search

In MS-CMES we use the MISO model as discussed in Appendix 9.1 and employ discrete Thompson sampling to sample \(f^*\). Moreover, we use optimize_acqf_mixed in BoTorch [44], using 3 restarts and 200 raw samples to maximise the acquisition function.

11.0.0.2 Vanilla Bayesian Optimisation

Here, we use the dimensionality-scaled lengthscale prior from [25] and the logarithmic version of Constrained Expected Improvements (logCEI), proposed in [27]. We use 64 MC samples and enable the sample_around_best option, embedded in BoTorch’s acquisition function optimiser. More information can be found in [26].

11.0.0.3 Scalable Constrained Bayesian Optimisation

We use the same parameters as in [23]. The method is implemented in BoTorch and can be found here: https://botorch.org/docs/tutorials/scalable_constrained_bo/

11.0.0.4 Feasibility-Driven Trust Region Bayesian Optimisation

FuRBO is directly built upon SCBO and differs only w.r.t. the trust region heuristic. Here again, we employ the same user-defined parameters as presented in [24]. The corresponding code was taken from https://anonymous.4open.science/r/FuRBO.

11.0.0.5 Constrained Max-value Entropy Search

Additionally, we compare against the original CMES-IBO methods, proposed in [22]. While the authors propose to define \(f^*\) as \[f^* := \begin{cases} \max\limits_{\mathbf{x} \in \mathcal{X}_f} f(\mathbf{x}) &\text{if } \mathcal{X}_f \neq \emptyset \\ -\infty & \text{else}, \end{cases}\] we found that our definition from Equation 5 yields better results as depicted in Figure 7 where we plot the mean and standard deviation over ten runs. For all benchmarks (this Section as well as Section 5.2) we use \(K=32\). Similar to before, we employ discrete Thompson sampling to sample \(f^*\) and use optimize_acqf in BoTorch [44], using 3 restarts and 200 raw samples.

11.0.0.6 Constrained Multi-Fidelity Bayesian Optimisation

We implemented our own version of the acquisition function in [14] in combination with the introduced MISO model. We use gradient-ascent to maximise the acquisition function. However, to obtain a similar ratio of auxiliary and target source evaluations we noted that we need to set the cost to \(c_L = 0\). Like before, we employ optimize_acqf in BoTorch [44], using 3 restarts and 200 raw samples.

12 Details on Cost Function↩︎

As the cost function needs to be set up such that it reflects the magnitude of the acquisition function scale, we found that in MS-CMES \(\lambda(\mathbf{x},\ell) = 1 + \frac{\ell}{10^{5}} c_\ell\) balances the number of target and auxiliary source evaluations well. However, for real world problems, the choice of the cost function can be problem specific and needs to reflect the sources it is trying to balance. We emphasise that the choice of cost function is arbitrary as it solely scales the utility value of the acquisition function.
 
To showcase when a target source evaluation is queried, we depict in Figure 8 the accumulated costs over the accumulated evaluations. While auxiliary sources provide cost-effective guidance in early iterations, our acquisition function also continues to select the target source throughout the optimisation.

13 Details on Benchmark Problems↩︎

We evaluate our method on several well-known constrained test problems to stress different aspects of performance:

• The Pressure Vessel benchmark was proposed by [39] with dimensionality \(d=4\) and number of constraints \(g=4\). The aim is to minimise a total cost of designing the pressure vessel, including the shell and head thicknesses, as well as the inner radius and length of the cylindrical section including some bounds.

• The benchmarks Different Powers, Rastrigin and Rotated Rastrigin (all \(d=40\) and \(g=9\)) were drawn from the BBOB-constrained suite [40], These functions are constructed by applying non-linear transformations (e.g.rotations, asymmetric scaling, oscillatory distortions) to the base functions, then overlaying constraint functions so that the feasible region becomes non-trivial. For instance, the Rastrigin and rotated Rastrigin variants combine the highly multimodal base with global rotations, while Different Powers applies coordinate-wise power scaling to create differing levels of smoothness along axes. This set covers both multimodal, non-separable and ill-conditioned landscapes with constraints, making it a rigorous benchmark for constrained multi-source BO.

• Additionally, we test our method on the Rosenbrock objective [38] augmented with two nonlinear constraints (\(d=100\) and \(g=2\)) defined in [23], which challenges the model in high dimensions and with narrow feasible channels.

14 Ablation and Parameter Study↩︎

References↩︎