1 //------------------------------------------------------------------------
2 // Copyright (C) 1996:
3 // Scientific computing department
4 // Sandia National Laboratories
5 //------------------------------------------------------------------------
6
7 #ifdef HAVE_CONFIG_H
8 #include "OPT++_config.h"
9 #endif
10
11 #ifdef HAVE_STD
12 #include <cstring>
13 #else
14 #include <string.h>
15 #endif
16
17 #include "OptBCQNewton.h"
18 #include "precisio.h"
19 #include "cblas.h"
20 #include "ioformat.h"
21
22 using NEWMAT::Real;
23 using NEWMAT::FloatingPointPrecision;
24 using NEWMAT::ColumnVector;
25 using NEWMAT::Matrix;
26 using NEWMAT::DiagonalMatrix;
27 using NEWMAT::SymmetricMatrix;
28 using NEWMAT::LowerTriangularMatrix;
29
30 namespace OPTPP {
31
32 static char* class_name = "OptBCQNewton";
33
34 //------------------------------------------------------------------------
35 // BCQNewton functions
36 // initHessian
37 // checkConvg
38 // checkDeriv
39 // initOpt
40 // printStatus
41 // stepTolNorm
42 // updateH
43 // computeSearch
44 // updateConstraints
45 // reset
46 //------------------------------------------------------------------------
initHessian()47 void OptBCQNewton::initHessian()
48 {
49 NLP1* nlp = nlprob();
50 int i,n = nlp->getDim();
51 Hessian.ReSize(n);
52 Hessian = 0.0;
53 for (i=1; i<=n; i++) Hessian(i,i) = 1.0;
54 return;
55 }
56
checkConvg()57 int OptBCQNewton::checkConvg() // Check convergence
58 {
59 NLP1* nlp = nlprob();
60 ColumnVector xc(nlp->getXc());
61 int i, n = nlp->getDim();
62
63 // Test 1. step tolerance
64 double step_tol = tol.getStepTol();
65 double snorm = stepTolNorm();
66 double xnorm = Norm2(xc);
67 double stol = step_tol*max(1.0,xnorm);
68 if (snorm <= stol) {
69 strcpy(mesg,"Step tolerance test passed");
70 *optout << "checkConvg: snorm = " << e(snorm,12,4)
71 << " stol = " << e(stol,12,4) << "\n";
72 return 1;
73 }
74
75 // Test 2. function tolerance
76 double ftol = tol.getFTol();
77 double fvalue = nlp->getF();
78 double rftol = ftol*max(1.0,fabs(fvalue));
79 Real deltaf = fprev - fvalue;
80 if (deltaf <= rftol) {
81 strcpy(mesg,"Function tolerance test passed");
82 *optout << "checkConvg: deltaf = " << e(deltaf,12,4)
83 << " ftol = " << e(ftol,12,4) << "\n";
84 return 2;
85 }
86
87 // Test 3. gradient tolerance
88 ColumnVector grad(nlp->getGrad());
89 double gtol = tol.getGTol();
90 double rgtol = gtol*max(1.0,fabs(fvalue));
91 for (i=1; i<=n; i++) if (work_set(i) == true) grad(i) = 0.0;
92 double gnorm = Norm2(grad);
93 if (gnorm <= rgtol) {
94 strcpy(mesg,"Gradient tolerance test passed");
95 *optout << "checkConvg: gnorm = " << e(gnorm,12,4)
96 << " gtol = " << e(rgtol, 12,4) << "\n";
97 return 3;
98 }
99
100 // Test 4. absolute gradient tolerance
101 if (gnorm <= gtol) {
102 strcpy(mesg,"Gradient tolerance test passed");
103 *optout << "checkConvg: gnorm = " << e(gnorm,12,4)
104 << " gtol = " << e(gtol, 12,4) << "\n";
105 return 4;
106 }
107
108 // Nothing to report
109 return 0;
110 }
111
checkDeriv()112 int OptBCQNewton::checkDeriv() // Check the analytic gradient with FD gradient
113 { return checkAnalyticFDGrad(); }
114
initOpt()115 void OptBCQNewton::initOpt()
116 {
117 OptBCNewtonLike::initOpt();
118 updateConstraints(0);
119 return;
120 }
121
printStatus(char * s)122 void OptBCQNewton::printStatus(char *s) // set Message
123 {
124 NLP1* nlp = nlprob();
125
126 if (debug_) *optout << class_name << "::printStatus: \n";
127
128 *optout << "\n\n========= " << s << " ===========\n\n";
129 *optout << "Optimization method = " << method << "\n";
130 *optout << "Dimension of the problem = " << nlp->getDim() << "\n";
131 *optout << "No. of bound constraints = " << nlp->getDim() << "\n";
132 *optout << "Return code = " << ret_code << " ("
133 << mesg << ")\n";
134 *optout << "No. iterations taken = " << iter_taken << "\n";
135 *optout << "No. function evaluations = " << nlp->getFevals() << "\n";
136 *optout << "No. gradient evaluations = " << nlp->getGevals() << "\n";
137
138 if (debug_) {
139 *optout << "Hessian \n";
140 Print(Hessian);
141 }
142
143 tol.printTol(optout);
144
145 nlp->fPrintState(optout, s);
146 }
147
stepTolNorm()148 real OptBCQNewton::stepTolNorm() const
149 {
150 NLP1* nlp = nlprob();
151 ColumnVector step(sx.AsDiagonal()*(nlp->getXc() - xprev));
152 return Norm2(step);
153 }
154
updateH(SymmetricMatrix & Hk,int k)155 SymmetricMatrix OptBCQNewton::updateH(SymmetricMatrix& Hk, int k)
156 {
157 Real mcheps = FloatingPointPrecision::Epsilon();
158 Real sqrteps = sqrt(mcheps);
159
160 NLP1* nlp = nlprob();
161 int i, nr = nlp->getDim();
162 ColumnVector grad(nr), xc;
163 xc = nlp->getXc();
164 grad = nlp->getGrad();
165
166 DiagonalMatrix D(nr);
167
168 // BFGS formula
169 if (k == 0) {
170 Hessian = 0.0;
171 Real typx, xmax, gnorm;
172
173 // Initialize xmax, typx and D to default values
174 xmax = -1.e30; typx = 1.0; D = 1.0;
175
176 gnorm = Norm2(grad);
177
178 for (i=1; i <= nr; i++) xmax = max(xmax,xc(i));
179 if(xmax != 0.0) typx = xmax;
180 if(gnorm!= 0.0) D = gnorm/typx;
181 if (debug_) {
182 *optout << "updateH: gnorm0 = " << gnorm
183 << "typx = " << typx << "\n";
184 }
185 for (i=1; i <= nr; i++) Hessian(i,i) = D(i);
186 return Hessian;
187 }
188
189 // update the portion of H corresponding to the free variable list only
190
191 ColumnVector yk(nr), sk(nr), Bsk(nr);
192 Matrix Htmp(nr,nr);
193
194 yk = grad - gprev;
195 sk = xc - xprev;
196 for (i=1; i<=nr; i++) {
197 if (work_set(i) == true) {
198 yk(i) = sk(i) = 0.0;
199 //*optout << "fixed variable = " << i << "\n";
200 }
201 }
202
203 Real gts = Dot(gprev,sk);
204 Real yts = Dot(yk,sk);
205
206 Real snorm = Norm2(sk);
207 Real ynorm = Norm2(yk);
208
209 if (debug_) {
210 *optout << "updateH: gts = " << gts
211 << " yts = " << yts << "\n";
212 *optout << "updateH: snorm = " << snorm
213 << " ynorm = " << ynorm << "\n";
214 }
215
216 if (yts <= sqrteps*snorm*ynorm) {
217 if (debug_) {
218 *optout << "updateH: <y,s> = " << e(yts,12,4) << " is too small\n";
219 *optout << "updateH: The BFGS update is skipped\n";
220 }
221 Hessian = Hk; return Hk;
222 }
223
224 ColumnVector res(nr);
225 res = yk - Hk*sk;
226 for (i=1; i<=nr; i++) if (work_set(i) == true) res(i) = 0.0;
227 if (res.NormInfinity() <= sqrteps) {
228 if (debug_) {
229 *optout << "updateH: <y,s> = " << e(yts,12,4) << " is too small\n";
230 *optout << "updateH: The BFGS update is skipped\n";
231 }
232 Hessian = Hk; return Hk;
233 }
234
235 Bsk = Hk*sk;
236 for (i=1; i<=nr; i++) if (work_set(i) == true) Bsk(i) = 0.0;
237 Real sBs = Dot(sk,Bsk);
238 Real etol = 1.e-8;
239
240 if (sBs <= etol*snorm*snorm) {
241 if (debug_) {
242 *optout << "updateH: <y,s> = " << e(yts,12,4) << " is too small\n";
243 *optout << "updateH: The BFGS update is skipped\n";
244 }
245 D = sx.AsDiagonal()*sx.AsDiagonal();
246 Hk = 0;
247 for (i=1; i <= nr; i++) Hk(i,i) = D(i);
248 Hessian = Hk; return Hk;
249 }
250
251 Htmp = - (Bsk * Bsk.t()) / sBs;
252 Htmp = Htmp + (yk * yk.t()) / yts;
253 Htmp = Hk + Htmp;
254 Hk << Htmp;
255 Htmp.Release();
256 ColumnVector Bgk(nr), ggrad(nr);
257 Bgk = Hk*grad;
258 ggrad = grad;
259 for (i=1; i<=nr; i++) if (work_set(i) == true) ggrad(i) = 0.0;
260 Real gBg = Dot(ggrad,Bgk);
261 Real gg = Dot(ggrad,ggrad);
262 Real ckp1= gBg/gg;
263 if (debug_) {
264 *optout << "\nupdateH: after update, k = " << k << "\n";
265 *optout << "updateH: sBs = " << sBs << "\n";
266 *optout << "updateH: ckp1 = " << ckp1 << "\n";
267 }
268 Hessian = Hk;
269 return Hk;
270 }
271
272 //----------------------------------------------------------------------------
273 // Compute the maximum step allowed along the search direction sk
274 // before we hit a constraint
275 //---------------------------------------------------------------------------
computeMaxStep(ColumnVector & sk)276 double OptBCQNewton::computeMaxStep(ColumnVector &sk)
277 {
278 NLP1* nlp = nlprob();
279 int i, n = nlp->getDim();
280 double gamma=FLT_MAX, delta;
281 ColumnVector lower = nlp->getConstraints()->getLower();
282 ColumnVector upper = nlp->getConstraints()->getUpper();
283 ColumnVector xc = nlp->getXc();
284
285 double snorm = Norm2(sk);
286 double feas_tol = 1.e-3;
287
288 for (i=1; i<=n; i++) {
289 if (work_set(i) == false) {
290 if (sk(i) > 0.0e0) {
291 delta = (upper(i)-xc(i)) / sk(i);
292 if (delta <= feas_tol) {
293 if (debug_)
294 *optout << "Hit an upper constraint for variable " << i << "\n";
295 }
296 }
297 else if (sk(i) < 0.0e0) {
298 delta = (lower(i)-xc(i)) / sk(i);
299 if (delta <= feas_tol) {
300 if (debug_)
301 *optout << "Hit a lower constraint for variable " << i << "\n";
302 }
303 }
304 gamma = min(gamma,delta);
305 }
306 }
307 if (debug_)
308 *optout << "computeMaxStep: maximum step allowed = " << gamma*snorm << "\n";
309 return gamma*snorm;
310 }
311
computeSearch(SymmetricMatrix & H)312 ColumnVector OptBCQNewton::computeSearch(SymmetricMatrix &H)
313 {
314 NLP1* nlp = nlprob();
315 int i, j, ncnt=0, *index_array, n = nlp->getDim();
316
317 ColumnVector gg(n), sk2(n), sk(n);
318 SymmetricMatrix H1;
319 LowerTriangularMatrix L;
320
321 // set up index_array to count the number of free variables
322
323 index_array = new int[n+1];
324 for (i=1; i<=n; i++) index_array[i] = 0;
325 for (i=1; i<=n; i++)
326 if (work_set(i) == false) index_array[i] = ++ncnt;
327 if (ncnt != (n-nactive)) {
328 *optout << "Number of fixed and free variables do not correspond. \n";
329 exit(-1);
330 }
331
332 // Form the projected Hessian
333
334 H1.ReSize(ncnt);
335 for (i=1; i<=n; i++) {
336 for (j=1; j<=n; j++)
337 if (index_array[i] != 0 && index_array[j] != 0)
338 H1(index_array[i],index_array[j]) = H(i,j);
339 }
340
341 // Form the projected gradient
342
343 gg.ReSize(ncnt,1);
344 for (i=1; i<=n; i++)
345 if (index_array[i] != 0) gg(index_array[i]) = gprev(i);
346
347 // Solve (H1 * sk2 = - gg) for projected search direction sk2
348
349 L.ReSize(ncnt);
350 sk2.ReSize(ncnt,1);
351 if (ncnt == 1) sk2(1) = - gg(1) / H1(1,1);
352 else {
353 L = MCholesky(H1);
354 sk2 = -(L.t().i()*(L.i()*gg));
355 }
356
357 // Form search direction sk from from projected search direction sk2
358
359 for (i=1; i<=n; i++) sk(i) = 0.0;
360 for (i=1; i<=n; i++)
361 if (index_array[i] != 0) sk(i) = sk2(index_array[i]);
362
363 // Sanitation and return
364
365 delete [] index_array;
366 return sk;
367 }
368
updateConstraints(int step_type)369 int OptBCQNewton::updateConstraints(int step_type)
370 {
371 NLP1* nlp = nlprob();
372 int n = nlp->getDim(), ret_flag=0;
373 int i, j, j2, k, *new_active, actcnt=0, notnew;
374 double reduced_grad_norm, feas_tol=1.0e-12;
375 ColumnVector lower(n), upper(n), xc(n), gg(n);
376
377 // initialization
378
379 lower = nlp->getConstraints()->getLower();
380 upper = nlp->getConstraints()->getUpper();
381 xc = nlp->getXc();
382 new_active = new int[n];
383
384 // Add variables to the working set
385
386 for (i=1; i<=n; i++) {
387 if ((fabs(upper(i)-xc(i))<feas_tol) || (fabs(lower(i)-xc(i))<feas_tol)) {
388 if (work_set(i) == false) {
389 new_active[actcnt++] = i; work_set(i) = true; nactive++;
390 *optout << "OptBCQNewton : variable added to working set : " << i << "\n";
391 }
392 }
393 }
394
395 // Delete variables from the active set
396 // First compute the norm of the reduced gradient
397
398 int jdel=0;
399 double ggdel=0;
400
401 gg = nlp->getGrad();
402 for (i=1; i<=n; i++) if(work_set(i) == true) gg(i) = 0.0;
403 reduced_grad_norm = Norm2(gg);
404 if (m_nconvgd > 0 || step_type < 0) {
405 gg = nlp->getGrad();
406 ret_flag = -1;
407 *optout << "OptBCQNewton : reduced_grad_norm = " << reduced_grad_norm << "\n";
408 for (i=1; i<=n; i++) {
409 notnew = true;
410 for (j=0; j<actcnt; j++) if (new_active[j] == i) notnew = false;
411 if (work_set(i) == true && notnew)
412 if (((fabs(upper(i)-xc(i))<feas_tol) && gg(i)>feas_tol) ||
413 ((fabs(lower(i)-xc(i))<feas_tol) && gg(i)<(-feas_tol))) {
414 if (fabs(gg(i)) > ggdel) {jdel = i; ggdel = fabs(gg(i)); }
415 }
416 }
417 if (jdel != 0) {
418 work_set(jdel) = false; nactive--;
419 *optout << "OptBCQNewton : variable deleted from working set : " << jdel << "\n";
420 ret_flag = 1;
421 }
422 }
423 if (nactive > 0) *optout << "OptbCNewton: Current working set \n";
424 k = 1;
425 for (i=1; i<=nactive; i+=10) {
426 *optout << " ----- variables : ";
427 j2 = min(i*10,nactive);
428 for (j=(i-1)*10+1; j<=j2; j++) {
429 while (work_set(k) == false) k++;
430 *optout << d(k,6); k++;
431 }
432 *optout << "\n ";
433 }
434 return ret_flag;
435 }
436
reset()437 void OptBCQNewton::reset()
438 {
439 NLP1* nlp = nlprob();
440 int n = nlp->getDim();
441 nlp->reset();
442 OptimizeClass::defaultReset(n);
443 nactive = 0;
444 work_set = false;
445 }
446
447 } // namespace OPTPP
448