A Unified Markov Chain Approach to Analysing Randomised
Search Heuristics

arXiv:1312.2368v1 [math.OC] 9 Dec 2013

Jun He

∗

Feidun He

Xin Yao

December 10, 2013

Abstract
The convergence, convergence rate and expected hitting time play fundamental roles in the analysis
of randomised search heuristics. This paper presents a unified Markov chain approach to studying them.
Using the approach, the sufficient and necessary conditions of convergence in distribution are established.
Then the average convergence rate is introduced to randomised search heuristics and its lower and upper
bounds are derived. Finally, novel average drift analysis and backward drift analysis are proposed for
bounding the expected hitting time. A computational study is also conducted to investigate the convergence, convergence rate and expected hitting time. The theoretical study belongs to a prior and general
study while the computational study belongs to a posterior and case study.

1

Introduction

Randomised search heuristics, such as evolutionary algorithms, have been widely applied to optimization
problems. Randomised search heuristics belong to iterative methods. As iterative methods, the following
three questions are fundamental in both theory and practice.
1. (Convergence) whether is a randomised search heuristic able to find an optimal solution eventually?
2. (Convergence rate) how fast does a randomised search heuristic converge to the optimal set per iteration?
3. (Hitting time) how many iterations are needed for obtaining an optimal solution?
Most randomised search heuristics satisfy the Markov property, that is, a population sequence (where a
population consists of one or more solutions) is generated subject to some probability distribution; and the
state of current population decides the state of next population in a probabilistic way. Hence Markov chain
theory provides a theoretical framework for analysing randomised search heuristics [31, 18],
Based on absorbing Markov chain theory [22, 13], a unified approach is used for studying the convergence,
convergence rate and expected hitting time of randomised search heuristics in this paper. The idea is described
as follows: the population sequence generated by a randomised search heuristic is modelled by an absorbing
Markov chain. Consider the probability distribution of a population in the non-optimal solution set and
represent it by a vector. Then the randomised search is equivalent to a matrix iteration. The vector 1-norm
is chosen to measure the distance between a population and the optimal solution set. This new feature makes
our current analysis different from previous work [30, 34, 15]. Using matrix iteration analysis [35, 24], all
theoretical results can be established in a unified manner.
The purpose of this paper is to seek new theoretical tools for analysing randomised search heuristics.
Indeed we have developed three new tools in the paper, which are novel average drift analysis and novel
backward drift analysis for bounding the expected hitting time; and the average convergence rate of randomised search heuristics. These new tools are seldom studied before.
This paper is organised as follows: literature review is given in Section 2. The Markov chain model
appears in Section 3. Convergence is analysed in Section 4. The average convergence rate of EAs is discussed
in Section 5. The expected hitting time is analysed in Section 6. Final conclusions are described in Section 7.
∗ Department

of Computer Science, Aberystwyth University, Aberystwyth, SY23 3DB, U.K.

1

2

Literature Review

In this section we review the work related to Markov chain analysis for randomised search heuristics, and
also show the difference between our work and previous ones.
Markov chains have been applied into analysing randomised search heuristics more than two decades [12,
10]. [31] gave a survey of Markov chain analysis of evolutionary algorithms up to 1998. [29] reviewed some
achievements after that year.
Markov chain theory provides a theoretical framework to model randomised search heuristics. For example, [27, 4, 5] modelled genetic algorithms by Markov chains. [2] presented Markov chain models of parallel
genetic algorithms. In general, randomised search heuristics for discrete optimisation may be modelled by
Markov chains, and randomised search heuristics for continuous optimisation by Markov processes.
Markov chain theory is widely applied to the limit behaviour of randomised search heuristics and their
convergence. For example, [11, 30, 19] proposed different convergent conditions of genetic algorithms. [32]
analysed convergence properties of some multi-objective evolutionary algorithms. [3] compared various selection algorithms using Markov chains.
The convergence rate of randomised search heuristics is a less studied topic. [34] analysed a simple GA by
evaluating the eigenvalues of the transition matrix of the Markov chain and computed its convergence rate.
[15] gave the convergence rates of general genetic algorithms by using the minorization condition. [9] studied
the convergence rates of gene expression programming by means of Markov chain and spectrum analysis.
The expected hitting time study has received more attentions recently. [16, 8] made two initial discussions
on evolutionary algorithms. According to absorbing Markov chain theory, the expected hitting time can be
calculated based on the fundamental matrix. Using this approach, [17] compared the expected hitting time
of (1 + 1) and (N + N) evolutionary algorithms. [18] gave a framework for analysing the expected hitting
time of evolutionary algorithms. [37] compared the runtime of three simple heuristic algorithms.
However, the fundamental matrix approach can only be suitable for simple algorithms and problems.
Thus [14, 16] introduced drift analysis to the expected hitting time study. Currently drift analysis becomes
a popular theoretical tool. Different variants have been developed, such as simplified drift analysis [28],
multiplicative drift analysis [7], adaptive drift analysis [6], and variable drift [25, 21]. Drift analysis have
been applied to both (1+1) EAs [36] and population-based EAs [23]. Drift theorems can be established using
either Markov chain theory or super-martingale theory [16, 26].
So far the convergence, convergence rate and expected hitting time are studied separately. Different from
existing work, we present a unified approach to bringing these three issues together.

3
3.1

Absorbing Makrov Chain Model
Theoretical Study

In this subsection, we describe absorbing Markov chain models of randomised search heuristics.
Consider a maximisation problem on a finite state space, that is
max f (x),

x ∈ D,

(1)

where x is a variable and D is its definition domain, a finite set. f (x) is called a fitness function.
A randomised search heuristic can be viewed as a randomised iteration process: initially construct a
population of solutions Φ0 ; based on Φ0 , then probabilistically generate a new population of solutions Φ1 ;
based on Φ1 , then probabilistically generate a new population of solutions Φ2 , and so on. This procedure is
repeated until a stopping criterion is satisfied. Then a sequence of populations is produced
Φ0 → Φ1 → Φ2 → · · · .
In order to guarantee that the best solution is always kept during the iteration, an extra archive is added
for recording the best found solution. The archive itself is not involved in generating new solutions. This
strategy is called elitist. A randomised search heuristic with an archive is described in Algorithm 1.
A population consists of one or more solutions while a solution is called an individual. The procedure of
generating new a population may include several steps, such as mutation, crossover and selection in a genetic
2

Algorithm 1 randomised search heuristic with an archive
1: set counter t to 0;
2: initialize a population of solutions Φ0 ;
3: archive Λ0 keeps the best solution in Φ0 ;
4: for t = 0, 1, 2, · · · do
5:
a new population of solutions Φt+1 is generated from Φt ;
6:
update the archive Λt+1 if the best solution in Φt+1 is better than Λt ;
7:
counter t is increased by 1;
8: end for
algorithm. For convenience of analysis, we only consider non-dynamical algorithms whose search operators
are not changed during the iteration. The algorithm runs for ever. This assumption is taken for convenience
of analysis of the hitting time. The fitness of a population at the t-th iteration is defined by the archive at
the t-th iteration, denoted by ft .
The sequence {Φt ; t = 0, 1, · · · } can be formulated by a Markov chain. Let X and Y be two populations.
The transition from X to Y happens with a probability P (X, Y ), that is,
P (X, Y ) := P (Φt+1 = Y | Φt = X),

X, Y ∈ S,

(2)

where S denotes the set of all populations. Both Φt and X represent a population, but Φt is a random
variable for representing the population at the t-th iteration; X its value, or called a state in the population
space.
For the sake of argument, we introduce an auxiliary Markov chain {Φ′t ; t = 0, 1, · · · } as follow. Let Φ′t = Φt
before the optimal solution is found for the first time. Once an optimal solution is found at the t-th iteration,
then assign Φ′s = Φt for any iteration s after t. This implies the optimal set is always absorbing in the new
chain {Φ′t ; t = 0, 1, · · · }. We don’t care about the behaviour of randomised search heuristics after the first
time to hit a optimal solution. To simplify notation, we still denote the new chain by {Φt ; t = 0, 1, · · · }.
As a result, the population sequence {Φt ; t = 0, 1, · · · } is modelled by a homogeneous Markov chain where
the optimal solution set is always absorbing.

3.2

Case Studies

In this subsection, we show how randomised heuristics can be easily modelled by Markov chains. For the
sake of illustration, we consider a simple maximisation problem
max f (x),

x ∈ {0, 1, · · · , 100}

(3)

Two randomised search heuristics are applied to the above problem. The first algorithm adopts random
walk with elitist selection, denoted by RSH-I (see Algorithm 2).
Algorithm 2 RSH-I
Random Walk: provided that Φt = x, then it walks to x − 1 (if x − 1 is not less than 0) with probability
0.01, or walks to x + 1 (if x + 1 is not more than 100) with probability 0.01. Denote the new position by
Φt+1/2 .
Elitist Selection: if f (Φt+1/2 ) > f (Φt ), then let Φt+1 ← Φt+1/2 ; otherwise Φt+1 ← Φt .
Th sequence {Φt ; t = 0, 1, · · · } is a Markov chain and its transition probabilities are given as follows: for
any x, y ∈ {0, 1, · · · , 100} which are not in the optimal solution set,


 0.01, if y = x − 1 and f (y) > f (x);

0.01, if y = x + 1 and f (y) > f (x);
P (x, y) =
p,
if y = x;



0,
otherwise.

where p = 1 − P (x, x − 1) − P (x, x + 1).

3

Algorithm 3 RSH-II
Random Walk: provided that Φt = x, then it walks to x − 1 (if x − 1 is not less than 0) with probability
0.01, or walks to x + 1 (if x + 1 is not more than 100) with probability 0.01. Denote the new position by
Φt+1/2 .
Non-elitist Selection: if f (Φt+1/2 ) > f (Φt ), then definitely let Φt+1 ← Φt+1/2 ; otherwise let Φt+1 ←
Φt+1/2 with probability 0.5.

The second algorithm adopts random walk with non-elitist selection, denoted by RSH-II (see Algorithm
3). Different from RSI-I, RSH-II allows a worse child to be accepted.
Th sequence {Φt ; t = 0, 1, · · · } is a Markov chain and its transition probabilities are given as follows: for
any x, y ∈ {0, 1, · · · , 100} which are not in the optimal solution set,

0.01, if y = x − 1 and f (y) > f (x);




0.005, if y = x − 1 and f (y) ≤ f (x);



0.01, if y = x + 1 and f (y) > f (x);
P (x, y) =
0.005, if y = x + 1 and f (y) ≤ f (x);




p,
if y = x;



0,
otherwise;
where p = 1 − P (x, x − 1) − P (x, x + 1).

4

Convergence in Distribution

4.1

Theory Study: Convergence Condition

In this subsection we define the convergence in distribution of randomised search heuristics and establish the
sufficient and necessary conditions of convergence in distribution.
Let Sopt denote the set of all populations which includes at least one optimal solution, and Snon the set
of all populations which doesn’t include any optimal solution.
Definition 1. A randomised search heuristic is called convergence in distribution if starting from any nonoptimal population, the probability of Φt in the optimal set goes towards 1 as t to the infinitely large. That
is
lim P (Φt ∈ Sopt ) = 1.
t→+∞

We consider the probability distribution of Φt in the non-optimal set. Let qt (X) denote the probability
of Φt at a non-optimal state X,
qt (X) := P (Φt = X).
Let (X1 , X2 , · · · ) represent all populations in the non-optimal set. Then the vector1 qt
qt := (qt (X1 ), qt (X2 ), · · · )T
denote the probability distribution of Φt over all non-optimal populations.
Notice that the vector 1-norm2 k qt k1 equals to the population Φt in the non-optimal set, i.e.,
k qt k1 = P (Φt ∈ Snon ).
Thus the convergence of a randomised search heuristic is rewritten in the norm form:
lim k qt k1 = 0.

t→+∞

v represents a column vector and vT the row column with
Pthe transpose operation.
Pn
an n × n matrix A = [ai,j ], its 1-norm is k A k1 = max1≤j≤n n
i=1 | aij |. Given a vector v, its 1-norm is
i=1 | vi |.

1 Notation
2 Given

4

The 1-norm k qt k1 plays the role of the distance between the population Φt and the optimal set. This
can be viewed as a special case of the analysis in [15], where a general vector norm has been used as the
distance.
In the following we draw the convergence condition of randomised search heuristics based on the absorbing
Markov chain theory.
Let P denote the transition matrix of the Markov chain associated with a randomised search heuristic,
whose entries are P (X, Y ), where X, Y ∈ S. Since a state in the optimal set is always absorbing, so the
transition matrix P can be written in the canonical form below,


I O
P=
,
(4)
R Q
where I is a unit matrix and O a zero matrix. Q is a matrix to denote probability transitions among nonoptimal populations. R is a matrix to represent probability transitions from non-optimal populations to the
optimal set.
According to absorbing Markov chain theory [13, Chapter 11],
X
P (Φt+1 = Y ) =
P (Φt = X)P (X, Y ),
X∈Snon

then the iteration Φt → Φt+1 is represented in an equivalent matrix iteration,
qTt+1 = qTt Q = qT0 Qt .

(5)

From the sufficient and necessary condition of convergence of iterative methods ([35, Theorem 1.10]), it
is straightforward to obtain the convergence condition for randomised search heuristics:
Lemma 1. A randomised search heuristic is convergent if and only if the spectral radius3 ρ(Q) < 1.
The above sufficient and necessary condition is less useful in practice since it is too difficult to calculate
the spectral radius of the transition matrix. Therefore we turn to find an equivalent condition which is much
easier to verify. The following lemma gives such a condition.
Lemma 2. A randomised search heuristic is convergent if and only if starting from any non-optimal population, it is possible to visit the optimal set after finite iterations. That is, there exists an integer k > 0 and
for any non-optimal state X and t ≥ 0,
P (Φt+k ∈ Sopt | Φt = X) > 0.

(6)

Proof. Since the chain {Φt ; t = 0, 1, · · · } is homogeneous, thus it is enough to prove the case of t = 0.
(i) The proof that the condition is sufficient.
Assume (6) holds. From
P (Φk ∈ Snon | Φ0 = X) + P (Φk ∈ Sopt | Φ0 = X) = 1,
then it gives k qk k1 < 1 for any k q0 k1 such that k q0 k1 = 1.
Then from the matrix iteration qk = (QT )k q0 , and the matrix 1-norm definition
k (QT )k k1 = max k (QT )k q0 k1 ,
kq0 k1 =1

it follows
k (QT )k k1 < 1.
Since the spectral radius of a matrix is never bigger than its consistent norm [24, Example 7.1.4], the
above inequality yields

ρ (QT )k ≤k (QT )k k1 < 1,
3 The spectral radius of a square matrix A, denoted by ρ(A), is the supremum among the absolute values of all eigenvalues
of A.

5

so that ρ(Q) = ρ(QT ) < 1. According to Lemma 1, this implies that the algorithm is convergent.
(ii) The proof that the condition is necessary.
Suppose the algorithm is convergent, then according to Lemma 1, ρ(Q) = ρ(QT ) < 1. From Gelfand’s
spectral radius formula4 , there exists an integer k > 0 such that
1/k

k (QT )k k1 < 1.

(7)

and then k (QT )k k1 < 1.
From qk = (QT )k q0 , it gives k qk k1 < 1.
From
P (Φk ∈ Snon | Φ0 = X) + P (Φk ∈ Sopt | Φ0 = X) = 1,
and
k qk k1 = P (Φk ∈ Snon | Φ0 = X),
It follows then that
P (Φk ∈ Sopt | Φ0 = X) > 0,
and this proves (6).
From the above lemma, we can easily draw the following sufficient and necessary condition, which is much
easier to verify in practice.
Theorem 1. A randomised search heuristic is convergent if and only if starting from any non-optimal state,
it is possible to reach a better state after finite iterations. That is, there exists an integer k > 0 and for any
non-optimal state X and t ≥ 0,
P (ft+k > ft | Φt = X) > 0.

4.2

(8)

Case Study

In this subsection, we show how the convergence conditions are applied to determine the convergence of a
randomised search heuristic.
Example 1. Consider RSH-I and RSH-II for solving the following maximisation problem,
max f (x),

x ∈ {0, 1, · · · , 100}.

(9)

RSH-I does not converge if there exists a state x such that f (x) > f (x − 1) and f (x) > f (x + 1). This
means f (x) is a multi-modal function. At the local optimum x, the algorithm cannot make any move.
RSH-II always converges. It is easy to verify that from any state X, the algorithm can make an improvement at most 100 iterations with a positive probability, thus according to Theorem 1, the algorithm is not
convergent.
Example 2. Consider RSH-I for solving the following maximisation problem
max x2 ,

x ∈ {0, 1, · · · , 100}.

(10)

It is a unimodal function with the optimum at 100.
The corresponding transition matrix Q is

0.99
0
···
0.01 0.99 · · ·


..
..
Q =  ...
.
.

 0
0
···
0
0
···

4 Gelfand’s

0
0
..
.

0
0
..
.






.

0.99
0 
0.01 0.99

(11)

spectral radius formula says that for any induced matrix norm: limt→+∞ k Ak k1/k = ρ(A) [24, Example 7.10.1].

6

We can prove the convergence of RSH-I using Lemma 1. Since ρ(Q) = 0.99 < 1, so RSH-I for maximising
f1 (x) is convergent.
We also can prove the convergence of RSH-I using Theorem 1 without calculating the spectral radius.
It is easy to see that the probability of obtaining a better child is 0.01, thus according to Theorem 1, the
algorithm is convergent.
Comparing with the two approaches, we see that Theorem 1 provides a faster approach to determine the
convergence.

4.3

Computational Study

In this subsection, we propose a computational approach to illustrating the convergence in distribution of
randomised search heuristics.
In the computation study, we run an algorithm for k times. Let n(Φt ∈ Sopt ) denote the number of Φt
(where t = 0, 1, · · · ) appearing in the optimal set for these k runs.
According to the law of large numbers, the probability P (Φt ∈ Sopt ) will be approximated by the relative
frequency as follows:
n(Φt ∈ Sopt )
, when k → +∞.
k
The above frequency is used as the probability P (Φt ∈ Sopt ) in the the computational study.
Example 3. Consider RSH-II for solving the following problem
max (x − 49)2 ,

x ∈ {0, 1, · · · , 100}.

(12)

It is a two-modal function, with one local optimum at 0 and one global optimum at 100.
We apply RSH-I to the problem. Run the algorithm for 100,000 times. The initial population is Φ0 = 20.
Figure 1 shows the probability of Φt in the optimal set is 0 when t ≤ 100, 000. In other words, no
convergence happens yet in 100, 000 iterations.
1

P (Φt ∈ Sopt )

0.5
0
-0.5
-1

0

20000 40000 60000 80000 100000
t

Figure 1: The probability P (Φt ∈ Sopt ) when RSH-I maximises (x − 49)2 .
Our approach is different from that of visualising the fitness value over t. The latter approach may be the
most popular used for illustrating the convergence of randomised search heuristics (for example, see Figures
1 to 3 in [1]). Figure 2 shows the fitness value ft ‘converges’ after about 5,000 iterations and thereafter no
change.
In Figure 2, RSH-I seems convergent, but this is only a kind of premature convergence to the local
optimum at 0, rather than the global optimum at 100. In fact, the solution becomes farther away from the
global optimum using the Euclidean distance: initially the distance between the solution and the optimum
is 100 − 20 = 80; then after 6000 iterations, the distance increases to 100 − 0 = 100.
The approach of using the convergence in distribution provides a more accurate description of converegnce
than that of visualising the fitness value does.

7

2600
2400
2200
2000
1800
1600
1400
1200
1000
800

ft

0

20000 40000 60000 80000 100000
t

Figure 2: The fitness value ft when RSH-I maximises (x − 49)2 . ft is the mean fitness value averaged over
100,000 runs.

5

Convergence Rate

5.1

Theoretical Study: Average Convergence Rate

In this subsection, we define the average convergence rate of randomised search heuristics and then present
lower and upper bounds on the average convergence rate. The convergence rate is how fast a randomised
search heuristic converges to the optimal set per iteration. It is an important measure of the performance of
randomised search heuristics, but less studied.
Since randomised search heuristics belong to iterative methods, we adopt the average convergence rate,
commonly used in iterative methods [35, Definition 3.1].
Definition 2. Assume the probability of the initial population Φ0 in the non-optimal set is larger than 0.
The average rate of convergence for t iterations is given by the following logarithmic reduction:
1 k qt k1
.
− ln
t k q0 k1

(13)

Since P (Φt ∈ Snon ) =k qt k1 , then the average convergence rate for t iterations can be rewritten as
1 P (Φt ∈ Snon )
− ln
.
t P (Φ0 ∈ Snon )

(14)

In the above definition, we don’t consider the case of the initial population in the optimal set with
probability 1. In this case the algorithm already converges and no need to discuss the convergence rate.
Notice that
k qt k1
k qt k1 k qt− k1
k q1 k1
=
···
,
k q0 k1
k qt−1 k1 k qt−2 k1
k q0 k1
and the average convergence rate is equal to the logarithmic mean

=
=

kqt k1
− 1t ln kq

 0 k1
kq1 k1
t k1
+ · · · + ln kq
− 1t ln kqkqt−1
k
1
0 k1


P (Φ1 ∈Snon )
t ∈Snon )
− 1t ln PP(Φ(Φt−1
+
·
·
·
+
ln
∈Snon )
P (Φ0 ∈Snon ) .

The last formula shows that the average convergence rate is the average reduction factor of the probability
of Φt in the non-optimal set per iteration in terms of the logarithmic mean.
In the following we estimate the lower bound and upper bound of the average convergence rate. The
following theorem gives a lower bound on the average convergence rate.

8

Theorem 2. If a randomised search heuristic is convergent, then the averaged convergence rate is lowerbounded by
1
1 k qt k1
≥ − ln k (QT )t k1 ,
ln
t k q0 k1
t
1 k qt k1
≥ − ln ρ(Q).
− lim
ln
t→+∞ t
k q0 k1
−

(15)
(16)

Proof. From the matrix iteration qt = (QT )t q0 , we get
1 k qt k1
1 k (QT )t q0 k1
1
= ln
≤ ln k (QT )t k1 .
ln
t k q0 k1
t
k q0 k1
t
Then the first conclusion is proven.
Let t → +∞ and apply Gelfand’s spectral radius formula, then
ln k (QT )t k1
= ln ρ(QT ) = ln ρ(Q).
t→+∞
t
lim

Thus we prove that
− lim

t→+∞

1 k qt k1
≥ − ln ρ(Q).
ln
t k q0 k1

We proves the second conclusion.
The following theorems gives an upper bound on the average convergence rate.
Theorem 3. If a randomised search heuristic is convergent, then the average convergence rate is upperbounded by
−1
1
1 k qt k1
≤ − ln k ((QT )−1 )t k1
ln
,
t k q0 k1
t
1 k qt k1
≤ ln ρ(Q−1 )
ln
− lim
t→+∞ t
k q0 k1

−

(17)
(18)

Proof. From qt = (QT )t q0 , we get
q0 = ((QT )−1 )t qt ,
Hence
k q0 k1 ≤k ((QT )−1 )t k1 k qt k1 ,
k q0 k1
≤k ((QT )−1 )t k1 ,
k qt k1
then
−1
1 k qt k1
1
− ln
≤ − ln k ((QT )−1 )t k1
.
t k q0 k1
t

Then the first conclusion is proven.
According to Gelfand’s spectral radius formula and the fact ρ((QT )−1 ) = ρ(Q−1 ), we get
1/t
lim k ((QT )−1 )t k1
= ρ((QT )−1 ) = ρ(Q−1 ).
t→+∞

Then

− lim

t→+∞

1 k qt k1
≤ ln ρ(Q−1 ).
ln
t k qt k1

which is the second conclusion.
From the theoretical viewpoint, the above two theorems show lower and upper bounds on the average
convergence rate. But in practice it is hard to apply the theoretical results since both spectral radii ρ(Q)
and ρ(Q−1 ) are too difficult to calculate in most cases.
9

5.2

Computational Study

In this subsection, we propose a computational approach to illustrating the average convergence rate of
randomised heuristics. Unlike the theoretical study, the calculation of the average convergence rate is rather
simple in the computational study.
We run a randomised search heuristic for k times. Let n(Φt ∈ Snon ) denote the number of Φt (where
t = 0, 1, · · · ) appearing in the non-optimal set for these k runs. Then in practice, we will take
1 n(Φt ∈ Snon )
− ln
t
k

(19)

as the average convergence rate.
Example 4. Consider RSH-I and RSH-II for solving the maximising problem
max x2 , x ∈ {0, · · · , 100}.
We run each algorithm for 100,000 times. The initial population is Φ0 = 20. If n(Φt ∈ Snon ) ≤ 10−5
happens, we don’t calculate the average convergence rate. It is due to the following reason: the event
n(Φt ∈ Snon ) ≤ 10−5 is a small probability event. Computer simulation in 100,000 runs is not enough from
the statistical viewpoint.
Figure 3 shows the average convergence rate of RSH-I is much higher than that of RSA-II. Initially the
average convergence rate of both algorithms is 0. Then the average convergence rate of RSH-I increases from
0 to about 0.0009 quickly, but the average rate of RSH-II increases from 0 to about 0.0004 slowly.
0.001
0.0009
0.0008
0.0007
0.0006
0.0005
0.0004
0.0003
0.0002
0.0001
0

RSH-I
RSH-II

0

10000

20000
t

30000

40000

Figure 3: The average convergence rate of RSH-I and RSH-II for maximising x2 .
The average convergence rate is different from the logarithmic progress rate, ln ft , used in some references
(for example, see Figures 8 and 9 in [33]). Such a logarithmic rate may provide an intuitive description of
the fitness change, but does not give a quantitative measure of the convergence rate itself. Let’s demonstrate
this by the following example.
Example 5. Consider RSH-II for solving the following two problems,
max x2 ,

max = 10x2 ,

x ∈ {0, 1, · · · , 100}.

We run the algorithm for 100,000 times. The initial population is Φ0 = 20. Figure 4 depicts that the
logarithmic progress rate of RSH-II on the function 10x2 is lager than that on the function x2 . It is caused
by the coefficient difference between the two fitnesses. The logarithmic progress rate is not a quantitative
measure of the convergence rate.

6
6.1

Expected Hitting Times
Theoretical Study 1:: Fundamental Matrix

In this subsection, we define the expected hitting time and fundamental matrix. The first hitting time is
the number of iterations to find an optimal solution for the first time, which is an important measure of the
10

12
11
10
9
8
7
6
5

x2
10x2

0

4000

8000

12000 16000 20000
t

Figure 4: The logarithmic progress rate ln ft of RSH-I for maximising x2 and 10x2 .
performance of randomised search heuristics. Its formal definition is given as below.
Definition 3. Suppose the initial population Φ0 = X. The mean number of iterations when a randomised
search heuristic encounters an optimal solution for the first time is called the expected hitting time, denoted
by h(X).
When we talk about the expected hitting time, we always assume that randomised search heuristics are
convergent. Otherwise the expected hitting time is infinite, and that is out of our interest.
Let (X1 , X2 , · · · ) represent all populations in the non-optimal set. Then the vector
h = (h(X1 ), h(X2 ), · · · )T ,
represents expected hitting times corresponding to all non-optimal populations.
In absorbing Markov chains, the fundamental matrix plays a crucial role which is defined as follows.
Definition 4. [13, Definition 11.3] For an absorbing Markov chain {Φt ; t = 0, 1, · · · } with the transition
matrix Q, the matrix N = (I − Q)−1 is called the fundamental matrix.
Now we explain the meaning of the entry N (X, Y ) of the fundamental matrix. Since the chain is convergent
(i.e., ρ(Q) < 1), then
N = (I − Q)−1 =

+∞
X

Qt .

t=0

Rewriting the above equality in the entry form, we get
N (X, Y ) =

+∞
X

P (Φt = Y | Φ0 = X).

t=0

Therefore N (X, Y ) is the mean number of the chain visiting Y when starting at X [13, Definition 11.3].
The hitting time vector h can be calculated by the fundamental matrix. According to Theorem 11.5 in
[13], for an absorbing Markov chain {Φt ; t = 0, 1, · · · }, its expected hitting times equal to
h = N1,

(20)

where 1 is a column vector all of whose entries are 1.
However, it is difficult to apply the above result to the analysis of expected hitting time since it is
impossible to calculate the fundamental matrix in most cases.

11

6.2

Theoretical Study 2: Average Drift Analysis

In this subsection, we present average drift analysis for bounding the expected hitting time, which is seldom
investigated before. The first work to use less point-wise drift was [20]. Recently average drift analysis was
applied to the runtime analysis of an EA for unimodal functions.5
It is too difficult to calculate the expected hitting time through the fundamental matrix. Instead it is more
realistic to obtain their lower and upper bounds on expected hitting time. Drift analysis was introduced in
bounding the expected hitting time of randomised search heuristics [14, 16]. In drift analysis, d(X) is called
a drift function if d(X) ≥ 0 for any non-optimal state X and d(X) = 0 for any optimal state X. Given a
drift function d(X), drift represents the progress rate of moving towards the optima per iteration.
Definition 5. Drift at point X is defined by
∆(X) := d(X) −

X

d(Y )P (X, Y ).

Y ∈Snon

Let (X1 , X2 , · · · ) represent all populations in the non-optimal set and the vector
d = (d(X1 ), d(X2 ), · · · )T
represents the drift function values corresponding to each non-optimal state. The vector
∆ = (∆(X1 ), ∆(X2 ), · · · )T
represents the drift value corresponding to each non-optimal state.
The drift ∆(X) is determined by a single state X. So it is called point-wise drift. Now we introduce
average drift which is the average of ∆(X) over the probability distribution of Φt = X.
Definition 6. Average drift at the t-iteration is
¯ t :=
∆

X

∆(X)

X∈Snon

P (Φt = X)
.
P (Φt ∈ Snon )

(21)

Let h(Φ0 ) denote the expected hitting time when the initial population is Φ0 , that is,
X
h(Φ0 ) =
h(X)P (Φ0 = X).
X∈Snon

Let d(Φ0 ) denote the expected drift function when the initial population is Φ0 , that is,
X
d(X)P (Φ0 = X).
d(Φ0 ) =
X∈Snon

The following average drift theorem is for upper-bounding the expected hitting time.
¯ t ≥ 1,
Theorem 4. Suppose a randomised search heuristic is convergent. If for any t ≥ 0, the average drift ∆
then the expected hitting time h(Φ0 ) is upper-bounded by d(Φ0 ).
Proof. Recall that the 1-norm equals to
k qt k1 = P (Φt ∈ Snon ),
then the average drift can be rewritten in an equivalent vector form:
¯t =
∆
5 Jun

qTt
(I − Q)d.
k qTt k1

He, Tianshi Chen, Xin Yao: Average Drift Analysis and its Application. CoRR abs/1308.3080 (2013)

12

(22)

¯ t ≥ 1 can be rewritten in an equivalent vector form,
The condition that ∆
qTt
(I − Q)d ≥ 1.
k qTt k1
It follows
qTt (I − Q)d ≥k qTt k1 ,
then from k qTt k1 = qTt 1, it follows
qTt (I − Q)d ≥ qTt 1.

(23)

From the matrix iteration qTt = qT0 Qt , it follows
qT0 Qt (I − Q)d ≥ qT0 Qt 1,
Equivalently
Now summing t from 0 to k, we get
k
X


qT0 Qt − Qt+1 d ≥ qT0 Qt 1,
qT0

t

Q −Q

t+1

t=0

and simplifying both sides, it follows
qT0 d

−

qT0 Qk+1 d



≥

d≥

k
X

qT0 Qt 1,

t=0

qT0

k
X

Q

t

t=0

!

1.

(24)

Due to ρ(Q) < 1, the following two limits exist,
lim Qk+1 = O,

k→+∞

lim

k→+∞

k
X

Qt = (I − Q)−1 .

t=0

Thus when k → +∞, (24) becomes
qT0 d ≥ qT0 (I − Q)−1 1.
Recalling that h = N1 = (I − Q)−1 1, we have
qT0 d ≥ qT0 h.
which proves the conclusion.
Similarly we can establish an average drift theorem for lower-bounding the expected hitting time. Its
proof is the almost the same as that for the above theorem. We omit the proof of theorem.
¯ t ≤ 1, then the expected hitting time h(Φ0 ) is
Theorem 5. If for any iteration t ≥ 0, the average drift ∆
lower-bounded by d(Φ0 ).
Previous point-wise drift theorems [18, Theorems 2 and 3] are direct corollaries of current average drift
theorems.
Corollary 1. If for any non-optimal population X, its drift ∆(X) ≥ 1, then the expected hitting time
h(X) ≤ d(X).
Corollary 2. If for any non-optimal population X, its drift ∆(X) ≤ 1, then the expected hitting time
h(X) ≥ d(X).
In point-wise drift theorems, its requirement is that the drift is not less than 1 (or not more than 1) for
all non-optimal states. In average drift theorems, the condition is replaced by that average drift is not less
than 1 (or not more than 1). Hence average drift analysis is more powerful than point-wise drift analysis.
13

6.3

Theoretical Study 3: Backward Drift Analysis

In this subsection, we present novel backward drift analysis, which was never discussed before. We call
the drift defined in the previous subsection forward drift in order to distinguish it from the backward drift
introduced in the current subsection. Forward and backward drift analysis can be regarded as a dual pair.
Starting from the fundamental matrix, we already know that the vector h = N1 represents all expected
hitting times in the non-optimal set. Similarly the vector
sT := 1T N
gives another type of important times for randomised search heuristics. Now we explain the intuitive meaning
of the vector s. Notice that the entry
X
s(Y ) =
N (X, Y ).
X∈Snon

and recall that N (X, Y ) is the expected number of that the Markov chain visits Y when starting at X, then
s(Y ) is the sum of the expected number of visiting state Y when starting from all non-optimal states. We
call s(Y ) the expected staying time in non-optimal state Y .
The expected hitting time and expected staying time have the following relationship.
Theorem 6. Let h(X) be the expected hitting time from non-optimal population X and s(Y ) the expected
staying time in non-optimal population Y . Then
X
X
s(Y )
h(X) =
Y ∈Snon

X∈Snon

Proof. From the equalities
X

h(X) = 1T N1,

X∈Snon

X

s(X) = 1T N1,

X∈Snon

we draw the conclusion.
The above theorem implies that the expected hitting time equals to the expected staying time when the
initial population is chosen at uniformly random.
Next we establish backward drift analysis for bounding the expected staying time. Like forward drift
analysis, a drift function is used in backward drift analysis. d(Y ) is called a drift function if d(Y ) ≥ 0 for
any non-optimal state Y and d(Y ) = 0 for any optimal state Y .
Definition 7. Let d(Y ) be a drift function. For a non-optimal population Y , the backward forward drift is
X
∇(Y ) := d(Y ) −
d(X)P (X, Y ).
X∈Snon

Backward drift is to measure the move from Y to X (backward). This is different from forward drift
X
d(Y )P (X, Y ).
∆(X) := d(X) −
Y ∈Snon

which is to measure the move from X to Y (forward).
The following backward drift theorem is used for upper-bounding the staying time.
Theorem 7. If for any non-optimal population Y , its backward drift ∇(Y ) ≥ 1, then the staying time
s(Y ) ≤ d(Y ).

14

Proof. Let (Y1 , Y2 , · · · ) represents all non-optimal populations and the vector
∇T = (∇(Y1 ), ∇(Y2 ), · · · )
represents their drift function values respectively. Then from the definition, ∇T can be rewritten in the
vector form as follows:
∇T := dT (I − Q).

(25)

The condition that ∇(X) ≥ 1 is rewritten in a vector form dT (I−Q) ≥ 1T , and it follows dT (I−Q)−1T ≥
T

0 .
Since the fundamental matrix N is non-negative, then

dT (I − Q) − 1T N ≥ 0T ,

Since N = (I − Q)−1 , so it follows dT − 1T N ≥ 0T , and sT ≤ dT which proves the conclusion.
Similarly we can establish a backward drift theorem for lower-bounding the staying time.
Theorem 8. If for any non-optimal population X, its backward drift ∇(X) ≤ 1, then the staying time
satisfies s(X) ≥ d(X).
Backward drift analysis provides an alternative way of bounding the expected hitting time when the initial
population is chosen at uniformly random. It is possible to establish average backward drift analysis similar
to average forward drift analysis. We will not discuss it in the detail.

6.4

Case Studies

In this subsection, we demonstrate how average drift analysis and backward drift analysis are applied to the
estimation of the expected hitting time. There are three steps when applying drift analysis. First, choose
an appropriate drift function; then estimate drift; finally, obtain a bound on the expected hitting time or
staying time.
The first example is to show that average drift theorems work well but point-wise drift theorems are not
applicable.
Example 6. Consider RSH-I for the maximizing problem
max x2 ,

x ∈ {0, 1, · · · , 100}.

Choose the drift function as follows
d(x) =



100×101
(100 − x),
99
d(x + 1),

if 1 ≤ x ≤ 99,
if x = 0.

Calculate the drift ∆(x). For 1 ≤ x ≤ 99,
∆(x) =d(x) − 0.99d(x) − 0.01d(x + 1) =

101
,
99

∆(0) =d(0) − 0.99d(0) − 0.01d(1) = 0.
Notice that ∆(0) = 0, thus point-wise drift theorems cannot be applied here.
However, average drift theorems work well. Assume that the initial population Φ0 is chosen at uniformly
random, that is, P (Φ0 = x) = 1/101. The average forward drift is
!
99
X
1
1
¯0 =
∆(0) +
∆(x) = 1.
∆
101
101 x=0

15

Since RSH-I adopts elitist selection, so that if initial population Φ0 = x is not at state 0, then for any
t ≥ 0, its offspring Φt never returns to state 0. Thus the average forward drift is not less than
!
99
X
1
1
¯t =
∆(0) +
∆(x) = 1.
∆
101
101 x=0
According to Theorem 4, the expected hitting time is not more than
99

1 X
d(x) = 5000.
101 x=0
The second example is to show that we can derive the same conclusion using backward drift analysis.
Example 7. Still consider RSH-I for the maximizing problem,
max x2 ,

x ∈ {0, 1, · · · , 100}.

Choose the drift function as follow
d(x) = 100(x + 1),

0 ≤ x ≤ 99,

Calculate the backward drift ∇(x). For 1 ≤ x ≤ 99,
∇(x) =d(x) − 0.99d(x) − 0.01d(x − 1) = 1,
∇(0) =d(0) − 0.99d(0) = 1.
According to Theorems 7 and 8, the staying time
s(x) = d(x) = 100(x + 1).
Furthermore if the initial population is chosen at uniformly random, then according to Theorem 6, the
expected hitting time equals to
99
99
1 X
1 X
s(x) =
100(x + 1) = 5000.
101 x=0
101 x=0

The result is the same as that in the first example by applying average drift analysis.

6.5

Computational Study

In this subsection, we illustrate a computational approach to the calculation of the expected hitting time.
In the computational study, we runs a randomised search heuristic for k times. Let τi (X) be the first
hitting time for the i-th run during these k runs.
From the law of large numbers, the expected hitting time h(x) is approximated by the average
k
1X
τi (X), when k → +∞.
k i=1

(26)

The above average value is taken as the expected hitting time h(X) in the computational study.
Sometimes it is easy to study the expected hitting time through the computational study.
Example 8. Consider RSH-I and RSH-II for solving the maximization problem
max x2 ,

x ∈ {0, 1, · · · , 100}.

We run each algorithm for 100,000 times. The initial population is Φ0 = 20. Figure 5 shows that the
expected hitting time of RSH-II is about 16000, which is twice as long as that of RSH-I, about 8000.
16

17000
16000
15000
14000
13000
12000
11000
10000
9000
8000
7000

RSH-II
RSH-I

0

50000
k

100000

Figure 5: Expected hitting times of RSH-I and RSH-II for maximising f2 (x). k is the number of runs.
However, sometimes the calculation of the expected hitting time needs an extremely long computation
time.
Example 9. Consider RSH-II for solving the maximization problem
max (x − 49)2 ,

x ∈ {0, 1, · · · , 100}.

We run each algorithm for 100,000 times. We found that RSH-II always got stuck at the local optimum
0. Even for the simple problem, it seems not easy to make a computational study of the expected hitting time.

7

Conclusions and Future Work

A unified Markov chain approach has been proposed for studying the convergence, convergence rate and
expected hitting time of randomised search heuristics in this paper. The core of the approach is to model
randomised search heuristics by absorbing Markov chains and then to study the chain based on matrix
iteration analysis. A novelty in the analysis is that the vector 1-norm is used to represent the probability of
a population in the non-optimal solution set. It plays the role of the distance between a population and the
optimal solution set. All theoretical results are proven in a unified manner.
The results of the paper are summarised as follow: First, Theoreom 1 establishes a sufficient and necessary
condition of convergence in distribution. The theorem states a randomised search heuristic is convergent if
and only if the algorithm can make an improvement in finite iterations suppose the current solution is not
optimal.
Then the average convergence rate is introduced, which refers to the average reduction factor of the
probability of a population in the non-optimal set per iteration in terms of the logarithmic mean. Theorems 2
and 3 provide lower and upper bounds on the average convergence rate.
Finally, two new types of drift analysis, average drift analysis and backward drift analysis, are proposed
for analysing the expected hitting time. Theorems 4 and 5 state that the expected hitting time can be
bounded by a drift function and related average drift. Theorem 6 reveals that the expected hitting time and
expected staying time are equal if the initial population is chosen at uniformly random. Theorems 7 and 8
state that the expected staying time can be bounded by a drift function and related backward drift.
Besides the theoretical study, computation approaches are also presented to study the convergence, average
convergence rate and expected hitting time.
Comparing computational and theoretical studies, we see none of them are perfect. The computational
study may provide an intuitive description of randomised search heuristics’ behaviour, but it belongs to a
posterior and case study. On the other hand, the theoretical study belongs to a prior and general study.
Theoretical results may provide some understanding of randomised search heuristics’ ability. But it is still
too hard to obtain an exact value of the expected hitting time or average convergence rate.
There is a gap between the theoretical and computational studies. Our future work is to improve theoretical tools, to apply them to different types of randomised search heuristics and to study their convergence,
convergence rate and expected hitting time.
17

References
[1] T. Bäck and H. Schwefel. An overview of evolutionary algorithms for parameter optimization. Evolutionary computation, 1(1):1–23, 1993.
[2] E. Cantú-Paz. Markov chain models of parallel genetic algorithms. Evolutionary Computation, IEEE
Transactions on, 4(3):216–226, 2000.
[3] U. K. Chakraborty, K. Deb, and M. Chakraborty. Analysis of selection algorithms: A markov chain
approach. Evolutionary Computation, 4(2):133–167, 1996.
[4] T. Davis and J. Principe. A markov chain framework for the simple genetic algorithm. Evolutionary
computation, 1(3):269–288, 1993.
[5] K. A. De Jong, W. M. Spears, and D. F. Gordon. Using markov chains to analyze gafos. In Foundations
of Genetic Algorithms 3, 1995.
[6] B. Doerr and L. A. Goldberg. Adaptive drift analysis. In Parallel Problem Solving from Nature, PPSN
XI, pages 32–41. Springer, 2010.
[7] B. Doerr, D. Johannsen, and C. Winzen. Multiplicative drift analysis. In Proceedings of the 12th Annual
Conference on Genetic and Evolutionary Computation, pages 1449–1456. ACM, 2010.
[8] S. Droste, T. Jansen, and I. Wegener. On the analysis of the (1+ 1) evolutionary algorithm. Theoretical
Computer Science, 276(1-2):51–81, 2002.
[9] X. Du and L. Ding. About the convergence rates of a class of gene expression programming. Science
China Information Sciences, 53(4):715–728, 2010.
[10] A. E. Eiben, E. H. Aarts, and K. M. Van Hee. Global convergence of genetic algorithms: A markov
chain analysis. In Parallel problem solving from nature, pages 3–12. Springer, 1991.
[11] D. Fogel. Asymptotic convergence properties of genetic algorithms and evolutionary programming:
analysis and experiments. Cybernetics and Systems, 25(3):389–407, 1994.
[12] D. E. Goldberg and P. Segrest. Finite markov chain analysis of genetic algorithms. In Proceedings of the
Second International Conference on Genetic Algorithms on Genetic Algorithms and Their Application,
pages 1–8, Hillsdale, NJ, USA, 1987. L. Erlbaum Associates Inc. ISBN 0-8058-0158-8.
[13] C. Grinstead and J. Snell. Introduction to Probability. American Mathematical Society, 1997.
[14] J. He. A study of the foundation of evolutionary computation. Postdoc techreport, Department of
Computer Science, Harbin Institute of Technology, 1998. (in Chinese).
[15] J. He and L. Kang. On the convergence rate of genetic algorithms. Theoretical Computer Science, 229
(1-2):23–39, 1999.
[16] J. He and X. Yao. Drift analysis and average time complexity of evolutionary algorithms. Artificial
Intelligence, 127(1):57–85, 2001.
[17] J. He and X. Yao. From an individual to a population: An analysis of the first hitting time of populationbased evolutionary algorithms. IEEE Transactions on Evolutionary Computation, 6(5):495–511, 2002.
[18] J. He and X. Yao. Towards an analytic framework for analysing the computation time of evolutionary
algorithms. Artificial Intelligence, 145(1-2):59–97, 2003.
[19] J. He and X. Yu. Conditions for the convergence of evolutionary algorithms. Journal of Systems
Architecture, 47(7):601–612, 2001.
[20] J. Jägersküpper. Combining Markov chain analysis and drift analysis: The (1+ 1) evolutionary algorithm
on linear functions reloaded. Algorithmica, 59(3):409–424, 2011.
18

[21] D. Johannsen. Random Combinatorial Structures and Randomized Search Heuristics. PhD thesis, Universität des Saarlandes, 2010.
[22] J. Kemeny and J. Snell. Finite Markov Chains. van Nostrand, 1960.
[23] P. K. Lehre. Negative drift in populations. In Parallel Problem Solving from Nature, PPSN XI, pages
244–253. Springer, 2010.
[24] C. Meyer. Matrix analysis and applied linear algebra: solutions manual. SIAM, 2000.
[25] B. Mitavskiy, J. Rowe, and C. Cannings. Theoretical analysis of local search strategies to optimize network communication subject to preserving the total number of links. International Journal of Intelligent
Computing and Cybernetics, 2(2):243–284, 2009.
[26] F. Neumann, D. Sudholt, and C. Witt. Analysis of different mmas aco algorithms on unimodal functions
and plateaus. Swarm Intelligence, 3(1):35–68, 2009.
[27] A. Nix and M. Vose. Modeling genetic algorithms with markov chains. Annals of mathematics and
artificial intelligence, 5(1):79–88, 1992.
[28] P. Oliveto and C. Witt. Simplified drift analysis for proving lower bounds in evolutionary computation.
Algorithmica, 59(3):369–386, 2011.
[29] P. Oliveto, J. He, and X. Yao. Time complexity of evolutionary algorithms for combinatorial optimization: A decade of results. International Journal of Automation and Computing, 4(3):281–293, 2007.
[30] G. Rudolph. Convergence analysis of canonical genetic algorithms. IEEE Transactions on Neural
Networks, 5(1):96–101, 1994.
[31] G. Rudolph. Finite Markov chain results in evolutionary computation: a tour d’horizon. Fundamenta
Informaticae, 35(1):67–89, 1998.
[32] G. Rudolph and A. Agapie. Convergence properties of some multi-objective evolutionary algorithms.
In Proceedings of the 2000 Congress on Evolutionary Computation, volume 2, pages 1010–1016. IEEE,
2000.
[33] R. Salomon. Evolutionary algorithms and gradient search: similarities and differences. IEEE Transactions on Evolutionary Computation, 2(2):45–55, 1998.
[34] J. Suzuki. A Markov chain analysis on simple genetic algorithms. IEEE Transactions on Systems, Man
and Cybernetics, 25(4):655–659, 1995.
[35] R. Varga. Matrix Iterative Analysis. Springer, 2009.
[36] C. Witt. Optimizing linear functions with randomized search heuristics-the robustness of mutation. In
29th Symposium on Theoretical Aspects of Computer Science (STACS12)., pages 420–431, 2012.
[37] Y. Zhou, J. He, and Q. Nie. A comparative runtime analysis of heuristic algorithms for satisfiability
problems. Artificial Intelligence, 173(2):240–257, 2009.

19

