xref: /petsc/src/tao/bound/impls/bnk/bnls.c (revision eb9107154965f6d42900b472f8cc894abc73f56d) 
1 #include <../src/tao/bound/impls/bnk/bnk.h>
2 #include <petscksp.h>
3 
4 /*
5  Implements Newton's Method with a line search approach for solving
6  unconstrained minimization problems.  A More'-Thuente line search
7  is used to guarantee that the bfgs preconditioner remains positive
8  definite.
9 
10  The method can shift the Hessian matrix.  The shifting procedure is
11  adapted from the PATH algorithm for solving complementarity
12  problems.
13 
14  The linear system solve should be done with a conjugate gradient
15  method, although any method can be used.
16 */
17 
18 static PetscErrorCode TaoSolve_BNLS(Tao tao)
19 {
20   PetscErrorCode               ierr;
21   TAO_BNK                      *bnk = (TAO_BNK *)tao->data;
22   TaoLineSearchConvergedReason ls_reason;
23 
24   PetscReal                    prered, actred, kappa;
25   PetscReal                    tau_1, tau_2, tau_max, tau_min;
26   PetscReal                    f_full, fold, gdx;
27   PetscReal                    step = 1.0;
28   PetscReal                    delta;
29   PetscReal                    norm_d = 0.0, e_min;
30 
31   PetscInt                     stepType;
32   PetscInt                     bfgsUpdates = 0;
33   PetscInt                     needH = 1;
34 
35   PetscFunctionBegin;
36   if (tao->XL || tao->XU || tao->ops->computebounds) {
37     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by nls algorithm\n");CHKERRQ(ierr);
38   }
39 
40   /* Check convergence criteria */
41   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, tao->gradient);CHKERRQ(ierr);
42   ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr);
43   if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
44 
45   tao->reason = TAO_CONTINUE_ITERATING;
46   ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
47   ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,step);CHKERRQ(ierr);
48   ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
49   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
50 
51   /* Initialize the preconditioner and trust radius */
52   ierr = TaoBNKInitialize(tao);CHKERRQ(ierr);
53 
54   /* Have not converged; continue with Newton method */
55   while (tao->reason == TAO_CONTINUE_ITERATING) {
56     ++tao->niter;
57     tao->ksp_its=0;
58 
59     /* Use the common BNK kernel to compute the step */
60     ierr = TaoBNKComputeStep(tao, &stepType);CHKERRQ(ierr);
61 
62     /* Perform the linesearch */
63     fold = bnk->f;
64     ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr);
65     ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr);
66 
67     ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, tao->gradient, bnk->D, &step, &ls_reason);CHKERRQ(ierr);
68     ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
69 
70     while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != BNK_GRADIENT) {
71       /* Linesearch failed, revert solution */
72       bnk->f = fold;
73       ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr);
74       ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr);
75 
76       switch(stepType) {
77       case BNK_NEWTON:
78         /* Failed to obtain acceptable iterate with Newton 1step
79            Update the perturbation for next time */
80         if (bnk->pert <= 0.0) {
81           /* Initialize the perturbation */
82           bnk->pert = PetscMin(bnk->imax, PetscMax(bnk->imin, bnk->imfac * bnk->gnorm));
83           if (bnk->is_gltr) {
84             ierr = KSPCGGLTRGetMinEig(tao->ksp,&e_min);CHKERRQ(ierr);
85             bnk->pert = PetscMax(bnk->pert, -e_min);
86           }
87         } else {
88           /* Increase the perturbation */
89           bnk->pert = PetscMin(bnk->pmax, PetscMax(bnk->pgfac * bnk->pert, bnk->pmgfac * bnk->gnorm));
90         }
91 
92         if (BNK_PC_BFGS != bnk->pc_type) {
93           /* We don't have the bfgs matrix around and being updated
94              Must use gradient direction in this case */
95           ierr = VecCopy(tao->gradient, bnk->D);CHKERRQ(ierr);
96           ++bnk->grad;
97           stepType = BNK_GRADIENT;
98         } else {
99           /* Attempt to use the BFGS direction */
100           ierr = MatLMVMSolve(bnk->M, tao->gradient, bnk->D);CHKERRQ(ierr);
101           /* Check for success (descent direction) */
102           ierr = VecDot(tao->solution, bnk->D, &gdx);CHKERRQ(ierr);
103           if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
104             /* BFGS direction is not descent or direction produced not a number
105                We can assert bfgsUpdates > 1 in this case
106                Use steepest descent direction (scaled) */
107 
108             if (bnk->f != 0.0) {
109               delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm);
110             } else {
111               delta = 2.0 / (bnk->gnorm*bnk->gnorm);
112             }
113             ierr = MatLMVMSetDelta(bnk->M, delta);CHKERRQ(ierr);
114             ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr);
115             ierr = MatLMVMUpdate(bnk->M, tao->solution, tao->gradient);CHKERRQ(ierr);
116             ierr = MatLMVMSolve(bnk->M, tao->gradient, bnk->D);CHKERRQ(ierr);
117 
118             bfgsUpdates = 1;
119             ++bnk->sgrad;
120             stepType = BNK_SCALED_GRADIENT;
121           } else {
122             ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr);
123             if (1 == bfgsUpdates) {
124               /* The first BFGS direction is always the scaled gradient */
125               ++bnk->sgrad;
126               stepType = BNK_SCALED_GRADIENT;
127             } else {
128               ++bnk->bfgs;
129               stepType = BNK_BFGS;
130             }
131           }
132         }
133         break;
134 
135       case BNK_BFGS:
136         /* Can only enter if pc_type == BNK_PC_BFGS
137            Failed to obtain acceptable iterate with BFGS step
138            Attempt to use the scaled gradient direction */
139 
140         if (bnk->f != 0.0) {
141           delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm);
142         } else {
143           delta = 2.0 / (bnk->gnorm*bnk->gnorm);
144         }
145         ierr = MatLMVMSetDelta(bnk->M, delta);CHKERRQ(ierr);
146         ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr);
147         ierr = MatLMVMUpdate(bnk->M, tao->solution, tao->gradient);CHKERRQ(ierr);
148         ierr = MatLMVMSolve(bnk->M, tao->gradient, bnk->D);CHKERRQ(ierr);
149 
150         bfgsUpdates = 1;
151         ++bnk->sgrad;
152         stepType = BNK_SCALED_GRADIENT;
153         break;
154 
155       case BNK_SCALED_GRADIENT:
156         /* Can only enter if pc_type == BNK_PC_BFGS
157            The scaled gradient step did not produce a new iterate;
158            attemp to use the gradient direction.
159            Need to make sure we are not using a different diagonal scaling */
160 
161         ierr = MatLMVMSetScale(bnk->M,0);CHKERRQ(ierr);
162         ierr = MatLMVMSetDelta(bnk->M,1.0);CHKERRQ(ierr);
163         ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr);
164         ierr = MatLMVMUpdate(bnk->M, tao->solution, tao->gradient);CHKERRQ(ierr);
165         ierr = MatLMVMSolve(bnk->M, tao->gradient, bnk->D);CHKERRQ(ierr);
166 
167         bfgsUpdates = 1;
168         ++bnk->grad;
169         stepType = BNK_GRADIENT;
170         break;
171       }
172       ierr = VecScale(bnk->D, -1.0);CHKERRQ(ierr);
173 
174       ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &bnk->f, tao->gradient, bnk->D, &step, &ls_reason);CHKERRQ(ierr);
175       ierr = TaoLineSearchGetFullStepObjective(tao->linesearch, &f_full);CHKERRQ(ierr);
176       ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
177     }
178 
179     if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
180       /* Failed to find an improving point */
181       bnk->f = fold;
182       ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr);
183       ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr);
184       step = 0.0;
185       tao->reason = TAO_DIVERGED_LS_FAILURE;
186       break;
187     }
188 
189     /* Update trust region radius */
190     if (bnk->is_nash || bnk->is_stcg || bnk->is_gltr) {
191       switch(bnk->update_type) {
192       case BNK_UPDATE_STEP:
193         if (stepType == BNK_NEWTON) {
194           if (step < bnk->nu1) {
195             /* Very bad step taken; reduce radius */
196             tao->trust = bnk->omega1 * PetscMin(norm_d, tao->trust);
197           } else if (step < bnk->nu2) {
198             /* Reasonably bad step taken; reduce radius */
199             tao->trust = bnk->omega2 * PetscMin(norm_d, tao->trust);
200           } else if (step < bnk->nu3) {
201             /*  Reasonable step was taken; leave radius alone */
202             if (bnk->omega3 < 1.0) {
203               tao->trust = bnk->omega3 * PetscMin(norm_d, tao->trust);
204             } else if (bnk->omega3 > 1.0) {
205               tao->trust = PetscMax(bnk->omega3 * norm_d, tao->trust);
206             }
207           } else if (step < bnk->nu4) {
208             /*  Full step taken; increase the radius */
209             tao->trust = PetscMax(bnk->omega4 * norm_d, tao->trust);
210           } else {
211             /*  More than full step taken; increase the radius */
212             tao->trust = PetscMax(bnk->omega5 * norm_d, tao->trust);
213           }
214         } else {
215           /*  Newton step was not good; reduce the radius */
216           tao->trust = bnk->omega1 * PetscMin(norm_d, tao->trust);
217         }
218         break;
219 
220       case BNK_UPDATE_REDUCTION:
221         if (stepType == BNK_NEWTON) {
222           /*  Get predicted reduction */
223           ierr = KSPCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
224           if (prered >= 0.0) {
225             /*  The predicted reduction has the wrong sign.  This cannot */
226             /*  happen in infinite precision arithmetic.  Step should */
227             /*  be rejected! */
228             tao->trust = bnk->alpha1 * PetscMin(tao->trust, norm_d);
229           } else {
230             if (PetscIsInfOrNanReal(f_full)) {
231               tao->trust = bnk->alpha1 * PetscMin(tao->trust, norm_d);
232             } else {
233               /*  Compute and actual reduction */
234               actred = fold - f_full;
235               prered = -prered;
236               if ((PetscAbsScalar(actred) <= bnk->epsilon) &&
237                   (PetscAbsScalar(prered) <= bnk->epsilon)) {
238                 kappa = 1.0;
239               } else {
240                 kappa = actred / prered;
241               }
242 
243               /*  Accept of reject the step and update radius */
244               if (kappa < bnk->eta1) {
245                 /*  Very bad step */
246                 tao->trust = bnk->alpha1 * PetscMin(tao->trust, norm_d);
247               } else if (kappa < bnk->eta2) {
248                 /*  Marginal bad step */
249                 tao->trust = bnk->alpha2 * PetscMin(tao->trust, norm_d);
250               } else if (kappa < bnk->eta3) {
251                 /*  Reasonable step */
252                 if (bnk->alpha3 < 1.0) {
253                   tao->trust = bnk->alpha3 * PetscMin(norm_d, tao->trust);
254                 } else if (bnk->alpha3 > 1.0) {
255                   tao->trust = PetscMax(bnk->alpha3 * norm_d, tao->trust);
256                 }
257               } else if (kappa < bnk->eta4) {
258                 /*  Good step */
259                 tao->trust = PetscMax(bnk->alpha4 * norm_d, tao->trust);
260               } else {
261                 /*  Very good step */
262                 tao->trust = PetscMax(bnk->alpha5 * norm_d, tao->trust);
263               }
264             }
265           }
266         } else {
267           /*  Newton step was not good; reduce the radius */
268           tao->trust = bnk->alpha1 * PetscMin(norm_d, tao->trust);
269         }
270         break;
271 
272       default:
273         if (stepType == BNK_NEWTON) {
274           ierr = KSPCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
275           if (prered >= 0.0) {
276             /*  The predicted reduction has the wrong sign.  This cannot */
277             /*  happen in infinite precision arithmetic.  Step should */
278             /*  be rejected! */
279             tao->trust = bnk->gamma1 * PetscMin(tao->trust, norm_d);
280           } else {
281             if (PetscIsInfOrNanReal(f_full)) {
282               tao->trust = bnk->gamma1 * PetscMin(tao->trust, norm_d);
283             } else {
284               actred = fold - f_full;
285               prered = -prered;
286               if ((PetscAbsScalar(actred) <= bnk->epsilon) && (PetscAbsScalar(prered) <= bnk->epsilon)) {
287                 kappa = 1.0;
288               } else {
289                 kappa = actred / prered;
290               }
291 
292               tau_1 = bnk->theta * gdx / (bnk->theta * gdx - (1.0 - bnk->theta) * prered + actred);
293               tau_2 = bnk->theta * gdx / (bnk->theta * gdx + (1.0 + bnk->theta) * prered - actred);
294               tau_min = PetscMin(tau_1, tau_2);
295               tau_max = PetscMax(tau_1, tau_2);
296 
297               if (kappa >= 1.0 - bnk->mu1) {
298                 /*  Great agreement */
299                 if (tau_max < 1.0) {
300                   tao->trust = PetscMax(tao->trust, bnk->gamma3 * norm_d);
301                 } else if (tau_max > bnk->gamma4) {
302                   tao->trust = PetscMax(tao->trust, bnk->gamma4 * norm_d);
303                 } else {
304                   tao->trust = PetscMax(tao->trust, tau_max * norm_d);
305                 }
306               } else if (kappa >= 1.0 - bnk->mu2) {
307                 /*  Good agreement */
308 
309                 if (tau_max < bnk->gamma2) {
310                   tao->trust = bnk->gamma2 * PetscMin(tao->trust, norm_d);
311                 } else if (tau_max > bnk->gamma3) {
312                   tao->trust = PetscMax(tao->trust, bnk->gamma3 * norm_d);
313                 } else if (tau_max < 1.0) {
314                   tao->trust = tau_max * PetscMin(tao->trust, norm_d);
315                 } else {
316                   tao->trust = PetscMax(tao->trust, tau_max * norm_d);
317                 }
318               } else {
319                 /*  Not good agreement */
320                 if (tau_min > 1.0) {
321                   tao->trust = bnk->gamma2 * PetscMin(tao->trust, norm_d);
322                 } else if (tau_max < bnk->gamma1) {
323                   tao->trust = bnk->gamma1 * PetscMin(tao->trust, norm_d);
324                 } else if ((tau_min < bnk->gamma1) && (tau_max >= 1.0)) {
325                   tao->trust = bnk->gamma1 * PetscMin(tao->trust, norm_d);
326                 } else if ((tau_1 >= bnk->gamma1) && (tau_1 < 1.0) && ((tau_2 < bnk->gamma1) || (tau_2 >= 1.0))) {
327                   tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
328                 } else if ((tau_2 >= bnk->gamma1) && (tau_2 < 1.0) && ((tau_1 < bnk->gamma1) || (tau_2 >= 1.0))) {
329                   tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
330                 } else {
331                   tao->trust = tau_max * PetscMin(tao->trust, norm_d);
332                 }
333               }
334             }
335           }
336         } else {
337           /*  Newton step was not good; reduce the radius */
338           tao->trust = bnk->gamma1 * PetscMin(norm_d, tao->trust);
339         }
340         break;
341       }
342 
343       /*  The radius may have been increased; modify if it is too large */
344       tao->trust = PetscMin(tao->trust, bnk->max_radius);
345     }
346 
347     /*  Check for termination */
348     ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr);
349     if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number");
350     needH = 1;
351     ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
352     ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,step);CHKERRQ(ierr);
353     ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
354   }
355   PetscFunctionReturn(0);
356 }
357 
358 PETSC_EXTERN PetscErrorCode TaoCreate_BNLS(Tao tao)
359 {
360   PetscErrorCode ierr;
361 
362   PetscFunctionBegin;
363   ierr = TaoCreate_BNK(tao);CHKERRQ(ierr);
364   tao->ops->solve=TaoSolve_BNLS;
365   PetscFunctionReturn(0);
366 }