Actual source code: taocg.c
petsc-3.5.4 2015-05-23
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"};
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;
26: if (tao->XL || tao->XU || tao->ops->computebounds) {
27: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");
28: }
30: /* Check convergence criteria */
31: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
32: VecNorm(tao->gradient,NORM_2,&gnorm);
33: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
35: TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);
36: if (reason != TAO_CONTINUE_ITERATING) return(0);
38: /* Set initial direction to -gradient */
39: VecCopy(tao->gradient, tao->stepdirection);
40: VecScale(tao->stepdirection, -1.0);
41: gnorm2 = gnorm*gnorm;
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;
57: while (1) {
58: /* Save the current gradient information */
59: f_old = f;
60: gnorm2_old = gnorm2;
61: VecCopy(tao->solution, cgP->X_old);
62: VecCopy(tao->gradient, cgP->G_old);
63: VecDot(tao->gradient, tao->stepdirection, &gd);
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: }
76: VecCopy(tao->gradient, tao->stepdirection);
77: VecScale(tao->stepdirection, -1.0);
78: }
80: /* Search direction for improving point */
81: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
82: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
83: TaoAddLineSearchCounts(tao);
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: VecCopy(cgP->X_old, tao->solution);
91: VecCopy(cgP->G_old, tao->gradient);
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: }
103: VecCopy(tao->gradient, tao->stepdirection);
104: VecScale(tao->stepdirection, -1.0);
106: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
107: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
108: TaoAddLineSearchCounts(tao);
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: VecCopy(cgP->X_old, tao->solution);
116: VecCopy(cgP->G_old, tao->gradient);
117: delta = 1.0;
118: VecCopy(tao->solution, tao->stepdirection);
119: VecScale(tao->stepdirection, -1.0);
121: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
122: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
123: TaoAddLineSearchCounts(tao);
124: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
126: /* Line search failed for last time -- give up */
127: f = f_old;
128: gnorm2 = gnorm2_old;
129: VecCopy(cgP->X_old, tao->solution);
130: VecCopy(cgP->G_old, tao->gradient);
131: step = 0.0;
132: reason = TAO_DIVERGED_LS_FAILURE;
133: tao->reason = TAO_DIVERGED_LS_FAILURE;
134: }
135: }
136: }
138: /* Check for bad value */
139: VecNorm(tao->gradient,NORM_2,&gnorm);
140: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN");
142: /* Check for termination */
143: gnorm2 =gnorm * gnorm;
144: iter++;
145: TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);
146: if (reason != TAO_CONTINUE_ITERATING) {
147: break;
148: }
150: /* Check for restart condition */
151: VecDot(tao->gradient, cgP->G_old, &ginner);
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;
162: case CG_PolakRibiere:
163: beta = (gnorm2 - ginner) / gnorm2_old;
164: break;
166: case CG_PolakRibierePlus:
167: beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
168: break;
170: case CG_HestenesStiefel:
171: VecDot(tao->gradient, tao->stepdirection, &gd);
172: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
173: beta = (gnorm2 - ginner) / (gd - gd_old);
174: break;
176: case CG_DaiYuan:
177: VecDot(tao->gradient, tao->stepdirection, &gd);
178: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
179: beta = gnorm2 / (gd - gd_old);
180: break;
182: default:
183: beta = 0.0;
184: break;
185: }
186: }
188: /* Compute the direction d=-g + beta*d */
189: VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);
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: return(0);
197: }
201: static PetscErrorCode TaoSetUp_CG(Tao tao)
202: {
203: TAO_CG *cgP = (TAO_CG*)tao->data;
207: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);}
208: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
209: if (!cgP->X_old) {VecDuplicate(tao->solution,&cgP->X_old);}
210: if (!cgP->G_old) {VecDuplicate(tao->gradient,&cgP->G_old); }
211: return(0);
212: }
216: static PetscErrorCode TaoDestroy_CG(Tao tao)
217: {
218: TAO_CG *cgP = (TAO_CG*) tao->data;
222: if (tao->setupcalled) {
223: VecDestroy(&cgP->X_old);
224: VecDestroy(&cgP->G_old);
225: }
226: TaoLineSearchDestroy(&tao->linesearch);
227: PetscFree(tao->data);
228: return(0);
229: }
233: static PetscErrorCode TaoSetFromOptions_CG(Tao tao)
234: {
235: TAO_CG *cgP = (TAO_CG*)tao->data;
239: TaoLineSearchSetFromOptions(tao->linesearch);
240: PetscOptionsHead("Nonlinear Conjugate Gradient method for unconstrained optimization");
241: PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,0);
242: PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, 0);
243: PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,0);
244: PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,0);
245: PetscOptionsTail();
246: return(0);
247: }
251: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
252: {
253: PetscBool isascii;
254: TAO_CG *cgP = (TAO_CG*)tao->data;
258: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
259: if (isascii) {
260: PetscViewerASCIIPushTab(viewer);
261: PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
262: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);
263: ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);
264: PetscViewerASCIIPopTab(viewer);
265: }
266: return(0);
267: }
269: /*MC
270: TAOCG - Nonlinear conjugate gradient method is an extension of the
271: nonlinear conjugate gradient solver for nonlinear optimization.
273: Options Database Keys:
274: + -tao_cg_eta <r> - restart tolerance
275: . -tao_cg_type <taocg_type> - cg formula
276: . -tao_cg_delta_min <r> - minimum delta value
277: - -tao_cg_delta_max <r> - maximum delta value
279: Notes:
280: CG formulas are:
281: "fr" - Fletcher-Reeves
282: "pr" - Polak-Ribiere
283: "prp" - Polak-Ribiere-Plus
284: "hs" - Hestenes-Steifel
285: "dy" - Dai-Yuan
286: Level: beginner
287: M*/
290: EXTERN_C_BEGIN
293: PetscErrorCode TaoCreate_CG(Tao tao)
294: {
295: TAO_CG *cgP;
296: const char *morethuente_type = TAOLINESEARCHMT;
300: tao->ops->setup = TaoSetUp_CG;
301: tao->ops->solve = TaoSolve_CG;
302: tao->ops->view = TaoView_CG;
303: tao->ops->setfromoptions = TaoSetFromOptions_CG;
304: tao->ops->destroy = TaoDestroy_CG;
306: tao->max_it = 2000;
307: tao->max_funcs = 4000;
308: tao->fatol = 1e-4;
309: tao->frtol = 1e-4;
311: /* Note: nondefault values should be used for nonlinear conjugate gradient */
312: /* method. In particular, gtol should be less that 0.5; the value used in */
313: /* Nocedal and Wright is 0.10. We use the default values for the */
314: /* linesearch because it seems to work better. */
315: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
316: TaoLineSearchSetType(tao->linesearch, morethuente_type);
317: TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
319: PetscNewLog(tao,&cgP);
320: tao->data = (void*)cgP;
321: cgP->eta = 0.1;
322: cgP->delta_min = 1e-7;
323: cgP->delta_max = 100;
324: cgP->cg_type = CG_PolakRibierePlus;
325: return(0);
326: }
327: EXTERN_C_END