Actual source code: taocg.c
petsc-3.6.1 2015-08-06
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;
25: if (tao->XL || tao->XU || tao->ops->computebounds) {
26: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");
27: }
29: /* Check convergence criteria */
30: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
31: VecNorm(tao->gradient,NORM_2,&gnorm);
32: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
34: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);
35: if (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: /* Save the current gradient information */
58: f_old = f;
59: gnorm2_old = gnorm2;
60: VecCopy(tao->solution, cgP->X_old);
61: VecCopy(tao->gradient, cgP->G_old);
62: VecDot(tao->gradient, tao->stepdirection, &gd);
63: if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
64: ++cgP->ngradsteps;
65: if (f != 0.0) {
66: delta = 2.0*PetscAbsScalar(f) / gnorm2;
67: delta = PetscMax(delta,cgP->delta_min);
68: delta = PetscMin(delta,cgP->delta_max);
69: } else {
70: delta = 2.0 / gnorm2;
71: delta = PetscMax(delta,cgP->delta_min);
72: delta = PetscMin(delta,cgP->delta_max);
73: }
75: VecCopy(tao->gradient, tao->stepdirection);
76: VecScale(tao->stepdirection, -1.0);
77: }
79: /* Search direction for improving point */
80: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
81: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
82: TaoAddLineSearchCounts(tao);
83: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
84: /* Linesearch failed */
85: /* Reset factors and use scaled gradient step */
86: ++cgP->nresetsteps;
87: f = f_old;
88: gnorm2 = gnorm2_old;
89: VecCopy(cgP->X_old, tao->solution);
90: VecCopy(cgP->G_old, tao->gradient);
92: if (f != 0.0) {
93: delta = 2.0*PetscAbsScalar(f) / gnorm2;
94: delta = PetscMax(delta,cgP->delta_min);
95: delta = PetscMin(delta,cgP->delta_max);
96: } else {
97: delta = 2.0 / gnorm2;
98: delta = PetscMax(delta,cgP->delta_min);
99: delta = PetscMin(delta,cgP->delta_max);
100: }
102: VecCopy(tao->gradient, tao->stepdirection);
103: VecScale(tao->stepdirection, -1.0);
105: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
106: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
107: TaoAddLineSearchCounts(tao);
109: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
110: /* Linesearch failed again */
111: /* switch to unscaled gradient */
112: f = f_old;
113: gnorm2 = gnorm2_old;
114: VecCopy(cgP->X_old, tao->solution);
115: VecCopy(cgP->G_old, tao->gradient);
116: delta = 1.0;
117: VecCopy(tao->solution, tao->stepdirection);
118: VecScale(tao->stepdirection, -1.0);
120: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
121: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
122: TaoAddLineSearchCounts(tao);
123: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
125: /* Line search failed for last time -- give up */
126: f = f_old;
127: gnorm2 = gnorm2_old;
128: VecCopy(cgP->X_old, tao->solution);
129: VecCopy(cgP->G_old, tao->gradient);
130: step = 0.0;
131: reason = TAO_DIVERGED_LS_FAILURE;
132: tao->reason = TAO_DIVERGED_LS_FAILURE;
133: }
134: }
135: }
137: /* Check for bad value */
138: VecNorm(tao->gradient,NORM_2,&gnorm);
139: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User-provided compute function generated Inf or NaN");
141: /* Check for termination */
142: gnorm2 =gnorm * gnorm;
143: tao->niter++;
144: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);
145: if (reason != TAO_CONTINUE_ITERATING) {
146: break;
147: }
149: /* Check for restart condition */
150: VecDot(tao->gradient, cgP->G_old, &ginner);
151: if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
152: /* Gradients far from orthognal; use steepest descent direction */
153: beta = 0.0;
154: } else {
155: /* Gradients close to orthogonal; use conjugate gradient formula */
156: switch (cgP->cg_type) {
157: case CG_FletcherReeves:
158: beta = gnorm2 / gnorm2_old;
159: break;
161: case CG_PolakRibiere:
162: beta = (gnorm2 - ginner) / gnorm2_old;
163: break;
165: case CG_PolakRibierePlus:
166: beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
167: break;
169: case CG_HestenesStiefel:
170: VecDot(tao->gradient, tao->stepdirection, &gd);
171: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
172: beta = (gnorm2 - ginner) / (gd - gd_old);
173: break;
175: case CG_DaiYuan:
176: VecDot(tao->gradient, tao->stepdirection, &gd);
177: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
178: beta = gnorm2 / (gd - gd_old);
179: break;
181: default:
182: beta = 0.0;
183: break;
184: }
185: }
187: /* Compute the direction d=-g + beta*d */
188: VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);
190: /* update initial steplength choice */
191: delta = 1.0;
192: delta = PetscMax(delta, cgP->delta_min);
193: delta = PetscMin(delta, cgP->delta_max);
194: }
195: return(0);
196: }
200: static PetscErrorCode TaoSetUp_CG(Tao tao)
201: {
202: TAO_CG *cgP = (TAO_CG*)tao->data;
206: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);}
207: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
208: if (!cgP->X_old) {VecDuplicate(tao->solution,&cgP->X_old);}
209: if (!cgP->G_old) {VecDuplicate(tao->gradient,&cgP->G_old); }
210: return(0);
211: }
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: }
232: static PetscErrorCode TaoSetFromOptions_CG(PetscOptions *PetscOptionsObject,Tao tao)
233: {
234: TAO_CG *cgP = (TAO_CG*)tao->data;
238: TaoLineSearchSetFromOptions(tao->linesearch);
239: PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");
240: PetscOptionsReal("-tao_cg_eta","restart tolerance", "", cgP->eta,&cgP->eta,NULL);
241: PetscOptionsEList("-tao_cg_type","cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type,NULL);
242: PetscOptionsReal("-tao_cg_delta_min","minimum delta value", "", cgP->delta_min,&cgP->delta_min,NULL);
243: PetscOptionsReal("-tao_cg_delta_max","maximum delta value", "", cgP->delta_max,&cgP->delta_max,NULL);
244: PetscOptionsTail();
245: return(0);
246: }
250: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
251: {
252: PetscBool isascii;
253: TAO_CG *cgP = (TAO_CG*)tao->data;
257: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
258: if (isascii) {
259: PetscViewerASCIIPushTab(viewer);
260: PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
261: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);
262: ierr= PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);
263: PetscViewerASCIIPopTab(viewer);
264: }
265: return(0);
266: }
268: /*MC
269: TAOCG - Nonlinear conjugate gradient method is an extension of the
270: nonlinear conjugate gradient solver for nonlinear optimization.
272: Options Database Keys:
273: + -tao_cg_eta <r> - restart tolerance
274: . -tao_cg_type <taocg_type> - cg formula
275: . -tao_cg_delta_min <r> - minimum delta value
276: - -tao_cg_delta_max <r> - maximum delta value
278: Notes:
279: CG formulas are:
280: "fr" - Fletcher-Reeves
281: "pr" - Polak-Ribiere
282: "prp" - Polak-Ribiere-Plus
283: "hs" - Hestenes-Steifel
284: "dy" - Dai-Yuan
285: Level: beginner
286: M*/
291: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
292: {
293: TAO_CG *cgP;
294: const char *morethuente_type = TAOLINESEARCHMT;
298: tao->ops->setup = TaoSetUp_CG;
299: tao->ops->solve = TaoSolve_CG;
300: tao->ops->view = TaoView_CG;
301: tao->ops->setfromoptions = TaoSetFromOptions_CG;
302: tao->ops->destroy = TaoDestroy_CG;
304: /* Override default settings (unless already changed) */
305: if (!tao->max_it_changed) tao->max_it = 2000;
306: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
307: if (!tao->fatol_changed) tao->fatol = 1e-4;
308: if (!tao->frtol_changed) tao->frtol = 1e-4;
310: /* Note: nondefault values should be used for nonlinear conjugate gradient */
311: /* method. In particular, gtol should be less that 0.5; the value used in */
312: /* Nocedal and Wright is 0.10. We use the default values for the */
313: /* linesearch because it seems to work better. */
314: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
315: TaoLineSearchSetType(tao->linesearch, morethuente_type);
316: TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
317: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
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: }