Actual source code: taocg.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/unconstrained/impls/cg/taocg.h>

  4: #define CG_FletcherReeves       0
  5: #define CG_PolakRibiere         1
  6: #define CG_PolakRibierePlus     2
  7: #define CG_HestenesStiefel      3
  8: #define CG_DaiYuan              4
  9: #define CG_Types                5

 11:  static const char *CG_Table[64] = {"fr", "pr", "prp", "hs", "dy"};

 15:  static PetscErrorCode TaoSolve_CG(Tao tao)
 16:  {
 17:    TAO_CG                       *cgP = (TAO_CG*)tao->data;
 18:    PetscErrorCode               ierr;
 19:    TaoConvergedReason           reason = TAO_CONTINUE_ITERATING;
 20:    TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
 21:    PetscReal                    step=1.0,f,gnorm,gnorm2,delta,gd,ginner,beta;
 22:    PetscReal                    gd_old,gnorm2_old,f_old;

 25:    if (tao->XL || tao->XU || tao->ops->computebounds) {
 26:      PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");
 27:    }

 29:    /*  Check convergence criteria */
 30:    TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 31:    VecNorm(tao->gradient,NORM_2,&gnorm);
 32:    if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");

 34:    TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);
 35:    if (reason != TAO_CONTINUE_ITERATING) return(0);

 37:    /*  Set initial direction to -gradient */
 38:    VecCopy(tao->gradient, tao->stepdirection);
 39:    VecScale(tao->stepdirection, -1.0);
 40:    gnorm2 = gnorm*gnorm;

 42:    /*  Set initial scaling for the function */
 43:    if (f != 0.0) {
 44:      delta = 2.0*PetscAbsScalar(f) / gnorm2;
 45:      delta = PetscMax(delta,cgP->delta_min);
 46:      delta = PetscMin(delta,cgP->delta_max);
 47:    } else {
 48:      delta = 2.0 / gnorm2;
 49:      delta = PetscMax(delta,cgP->delta_min);
 50:      delta = PetscMin(delta,cgP->delta_max);
 51:    }
 52:    /*  Set counter for gradient and reset steps */
 53:    cgP->ngradsteps = 0;
 54:    cgP->nresetsteps = 0;

 56:    while (1) {
 57:      /*  Save the current gradient information */
 58:      f_old = f;
 59:      gnorm2_old = gnorm2;
 60:      VecCopy(tao->solution, cgP->X_old);
 61:      VecCopy(tao->gradient, cgP->G_old);
 62:      VecDot(tao->gradient, tao->stepdirection, &gd);
 63:      if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
 64:        ++cgP->ngradsteps;
 65:        if (f != 0.0) {
 66:          delta = 2.0*PetscAbsScalar(f) / gnorm2;
 67:          delta = PetscMax(delta,cgP->delta_min);
 68:          delta = PetscMin(delta,cgP->delta_max);
 69:        } else {
 70:          delta = 2.0 / gnorm2;
 71:          delta = PetscMax(delta,cgP->delta_min);
 72:          delta = PetscMin(delta,cgP->delta_max);
 73:        }

 75:        VecCopy(tao->gradient, tao->stepdirection);
 76:        VecScale(tao->stepdirection, -1.0);
 77:      }

 79:      /*  Search direction for improving point */
 80:      TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
 81:      TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
 82:      TaoAddLineSearchCounts(tao);
 83:      if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
 84:        /*  Linesearch failed */
 85:        /*  Reset factors and use scaled gradient step */
 86:        ++cgP->nresetsteps;
 87:        f = f_old;
 88:        gnorm2 = gnorm2_old;
 89:        VecCopy(cgP->X_old, tao->solution);
 90:        VecCopy(cgP->G_old, tao->gradient);

 92:        if (f != 0.0) {
 93:          delta = 2.0*PetscAbsScalar(f) / gnorm2;
 94:          delta = PetscMax(delta,cgP->delta_min);
 95:          delta = PetscMin(delta,cgP->delta_max);
 96:        } else {
 97:          delta = 2.0 / gnorm2;
 98:          delta = PetscMax(delta,cgP->delta_min);
 99:          delta = PetscMin(delta,cgP->delta_max);
100:        }

102:        VecCopy(tao->gradient, tao->stepdirection);
103:        VecScale(tao->stepdirection, -1.0);

105:        TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
106:        TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
107:        TaoAddLineSearchCounts(tao);

109:        if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
110:          /*  Linesearch failed again */
111:          /*  switch to unscaled gradient */
112:          f = f_old;
113:          gnorm2 = gnorm2_old;
114:          VecCopy(cgP->X_old, tao->solution);
115:          VecCopy(cgP->G_old, tao->gradient);
116:          delta = 1.0;
117:          VecCopy(tao->solution, tao->stepdirection);
118:          VecScale(tao->stepdirection, -1.0);

120:          TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
121:          TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
122:          TaoAddLineSearchCounts(tao);
123:          if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {

125:            /*  Line search failed for last time -- give up */
126:            f = f_old;
127:            gnorm2 = gnorm2_old;
128:            VecCopy(cgP->X_old, tao->solution);
129:            VecCopy(cgP->G_old, tao->gradient);
130:            step = 0.0;
131:            reason = TAO_DIVERGED_LS_FAILURE;
132:            tao->reason = TAO_DIVERGED_LS_FAILURE;
133:          }
134:        }
135:      }

137:      /*  Check for bad value */
138:      VecNorm(tao->gradient,NORM_2,&gnorm);
139:      if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN");

141:      /*  Check for termination */
142:      gnorm2 =gnorm * gnorm;
143:      tao->niter++;
144:      TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);
145:      if (reason != TAO_CONTINUE_ITERATING) {
146:        break;
147:      }

149:      /*  Check for restart condition */
150:      VecDot(tao->gradient, cgP->G_old, &ginner);
151:      if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
152:        /*  Gradients far from orthognal; use steepest descent direction */
153:        beta = 0.0;
154:      } else {
155:        /*  Gradients close to orthogonal; use conjugate gradient formula */
156:        switch (cgP->cg_type) {
157:        case CG_FletcherReeves:
158:          beta = gnorm2 / gnorm2_old;
159:          break;

161:        case CG_PolakRibiere:
162:          beta = (gnorm2 - ginner) / gnorm2_old;
163:          break;

165:        case CG_PolakRibierePlus:
166:          beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
167:          break;

169:        case CG_HestenesStiefel:
170:          VecDot(tao->gradient, tao->stepdirection, &gd);
171:          VecDot(cgP->G_old, tao->stepdirection, &gd_old);
172:          beta = (gnorm2 - ginner) / (gd - gd_old);
173:          break;

175:        case CG_DaiYuan:
176:          VecDot(tao->gradient, tao->stepdirection, &gd);
177:          VecDot(cgP->G_old, tao->stepdirection, &gd_old);
178:          beta = gnorm2 / (gd - gd_old);
179:          break;

181:        default:
182:          beta = 0.0;
183:          break;
184:        }
185:      }

187:      /*  Compute the direction d=-g + beta*d */
188:      VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);

190:      /*  update initial steplength choice */
191:      delta = 1.0;
192:      delta = PetscMax(delta, cgP->delta_min);
193:      delta = PetscMin(delta, cgP->delta_max);
194:    }
195:    return(0);
196:  }

200:  static PetscErrorCode TaoSetUp_CG(Tao tao)
201:  {
202:    TAO_CG         *cgP = (TAO_CG*)tao->data;

206:    if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);}
207:    if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
208:    if (!cgP->X_old) {VecDuplicate(tao->solution,&cgP->X_old);}
209:    if (!cgP->G_old) {VecDuplicate(tao->gradient,&cgP->G_old); }
210:     return(0);
211:  }

215:  static PetscErrorCode TaoDestroy_CG(Tao tao)
216:  {
217:    TAO_CG         *cgP = (TAO_CG*) tao->data;

221:    if (tao->setupcalled) {
222:      VecDestroy(&cgP->X_old);
223:      VecDestroy(&cgP->G_old);
224:    }
225:    TaoLineSearchDestroy(&tao->linesearch);
226:    PetscFree(tao->data);
227:    return(0);
228:  }

232: static PetscErrorCode TaoSetFromOptions_CG(PetscOptions *PetscOptionsObject,Tao tao)
233:  {
234:     TAO_CG         *cgP = (TAO_CG*)tao->data;

238:     TaoLineSearchSetFromOptions(tao->linesearch);
239:     PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");
240:     PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,NULL);
241:     PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type,NULL);
242:     PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,NULL);
243:     PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,NULL);
244:    PetscOptionsTail();
245:    return(0);
246: }

250: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
251: {
252:   PetscBool      isascii;
253:   TAO_CG         *cgP = (TAO_CG*)tao->data;

257:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
258:   if (isascii) {
259:     PetscViewerASCIIPushTab(viewer);
260:     PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
261:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);
262:     ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);
263:     PetscViewerASCIIPopTab(viewer);
264:   }
265:   return(0);
266: }

268: /*MC
269:      TAOCG -   Nonlinear conjugate gradient method is an extension of the
270: nonlinear conjugate gradient solver for nonlinear optimization.

272:    Options Database Keys:
273: +      -tao_cg_eta <r> - restart tolerance
274: .      -tao_cg_type <taocg_type> - cg formula
275: .      -tao_cg_delta_min <r> - minimum delta value
276: -      -tao_cg_delta_max <r> - maximum delta value

278:   Notes:
279:      CG formulas are:
280:          "fr" - Fletcher-Reeves
281:          "pr" - Polak-Ribiere
282:          "prp" - Polak-Ribiere-Plus
283:          "hs" - Hestenes-Steifel
284:          "dy" - Dai-Yuan
285:   Level: beginner
286: M*/


291: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
292: {
293:   TAO_CG         *cgP;
294:   const char     *morethuente_type = TAOLINESEARCHMT;

298:   tao->ops->setup = TaoSetUp_CG;
299:   tao->ops->solve = TaoSolve_CG;
300:   tao->ops->view = TaoView_CG;
301:   tao->ops->setfromoptions = TaoSetFromOptions_CG;
302:   tao->ops->destroy = TaoDestroy_CG;

304:   /* Override default settings (unless already changed) */
305:   if (!tao->max_it_changed) tao->max_it = 2000;
306:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
307:   if (!tao->fatol_changed) tao->fatol = 1e-4;
308:   if (!tao->frtol_changed) tao->frtol = 1e-4;

310:   /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
311:   /*  method.  In particular, gtol should be less that 0.5; the value used in  */
312:   /*  Nocedal and Wright is 0.10.  We use the default values for the  */
313:   /*  linesearch because it seems to work better. */
314:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
315:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
316:   TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
317:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

319:   PetscNewLog(tao,&cgP);
320:   tao->data = (void*)cgP;
321:   cgP->eta = 0.1;
322:   cgP->delta_min = 1e-7;
323:   cgP->delta_max = 100;
324:   cgP->cg_type = CG_PolakRibierePlus;
325:   return(0);
326: }