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