xref: /petsc/src/tao/bound/impls/bncg/bncg.c (revision feff33ee0b5b037fa8f9f294dede656a2f85cc47) 
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 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;
64 
65   PetscFunctionBegin;
66   /*   Project the current point onto the feasible set */
67   ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr);
68   ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr);
69 
70   /* Project the initial point onto the feasible region */
71   ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr);
72 
73   /*  Compute the objective function and criteria */
74   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, cg->unprojected_gradient);CHKERRQ(ierr);
75   ierr = VecNorm(cg->unprojected_gradient,NORM_2,&gnorm);CHKERRQ(ierr);
76   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
77 
78   /* Project the gradient and calculate the norm */
79   ierr = VecBoundGradientProjection(cg->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr);
80   ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
81   gnorm2 = gnorm*gnorm;
82 
83   /* Convergence check */
84   tao->reason = TAO_CONTINUE_ITERATING;
85   ierr = TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);CHKERRQ(ierr);
86   ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);CHKERRQ(ierr);
87   ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
88   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
89 
90   /* Start optimization iterations */
91   f_old = f;
92   gnorm2_old = gnorm2;
93   ierr = VecCopy(tao->solution, cg->X_old);CHKERRQ(ierr);
94   ierr = VecCopy(tao->gradient, cg->G_old);CHKERRQ(ierr);
95   ierr = VecCopy(cg->unprojected_gradient, cg->unprojected_gradient_old);CHKERRQ(ierr);
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     ierr = VecDot(tao->gradient, cg->G_old, &ginner);CHKERRQ(ierr);
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;
120 
121       case CG_PolakRibiere:
122         beta = (gnorm2 - ginner) / gnorm2_old;
123         break;
124 
125       case CG_PolakRibierePlus:
126         beta = PetscMax((gnorm2-ginner)/gnorm2_old, 0.0);
127         break;
128 
129       case CG_HestenesStiefel:
130         ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
131         ierr = VecDot(cg->G_old, tao->stepdirection, &gd_old);CHKERRQ(ierr);
132         beta = (gnorm2 - ginner) / (gd - gd_old);
133         break;
134 
135       case CG_DaiYuan:
136         ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
137         ierr = VecDot(cg->G_old, tao->stepdirection, &gd_old);CHKERRQ(ierr);
138         beta = gnorm2 / (gd - gd_old);
139         break;
140 
141       default:
142         beta = 0.0;
143         break;
144       }
145     }
146 
147     /*  Compute the direction d=-g + beta*d */
148     ierr = VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient);CHKERRQ(ierr);
149     ierr = TaoBNCGResetStepForNewInactives(tao, tao->stepdirection);CHKERRQ(ierr);
150 
151     /* Verify that this is a descent direction */
152     ierr = VecDot(tao->gradient, tao->stepdirection, &gd);CHKERRQ(ierr);
153     ierr = 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       ierr = VecAXPBY(tao->stepdirection, -1.0, 0.0, tao->gradient);CHKERRQ(ierr);
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     ierr = VecCopy(tao->solution, cg->X_old);CHKERRQ(ierr);
168     ierr = VecCopy(tao->gradient, cg->G_old);CHKERRQ(ierr);
169     ierr = VecCopy(cg->unprojected_gradient, cg->unprojected_gradient_old);CHKERRQ(ierr);
170     gnorm2_old = gnorm2;
171     f_old = f;
172 
173     /* Perform bounded line search */
174     ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);CHKERRQ(ierr);
175     ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
176     ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
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       ierr = VecCopy(cg->X_old, tao->solution);CHKERRQ(ierr);
185       ierr = VecCopy(cg->G_old, tao->gradient);CHKERRQ(ierr);
186       ierr = VecCopy(cg->unprojected_gradient_old, cg->unprojected_gradient);CHKERRQ(ierr);
187 
188       /* Fall back on the unscaled gradient step */
189       delta = 1.0;
190       ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr);
191       ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
192 
193       ierr = TaoLineSearchSetInitialStepLength(tao->linesearch,delta);CHKERRQ(ierr);
194       ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status);CHKERRQ(ierr);
195       ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
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         ierr = VecCopy(cg->X_old, tao->solution);CHKERRQ(ierr);
203         ierr = VecCopy(cg->G_old, tao->gradient);CHKERRQ(ierr);
204         ierr = VecCopy(cg->unprojected_gradient_old, cg->unprojected_gradient);CHKERRQ(ierr);
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     }
211 
212     /* Compute the projected gradient and its norm */
213     ierr = VecBoundGradientProjection(cg->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr);
214     ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
215     gnorm2 = gnorm*gnorm;
216 
217     /* Convergence test */
218     tao->niter++;
219     ierr = TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);CHKERRQ(ierr);
220     ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);CHKERRQ(ierr);
221     ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
222   }
223   PetscFunctionReturn(0);
224 }
225 
226 static PetscErrorCode TaoSetUp_BNCG(Tao tao)
227 {
228   TAO_BNCG         *cg = (TAO_BNCG*)tao->data;
229   PetscErrorCode ierr;
230 
231   PetscFunctionBegin;
232   if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);}
233   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr); }
234   if (!cg->X_old) {ierr = VecDuplicate(tao->solution,&cg->X_old);CHKERRQ(ierr);}
235   if (!cg->G_old) {ierr = VecDuplicate(tao->gradient,&cg->G_old);CHKERRQ(ierr); }
236   if (!cg->unprojected_gradient) {ierr = VecDuplicate(tao->gradient,&cg->unprojected_gradient);CHKERRQ(ierr);}
237   if (!cg->unprojected_gradient_old) {ierr = VecDuplicate(tao->gradient,&cg->unprojected_gradient_old);CHKERRQ(ierr);}
238   PetscFunctionReturn(0);
239 }
240 
241 static PetscErrorCode TaoDestroy_BNCG(Tao tao)
242 {
243   TAO_BNCG       *cg = (TAO_BNCG*) tao->data;
244   PetscErrorCode ierr;
245 
246   PetscFunctionBegin;
247   if (tao->setupcalled) {
248     ierr = VecDestroy(&cg->X_old);CHKERRQ(ierr);
249     ierr = VecDestroy(&cg->G_old);CHKERRQ(ierr);
250     ierr = VecDestroy(&cg->unprojected_gradient);CHKERRQ(ierr);
251     ierr = VecDestroy(&cg->unprojected_gradient_old);CHKERRQ(ierr);
252   }
253   ierr = TaoLineSearchDestroy(&tao->linesearch);CHKERRQ(ierr);
254   ierr = PetscFree(tao->data);CHKERRQ(ierr);
255   PetscFunctionReturn(0);
256 }
257 
258 static PetscErrorCode TaoSetFromOptions_BNCG(PetscOptionItems *PetscOptionsObject,Tao tao)
259  {
260     TAO_BNCG       *cg = (TAO_BNCG*)tao->data;
261     PetscErrorCode ierr;
262 
263     PetscFunctionBegin;
264     ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
265     ierr = PetscOptionsHead(PetscOptionsObject,"Nonlinear Conjugate Gradient method for unconstrained optimization");CHKERRQ(ierr);
266     ierr = PetscOptionsReal("-tao_BNCG_eta","restart tolerance", "", cg->eta,&cg->eta,NULL);CHKERRQ(ierr);
267     ierr = PetscOptionsReal("-tao_BNCG_rho","descent direction tolerance", "", cg->rho,&cg->rho,NULL);CHKERRQ(ierr);
268     ierr = PetscOptionsReal("-tao_BNCG_pow","descent direction exponent", "", cg->pow,&cg->pow,NULL);CHKERRQ(ierr);
269     ierr = PetscOptionsEList("-tao_BNCG_type","cg formula", "", CG_Table, CG_Types, CG_Table[cg->cg_type], &cg->cg_type,NULL);CHKERRQ(ierr);
270     ierr = PetscOptionsReal("-tao_BNCG_delta_min","minimum delta value", "", cg->delta_min,&cg->delta_min,NULL);CHKERRQ(ierr);
271     ierr = PetscOptionsReal("-tao_BNCG_delta_max","maximum delta value", "", cg->delta_max,&cg->delta_max,NULL);CHKERRQ(ierr);
272    ierr = PetscOptionsTail();CHKERRQ(ierr);
273    PetscFunctionReturn(0);
274 }
275 
276 static PetscErrorCode TaoView_BNCG(Tao tao, PetscViewer viewer)
277 {
278   PetscBool      isascii;
279   TAO_BNCG       *cg = (TAO_BNCG*)tao->data;
280   PetscErrorCode ierr;
281 
282   PetscFunctionBegin;
283   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
284   if (isascii) {
285     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
286     ierr = PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cg->cg_type]);CHKERRQ(ierr);
287     ierr = PetscViewerASCIIPrintf(viewer, "Resets: %i\n", cg->resets);CHKERRQ(ierr);
288     ierr = PetscViewerASCIIPrintf(viewer, "  Broken ortho: %i\n", cg->broken_ortho);CHKERRQ(ierr);
289     ierr = PetscViewerASCIIPrintf(viewer, "  Not a descent dir.: %i\n", cg->descent_error);CHKERRQ(ierr);
290     ierr = PetscViewerASCIIPrintf(viewer, "Line search fails: %i\n", cg->ls_fails);CHKERRQ(ierr);
291     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
292   }
293   PetscFunctionReturn(0);
294 }
295 
296 /*MC
297      TAOBNCG -   Bound-constrained Nonlinear Conjugate Gradient method.
298 
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
304 
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*/
314 
315 
316 PETSC_EXTERN PetscErrorCode TaoCreate_BNCG(Tao tao)
317 {
318   TAO_BNCG       *cg;
319   const char     *morethuente_type = TAOLINESEARCHMT;
320   PetscErrorCode ierr;
321 
322   PetscFunctionBegin;
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;
328 
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;
332 
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   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);CHKERRQ(ierr);
338   ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr);
339   ierr = TaoLineSearchSetType(tao->linesearch, morethuente_type);CHKERRQ(ierr);
340   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch, tao);CHKERRQ(ierr);
341   ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
342 
343   ierr = PetscNewLog(tao,&cg);CHKERRQ(ierr);
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   PetscFunctionReturn(0);
352 }
353