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