ROL
ROL_Brents.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_BRENTS_H
45#define ROL_BRENTS_H
46
51#include "ROL_LineSearch.hpp"
52#include "ROL_BackTracking.hpp"
53
54namespace ROL {
55
56template<class Real>
57class Brents : public LineSearch<Real> {
58private:
59 Real tol_;
60 int niter_;
61 bool test_;
62
63 ROL::Ptr<Vector<Real> > xnew_;
64// ROL::Ptr<LineSearch<Real> > btls_;
65
66public:
67
68 virtual ~Brents() {}
69
70 // Constructor
71 Brents( ROL::ParameterList &parlist ) : LineSearch<Real>(parlist) {
72 Real oem10(1.e-10);
73 ROL::ParameterList &list
74 = parlist.sublist("Step").sublist("Line Search").sublist("Line-Search Method").sublist("Brent's");
75 tol_ = list.get("Tolerance",oem10);
76 niter_ = list.get("Iteration Limit",1000);
77 test_ = list.get("Run Test Upon Initialization",true);
78// tol_ = parlist.sublist("Step").sublist("Line Search").sublist("Line-Search Method").get("Bracketing Tolerance",1.e-8);
79// btls_ = ROL::makePtr<BackTracking<Real>>(parlist);
80 }
81
82 void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g,
84 LineSearch<Real>::initialize(x,s,g,obj,con);
85 xnew_ = x.clone();
86// btls_->initialize(x,s,g,obj,con);
87
88 if ( test_ ) {
89 if ( test_brents() ) {
90 std::cout << "Brent's Test Passed!\n";
91 }
92 else {
93 std::cout << "Brent's Test Failed!\n";
94 }
95 }
96 }
97
98 // Find the minimum of phi(alpha) = f(x + alpha*s) using Brent's method
99 void run( Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad,
100 const Real &gs, const Vector<Real> &s, const Vector<Real> &x,
102 ls_neval = 0; ls_ngrad = 0;
103
104 // Get initial line search parameter
105 alpha = LineSearch<Real>::getInitialAlpha(ls_neval,ls_ngrad,fval,gs,x,s,obj,con);
106
107 // TODO: Bracketing
108
109 // Run Brents
110 ROL::Ptr<typename LineSearch<Real>::ScalarFunction> phi
111 = ROL::makePtr<typename LineSearch<Real>::Phi>(*xnew_,x,s,obj,con);
112 int neval = 0;
113 Real A(0), B = alpha;
114 run_brents(neval, fval, alpha, *phi, A, B);
115 ls_neval += neval;
116 }
117
118private:
119 void run_brents(int &neval, Real &fval, Real &alpha,
121 const Real A, const Real B) const {
122 neval = 0;
123 // ---> Set algorithmic constants
124 const Real zero(0), half(0.5), one(1), two(2), three(3), five(5);
125 const Real c = half*(three - std::sqrt(five));
126 const Real eps = std::sqrt(ROL_EPSILON<Real>());
127 // ---> Set end points and initial guess
128 Real a = A, b = B;
129 alpha = a + c*(b-a);
130 // ---> Evaluate function
131 Real fx = phi.value(alpha);
132 neval++;
133 // ---> Initialize algorithm storage
134 Real v = alpha, w = v, u(0), fu(0);
135 Real p(0), q(0), r(0), d(0), e(0);
136 Real fv = fx, fw = fx, tol(0), t2(0), m(0);
137 for (int i = 0; i < niter_; i++) {
138 m = half*(a+b);
139 tol = eps*std::abs(alpha) + tol_; t2 = two*tol;
140 // Check stopping criterion
141 if (std::abs(alpha-m) <= t2 - half*(b-a)) {
142 break;
143 }
144 p = zero; q = zero; r = zero;
145 if ( std::abs(e) > tol ) {
146 // Fit parabola
147 r = (alpha-w)*(fx-fv); q = (alpha-v)*(fx-fw);
148 p = (alpha-v)*q - (alpha-w)*r; q = two*(q-r);
149 if ( q > zero ) {
150 p *= -one;
151 }
152 q = std::abs(q);
153 r = e; e = d;
154 }
155 if ( std::abs(p) < std::abs(half*q*r) && p > q*(a-alpha) && p < q*(b-alpha) ) {
156 // A parabolic interpolation step
157 d = p/q; u = alpha + d;
158 // f must not be evaluated too close to a or b
159 if ( (u - a) < t2 || (b - u) < t2 ) {
160 d = (alpha < m) ? tol : -tol;
161 }
162 }
163 else {
164 // A golden section step
165 e = ((alpha < m) ? b : a) - alpha; d = c*e;
166 }
167 // f must not be evaluated too close to alpha
168 u = alpha + ((std::abs(d) >= tol) ? d : ((d > zero) ? tol : -tol));
169 fu = phi.value(u);
170 neval++;
171 // Update a, b, v, w, and alpha
172 if ( fu <= fx ) {
173 if ( u < alpha ) {
174 b = alpha;
175 }
176 else {
177 a = alpha;
178 }
179 v = w; fv = fw; w = alpha; fw = fx; alpha = u; fx = fu;
180 }
181 else {
182 if ( u < alpha ) {
183 a = u;
184 }
185 else {
186 b = u;
187 }
188 if ( fu <= fw || w == alpha ) {
189 v = w; fv = fw; w = u; fw = fu;
190 }
191 else if ( fu <= fv || v == alpha || v == w ) {
192 v = u; fv = fu;
193 }
194 }
195 }
196 fval = fx;
197 }
198
199 class testFunction : public LineSearch<Real>::ScalarFunction {
200 public:
201 Real value(const Real x) {
202 Real val(0), I(0), two(2), five(5);
203 for (int i = 0; i < 20; i++) {
204 I = (Real)(i+1);
205 val += std::pow((two*I - five)/(x-(I*I)),two);
206 }
207 return val;
208 }
209 };
210
211 bool test_brents(void) const {
212 ROL::Ptr<typename LineSearch<Real>::ScalarFunction> phi
213 = ROL::makePtr<testFunction>();
214 Real A(0), B(0), alpha(0), fval(0);
215 Real error(0), error_i(0);
216 Real zero(0), two(2), three(3);
217 int neval = 0;
218 std::vector<Real> fvector(19,zero), avector(19,zero);
219 fvector[0] = 3.6766990169; avector[0] = 3.0229153;
220 fvector[1] = 1.1118500100; avector[1] = 6.6837536;
221 fvector[2] = 1.2182217637; avector[2] = 11.2387017;
222 fvector[3] = 2.1621103109; avector[3] = 19.6760001;
223 fvector[4] = 3.0322905193; avector[4] = 29.8282273;
224 fvector[5] = 3.7583856477; avector[5] = 41.9061162;
225 fvector[6] = 4.3554103836; avector[6] = 55.9535958;
226 fvector[7] = 4.8482959563; avector[7] = 71.9856656;
227 fvector[8] = 5.2587585400; avector[8] = 90.0088685;
228 fvector[9] = 5.6036524295; avector[9] = 110.0265327;
229 fvector[10] = 5.8956037976; avector[10] = 132.0405517;
230 fvector[11] = 6.1438861542; avector[11] = 156.0521144;
231 fvector[12] = 6.3550764593; avector[12] = 182.0620604;
232 fvector[13] = 6.5333662003; avector[13] = 210.0711010;
233 fvector[14] = 6.6803639849; avector[14] = 240.0800483;
234 fvector[15] = 6.7938538365; avector[15] = 272.0902669;
235 fvector[16] = 6.8634981053; avector[16] = 306.1051233;
236 fvector[17] = 6.8539024631; avector[17] = 342.1369454;
237 fvector[18] = 6.6008470481; avector[18] = 380.2687097;
238 for ( int i = 0; i < 19; i++ ) {
239 A = std::pow((Real)(i+1),two);
240 B = std::pow((Real)(i+2),two);
241 run_brents(neval, fval, alpha, *phi, A, B);
242 error_i = std::max(std::abs(fvector[i]-fval)/fvector[i],
243 std::abs(avector[i]-alpha)/avector[i]);
244 error = std::max(error,error_i);
245 }
246 return (error < three*(std::sqrt(ROL_EPSILON<Real>())*avector[18]+tol_)) ? true : false;
247 }
248
249};
250
251}
252
253#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the interface to apply upper and lower bound constraints.
Real value(const Real x)
Implements a Brent's method line search.
void run_brents(int &neval, Real &fval, Real &alpha, typename LineSearch< Real >::ScalarFunction &phi, const Real A, const Real B) const
Brents(ROL::ParameterList &parlist)
void run(Real &alpha, Real &fval, int &ls_neval, int &ls_ngrad, const Real &gs, const Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con)
ROL::Ptr< Vector< Real > > xnew_
bool test_brents(void) const
virtual ~Brents()
void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &con)
Provides interface for and implements line searches.
virtual void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &con)
virtual Real getInitialAlpha(int &ls_neval, int &ls_ngrad, const Real fval, const Real gs, const Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con)
Provides the interface to evaluate objective functions.
Defines the linear algebra or vector space interface.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.