1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2006 Ferdinando Ametrano
5 Copyright (C) 2009 Frédéric Degraeve
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
21#include <ql/math/optimization/armijo.hpp>
22#include <ql/math/optimization/linesearch.hpp>
23#include <ql/math/optimization/linesearchbasedmethod.hpp>
24#include <ql/math/optimization/problem.hpp>
25#include <utility>
26
27namespace QuantLib {
28
29 LineSearchBasedMethod::LineSearchBasedMethod(ext::shared_ptr<LineSearch> lineSearch)
30 : lineSearch_(std::move(lineSearch)) {
31 if (!lineSearch_)
32 lineSearch_ = ext::shared_ptr<LineSearch>(new ArmijoLineSearch);
33 }
34
35 EndCriteria::Type
36 LineSearchBasedMethod::minimize(Problem& P,
37 const EndCriteria& endCriteria) {
38 // Initializations
39 Real ftol = endCriteria.functionEpsilon();
40 Size maxStationaryStateIterations_
41 = endCriteria.maxStationaryStateIterations();
42 EndCriteria::Type ecType = EndCriteria::None; // reset end criteria
43 P.reset(); // reset problem
44 Array x_ = P.currentValue(); // store the starting point
45 Size iterationNumber_ = 0;
46 // dimension line search
47 lineSearch_->searchDirection() = Array(x_.size());
48 bool done = false;
49
50 // function and squared norm of gradient values;
51 Real fnew, fold, gold2;
52 Real fdiff;
53 // classical initial value for line-search step
54 Real t = 1.0;
55 // Set gradient g at the size of the optimization problem
56 // search direction
57 Size sz = lineSearch_->searchDirection().size();
58 Array prevGradient(sz), d(sz), sddiff(sz), direction(sz);
59 // Initialize cost function, gradient prevGradient and search direction
60 P.setFunctionValue(P.valueAndGradient(grad_f&: prevGradient, x: x_));
61 P.setGradientNormValue(DotProduct(v1: prevGradient, v2: prevGradient));
62 lineSearch_->searchDirection() = -prevGradient;
63
64 bool first_time = true;
65 // Loop over iterations
66 do {
67 // Linesearch
68 if (!first_time)
69 prevGradient = lineSearch_->lastGradient();
70 t = (*lineSearch_)(P, ecType, endCriteria, t);
71 // don't throw: it can fail just because maxIterations exceeded
72 //QL_REQUIRE(lineSearch_->succeed(), "line-search failed!");
73 if (lineSearch_->succeed())
74 {
75 // Updates
76
77 // New point
78 x_ = lineSearch_->lastX();
79 // New function value
80 fold = P.functionValue();
81 P.setFunctionValue(lineSearch_->lastFunctionValue());
82 // New gradient and search direction vectors
83
84 // orthogonalization coef
85 gold2 = P.gradientNormValue();
86 P.setGradientNormValue(lineSearch_->lastGradientNorm2());
87
88 // conjugate gradient search direction
89 direction = getUpdatedDirection(P, gold2, gradient: prevGradient);
90
91 sddiff = direction - lineSearch_->searchDirection();
92 lineSearch_->searchDirection() = direction;
93 // Now compute accuracy and check end criteria
94 // Numerical Recipes exit strategy on fx (see NR in C++, p.423)
95 fnew = P.functionValue();
96 fdiff = 2.0*std::fabs(x: fnew-fold) /
97 (std::fabs(x: fnew) + std::fabs(x: fold) + QL_EPSILON);
98 if (fdiff < ftol ||
99 endCriteria.checkMaxIterations(iteration: iterationNumber_, ecType)) {
100 endCriteria.checkStationaryFunctionValue(fxOld: 0.0, fxNew: 0.0,
101 statStateIterations&: maxStationaryStateIterations_, ecType);
102 endCriteria.checkMaxIterations(iteration: iterationNumber_, ecType);
103 return ecType;
104 }
105 P.setCurrentValue(x_); // update problem current value
106 ++iterationNumber_; // Increase iteration number
107 first_time = false;
108 } else {
109 done = true;
110 }
111 } while (!done);
112 P.setCurrentValue(x_);
113 return ecType;
114 }
115
116}
117

source code of quantlib/ql/math/optimization/linesearchbasedmethod.cpp