Do The Evolution

| 0 Comments

In the last few posts[1][2] we have taken a look at genetic algorithms, which use simple models of biological evolution to search for global maxima of functions, being those points at which they return their greatest possible values.
These models typically represent the arguments of the function as genes within the binary chromosomes of individuals whose fitnesses are the values of the function for those arguments, exchange genetic information between them with a crossover operator, make small random changes to them with a mutation operator and, most importantly, favour the fitter individuals in the population for reproduction into the next generation with a selection operator.
We used a theoretical analysis of a simple genetic algorithm[3] to suggest improved versions of the crossover operator, as well as proposing more robust schemes for selection and the genetic encoding of the parameters.
In this post we shall use some of them to implement a genetic algorithm for the ak library.

Firstly, we shall use the uniform crossover operator which swaps the bits between a pair of chromosomes at each position, or locus, with a fifty percent probability.
Secondly, we shall use tournament selection which picks two individuals at random and copies the better of them into the next generation, neatly avoiding any direct dependence upon the actual values of their fitnesses.
Finally, we shall use Gray code[4], in which consecutive integers differ at only one bit, to map a chromosome's genes to the function's parameters, ensuring that the mutation operator has a realistic chance of making a small change to them.

The Chromosome And Population

So far we have been representing chromosomes with arrays of bits which, given that we have 32 bits to play with when using bitwise operators on JavaScript's numbers, is a shocking waste of memory. Much better then to treat those numbers as "arrays" of 32 bits and represent the chromosome as an array of arrays of bits. We shall, of course, need to use bitwise operators to access the individual bits, but that's not such a big deal really, as illustrated by listing 1 which uses a random bit generator like our ak.wardMoltenoRnd, for example, to fill the bits of such a chromosome.

Listing 1: The Chromosome
function chromosome(length, bitRnd) {
 var n = ak.ceil(length/32);
 var bits = new Array(n);
 var pos = 0;
 var bit = 1;
 var i;

 for(i=0;i<n;++i) bits[i] = 0;

 while(length-->0) {
  if(bitRnd()!==0) bits[pos] |= bit;
  if((bit<<=1)===0) {bit=1; ++pos;}
 }
 return {bits: bits, fitness: undefined};
}

Here pos indicates which array of bits we're currently initialising and bit indicates which bit within it we are. Note that we're relying upon the left shift operator throwing away the most significant bit and setting the least significant bit to zero so that bit will equal zero once we have shifted it 32 times from its initial value of one, indicating that we must reset it to one and increment pos.
Finally, we initialise the fitness with undefined to indicate that it hasn't been evaluated.

Now this does mean that it's possible for a gene to be comprised of bits contained in more than one number, although if we don't allow them to have more than 32 bits each they can be located in at most two numbers. Even so, extracting them is a slightly fiddly business, as demonstrated by listing 2.

Listing 2: Extracting A Gene And Converting It To An Argument
function uniArg(bits, lb, ub, locus, length) {
 var pos, from, to, u, mask;

 if(length===0) return lb;

 pos = ak.floor(locus/32);
 from = locus%32;
 to = (locus+length-1)%32;

 if(from+length<=32) {
  u = (bits[pos]<<(31-to))>>>(32-length);
 }
 else {
  u = (bits[pos]>>>from) | ((bits[pos+1]&(0xffffffff>>>(31-to)))<<(32-from));
 }

 u ^= u>>>1;
 u ^= u>>>2;
 u ^= u>>>4;
 u ^= u>>>8;
 u ^= u>>>16;
 if(u<0) u += 0x100000000;

 u /= Math.pow(2, length);
 return lb*(1-u) + ub*u;
}

If a gene resides entirely within a single number, we first shift any bits after its last out past the most significant bit, discarding them, after which we shift it back so that its first bit is the least significant bit.
Otherwise, we shift the gene's first part down so that its first bit is the least significant, then mask off any bits after the end of its second part before shifting it up to the bit after the last in the first part and combining them with a bitwise or.
That minor complication dealt with, we convert from Gray code to binary using the 32 bit scheme that we proved works in the last post, scale the result to a number greater than or equal to zero and strictly less than one and finally scale it again onto one between the lower and upper bounds lb and ub.

Listing 3 gives the multiArg function which converts a set of genes into a vector argument for a multivariate function, represented by our ak.vector type.

Listing 3: Converting All Of The Genes To A Vector Argument
function multiArg(bits, lb, ub, lengths) {
 var n = lb.length;
 var locus = 0;
 var arg = function(i) {
  var a = uniArg(bits, lb[i], ub[i], locus, lengths[i]);
  locus += lengths[i];
  return a;
 };
 return ak.vector(n, arg);
}

To initialise a population we simply copy freshly created individuals into an array, as shown by listing 4.

Listing 4: Initialising The Population
function populate(pop, length) {
 var n = pop.length;
 var bitRnd = ak.wardMoltenoRnd();
 var i;

 for(i=0;i<n;++i) pop[i] = chromosome(length, bitRnd);
}

When it comes to evaluating the fitness of an individual, we must distinguish between univariate fitness functions, which take a number as an argument, and multivariate fitness functions, which take a vector.
The uniEval and multiEval functions, given in listing 5, create functions that convert individuals' chromosomes into number or vector arguments, as appropriate, and sets their fitnesses to the result of the fitness function f for those arguments.

Listing 5: Evaluating Individuals
function uniEval(f, lb, ub, length) {
 return function(c) {
  var arg = uniArg(c.bits, lb, ub, 0, length);
  c.fitness = f(arg);
 };
}

function multiEval(f, lb, ub, lengths) {
 return function(c) {
  var arg = multiArg(c.bits, lb, ub, lengths);
  c.fitness = f(arg);
 };
}

The Crossover And Mutation Operators

To implement the crossover and mutation operators we'll need to swap bits between and flip bits within chromosomes, which the swapBit and flipBit functions given in listing 6 do with some relatively simple bit twiddling.

Listing 6: Swapping And Flipping Bits
function swapBit(bits0, bits1, locus) {
 var pos = ak.floor(locus/32);
 var bit = 1 << locus%32;
 var tmp = bits0[pos];

 bits0[pos] &= ~bit;
 bits0[pos] |= bits1[pos] & bit;
 bits1[pos] &= ~bit;
 bits1[pos] |= tmp & bit;
}

function flipBit(bits, locus) {
 var pos = ak.floor(locus/32);
 var bit = 1 << locus%32;
 bits[pos] ^= bit;
}

In the crossover operator we simply iterate over pairs of individuals in the proportion of the population chosen for crossover and swap each of their bits with a probability of one half. Note that we don't need to randomise the population since both its initialisation and selection into subsequent generations will already have done so. This is implemented by the crossover function given in listing 7 for a proportion rate of the population.

Listing 7: Crossover
function crossover(pop, length, rate, rnd) {
 var n = pop.length;
 var count = ak.floor(n*rate);
 var i, j;

 for(i=1;i<count;i+=2) {
  for(j=0;j<length;++j) {
   if(rnd()<0.5) swapBit(pop[i-1].bits, pop[i].bits, j);
  }
  pop[i-1].fitness = undefined;
  pop[i].fitness = undefined;
 }
}

Listing 8 gives an implementation of the mutation operator which simply flips a proportion rate of the bits in the population, choosing each at random.

Listing 8: Mutation
function mutation(pop, length, rate, rnd) {
 var n = pop.length;
 var count = ak.floor(n*length*rate);
 var i, j, locus;

 for(i=0;i<count;++i) {
  j = ak.floor(rnd()*n);
  locus = ak.floor(rnd()*length);
  flipBit(pop[j].bits, locus);
  pop[j].fitness = undefined;
 }
}

Note that in both of these function, the affected individuals have their fitnesses set to undefined to indicate they they need to by re-evaluated.

The Evaluation And Selection Operators

When evaluating the population we shall also keep track of the best individual that we have seen throughout the algorithm's operation. Since we shall be using tournament selection, which simply compares the fitnesses of two individuals to decide which should be copied into the next generation, we can make our genetic algorithm a function minimiser by defining better to mean lesser, as done in listing 9.

Listing 9: Evaluation
function evaluation(pop, best, evaluate) {
 var n = pop.length;
 var i;

 for(i=0;i<n;++i) {
  if(ak.nativeType(pop[i].fitness)===ak.UNDEFINED_T) {
   evaluate(pop[i]);
   if(isNaN(pop[i].fitness)) pop[i].fitness = ak.INFINITY;

   if(pop[i].fitness<best.fitness) {
    best.bits = pop[i].bits.slice(0);
    best.fitness = pop[i].fitness;
   }
  }
 }
}

Here the evaluate argument will be a function created by either uniEval or multiEval, depending on whether the fitness function takes a single argument of a vector of them.
Note that we ensure that we shall aggressively select against individuals with NaN fitnesses by setting them to infinity.

To give us some control over the rate at which tournament selection drives the algorithm to convergence we shall introduce a measure of selection pressure by giving the worse of the pair of individuals a chance of being selected over the better of them, as done in listing 10, allowing the algorithm to spend more time exploring the search space.

Listing 10: Selection
function selection(pop, pressure, rnd) {
 var n = pop.length;
 var tmp = new Array(n);
 var i, i0, i1, tmp;

 for(i=0;i<n;++i) {
  i0 = ak.floor(rnd()*n);
  i1 = ak.floor(rnd()*n);

  if(rnd()<pressure) {
   tmp[i] = pop[i0].fitness<pop[i1].fitness ? pop[i0] : pop[i1];
  }
  else {
   tmp[i] = pop[i0].fitness<pop[i1].fitness ? pop[i1] : pop[i0];
  }
 }

 for(i=0;i<n;++i) {
  pop[i].bits = tmp[i].bits.slice(0);
  pop[i].fitness = tmp[i].fitness;
 }
}

Note that if pressure equals one half then the individuals will be chosen entirely at random and so there will be no selection pressure at all, whereas if it equals one then we shall always choose the better of them.

Finally, the initialise and generation functions given in listing 11 set up and evaluate the first population and apply the crossover, mutation, evaluation and selection operators to create subsequent populations respectively.

Listing 11: Initialisation And Generation
function initialise(pop, best, length, evaluate) {
 populate(pop, length);
 evaluate(pop[0]);
 if(isNaN(pop[0].fitness)) pop[0].fitness = ak.INFINITY;

 best.bits = pop[0].bits.slice(0);
 best.fitness = pop[0].fitness;
 evaluation(pop, best, evaluate);
}

function generation(pop, best, length, cross, mutate, pressure, evaluate, rnd) {
 crossover(pop, length, cross, rnd);
 mutation(pop, length, mutate, rnd);
 evaluation(pop, best, evaluate);
 selection(pop, pressure, rnd);
}

ak.geneticMinimum

We're finally ready to implement our genetic algorithm for function minimisation, given in listing 12.

Listing 12: ak.geneticMinimum
ak.geneticMinimum = function(f, size, cross, mutate, pressure, steps, rnd) {
 if(ak.nativeType(f)!==ak.FUNCTION_T) {
  throw new Error('invalid function in ak.geneticMinimum');
 }

 size = ak.nativeType(size)===ak.UNDEFINED_T ? 50 : ak.floor(Math.abs(size));
 if(!isFinite(size)) {
  throw new Error('invalid population size in ak.geneticMinimum');
 }

 if(ak.nativeType(cross)===ak.UNDEFINED_T) {
  cross = 0.6;
 }
 else if(ak.nativeType(cross)!==ak.NUMBER_T || !(cross>=0 && cross<=1)) {
  throw new Error('invalid crossover rate in ak.geneticMinimum');
 }

 if(ak.nativeType(mutate)===ak.UNDEFINED_T) {
  mutate = 0.01;
 }
 else if(ak.nativeType(mutate)!==ak.NUMBER_T || !(mutate>=0 && mutate<=1)) {
  throw new Error('invalid mutation rate in ak.geneticMinimum');
 }

 if(ak.nativeType(pressure)===ak.UNDEFINED_T) {
  pressure = 0.75;
 }
 else if(ak.nativeType(pressure)!==ak.NUMBER_T || !(pressure>=0 && pressure<=1)) {
  throw new Error('invalid selection pressure in ak.geneticMinimum');
 }

 steps = ak.nativeType(steps)===ak.UNDEFINED_T ? 100 : ak.floor(Math.abs(steps));
 if(!isFinite(steps)) {
  throw new Error('invalid number of steps in ak.geneticMinimum');
 }

 if(ak.nativeType(rnd)===ak.UNDEFINED_T) {
  rnd = Math.random;
 }
 else if(ak.nativeType(rnd)!==ak.FUNCTION_T) {
  throw new Error('invalid random number generator in ak.geneticMinimum');
 }

 pressure = 0.5*(1+pressure);

 return function(lb, ub, lengths) {
  var minimum = ak.nativeType(lb)===ak.NUMBER_T ? uniMinimum : multiMinimum;
  return minimum(f, lb, ub, lengths, size, cross, mutate, pressure, steps, rnd);
 };
};

Most of the code is concerned with providing default values for the algorithm's parameters and checking that those that were specified are valid.
Specifically, we have a default population size of fifty, a crossover rate of sixty percent, a mutation rate of one percent, a selection pressure of three quarters and one hundred steps, using Math.random if the user hasn't provided a random number generator of their own.
Note that we transform the pressure argument so that values between zero and one are translated to probabilities between one half and one. This means that zero equates to no selection pressure and one equates to maximal selection pressure.
These details out of the way, we return a function that takes lower bounds, upper bounds and lengths for the individuals' genes and returns the result of either uniMinimum or multiMinimum depending upon whether the lower bound is a number or not.

The uniMinimum function, given in listing 13, first checks that both the lower and upper bounds are finite numbers and that length of the individuals' single gene is a number no greater than 32, simply returning the lower bound if either the size of the population or the length of the gene equal zero. Finally, it creates an array for the population, an object to keep track of the best individual and uses uniEval to create the function to evaluate individuals before initialising the population and executing the generation once per step.

Listing 13: uniMinimum
function uniMinimum(f, lb, ub, length, size, cross, mutate, pressure, steps, rnd) {
 var pop, best, evaluate;

 length = ak.floor(Math.abs(length));

 if(!isFinite(lb)) {
  throw new Error('invalid lower bound in ak.geneticMinimum');
 }

 if(ak.nativeType(ub)!==ak.NUMBER_T || !isFinite(ub)) {
  throw new Error('invalid upper bound in ak.geneticMinimum');
 }

 if(!(length<=32)) {
  throw new Error('invalid gene length in ak.geneticMinimum');
 }

 if(size===0 || length===0) return lb;

 pop = new Array(size);
 best = {};
 evaluate = uniEval(f, lb, ub, length);

 initialise(pop, best, length, evaluate);
 while(steps-->0) {
  generation(pop, best, length, cross, mutate, pressure, evaluate, rnd);
 }
 
 if(!isFinite(best.fitness)) evaluate(best);
 return !isNaN(best.fitness) ? uniArg(best.bits, lb, ub, 0, length) : ak.NaN;
}

Note that if, at the end of the algorithm's operation, the best individual has an infinite fitness we must re-evaluate it to make sure that it wasn't really NaN, converting its bits to a number if it wasn't and returning NaN if it was.

The multiMinimum function, given in listing 14, proceeds in much the same manner, albeit checking that the lower and upper bounds are ak.vector objects having the same non-zero number of dimensions and containing finite values, that the lengths of the genes are contained in an array of the same length and are each no greater than 32, adding them together to yield the length of each individual's chromosome, and using multiEval to create the function to evaluate them.

Listing 14: multiMinimum
function multiMinimum(f, lb, ub, lengths, size, cross, mutate, pressure, steps, rnd) {
 var n, i, length, pop, best, evaluate;

 if(ak.type(lb)!==ak.VECTOR_T || lb.dims()===0) {
  throw new Error('invalid lower bounds in ak.geneticMinimum');
 }

 if(ak.type(ub)!==ak.VECTOR_T || ub.dims()!==lb.dims()) {
  throw new Error('invalid upper bounds in ak.geneticMinimum');
 }

 lb = lb.toArray();
 ub = ub.toArray();

 n = lb.length;

 if(ak.nativeType(lengths)!==ak.ARRAY_T) {
  length = lengths;
  lengths = new Array(n);
  for(i=0;i<n;++i) lengths[i] = length;
 }
 else if(lengths.length!==n) {
  throw new Error('invalid gene lengths in ak.geneticMinimum');
 }

 length = 0;
 for(i=0;i<n;++i) {
  lengths[i] = ak.floor(Math.abs(lengths[i]));
  if(!isFinite(lb[i])) throw new Error('invalid lower bound in ak.geneticMinimum');
  if(!isFinite(ub[i])) throw new Error('invalid upper bound in ak.geneticMinimum');
  if(!(lengths[i]<=32)) throw new Error('invalid gene length in ak.geneticMinimum');
  length += lengths[i];
 }

 if(size===0 || length===0) return ak.vector(lb);

 pop = new Array(size);
 best = {};
 evaluate = multiEval(f, lb, ub, lengths);

 initialise(pop, best, length, evaluate);
 while(steps-->0) {
  generation(pop, best, length, cross, mutate, pressure, evaluate, rnd);
 }

 if(!isFinite(best.fitness)) evaluate(best);
 return !isNaN(best.fitness) ? multiArg(best.bits, lb, ub, lengths)
                             : ak.vector(n, ak.NaN);
}

Finally, this implementation can be found in GeneticMinimum.js.
Figure 1: Rosenbrock's Function

Testing The Implementation

To test our implementation we shall first take a look of its results for Rosenbrock's function, defined as
\[ f(x, y) = 100 (y - x^2)^2 + (1-x)^2 \]
and having a long curved valley that gently slopes to a global minimum of zero when both arguments are equal to one, as illustrated by figure 1.
Recall that our original version of the genetic algorithm would converge somewhere on the valley floor, but frequently some way away from its minimum which wasn't entirely surprising given that it was designed to test local, rather than global, minimisation algorithms.
Program 1 demonstrates how our new version does by plotting the results of running it twenty times with red crosses on the contour plot of Rosenbrock's function.

Program 1: Minimising Rosenbrock's Function

Figure 2: De Jong's f5
The results are slightly better than they were before, concentrated around the minimum and less likely to be very far from it. This is primarily because Gray coding makes it easier for the mutation operator to make small incremental changes to the arguments associated with the individuals' genes.

De Jong's \(f_5\)[5] is a more appropriate test for global optimisation algorithms since it consists of a gently sloping plane of values in the region of 500 with 25 narrow local minima ranging from a value of approximately 24 at the top right to a value of approximately one at the bottom left, as illustrated by figure 2.
We can define it in terms of the constants
\[ \begin{align*} a_{0i} &= -32 + \left(i \,\%\, 5\right) \times 16\\ a_{1i} &= -32 + \bigg\lfloor\frac{i}{5}\bigg\rfloor \times 16 \end{align*} \]
where \(\%\) is the remainder of the term to its left after integer division by the term to its right and the oddly shaped braces stand for the greatest integer less than or equal to the value between them, with
\[ f_5(x_0, x_1) = \left(0.002 + \sum_{i=0}^{24} \frac{1}{i+1 + \left(x_0-a_{0i}\right)^6 + \left(x_1-a_{1i}\right)^6}\right)^{-1} \]
Program 2 demonstrates our genetic algorithm's results for De Jong's \(f_5\) in the same way that program 1 did for Rosenbrock's function.

Program 2: Minimising De Jong's f5

Well it certainly converges upon the global minimum more consistently that our original version of the algorithm did, lending a little weight to my claim that the new versions of the genetic encoding and operators are indeed improvements.

Afterword

Now as much fun as genetic algorithms are to implement, it must be noted that they have even more controlling parameters than simulated annealing[6] which are, if anything, even less directly related to the balance between exploring the search space and concentrating on regions containing the best arguments that it finds.
I would therefore suggest that you turn to our ak.annealMinimum function minimiser in preference to them when searching for global minima.

References

[1] It's All In The Genes, www.thusspakeak.com, 2017.

[2] The Best Laid Schemata, www.thusspakeak.com, 2017.

[3] Holland, J., Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence, University of Michigan, 1975.

[4] Gray, F., Pulse Code Communication, US2632058A, United States Patent and Trademark Office, 1953.

[5] De Jong, K., An analysis of the behaviour of a class of genetic adaptive systems, Ph.D. thesis, University of Michigan, 1975.

[6] Annealing Down, www.thusspakeak.com, 2017.

Leave a comment