xref: /petsc/src/tao/unconstrained/impls/cg/taocg.c (revision bebe2cf65d55febe21a5af8db2bd2e168caaa2e7) 
1 #include <petsctaolinesearch.h>
2 #include <../src/tao/unconstrained/impls/cg/taocg.h>
3 
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
10 
11  static const char *CG_Table[64] = {"fr", "pr", "prp", "hs", "dy"};
12 
13  #undef __FUNCT__
14  #define __FUNCT__ "TaoSolve_CG"
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;
23 
24    PetscFunctionBegin;
25    if (tao->XL || tao->XU || tao->ops->computebounds) {
26      ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");CHKERRQ(ierr);
27    }
28 
29    /*  Check convergence criteria */
30    ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
31    ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
32    if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
33 
34    ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);CHKERRQ(ierr);
35    if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
36 
37    /*  Set initial direction to -gradient */
38    ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
39    ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
40    gnorm2 = gnorm*gnorm;
41 
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;
55 
56    while (1) {
57      /*  Save the current gradient information */
58      f_old = f;
59      gnorm2_old = gnorm2;
60      ierr = VecCopy(tao->solution, cgP->X_old);CHKERRQ(ierr);
61      ierr = VecCopy(tao->gradient, cgP->G_old);CHKERRQ(ierr);
62      ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
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        }
74 
75        ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
76        ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
77      }
78 
79      /*  Search direction for improving point */
80      ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);CHKERRQ(ierr);
81      ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
82      ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
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        ierr = VecCopy(cgP->X_old, tao->solution);CHKERRQ(ierr);
90        ierr = VecCopy(cgP->G_old, tao->gradient);CHKERRQ(ierr);
91 
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        }
101 
102        ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr);
103        ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
104 
105        ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);CHKERRQ(ierr);
106        ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
107        ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
108 
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          gnorm2 = gnorm2_old;
114          ierr = VecCopy(cgP->X_old, tao->solution);CHKERRQ(ierr);
115          ierr = VecCopy(cgP->G_old, tao->gradient);CHKERRQ(ierr);
116          delta = 1.0;
117          ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr);
118          ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
119 
120          ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);CHKERRQ(ierr);
121          ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
122          ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
123          if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
124 
125            /*  Line search failed for last time -- give up */
126            f = f_old;
127            gnorm2 = gnorm2_old;
128            ierr = VecCopy(cgP->X_old, tao->solution);CHKERRQ(ierr);
129            ierr = VecCopy(cgP->G_old, tao->gradient);CHKERRQ(ierr);
130            step = 0.0;
131            reason = TAO_DIVERGED_LS_FAILURE;
132            tao->reason = TAO_DIVERGED_LS_FAILURE;
133          }
134        }
135      }
136 
137      /*  Check for bad value */
138      ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
139      if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN");
140 
141      /*  Check for termination */
142      gnorm2 =gnorm * gnorm;
143      tao->niter++;
144      ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);CHKERRQ(ierr);
145      if (reason != TAO_CONTINUE_ITERATING) {
146        break;
147      }
148 
149      /*  Check for restart condition */
150      ierr = VecDot(tao->gradient, cgP->G_old, &ginner);CHKERRQ(ierr);
151      if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
152        /*  Gradients far from orthognal; use steepest descent direction */
153        beta = 0.0;
154      } else {
155        /*  Gradients close to orthogonal; use conjugate gradient formula */
156        switch (cgP->cg_type) {
157        case CG_FletcherReeves:
158          beta = gnorm2 / gnorm2_old;
159          break;
160 
161        case CG_PolakRibiere:
162          beta = (gnorm2 - ginner) / gnorm2_old;
163          break;
164 
165        case CG_PolakRibierePlus:
166          beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
167          break;
168 
169        case CG_HestenesStiefel:
170          ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
171          ierr = VecDot(cgP->G_old, tao->stepdirection, &gd_old);CHKERRQ(ierr);
172          beta = (gnorm2 - ginner) / (gd - gd_old);
173          break;
174 
175        case CG_DaiYuan:
176          ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
177          ierr = VecDot(cgP->G_old, tao->stepdirection, &gd_old);CHKERRQ(ierr);
178          beta = gnorm2 / (gd - gd_old);
179          break;
180 
181        default:
182          beta = 0.0;
183          break;
184        }
185      }
186 
187      /*  Compute the direction d=-g + beta*d */
188      ierr = VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);CHKERRQ(ierr);
189 
190      /*  update initial steplength choice */
191      delta = 1.0;
192      delta = PetscMax(delta, cgP->delta_min);
193      delta = PetscMin(delta, cgP->delta_max);
194    }
195    PetscFunctionReturn(0);
196  }
197 
198  #undef __FUNCT__
199  #define __FUNCT__ "TaoSetUp_CG"
200  static PetscErrorCode TaoSetUp_CG(Tao tao)
201  {
202    TAO_CG         *cgP = (TAO_CG*)tao->data;
203    PetscErrorCode ierr;
204 
205    PetscFunctionBegin;
206    if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);}
207    if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr); }
208    if (!cgP->X_old) {ierr = VecDuplicate(tao->solution,&cgP->X_old);CHKERRQ(ierr);}
209    if (!cgP->G_old) {ierr = VecDuplicate(tao->gradient,&cgP->G_old);CHKERRQ(ierr); }
210     PetscFunctionReturn(0);
211  }
212 
213  #undef __FUNCT__
214  #define __FUNCT__ "TaoDestroy_CG"
215  static PetscErrorCode TaoDestroy_CG(Tao tao)
216  {
217    TAO_CG         *cgP = (TAO_CG*) tao->data;
218    PetscErrorCode ierr;
219 
220    PetscFunctionBegin;
221    if (tao->setupcalled) {
222      ierr = VecDestroy(&cgP->X_old);CHKERRQ(ierr);
223      ierr = VecDestroy(&cgP->G_old);CHKERRQ(ierr);
224    }
225    ierr = TaoLineSearchDestroy(&tao->linesearch);CHKERRQ(ierr);
226    ierr = PetscFree(tao->data);CHKERRQ(ierr);
227    PetscFunctionReturn(0);
228  }
229 
230  #undef __FUNCT__
231  #define __FUNCT__ "TaoSetFromOptions_CG"
232 static PetscErrorCode TaoSetFromOptions_CG(PetscOptions *PetscOptionsObject,Tao tao)
233  {
234     TAO_CG         *cgP = (TAO_CG*)tao->data;
235     PetscErrorCode ierr;
236 
237     PetscFunctionBegin;
238     ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
239     ierr = PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");CHKERRQ(ierr);
240     ierr = PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,NULL);CHKERRQ(ierr);
241     ierr = PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type,NULL);CHKERRQ(ierr);
242     ierr = PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,NULL);CHKERRQ(ierr);
243     ierr = PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,NULL);CHKERRQ(ierr);
244    ierr = PetscOptionsTail();CHKERRQ(ierr);
245    PetscFunctionReturn(0);
246 }
247 
248 #undef __FUNCT__
249 #define __FUNCT__ "TaoView_CG"
250 static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
251 {
252   PetscBool      isascii;
253   TAO_CG         *cgP = (TAO_CG*)tao->data;
254   PetscErrorCode ierr;
255 
256   PetscFunctionBegin;
257   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
258   if (isascii) {
259     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
260     ierr = PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);CHKERRQ(ierr);
261     ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);CHKERRQ(ierr);
262     ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);CHKERRQ(ierr);
263     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
264   }
265   PetscFunctionReturn(0);
266 }
267 
268 /*MC
269      TAOCG -   Nonlinear conjugate gradient method is an extension of the
270 nonlinear conjugate gradient solver for nonlinear optimization.
271 
272    Options Database Keys:
273 +      -tao_cg_eta <r> - restart tolerance
274 .      -tao_cg_type <taocg_type> - cg formula
275 .      -tao_cg_delta_min <r> - minimum delta value
276 -      -tao_cg_delta_max <r> - maximum delta value
277 
278   Notes:
279      CG formulas are:
280          "fr" - Fletcher-Reeves
281          "pr" - Polak-Ribiere
282          "prp" - Polak-Ribiere-Plus
283          "hs" - Hestenes-Steifel
284          "dy" - Dai-Yuan
285   Level: beginner
286 M*/
287 
288 
289 #undef __FUNCT__
290 #define __FUNCT__ "TaoCreate_CG"
291 PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
292 {
293   TAO_CG         *cgP;
294   const char     *morethuente_type = TAOLINESEARCHMT;
295   PetscErrorCode ierr;
296 
297   PetscFunctionBegin;
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   /* Override default settings (unless already changed) */
305   if (!tao->max_it_changed) tao->max_it = 2000;
306   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
307   if (!tao->fatol_changed) tao->fatol = 1e-4;
308   if (!tao->frtol_changed) tao->frtol = 1e-4;
309 
310   /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
311   /*  method.  In particular, gtol should be less that 0.5; the value used in  */
312   /*  Nocedal and Wright is 0.10.  We use the default values for the  */
313   /*  linesearch because it seems to work better. */
314   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
315   ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr);
316   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch, tao);CHKERRQ(ierr);
317   ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
318 
319   ierr = PetscNewLog(tao,&cgP);CHKERRQ(ierr);
320   tao->data = (void*)cgP;
321   cgP->eta = 0.1;
322   cgP->delta_min = 1e-7;
323   cgP->delta_max = 100;
324   cgP->cg_type = CG_PolakRibierePlus;
325   PetscFunctionReturn(0);
326 }
327