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: PetscFunctionBegin;
21: if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by cg algorithm\n"));
23: /* Check convergence criteria */
24: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
25: PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
26: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
28: tao->reason = TAO_CONTINUE_ITERATING;
29: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
30: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
31: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
32: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
34: /* Set initial direction to -gradient */
35: PetscCall(VecCopy(tao->gradient, tao->stepdirection));
36: PetscCall(VecScale(tao->stepdirection, -1.0));
37: gnorm2 = gnorm * gnorm;
39: /* Set initial scaling for the function */
40: if (f != 0.0) {
41: delta = 2.0 * PetscAbsScalar(f) / gnorm2;
42: delta = PetscMax(delta, cgP->delta_min);
43: delta = PetscMin(delta, cgP->delta_max);
44: } else {
45: delta = 2.0 / gnorm2;
46: delta = PetscMax(delta, cgP->delta_min);
47: delta = PetscMin(delta, cgP->delta_max);
48: }
49: /* Set counter for gradient and reset steps */
50: cgP->ngradsteps = 0;
51: cgP->nresetsteps = 0;
53: while (1) {
54: /* Call general purpose update function */
55: if (tao->ops->update) {
56: PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
57: PetscCall(TaoComputeObjective(tao, tao->solution, &f));
58: }
59: /* Save the current gradient information */
60: f_old = f;
61: gnorm2_old = gnorm2;
62: PetscCall(VecCopy(tao->solution, cgP->X_old));
63: PetscCall(VecCopy(tao->gradient, cgP->G_old));
64: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
65: if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
66: ++cgP->ngradsteps;
67: if (f != 0.0) {
68: delta = 2.0 * PetscAbsScalar(f) / gnorm2;
69: delta = PetscMax(delta, cgP->delta_min);
70: delta = PetscMin(delta, cgP->delta_max);
71: } else {
72: delta = 2.0 / gnorm2;
73: delta = PetscMax(delta, cgP->delta_min);
74: delta = PetscMin(delta, cgP->delta_max);
75: }
77: PetscCall(VecCopy(tao->gradient, tao->stepdirection));
78: PetscCall(VecScale(tao->stepdirection, -1.0));
79: }
81: /* Search direction for improving point */
82: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta));
83: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status));
84: PetscCall(TaoAddLineSearchCounts(tao));
85: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
86: /* Linesearch failed */
87: /* Reset factors and use scaled gradient step */
88: ++cgP->nresetsteps;
89: f = f_old;
90: gnorm2 = gnorm2_old;
91: PetscCall(VecCopy(cgP->X_old, tao->solution));
92: PetscCall(VecCopy(cgP->G_old, tao->gradient));
94: if (f != 0.0) {
95: delta = 2.0 * PetscAbsScalar(f) / gnorm2;
96: delta = PetscMax(delta, cgP->delta_min);
97: delta = PetscMin(delta, cgP->delta_max);
98: } else {
99: delta = 2.0 / gnorm2;
100: delta = PetscMax(delta, cgP->delta_min);
101: delta = PetscMin(delta, cgP->delta_max);
102: }
104: PetscCall(VecCopy(tao->gradient, tao->stepdirection));
105: PetscCall(VecScale(tao->stepdirection, -1.0));
107: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta));
108: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status));
109: PetscCall(TaoAddLineSearchCounts(tao));
111: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
112: /* Linesearch failed again */
113: /* switch to unscaled gradient */
114: f = f_old;
115: PetscCall(VecCopy(cgP->X_old, tao->solution));
116: PetscCall(VecCopy(cgP->G_old, tao->gradient));
117: delta = 1.0;
118: PetscCall(VecCopy(tao->solution, tao->stepdirection));
119: PetscCall(VecScale(tao->stepdirection, -1.0));
121: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta));
122: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status));
123: PetscCall(TaoAddLineSearchCounts(tao));
124: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
125: /* Line search failed for last time -- give up */
126: f = f_old;
127: PetscCall(VecCopy(cgP->X_old, tao->solution));
128: PetscCall(VecCopy(cgP->G_old, tao->gradient));
129: step = 0.0;
130: tao->reason = TAO_DIVERGED_LS_FAILURE;
131: }
132: }
133: }
135: /* Check for bad value */
136: PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
137: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User-provided compute function generated Inf or NaN");
139: /* Check for termination */
140: gnorm2 = gnorm * gnorm;
141: tao->niter++;
142: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
143: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
144: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
145: if (tao->reason != TAO_CONTINUE_ITERATING) break;
147: /* Check for restart condition */
148: PetscCall(VecDot(tao->gradient, cgP->G_old, &ginner));
149: if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
150: /* Gradients far from orthogonal; use steepest descent direction */
151: beta = 0.0;
152: } else {
153: /* Gradients close to orthogonal; use conjugate gradient formula */
154: switch (cgP->cg_type) {
155: case CG_FletcherReeves:
156: beta = gnorm2 / gnorm2_old;
157: break;
159: case CG_PolakRibiere:
160: beta = (gnorm2 - ginner) / gnorm2_old;
161: break;
163: case CG_PolakRibierePlus:
164: beta = PetscMax((gnorm2 - ginner) / gnorm2_old, 0.0);
165: break;
167: case CG_HestenesStiefel:
168: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
169: PetscCall(VecDot(cgP->G_old, tao->stepdirection, &gd_old));
170: beta = (gnorm2 - ginner) / (gd - gd_old);
171: break;
173: case CG_DaiYuan:
174: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
175: PetscCall(VecDot(cgP->G_old, tao->stepdirection, &gd_old));
176: beta = gnorm2 / (gd - gd_old);
177: break;
179: default:
180: beta = 0.0;
181: break;
182: }
183: }
185: /* Compute the direction d=-g + beta*d */
186: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient));
188: /* update initial steplength choice */
189: delta = 1.0;
190: delta = PetscMax(delta, cgP->delta_min);
191: delta = PetscMin(delta, cgP->delta_max);
192: }
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: static PetscErrorCode TaoSetUp_CG(Tao tao)
197: {
198: TAO_CG *cgP = (TAO_CG *)tao->data;
200: PetscFunctionBegin;
201: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
202: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
203: if (!cgP->X_old) PetscCall(VecDuplicate(tao->solution, &cgP->X_old));
204: if (!cgP->G_old) PetscCall(VecDuplicate(tao->gradient, &cgP->G_old));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: static PetscErrorCode TaoDestroy_CG(Tao tao)
209: {
210: TAO_CG *cgP = (TAO_CG *)tao->data;
212: PetscFunctionBegin;
213: if (tao->setupcalled) {
214: PetscCall(VecDestroy(&cgP->X_old));
215: PetscCall(VecDestroy(&cgP->G_old));
216: }
217: PetscCall(TaoLineSearchDestroy(&tao->linesearch));
218: PetscCall(PetscFree(tao->data));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode TaoSetFromOptions_CG(Tao tao, PetscOptionItems *PetscOptionsObject)
223: {
224: TAO_CG *cgP = (TAO_CG *)tao->data;
226: PetscFunctionBegin;
227: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
228: PetscOptionsHeadBegin(PetscOptionsObject, "Nonlinear Conjugate Gradient method for unconstrained optimization");
229: PetscCall(PetscOptionsReal("-tao_cg_eta", "restart tolerance", "", cgP->eta, &cgP->eta, NULL));
230: PetscCall(PetscOptionsEList("-tao_cg_type", "cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, NULL));
231: PetscCall(PetscOptionsReal("-tao_cg_delta_min", "minimum delta value", "", cgP->delta_min, &cgP->delta_min, NULL));
232: PetscCall(PetscOptionsReal("-tao_cg_delta_max", "maximum delta value", "", cgP->delta_max, &cgP->delta_max, NULL));
233: PetscOptionsHeadEnd();
234: PetscFunctionReturn(PETSC_SUCCESS);
235: }
237: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
238: {
239: PetscBool isascii;
240: TAO_CG *cgP = (TAO_CG *)tao->data;
242: PetscFunctionBegin;
243: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
244: if (isascii) {
245: PetscCall(PetscViewerASCIIPushTab(viewer));
246: PetscCall(PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]));
247: PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", cgP->ngradsteps));
248: PetscCall(PetscViewerASCIIPrintf(viewer, "Reset steps: %" PetscInt_FMT "\n", cgP->nresetsteps));
249: PetscCall(PetscViewerASCIIPopTab(viewer));
250: }
251: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
280: tao->ops->setup = TaoSetUp_CG;
281: tao->ops->solve = TaoSolve_CG;
282: tao->ops->view = TaoView_CG;
283: tao->ops->setfromoptions = TaoSetFromOptions_CG;
284: tao->ops->destroy = TaoDestroy_CG;
286: /* Override default settings (unless already changed) */
287: PetscCall(TaoParametersInitialize(tao));
288: PetscObjectParameterSetDefault(tao, max_it, 2000);
289: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
291: /* Note: nondefault values should be used for nonlinear conjugate gradient */
292: /* method. In particular, gtol should be less that 0.5; the value used in */
293: /* Nocedal and Wright is 0.10. We use the default values for the */
294: /* linesearch because it seems to work better. */
295: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
296: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
297: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
298: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
299: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
301: PetscCall(PetscNew(&cgP));
302: tao->data = (void *)cgP;
303: cgP->eta = 0.1;
304: cgP->delta_min = 1e-7;
305: cgP->delta_max = 100;
306: cgP->cg_type = CG_PolakRibierePlus;
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }