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