QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
particleswarmoptimization.hpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4Copyright (C) 2015 Andres Hernandez
5
6This file is part of QuantLib, a free-software/open-source library
7for financial quantitative analysts and developers - http://quantlib.org/
8
9QuantLib is free software: you can redistribute it and/or modify it
10under the terms of the QuantLib license. You should have received a
11copy of the license along with this program; if not, please email
12<quantlib-dev@lists.sf.net>. The license is also available online at
13<http://quantlib.org/license.shtml>.
14
15This program is distributed in the hope that it will be useful, but WITHOUT
16ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file particleswarmoptimization.hpp
21\brief Implementation based on:
22Clerc, M., Kennedy, J. (2002) The particle swarm-explosion, stability and
23convergence in a multidimensional complex space. IEEE Transactions on Evolutionary
24Computation, 6(2): 58–73.
25*/
26
27#ifndef quantlib_optimization_particleswarmoptimization_hpp
28#define quantlib_optimization_particleswarmoptimization_hpp
29
36
37#include <random>
38
39namespace QuantLib {
40
41 /*! The process is as follows:
42 M individuals are used to explore the N-dimensional parameter space:
43 \f$ X_{i}^k = (X_{i, 1}^k, X_{i, 2}^k, \ldots, X_{i, N}^k) \f$ is the kth-iteration for the ith-individual.
44
45 X is updated via the rule
46 \f[
47 X_{i, j}^{k+1} = X_{i, j}^k + V_{i, j}^{k+1}
48 \f]
49 with V being the "velocity" that updates the position:
50 \f[
51 V_{i, j}^{k+1} = \chi\left(V_{i, j}^k + c_1 r_{i, j}^k (P_{i, j}^k - X_{i, j}^k)
52 + c_2 R_{i, j}^k (G_{i, j}^k - X_{i, j}^k)\right)
53 \f]
54 where c are constants, r and R are uniformly distributed random numbers in the range [0, 1], and
55 \f$ P_{i, j} \f$ is the personal best parameter set for individual i up to iteration k
56 \f$ G_{i, j} \f$ is the global best parameter set for the swarm up to iteration k.
57 \f$ c_1 \f$ is the self recognition coefficient
58 \f$ c_2 \f$ is the social recognition coefficient
59
60 This version is known as the PSO with constriction factor (PSO-Co).
61 PSO with inertia factor (PSO-In) updates the velocity according to:
62 \f[
63 V_{i, j}^{k+1} = \omega V_{i, j}^k + \hat{c}_1 r_{i, j}^k (P_{i, j}^k - X_{i, j}^k)
64 + \hat{c}_2 R_{i, j}^k (G_{i, j}^k - X_{i, j}^k)
65 \f]
66 and is accessible from PSO-Co by setting \f$ \omega = \chi \f$,
67 and \f$ \hat{c}_{1,2} = \chi c_{1,2} \f$.
68
69 These two versions of PSO are normally referred to as canonical PSO.
70
71 Convergence of PSO-Co is improved if \f$ \chi \f$ is chosen as
72 \f$ \chi = \frac{2}{\vert 2-\phi-\sqrt{\phi^2 - 4\phi}\vert} \f$,
73 with \f$ \phi = c_1 + c_2 \f$.
74 Stable convergence is achieved if \f$ \phi >= 4 \f$. Clerc and Kennedy recommend
75 \f$ c_1 = c_2 = 2.05 \f$ and \f$ \phi = 4.1 \f$.
76
77 Different topologies can be chosen for G, e.g. instead of it being the best
78 of the swarm, it is the best of the nearest neighbours, or some other form.
79
80 In the canonical PSO, the inertia function is trivial. It is simply a
81 constant (the inertia) multiplying the previous iteration's velocity. The
82 value of the inertia constant determines the weight of a global search over
83 local search. Like in the case of the topology, other possibilities for the
84 inertia function are also possible, e.g. a function that interpolates between a
85 high inertia at the beginning of the optimization (hence prioritizing a global
86 search) and a low inertia towards the end of the optimization (hence prioritizing
87 a local search).
88
89 The optimization stops either because the number of iterations has been reached
90 or because the stationary function value limit has been reached.
91 */
93 public:
94 class Inertia;
95 class Topology;
97 ext::shared_ptr<Topology> topology,
98 ext::shared_ptr<Inertia> inertia,
99 Real c1 = 2.05,
100 Real c2 = 2.05,
101 unsigned long seed = SeedGenerator::instance().get());
103 ext::shared_ptr<Topology> topology,
104 ext::shared_ptr<Inertia> inertia,
105 Real omega,
106 Real c1,
107 Real c2,
108 unsigned long seed = SeedGenerator::instance().get());
109 void startState(Problem &P, const EndCriteria &endCriteria);
110 EndCriteria::Type minimize(Problem& P, const EndCriteria& endCriteria) override;
111
112 protected:
113 std::vector<Array> X_, V_, pBX_, gBX_;
119 ext::shared_ptr<Topology> topology_;
120 ext::shared_ptr<Inertia> inertia_;
121 };
122
123 //! Base inertia class used to alter the PSO state
124 /*! This pure virtual base class provides the access to the PSO state
125 which the particular inertia algorithm will change upon each iteration.
126 */
129 public:
130 virtual ~Inertia() = default;
131 //! initialize state for current problem
132 virtual void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria) = 0;
133 //! produce changes to PSO state for current iteration
134 virtual void setValues() = 0;
135 protected:
137 std::vector<Array> *X_, *V_, *pBX_, *gBX_;
140
141 virtual void init(ParticleSwarmOptimization *pso) {
142 pso_ = pso;
143 X_ = &pso_->X_;
144 V_ = &pso_->V_;
145 pBX_ = &pso_->pBX_;
146 gBX_ = &pso_->gBX_;
147 pBF_ = &pso_->pBF_;
148 gBF_ = &pso_->gBF_;
149 lX_ = &pso_->lX_;
150 uX_ = &pso_->uX_;
151 }
152 };
153
154 //! Trivial Inertia
155 /* Inertia is a static value
156 */
158 public:
159 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
160 c0_ = c0;
161 M_ = M;
162 }
163 inline void setValues() override {
164 for (Size i = 0; i < M_; i++) {
165 (*V_)[i] *= c0_;
166 }
167 }
168
169 private:
172 };
173
174 //! Simple Random Inertia
175 /* Inertia value gets multiplied with a random number
176 between (threshold, 1)
177 */
179 public:
180 SimpleRandomInertia(Real threshold = 0.5, unsigned long seed = SeedGenerator::instance().get())
181 : threshold_(threshold), rng_(seed) {
182 QL_REQUIRE(threshold_ >= 0.0 && threshold_ < 1.0, "Threshold must be a Real in [0, 1)");
183 }
184 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
185 M_ = M;
186 c0_ = c0;
187 }
188 inline void setValues() override {
189 for (Size i = 0; i < M_; i++) {
190 Real val = c0_*(threshold_ + (1.0 - threshold_)*rng_.nextReal());
191 (*V_)[i] *= val;
192 }
193 }
194
195 private:
199 };
200
201 //! Decreasing Inertia
202 /* Inertia value gets decreased every iteration until it reaches
203 a value of threshold when iteration reaches the maximum level
204 */
206 public:
207 DecreasingInertia(Real threshold = 0.5)
208 : threshold_(threshold) {
209 QL_REQUIRE(threshold_ >= 0.0 && threshold_ < 1.0, "Threshold must be a Real in [0, 1)");
210 }
211 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
212 N_ = N;
213 c0_ = c0;
214 iteration_ = 0;
215 maxIterations_ = endCriteria.maxIterations();
216 }
217 inline void setValues() override {
219 for (Size i = 0; i < M_; i++) {
220 (*V_)[i] *= c0;
221 }
222 }
223
224 private:
227 };
228
229 //! AdaptiveInertia
230 /* Alen Lukic, Approximating Kinetic Parameters Using Particle
231 Swarm Optimization.
232 */
234 public:
235 AdaptiveInertia(Real minInertia, Real maxInertia, Size sh = 5, Size sl = 2)
236 :minInertia_(minInertia), maxInertia_(maxInertia),
237 sh_(sh), sl_(sl) {};
238 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
239 M_ = M;
240 c0_ = c0;
241 adaptiveCounter = 0;
243 started_ = false;
244 }
245 void setValues() override;
246
247 private:
254 };
255
256 //! Levy Flight Inertia
257 /* As long as the particle keeps getting frequent updates to its
258 personal best value, the inertia behaves like a SimpleRandomInertia,
259 but after a number of iterations without improvement, the behaviour
260 changes to that of a Levy flight ~ u^{-1/\alpha}
261 */
263 public:
265 unsigned long seed = SeedGenerator::instance().get())
267 1, Array(1, 1.0), seed),
268 threshold_(threshold) {};
269 inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
270 M_ = M;
271 N_ = N;
272 c0_ = c0;
273 adaptiveCounter_ = std::vector<Size>(M_, 0);
274 }
275 inline void setValues() override {
276 for (Size i = 0; i < M_; i++) {
277 if ((*pBF_)[i] < personalBestF_[i]) {
278 personalBestF_[i] = (*pBF_)[i];
279 adaptiveCounter_[i] = 0;
280 }
281 else {
282 adaptiveCounter_[i]++;
283 }
284 if (adaptiveCounter_[i] <= threshold_) {
285 //Simple Random Inertia
286 (*V_)[i] *= c0_*(0.5 + 0.5*rng_.nextReal());
287 }
288 else {
289 //If particle has not found a new personal best after threshold_ iterations
290 //then trigger a Levy flight pattern for the speed
291 flight_.nextReal<Real *>(&(*V_)[i][0]);
292 }
293 }
294 }
295
296 protected:
297 void init(ParticleSwarmOptimization* pso) override {
300 flight_.setDimension(N_, *lX_, *uX_);
301 }
302
303 private:
305 std::mt19937 generator_;
308 std::vector<Size> adaptiveCounter_;
312 };
313
314 //! Base topology class used to determine the personal and global best
315 /*! This pure virtual base class provides the access to the PSO state
316 which the particular topology algorithm will change upon each iteration.
317 */
320 public:
321 virtual ~Topology() = default;
322 //! initialize state for current problem
323 virtual void setSize(Size M) = 0;
324 //! produce changes to PSO state for current iteration
325 virtual void findSocialBest() = 0;
326 protected:
328 std::vector<Array> *X_, *V_, *pBX_, *gBX_;
330 private:
332 pso_ = pso;
333 X_ = &pso_->X_;
334 V_ = &pso_->V_;
335 pBX_ = &pso_->pBX_;
336 gBX_ = &pso_->gBX_;
337 pBF_ = &pso_->pBF_;
338 gBF_ = &pso_->gBF_;
339 }
340 };
341
342 //! Global Topology
343 /* The global best as seen by each particle is the best from amongst
344 all particles
345 */
347 public:
348 inline void setSize(Size M) override { M_ = M; }
349 inline void findSocialBest() override {
350 Real bestF = (*pBF_)[0];
351 Size bestP = 0;
352 for (Size i = 1; i < M_; i++) {
353 if (bestF < (*pBF_)[i]) {
354 bestF = (*pBF_)[i];
355 bestP = i;
356 }
357 }
358 Array& x = (*pBX_)[bestP];
359 for (Size i = 0; i < M_; i++) {
360 if (i != bestP) {
361 (*gBX_)[i] = x;
362 (*gBF_)[i] = bestF;
363 }
364 }
365 }
366
367 private:
369 };
370
371 //! K-Neighbor Topology
372 /* The global best as seen by each particle is the best from amongst
373 the previous K and next K neighbors. For particle I, the best is
374 then taken from amongst the [I - K, I + K] particles.
375 */
377 public:
378 KNeighbors(Size K = 1) :K_(K) {
379 QL_REQUIRE(K > 0, "Neighbors need to be larger than 0");
380 }
381 inline void setSize(Size M) override {
382 M_ = M;
383 QL_ENSURE(K_ < M, "Number of neighbors need to be smaller than total particles in swarm");
384 }
385 void findSocialBest() override;
386
387 private:
389 };
390
391 //! Clubs Topology
392 /* H.M. Emara, Adaptive Clubs-based Particle Swarm Optimization
393 Each particle is originally assigned to a default number of clubs
394 from among the total set. The best as seen by each particle is the
395 best from amongst the clubs to which the particle belongs.
396 Underperforming particles join more clubs randomly (up to a maximum
397 number) to widen the particles that influence them, while
398 overperforming particles leave clubs randomly (down to a minimum
399 number) to avoid early convergence to local minima.
400 */
402 public:
403 ClubsTopology(Size defaultClubs, Size totalClubs,
404 Size maxClubs, Size minClubs,
405 Size resetIteration, unsigned long seed = SeedGenerator::instance().get());
406 void setSize(Size M) override;
407 void findSocialBest() override;
408
409 private:
413 std::vector<std::vector<bool> > clubs4particles_;
414 std::vector<std::vector<bool> > particles4clubs_;
415 std::vector<Size> bestByClub_;
416 std::vector<Size> worstByClub_;
417 std::mt19937 generator_;
418 std::uniform_int_distribution<QuantLib::Size> distribution_;
420
421 void leaveRandomClub(Size particle, Size currentClubs);
422 void joinRandomClub(Size particle, Size currentClubs);
423 };
424
425}
426
427#endif
void setValues() override
produce changes to PSO state for current iteration
void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria) override
initialize state for current problem
AdaptiveInertia(Real minInertia, Real maxInertia, Size sh=5, Size sl=2)
1-D array used in linear algebra.
Definition: array.hpp:52
std::uniform_int_distribution< QuantLib::Size > distribution_
std::vector< std::vector< bool > > particles4clubs_
std::vector< std::vector< bool > > clubs4particles_
void setSize(Size M) override
initialize state for current problem
decltype(distribution_)::param_type param_type
void findSocialBest() override
produce changes to PSO state for current iteration
void joinRandomClub(Size particle, Size currentClubs)
void leaveRandomClub(Size particle, Size currentClubs)
void setValues() override
produce changes to PSO state for current iteration
void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria) override
initialize state for current problem
Criteria to end optimization process:
Definition: endcriteria.hpp:40
Size maxIterations() const
void setSize(Size M) override
initialize state for current problem
void findSocialBest() override
produce changes to PSO state for current iteration
void setSize(Size M) override
initialize state for current problem
void findSocialBest() override
produce changes to PSO state for current iteration
void setValues() override
produce changes to PSO state for current iteration
IsotropicRandomWalk< LevyFlightDistribution, std::mt19937 > flight_
LevyFlightInertia(Real alpha, Size threshold, unsigned long seed=SeedGenerator::instance().get())
void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria) override
initialize state for current problem
void init(ParticleSwarmOptimization *pso) override
Uniform random number generator.
Real nextReal() const
return a random number in the (0.0, 1.0)-interval
Abstract class for constrained optimization method.
Definition: method.hpp:36
Base inertia class used to alter the PSO state.
virtual void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria)=0
initialize state for current problem
virtual void init(ParticleSwarmOptimization *pso)
virtual void setValues()=0
produce changes to PSO state for current iteration
Base topology class used to determine the personal and global best.
virtual void setSize(Size M)=0
initialize state for current problem
virtual void findSocialBest()=0
produce changes to PSO state for current iteration
void startState(Problem &P, const EndCriteria &endCriteria)
EndCriteria::Type minimize(Problem &P, const EndCriteria &endCriteria) override
minimize the optimization problem P
Constrained optimization problem.
Definition: problem.hpp:42
void setValues() override
produce changes to PSO state for current iteration
SimpleRandomInertia(Real threshold=0.5, unsigned long seed=SeedGenerator::instance().get())
void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria) override
initialize state for current problem
static SeedGenerator & instance()
access to the unique instance
Definition: singleton.hpp:104
void setValues() override
produce changes to PSO state for current iteration
void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria) override
initialize state for current problem
Abstract constraint class.
#define QL_ENSURE(condition, message)
throw an error if the given post-condition is not verified
Definition: errors.hpp:130
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
#define QL_MAX_REAL
Definition: qldefines.hpp:176
QL_REAL Real
real number
Definition: types.hpp:50
std::size_t Size
size of a container
Definition: types.hpp:58
Isotropic random walk.
Levy Flight, aka Pareto Type I, distribution.
Mersenne Twister uniform random number generator.
Definition: any.hpp:35
Abstract optimization problem class.
Real alpha
Definition: sabr.cpp:200
Random seed generator.