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) PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by cg algorithm\n");
22: /* Check convergence criteria */
23: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
24: VecNorm(tao->gradient, NORM_2, &gnorm);
27: tao->reason = TAO_CONTINUE_ITERATING;
28: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
29: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
30: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
31: if (tao->reason != TAO_CONTINUE_ITERATING) return 0;
33: /* Set initial direction to -gradient */
34: VecCopy(tao->gradient, tao->stepdirection);
35: VecScale(tao->stepdirection, -1.0);
36: gnorm2 = gnorm * gnorm;
38: /* Set initial scaling for the function */
39: if (f != 0.0) {
40: delta = 2.0 * PetscAbsScalar(f) / gnorm2;
41: delta = PetscMax(delta, cgP->delta_min);
42: delta = PetscMin(delta, cgP->delta_max);
43: } else {
44: delta = 2.0 / gnorm2;
45: delta = PetscMax(delta, cgP->delta_min);
46: delta = PetscMin(delta, cgP->delta_max);
47: }
48: /* Set counter for gradient and reset steps */
49: cgP->ngradsteps = 0;
50: cgP->nresetsteps = 0;
52: while (1) {
53: /* Call general purpose update function */
54: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
56: /* Save the current gradient information */
57: f_old = f;
58: gnorm2_old = gnorm2;
59: VecCopy(tao->solution, cgP->X_old);
60: VecCopy(tao->gradient, cgP->G_old);
61: VecDot(tao->gradient, tao->stepdirection, &gd);
62: if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
63: ++cgP->ngradsteps;
64: if (f != 0.0) {
65: delta = 2.0 * PetscAbsScalar(f) / gnorm2;
66: delta = PetscMax(delta, cgP->delta_min);
67: delta = PetscMin(delta, cgP->delta_max);
68: } else {
69: delta = 2.0 / gnorm2;
70: delta = PetscMax(delta, cgP->delta_min);
71: delta = PetscMin(delta, cgP->delta_max);
72: }
74: VecCopy(tao->gradient, tao->stepdirection);
75: VecScale(tao->stepdirection, -1.0);
76: }
78: /* Search direction for improving point */
79: TaoLineSearchSetInitialStepLength(tao->linesearch, delta);
80: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
81: TaoAddLineSearchCounts(tao);
82: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
83: /* Linesearch failed */
84: /* Reset factors and use scaled gradient step */
85: ++cgP->nresetsteps;
86: f = f_old;
87: gnorm2 = gnorm2_old;
88: VecCopy(cgP->X_old, tao->solution);
89: VecCopy(cgP->G_old, tao->gradient);
91: if (f != 0.0) {
92: delta = 2.0 * PetscAbsScalar(f) / gnorm2;
93: delta = PetscMax(delta, cgP->delta_min);
94: delta = PetscMin(delta, cgP->delta_max);
95: } else {
96: delta = 2.0 / gnorm2;
97: delta = PetscMax(delta, cgP->delta_min);
98: delta = PetscMin(delta, cgP->delta_max);
99: }
101: VecCopy(tao->gradient, tao->stepdirection);
102: VecScale(tao->stepdirection, -1.0);
104: TaoLineSearchSetInitialStepLength(tao->linesearch, delta);
105: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
106: TaoAddLineSearchCounts(tao);
108: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
109: /* Linesearch failed again */
110: /* switch to unscaled gradient */
111: f = f_old;
112: VecCopy(cgP->X_old, tao->solution);
113: VecCopy(cgP->G_old, tao->gradient);
114: delta = 1.0;
115: VecCopy(tao->solution, tao->stepdirection);
116: VecScale(tao->stepdirection, -1.0);
118: TaoLineSearchSetInitialStepLength(tao->linesearch, delta);
119: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status);
120: TaoAddLineSearchCounts(tao);
121: 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: tao->reason = TAO_DIVERGED_LS_FAILURE;
128: }
129: }
130: }
132: /* Check for bad value */
133: VecNorm(tao->gradient, NORM_2, &gnorm);
136: /* Check for termination */
137: gnorm2 = gnorm * gnorm;
138: tao->niter++;
139: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
140: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
141: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
142: if (tao->reason != TAO_CONTINUE_ITERATING) break;
144: /* Check for restart condition */
145: VecDot(tao->gradient, cgP->G_old, &ginner);
146: if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
147: /* Gradients far from orthogonal; use steepest descent direction */
148: beta = 0.0;
149: } else {
150: /* Gradients close to orthogonal; use conjugate gradient formula */
151: switch (cgP->cg_type) {
152: case CG_FletcherReeves:
153: beta = gnorm2 / gnorm2_old;
154: break;
156: case CG_PolakRibiere:
157: beta = (gnorm2 - ginner) / gnorm2_old;
158: break;
160: case CG_PolakRibierePlus:
161: beta = PetscMax((gnorm2 - ginner) / gnorm2_old, 0.0);
162: break;
164: case CG_HestenesStiefel:
165: VecDot(tao->gradient, tao->stepdirection, &gd);
166: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
167: beta = (gnorm2 - ginner) / (gd - gd_old);
168: break;
170: case CG_DaiYuan:
171: VecDot(tao->gradient, tao->stepdirection, &gd);
172: VecDot(cgP->G_old, tao->stepdirection, &gd_old);
173: beta = gnorm2 / (gd - gd_old);
174: break;
176: default:
177: beta = 0.0;
178: break;
179: }
180: }
182: /* Compute the direction d=-g + beta*d */
183: VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);
185: /* update initial steplength choice */
186: delta = 1.0;
187: delta = PetscMax(delta, cgP->delta_min);
188: delta = PetscMin(delta, cgP->delta_max);
189: }
190: return 0;
191: }
193: static PetscErrorCode TaoSetUp_CG(Tao tao)
194: {
195: TAO_CG *cgP = (TAO_CG *)tao->data;
197: if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
198: if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
199: if (!cgP->X_old) VecDuplicate(tao->solution, &cgP->X_old);
200: if (!cgP->G_old) VecDuplicate(tao->gradient, &cgP->G_old);
201: return 0;
202: }
204: static PetscErrorCode TaoDestroy_CG(Tao tao)
205: {
206: TAO_CG *cgP = (TAO_CG *)tao->data;
208: if (tao->setupcalled) {
209: VecDestroy(&cgP->X_old);
210: VecDestroy(&cgP->G_old);
211: }
212: TaoLineSearchDestroy(&tao->linesearch);
213: PetscFree(tao->data);
214: return 0;
215: }
217: static PetscErrorCode TaoSetFromOptions_CG(Tao tao, PetscOptionItems *PetscOptionsObject)
218: {
219: TAO_CG *cgP = (TAO_CG *)tao->data;
221: TaoLineSearchSetFromOptions(tao->linesearch);
222: PetscOptionsHeadBegin(PetscOptionsObject, "Nonlinear Conjugate Gradient method for unconstrained optimization");
223: PetscOptionsReal("-tao_cg_eta", "restart tolerance", "", cgP->eta, &cgP->eta, NULL);
224: PetscOptionsEList("-tao_cg_type", "cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, NULL);
225: PetscOptionsReal("-tao_cg_delta_min", "minimum delta value", "", cgP->delta_min, &cgP->delta_min, NULL);
226: PetscOptionsReal("-tao_cg_delta_max", "maximum delta value", "", cgP->delta_max, &cgP->delta_max, NULL);
227: PetscOptionsHeadEnd();
228: return 0;
229: }
231: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
232: {
233: PetscBool isascii;
234: TAO_CG *cgP = (TAO_CG *)tao->data;
236: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
237: if (isascii) {
238: PetscViewerASCIIPushTab(viewer);
239: PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]);
240: PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", cgP->ngradsteps);
241: PetscViewerASCIIPrintf(viewer, "Reset steps: %" PetscInt_FMT "\n", cgP->nresetsteps);
242: PetscViewerASCIIPopTab(viewer);
243: }
244: return 0;
245: }
247: /*MC
248: TAOCG - Nonlinear conjugate gradient method is an extension of the
249: nonlinear conjugate gradient solver for nonlinear optimization.
251: Options Database Keys:
252: + -tao_cg_eta <r> - restart tolerance
253: . -tao_cg_type <taocg_type> - cg formula
254: . -tao_cg_delta_min <r> - minimum delta value
255: - -tao_cg_delta_max <r> - maximum delta value
257: Notes:
258: CG formulas are:
259: "fr" - Fletcher-Reeves
260: "pr" - Polak-Ribiere
261: "prp" - Polak-Ribiere-Plus
262: "hs" - Hestenes-Steifel
263: "dy" - Dai-Yuan
264: Level: beginner
265: M*/
267: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
268: {
269: TAO_CG *cgP;
270: const char *morethuente_type = TAOLINESEARCHMT;
272: tao->ops->setup = TaoSetUp_CG;
273: tao->ops->solve = TaoSolve_CG;
274: tao->ops->view = TaoView_CG;
275: tao->ops->setfromoptions = TaoSetFromOptions_CG;
276: tao->ops->destroy = TaoDestroy_CG;
278: /* Override default settings (unless already changed) */
279: if (!tao->max_it_changed) tao->max_it = 2000;
280: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
282: /* Note: nondefault values should be used for nonlinear conjugate gradient */
283: /* method. In particular, gtol should be less that 0.5; the value used in */
284: /* Nocedal and Wright is 0.10. We use the default values for the */
285: /* linesearch because it seems to work better. */
286: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
287: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
288: TaoLineSearchSetType(tao->linesearch, morethuente_type);
289: TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
290: TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);
292: PetscNew(&cgP);
293: tao->data = (void *)cgP;
294: cgP->eta = 0.1;
295: cgP->delta_min = 1e-7;
296: cgP->delta_max = 100;
297: cgP->cg_type = CG_PolakRibierePlus;
298: return 0;
299: }