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