Random Sampling from a Discrete Distribution

I can't recommend enough the article "Darts, Dice, and Coins: Sampling from a Discrete Distribution". It develops Vose Alias Algorithm step by step from a series of algorithms with increasing performance, proving the correctness of each algorithm along the way. I feel like I have a much better intuition about this awesome algorithm now, and the pictures and examples went a long way in helping along the way.

Good stuff.

Power vs Safety

I heard one time the idea that one should use the least powerful tool powerful enough to solve a problem. I didn't understand this at the time, but I've recently come to an understanding (not necessarily the meaning of the person who wrote it).

My first interpretation of this idea was that one should use the tool (I'm thinking in particular about programming languages, but this is more general then just that) with the fewest features. I thought I would rather use the *most* powerful tool I could find. I can understand, for example, not embedding a full programming language into a program if all you want is to let users enter arithmetic expressions, but not, for example, using BASIC rather than Haskell if its an option. You can benefit from simplicity, but I don't see an advantage to using languages without useful features just because you don't strictly *need* them.

Recently I came to another interpretation of this idea. It is common in computer science for there to be a series of systems of varying power that can solve increasingly difficult problems. The catch is that the more powerful the systems the weaker the statements one can make about them. This is true in automata theory (FSM < DFM < Pushdown Automata < Nondeterministic Pushdown Automata < Turing Machine, (although DFM < FSM as well)), the typed lambda calculus (the hierarchy is less obvious (at least to me), but much research goes into coming up with powerful systems (the calculus of constructions, for example) that are still weaker than the untyped lambda calculus to retain properties such as strong normalization or confluence), optimization theory (where you classify problems by the weakest form of optimization that can solve them), recursion theory (primitive recursion < mu-recursion (that about all I know about recursion theory)), abstract grammars (regular grammars < context free grammars < context sensitive grammar < universal grammar). Each of these situations (and more, I'm sure) have many more levels between the ones I happen to have heard of.

Taking automata as an example, each level in that hierarchy can solve all problems solvable by ones lower down, but can't solve all problems solvable by ones at higher levels. However, at each level we lose some assurances about termination or resource use. Finite automata are guaranteed to terminate, but the more powerful Turing Machines are not. Expressions typeable in the polymorphic lambda calculus reduce to a normal form, where untyped expressions do not always have such a form. Equality of regular grammars is decidable, but not so for context free grammars in general. We won't have strange optimization techniques like Genetic Algorithms if every problem were linear.

This is really a very common situation, where we would like to determine the least powerful system powerful enough to solve our problem simply so that we can make as many assurances about the solution as possible. This is also one of the ways that things like closure operators and adjunctions come into computer science ("the least thing great enough to" or "the greatest thing less then"). I think it may even be a fundamental part of computation in some way that I don't know how to formalize. Regardless, it is certainly true that a great deal of effort goes into doing exactly this- using the least powerful tool powerful enough to solve a problem.

The R Language for Machine Learning

The other day I bought a book called "Machine Learning for Hackers," which uses R to show off some machine learning algorithms. There is a lot I don't know about machine learning, as I've only studied a subfield and not really the main statistical methods that are useful for the vast majority of problems. One reason I bought this book (in addition to having lots of gift cards) is that a friend of mine mentioned learning R recently, and I would like to know more about it. I don't know almost anything about R, so beware this post- its my initial impressions and misunderstandings.

My understanding is that R has a hodgepodge of properties that make it a bit strange. The basic data type is a vector (which is pretty cool and probably allows an efficient implementation of a lot of primitive operators on vectors), all data seems to have metadata attached, it is dynamically typed (properties of data are determined by metadata modifiable at runtime), there is an object system build on that metadata, it has higher-order function (terribly important, if you ask me), it is lazy (although common operations seem to force evaluation), its call-by-value in the sense of copying (although its looks like copying is done only when necessary), and it allows mutation of structures (although it limits propagation of effects on data through copying).

Overall its probably a good language for its domain- it got a lot attention to convenience, its designed to manipulate lots of complex data, and it supports some functional programming. I expect that it could get very messy with edge cases, odd interactions between mechanisms, possibly clumsy support for more regular programming tasks (I don't know either way about this one), and undefined semantics. The last one is a shame- languages should really know what their semantics or its programmings won't know what their programs really mean, and its hard to come up with variations of the language.

I'm pretty excited to learn these techniques, although this book does not go deep into theory. Its more about getting interested and involved in these algorithms with case studies and practical concerns. I think this is perfect- I'll learn what theory I want to know somewhere else. For now, I'm just glad to get a feel for the broader reach of this field that I like so much.

Genetic Algorithm and RGEP Complexity

I've been looking at efficient implementations of genetic operators lately. Obviously the actual complexity depends on the fitness function (sometimes more then anything else), the choice of genetic operators, the implementation of those operators, and the representation of the individuals. However, I'm really interested in the ones used in RGEP, so point mutation, one and two point crossover, rotation, and tournament selection.

I'm slowly converging on a set of structures and algorithms that seem to give a complexity on the order of O(g*p*log(n)) with g the number of generations, p the population size, and n the size of the individuals. This also applies to RGEP, as it is pretty much just a genetic algorithm with a Rotation operator. A genetic algorithm would normally have O(g*p*n) instead, assuming the usual genetic operators, the usual array of arrays representation, and a reasonable implementation of the genetic operators.

To get this time complexity, I plan on using a Vector (the Haskell library Data.Vector) for the population (assuming the usual multiset population) and a Sequence (the Haskell library Data.Sequence, which implements 2-3 finger trees) for the individuals. These are chosen because Vectors have efficient indexing and a nice function "backpermute" which can be used to select the individuals for a new population (given their indices from a selection operator) in O(n) time. Sequences have an amortized O(min(log(i, n-i))) split operator where i is the index to split on and n is the size of the sequence. They also have a O(min(n,m)) concatenation operation (n and m are the sizes of the two sequences). This is important as all of point mutation, crossover, and rotation are really not about copying, but about splitting and concatenating, and should really take less then linear time.

I have learned that roulette wheel selection can be done in O(n) time, sampling n times from a population where each sample can take O(1) time using either the algorithm in "Roulette-Wheel Selection via Stochastic Acceptance", or in my case maybe the one in "A Linear Algorithm for Generating Random Numbers with a Given Distribution". Even using Tournament Selection, which is my personal preference, we need only linear time in the population size to do selection.

I've talked before about point mutation by sampling from a geometric distribution, and given that we need only set up the distribution once and each sample will take O(1) time depends only on the number of mutations in the population. This is usually chosen so that each individual is mutated only a couple of times, so it occurs with about the same frequency regardless of the individuals size (and therefore shouldn't depend on that size). It is unfortunate that this is the one place where a vector would be better than a sequence (assuming it is a mutable vector), but I'm not terribly worried about this. If necessary, all the points can be either mutated in a single pass down the whole population, only changing those points that were determined to be mutation points, or by a series of applications of the function "adjust", which is where I'm getting what I'm going to call p*logn complexity because the number of mutations is determined by the size of the population and the complexity of splitting and rejoining an individual at a given point.

Given these observations, the time complexity is O(g * (mut + cross + rot + select)) where the complex of mut, cross, rot, and select are plogn, plogn, plogn, and n respectively. This makes the total complexity only O(g*p*logn).

This may not be the best we can get (I'm not sure, but I'm looking in to it) but its pretty dang good compared to my usually naive implementations.

Efficient Roulette Wheel Selection

I ran into a pretty nice way of doing roulette wheel selection the other day. It requires O(n) set up time, and each selection takes O(log(n)) time, for a total of O(n+nlog(n)) or O(nlog(n)) time for the whole selection process. This is a lot better than the completely naive way (generating the "spin" of the wheel, and then running through the individuals looking for the one that was "hit") which takes O(n^2), and which I've been guilty of in the past.

First off is the setup. Given a population of individuals, paired with their fitnesses. Lets assume that this is in some linear structure like a linked list or vector. Go through each individual, setting their fitness value to the total sum of all previous individuals plus their own fitness. The summing of the fitnesses takes time linear in the population size.

Next is doing the actual selection. This is amazingly easy- generate a value between 0 and the sum of all fitness values. Then do a binary search through the population with the new (summed) fitness values. Each selection takes log(n) time (to do the binary search) and you need to do it n times, so the second part of this selection technique takes nlog(n) time.

So there you have it- an easy and quick O(nlog(n)) time implementation of roulette wheel selection. I will probably end up using this soon enough.

Haskell Evolutionary Algorithm Library, take 2

I've been working on rewriting the program I did for my thesis- an evolutionary algorithm library in Haskell. It should be more general then just evolutionary algorithms, but I can't think of a better name involving Haskell and Machine Learning. This post is a progress report.

So far its going great. The structure of the library is a mixture of functional lenses (from Data-Lens), pipes (from Conduits), random variables/probability distributions (from RVar), and a little monad for keeping state and doing things randomly while allowing mutable references. I'm not completely sold up on the mutable references thing (ST monad) for this situation but I would like to play around with using it for all operators and seeing if the efficiency is worth the extra complexity and type variable in all my definitions.

The plan is to provide a couple of implementations of the usual genetic operators- a very general one that should work with different structure nicely and can be easily used, a very efficient one that can be used on fewer structures (probably vectors of unboxed types) but is hopefully insanely fast compared to what I've done before (see my previous posts), and one thats mostly for fun where everything is done entirely in terms of pipes. The last one would probably be slow, but it would be interesting to see a GA defined not in terms of populations and individuals, but completely in terms of pipes of elements (bits for example).

One of my major ease-of-use goals is that the user should be able to run a Genetic Algorithm simply by providing a fitness function. It should be easy to change parameters and all that, and slightly less easy to add operators or define a new algorithm, but for a quick run the fitness function should be it. I don't just mean a fitness function of a fixed type- I would like the generic GA to be so generic that you can give it a fitness function that expects lists of Bools and get lists of Bools, or if you want vectors of Ints, you get vectors of Ints, etc.

The current next step is finishing the implementation of point mutation and crossover, and then starting on some selections. I've notice that the pipe concept (and the nice implementation from Conduits) vastly simplifies this sort of algorithm. I really believe it is the Right Thing for the job.

Efficient Crossover and Selection for Genetic Algorithms

In the last post I discussed an efficient implementation of point mutation for a Genetic Algorithm. This time we look at crossover and selection (and rotation, but thats for PGEP and RGEP, not a GA). These operators take linear time each, and we do no better than that here. However, we can reduce the coefficient nicely, which is important in practice.

First off, lets look at crossover- especially 1 or 2 point crossover (or both). Crossover is not a terribly expensive operator- we simply need to generate a random value between 0 and 1 (to determine if we will cross a pair of individuals), and then a cross point between 0 and the length of the individual. The only somewhat expensive thing here is moving pieces of the individuals around, so lets not do this yet. Instead, we just make a list/array/vector of the ranges each individual occupies. One way to do this would be to consider the population a single large vector (or a 2 dimensional vector) and each individual a list of ranges from their start to end. After 1 point crossover an individual is a 2 element list of pairs- start to cross point, other side of crosspoint to end. In other words, we record what bits in the population make up each individual without actually changing the population at all. If we do 2 point crossover, or multiple crossovers in a single generation, these ranges may consist of more pairs of locations, but are still basically the same. We don't want the individuals to become too fragmented like this, as with completely fragmented individual this would be less efficient then the usual implementation, but for a single generation fragmentation shouldn't be too bad.

Selection can be similarly performed without changing the population. We just record the winner of each tournament or the result of spins of the roulette wheel (or whatever selection one desires). One way to do this is an initially all 0 vector, counting the number of occurrences of each individual in the new population. We add one to a location if the individual corresponding to that location is selected.

While this post is mainly about Genetic Algorithms, I want to mention that the rotation operator from PGEP and RGEP can be performed in a similar fashion. We either record the number of locations to rotate for each individual (with 0 for individuals that aren't rotated that generation) or modify the ranges we created during crossover described above.

So, we have an untouched population, a vector of ranges for each individual that describes what parts of the population they occupy, and a count saying how many times they occur in the next generation. It might be nice to be able to do the operators this way several times before touching the population (except mutation, which can be done independently however we desire), but since we need the full individual for fitness evaluation anyway, we might as well perform the crossover and selection we recorded at the end of each generation.

Now all thats left is to copy each individual from the old generation to the new one in one single pass over the population. This involves inspecting each individual and copying them the number of times they were selected into successive places in the new population (I'm thinking a double buffering style scheme where you copy back and forth between two equal sized populations). Each individual consists of a series of ranges, so getting the individual to the new population means copying each range one after the next to successive locations. An individual [(0, 10), (103, 107), (23, 26)] would copy the values between the locations (assuming the population is considered one big vector) 0 and 10, 103 and 107, and 23 to 26 to a single vector from, say, 0 to 20.

The advantage here is that while each copy is somewhat complex, we end up doing a small amount of work to determine what actions we need to take are (recording how to do each operator instead of doing it) which is sparse data about rearranging the individuals genetic material, and then actually copying each location only once per generation. I'm hoping to try this out when I get around to actually implementing this evolutionary algorithm library, although I will provide the simpler implementations as well.