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:   PetscFunctionBegin;
 21:   if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by cg algorithm\n"));

 23:   /*  Check convergence criteria */
 24:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
 25:   PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
 26:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

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

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

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

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

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

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

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

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

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

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

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

135:     /*  Check for bad value */
136:     PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
137:     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User-provided compute function generated Inf or NaN");

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

147:     /*  Check for restart condition */
148:     PetscCall(VecDot(tao->gradient, cgP->G_old, &ginner));
149:     if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
150:       /*  Gradients far from orthogonal; 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:         PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
169:         PetscCall(VecDot(cgP->G_old, tao->stepdirection, &gd_old));
170:         beta = (gnorm2 - ginner) / (gd - gd_old);
171:         break;

173:       case CG_DaiYuan:
174:         PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
175:         PetscCall(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:     PetscCall(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:   PetscFunctionReturn(PETSC_SUCCESS);
194: }

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

200:   PetscFunctionBegin;
201:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
202:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
203:   if (!cgP->X_old) PetscCall(VecDuplicate(tao->solution, &cgP->X_old));
204:   if (!cgP->G_old) PetscCall(VecDuplicate(tao->gradient, &cgP->G_old));
205:   PetscFunctionReturn(PETSC_SUCCESS);
206: }

208: static PetscErrorCode TaoDestroy_CG(Tao tao)
209: {
210:   TAO_CG *cgP = (TAO_CG *)tao->data;

212:   PetscFunctionBegin;
213:   if (tao->setupcalled) {
214:     PetscCall(VecDestroy(&cgP->X_old));
215:     PetscCall(VecDestroy(&cgP->G_old));
216:   }
217:   PetscCall(TaoLineSearchDestroy(&tao->linesearch));
218:   PetscCall(PetscFree(tao->data));
219:   PetscFunctionReturn(PETSC_SUCCESS);
220: }

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

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

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

242:   PetscFunctionBegin;
243:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
244:   if (isascii) {
245:     PetscCall(PetscViewerASCIIPushTab(viewer));
246:     PetscCall(PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]));
247:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", cgP->ngradsteps));
248:     PetscCall(PetscViewerASCIIPrintf(viewer, "Reset steps: %" PetscInt_FMT "\n", cgP->nresetsteps));
249:     PetscCall(PetscViewerASCIIPopTab(viewer));
250:   }
251:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
280:   tao->ops->setup          = TaoSetUp_CG;
281:   tao->ops->solve          = TaoSolve_CG;
282:   tao->ops->view           = TaoView_CG;
283:   tao->ops->setfromoptions = TaoSetFromOptions_CG;
284:   tao->ops->destroy        = TaoDestroy_CG;

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

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

300:   PetscCall(PetscNew(&cgP));
301:   tao->data      = (void *)cgP;
302:   cgP->eta       = 0.1;
303:   cgP->delta_min = 1e-7;
304:   cgP->delta_max = 100;
305:   cgP->cg_type   = CG_PolakRibierePlus;
306:   PetscFunctionReturn(PETSC_SUCCESS);
307: }