xref: /petsc/src/tao/unconstrained/impls/cg/taocg.c (revision b8a3fda2eec00bbbfad0cf0db17d1d1f04b07f29) 
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          ierr = VecCopy(cgP->X_old, tao->solution);CHKERRQ(ierr);
114          ierr = VecCopy(cgP->G_old, tao->gradient);CHKERRQ(ierr);
115          delta = 1.0;
116          ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr);
117          ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
118 
119          ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);CHKERRQ(ierr);
120          ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
121          ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
122          if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
123 
124            /*  Line search failed for last time -- give up */
125            f = f_old;
126            ierr = VecCopy(cgP->X_old, tao->solution);CHKERRQ(ierr);
127            ierr = VecCopy(cgP->G_old, tao->gradient);CHKERRQ(ierr);
128            step = 0.0;
129            reason = TAO_DIVERGED_LS_FAILURE;
130            tao->reason = TAO_DIVERGED_LS_FAILURE;
131          }
132        }
133      }
134 
135      /*  Check for bad value */
136      ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
137      if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN");
138 
139      /*  Check for termination */
140      gnorm2 =gnorm * gnorm;
141      tao->niter++;
142      ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);CHKERRQ(ierr);
143      if (reason != TAO_CONTINUE_ITERATING) {
144        break;
145      }
146 
147      /*  Check for restart condition */
148      ierr = VecDot(tao->gradient, cgP->G_old, &ginner);CHKERRQ(ierr);
149      if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
150        /*  Gradients far from orthognal; 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;
158 
159        case CG_PolakRibiere:
160          beta = (gnorm2 - ginner) / gnorm2_old;
161          break;
162 
163        case CG_PolakRibierePlus:
164          beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
165          break;
166 
167        case CG_HestenesStiefel:
168          ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
169          ierr = VecDot(cgP->G_old, tao->stepdirection, &gd_old);CHKERRQ(ierr);
170          beta = (gnorm2 - ginner) / (gd - gd_old);
171          break;
172 
173        case CG_DaiYuan:
174          ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
175          ierr = VecDot(cgP->G_old, tao->stepdirection, &gd_old);CHKERRQ(ierr);
176          beta = gnorm2 / (gd - gd_old);
177          break;
178 
179        default:
180          beta = 0.0;
181          break;
182        }
183      }
184 
185      /*  Compute the direction d=-g + beta*d */
186      ierr = VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);CHKERRQ(ierr);
187 
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(0);
194  }
195 
196  #undef __FUNCT__
197  #define __FUNCT__ "TaoSetUp_CG"
198  static PetscErrorCode TaoSetUp_CG(Tao tao)
199  {
200    TAO_CG         *cgP = (TAO_CG*)tao->data;
201    PetscErrorCode ierr;
202 
203    PetscFunctionBegin;
204    if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);}
205    if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr); }
206    if (!cgP->X_old) {ierr = VecDuplicate(tao->solution,&cgP->X_old);CHKERRQ(ierr);}
207    if (!cgP->G_old) {ierr = VecDuplicate(tao->gradient,&cgP->G_old);CHKERRQ(ierr); }
208     PetscFunctionReturn(0);
209  }
210 
211  #undef __FUNCT__
212  #define __FUNCT__ "TaoDestroy_CG"
213  static PetscErrorCode TaoDestroy_CG(Tao tao)
214  {
215    TAO_CG         *cgP = (TAO_CG*) tao->data;
216    PetscErrorCode ierr;
217 
218    PetscFunctionBegin;
219    if (tao->setupcalled) {
220      ierr = VecDestroy(&cgP->X_old);CHKERRQ(ierr);
221      ierr = VecDestroy(&cgP->G_old);CHKERRQ(ierr);
222    }
223    ierr = TaoLineSearchDestroy(&tao->linesearch);CHKERRQ(ierr);
224    ierr = PetscFree(tao->data);CHKERRQ(ierr);
225    PetscFunctionReturn(0);
226  }
227 
228  #undef __FUNCT__
229  #define __FUNCT__ "TaoSetFromOptions_CG"
230 static PetscErrorCode TaoSetFromOptions_CG(PetscOptionItems *PetscOptionsObject,Tao tao)
231  {
232     TAO_CG         *cgP = (TAO_CG*)tao->data;
233     PetscErrorCode ierr;
234 
235     PetscFunctionBegin;
236     ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
237     ierr = PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");CHKERRQ(ierr);
238     ierr = PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,NULL);CHKERRQ(ierr);
239     ierr = PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type,NULL);CHKERRQ(ierr);
240     ierr = PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,NULL);CHKERRQ(ierr);
241     ierr = PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,NULL);CHKERRQ(ierr);
242    ierr = PetscOptionsTail();CHKERRQ(ierr);
243    PetscFunctionReturn(0);
244 }
245 
246 #undef __FUNCT__
247 #define __FUNCT__ "TaoView_CG"
248 static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
249 {
250   PetscBool      isascii;
251   TAO_CG         *cgP = (TAO_CG*)tao->data;
252   PetscErrorCode ierr;
253 
254   PetscFunctionBegin;
255   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
256   if (isascii) {
257     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
258     ierr = PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);CHKERRQ(ierr);
259     ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);CHKERRQ(ierr);
260     ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);CHKERRQ(ierr);
261     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
262   }
263   PetscFunctionReturn(0);
264 }
265 
266 /*MC
267      TAOCG -   Nonlinear conjugate gradient method is an extension of the
268 nonlinear conjugate gradient solver for nonlinear optimization.
269 
270    Options Database Keys:
271 +      -tao_cg_eta <r> - restart tolerance
272 .      -tao_cg_type <taocg_type> - cg formula
273 .      -tao_cg_delta_min <r> - minimum delta value
274 -      -tao_cg_delta_max <r> - maximum delta value
275 
276   Notes:
277      CG formulas are:
278          "fr" - Fletcher-Reeves
279          "pr" - Polak-Ribiere
280          "prp" - Polak-Ribiere-Plus
281          "hs" - Hestenes-Steifel
282          "dy" - Dai-Yuan
283   Level: beginner
284 M*/
285 
286 
287 #undef __FUNCT__
288 #define __FUNCT__ "TaoCreate_CG"
289 PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
290 {
291   TAO_CG         *cgP;
292   const char     *morethuente_type = TAOLINESEARCHMT;
293   PetscErrorCode ierr;
294 
295   PetscFunctionBegin;
296   tao->ops->setup = TaoSetUp_CG;
297   tao->ops->solve = TaoSolve_CG;
298   tao->ops->view = TaoView_CG;
299   tao->ops->setfromoptions = TaoSetFromOptions_CG;
300   tao->ops->destroy = TaoDestroy_CG;
301 
302   /* Override default settings (unless already changed) */
303   if (!tao->max_it_changed) tao->max_it = 2000;
304   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
305 
306   /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
307   /*  method.  In particular, gtol should be less that 0.5; the value used in  */
308   /*  Nocedal and Wright is 0.10.  We use the default values for the  */
309   /*  linesearch because it seems to work better. */
310   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
311   ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr);
312   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch, tao);CHKERRQ(ierr);
313   ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
314 
315   ierr = PetscNewLog(tao,&cgP);CHKERRQ(ierr);
316   tao->data = (void*)cgP;
317   cgP->eta = 0.1;
318   cgP->delta_min = 1e-7;
319   cgP->delta_max = 100;
320   cgP->cg_type = CG_PolakRibierePlus;
321   PetscFunctionReturn(0);
322 }
323