xref: /petsc/src/tao/unconstrained/impls/cg/taocg.c (revision a3fa217bfefef3d54995bb82b2b5ea8a3c3e9019) 
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,0);CHKERRQ(ierr);
242    ierr = PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, 0);CHKERRQ(ierr);
243    ierr = PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,0);CHKERRQ(ierr);
244    ierr = PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,0);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 
270 EXTERN_C_BEGIN
271 #undef __FUNCT__
272 #define __FUNCT__ "TaoCreate_CG"
273 PetscErrorCode TaoCreate_CG(Tao tao)
274 {
275   TAO_CG         *cgP;
276   const char     *morethuente_type = TAOLINESEARCHMT;
277   PetscErrorCode ierr;
278 
279   PetscFunctionBegin;
280   tao->ops->setup = TaoSetUp_CG;
281   tao->ops->solve = TaoSolve_CG;
282   tao->ops->view = TaoView_CG;
283   tao->ops->setfromoptions = TaoSetFromOptions_CG;
284   tao->ops->destroy = TaoDestroy_CG;
285 
286   tao->max_it = 2000;
287   tao->max_funcs = 4000;
288   tao->fatol = 1e-4;
289   tao->frtol = 1e-4;
290 
291   /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
292   /*  method.  In particular, gtol should be less that 0.5; the value used in  */
293   /*  Nocedal and Wright is 0.10.  We use the default values for the  */
294   /*  linesearch because it seems to work better. */
295   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
296   ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr);
297   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch, tao);CHKERRQ(ierr);
298 
299   ierr = PetscNewLog(tao,&cgP);CHKERRQ(ierr);
300   tao->data = (void*)cgP;
301   cgP->eta = 0.1;
302   cgP->delta_min = 1e-7;
303   cgP->delta_max = 100;
304   cgP->cg_type = CG_PolakRibierePlus;
305   PetscFunctionReturn(0);
306 }
307 EXTERN_C_END
308