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: */