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(TaoLineSearch ls, PetscOptionItems *PetscOptionsObject)
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, ostepmin, ostepmax;
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: ostepmax = ls->stepmax;
63: ostepmin = ls->stepmin;
65: if (ls->bounded) {
66: /* Compute step length needed to make all variables equal a bound */
67: /* Compute the smallest steplength that will make one nonbinding variable
68: equal the bound */
69: VecGetLocalSize(ls->upper, &n1);
70: VecGetLocalSize(mt->x, &n2);
71: VecGetSize(ls->upper, &nn1);
72: VecGetSize(mt->x, &nn2);
74: VecScale(s, -1.0);
75: VecBoundGradientProjection(s, x, ls->lower, ls->upper, s);
76: VecScale(s, -1.0);
77: VecStepBoundInfo(x, s, ls->lower, ls->upper, &bstepmin1, &bstepmin2, &bstepmax);
78: ls->stepmax = PetscMin(bstepmax, ls->stepmax);
79: }
81: VecDot(g, s, &dginit);
82: if (PetscIsInfOrNanReal(dginit)) {
83: PetscInfo(ls, "Initial Line Search step * g is Inf or Nan (%g)\n", (double)dginit);
84: ls->reason = TAOLINESEARCH_FAILED_INFORNAN;
85: return 0;
86: }
87: if (dginit >= 0.0) {
88: PetscInfo(ls, "Initial Line Search step * g is not descent direction (%g)\n", (double)dginit);
89: ls->reason = TAOLINESEARCH_FAILED_ASCENT;
90: return 0;
91: }
93: /* Initialization */
94: mt->bracket = 0;
95: stage1 = 1;
96: finit = *f;
97: dgtest = ls->ftol * dginit;
98: width = ls->stepmax - ls->stepmin;
99: width1 = width * 2.0;
100: VecCopy(x, mt->work);
101: /* Variable dictionary:
102: stx, fx, dgx - the step, function, and derivative at the best step
103: sty, fy, dgy - the step, function, and derivative at the other endpoint
104: of the interval of uncertainty
105: step, f, dg - the step, function, and derivative at the current step */
107: stx = 0.0;
108: fx = finit;
109: dgx = dginit;
110: sty = 0.0;
111: fy = finit;
112: dgy = dginit;
114: ls->step = ls->initstep;
115: for (i = 0; i < ls->max_funcs; i++) {
116: /* Set min and max steps to correspond to the interval of uncertainty */
117: if (mt->bracket) {
118: ls->stepmin = PetscMin(stx, sty);
119: ls->stepmax = PetscMax(stx, sty);
120: } else {
121: ls->stepmin = stx;
122: ls->stepmax = ls->step + xtrapf * (ls->step - stx);
123: }
125: /* Force the step to be within the bounds */
126: ls->step = PetscMax(ls->step, ls->stepmin);
127: ls->step = PetscMin(ls->step, ls->stepmax);
129: /* If an unusual termination is to occur, then let step be the lowest
130: point obtained thus far */
131: if (stx != 0 && ((mt->bracket && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || (mt->bracket && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) || (ls->nfeval + ls->nfgeval >= ls->max_funcs - 1) || mt->infoc == 0))
132: ls->step = stx;
134: VecWAXPY(mt->work, ls->step, s, x); /* W = X + step*S */
136: if (ls->bounded) VecMedian(ls->lower, mt->work, ls->upper, mt->work);
137: if (ls->usegts) {
138: TaoLineSearchComputeObjectiveAndGTS(ls, mt->work, f, &dg);
139: g_computed = PETSC_FALSE;
140: } else {
141: TaoLineSearchComputeObjectiveAndGradient(ls, mt->work, f, g);
142: g_computed = PETSC_TRUE;
143: if (ls->bounded) {
144: VecDot(g, x, &dg);
145: VecDot(g, mt->work, &dg2);
146: dg = (dg2 - dg) / ls->step;
147: } else {
148: VecDot(g, s, &dg);
149: }
150: }
152: /* update bracketing parameters in the MT context for printouts in monitor */
153: mt->stx = stx;
154: mt->fx = fx;
155: mt->dgx = dgx;
156: mt->sty = sty;
157: mt->fy = fy;
158: mt->dgy = dgy;
159: TaoLineSearchMonitor(ls, i + 1, *f, ls->step);
161: if (i == 0) ls->f_fullstep = *f;
163: if (PetscIsInfOrNanReal(*f) || PetscIsInfOrNanReal(dg)) {
164: /* User provided compute function generated Not-a-Number, assume
165: domain violation and set function value and directional
166: derivative to infinity. */
167: *f = PETSC_INFINITY;
168: dg = PETSC_INFINITY;
169: }
171: ftest1 = finit + ls->step * dgtest;
172: if (ls->bounded) ftest2 = finit + ls->step * dgtest * ls->ftol;
174: /* Convergence testing */
175: if ((*f - ftest1 <= PETSC_SMALL * PetscAbsReal(finit)) && (PetscAbsReal(dg) + ls->gtol * dginit <= 0.0)) {
176: PetscInfo(ls, "Line search success: Sufficient decrease and directional deriv conditions hold\n");
177: ls->reason = TAOLINESEARCH_SUCCESS;
178: break;
179: }
181: /* Check Armijo if beyond the first breakpoint */
182: if (ls->bounded && *f <= ftest2 && ls->step >= bstepmin2) {
183: PetscInfo(ls, "Line search success: Sufficient decrease.\n");
184: ls->reason = TAOLINESEARCH_SUCCESS;
185: break;
186: }
188: /* Checks for bad cases */
189: if ((mt->bracket && (ls->step <= ls->stepmin || ls->step >= ls->stepmax)) || !mt->infoc) {
190: PetscInfo(ls, "Rounding errors may prevent further progress. May not be a step satisfying\nsufficient decrease and curvature conditions. Tolerances may be too small.\n");
191: ls->reason = TAOLINESEARCH_HALTED_OTHER;
192: break;
193: }
194: if (ls->step == ls->stepmax && *f <= ftest1 && dg <= dgtest) {
195: PetscInfo(ls, "Step is at the upper bound, stepmax (%g)\n", (double)ls->stepmax);
196: ls->reason = TAOLINESEARCH_HALTED_UPPERBOUND;
197: break;
198: }
199: if (ls->step == ls->stepmin && *f >= ftest1 && dg >= dgtest) {
200: PetscInfo(ls, "Step is at the lower bound, stepmin (%g)\n", (double)ls->stepmin);
201: ls->reason = TAOLINESEARCH_HALTED_LOWERBOUND;
202: break;
203: }
204: if (mt->bracket && (ls->stepmax - ls->stepmin <= ls->rtol * ls->stepmax)) {
205: PetscInfo(ls, "Relative width of interval of uncertainty is at most rtol (%g)\n", (double)ls->rtol);
206: ls->reason = TAOLINESEARCH_HALTED_RTOL;
207: break;
208: }
210: /* In the first stage, we seek a step for which the modified function
211: has a nonpositive value and nonnegative derivative */
212: if (stage1 && *f <= ftest1 && dg >= dginit * PetscMin(ls->ftol, ls->gtol)) stage1 = 0;
214: /* A modified function is used to predict the step only if we
215: have not obtained a step for which the modified function has a
216: nonpositive function value and nonnegative derivative, and if a
217: lower function value has been obtained but the decrease is not
218: sufficient */
220: if (stage1 && *f <= fx && *f > ftest1) {
221: fm = *f - ls->step * dgtest; /* Define modified function */
222: fxm = fx - stx * dgtest; /* and derivatives */
223: fym = fy - sty * dgtest;
224: dgm = dg - dgtest;
225: dgxm = dgx - dgtest;
226: dgym = dgy - dgtest;
228: /* if (dgxm * (ls->step - stx) >= 0.0) */
229: /* Update the interval of uncertainty and compute the new step */
230: Tao_mcstep(ls, &stx, &fxm, &dgxm, &sty, &fym, &dgym, &ls->step, &fm, &dgm);
232: fx = fxm + stx * dgtest; /* Reset the function and */
233: fy = fym + sty * dgtest; /* gradient values */
234: dgx = dgxm + dgtest;
235: dgy = dgym + dgtest;
236: } else {
237: /* Update the interval of uncertainty and compute the new step */
238: Tao_mcstep(ls, &stx, &fx, &dgx, &sty, &fy, &dgy, &ls->step, f, &dg);
239: }
241: /* Force a sufficient decrease in the interval of uncertainty */
242: if (mt->bracket) {
243: if (PetscAbsReal(sty - stx) >= 0.66 * width1) ls->step = stx + 0.5 * (sty - stx);
244: width1 = width;
245: width = PetscAbsReal(sty - stx);
246: }
247: }
248: if (ls->nfeval + ls->nfgeval > ls->max_funcs) {
249: PetscInfo(ls, "Number of line search function evals (%" PetscInt_FMT ") > maximum (%" PetscInt_FMT ")\n", ls->nfeval + ls->nfgeval, ls->max_funcs);
250: ls->reason = TAOLINESEARCH_HALTED_MAXFCN;
251: }
252: ls->stepmax = ostepmax;
253: ls->stepmin = ostepmin;
255: /* Finish computations */
256: PetscInfo(ls, "%" PetscInt_FMT " function evals in line search, step = %g\n", ls->nfeval + ls->nfgeval, (double)ls->step);
258: /* Set new solution vector and compute gradient if needed */
259: VecCopy(mt->work, x);
260: if (!g_computed) TaoLineSearchComputeGradient(ls, mt->work, g);
261: return 0;
262: }
264: /*MC
265: TAOLINESEARCHMT - Line-search type with cubic interpolation that satisfies both the sufficient decrease and
266: curvature conditions. This method can take step lengths greater than 1.
268: More-Thuente line-search can be selected with "-tao_ls_type more-thuente".
270: References:
271: . * - JORGE J. MORE AND DAVID J. THUENTE, LINE SEARCH ALGORITHMS WITH GUARANTEED SUFFICIENT DECREASE.
272: ACM Trans. Math. Software 20, no. 3 (1994): 286-307.
274: Level: developer
276: .seealso: `TaoLineSearchCreate()`, `TaoLineSearchSetType()`, `TaoLineSearchApply()`
278: .keywords: Tao, linesearch
279: M*/
280: PETSC_EXTERN PetscErrorCode TaoLineSearchCreate_MT(TaoLineSearch ls)
281: {
282: TaoLineSearch_MT *ctx;
285: PetscNew(&ctx);
286: ctx->bracket = 0;
287: ctx->infoc = 1;
288: ls->data = (void *)ctx;
289: ls->initstep = 1.0;
290: ls->ops->setup = NULL;
291: ls->ops->reset = NULL;
292: ls->ops->apply = TaoLineSearchApply_MT;
293: ls->ops->destroy = TaoLineSearchDestroy_MT;
294: ls->ops->setfromoptions = TaoLineSearchSetFromOptions_MT;
295: ls->ops->monitor = TaoLineSearchMonitor_MT;
296: return 0;
297: }
299: /*
300: The subroutine mcstep is taken from the work of Jorge Nocedal.
301: this is a variant of More' and Thuente's routine.
303: subroutine mcstep
305: the purpose of mcstep is to compute a safeguarded step for
306: a linesearch and to update an interval of uncertainty for
307: a minimizer of the function.
309: the parameter stx contains the step with the least function
310: value. the parameter stp contains the current step. it is
311: assumed that the derivative at stx is negative in the
312: direction of the step. if bracket is set true then a
313: minimizer has been bracketed in an interval of uncertainty
314: with endpoints stx and sty.
316: the subroutine statement is
318: subroutine mcstep(stx,fx,dx,sty,fy,dy,stp,fp,dp,bracket,
319: stpmin,stpmax,info)
321: where
323: stx, fx, and dx are variables which specify the step,
324: the function, and the derivative at the best step obtained
325: so far. The derivative must be negative in the direction
326: of the step, that is, dx and stp-stx must have opposite
327: signs. On output these parameters are updated appropriately.
329: sty, fy, and dy are variables which specify the step,
330: the function, and the derivative at the other endpoint of
331: the interval of uncertainty. On output these parameters are
332: updated appropriately.
334: stp, fp, and dp are variables which specify the step,
335: the function, and the derivative at the current step.
336: If bracket is set true then on input stp must be
337: between stx and sty. On output stp is set to the new step.
339: bracket is a logical variable which specifies if a minimizer
340: has been bracketed. If the minimizer has not been bracketed
341: then on input bracket must be set false. If the minimizer
342: is bracketed then on output bracket is set true.
344: stpmin and stpmax are input variables which specify lower
345: and upper bounds for the step.
347: info is an integer output variable set as follows:
348: if info = 1,2,3,4,5, then the step has been computed
349: according to one of the five cases below. otherwise
350: info = 0, and this indicates improper input parameters.
352: subprograms called
354: fortran-supplied ... abs,max,min,sqrt
356: argonne national laboratory. minpack project. june 1983
357: jorge j. more', david j. thuente
359: */
361: static PetscErrorCode Tao_mcstep(TaoLineSearch ls, PetscReal *stx, PetscReal *fx, PetscReal *dx, PetscReal *sty, PetscReal *fy, PetscReal *dy, PetscReal *stp, PetscReal *fp, PetscReal *dp)
362: {
363: TaoLineSearch_MT *mtP = (TaoLineSearch_MT *)ls->data;
364: PetscReal gamma1, p, q, r, s, sgnd, stpc, stpf, stpq, theta;
365: PetscInt bound;
367: /* Check the input parameters for errors */
368: mtP->infoc = 0;
373: /* Determine if the derivatives have opposite sign */
374: sgnd = *dp * (*dx / PetscAbsReal(*dx));
376: if (*fp > *fx) {
377: /* Case 1: a higher function value.
378: The minimum is bracketed. If the cubic step is closer
379: to stx than the quadratic step, the cubic step is taken,
380: else the average of the cubic and quadratic steps is taken. */
382: mtP->infoc = 1;
383: bound = 1;
384: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
385: s = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
386: s = PetscMax(s, PetscAbsReal(*dp));
387: gamma1 = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s));
388: if (*stp < *stx) gamma1 = -gamma1;
389: /* Can p be 0? Check */
390: p = (gamma1 - *dx) + theta;
391: q = ((gamma1 - *dx) + gamma1) + *dp;
392: r = p / q;
393: stpc = *stx + r * (*stp - *stx);
394: stpq = *stx + ((*dx / ((*fx - *fp) / (*stp - *stx) + *dx)) * 0.5) * (*stp - *stx);
396: if (PetscAbsReal(stpc - *stx) < PetscAbsReal(stpq - *stx)) stpf = stpc;
397: else stpf = stpc + 0.5 * (stpq - stpc);
398: mtP->bracket = 1;
399: } else if (sgnd < 0.0) {
400: /* Case 2: A lower function value and derivatives of
401: opposite sign. The minimum is bracketed. If the cubic
402: step is closer to stx than the quadratic (secant) step,
403: the cubic step is taken, else the quadratic step is taken. */
405: mtP->infoc = 2;
406: bound = 0;
407: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
408: s = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
409: s = PetscMax(s, PetscAbsReal(*dp));
410: gamma1 = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s));
411: if (*stp > *stx) gamma1 = -gamma1;
412: p = (gamma1 - *dp) + theta;
413: q = ((gamma1 - *dp) + gamma1) + *dx;
414: r = p / q;
415: stpc = *stp + r * (*stx - *stp);
416: stpq = *stp + (*dp / (*dp - *dx)) * (*stx - *stp);
418: if (PetscAbsReal(stpc - *stp) > PetscAbsReal(stpq - *stp)) stpf = stpc;
419: else stpf = stpq;
420: mtP->bracket = 1;
421: } else if (PetscAbsReal(*dp) < PetscAbsReal(*dx)) {
422: /* Case 3: A lower function value, derivatives of the
423: same sign, and the magnitude of the derivative decreases.
424: The cubic step is only used if the cubic tends to infinity
425: in the direction of the step or if the minimum of the cubic
426: is beyond stp. Otherwise the cubic step is defined to be
427: either stepmin or stepmax. The quadratic (secant) step is also
428: computed and if the minimum is bracketed then the step
429: closest to stx is taken, else the step farthest away is taken. */
431: mtP->infoc = 3;
432: bound = 1;
433: theta = 3 * (*fx - *fp) / (*stp - *stx) + *dx + *dp;
434: s = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dx));
435: s = PetscMax(s, PetscAbsReal(*dp));
437: /* The case gamma1 = 0 only arises if the cubic does not tend
438: to infinity in the direction of the step. */
439: gamma1 = s * PetscSqrtScalar(PetscMax(0.0, PetscPowScalar(theta / s, 2.0) - (*dx / s) * (*dp / s)));
440: if (*stp > *stx) gamma1 = -gamma1;
441: p = (gamma1 - *dp) + theta;
442: q = (gamma1 + (*dx - *dp)) + gamma1;
443: r = p / q;
444: if (r < 0.0 && gamma1 != 0.0) stpc = *stp + r * (*stx - *stp);
445: else if (*stp > *stx) stpc = ls->stepmax;
446: else stpc = ls->stepmin;
447: stpq = *stp + (*dp / (*dp - *dx)) * (*stx - *stp);
449: if (mtP->bracket) {
450: if (PetscAbsReal(*stp - stpc) < PetscAbsReal(*stp - stpq)) stpf = stpc;
451: else stpf = stpq;
452: } else {
453: if (PetscAbsReal(*stp - stpc) > PetscAbsReal(*stp - stpq)) stpf = stpc;
454: else stpf = stpq;
455: }
456: } else {
457: /* Case 4: A lower function value, derivatives of the
458: same sign, and the magnitude of the derivative does
459: not decrease. If the minimum is not bracketed, the step
460: is either stpmin or stpmax, else the cubic step is taken. */
462: mtP->infoc = 4;
463: bound = 0;
464: if (mtP->bracket) {
465: theta = 3 * (*fp - *fy) / (*sty - *stp) + *dy + *dp;
466: s = PetscMax(PetscAbsReal(theta), PetscAbsReal(*dy));
467: s = PetscMax(s, PetscAbsReal(*dp));
468: gamma1 = s * PetscSqrtScalar(PetscPowScalar(theta / s, 2.0) - (*dy / s) * (*dp / s));
469: if (*stp > *sty) gamma1 = -gamma1;
470: p = (gamma1 - *dp) + theta;
471: q = ((gamma1 - *dp) + gamma1) + *dy;
472: r = p / q;
473: stpc = *stp + r * (*sty - *stp);
474: stpf = stpc;
475: } else if (*stp > *stx) {
476: stpf = ls->stepmax;
477: } else {
478: stpf = ls->stepmin;
479: }
480: }
482: /* Update the interval of uncertainty. This update does not
483: depend on the new step or the case analysis above. */
485: if (*fp > *fx) {
486: *sty = *stp;
487: *fy = *fp;
488: *dy = *dp;
489: } else {
490: if (sgnd < 0.0) {
491: *sty = *stx;
492: *fy = *fx;
493: *dy = *dx;
494: }
495: *stx = *stp;
496: *fx = *fp;
497: *dx = *dp;
498: }
500: /* Compute the new step and safeguard it. */
501: stpf = PetscMin(ls->stepmax, stpf);
502: stpf = PetscMax(ls->stepmin, stpf);
503: *stp = stpf;
504: if (mtP->bracket && bound) {
505: if (*sty > *stx) *stp = PetscMin(*stx + 0.66 * (*sty - *stx), *stp);
506: else *stp = PetscMax(*stx + 0.66 * (*sty - *stx), *stp);
507: }
508: return 0;
509: }