16 Surrogate Optimization
The previous chapter explained how to use a probabilistic surrogate model, inparticular a Gaussian process, to infer probability distributions over the true
objective function. These distributions can be used to guide an optimization
process toward better design points.1 This chapter outlines several common
techniques for choosing which design point to evaluate next. The techniques
we discuss here greedily optimize various metrics.2 We will also discuss how
surrogate models can be used to optimize an objective measure in a safe manner.

16.1

Prediction-Based Exploration

In prediction-based exploration, we select the minimizer of the surrogate function.
An example of this approach is the quadratic fit search that we discussed earlier
in section 3.5. With quadratic fit search, we use a quadratic surrogate model to fit
the last three bracketing points and then select the point at the minimum of thequadratic function.
If we use a Gaussian process surrogate model, prediction-based optimization
has us select the minimizer of the mean function
x(m+1) = arg min µ̂(x)x∈X(16.1)

where µ̂(x) is the predicted mean of a Gaussian process at a design point x basedon the previous m design points. The process is illustrated in figure 16.1.
Prediction-based optimization does not take uncertainty into account, and new
samples can be generated very close to existing samples. Sampling at locations
where we are already confident in the objective value is a waste of function
evaluations.

A. Forrester, A. Sobester, and A.
Keane, Engineering Design via Surrogate Modelling: A Practical Guide.
Wiley, 2008.
1

An alternative to greedy optimization is to frame the problem
as a partially observable Markov decision process and plan ahead somenumber of steps as outlined by M.
Toussaint, “The Bayesian SearchGame,” in Theory and Principled
Methods for the Design of Metaheuristics, Y. Borenstein and A. Moraglio,
eds. Springer, 2014, pp. 129–144.
See also R. Lam, K. Willcox, andD. H. Wolpert, “Bayesian Optimization with a Finite Budget: An Approximate Dynamic Programming
Approach,” in Advances in Neural Information Processing Systems(NIPS), 2016.
2

c ha p te r 16 . su rrog a te op ti mi z a ti on

true
predicted
fit points
next sample

y

292

x

x

x

x
Figure 16.1. Prediction-based optimization selects the point that minimizes the mean of the objectivefunction.

16.2 Error-Based Exploration
Error-based exploration seeks to increase confidence in the true function. A Gaussian
process can tell us both the mean and standard deviation at every point. A large
standard deviation indicates low confidence, so error-based exploration samples
at design points with maximum uncertainty.
The next sample point is:(16.2)x (m+1) = arg max σ̂(x)x∈X

where σ̂ (x) is the standard deviation of a Gaussian process at a design point x
based on the previous m design points. The process is illustrated in figure 16.2.

yyconfidence regiontrue
predicted
fit points
sampled

x

x

x

x

Gaussian processes are often defined over all of R n . Optimization problems

with unbounded feasible sets will always have high uncertainty far away fromsampled points, making it impossible to become confident in the true underlying
function over the entire domain. Error-based exploration must thus be constrained
to a closed region.

Figure 16.2. Error-based exploration selects a point with maximal
uncertainty.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

16 . 3. l ow e r c on fi d e n c e bou n d e x p l ora ti on 293

16.3 Lower Confidence Bound Exploration
While error-based exploration reduces the uncertainty in the objective function
overall, its samples are often in regions that are unlikely to contain a global minimum. Lower confidence bound exploration trades off between greedy minimization
employed by prediction-based optimization and uncertainty reduction employedby error-based exploration. The next sample minimizes the lower confidence boundof the objective function
LB(x) = µ̂(x) − ασ̂(x)(16.3)

where α ≥ 0 is a constant that controls the trade-off between exploration and exploitation. Exploration involves minimizing uncertainty, and exploitation involves
minimizing the predicted mean. We have prediction-based optimization with
α = 0, and we have error-based exploration as α approaches ∞. The process is
illustrated in figure 16.3.

yyconfidence regiontruelower boundpredicted
fit points
sampled

x

x

x

x

16.4 Probability of Improvement Exploration
We can sometimes obtain faster convergence by selecting the design point thatmaximizes the chance that the new point will be better than the samples we haveseen so far. The improvement for a function sampled at x producing y = f (x) isymin − y if y < ymin
(16.4)
I (y) =
0otherwise

Figure 16.3. Lower confidencebound exploration trades off between minimizing uncertainty andminimizing the predicted function.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

294 c ha p te r 16 . su rrog a te op ti mi z a ti on

where ymin is the minimum value sampled so far.
The probability of improvement at points where σ̂ > 0 is
Z ymin(16.5)N (y | µ̂, σ̂2 )dyymin − µ̂=Φσ̂P(y < ymin ) =−∞(16.6)

where Φ is the standard normal cumulative distribution function (see appendix C.7).
This calculation (algorithm 16.1) is shown in figure 16.4. Figure 16.5 illustrates
this process. When σ̂ = 0, which occurs at points where we have noiseless measurements, the probability of improvement is zero.
true
predicted
fit points

Figure 16.4. The probability of improvement is the probability that
evaluating a particular point will
yield a better result than the best
so far. This figure shows the probability density function predicted
at a query point, with the shaded
region below ymin corresponding
to the probability of improvement.

prob_of_improvement(y_min, μ, σ) = cdf(Normal(μ, σ), y_min)

Algorithm 16.1. Computing theprobability of improvement for agiven best y value y_min, mean μ,
and variance ν.

ybest so far

probability of improvement
query pointx

16.5

Expected Improvement Exploration

Optimization is concerned with finding the minimum of the objective function.
While maximizing the probability of improvement will tend to decrease theobjective function over time, it does not improve very much with each iteration.
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

16 . 5 . e x p e c te d i mp rove me n t e x p l ora ti on 295

yconfidence regiontrue
predicted
fit points
sampledP[ f (x) < ymin ]

1

0

x

x

x

x

We can focus our exploration of points that maximize our expected improvement
over the current best function value.
Through a substitutionz=y − µ̂
σ̂

′
ymin
=

ymin − µ̂σ̂Figure 16.5. Maximizing the probability of improvement selects samples most likely to produce lower
objective point values.

(16.7)

we can write the improvement in equation (16.4) asσ̂(y′ − z) if z < y′ and σ̂ > 0minmin
I (y) =
0otherwise(16.8)

where µ̂ and σ̂ are the predicted mean and standard deviation at the sample pointx.
We can calculate the expected improvement using the distribution predictedby the Gaussian process:E [ I (y)] = σ̂

Z y′

min

−∞

(16.9)

′ymin− z N (z | 0, 1) dz



Z y′
Z y′
min
min
′
= σ̂ ymin
N (z | 0, 1) dz −z N (z | 0, 1) dz
−∞
−∞



′
′
′
P(z ≤ ymin
) + N (ymin| 0, 1) − N (−∞ | 0, 1)
= σ̂ymin
|
{z
}
=0

= (ymin − µ̂) P(y ≤ ymin ) + σ̂2 N (ymin | µ̂, σ̂2 )(16.10)(16.11)(16.12)

Figure 16.6 illustrates this process using algorithm 16.2.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

296 c ha p te r 16 . su rrog a te op ti mi z a ti on

Algorithm 16.2. Computing the expected improvement for a given
best y value y_min, mean μ, andstandard deviation σ.

function expected_improvement(y_min, μ, σ)p_imp = prob_of_improvement(y_min, μ, σ)p_ymin = pdf(Normal(μ, σ), y_min)
return (y_min - μ)*p_imp + σ^2*p_ymin
endyconfidence regiontrue
predicted
fit points
sampledE [ I (y)]

1

x

x

x

x

16.6 Safe Optimization
In some contexts, it may be costly to evaluate points that are deemed unsafe,
which may correspond to low performing or infeasible points. Problems such as
the in-flight tuning of the controller of a drone or safe movie recommendations
require safe exploration—searching for an optimal design point while carefully
avoiding sampling an unsafe design.
This section outlines the SafeOpt algorithm,3 which addresses a class of safe
exploration problems. We sample a series of design points x(1) , . . . , x(m) in pursuitof a minimum but without f (x(i) ) exceeding a critical safety threshold ymax . In
addition, we receive only noisy measurements of the objective function, where thenoise is zero-mean with variance ν. Such an objective function and its associatedsafe regions are shown in figure 16.7.
The SafeOpt algorithm uses Gaussian process surrogate models for prediction.
At each iteration, we fit a Gaussian process to the noisy samples from f . After the

Figure 16.6. Maximizing the expected improvement selects samples which are likely to improve thelower bound by as much as possible.
ymax

y

0

x∗
x
Figure 16.7. SafeOpt solves safe exploration problems that minimize
f while remaining within safe regions defined by maximum objective function values.
3
Y. Sui, A. Gotovos, J. Burdick,and A. Krause, “Safe Exploration
for Optimization with GaussianProcesses,” in International Conference on Machine Learning (ICML),
vol. 37, 2015.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

16 . 6 . sa fe op ti mi z a ti on 297

ith sample, SafeOpt calculates the upper and lower confidence bounds:q(16.13)ui (x) = µ̂i−1 (x) + βν̂i−1 (x)qℓi (x) = µ̂i−1 (x) − βν̂i−1 (x)(16.14)

where larger values of β yield wider confidence regions. Such bounds are shownin figure 16.8.
confidence regiontrue
predictedℓ( x )u( x )
fit points

Figure 16.8. An illustration of functions based on the predictions of a
Gaussian process used by SafeOpt.

x

The Gaussian process predicts a distribution over f (x) for any design point.
Being Gaussian, these predictions can provide only a probabilistic guarantee ofsafety up to an arbitrary factor:4!
ymax − µ̂(x)pP( f (x) ≤ ymax ) = Φ≥ Psafe(16.15)ν̂(x)The predicted safe region S consists of the design points that provide a probabilityof safety greater than the required level Psafe , as illustrated in figure 16.9. The safe
region can also be defined in terms of Lipschitz upper bounds constructed fromupper bounds evaluated at previously sampled points.
SafeOpt chooses a safe sample point that balances the desire to localize areachable minimizer of f and to expand the safe region. The set of potential
minimizers of f is denoted M (figure 16.10), and the set of points that will
potentially lead to the expansion of the safe regions is denoted E (figure 16.11).
To trade off exploration and exploitation, we choose the design point x with thelargest predictive variance among both sets M and E .5
The set of potential minimizers consists of the safe points whose lower confidence bound is lower than the lowest upper bound:′
Mi = x ∈ Si | ℓi (x) ≤ min ui (x )(16.16)x′ ∈Si

Note the similarity to the probability of improvement.

4

5
For a variation of this algorithm,
see F. Berkenkamp, A. P. Schoellig,and A. Krause, “Safe Controller
Optimization for Quadrotors with
Gaussian Processes,” in IEEE International Conference on Robotics andAutomation (ICRA), 2016.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

y298 c ha p te r 16 . su rrog a te op ti mi z a ti on

confidence interval
objective function ffit points
safety threshold
estimated safe region S

Figure 16.9. The safety regions(green) predicted by a Gaussian
process.

estimated safe region S

Figure 16.10. The potential minimizers are the safe points whose
lower bounds are lower than thebest, safe upper bound.

safety probability

x
1

Psafe

0

x

potential minimizers Myconfidence region
objective function
fit points
safety threshold
best upper boundx© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

16 . 6 . sa fe op ti mi z a ti on 299

At step i, the set of potential expanders Ei consists of the safe points that, if
added to the Gaussian process, optimistically assuming the lower bound, produce
a posterior distribution with a larger safe set. The potential expanders naturally
lie near the boundary of the safe region.
estimated safe region S
expected expanders E

Figure 16.11. The set of potential
expanders.

yconfidence region
objective function
fit points
safety thresholdx

Given an initial safe point6 x(1) , SafeOpt chooses the design point among sets
M and E with the greatest uncertainty, as quantified by the width wi ( x ) =u( x ) − ℓ( x ):wi (x)(16.17)x (i) = arg max
Sx∈Mi

SafeOpt cannot guarantee safety
if it is not initialized with at least
one point that it knows is safe.

6

Ei

SafeOpt proceeds until a termination condition is met. It is common to run thealgorithm for a fixed number of iterations or until the maximum width is lessthan a set threshold.
Maintaining sets in multidimensional spaces can be computationally challenging. SafeOpt assumes a finite design space X that can be obtained with a sampling
method applied over the continuous search domain. Increasing the density of thefinite design space leads to more accurate results with respect to the continuous
space, but it takes longer per iteration.
SafeOpt is implemented in algorithm 16.3, and calls algorithm 16.4 to update
the predicted confidence intervals; algorithm 16.5 to compute the safe, minimizer,
and expander regions; and algorithm 16.6 to select a query point. The progressionof SafeOpt is shown for one dimension in figure 16.12, and for two dimensions infigure 16.13.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

300 c ha p te r 16 . su rrog a te op ti mi z a ti on

function safe_opt(GP, X, i, f, y_max; β=3.0, k_max=10)
push!(GP, X[i], f(X[i])) # make first observationm = length(X)u, l = fill(Inf, m), fill(-Inf, m)
S, M, E = falses(m), falses(m), falses(m)for k in 1 : k_max
update_confidence_intervals!(GP, X, u, l, β)
compute_sets!(GP, S, M, E, X, u, l, y_max, β)
i = get_new_query_point(M, E, u, l)
i != 0 || break
push!(GP, X[i], f(X[i]))
end

end

Algorithm 16.3. The SafeOpt algorithm applied to an empty Gaussian process GP, a finite design
space X, index of initial safe pointi, objective function f, and safety
threshold y_max. The optional parameters are the confidence scalar
β and the number of iterations
k_max. A tuple containing the bestsafe upper bound and its index in
X is returned.

# return the best point
update_confidence_intervals!(GP, X, u, l, β)S[:] = u .≤ y_max
if any(S)u_best, i_best = findmin(u[S])i_best = findfirst(isequal(i_best), cumsum(S))
return (u_best, i_best)
else
return (NaN,0)
end

function update_confidence_intervals!(GP, X, u, l, β)
μₚ, νₚ = predict(GP, X)
u[:] = μₚ + sqrt.(β*νₚ)l[:] = μₚ - sqrt.(β*νₚ)return (u, l)end

Algorithm 16.4. A method for
updating the lower and upperbounds used in SafeOpt, whichtakes the Gaussian process GP, thefinite search space X, the upper andlower-bound vectors u and l, andthe confidence scalar β.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

16 . 6 . sa fe op ti mi z a ti on 301

function compute_sets!(GP, S, M, E, X, u, l, y_max, β)
fill!(M, false)
fill!(E, false)# safe set
S[:] = u .≤ y_max
if any(S)# potential minimizers
M[S] = l[S] .< minimum(u[S])Algorithm 16.5. A method for updating the safe S, minimizer M, andexpander E sets used in SafeOpt.
The sets are all Boolean vectors indicating whether the corresponding design point in X is in the set.
The method also takes the Gaussian process GP, the upper andlower bounds u and l, respectively,
the safety threshold y_max, and theconfidence scalar β.

# maximum width (in M)
w_max = maximum(u[M] - l[M])

end
end

# expanders - skip values in M or those with w ≤ w_max
E[:] = S .& .~M # skip points in M
if any(E)E[E] .= maximum(u[E] - l[E]) .> w_max
for (i,e) in enumerate(E)if e && u[i] - l[i] > w_max
push!(GP, X[i], l[i])μₚ, νₚ = predict(GP, X[.~S])
pop!(GP)E[i] = any(μₚ + sqrt.(β*νₚ) .≥ y_max)
if E[i]; w_max = u[i] - l[i]; end
end
end
end

return (S,M,E)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

302 c ha p te r 16 . su rrog a te op ti mi z a ti on

function get_new_query_point(M, E, u, l)ME = M .| E
if any(ME)v = argmax(u[ME] - l[ME])return findfirst(isequal(v), cumsum(ME))
else
return 0
end
end

Algorithm 16.6. A method for obtaining the next query point in
SafeOpt. The index of the point in X
with the greatest width is returned.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

16 . 6 . sa fe op ti mi z a ti on 303

Figure 16.12. The first eight iterations of SafeOpt on a univariate
function. SafeOpt can never reachthe global optimum on the righthand side because it requires crossing an unsafe region. We can only
hope to find the global minima in
our locally reachable safe region.

y

y

y

y

y

y

y

y

y

estimated safe region
potential minimizers
potential expandersx© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

304 c ha p te r 16 . su rrog a te op ti mi z a ti on

Figure 16.13. SafeOpt applied to
the flower function (appendix B.4)
with ymax = 2, a Gaussian process mean of µ(x) = 2.5, varianceν = 0.01, β = 10, a 51 × 51 uniform
grid over the search space, and an
initial point x (1) = [−2.04, 0.96].
The color indicates the value of theupper bound, the cross indicatesthe safe point with the lowest upper bound, and the white contour
line is the estimated safe region.
The objective function with thetrue safe region outlined in white:

x2

2
0

−2

x2

2
0

−2

2
x2

x2

2
0

0

−2

−2

−2

x2

2
0

−2
−2

0
x1

2

−2

0
x1

2

−2

0
x1

2

−2

0
x1

2

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

0
x1

2

16 . 7. su mma ry

305

16.7 Summary• Gaussian processes can be used to guide the optimization process using avariety of strategies that use estimates of quantities such as the lower confidencebound, probability of improvement, and expected improvement.
• Some problems do not allow for the evaluation of unsafe designs, in whichcase we can use safe exploration strategies that rely on Gaussian processes.

16.8

Exercises

Exercise 16.1. Give an example in which prediction-based optimization fails.
Exercise 16.2. What is the main difference between lower confidence boundexploration and error-based exploration in the context of optimization?
Exercise 16.3. We have a function f ( x ) = ( x − 2)2 /40 − 0.5 with x ∈ [−5, 5],
and we have evaluation points at −1 and 1. Assume we use a Gaussian processsurrogate model with a zero-mean function, and a squared exponential kernelexp(−r2 /2), where r is the Euclidean distance between two points. Which value
for x would we evaluate next if we were maximizing probability of improvement?
Which value for x would we evaluate next if we were maximizing expected
improvement?

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com