Bayesian Methods for Efficient Genetic Programming.
A Bayesian framework for genetic programming (GP) is presented. This is motivated by the observation that genetic programming iteratively searches populations of fitter programs and thus the information gained in the previous generation can be used in the next generation. The Bayesian GP makes use of Bayes theorem to estimate the posterior distribution of programs from their prior distribution and likelihood for the fitness data observed. Offspring programs are then generated by sampling from the posterior distribution by genetic variation operators. We present two GP algorithms derived from the Bayesian GP framework. One is the genetic programming with the adaptive Occam's razor (AOR) designed to evolve parsimonious programs. The other is the genetic programming with incremental data inheritance (IDI) designed to accelerate evolution by active selection of fitness cases. A multiagent learning task is used to demonstrate the effectiveness of the presented methods. In a series of experiments, AOR reduced solution complexity by 20% and IDI doubled evolution speed, both without loss of solution accuracy.
Table

1. A summary of previous efforts for upscaling genetic programming
Subjects Selected references
Angeline w3x Banzhaf w7x
Bloat Landgon w22x Nordin w27x
Wu w39x Zhang w44x
Blickle w10x Iba w18x
Parsimony Kinnear w19x Rosca w31x
Soule w36x Zhang w44, 45x
Andre w1x Angeline w2x
Modularity Koza w21x O'Reilly w28x
Spector w37x Rosca w32x
Angeline w4x Banzhaf w7x
Operators Chellapilla w11x Koza w20x
Luke w24x Poli w29x
Gathercole w13, 14x Hillis w17x
Subset Schoenauer w33x Siegel w34x
Teller w38x Zhang w40, 42x
This list contains references which are directly related to the present work,
and thus it is far from complete.
produced by the solution. Soule et al. w35xobserve that removing non-functional
codes at every generation does not halt the program's growth. Instead, the
programs generate code which, while functional, is never actually executed.
Similar observations have been made in "biological" evolution. That is, introns
or non-coding segments emerge in DNA as well. Wu and Lindsay w39xgive a recent
review of biological introns. Banzhaf et al. w7xcharacterize introns as having two
salient features: An intron is a segment of the genotype that emerges from the
process of the evolution of variable length structures, and an intron does not affect
the survivability of the individual directly. Introns in GP turned out to be a mixed
blessing. On one hand, they may benefit evolution since they enable a genetic
program to protect itself against the destructive effect of crossover, and allow the
population to preserve highly-fit building blocks w27x. On the other hand, from the
practical point of view, introns usually result in run stagnation, poor results, and a
heavy drain on memory and CPU time w8, 35x.
An explanation for the cause of bloat in GP was provided by Zhang and
Muhlenbein w44x. Based on statistical learning theory, they show that the total error
of genetic programs can be decomposed into two terms attributed to bias and
variance. Since complex models are more expressive and thus better at reducing
the bias error than simple models, GP programs tend to grow until they fit the
fitness data perfectly unless complex models are penalized to reduce the variance
error. Empirical evidence supporting the explanation was given by an analysis of
the distributions of fitness vs. complexity of genetic programs generated by a large
number of GP runs for solving 7-parity problems.
Recently, Langdon and Poli w22xprovide a similar but more general explanation
for bloat. They argue that any stochastic search techniques, including GP, will tend
to find the most common programs in the search space of the current best fitness.
Since in general there are more of these which are long than there are which are
short but GP starts with the shorter ones.the population tends to be filled with
longer and longer programs. Based on this argument, Langdon w23xpresents new
crossover operators that carefully control variation in size and produce much less
bloat.
For the purposes of upscaling genetic programming, there are several reasons for
generally preferring parsimonious programs to complex ones. Parsimonious programs
typically require less time and less space to run. This is particularly
important during the GP process which may need to store and evaluate populations
of hundreds or thousands of programs. When hardware implementation of
GP solutions is in mind, or when GP evolves hardware circuits online "evolvable
hardware"., simple circuits require less hardware resource and execution time than
complex circuits.
In addition, statistical theory says that simpler models are likely to generalize
better on unseen fitness cases. As shown in w44x, when biases are equal, a simple
program has less variance in average than a complex program, resulting in smaller
total error. This is the statistical background behind the principle of Occam's
Razor w43xand the necessity for parsimony pressure in genetic programming.
However, an effective control of program growth is a difficult task, since too much
pressure on parsimony that reduces variance error.may lead to loss of diversity
and thus result in low accuracy greater bias error. The adaptive Occam method
w44xwas presented as a method for striking a balance between accuracy and
parsimony of GP solutions. It adapts the complexity penalty in the fitness function
during a run so that programs of minimal complexity are evolved without loss of
their accuracy. Several researchers support that parsimony pressure is one of the
easiest and effective methods for avoiding bloat w10, 18, 19, 31, 35, 36x.
Controlling program growth is one method for reducing time complexity as well
as space complexity of genetic programming. Another approach to upscaling GP is
by increasing modularity and reusability of programs. Koza w21xintroduces the
automatically defined functions or ADFs. An ADF is a subroutine that is evolved
during a run of genetic programming and which may be called by the main
program that is being simultaneously evolved during the same run. He reports that
GP with ADFs produced more parsimonious solutions and required fewer fitness
evaluation. Angeline and Pollack w2xdevelop an alternative method called module
acquisition MA. Rosca and Ballard w32xpresent mechanisms for adaptive representation
learning ARL.that creates new subroutines through discovery and
generalization of blocks of code. Spector w37xpresents a method for evolving a
collection of automatically defined macros ADMs. He shows that ADMs sometimes
provide a greater benefit than ADFs.
The search process of GP can be made more efficient by designing novel genetic
operators. Traditionally, crossover has been considered as the primary operator
and mutation as the secondary. This view is consistent with the notion that GP is
best made by combining building blocks, as hypothesized in bitstring genetic
algorithms. Koza w20xhas argued that mutations have little utility in GP because of
the position-independence of GP subtrees. However, the roles of building blocks
and of crossover have become increasingly controversial in recent years. Luke and
BAYESIAN GENETIC PROGRAMMING 221
Spector w24xexperimentally demonstrate that mutation can in fact have utility and
that crossover does not consistently have a considerable advantage over mutation.
Banzhaf et al. w6xalso found that increasing the mutation rate can significantly
improve the generalization performance of genetic programming. Angeline w4x
shows that mutation operations can perform on par with subtree crossover and
suggests that the building block hypothesis may not accurately describe the
operational chracteristics of subtree crossover. Chellapilla w11xproposes a method
for evolving computer programs without crossover. Poli and Langdon w29xpropose
various crossover operators and compare their search properties. Haynes w16x
how small changes in representation, decoding, and evaluation of
genetic programs can increase the probability of destructive crossover and mutation
while not changing the search space.
More recently, several authors have shown that fitness evaluation in genetic
programming can be significantly accelerated by selecting a subset of fitness cases.
The basic idea is that, other things being equal, the evolution time can be
minimized by reducing the effective data size for each generation. Siegel w34x, for
example, describes a GP method for evolving decision trees using subsets of given
training cases. Each fitness case has a fitness measure and the training cases which
tend to be incorrectly classified by decision trees become more fit, and therefore
selected more frequently. This method is motivated by competitive selection of
fitness cases in the host-parasite model of Hillis w17x. Similar ideas of competitive
selection of fitness cases have been refined by Teller and Andre w38xand by
Gathercole and Ross w14x. Here, the fitness of a training case is evaluated only
when the cost of evaluating another fitness case is outweighed by the expected
utility that the new information will provide.
Another class of methods for fitness case selection is based on incremental
learning. Schoenauer et al. w33xpropose the successive optimization scheme that
gradually refines the fitness. Evolution typically starts by considering a unique
fitness case, and additional fitness cases are "gradually" taken into account when
the current population meets some performance criterion. Though in a slightly
different context, Zhang w40xpresents a selective incremental learning SEL.
method in which an "incrementally" chosen subset of given training data is used
for fitness evaluation of the model. The data sets are divided into two disjoint sets
of candidate and training sets. Each candidate data point has a fitness value, called
criticality, defined as the error made by the model. Fitter candidates are then
incrementally selected into the training set as learning proceeds. SEL evolves a
single data set for a single model. The method presented in Section 5 generalizes
this to populations of multiple models and multiple data sets.
2.2. Eoling multiagent cooperation strategies using GP
A multiagent learning task is used as a testbed for the Bayesian genetic programming
methods. In an single table and four robotic agents are
placed at random positions, as shown in Figure 1. A specific location is designated
as the destination. The goal of the robots is to transport the table to the
222 ZHANG

Figure

1. The environment for multiagent learning. On a grid world of there are four robots,
obstacles, and a table. The task of the robots is to transport the table in group motion to the
designated target position G.
destination in group motion. The robots need to move in herd since the table is too
heavy and large to be transported by single robots. The robots share a common
control program A .
best
The objective of a GP run is to find a multi-robot algorithm A that, when
best
executed by the robots in parallel, causes efficient table transport behavior in
group. To evolve the programs, the robots activate each candidate program A for
in parallel to run a team trial. At the beginning of the trial, the robot
locations are chosen at random in the arena. They have different positions and
orientations. During a trial, each robot is granted a total of S elementary
movements. The robot is allowed to stop in less than S steps if it reaches the
goal. At the end of the trial, each robot i gets a fitness value which was measured
by summing the contributions from various factors.
The terminal and function symbols used for building GP trees to solve this
problem are listed in Tables 2 and 3. The terminal set consists of six primitive
actions: FORWARD, AVOID, RANDOM- MOVE, TURN- TABLE, TURN- GOAL, and STOP.
The function set consists of six primitives: IF- OBSTACLE, IF- ROBOT, IF- TABLE,
IF- GOAL, PROG2, and PROG3. Each fitness case represents a world of 32 by
grid on which there are four robots, 64 obstacles, and the table to be transported
see

Figure

1 for an example of fitness cases. A set of 200 training cases that are
generated at random is used for evolving the programs. An independent set of 200
random cases is used to test the performance of the evolved programs.
All the robots use the same control program. To evaluate the fitness of robots,
we made a complete run of the program for one robot before the fitness of another
is measured. The raw fitness value, eg., of individual i at generation g against
case c is computed by considering various factors. These include the distance
between the target and the robot, the number of steps moved by the robot, the
number of collisions made by the robot, the distance between starting and final
position of the robot, and the penalty for moving away from other robots.
BAYESIAN GENETIC PROGRAMMING 223

Table

2. GP terminal symbols for the multiagent task
Symbol Description
FORWARD Move one step forward in the current
direction
AVOID Check clockwise and make one step in
the first direction that avoids collision
RANDOM- MOVE Move one step in the random
direction
TURN-


Make a clockwise turn to the nearest
direction of the table
TURN- GOAL Make a clockwise turn to the nearest
direction of the goal
STOP Stay at the same position

Table

3. GP function symbols for the multiagent task
Symbol Description
IF- OBSTACLE Check collision with obstacles
IF- ROBOT Check collision with other robots
IF-


Check if the table is nearby
IF- GOAL Check if the table is nearby
PROG2, PROG3 Evalute two or three.subtrees in
sequence
3. Bayesian genetic programming
In this section we present a theory of genetic programing that is based on Bayesian
inference. The general theory is then applied to addressing two important issues in
genetic programming, i.e., control of program growth and acceleration of fitness
evaluation. This section aims to provide an outline of the Bayesian GP approach

Figure

2. An example genetic program for the multiagent learning task. Non-terminal nodes denote
checking sensor inputs of the robots, and terminals indicate actions to be taken. The meaning of the
symbols is described in Tables 2 and 3.
and some distinguishing features in different implementations. The algorithms are
detailed in the following two sections.
3.1. Bayesian formulation of genetic programming
Genetic programming works by initializing a population of programs and iteratively
producing the next generation of fitter programs, i.e. Ag.s Ag4 M, where Ag
denotes the ith program at generation g, and M is the population size. Genetic
operators such as mutation and crossover are used to produce offspring programs
from the parent programs. New generations are produced repeatedly until the
maximum number of generations g is reached or some other termination
condition is satisfied. The goodness or fitness of a program is measured in terms of
a set D of fitness cases or training data and the programs can be considered as a
model of the unknown process f generating the data.
In the Bayesian GP, the best program is defined as the most probable model,
given the data D plus the prior knowledge on the problem domain. Bayes theorem
provides a direct method for calculating such probabilities w15x. It states that the
posterior i.e. after observing the data D.probability of a program A is
PDN A.PA.
PAN D.s
PD.
PDN A.PA.
s ,1.
A
where A is the space of all possible programs in case of A taking discrete values,
the integral will be replaced by summation. Here PA.is the prior i.e. before
observing the data.probability distribution for the programs, and PDN A.is the
likelihood of the program for the data.
The relationship between the states of models and their probabilities is established
by borrowing the concept of "energy" from statistical physics w26x. We regard
the GP system as a thermodynamic system. Every possible state s of the system has
some definite energy Es. The system energy can fluctuate and it is assumed that
the probability of a thermodynamic.system being in state s, given that the
temperature is T is given asPs.s expyEs.rT .,2.
Z
where Es.is the energy of the system and Z is the normalization constant needed
to make the distribution integrate or sum.to one. This distribution is known as
the canonical or Boltzmann.distribution.
BAYESIAN GENETIC PROGRAMMING 225
Under the canonical distribution, the prior distribution of programs A can be
expressed asPA.s expyaFA .,3.
ZAa .
where Z a.is a normalizing constant, a s 1rT, and F is the "energy" of model
AA
A in the "equilibrium" state. For example, the energy function, F , can be chosen
A
as the total number of nodes in genetic program A. This choice of prior distribution
says that we expect the program size to be small rather than large, thus
implementing a parsimony pressure. Similarly, the likelihood of models A for the
fitness data D can be expressed asPDN A.s expybFD .,4.
ZD b .
where F is an error function for the fitness set D, b controls the variance of the
noise, and Z b .is a normalization factor. The likelihood factor gives preference
to programs that fit better to have less error for.the fitness cases.
Initially, the shape of the prior probability distribution of programs PA.is flat
to reflect the fact that little is known in advance. Evolution is considered as an
iterative process of revising the posterior distribution of models PAN D.by
combining the prior PA.with the likelihood PDN A. In each generation, Bayes
theorem 1.is used to estimate the posterior fitness of individuals from their prior
fitness values. The posterior distribution PAN D.is then used to generate its
offspring.
The objective of Bayesian genetic programming Figure 3.is to find a program
Ag that maximizes the posterior probability:
best
est s min arg max PgAig N D. 5.
AiggAg.

Figure

3. Outline of the Bayesian genetic programming procedure.
226 ZHANG
where g is the maximum number of generations and Ag.is the populations of
size M:
Ag.s Aig 4Mis1.6.
The posterior probability P Ag N D. of program Ag is computed with respect to
the gth population:
PDN Aig .PAig .
PgAig N D.s Mg,7.
where P DN Ag .is the likelihood and PAg .is the prior probability of or degree
ii
of belief in. Ag. Note that the posterior probability is approximated by a fixed-size
population Ag.which is typically a small subset of the entire program space A.
Genetic operators are applied to generate L offspring AX , k s 1,.,L. For-
mally, this proceeds in two steps. First, candidates are generated by sampling from
the proposal distribution:
QgAXkN Aig .8.
The specific form of Q ?N ?.is determined by variation operators. For example,
for subtree crossover from two parents it is given as:
QCAXkN Ai.s PSAjN Ai.PRAXkN Ai, Aj.9.
AjgAg.
Here PAN A .is the probability of individual A in Ag.being selected as a
mate for A and P AX N A , A . is the probability of A and A producing AX by
crossover.
Then, each candidate generated by genetic operators is accepted with probability
PgAXkN D.
agAXkN Aig .s min 1, PgAig N D.5, 10.
where P Ag N D.is computed by 7.and PAX N D. is the posterior probability of
A estimated with respect to the current population:
PDN AXk.PAXk.
PgAXkN D.s Mg.1.
If AX is rejected in 10., then Ag is retained, i.e., AX / Ag. Note that this
kiki
acceptance function does not exclude the case that A is generated by crossover
from Agand another parent Ag Ag., j / i.
BAYESIAN GENETIC PROGRAMMING 227
After L offspring AX , k s 1,.,L, are generated, M of them are selected to
build the new population:
Ag q 1.s Aigq14is1.12.
This defines the posterior distribution P Ag N D. at the next generation. It
should be mentioned that this formulation of offspring selection is intentionally
very general so that it can accommodate various forms of existing selection
schemes, such as m, l.selection w5, 25x.
In effect, the evolutionary inference step from generation g to g q 1 is considered
to induce a new fitness distribution P Ag N D.from priors PAg . through
posterior distribution P Ag N D. following Bayes formula, using genetic operators.
Based on this theoretical framework we present in the following subsections two
examples of Bayesian genetic programming that employ specific techniques for
effective control of evolutionary dynamics. Detailed procedures and experimental
results are described in Sections 4 and 5.
3.2. GP with the adaptie Occam's razor
The first Bayesian GP focuses on the fact that genetic programming iteratively
searches populations of more probable more likely in terms of the data and the
priors. programs and thus the information gained in the previous generation can
be used in the next generation to revise the prior before seeing the data.belief in
true programs. Thus, the posterior distribution can be written in the form
PDN Aig .Pgy1Aig .
PgAig N D.s ,13.
Mjs1 PDN Ajg.Pgy1Ajg.
where the priors P Ag . are now expressed explicitly as a function of generation
rather than the fixed prior P Ag .as in 7. After computing the posterior
probabilities for Ag, the priors PA.are revised to PA.by a belief update
function u?, ? .:
PgA.s uPgy1A., PgAig N D.14.
An implementation of u?, ? .will be described in the next section.
Basically, the Bayesian GP starts with generating programs according to an
initial prior distribution PA.on the program sizes. Typically, the prior distribution
is given as uniform, reflecting the fact that little is known in advance about the
optimal size of genetic programs. Then, the fitness Fg.of the programs is
measured on the training cases, which results in the estimation of the likelihood
programs more details in Section 4. Combining the prior and the
likelihood by Bayes formula, we get the posterior probabilities P Ag N D. The
current prior is then updated to reflect the new information gained from the
posterior distribution.
Note that we assign prior distributions on the complexity of models programs.
In addition, information theory w12xsays that the probability and code length
complexity.of models are related as LA.sylog PA. By making use of this,
the complexity of programs can be controlled to evolve parsimonious and accurate
programs. In fact, we show in Section 4 that the adaptive Occam method presented
in w44xis derived from the Bayesian genetic programming framework.
3.3. GP with incremental data inheritance
In the second example of the Bayesian approach to GP, we make use of the fact
that the Bayes formula suggests an incremental, evolutionary learning rule: infer
programs of higher posterior probability from the existing programs by observing
new fitness cases. This leads to writing the posterior distribution as
PDig N Aig .Pgy1Aig N Digy1 .
PgAigN Dig.s Mg ,15.
where the data set Dg is now a variable of generation number g and specific to a
single program Ag. The collection of Dgconstitutes the data population Dg.
ii
The observation of new data will lead to update of the prior distribution:
PgAN Dig .s uPgy1AN Digy1 ., PgAig N Dig .16.
The revision of prior distribution is the same as the process for equation 14.,
except that the data Dgy1 for the estimation of likelihood and thus the posterior
probabilities of programs. is now a function of generation rather than fixed as D.
The genetic programming with the incremental data inheritance IDI.method
w42xis an example of this approach, where Dg specifically satisfies the following
conditions:
where D is the entire data set given for training. In this method, the fitness of
programs is estimated on incrementally chosen data subsets Dg, rather than on the
whole data set D, and thus the evolution is accelerated by reducing the effective
number of fitness evaluations. More details are described in Section 5.
4. Bayesian GP for parsimonious solutions
Statistical theories suggest that models which are too simple lack sufficient learning
capability while models which are too complex may generalize poorly on unseen
BAYESIAN GENETIC PROGRAMMING 229
data. As reviewed in Section 2.1, several researchers have observed that the
program size tends to grow without bound or "bloat" see, for example, w22x. In
this section we describe a Bayesian method for evolving parsimonious programs.
4.1. Algorithm description
The GP algorithm for the adaptive Occam method Figure 4.is the same as that of
the Bayesian GP procedure described in Figure 3, except three differences. The
first is the raw fitness calculated in step 2. The derivation of the fitness function is
given in Section 4.2. The second difference is in step 3, where the posterior
probability 13.is substituted for 7.; in 13.the prior is revised each generation
while in 7.it is constant. The third is the additional step 6 for the revision of
priors.
4.2. Fitness ealuation
For a convenient implementation of the Bayesian evolutionary algorithm we take
the negative logarithm of the posterior probability P Ag N D. and use it as the
fitness function
Fig.sylog PgAig N D.,18.

Figure

4. Outline of the Bayesian genetic programming with the adaptive Occam's razor AOR.
where Ag g Ag. Then the evolutionary process is reformulated as a minimization
process
est s min arg min Fig.,19.
AiggAg.
where the fitness function is expressed as
Fig.sylog PDN Aig .y log PAig .20.
As described in the previous section, we can write the likelihood function in
Bayes' theorem 1.in the form w9xPDN A.s expybFD . 21.
ZD b .
where F is an error function, b controls the variance of the noise, and Z b .is a
normalization factor. If we assume that the data has additive zero-mean Gaussian
noise, then the probability of observing a data value y for a given input vector x
would be
where Z b .is a normalizing constant. Provided the data points are drawn
independently from this distribution, we have
PDN Aig.s PycN xc, Aig. 23.
cs1s expybFD .24.
ZD b .
where x , y .g D are training cases and F is given as
cc D
FDs EDN Aig.s f xc; Aig.y yc.2.25.cs1
If we also assume that a Gaussian prior on the architecture of program Ag,we
ZAa .
BAYESIAN GENETIC PROGRAMMING 231
where Z a.is a normalizing constant. For example, F can be chosen in the form
AA
FAs CAig .s uk2 27.ks1
where u are the parameters defining the program Ag. This choice of prior
distribution says that we expect the complexity parameters to be small rather than
large, thus implementing a parsimony pressure.
Substituting 24.and 26.into 20., the fitness function can be expressed as
Fig.s bFDq aFA 28.
s bEDN Aig .q aCAig .,29.
where the first term reflects the error and the second the model complexity.
The exact calculation of the constants b and a in equation 29.requires the
true probability distribution of underlying data structure, which is in most real
situations unknown. Instead, we define an adaptive fitness function in its most
general form as
Fig.s Eig.q a g.Cig.,30.
where Eg.and Cg.are the measures for error and complexity of the program,
ii
and the parameter b is absorbed into the adaptive parameter ag.which balances
the error and complexity factors as follows:
if Ebestgy 1.) e
a g.s~ N 2 C^bestg. 31.N 2 Ebestgy 1.? C^bestg. otherwise.
This is the adaptive Occam method w44x. User-defined constant e specifies the
maximum training error allowed for the run. Egy 1.is the error of the best
best
program of generation g y 1. C^ g.is the size of the best program at generation
best
g estimated at generation g y 1. These are used to balance the error and
complexity terms to obtain programs as parsimonious as possible while not sacrificing
their accuracy. The procedure for estimating C^ g.is given in w44x.
best
4.3. Discussion
The necessity and difficulty of non-coding segments or introns in genetic programming
have been studied by many authors w22, 27, 39x. The adaptive Occam method
deals with the non-coding segment problem using a two-phase strategy by means of
an adaptive fitness function. In the first stage during Egy 1.) e., the
best
growth of non-coding segments is encouraged to increase the diversity of partial
solutions. In the second stage during Egy 1.F e., a strong parsimony
best
pressure is enforced to prefer compact solutions. The transfer from the first stage
to the second is controlled by Bayesian inference under the constraint of the
user-defined parameter e.
Note that the posterior probability P Ag N D. can be computed from the fitness
values Fg.by taking its exponential function:
PDN Aig .Pgy1Aig .
Mjs1 PDN Ajg.Pgy1Ajg.
expyEig.y a g.Cig.
s 32.
Mjs1 expyEjg.y a g.Cjg.
This shows that the revision of priors P Ag . is implemented as the update of
ag.in the adaptive Occam method since the priors are reflected in C^ g.which
best
leads to revision of ag.
We also note that minimization of Fg.is equivalent to the minimum description
length MDL.principle w18, 30x: the best model is a model whose total code
length for model description L AgN Dg.and data description LDgN Ag. are
ii ii
minimal. In information theory w12xthe optimal description length for a model is
given as the negative logarithm of probability of the model:
LAig .sylog PAig .3.
Similarly, the code length for the data given the program Ag is given as:
LDN Aig .sylog PDN Aig .4.34.
.correspond to the two terms in 20.
ii
4.4. Experimental results
The adaptive Occam method for complexity control was applied to the multiagent
learning task. Experiments have been performed using the parameter values listed
in

Table

4. The terminal set and function set consist of six primitives, respectively,
as given in Tables 2 and 3. A set of 200 training cases was used for evolving the
programs. An independent set of 200 cases was used for evaluating the generalization
performance of evolved programs.
The Eg.-values of program i at generation g are measured as the average of
its raw fitness values eg.less is better.for the cases c:
Eig.s ei, cg.,35.
BAYESIAN GENETIC PROGRAMMING 233

Table

4. Parameters used in the experiments for GP with the adaptive Occam's razor
Parameter Value
Population size 100
Max generation
Crossover rate 0.9
Mutation rate 0.1
Training cases 200
Test cases 200
where N is the number of fitness cases in the training set. Each eg.-value was
computed by considering such factors as the distance between the target and the
robot, the number of collisions made by the robot, and the distance between
starting and final position of the robot. Note that in the multiagent task there is no
target output given for the training case. The goodness of a program is measured
by scores or penalties.it collects by running the robots.
The complexity of a GP tree is defined as the size and depth of the tree:
Cig.s
where k is a constant. We used k s 2.

Figure

5 compares the fitness values for GP with and without AOR, averaged
runs. Shown are the Eg.-values of the best individuals in each generation
for each method. A tendency can be observed that the GP with the adaptive
Occam's razor converges slightly faster than the baseline GP, though there is no
significant difference in their final fitness. Figure 6 compares the typical changes of
program complexity for both methods. The tree complexity was measured in terms
of the tree size plus

Figure

5. Comparison of fitness values Eig.-component only.of the genetic programs for the
multiagent task. A tendency can be observed that the GP with the adaptive Occam's razor converges
slightly faster than the baseline GP, though there is no significant difference in their final fitness.

Figure

6. Evolution of the complexity of GP trees for the multiagent task. The GP with the adaptive
Occam's razor AOR.promotes the trees to grow when significant fitness error.reduction is required,
while it prefers smaller trees to larger ones when their fitness error.is comparable. In contrast, the GP
without the Occam factor tends to grow as generation goes on.
More detailed results are summarized in Table 5. The GP with AOR achieved,
on average, a 20% reduction of program complexity without loss of solution
accuracy in fact, with a slight improvement both in training and test performances.
It can be concluded that the adaptive Occam method evolves smaller programs
without loss of generalization capability of the programs.
5. Bayesian GP for accelerated evolution
5.1. Algorithm description
The algorithm for incremental data inheritance Figure 7.is the same as the GP
with AOR, except that the training set Dg now increases with generation. The
additional step for this modification is step 5, where the training set grows by
inheritance. The next section details the data inheritance procedure.

Table

5. Effects of the adaptive Occam's razor
GP programs for the multiagent learning task
Method Complexity
Baseline 31.4 " 6.5
AOR.on the complexity and average fitness values of
Average Fitness
Training Test
The sizes of training and test sets were 200, respectively. The values are averaged over ten runs.
BAYESIAN GENETIC PROGRAMMING 235

Figure

7. Outline of the Bayesian genetic programming with incremental data inheritance IDI.
5.2. Data inheritance procedure
The basic idea in genetic programming with incremental data inheritance is that
programs and their data are evolved at the same time. With each program is
associated a separate data set. We describe a variant of uniform crossover that we
call uniform data crossover. A simplified example for illustrating this process is
given in Figure 8.
First, two parent data sets, Dg and Dg, are crossed to inherit their subsets to two
offspring data sets, Dgq1 and Dgq1. Second, the data of parents' are mixed into a
union set
which is then redistributed to two offspring Dgq1 and Dgq1, where the size of
offspring data sets is equal to N s N q l, where l G 1 is the data increment
size. Thus, the size of data sets monotonically increases as generation goes on.
To maintain the diversity of the training data during inheritance we import some
portion of data from the base set. The import rate r is given as
ris r ? 1 y di.,0F r F 1. 38.

Figure

8. An illustrative example for data inheritance in IDI. Two parent data sets, Dig and Dig, are
merged to form the union set, Digqj, which is then inherited to two offspring data sets, Digq1 and Djgq1.
The off-spring data points with shaded circles are the examples imported from the base data set D to
maintain the diversity of data sets.
where r is a constant for import strength. The diversity d is measured as the ratio
of distinctive examples in the union set:
d sy1, 0F d F 1. 39.

Figure

8 illustrates the process of data inheritance, where two parent data sets of
size 6 each are unioned to form the genetic pool of size 10, from which two
offspring data sets of size 8 each are inherited. Marked are the data imported from
the basis data set to maintain the diversity.
5.3. Fitness ealuation
As in the implementation of the Bayesian GP with the adaptive Occam's razor, we
take the negative logarithm of P Ag N D. and use it as the fitness function
Fig.sylog PgAig N Dig .,40.
where Ag g Ag., Dg ; D, and P AgN Dg.are defined by 15.:
ii gi
gg PDig N Aig .Pgy1Aig N Digy1 .
PgAiN Di.s Mg gy1 .41.
Then the evolutionary process is reformulated as a minimization process
est s min arg min Fig.,42.
gFgmam Aig, Dig
where the fitness function is expressed as
Fig.sylog PDig N Aig .y log Pgy1Aig N Digy1 .43.
BAYESIAN GENETIC PROGRAMMING 237
Using the similar arguments as in Section 4.2 we can write the likelihood
function and the prior distribution of programs in the formPDigN Aig.s expybFD. 44.
ZAa .
where Z b .and Z a.are normalizing constants.
D A
We define an adaptive fitness function in its most general form as
Fig.s Eig.q a g.Cig.,46.
where Eg.and Cg.represents the error factor F and complexity factor F
ii D A
with a and b absorbed into ag. The parameter ag.balances these two factors
as follows
if Ebestgy 1.) e
a g.s~ Ng2 C^bestg. 47.Ng2 Ebestgy 1.? C^bestg. otherwise.
This is a generalization of the adaptive Occam method in that N is now a variable,
rather than fixed value N, and scheduled to increase monotonically as a function of
generation g.
5.4. Discussion
The incremental data inheritance method has interesting properties that characterize
the stability of the algorithm as generation goes on. It maintains a growing
subset of given fitness cases for each program. The monotonic growth of the fitness
subset ensures that the IDI method eventually achieves the same level of performance
as that of the baseline algorithm. This is contrasted with many of existing
methods for fitness case selection, including LEF w14xand RAT w38x. Maintaining a
separate subset for each program allows the learned component of the program
can be retained during evolution. At the same time, the import of new fitness cases
from the base set allows the programs to learn new situations, thus leading to
gradual improvement in performance.
5.5. Experimental results
We compare the performance of the GP with incremental data inheritance IDI.to
the baseline GP, i.e. one that uses the complete training cases from the outset.
Experiments have been performed using the parameter values listed in Table 6. GP

Table

6. Parameters used in the experiments for GP with
incremental data inheritance
Parameter Value
Population size 100
Max generation
Crossover rate 0.9
Mutation rate 0.1
Training cases 200
Test cases 200
Initial data size N . 20Data increment size l. 6
runs with IDI used 20 q 6g examples for each generation g selected from the
given data set, i.e. N s 20, l s 6, for fitness evaluation. For all methods, a totalof 200 training cases were used for training and an independent set of 200 test
cases was used for evaluating the generalization performance of evolved programs.
Fitness factors Eg.and Cg.are the same as in the experiments for GP with the
ii
adaptive Occam's razor, except that the training set size for Eg.is now N
ig
instead of N.
Results are shown in two different forms. Figure 9 shows the evolution of fitness
as a function of the generation number, in which there is no significant difference
in the performance. Since the GP with IDI uses variable data size, computing time
at each generation should be measured by a product of the population size and the
data size. This result is shown in Figure 10, where IDI achieved a speed-up factor
of approximately four compared with the baseline GP. More detailed results are
summarized in Table 7. The incremental data inheritance IDI.method used only
50% of the time required for the baseline GP. It is interesting to see that, despite

Figure

9. Comparison of fitness values as a function of the generation number. The curves are mean
values over ten runs. The GP with IDI shows no loss of fitness values compared to that of the baseline
GP algorithm.
BAYESIAN GENETIC PROGRAMMING 239

Figure

10. Comparison of fitness values as a function of the number of function evaluations. The
curves are mean values over ten runs. The GP with IDI converges much faster than the baseline GP
algorithm.

Table

7. Effects of incremental data inheritance IDI.on the time and average fitness
values lower is better.of GP programs for the multiagent learning task. The sizes of
training and test sets were 200, respectively. The values are averaged over ten runs. Time is
measured as the total number of fitness evaluations. Also shown are the standard deviations.
Average Fitness
Method Time Training Test
Baseline 1,220,000 215.8 " 8.3                                     221.3 " 7.6
IDI
the reduced data size, the generalization performance of the GP with IDI was
slightly better than that of the baseline GP.
6. Concluding remarks
We have presented a Bayesian framework for genetic programming. Two specific
GP algorithms were derived from this framework for reducing the time and space
complexity of genetic programming. Applied to the multiagent learning task, the
first method, i.e. GP with the adaptive Occam's razor AOR., achieved approximately
20% reduction of program complexity without any loss in fitness values. The
second Bayesian approach to genetic programming, i.e. the incremental data
inheritance IDI.method, used only 50% of the time required for the baseline GP
to achieve the same or a little better fitness level.
Though this improvement is significant, it should be mentioned that there is still
much room for further reduction of spacertime complexity of genetic program-
ming. Future work should address the following issues, among others. One is
incorporating better genetic operators. We have focused in the present work on
dynamics at the phenotypic level. Further improvement can be achieved by finding
more "intelligent" variation operators that adapt to the dynamics at the genotypic
level. In terms of Bayesian genetic programming, this involves adapting proposal
functions. The second issue is to improve modularity of genetic programs by
automatically designing and adapting reusable submodules such as ADFs or
libraries. This, combined with the Occam's razor, will further improve the comprehensibility
and reusability of genetic programs as well as speed up the GP process.
From the theoretical point of view, the Bayesian framework for genetic programming
provides a number of important features. One is that, by formulating the GP
process as Bayesian inference, principled techniques for driving evolutionary dynamics
of conventional GPs can be developed. In addition to the methods presented
in this paper, one can think of other methods within the framework. For
example, decision-making can be made more robust by combining multiple programs
instead of a single GP tree. The Bayesian GP framework provides a
principled way to combine multiple programs to build a committee machine.
Another important feature of Bayesian inference is that it allows background
knowledge in the problem domain to be incorporated in a formal way. For
instance, if we can guess the distribution of specific function symbols for good GP
trees, this knowledge can be reflected in the application probabilities of genetic
operators. Background knowledge is important especially for solving real-life
problems of practical interest.
Finally, the Bayesian analysis of genetic programming appears to be a useful tool
for incorporating various genetic programming procedures into a uniform frame-
work. A theoretical framework is crucial for the design and comparative analysis of
various genetic programming algorithms.

Acknowledgments

This research was supported by the Korea Science and Engineering Foundation
KOSEF.Granta 981-0920-350-2, by the Korea Research Foundation KRF.
Granta 1998-001-E01025, and by the Korea Ministry of Science and Technology
through KISTEP Granta BR-2-1-G-06. Thanks to Wolfgang Banzhaf and three
anonymous reviewers for helpful comments that improved the readability of the
paper.


--R

"Automatically defined features: The simultaneous evolution of 2-dimensional feature detectors and an algorithm for using them,"

"Genetic programming and emergent intelligence,"
An Introduction
"Evolving compact solutions in genetic programming: A case study,"
"Evolutionary programming with tree mutations: Evolving computer programs without crossover,"
Elements of Information Theory
"Dynamic training subset selection for supervised learning in genetic programming,"
"Small populations over many generations can beat large populations over few generations in genetic programming,"
Bayesian Data Analysis
"Perturbing the representation, decoding, and evaluation of chromosomes,"
"Co-evolving parasites improves simulated evolution as an optimization procedure,"
"Genetic programming using a minimum description length principle,"
"Generality and difficulty in genetic programming: Evolving a sort,"
On the Programming of Computers by Means of Natural Selection
"Scalable learning in genetic programming using automatic function definition,"
"Fitness causes bloat: Mutation,"
"Size fair and homologous tree crossovers,"
"A comparison of crossover and mutation in genetic programming,"
"The science of breeding and its application to the breeder genetic algorithm,"
"Probabilistic inference using Markov chain Monte Carlo methods,"
"Explicitly defined introns and destructive crossover in genetic programming,"

"On the search properties of different crossover operators in genetic programming,"
"Stochastic complexity and modeling,"
"Analysis of complexity drift in genetic programming,"
"Discovery of subroutines in genetic programming,"
"Evolutionary identification of macro-mechanical models,"
"Competitively evolving decision trees against fixed training cases for natural language processing,"
"Code growth in genetic programming,"
"Effects of code growth and parsimony pressure on populations in genetic programming,"
"Simultaneous evolution of programs and their control structures,"
"Automatically choosing the number of fitness cases: The rational allocation of trials,"
"Empirical studies of the genetic algorithm with noncoding segments,"
"Accelerated learning by active example selection,"
"A Bayesian framework for evolutionary computation,"
"Genetic programming with incremental data inheritance,"
"Genetic programming of minimal neural nets using Occam's razor,"
"Balancing accuracy and parsimony in genetic programming,"
"Evolutionary induction of sparse neural trees,"
--TR

--CTR
Byoung-Tak Zhang, A unified Bayesian framework for evolutionary learning and optimization, Advances in evolutionary computing: theory and applications, Springer-Verlag New York, Inc., New York, NY,
Sean Luke, Modification point depth and genome growth in genetic programming, Evolutionary Computation, v.11 n.1, p.67-106, Spring
bayesian genetic programming;incremental data inheritance;adaptive Occam's razor;probabilistic evolution;parsimony pressure;data subset selection