# Dynamics and Stability of the Sampling Distribution  of

Authors: Riccardo Poli
Journal: Journal of Artificial Evolution and Applications
 (2008)

---

## Abstract

For stochastic optimisation algorithms, knowing the probability distribution with which
an algorithm allocates new samples in the search space is very important, since this explains
how the algorithm really works and is a prerequisite to being able to match algorithms to
problems. This is the only way to beat the limitations highlighted by the no-free lunch theory.
Yet, the sampling distribution for velocity-based particle swarm optimisers has remained a
mystery for the whole of the first decade of PSO research. In this paper, a method is presented
that allows one to exactly determine all the characteristics of a PSO's sampling distribution
and explain how it changes over time during stagnation (i.e., while particles are in search for
a better personal best) for a large class of PSO's.


---

## Body

## 1. Introduction

Let us consider the basic form of PSO with inertia
weight, which is controlled by the following two equations:(1)vt+1i=wvti+ϕ1⊗(yi−xti)+ϕ2⊗(y^−xti),xt+1i=xti+vt+1i,  where ϕ1 and ϕ2 represent
vectors of random numbers uniformly distributed in [0,c1] and [0,c2], respectively, (c1 and c2 are nonnegative constants), ⊗ is
component wise multiplication, yi is the position
of the ith particle's
personal best, y^ is the position
of the particle's neighbourhood best, and w is a constant
in [0,1]. Following Kennedy's graphical examinations of the
trajectories of individual particles and their responses to variations in key
parameters [1], traditionally, this model (and the related PSO with
constriction [2], which is, provably, mathematically equivalent) has been
studied theoretically under strong simplifying assumptions such as isolated
single individuals, search stagnation (i.e., no improved solutions are found)
and absence of randomness [2–12]. While these assumptions make it possible to make
theoretical progress, they also make the results of the analyses difficult to
carry over to the real PSO (where stochasticity is present).Only very few attempts to understand the behaviour of
the PSO in thepresence of
stochasticity have been made. For example, Clerc [13] analysed the
distribution of velocities of one particle controlled by the standard PSO
update rule with inertia and stochastic forces and was able to show that a
particle's new velocity is the sum of three components: a forward force, a
backward force, and noise. Kadirkamanathan et al. [14] were able to study the
stability of particles in the presence of stochasticity by using Lyapunov
stability analysis.These are good and important first steps in the
direction of modelling the real PSO. However, they only marginally scratch the
surface in relation to one of the most important open problems in PSO research:
the characterisation of the PSO's sampling distribution and its changes over
time. Why is this question so fundamental? For two reasons. Firstly, this
scientific knowledge gives us much greater understanding of an algorithm's
behaviour. Secondly, on a more practical level, the only way to beat the
limitations implied by NFL [15] is to match algorithms to problems. Naturally,
at least in part, this can be achieved by a trial-and-error approach involving
testing different algorithms and parameter settings on a particular problem or
class of problems until a suitable match is found. A more analytic approach to
this problem, however, requires two things: (a) knowing the search distribution
of an algorithm and (b) checking whether this is well matched to the features
of a problem's fitness landscape. Unfortunately, the sampling distribution of
the PSO has been unavailable during the whole first decade (and beyond) of PSO
research.In recent work [16], we introduced a novel method,
which allowed us to exactly determine the mean and standard deviation of the
sampling distribution of the canonical PSO as well as their changes over any
number of generations. The only assumption we made was stagnation, that is, we
studied the sampling distribution produced by particles in search for a better
personal best. In this paper, we generalise the technique and show how it can
be used to determine any number of moments of the sampling distribution of
PSO's.The paper is organised as follows. In Section2, we
provide a summary of the results presented in [16]. In Section 3, we
generalised the method to moments of any order. In Section 4, we apply the
technique to compute the skewness and kurtosis of the canonical PSO's sampling
distribution. In Section 5, we further extend the generality of our technique
so as to cover a variety of PSO's, not just the one in (1). This allows us
to compare, on a theoretical basis, different forms of PSO. In Section 6, we
then use the moments of the PSO's sampling distribution to approximately
reconstruct the distribution itself, making it possible, for the first time, to
understand the search strategy implemented by the canonical PSO, while
searching for new improvements. We discuss the results and conclude in Section
7.


## 2. Dynamics of First- and Second Order Moments of the PSO's Sampling Distribution

### 2.1. Derivation of Dynamic Equations for
the Moments

During stagnation, each particle behaves
independently. Also, each dimension is treated independently. So we can analyse
each particle's behaviour in isolation. Therefore, we can drop the superscripti in (1),
and rewrite them as(2)xt+1=xt(1+w)−xt(ϕ1+ϕ2)−wxt−1+ϕ1y+ϕ2y^,  where we used
the relation vt=xt−xt−1. Note that on the r.h.s. of this equation, while w, y, and y^ are constants, xt, ϕ1, ϕ2, and xt−1 are stochastic
variables. As a consequence, xt+1 on the
l.h.s. of the equation is a
stochastic variable too (in fact, it is a stochastic function of the stochastic
variables appearing in the r.h.s.).In [16], we applied the expectation operator to both
sides of the equation obtaining(3)E[xt+1]=E[xt](1+w)−E[xt](E[ϕ1]+E[ϕ2])−wE[xt−1]+E[ϕ1]y+E[ϕ2]y^,  where we
performed the substitution E[xtϕi]=E[xt]E[ϕi] because of the
statistical independence between ϕi and xt. Assuming that c1=c2=c, that is, the ϕ1 and ϕ2 are both
uniformly distributed in [0,c], we have E[ϕi]=c/2, and so(4)E[xt+1]=E[xt](ρ−c)−wE[xt−1]+cy+y^2,  where we
renamed (1+w)=ρ.If we square both sides of (2), we obtain(5)xt+12=(xtρ−xtϕ1−xtϕ2−wxt−1+ϕ1y+ϕ2y^)2.  Expanding the
r.h.s., we obtain an equation expressing the stochastic variable xt+12 as a function
of other stochastic variables, that is, xt2, xt−1xt, ϕ12, and so forth. By applying the expectation operator
to both sides of this equation, one obtains(6)E[xt+12]=E[xt2](ρ2−4μρ+2ν+2μ2)+E[xt−1xt](−2wρ+4wμ)+E[xt−12](w2)+E[xt](2μyρ+2μy^ρ−2νy−2μ2y^−2μ2y−2νy^)+E[xt−1](−2μyw−2μy^w)+νy2+2μ2yy^+νy^2,  where we set μ=E[ϕi]=c/2 and ν=E[ϕi2]=c2/3, for brevity.If, instead, we multiply both sides of (2) by xt and apply the
expectation operator, we obtain(7)E[xt+1xt]=E[xt2](ρ−c)−wE[xtxt−1]+E[xt]cy+y^2.Equations (4), (6), and (7) form a system of coupled
second-order difference equations using which, given appropriate initial
conditions, one can compute E[xt], E[xt2] and E[xtxt−1] for any t.


### 2.2. Analysis of the Dynamics

Equations (4), (6), and (7) can also be written in
matrix notation as the extended first order system(8)z(t+1)=Mz(t)+b,  where(9)z(t)=(E[xt]E[xt−1]E[xt2]E[xt−12]E[xtxt−1])T,M=(ρ−c−w00010000A1A2A3w22w(2μ−ρ)00100cp0ρ−c0−w),b=(cp0νy2+2μ2yy^+νy^200),where A1=4p(μρ−ν−μ2), A2=−4μwp, A3=ρ2−4μρ+2ν+2μ2.It is then trivial to verify under what conditionsE[xt], E[xt2], and E[xtxt−1] will converge
to stable fixed points. We need to have that all eigenvalues of M must be within
the unit circle, that is, Λm=maxi|λi|<1. When this happens, we will say that the PSO is order-2 stable.Naturally, whenΛm<1, in principle, we could derive (either symbolically
or numerically) the fixed point for the system, which we will denote as z*. This would be simply given by(10)z*=(I−M)−1b. When the system is order-2 stable, by the simple
change of variables z˜(t)=z(t)−z*, we can represent the dynamics of the system via the
following linear homogeneous equation:(11)z˜(t+1)=Mz˜(t), which can trivially be integrated to obtain the
following explicit solution for the dynamics of the first two moments:(12)z˜(t)=Mtz˜(0).


## 2.1. Derivation of Dynamic Equations for
the Moments

During stagnation, each particle behaves
independently. Also, each dimension is treated independently. So we can analyse
each particle's behaviour in isolation. Therefore, we can drop the superscripti in (1),
and rewrite them as(2)xt+1=xt(1+w)−xt(ϕ1+ϕ2)−wxt−1+ϕ1y+ϕ2y^,  where we used
the relation vt=xt−xt−1. Note that on the r.h.s. of this equation, while w, y, and y^ are constants, xt, ϕ1, ϕ2, and xt−1 are stochastic
variables. As a consequence, xt+1 on the
l.h.s. of the equation is a
stochastic variable too (in fact, it is a stochastic function of the stochastic
variables appearing in the r.h.s.).In [16], we applied the expectation operator to both
sides of the equation obtaining(3)E[xt+1]=E[xt](1+w)−E[xt](E[ϕ1]+E[ϕ2])−wE[xt−1]+E[ϕ1]y+E[ϕ2]y^,  where we
performed the substitution E[xtϕi]=E[xt]E[ϕi] because of the
statistical independence between ϕi and xt. Assuming that c1=c2=c, that is, the ϕ1 and ϕ2 are both
uniformly distributed in [0,c], we have E[ϕi]=c/2, and so(4)E[xt+1]=E[xt](ρ−c)−wE[xt−1]+cy+y^2,  where we
renamed (1+w)=ρ.If we square both sides of (2), we obtain(5)xt+12=(xtρ−xtϕ1−xtϕ2−wxt−1+ϕ1y+ϕ2y^)2.  Expanding the
r.h.s., we obtain an equation expressing the stochastic variable xt+12 as a function
of other stochastic variables, that is, xt2, xt−1xt, ϕ12, and so forth. By applying the expectation operator
to both sides of this equation, one obtains(6)E[xt+12]=E[xt2](ρ2−4μρ+2ν+2μ2)+E[xt−1xt](−2wρ+4wμ)+E[xt−12](w2)+E[xt](2μyρ+2μy^ρ−2νy−2μ2y^−2μ2y−2νy^)+E[xt−1](−2μyw−2μy^w)+νy2+2μ2yy^+νy^2,  where we set μ=E[ϕi]=c/2 and ν=E[ϕi2]=c2/3, for brevity.If, instead, we multiply both sides of (2) by xt and apply the
expectation operator, we obtain(7)E[xt+1xt]=E[xt2](ρ−c)−wE[xtxt−1]+E[xt]cy+y^2.Equations (4), (6), and (7) form a system of coupled
second-order difference equations using which, given appropriate initial
conditions, one can compute E[xt], E[xt2] and E[xtxt−1] for any t.


## 2.2. Analysis of the Dynamics

Equations (4), (6), and (7) can also be written in
matrix notation as the extended first order system(8)z(t+1)=Mz(t)+b,  where(9)z(t)=(E[xt]E[xt−1]E[xt2]E[xt−12]E[xtxt−1])T,M=(ρ−c−w00010000A1A2A3w22w(2μ−ρ)00100cp0ρ−c0−w),b=(cp0νy2+2μ2yy^+νy^200),where A1=4p(μρ−ν−μ2), A2=−4μwp, A3=ρ2−4μρ+2ν+2μ2.It is then trivial to verify under what conditionsE[xt], E[xt2], and E[xtxt−1] will converge
to stable fixed points. We need to have that all eigenvalues of M must be within
the unit circle, that is, Λm=maxi|λi|<1. When this happens, we will say that the PSO is order-2 stable.Naturally, whenΛm<1, in principle, we could derive (either symbolically
or numerically) the fixed point for the system, which we will denote as z*. This would be simply given by(10)z*=(I−M)−1b. When the system is order-2 stable, by the simple
change of variables z˜(t)=z(t)−z*, we can represent the dynamics of the system via the
following linear homogeneous equation:(11)z˜(t+1)=Mz˜(t), which can trivially be integrated to obtain the
following explicit solution for the dynamics of the first two moments:(12)z˜(t)=Mtz˜(0).


## 3. Higher-Order Moments

In [16], we obtained recursions which describe the
dynamics of first- and second-order moments of the sampling distribution of a
standard PSO during stagnation, as summarised in the previous section. One may
then wonder whether it would be possible to follow a similar approach to study
the dynamics of higher-order moments.


### 3.1. Derivation of Dynamic Equations Forhigher-Order Moments

The fundamental question is “what quantities
would we have to deal with if we took higher powers of both sides of (2) as we
did to derive (5)?” Generally, the r.h.s. would be a sum of terms of the
form(13)a0xta1xt−1a2wa3ϕ1a4ϕ2a5ya6y^a7,  where ak are suitable
constants. Naturally, taking powers of (2) and then multiplying both sides by
some power of xt would also lead
to equations involving terms such as those in (13). That is, for any choice of b1∈ℕ and b2∈ℕ,(14)xt+1b1xtb2=∑ia0ixta1ixt−1a2iwa3iϕ1a4iϕ2a5iya6iy^a7i,  where aki are suitable
constants. If we then take expectations for both sides, we obtain(15)E[xt+1b1xtb2]=∑ia0iwa3iya6iy^a7iE[ϕ1a4i]E[ϕ2a5i]E[xta1ixt−1a2i],  where we used
the independence of ϕ1, ϕ2, xtxt−1, and, of course, their powers. If ϕ1 and ϕ2 are uniformly
distributed in the same range, namely, [0,c], it is easy to verify that(16)E[ϕjn]=cnn+1.  So(17)E[xt+1b1xtb2]=∑iωiE[xta1ixt−1a2i],  where(18)ωi=(a0iwa3ica4i+a5iya6iy^a7i(1+a4i)(1+a5i)).  It is important
to note here that because (2) is linear in xt and xt−1, all the terms on the r.h.s. of (14) and (17) respect
the relation a1i+a2i≤b1+b2. This implies that it is possible to construct recursions for
moments of arbitrary order.


### 3.2. Complexity of the Calculations

If one wanted to push the analysis developed in
Section2 up to order 3, one would need to instantiate (17) for E[xt+13], E[xt+12xt], E[xt+1xt2] and add the
resulting equations to (4), (6), and (7). If one wanted to go to order 4, an
additional set of four equations (for E[xt+14], E[xt+13xt], E[xt+12xt2], and E[xt+1xt3]) would be
needed, bringing the total to 10.More generally, in order to compute statistics of
ordern, one needs to construct and iterate(19)Q(n)=n×(n+1)2 second-order
difference equations. Since, after expansion, the r.h.s. of (2) contains 7
atomic terms of the form in (13), the r.h.s. of (17) contains 7b1 such
terms (Note that the exponent b2 does not influence
the number of terms because the recursion for E[xt+1b1xtb2] is
obtained as follows: (a) we compute xt+1b1, which is given by an expression containing 7b1 terms; (b) we
multiply each term by xtb2, which changes the exponents a1i but does not
alter the number of terms; (c) we apply the expectation operator, which again
does not modify the number of terms.) So the total number of terms one needs to
compute to construct the equations for order- n statistics is 7 (for the order
1 equations) plus (72+7) (for the order
2 equations) plus (73+72+7) (for the order
3 equations), and so forth. This gives us a total of(20)T(n)=∑i=1n(n−i+1)×7i terms. Note that T(n) grows
exponentially approximately as 1.36×7n. So although the number of equations one needs to
deal with grows quadratically, the computational effort required to instantiate
them is exponential. For instance, T(1)=7, T(2)=63, T(3)=462, T(4)=3262, and T(5)=22869. The growth in number of terms can be reduced if one
makes explicit use of ρ (i.e., by
adding the factor ρa8 in (13)). Then T(n)=O(6n). Either way, manually deriving equations for moments
of order 3 is already very laborious. The process, however, is clearly
mechanisable. This can be done using computer algebra systems, or by explicitly
representing and manipulating the ωi's in each
equation (this is what we did). As a result of mechanisation, computing the
equations for up to order 6 or 7 is feasible with an ordinary personal
computer.Some of theωi's in (17)
present the same pattern of exponents for w, c, y, and y^, so terms can be collected leading to more compact
equations (e.g., compare (5) and (6)). Also, given their size, one will
normally want to study (e.g., integrate) (17) numerically. In this case, w, c, y, and y^ are all
numerical parameters. So the ωi's become
constants and, after collecting terms, each equation contains at most Q(n) terms, which,
as we know, is quadratic in the order n. As a result, although the complexity of the
construction of the motion equations for the moments is exponential in the
order of the moments, their numerical integration is only of order O(n4).Naturally, the system ofQ(n) second-order
difference equations necessary to predict the dynamics of moments of order 1 to n can be turned
into a system of order 1 of the form in (8), via the choice(21)z(t)=(E[xt]E[xt−1]E[xt2]E[xtxt−1]E[xt−12]E[xt3]E[xt2xt−1]E[xtxt−12]E[xt−13]⋮E[xt−1n]).This effectively means adding artificial update
equations of the formE[xtk]=E[xtk] for k=1,…,n, bringing the total to Q′(n)=Q(n)+n. The transition matrix for the system is, therefore,
of size Q′(n)×Q′(n). We will denote this with Mn. For example, for n=4, which would allow one to study the mean, variance,
skewness, and kurtosis of the sampling distribution as a function of t, M4 is merely a 14×14 matrix.Interestingly,Q′(n) grows so slowly
that one can perform an eigenvalue analysis for any Mn that one is
able to compute. That is, the expensive part of the process is the construction
of Mn. Once this is done, iterating the system,
establishing its stability, or finding its fixed points is a trivial matter.In the next section, we provide results for statistics
of order 3 and 4, that is,n=4, a value of n for which
computing Mn takes only a
few seconds. However, before we do this, we need to consider the initial
conditions for the system.


### 3.3. Initial Conditions

In order to provide initial conditions for the
moments' dynamics, we need to computeE[x0k] and E[x1kx0l] for generic k>0 and l≥0.If we are at the very first iteration of the PSO
algorithm, the initial conditions for the dynamics of the moments are related
to the ranges used to initialise positions and velocities in the PSO. Under the
assumption that a particle's initial position,x0, is chosen uniformly at random in a symmetric range [−Ω,Ω], we have(22)E[x0k]={0ifk⁢is odd,Ωkk+1otherwise.  In order to
compute E[x1kx0l], we need to consider the equation(23)x1=x0+wv0−x0(ϕ1+ϕ2)+ϕ1y+ϕ2y^,  where a
particle's initial velocity, v0, is a stochastic variable uniformly distributed the
range [−Ωv,Ωv] (often Ωv=Ω). By taking
the kth power of
both sides of the equation, multiplying by x0l, and taking expectations, as we did to construct
(17), one obtains the desired expressions for E[x1kx0l]. Like for (17), these expressions contain a number of
terms that grows exponentially with n. However, this process, too, can be trivially
mechanised.If the PSO is not at its first iteration, and a new
improvement has just been found to a particle's personal best or neighbourhood
best, thenx0 is not a
stochastic variable, but a constant determined by the position where the
particle found the last improvement.


## 3.1. Derivation of Dynamic Equations Forhigher-Order Moments

The fundamental question is “what quantities
would we have to deal with if we took higher powers of both sides of (2) as we
did to derive (5)?” Generally, the r.h.s. would be a sum of terms of the
form(13)a0xta1xt−1a2wa3ϕ1a4ϕ2a5ya6y^a7,  where ak are suitable
constants. Naturally, taking powers of (2) and then multiplying both sides by
some power of xt would also lead
to equations involving terms such as those in (13). That is, for any choice of b1∈ℕ and b2∈ℕ,(14)xt+1b1xtb2=∑ia0ixta1ixt−1a2iwa3iϕ1a4iϕ2a5iya6iy^a7i,  where aki are suitable
constants. If we then take expectations for both sides, we obtain(15)E[xt+1b1xtb2]=∑ia0iwa3iya6iy^a7iE[ϕ1a4i]E[ϕ2a5i]E[xta1ixt−1a2i],  where we used
the independence of ϕ1, ϕ2, xtxt−1, and, of course, their powers. If ϕ1 and ϕ2 are uniformly
distributed in the same range, namely, [0,c], it is easy to verify that(16)E[ϕjn]=cnn+1.  So(17)E[xt+1b1xtb2]=∑iωiE[xta1ixt−1a2i],  where(18)ωi=(a0iwa3ica4i+a5iya6iy^a7i(1+a4i)(1+a5i)).  It is important
to note here that because (2) is linear in xt and xt−1, all the terms on the r.h.s. of (14) and (17) respect
the relation a1i+a2i≤b1+b2. This implies that it is possible to construct recursions for
moments of arbitrary order.


## 3.2. Complexity of the Calculations

If one wanted to push the analysis developed in
Section2 up to order 3, one would need to instantiate (17) for E[xt+13], E[xt+12xt], E[xt+1xt2] and add the
resulting equations to (4), (6), and (7). If one wanted to go to order 4, an
additional set of four equations (for E[xt+14], E[xt+13xt], E[xt+12xt2], and E[xt+1xt3]) would be
needed, bringing the total to 10.More generally, in order to compute statistics of
ordern, one needs to construct and iterate(19)Q(n)=n×(n+1)2 second-order
difference equations. Since, after expansion, the r.h.s. of (2) contains 7
atomic terms of the form in (13), the r.h.s. of (17) contains 7b1 such
terms (Note that the exponent b2 does not influence
the number of terms because the recursion for E[xt+1b1xtb2] is
obtained as follows: (a) we compute xt+1b1, which is given by an expression containing 7b1 terms; (b) we
multiply each term by xtb2, which changes the exponents a1i but does not
alter the number of terms; (c) we apply the expectation operator, which again
does not modify the number of terms.) So the total number of terms one needs to
compute to construct the equations for order- n statistics is 7 (for the order
1 equations) plus (72+7) (for the order
2 equations) plus (73+72+7) (for the order
3 equations), and so forth. This gives us a total of(20)T(n)=∑i=1n(n−i+1)×7i terms. Note that T(n) grows
exponentially approximately as 1.36×7n. So although the number of equations one needs to
deal with grows quadratically, the computational effort required to instantiate
them is exponential. For instance, T(1)=7, T(2)=63, T(3)=462, T(4)=3262, and T(5)=22869. The growth in number of terms can be reduced if one
makes explicit use of ρ (i.e., by
adding the factor ρa8 in (13)). Then T(n)=O(6n). Either way, manually deriving equations for moments
of order 3 is already very laborious. The process, however, is clearly
mechanisable. This can be done using computer algebra systems, or by explicitly
representing and manipulating the ωi's in each
equation (this is what we did). As a result of mechanisation, computing the
equations for up to order 6 or 7 is feasible with an ordinary personal
computer.Some of theωi's in (17)
present the same pattern of exponents for w, c, y, and y^, so terms can be collected leading to more compact
equations (e.g., compare (5) and (6)). Also, given their size, one will
normally want to study (e.g., integrate) (17) numerically. In this case, w, c, y, and y^ are all
numerical parameters. So the ωi's become
constants and, after collecting terms, each equation contains at most Q(n) terms, which,
as we know, is quadratic in the order n. As a result, although the complexity of the
construction of the motion equations for the moments is exponential in the
order of the moments, their numerical integration is only of order O(n4).Naturally, the system ofQ(n) second-order
difference equations necessary to predict the dynamics of moments of order 1 to n can be turned
into a system of order 1 of the form in (8), via the choice(21)z(t)=(E[xt]E[xt−1]E[xt2]E[xtxt−1]E[xt−12]E[xt3]E[xt2xt−1]E[xtxt−12]E[xt−13]⋮E[xt−1n]).This effectively means adding artificial update
equations of the formE[xtk]=E[xtk] for k=1,…,n, bringing the total to Q′(n)=Q(n)+n. The transition matrix for the system is, therefore,
of size Q′(n)×Q′(n). We will denote this with Mn. For example, for n=4, which would allow one to study the mean, variance,
skewness, and kurtosis of the sampling distribution as a function of t, M4 is merely a 14×14 matrix.Interestingly,Q′(n) grows so slowly
that one can perform an eigenvalue analysis for any Mn that one is
able to compute. That is, the expensive part of the process is the construction
of Mn. Once this is done, iterating the system,
establishing its stability, or finding its fixed points is a trivial matter.In the next section, we provide results for statistics
of order 3 and 4, that is,n=4, a value of n for which
computing Mn takes only a
few seconds. However, before we do this, we need to consider the initial
conditions for the system.


## 3.3. Initial Conditions

In order to provide initial conditions for the
moments' dynamics, we need to computeE[x0k] and E[x1kx0l] for generic k>0 and l≥0.If we are at the very first iteration of the PSO
algorithm, the initial conditions for the dynamics of the moments are related
to the ranges used to initialise positions and velocities in the PSO. Under the
assumption that a particle's initial position,x0, is chosen uniformly at random in a symmetric range [−Ω,Ω], we have(22)E[x0k]={0ifk⁢is odd,Ωkk+1otherwise.  In order to
compute E[x1kx0l], we need to consider the equation(23)x1=x0+wv0−x0(ϕ1+ϕ2)+ϕ1y+ϕ2y^,  where a
particle's initial velocity, v0, is a stochastic variable uniformly distributed the
range [−Ωv,Ωv] (often Ωv=Ω). By taking
the kth power of
both sides of the equation, multiplying by x0l, and taking expectations, as we did to construct
(17), one obtains the desired expressions for E[x1kx0l]. Like for (17), these expressions contain a number of
terms that grows exponentially with n. However, this process, too, can be trivially
mechanised.If the PSO is not at its first iteration, and a new
improvement has just been found to a particle's personal best or neighbourhood
best, thenx0 is not a
stochastic variable, but a constant determined by the position where the
particle found the last improvement.


## 4. Skewness and Kurtosis of the PSO's Sampling Distribution

We constructed the recursions for moments of up to
order 4 as described in the previous section for the canonical PSO. Naturally,
for higher-order moments, we can do exactly the same type of eigenvalue
analysis we did in [16] for the mean and standard deviation of the sampling
distribution. For example, Figure 1 shows the magnitude of the largest
eigenvalue, Λm, of M4 as a function
of the parameters w and c. The system is order-4 stable, that is, mean,
variance, skewness, and kurtosis have a stable fixed point, whenever Λm<1, that is, within the curved region enclosed by the
contour shown in the figure.Figure 1
Magnitude of the largest
eigenvalue ofM4 as a function
of the parameters w and c. The curved line on the surface encloses the order-4
stable region.For easier comparison, we show the lines whereΛm=1 for M1, M2, M3, and M4 in Figure 2
(ordered from top to bottom, resp.). The regions of order-1, -2, -3, and
-4 stability are nested. Note how the Λm=1 lines for M2 and M3 coincide for
many values of w. Note also that the standard setting, c=1.49618 and w=0.7298, lays within the narrow region of order-3 stability.
This implies that while mean, variance, and skewness of the standard PSO tend
to a fixed point, kurtosis is unstable and will tend to grow indefinitely.
Interestingly, a growth in the kurtosis of samples was observed by Kennedy
[17], although this was effectively computed under the assumption that the
sampling distribution is time independent. So the values of xt recorded in a
run at t grows were
treated as different samples from the same distribution, while we now know this
is incorrect.Figure 2
Plot of the regions
of order-1, -2, -3, and -4 stability for the canonical PSO.That the predictions of the model are exact is also
confirmed by the comparison of the dynamics of predicted and recorded
higher-order moments. Figure3(a) shows a comparison between the skewness E[(xt−μt)3]/σt3 computed using
our model and the average positions of the particle recorded in one billion (1
000 000 000) real runs in the first 30 iterations for the case c=1.49618, w=0.7298, y=0, y^=1, and Ω=5. As one can see, there is a very good match between
the model's predictions and the stagnation behaviour of particles in real runs.
Only after about 20–25 generations, the sampling errors start accumulating
enough to show significant differences. We know that these differences are due
to sampling errors because before performing 1 000 000 000 runs, we attempted
to use much smaller sets of runs. With fewer runs, we observed that errors were
visible well before generation 25. For example, with 1 000 000 runs, there is a
good match up to around iteration 15.Comparison between predicted
and experimental skewness and (excess) kurtosis of the PSO sampling
distribution forc=1.49618, w=0.7298, y=0, y^=1, and Ω=5. Kurtosis grows exponentially, and so it is plotted
on a logarithmic scale. The first point is not plotted because the excess
kurtosis was negative (−1.2).
(a)(b)As shown in Figure3(b), the model also predicts
very well the behaviour of the (excess) kurtosis E[(xt−μt)4]/σt2−3 of the sampling
distribution. (Following standard practice, in this paper whenever we use the
term “kurtosis,” we will refer to the excess kurtosis E[(xt−μt)4]/σt2−3. The excess kurtosis of the normal distribution to
0.) Note that for c=1.49618 and w=0.7298, the system is order-3 stable, and so, although the
oscillations of the skewness shown in Figure 3(a) appear to grow bigger and
bigger, suggesting instability, this is actually only a transient effect, as
shown in Figure 4 where we integrate the equations over 200 generations instead
of 30.Figure 4
Predicted dynamics of
the skewness of the PSO sampling distribution over 200 generations forc=1.49618, w=0.7298, y=0, y^=1, and Ω=5.


## 5. Comparison between PSO's

In the previous sections, we studied the canonical PSO
with the restriction that the acceleration coefficients,c1 and c2, were identical: c1=c2=c. One may wonder, however, whether allowing such
coefficients to differ would produce qualitatively very different dynamics. For
example, what would happen if we set one of the ci's
to
zero as in a purely cognitive or purely social PSO model? This effectively
would reduce to one the sources of random influences on a particle's dynamics.
Conversely, one might wonder what would happen if we increased the number of
such sources of influence, as is done, for example, in the fully informed
particle swarm (FIPS) [18, 19].


### 5.1. Moments of Generalised PSO's

To answer these (and other) important questions on the
sampling distribution of different PSO models, we adopt a FIPS-like general
class of PSO's described by the following difference equation:(24)xt+1=xt+wvt+∑i=1mϕi(y^i−xt),  where the ϕi's are
stochastic variables uniformly distributed in the range [0,ci], ci are constants, and
the y^i's are the
personal best positions of neighbours of the particle (the particle itself may
be included in its own neighbourhood). Naturally, this equation can be
converted into the following:(25)xt+1=xt(1+w)−wxt−1−∑i=1mϕixt+∑iϕiy^i,  which is a generalisation
of (2). All of the steps we performed in Section 3 can be repeated for (25).
These lead to recursions of the form in (17) with the only difference that the
coefficients ωi take the more
general form(26)ωi=(a0iwawic1ac1i⋯cmacmiy^1ay1i⋯y^maymi∏j=1m(1+acji)),  where a0i, awi, ac1i,…,acmi, ay1i,…,aymi are appropriate
constants.Because (25) contains 3+2×m terms, the
complexity of the expansion now grows as O((3+2×m)n), where n is the order of
the moments we are interested in. So the larger m, the heavier the computation load required to compute Mn. Once the transition matrices Mn are computed,
however, they are exactly of the same size for all PSO models within the class
defined by (24). Initial conditions can be found following the approach
described in Section 3. Calculations are expensive but can be mechanised. We
did this for the examples described below.


### 5.2. Three PSO's We Want to Compare

An extensive comparison of different PSO's is beyond
the scope of this paper. However, as an example of the kind of comparisons one
can make using our approach, we considered the PSO's in (24) with N=3. Within this class of PSO's, we considered three
variants as follows:(a)
a purely social
variant of PSO, which we will callsocial PSO for brevity, where c1>0 and c2=c3=0 (due to
symmetries, the behaviour of a purely cognitive PSO where c2>0 and c1=c3=0 is effectively
identical to that of this social PSO);(b)
thecanonical PSO we have studied so far in
the paper, which is obtained by setting c1=c2>0 and c3=0;(c)
the simplest
version of FIPS with a neighbourhood of three individuals, for example,
obtained using anlbest topology
and an interaction radius of 1, where c1=c2=c3>0. We will call this version FIPS3.In order to perform a fair comparison of the stability
properties of these PSO variants, we study them in conditions where the sum of
the amplitudes of the random components,ϕi, is identical across models. That is, we set c=∑ici, and compare models with
the same c value. Again,
we analyse eigenvalues.


### 5.3. Results of the Comparison

Figures5–8 show the lines in the (w,c) plane where the
magnitude of the largest eigenvalue of Mn , Λm is 1 for n=1,2,3,4 and for the
three PSO variants mentioned above. Each figure represents a different moment.
In the plot for the nth moment, the
region below each line represents the set of all (w,c) pairs for which
the nth moment of
the distribution of the PSO corresponding to that line is stable. Let us
analyse these figures in detail.Figure 5
Lines below which the
mean of the sampling distributions for a social PSO, a canonical PSO, and FIPS
with a neighbourhood of three individuals have a fixed point (order-1
stability). The three lines coincide.Figure 6
Lines below which
the variance of the sampling distributions for a social PSO (bottom), a
canonical PSO (middle), and FIPS with a neighbourhood of three individuals
(top) have a fixed point (order-2 stability).Figure 7
Lines below which a
social PSO, a canonical PSO, and FIPS3 are order-3 stable (from bottom to top,
resp.).Figure 8
Lines below which
the kurtosis of the sampling distributions for a social PSO (bottom), a
canonical PSO (middle), and FIPS3 (top) have a fixed point (order-4 stability).Firstly, we should note that the regions of order-1
stability for the three models are identical. This is because the dynamics of
the mean of the three models is governed by equations of the same form, namely,(27)E[xt+1]=E[xt](1+w−c2)−wE[xt−1]+constant,  where the
constant term may differ in different PSO variants. (This is irrelevant for the
stability of the system, since stability is determined by the homogeneous part
of the equation.) Note also that the rightmost point in each plot is an
artifact due to the fact that, at w=1, Λm=1 for M1 irrespective of
the value of c.The regions where the variance is stable for the three
models, instead, are different, with FIPS3 having the largest region of order-2
stability, followed by the canonical PSO, and finally, by the social PSO.
Exactly the same happens with skewness (Figure7) and kurtosis (Figure 8), with
the order-3 stability region largely coinciding with the order-2 ones also for
FIPS3 and the social PSO.These results are counter intuitive. One would expect
that the more sources of randomness, theϕi's, there are,
the more a PSO should be unstable. However, the exact opposite happens. The
social PSO, where the only influence is ϕ1, is the least stable of all models, while FIPS3,
which has three sources of randomness, is the most stable. What are the reasons
for this behaviour?We can understand this by rewriting (25) as follows:(28)xt+1=xt(1+w)−wxt−1−xtΦm+Ψm,  where Φm=∑i=1mϕi and Ψm=∑iϕiy^i. Both Φm and Ψm are the sum of
independent and uniformly distributed variables: the variables ϕi in the case of Φm and the
variables y^iϕi in the case of Ψm.We know that∑ici=c. To simplify our treatment, let us further assume
that the ϕi's are i.i.d.,
that is, that all ci are identical,
and so ci=c/m. We can then apply the central limit theorem to Φm. This predicts that for sufficiently large m, the distribution of Φm is
approximately Gaussian with mean ∑ici/2=c/2 and variance ∑i(ci2/3−(ci/2)2)=∑ici2/12=c2/(12m). So the larger m, the smaller the variance of Φm and the
stochasticity of (28).In the case of the stochastic variableΨm, the quantities ϕiy^i are not
identically distributed even if all ci are identical.
This is because, in principle, each y^i may be
different. This prevents the use of the standard central limit theorem. We can,
however, apply Lyapunov's central limit theorem to Ψm. The conditions for its application are as
follows:(1)
the variablesϕiy^i must have
finite mean, which is the case since μi=E[ϕiy^i]=ciy^i/2=cy^i/(2m),(2)
the variablesϕiy^i must have
finite variance, which, again, is the case since σi2=E[(ϕiy^i−μi))2]=(ciy^i)2/12=(cy^i/m)2/12,(3)
the variablesϕiy^i must have
finite third central moment, which is satisfied since ri3=E[(ϕiy^i−μi))3]=0, and finally,(4)
the Lyapunov
condition,lim⁡m→∞((∑i=1mri3)1/3/(∑i=1mσi2)1/2)=0, must be satisfied, which, again, is the case since
all ri3=0.Then for
sufficiently largem, also the distribution of Ψm is
approximately Gaussian with mean ∑icy^i/(2m)=(c/2)×(∑iy^i/m) and variance (c2/12m)×(∑iy^i2/m). Note that ∑iy^i/m and ∑iy^i2/m are the mean y^i and the mean y^i2, respectively. So these are finite quantities if, as
is normally the case, all y^i are finite.
(PSO search is normally confined to a prefixed, finite region of ℝN, and so all y^i must be
finite.) So like for Φm, the larger m, the smaller the variance of Ψm, and consequently, the less the stochasticity of
(28).Effectively, the largerm, the more Φm and Ψm become
deterministic and approach constant values. This explains why adding more and
more sources of randomness—while keeping c constant—produces progressively more and more stable PSO's.


## 5.1. Moments of Generalised PSO's

To answer these (and other) important questions on the
sampling distribution of different PSO models, we adopt a FIPS-like general
class of PSO's described by the following difference equation:(24)xt+1=xt+wvt+∑i=1mϕi(y^i−xt),  where the ϕi's are
stochastic variables uniformly distributed in the range [0,ci], ci are constants, and
the y^i's are the
personal best positions of neighbours of the particle (the particle itself may
be included in its own neighbourhood). Naturally, this equation can be
converted into the following:(25)xt+1=xt(1+w)−wxt−1−∑i=1mϕixt+∑iϕiy^i,  which is a generalisation
of (2). All of the steps we performed in Section 3 can be repeated for (25).
These lead to recursions of the form in (17) with the only difference that the
coefficients ωi take the more
general form(26)ωi=(a0iwawic1ac1i⋯cmacmiy^1ay1i⋯y^maymi∏j=1m(1+acji)),  where a0i, awi, ac1i,…,acmi, ay1i,…,aymi are appropriate
constants.Because (25) contains 3+2×m terms, the
complexity of the expansion now grows as O((3+2×m)n), where n is the order of
the moments we are interested in. So the larger m, the heavier the computation load required to compute Mn. Once the transition matrices Mn are computed,
however, they are exactly of the same size for all PSO models within the class
defined by (24). Initial conditions can be found following the approach
described in Section 3. Calculations are expensive but can be mechanised. We
did this for the examples described below.


## 5.2. Three PSO's We Want to Compare

An extensive comparison of different PSO's is beyond
the scope of this paper. However, as an example of the kind of comparisons one
can make using our approach, we considered the PSO's in (24) with N=3. Within this class of PSO's, we considered three
variants as follows:(a)
a purely social
variant of PSO, which we will callsocial PSO for brevity, where c1>0 and c2=c3=0 (due to
symmetries, the behaviour of a purely cognitive PSO where c2>0 and c1=c3=0 is effectively
identical to that of this social PSO);(b)
thecanonical PSO we have studied so far in
the paper, which is obtained by setting c1=c2>0 and c3=0;(c)
the simplest
version of FIPS with a neighbourhood of three individuals, for example,
obtained using anlbest topology
and an interaction radius of 1, where c1=c2=c3>0. We will call this version FIPS3.In order to perform a fair comparison of the stability
properties of these PSO variants, we study them in conditions where the sum of
the amplitudes of the random components,ϕi, is identical across models. That is, we set c=∑ici, and compare models with
the same c value. Again,
we analyse eigenvalues.


## 5.3. Results of the Comparison

Figures5–8 show the lines in the (w,c) plane where the
magnitude of the largest eigenvalue of Mn , Λm is 1 for n=1,2,3,4 and for the
three PSO variants mentioned above. Each figure represents a different moment.
In the plot for the nth moment, the
region below each line represents the set of all (w,c) pairs for which
the nth moment of
the distribution of the PSO corresponding to that line is stable. Let us
analyse these figures in detail.Figure 5
Lines below which the
mean of the sampling distributions for a social PSO, a canonical PSO, and FIPS
with a neighbourhood of three individuals have a fixed point (order-1
stability). The three lines coincide.Figure 6
Lines below which
the variance of the sampling distributions for a social PSO (bottom), a
canonical PSO (middle), and FIPS with a neighbourhood of three individuals
(top) have a fixed point (order-2 stability).Figure 7
Lines below which a
social PSO, a canonical PSO, and FIPS3 are order-3 stable (from bottom to top,
resp.).Figure 8
Lines below which
the kurtosis of the sampling distributions for a social PSO (bottom), a
canonical PSO (middle), and FIPS3 (top) have a fixed point (order-4 stability).Firstly, we should note that the regions of order-1
stability for the three models are identical. This is because the dynamics of
the mean of the three models is governed by equations of the same form, namely,(27)E[xt+1]=E[xt](1+w−c2)−wE[xt−1]+constant,  where the
constant term may differ in different PSO variants. (This is irrelevant for the
stability of the system, since stability is determined by the homogeneous part
of the equation.) Note also that the rightmost point in each plot is an
artifact due to the fact that, at w=1, Λm=1 for M1 irrespective of
the value of c.The regions where the variance is stable for the three
models, instead, are different, with FIPS3 having the largest region of order-2
stability, followed by the canonical PSO, and finally, by the social PSO.
Exactly the same happens with skewness (Figure7) and kurtosis (Figure 8), with
the order-3 stability region largely coinciding with the order-2 ones also for
FIPS3 and the social PSO.These results are counter intuitive. One would expect
that the more sources of randomness, theϕi's, there are,
the more a PSO should be unstable. However, the exact opposite happens. The
social PSO, where the only influence is ϕ1, is the least stable of all models, while FIPS3,
which has three sources of randomness, is the most stable. What are the reasons
for this behaviour?We can understand this by rewriting (25) as follows:(28)xt+1=xt(1+w)−wxt−1−xtΦm+Ψm,  where Φm=∑i=1mϕi and Ψm=∑iϕiy^i. Both Φm and Ψm are the sum of
independent and uniformly distributed variables: the variables ϕi in the case of Φm and the
variables y^iϕi in the case of Ψm.We know that∑ici=c. To simplify our treatment, let us further assume
that the ϕi's are i.i.d.,
that is, that all ci are identical,
and so ci=c/m. We can then apply the central limit theorem to Φm. This predicts that for sufficiently large m, the distribution of Φm is
approximately Gaussian with mean ∑ici/2=c/2 and variance ∑i(ci2/3−(ci/2)2)=∑ici2/12=c2/(12m). So the larger m, the smaller the variance of Φm and the
stochasticity of (28).In the case of the stochastic variableΨm, the quantities ϕiy^i are not
identically distributed even if all ci are identical.
This is because, in principle, each y^i may be
different. This prevents the use of the standard central limit theorem. We can,
however, apply Lyapunov's central limit theorem to Ψm. The conditions for its application are as
follows:(1)
the variablesϕiy^i must have
finite mean, which is the case since μi=E[ϕiy^i]=ciy^i/2=cy^i/(2m),(2)
the variablesϕiy^i must have
finite variance, which, again, is the case since σi2=E[(ϕiy^i−μi))2]=(ciy^i)2/12=(cy^i/m)2/12,(3)
the variablesϕiy^i must have
finite third central moment, which is satisfied since ri3=E[(ϕiy^i−μi))3]=0, and finally,(4)
the Lyapunov
condition,lim⁡m→∞((∑i=1mri3)1/3/(∑i=1mσi2)1/2)=0, must be satisfied, which, again, is the case since
all ri3=0.Then for
sufficiently largem, also the distribution of Ψm is
approximately Gaussian with mean ∑icy^i/(2m)=(c/2)×(∑iy^i/m) and variance (c2/12m)×(∑iy^i2/m). Note that ∑iy^i/m and ∑iy^i2/m are the mean y^i and the mean y^i2, respectively. So these are finite quantities if, as
is normally the case, all y^i are finite.
(PSO search is normally confined to a prefixed, finite region of ℝN, and so all y^i must be
finite.) So like for Φm, the larger m, the smaller the variance of Ψm, and consequently, the less the stochasticity of
(28).Effectively, the largerm, the more Φm and Ψm become
deterministic and approach constant values. This explains why adding more and
more sources of randomness—while keeping c constant—produces progressively more and more stable PSO's.


## 6. The Density Function of the Canonical PSO's Sampling Distribution

The technique described in this paper, in principle,
would allow one to determine all the moments of the sampling distribution of
the PSO at all times. The question then is “could we derive the PSO
sampling distribution itself ?” The answer is of course in the positive
since knowing all the moments of a distribution implies knowing its moment
generating function. This, in turn, allows one to obtain the density function
of the distribution via inverse Laplace transform.In practice, however, it is impossible to compute all
the moments of the PSO sampling distribution. This is for two reasons. Firstly,
there are infinitely many such moments. Secondly, as we have seen in the
previous sections, the cost of computing moments is exponential in the order of
the moments. The next question is then “to what extent can we still
reconstruct the PSO's density function from a finite number of moments?”
This is an instance of the well-knowntruncated moment problem, a difficult,
inverse ill-posed problem for which many approaches have been proposed. Here,
we consider only one such approach.A particularly simple idea is to consider a family of
density functionsf(x;λ1,λ2,…) with parameters λ1, λ2, and so forth,
with sufficient expressive power to represent distributions with widely
different shapes, with more or less asymmetries, with tails of different
characteristics, and so on. Then one can use optimisation techniques to
identify the parameters of the distribution f which minimise
the difference between the moments of f and the moments
of the PSO's sampling distribution. This is called the moment-matching method. Once the
parameters λ1, λ2, and so forth, are identified, f can be used as
an approximation of the true PSO sampling distribution. This approach to
reconstructing probability distributions from moments was proposed [20] (see
also [21, 22]), where a generalised
lambda distribution (GLD) was used. We adopt this same approach here.GLD is a four-parameter distribution defined via its
quantile function as follows:(29)Q(u)=λ1+1λ2(uλ3−1λ3+1−(1−u)λ4λ4),  where u∈[0,1]. Its density function is given by(30)f(x;λ1,λ2,λ3,λ4)=(dQ(u)du)−1=λ2u(λ3−1)+(1−u)(λ4−1), where u=Q−1(x).The GLD is enormously flexible in terms of the shape
of the distribution. For example, as shown in Figure9, the uniform, Gaussian,
exponential, and Gamma distributions are all special
cases of GLD. Effectively, λ1 determines the
location of the distribution, λ2 determines its
scale, while λ3 and λ4 determine other
shape characteristics. In particular, only if λ3=λ4, the distribution is symmetric.Figure 9
Examples of GLD density functions.Because GLD has 4 parameters, all we need
is four moments—the mean, variance, skewness,
and kurtosis—of the PSO's sampling distribution in order to identify such
parameters with the moment-matching method described above.As an illustration, we apply this technique to
reconstruct the sampling distribution during stagnation of a canonical PSO with
parametersc=1.49618, w=0.7298, y=0, y^=10, and Ω=5. With these parameters, lim⁡t→∞E[xt]=5. In Figure 10, we show snapshots at times t=0,2,4,12, and 24 of the theoretical sampling
distribution together with estimates of the distribution based on 1 000 000
actual runs. In all cases, the match between the moments of the GLD and those
of the PSO sampling distribution was exact (within experimental errors). Also,
there is considerable agreement between the theoretical lines and histograms
obtained in real runs. Note how widely the mean of the density function
oscillates around 5 in the first few generations. Also note the asymmetry in
the distributions due to the oscillations of the skewness.Figure 10
Estimates of the
sampling distribution of a canonical PSO with parametersc=1.49618, w=0.7298, y=0, y^=10, and Ω=5 during
stagnation, reconstructed via GLD best fitting vers. histograms over 1 000 000
real runs. Snapshots at times t=0, 2, 4, 12, and 24 are shown. For each theoretical
sampling distribution, we report the parameters of the corresponding GLD.


## 7. Discussion and Conclusions

Several theoretical analyses of the dynamics of
particle swarms have been offered in the literature over the last decade. These
have been very illuminating. However, virtually, all have relied on substantial
simplifications and on the assumption that the particles are deterministic.
Naturally, these simplifications make it impossible to derive an exact characterisation
of the sampling distribution of the PSO.In previous work [16], by using surprisingly simple
techniques, we started by exactly determining, perhaps, the most important
characteristic of a PSO's sampling distribution, its variance, and we have been
able to explain how it changes over any number of generations. The only
assumption we made is stagnation. Here, we extended this technique to the study
of higher-order statistics. In particular, we analysed, in detail, the skewness
and kurtosis of the distribution. Because of the complexity of the calculations
involved, this required mechanising the derivation of the recursions for these
moments.We applied the analysis to the PSO with inertia
weight, but the analysis is also valid for the PSO with constriction, because
of the well-known equivalence of these two models (via a simple parameter
mapping). We also generalised our model so as to include FIPS. This made it
possible to explicitly compare the stability of different forms of PSO, leading
to a deeper understanding of their properties. In particular, we showed that
while FIPS and standard forms of PSO present exactly the same order-1
stability, in FIPS, higher-order moments are more stable than in the other
PSO's, and we were able to explain why this is the case, using two forms of
central limit theorem (Section5). Elsewhere [23], we also applied the same
type of analysis to derive the stability regions of a simpler form of
velocity-based PSO.Finally, with all these tools in hand, we went in
search for the “holy grail”—the actual PSO sampling density function. We
treated the problem as an ill-posed inverse problem, which we regularised
thanks to the use (and best fit) of a family of distributions—the
generalised lambda distribution (GLD). All empirical evidence we have suggests
that this distribution approximates very closely the sampling behaviour of
PSO's (naturally, with parameters that are functions of time, i.e.,λi=λi(t)). We will
attempt to characterise the PSO's sampling distribution in more depth in future
research.Whether or not GLD is the exact PSO sampling
distribution or just a very good approximation, if one could determine (again,
either exactly or approximately) how theλi(t)'s are affected
by the parameters c, w, y, and y^ and by the
initial conditions x0 and v0, it would then be possible to accurately simulate the
behaviour of the PSO by sampling from f(x;λ1(t),λ2(t),λ3(t),λ4(t)). This is easily done since GLD deviates can trivially
be produced by picking u uniformly at
random in the interval [0,1] and applying
(29), that is, Q(U[0,1]) is generalised
lambda distributed. We will study this more sophisticated form of bare-bones
PSO [17] in future research.Bare-bone PSO's have some resemblance with estimation
of distribution algorithms (EDA's). EDA's are powerful population-based searchers
where the variation operations traditionally implemented via crossover and
mutation in evolutionary algorithms are replaced by the process of sampling
from a distribution (a review of EDA's is available in [24]). The distribution
is modified generation after generation, by using information obtained from the
more fit individuals in the population. The objective of these changes is to
increase the probability of generating individuals with high fitness. As one
can clearly see from this description, a bare-bone PSO performs a very similar
form of search to an EDA. However, there is a difference between the two. While
in a typical EDA the whole new population is generated from a single
distribution, in a bare-bone PSO, each individual has a personal sampling
distribution from which only one sample is drawn at each iteration. (The
sampling distributions of different particles may share some characteristics,
since the sampling distribution of a particle is determined by both personal
best and neighbourhood best, and the neighbourhood best can easily be the same
for different particles.) Nonetheless, one could still interpret a bare-bone
PSO as an unusual type of EDE in which multiple distributions are used and
updated based on fitness during the search.



---

*Source: 761459-2008-03-30.xml*