QuantLib: a free/open-source library for quantitative finance
fully annotated source code - version 1.34
Loading...
Searching...
No Matches
levenbergmarquardt.cpp
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) 2006 Klaus Spanderen
5 Copyright (C) 2015 Peter Caspers
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
24#include <ql/functional.hpp>
25#include <memory>
26
27namespace QuantLib {
28
30 Real xtol,
31 Real gtol,
32 bool useCostFunctionsJacobian)
33 : epsfcn_(epsfcn), xtol_(xtol), gtol_(gtol),
34 useCostFunctionsJacobian_(useCostFunctionsJacobian) {}
35
37 return info_;
38 }
39
41 const EndCriteria& endCriteria) {
43 P.reset();
44 Array x_ = P.currentValue();
45 currentProblem_ = &P;
47 int m = initCostValues_.size();
48 int n = x_.size();
52 }
53 std::unique_ptr<Real[]> xx(new Real[n]);
54 std::copy(x_.begin(), x_.end(), xx.get());
55 std::unique_ptr<Real[]> fvec(new Real[m]);
56 std::unique_ptr<Real[]> diag(new Real[n]);
57 int mode = 1;
58 Real factor = 1;
59 int nprint = 0;
60 int info = 0;
61 int nfev =0;
62 std::unique_ptr<Real[]> fjac(new Real[m*n]);
63 int ldfjac = m;
64 std::unique_ptr<int[]> ipvt(new int[n]);
65 std::unique_ptr<Real[]> qtf(new Real[n]);
66 std::unique_ptr<Real[]> wa1(new Real[n]);
67 std::unique_ptr<Real[]> wa2(new Real[n]);
68 std::unique_ptr<Real[]> wa3(new Real[n]);
69 std::unique_ptr<Real[]> wa4(new Real[m]);
70 // requirements; check here to get more detailed error messages.
71 QL_REQUIRE(n > 0, "no variables given");
72 QL_REQUIRE(m >= n,
73 "less functions (" << m <<
74 ") than available variables (" << n << ")");
75 QL_REQUIRE(endCriteria.functionEpsilon() >= 0.0,
76 "negative f tolerance");
77 QL_REQUIRE(xtol_ >= 0.0, "negative x tolerance");
78 QL_REQUIRE(gtol_ >= 0.0, "negative g tolerance");
79 QL_REQUIRE(endCriteria.maxIterations() > 0,
80 "null number of evaluations");
81
82 // call lmdif to minimize the sum of the squares of m functions
83 // in n variables by the Levenberg-Marquardt algorithm.
84 MINPACK::LmdifCostFunction lmdifCostFunction =
85 [this](const auto m, const auto n, const auto x, const auto fvec, const auto iflag) {
86 this->fcn(m, n, x, fvec, iflag);
87 };
88 MINPACK::LmdifCostFunction lmdifJacFunction =
90 ? [this](const auto m, const auto n, const auto x, const auto fjac, const auto iflag) {
91 this->jacFcn(m, n, x, fjac, iflag);
92 }
94 MINPACK::lmdif(m, n, xx.get(), fvec.get(),
95 endCriteria.functionEpsilon(),
96 xtol_,
97 gtol_,
98 endCriteria.maxIterations(),
99 epsfcn_,
100 diag.get(), mode, factor,
101 nprint, &info, &nfev, fjac.get(),
102 ldfjac, ipvt.get(), qtf.get(),
103 wa1.get(), wa2.get(), wa3.get(), wa4.get(),
104 lmdifCostFunction,
105 lmdifJacFunction);
106 info_ = info;
107 // check requirements & endCriteria evaluation
108 QL_REQUIRE(info != 0, "MINPACK: improper input parameters");
109 //QL_REQUIRE(info != 6, "MINPACK: ftol is too small. no further "
110 // "reduction in the sum of squares "
111 // "is possible.");
113 //QL_REQUIRE(info != 5, "MINPACK: number of calls to fcn has "
114 // "reached or exceeded maxfev.");
115 endCriteria.checkMaxIterations(nfev, ecType);
116 QL_REQUIRE(info != 7, "MINPACK: xtol is too small. no further "
117 "improvement in the approximate "
118 "solution x is possible.");
119 QL_REQUIRE(info != 8, "MINPACK: gtol is too small. fvec is "
120 "orthogonal to the columns of the "
121 "jacobian to machine precision.");
122 // set problem
123 std::copy(xx.get(), xx.get()+n, x_.begin());
126
127 return ecType;
128 }
129
130 void LevenbergMarquardt::fcn(int, int n, Real* x, Real* fvec, int*) {
131 Array xt(n);
132 std::copy(x, x+n, xt.begin());
133 // constraint handling needs some improvement in the future:
134 // starting point should not be close to a constraint violation
135 if (currentProblem_->constraint().test(xt)) {
136 const Array& tmp = currentProblem_->values(xt);
137 std::copy(tmp.begin(), tmp.end(), fvec);
138 } else {
139 std::copy(initCostValues_.begin(), initCostValues_.end(), fvec);
140 }
141 }
142
143 void LevenbergMarquardt::jacFcn(int m, int n, Real* x, Real* fjac, int*) {
144 Array xt(n);
145 std::copy(x, x+n, xt.begin());
146 // constraint handling needs some improvement in the future:
147 // starting point should not be close to a constraint violation
148 if (currentProblem_->constraint().test(xt)) {
149 Matrix tmp(m,n);
151 Matrix tmpT = transpose(tmp);
152 std::copy(tmpT.begin(), tmpT.end(), fjac);
153 } else {
155 std::copy(tmpT.begin(), tmpT.end(), fjac);
156 }
157 }
158
159}
1-D array used in linear algebra.
Definition: array.hpp:52
const_iterator end() const
Definition: array.hpp:511
Size size() const
dimension of the array
Definition: array.hpp:495
const_iterator begin() const
Definition: array.hpp:503
bool test(const Array &p) const
Definition: constraint.hpp:57
virtual Array values(const Array &x) const =0
method to overload to compute the cost function values in x
virtual Real value(const Array &x) const
method to overload to compute the cost function value in x
virtual void jacobian(Matrix &jac, const Array &x) const
method to overload to compute J_f, the jacobian of
Criteria to end optimization process:
Definition: endcriteria.hpp:40
Real functionEpsilon() const
Size maxIterations() const
bool checkMaxIterations(Size iteration, EndCriteria::Type &ecType) const
Definition: endcriteria.cpp:56
void fcn(int m, int n, Real *x, Real *fvec, int *iflag)
virtual Integer getInfo() const
LevenbergMarquardt(Real epsfcn=1.0e-8, Real xtol=1.0e-8, Real gtol=1.0e-8, bool useCostFunctionsJacobian=false)
void jacFcn(int m, int n, Real *x, Real *fjac, int *iflag)
EndCriteria::Type minimize(Problem &P, const EndCriteria &endCriteria) override
minimize the optimization problem P
Matrix used in linear algebra.
Definition: matrix.hpp:41
const_iterator begin() const
Definition: matrix.hpp:327
const_iterator end() const
Definition: matrix.hpp:335
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
Array values(const Array &x)
call cost values computation and increment evaluation counter
Definition: problem.hpp:121
void setFunctionValue(Real functionValue)
Definition: problem.hpp:83
CostFunction & costFunction() const
Cost function.
Definition: problem.hpp:74
void setCurrentValue(const Array &currentValue)
Definition: problem.hpp:76
Abstract constraint class.
#define QL_REQUIRE(condition, message)
throw an error if the given pre-condition is not verified
Definition: errors.hpp:117
Maps function, bind and cref to either the boost or std implementation.
QL_REAL Real
real number
Definition: types.hpp:50
QL_INTEGER Integer
integer number
Definition: types.hpp:35
Levenberg-Marquardt optimization method.
wrapper for MINPACK minimization routine
void lmdif(int m, int n, Real *x, Real *fvec, Real ftol, Real xtol, Real gtol, int maxfev, Real epsfcn, Real *diag, int mode, Real factor, int nprint, int *info, int *nfev, Real *fjac, int ldfjac, int *ipvt, Real *qtf, Real *wa1, Real *wa2, Real *wa3, Real *wa4, const QuantLib::MINPACK::LmdifCostFunction &fcn, const QuantLib::MINPACK::LmdifCostFunction &jacFcn)
Definition: lmdif.cpp:1131
ext::function< void(int, int, Real *, Real *, int *)> LmdifCostFunction
Definition: lmdif.hpp:38
Definition: any.hpp:35
Matrix transpose(const Matrix &m)
Definition: matrix.hpp:700
std::uint64_t x_