| 1 | /* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ |
| 2 | |
| 3 | /* |
| 4 | Copyright (C) 2015 Cheng Li |
| 5 | This file is part of QuantLib, a free-software/open-source library |
| 6 | for financial quantitative analysts and developers - http://quantlib.org/ |
| 7 | |
| 8 | QuantLib is free software: you can redistribute it and/or modify it |
| 9 | under the terms of the QuantLib license. You should have received a |
| 10 | copy of the license along with this program; if not, please email |
| 11 | <quantlib-dev@lists.sf.net>. The license is also available online at |
| 12 | <http://quantlib.org/license.shtml>. |
| 13 | |
| 14 | This program is distributed in the hope that it will be useful, but WITHOUT |
| 15 | ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS |
| 16 | FOR A PARTICULAR PURPOSE. See the license for more details. |
| 17 | */ |
| 18 | |
| 19 | #include <ql/math/optimization/goldstein.hpp> |
| 20 | #include <ql/math/optimization/method.hpp> |
| 21 | #include <ql/math/optimization/problem.hpp> |
| 22 | #include <ql/math/comparison.hpp> |
| 23 | |
| 24 | namespace QuantLib { |
| 25 | |
| 26 | Real GoldsteinLineSearch::operator()(Problem& P, |
| 27 | EndCriteria::Type& ecType, |
| 28 | const EndCriteria& endCriteria, |
| 29 | const Real t_ini) |
| 30 | { |
| 31 | Constraint& constraint = P.constraint(); |
| 32 | succeed_=true; |
| 33 | bool maxIter = false; |
| 34 | Real t = t_ini; |
| 35 | Size loopNumber = 0; |
| 36 | |
| 37 | Real q0 = P.functionValue(); |
| 38 | Real qp0 = P.gradientNormValue(); |
| 39 | |
| 40 | Real tl = 0.0; |
| 41 | Real tr = 0.0; |
| 42 | |
| 43 | qt_ = q0; |
| 44 | qpt_ = (gradient_.empty()) ? qp0 : -DotProduct(v1: gradient_,v2: searchDirection_); |
| 45 | |
| 46 | // Initialize gradient |
| 47 | gradient_ = Array(P.currentValue().size()); |
| 48 | // Compute new point |
| 49 | xtd_ = P.currentValue(); |
| 50 | t = update(params&: xtd_, direction: searchDirection_, beta: t, constraint); |
| 51 | // Compute function value at the new point |
| 52 | qt_ = P.value (x: xtd_); |
| 53 | |
| 54 | while ((qt_ - q0) < -beta_*t*qpt_ || (qt_ - q0) > -alpha_*t*qpt_) { |
| 55 | if ((qt_ - q0) > -alpha_*t*qpt_) |
| 56 | tr = t; |
| 57 | else |
| 58 | tl = t; |
| 59 | ++loopNumber; |
| 60 | |
| 61 | // calculate the new step |
| 62 | if (close_enough(x: tr, y: 0.0)) |
| 63 | t *= extrapolation_; |
| 64 | else |
| 65 | t = (tl + tr) / 2.0; |
| 66 | |
| 67 | // New point value |
| 68 | xtd_ = P.currentValue(); |
| 69 | t = update(params&: xtd_, direction: searchDirection_, beta: t, constraint); |
| 70 | |
| 71 | // Compute function value at the new point |
| 72 | qt_ = P.value (x: xtd_); |
| 73 | P.gradient (grad_f&: gradient_, x: xtd_); |
| 74 | // and it squared norm |
| 75 | maxIter = endCriteria.checkMaxIterations(iteration: loopNumber, ecType); |
| 76 | |
| 77 | if (maxIter) |
| 78 | break; |
| 79 | } |
| 80 | |
| 81 | if (maxIter) |
| 82 | succeed_ = false; |
| 83 | |
| 84 | // Compute new gradient |
| 85 | P.gradient(grad_f&: gradient_, x: xtd_); |
| 86 | // and it squared norm |
| 87 | qpt_ = DotProduct(v1: gradient_, v2: gradient_); |
| 88 | |
| 89 | // Return new step value |
| 90 | return t; |
| 91 | } |
| 92 | |
| 93 | } |
| 94 | |