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) {
 21:     PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");
 22:   }

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

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

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

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

 54:   while (1) {
 55:     /* Call general purpose update function */
 56:     if (tao->ops->update) {
 57:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
 58:     }

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

 78:       VecCopy(tao->gradient, tao->stepdirection);
 79:       VecScale(tao->stepdirection, -1.0);
 80:     }

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

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

105:       VecCopy(tao->gradient, tao->stepdirection);
106:       VecScale(tao->stepdirection, -1.0);

108:       TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
109:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
110:       TaoAddLineSearchCounts(tao);

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

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

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

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

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

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

163:       case CG_PolakRibiere:
164:         beta = (gnorm2 - ginner) / gnorm2_old;
165:         break;

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

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

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

183:       default:
184:         beta = 0.0;
185:         break;
186:       }
187:     }

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

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

200: static PetscErrorCode TaoSetUp_CG(Tao tao)
201: {
202:   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: }

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

215:   if (tao->setupcalled) {
216:     VecDestroy(&cgP->X_old);
217:     VecDestroy(&cgP->G_old);
218:   }
219:   TaoLineSearchDestroy(&tao->linesearch);
220:   PetscFree(tao->data);
221:   return 0;
222: }

224: static PetscErrorCode TaoSetFromOptions_CG(PetscOptionItems *PetscOptionsObject,Tao tao)
225: {
226:   TAO_CG         *cgP = (TAO_CG*)tao->data;

228:   TaoLineSearchSetFromOptions(tao->linesearch);
229:   PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");
230:   PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,NULL);
231:   PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type,NULL);
232:   PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,NULL);
233:   PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,NULL);
234:   PetscOptionsTail();
235:   return 0;
236: }

238: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
239: {
240:   PetscBool      isascii;
241:   TAO_CG         *cgP = (TAO_CG*)tao->data;

243:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
244:   if (isascii) {
245:     PetscViewerASCIIPushTab(viewer);
246:     PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
247:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);
248:     PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);
249:     PetscViewerASCIIPopTab(viewer);
250:   }
251:   return 0;
252: }

254: /*MC
255:      TAOCG -   Nonlinear conjugate gradient method is an extension of the
256: nonlinear conjugate gradient solver for nonlinear optimization.

258:    Options Database Keys:
259: +      -tao_cg_eta <r> - restart tolerance
260: .      -tao_cg_type <taocg_type> - cg formula
261: .      -tao_cg_delta_min <r> - minimum delta value
262: -      -tao_cg_delta_max <r> - maximum delta value

264:   Notes:
265:      CG formulas are:
266:          "fr" - Fletcher-Reeves
267:          "pr" - Polak-Ribiere
268:          "prp" - Polak-Ribiere-Plus
269:          "hs" - Hestenes-Steifel
270:          "dy" - Dai-Yuan
271:   Level: beginner
272: M*/

274: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
275: {
276:   TAO_CG         *cgP;
277:   const char     *morethuente_type = TAOLINESEARCHMT;

279:   tao->ops->setup = TaoSetUp_CG;
280:   tao->ops->solve = TaoSolve_CG;
281:   tao->ops->view = TaoView_CG;
282:   tao->ops->setfromoptions = TaoSetFromOptions_CG;
283:   tao->ops->destroy = TaoDestroy_CG;

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

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

299:   PetscNewLog(tao,&cgP);
300:   tao->data = (void*)cgP;
301:   cgP->eta = 0.1;
302:   cgP->delta_min = 1e-7;
303:   cgP->delta_max = 100;
304:   cgP->cg_type = CG_PolakRibierePlus;
305:   return 0;
306: }