Evolution in a Computer, pt. 3

After months of waiting, I welcome you to the final part of my humble introductory mini-series about evolutionary computation. And today's topic is Genetic Programming

Discovering functions

Recall the example of finding the minimum of the function from the first part of this mini-series? We had a function we didn't know anything about (actually, we did, but that was just for the sake of showing what we can do when we know) and we wanted to find a minimum of that function. Now, let's look at a different problem: we don't have the function at all, only a handful of \(x, y\) pairs, and we want to find out what the function actually looks like. How to do it?

If we restrict ourselves to a specific set of functions, e.g. linear functions, we can use specialized methods like (multiple) linear regression that will just throw the function at us given the data. Period. However, this has the problem of restricting to only a specific set of functions. If the "secret" function was actually a linear function, it would be the best we could do but if it wasn't, we would get errorneous function.

We could take a different approach and use more complicated machine learning methods like Artificial neural networks (ANN) that have the universal approximation property (simplified: given enough neurons they are, in principle, capable of representing almost any function with arbitrarily low error). However, their problem is that they can be immensely complex and if one wanted to use it as a mathematical expression, it would be mess.

As you might already have guessed, we are going to use.. Evolutionary algoritm! Well, we know how to evolve binary strings and real numbers... how do we use it to find a function? Well, it is very complicated. You could design some algorithm that would translate binary string to a function but it would probably be very messy and not efficient, not speaking about the fixed length of the binary strings which would be a major problem. However, there is a representation that can be used... trees!

Tree representation

A mathematical expression can be described like a tree. Let's take some made-up expression like $$ (x + \sin{y})^2 - \cos{3x} $$ Where is the tree in there? It's easy - the operators (i.e. addition, mutiplication...) and functions (i.e. sin, cos...) are inner nodes and the variables (i.e. x, y...) and constants (i.e. numbers) are leaves of the tree. The corresponding tree for the expression above looks like this:

Evolution of trees

So now that we have a representation for math expressions, how do we actually evolve it? Well, pretty much in the same way as in classical GA - by selection, recombination and mutation. Selection is identical - the better solutions (in this case those that have lower error) are preferred more than the worse ones. Recombination, i.e. crossover, and muation are the more interesting ones.

Crossover

There are more types of crossover but I'm going to speak about the most basic type - the subtree crossover. This crossover works in the following way:

1. Pick a random node \( N_1 \) in the first parent
2. Pick a random node \( N_2 \) in the second parent
3. Exchange the subtrees rooting in \( N_1 \) and \( N_2 \) between the two parents

Look at this example: And that's it! It's quite simple. Sometimes some tweaks are used, e.g. giving more preference to nodes that are deeper in the tree, but in essence, the basics are the same - just swapping subtrees.

Mutation

As there is subtree crossover, there is subtree mutation too. Again, it is based around mangling subtrees. However, contrary to the crossover, in mutation we don't exchange subtrees between multiple solutions but we modify only one. The most basic subtree mutation looks like this:

1. Randomly pick a node
2. Throw away that node (and its subtree) and replace it with a new randomly generated subtree

No illustration here - it's really that simple.

Demo

This is a little demo of GP in action. The task is to find the function based on a bunch of samples (the red points; the sought function is the dashed curve). The actual function that we want the GP to find is $$ f(x) = 5 \sin{x^2} $$ ATTENTION: you might find that the algorithm "runs away" with bigger and bigger and messier and messier functions. I suggest you stop the algorithm and start again. This algorithm is really not tuned very well, it's just a demo.

Clicking Start starts the evolution. Clicking Cancel during the evolution terminates it. You can play with the population size (must be positive even number), crossover and mutation probs (must be between 0 and 1) and the tournament size (must be between 2 and size of population). The evolution runs indefinitely (i.e. until you stop it). You can also play with the number of available data points (the more the easier it should be for the algorithm to find the solution). This thing is written in Dart and you can find the source here.

Not just functions

We can evolve functions. But one could say “Fine, but where's the programming in that?” and she'll be totally right. So far we can make mathematical expressions. But the tree representation goes beyond that. All you need to do is to represent a program like a tree. And it is quite simple. Let's imagine that we control an ant in a 2D grid world. The ant can turn left, turn right and decide what to do based on whether there is a piece of food just ahead of it. And here is what a program in this domain might look like: As you can see, the only difference is that the branching point (or non-terminal) is not an operator but a decision-making point and the leaves are not variables and constants but actions. In fact, there could be many more non-terminals, e.g. cycles, function calls etc.

Closure property

All the Genetic Programming I introduced above is called Koza-style GP because it is the idea of John R. Koza which he published in his book [1]. For this type of GP to work a so called closure property must hold. Closure property holds if an output of every terminal and non-terminal is acceptable by every non-terminal as an input. That menas that whatever our nodes output, it must be possible to chain it with an arbitrary node.

If our non-terminal set is made of plus, minus and multiplication and our terminal set of variables and numbers, closure property holds - however we add, subtract or multiply numbers, we always get a number that can be furter added to, subtracted from or multiplied. However, if we extended our non-terminal set by a non-terminal such as if that would take three children - a condition, a result when condition is true and a result when condition is false - then the closure property does not hold - the condition requires a boolean value but that is different from numbers. We would have to employ some rule that would define true and false in terms of numbers.

In fact, it would be enough to just add division to the non-terminal set. Yes, it would still produce only numbers but what if the second argument to a division resulted to zero? We can't divide by zero so division cannot accept zero as its second argument. You see, it is sufficient that a non-terminal doesn't accept some of the values even though they are of the same type.

Summary

In this long-awaited (I hope :)) post I described a little bit of what is known as Genetic Programming. It is basically a way of evolving mathematical expressions, computer programs or basically any tree-like structure using the principles of genetic algorithms. It might (or might not) seem as a brilliant technique able to solve anything... but don't let yourself be fooled. Finding functions and programs is really hard and most of the time this vanilla GP doesn't work very well. But it's a start and it can solve at least some easy problems.

Regarding further posts I'm going to publish here, I don't have any clear plan. I suppose I might dive into specific algorithms or things I do or something totally different. We'll see :).

---

^[1] Koza, John R. Genetic programming: on the programming of computers by means of natural selection. Vol. 1. MIT press, 1992.

Evolution in a Computer, Pt. 2

After a long waiting, it's here! To seamlessly continue from the previous part, I prepared a little demo of a genetic algorithm. What we have here is a so-called Facility location problem, or a variation of it.

There is a bunch of cities (50 in our case) and we want to supply the cities with goods so we must open some factories. We must decide how many factories we are going to open and in which cities. We won't explicitly assign which cities are to be supplied by which factories but a city will be assigned to the nearest factory (imagine the factories have unlmitted capacity). The catch: we want to minimize the total cost:

• every opened factory costs us some fixed amount of money
• delivering goods from a factory to a city costs us an amount of money equal to the (euclidean) distance from the factory to the city

Clearly, we cannot try every combination of cities with a factory as it's \(2^{50} \approx 1.13 \cdot 10^{15}\) combinations. That would took almost 36 years if checking one solution took one microsecond. Yes, this problem is NP-hard. Instead, we are going to use a... Genetic Algorithm!

So how do we represent a solution to our problem? Well, we'll use a binary string, as in proper GA, with the number of bits equal to the number of cities. Then \(1\) in position \(n\) means that the \(n\)-th city has a factory and vice versa, e.g. a string like 001000110 would mean that the third, seventh and eighth city would have a factory while the first, second, fourth, fifth, sixth and ninth wouldn't. Now, give it a try!

Clicking Start starts the genetic algorithm from beginning with the current map. Clicking New map generates a new (random) map. Clicking Stop during the optimization stops it. You can edit the cost of opening a facility, size of population (must be positive even number), probability of crossover and mutation (must be between 0 and 1) and size of tournament (must be between 2 and size of population). The evolution runs indefinitely (i.e. until you stop it). The world is 100 by 50 units. This thing is written in Dart and you can find the source here.

Real-valued Evolution

So far we have spoken about a genetic algorithm that uses binary strings. This means that our solutions are inherently discrete. But what if we wanted to really find the minimum of the function from the part one?

The „standard“ way

The first way of dealing with real number is to just express them using a binary standard, e.g. the IEEE 754 single-precision floating-point format. But this has problems. Imagine numbers 1.763 and -2.611. These would look like this (e.g. using this tool):

 1.763 = 00111111111000011010100111111100
-2.611 = 11000000001001110001101010100000

And now imagine what would happen when we crossed them over using single-point crossover after the 5th bit:

00111111111000011010100111111100
11000000001001110001101010100000
turns into
00111000001001110001101010100000 = 3.9840699e-05
11000111111000011010100111111100 = -115539.96875

One could say "Well, that escalated quickly!" and indeed, these problems arise all the time. If we explore all the possible crossovers (there are 32 of them resulting in 64 numbers) they include gems like -3.252161001305617e+19, 3.591135600494495e-39 and NaN. Bit-flip mutation is going to have similar problems.

The good way

Much better way than sticking to the binary representation is just use the real numbers. Then, if we have a one-dimensional problem (like in searching for the minimum of that function), the genotype is just one real number. This brings a lot of changes in the crossover and mutation operators.

As we have seen, just crossing over the binary representation is bad. Therefore we need to find a better way of combining two individuals. We would like the offspring(s) to be „close“ to the parents, but be different from them. Since we are in the real numbers, we can use the power of probability and random sampling! One of the simplest real-valued crossover is a so-called arithmetic crossover, which is basically a weighted mean of two (or more!) parents with randomly chosen weights. Or more formally:

$$ \mathbf{x}^{Offs} = \sum_{i=1}^{N_{Pa}} \alpha_i \mathbf{x}_i $$ where \(\sum_{i=1}^{N_{Pa}} \alpha_i = 1, \alpha_i \in (0, 1)\); \(\alpha_i\) are picked randomly for each offspring.

And there are many other crossover operators, e.g. BLX-α, SBX, PCX...

The mutation is easier. One of the most used way is to just add a randomly generated number (or vector if our genotype is more numbers) from a distribution of our choice. The most used distributions are the uniform distribution (where we need to specify the bounds) and the normal distribution (where we need to specify the variance).

Next...

The final part of this introductory mini-series, which will finally get us to what I'm actually doing, is going to be big. It is going to be about Genetic Programming. So stay tuned!

Evolution in a Computer

As I already mentioned in The First post, my main interest is the evolutionary computation, so I decided I'll write about some basics of it.

Optimisation

Because evolution is in fact a sort of optimisation, what that one actually is? Generally, optimisation is a process of developing or tuning a solution to some problem such that the solution is as good as we can get.

Example: suppose we want to find the minimum (e.g. because it is the money we have to spend) of the following function:

$$ f(x) = \mathrm{e}^x\left(x^3 - 3x^2 + 5x − 5\right) + 6 $$

In this case it is super-easy - we can just compute the derivative

$$ f'(x) = \mathrm{e}^x x \left(x + 1\right) \left(x - 1\right) $$

and then find the points where it equals zero

$$ x = \left\lbrace \begin{matrix} -1 \\ 0 \\ 1 \end{matrix} \right. $$

and from them pick the one with the lowest value

$$ f(-1) \approx 0.849688 \quad f(0) = 1 \quad f(1) \approx 0.563436 $$

... and the winner is... \(x = 1\)!

Black-Box Optimisation

The previous example was nice, because we had the analytical description of the function (and it was very "user-friendly"). But what if \(f(x)\) is not differentiable (e.g. it's discrete)? Or what if \(f(x)\) is a black-box and we know nothing about it?

We might want to systematically try some points, e.g. in a grid, or pick them randomly and try luck... but what if the process of getting \(f(x)\) is expensive? What if \(x\) is a set of parameters of an aerodynamical body and evaluating \(f(x)\) means physically manufacturing the body and testing it in a wind tunnel? It is clear that we have to find some clever way of exploring the space of \(x\).

Evolutionary Algorithms

Believe it or not, natural evolution is a process of optimisation. What is the \(x\) and \(f(x)\) then? Well, in nature, \(x\) is the DNA that determines the organism (very simplified) and \(f(x)\) is the success of the particular organism of surviving and reproducing in the given environment. So natural evolution optimizes the organisms to surivive and reproduce well in the environment; it seeks the "solution" to the life in the environment.

Nature is, generally, quite successful in finding good solutions. We mimick nature all the time (insect-like robots, bat-like sonar...) and evolutionary algorithms (EAs) are no exception.

The general structure of an EA is following:

initialize a population of solutions \(P\)
while not stopping condition
    \(\forall x \in P:\) compute \(f(x)\)
    select parents from \(P\)
    generate children from parents
    mutate children
    put children back in \(P\)
    possibly remove (kill) some individuals in \(P\)
return the individual with the best \(f(x)\) encountered

Genetic Algorithm

Genetic Algorithm (GA) is one particular type of EA. There are others but GA is probably the most famous and the one being extended and based upon the most. So how does the GA fit to the EA structure described above?

In GA, the genotype of a solution, or \(x\), is a string of 0s and 1s of predefined fixed length (depends on the problem we are solving).

Genotype, phenotype

Genotype is the representation of the solution in the evolution. In natural evolution it's the DNA. In GA it's the binary string.

Phenotype is the "thing" encoded by the genotype, the solution itself. In natural evolution it's the organism. In GA it's some solution to the problem we are solving, e.g. that aerodynamical body I mentioned earlier.

Genotype-phenotype mapping is the process of translating the genotype to a phenotype. Don't ask me to explain what does this mean in nature, I'm not a biologist :). Properties and aspects of this mapping may greatly influence the performance of GA (or any EA, in fact).

The population is a fixed-sized collection of solutions (the right word would be genotype but most of the time the difference is not important so I'll use solution or individual).

Selection is the crucial part of the algorithm. That's the part "where it all happens". Selection's job is to pick such individuals from \(P\) which are better than the others but in a non-greedy way. One of the most widely used selection strategies is the Tournament selection.

Generating children from parets is done by crossover. In general, crossover is a process of combining the genotypes of two (or more) parents to create a genotype of one (or more) child. One of the most simple (and also the most widely used) crossover strategy is the One-point (or single-point) crossover.

Mutation is the second crutial part of the algorithm (together with crossover). It's the part where new genetic material is created. Mutation consists of a small random change in the genotype. The easiest and the most simple mutation is just to randomly pick one (or more) bits and flip their values.

The last two steps, that are putting children back into \(P\) and killing some individuals in \(P\), are usually bound together. In a so-called generational GA, the process of selecting parents, crossing them over to children and mutating the children is repeated until there is a number of children equal to the size of \(P\) and then the \(P\) is discarded and replaced with the generated children. On the other hand, in a so-called steady-state GA, only a few parents are selected (usualy just enough so that crossover can be performed) producing only a few children (defined by the crossover method). These few children are mutated and then they are put back into \(P\) such that the size of \(P\) remains the same. This means that some individuals must be killed. The strategy of chosing those individuals is called replacement strategy and has to be chosen. A very simple strategy might be just killing the worst individual(s).

Both the crossover and mutation are usually applied probabilistically, i.e. there is a predefined probability (chosen by the user) of two parents being subject to crossover (if they are not, they are usually just copied so that the children are identical to their parents) and another probability (also chosen by the user) of a child being subject to mutation. Crossover probability is usually rather high (0.8 is a usual value) while the mutation probability is usually very low (e.g. 0.01 per bit).

Final words

GA is definitely the only EA in town. I also only scratched the surface. The Evolution in a Computer is probably going to become a little mini-series covering some aspects of evolutionary computation. In the next post I'll write about other EAs which are more suitable for continuous problems. And also await a post about evolutionary programming! But it is really too early for that one.

---

This post uses JSXGraph for rendering the function graph, MathJax for displaying math, my fork of David Stutz's jquery-pseudocode for rendering (pseudo)code and some my CSS and JavaScript snippets which you can also find on GitHub.