8

Stochastic Methods

This chapter presents a variety of stochastic methods that use randomization strategically to help explore the design space for an optimum. Randomness can help escape local optima and increase the chances of finding a global optimum. Stochastic
methods typically use pseudo-random number generators to ensure repeatability.1
A large amount of randomness is generally ineffective because it prevents us from
effectively using previous evaluation points to help guide the search. This chapter
discusses a variety of ways to control the degree of randomness in our search.

8.1 Noisy Descent
Adding stochasticity to gradient descent can be beneficial in large nonlinear
optimization problems. Saddle points, where the gradient is very close to zero,can cause descent methods to select step sizes that are too small to be useful. One
approach is to add Gaussian noise at each descent step:2
x(k+1) ← x(k) + αg(k) + ǫ(k)(8.1)

where ǫ(k) is zero-mean Gaussian noise with standard deviation σ. The amountof noise is typically reduced over time. The standard deviation of the noise istypically a decreasing sequence σ(k) such as 1/k.3 Algorithm 8.1 provides animplementation of this method. Figure 8.1 compares descent with and without
noise on a saddle function.
A common approach for training neural networks is stochastic gradient descent,which uses a noisy gradient approximation. In addition to helping traverse pastsaddle points, evaluating noisy gradients using randomly chosen subsets of the
training data4 is significantly less expensive computationally than calculating thetrue gradient at every iteration.

Although pseudo-random number generators produce numbersthat appear random, they are actually a result of a deterministic process. Pseudo-random numbers can be produced through calls
to the rand function. The process
can be reset to an initial state using the seed! function from theRandom.jl package.

1

G. Hinton and S. Roweis,“Stochastic Neighbor Embedding,”
in Advances in Neural Information
Processing Systems (NIPS), 2003.
2

The Hinton and Roweis paper
used a fixed standard deviation for
the first 3,500 iterations and set thestandard deviation to zero thereafter.
3

4

These subsets are called batches.

126 c ha p te r 8 . stoc ha sti c me thod s

Algorithm 8.1. A noisy descentmethod, which augments anotherdescent method with additive
Gaussian noise. The method
takes another DescentMethod
submethod, a noise sequence σ,
and stores the iteration count k.

mutable struct NoisyDescent <: DescentMethod
submethodσkendfunction init!(M::NoisyDescent, f, ∇f, x)
init!(M.submethod, f, ∇f, x)
M.k = 1
return M
endfunction step!(M::NoisyDescent, f, ∇f, x)
x = step!(M.submethod, f, ∇f, x)
σ = M.σ(M.k)
x += σ.*randn(length(x))
M.k += 1
return x
end

∞

∞ 

∑ α(k) = ∞∑ α(k)

k =1

k =1

2

<∞

x2

Convergence guarantees for stochastic gradient descent require that the positive
step sizes be chosen such that:(8.2)

These conditions ensure that the step sizes decrease and allow the method toconverge, but not too quickly so as to become stuck away from a local minimum.

8.2 Mesh Adaptive Direct Search
The generalized pattern search methods covered in section 7.4 restricted local
exploration to a fixed set of directions. In contrast, mesh adaptive direct search uses
random positive spanning directions.5
The procedure used to sample positive spanning sets (see example 8.1) begins
by constructing an initial linearly spanning set in the form
√ of a lower triangularmatrix L. The diagonal terms in L are sampled from ±1/ α(k) , where α(k) is thestep size at iteration k. The lower components of L are sampled from
n
o
p
p
p
−1/ α(k) + 1, −1/ α(k) + 2, . . . , 1/ α(k) − 1
(8.3)

x1
stochastic gradient descentsteepest descentFigure 8.1. Adding stochasticity
to a descent method helps withtraversing saddle points such as
f (x) = x12 − x22 shown here. Dueto the initialization, the steepest
descent method converges to thesaddle point where the gradient is
zero.
5This section follows the lower
triangular mesh adaptive direct
search given by C. Audet and J. E.
Dennis Jr., “Mesh Adaptive Direct
Search Algorithms for Constrained
Optimization,” SIAM Journal on Optimization, vol. 17, no. 1, pp. 188–
217, 2006.

© 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

8 . 2. me sh a d a p ti ve d i re c t se a rc h 127

Consider positive spanning sets constructed from the nonzero directions
d1 , d2 ∈ {−1, 0, 1}. There are 8 positive spanning sets with 3 elements thatcan be constructed from these directions:

Example 8.1. Positive spanningsets for R2 . Note that the lower
triangular generation strategy can
only generate the first two columnsof spanning sets.

The rows and columns of L are then randomly permuted to obtain a matrix Dnwhose columns correspond to n directions√ that linearly span R . The maximum
magnitude among these directions is 1/ α(k) .
Two common methods for obtaining a positive spanning set from the linearlyspanning set are to add one additional direction d(n+1) = − ∑in=1 d(i) and to add
n additional directions d(n+ j) = −d( j) for j in {1, . . . , n}. We use the first method
in algorithm 8.2.
The step size α starts at 1, is always a power of 4, and never exceeds 1. Using a
power of 4 causes the maximum possible step size taken in each iteration to be√√scaled by a factor of 2, as the maximum step size α/ α has length 4m / 4m = 2m
for integer m < 1. The step size is updated according to:α(k) /4
if no improvement was found in this iteration
α ( k +1) ←min(1, 4α(k) ) otherwise(8.4)
Mesh adaptive direct search is opportunistic but cannot support dynamicordering6 since, after a successful iteration, the step size is increased, and anotherstep in the successful direction would lie outside of the mesh. The algorithm
queries a new design point along the accepted descent direction. If f (x(k) =x(k−1) + αd) < f (x(k−1) ), then the queried point is x(k−1) + 4αd = x(k) + 3αd. The

6

See section 7.4.

© 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

128

c ha p te r 8 . stoc ha sti c me thod sfunction rand_positive_spanning_set(α, n)δ = round(Int, 1/sqrt(α))
L = Matrix(Diagonal(δ*rand([1,-1], n)))for i in 1 : n-1
for j in 1:i-1
L[i,j] = rand(-δ+1:δ-1)
end
end
D = L[randperm(n),:]
D = D[:,randperm(n)]D = hcat(D, -sum(D,dims=2))
return [D[:,i] for i in 1 : n+1]
end

Algorithm 8.2. Randomly sampling a positive spanning set ofn + 1 directions according to mesh
adaptive direct search with step
size α and number of dimensions
n.

procedure is outlined in algorithm 8.3. Figure 8.2 illustrates how this algorithm
explores the search space.

8.3 Simulated Annealing
Simulated annealing7 borrows inspiration from metallurgy.8 Temperature is used
to control the degree of stochasticity during the randomized search. The temperature starts high, allowing the process to freely move about the search space,with the hope that in this phase the process will find a good region with thebest local minimum. The temperature is then slowly brought down, reducing
the stochasticity and forcing the search to converge to a minimum. Simulated
annealing is often used on functions with many local minima due to its ability toescape local minima.
At every iteration, a candidate transition from x to x′ is sampled from a transition distribution T and is accepted with probability1
if ∆y ≤ 0(8.5)
e−∆y/t if ∆y > 0

7
S. Kirkpatrick, C. D. Gelatt Jr.,and M. P. Vecchi, “Optimization
by Simulated Annealing,” Science,
vol. 220, no. 4598, pp. 671–680,
1983.
8
Annealing is a process in which amaterial is heated and then cooled,making it more workable. Whenhot, the atoms in the material are
more free to move around, and,through random motion, tend tosettle into better positions. A slow
cooling brings the material to anordered, crystalline state. A fast,
abrupt quenching causes defects
because the material is forced tosettle in its current condition.

where ∆y = f (x′ ) − f (x) is the difference in the objective and t is the temperature.
It is this acceptance probability, known as the Metropolis criterion, that allows thealgorithm to escape from local minima when the temperature is high.

© 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

8 . 3. si mu l a te d a n n e a l i n g 129

function mesh_adaptive_direct_search(f, x, ϵ)
α, y, n = 1, f(x), length(x)while α > ϵimproved = falsefor (i,d) in enumerate(rand_positive_spanning_set(α, n))x′ = x + α*dy′ = f(x′)
if y′ < yx, y, improved = x′, y′, truex′ = x + 3α*dy′ = f(x′)
if y′ < yx, y = x′, y′
end
break
end
end
α = improved ? min(4α, 1) : α/4
end
return x
end

Algorithm 8.3. Mesh adaptive direct search for an objective function
f , an initial design x, and a tolerance ϵ.

Figure 8.2. Mesh adaptive direct
search proceeding left to right and
top to bottom.
© 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

130 c ha p te r 8 . stoc ha sti c me thod st(k+1) = γt(k)10temperature

The temperature parameter t controls the acceptance probability. An annealingschedule is used to slowly bring down the temperature as the algorithm progresses, as illustrated by figure 8.3. The temperature must be brought down toensure convergence. If it is brought down too quickly, the search method may notcover the portion of the search space containing the global minimum.
It can be shown that a logarithmic annealing schedule of t(k) = t(1) ln(2)/ ln(k + 1)
for the kth iteration is guaranteed to asymptotically reach the global optimum
under certain conditions,9 but it can be very slow in practice. The exponential
annealing schedule, which is more common, uses a simple decay factor:

0100function simulated_annealing(f, x, T, t, k_max)
y = f(x)x_best, y_best = x, y
for k in 1 : k_max
x′ = x + rand(T)y′ = f(x′)Δy = y′ - y
if Δy ≤ 0 || rand() < exp(-Δy/t(k))
x, y = x′, y′end
if y′ < y_best
x_best, y_best = x′, y′end
end
return x_bestend

A more sophisticated algorithm was introduced by Corana et al. in 1987 thatallows for the step size to change during the search.11 Rather than using a fixed102
iteration

104

logarithmic
exponential, γ = 1/4
exponential, γ = 1/2
exponential, γ = 3/4
fast

(8.6)

for some γ ∈ (0, 1). Another common annealing schedule, fast annealing,10 uses atemperature of
t (1)
(8.7)
t(k) =k
A basic implementation of simulated annealing is provided by algorithm 8.4.
Example 8.2 shows the effect different transition distributions and annealingschedules have on the optimization process.

5

Figure 8.3. Several annealing schedules commonly used in simulated
annealing. The schedules have an
initial temperature of 10.
9
B. Hajek, “Cooling Schedules for
Optimal Annealing,” Mathematics
of Operations Research, vol. 13, no. 2,
pp. 311–329, 1988.

Algorithm 8.4. Simulated annealing, which takes as input an objective function f, an initial point x, atransition distribution T, an annealing schedule t, and the number of
iterations k_max.
H. Szu and R. Hartley, “Fast Simulated Annealing,” Physics Letters
A, vol. 122, no. 3-4, pp. 157–162,
1987.
10

A. Corana, M. Marchesi, C. Martini, and S. Ridella, “Minimizing
Multimodal Functions of Continuous Variables with the ‘Simulated Annealing’ Algorithm,” ACM
Transactions on Mathematical Software, vol. 13, no. 3, pp. 262–280,
1987.
11

© 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

8 . 3. si mu l a te d a n n e a l i n g 131

We can use simulated annealing to optimize Ackley’s function, appendix B.1.
Ackley’s function has many local minima, making it easy for gradient-based
methods to get stuck.
Suppose we start at x(1) = [15, 15] and run 100 iterations. Below we show
the distribution over iterations for multiple runs with different combinationsof three zero-mean, diagonal covariance (σI) Gaussian transition distributions, and three different temperature schedules t(k) = t(1) /k.
σ = 1, t(1) = 1

σ = 1, t(1) = 10

σ = 1, t(1) = 25

σ = 5, t(1) = 1

σ = 5, t(1) = 10

σ = 5, t(1) = 25

σ = 25, t(1) = 1

σ = 25, t(1) = 10

σ = 25, t(1) = 25

Example 8.2. Exploring the effectof distribution variance and temperature on the performance of
simulated annealing. The blue regions indicate the 5% to 95% and
25% to 75% empirical Gaussianquantiles of the objective function
value.

30

y

20
10
0

30

y

20
10
0

30

y

20
10
0
20 40 60 80 100
iteration

20 40 60 80 100
iteration

20 40 60 80 100
iteration

In this case, the spread of the transition distribution has the greatest impact
on performance.

© 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

132 c ha p te r 8 . stoc ha sti c me thod stransition distribution, this adaptive simulated annealing method keeps track ofa separate step size v for each coordinate direction. For a given point x, a cycle ofrandom moves is performed in each coordinate direction i according to:(8.8)

where r is drawn uniformly at random from [−1, 1] and vi is the maximum step
size in the ith coordinate direction. Each new point is accepted according to theMetropolis criterion. The number of accepted points in each coordinate direction
is stored in a vector a.
After ns cycles, the step sizes are adjusted with the aim to maintain an approximately equal number of accepted and rejected designs with an average acceptance
rate near one-half. Rejecting too many moves is a waste of computational effort,
while accepting too many moves indicates that the configuration is evolving too
slowly because candidate points are too similar to the current location. The update
formula used by Corana et al. is: 
ai /ns −0.6

if ai > 0.6ns
v
1
+
c

i
0.4

 i
 −1
vi = vi 1 + ci 0.4− ai /ns(8.9)
if ai < 0.4ns
0.4



v
otherwisei
The ci parameter controls the step variation along each direction and is typically
set to 2 as shown in figure 8.4. An implementation of the update is shown inalgorithm 8.5.
function corana_update!(v, a, c, ns)
for i in 1 : length(v)ai, ci = a[i], c[i]if ai > 0.6nsv[i] *= (1 + ci*(ai/ns - 0.6)/0.4)
elseif ai < 0.4nsv[i] /= (1 + ci*(0.4-ai/ns)/0.4)
end
end
return v
end

Step Multiplication Factor

x′ = x + rvi ei

1+c

1
1
1+ c

0

0.4 0.6
1
Acceptance Ratio ai /ns

Figure 8.4. The step multiplication
factor as a function of acceptance
rate for c = 2.
Algorithm 8.5. The update formula used by Corana et al. in adaptive simulated annealing, where v
is a vector of coordinate step sizes,
a is a vector of the number of accepted steps in each coordinate direction, c is a vector of step scaling
factors for each coordinate direction, and ns is the number of cycles
before running the step size adjustment.

© 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

8 . 4 . c ross-e n trop y me thod

133

Temperature reduction occurs every nt step adjustments, which is every ns · nt
cycles. The original implementation simply multiplies the temperature by a reduction factor.
The process is terminated when the temperature sinks low enough such that
improvement can no longer be expected. Termination occurs when the most recent
function value is no farther than ǫ from the previous nǫ iterations and the best
function value obtained over the course of execution. Algorithm 8.6 provides animplementation and the algorithm is visualized in figure 8.5.

8.4 Cross-Entropy Method
The cross-entropy method,12 in contrast with the methods we have discussed so
far in this chapter, maintains an explicit probability distribution over the design
space.13 This probability distribution, often called a proposal distribution, is used topropose new samples for the next iteration. At each iteration, we sample from the
proposal distribution and then update the proposal distribution to fit a collectionof the best samples. The aim at convergence is for the proposal distribution to
focus on the global optima. Algorithm 8.7 provides an implementation.
The cross-entropy method requires choosing a family of distributions parameterized by θ. One common choice is the family of multivariate normal distributions
parameterized by a mean vector and a covariance matrix. The algorithm also
requires us to specify the number of elite samples, melite , to use when fitting theparameters for the next iteration.
Depending on the choice of distribution family, the process of fitting the distribution to the elite samples can be done analytically. In the case of the multivariate
normal distribution, the parameters are updated according to the maximum
likelihood estimate:µ( k +1) =Σ ( k +1) =1

melite1melitemelite i∑=1melite i∑=1x( i )(8.10)(x(i) − µ(k+1) )(x(i) − µ(k+1) )⊤(8.11)

12
R. Y. Rubinstein and D. P. Kroese,
The Cross-Entropy Method: A Unified
Approach to Combinatorial Optimization, Monte-Carlo Simulation, andMachine Learning. Springer, 2004.

The name of this method comesfrom the fact that the process offitting the distribution involves
minimizing cross-entropy, which isalso called the Kullback–Leibler divergence. Under certain conditions,minimizing the cross-entropy corresponds to finding the maximum
likelihood estimate of the parameters of the distribution.

13

Example 8.3 applies the cross-entropy method to a simple function. Figure 8.6
shows several iterations on a more complex function. Example 8.4 shows thepotential limitation of using a multivariate normal distribution for fitting elite
samples.
© 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

134 c ha p te r 8 . stoc ha sti c me thod sfunction adaptive_simulated_annealing(f, x, v, t, ϵ;
ns=20, nϵ=4, nt=max(100,5length(x)),γ=0.85, c=fill(2,length(x)) )y = f(x)x_best, y_best = x, yy_arr, n, U = [], length(x), Uniform(-1.0,1.0)a,counts_cycles,counts_resets = zeros(n), 0, 0
while true
for i in 1:nx′ = x + basis(i,n)*rand(U)*v[i]y′ = f(x′)Δy = y′ - y
if Δy < 0 || rand() < exp(-Δy/t)x, y = x′, y′a[i] += 1
if y′ < y_best; x_best, y_best = x′, y′; end
end
end
counts_cycles += 1
counts_cycles ≥ ns || continue
counts_cycles = 0
corana_update!(v, a, c, ns)
fill!(a, 0)counts_resets += 1
counts_resets ≥ nt || continue

Algorithm 8.6. The adaptive simulated annealing algorithm, where
f is the multivariate objective function, x is the starting point, v isstarting step vector, t is the starting temperature, and ϵ is the termination criterion parameter. The optional parameters are the numberof cycles before running the step
size adjustment ns, the number ofcycles before reducing the temperature nt, the number of successive temperature reductions to test
for termination nϵ, the temperature reduction coefficient γ, and the
direction-wise varying criterion c.
Below is a flowchart for theadaptive simulated annealing algorithm as presented in the original
paper.
Initialize Parameters
Perform a cycle of random moves,each along one coordinate direction. Accept or reject each pointaccording to the Metropolis criterion. Record the optimum point
reached so far.

t *= γcounts_resets = 0
push!(y_arr, y)
if !(length(y_arr) > nϵ && y_arr[end] - y_best ≤ ϵ &&all(abs(y_arr[end]-y_arr[end-u]) ≤ ϵ for u in 1:nϵ))x, y = x_best, y_best
else
break
end

end

end
return x_best© 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

No. cycles ≥ ns ?
yesnoAdjust step vector v.
Reset no. cycles to 0.
No. step adjustments ≥ nt ?
yesnoReduce temperature.
Reset no. adjustments to 0.
Add point to optimum list.
Stopping criterion satisfied?
yes
End

no

135

y

y

8 . 4 . c ross-e n trop y me thod

x

x

x

x
Figure 8.5. Simulated annealing
with an exponentially decaying
temperature, where the histograms
indicate the probability of simulated annealing being at a particular position at that iteration.

using Distributions
function cross_entropy_method(f, P, k_max, m=100, m_elite=10)
for k in 1 : k_max
samples = rand(P, m)order = sortperm([f(samples[:,i]) for i in 1:m])
P = fit(typeof(P), samples[:,order[1:m_elite]])end
return P
end

Algorithm 8.7. The cross-entropy
method, which takes an objectivefunction f to be minimized, a proposal distribution P, an iteration
count k_max, a sample size m, andthe number of samples to use when
refitting the distribution m_elite.
It returns the updated distribution
over where the global minimum islikely to exist.

© 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

136

c ha p te r 8 . stoc ha sti c me thod s

We can use Distributions.jl to represent, sample from, and fit proposal
distributions. The parameter vector θ is replaced by a distribution P. Callingrand(P,m) will produce an n × m matrix corresponding to m samples ofn-dimensional samples from P, and calling fit will fit a new distribution of
the given input type.

Example 8.3. An example of using
the cross-entropy method.

import Random: seed!
import LinearAlgebra: norm
seed!(0) # set random seed for reproducible results
f = x->norm(x)μ = [0.5, 1.5]
Σ = [1.0 0.2; 0.2 2.0]
P = MvNormal(μ, Σ)
k_max = 10
P = cross_entropy_method(f, P, k_max)@show P.μx2

P.μ = [-6.136226751273704e-7, -1.372160422208573e-6]

x1

x1

x1

x1

Figure 8.6. The cross-entropy
method with m = 40 applied to
the Branin function (appendix B.3)
using a multivariate Gaussian proposal distribution. The 10 elite samples in each iteration are in red.
© 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

8 . 5 . n a tu ra l e vol u ti on stra te g i e s 137

The distribution family should be flexible enough to capture the relevant
features of the objective function. Here we show the limitations of using
a normal distribution on a multimodal objective function, which assigns
greater density in between the two minima. A mixture model is able to center
itself over each minimum.
0.3

Example 8.4. The normal distribution is unable to capture multiple
local minima, in contrast to mixture models which can maintainseveral.

f
mixture model fit
normal fit
samples

0.2

0.1

0

8.5

−5

0
x

5

Natural Evolution Strategies

Like the cross-entropy method, natural evolution strategies14 optimize a proposal
distribution parameterized by θ. We have to specify the proposal distribution
family and the number of samples. The aim is to minimize the expectation
Ex∼ p(·|θ) [ f (x)]. Instead of fitting elite samples, evolution strategies apply gradientdescent. The gradient is estimated from the samples:15∇θ Ex∼ p(·|θ) [ f (x)] ===Z∇θ p(x | θ) f (x) dx(8.12)(8.13)Zp (x | θ)∇ p(x | θ) f (x) dx
p (x | θ) θp(x | θ)∇θ log p(x | θ) f (x) dx

Z(8.14)= Ex∼ p(·|θ) [ f (x)∇θ log p(x | θ)](8.15)≈(8.16)m1
f (x(i) )∇θ log p(x(i) | θ)
m i∑
=1

14
I. Rechenberg, Evolutionsstrategie Optimierung technischer Systeme
nach Prinzipien der biologischen Evolution. Frommann-Holzboog, 1973.

This gradient estimation has recently been successfully applied
to proposal distributions represented by deep neural networks.
T. Salimans, J. Ho, X. Chen, andI. Sutskever, “Evolution Strategies as a Scalable Alternative toReinforcement Learning,” ArXiv,
no. 1703.03864, 2017.
15

Although we do not need the gradient of the objective function, we do need thegradient of the log likelihood, log p(x | θ). Example 8.5 shows how to compute© 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

138

c ha p te r 8 . stoc ha sti c me thod s

the gradient of the log likelihood for the multivariate normal distribution. The
estimated gradient can be used along with any of the descent methods discussed
in previous chapters to improve θ. Algorithm 8.8 uses gradient descent with afixed step size. Figure 8.7 shows a few iterations of the algorithm.
using Distributions
function natural_evolution_strategies(f, θ, k_max; m=100, α=0.01)
for k in 1 : k_max
samples = [rand(θ) for i in 1 : m]θ -= α*sum(f(x)*∇logp(x, θ) for x in samples)/m
end
return θ
end

Algorithm 8.8. The natural evolution strategies method, whichtakes an objective function f to
be minimized, an initial distribution parameter vector θ, an iteration count k_max, a sample size m,
and a step factor α. An optimized
parameter vector is returned. The
method rand(θ) should sample
from the distribution parameterized by θ, and ∇logp(x, θ) should
return the log likelihood gradient.

8.6 Covariance Matrix Adaptation
Another popular method is covariance matrix adaptation,16 which is also referred toas CMA-ES for covariance matrix adaptation evolutionary strategy. It has similarities
with natural evolution strategies from section 8.5, but the two should not be
confused. This method maintains a covariance matrix and is robust and sample
efficient. Like the cross-entropy method and natural evolution strategies, a distribution is improved over time based on samples. Covariance matrix adaptation
uses multivariate Gaussian distributions.17
Covariance matrix adaptation maintains a mean vector µ, a covariance matrix
Σ, and an additional step-size scalar σ. The covariance matrix only increases ordecreases in a single direction with every iteration, whereas the step-size scalaris adapted to control the overall spread of the distribution. At every iteration, mdesigns are sampled from the multivariate Gaussian:18
x ∼ N (µ, σ2 Σ)(8.17)The designs are then sorted according to their objective function values suchthat f (x(1) ) ≤ f (x(2) ) ≤ · · · ≤ f (x(m) ). A new mean vector µ(k+1) is formed
using a weighted average of the first melite sampled designs:
meliteµ( k +1) ← ∑ wi x( i )It is common to use the phrase
evolution strategies to refer specifically to covariance matrix adaptation.

16

N. Hansen, “The CMA Evolution Strategy: A Tutorial,” ArXiv,
no. 1604.00772, 2016.

17

For optimization in R n , it is recommended to use at least m =4 + ⌊3 ln n⌋ samples per iteration,and melite = ⌊m/2⌋ elite samples.

18(8.18)

i =1© 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

8 . 6 . c ova ri a n c e ma tri x a d a p ta ti on

The multivariate normal distribution N (µ, Σ) with mean µ and covarianceΣ is a popular distribution family due to having analytic solutions. The
likelihood in d dimensions has the form


1
− 12
⊤ −1
− 2d
p(x | µ, Σ) = (2π ) |Σ| exp − (x − µ) Σ (x − µ)2where |Σ| is the determinant of Σ. The log likelihood isd1
1
log p(x | µ, Σ) = − log(2π ) − log |Σ| − (x − µ)⊤ Σ−1 (x − µ)
2
2
2

139

Example 8.5. A derivation of thelog likelihood gradient equations
for the multivariate Gaussian distribution. For the original derivation and several more sophisticatedsolutions for handling the positive
definite covariance matrix, see D.
Wierstra, T. Schaul, T. Glasmachers,
Y. Sun, and J. Schmidhuber, “Natural Evolution Strategies,” ArXiv,
no. 1106.4487, 2011.

The parameters can be updated using their log likelihood gradients:∇(µ) log p(x | µ, Σ) = Σ−1 (x − µ)11∇(Σ) log p(x | µ, Σ) = Σ−1 (x − µ)(x − µ)⊤ Σ−1 − Σ−122

Figure 8.7. Natural evolution strategies using multivariate Gaussian
distributions applied to Wheeler’s
Ridge, appendix B.7.

x2

The term ∇(Σ) contains the partial derivative of each entry of Σ with respect
to the log likelihood.
Directly updating Σ may not result in a positive definite matrix, as is required for covariance matrices. One solution is to represent Σ as a productA⊤ A, which guarantees that Σ remains positive semidefinite, and then update A instead. Replacing Σ by A⊤ A and taking the gradient with respect toA yields:hi∇(A) log p(x | µ, A) = A ∇(Σ) log p(x | µ, Σ) + ∇(Σ) log p(x | µ, Σ)⊤

x1

x1

x1

x1

© 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

140 c ha p te r 8 . stoc ha sti c me thod s

where the first melite weights sum to 1, and all the weights approximately sum to
0 and are ordered largest to smallest:19
melitem∑ wi ≈ 0∑ wi = 1

w1 ≥ w2 ≥ · · · ≥ w m

i =1

i =1

(8.19)

We can approximate the mean update in the cross-entropy method by setting
the first melite weights to 1/melite , and setting the remaining weights to zero. Covariance matrix adaptation instead distributes decreasing weight to all m designs,
including some negative weights. The recommended weighting is obtained bynormalizing
m+1
wi′ = ln
− ln i for i ∈ {1, . . . , m}
(8.20)
2
to obtain w. The positive and negative weights are normalized separately. Figure 8.8 compares the mean updates for covariance matrix adaptation and thecross-entropy method.
The step size is updated using a cumulative variable pσ that tracks steps over
time:
(1)
pσ = 0q(8.21)( k +1)← (1 − cσ )pσ + cσ (2 − cσ )µeff (Σ(k) )−1/2 δw
pσ
where cσ < 1 controls the rate of decay and the right hand term determines
whether the step size should be increased or decreased based on the observed
samples with respect to the present scale of the distribution. The variance effective
selection mass µeff has the form
µeff =

1

In the recommended weighting,the first melite samples are positive, and the remaining samples
are non-positive.

19

Figure 8.8. Shown is an initial
proposal distribution (white contours), six samples (white dots),and the new updated means forboth covariance matrix adaptation(blue dot) and the cross-entropy
method (red dot) using three elite
samples. Covariance matrix adaptation tends to update the mean
more aggressively than the crossentropy method (red dot), as it assigns higher weight to better sampled designs, and negative weight
to worse sampled designs.

(8.22)m2∑i=elite
1 wiand δw is computed from the sampled deviations:meliteδw = ∑ wi δ(i) for δ(i) =i =1x( i ) − µ( k )σ (k)The new step size is obtained according to
 

cσ
k pσ kσ(k+1) ← σ(k) exp−1dσ E kN (0, I)k

where

√ Γ
E kN (0, I)k = 2



n +1
2

n
Γ 2



≈

√



1
1
n 1−
+
4n 21n2(8.23)(8.24)(8.25)© 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

8 . 6 . c ova ri a n c e ma tri x a d a p ta ti on 141

is the expected length of a vector drawn from a Gaussian distribution. Comparing
the length of pσ to its expected length under random selection provides themechanism by which σ is increased or decreased. The constants cσ and dσ haverecommended values:cσ = (µeff + 2)/(n + µeff + 5) q
dσ = 1 + 2 max 0, (µeff − 1)/(n + 1) − 1 + cσ

(8.26)

The covariance matrix is also updated using a cumulative vector:(1)pΣ = 0( k +1)

pΣ
where
hσ =

(k)

← (1 − c Σ )pΣ + h σ1 if qkpσ k1−(1−cσ0 otherwise)2( k +1)q

(8.27)cΣ (2 − cΣ )µeff δw

2< 1.4 + n+
1 E kN (0, I)k(8.28)The hσ stalls the update of pΣ if kpΣ k is too large, thereby preventing excessive
increases in Σ when the step size is too small.
The update requires the adjusted weights w◦ :wi
if wi ≥ 0◦wi =(8.29)nwi Σ−1/2 δ(i) 2 otherwise
The covariance update is then






Σ(k+1) ← 1 + c1 cΣ (1 − hσ )(2 − cΣ ) − c1 − cµ Σ(k) +|{z}typically zero

The constants cΣ , c1 and cµ have recommended values

4 + µeff /n
n + 4 + 2µeff /n2
c1 =(n + 1.3)2 + µeff
µ − 2 + 1/µeff
cµ = min 1 − c1 , 2 eff(n + 2)2 + µeff

c 1 pΣ p⊤| {z Σ}rank-one updateµ⊤+ cµ ∑ wi◦ δ(i) δ(i)

|

i =1

{z

rank-µ update(8.30)}cΣ =(8.31)© 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

142 c ha p te r 8 . stoc ha sti c me thod s

The covariance update consists of three components: the previous covariancematrix Σ(k) , a rank-one update, and a rank-µ update. The rank-one update gets
its name from the fact that pΣ p⊤
Σ has rank one; it has only one eigenvector alongpΣ . Rank-one updates using the cumulation vector allow for correlations between
consecutive steps to be exploited, permitting the covariance matrix to elongate
itself more quickly along a favorable axis.
⊤µThe rank-µ update gets its name from the fact that ∑i=1 wi◦ δ(i) δ(i)has rank
min(µ, n). One important difference between the empirical covariance matrix
update used by the cross-entropy method and the rank-µ update is that theformer estimates the covariance about the new mean µ(k+1) , whereas the latter
estimates the covariance about the original mean µ(k) . The δ(i) values thus help
estimate the variances of the sampled steps rather than the variance within thesampled designs.
Covariance matrix adaptation is depicted in figure 8.9.

8.7 Summary• Stochastic methods employ random numbers during the optimization process.
• Simulated annealing uses a temperature that controls random exploration andwhich is reduced over time to converge on a local minimum.
• The cross-entropy method and evolution strategies maintain proposal distributions from which they sample in order to inform updates.
• Natural evolution strategies uses gradient descent with respect to the log
likelihood to update its proposal distribution.
• Covariance matrix adaptation is a robust and sample-efficient optimizer thatmaintains a multivariate Gaussian proposal distribution with a full covariancematrix.

8.8 Exercises
Exercise 8.1. We have shown that mixture proposal distributions can better
capture multiple minima. Why might their use in the cross-entropy method belimited?
© 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

8 . 8 . e x e rc i se s 143

function covariance_matrix_adaptation(f, x, k_max;
σ = 1.0,m = 4 + floor(Int, 3*log(length(x))),m_elite = div(m,2))μ, n = copy(x), length(x)ws = log((m+1)/2) .- log.(1:m)
ws[1:m_elite] ./= sum(ws[1:m_elite])
μ_eff = 1 / sum(ws[1:m_elite].^2)
cσ = (μ_eff + 2)/(n + μ_eff + 5)
dσ = 1 + 2max(0, sqrt((μ_eff-1)/(n+1))-1) + cσcΣ = (4 + μ_eff/n)/(n + 4 + 2μ_eff/n)
c1 = 2/((n+1.3)^2 + μ_eff)cμ = min(1-c1, 2*(μ_eff-2+1/μ_eff)/((n+2)^2 + μ_eff))
ws[m_elite+1:end] .*= -(1 + c1/cμ)/sum(ws[m_elite+1:end])
E = n^0.5*(1-1/(4n)+1/(21*n^2))pσ, pΣ, Σ = zeros(n), zeros(n), Matrix(1.0I, n, n)
for k in 1 : k_max
P = MvNormal(μ, σ^2*Σ)
xs = [rand(P) for i in 1 : m]
ys = [f(x) for x in xs]
is = sortperm(ys) # best to worst

Algorithm 8.9. Covariance matrix
adaptation, which takes an objective function f to be minimized, an
initial design point x, and an iteration count k_max. One can optionally specify the step-size scalar σ,
the sample size m, and the numberof elite samples m_elite.
The best candidate design pointis returned, which is the mean of
the final sample distribution.
The covariance matrix undergoes an additional operation to ensure that it remains symmetric; otherwise small numerical inconsistencies can cause the matrix nolonger to be positive definite.
This implementation uses a
simplified normalization strategy
for the negative weights. The original can be found in Equations
50–53 of N. Hansen, “The CMA
Evolution Strategy: A Tutorial,”
ArXiv, no. 1604.00772, 2016.

# selection and mean update
δs = [(x - μ)/σ for x in xs]
δw = sum(ws[i]*δs[is[i]] for i in 1 : m_elite)
μ += σ*δw# step-size control
C = Σ^-0.5
pσ = (1-cσ)*pσ + sqrt(cσ*(2-cσ)*μ_eff)*C*δw
σ *= exp(cσ/dσ * (norm(pσ)/E - 1))# covariance adaptation
hσ = Int(norm(pσ)/sqrt(1-(1-cσ)^(2k)) < (1.4+2/(n+1))*E)pΣ = (1-cΣ)*pΣ + hσ*sqrt(cΣ*(2-cΣ)*μ_eff)*δw
w0 = [ws[i]≥0 ? ws[i] : n*ws[i]/norm(C*δs[is[i]])^2
for i in 1:m]
Σ = (1-c1-cμ) * Σ +c1*(pΣ*pΣ' + (1-hσ) * cΣ*(2-cΣ) * Σ) +cμ*sum(w0[i]*δs[is[i]]*δs[is[i]]' for i in 1 : m)
Σ = triu(Σ)+triu(Σ,1)' # enforce symmetry

end

end
return μ© 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

c ha p te r 8 . stoc ha sti c me thod s

x2

x2

x2

x2

144

x1

x1

x1

x1
Figure 8.9. Covariance matrix adaptation using multivariate Gaussian
distributions applied to the flower
function, appendix B.4.

© 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

8 . 8 . e x e rc i se s 145

Exercise 8.2. In the cross-entropy method, what is a potential effect of using anelite sample size that is very close to the total sample size?
Exercise 8.3. The log-likelihood of a value sampled from a Gaussian distribution
with mean µ and variance ν is:1
1( x − µ )2ℓ( x | µ, ν) = − ln 2π − ln ν −222ν
Show why evolution strategies using Gaussian distributions may encounter difficulties while applying a descent update on the variance when the mean is on the
optimum, µ = x ∗ .
Exercise 8.4. Derive the maximum likelihood estimate for the cross-entropy
method using multivariate normal distributions:µ( k +1) =1 m (i )xm i∑=1

Σ ( k +1) =1 m (i )(x − µ(k+1) )(x(i) − µ(k+1) )⊤m i∑=1where the maximum likelihood estimates aren the parameterovalues that maximize
the likelihood of sampling the individuals x(1) , · · · , x(m) .

© 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