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: }