Actual source code: bncg.c
petsc-3.9.4 2018-09-11
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/bound/impls/bncg/bncg.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: PetscErrorCode TaoBNCGResetStepForNewInactives(Tao tao, Vec step)
14: {
15: TAO_BNCG *cg = (TAO_BNCG*)tao->data;
16: PetscErrorCode ierr;
17: const PetscScalar *xl, *xo, *xn, *xu, *gn, *go;
18: PetscInt size, i;
19: PetscScalar *s;
22: VecGetLocalSize(tao->solution, &size);
23: VecGetArrayRead(cg->unprojected_gradient_old, &go);
24: VecGetArrayRead(cg->unprojected_gradient, &gn);
25: VecGetArrayRead(cg->X_old, &xo);
26: VecGetArrayRead(tao->solution, &xn);
27: VecGetArrayRead(tao->XL, &xl);
28: VecGetArrayRead(tao->XU, &xu);
29: VecGetArray(step, &s);
30: for (i=0; i<size; i++) {
31: if (xl[i] == xu[i]) {
32: s[i] = 0.0;
33: } else {
34: if (xl[i] > PETSC_NINFINITY) {
35: if ((xn[i] == xl[i] && gn[i] < 0.0) && (xo[i] == xl[i] && go[i] >= 0.0)) {
36: s[i] = -gn[i];
37: }
38: }
39: if (xu[i] < PETSC_NINFINITY) {
40: if ((xn[i] == xu[i] && gn[i] > 0.0) && (xo[i] == xu[i] && go[i] <= 0.0)) {
41: s[i] = -gn[i];
42: }
43: }
44: }
45: }
46: VecRestoreArrayRead(cg->unprojected_gradient_old, &go);
47: VecRestoreArrayRead(cg->unprojected_gradient, &gn);
48: VecRestoreArrayRead(cg->X_old, &xo);
49: VecRestoreArrayRead(tao->solution, &xn);
50: VecRestoreArrayRead(tao->XL, &xl);
51: VecRestoreArrayRead(tao->XU, &xu);
52: VecRestoreArray(step, &s);
53: return(0);
54: }
56: static PetscErrorCode TaoSolve_BNCG(Tao tao)
57: {
58: TAO_BNCG *cg = (TAO_BNCG*)tao->data;
59: PetscErrorCode ierr;
60: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
61: PetscReal step=1.0,f,gnorm,gnorm2,delta,gd,ginner,beta,dnorm;
62: PetscReal gd_old,gnorm2_old,f_old;
63: PetscBool cg_restart;
66: /* Project the current point onto the feasible set */
67: TaoComputeVariableBounds(tao);
68: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
69:
70: /* Project the initial point onto the feasible region */
71: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
73: /* Compute the objective function and criteria */
74: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, cg->unprojected_gradient);
75: VecNorm(cg->unprojected_gradient,NORM_2,&gnorm);
76: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
78: /* Project the gradient and calculate the norm */
79: VecBoundGradientProjection(cg->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
80: VecNorm(tao->gradient,NORM_2,&gnorm);
81: gnorm2 = gnorm*gnorm;
82:
83: /* Convergence check */
84: tao->reason = TAO_CONTINUE_ITERATING;
85: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
86: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
87: (*tao->ops->convergencetest)(tao,tao->cnvP);
88: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
89:
90: /* Start optimization iterations */
91: f_old = f;
92: gnorm2_old = gnorm2;
93: VecCopy(tao->solution, cg->X_old);
94: VecCopy(tao->gradient, cg->G_old);
95: VecCopy(cg->unprojected_gradient, cg->unprojected_gradient_old);
96: tao->niter = cg->ls_fails = cg->broken_ortho = cg->descent_error = 0;
97: cg->resets = -1;
98: while (tao->reason == TAO_CONTINUE_ITERATING) {
99: /* Check restart conditions for using steepest descent */
100: cg_restart = PETSC_FALSE;
101: VecDot(tao->gradient, cg->G_old, &ginner);
102: if (tao->niter == 0) {
103: /* 1) First iteration */
104: cg_restart = PETSC_TRUE;
105: } else if (PetscAbsScalar(ginner) >= cg->eta * gnorm2) {
106: /* 2) Gradients are far from orthogonal */
107: cg_restart = PETSC_TRUE;
108: cg->broken_ortho++;
109: }
110:
111: /* Compute CG step */
112: if (cg_restart) {
113: beta = 0.0;
114: cg->resets++;
115: } else {
116: switch (cg->cg_type) {
117: case CG_FletcherReeves:
118: beta = gnorm2 / gnorm2_old;
119: break;
121: case CG_PolakRibiere:
122: beta = (gnorm2 - ginner) / gnorm2_old;
123: break;
125: case CG_PolakRibierePlus:
126: beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
127: break;
129: case CG_HestenesStiefel:
130: VecDot(tao->gradient, tao->stepdirection, &gd);
131: VecDot(cg->G_old, tao->stepdirection, &gd_old);
132: beta = (gnorm2 - ginner) / (gd - gd_old);
133: break;
135: case CG_DaiYuan:
136: VecDot(tao->gradient, tao->stepdirection, &gd);
137: VecDot(cg->G_old, tao->stepdirection, &gd_old);
138: beta = gnorm2 / (gd - gd_old);
139: break;
141: default:
142: beta = 0.0;
143: break;
144: }
145: }
146:
147: /* Compute the direction d=-g + beta*d */
148: VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);
149: TaoBNCGResetStepForNewInactives(tao, tao->stepdirection);
150:
151: /* Verify that this is a descent direction */
152: VecDot(tao->gradient, tao->stepdirection, &gd);
153: VecNorm(tao->stepdirection, NORM_2, &dnorm);
154: if (gd > -cg->rho*PetscPowReal(dnorm, cg->pow)) {
155: /* Not a descent direction, so we reset back to projected gradient descent */
156: VecAXPBY(tao->stepdirection, -1.0, 0.0, tao->gradient);
157: cg->resets++;
158: cg->descent_error++;
159: }
160:
161: /* update initial steplength choice */
162: delta = 1.0;
163: delta = PetscMax(delta, cg->delta_min);
164: delta = PetscMin(delta, cg->delta_max);
165:
166: /* Store solution and gradient info before it changes */
167: VecCopy(tao->solution, cg->X_old);
168: VecCopy(tao->gradient, cg->G_old);
169: VecCopy(cg->unprojected_gradient, cg->unprojected_gradient_old);
170: gnorm2_old = gnorm2;
171: f_old = f;
172:
173: /* Perform bounded line search */
174: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
175: TaoLineSearchApply(tao->linesearch, tao->solution, &f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status);
176: TaoAddLineSearchCounts(tao);
177:
178: /* Check linesearch failure */
179: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
180: cg->ls_fails++;
181: /* Restore previous point */
182: gnorm2 = gnorm2_old;
183: f = f_old;
184: VecCopy(cg->X_old, tao->solution);
185: VecCopy(cg->G_old, tao->gradient);
186: VecCopy(cg->unprojected_gradient_old, cg->unprojected_gradient);
187:
188: /* Fall back on the unscaled gradient step */
189: delta = 1.0;
190: VecCopy(tao->solution, tao->stepdirection);
191: VecScale(tao->stepdirection, -1.0);
192:
193: TaoLineSearchSetInitialStepLength(tao->linesearch,delta);
194: TaoLineSearchApply(tao->linesearch, tao->solution, &f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status);
195: TaoAddLineSearchCounts(tao);
196:
197: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER){
198: cg->ls_fails++;
199: /* Restore previous point */
200: gnorm2 = gnorm2_old;
201: f = f_old;
202: VecCopy(cg->X_old, tao->solution);
203: VecCopy(cg->G_old, tao->gradient);
204: VecCopy(cg->unprojected_gradient_old, cg->unprojected_gradient);
205:
206: /* Nothing left to do but fail out of the optimization */
207: step = 0.0;
208: tao->reason = TAO_DIVERGED_LS_FAILURE;
209: }
210: }
212: /* Compute the projected gradient and its norm */
213: VecBoundGradientProjection(cg->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
214: VecNorm(tao->gradient,NORM_2,&gnorm);
215: gnorm2 = gnorm*gnorm;
216:
217: /* Convergence test */
218: tao->niter++;
219: TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
220: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
221: (*tao->ops->convergencetest)(tao,tao->cnvP);
222: }
223: return(0);
224: }
226: static PetscErrorCode TaoSetUp_BNCG(Tao tao)
227: {
228: TAO_BNCG *cg = (TAO_BNCG*)tao->data;
232: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);}
233: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
234: if (!cg->X_old) {VecDuplicate(tao->solution,&cg->X_old);}
235: if (!cg->G_old) {VecDuplicate(tao->gradient,&cg->G_old); }
236: if (!cg->unprojected_gradient) {VecDuplicate(tao->gradient,&cg->unprojected_gradient);}
237: if (!cg->unprojected_gradient_old) {VecDuplicate(tao->gradient,&cg->unprojected_gradient_old);}
238: return(0);
239: }
241: static PetscErrorCode TaoDestroy_BNCG(Tao tao)
242: {
243: TAO_BNCG *cg = (TAO_BNCG*) tao->data;
247: if (tao->setupcalled) {
248: VecDestroy(&cg->X_old);
249: VecDestroy(&cg->G_old);
250: VecDestroy(&cg->unprojected_gradient);
251: VecDestroy(&cg->unprojected_gradient_old);
252: }
253: TaoLineSearchDestroy(&tao->linesearch);
254: PetscFree(tao->data);
255: return(0);
256: }
258: static PetscErrorCode TaoSetFromOptions_BNCG(PetscOptionItems *PetscOptionsObject,Tao tao)
259: {
260: TAO_BNCG *cg = (TAO_BNCG*)tao->data;
264: TaoLineSearchSetFromOptions(tao->linesearch);
265: PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");
266: PetscOptionsReal("-tao_BNCG_eta","restart tolerance", "", cg->eta,&cg->eta,NULL);
267: PetscOptionsReal("-tao_BNCG_rho","descent direction tolerance", "", cg->rho,&cg->rho,NULL);
268: PetscOptionsReal("-tao_BNCG_pow","descent direction exponent", "", cg->pow,&cg->pow,NULL);
269: PetscOptionsEList("-tao_BNCG_type","cg formula", "", CG_Table, CG_Types, CG_Table[cg->cg_type], &cg->cg_type,NULL);
270: PetscOptionsReal("-tao_BNCG_delta_min","minimum delta value", "", cg->delta_min,&cg->delta_min,NULL);
271: PetscOptionsReal("-tao_BNCG_delta_max","maximum delta value", "", cg->delta_max,&cg->delta_max,NULL);
272: PetscOptionsTail();
273: return(0);
274: }
276: static PetscErrorCode TaoView_BNCG(Tao tao, PetscViewer viewer)
277: {
278: PetscBool isascii;
279: TAO_BNCG *cg = (TAO_BNCG*)tao->data;
283: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
284: if (isascii) {
285: PetscViewerASCIIPushTab(viewer);
286: PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cg->cg_type]);
287: PetscViewerASCIIPrintf(viewer, "Resets: %i\n", cg->resets);
288: PetscViewerASCIIPrintf(viewer, " Broken ortho: %i\n", cg->broken_ortho);
289: PetscViewerASCIIPrintf(viewer, " Not a descent dir.: %i\n", cg->descent_error);
290: PetscViewerASCIIPrintf(viewer, "Line search fails: %i\n", cg->ls_fails);
291: PetscViewerASCIIPopTab(viewer);
292: }
293: return(0);
294: }
296: /*MC
297: TAOBNCG - Bound-constrained Nonlinear Conjugate Gradient method.
299: Options Database Keys:
300: + -tao_BNCG_eta <r> - restart tolerance
301: . -tao_BNCG_type <taocg_type> - cg formula
302: . -tao_BNCG_delta_min <r> - minimum delta value
303: - -tao_BNCG_delta_max <r> - maximum delta value
305: Notes:
306: CG formulas are:
307: "fr" - Fletcher-Reeves
308: "pr" - Polak-Ribiere
309: "prp" - Polak-Ribiere-Plus
310: "hs" - Hestenes-Steifel
311: "dy" - Dai-Yuan
312: Level: beginner
313: M*/
316: PETSC_EXTERN PetscErrorCode TaoCreate_BNCG(Tao tao)
317: {
318: TAO_BNCG *cg;
319: const char *morethuente_type = TAOLINESEARCHMT;
323: tao->ops->setup = TaoSetUp_BNCG;
324: tao->ops->solve = TaoSolve_BNCG;
325: tao->ops->view = TaoView_BNCG;
326: tao->ops->setfromoptions = TaoSetFromOptions_BNCG;
327: tao->ops->destroy = TaoDestroy_BNCG;
329: /* Override default settings (unless already changed) */
330: if (!tao->max_it_changed) tao->max_it = 2000;
331: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
333: /* Note: nondefault values should be used for nonlinear conjugate gradient */
334: /* method. In particular, gtol should be less that 0.5; the value used in */
335: /* Nocedal and Wright is 0.10. We use the default values for the */
336: /* linesearch because it seems to work better. */
337: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
338: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
339: TaoLineSearchSetType(tao->linesearch, morethuente_type);
340: TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
341: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
343: PetscNewLog(tao,&cg);
344: tao->data = (void*)cg;
345: cg->rho = 1e-4;
346: cg->pow = 2.1;
347: cg->eta = 0.5;
348: cg->delta_min = 1e-7;
349: cg->delta_max = 100;
350: cg->cg_type = CG_DaiYuan;
351: return(0);
352: }