Actual source code: cg.c

  1: #include "taolinesearch.h"
  2: #include "cg.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] = {
 12:   "fr", "pr", "prp", "hs", "dy"
 13: };


 18: static PetscErrorCode TaoSolve_CG(TaoSolver tao)
 19: {
 20:     TAO_CG *cgP = (TAO_CG*)tao->data;
 22:     TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING;
 23:     TaoLineSearchTerminationReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
 24:     PetscReal step=1.0,f,gnorm,gnorm2,delta,gd,ginner,beta;
 25:     PetscReal gd_old,gnorm2_old,f_old;
 26:     PetscInt iter=0;

 28:     

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

 35:     /*  Check convergence criteria */
 36:     TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient); 
 37:     VecNorm(tao->gradient,NORM_2,&gnorm); 
 38:     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
 39:         SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
 40:     }
 41:     
 42:     TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); 
 43:     if (reason != TAO_CONTINUE_ITERATING) {
 44:         return(0);
 45:     }

 47:     
 48:     /*  Set initial direction to -gradient */
 49:     VecCopy(tao->gradient, tao->stepdirection); 
 50:     VecScale(tao->stepdirection, -1.0); 
 51:     gnorm2 = gnorm*gnorm;
 52:     
 53:     /*  Set initial scaling for the function */
 54:     if (f != 0.0) {
 55:         delta = 2.0*PetscAbsScalar(f) / gnorm2;
 56:         delta = PetscMax(delta,cgP->delta_min);
 57:         delta = PetscMin(delta,cgP->delta_max);
 58:     } else {
 59:         delta = 2.0 / gnorm2;
 60:         delta = PetscMax(delta,cgP->delta_min);
 61:         delta = PetscMin(delta,cgP->delta_max);
 62:     }
 63:     /*  Set counter for gradient and reset steps */
 64:     cgP->ngradsteps = 0;
 65:     cgP->nresetsteps = 0;
 66:     
 67:     while (1) {
 68:         /*  Save the current gradient information */
 69:         f_old = f;
 70:         gnorm2_old = gnorm2;
 71:         VecCopy(tao->solution, cgP->X_old); 
 72:         VecCopy(tao->gradient, cgP->G_old); 
 73:         VecDot(tao->gradient, tao->stepdirection, &gd); 
 74:         if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
 75:             ++cgP->ngradsteps;
 76:             if (f != 0.0) {
 77:                 delta = 2.0*PetscAbsScalar(f) / gnorm2;
 78:                 delta = PetscMax(delta,cgP->delta_min);
 79:                 delta = PetscMin(delta,cgP->delta_max);
 80:             } else {
 81:                 delta = 2.0 / gnorm2;
 82:                 delta = PetscMax(delta,cgP->delta_min);
 83:                 delta = PetscMin(delta,cgP->delta_max);
 84:             }
 85:             
 86:             VecCopy(tao->gradient, tao->stepdirection); 
 87:             VecScale(tao->stepdirection, -1.0); 
 88:         }
 89:         
 90:         /*  Search direction for improving point */
 91:         TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
 92:         TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status); 
 93:         TaoAddLineSearchCounts(tao); 
 94:         if (ls_status < 0) {
 95:             /*  Linesearch failed */
 96:             /*  Reset factors and use scaled gradient step */
 97:             ++cgP->nresetsteps;
 98:             f = f_old;
 99:             gnorm2 = gnorm2_old;
100:             VecCopy(cgP->X_old, tao->solution); 
101:             VecCopy(cgP->G_old, tao->gradient); 
102:             
103:             if (f != 0.0) {
104:                 delta = 2.0*PetscAbsScalar(f) / gnorm2;
105:                 delta = PetscMax(delta,cgP->delta_min);
106:                 delta = PetscMin(delta,cgP->delta_max);
107:             } else {
108:                 delta = 2.0 / gnorm2;
109:                 delta = PetscMax(delta,cgP->delta_min);
110:                 delta = PetscMin(delta,cgP->delta_max);
111:             }
112:             
113:             VecCopy(tao->gradient, tao->stepdirection); 
114:             VecScale(tao->stepdirection, -1.0); 

116:             TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
117:             TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status); 
118:             TaoAddLineSearchCounts(tao); 
119:             
120:             if (ls_status < 0) {
121:                 /*  Linesearch failed again */
122:                 /*  switch to unscaled gradient */
123:                 f = f_old;
124:                 gnorm2 = gnorm2_old;
125:                 VecCopy(cgP->X_old, tao->solution); 
126:                 VecCopy(cgP->G_old, tao->gradient); 
127:                 delta = 1.0;
128:                 VecCopy(tao->solution, tao->stepdirection); 
129:                 VecScale(tao->stepdirection, -1.0); 

131:                 TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
132:                 TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status); 
133:                 TaoAddLineSearchCounts(tao); 
134:                 if (ls_status < 0) {
135:                     /*  Line search failed for last time -- give up */
136:                     f = f_old;
137:                     gnorm2 = gnorm2_old;
138:                     VecCopy(cgP->X_old, tao->solution); 
139:                     VecCopy(cgP->G_old, tao->gradient); 
140:                     step = 0.0;
141:                 }
142:             }
143:         }
144:         
145:         /*  Check for bad value */
146:         VecNorm(tao->gradient,NORM_2,&gnorm); 
147:         if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
148:             SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN");
149:         }
150:         
151:         /*  Check for termination */
152:         gnorm2 =gnorm * gnorm;
153:         iter++;
154:         TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason); 
155:         if (reason != TAO_CONTINUE_ITERATING) {
156:             break;
157:         }

159:         /*  Check for restart condition */
160:         VecDot(tao->gradient, cgP->G_old, &ginner); 
161:         if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
162:             /*  Gradients far from orthognal; use steepest descent direction */
163:             beta = 0.0;
164:         } else {
165:             /*  Gradients close to orthogonal; use conjugate gradient formula */
166:             switch (cgP->cg_type) {
167:                 case CG_FletcherReeves:
168:                     beta = gnorm2 / gnorm2_old;
169:                     break;
170:                     
171:                 case CG_PolakRibiere:
172:                     beta = (gnorm2 - ginner) / gnorm2_old;
173:                     break;
174:                 
175:                 case CG_PolakRibierePlus:
176:                     beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
177:                     break;

179:                 case CG_HestenesStiefel:
180:                     VecDot(tao->gradient, tao->stepdirection, &gd); 
181:                     VecDot(cgP->G_old, tao->stepdirection, &gd_old); 
182:                     beta = (gnorm2 - ginner) / (gd - gd_old);
183:                     break;

185:                 case CG_DaiYuan:
186:                     VecDot(tao->gradient, tao->stepdirection, &gd); 
187:                     VecDot(cgP->G_old, tao->stepdirection, &gd_old); 
188:                     beta = gnorm2 / (gd - gd_old);
189:                     break;

191:                 default:
192:                     beta = 0.0;
193:                     break;
194:             }
195:         }
196:         
197:         /*  Compute the direction d=-g + beta*d */
198:         VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient); 
199:         
200:         /*  update initial steplength choice */
201:         delta = 1.0;
202:         delta = PetscMax(delta, cgP->delta_min);
203:         delta = PetscMin(delta, cgP->delta_max);
204:     }
205:     return(0);
206: }

210: static PetscErrorCode TaoSetUp_CG(TaoSolver tao)
211: {
212:     TAO_CG *cgP = (TAO_CG*)tao->data;
214:     
216:     if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
217:     if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection);  }
218:     if (!cgP->X_old) {VecDuplicate(tao->solution,&cgP->X_old); }
219:     if (!cgP->G_old) {VecDuplicate(tao->gradient,&cgP->G_old);  }

221:     return(0);
222: }


227: static PetscErrorCode TaoDestroy_CG(TaoSolver tao)
228: {
229:     TAO_CG *cgP = (TAO_CG*) tao->data;


234:     if (tao->setupcalled) {
235:       VecDestroy(&cgP->X_old); 
236:       VecDestroy(&cgP->G_old); 
237:     }
238:     TaoLineSearchDestroy(&tao->linesearch); 

240:     PetscFree(tao->data); 
241:     tao->data = PETSC_NULL;
242:         
243:     return(0);
244: }

248: static PetscErrorCode TaoSetFromOptions_CG(TaoSolver tao)
249: {
250:     TAO_CG *cgP = (TAO_CG*)tao->data;
252:     
254:     TaoLineSearchSetFromOptions(tao->linesearch); 
255:     PetscOptionsHead("Nonlinear Conjugate Gradient method for unconstrained optimization"); 
256:     PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,0);
257:     PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, 0); 
258:     PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->
259: delta_min,&cgP->delta_min,0); 
260:     PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->
261: delta_max,&cgP->delta_max,0); 
262:     PetscOptionsTail(); 
263:     return(0);
264: }       

268: static PetscErrorCode TaoView_CG(TaoSolver tao, PetscViewer viewer)
269: {
270:     PetscBool isascii;
271:     TAO_CG *cgP = (TAO_CG*)tao->data;

275:     PetscTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii); 
276:     if (isascii) {
277:         PetscViewerASCIIPushTab(viewer); 
278:         PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]); 
279:         PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps); 
280:         ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps); 
281:         PetscViewerASCIIPopTab(viewer); 
282:     } else {
283:       SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO CG",((PetscObject)viewer)->type_name);
284:     }
285:     return(0);
286: }


292: PetscErrorCode TaoCreate_CG(TaoSolver tao)
293: {
294:     TAO_CG *cgP;
295:     const char *morethuente_type = TAOLINESEARCH_MT;
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;
303:     
304:     tao->max_it = 2000;
305:     tao->max_funcs = 4000;
306:     tao->fatol = 1e-4;
307:     tao->frtol = 1e-4;

309:     /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
310:     /*  method.  In particular, gtol should be less that 0.5; the value used in  */
311:     /*  Nocedal and Wright is 0.10.  We use the default values for the  */
312:     /*  linesearch because it seems to work better. */
313:     TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); 
314:     TaoLineSearchSetType(tao->linesearch, morethuente_type); 
315:     TaoLineSearchUseTaoSolverRoutines(tao->linesearch, tao); 
316:     
317:     
318:     PetscNewLog(tao,TAO_CG,&cgP); 
319:     tao->data = (void*)cgP;
320:     cgP->eta = 0.1;
321:     cgP->delta_min = 1e-7;
322:     cgP->delta_max = 100;
323:     cgP->cg_type = CG_PolakRibierePlus;

325:     return(0);
326: }


330:             

332:                     
333: