This is about estimating the *mean first passage time* (MFPT) using
random walks.

I’ve been using random walks to model behaviour of evolutionary
algorithms. Each individual is a state, and the probability of
mutating from individual `u`

to individual `v`

is the transition
probability `p(v|u)`

. Given a particular evolutionary operator, say
subtree mutation, transition probabilities can be calculated exactly
and collected in a transition matrix for all `u`

and `v`

.

One of the properties of most interest to me is the mean first passage
time, also known as the hitting time, from `u`

to `v`

. It’s the mean
number of steps you expect to have to wait, if you start at `u`

and
random-walk until you hit `v`

. It’s interesting because it’s a natural
way of quantifying the distance, in the fitness landscape, between
evolutionary individuals.

It’s kind of amazing that, given a transition matrix, you can
calculate the mean first passage time matrix efficiently in closed
form. Kemeny and Snell gave a method in *Finite Markov Chains*, 1976.
A Python implementation (not written by me) is
here.

However, that can only work where the number of states `n < 4000`

or
so – depending on your available RAM – simply because you need to
store the `n x n`

transition matrix in memory, together with a few
other matrices, and finally the MFPT matrix of the same size. In my
current work, I’m looking at spaces with many billions of elements, so
that is completely out. We might not even know the names of all the
states, or how many there are. Maybe the space could even be infinite?
But don’t hold me to that, because maybe MFPT is undefined then.
Anyway, there is a subset of states in which we’re interested, and for
all pairs `(u, v)`

in that subset, I want to be able to estimate the
MFPT.

So, the natural thing to do is to simulate some random walks and see
what happens. This is possible because instead of storing a transition
matrix, we just have to have a transition function – given a state,
it will transition to a new state. In an evolutionary algorithm, the
transition function is just mutation. In many interesting algorithms,
such as tree-based genetic programming, mutation is non-uniform – we
are more likely to mutate to certain individuals than others. We
iterate this function, counting the number of steps: if we start at
`u`

and stop when we hit `v`

, the number of steps taken is a sample of
`MFPT(u, v)`

. We can then repeat the whole thing to get another
sample. Writing this algorithm is very easy.

However, running it might take a long time. Each walk might take
billions of time-steps, and we’ll need to collect quite a few samples
of the length in order for our estimate to be reliable. Remember we
are interested in multiple pairs `(u, v)`

. If we only think about two
states at a time, we can’t collect information on all those other
states of interest even though we might be going through them. So this
way of running the random walk is very inefficient.

We can do a lot better if use a single long walk and some book-keeping to collect samples for many pairs at the same time. We start walking, and every time we hit one of those states of interest, we record how long it took to get here from the previous sighting of the other states of interest. It’s like running multiple interleaved walks at once. This is much more efficient, and it sounds pretty simple: just record the time-step of the last occurrence of each state of interest, and when you hit a state of interest, subtract the last occurrence..?

Actually, it’s not that simple. There are a couple of tricky issues which could cause bad estimates, and that’s why I want to document my algorithm here.

First, consider this situation. We start out, hit `u`

at time-step 100,
then walk for a while, hitting `u`

again at time-step 200, and so on.

`others -> u (100) -> others -> u (200) -> others -> v (300) -> others -> v (400)`

What happens when we hit the `v`

at time-step 300? Should we record a
sample for `MFPT(u, v)`

of 100? 200? Or both? The answer is 200, only.
We started at `u`

at time-step 100, and even though we arrive there
again at time-step 200, for the purposes of sampling `MFPT(u, v)`

that
is not the start of a new walk. So instead of recording the last
sighting of `u`

, we need to record when the walk from `u`

to `v`

started. In other words, we need a variable `rw_started(u, v)`

, which
records the time-step at which the current walk from `u`

to `v`

started. If -1, that means we are not currently in a walk from `u`

to
`v`

. When we hit `v`

, we regard it as the end of a walk from `u`

to
`v`

(and thus save a sample) only if `rw_started(u, v)`

is not -1. We
then update `rw_started(u, v)`

to -1, and `rw_started(v, u)`

to the
current time-step, to say that a new walk from `v`

to `u`

is now
starting.

Similarly, what should we do when we hit `v`

at time-step 400? Save a
sample 200? 300? Both? This time the answer is to save *neither*. At
time-step 400, we are in the middle of a walk from `v`

to `u`

. We’re
not in a walk from `u`

to `v`

.

However, at time-step 400 we *do* save a sample of 100 for ```
MFPT(v,
v)
```

, ie the random walk from `v`

to itself.

So here is what should happen at each step:

Time step 100: set `rw_started(u, v) = 100`

, set `rw_started(u, u) = 100`

Time step 200: save a sample `MFPT(u, u) = 100`

, set `rw_started(u, u) = 200`

Time step 300: save a sample `MFPT(u, v) = 200`

, set
`rw_started(v, u) = 300`

, set `rw_started(u, v) = -1`

, set
`rw_started(v, v) = 300`

Time step 400: save a sample `MFPT(v, v) = 100`

, set `rw_started(v, v) = 400`

One more issue. After running for a while, hopefully we’ll have saved
some samples. Can we trust them? There is a danger of a systematic
under-estimation. Why? Suppose our walk is of length, say, 1 billion
steps, and there are, say, 10 billion states. Then there must be some
pairs whose true MFPT is longer than the length of our walk. For those
pairs, the only samples we *can* collect are those which are
underestimates of the MFPT. If we get lucky and collect some of those
samples, they can only be underestimates. So to be safe, it’s best to
use the results only for pairs where there are quite a few samples. We
then hope that this happened because these are pairs whose true MFPT
is less than the length of the walk.

The code below carries out this algorithm. Given about 50 steps, it’s already getting pretty good estimates for MFPT between all pairs.

Download this Python code here. A Java version, somewhat coupled to the surrounding code, is available here.