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/linesearchbasedmethod.hpp>
22 #include <ql/math/optimization/problem.hpp>
23 #include <ql/math/optimization/linesearch.hpp>
24 #include <ql/math/optimization/armijo.hpp>
25 
26 namespace QuantLib {
27 
LineSearchBasedMethod(const ext::shared_ptr<LineSearch> & lineSearch)28     LineSearchBasedMethod::LineSearchBasedMethod(
29                            const ext::shared_ptr<LineSearch>& lineSearch)
30     : lineSearch_(lineSearch) {
31         if (!lineSearch_)
32            lineSearch_ = ext::shared_ptr<LineSearch>(new ArmijoLineSearch);
33     }
34 
35     EndCriteria::Type
minimize(Problem & P,const EndCriteria & endCriteria)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(prevGradient, x_));
61         P.setGradientNormValue(DotProduct(prevGradient, 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, 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(fnew-fold) /
97                         (std::fabs(fnew) + std::fabs(fold) + QL_EPSILON);
98                 if (fdiff < ftol ||
99                     endCriteria.checkMaxIterations(iterationNumber_, ecType)) {
100                     endCriteria.checkStationaryFunctionValue(0.0, 0.0,
101                         maxStationaryStateIterations_, ecType);
102                     endCriteria.checkMaxIterations(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