QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
simulatedannealing.hpp
Go to the documentation of this file.
1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2013 Peter Caspers
5
6 This file is part of QuantLib, a free-software/open-source library
7 for financial quantitative analysts and developers - http://quantlib.org/
8
9 QuantLib is free software: you can redistribute it and/or modify it
10 under the terms of the QuantLib license. You should have received a
11 copy 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
15 This program is distributed in the hope that it will be useful, but WITHOUT
16 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17 FOR A PARTICULAR PURPOSE. See the license for more details.
18*/
19
20/*! \file simulatedannealing.hpp
21 \brief Numerical Recipes in C (second edition), Chapter 10.9,
22 with the original exit criterion in f(x) replaced by one
23 in x (see simplex.cpp for a reference to GSL concerning this)
24*/
25
26#ifndef quantlib_optimization_simulatedannealing_hpp
27#define quantlib_optimization_simulatedannealing_hpp
28
32#include <cmath>
33
34namespace QuantLib {
35
36 /*! Class RNG must implement the following interface:
37 \code
38 RNG::sample_type RNG::next() const;
39 \endcode
40
41 \ingroup optimizers
42 */
43
44 //! Simulated Annealing
45 template <class RNG = MersenneTwisterUniformRng>
47
48 public:
49
50 enum Scheme {
53 };
54
55 /*! reduce temperature T by a factor of \f$ (1-\epsilon) \f$ after m moves */
56 SimulatedAnnealing(const Real lambda, const Real T0,
57 const Real epsilon, const Size m,
58 const RNG &rng = RNG())
59 : scheme_(ConstantFactor), lambda_(lambda), T0_(T0),
60 epsilon_(epsilon), alpha_(0.0), K_(0), rng_(rng), m_(m) {}
61
62 /*! budget a total of K moves, set temperature T to the initial
63 temperature times \f$ ( 1 - k/K )^\alpha \f$ with k being the total number
64 of moves so far. After K moves the temperature is guaranteed to be
65 zero, after that the optimization runs like a deterministic simplex
66 algorithm.
67 */
68 SimulatedAnnealing(const Real lambda, const Real T0, const Size K,
69 const Real alpha, const RNG &rng = RNG())
70 : scheme_(ConstantBudget), lambda_(lambda), T0_(T0), epsilon_(0.0),
71 alpha_(alpha), K_(K), rng_(rng) {}
72
73 EndCriteria::Type minimize(Problem& P, const EndCriteria& ec) override;
74
75 private:
76
79 const Size K_;
80 const RNG rng_;
81
83 void amotsa(Problem &, Real);
84
86 std::vector<Array> vertices_;
93 };
94
95 template <class RNG>
97 // simplex.cpp
98 Array center(vertices_.front().size(), 0);
99 for (auto& vertice : vertices_)
100 center += vertice;
101 center *= 1 / Real(vertices_.size());
102 Real result = 0;
103 for (auto& vertice : vertices_) {
104 Array temp = vertice - center;
105 result += Norm2(temp);
106 }
107 return result / Real(vertices_.size());
108 }
109
110 template <class RNG>
112 fac1_ = (1.0 - fac) / ((Real)n_);
113 fac2_ = fac1_ - fac;
114 for (j_ = 0; j_ < n_; j_++) {
115 ptry_[j_] = sum_[j_] * fac1_ - vertices_[ihi_][j_] * fac2_;
116 }
117 if (!P.constraint().test(ptry_))
118 ytry_ = QL_MAX_REAL;
119 else
120 ytry_ = P.value(ptry_);
121 if (std::isnan(ytry_)) {
122 ytry_ = QL_MAX_REAL;
123 }
124 if (ytry_ <= yb_) {
125 yb_ = ytry_;
126 pb_ = ptry_;
127 }
128 yflu_ = ytry_ - tt_ * std::log(rng_.next().value);
129 if (yflu_ < yhi_) {
130 values_[ihi_] = ytry_;
131 yhi_ = yflu_;
132 for (j_ = 0; j_ < n_; j_++) {
133 sum_[j_] += ptry_[j_] - vertices_[ihi_][j_];
134 vertices_[ihi_][j_] = ptry_[j_];
135 }
136 }
137 ytry_ = yflu_;
138 }
139
140 template <class RNG>
142 const EndCriteria &ec) {
143
144 Size stationaryStateIterations_ = 0;
146 P.reset();
147 Array x = P.currentValue();
148 iteration_ = 0;
149 n_ = x.size();
150 ptry_ = Array(n_, 0.0);
151
152 // build vertices
153
154 vertices_ = std::vector<Array>(n_ + 1, x);
155 for (i_ = 0; i_ < n_; i_++) {
156 Array direction(n_, 0.0);
157 direction[i_] = 1.0;
158 P.constraint().update(vertices_[i_ + 1], direction, lambda_);
159 }
160 values_ = Array(n_ + 1, 0.0);
161 for (i_ = 0; i_ <= n_; i_++) {
162 if (!P.constraint().test(vertices_[i_]))
163 values_[i_] = QL_MAX_REAL;
164 else
165 values_[i_] = P.value(vertices_[i_]);
166 if (std::isnan(ytry_)) { // handle NAN
167 values_[i_] = QL_MAX_REAL;
168 }
169 }
170
171 // minimize
172
173 T_ = T0_;
174 yb_ = QL_MAX_REAL;
175 pb_ = Array(n_, 0.0);
176 do {
177 iterationT_ = iteration_;
178 do {
179 sum_ = Array(n_, 0.0);
180 for (i_ = 0; i_ <= n_; i_++)
181 sum_ += vertices_[i_];
182 tt_ = -T_;
183 ilo_ = 0;
184 ihi_ = 1;
185 ynhi_ = values_[0] + tt_ * std::log(rng_.next().value);
186 ylo_ = ynhi_;
187 yhi_ = values_[1] + tt_ * std::log(rng_.next().value);
188 if (ylo_ > yhi_) {
189 ihi_ = 0;
190 ilo_ = 1;
191 ynhi_ = yhi_;
192 yhi_ = ylo_;
193 ylo_ = ynhi_;
194 }
195 for (i_ = 2; i_ < n_ + 1; i_++) {
196 yt_ = values_[i_] + tt_ * std::log(rng_.next().value);
197 if (yt_ <= ylo_) {
198 ilo_ = i_;
199 ylo_ = yt_;
200 }
201 if (yt_ > yhi_) {
202 ynhi_ = yhi_;
203 ihi_ = i_;
204 yhi_ = yt_;
205 } else {
206 if (yt_ > ynhi_) {
207 ynhi_ = yt_;
208 }
209 }
210 }
211
212 // rtol_ = 2.0 * std::fabs(yhi_ - ylo_) /
213 // (std::fabs(yhi_) + std::fabs(ylo_));
214 // check rtol against some ftol... // NR end criterion in f(x)
215
216 // GSL end criterion in x (cf. above)
217 if (ec.checkStationaryPoint(simplexSize(), 0.0,
218 stationaryStateIterations_,
219 ecType) ||
220 ec.checkMaxIterations(iteration_, ecType)) {
221 // no matter what, we return the best ever point !
222 P.setCurrentValue(pb_);
223 P.setFunctionValue(yb_);
224 return ecType;
225 }
226
227 iteration_ += 2;
228 amotsa(P, -1.0);
229 if (ytry_ <= ylo_) {
230 amotsa(P, 2.0);
231 } else {
232 if (ytry_ >= ynhi_) {
233 ysave_ = yhi_;
234 amotsa(P, 0.5);
235 if (ytry_ >= ysave_) {
236 for (i_ = 0; i_ < n_ + 1; i_++) {
237 if (i_ != ilo_) {
238 for (j_ = 0; j_ < n_; j_++) {
239 sum_[j_] = 0.5 * (vertices_[i_][j_] +
240 vertices_[ilo_][j_]);
241 vertices_[i_][j_] = sum_[j_];
242 }
243 values_[i_] = P.value(sum_);
244 }
245 }
246 iteration_ += n_;
247 for (i_ = 0; i_ < n_; i_++)
248 sum_[i_] = 0.0;
249 for (i_ = 0; i_ <= n_; i_++)
250 sum_ += vertices_[i_];
251 }
252 } else {
253 iteration_ += 1;
254 }
255 }
256 } while (iteration_ <
257 iterationT_ + (scheme_ == ConstantFactor ? m_ : 1));
258
259 switch (scheme_) {
260 case ConstantFactor:
261 T_ *= (1.0 - epsilon_);
262 break;
263 case ConstantBudget:
264 if (iteration_ <= K_)
265 T_ = T0_ *
266 std::pow(1.0 - (Real)iteration_ / (Real)K_, alpha_);
267 else
268 T_ = 0.0;
269 break;
270 }
271
272 } while (true);
273 }
274}
275
276#endif
1-D array used in linear algebra.
Definition: array.hpp:52
Size size() const
dimension of the array
Definition: array.hpp:495
bool test(const Array &p) const
Definition: constraint.hpp:57
Real update(Array &p, const Array &direction, Real beta) const
Definition: constraint.cpp:27
Criteria to end optimization process:
Definition: endcriteria.hpp:40
bool checkStationaryPoint(Real xOld, Real xNew, Size &statStateIterations, EndCriteria::Type &ecType) const
Definition: endcriteria.cpp:64
bool checkMaxIterations(Size iteration, EndCriteria::Type &ecType) const
Definition: endcriteria.cpp:56
Abstract class for constrained optimization method.
Definition: method.hpp:36
Constrained optimization problem.
Definition: problem.hpp:42
const Array & currentValue() const
current value of the local minimum
Definition: problem.hpp:81
Constraint & constraint() const
Constraint.
Definition: problem.hpp:71
Real value(const Array &x)
call cost function computation and increment evaluation counter
Definition: problem.hpp:116
void setFunctionValue(Real functionValue)
Definition: problem.hpp:83
void setCurrentValue(const Array &currentValue)
Definition: problem.hpp:76
EndCriteria::Type minimize(Problem &P, const EndCriteria &ec) override
minimize the optimization problem P
SimulatedAnnealing(const Real lambda, const Real T0, const Real epsilon, const Size m, const RNG &rng=RNG())
SimulatedAnnealing(const Real lambda, const Real T0, const Size K, const Real alpha, const RNG &rng=RNG())
Abstract constraint class.
const Matrix m_
Definition: expm.cpp:49
#define QL_MAX_REAL
Definition: qldefines.hpp:176
QL_REAL Real
real number
Definition: types.hpp:50
QL_INTEGER Integer
integer number
Definition: types.hpp:35
std::size_t Size
size of a container
Definition: types.hpp:58
const std::unique_ptr< T > scheme_
Mersenne Twister uniform random number generator.
Definition: any.hpp:35
Real Norm2(const Array &v)
Definition: array.hpp:560
Abstract optimization problem class.
Real alpha
Definition: sabr.cpp:200