Actual source code: morethuente.c
1: /*$Id$*/
3: #include "morethuente.h" /*I "tao_solver.h" I*/
4: #include <stdlib.h>
6: /*
7: The subroutine mcstep is taken from the work of Jorge Nocedal.
8: this is a variant of More' and Thuente's routine.
10: subroutine mcstep
12: the purpose of mcstep is to compute a safeguarded step for
13: a linesearch and to update an interval of uncertainty for
14: a minimizer of the function.
16: the parameter stx contains the step with the least function
17: value. the parameter stp contains the current step. it is
18: assumed that the derivative at stx is negative in the
19: direction of the step. if bracket is set true then a
20: minimizer has been bracketed in an interval of uncertainty
21: with endpoints stx and sty.
23: the subroutine statement is
25: subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
26: stpmin,stpmax,info)
28: where
30: stx, fx, and dx are variables which specify the step,
31: the function, and the derivative at the best step obtained
32: so far. The derivative must be negative in the direction
33: of the step, that is, dx and stp-stx must have opposite
34: signs. On output these parameters are updated appropriately.
36: sty, fy, and dy are variables which specify the step,
37: the function, and the derivative at the other endpoint of
38: the interval of uncertainty. On output these parameters are
39: updated appropriately.
41: stp, fp, and dp are variables which specify the step,
42: the function, and the derivative at the current step.
43: If bracket is set true then on input stp must be
44: between stx and sty. On output stp is set to the new step.
46: bracket is a logical variable which specifies if a minimizer
47: has been bracketed. If the minimizer has not been bracketed
48: then on input bracket must be set false. If the minimizer
49: is bracketed then on output bracket is set true.
51: stpmin and stpmax are input variables which specify lower
52: and upper bounds for the step.
54: info is an integer output variable set as follows:
55: if info = 1,2,3,4,5, then the step has been computed
56: according to one of the five cases below. otherwise
57: info = 0, and this indicates improper input parameters.
59: subprograms called
61: fortran-supplied ... abs,max,min,sqrt
63: argonne national laboratory. minpack project. june 1983
64: jorge j. more', david j. thuente
66: */
70: static int TaoStep_LineSearch(TAO_SOLVER tao,
71: double *stx, double *fx, double *dx,
72: double *sty, double *fy, double *dy,
73: double *stp, double *fp, double *dp)
74: {
75: TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
76: double gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
77: TaoInt bound;
79: TaoFunctionBegin;
81: // Check the input parameters for errors
82: neP->infoc = 0;
83: if (neP->bracket && (*stp <= TaoMin(*stx,*sty) || (*stp >= TaoMax(*stx,*sty)))) SETERRQ(1,"bad stp in bracket");
84: if (*dx * (*stp-*stx) >= 0.0) SETERRQ(1,"dx * (stp-stx) >= 0.0");
85: if (neP->stepmax < neP->stepmin) SETERRQ(1,"stepmax > stepmin");
87: // Determine if the derivatives have opposite sign */
88: sgnd = *dp * (*dx / TaoAbsDouble(*dx));
90: if (*fp > *fx) {
91: // Case 1: a higher function value.
92: // The minimum is bracketed. If the cubic step is closer
93: // to stx than the quadratic step, the cubic step is taken,
94: // else the average of the cubic and quadratic steps is taken.
96: neP->infoc = 1;
97: bound = 1;
98: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
99: s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dx));
100: s = TaoMax(s,TaoAbsDouble(*dp));
101: gamma1 = s*sqrt(pow(theta/s,2.0) - (*dx/s)*(*dp/s));
102: if (*stp < *stx) gamma1 = -gamma1;
103: /* Can p be 0? Check */
104: p = (gamma1 - *dx) + theta;
105: q = ((gamma1 - *dx) + gamma1) + *dp;
106: r = p/q;
107: stpc = *stx + r*(*stp - *stx);
108: stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
110: if (TaoAbsDouble(stpc-*stx) < TaoAbsDouble(stpq-*stx)) {
111: stpf = stpc;
112: }
113: else {
114: stpf = stpc + 0.5*(stpq - stpc);
115: }
116: neP->bracket = 1;
117: }
118: else if (sgnd < 0.0) {
119: // Case 2: A lower function value and derivatives of
120: // opposite sign. The minimum is bracketed. If the cubic
121: // step is closer to stx than the quadratic (secant) step,
122: // the cubic step is taken, else the quadratic step is taken.
124: neP->infoc = 2;
125: bound = 0;
126: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
127: s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dx));
128: s = TaoMax(s,TaoAbsDouble(*dp));
129: gamma1 = s*sqrt(pow(theta/s,2.0) - (*dx/s)*(*dp/s));
130: if (*stp > *stx) gamma1 = -gamma1;
131: p = (gamma1 - *dp) + theta;
132: q = ((gamma1 - *dp) + gamma1) + *dx;
133: r = p/q;
134: stpc = *stp + r*(*stx - *stp);
135: stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
137: if (TaoAbsDouble(stpc-*stp) > TaoAbsDouble(stpq-*stp)) {
138: stpf = stpc;
139: }
140: else {
141: stpf = stpq;
142: }
143: neP->bracket = 1;
144: }
145: else if (TaoAbsDouble(*dp) < TaoAbsDouble(*dx)) {
146: // Case 3: A lower function value, derivatives of the
147: // same sign, and the magnitude of the derivative decreases.
148: // The cubic step is only used if the cubic tends to infinity
149: // in the direction of the step or if the minimum of the cubic
150: // is beyond stp. Otherwise the cubic step is defined to be
151: // either stepmin or stepmax. The quadratic (secant) step is also
152: // computed and if the minimum is bracketed then the the step
153: // closest to stx is taken, else the step farthest away is taken.
155: neP->infoc = 3;
156: bound = 1;
157: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
158: s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dx));
159: s = TaoMax(s,TaoAbsDouble(*dp));
161: // The case gamma1 = 0 only arises if the cubic does not tend
162: // to infinity in the direction of the step.
163: gamma1 = s*sqrt(TaoMax(0.0,pow(theta/s,2.0) - (*dx/s)*(*dp/s)));
164: if (*stp > *stx) gamma1 = -gamma1;
165: p = (gamma1 - *dp) + theta;
166: q = (gamma1 + (*dx - *dp)) + gamma1;
167: r = p/q;
168: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
169: else if (*stp > *stx) stpc = neP->stepmax;
170: else stpc = neP->stepmin;
171: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
173: if (neP->bracket) {
174: if (TaoAbsDouble(*stp-stpc) < TaoAbsDouble(*stp-stpq)) {
175: stpf = stpc;
176: }
177: else {
178: stpf = stpq;
179: }
180: }
181: else {
182: if (TaoAbsDouble(*stp-stpc) > TaoAbsDouble(*stp-stpq)) {
183: stpf = stpc;
184: }
185: else {
186: stpf = stpq;
187: }
188: }
189: }
190: else {
191: // Case 4: A lower function value, derivatives of the
192: // same sign, and the magnitude of the derivative does
193: // not decrease. If the minimum is not bracketed, the step
194: // is either stpmin or stpmax, else the cubic step is taken.
196: neP->infoc = 4;
197: bound = 0;
198: if (neP->bracket) {
199: theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
200: s = TaoMax(TaoAbsDouble(theta),TaoAbsDouble(*dy));
201: s = TaoMax(s,TaoAbsDouble(*dp));
202: gamma1 = s*sqrt(pow(theta/s,2.0) - (*dy/s)*(*dp/s));
203: if (*stp > *sty) gamma1 = -gamma1;
204: p = (gamma1 - *dp) + theta;
205: q = ((gamma1 - *dp) + gamma1) + *dy;
206: r = p/q;
207: stpc = *stp + r*(*sty - *stp);
208: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
210: stpf = stpc;
211: }
212: else if (*stp > *stx) {
213: stpf = neP->stepmax;
214: }
215: else {
216: stpf = neP->stepmin;
217: }
218: }
219:
220: // Update the interval of uncertainty. This update does not
221: // depend on the new step or the case analysis above.
223: if (*fp > *fx) {
224: *sty = *stp;
225: *fy = *fp;
226: *dy = *dp;
227: }
228: else {
229: if (sgnd < 0.0) {
230: *sty = *stx;
231: *fy = *fx;
232: *dy = *dx;
233: }
234: *stx = *stp;
235: *fx = *fp;
236: *dx = *dp;
237: }
238:
239: // Compute the new step and safeguard it.
240: stpf = TaoMin(neP->stepmax,stpf);
241: stpf = TaoMax(neP->stepmin,stpf);
242: *stp = stpf;
243: if (neP->bracket && bound) {
244: if (*sty > *stx) {
245: *stp = TaoMin(*stx+0.66*(*sty-*stx),*stp);
246: }
247: else {
248: *stp = TaoMax(*stx+0.66*(*sty-*stx),*stp);
249: }
250: }
251: TaoFunctionReturn(0);
252: }
256: /* @
257: TaoApply_LineSearch - This routine performs a line search algorithm.
259: Input Parameters:
260: + tao - TAO_SOLVER context
261: . X - current iterate (on output X contains new iterate, X + step*S)
262: . S - search direction
263: . f - objective function evaluated at X
264: . G - gradient evaluated at X
265: . W - work vector
266: - step - initial estimate of step length
268: Output parameters:
269: + f - objective function evaluated at new iterate, X + step*S
270: . G - gradient evaluated at new iterate, X + step*S
271: . X - new iterate
272: . gnorm - 2-norm of G
273: - step - final step length
276: Info is set to one of:
277: + 0 - the line search succeeds; the sufficient decrease
278: condition and the directional derivative condition hold
280: negative number if an input parameter is invalid
281: . -1 - step < 0
282: . -2 - ftol < 0
283: . -3 - rtol < 0
284: . -4 - gtol < 0
285: . -5 - stepmin < 0
286: . -6 - stepmax < stepmin
287: - -7 - maxfev < 0
289: positive number > 1 if the line search otherwise terminates
290: + 2 - Relative width of the interval of uncertainty is
291: at most rtol.
292: . 3 - Maximum number of function evaluations (maxfev) has
293: been reached.
294: . 4 - Step is at the lower bound, stepmin.
295: . 5 - Step is at the upper bound, stepmax.
296: . 6 - Rounding errors may prevent further progress.
297: There may not be a step that satisfies the
298: sufficient decrease and curvature conditions.
299: Tolerances may be too small.
300: - 7 - Search direction is not a descent direction.
302: Notes:
303: This algorithm is taken from More' and Thuente, "Line search algorithms
304: with guaranteed sufficient decrease", Argonne National Laboratory,
305: Technical Report MCS-P330-1092.
307: Notes:
308: This routine is used within the following TAO unconstrained minimization
309: solvers: Newton linesearch (tao_nls), limited memory variable metric
310: (tao_lmvm), and conjugate gradient (tao_cg).
312: Level: advanced
314: .keywords: TAO_SOLVER, linesearch
315: @ */
316: static int TaoApply_LineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f, double *f_full,
317: double *step,TaoInt *info2,void*ctx)
318: {
319: TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
320: double xtrapf = 4.0;
321: double finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
322: double dgx, dgy, dg, fx, fy, stx, sty, dgtest, ftest1=0.0;
323: int info;
324: TaoInt i, stage1;
326: #if defined(PETSC_USE_COMPLEX)
327: PetscScalar cdginit, cdg, cstep = 0.0;
328: #endif
330: TaoFunctionBegin;
331: /* neP->stepmin - lower bound for step */
332: /* neP->stepmax - upper bound for step */
333: /* neP->rtol - relative tolerance for an acceptable step */
334: /* neP->ftol - tolerance for sufficient decrease condition */
335: /* neP->gtol - tolerance for curvature condition */
336: /* neP->nfev - number of function evaluations */
337: /* neP->maxfev - maximum number of function evaluations */
339: *info2 = 0;
341: /* Check input parameters for errors */
342: if (*step < 0.0) {
343: info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: step (%g) < 0\n",*step); CHKERRQ(info);
344: *info2 = -1;
345: }
346:
347: if (neP->ftol < 0.0) {
348: info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: ftol (%g) < 0\n",neP->ftol); CHKERRQ(info);
349: *info2 = -2;
350: }
351:
352: if (neP->rtol < 0.0) {
353: info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: rtol (%g) < 0\n",neP->rtol); CHKERRQ(info);
354: *info2 = -3;
355: }
357: if (neP->gtol < 0.0) {
358: info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: gtol (%g) < 0\n",neP->gtol); CHKERRQ(info);
359: *info2 = -4;
360: }
362: if (neP->stepmin < 0.0) {
363: info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: stepmin (%g) < 0\n",neP->stepmin); CHKERRQ(info);
364: *info2 = -5;
365: }
367: if (neP->stepmax < neP->stepmin) {
368: info = PetscInfo2(tao, "TaoApply_LineSearch: Line search error: stepmax (%g) < stepmin (%g)\n", neP->stepmax,neP->stepmin); CHKERRQ(info);
369: *info2 = -6;
370: }
371:
372: if (neP->maxfev < 0) {
373: info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: maxfev (%d) < 0\n",neP->maxfev); CHKERRQ(info);
374: *info2 = -7;
375: }
377: if (TaoInfOrNaN(*f)) {
378: info = PetscInfo1(tao, "TaoApply_LineSearch: Line search error: function (%g) inf or nan\n", *f); CHKERRQ(info);
379: *info2 = -8;
380: }
382: /* Check that search direction is a descent direction */
384: #if defined(PETSC_USE_COMPLEX)
385: info = G->Dot(S,&cdginit);CHKERRQ(info); dginit = TaoReal(cdginit);
386: #else
387: info = G->Dot(S,&dginit);CHKERRQ(info);
388: #endif
390: if (TaoInfOrNaN(dginit)) {
391: info = PetscInfo1(tao,"TaoApply_LineSearch: Line search error: dginit (%g) inf or nan\n", dginit); CHKERRQ(info);
392: *info2 = -9;
393: }
395: if (dginit >= 0.0) {
396: info = PetscInfo(tao,"TaoApply_LineSearch:Search direction not a descent direction\n"); CHKERRQ(info);
397: *info2 = -10;
398: }
400: if (*info2) {
401: TaoFunctionReturn(0);
402: }
404: /* Initialization */
405: neP->bracket = 0;
406: *info2 = 0;
407: stage1 = 1;
408: finit = *f;
409: dgtest = neP->ftol * dginit;
410: width = neP->stepmax - neP->stepmin;
411: width1 = width * 2.0;
412: info = W->CopyFrom(X);CHKERRQ(info);
413: /* Variable dictionary:
414: stx, fx, dgx - the step, function, and derivative at the best step
415: sty, fy, dgy - the step, function, and derivative at the other endpoint
416: of the interval of uncertainty
417: step, f, dg - the step, function, and derivative at the current step */
419: stx = 0.0;
420: fx = finit;
421: dgx = dginit;
422: sty = 0.0;
423: fy = finit;
424: dgy = dginit;
425:
426: neP->nfev = 0;
427: tao->new_search=TAO_TRUE;
428: for (i=0; i< neP->maxfev; i++) {
429: /* Set min and max steps to correspond to the interval of uncertainty */
430: if (neP->bracket) {
431: neP->stepmin = TaoMin(stx,sty);
432: neP->stepmax = TaoMax(stx,sty);
433: }
434: else {
435: neP->stepmin = stx;
436: neP->stepmax = *step + xtrapf * (*step - stx);
437: }
439: /* Force the step to be within the bounds */
440: *step = TaoMax(*step,neP->stepmin);
441: *step = TaoMin(*step,neP->stepmax);
442:
443: /* If an unusual termination is to occur, then let step be the lowest
444: point obtained thus far */
445: if (((neP->bracket) && (*step <= neP->stepmin || *step >= neP->stepmax)) ||
446: ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol * neP->stepmax)) ||
447: (neP->nfev >= neP->maxfev - 1) || (neP->infoc == 0)) {
448: *step = stx;
449: }
451: #if defined(PETSC_USE_COMPLEX)
452: cstep = *step;
453: info = X->Waxpby(cstep,S,1.0,W);CHKERRQ(info);
454: #else
455: info = X->Waxpby(*step,S,1.0,W);CHKERRQ(info); /* X = W + step*S */
456: #endif
457: tao->current_step=*step;
458: info = TaoComputeMeritFunctionGradient(tao,X,f,G);CHKERRQ(info);
459: tao->new_search=TAO_FALSE;
460: if (0 == i) {
461: *f_full = *f;
462: }
464: neP->nfev++;
465: #if defined(PETSC_USE_COMPLEX)
466: info = G->Dot(S,&cdg);CHKERRQ(info); dg = TaoReal(cdg);
467: #else
468: info = G->Dot(S,&dg);CHKERRQ(info); /* dg = G^T S */
469: #endif
471: if (TaoInfOrNaN(*f) || TaoInfOrNaN(dg)) {
472: // User provided compute function generated Not-a-Number, assume
473: // domain violation and set function value and directional
474: // derivative to infinity.
475: *f = TAO_INFINITY;
476: dg = TAO_INFINITY;
477: }
479: ftest1 = finit + *step * dgtest;
481: /* Convergence testing */
482: if (((*f - ftest1 <= 1.0e-10 * fabs(finit)) &&
483: (TaoAbsDouble(dg) + neP->gtol*dginit <= 0.0))) {
484: info = PetscInfo(tao, "TaoApply_LineSearch:Line search success: Sufficient decrease and directional deriv conditions hold\n"); CHKERRQ(info);
485: *info2 = 0;
486: break;
487: }
489: /* Checks for bad cases */
490: if (((neP->bracket) && (*step <= neP->stepmin||*step >= neP->stepmax)) || (!neP->infoc)) {
491: info = PetscInfo(tao,"TaoApply_LineSearch:Rounding errors may prevent further progress. May not be a step satisfying\n"); CHKERRQ(info);
492: info = PetscInfo(tao,"TaoApply_LineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n"); CHKERRQ(info);
493: *info2 = 6; break;
494: }
495: if ((*step == neP->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
496: info = PetscInfo1(tao,"TaoApply_LineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax); CHKERRQ(info);
497: *info2 = 5; break;
498: }
499: if ((*step == neP->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
500: info = PetscInfo1(tao,"TaoApply_LineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin); CHKERRQ(info);
501: *info2 = 4; break;
502: }
503: if (neP->nfev >= neP->maxfev) {
504: info = PetscInfo2(tao,"TaoApply_LineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev); CHKERRQ(info);
505: *info2 = 3; break;
506: }
507: if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
508: info = PetscInfo1(tao,"TaoApply_LineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol); CHKERRQ(info);
509: *info2 = 2; break;
510: }
512: /* In the first stage, we seek a step for which the modified function
513: has a nonpositive value and nonnegative derivative */
514: if ((stage1) && (*f <= ftest1) && (dg >= dginit * TaoMin(neP->ftol, neP->gtol))) {
515: stage1 = 0;
516: }
518: /* A modified function is used to predict the step only if we
519: have not obtained a step for which the modified function has a
520: nonpositive function value and nonnegative derivative, and if a
521: lower function value has been obtained but the decrease is not
522: sufficient */
524: if ((stage1) && (*f <= fx) && (*f > ftest1)) {
525: fm = *f - *step * dgtest; /* Define modified function */
526: fxm = fx - stx * dgtest; /* and derivatives */
527: fym = fy - sty * dgtest;
528: dgm = dg - dgtest;
529: dgxm = dgx - dgtest;
530: dgym = dgy - dgtest;
531:
532: /* Update the interval of uncertainty and compute the new step */
533: info = TaoStep_LineSearch(tao,&stx,&fxm,&dgxm,&sty,&fym,&dgym,step,&fm,&dgm);CHKERRQ(info);
534:
535: fx = fxm + stx * dgtest; /* Reset the function and */
536: fy = fym + sty * dgtest; /* gradient values */
537: dgx = dgxm + dgtest;
538: dgy = dgym + dgtest;
539: }
540: else {
541: /* Update the interval of uncertainty and compute the new step */
542: info = TaoStep_LineSearch(tao,&stx,&fx,&dgx,&sty,&fy,&dgy,step,f,&dg);CHKERRQ(info);
543: }
544:
545: /* Force a sufficient decrease in the interval of uncertainty */
546: if (neP->bracket) {
547: if (TaoAbsDouble(sty - stx) >= 0.66 * width1) *step = stx + 0.5*(sty - stx);
548: width1 = width;
549: width = TaoAbsDouble(sty - stx);
550: }
551: }
552:
553: /* Finish computations */
554: info = PetscInfo2(tao,"TaoApply_LineSearch:%d function evals in line search, step = %10.4f\n",neP->nfev,*step); CHKERRQ(info);
555: TaoFunctionReturn(0);
556: }
560: /* @
561: TaoApply_BoundLineSearch - This routine performs a line search algorithm.
563: Input Parameters:
564: + tao - TAO_SOLVER context
565: . X - current iterate (on output X contains new iterate, X + step*S)
566: . S - search direction
567: . f - objective function evaluated at X
568: . G - gradient evaluated at X
569: . W - work vector
570: - step - initial estimate of step length
572: Output parameters:
573: + f - objective function evaluated at new iterate, X + step*S
574: . G - gradient evaluated at new iterate, X + step*S
575: . X - new iterate
576: . gnorm - 2-norm of G
577: - step - final step length
580: Info is set to one of:
581: + 0 - the line search succeeds; the sufficient decrease
582: condition and the directional derivative condition hold
584: negative number if an input parameter is invalid
585: . -1 - step < 0
586: . -2 - ftol < 0
587: . -3 - rtol < 0
588: . -4 - gtol < 0
589: . -5 - stepmin < 0
590: . -6 - stepmax < stepmin
591: - -7 - maxfev < 0
593: positive number > 1 if the line search otherwise terminates
594: + 2 - Relative width of the interval of uncertainty is
595: at most rtol.
596: . 3 - Maximum number of function evaluations (maxfev) has
597: been reached.
598: . 4 - Step is at the lower bound, stepmin.
599: . 5 - Step is at the upper bound, stepmax.
600: . 6 - Rounding errors may prevent further progress.
601: There may not be a step that satisfies the
602: sufficient decrease and curvature conditions.
603: Tolerances may be too small.
604: - 7 - Search direction is not a descent direction.
606: Notes:
607: This algorithm is is a modification of the algorithm by More' and
608: Thuente. The modifications concern bounds. This algorithm steps
609: in the direction passed into this routine. This point get
610: projected back into the feasible set. In the context of bound
611: constrained optimization, there may not be a point in the piecewise
612: linear path that satisfies the Wolfe conditions. When the active
613: set is changing, decrease in the objective function may be
614: sufficient to terminate this line search.
616: Note: Much of this coded is identical to the More' Thuente line search.
617: Variations to the code are commented.
619: Notes:
620: This routine is used within the following TAO bound constrained minimization
621: solvers: Newton linesearch (tao_tron) and limited memory variable metric
622: (tao_blmvm).
624: Level: advanced
626: .keywords: TAO_SOLVER, linesearch
627: @ */
628: static int TaoApply_BoundLineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f, double *f_full,
629: double *step,TaoInt *info2,void*ctx)
630: {
631: TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
632: TaoVec *XL,*XU;
633: double xtrapf = 4.0;
634: double finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
635: double dgx, dgy, dg, fx, fy, stx, sty, dgtest, ftest1=0.0, ftest2=0.0;
636: double bstepmin1, bstepmin2, bstepmax;
637: double dg1, dg2;
638: int info;
639: TaoInt i, stage1;
641: #if defined(PETSC_USE_COMPLEX)
642: PetscScalar cdginit, cdg, cstep = 0.0;
643: #endif
645: TaoFunctionBegin;
646: /* neP->stepmin - lower bound for step */
647: /* neP->stepmax - upper bound for step */
648: /* neP->rtol - relative tolerance for an acceptable step */
649: /* neP->ftol - tolerance for sufficient decrease condition */
650: /* neP->gtol - tolerance for curvature condition */
651: /* neP->nfev - number of function evaluations */
652: /* neP->maxfev - maximum number of function evaluations */
654: /* Check input parameters for errors */
655: if (*step < 0.0) {
656: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: step (%g) < 0\n",*step); CHKERRQ(info);
657: *info2 = -1; TaoFunctionReturn(0);
658: }
659: else if (neP->ftol < 0.0) {
660: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: ftol (%g) < 0\n",neP->ftol); CHKERRQ(info);
661: *info2 = -2; TaoFunctionReturn(0);
662: }
663: else if (neP->rtol < 0.0) {
664: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: rtol (%g) < 0\n",neP->rtol); CHKERRQ(info);
665: *info2 = -3; TaoFunctionReturn(0);
666: }
667: else if (neP->gtol < 0.0) {
668: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: gtol (%g) < 0\n",neP->gtol); CHKERRQ(info);
669: *info2 = -4; TaoFunctionReturn(0);
670: }
671: else if (neP->stepmin < 0.0) {
672: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: stepmin (%g) < 0\n",neP->stepmin); CHKERRQ(info);
673: *info2 = -5; TaoFunctionReturn(0);
674: }
675: else if (neP->stepmax < neP->stepmin) {
676: info = PetscInfo2(tao,"TaoApply_BoundLineSearch:Line search error: stepmax (%g) < stepmin (%g)\n", neP->stepmax,neP->stepmin); CHKERRQ(info);
677: *info2 = -6; TaoFunctionReturn(0);
678: }
679: else if (neP->maxfev < 0.0) {
680: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Line search error: maxfev (%d) < 0\n",neP->maxfev); CHKERRQ(info);
681: *info2 = -7; TaoFunctionReturn(0);
682: }
684: /* Compute step length needed to make all variables equal a bound */
685: /* Compute the smallest steplength that will make one nonbinding */
686: /* variable equal the bound */
687: info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
688: info = S->Negate(); CHKERRQ(info);
689: info = S->BoundGradientProjection(S,XL,X,XU); CHKERRQ(info);
690: info = S->Negate(); CHKERRQ(info);
691: info = X->StepBoundInfo(XL,XU,S,&bstepmin1,&bstepmin2,&bstepmax); CHKERRQ(info);
692: neP->stepmax=TaoMin(bstepmax,1.0e+15);
694: /* Check that search direction is a descent direction */
696: #if defined(PETSC_USE_COMPLEX)
697: info = G->Dot(S,&cdginit);CHKERRQ(info); dginit = TaoReal(cdginit);
698: #else
699: info = G->Dot(S,&dginit);CHKERRQ(info);
700: #endif
702: if (dginit >= 0.0) {
703: info = PetscInfo(tao,"TaoApply_BoundLineSearch:Search direction not a descent direction\n"); CHKERRQ(info);
704: *info2 = 7; TaoFunctionReturn(0);
705: }
707: /* Initialization */
708: neP->bracket = 0;
709: *info2 = 0;
710: stage1 = 1;
711: finit = *f;
712: dgtest = neP->ftol * dginit;
713: width = neP->stepmax - neP->stepmin;
714: width1 = width * 2.0;
715: info = W->CopyFrom(X);CHKERRQ(info);
716: /* Variable dictionary:
717: stx, fx, dgx - the step, function, and derivative at the best step
718: sty, fy, dgy - the step, function, and derivative at the other endpoint
719: of the interval of uncertainty
720: step, f, dg - the step, function, and derivative at the current step */
722: stx = 0.0;
723: fx = finit;
724: dgx = dginit;
725: sty = 0.0;
726: fy = finit;
727: dgy = dginit;
728:
729: neP->nfev = 0;
730: tao->new_search=TAO_TRUE;
731: for (i=0; i< neP->maxfev; i++) {
732: /* Set min and max steps to correspond to the interval of uncertainty */
733: if (neP->bracket) {
734: neP->stepmin = TaoMin(stx,sty);
735: neP->stepmax = TaoMax(stx,sty);
736: } else {
737: neP->stepmin = stx;
738: neP->stepmax = *step + xtrapf * (*step - stx);
739: }
741: /* Force the step to be within the bounds */
742: *step = TaoMax(*step,neP->stepmin);
743: *step = TaoMin(*step,neP->stepmax);
745: /* If an unusual termination is to occur, then let step be the lowest
746: point obtained thus far */
747: if (((neP->bracket) && (*step <= neP->stepmin || *step >= neP->stepmax)) ||
748: ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol * neP->stepmax)) ||
749: (neP->nfev >= neP->maxfev - 1) || (neP->infoc == 0)) {
750: *step = stx;
751: }
752:
753: #if defined(PETSC_USE_COMPLEX)
754: cstep = *step;
755: info = X->Waxpby(cstep,S,1.0,W);CHKERRQ(info);
756: #else
757: info = X->Waxpby(*step,S,1.0,W);CHKERRQ(info); /* X = W + step*S */
758: #endif
760: info=X->Median(XL,X,XU);CHKERRQ(info);
761: tao->current_step=*step;
762: info = TaoComputeMeritFunctionGradient(tao,X,f,G);CHKERRQ(info);
763: tao->new_search=TAO_FALSE;
764: if (0 == i) {
765: *f_full = *f;
766: }
768: neP->nfev++;
769: #if defined(PETSC_USE_COMPLEX)
770: info = G->Dot(S,&cdg);CHKERRQ(info); dg = TaoReal(cdg);
771: #else
772: info = G->Dot(W,&dg1);CHKERRQ(info); /* dg = G^T S */
773: info = G->Dot(X,&dg2);CHKERRQ(info); /* dg = G^T S */
774: dg = (dg2-dg1) / (*step);
775: #endif
777: if ((*f != *f) || (dg1 != dg1) || (dg2 != dg2) || (dg != dg)) {
778: // User provided compute function generated Not-a-Number, assume
779: // domain violation and set function value and directional
780: // derivative to infinity.
781: *f = TAO_INFINITY;
782: dg = TAO_INFINITY;
783: dg1 = TAO_INFINITY;
784: dg2 = TAO_INFINITY;
785: }
787: ftest1 = finit + (*step) * dgtest;
788: ftest2 = finit + (*step) * dgtest * neP->ftol; // Armijo
790: /* Convergence testing */
791: if ((*f <= ftest1) && (TaoAbsDouble(dg) <= neP->gtol*(-dginit))) {
792: info = PetscInfo(tao, "TaoApply_BoundLineSearch:Line search success: Sufficient decrease and directional deriv conditions hold\n"); CHKERRQ(info);
793: *info2 = 0;
794: break;
795: }
797: /* Check Armijo if beyond the first breakpoint */
798: if ((*f <= ftest2) && (*step >= bstepmin2)) {
799: info = PetscInfo(tao,"TaoApply_BoundLineSearch:Line search success: Sufficient decrease\n"); CHKERRQ(info);
800: *info2 = 0;
801: break;
802: }
804: /* Checks for bad cases */
805: if (((neP->bracket) && (*step <= neP->stepmin||*step >= neP->stepmax)) || (!neP->infoc)) {
806: info = PetscInfo(tao,"TaoApply_LineSearch:Rounding errors may prevent further progress. May not be a step satisfying\n"); CHKERRQ(info);
807: info = PetscInfo(tao,"TaoApply_BoundLineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n"); CHKERRQ(info);
808: *info2 = 6; break;
809: }
810: if ((*step == neP->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
811: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax); CHKERRQ(info);
812: *info2 = 5; break;
813: }
814: if ((*step == neP->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
815: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin); CHKERRQ(info);
816: *info2 = 4; break;
817: }
818: if (neP->nfev >= neP->maxfev) {
819: info = PetscInfo2(tao,"TaoApply_BoundLineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev); CHKERRQ(info);
820: *info2 = 3; break;
821: }
822: if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
823: info = PetscInfo1(tao,"TaoApply_BoundLineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol); CHKERRQ(info);
824: *info2 = 2; break;
825: }
827: /* In the first stage, we seek a step for which the modified function
828: has a nonpositive value and nonnegative derivative */
829: if ((stage1) && (*f <= ftest1) && (dg >= dginit * TaoMin(neP->ftol, neP->gtol)))
830: stage1 = 0;
832: /* A modified function is used to predict the step only if we
833: have not obtained a step for which the modified function has a
834: nonpositive function value and nonnegative derivative, and if a
835: lower function value has been obtained but the decrease is not
836: sufficient */
838: if ((stage1) && (*f <= fx) && (*f > ftest1)) {
839: fm = *f - *step * dgtest; /* Define modified function */
840: fxm = fx - stx * dgtest; /* and derivatives */
841: fym = fy - sty * dgtest;
842: dgm = dg - dgtest;
843: dgxm = dgx - dgtest;
844: dgym = dgy - dgtest;
846: /* Update the interval of uncertainty and compute the new step */
847: info = TaoStep_LineSearch(tao,&stx,&fxm,&dgxm,&sty,&fym,&dgym,step,&fm,&dgm);CHKERRQ(info);
849: fx = fxm + stx * dgtest; /* Reset the function and */
850: fy = fym + sty * dgtest; /* gradient values */
851: dgx = dgxm + dgtest;
852: dgy = dgym + dgtest;
853: } else {
854: /* Update the interval of uncertainty and compute the new step */
855: info = TaoStep_LineSearch(tao,&stx,&fx,&dgx,&sty,&fy,&dgy,step,f,&dg);CHKERRQ(info);
856: }
858: /* Force a sufficient decrease in the interval of uncertainty */
859: if (neP->bracket) {
860: if (TaoAbsDouble(sty - stx) >= 0.66 * width1) *step = stx + 0.5*(sty - stx);
861: width1 = width;
862: width = TaoAbsDouble(sty - stx);
863: }
864: }
866: /* Finish computations */
867: info = PetscInfo2(tao,"TaoApply_BoundLineSearch:%d function evals in line search, step = %10.4e\n",neP->nfev,*step); CHKERRQ(info);
868: TaoFunctionReturn(0);
869: }
873: static int TaoDestroy_LineSearch(TAO_SOLVER tao, void *ctx)
874: {
875: int info;
877: TaoFunctionBegin;
878: info = TaoFree(tao->linectx);CHKERRQ(info);
879: TaoFunctionReturn(0);
880: }
884: static int TaoSetOptions_LineSearch(TAO_SOLVER tao, void*linectx)
885: {
886: TAO_LINESEARCH *ctx = (TAO_LINESEARCH *)linectx;
887: int info;
889: TaoFunctionBegin;
890: info = TaoOptionsHead("More-Thuente line line search options for minimization");CHKERRQ(info);
891: info = TaoOptionInt("-tao_ls_maxfev","max function evals in line search","",ctx->maxfev,&ctx->maxfev,0);CHKERRQ(info);
892: info = TaoOptionDouble("-tao_ls_ftol","tol for sufficient decrease","",ctx->ftol,&ctx->ftol,0);CHKERRQ(info);
893: info = TaoOptionDouble("-tao_ls_gtol","tol for curvature condition","",ctx->gtol,&ctx->gtol,0);CHKERRQ(info);
894: info = TaoOptionDouble("-tao_ls_rtol","relative tol for acceptable step","",ctx->rtol,&ctx->rtol,0);CHKERRQ(info);
895: info = TaoOptionDouble("-tao_ls_stepmin","lower bound for step","",ctx->stepmin,&ctx->stepmin,0);CHKERRQ(info);
896: info = TaoOptionDouble("-tao_ls_stepmax","upper bound for step","",ctx->stepmax,&ctx->stepmax,0);CHKERRQ(info);
897: info = TaoOptionsTail();CHKERRQ(info);
898: TaoFunctionReturn(0);
899: }
903: static int TaoView_LineSearch(TAO_SOLVER tao, void*ctx)
904: {
905: TAO_LINESEARCH *ls = (TAO_LINESEARCH *)ctx;
906: int info;
908: TaoFunctionBegin;
909: info=TaoPrintInt(tao," More'-Thuente Line search: maxf=%d,",ls->maxfev);CHKERRQ(info);
910: info=TaoPrintDouble(tao," ftol=%g,",ls->ftol);CHKERRQ(info);
911: info=TaoPrintDouble(tao," rtol=%g,",ls->rtol);CHKERRQ(info);
912: info=TaoPrintDouble(tao," gtol=%g\n",ls->gtol);CHKERRQ(info);
913: TaoFunctionReturn(0);
914: }
918: /*@C
919: TaoCreateMoreThuenteLineSearch - Create a line search
921: Input Parameters:
922: + tao - TAO_SOLVER context
923: . fftol - the sufficient descent parameter , greater than 0.
924: - ggtol - the curvature tolerance, greater than 0, less than 1.
927: Note:
928: If either fftol or ggtol is 0, default parameters will be used.
930: Note:
931: This algorithm is taken from More' and Thuente, "Line search algorithms
932: with guaranteed sufficient decrease", Argonne National Laboratory,
933: Technical Report MCS-P330-1092.
935: Note:
936: This line search enforces the strong Wolfe conditions for unconstrained
937: optimization. This routine is used within the following TAO unconstrained
938: minimization solvers: Newton linesearch (tao_nls), limited memory variable
939: metric (tao_lmvm), and nonlinear conjugate gradient methods.
941: Level: developer
943: .keywords: TAO_SOLVER, linesearch
944: @*/
945: int TaoCreateMoreThuenteLineSearch(TAO_SOLVER tao, double fftol, double ggtol)
946: {
947: int info;
948: TAO_LINESEARCH *neP;
950: TaoFunctionBegin;
952: info = TaoNew(TAO_LINESEARCH,&neP); CHKERRQ(info);
953: info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH)); CHKERRQ(info);
955: if (fftol<=0) fftol=0.0001;
956: if (ggtol<=0) ggtol=0.9;
958: neP->ftol = fftol;
959: neP->rtol = 1.0e-10;
960: neP->gtol = ggtol;
961: neP->stepmin = 1.0e-20;
962: neP->stepmax = 1.0e+20;
963: neP->nfev = 0;
964: neP->bracket = 0;
965: neP->infoc = 1;
966: neP->maxfev = 30;
968: info = TaoSetLineSearch(tao,0,
969: TaoSetOptions_LineSearch,
970: TaoApply_LineSearch,
971: TaoView_LineSearch,
972: TaoDestroy_LineSearch,
973: (void *) neP);CHKERRQ(info);
975: TaoFunctionReturn(0);
976: }
980: /*@C
981: TaoCreateMoreThuenteBoundLineSearch - Create a line search
983: Input Parameters:
984: + tao - TAO_SOLVER context
985: . fftol - the sufficient descent parameter , greater than 0.
986: - ggtol - the curvature tolerance, greater than 0, less than 1.
989: Note:
990: If either fftol or ggtol is 0, default parameters will be used.
992: Note:
993: This algorithm is is a modification of the algorithm by More' and Thuente.
994: The modifications concern bounds. This algorithm steps in the direction
995: passed into this routine. This point get projected back into the feasible set.
996: It tries to satisfy the Wolfe conditions, but in the context of bound constrained
997: optimization, there may not be a point in the piecewise linear
998: path that satisfies the Wolfe conditions. When the active set
999: is changing, decrease in the objective function may be sufficient
1000: to terminate this line search.
1002: Note:
1003: This routine is used within the following TAO bound constrained
1004: minimization solvers: Newton trust region (tao_tron) and limited memory variable
1005: metric (tao_blmvm).
1007: Level: developer
1009: .keywords: TAO_SOLVER, linesearch
1010: @*/
1011: int TaoCreateMoreThuenteBoundLineSearch(TAO_SOLVER tao, double fftol, double ggtol)
1012: {
1013: int info;
1014: TAO_LINESEARCH *neP;
1016: TaoFunctionBegin;
1018: info = TaoNew(TAO_LINESEARCH,&neP); CHKERRQ(info);
1019: info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH)); CHKERRQ(info);
1021: if (fftol<=0) fftol=0.0001;
1022: if (ggtol<=0) ggtol=0.9;
1024: neP->ftol = fftol;
1025: neP->rtol = 1.0e-10;
1026: neP->gtol = ggtol;
1027: neP->stepmin = 1.0e-20;
1028: neP->stepmax = 1.0e+20;
1029: neP->nfev = 0;
1030: neP->bracket = 0;
1031: neP->infoc = 1;
1032: neP->maxfev = 30;
1034: info = TaoSetLineSearch(tao,0,
1035: TaoSetOptions_LineSearch,
1036: TaoApply_BoundLineSearch,
1037: TaoView_LineSearch,
1038: TaoDestroy_LineSearch,
1039: (void *) neP);CHKERRQ(info);
1041: TaoFunctionReturn(0);
1042: }