1 //------------------------------------------------------------------------
2 // Copyright (C) 1993:
3 // J.C. Meza
4 // Sandia National Laboratories
5 // meza@california.sandia.gov
6 //------------------------------------------------------------------------
7
8 #ifdef HAVE_CONFIG_H
9 #include "OPT++_config.h"
10 #endif
11
12 #ifdef HAVE_STD
13 #include <cstring>
14 #include <ctime>
15 #else
16 #include <string.h>
17 #include <time.h>
18 #endif
19
20 #include "OptCG.h"
21 #include "cblas.h"
22 #include "ioformat.h"
23
24 using namespace std;
25
26 using NEWMAT::Real;
27 using NEWMAT::ColumnVector;
28
29 namespace OPTPP {
30
checkDeriv()31 int OptCGLike::checkDeriv() // check the analytic gradient with FD gradient
32 {return GOOD;}
33
checkConvg()34 int OptCGLike::checkConvg() // check convergence
35 {
36 NLP1* nlp = nlprob();
37 ColumnVector xc(nlp->getXc());
38
39 // Test 1. step tolerance
40
41 double step_tol = tol.getStepTol();
42 double snorm = stepTolNorm();
43 double xnorm = Norm2(xc);
44 double stol = step_tol*max(1.0,xnorm);
45 if (snorm <= stol) {
46 strcpy(mesg,"Step tolerance test passed");
47 *optout << "checkConvg: snorm = " << e(snorm,12,4)
48 << " stol = " << e(stol,12,4) << "\n";
49 return 1;
50 }
51
52 // Test 2. function tolerance
53 double ftol = tol.getFTol();
54 double fvalue = nlp->getF();
55 double rftol = ftol*max(1.0,fabs(fvalue));
56 Real deltaf = fprev - fvalue;
57 if (deltaf <= rftol) {
58 strcpy(mesg,"Function tolerance test passed");
59 *optout << "checkConvg: deltaf = " << e(deltaf,12,4)
60 << " ftol = " << e(ftol,12,4) << "\n";
61 return 2;
62 }
63
64
65 // Test 3. gradient tolerance
66
67 ColumnVector grad(nlp->getGrad());
68 double gtol = tol.getGTol();
69 double rgtol = gtol*max(1.0,fabs(fvalue));
70 double gnorm = Norm2(grad);
71 if (gnorm <= rgtol) {
72 strcpy(mesg,"Gradient tolerance test passed");
73 *optout << "checkConvg: gnorm = " << e(gnorm,12,4)
74 << " gtol = " << e(rgtol, 12,4) << "\n";
75 return 3;
76 }
77
78
79 // Test 4. absolute gradient tolerance
80
81 if (gnorm <= gtol) {
82 strcpy(mesg,"Gradient tolerance test passed");
83 *optout << "checkConvg: gnorm = " << e(gnorm,12,4)
84 << " gtol = " << e(gtol, 12,4) << "\n";
85 return 4;
86 }
87
88 // Nothing to report
89
90 return 0;
91
92 }
93
printStatus(char * s)94 void OptCG::printStatus(char *s) // set Message
95 {
96
97 *optout << "\n\n========= " << s << " ===========\n\n";
98 *optout << "Optimization method = " << method << "\n";
99 *optout << "Dimension of the problem = " << dim << "\n";
100 *optout << "Return code = " << ret_code << " ("
101 << mesg << ")\n";
102 *optout << "No. iterations taken = " << iter_taken << "\n";
103 *optout << "No. function evaluations = " << fcn_evals << "\n";
104 *optout << "No. gradient evaluations = " << grad_evals << "\n";
105
106 tol.printTol(optout);
107
108 nlp->fPrintState(optout, s);
109 }
110
reset()111 void OptCG::reset() // Reset parameters
112 {
113 NLP1* nlp = nlprob();
114 int n = nlp->getDim();
115 nlp->reset();
116 OptimizeClass::defaultReset(n);
117 grad_evals = 0;
118 }
119
checkDeriv()120 int OptCG::checkDeriv() // check the analytic gradient with FD gradient
121 {return GOOD;}
122
stepTolNorm()123 real OptCG::stepTolNorm() const
124 {
125 return Norm2(nlp->getXc()-xprev);
126 }
127
computeStep(ColumnVector sk)128 int OptCG::computeStep(ColumnVector sk)
129 //----------------------------------------------------------------------------
130 //
131 // compute a step along the direction sk using either a backtrack line search
132 // or a More-Thuente search
133 //
134 //----------------------------------------------------------------------------
135 {
136 int step_type;
137 int itnmax = tol.getMaxBacktrackIter();
138 real stp_length = 1.0;
139 real stpmax = tol.getMaxStep();
140 real stpmin = tol.getMinStep(); // 1.e-9;
141 real ftol = 5.e-1;
142 real xtol = 2.2e-16;
143 real gtol = 5.e-1;
144
145 step_type = linesearch(nlp, optout, sk, sx, &stp_length, stpmax, stpmin,
146 itnmax, ftol, xtol, gtol);
147 if (step_type < 0) {
148 setMesg("OptCG: Step does not satisfy sufficient decrease condition");
149 ret_code = -1;
150 setReturnCode(ret_code);
151 return(-1);
152 }
153 fcn_evals = nlp->getFevals();
154 grad_evals = nlp->getGevals();
155 step_length = stp_length;
156 return(step_type);
157 }
158
initOpt()159 void OptCG::initOpt()
160 {
161 time_t t;
162 char *c;
163
164 // get date and print out header
165
166 t = time(NULL);
167 c = asctime(localtime(&t));
168
169 *optout << "************************************************************\n";
170 *optout << "OPT++ version " << OPT_GLOBALS::OPT_VERSION << "\n";
171 *optout << "Job run at " << c << "\n";
172 copyright();
173 *optout << "************************************************************\n";
174
175 nlp->initFcn();
176 ret_code = 0;
177
178 if(nlp->hasConstraints()){
179 CompoundConstraint* constraints = nlp->getConstraints();
180 ColumnVector xstart = nlp->getXc();
181 double feas_tol = tol.getCTol();
182 bool feasible = constraints->amIFeasible(xstart, feas_tol);
183 if (!feasible) {
184 *optout << "OptCG WARNING: Initial guess not feasible.\n"
185 << "CG may be unable to make progress." << endl;
186 }
187 //cerr << "Error: OptCG does not support bound, linear, or nonlinear "
188 // << "constraints.\n Please select a different method for "
189 // << "constrained problems." << endl;
190 //abort_handler(-1);
191 }
192
193 if (ret_code == 0) {
194 double fvalue, gnorm;
195
196 int n = nlp->getDim();
197 nlp->eval();
198
199 fvalue = nlp->getF();
200 fprev = fvalue;
201 xprev = nlp->getXc();
202 gprev = nlp->getGrad(); gnorm = Norm2(gprev);
203
204
205 *optout << "\n\t\t\t\tNonlinear CG"
206 << "\n Iter F(x) ||grad|| "
207 << "||step|| beta gtp fcn\n\n"
208 << d(0,5) << " " << e(fvalue,12,4) << " " << e(gnorm,12,4) << endl;
209
210 if (debug_) {
211 nlp->fPrintState(optout, "qnewton: Initial Guess");
212 *optout << "xc, grad, step\n";
213 for(int i=1; i<=n; i++)
214 *optout << d(i,6) << e(xprev(i),24,16) << e(gprev(i),24,16) << "\n";
215 }
216 }
217 }
218
optimize()219 void OptCG::optimize()
220 //------------------------------------------------------------------------
221 // Nonlinear Preconditioned Conjugate Gradient
222 //
223 // Given a nonlinear operator objfcn find the minimizer using a
224 // nonlinear conjugate gradient method
225 // This version uses the Polak-Ribiere formula.
226 // and a line search routine due to More and Thuente as implemented
227 // in the routines mcsrch and mcstep
228 //
229 // Notes: The parameters ftol and gtol should be set so that
230 // 0 < ftol < gtol < 0.5
231 // Default values: ftol = 1.e-1, gtol = 5.e-1
232 // This results in a fairly accurate line search
233 //
234 // Here is the mathematical description of the algorithm
235 // (g = grad f).
236 // -1
237 // 1. set z = M g, search = -z;
238 //
239 // 2. for i=0 until convergence
240 //
241 // find alpha that minimizes f(x + alpha*search)
242 // subject to the strong Wolfe conditions
243 //
244 // Test for convergence
245 //
246 //
247 // beta = -( g , (z - z ) ) / ( g , z )
248 // i+1 i + 1 i i i
249 //
250 // search = - z + beta * search
251 // i+1 i+1 i
252 //
253 //----------------------------------------------------------------------------
254
255 {
256
257 int i, nlcg_iter;
258 int convgd = 0;
259 int step_type;
260
261 double beta;
262 double delta, delta_old, delta_mid, delta_new;
263 double slope, gnorm;
264
265 double step;
266 double zero = 0.;
267
268 // Allocate local vectors
269
270 int n = dim;
271 int maxiter;
272 double fvalue;
273 ColumnVector search(n), grad(n), z(n), diag(n), xc(n);
274
275 // Initialize iteration
276
277 maxiter = tol.getMaxIter();
278
279 initOpt();
280
281 if (ret_code == 0) {
282 // compute preconditioned gradient
283
284 diag = getFcnScale();
285 grad = nlp->getGrad();
286 for (i=1; i<=n; i++) z(i) = grad(i)/diag(i);
287
288 search = -z;
289 delta_old = delta_new = Dot(grad,z);
290 gnorm = sqrt(Dot(grad,grad));
291
292 step = 1.0/gnorm;
293
294 //---------------------------------------------------------------------------
295 //
296 //
297 //
298 //
299 //---------------------------------------------------------------------------
300
301 for (nlcg_iter=1; nlcg_iter <= maxiter; nlcg_iter++) {
302
303 iter_taken = nlcg_iter;
304
305 // compute a step along the direction search
306
307 if ((step_type = computeStep(search)) < 0) {
308 setMesg("nlcg: Step does not satisfy sufficient decrease condition");
309 ret_code = step_type;
310 setReturnCode(ret_code);
311 return;
312 }
313
314 // Accept this step and update the nonlinear model
315
316 acceptStep(nlcg_iter, step_type);
317 updateModel(nlcg_iter, n, xprev);
318
319 xc = nlp->getXc();
320 mem_step = xc - xprev;
321 step = Norm2(mem_step);
322
323 fvalue = nlp->getF();
324 grad = nlp->getGrad();
325 gnorm = sqrt(Dot(grad,grad));
326 slope = Dot(grad,search);
327
328 // Test for Convergence
329
330 convgd = checkConvg();
331 if (convgd > 0) {
332 ret_code = convgd;
333 setReturnCode(ret_code);
334 setMesg("OptCG: Algorithm converged");
335 *optout << d(nlcg_iter,5) << " " << e(fvalue,12,4) << " "
336 << e(gnorm,12,4) << e(step,12,4) << "\n";
337 return;
338 }
339
340 //
341 // compute a new search direction
342 // 1. compute preconditioned gradient, z = grad;
343 // 2. beta is computed using Polak-Ribiere Formula constrained
344 // so that beta > 0
345 // 3 Update search direction and norms
346
347 delta_old = delta_new; delta_mid = Dot(grad,z);
348
349 for (i=1; i<=n; i++) z(i) = grad(i)/diag(i);
350
351 delta_new = Dot(grad,z);
352 delta = delta_new - delta_mid;
353 beta = max(zero,delta/delta_old);
354
355 search = -z + search*beta;
356
357 xprev = nlp->getXc();
358 fprev = fvalue;
359 gprev = grad;
360
361 *optout
362 << d(nlcg_iter,5) << " " << e(fvalue,12,4) << " " << e(gnorm,12,4)
363 << e(step,12,4) << " " << e(beta,12,4) << " " << e(slope,12,4)
364 << d(fcn_evals,4) << " " << d(grad_evals,4) << endl;
365 }
366
367 setMesg("Maximum number of iterations in nlcg");
368 ret_code = 4;
369 setReturnCode(ret_code);
370 }
371 }
372
373 } // namespace OPTPP
374