Actual source code: taocg.c
petsc-3.12.5 2020-03-29
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/unconstrained/impls/cg/taocg.h>
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
11: static const char *CG_Table[64] = {"fr", "pr", "prp", "hs", "dy"};
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;
22: if (tao->XL || tao->XU || tao->ops->computebounds) {
23: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");
24: }
26: /* Check convergence criteria */
27: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
28: VecNorm(tao->gradient,NORM_2,&gnorm);
29: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
30:
31: tao->reason = TAO_CONTINUE_ITERATING;
32: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
33: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
34: (*tao->ops->convergencetest)(tao,tao->cnvP);
35: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
37: /* Set initial direction to -gradient */
38: VecCopy(tao->gradient, tao->stepdirection);
39: VecScale(tao->stepdirection, -1.0);
40: gnorm2 = gnorm*gnorm;
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;
56: while (1) {
57: /* Call general purpose update function */
58: if (tao->ops->update) {
59: (*tao->ops->update)(tao, tao->niter, tao->user_update);
60: }
61:
62: /* Save the current gradient information */
63: f_old = f;
64: gnorm2_old = gnorm2;
65: VecCopy(tao->solution, cgP->X_old);
66: VecCopy(tao->gradient, cgP->G_old);
67: VecDot(tao->gradient, tao->stepdirection, &gd);
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: }
80: VecCopy(tao->gradient, tao->stepdirection);
81: VecScale(tao->stepdirection, -1.0);
82: }
84: /* Search direction for improving point */
85: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
86: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
87: TaoAddLineSearchCounts(tao);
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: VecCopy(cgP->X_old, tao->solution);
95: VecCopy(cgP->G_old, tao->gradient);
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: }
107: VecCopy(tao->gradient, tao->stepdirection);
108: VecScale(tao->stepdirection, -1.0);
110: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
111: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
112: TaoAddLineSearchCounts(tao);
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: VecCopy(cgP->X_old, tao->solution);
119: VecCopy(cgP->G_old, tao->gradient);
120: delta = 1.0;
121: VecCopy(tao->solution, tao->stepdirection);
122: VecScale(tao->stepdirection, -1.0);
124: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
125: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
126: TaoAddLineSearchCounts(tao);
127: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
129: /* Line search failed for last time -- give up */
130: f = f_old;
131: VecCopy(cgP->X_old, tao->solution);
132: VecCopy(cgP->G_old, tao->gradient);
133: step = 0.0;
134: tao->reason = TAO_DIVERGED_LS_FAILURE;
135: }
136: }
137: }
139: /* Check for bad value */
140: VecNorm(tao->gradient,NORM_2,&gnorm);
141: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN");
143: /* Check for termination */
144: gnorm2 =gnorm * gnorm;
145: tao->niter++;
146: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
147: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
148: (*tao->ops->convergencetest)(tao,tao->cnvP);
149: if (tao->reason != TAO_CONTINUE_ITERATING) {
150: break;
151: }
153: /* Check for restart condition */
154: VecDot(tao->gradient, cgP->G_old, &ginner);
155: if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
156: /* Gradients far from orthognal; 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;
165: case CG_PolakRibiere:
166: beta = (gnorm2 - ginner) / gnorm2_old;
167: break;
169: case CG_PolakRibierePlus:
170: beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
171: break;
173: case CG_HestenesStiefel:
174: VecDot(tao->gradient, tao->stepdirection, &gd);
175: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
176: beta = (gnorm2 - ginner) / (gd - gd_old);
177: break;
179: case CG_DaiYuan:
180: VecDot(tao->gradient, tao->stepdirection, &gd);
181: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
182: beta = gnorm2 / (gd - gd_old);
183: break;
185: default:
186: beta = 0.0;
187: break;
188: }
189: }
191: /* Compute the direction d=-g + beta*d */
192: VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);
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: return(0);
200: }
202: static PetscErrorCode TaoSetUp_CG(Tao tao)
203: {
204: TAO_CG *cgP = (TAO_CG*)tao->data;
208: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);}
209: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
210: if (!cgP->X_old) {VecDuplicate(tao->solution,&cgP->X_old);}
211: if (!cgP->G_old) {VecDuplicate(tao->gradient,&cgP->G_old); }
212: return(0);
213: }
215: static PetscErrorCode TaoDestroy_CG(Tao tao)
216: {
217: TAO_CG *cgP = (TAO_CG*) tao->data;
221: if (tao->setupcalled) {
222: VecDestroy(&cgP->X_old);
223: VecDestroy(&cgP->G_old);
224: }
225: TaoLineSearchDestroy(&tao->linesearch);
226: PetscFree(tao->data);
227: return(0);
228: }
230: static PetscErrorCode TaoSetFromOptions_CG(PetscOptionItems *PetscOptionsObject,Tao tao)
231: {
232: TAO_CG *cgP = (TAO_CG*)tao->data;
236: TaoLineSearchSetFromOptions(tao->linesearch);
237: PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");
238: PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,NULL);
239: PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type,NULL);
240: PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,NULL);
241: PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,NULL);
242: PetscOptionsTail();
243: return(0);
244: }
246: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
247: {
248: PetscBool isascii;
249: TAO_CG *cgP = (TAO_CG*)tao->data;
253: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
254: if (isascii) {
255: PetscViewerASCIIPushTab(viewer);
256: PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
257: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);
258: ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);
259: PetscViewerASCIIPopTab(viewer);
260: }
261: return(0);
262: }
264: /*MC
265: TAOCG - Nonlinear conjugate gradient method is an extension of the
266: nonlinear conjugate gradient solver for nonlinear optimization.
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
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*/
285: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
286: {
287: TAO_CG *cgP;
288: const char *morethuente_type = TAOLINESEARCHMT;
292: tao->ops->setup = TaoSetUp_CG;
293: tao->ops->solve = TaoSolve_CG;
294: tao->ops->view = TaoView_CG;
295: tao->ops->setfromoptions = TaoSetFromOptions_CG;
296: tao->ops->destroy = TaoDestroy_CG;
298: /* Override default settings (unless already changed) */
299: if (!tao->max_it_changed) tao->max_it = 2000;
300: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
302: /* Note: nondefault values should be used for nonlinear conjugate gradient */
303: /* method. In particular, gtol should be less that 0.5; the value used in */
304: /* Nocedal and Wright is 0.10. We use the default values for the */
305: /* linesearch because it seems to work better. */
306: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
307: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
308: TaoLineSearchSetType(tao->linesearch, morethuente_type);
309: TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
310: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
312: PetscNewLog(tao,&cgP);
313: tao->data = (void*)cgP;
314: cgP->eta = 0.1;
315: cgP->delta_min = 1e-7;
316: cgP->delta_max = 100;
317: cgP->cg_type = CG_PolakRibierePlus;
318: return(0);
319: }