10 Policy Search

Policy search involves searching the space of policies without directly computing a
value function. The policy space is often lower-dimensional than the state space
and can often be searched more efficiently. Policy optimization optimizes the
parameters in a parameterized policy in order to maximize utility. This parameter-
ized policy can take many forms, such as neural networks, decision trees, and
computer programs. This chapter begins by discussing a way to estimate the value
of a policy given an initial state distribution. We will then discuss search methods
that do not use estimates of the gradient of the policy, saving gradient methods
for the next chapter. Although local search can be quite effective in practice, we
will also discuss a few alternative optimization approaches that can avoid local
optima.1 1
There are many other optimiza-
tion approaches, as discussed
by M. J. Kochenderfer and
10.1 Approximate Policy Evaluation T. A. Wheeler, Algorithms for
Optimization. MIT Press, 2019.

As introduced in section 7.2, we can compute the expected discounted return

when following a policy π from a state s. This expected discounted return U π (s)
can be computed iteratively (algorithm 7.3) or through matrix operations (algo-
rithm 7.4) when the state space is discrete and relatively small. We can use these
results to compute the expected discounted return of π:

U (π ) = ∑ b ( s )U π ( s ) (10.1)
s

assuming an initial state distribution b(s).

We will use this definition of U (π ) throughout this chapter. However, we often
cannot compute U (π ) exactly when the state space is large or continuous. Instead,
we can approximate U (π ) by sampling trajectories, consisting of states, actions,
214 c ha p te r 10 . p ol i c y se a rc h

and rewards when following π. The definition of U (π ) can be rewritten as

Z
U (π ) = E τ [ R(τ )] = pπ (τ ) R(τ ) dτ (10.2)
U (π )
where pπ (τ ) is the probability density associated with trajectory τ when following
policy π, starting from initial state distribution b. The trajectory reward R(τ ) is the
discounted return associated with τ. Figure 10.1 illustrates the computation of
U (π ) in terms of trajectories sampled from an initial state distribution.
Figure 10.1. The utility associated
Monte Carlo policy evaluation (algorithm 10.1) involves approximating equa- with a policy from an initial state
tion (10.2) with m trajectory rollouts of π: distribution is computed from the
return associated with all possible
1 m trajectories under the given policy,
U (π ) ≈ ∑ R ( τ (i ) ) (10.3) weighted according to their likeli-
m i =1 hood.

where τ (i) is the ith trajectory sample.

struct MonteCarloPolicyEvaluation Algorithm 10.1. Monte Carlo pol-

𝒫 # problem icy evaluation of a policy π. The
b # initial state distribution method runs m rollouts to depth d
d # depth according to the dynamics speci-
m # number of samples fied by the problem 𝒫. Each rollout
end is run from an initial state sampled
from state distribution b. The final
function (U::MonteCarloPolicyEvaluation)(π) line in this algorithm block evalu-
R(π) = rollout(U.𝒫, rand(U.b), π, U.d) ates a policy π parameterized by
return mean(R(π) for i = 1:U.m) θ, which will be useful in the algo-
end rithms in this chapter that attempt
to find a value of θ that maximizes
(U::MonteCarloPolicyEvaluation)(π, θ) = U(s->π(θ, s)) U.

50
U (π )

Monte Carlo policy evaluation is stochastic. Multiple evaluations of equa-

40
30
tion (10.1) with the same policy can give different estimates. Increasing the 20
number of rollouts decreases the variance of the evaluation, as demonstrated in 100 101 102
figure 10.2. number of samples

We will use πθ to denote a policy parameterized by θ. For convenience, we Figure 10.2. The effect of the depth
will use U (θ ) as shorthand for U (πθ ) in cases where it is not ambiguous. The and sample count for Monte Carlo
parameter θ may be a vector or some other more complex representation. For policy evaluation of a uniform ran-
dom policy on the cart-pole prob-
example, we may want to represent our policy using a neural network with a lem (appendix F.3). The variance
particular structure. We would use θ to represent the weights in the network. decreases as the number of sam-
ples increases. The blue regions in-
Many optimization algorithms assume that θ is a vector with a fixed number of dicate the 5 % to 95 % and 25 % to
75 % empirical quantiles of U (π ).
© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
 10 . 2. l oc a l se a rc h 215

components. Other optimization algorithms allow more flexible representations,

including representations like decision trees or computational expressions.2 2
We will not be discussing
those representations here,
but some are implemented in
10.2 Local Search [Link].

A common approach to optimization is local search, where we begin with an

initial parameterization and incrementally move from neighbor to neighbor in
the search space until convergence occurs. We discussed this type of approach in
chapter 5, in the context of optimizing Bayesian network structures with respect
to the Bayesian score. Here, we are optimizing policies parameterized by θ. We
are trying to find a value of θ that maximizes U (θ ).
There are many local search algorithms, but this section will focus on the
Hooke-Jeeves method (algorithm 10.2).3 This algorithm assumes that our policy 3
R. Hooke and T. A. Jeeves, “Direct
is parameterized by an n-dimensional vector θ. The algorithm takes a step of Search Solution of Numerical and
Statistical Problems,” Journal of the
size ±α in each of the coordinate directions from the current θ. These 2n points ACM (JACM), vol. 8, no. 2, pp. 212–
correspond to the neighborhood of θ. If no improvements to the policy are found, 229, 1961.
then the step size α is decreased by some factor. If an improvement is found, it
moves to the best point. The process continues until α drops below some threshold
ǫ > 0. An example involving policy optimization is provided in example 10.1,
and figure 10.3 illustrates this process.

10.3 Genetic Algorithms

4
D. E. Goldberg, Genetic Algorithms
in Search, Optimization, and Machine
A potential issue with local search algorithms like the Hooke-Jeeves method is Learning. Addison-Wesley, 1989.
that the optimization can get stuck in a local optimum. There are a wide variety of
approaches that involve maintaining a population consisting of samples of points 5
F. P. Such, V. Madhavan, E. Conti,
J. Lehman, K. O. Stanley, and J.
in the parameter space, evaluating them in parallel with respect to our objective, Clune, “Deep Neuroevolution: Ge-
and then recombining them in some way to drive the population toward a global netic Algorithms Are a Competi-
optimum. A genetic algorithm4 is one such approach, which derives inspiration tive Alternative for Training Deep
Neural Networks for Reinforce-
from biological evolution. It is a general optimization method, but it has been ment Learning,” 2017. arXiv: 171
successful in the context of optimizing policies. For example, this approach has 2.06567v3. The implementation in
this section follows their relatively
been used to optimize policies for Atari video games, where the policy parameters simple formulation. Their formu-
correspond to weights in a neural network.5 lation does not include crossover,
A simple version of this approach (algorithm 10.3) begins with a population of which is typically used to mix pa-
rameterizations across a popula-
m random parameterizations, θ(1) , . . . , θ(m) . We compute U (θ(i) ) for each sample tion.

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
216 c ha p te r 10 . p ol i c y se a rc h

struct HookeJeevesPolicySearch Algorithm 10.2. Policy search

θ # initial parameterization using the Hooke-Jeeves method,
α # step size which returns a θ that has been op-
c # step size reduction factor timized with respect to U. The pol-
ϵ # termination step size icy π takes as input a parameter θ
end and state s. This implementation
starts with an initial value of θ. The
function optimize(M::HookeJeevesPolicySearch, π, U) step size α is reduced by a factor of
θ, θ′, α, c, ϵ = copy(M.θ), similar(M.θ), M.α, M.c, M.ϵ c if no neighbor improves the ob-
u, n = U(π, θ), length(θ) jective. Iterations are run until the
while α > ϵ step size is less than ϵ.
copyto!(θ′, θ)
best = (i=0, sgn=0, u=u)
for i in 1:n
for sgn in (-1,1)
θ′[i] = θ[i] + sgn*α
u′ = U(π, θ′)
if u′ > best.u
best = (i=i, sgn=sgn, u=u′)
end
end
θ′[i] = θ[i]
end
if best.i != 0
θ[best.i] += [Link]*α
u = best.u
else
α *= c
end
end
return θ
end

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
 10 . 3. g e n e ti c a l g ori thms 217

Suppose we want to optimize a policy for the simple regulator problem Example 10.1. Using a policy opti-
mization algorithm to optimize the
described in appendix F.5. We define a stochastic policy π parameterized by parameters of a stochastic policy.
θ such that the action is generated according to

a ∼ N (θ1 s, (|θ2 | + 10−5 )2 ) (10.4)

The following code defines the parameterized stochastic policy π, evaluation

function U, and method M. It then calls optimize(M, π, U), which returns
an optimized value for θ. In this case, we use the Hooke-Jeeves method, but
the other methods discussed in this chapter can be passed in as M instead:
function π(θ, s)
return rand(Normal(θ[1]*s, abs(θ[2]) + 0.00001))
end
b, d, n_rollouts = Normal(0.3,0.1), 10, 3
U = MonteCarloPolicyEvaluation(𝒫, b, d, n_rollouts)
θ, α, c, ϵ = [0.0,1.0], 0.75, 0.75, 0.01
M = HookeJeevesPolicySearch(θ, α, c, ϵ)
θ = optimize(M, π, U)

Figure 10.3. The Hooke-Jeeves

5
method applied to optimizing a
policy in the simple regulator prob-
3 lem discussed in example 10.1. The
θ2

evaluations at each iteration are

1 shown as white points. Iterations
proceed left to right and top to bot-
tom, and the background is colored
according to the expected utility,
5 with yellow indicating lower util-
ity and dark blue indicating higher
3 utility.
θ2

−2 0 2 −2 0 2 −2 0 2 −2 0 2
θ1 θ1 θ1 θ1

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
218 c ha p te r 10 . p ol i c y se a rc h

i in the population. Since these evaluations potentially involve many rollout

simulations and are therefore computationally expensive, they are often run in
parallel. These evaluations help us identify the elite samples, which are the top
melite samples according to U.
The population at the next iteration is generated by producing m − 1 new pa-
rameterizations by repeatedly selecting a random elite sample θ and perturbing
it with isotropic Gaussian noise, θ + σǫ, where ǫ ∼ N (0, I). The best parame-
terization, unperturbed, is included as the mth sample. Because the evaluations
involve stochastic rollouts, a variation of this algorithm could involve running
additional rollouts to help identify which of the elite samples is truly the best.
Figure 10.4 shows several iterations, or generations, of this approach in a sample
problem.

struct GeneticPolicySearch Algorithm 10.3. A genetic policy

θs # initial population search method for iteratively up-
σ # initial standard deviation dating a population of policy pa-
m_elite # number of elite samples rameterizations θs, which takes a
k_max # number of iterations policy evaluation function U, a pol-
end icy π(θ, s), a perturbation stan-
dard deviation σ, an elite sample
function optimize(M::GeneticPolicySearch, π, U) count m_elite, and an iteration
θs, σ = M.θs, M.σ count k_max. The best m_elite sam-
n, m = length(first(θs)), length(θs) ples from each iteration are used to
for k in 1:M.k_max generate the samples for the subse-
us = [U(π, θ) for θ in θs] quent iteration.
sp = sortperm(us, rev=true)
θ_best = θs[sp[1]]
rand_elite() = θs[sp[rand(1:M.m_elite)]]
θs = [rand_elite() + σ.*randn(n) for i in 1:(m-1)]
push!(θs, θ_best)
end
return last(θs)
end

6
S. Mannor, R. Y. Rubinstein, and Y.
Gat, “The Cross Entropy Method
10.4 Cross Entropy Method for Fast Policy Search,” in Interna-
tional Conference on Machine Learn-
The cross entropy method (algorithm 10.4) involves updating a search distribution ing (ICML), 2003.

over the parameterized space of policies at each iteration.6 We parameterize this

search distribution p(θ | ψ) with ψ.7 This distribution can belong to any family, 7
Often, θ and ψ are vectors, but
because this assumption is not re-
but a Gaussian distribution is a common choice, where ψ represents the mean and
quired for this method, we will not
bold them in this section.
© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
 10 . 5 . e vol u ti on stra te g i e s 219

Figure 10.4. Genetic policy search

5
with σ = 0.25 applied to the simple
regulator problem using 25 sam-
3 ples per iteration. The five elite
θ2

samples in each generation are

1 shown in red, with the best sam-
ple indicated by a larger dot.
−2 0 2 −2 0 2 −2 0 2 −2 0 2
θ1 θ1 θ1 θ1

covariance of the distribution. The objective is to find a value of ψ∗ that maximizes

the expectation of U (θ ) when θ is drawn from the search distribution:
Z
ψ∗ = arg max E [U (θ )] = arg max U (θ ) p(θ | ψ) dθ (10.5)
ψ θ ∼ p(·|ψ) ψ

Directly maximizing equation (10.5) is typically computationally infeasible.

The approach taken in the cross entropy method is to start with an initial value of
ψ, typically chosen so that the distribution is spread over the relevant parameter
space. At each iteration, we draw m samples from the associated distribution and
then update ψ to fit the elite samples. For the fit, we typically use the maximum 8
The maximum likelihood esti-
mate corresponds to the choice of
likelihood estimate (section 4.1).8 We stop after a fixed number of iterations, or
ψ that minimizes the cross entropy
until the search distribution becomes highly focused on an optimum. Figure 10.5 (see appendix A.9) between the
demonstrates the algorithm on a simple problem. search distribution and the elite
samples.

10.5 Evolution Strategies

Evolution strategies9 update a search distribution parameterized by a vector ψ at 9

D. Wierstra, T. Schaul, T. Glas-
each iteration. However, instead of fitting the distribution to a set of elite samples, machers, Y. Sun, J. Peters, and
J. Schmidhuber, “Natural Evolu-
they update the distribution by taking a step in the direction of the gradient.10 tion Strategies,” Journal of Machine
The gradient of the objective in equation (10.5) can be computed as follows:11 Learning Research, vol. 15, pp. 949–
980, 2014.
10
We are effectively doing gradi-
ent ascent, which is reviewed in ap-
pendix A.11.
11
The policy parameter θ is not
bolded here because it is not re-
quired to be a vector. However, ψ
is in bold because we require it to
be a vector when we work with the
gradient of the objective.

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
220 c ha p te r 10 . p ol i c y se a rc h

struct CrossEntropyPolicySearch Algorithm 10.4. Cross entropy pol-

p # initial distribution icy search, which iteratively im-
m # number of samples proves a search distribution ini-
m_elite # number of elite samples tially set to p. This algorithm takes
k_max # number of iterations as input a parameterized policy
end π(θ, s) and a policy evaluation
function U. In each iteration, m
function optimize_dist(M::CrossEntropyPolicySearch, π, U) samples are drawn and the top
p, m, m_elite, k_max = M.p, M.m, M.m_elite, M.k_max m_elite are used to refit the dis-
for k in 1:k_max tribution. The algorithm termi-
θs = rand(p, m) nates after k_max iterations. The
us = [U(π, θs[:,i]) for i in 1:m] distribution p can be defined us-
θ_elite = θs[:,sortperm(us)[(m-m_elite+1):m]] ing the [Link] pack-
p = [Link](typeof(p), θ_elite) age. For example, we might define
end μ = [0.0,1.0]
return p
Σ = [1.0 0.0; 0.0 1.0]
end
p = MvNormal(μ,Σ)

function optimize(M, π, U)
return [Link](optimize_dist(M, π, U))
end

Figure 10.5. The cross entropy

5
method applied to the simple regu-
lator problem using a multivariate
3 Gaussian search distribution. The
θ2

five elite samples in each iteration

1 are shown in red. The initial distri-
bution is set to N ([0, 3], 2I).
−2 0 2 −2 0 2 −2 0 2 −2 0 2
θ1 θ1 θ1 θ1

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
 10 . 5 . e vol u ti on stra te g i e s 221

Z
∇ψ E [U (θ )] = ∇ψ U (θ ) p(θ | ψ) dθ (10.6)
θ ∼ p(·|ψ)
Z
= U (θ )∇ψ p(θ | ψ) dθ (10.7)
p ( θ | ψ)
Z
= U (θ )∇ψ p(θ | ψ) dθ (10.8)
p ( θ | ψ)
Z
(10.9)


= U (θ )∇ψ log p(θ | ψ) p(θ | ψ) dθ

(10.10)
 
= E U (θ )∇ψ log p(θ | ψ)
θ ∼ p(·|ψ)

The introduction of the logarithm above comes from what is called the log deriva-
tive trick, which observes that ∇ψ log p(θ | ψ) = ∇ψ p(θ | ψ)/p(θ | ψ). This
computation requires knowing ∇ψ log p(θ | ψ), but we can often compute this
analytically, as discussed in example 10.2.
The search gradient can be estimated from m samples: θ (1) , . . . , θ (m) ∼ p(· | ψ):
m
1
∇ψ E [U (θ )] ≈ ∑ U (θ (i) )∇ψ log p(θ (i) | ψ) (10.11)
θ ∼ p(·|ψ) m i =1

This estimate depends on the evaluated expected utility, which itself can vary
widely. We can make our gradient estimate more resilient with rank shaping,
which replaces the utility values with weights based on the relative performance
of each sample to the other samples in its iteration. The m samples are sorted in
descending order of expected utility. Weight w(i) is assigned to sample i according
to a weighting scheme with w(1) ≥ · · · ≥ w(m) . The search gradient becomes
m
∇ψ E [U (θ )] ≈ ∑ w(i) ∇ψ log p(θ (i) | ψ) (10.12)
θ ∼ p(·|ψ) i =1

A common weighting scheme is12 12

N. Hansen and A. Ostermeier,
“Adapting Arbitrary Normal Mu-
max 0, log m2 + 1 − log(i ) tation Distributions in Evolution



1
w (i )
= m m


− (10.13) Strategies: The Covariance Matrix
∑ j=1 max 0, log 2 + 1 − log( j) m Adaptation,” in IEEE International
Conference on Evolutionary Computa-
These weights, shown in figure 10.6, favor better samples and give most samples tion, 1996.
a small negative weight. Rank-shaping reduces the influence of outliers.
Algorithm 10.5 provides an implementation of the evolution strategies method.
Figure 10.7 shows an example of a search progression.

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
222 c ha p te r 10 . p ol i c y se a rc h

The multivariate normal distribution N (µ, Σ), with mean µ and covariance Example 10.2. A derivation of the
log likelihood gradient equations
Σ, is a common distribution family. The likelihood in d dimensions takes the for the multivariate Gaussian dis-
form tribution. For the original deriva-
tion and several more sophisticated
solutions for handling the positive
 
d 1 1
p(x | µ, Σ) = (2π )− 2 |Σ|− 2 exp − (x − µ)⊤ Σ−1 (x − µ) definite covariance matrix, see D.
2
Wierstra, T. Schaul, T. Glasmach-
ers, Y. Sun, J. Peters, and J. Schmid-
where |Σ| is the determinant of Σ. The log likelihood is huber, “Natural Evolution Strate-
gies,” Journal of Machine Learning
d 1 1 Research, vol. 15, pp. 949–980, 2014.
log p(x | µ, Σ) = − log(2π ) − log |Σ| − (x − µ)⊤ Σ−1 (x − µ)
2 2 2
The parameters can be updated using their log likelihood gradients:

∇µ log p(x | µ, Σ) = Σ−1 (x − µ)

1 1
∇Σ log p(x | µ, Σ) = Σ−1 (x − µ)(x − µ)⊤ Σ−1 − Σ−1
2 2
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 product
A⊤ A, which guarantees that Σ remains positive semidefinite, and then to
update A instead. Replacing Σ by A⊤ A and taking the gradient with respect
to A yields
h i
∇(A) log p(x | µ, A) = A ∇Σ log p(x | µ, Σ) + ∇Σ log p(x | µ, Σ)⊤

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
 10 . 5 . e vol u ti on stra te g i e s 223

m =2 Figure 10.6. Several weight-

0.5
ings constructed using equa-
m =4
tion (10.13).
m =7
m = 10
weight, w(i)

−0.5

1 2 3 4 5 6 7 8 9 10
index in order by decreasing expected utility, i

struct EvolutionStrategies Algorithm 10.5. An evolution

D # distribution constructor strategies method for updating a
ψ # initial distribution parameterization search distribution D(ψ) over pol-
∇logp # log search likelihood gradient icy parameterizations for policy
m # number of samples π(θ, s). This implementation also
α # step factor takes an initial search distribution
k_max # number of iterations parameterization ψ, the log search
end likelihood gradient ∇logp(ψ, θ),
a policy evaluation function U,
function evolution_strategy_weights(m) and an iteration count k_max. In
ws = [max(0, log(m/2+1) - log(i)) for i in 1:m] each iteration, m parameteriza-
ws ./= sum(ws) tion samples are drawn and are
ws .-= 1/m used to estimate the search gradi-
return ws ent. This gradient is then applied
end with a step factor α. We can use
[Link] to define D(ψ).
function optimize_dist(M::EvolutionStrategies, π, U)
For example, if we want to define
D, ψ, m, ∇logp, α = M.D, M.ψ, M.m, M.∇logp, M.α
D to construct a Gaussian with a
ws = evolution_strategy_weights(m)
given mean ψ and fixed covariance
for k in 1:M.k_max
Σ, we can use
θs = rand(D(ψ), m)
D(ψ) = MvNormal(ψ, Σ).
us = [U(π, θs[:,i]) for i in 1:m]
sp = sortperm(us, rev=true)
∇ = sum(w.*∇logp(ψ, θs[:,i]) for (w,i) in zip(ws,sp))
ψ += α.*∇
end
return D(ψ)
end

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
224 c ha p te r 10 . p ol i c y se a rc h

Figure 10.7. Evolution strategies

5
(algorithm 10.5) applied to the sim-
ple regulator problem using a mul-
3 tivariate Gaussian search distribu-
θ2

tion. Samples are shown in white,

1 along with their search gradient
contributions, w∇ log p.
−2 0 2 −2 0 2 −2 0 2 −2 0 2
θ1 θ1 θ1 θ1

10.6 Isotropic Evolutionary Strategies

The previous section introduced evolutionary strategies that can work with gen- 13
An example of this approach ap-
eral search distributions. This section will make the assumption that the search plied to policy search is explored
distribution is a spherical or isotropic Gaussian, where the covariance matrix takes by T. Salimans, J. Ho, X. Chen, S.
Sidor, and I. Sutskever, “Evolution
the form σ2 I.13 Under this assumption, the expected utility of the distribution Strategies as a Scalable Alternative
introduced in equation (10.5) simplifies to14 to Reinforcement Learning,” 2017.
arXiv: 1703.03864v2.
E [U (θ)] = E [U (ψ + σǫ)] (10.14) 14
In general, if A⊤ A = Σ, then θ =
θ∼N (ψ,σ2 I) ǫ∼N (0,I) µ + A⊤ ǫ transforms ǫ ∼ N (0, I)
into a sample θ ∼ N (µ, Σ).
The search gradient reduces to

(10.15)
h i
∇ψ E [U (θ)] = E U (θ)∇ψ log p(θ | ψ, σ2 I)
θ∼N (ψ,σ2 I) θ∼N (ψ,σ2 I)
  15
D. Brockhoff, A. Auger, N.
1 Hansen, D. Arnold, and T.
= E U (θ) 2 (θ − ψ) (10.16)
2
θ∼N (ψ,σ I) σ Hohm, “Mirrored Sampling and
  Sequential Selection for Evolution
1 Strategies,” in International Confer-
= E U (ψ + σǫ) 2 (σǫ) (10.17)
ǫ∼N (0,I) σ ence on Parallel Problem Solving from
Nature, 2010.
1
= E [U (ψ + σǫ)ǫ] (10.18)
σ ǫ∼N (0,I) 16
This technique was implemented
by T. Salimans, J. Ho, X. Chen, S.
Algorithm 10.6 provides an implementation of this strategy. This implemen-
Sidor, and I. Sutskever, “Evolution
tation incorporates mirrored sampling.15 We sample m/2 values from the search Strategies as a Scalable Alterna-
distribution and then generate the other m/2 samples by mirroring them about tive to Reinforcement Learning,”
2017. arXiv: 1703 . 03864v2. They
the mean. Mirrored samples reduce the variance of the gradient estimate.16 The included other techniques as well,
benefit of using this technique is shown in figure 10.8. including weight decay.

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
 10 . 6 . i sotrop i c e vol u ti on a ry stra te g i e s 225

struct IsotropicEvolutionStrategies Algorithm 10.6. An evolution

ψ # initial mean strategies method for updating
σ # initial standard deviation an isotropic multivariate Gaus-
m # number of samples sian search distribution with mean
α # step factor ψ and covariance σ2 I over pol-
k_max # number of iterations icy parameterizations for a policy
end π(θ, s). This implementation also
takes a policy evaluation function
function optimize_dist(M::IsotropicEvolutionStrategies, π, U) U, a step factor α, and an itera-
ψ, σ, m, α, k_max = M.ψ, M.σ, M.m, M.α, M.k_max tion count k_max. In each iteration,
n = length(ψ) m/2 parameterization samples are
ws = evolution_strategy_weights(2*div(m,2)) drawn and mirrored and are then
for k in 1:k_max used to estimate the search gradi-
ϵs = [randn(n) for i in 1:div(m,2)] ent.
append!(ϵs, -ϵs) # weight mirroring
us = [U(π, ψ + σ.*ϵ) for ϵ in ϵs]
sp = sortperm(us, rev=true)
∇ = sum(w.*ϵs[i] for (w,i) in zip(ws,sp)) / σ
ψ += α.*∇
end
return MvNormal(ψ, σ)
end

800 Figure 10.8. A demonstration of

the effect that mirrored sampling
has on isotropic evolution strate-
600 gies. Two-layer neural network
expected utility

policies were trained on the cart-

pole problem (appendix F.3) us-
400 ing m = 10, and σ = 0.25, with six
rollouts per evaluation. Mirrored
200 sampling significantly speeds and
stabilizes learning.
without mirroring
0 with mirroring

0 10 20 30 40 50 60
iteration

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
226 c ha p te r 10 . p ol i c y se a rc h

10.7 Summary

• Monte Carlo policy evaluation involves computing the expected utility associ-
ated with a policy using a large number of rollouts from states sampled from
an initial state distribution.

• Local search methods, such as the Hooke-Jeeves method, improve a policy

based on small, local changes.

• Genetic algorithms maintain a population of points in the parameter space,

recombining them in different ways in attempt to drive the population toward
a global optimum.

• The cross entropy method iteratively improves a search distribution over policy
parameters by refitting the distribution to elite samples at each iteration.

• Evolutionary strategies attempt to improve the search distribution using gradi-

ent information from samples from that distribution.

• Isotropic evolutionary strategies make the assumption that the search distribu-
tion is an isotropic Gaussian.

10.8 Exercises
Exercise 10.1. In Monte Carlo policy evaluation, how is the variance of the utility estimate
affected by the number of samples?

Solution: The variance of Monte Carlo policy evaluation is the variance of the mean of m
samples. These samples are assumed to be independent, and so the variance of the mean
is the variance of a single rollout evaluation divided by the sample size:
" #
m
(i ) 1
Var[Û (π )] = Var ∑ R(τ ) = Varτ [ R(τ )]
i =1
m

where Û (π ) is the utility from Monte Carlo policy evaluation and R(τ ) is the trajectory
reward for a sampled trajectory τ. The sample variance, therefore, decreases with 1/m.

Exercise 10.2. What effect does varying the number of samples m and the number of elite
samples melite have on cross entropy policy search?

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
 10 . 8 . e x e rc i se s 227

Solution: The computational cost per iteration scales linearly with the number of samples.
More samples will better cover the search space, resulting in a better chance of identifying
better elite samples to improve the policy. The number of elite samples also has an effect.
Making all samples elite provides no feedback to the improvement process. Having too
few elite samples can lead to premature convergence to a suboptimal solution.

Exercise 10.3. Consider using evolution strategies with a univariate Gaussian distribution,
θ ∼ N (µ, ν). What is the search gradient with respect to the variance ν? What issue arises
as the variance becomes small?

Solution: The search gradient is the gradient of the log-likelihood:

!
∂ ∂ 1 ( x − µ )2
log p( x | µ, ν) = log √ exp −
∂ν ∂ν 2πν 2ν
!
∂ 1 1 ( x − µ )2
= − log(2π ) − log(ν) −
∂ν 2 2 2ν

1 ( x − µ )2
=− +
2ν 2ν2
log p(1 | 0, ν)

0.4

0.2

0
∂ν
∂

1 2 3 4 5
ν

We find that the gradient goes to infinity as the variance approaches zero. This is a
problem because the variance should be small when the search distribution converges.
Very large gradients can cause simple ascent methods to overshoot optima.

Exercise 10.4. Equation (10.14) defines the objective in terms of a search distribution
θ ∼ N (ψ, Σ). What advantage does this objective have over directly optimizing θ using
the expected utility objective in equation (10.1)?

Solution: The added Gaussian noise around the policy parameters can smooth discontinu-
ities in the original objective, which can make optimization more reliable.

Exercise 10.5. Which of the methods in this chapter are best suited to the fact that multiple
types of policies could perform well in a given problem?

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
228 c ha p te r 10 . p ol i c y se a rc h

Solution: The Hooke-Jeeves method improves a single policy parameterization, so it cannot

retain multiple policies. Both the cross entropy method and evolution strategies use search
distributions. In order to successfully represent multiple types of policies, a multimodal
distribution would have to be used. One common multimodal distribution is a mixture of
Gaussians. A mixture of Gaussians cannot be fit analytically, but they can be reliably fit
using expectation maximization (EM), as demonstrated in example 4.4. Genetic algorithms
can retain multiple policies if the population size is sufficiently large.

Exercise 10.6. Suppose we have a parameterized policy πθ that we would like to optimize
using the Hooke-Jeeves method. If we initialize our parameter θ = 0 and the utility function
is U (θ ) = −3θ 2 + 4θ + 1, what is the largest step size α that would still guarantee policy
improvement in the first iteration of the Hooke-Jeeves method?

Solution: The Hooke-Jeeves method evaluates the objective function at the center point ±α
along each coordinate direction. In order to guarantee improvement in the first iteration of
Hooke-Jeeves search, at least one of the objective function values at the new points must
improve the objective function value. For our policy optimization problem, this means
that we are searching for the largest step size α such that either U (θ + α) or U (θ − α) is
greater than U (θ ).
Since the underlying utility function is parabolic and concave, the largest step size that
would still lead to improvement is slightly less than the width of the parabola at the current
point. Thus, we compute the point on the parabola opposite the current point, θ ′ at which
U ( θ ′ ) = U ( θ ):

U (θ ) = −3θ 2 + 4θ + 1 = −3(0)2 + 4(0) + 1 = 1

U (θ ) = U (θ ′ )
1 = −3θ ′2 + 4θ ′ + 1
0 = −3θ ′2 + 4θ ′ + 0
p
′ −4 ± 42 − 4(−3)(0) −4 ± 4 2±2 n o
θ = = = = 0, 34
2(−3) −6 3

The point on the parabola opposite the current point is thus θ ′ = 34 . The distance be-
tween θ and θ ′ is 43 − 0 = 43 . Thus, the maximal step size we can take and still guarantee
improvement in the first iteration is just under 43 .

Exercise 10.7. Suppose we have a policy parameterized by a single parameter θ. We

take an evolution strategies approach with a search distribution that follows a Bernoulli
distribution p(θ | ψ) = ψθ (1 − ψ)1−θ . Compute the log-likelihood gradient ∇ψ log p(θ |
ψ ).

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]
 10 . 8 . e x e rc i se s 229

Solution: The log-likelihood gradient can be computed as follows:

p ( θ | ψ ) = ψ θ (1 − ψ )1− θ
 
log p(θ | ψ) = log ψθ (1 − ψ)1−θ

log p(θ | ψ) = θ log ψ + (1 − θ ) log(1 − ψ)

d
∇ψ log p(θ | ψ) = [θ log ψ + (1 − θ ) log(1 − ψ)]
dψ
θ 1−θ
∇ψ log p(θ | ψ) = −
ψ 1−ψ

Exercise 10.8. Compute the sample weights for search gradient estimation with rank
shaping given m = 3 samples.

Solution: We first compute the numerator of the first term from equation (10.13), for all i:
   
i=1 max 0, log 32 + 1 − log 1 = log 25
   
i=2 max 0, log 23 + 1 − log 2 = log 45
   
i=3 max 0, log 23 + 1 − log 3 = 0

Now, we compute the weights:

log 25 1
w (1) = 5 5
− = 0.47
log 2 + log + 04
3
log 45 1
w (2) = 5
− = −0.14
log + log 45 + 0 3
2
0 1
w (3) = 5 5
− = −0.33
log 2 + log 4 + 0 3

© 2022 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-07 21:55:16-07:00, comments to bugs@[Link]