xref: /petsc/src/tao/unconstrained/impls/ntl/ntl.c (revision 03047865b8d8757cf1cf9cda45785c1537b01dc1) 
1 #include <../src/tao/unconstrained/impls/ntl/ntlimpl.h>
2 
3 #include <petscksp.h>
4 
5 #define NTL_INIT_CONSTANT      0
6 #define NTL_INIT_DIRECTION     1
7 #define NTL_INIT_INTERPOLATION 2
8 #define NTL_INIT_TYPES         3
9 
10 #define NTL_UPDATE_REDUCTION     0
11 #define NTL_UPDATE_INTERPOLATION 1
12 #define NTL_UPDATE_TYPES         2
13 
14 static const char *NTL_INIT[64] = {"constant", "direction", "interpolation"};
15 
16 static const char *NTL_UPDATE[64] = {"reduction", "interpolation"};
17 
18 /* Implements Newton's Method with a trust-region, line-search approach for
19    solving unconstrained minimization problems.  A More'-Thuente line search
20    is used to guarantee that the bfgs preconditioner remains positive
21    definite. */
22 
23 #define NTL_NEWTON          0
24 #define NTL_BFGS            1
25 #define NTL_SCALED_GRADIENT 2
26 #define NTL_GRADIENT        3
27 
TaoSolve_NTL(Tao tao)28 static PetscErrorCode TaoSolve_NTL(Tao tao)
29 {
30   TAO_NTL                     *tl = (TAO_NTL *)tao->data;
31   KSPType                      ksp_type;
32   PetscBool                    is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set;
33   KSPConvergedReason           ksp_reason;
34   PC                           pc;
35   TaoLineSearchConvergedReason ls_reason;
36 
37   PetscReal fmin, ftrial, prered, actred, kappa, sigma;
38   PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
39   PetscReal f, fold, gdx, gnorm;
40   PetscReal step = 1.0;
41 
42   PetscReal norm_d = 0.0;
43   PetscInt  stepType;
44   PetscInt  its;
45 
46   PetscInt bfgsUpdates = 0;
47   PetscInt needH;
48 
49   PetscInt i_max = 5;
50   PetscInt j_max = 1;
51   PetscInt i, j, n, N;
52 
53   PetscInt tr_reject;
54 
55   PetscFunctionBegin;
56   if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n"));
57 
58   PetscCall(KSPGetType(tao->ksp, &ksp_type));
59   PetscCall(PetscStrcmp(ksp_type, KSPNASH, &is_nash));
60   PetscCall(PetscStrcmp(ksp_type, KSPSTCG, &is_stcg));
61   PetscCall(PetscStrcmp(ksp_type, KSPGLTR, &is_gltr));
62   PetscCheck(is_nash || is_stcg || is_gltr, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAO_NTR requires nash, stcg, or gltr for the KSP");
63 
64   /* Initialize the radius and modify if it is too large or small */
65   tao->trust = tao->trust0;
66   tao->trust = PetscMax(tao->trust, tl->min_radius);
67   tao->trust = PetscMin(tao->trust, tl->max_radius);
68 
69   /* Allocate the vectors needed for the BFGS approximation */
70   PetscCall(KSPGetPC(tao->ksp, &pc));
71   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs));
72   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi));
73   if (is_bfgs) {
74     tl->bfgs_pre = pc;
75     PetscCall(PCLMVMGetMatLMVM(tl->bfgs_pre, &tl->M));
76     PetscCall(VecGetLocalSize(tao->solution, &n));
77     PetscCall(VecGetSize(tao->solution, &N));
78     PetscCall(MatSetSizes(tl->M, n, n, N, N));
79     PetscCall(MatLMVMAllocate(tl->M, tao->solution, tao->gradient));
80     PetscCall(MatIsSymmetricKnown(tl->M, &sym_set, &is_symmetric));
81     PetscCheck(sym_set && is_symmetric, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
82   } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc, PETSC_TRUE));
83 
84   /* Check convergence criteria */
85   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
86   PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
87   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
88   needH = 1;
89 
90   tao->reason = TAO_CONTINUE_ITERATING;
91   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
92   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
93   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
94   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
95 
96   /* Initialize trust-region radius */
97   switch (tl->init_type) {
98   case NTL_INIT_CONSTANT:
99     /* Use the initial radius specified */
100     break;
101 
102   case NTL_INIT_INTERPOLATION:
103     /* Use the initial radius specified */
104     max_radius = 0.0;
105 
106     for (j = 0; j < j_max; ++j) {
107       fmin  = f;
108       sigma = 0.0;
109 
110       if (needH) {
111         PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
112         needH = 0;
113       }
114 
115       for (i = 0; i < i_max; ++i) {
116         PetscCall(VecCopy(tao->solution, tl->W));
117         PetscCall(VecAXPY(tl->W, -tao->trust / gnorm, tao->gradient));
118 
119         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
120         if (PetscIsInfOrNanReal(ftrial)) {
121           tau = tl->gamma1_i;
122         } else {
123           if (ftrial < fmin) {
124             fmin  = ftrial;
125             sigma = -tao->trust / gnorm;
126           }
127 
128           PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection));
129           PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered));
130 
131           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
132           actred = f - ftrial;
133           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
134             kappa = 1.0;
135           } else {
136             kappa = actred / prered;
137           }
138 
139           tau_1   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
140           tau_2   = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
141           tau_min = PetscMin(tau_1, tau_2);
142           tau_max = PetscMax(tau_1, tau_2);
143 
144           if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu1_i) {
145             /* Great agreement */
146             max_radius = PetscMax(max_radius, tao->trust);
147 
148             if (tau_max < 1.0) {
149               tau = tl->gamma3_i;
150             } else if (tau_max > tl->gamma4_i) {
151               tau = tl->gamma4_i;
152             } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
153               tau = tau_1;
154             } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
155               tau = tau_2;
156             } else {
157               tau = tau_max;
158             }
159           } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu2_i) {
160             /* Good agreement */
161             max_radius = PetscMax(max_radius, tao->trust);
162 
163             if (tau_max < tl->gamma2_i) {
164               tau = tl->gamma2_i;
165             } else if (tau_max > tl->gamma3_i) {
166               tau = tl->gamma3_i;
167             } else {
168               tau = tau_max;
169             }
170           } else {
171             /* Not good agreement */
172             if (tau_min > 1.0) {
173               tau = tl->gamma2_i;
174             } else if (tau_max < tl->gamma1_i) {
175               tau = tl->gamma1_i;
176             } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
177               tau = tl->gamma1_i;
178             } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
179               tau = tau_1;
180             } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
181               tau = tau_2;
182             } else {
183               tau = tau_max;
184             }
185           }
186         }
187         tao->trust = tau * tao->trust;
188       }
189 
190       if (fmin < f) {
191         f = fmin;
192         PetscCall(VecAXPY(tao->solution, sigma, tao->gradient));
193         PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
194 
195         PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
196         PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
197         needH = 1;
198 
199         PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
200         PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
201         PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
202         if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
203       }
204     }
205     tao->trust = PetscMax(tao->trust, max_radius);
206 
207     /* Modify the radius if it is too large or small */
208     tao->trust = PetscMax(tao->trust, tl->min_radius);
209     tao->trust = PetscMin(tao->trust, tl->max_radius);
210     break;
211 
212   default:
213     /* Norm of the first direction will initialize radius */
214     tao->trust = 0.0;
215     break;
216   }
217 
218   /* Set counter for gradient/reset steps */
219   tl->ntrust = 0;
220   tl->newt   = 0;
221   tl->bfgs   = 0;
222   tl->grad   = 0;
223 
224   /* Have not converged; continue with Newton method */
225   while (tao->reason == TAO_CONTINUE_ITERATING) {
226     /* Call general purpose update function */
227     if (tao->ops->update) {
228       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
229       PetscCall(TaoComputeObjective(tao, tao->solution, &f));
230     }
231     ++tao->niter;
232     tao->ksp_its = 0;
233     /* Compute the Hessian */
234     if (needH) PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
235 
236     if (tl->bfgs_pre) {
237       /* Update the limited memory preconditioner */
238       PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
239       ++bfgsUpdates;
240     }
241     PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
242     /* Solve the Newton system of equations */
243     PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
244     PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
245     PetscCall(KSPGetIterationNumber(tao->ksp, &its));
246     tao->ksp_its += its;
247     tao->ksp_tot_its += its;
248     PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));
249 
250     if (0.0 == tao->trust) {
251       /* Radius was uninitialized; use the norm of the direction */
252       if (norm_d > 0.0) {
253         tao->trust = norm_d;
254 
255         /* Modify the radius if it is too large or small */
256         tao->trust = PetscMax(tao->trust, tl->min_radius);
257         tao->trust = PetscMin(tao->trust, tl->max_radius);
258       } else {
259         /* The direction was bad; set radius to default value and re-solve
260            the trust-region subproblem to get a direction */
261         tao->trust = tao->trust0;
262 
263         /* Modify the radius if it is too large or small */
264         tao->trust = PetscMax(tao->trust, tl->min_radius);
265         tao->trust = PetscMin(tao->trust, tl->max_radius);
266 
267         PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
268         PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
269         PetscCall(KSPGetIterationNumber(tao->ksp, &its));
270         tao->ksp_its += its;
271         tao->ksp_tot_its += its;
272         PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));
273 
274         PetscCheck(norm_d != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero");
275       }
276     }
277 
278     PetscCall(VecScale(tao->stepdirection, -1.0));
279     PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason));
280     if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tl->bfgs_pre)) {
281       /* Preconditioner is numerically indefinite; reset the
282          approximate if using BFGS preconditioning. */
283       PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
284       PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
285       bfgsUpdates = 1;
286     }
287 
288     /* Check trust-region reduction conditions */
289     tr_reject = 0;
290     if (NTL_UPDATE_REDUCTION == tl->update_type) {
291       /* Get predicted reduction */
292       PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
293       if (prered >= 0.0) {
294         /* The predicted reduction has the wrong sign.  This cannot
295            happen in infinite precision arithmetic.  Step should
296            be rejected! */
297         tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
298         tr_reject  = 1;
299       } else {
300         /* Compute trial step and function value */
301         PetscCall(VecCopy(tao->solution, tl->W));
302         PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
303         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
304 
305         if (PetscIsInfOrNanReal(ftrial)) {
306           tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
307           tr_reject  = 1;
308         } else {
309           /* Compute and actual reduction */
310           actred = f - ftrial;
311           prered = -prered;
312           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
313             kappa = 1.0;
314           } else {
315             kappa = actred / prered;
316           }
317 
318           /* Accept of reject the step and update radius */
319           if (kappa < tl->eta1) {
320             /* Reject the step */
321             tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
322             tr_reject  = 1;
323           } else {
324             /* Accept the step */
325             if (kappa < tl->eta2) {
326               /* Marginal bad step */
327               tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
328             } else if (kappa < tl->eta3) {
329               /* Reasonable step */
330               tao->trust = tl->alpha3 * tao->trust;
331             } else if (kappa < tl->eta4) {
332               /* Good step */
333               tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
334             } else {
335               /* Very good step */
336               tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
337             }
338           }
339         }
340       }
341     } else {
342       /* Get predicted reduction */
343       PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
344       if (prered >= 0.0) {
345         /* The predicted reduction has the wrong sign.  This cannot
346            happen in infinite precision arithmetic.  Step should
347            be rejected! */
348         tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
349         tr_reject  = 1;
350       } else {
351         PetscCall(VecCopy(tao->solution, tl->W));
352         PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
353         PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
354         if (PetscIsInfOrNanReal(ftrial)) {
355           tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
356           tr_reject  = 1;
357         } else {
358           PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx));
359 
360           actred = f - ftrial;
361           prered = -prered;
362           if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
363             kappa = 1.0;
364           } else {
365             kappa = actred / prered;
366           }
367 
368           tau_1   = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
369           tau_2   = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
370           tau_min = PetscMin(tau_1, tau_2);
371           tau_max = PetscMax(tau_1, tau_2);
372 
373           if (kappa >= 1.0 - tl->mu1) {
374             /* Great agreement; accept step and update radius */
375             if (tau_max < 1.0) {
376               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
377             } else if (tau_max > tl->gamma4) {
378               tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
379             } else {
380               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
381             }
382           } else if (kappa >= 1.0 - tl->mu2) {
383             /* Good agreement */
384 
385             if (tau_max < tl->gamma2) {
386               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
387             } else if (tau_max > tl->gamma3) {
388               tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
389             } else if (tau_max < 1.0) {
390               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
391             } else {
392               tao->trust = PetscMax(tao->trust, tau_max * norm_d);
393             }
394           } else {
395             /* Not good agreement */
396             if (tau_min > 1.0) {
397               tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
398             } else if (tau_max < tl->gamma1) {
399               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
400             } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
401               tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
402             } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
403               tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
404             } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
405               tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
406             } else {
407               tao->trust = tau_max * PetscMin(tao->trust, norm_d);
408             }
409             tr_reject = 1;
410           }
411         }
412       }
413     }
414 
415     if (tr_reject) {
416       /* The trust-region constraints rejected the step.  Apply a linesearch.
417          Check for descent direction. */
418       PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
419       if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
420         /* Newton step is not descent or direction produced infinity or NaN */
421 
422         if (!tl->bfgs_pre) {
423           /* We don't have the bfgs matrix around and updated
424              Must use gradient direction in this case */
425           PetscCall(VecCopy(tao->gradient, tao->stepdirection));
426           PetscCall(VecScale(tao->stepdirection, -1.0));
427           ++tl->grad;
428           stepType = NTL_GRADIENT;
429         } else {
430           /* Attempt to use the BFGS direction */
431           PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
432           PetscCall(VecScale(tao->stepdirection, -1.0));
433 
434           /* Check for success (descent direction) */
435           PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
436           if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
437             /* BFGS direction is not descent or direction produced not a number
438                We can assert bfgsUpdates > 1 in this case because
439                the first solve produces the scaled gradient direction,
440                which is guaranteed to be descent */
441 
442             /* Use steepest descent direction (scaled) */
443             PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
444             PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
445             PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
446             PetscCall(VecScale(tao->stepdirection, -1.0));
447 
448             bfgsUpdates = 1;
449             ++tl->grad;
450             stepType = NTL_GRADIENT;
451           } else {
452             PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
453             if (1 == bfgsUpdates) {
454               /* The first BFGS direction is always the scaled gradient */
455               ++tl->grad;
456               stepType = NTL_GRADIENT;
457             } else {
458               ++tl->bfgs;
459               stepType = NTL_BFGS;
460             }
461           }
462         }
463       } else {
464         /* Computed Newton step is descent */
465         ++tl->newt;
466         stepType = NTL_NEWTON;
467       }
468 
469       /* Perform the linesearch */
470       fold = f;
471       PetscCall(VecCopy(tao->solution, tl->Xold));
472       PetscCall(VecCopy(tao->gradient, tl->Gold));
473 
474       step = 1.0;
475       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
476       PetscCall(TaoAddLineSearchCounts(tao));
477 
478       while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */
479         /* Linesearch failed */
480         f = fold;
481         PetscCall(VecCopy(tl->Xold, tao->solution));
482         PetscCall(VecCopy(tl->Gold, tao->gradient));
483 
484         switch (stepType) {
485         case NTL_NEWTON:
486           /* Failed to obtain acceptable iterate with Newton step */
487 
488           if (tl->bfgs_pre) {
489             /* We don't have the bfgs matrix around and being updated
490                Must use gradient direction in this case */
491             PetscCall(VecCopy(tao->gradient, tao->stepdirection));
492             ++tl->grad;
493             stepType = NTL_GRADIENT;
494           } else {
495             /* Attempt to use the BFGS direction */
496             PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
497 
498             /* Check for success (descent direction) */
499             PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
500             if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
501               /* BFGS direction is not descent or direction produced
502                  not a number.  We can assert bfgsUpdates > 1 in this case
503                  Use steepest descent direction (scaled) */
504               PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
505               PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
506               PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
507 
508               bfgsUpdates = 1;
509               ++tl->grad;
510               stepType = NTL_GRADIENT;
511             } else {
512               PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
513               if (1 == bfgsUpdates) {
514                 /* The first BFGS direction is always the scaled gradient */
515                 ++tl->grad;
516                 stepType = NTL_GRADIENT;
517               } else {
518                 ++tl->bfgs;
519                 stepType = NTL_BFGS;
520               }
521             }
522           }
523           break;
524 
525         case NTL_BFGS:
526           /* Can only enter if pc_type == NTL_PC_BFGS
527              Failed to obtain acceptable iterate with BFGS step
528              Attempt to use the scaled gradient direction */
529           PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
530           PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
531           PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
532 
533           bfgsUpdates = 1;
534           ++tl->grad;
535           stepType = NTL_GRADIENT;
536           break;
537         }
538         PetscCall(VecScale(tao->stepdirection, -1.0));
539 
540         /* This may be incorrect; linesearch has values for stepmax and stepmin
541            that should be reset. */
542         step = 1.0;
543         PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
544         PetscCall(TaoAddLineSearchCounts(tao));
545       }
546 
547       if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
548         /* Failed to find an improving point */
549         f = fold;
550         PetscCall(VecCopy(tl->Xold, tao->solution));
551         PetscCall(VecCopy(tl->Gold, tao->gradient));
552         tao->trust  = 0.0;
553         step        = 0.0;
554         tao->reason = TAO_DIVERGED_LS_FAILURE;
555         break;
556       } else if (stepType == NTL_NEWTON) {
557         if (step < tl->nu1) {
558           /* Very bad step taken; reduce radius */
559           tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
560         } else if (step < tl->nu2) {
561           /* Reasonably bad step taken; reduce radius */
562           tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
563         } else if (step < tl->nu3) {
564           /* Reasonable step was taken; leave radius alone */
565           if (tl->omega3 < 1.0) {
566             tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
567           } else if (tl->omega3 > 1.0) {
568             tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
569           }
570         } else if (step < tl->nu4) {
571           /* Full step taken; increase the radius */
572           tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
573         } else {
574           /* More than full step taken; increase the radius */
575           tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
576         }
577       } else {
578         /* Newton step was not good; reduce the radius */
579         tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
580       }
581     } else {
582       /* Trust-region step is accepted */
583       PetscCall(VecCopy(tl->W, tao->solution));
584       f = ftrial;
585       PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
586       ++tl->ntrust;
587     }
588 
589     /* The radius may have been increased; modify if it is too large */
590     tao->trust = PetscMin(tao->trust, tl->max_radius);
591 
592     /* Check for converged */
593     PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
594     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
595     needH = 1;
596 
597     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
598     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
599     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
600   }
601   PetscFunctionReturn(PETSC_SUCCESS);
602 }
603 
TaoSetUp_NTL(Tao tao)604 static PetscErrorCode TaoSetUp_NTL(Tao tao)
605 {
606   TAO_NTL *tl = (TAO_NTL *)tao->data;
607 
608   PetscFunctionBegin;
609   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
610   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
611   if (!tl->W) PetscCall(VecDuplicate(tao->solution, &tl->W));
612   if (!tl->Xold) PetscCall(VecDuplicate(tao->solution, &tl->Xold));
613   if (!tl->Gold) PetscCall(VecDuplicate(tao->solution, &tl->Gold));
614   tl->bfgs_pre = NULL;
615   tl->M        = NULL;
616   PetscFunctionReturn(PETSC_SUCCESS);
617 }
618 
TaoDestroy_NTL(Tao tao)619 static PetscErrorCode TaoDestroy_NTL(Tao tao)
620 {
621   TAO_NTL *tl = (TAO_NTL *)tao->data;
622 
623   PetscFunctionBegin;
624   if (tao->setupcalled) {
625     PetscCall(VecDestroy(&tl->W));
626     PetscCall(VecDestroy(&tl->Xold));
627     PetscCall(VecDestroy(&tl->Gold));
628   }
629   PetscCall(KSPDestroy(&tao->ksp));
630   PetscCall(PetscFree(tao->data));
631   PetscFunctionReturn(PETSC_SUCCESS);
632 }
633 
TaoSetFromOptions_NTL(Tao tao,PetscOptionItems PetscOptionsObject)634 static PetscErrorCode TaoSetFromOptions_NTL(Tao tao, PetscOptionItems PetscOptionsObject)
635 {
636   TAO_NTL *tl = (TAO_NTL *)tao->data;
637 
638   PetscFunctionBegin;
639   PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region with line search method for unconstrained optimization");
640   PetscCall(PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, NULL));
641   PetscCall(PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, NULL));
642   PetscCall(PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, NULL));
643   PetscCall(PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, NULL));
644   PetscCall(PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, NULL));
645   PetscCall(PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, NULL));
646   PetscCall(PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, NULL));
647   PetscCall(PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, NULL));
648   PetscCall(PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, NULL));
649   PetscCall(PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, NULL));
650   PetscCall(PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, NULL));
651   PetscCall(PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, NULL));
652   PetscCall(PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, NULL));
653   PetscCall(PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, NULL));
654   PetscCall(PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, NULL));
655   PetscCall(PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, NULL));
656   PetscCall(PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, NULL));
657   PetscCall(PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, NULL));
658   PetscCall(PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, NULL));
659   PetscCall(PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, NULL));
660   PetscCall(PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, NULL));
661   PetscCall(PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, NULL));
662   PetscCall(PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, NULL));
663   PetscCall(PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, NULL));
664   PetscCall(PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, NULL));
665   PetscCall(PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, NULL));
666   PetscCall(PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, NULL));
667   PetscCall(PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, NULL));
668   PetscCall(PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, NULL));
669   PetscCall(PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, NULL));
670   PetscCall(PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, NULL));
671   PetscCall(PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, NULL));
672   PetscCall(PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, NULL));
673   PetscCall(PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, NULL));
674   PetscCall(PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, NULL));
675   PetscCall(PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, NULL));
676   PetscCall(PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, NULL));
677   PetscOptionsHeadEnd();
678   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
679   PetscCall(KSPSetFromOptions(tao->ksp));
680   PetscFunctionReturn(PETSC_SUCCESS);
681 }
682 
TaoView_NTL(Tao tao,PetscViewer viewer)683 static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer)
684 {
685   TAO_NTL  *tl = (TAO_NTL *)tao->data;
686   PetscBool isascii;
687 
688   PetscFunctionBegin;
689   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
690   if (isascii) {
691     PetscCall(PetscViewerASCIIPushTab(viewer));
692     PetscCall(PetscViewerASCIIPrintf(viewer, "Trust-region steps: %" PetscInt_FMT "\n", tl->ntrust));
693     PetscCall(PetscViewerASCIIPrintf(viewer, "Newton search steps: %" PetscInt_FMT "\n", tl->newt));
694     PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS search steps: %" PetscInt_FMT "\n", tl->bfgs));
695     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient search steps: %" PetscInt_FMT "\n", tl->grad));
696     PetscCall(PetscViewerASCIIPopTab(viewer));
697   }
698   PetscFunctionReturn(PETSC_SUCCESS);
699 }
700 
701 /*MC
702   TAONTL - Newton's method with trust region globalization and line search fallback.
703   At each iteration, the Newton trust region method solves the system for d
704   and performs a line search in the d direction:
705 
706             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k
707 
708   Options Database Keys:
709 + -tao_ntl_init_type - "constant","direction","interpolation"
710 . -tao_ntl_update_type - "reduction","interpolation"
711 . -tao_ntl_min_radius - lower bound on trust region radius
712 . -tao_ntl_max_radius - upper bound on trust region radius
713 . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
714 . -tao_ntl_mu1_i - mu1 interpolation init factor
715 . -tao_ntl_mu2_i - mu2 interpolation init factor
716 . -tao_ntl_gamma1_i - gamma1 interpolation init factor
717 . -tao_ntl_gamma2_i - gamma2 interpolation init factor
718 . -tao_ntl_gamma3_i - gamma3 interpolation init factor
719 . -tao_ntl_gamma4_i - gamma4 interpolation init factor
720 . -tao_ntl_theta_i - theta1 interpolation init factor
721 . -tao_ntl_eta1 - eta1 reduction update factor
722 . -tao_ntl_eta2 - eta2 reduction update factor
723 . -tao_ntl_eta3 - eta3 reduction update factor
724 . -tao_ntl_eta4 - eta4 reduction update factor
725 . -tao_ntl_alpha1 - alpha1 reduction update factor
726 . -tao_ntl_alpha2 - alpha2 reduction update factor
727 . -tao_ntl_alpha3 - alpha3 reduction update factor
728 . -tao_ntl_alpha4 - alpha4 reduction update factor
729 . -tao_ntl_alpha4 - alpha4 reduction update factor
730 . -tao_ntl_mu1 - mu1 interpolation update
731 . -tao_ntl_mu2 - mu2 interpolation update
732 . -tao_ntl_gamma1 - gamma1 interpolcation update
733 . -tao_ntl_gamma2 - gamma2 interpolcation update
734 . -tao_ntl_gamma3 - gamma3 interpolcation update
735 . -tao_ntl_gamma4 - gamma4 interpolation update
736 - -tao_ntl_theta - theta1 interpolation update
737 
738   Level: beginner
739 M*/
TaoCreate_NTL(Tao tao)740 PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao)
741 {
742   TAO_NTL    *tl;
743   const char *morethuente_type = TAOLINESEARCHMT;
744 
745   PetscFunctionBegin;
746   PetscCall(PetscNew(&tl));
747   tao->ops->setup          = TaoSetUp_NTL;
748   tao->ops->solve          = TaoSolve_NTL;
749   tao->ops->view           = TaoView_NTL;
750   tao->ops->setfromoptions = TaoSetFromOptions_NTL;
751   tao->ops->destroy        = TaoDestroy_NTL;
752 
753   /* Override default settings (unless already changed) */
754   PetscCall(TaoParametersInitialize(tao));
755   PetscObjectParameterSetDefault(tao, max_it, 50);
756   PetscObjectParameterSetDefault(tao, trust0, 100.0);
757 
758   tao->data = (void *)tl;
759 
760   /* Default values for trust-region radius update based on steplength */
761   tl->nu1 = 0.25;
762   tl->nu2 = 0.50;
763   tl->nu3 = 1.00;
764   tl->nu4 = 1.25;
765 
766   tl->omega1 = 0.25;
767   tl->omega2 = 0.50;
768   tl->omega3 = 1.00;
769   tl->omega4 = 2.00;
770   tl->omega5 = 4.00;
771 
772   /* Default values for trust-region radius update based on reduction */
773   tl->eta1 = 1.0e-4;
774   tl->eta2 = 0.25;
775   tl->eta3 = 0.50;
776   tl->eta4 = 0.90;
777 
778   tl->alpha1 = 0.25;
779   tl->alpha2 = 0.50;
780   tl->alpha3 = 1.00;
781   tl->alpha4 = 2.00;
782   tl->alpha5 = 4.00;
783 
784   /* Default values for trust-region radius update based on interpolation */
785   tl->mu1 = 0.10;
786   tl->mu2 = 0.50;
787 
788   tl->gamma1 = 0.25;
789   tl->gamma2 = 0.50;
790   tl->gamma3 = 2.00;
791   tl->gamma4 = 4.00;
792 
793   tl->theta = 0.05;
794 
795   /* Default values for trust region initialization based on interpolation */
796   tl->mu1_i = 0.35;
797   tl->mu2_i = 0.50;
798 
799   tl->gamma1_i = 0.0625;
800   tl->gamma2_i = 0.5;
801   tl->gamma3_i = 2.0;
802   tl->gamma4_i = 5.0;
803 
804   tl->theta_i = 0.25;
805 
806   /* Remaining parameters */
807   tl->min_radius = 1.0e-10;
808   tl->max_radius = 1.0e10;
809   tl->epsilon    = 1.0e-6;
810 
811   tl->init_type   = NTL_INIT_INTERPOLATION;
812   tl->update_type = NTL_UPDATE_REDUCTION;
813 
814   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
815   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
816   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
817   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
818   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
819   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
820   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
821   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
822   PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_ntl_"));
823   PetscCall(KSPSetType(tao->ksp, KSPSTCG));
824   PetscFunctionReturn(PETSC_SUCCESS);
825 }
826