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