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