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