Constrained Efficient Global Optimization of Expensive Black-box Functions

Donald R. JonesMatthias SchonlauW. Welch

article2022ICML1,867 citations

Introduces the Efficient Global Optimization (EGO) algorithm using Kriging-based stochastic process models and expected improvement to find global optima of expensive black-box functions with minimal evaluations.

Listen

Engineering and industrial design increasingly rely on high-fidelity computer simulations to test alternative designs and eliminate expensive physical prototypes. However, these simulations often take hours or days per run, severely limiting the number of evaluations teams can perform. Conventional global optimization methods typically require hundreds or thousands of evaluations, making them impractical for expensive engineering functions.

The article develops and demonstrates an automated global optimization algorithm called Efficient Global Optimization (EGO), paired with a statistical response surface framework, to locate optimal designs with very few function evaluations while providing intuitive tools for design exploration.

The authors use a stochastic process framework known as the Design and Analysis of Computer Experiments (DACE), also known as kriging. The approach begins by evaluating a modest, space-filling initial sample (typically about 10 points per input dimension). It fits an approximating surface that models both the predicted values and the prediction uncertainty across the design space. Before optimizing, the framework validates model credibility using cross-validation diagnostic tests. The algorithm then guides subsequent search by maximizing an "expected improvement" metric via branch-and-bound optimization, systematically balancing local exploitation (sampling where the model predicts low values) with global exploration (sampling where uncertainty is high).

The article demonstrates five major findings. First, EGO locates near-global optima with exceptional sample efficiency; on standard benchmark problems ranging from two to six dimensions, the algorithm met a 1% stopping criterion in only 28 to 84 total evaluations. Second, the expected improvement criterion provides a credible, objective stopping rule that reflects the marginal potential gain from continued searching. Third, cross-validation diagnostic plots reliably identify poor model fits, allowing users to apply transformations (such as logarithmic scaling) to restore predictive accuracy. Fourth, global sensitivity analysis allows complex models to be simplified; in a 36-variable integrated-circuit case study, just two variables accounted for 66.4% of total output variation, enabling engineers to focus on roughly five key drivers. Fifth, the surrogate surfaces allow rapid multi-objective tradeoff analysis, as demonstrated in an automotive material formulation where engineers uncovered a previously unknown design region that simultaneously improved viscosity and yield stress.

These findings mean engineering teams can dramatically compress development timelines, reduce compute costs, and mitigate risk when optimizing computationally expensive products. Rather than relying on manual trial-and-error, organizations can automate global search with confidence that the algorithm will avoid getting trapped in local optima while minimizing total simulation calls.

To adopt this methodology, engineering organizations should implement the initial space-filling design and validate model fits using standardized cross-validation residuals before optimizing. When the stopping rule slightly understates remaining uncertainty, users should set tighter stopping thresholds (such as 0.1% or requiring criteria satisfaction across consecutive iterations). Practitioners facing high-dimensional problems should run variance decompositions to eliminate non-influential parameters. Promising future opportunities include integrating both low- and high-fidelity simulations into a unified multi-fidelity model and utilizing available gradient data.

Key limitations include computational scaling and numerical sensitivity. The underlying correlation matrices can become ill-conditioned when sample points cluster closely, requiring numerical stabilizing techniques such as singular value decomposition. Furthermore, solving the branch-and-bound subproblem becomes computationally heavy in higher dimensions, requiring limited-memory approximations. Within these identified boundary conditions, the findings provide high confidence that EGO is a rigorous and highly efficient methodology for optimizing expensive engineering systems.

General Motors Research Laboratories National Institute of Statistical Sciences University of Waterloo engineering optimization optimization problems efficient global optimization expensive black-box functions design optimization

arXiv: 2211.00162 JackieXuw/CONFIG

Paper: A Tutorial on Bayesian Optimization, Peter I. Frazier (2018). Provides a comprehensive foundational tutorial on Bayesian optimization, Gaussian process regression, and the expected improvement acquisition criterion foundational to Efficient Global Optimization.
Paper: A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning, Eric Brochu et al. (2010). Explains core Bayesian optimization principles, covariance modeling, and the exploration-exploitation balance via expected improvement on expensive black-box objective functions.
Paper: Practical Bayesian Optimization of Machine Learning Algorithms, Jasper Snoek et al. (2012). Demonstrates practical sequential model-based Bayesian optimization with Gaussian processes and expected improvement acquisition functions for optimizing computationally demanding systems.
Paper: Algorithms for Hyper-Parameter Optimization, James Bergstra et al. (2011). Introduces sequential model-based optimization methodologies using expected improvement criteria to navigate complex, expensive objective spaces.
Paper: On Hyperparameter Optimization of Machine Learning Algorithms: Theory and Practice, Li Yang et al. (2020). Surveys the theoretical and practical foundations of black-box optimization techniques, including surrogate-assisted and Bayesian optimization frameworks.
Paper: Active Learning with Statistical Models, D. Cohn et al. (1996). Establishes variance-reduction principles and statistical active learning criteria that underpin uncertainty-guided exploration in surrogate-based experiment design.
Paper: Optuna: A Next-generation Hyperparameter Optimization Framework, Takuya Akiba et al. (2019). Provides practical context for implementing automated black-box optimization and surrogate sampling strategies across multidimensional parameter spaces.

Paper: Model Discovery Agent: LLM-assisted Bayesian experiment design for data-efficient discovery of mechanistic world models, Kevin Murphy (2026). Extends data-efficient Bayesian experimental design by combining probabilistic inference with large language model agents to discover complex mechanistic world models from minimal physical trials.

Abstract

In many engineering optimization problems, the number of function evaluations is severely limited by time or cost. These problems pose a special challenge to the field of global optimization, since existing methods often require more function evaluations than can be comfortably afforded. One way to address this challenge is to fit response surfaces to data collected by evaluating the objective and constraint functions at a few points. These surfaces can then be used for visualization, tradeoff analysis, and optimization. In this paper, we introduce the reader to a response surface methodology that is especially good at modeling the nonlinear, multimodal functions that often occur in engineering. We then show how these approximating functions can be used to construct an efficient global optimization algorithm with a credible stopping rule. The key to using response surfaces for global optimization lies in balancing the need to exploit the approximating surface (by sampling where it is minimized) with the need to improve the approximation (by sampling where prediction error may be high). Striking this balance requires solving certain auxiliary problems which have previously been considered intractable, but we show how these computational obstacles can be overcome.

Topics

engineering optimization optimization problems visualization trade-off analysis

Keywords

efficient global optimization expensive black-box functions design optimization constraint functions improvement function maximum likelihood estimation convex relaxation multimodal functions high-fidelity models

1 Introduction
2 Notation and Preliminaries
2.1 Problem statement
2.2 Gaussian process surrogates
2.3 Performance metric
3 Algorithm
4 Analysis
5 Infeasibility Declaration
6 Experiments
6.1 Numerical Results over Sampled Instances from Gaussian Process
6.2 Artificial Numerical Instances
6.3 Tuning The P Controller of a Building
7 Conclusion
References
A Examples of kernel functions
B Proof of Corollary
C Proof of Lemma
D Proof of Theorem
E Proof of Theorem
F Proof of Theorem
G Explicit functions for the artificial numerical instances

Knowls

Knowl 1 — Efficient Global Optimization (EGO) Algorithm
algorithm
The Efficient Global Optimization (EGO) algorithm optimizes expensive, multimodal, deterministic black-box functions $f: \mathcal{D} \to \mathbb{R}$ over a bounded hypercube $\mathcal{D} = \prod_{h=1}^k [\ell_h, u_h]$ by balancing local exploitation and global exploration using the Expected Improvement ( $EI$ ) criterion.
Input: Objective function $f$ , design bounds $[\boldsymbol{\ell}, \mathbf{u}] \subset \mathbb{R}^k$ , stopping relative tolerance $\epsilon_{\text{tol}} = 0.01$
Output: Estimated global minimum $\mathbf{x}^*$ and optimal value $f_{\min}$
Select initial space-filling design $\mathbf{X}_0 = \{\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(n_0)}\}$ (with $n_0 \approx 10k$ points using Latin hypercube sampling)
Evaluate initial sample $y^{(i)} = f(\mathbf{x}^{(i)})$ for $i = 1, \dots, n_0$
Set data set $\mathcal{S} \leftarrow \{(\mathbf{x}^{(i)}, y^{(i)})\}_{i=1}^{n_0}$
Fit DACE stochastic process model parameters $(\mu, \sigma^2, \boldsymbol{\theta}, \mathbf{p})$ on $\mathcal{S}$ via maximum likelihood estimation
Perform cross-validation diagnostics on $\mathcal{S}$ ; if standardized residuals violate $[-3, 3]$ or show severe bias, apply a transformation (e.g., $\ln(y)$ or $-1/y$ ) to response values
loop
    Compute current minimum $f_{\min} = \min \{y : (\mathbf{x}, y) \in \mathcal{S}\}$
    Find point $\mathbf{x}^* = \arg\max_{\mathbf{x} \in \mathcal{D}} E[I(\mathbf{x})]$ using branch-and-bound global search
    if $E[I(\mathbf{x}^*)] < \epsilon_{\text{tol}} \cdot |f_{\min}|$ (or $E[I(\mathbf{x}^*)] < \epsilon_{\text{tol}}$ on log scale) then
        break
    end if
    Evaluate function at new point: $y^* = f(\mathbf{x}^*)$
    Augment dataset: $\mathcal{S} \leftarrow \mathcal{S} \cup \{(\mathbf{x}^*, y^*)\}$
    Re-estimate DACE model parameters $(\mu, \sigma^2, \boldsymbol{\theta})$ on $\mathcal{S}$ via maximum likelihood estimation
end loop
return point $\mathbf{x}^*$ corresponding to $f_{\min}$ and value $f_{\min}$
In high dimensions (such as $k=6$ ), the branch-and-bound optimization of $E[I(\mathbf{x})]$ can be accelerated using a limited-memory variant that retains only the best 50 subregion nodes across a maximum of 500 iterations, followed by local optimization.
Knowl 2 — Closed-Form Expected Improvement Criterion and Monotonicity Properties
equation
Let $f_{\min} = \min(y^{(1)}, \dots, y^{(n)})$ be the best objective function value observed so far. At an unobserved point $\mathbf{x}$ , the objective value is modeled as a normal random variable $Y(\mathbf{x}) \sim \mathcal{N}(\hat{y}(\mathbf{x}), s^2(\mathbf{x}))$ , where $\hat{y}(\mathbf{x})$ is the DACE/kriging predictor and $s(\mathbf{x})$ is its standard error of prediction. The improvement random variable is defined as $I(\mathbf{x}) = \max(f_{\min} - Y(\mathbf{x}), 0)$ .

The expected improvement $E[I(\mathbf{x})] \equiv E[\max(f_{\min} - Y(\mathbf{x}), 0)]$ evaluates analytically to:

$E[I(\mathbf{x})] = (f_{\min} - \hat{y}(\mathbf{x})) \Phi\left(\frac{f_{\min} - \hat{y}(\mathbf{x})}{s(\mathbf{x})}\right) + s(\mathbf{x}) \phi\left(\frac{f_{\min} - \hat{y}(\mathbf{x})}{s(\mathbf{x})}\right)$

where $\phi(\cdot)$ and $\Phi(\cdot)$ denote the standard normal probability density function and cumulative distribution function, respectively.

The partial derivatives of $E[I(\mathbf{x})]$ with respect to $\hat{y}$ and $s$ satisfy:

$\frac{\partial E[I]}{\partial \hat{y}} = -\Phi\left(\frac{f_{\min} - \hat{y}}{s}\right) < 0$

$\frac{\partial E[I]}{\partial s} = \phi\left(\frac{f_{\min} - \hat{y}}{s}\right) > 0$

Because $E[I]$ is strictly decreasing in $\hat{y}$ and strictly increasing in $s$ , an upper bound on $E[I(\mathbf{x})]$ over any subregion $B$ can be computed by substituting a lower bound $y^L \le \min_{\mathbf{x} \in B} \hat{y}(\mathbf{x})$ and an upper bound $s^U \ge \max_{\mathbf{x} \in B} s(\mathbf{x})$ into the closed-form expected improvement formula:

$\max_{\mathbf{x} \in B} E[I(\mathbf{x})] \le (f_{\min} - y^L) \Phi\left(\frac{f_{\min} - y^L}{s^U}\right) + s^U \phi\left(\frac{f_{\min} - y^L}{s^U}\right)$
Knowl 3 — DACE Kriging Stochastic Process Model and Predictor
model/method
Let a deterministic function of $k$ variables be evaluated at $n$ sample points $\mathbf{x}^{(1)}, \dots, \mathbf{x}^{(n)} \in \mathbb{R}^k$ with responses $\mathbf{y} = (y(x^{(1)}), \dots, y(x^{(n)}))'$ . The DACE (Design and Analysis of Computer Experiments) framework models the response as:

$y(\mathbf{x}^{(i)}) = \mu + \epsilon(\mathbf{x}^{(i)}) \quad (i = 1, \dots, n)$

where $\mu$ is an unknown constant mean, and $\epsilon(\mathbf{x})$ is a stationary Gaussian stochastic process with mean zero, variance $\sigma^2$ , and parameterized spatial correlation:

$\text{Corr}(\epsilon(\mathbf{x}^{(i)}), \epsilon(\mathbf{x}^{(j)})) = \exp\left( -d(\mathbf{x}^{(i)}, \mathbf{x}^{(j)}) \right) = \exp\left( -\sum_{h=1}^k \theta_h |x_h^{(i)} - x_h^{(j)}|^{p_h} \right)$

with activity parameters $\theta_h \ge 0$ and smoothness parameters $p_h \in [1, 2]$ . Let $\mathbf{R}$ denote the $n \times n$ correlation matrix with entries $R_{ij} = \text{Corr}(\epsilon(\mathbf{x}^{(i)}), \epsilon(\mathbf{x}^{(j)}))$ , and let $\mathbf{1}$ be an $n$ -vector of ones. Given $\boldsymbol{\theta}$ and $\mathbf{p}$ , the maximum likelihood estimates for $\mu$ and $\sigma^2$ are:

$\hat{\mu} = \frac{\mathbf{1}'\mathbf{R}^{-1}\mathbf{y}}{\mathbf{1}'\mathbf{R}^{-1}\mathbf{1}}$

$\hat{\sigma}^2 = \frac{(\mathbf{y} - \mathbf{1}\hat{\mu})'\mathbf{R}^{-1}(\mathbf{y} - \mathbf{1}\hat{\mu})}{n}$

Substituting $\hat{\mu}$ and $\hat{\sigma}^2$ into the likelihood yields the concentrated likelihood function optimized to estimate $\boldsymbol{\theta}$ and $\mathbf{p}$ .

For an unobserved point $\mathbf{x}^*$ , let $\mathbf{r}(\mathbf{x}^*) = (\text{Corr}(\epsilon(\mathbf{x}^*), \epsilon(\mathbf{x}^{(1)})), \dots, \text{Corr}(\epsilon(\mathbf{x}^*), \epsilon(\mathbf{x}^{(n)})))'$ . The Best Linear Unbiased Predictor (BLUP) $\hat{y}(\mathbf{x}^*)$ and its mean squared error $s^2(\mathbf{x}^*)$ are given by:

$\hat{y}(\mathbf{x}^*) = \hat{\mu} + \mathbf{r}'\mathbf{R}^{-1}(\mathbf{y} - \mathbf{1}\hat{\mu})$

$s^2(\mathbf{x}^*) = \sigma^2 \left[ 1 - \mathbf{r}'\mathbf{R}^{-1}\mathbf{r} + \frac{(1 - \mathbf{1}'\mathbf{R}^{-1}\mathbf{r})^2}{\mathbf{1}'\mathbf{R}^{-1}\mathbf{1}} \right]$

The predictor interpolates the data ( $\hat{y}(\mathbf{x}^{(i)}) = y(\mathbf{x}^{(i)})$ ) with zero uncertainty ( $s^2(\mathbf{x}^{(i)}) = 0$ ) at sampled points.
Knowl 4 — Bounding Kriging Mean Squared Error via Convex Relaxation
model/method
To find an upper bound $s^U$ on the kriging standard error $s(\mathbf{x}) = \sqrt{s^2(\mathbf{x})}$ over a rectangular subregion $[\boldsymbol{\ell}, \mathbf{u}] = \prod_{h=1}^k [\ell_h, u_h]$ , the maximization of $s^2(\mathbf{x})$ subject to correlation constraints is reformulated as the minimization of $-s^2(\mathbf{r})$ :

$\min_{\mathbf{r}, \mathbf{x}} -\sigma^2 \left[ 1 - \mathbf{r}'\mathbf{R}^{-1}\mathbf{r} + \frac{(1 - \mathbf{1}'\mathbf{R}^{-1}\mathbf{r})^2}{\mathbf{1}'\mathbf{R}^{-1}\mathbf{1}} \right]$

subject to $\ln(r_i) + \sum_{h=1}^k \theta_h |x_h - x_h^{(i)}|^{p_h} \le 0$ , $-\ln(r_i) - \sum_{h=1}^k \theta_h |x_h - x_h^{(i)}|^{p_h} \le 0$ , $\ell_h \le x_h \le u_h$ , and induced bounds $r_i^L \le r_i \le r_i^U$ obtained via interval analysis.

The Hessian of this objective function with respect to $\mathbf{r}$ is:

$\mathbf{H}_r = 2\sigma^2 \left[ \mathbf{R}^{-1} - \frac{(\mathbf{R}^{-1}\mathbf{1})(\mathbf{R}^{-1}\mathbf{1})'}{\mathbf{1}'\mathbf{R}^{-1}\mathbf{1}} \right]$

Using the $\alpha$ BB convexification technique, if $\lambda_{\min}$ is the minimum eigenvalue of $\mathbf{H}_r$ , adding the term $\alpha \sum_{i=1}^n (r_i - r_i^L)(r_i - r_i^U)$ with $\alpha = \max(0, -\lambda_{\min}/2)$ ensures the modified objective function is convex while underestimating the original objective on $[\mathbf{r}^L, \mathbf{r}^U]$ .

The nonconvex terms $\ln(r_i)$ , $-\ln(r_i)$ , $|x_h - x_h^{(i)}|^{p_h}$ , and $-|x_h - x_h^{(i)}|^{p_h}$ within the constraints are replaced by linear underestimators (chords and tangent lines over the coordinate intervals). The resulting relaxed problem has a convex objective function with linear constraints and simple bounds, which is solved by standard local optimization to produce a rigorous lower bound to the minimization problem (and hence an upper bound on $s^2(\mathbf{x})$ ).
Knowl 5 — Lower Bounding the Kriging Predictor via Nonconvex Separable Relaxation
model/method
When smoothness parameters $p_h = 2$ for all $h = 1, \dots, k$ , the DACE predictor over a box subregion $[\boldsymbol{\ell}, \mathbf{u}] = \prod_{h=1}^k [\ell_h, u_h]$ can be written as:

$\hat{y}(\mathbf{x}) = \hat{\mu} + \sum_{i=1}^n c_i \exp[-z_i(\mathbf{x})]$

where $\mathbf{c} = \mathbf{R}^{-1}(\mathbf{y} - \mathbf{1}\hat{\mu})$ and $z_i(\mathbf{x}) = \sum_{h=1}^k \theta_h (x_h - x_h^{(i)})^2$ . Interval analysis on $[\boldsymbol{\ell}, \mathbf{u}]$ provides lower and upper bounds $[z_i^L, z_i^U]$ for each $z_i(\mathbf{x})$ .

For each term $c_i \exp(-z_i)$ , a univariate linear underestimator $a_i + b_i z_i$ is constructed over $[z_i^L, z_i^U]$ : if $c_i \ge 0$ , the underestimator is the tangent line to $c_i \exp(-z_i)$ at the interval midpoint; if $c_i < 0$ , it is the secant chord connecting the endpoints $(z_i^L, c_i \exp(-z_i^L))$ and $(z_i^U, c_i \exp(-z_i^U))$ .

Substituting these linear underestimators yields a separable lower bound function:

$\hat{y}(\mathbf{x}) \ge \hat{\mu} + \sum_{i=1}^n a_i + \sum_{h=1}^k \theta_h \left( \sum_{i=1}^n b_i (x_h - x_h^{(i)})^2 \right)$

A rigorous lower bound $y^L \le \min_{\mathbf{x} \in [\boldsymbol{\ell}, \mathbf{u}]} \hat{y}(\mathbf{x})$ is obtained by minimizing each univariate quadratic independently over $[\ell_h, u_h]$ :

$y^L = \hat{\mu} + \sum_{i=1}^n a_i + \sum_{h=1}^k \theta_h \min_{\ell_h \le x_h \le u_h} \left[ \sum_{i=1}^n b_i (x_h - x_h^{(i)})^2 \right]$

This one-dimensional minimization is solved in closed form for each coordinate $h$ regardless of whether the quadratic is convex or concave.
Knowl 6 — Cross-Validation Diagnostics and Transformations for DACE Models
model/method
Before using a DACE surrogate to guide optimization, its predictive validity must be verified using leave-one-out cross-validation. For each sampled point $i \in \{1, \dots, n\}$ , the point $(\mathbf{x}^{(i)}, y(\mathbf{x}^{(i)}))$ is held out, and the model predicts the cross-validated mean $\hat{y}_{-i}(\mathbf{x}^{(i)})$ and cross-validated standard error $s_{-i}(\mathbf{x}^{(i)})$ using the remaining $n - 1$ points while keeping estimated correlation parameters fixed.

The standardized cross-validated residual for point $i$ is defined as:

$e_i = \frac{y(\mathbf{x}^{(i)}) - \hat{y}_{-i}(\mathbf{x}^{(i)})}{s_{-i}(\mathbf{x}^{(i)})}$

A DACE model is considered well-calibrated and valid if:
1. All standardized residuals $e_i$ fall approximately within the interval $[-3, +3]$ .
2. Plots of $e_i$ versus $\hat{y}_{-i}(\mathbf{x}^{(i)})$ show no systematic trends or heteroscedasticity.
3. A normal Q-Q plot of ordered $e_i$ values against standard normal quantiles lies along the $45^\circ$ line.
If the model fails these diagnostics (e.g., presence of large standardized residuals $> 4$ or severe non-normality), variance-stabilizing monotonic transformations of the response variable—such as $y \to \ln(y)$ or $y \to -1/y$ (or $-\ln(-y)$ for negative objectives)—should be applied before re-fitting the surrogate.
Knowl 7 — Functional ANOVA Variance Decomposition for DACE Models
model/method
The DACE predictor $\hat{y}(\mathbf{x})$ over a hypercube domain $[0, 1]^k$ can be decomposed into an overall mean, main effects, and interaction effects of increasing order:

$\hat{y}(\mathbf{x}) = a_0 + \sum_{h=1}^k [a_h(x_h) - a_0] + \sum_{h < j} [a_{hj}(x_h, x_j) - a_h(x_h) - a_j(x_j) + a_0] + \dots$

where the terms are defined by integrating over unindexed coordinates:

$a_0 = \int_0^1 \dots \int_0^1 \hat{y}(\mathbf{x}) \, dx_1 \dots dx_k$

$a_h(x_h) = \int_0^1 \dots \int_0^1 \hat{y}(\mathbf{x}) \prod_{j \ne h} dx_j$

$a_{hj}(x_h, x_j) = \int_0^1 \dots \int_0^1 \hat{y}(\mathbf{x}) \prod_{m \ne h, j} dx_m$

Because the summands are mutually orthogonal, the total variance of the predictor decomposes additively:

$\int_{[0, 1]^k} [\hat{y}(\mathbf{x}) - a_0]^2 \, d\mathbf{x} = \sum_{h=1}^k \int_0^1 [a_h(x_h) - a_0]^2 \, dx_h + \sum_{h < j} \int_0^1 \int_0^1 [a_{hj}(x_h, x_j) - a_h(x_h) - a_j(x_j) + a_0]^2 \, dx_h dx_j + \dots$

Dividing individual variance components by total variance yields the percentage of variation attributable to each factor and interaction. Because the product correlation function $\prod_{h=1}^k \exp(-\theta_h |x_h^{(i)} - x_h^{(j)}|^{p_h})$ is separable across dimensions, all multidimensional integrals reduce to products of 1D integrals.
Knowl 8 — Global Convergence Density Conjecture for EGO
theoretical result
Let $f: \mathcal{D} \to \mathbb{R}$ be a continuous objective function defined over a compact hyperbox $\mathcal{D} \subset \mathbb{R}^k$ with finite lower and upper bounds. If the correlation parameters $\theta_h$ in the DACE model are strictly bounded away from zero ( $\theta_h > 0$ for all $h = 1, \dots, k$ ), then the sequence of search points $\{\mathbf{x}^{(1)}, \mathbf{x}^{(2)}, \dots\}$ generated by running the EGO algorithm indefinitely without stopping is conjectured to form a dense subset of the feasible domain $\mathcal{D}$ .

Knowl 9 — Empirical Performance of EGO on Benchmark Global Optimization Functions

data/table

The performance of the EGO algorithm was evaluated on four standard global optimization test problems: Branin ( $k=2$ ), Goldstein–Price ( $k=2$ , modeled on $\ln(y)$ ), Hartman 3 ( $k=3$ ), and Hartman 6 ( $k=6$ , modeled on $-\ln(-y)$ ). For all tests, $p_h = 2$ was fixed in the correlation function.

Test	Evaluations to meet	Actual error	Evaluations required
problem	stopping criterion	when stopped	for 1% accuracy
Branin	28	0.2%	28
Goldstein–Price	32	0.1%	32
Hartman 3	34	1.7%	35
Hartman 6	84	1.9%	121

The algorithm terminates when the expected improvement drops below 1% of the current best value (or $< 0.01$ absolute on the log scale). On Hartman 3 and Hartman 6, the actual error at stopping slightly exceeds the 1% target (1.7% and 1.9%) due to two factors: the prediction standard error ignores parameter estimation uncertainty in $\boldsymbol{\theta}$ and $\sigma^2$ , and the 1-step expected improvement metric understates the potential gain of all remaining evaluations. Requiring the stopping criterion to be met for two consecutive iterations reduced the actual relative error to 0.5% for Hartman 3 and 0.6% for Hartman 6.

Knowl 10 — Multi-Fidelity DACE Surrogate Modeling
model/method
When expensive high-fidelity computer models can be approximated by faster low-fidelity models, DACE surrogate accuracy can be improved without substantial additional computational cost via two approaches:
1. Variable Augmentation: The output of the low-fidelity model evaluated at $\mathbf{x}$ is treated as an additional explanatory input coordinate $x_{k+1} = y_{\text{low}}(\mathbf{x})$ in a standard DACE model fitted to the sparse high-fidelity evaluations $(\mathbf{x}^{(i)}, y_{\text{high}}(\mathbf{x}^{(i)}))$ . Predictions at a new design $\mathbf{x}^*$ are obtained by first evaluating $y_{\text{low}}(\mathbf{x}^*)$ and feeding $(\mathbf{x}^*, y_{\text{low}}(\mathbf{x}^*))$ into the $(k+1)$ -dimensional DACE predictor.
2. Co-Kriging: The low-fidelity and high-fidelity model responses are jointly treated as correlated dependent outputs in a multivariate Gaussian process framework.

Coverage note — Omitted the proprietary automotive multi-objective case study details (which illustrated Pareto tradeoff tracing rather than introducing a new algorithm) and Appendix 2's proposed max-min formulation for MSE bounding (which was explicitly presented as an unworked future concept).

References

1.Alexandrov, N.M., Dennis, Jr., J.E., Lewis, R.M. and Torczon, V. (1998), A trust-region framework for managing the use of approximation models in optimization, Structural Optimization 15: 16–23.
2.Androulakis, I.P., Maranas, C.D. and Floudas, C.A. (1995), $\alpha$ BB: a global optimization method for general constrained nonconvex problems, Journal of Global Optimization 7: 337–363.
3.Aslett, R., Buck, R.J., Duvall, S.G., Sacks, J. and Welch, W.J. (1998), Circuit optimization via sequential computer experiments: design of an output buffer, Applied Statistics 47: 31–48.
4.Betro, B. (1991), Bayesian methods in global optimization, Journal of Global Optimization 1: 1–14.
5.Boender, C.G.E. and Romeijn, H.E. (1995), Stochastic methods, in R. Horst and P.M. Pardalos, eds., Handbook of Global Optimization, pp. 829–869, Kluwer Academic Publishers, Dordrecht/Boston/London.
6.Booker, A.J., Conn, A.R., Dennis, J.E., Frank, P.D., Trossett, M. and Torczon, V. (1995), Global modeling for optimization, Boeing Information and Support Services, Technical Report ISSTECH-95-032.
7.Booker, A.J., Dennis, J.E., Frank, P.D., Serafini, D.B. and Torczon, V. (1997), Optimization using surrogate objectives on a helicopter test example, Boeing Shared Services Group, Technical Report SSGTECH-97-027.
8.Box, G.E.P., Hunter, W.G. and Hunter, J.S. (1978), Statistics for Experimenters, John Wiley, New York.
9.Cox, D.D. and John. S. (1997), SDO: A statistical method for global optimization, in N. Alexandrov and M.Y. Hussaini, eds., Multidisciplinary Design Optimization: State of the Art, pp. 315–329, SIAM, Philadelphia.
10.Cressie. N. (1989), Geostatistics, The American Statistician 43: 197–202. See also the comment on this article by G. Wahba and the reply by N. Cressie (1990), in The American Statistician 44: 255–258.
11.Cressie, N. (1990), The origins of kriging, Mathematical Geology 22: 239–252.
12.Cressie, N. (1993), Statistics for Spatial Data, John Wiley, New York.
13.Currin, C., Mitchell, T., Morris, M. and Ylvisaker, D. (1991), Bayesian prediction of deterministic functions, with applications to the design and analysis of computer experiments, Journal of the American Statistical Association 86: 953–963.
14.Dixon, L.C.W. and Szego, G.P. (1978), The global optimisation problem: an introduction, in L.C.W. Dixon and G.P. Szego (eds.), Towards Global Optimisation, Vol. 2, pp. 1–15. North Holland, Amsterdam.
15.Eby, D., Averill, R.C., Punch III, W.F. and Goodman, E.D. (1998), Evaluation of injection island GA performance on flywheel design optimization, in I.C. Parmee (ed.), Adaptive Computing in Design and Manufacture, Springer Verlag.
16.Elder IV, J.F. (1992), Global $\mathrm{R}^{d}$ optimization when probes are expensive: the GROPE algorithm. Proceedings of the 1992 IEEE International Conference on Systems, Man, and Cybernetics, Vol. 1, pp. 577–582, Chicago.
17.Gao, F., Sacks, J. and Welch, W.J. (1996), Predicting urban ozone levels and trends with semiparametric modeling, Journal of Agricultural, Biological, and Environmental Statistics 1: 404–425.
18.Koehler, J. and Owen, A. (1996), Computer experiments, in S. Ghosh and C.R. Rao (eds.), Handbook of Statistics, 13: Design and Analysis of Experiments, pp. 261–308, Elsevier, Amsterdam.
19.Kushner, H.J. (1964), A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise, Journal of Basic Engineering 86: 97–106.
20.Locatelli, M. (1997), Bayesian algorithms for one-dimensional global optimization, Journal of Global Optimization 10: 57–76.
21.Matheron, G. (1963), Principles of geostatistics, Economic Geology 58: 1246–1266.
22.McKay, M.D., Conover, W.J. and Beckman, R.J. (1979), A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 21: 239–245.
23.Mockus, J. (1994), Application of Bayesian approach to numerical methods of global and stochastic optimization, Journal of Global Optimization 4: 347–365.
24.Mockus, J., Tiesis, V. and Zilinskas, A. (1978), The application of Bayesian methods for seeking the extremum, in L.C.W. Dixon and G.P. Szego (eds.), Towards Global Optimisation, Vol.2, pp. 117–129. North Holland, Amsterdam.
25.Morris, M.D., Mitchell, T.J. and Ylvisaker, D. (1993), Bayesian design and analysis of computer experiments: use of derivatives in surface prediction, Technometrics 35: 243–255.
26.Parzen, E. (1963), A new approach to the synthesis of optimal smoothing and prediction systems, in R. Bellman, (ed.), Mathematical Optimization Techniques, pp. 75–108, University of California Press, Berkeley.
27.Perttunen, C. (1991), A computational geometric approach to feasible region division in constrained global optimization, Proceedings of the 1991 IEEE Conference on Systems, Man, and Cybernetics, Vol. 1, pp. 585–590.
28.Press, W.H., Flannery, B.P., Teukolsky, S.A. and Vetterling, W.T. (1993), Numerical Recipes in FORTRAN, Cambridge University Press.
29.Sacks, J., Welch, W.J., Mitchell, T.J. and Wynn, H.P. (1989), Design and analysis of computer experiments (with discussion), Statistical Science 4: 409–435.
30.Sacks, J. and Ylvisaker, D. (1970), Statistical designs and integral approximation, in R. Pyke (ed.), Proceedings of the Twelfth Biennial Seminar of the Canadian Mathematical Congress, pp. 115–136, Canadian Mathematical Congress, Montreal.
31.Schonlau, M. (1997), Computer experiments and global optimization, Ph.D. Thesis, University of Waterloo, Waterloo, Ontario, Canada.
32.Schonlau, M., Welch, W.J. and Jones, D.R. (1998), Global versus local search in constrained optimization of computer models, to appear in N. Flournoy, W.F. Rosenberger and W.K. Wong (eds.), New Developments and Applications in Experimental Design, Institute of Mathematical Statistics. Also available as Technical Report RR-97-11, Institute for Improvement in Quality and Productivity, University of Waterloo, Waterloo, Ontario, Canada, December 1997.
33.Stuckman, B.E. (1988), A global search method for optimizing nonlinear systems, IEEE Transactions on Systems, Man, and Cybernetics 18: 965–977.
34.Theil, H. (1971), Principles of Econometrics, John Wiley, New York.
35.Ver Hoef, J.M. and Cressie, N. (1993), Multivariable spatial prediction, Mathematical Geology 25: 219–240.
36.Welch, W.J., Buck, R.J., Sacks, J., Wynn, H.P., Mitchell, T.J. and Morris, M.D. (1992), Screening, predicting, and computer experiments, Technometrics 34: 15–25.
37.Zilinskas, A. (1992), A review of statistical models for global optimization, Journal of Global Optimization 2: 145–153.

Citation

MLA

Xu, W., et al. “Constrained Efficient Global Optimization of Expensive Black-box Functions”. arXiv, 2022, http://arxiv.org/abs/2211.00162v4.

APA

Xu, W., Jiang, Y., Svetozarevic, B., & Jones, C. N. (2022). Constrained Efficient Global Optimization of Expensive Black-box Functions. arXiv. http://arxiv.org/abs/2211.00162v4

Chicago

Xu, W., Y. Jiang, B. Svetozarevic, and C. N. Jones. 2022. “Constrained Efficient Global Optimization of Expensive Black-box Functions”. arXiv. http://arxiv.org/abs/2211.00162v4.

Harvard

Xu, W. et al. (2022) “Constrained Efficient Global Optimization of Expensive Black-box Functions”, arXiv [Preprint]. Available at: http://arxiv.org/abs/2211.00162v4.

Vancouver

1. Xu W, Jiang Y, Svetozarevic B, Jones CN (2022) Constrained Efficient Global Optimization of Expensive Black-box Functions. arXiv

BibTeX

@article{xu2022constrained,
  title = {Constrained Efficient Global Optimization of Expensive Black-box Functions},
  author = {Xu, Wenjie and Jiang, Yuning and Svetozarevic, Bratislav and Jones, Colin N.},
  year = {2022},
  journal = {arXiv},
  url = {http://arxiv.org/abs/2211.00162v4},
  eprint = {2211.00162}
}

Metadata:arXiv

Source Code

This paper has an official code repository available. Click below to access the source code.

View Repository

Access the Paper

This paper is available from its original source. Click below to access the PDF.

Open PDF

License: https://creativecommons.org/licenses/by/4.0/

Constrained Efficient Global Optimization of Expensive Black-box Functions

Abstract

Topics

Keywords

Table of Contents

Knowls

Knowl 1 — Efficient Global Optimization (EGO) Algorithm

Knowl 2 — Closed-Form Expected Improvement Criterion and Monotonicity Properties

Knowl 3 — DACE Kriging Stochastic Process Model and Predictor

Knowl 4 — Bounding Kriging Mean Squared Error via Convex Relaxation

Knowl 5 — Lower Bounding the Kriging Predictor via Nonconvex Separable Relaxation

Knowl 6 — Cross-Validation Diagnostics and Transformations for DACE Models

Knowl 7 — Functional ANOVA Variance Decomposition for DACE Models

Knowl 8 — Global Convergence Density Conjecture for EGO