20 Expression Optimization
Previous chapters discussed optimization over a fixed set of design variables. For
many problems, the number of variables is unknown, such as in the optimizationof graphical structures or computer programs. Designs in these contexts can be
represented by expressions that belong to a grammar. This chapter discusses
ways to make the search of optimal designs more efficient by accounting for thegrammatical structure of the design space.

20.1 Grammars
An expression can be represented by a tree of symbols. For example, the mathematical expression x + ln 2 can be represented using the tree in figure 20.1 consisting of
the symbols +, x, ln, and 2. Grammars specify constraints on the space of possible
expressions.
A grammar is represented by a set of production rules. These rules involvesymbols as well as types. A type can be interpreted as a set of expression trees. A
production rule represents a possible expansion of type to an expression involving
symbols or types. If a rule expands only to symbols, then it is called terminal
because it cannot be expanded further. An example of a nonterminal rule is
R 7→ R + R, which means that the type R can consist of elements of the set R
added to elements in the set R.1
We can generate an expression from a grammar by starting with a start type
and then recursively applying different production rules. We stop when thetree contains only symbols. Figure 20.2 illustrates this process for the expression
x + ln 2. An application to natural language expressions is shown in example 20.1.
The number of possible expressions allowed by a grammar can be infinite.
Example 20.2 shows a grammar that allows for infinitely many valid expressions.

+
x

ln

2
Figure 20.1. The expression x +ln 2 represented as a tree.

This chapter focuses on contextfree grammars, but other forms exist.
See L. Kallmeyer, Parsing Beyond
Context-Free Grammars. Springer,2010.
1362 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

+

R
R

+
x

R

Figure 20.2. Using the production
rules
R 7→ R + R

+
ln

x

R 7→ xln

R 7→ ln(R )R

2Consider a grammar that allows for the generation of simple English statements:
S 7→ N V
V 7→ V A

A 7→ rapidly | efficiently

N 7→ Alice | Bob | Mykel | Tim
V 7→ runs | reads | writes

R 7→ 2

to generate x + ln 2. Blue nodes are
unexpanded types.

Example 20.1. A grammar for producing simple English statements.
Using | on the right-hand side ofan expression is shorthand for or.
Thus, the rule
A 7→ rapidly | efficiently

is equivalent to having two rules,
A 7→ rapidly and A 7→ efficiently.

The types S, N, V, and A correspond to statements, nouns, verbs, andadverbs, respectively. An expression is generated by starting with the type S
and iteratively replacing types:
S
NV
Mykel V A
Mykel writes rapidly
Not all terminal symbol categories must be used. For instance, the statement‘‘Alice runs’’ can also be generated.

© 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

20 . 1. g ra mma rs

363

Expression optimization often constrains expressions to a maximum depth or
penalizes expressions based on their depth or node count. Even if the grammar
allows a finite number of expressions, the space is often too vast to search exhaustively. Hence, there is a need for algorithms that efficiently search the space ofpossible expressions for one that optimizes an objective function.
Consider a four-function calculator grammar that applies addition, subtraction, multiplication, and division to the ten digits:
R 7→ R + R

Example 20.2. Some of the challenges associated with grammars,
as illustrated with a four-function
calculator grammar.

R 7→ R − R
R 7→ R × R
R 7→ R / R

R 7→ 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
An infinite number of expressions can be generated because the nonterminal
R can always be expanded into one of the calculator operations.
Many expressions will produce the same value. Addition and multiplication operators are commutative, meaning that the order does not matter.
For example, a + b is the same as b + a. These operations are also associative,
meaning the order in which multiple operations of the same type occur do
not matter. For example, a × b × c is the same as c × b × a. Other operations
preserve values, like adding zero or multiplying by one.
Not all expressions under this grammar are mathematically valid. For
example, division by zero is undefined. Removing zero as a terminal symbol
is insufficient to prevent this error because zero can be constructed using
other operations, such as 1 − 1. Such exceptions are often handled by theobjective function, which can catch exceptions and penalize them.
The expression optimization routines covered in the chapter use ExprRules.jl.
Grammars can be defined using the grammar macro by listing the production
rules, as shown in example 20.3.
Many of the expression optimization algorithms involve manipulating components of an expression tree in a way that preserves the way the types wereexpanded. A RuleNode object tracks which production rules were applied when
© 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

364 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

We may define a grammar using the grammar macro. The nonterminals areon the left of the equal sign, and the expressions with terminals and nonterminals are the on the right. The package includes some syntax to represent
grammars more compactly.

Example 20.3. Example ofdefining a grammar using the
ExprRules.jl package.

using ExprRules
grammar = @grammar begin
R = x# reference a variable
R = R * A# multiple children
R = f(R)# call a function
R = _(randn()) # random variable generated on node creation
R = 1 | 2 | 3 # equivalent to R = 1, R = 2, and R = 3
R = |(4:6)# equivalent to R = 4, R = 5, and R = 6
A = 7# rules for different return types
end;doing an expansion. Calling rand with a specified starting type will generate
a random expression represented by a RuleNode. Calling sample will select arandom RuleNode from an existing RuleNode tree. Nodes are evaluated usingeval.
The method return_type returns the node’s return type as a symbol, isterminal
returns whether the symbol is terminal, child_types returns the list of nonterminal symbols associated with the node’s production rule, and nchildren returns
the number of children. These four methods each take as input the grammar
and the node. The number of nodes in an expression tree is obtained using
length(node), and the depth is obtained using depth(node).
A third type, NodeLoc, is used to refer to a node’s location in the expression
tree. Subtrees manipulation often requires NodeLocs.
loc = sample(NodeLoc, node); # uniformly sample a node loc
loc = sample(NodeLoc, node, :R, grammar); # sample a node loc of type R
subtree = get(node, loc);

20.2 Genetic Programming
Genetic algorithms (see chapter 9) use chromosomes that encode design points
in a sequential format. Genetic programming2 represents individuals using trees

2
J. R. Koza, Genetic Programming:
On the Programming of Computers
by Means of Natural Selection. MIT
Press, 1992.

© 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

20 . 2. g e n e ti c p rog ra mmi n g 365

instead (figure 20.3), which are better at representing mathematical functions,programs, decision trees, and other hierarchical structures.
Figure 20.3. A tree representationof the Julia method:

if

×

≤
x

2

x ≤ 2 ? x*(y+3) : 0

0

+

x
y

3

Similar to genetic algorithms, genetic programs are initialized randomly andsupport crossover and mutation. In tree crossover (figure 20.4), two parent trees
are mixed to form a child tree. A random node is chosen in each parent, and thesubtree at the chosen node in the first parent is replaced with the subtree at thechosen node of the second parent. Tree crossover works on parents with different
sizes and shapes, allowing arbitrary trees to mix. In some cases one must ensurethat replacement nodes have certain types, such as Boolean values input into thecondition of an if statement.3 Tree crossover is implemented in algorithm 20.1.

This book focuses only on genetic operations that adhere to
the constraints of the grammar.
Sometimes genetic programmingwith this restriction is referred toas strongly typed genetic programming, as discussed in D. J. Montana, “Strongly Typed Genetic Programming,” Evolutionary Computation, vol. 3, no. 2, pp. 199–230, 1995.
3

Figure 20.4. Tree crossover is used
to combine two parent trees to produce a child tree.

Parent A

Parent B

Child

Tree crossover tends to produce trees with greater depth than the parent trees.
Each generation tends to increase in complexity, which often results in overly
complicated solutions and slower runtimes. We encourage parsimony, or simplicity,in the solution, by introducing a small bias in the objective function value basedon a tree’s depth or node count.
Applying tree mutation (figure 20.5) starts by choosing a random node in thetree. The subtree rooted at that node is deleted, and a new random subtree isgenerated to replace the old subtree. In contrast to mutation in binary chromo© 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

366 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

struct TreeCrossover <: CrossoverMethod
grammar
max_depthendfunction crossover(C::TreeCrossover, a, b)
child = deepcopy(a)crosspoint = sample(b)typ = return_type(C.grammar, crosspoint.ind)
d_subtree = depth(crosspoint)
d_max = C.max_depth + 1 - d_subtree
if d_max > 0 && contains_returntype(child,C.grammar,typ,d_max)
loc = sample(NodeLoc, child, typ, C.grammar, d_max)
insert!(child, loc, deepcopy(crosspoint))
end
child
end

Algorithm 20.1. Tree crossover
implemented for a and b of type
RuleNode from ExprRules.jl. The
TreeCrossover struct contains arule set grammar and a maximum
depth max_depth.

somes, tree mutation can typically occur at most once, often with a low probability
around 1%. Tree mutation is implemented in algorithm 20.2.
Figure 20.5. Tree mutation deletes
a random subtree and generates a
new one to replace it.

BeforeAfter© 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

20 . 2. g e n e ti c p rog ra mmi n g 367

struct TreeMutation <: MutationMethod
grammarpendfunction mutate(M::TreeMutation, a)
child = deepcopy(a)if rand() < M.p
loc = sample(NodeLoc, child)
typ = return_type(M.grammar, get(child, loc).ind)
subtree = rand(RuleNode, M.grammar, typ)insert!(child, loc, subtree)
end
return child
end

Algorithm 20.2. Tree mutation implemented for an individual a oftype RuleNode from ExprRules.jl.
The TreeMutation struct contains
a rule set grammar and a mutation
probability p.

Tree permutation (figure 20.6) is a second form of genetic mutation. Here, thechildren of a randomly chosen node are randomly permuted. Tree permutation
alone is typically not sufficient to introduce new genetic material into a populationand is often combined with tree mutation. Tree permutation is implemented inalgorithm 20.3.
Figure 20.6. Tree permutation permutes the children of a randomly
chosen node.

BeforeAfter

The implementation of genetic programming is otherwise identical to thatof genetic algorithms. More care must typically be taken in implementing thecrossover and mutation routines, particularly when determining what sorts ofnodes can be generated and that only syntactically correct trees are produced.
Genetic programming is used to generate an expression that approximates π in
example 20.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

368 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

struct TreePermutation <: MutationMethod
grammarpendfunction mutate(M::TreePermutation, a)
child = deepcopy(a)if rand() < M.p
node = sample(child)n = length(node.children)
types = child_types(M.grammar, node)
for i in 1 : n-1
c = 1
for k in i+1 : n
if types[k] == types[i] &&rand() < 1/(c+=1)

end

end

end

Algorithm 20.3. Tree permutation
implemented for an individual a oftype RuleNode from ExprRules.jl,
where p is the mutation probability.

node.children[i], node.children[k] =node.children[k], node.children[i]endend
return child© 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

20 . 2. g e n e ti c p rog ra mmi n g 369

Consider approximating π using only operations on a four-function calculator. We can solve this problem using genetic programming where nodes canbe any of the elementary operations: add, subtract, multiply, and divide, andthe digits 1 − 9.
We use ExprRules.jl to specify our grammar:
grammar = @grammar begin
R = |(1:9)
R = R + R
R = R - R
R = R / R
R = R * Rend

We construct an objective function and penalize large trees:
function f(node)
value = eval(node, grammar)
if isinf(value) || isnan(value)
return Inf
end
Δ = abs(value - π)return log(Δ) + length(node)/1e3
end

Example 20.4. Using genetic programming to estimate π using only
digits and the four principle arithmetic operations.
/
+

4
+

3
/

9
+

7

+

7
/

5
+

3
/
3

8
9

We finally run our genetic program using a call to the genetic_algorithm
function from section 9.2:
population = [rand(RuleNode, grammar, :R) for i in 1:1000]
best_tree = genetic_algorithm(f, population, 30,
TruncationSelection(50),
TreeCrossover(grammar, 10),
TreeMutation(grammar, 0.25))

The best performing tree is shown on the right. It evaluates to 3.141586, whichmatches π to four decimal places.

© 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

370 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

20.3

Grammatical Evolution

Grammatical evolution4 operates on an integer array instead of a tree, allowing thesame techniques employed in genetic algorithms to be applied. Unlike genetic
algorithms, the chromosomes in grammatical evolution encode expressions basedon a grammar. Grammatical evolution was inspired by genetic material, which isinherently serial like the chromosomes used in genetic algorithms. 5
In grammatical evolution, designs are integer arrays much like the chromosomes used in genetic algorithms. Each integer is unbounded because indexing
is performed using modular arithmetic. The integer array can be translated into
an expression tree by parsing it from left to right.
We begin with a starting symbol and a grammar. Suppose n rules in thegrammar can be applied to the starting symbol. The jth rule is applied, wherej = i mod1 n and i is the first integer in the integer array.6
We then consider the rules applicable to the resulting expression and use
similar modular arithmetic based on the second integer in the array to selectwhich rule to apply. This process is repeated until no rules can be applied and thephenotype is complete.7 The decoding process is implemented in algorithm 20.4
and is worked through in example 20.5.
It is possible for the integer array to be too short, thereby causing the translation
process to run past the length of the array. Rather than producing an invalid
individual and penalizing it in the objective function, the process wraps aroundto the beginning of the array instead. This wrap-around effect means that the samedecision can be read several times during the transcription process. Transcription
can result in infinite loops, which can be prevented by a maximum depth.
Genetic operations work directly on the integer design array. We can adopt
the operations used on real-valued chromosomes and apply them to the integervalued chromosomes. The only change is that mutation must preserve real values.
A mutation method for integer-valued chromosomes using zero-mean Gaussian
perturbations is implemented in algorithm 20.5.
Grammatical evolution uses two additional genetic operators. The first, gene
duplication, occurs naturally as an error in DNA replication and repair. Geneduplication can allow new genetic material to be generated and can store a second
copy of a useful gene to reduce the chance of a lethal mutation removing the genefrom the gene pool. Gene duplication chooses a random interval of genes in the

C. Ryan, J. J. Collins, and M. O.
Neill, “Grammatical Evolution:
Evolving Programs for an Arbitrary Language,” in European
Conference on Genetic Programming,
1998.
5
Our serial DNA is read and used
to construct complicated protein
structures. DNA is often referred
to as the genotype—the object
on which genetic operations are
performed. The protein structureis the phenotype—the object encoded by the genotype whose performance is evaluated. The grammatical evolution literature often
refers to the integer design vector
as the genotype and the resulting
expression as the phenotype.
4

We use x mod1 n to refer to the 1index modulus:

6(( x − 1) mod n) + 1

This type of modulus is useful with
1-based indexing. The corresponding Julia function is mod1.
7
No genetic information is readwhen there is only a single applicable rule.

© 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

20 . 3. g ra mma ti c a l e vol u ti on

struct DecodedExpression
node
n_rules_applied
end
function decode(x, grammar, sym, c_max=1000, c=0)
node, c = _decode(x, grammar, sym, c_max, c)
DecodedExpression(node, c)endfunction _decode(x, grammar, typ, c_max, c)
types = grammar[typ]if length(types) > 1
g = x[mod1(c+=1, length(x))]
rule = types[mod1(g, length(types))]else
rule = types[1]endnode = RuleNode(rule)childtypes = child_types(grammar, node)
if !isempty(childtypes) && c < c_max
for ctyp in childtypes
cnode, c = _decode(x, grammar, ctyp, c_max, c)
push!(node.children, cnode)end
end
return (node, c)
end

371

Algorithm 20.4. A method for decoding an integer design vector
to produce an expression, where
x is a vector of integers, grammar
is a Grammar, and sym is the root
symbol. The counter c is usedduring the recursion process andthe parameter c_max is an upper
limit on the maximum number ofrule applications, to prevent an infinite loop. The method returns
a DecodedExpression, which contains the expression tree and thenumber of rules applied during thedecoding process.

© 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

372 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

Consider a grammar for real-valued strings:
R 7→ D D ′ P E

D ′ 7→ D D ′ | ǫ

P 7→ . D D ′ | ǫ

E 7→ E S D D ′ | ǫ
S 7→ + | − | ǫ

Example 20.5. The process bywhich an integer design vector ingrammatical evolution is decoded
into an expression.
Our implementation is depthfirst. If 52 were instead 51, the ruleD ′ 7→ D D ′ would be applied,followed by selecting a rule forthe new D, eventually resulting in
43950.950E+8.

D 7→ 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9
where R is a real value, D is a terminal decimal, D ′ is a nonterminal decimal,
P is the decimal part, E is the exponent, and S is the sign. Any ǫ values
produce empty strings.
Suppose our design is [205, 52, 4, 27, 10, 59, 6] and we have the starting
symbol R. There is only one applicable rule, so we do not use any genetic
information and we replace R with DD ′ PE.
Next we must replace D. There are 10 options. We select 205 mod1 10 = 5,
and thus obtain 4D ′ PE
Next we replace D ′ , which has two options. We select index 52 mod1 2 = 2,
which corresponds to ǫ.
Continuing in this manner we produce the string 4E+8.
The above grammar can be implemented in ExprRules using:grammar = @grammar begin
R = D * De * P * E
De = D * De | ""P = "." * D * De | ""E = "E" * S * D * De | ""S = "+"|"-"|""
D = "0"|"1"|"2"|"3"|"4"|"5"|"6"|"7"|"8"|"9"end

and can be evaluated using:x = [205, 52, 4, 27, 10, 59, 6]
str = eval(decode(x, grammar, :R).node, grammar)© 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

20 . 3. g ra mma ti c a l e vol u ti on

struct IntegerGaussianMutation <: MutationMethodσendfunction mutate(M::IntegerGaussianMutation, child)return child + round.(Int, randn(length(child)).*M.σ)
end

373

Algorithm 20.5. The Gaussian mutation method modified to preserve integer values for integervalued chromosomes. Each value
is perturbed by a zero-mean Gaussian random value with standard
deviation σ and then rounded to
the nearest integer.

chromosome to duplicate. A copy of the selected interval is appended to the backof the chromosome. Duplication is implemented in algorithm 20.6.
struct GeneDuplication <: MutationMethod
endfunction mutate(M::GeneDuplication, child)n = length(child)i, j = rand(1:n), rand(1:n)
interval = min(i,j) : max(i,j)
return vcat(child, deepcopy(child[interval]))end

The second genetic operation, pruning, tackles a problem encountered duringcrossover. As illustrated in figure 20.7, crossover will select a crossover point at
random in each chromosome and construct a new chromosome using the left
side of the first and the right side of the second chromosome. Unlike genetic
algorithms, the trailing entries in chromosomes of grammatical evolution may
not be used; during parsing, once the tree is complete, the remaining entries
are ignored. The more unused entries, the more likely it is that the crossover
point lies in the inactive region, thus not providing new beneficial material. An
individual is pruned with a specified probability, and, if pruned, its chromosome
is truncated such that only active genes are retained. Pruning is implemented inalgorithm 20.7 and is visualized in figure 20.8.
Like genetic programming, grammatical evolution can use the genetic algorithm method.8 We can construct a MutationMethod that applies multiple mutation methods in order to use pruning, duplication, and standard mutation
approaches in tandem. Such a method is implemented in algorithm 20.8.
Grammatical evolution suffers from two primary drawbacks. First, it is difficult
to tell whether the chromosome is feasible without decoding it into an expression.

Algorithm 20.6. The gene duplication method used in grammaticalevolution.

8
Genotype to phenotype mapping
would occur in the objective 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

374 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

parent aparent b
child

struct GenePruning <: MutationMethodpgrammartypendfunction mutate(M::GenePruning, child)
if rand() < M.p
c = decode(child, M.grammar, M.typ).n_rules_applied
if c < length(child)child = child[1:c]end
end
return childend

beforeafter

struct MultiMutate <: MutationMethodMs
endfunction mutate(M::MultiMutate, child)for m in M.Ms
child = mutate(m, child)end
return childend

Figure 20.7. Crossover applied tochromosomes in grammatical evolution may not affect the active
genes in the front of the chromosome. The child shown here inherits all of the active shaded genes
from parent a so it will effectively
act as an identical expression. Pruning was developed to overcome
this issue.
Algorithm 20.7. The gene pruning
method used in grammatical evolution.

Figure 20.8. Pruning truncates thechromosome such that only active
genes remain.

Algorithm

20.8.
A
for applying
all mutation methods stored in thevector Ms.
MutationMethod© 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

20 . 4 . p roba bi l i sti c g ra mma rs 375

Second, a small change in the chromosome may produce a large change in thecorresponding expression.

20.4 Probabilistic Grammars
A probabilistic grammar9 adds a weight to each rule in the genetic program’s grammar. When sampling from all applicable rules for a given node, we select a rulestochastically according to the relative weights. The probability of an expression is
the product of the probabilities of sampling each rule. Algorithm 20.9 implements
the probability calculation. Example 20.6 demonstrates sampling an expression
from a probabilistic grammar and computes its likelihood.
struct ProbabilisticGrammar
grammarws
endfunction probability(probgram, node)
typ = return_type(probgram.grammar, node)i = findfirst(isequal(node.ind), probgram.grammar[typ])prob = probgram.ws[typ][i] / sum(probgram.ws[typ])for (i,c) in enumerate(node.children)
prob *= probability(probgram, c)
end
return prob
end

Optimization using a probabilistic grammar improves its weights with each
iteration using elite samples from a population. At each iteration, a population of
expressions is sampled and their objective function values are computed. Some
number of the best expressions are considered the elite samples and can be used
to update the weights. A new set of weights is generated for the probabilistic
grammar, where the weight wiT for the ith production rule applicable to return
type T is set to the number of times the production rule was used in generating
the elite samples. This update procedure is implemented in algorithm 20.10.
The probabilistic grammars above can be extended to more complicated probability distributions that consider other factors, such as the depth in the expression
or local dependencies among siblings in subtrees. One approach is to use Bayesian
networks.10

T. L. Booth and R. A. Thompson,“Applying Probability Measures to
Abstract Languages,” IEEE Transactions on Computers, vol. C-22, no. 5,
pp. 442–450, 1973.
9

Algorithm 20.9. A method for computing the probability of an expression based on a probabilistic grammar, where probgram is a probabilistic grammar consisting of agrammar grammar and a mappingof types to weights for all applicable rules ws, and node is a RuleNode
expression.

P. K. Wong, L. Y. Lo, M. L. Wong,and K. S. Leung, “GrammarBased Genetic Programmingwith Bayesian Network,” inIEEE Congress on Evolutionary
Computation (CEC), 2014.
10© 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

376 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

Consider a probabilistic grammar for strings composed entirely of ‘‘a’’s:
A 7→ a A

w1A = 17→ ǫ

w3A = 27→ ǫw2B = 17→ a B a A

B 7→ a B

Example 20.6. Sampling an expression from a probabilistic grammar
and computing the expression’s
likelihood.

w2A = 3
w1B = 4

where we have a set of weights w for each parent type and ǫ is an empty
string.
Suppose we generate an expression starting with the type A. The probability distribution over the three possible rules is:P(A 7→ a A ) = 1/(1 + 3 + 2)

P(A 7→ a B a A ) = 3/(1 + 3 + 2)

P(A 7→ ǫ) = 2/(1 + 3 + 2)

= 1/6
= 1/2
= 1/3

Suppose we sample the second rule and obtain a B a A.
Next we sample a rule to apply to B. The probability distribution over thetwo possible rules is:P(B 7→ a B ) = 4/(4 + 1)

P(B 7→ ǫ) = 1/(4 + 1)

= 4/5
= 1/5

Suppose we sample the second rule and obtain a ǫ a A.
Next we sample a rule to apply to A. Suppose we sample A 7→ ǫ to obtaina ǫ a ǫ, which produces the ‘‘a’’-string ‘‘aa’’. The probability of the sequenceof rules applied to produce ‘‘aa’’ under the probabilistic grammar is:
P (A 7 → a B a A ) P (B 7 → ǫ ) P (A 7 → ǫ ) =

1
1 1 1
· · =
2 5 3
30

Note that this is not the same as the probability of obtaining ‘‘aa’’, as othersequences of production rules could also have produced it.

© 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

20 . 5 . p roba bi l i sti c p rototy p e tre e sfunction _update!(probgram, x)
grammar = probgram.grammar
typ = return_type(grammar, x)
i = findfirst(isequal(x.ind), grammar[typ])probgram.ws[typ][i] += 1
for c in x.children_update!(probgram, c)
end
return probgram
end
function update!(probgram, Xs)
for w in values(probgram.ws)
fill!(w,0)
end
for x in Xs_update!(probgram, x)
end
return probgram
end

20.5

377

Algorithm 20.10. A method for applying a learning update to a probabilistic grammar probgram basedon an elite sample of expressions
Xs.

Probabilistic Prototype Trees

The probabilistic prototype tree11 is a different approach that learns a distribution
for every node in the expression tree. Each node in a probabilistic prototype
tree contains a probability vector representing a categorical distribution overthe grammar’s production rules. The probability vectors are updated to reflectknowledge gained from successive generations of expressions. The maximum
number of children for a node is the maximum number of nonterminals among
rules in the grammar.12
Probability vectors are randomly initialized when a node is created. Random
probability vectors can be drawn from a Dirichlet distribution.13 The original implementation initializes terminals to a scalar value of 0.6 and nonterminals to 0.4.
In order to handle strongly-typed grammars we maintain a probability vector for
applicable rules to each parent type. Algorithm 20.11 defines a node type andimplements this initialization method.
Expressions are sampled using the probability vectors in the probabilistic
prototype tree. A rule in a node is drawn from the categorical distribution definedby the node’s probability vector for the required return type, normalizing theR. Salustowicz and J. Schmidhuber, “Probabilistic Incremental Program Evolution,” Evolutionary Computation, vol. 5, no. 2, pp. 123–141,
1997.
11

The arity of a function is thenumber of arguments. The arityof a grammar rule, which can beviewed as a function, is the number of nonterminals in the rule.
13
The Dirichlet distribution is often used to represent a distribution over discrete distributions. D.
Barber, Bayesian Reasoning and Machine Learning. Cambridge University Press, 2012.
12© 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

378 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

Algorithm 20.11. A probabilistic
prototype tree node type and associated initialization function where
ps is a dictionary mapping a symbol corresponding to a return typeto a probability vector over applicable rules, and children is a list of
PPTNodes. The method get_childwill automatically expand the tree
when attempting to access a nonexistent child.

struct PPTNodepschildrenendfunction PPTNode(grammar;
w_terminal = 0.6,
w_nonterm = 1-w_terminal,)
ps = Dict(typ => normalize!([isterminal(grammar, i) ?
w_terminal : w_nonterm
for i in grammar[typ]], 1)
for typ in nonterminals(grammar))
PPTNode(ps, PPTNode[])end
function get_child(ppt::PPTNode, grammar, i)
if i > length(ppt.children)
push!(ppt.children, PPTNode(grammar))end
return ppt.children[i]endassociated probability vector values to obtain a valid probability distribution. The
tree is traversed in depth-first order. This sampling procedure is implemented inalgorithm 20.12 and visualized in figure 20.9.
Learning can use information either from an entire sampled population orfrom elite samples. Let the best expression in the current generation be xbest andthe best expression found so far be xelite . The node probabilities are updated toincrease the likelihood of generating xbest .14
The probability of generating xbest is the product of the probabilities of choosing
each rule in xbest when traversing through the probabilistic prototype tree. We
compute a target probability for P( xbest ):Ptarget = P( xbest ) + (1 − P( xbest )) · α ·ǫ − yeliteǫ − ybest

The original probabilistic prototype tree implementation will periodically increase the likelihood ofgenerating xelite .

14(20.1)

where α and ǫ are positive constants. The fraction on the right-hand side produces
larger steps toward expressions with better objective function values. The target
probability can be calculated using algorithm 20.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

20 . 5 . p roba bi l i sti c p rototy p e tre e s 379

Algorithm 20.12. A method forsampling an expression from aprobabilistic prototype tree. The
tree is expanded as needed.

function rand(ppt, grammar, typ)rules = grammar[typ]rule_index = sample(rules, Weights(ppt.ps[typ]))
ctypes = child_types(grammar, rule_index)
arr = Vector{RuleNode}(undef, length(ctypes))
node = iseval(grammar, rule_index) ?
RuleNode(rule_index, eval(grammar, rule_index), arr) :RuleNode(rule_index, arr)endfor (i,typ) in enumerate(ctypes)
node.children[i] =rand(get_child(ppt, grammar, i), grammar, typ)end
return node

sample 1st expression
×
ln
p1

x

2
p11
p111

Figure 20.9. A probabilistic prototype tree initially contains only
the root node but expands as additional nodes are needed during
expression generation.

sample 2nd expression
×
p1

−
p12

3

p1

+
x

2

7

p11
p111

p12
p112

function probability(ppt, grammar, expr)
typ = return_type(grammar, expr)
i = findfirst(isequal(expr.ind), grammar[typ])p = ppt.ps[typ][i]for (i,c) in enumerate(expr.children)p *= probability(get_child(ppt, grammar, i), grammar, c)
end
return pendfunction p_target(ppt, grammar, x_best, y_best, y_elite, α, ϵ)p_best = probability(ppt, grammar, x_best)
return p_best + (1-p_best)*α*(ϵ - y_elite)/(ϵ - y_best)
end

p121

p122

Algorithm 20.13. Methods for computing the probability of an expression and the target probability, where ppt is the root nodeof the probabilistic prototype tree,grammar is the grammar, expr andx_best are RuleNode expressions,y_best and y_elite are scalar objective function values, and α andϵ are scalar parameters.

© 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

380 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

The target probability is used to adjust the probability vectors in the probabilistic prototype tree. The probabilities associated with the chosen nodes are
increased iteratively until the target probability is exceeded:







(i )
(i )
(i )
for all i in {1, 2, . . .}(20.2)
P xbest ← P xbest + c · α · 1 − P xbest

where xbest is the ith rule applied in expression xbest and c is a scalar.15
The adapted probability vectors are then renormalized to 1 by downscaling thevalues of all nonincreased vector components proportionally to their current value.
The probability vector p, where the ith component was increased, is adjusted
according to:
1 − pipj ← pj
for j 6= i
(20.3)
kpk1 − p i
(i )

15

A recommended value is c = 0.1.

The learning update is implemented in algorithm 20.14.

function _update!(ppt, grammar, x, c, α)
typ = return_type(grammar, x)
i = findfirst(isequal(x.ind), grammar[typ])p = ppt.ps[typ]p[i] += c*α*(1-p[i])psum = sum(p)for j in 1 : length(p)if j != ip[j] *= (1- p[i])/(psum - p[i])endendfor (pptchild,xchild) in zip(ppt.children, x.children)_update!(pptchild, grammar, xchild, c, α)
end
return ppt
endfunction update!(ppt, grammar, x_best, y_best, y_elite, α, c, ϵ)
p_targ = p_target(ppt, grammar, x_best, y_best, y_elite, α, ϵ)while probability(ppt, grammar, x_best) < p_targ_update!(ppt, grammar, x_best, c, α)
end
return ppt
end

Algorithm 20.14. A method for
applying a learning update to aprobabilistic prototype tree withroot ppt, grammar grammar, best
expression x_best with objectivefunction value y_best, elite objective function value y_elite, learning rate α, learning rate multiplierc, and parameter ϵ.

In addition to population-based learning, probabilistic prototype trees can also
explore the design space via mutations. The tree is mutated to explore the region© 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

20 . 5 . p roba bi l i sti c p rototy p e tre e s 381

around xbest . Let p be a probability vector in a node that was accessed whengenerating xbest . Each component in p is mutated with a probability proportional
to the problem size:
pmutation√(20.4)#p #xbestwhere pmutation is a mutation parameter, #p is the number of components in p,
and #xbest is the number of rules applied in xbest . A component i, selected formutation, is adjusted according to:p i ← p i + β · (1 − p i )(20.5)

where β controls the amount of mutation. Small probabilities undergo larger
mutations than do larger probabilities. All mutated probability vectors must be
renormalized. Mutation is implemented in algorithm 20.15 and visualized infigure 20.10.
function mutate!(ppt, grammar, x_best, p_mutation, β;
sqrtlen = sqrt(length(x_best)),)typ = return_type(grammar, x_best)p = ppt.ps[typ]prob = p_mutation/(length(p)*sqrtlen)
for i in 1 : length(p)if rand() < prob
p[i] += β*(1-p[i])endendnormalize!(p, 1)
for (pptchild,xchild) in zip(ppt.children, x_best.children)
mutate!(pptchild, grammar, xchild, p_mutation, β,
sqrtlen=sqrtlen)
end
return ppt
end

Algorithm 20.15. A method for mutating a probabilistic prototype treewith root ppt, grammar grammar,
best expression x_best, mutation
parameter p_mutation, and mutation rate β.

Finally, subtrees in the probabilistic prototype tree are pruned in order toremove stale portions of the tree. A child node is removed if its parent contains
a probability component above a specified threshold such that, when chosen,
causes the child to be irrelevant. This is always the case for a terminal and may© 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

382 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

probability

10increaseprobability vector entries

normalize

probability vector entries

be the case for a nonterminal. Pruning is implemented in algorithm 20.16 anddemonstrated in example 20.7.

20.6 Summary

probability vector entries
Figure 20.10. Mutating a probability vector in a probabilistic
prototype tree with β = 0.5.
The mutated components wereincreased according to equation (20.5) and the resulting
probability vector was renormalized. Notice how smaller
probabilities
receivegreaterincreases.

• Expression optimization allows for optimizing tree structures that, under agrammar, can express sophisticated programs, structures, and other designs
lacking a fixed size.
• Grammars define the rules used to construct expressions.
• Genetic programming adapts genetic algorithms to perform mutation andcrossover on expression trees.
• Grammatical evolution operates on an integer array that can be decoded into
an expression tree.
• Probabilistic grammars learn which rules are best to generate, and probabilistic
prototype trees learn probabilities for every iteration of the expression rule
generation process.

© 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

20 . 6 . su mma ry 383

function prune!(ppt, grammar; p_threshold=0.99)
kmax, pmax = :None, 0.0
for (k, p) in ppt.ps
pmax′ = maximum(p)if pmax′ > pmax
kmax, pmax = k, pmax′
endend
if pmax > p_thresholdi = argmax(ppt.ps[kmax])
if isterminal(grammar, i)empty!(ppt.children[kmax])else
max_arity_for_rule = maximum(nchildren(grammar, r) forr in grammar[kmax])
while length(ppt.children) > max_arity_for_rule
pop!(ppt.children)
end
end
end
for child in ppt.children
prune!(child, grammar, p_threshold=threshold)end
return ppt
end

Algorithm 20.16. A method forpruning a probabilistic prototype
tree with root ppt, grammar
grammar, and pruning probability
threshold p_treshold.

Consider a node with a probability vector over the rule set:

Example 20.7. An example ofwhen pruning for probabilistic prototype trees should be applied.

R 7→ R + R
R 7→ ln(R )R 7→ 2 | x

R 7→ S

If the probability of selecting 2 or x grows large, then any children in theprobabilistic prototype tree are unlikely to be needed and can be pruned.
Similarly, if the probability of choosing S grows large, any children withreturn type R are unneeded and can be pruned.

© 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

384 c ha p te r 20 . e x p re ssi on op ti mi z a ti on

20.7

Exercises

Exercise 20.1. How many expression trees can be generated using the following
grammar and the starting set {R, I, F }?
R 7→ I | F
I 7→ 1 | 2

F 7→ πExercise 20.2. The number of expression trees up to height h that can be generated
under a grammar grows super-exponentially. As a reference, calculate the numberof expressions of height h can be generated using the grammar:16
N 7→ {N, N } | {N, } | {, N } | {}(20.6)

Let an empty expression haveheight 0, the expression {} haveheight 1, and so on.

16

Exercise 20.3. Define a grammar which can generate any nonnegative integer.
Exercise 20.4. How do expression optimization methods handle divide-by-zerovalues or other exceptions encountered when generating random subtrees?
Exercise 20.5. Consider an arithmetic grammar such as:
R 7→ x | y | z | R + R | R − R | R × R | R/R | ln R | sin R | cos R
Suppose the variables x, y, and z each have units, and the output is expected to
be in particular units. How might such a grammar be modified to respect units?
Exercise 20.6. Consider the grammar
S 7→ NP VP

NP 7→ ADJ NP | ADJ N
VP 7→ V ADVADJ 7→ a | the | big | little | blue | red
N 7→ mouse | cat | dog | pony
V 7→ ran | sat | slept | ate

ADV 7→ quietly | quickly | soundly | happily

What is the phenotype corresponding to the genotype [2, 10, 19, 0, 6] and thestarting symbol S?
© 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

20 . 7. e x e rc i se s 385
rpExercise 20.7. Use genetic programming to evolve the gear ratios for a clock. Assume all gears are restricted to have radii selected from R = {10, 25, 30, 50, 60, 100}.
Each gear can either be attached to its parent’s axle, thereby sharing the samerotation period, or be interlocked on its parent’s rim, thereby having a rotation
period depending on the parent’s rotation period and on the gear ratio as shownin figure 20.11.
The clock can also contain hands, which are mounted on the axle of a parent
gear. Assume the root gear turns with a period of troot = 0.1 s and has a radius of
25. The objective is to produce a clock with a second, minute, and hour hand.
Score each individual according to:minimize (1 − thand )
hands

2





+ minimize (60 − thand )
hands

2





rc

tc = t p rrpc
Figure 20.11. The rotation period
tc of a child gear attached to therim of a parent gear depends on
the rotation period of the parent
gear, t p and the ratio of the gears’
radii.

+ minimize (3600 − thand )
hands

2



+ #nodes · 10−3

where thand is the rotation period of a particular hand in seconds and #nodes isthe number of nodes in the expression tree. Ignore rotation direction.
Exercise 20.8. The four 4s puzzle17 is a mathematical challenge in which we use
four 4 digits and mathematical operations to generate expressions for each of theintegers from 0 to 100. For example, the first two integers can be produced by
4 + 4 − 4 − 4 and 44/44, respectively. Complete the four 4s puzzle.

W. W. R. Ball, Mathematical Recreations and Essays. Macmillan, 1892.

17

Exercise 20.9. Consider the probabilistic grammar
R 7→ R + R | R × R | F | I

wR = [1, 1, 5, 5]

I 7→ 1 | 2 | 3

wI = [1, 1, 1]

F 7→ 1.5 | ∞wF = [4, 3]

What is the generation probability of the expression 1.5 + 2?
Exercise 20.10. What is the probabilistic grammar from the previous question
after clearing the counts and applying a learning update on 1.5 + 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