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