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