Actual source code: morethuente.c
1: #include <petsc/private/taolinesearchimpl.h>
2: #include <../src/tao/linesearch/impls/morethuente/morethuente.h>
4: /*
5: This algorithm is taken from More' and Thuente, "Line search algorithms
6: with guaranteed sufficient decrease", Argonne National Laboratory,
7: Technical Report MCS-P330-1092.
8: */
10: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp);
12: static PetscErrorCode TaoLineSearchDestroy_MT(TaoLineSearch ls)
13: {
14: TaoLineSearch_MT *mt = (TaoLineSearch_MT*)(ls->data);
16: PetscObjectDereference((PetscObject)mt->x);
17: VecDestroy(&mt->work);
18: PetscFree(ls->data);
19: return 0;
20: }
22: static PetscErrorCode TaoLineSearchSetFromOptions_MT(PetscOptionItems *PetscOptionsObject,TaoLineSearch ls)
23: {
24: return 0;
25: }
27: static PetscErrorCode TaoLineSearchMonitor_MT(TaoLineSearch ls)
28: {
29: TaoLineSearch_MT *mt = (TaoLineSearch_MT*)ls->data;
31: PetscViewerASCIIPrintf(ls->viewer, "stx: %g, fx: %g, dgx: %g\n", (double)mt->stx, (double)mt->fx, (double)mt->dgx);
32: PetscViewerASCIIPrintf(ls->viewer, "sty: %g, fy: %g, dgy: %g\n", (double)mt->sty, (double)mt->fy, (double)mt->dgy);
33: return 0;
34: }
36: static PetscErrorCode TaoLineSearchApply_MT(TaoLineSearch ls, Vec x, PetscReal *f, Vec g, Vec s)
37: {
38: TaoLineSearch_MT *mt = (TaoLineSearch_MT*)(ls->data);
39: PetscReal xtrapf = 4.0;
40: PetscReal finit, width, width1, dginit, fm, fxm, fym, dgm, dgxm, dgym;
41: PetscReal dgx, dgy, dg, dg2, fx, fy, stx, sty, dgtest;
42: PetscReal ftest1=0.0, ftest2=0.0;
43: PetscInt i, stage1,n1,n2,nn1,nn2;
44: PetscReal bstepmin1, bstepmin2, bstepmax;
45: PetscBool g_computed = PETSC_FALSE; /* to prevent extra gradient computation */
47: ls->reason = TAOLINESEARCH_CONTINUE_ITERATING;
48: TaoLineSearchMonitor(ls, 0, *f, 0.0);
49: /* Check work vector */
50: if (!mt->work) {
51: VecDuplicate(x,&mt->work);
52: mt->x = x;
53: PetscObjectReference((PetscObject)mt->x);
54: } else if (x != mt->x) {
55: VecDestroy(&mt->work);
56: VecDuplicate(x,&mt->work);
57: PetscObjectDereference((PetscObject)mt->x);
58: mt->x = x;
59: PetscObjectReference((PetscObject)mt->x);
60: }
62: if (ls->bounded) {
63: /* Compute step length needed to make all variables equal a bound */
64: /* Compute the smallest steplength that will make one nonbinding variable
65: equal the bound */
66: VecGetLocalSize(ls->upper,&n1);
67: VecGetLocalSize(mt->x, &n2);
68: VecGetSize(ls->upper,&nn1);
69: VecGetSize(mt->x,&nn2);
71: VecScale(s,-1.0);
72: VecBoundGradientProjection(s,x,ls->lower,ls->upper,s);
73: VecScale(s,-1.0);
74: VecStepBoundInfo(x,s,ls->lower,ls->upper,&bstepmin1,&bstepmin2,&bstepmax);
75: ls->stepmax = PetscMin(bstepmax,1.0e15);
76: }
78: VecDot(g,s,&dginit);
79: if (PetscIsInfOrNanReal(dginit)) {
80: PetscInfo(ls,"Initial Line Search step * g is Inf or Nan (%g)\n",(double)dginit);
81: ls->reason = TAOLINESEARCH_FAILED_INFORNAN;
82: return 0;
83: }
84: if (dginit >= 0.0) {
85: PetscInfo(ls,"Initial Line Search step * g is not descent direction (%g)\n",(double)dginit);
86: ls->reason = TAOLINESEARCH_FAILED_ASCENT;
87: return 0;
88: }
90: /* Initialization */
91: mt->bracket = 0;
92: stage1 = 1;
93: finit = *f;
94: dgtest = ls->ftol * dginit;
95: width = ls->stepmax - ls->stepmin;
96: width1 = width * 2.0;
97: VecCopy(x,mt->work);
98: /* Variable dictionary:
99: stx, fx, dgx - the step, function, and derivative at the best step
100: sty, fy, dgy - the step, function, and derivative at the other endpoint
101: of the interval of uncertainty
102: step, f, dg - the step, function, and derivative at the current step */
104: stx = 0.0;
105: fx = finit;
106: dgx = dginit;
107: sty = 0.0;
108: fy = finit;
109: dgy = dginit;
111: ls->step = ls->initstep;
112: for (i=0; i< ls->max_funcs; i++) {
113: /* Set min and max steps to correspond to the interval of uncertainty */
114: if (mt->bracket) {
115: ls->stepmin = PetscMin(stx,sty);
116: ls->stepmax = PetscMax(stx,sty);
117: } else {
118: ls->stepmin = stx;
119: ls->stepmax = ls->step + xtrapf * (ls->step - stx);
120: }
122: /* Force the step to be within the bounds */
123: ls->step = PetscMax(ls->step,ls->stepmin);
124: ls->step = PetscMin(ls->step,ls->stepmax);
126: /* If an unusual termination is to occur, then let step be the lowest
127: point obtained thus far */
128: if ((stx!=0) && (((mt->bracket) && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) ||
129: ((ls->nfeval+ls->nfgeval) >= ls->max_funcs - 1) || (mt->infoc == 0))) {
130: ls->step = stx;
131: }
133: VecCopy(x,mt->work);
134: VecAXPY(mt->work,ls->step,s); /* W = X + step*S */
136: if (ls->bounded) {
137: VecMedian(ls->lower, mt->work, ls->upper, mt->work);
138: }
139: if (ls->usegts) {
140: TaoLineSearchComputeObjectiveAndGTS(ls,mt->work,f,&dg);
141: g_computed = PETSC_FALSE;
142: } else {
143: TaoLineSearchComputeObjectiveAndGradient(ls,mt->work,f,g);
144: g_computed = PETSC_TRUE;
145: if (ls->bounded) {
146: VecDot(g,x,&dg);
147: VecDot(g,mt->work,&dg2);
148: dg = (dg2 - dg)/ls->step;
149: } else {
150: VecDot(g,s,&dg);
151: }
152: }
154: /* update bracketing parameters in the MT context for printouts in monitor */
155: mt->stx = stx;
156: mt->fx = fx;
157: mt->dgx = dgx;
158: mt->sty = sty;
159: mt->fy = fy;
160: mt->dgy = dgy;
161: TaoLineSearchMonitor(ls, i+1, *f, ls->step);
163: if (i == 0) ls->f_fullstep=*f;
165: if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
166: /* User provided compute function generated Not-a-Number, assume
167: domain violation and set function value and directional
168: derivative to infinity. */
169: *f = PETSC_INFINITY;
170: dg = PETSC_INFINITY;
171: }
173: ftest1 = finit + ls->step * dgtest;
174: if (ls->bounded) ftest2 = finit + ls->step * dgtest * ls->ftol;
176: /* Convergence testing */
177: if (((*f - ftest1 <= 1.0e-10 * PetscAbsReal(finit)) && (PetscAbsReal(dg) + ls->gtol*dginit <= 0.0))) {
178: PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n");
179: ls->reason = TAOLINESEARCH_SUCCESS;
180: break;
181: }
183: /* Check Armijo if beyond the first breakpoint */
184: if (ls->bounded && (*f <= ftest2) && (ls->step >= bstepmin2)) {
185: PetscInfo(ls,"Line search success: Sufficient decrease.\n");
186: ls->reason = TAOLINESEARCH_SUCCESS;
187: break;
188: }
190: /* Checks for bad cases */
191: if (((mt->bracket) && (ls->step <= ls->stepmin||ls->step >= ls->stepmax)) || (!mt->infoc)) {
192: PetscInfo(ls,"Rounding errors may prevent further progress. May not be a step satisfying\n");
193: PetscInfo(ls,"sufficient decrease and curvature conditions. Tolerances may be too small.\n");
194: ls->reason = TAOLINESEARCH_HALTED_OTHER;
195: break;
196: }
197: if ((ls->step == ls->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
198: PetscInfo(ls,"Step is at the upper bound, stepmax (%g)\n",(double)ls->stepmax);
199: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
200: break;
201: }
202: if ((ls->step == ls->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
203: PetscInfo(ls,"Step is at the lower bound, stepmin (%g)\n",(double)ls->stepmin);
204: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
205: break;
206: }
207: if ((mt->bracket) && (ls->stepmax - ls->stepmin <= ls->rtol*ls->stepmax)) {
208: PetscInfo(ls,"Relative width of interval of uncertainty is at most rtol (%g)\n",(double)ls->rtol);
209: ls->reason = TAOLINESEARCH_HALTED_RTOL;
210: break;
211: }
213: /* In the first stage, we seek a step for which the modified function
214: has a nonpositive value and nonnegative derivative */
215: if ((stage1) && (*f <= ftest1) && (dg >= dginit * PetscMin(ls->ftol, ls->gtol))) stage1 = 0;
217: /* A modified function is used to predict the step only if we
218: have not obtained a step for which the modified function has a
219: nonpositive function value and nonnegative derivative, and if a
220: lower function value has been obtained but the decrease is not
221: sufficient */
223: if ((stage1) && (*f <= fx) && (*f > ftest1)) {
224: fm = *f - ls->step * dgtest; /* Define modified function */
225: fxm = fx - stx * dgtest; /* and derivatives */
226: fym = fy - sty * dgtest;
227: dgm = dg - dgtest;
228: dgxm = dgx - dgtest;
229: dgym = dgy - dgtest;
231: /* if (dgxm * (ls->step - stx) >= 0.0) */
232: /* Update the interval of uncertainty and compute the new step */
233: Tao_mcstep(ls,&stx,&fxm,&dgxm,&sty,&fym,&dgym,&ls->step,&fm,&dgm);
235: fx = fxm + stx * dgtest; /* Reset the function and */
236: fy = fym + sty * dgtest; /* gradient values */
237: dgx = dgxm + dgtest;
238: dgy = dgym + dgtest;
239: } else {
240: /* Update the interval of uncertainty and compute the new step */
241: Tao_mcstep(ls,&stx,&fx,&dgx,&sty,&fy,&dgy,&ls->step,f,&dg);
242: }
244: /* Force a sufficient decrease in the interval of uncertainty */
245: if (mt->bracket) {
246: if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5*(sty - stx);
247: width1 = width;
248: width = PetscAbsReal(sty - stx);
249: }
250: }
251: if ((ls->nfeval+ls->nfgeval) > ls->max_funcs) {
252: PetscInfo(ls,"Number of line search function evals (%D) > maximum (%D)\n",(ls->nfeval+ls->nfgeval),ls->max_funcs);
253: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
254: }
256: /* Finish computations */
257: PetscInfo(ls,"%D function evals in line search, step = %g\n",(ls->nfeval+ls->nfgeval),(double)ls->step);
259: /* Set new solution vector and compute gradient if needed */
260: VecCopy(mt->work,x);
261: if (!g_computed) {
262: TaoLineSearchComputeGradient(ls,mt->work,g);
263: }
264: return 0;
265: }
267: /*MC
268: TAOLINESEARCHMT - Line-search type with cubic interpolation that satisfies both the sufficient decrease and
269: curvature conditions. This method can take step lengths greater than 1.
271: More-Thuente line-search can be selected with "-tao_ls_type more-thuente".
273: References:
274: . * - JORGE J. MORE AND DAVID J. THUENTE, LINE SEARCH ALGORITHMS WITH GUARANTEED SUFFICIENT DECREASE.
275: ACM Trans. Math. Software 20, no. 3 (1994): 286-307.
277: Level: developer
279: .seealso: TaoLineSearchCreate(), TaoLineSearchSetType(), TaoLineSearchApply()
281: .keywords: Tao, linesearch
282: M*/
283: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)
284: {
285: TaoLineSearch_MT *ctx;
288: PetscNewLog(ls,&ctx);
289: ctx->bracket = 0;
290: ctx->infoc = 1;
291: ls->data = (void*)ctx;
292: ls->initstep = 1.0;
293: ls->ops->setup = NULL;
294: ls->ops->reset = NULL;
295: ls->ops->apply = TaoLineSearchApply_MT;
296: ls->ops->destroy = TaoLineSearchDestroy_MT;
297: ls->ops->setfromoptions = TaoLineSearchSetFromOptions_MT;
298: ls->ops->monitor = TaoLineSearchMonitor_MT;
299: return 0;
300: }
302: /*
303: The subroutine mcstep is taken from the work of Jorge Nocedal.
304: this is a variant of More' and Thuente's routine.
306: subroutine mcstep
308: the purpose of mcstep is to compute a safeguarded step for
309: a linesearch and to update an interval of uncertainty for
310: a minimizer of the function.
312: the parameter stx contains the step with the least function
313: value. the parameter stp contains the current step. it is
314: assumed that the derivative at stx is negative in the
315: direction of the step. if bracket is set true then a
316: minimizer has been bracketed in an interval of uncertainty
317: with endpoints stx and sty.
319: the subroutine statement is
321: subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
322: stpmin,stpmax,info)
324: where
326: stx, fx, and dx are variables which specify the step,
327: the function, and the derivative at the best step obtained
328: so far. The derivative must be negative in the direction
329: of the step, that is, dx and stp-stx must have opposite
330: signs. On output these parameters are updated appropriately.
332: sty, fy, and dy are variables which specify the step,
333: the function, and the derivative at the other endpoint of
334: the interval of uncertainty. On output these parameters are
335: updated appropriately.
337: stp, fp, and dp are variables which specify the step,
338: the function, and the derivative at the current step.
339: If bracket is set true then on input stp must be
340: between stx and sty. On output stp is set to the new step.
342: bracket is a logical variable which specifies if a minimizer
343: has been bracketed. If the minimizer has not been bracketed
344: then on input bracket must be set false. If the minimizer
345: is bracketed then on output bracket is set true.
347: stpmin and stpmax are input variables which specify lower
348: and upper bounds for the step.
350: info is an integer output variable set as follows:
351: if info = 1,2,3,4,5, then the step has been computed
352: according to one of the five cases below. otherwise
353: info = 0, and this indicates improper input parameters.
355: subprograms called
357: fortran-supplied ... abs,max,min,sqrt
359: argonne national laboratory. minpack project. june 1983
360: jorge j. more', david j. thuente
362: */
364: static PetscErrorCode Tao_mcstep(TaoLineSearch ls,PetscReal *stx,PetscReal *fx,PetscReal *dx,PetscReal *sty,PetscReal *fy,PetscReal *dy,PetscReal *stp,PetscReal *fp,PetscReal *dp)
365: {
366: TaoLineSearch_MT *mtP = (TaoLineSearch_MT *) ls->data;
367: PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
368: PetscInt bound;
370: /* Check the input parameters for errors */
371: mtP->infoc = 0;
376: /* Determine if the derivatives have opposite sign */
377: sgnd = *dp * (*dx / PetscAbsReal(*dx));
379: if (*fp > *fx) {
380: /* Case 1: a higher function value.
381: The minimum is bracketed. If the cubic step is closer
382: to stx than the quadratic step, the cubic step is taken,
383: else the average of the cubic and quadratic steps is taken. */
385: mtP->infoc = 1;
386: bound = 1;
387: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
388: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
389: s = PetscMax(s,PetscAbsReal(*dp));
390: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
391: if (*stp < *stx) gamma1 = -gamma1;
392: /* Can p be 0? Check */
393: p = (gamma1 - *dx) + theta;
394: q = ((gamma1 - *dx) + gamma1) + *dp;
395: r = p/q;
396: stpc = *stx + r*(*stp - *stx);
397: stpq = *stx + ((*dx/((*fx-*fp)/(*stp-*stx)+*dx))*0.5) * (*stp - *stx);
399: if (PetscAbsReal(stpc-*stx) < PetscAbsReal(stpq-*stx)) {
400: stpf = stpc;
401: } else {
402: stpf = stpc + 0.5*(stpq - stpc);
403: }
404: mtP->bracket = 1;
405: } else if (sgnd < 0.0) {
406: /* Case 2: A lower function value and derivatives of
407: opposite sign. The minimum is bracketed. If the cubic
408: step is closer to stx than the quadratic (secant) step,
409: the cubic step is taken, else the quadratic step is taken. */
411: mtP->infoc = 2;
412: bound = 0;
413: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
414: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
415: s = PetscMax(s,PetscAbsReal(*dp));
416: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s));
417: if (*stp > *stx) gamma1 = -gamma1;
418: p = (gamma1 - *dp) + theta;
419: q = ((gamma1 - *dp) + gamma1) + *dx;
420: r = p/q;
421: stpc = *stp + r*(*stx - *stp);
422: stpq = *stp + (*dp/(*dp-*dx))*(*stx - *stp);
424: if (PetscAbsReal(stpc-*stp) > PetscAbsReal(stpq-*stp)) {
425: stpf = stpc;
426: } else {
427: stpf = stpq;
428: }
429: mtP->bracket = 1;
430: } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
431: /* Case 3: A lower function value, derivatives of the
432: same sign, and the magnitude of the derivative decreases.
433: The cubic step is only used if the cubic tends to infinity
434: in the direction of the step or if the minimum of the cubic
435: is beyond stp. Otherwise the cubic step is defined to be
436: either stepmin or stepmax. The quadratic (secant) step is also
437: computed and if the minimum is bracketed then the step
438: closest to stx is taken, else the step farthest away is taken. */
440: mtP->infoc = 3;
441: bound = 1;
442: theta = 3*(*fx - *fp)/(*stp - *stx) + *dx + *dp;
443: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dx));
444: s = PetscMax(s,PetscAbsReal(*dp));
446: /* The case gamma1 = 0 only arises if the cubic does not tend
447: to infinity in the direction of the step. */
448: gamma1 = s*PetscSqrtScalar(PetscMax(0.0,PetscPowScalar(theta/s,2.0) - (*dx/s)*(*dp/s)));
449: if (*stp > *stx) gamma1 = -gamma1;
450: p = (gamma1 - *dp) + theta;
451: q = (gamma1 + (*dx - *dp)) + gamma1;
452: r = p/q;
453: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r*(*stx - *stp);
454: else if (*stp > *stx) stpc = ls->stepmax;
455: else stpc = ls->stepmin;
456: stpq = *stp + (*dp/(*dp-*dx)) * (*stx - *stp);
458: if (mtP->bracket) {
459: if (PetscAbsReal(*stp-stpc) < PetscAbsReal(*stp-stpq)) {
460: stpf = stpc;
461: } else {
462: stpf = stpq;
463: }
464: } else {
465: if (PetscAbsReal(*stp-stpc) > PetscAbsReal(*stp-stpq)) {
466: stpf = stpc;
467: } else {
468: stpf = stpq;
469: }
470: }
471: } else {
472: /* Case 4: A lower function value, derivatives of the
473: same sign, and the magnitude of the derivative does
474: not decrease. If the minimum is not bracketed, the step
475: is either stpmin or stpmax, else the cubic step is taken. */
477: mtP->infoc = 4;
478: bound = 0;
479: if (mtP->bracket) {
480: theta = 3*(*fp - *fy)/(*sty - *stp) + *dy + *dp;
481: s = PetscMax(PetscAbsReal(theta),PetscAbsReal(*dy));
482: s = PetscMax(s,PetscAbsReal(*dp));
483: gamma1 = s*PetscSqrtScalar(PetscPowScalar(theta/s,2.0) - (*dy/s)*(*dp/s));
484: if (*stp > *sty) gamma1 = -gamma1;
485: p = (gamma1 - *dp) + theta;
486: q = ((gamma1 - *dp) + gamma1) + *dy;
487: r = p/q;
488: stpc = *stp + r*(*sty - *stp);
489: stpf = stpc;
490: } else if (*stp > *stx) {
491: stpf = ls->stepmax;
492: } else {
493: stpf = ls->stepmin;
494: }
495: }
497: /* Update the interval of uncertainty. This update does not
498: depend on the new step or the case analysis above. */
500: if (*fp > *fx) {
501: *sty = *stp;
502: *fy = *fp;
503: *dy = *dp;
504: } else {
505: if (sgnd < 0.0) {
506: *sty = *stx;
507: *fy = *fx;
508: *dy = *dx;
509: }
510: *stx = *stp;
511: *fx = *fp;
512: *dx = *dp;
513: }
515: /* Compute the new step and safeguard it. */
516: stpf = PetscMin(ls->stepmax,stpf);
517: stpf = PetscMax(ls->stepmin,stpf);
518: *stp = stpf;
519: if (mtP->bracket && bound) {
520: if (*sty > *stx) {
521: *stp = PetscMin(*stx+0.66*(*sty-*stx),*stp);
522: } else {
523: *stp = PetscMax(*stx+0.66*(*sty-*stx),*stp);
524: }
525: }
526: return 0;
527: }