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