Actual source code: taocg.c

petsc-3.7.3 2016-08-01
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:          VecCopy(cgP->X_old, tao->solution);
114:          VecCopy(cgP->G_old, tao->gradient);
115:          delta = 1.0;
116:          VecCopy(tao->solution, tao->stepdirection);
117:          VecScale(tao->stepdirection, -1.0);

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

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

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

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

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

159:        case CG_PolakRibiere:
160:          beta = (gnorm2 - ginner) / gnorm2_old;
161:          break;

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

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

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

179:        default:
180:          beta = 0.0;
181:          break;
182:        }
183:      }

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

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

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

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

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

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

230: static PetscErrorCode TaoSetFromOptions_CG(PetscOptionItems *PetscOptionsObject,Tao tao)
231:  {
232:     TAO_CG         *cgP = (TAO_CG*)tao->data;

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

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

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

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

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

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


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

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

302:   /* Override default settings (unless already changed) */
303:   if (!tao->max_it_changed) tao->max_it = 2000;
304:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

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

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