Actual source code: taocg.c

  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"};

 13: static PetscErrorCode TaoSolve_CG(Tao tao)
 14: {
 15:   TAO_CG                      *cgP       = (TAO_CG *)tao->data;
 16:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
 17:   PetscReal                    step      = 1.0, f, gnorm, gnorm2, delta, gd, ginner, beta;
 18:   PetscReal                    gd_old, gnorm2_old, f_old;

 20:   if (tao->XL || tao->XU || tao->ops->computebounds) PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");

 22:   /*  Check convergence criteria */
 23:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 24:   VecNorm(tao->gradient, NORM_2, &gnorm);

 27:   tao->reason = TAO_CONTINUE_ITERATING;
 28:   TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
 29:   TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
 30:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 31:   if (tao->reason != TAO_CONTINUE_ITERATING) return 0;

 33:   /*  Set initial direction to -gradient */
 34:   VecCopy(tao->gradient, tao->stepdirection);
 35:   VecScale(tao->stepdirection, -1.0);
 36:   gnorm2 = gnorm * gnorm;

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

 52:   while (1) {
 53:     /* Call general purpose update function */
 54:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);

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

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

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

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

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

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

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

118:         TaoLineSearchSetInitialStepLength(tao->linesearch, delta);
119:         TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
120:         TaoAddLineSearchCounts(tao);
121:         if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
122:           /*  Line search failed for last time -- give up */
123:           f = f_old;
124:           VecCopy(cgP->X_old, tao->solution);
125:           VecCopy(cgP->G_old, tao->gradient);
126:           step        = 0.0;
127:           tao->reason = TAO_DIVERGED_LS_FAILURE;
128:         }
129:       }
130:     }

132:     /*  Check for bad value */
133:     VecNorm(tao->gradient, NORM_2, &gnorm);

136:     /*  Check for termination */
137:     gnorm2 = gnorm * gnorm;
138:     tao->niter++;
139:     TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
140:     TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
141:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
142:     if (tao->reason != TAO_CONTINUE_ITERATING) break;

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

156:       case CG_PolakRibiere:
157:         beta = (gnorm2 - ginner) / gnorm2_old;
158:         break;

160:       case CG_PolakRibierePlus:
161:         beta = PetscMax((gnorm2 - ginner) / gnorm2_old, 0.0);
162:         break;

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

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

176:       default:
177:         beta = 0.0;
178:         break;
179:       }
180:     }

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

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

193: static PetscErrorCode TaoSetUp_CG(Tao tao)
194: {
195:   TAO_CG *cgP = (TAO_CG *)tao->data;

197:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
198:   if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
199:   if (!cgP->X_old) VecDuplicate(tao->solution, &cgP->X_old);
200:   if (!cgP->G_old) VecDuplicate(tao->gradient, &cgP->G_old);
201:   return 0;
202: }

204: static PetscErrorCode TaoDestroy_CG(Tao tao)
205: {
206:   TAO_CG *cgP = (TAO_CG *)tao->data;

208:   if (tao->setupcalled) {
209:     VecDestroy(&cgP->X_old);
210:     VecDestroy(&cgP->G_old);
211:   }
212:   TaoLineSearchDestroy(&tao->linesearch);
213:   PetscFree(tao->data);
214:   return 0;
215: }

217: static PetscErrorCode TaoSetFromOptions_CG(Tao tao, PetscOptionItems *PetscOptionsObject)
218: {
219:   TAO_CG *cgP = (TAO_CG *)tao->data;

221:   TaoLineSearchSetFromOptions(tao->linesearch);
222:   PetscOptionsHeadBegin(PetscOptionsObject, "Nonlinear Conjugate Gradient method for unconstrained optimization");
223:   PetscOptionsReal("-tao_cg_eta", "restart tolerance", "", cgP->eta, &cgP->eta, NULL);
224:   PetscOptionsEList("-tao_cg_type", "cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, NULL);
225:   PetscOptionsReal("-tao_cg_delta_min", "minimum delta value", "", cgP->delta_min, &cgP->delta_min, NULL);
226:   PetscOptionsReal("-tao_cg_delta_max", "maximum delta value", "", cgP->delta_max, &cgP->delta_max, NULL);
227:   PetscOptionsHeadEnd();
228:   return 0;
229: }

231: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
232: {
233:   PetscBool isascii;
234:   TAO_CG   *cgP = (TAO_CG *)tao->data;

236:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
237:   if (isascii) {
238:     PetscViewerASCIIPushTab(viewer);
239:     PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
240:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", cgP->ngradsteps);
241:     PetscViewerASCIIPrintf(viewer, "Reset steps: %" PetscInt_FMT "\n", cgP->nresetsteps);
242:     PetscViewerASCIIPopTab(viewer);
243:   }
244:   return 0;
245: }

247: /*MC
248:      TAOCG -   Nonlinear conjugate gradient method is an extension of the
249: nonlinear conjugate gradient solver for nonlinear optimization.

251:    Options Database Keys:
252: +      -tao_cg_eta <r> - restart tolerance
253: .      -tao_cg_type <taocg_type> - cg formula
254: .      -tao_cg_delta_min <r> - minimum delta value
255: -      -tao_cg_delta_max <r> - maximum delta value

257:   Notes:
258:      CG formulas are:
259:          "fr" - Fletcher-Reeves
260:          "pr" - Polak-Ribiere
261:          "prp" - Polak-Ribiere-Plus
262:          "hs" - Hestenes-Steifel
263:          "dy" - Dai-Yuan
264:   Level: beginner
265: M*/

267: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
268: {
269:   TAO_CG     *cgP;
270:   const char *morethuente_type = TAOLINESEARCHMT;

272:   tao->ops->setup          = TaoSetUp_CG;
273:   tao->ops->solve          = TaoSolve_CG;
274:   tao->ops->view           = TaoView_CG;
275:   tao->ops->setfromoptions = TaoSetFromOptions_CG;
276:   tao->ops->destroy        = TaoDestroy_CG;

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

282:   /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
283:   /*  method.  In particular, gtol should be less that 0.5; the value used in  */
284:   /*  Nocedal and Wright is 0.10.  We use the default values for the  */
285:   /*  linesearch because it seems to work better. */
286:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
287:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
288:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
289:   TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
290:   TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);

292:   PetscNew(&cgP);
293:   tao->data      = (void *)cgP;
294:   cgP->eta       = 0.1;
295:   cgP->delta_min = 1e-7;
296:   cgP->delta_max = 100;
297:   cgP->cg_type   = CG_PolakRibierePlus;
298:   return 0;
299: }