Actual source code: boundmorethuente.c

      1: /*$Id: morethuente.c,v 1.1 2002/08/11 01:15:06 benson Exp $*/

      3: #include "morethuente.h"  /*I "tao_solver.h
      4: " I*/

      6: static int TaoSetOptions_BoundLineSearch(TAO_SOLVER,void*);
      7: static int TaoApply_BoundLineSearch(TAO_SOLVER,TaoVec*,TaoVec*,TaoVec*,TaoVec*,
      8:                                double*,double*,int*,void*);

     10: /*------------------------------------------------------------*/
     13: static int TaoSetOptions_BoundLineSearch(TAO_SOLVER tao, void*linectx)
     14: {
     15:   TAO_LINESEARCH *ctx = (TAO_LINESEARCH *)linectx;
     16:   int            info;

     18:   TaoFunctionBegin;
     19:   info = TaoOptionsHead("Modified More-Thuente line line search options for bound constrained minimization");CHKERRQ(info);
     20:   info = TaoOptionsTail();CHKERRQ(info);
     21:   info = TaoSetOptions_LineSearch(tao,ctx);
     22:   TaoFunctionReturn(0);
     23: }
     24: /*------------------------------------------------------------*/

     26: /* ---------------------------------------------------------- */
     29: /* @
     30:    TaoApply_BoundLineSearch - This routine performs a line search algorithm.

     32:    Input Parameters:
     33: +  tao - TAO_SOLVER context
     34: .  X - current iterate (on output X contains new iterate, X + step*S)
     35: .  S - search direction
     36: .  f - objective function evaluated at X
     37: .  G - gradient evaluated at X
     38: .  W - work vector
     39: -  step - initial estimate of step length

     41:    Output parameters:
     42: +  f - objective function evaluated at new iterate, X + step*S
     43: .  G - gradient evaluated at new iterate, X + step*S
     44: .  X - new iterate
     45: .  gnorm - 2-norm of G
     46: -  step - final step length

     49:    Info is set to one of:
     50: +   0 - the line search succeeds; the sufficient decrease
     51:    condition and the directional derivative condition hold

     53:    negative number if an input parameter is invalid
     54: .   -1 -  step < 0
     55: .   -2 -  ftol < 0
     56: .   -3 -  rtol < 0
     57: .   -4 -  gtol < 0
     58: .   -5 -  stepmin < 0
     59: .   -6 -  stepmax < stepmin
     60: -   -7 -  maxfev < 0

     62:    positive number > 1 if the line search otherwise terminates
     63: +    2 -  Relative width of the interval of uncertainty is
     64:          at most rtol.
     65: .    3 -  Maximum number of function evaluations (maxfev) has
     66:          been reached.
     67: .    4 -  Step is at the lower bound, stepmin.
     68: .    5 -  Step is at the upper bound, stepmax.
     69: .    6 -  Rounding errors may prevent further progress.
     70:          There may not be a step that satisfies the
     71:          sufficient decrease and curvature conditions.
     72:          Tolerances may be too small.
     73: -    7 -  Search direction is not a descent direction.

     75:    Notes:
     76:    This algorithm is is a modification of the algorithm by More' and Thuente.
     77:    The modifications concern bounds.  This algorithm steps in the direction
     78:    passed into this routine.  This point get projected back into the feasible set.
     79:    In the context of bound constrained
     80:    optimization, there may not be a point in the piecewise linear
     81:    path that satisfies the Wolfe conditions.  When the active set
     82:    is changing, decrease in the objective function may be sufficient
     83:    to terminate this line search.

     85:    Note: Much of this coded is identical to the More' Thuente line search.
     86:    Variations to the code are commented.

     88:    Notes:
     89:    This routine is used within the following TAO boiund constrained minimization
     90:    solvers: Newton linesearch (tao_tron) and limited memory variable metric
     91:    (tao_blmvm).

     93:    Level: advanced

     95: .keywords: TAO_SOLVER, linesearch
     96: @ */
     97: static int TaoApply_BoundLineSearch(TAO_SOLVER tao,TaoVec* X,TaoVec* G,TaoVec* S,TaoVec* W,double *f,
     98:                         double *step,int *info2,void*ctx)
     99: {
    100:   TAO_LINESEARCH *neP = (TAO_LINESEARCH *) tao->linectx;
    101:   TaoVec *XL,*XU;
    102:   double    zero = 0.0, two = 2.0, p5 = 0.5, p66 = 0.66, xtrapf = 4.0;
    103:   double    finit, width, width1, dginit,fm, fxm, fym, dgm, dgxm, dgym;
    104:   double    dgx, dgy, dg, fx, fy, stx, sty, dgtest, ftest1=0.0,ftest2=0.0;
    105:   double    bstepmin1,bstepmin2,bstepmax;
    106:   double    dg1,dg2;
    107:   int       info, i, stage1;

    109: #if defined(PETSC_USE_COMPLEX)
    110:   PetscScalar    cdginit, cdg, cstep = 0.0;
    111: #endif

    113:   TaoFunctionBegin;
    114:   /* neP->stepmin - lower bound for step */
    115:   /* neP->stepmax - upper bound for step */
    116:   /* neP->rtol           - relative tolerance for an acceptable step */
    117:   /* neP->ftol           - tolerance for sufficient decrease condition */
    118:   /* neP->gtol           - tolerance for curvature condition */
    119:   /* neP->nfev           - number of function evaluations */
    120:   /* neP->maxfev  - maximum number of function evaluations */

    122:   /* Check input parameters for errors */
    123:   if (*step < zero) {
    124:     info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: step (%g) < 0\n",*step)); CHKERRQ(info);
    125:     *info2 = -1; TaoFunctionReturn(0);
    126:   } else if (neP->ftol < zero) {
    127:     info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: ftol (%g) < 0\n",neP->ftol)); CHKERRQ(info);
    128:     *info2 = -2; TaoFunctionReturn(0);
    129:   } else if (neP->rtol < zero) {
    130:     info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: rtol (%g) < 0\n",neP->rtol)); CHKERRQ(info);
    131:     *info2 = -3; TaoFunctionReturn(0);
    132:   } else if (neP->gtol < zero) {
    133:     info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: gtol (%g) < 0\n",neP->gtol)); CHKERRQ(info);
    134:     *info2 = -4; TaoFunctionReturn(0);
    135:   } else if (neP->stepmin < zero) {
    136:     info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: stepmin (%g) < 0\n",neP->stepmin)); CHKERRQ(info);
    137:     *info2 = -5; TaoFunctionReturn(0);
    138:   } else if (neP->stepmax < neP->stepmin) {
    139:     info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: stepmax (%g) < stepmin (%g)\n",
    140:                          neP->stepmax,neP->stepmin)); CHKERRQ(info);
    141:     *info2 = -6; TaoFunctionReturn(0);
    142:   } else if (neP->maxfev < zero) {
    143:     info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search error: maxfev (%d) < 0\n",neP->maxfev)); CHKERRQ(info);
    144:     *info2 = -7; TaoFunctionReturn(0);
    145:   }

    147:   /* Next 6 lines added by S. Benson to accomodate bounds on the variables */
    148:   /* Compute step length needed to make all variables equal a bound */
    149:   /* Compute the smallest steplength that will make one nonbinding variable equal the bound */
    150:   info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
    151:   info = S->Negate(); CHKERRQ(info);
    152:   info = S->BoundGradientProjection(S,XL,X,XU); CHKERRQ(info);
    153:   info = S->Negate(); CHKERRQ(info);
    154:   info = X->StepBoundInfo(XL,XU,S,&bstepmin1,&bstepmin2,&bstepmax); CHKERRQ(info);
    155:   neP->stepmax=TaoMin(bstepmax,1.0e+15);

    157:   /* Check that search direction is a descent direction */

    159: #if defined(PETSC_USE_COMPLEX)
    160:   info = G->Dot(S,&cdginit);CHKERRQ(info); dginit = TaoReal(cdginit);
    161: #else
    162:   info = G->Dot(S,&dginit);CHKERRQ(info);
    163: #endif

    166:   if (dginit >= zero) {
    167:     info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Search direction not a descent direction\n")); CHKERRQ(info);
    168:     *info2 = 7; TaoFunctionReturn(0);
    169:   }

    171:   /* Initialization */
    172:   neP->bracket = 0;
    173:   *info2          = 0;
    174:   stage1  = 1;
    175:   finit   = *f;
    176:   dgtest  = neP->ftol * dginit;
    177:   width   = neP->stepmax - neP->stepmin;
    178:   width1  = width * two;
    179:   info = W->CopyFrom(X);CHKERRQ(info);
    180:   /* Variable dictionary:
    181:      stx, fx, dgx - the step, function, and derivative at the best step
    182:      sty, fy, dgy - the step, function, and derivative at the other endpoint
    183:                    of the interval of uncertainty
    184:      step, f, dg - the step, function, and derivative at the current step */

    186:   stx = zero;
    187:   fx  = finit;
    188:   dgx = dginit;
    189:   sty = zero;
    190:   fy  = finit;
    191:   dgy = dginit;
    192:
    193:   neP->nfev = 0;
    194:   for (i=0; i< neP->maxfev; i++) {
    195:     /* Set min and max steps to correspond to the interval of uncertainty */
    196:     if (neP->bracket) {
    197:       neP->stepmin = TaoMin(stx,sty);
    198:       neP->stepmax = TaoMax(stx,sty);
    199:     } else {
    200:       neP->stepmin = stx;
    201:       neP->stepmax = *step + xtrapf * (*step - stx);
    202:     }

    204:     /* Force the step to be within the bounds */
    205:     *step = TaoMax(*step,neP->stepmin);
    206:     *step = TaoMin(*step,neP->stepmax);

    208:     /* If an unusual termination is to occur, then let step be the lowest
    209:        point obtained thus far */
    210:     if (((neP->bracket) && (*step <= neP->stepmin || *step >= neP->stepmax)) ||
    211:         ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol * neP->stepmax)) ||
    212:         (neP->nfev >= neP->maxfev - 1) || (neP->infoc == 0))
    213:       *step = stx;
    214:
    215:     info = PetscLogInfo((tao,"TaoApply_LineSearch: Try step size: %6.4e\n",*step)); CHKERRQ(info);
    216: #if defined(PETSC_USE_COMPLEX)
    217:     cstep = *step;
    218:     info = X->Waxpby(cstep,S,1.0,W);CHKERRQ(info);
    219: #else
    220:     info = X->Waxpby(*step,S,1.0,W);CHKERRQ(info);         /* X = W + step*S */
    221: #endif
    222:     /* This if statement added by S. Benson.   --  Project the solution, if necessary, to keep feasible */
    223:     if (*step >= bstepmin1){
    224:       info=X->Median(XL,X,XU);CHKERRQ(info);
    225:     }
    226:     info = TaoComputeMeritFunctionGradient(tao,X,f,G);CHKERRQ(info);
    227:     neP->nfev++;
    228: #if defined(PETSC_USE_COMPLEX)
    229:     info = G->Dot(S,&cdg);CHKERRQ(info); dg = TaoReal(cdg);
    230:     return 1;
    231: #else
    232:     //    info = G->Dot(S,&dg);CHKERRQ(info);                /* dg = G^T S */
    233:      /* Next 3 lines  added by S. Benson.   to compute an estimate of decrease in the actual direction and curvature */
    234:     info = G->Dot(W,&dg1);CHKERRQ(info);                /* dg = G^T S */
    235:     info = G->Dot(X,&dg2);CHKERRQ(info);                /* dg = G^T S */
    236:     dg=dg2-dg1;
    237: #endif
    238:     ftest1 = finit + (*step) * dgtest;

    240:     /* This test added by S. Benson.   --  Be satisfied with Armijo condition if the Projection X->Median(XL,X,XU) changed X. */
    241:     ftest2 = finit + (*step) * dgtest * neP->ftol;
    242:

    244:     /* Convergence testing */
    245:     if (((neP->bracket) && (*step <= neP->stepmin||*step >= neP->stepmax)) || (!neP->infoc)) {
    246:       info = PetscLogInfo((tao,"TaoApply_LineSearch:Rounding errors may prevent further progress.  May not be a step satisfying\n")); CHKERRQ(info);
    247:       info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:sufficient decrease and curvature conditions. Tolerances may be too small.\n")); CHKERRQ(info);
    248:       *info2 = 6; break;
    249:     }
    250:     if ((*step == neP->stepmax) && (*f <= ftest1) && (dg <= dgtest)) {
    251:       info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Step is at the upper bound, stepmax (%g)\n",neP->stepmax)); CHKERRQ(info);
    252:       *info2 = 5; break;
    253:     }
    254:     if ((*step == neP->stepmin) && (*f >= ftest1) && (dg >= dgtest)) {
    255:       info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Step is at the lower bound, stepmin (%g)\n",neP->stepmin)); CHKERRQ(info);
    256:       *info2 = 4; break;
    257:     }
    258:     if (neP->nfev >= neP->maxfev) {
    259:       info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Number of line search function evals (%d) > maximum (%d)\n",neP->nfev,neP->maxfev)); CHKERRQ(info);
    260:       *info2 = 3; break;
    261:     }
    262:     if ((neP->bracket) && (neP->stepmax - neP->stepmin <= neP->rtol*neP->stepmax)){
    263:       info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Relative width of interval of uncertainty is at most rtol (%g)\n",neP->rtol)); CHKERRQ(info);
    264:       *info2 = 2; break;
    265:     }
    266:     if ((*f <= ftest1) && (TaoAbsDouble(dg) <= neP->gtol*(-dginit))) {
    267:       info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search success: Sufficient decrease and directional deriv conditions hold, %6.4e, %6.4e\n",TaoAbsDouble(dg), -dginit )); CHKERRQ(info);
    268:       *info2 = 0; break;
    269:     }

    271:     /* This test added by S. Benson.   --  Be satisfied with Armijo condition if the Projection X->Median(XL,X,XU) changed X. */
    272:     if ( (*f <= ftest2) && ((*step) < bstepmin2) ) {
    273:       info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:Line search success: Sufficient decrease new bounds apply\n")); CHKERRQ(info);
    274:       *info2 = 0; break;
    275:     }

    277:     /* In the first stage, we seek a step for which the modified function
    278:         has a nonpositive value and nonnegative derivative */
    279:     if ((stage1) && (*f <= ftest1) && (dg >= dginit * TaoMin(neP->ftol, neP->gtol)))
    280:       stage1 = 0;

    282:     /* A modified function is used to predict the step only if we
    283:        have not obtained a step for which the modified function has a
    284:        nonpositive function value and nonnegative derivative, and if a
    285:        lower function value has been obtained but the decrease is not
    286:        sufficient */

    288:     if ((stage1) && (*f <= fx) && (*f > ftest1)) {
    289:       fm   = *f - *step * dgtest;        /* Define modified function */
    290:       fxm  = fx - stx * dgtest;                /* and derivatives */
    291:       fym  = fy - sty * dgtest;
    292:       dgm  = dg - dgtest;
    293:       dgxm = dgx - dgtest;
    294:       dgym = dgy - dgtest;

    296:       /* Update the interval of uncertainty and compute the new step */
    297:       info = TaoStep_LineSearch(tao,&stx,&fxm,&dgxm,&sty,&fym,&dgym,step,&fm,&dgm);CHKERRQ(info);

    299:       fx  = fxm + stx * dgtest;        /* Reset the function and */
    300:       fy  = fym + sty * dgtest;        /* gradient values */
    301:       dgx = dgxm + dgtest;
    302:       dgy = dgym + dgtest;
    303:     } else {
    304:       /* Update the interval of uncertainty and compute the new step */
    305:       info = TaoStep_LineSearch(tao,&stx,&fx,&dgx,&sty,&fy,&dgy,step,f,&dg);CHKERRQ(info);
    306:     }

    308:    /* Force a sufficient decrease in the interval of uncertainty */
    309:    if (neP->bracket) {
    310:      if (TaoAbsDouble(sty - stx) >= p66 * width1) *step = stx + p5*(sty - stx);
    311:        width1 = width;
    312:        width = TaoAbsDouble(sty - stx);
    313:      }
    314:    }

    316:   /* Finish computations */
    317:   info = PetscLogInfo((tao,"TaoApply_BoundLineSearch:%d function evals in line search, step = %10.4e\n",neP->nfev,*step)); CHKERRQ(info);
    318:   /*
    319:   info = G->Norm2(dginitt);CHKERRQ(info);
    320:   */
    321:   TaoFunctionReturn(0);
    322: }

    324: /* ---------------------------------------------------------- */

    328: /*@C
    329:    TaoCreateMoreThuenteBoundLineSearch - Create a line search

    331:    Input Parameters:
    332: +  tao - TAO_SOLVER context
    333: .  fftol - the sufficient descent parameter , greater than 0.
    334: -  ggtol - the curvature tolerance, greater than 0, less than 1.

    337:    Note:
    338:    If either fftol or ggtol is 0, default parameters will be used.

    340:    Note:
    341:    This algorithm is is a modification of the algorithm by More' and Thuente.
    342:    The modifications concern bounds.  This algorithm steps in the direction
    343:    passed into this routine.  This point get projected back into the feasible set.
    344:    It tries to satisfy the Wolfe conditions, but in the context of bound constrained
    345:    optimization, there may not be a point in the piecewise linear
    346:    path that satisfies the Wolfe conditions.  When the active set
    347:    is changing, decrease in the objective function may be sufficient
    348:    to terminate this line search.

    350:    Note:
    351:    This routine is used within the following TAO bound constrained
    352:    minimization solvers: Newton trust region (tao_tron) and limited memory variable
    353:    metric (tao_blmvm).

    355:    Level: developer

    357: .keywords: TAO_SOLVER, linesearch
    358: @*/
    359: int TaoCreateMoreThuenteBoundLineSearch(TAO_SOLVER tao, double fftol, double ggtol)
    360: {
    361:   int info;
    362:   TAO_LINESEARCH *neP;

    364:   TaoFunctionBegin;

    366:   info = TaoNew(TAO_LINESEARCH,&neP); CHKERRQ(info);
    367:   info = PetscLogObjectMemory(tao,sizeof(TAO_LINESEARCH)); CHKERRQ(info);

    369:   if (fftol<=0) fftol=0.001;
    370:   if (ggtol<=0) ggtol=0.9;

    372:   neP->ftol                  = fftol;
    373:   neP->rtol                  = 1.0e-10;
    374:   neP->gtol                  = ggtol;
    375:   neP->stepmin                  = 1.0e-30;
    376:   neP->stepmax                  = 1.0e+20;
    377:   neP->nfev                  = 0;
    378:   neP->bracket                  = 0;
    379:   neP->infoc              = 1;
    380:   neP->maxfev                  = 30;

    382:   info = TaoSetLineSearch(tao,0,
    383:                           TaoSetOptions_BoundLineSearch,
    384:                           TaoApply_BoundLineSearch,
    385:                           TaoView_LineSearch,
    386:                           TaoDestroy_LineSearch,
    387:                           (void *) neP);CHKERRQ(info);

    389:   TaoFunctionReturn(0);
    390: }

    395: /*

    398: */