Actual source code: taocg.c
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: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
17: PetscReal step=1.0,f,gnorm,gnorm2,delta,gd,ginner,beta;
18: PetscReal gd_old,gnorm2_old,f_old;
20: if (tao->XL || tao->XU || tao->ops->computebounds) {
21: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");
22: }
24: /* Check convergence criteria */
25: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
26: VecNorm(tao->gradient,NORM_2,&gnorm);
29: tao->reason = TAO_CONTINUE_ITERATING;
30: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
31: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
32: (*tao->ops->convergencetest)(tao,tao->cnvP);
33: if (tao->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: /* Call general purpose update function */
56: if (tao->ops->update) {
57: (*tao->ops->update)(tao, tao->niter, tao->user_update);
58: }
60: /* Save the current gradient information */
61: f_old = f;
62: gnorm2_old = gnorm2;
63: VecCopy(tao->solution, cgP->X_old);
64: VecCopy(tao->gradient, cgP->G_old);
65: VecDot(tao->gradient, tao->stepdirection, &gd);
66: if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
67: ++cgP->ngradsteps;
68: if (f != 0.0) {
69: delta = 2.0*PetscAbsScalar(f) / gnorm2;
70: delta = PetscMax(delta,cgP->delta_min);
71: delta = PetscMin(delta,cgP->delta_max);
72: } else {
73: delta = 2.0 / gnorm2;
74: delta = PetscMax(delta,cgP->delta_min);
75: delta = PetscMin(delta,cgP->delta_max);
76: }
78: VecCopy(tao->gradient, tao->stepdirection);
79: VecScale(tao->stepdirection, -1.0);
80: }
82: /* Search direction for improving point */
83: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
84: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
85: TaoAddLineSearchCounts(tao);
86: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
87: /* Linesearch failed */
88: /* Reset factors and use scaled gradient step */
89: ++cgP->nresetsteps;
90: f = f_old;
91: gnorm2 = gnorm2_old;
92: VecCopy(cgP->X_old, tao->solution);
93: VecCopy(cgP->G_old, tao->gradient);
95: if (f != 0.0) {
96: delta = 2.0*PetscAbsScalar(f) / gnorm2;
97: delta = PetscMax(delta,cgP->delta_min);
98: delta = PetscMin(delta,cgP->delta_max);
99: } else {
100: delta = 2.0 / gnorm2;
101: delta = PetscMax(delta,cgP->delta_min);
102: delta = PetscMin(delta,cgP->delta_max);
103: }
105: VecCopy(tao->gradient, tao->stepdirection);
106: VecScale(tao->stepdirection, -1.0);
108: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
109: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
110: TaoAddLineSearchCounts(tao);
112: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
113: /* Linesearch failed again */
114: /* switch to unscaled gradient */
115: f = f_old;
116: VecCopy(cgP->X_old, tao->solution);
117: VecCopy(cgP->G_old, tao->gradient);
118: delta = 1.0;
119: VecCopy(tao->solution, tao->stepdirection);
120: VecScale(tao->stepdirection, -1.0);
122: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
123: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
124: TaoAddLineSearchCounts(tao);
125: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
127: /* Line search failed for last time -- give up */
128: f = f_old;
129: VecCopy(cgP->X_old, tao->solution);
130: VecCopy(cgP->G_old, tao->gradient);
131: step = 0.0;
132: tao->reason = TAO_DIVERGED_LS_FAILURE;
133: }
134: }
135: }
137: /* Check for bad value */
138: VecNorm(tao->gradient,NORM_2,&gnorm);
141: /* Check for termination */
142: gnorm2 =gnorm * gnorm;
143: tao->niter++;
144: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
145: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
146: (*tao->ops->convergencetest)(tao,tao->cnvP);
147: if (tao->reason != TAO_CONTINUE_ITERATING) {
148: break;
149: }
151: /* Check for restart condition */
152: VecDot(tao->gradient, cgP->G_old, &ginner);
153: if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
154: /* Gradients far from orthogonal; use steepest descent direction */
155: beta = 0.0;
156: } else {
157: /* Gradients close to orthogonal; use conjugate gradient formula */
158: switch (cgP->cg_type) {
159: case CG_FletcherReeves:
160: beta = gnorm2 / gnorm2_old;
161: break;
163: case CG_PolakRibiere:
164: beta = (gnorm2 - ginner) / gnorm2_old;
165: break;
167: case CG_PolakRibierePlus:
168: beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
169: break;
171: case CG_HestenesStiefel:
172: VecDot(tao->gradient, tao->stepdirection, &gd);
173: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
174: beta = (gnorm2 - ginner) / (gd - gd_old);
175: break;
177: case CG_DaiYuan:
178: VecDot(tao->gradient, tao->stepdirection, &gd);
179: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
180: beta = gnorm2 / (gd - gd_old);
181: break;
183: default:
184: beta = 0.0;
185: break;
186: }
187: }
189: /* Compute the direction d=-g + beta*d */
190: VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);
192: /* update initial steplength choice */
193: delta = 1.0;
194: delta = PetscMax(delta, cgP->delta_min);
195: delta = PetscMin(delta, cgP->delta_max);
196: }
197: return 0;
198: }
200: static PetscErrorCode TaoSetUp_CG(Tao tao)
201: {
202: TAO_CG *cgP = (TAO_CG*)tao->data;
204: if (!tao->gradient) VecDuplicate(tao->solution,&tao->gradient);
205: if (!tao->stepdirection) VecDuplicate(tao->solution,&tao->stepdirection);
206: if (!cgP->X_old) VecDuplicate(tao->solution,&cgP->X_old);
207: if (!cgP->G_old) VecDuplicate(tao->gradient,&cgP->G_old);
208: return 0;
209: }
211: static PetscErrorCode TaoDestroy_CG(Tao tao)
212: {
213: TAO_CG *cgP = (TAO_CG*) tao->data;
215: if (tao->setupcalled) {
216: VecDestroy(&cgP->X_old);
217: VecDestroy(&cgP->G_old);
218: }
219: TaoLineSearchDestroy(&tao->linesearch);
220: PetscFree(tao->data);
221: return 0;
222: }
224: static PetscErrorCode TaoSetFromOptions_CG(PetscOptionItems *PetscOptionsObject,Tao tao)
225: {
226: 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;
243: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
244: if (isascii) {
245: PetscViewerASCIIPushTab(viewer);
246: PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
247: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", cgP->ngradsteps);
248: PetscViewerASCIIPrintf(viewer, "Reset steps: %D\n", cgP->nresetsteps);
249: PetscViewerASCIIPopTab(viewer);
250: }
251: return 0;
252: }
254: /*MC
255: TAOCG - Nonlinear conjugate gradient method is an extension of the
256: nonlinear conjugate gradient solver for nonlinear optimization.
258: Options Database Keys:
259: + -tao_cg_eta <r> - restart tolerance
260: . -tao_cg_type <taocg_type> - cg formula
261: . -tao_cg_delta_min <r> - minimum delta value
262: - -tao_cg_delta_max <r> - maximum delta value
264: Notes:
265: CG formulas are:
266: "fr" - Fletcher-Reeves
267: "pr" - Polak-Ribiere
268: "prp" - Polak-Ribiere-Plus
269: "hs" - Hestenes-Steifel
270: "dy" - Dai-Yuan
271: Level: beginner
272: M*/
274: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
275: {
276: TAO_CG *cgP;
277: const char *morethuente_type = TAOLINESEARCHMT;
279: tao->ops->setup = TaoSetUp_CG;
280: tao->ops->solve = TaoSolve_CG;
281: tao->ops->view = TaoView_CG;
282: tao->ops->setfromoptions = TaoSetFromOptions_CG;
283: tao->ops->destroy = TaoDestroy_CG;
285: /* Override default settings (unless already changed) */
286: if (!tao->max_it_changed) tao->max_it = 2000;
287: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
289: /* Note: nondefault values should be used for nonlinear conjugate gradient */
290: /* method. In particular, gtol should be less that 0.5; the value used in */
291: /* Nocedal and Wright is 0.10. We use the default values for the */
292: /* linesearch because it seems to work better. */
293: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
294: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
295: TaoLineSearchSetType(tao->linesearch, morethuente_type);
296: TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
297: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
299: PetscNewLog(tao,&cgP);
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: return 0;
306: }