xref: /petsc/src/tao/unconstrained/impls/ntr/ntr.c (revision 68ed366e2c7210eabb52d7e38c0e86223c5dc07d) 
1 #include <../src/tao/unconstrained/impls/ntr/ntrimpl.h>
2 
3 #include <petscksp.h>
4 
5 #define NTR_INIT_CONSTANT         0
6 #define NTR_INIT_DIRECTION        1
7 #define NTR_INIT_INTERPOLATION    2
8 #define NTR_INIT_TYPES            3
9 
10 #define NTR_UPDATE_REDUCTION      0
11 #define NTR_UPDATE_INTERPOLATION  1
12 #define NTR_UPDATE_TYPES          2
13 
14 static const char *NTR_INIT[64] = {"constant","direction","interpolation"};
15 
16 static const char *NTR_UPDATE[64] = {"reduction","interpolation"};
17 
18 /*
19    TaoSolve_NTR - Implements Newton's Method with a trust region approach
20    for solving unconstrained minimization problems.
21 
22    The basic algorithm is taken from MINPACK-2 (dstrn).
23 
24    TaoSolve_NTR computes a local minimizer of a twice differentiable function
25    f by applying a trust region variant of Newton's method.  At each stage
26    of the algorithm, we use the prconditioned conjugate gradient method to
27    determine an approximate minimizer of the quadratic equation
28 
29         q(s) = <s, Hs + g>
30 
31    subject to the trust region constraint
32 
33         || s ||_M <= radius,
34 
35    where radius is the trust region radius and M is a symmetric positive
36    definite matrix (the preconditioner).  Here g is the gradient and H
37    is the Hessian matrix.
38 
39    Note:  TaoSolve_NTR MUST use the iterative solver KSPCGNASH, KSPCGSTCG,
40           or KSPCGGLTR.  Thus, we set KSPCGNASH, KSPCGSTCG, or KSPCGGLTR in this
41           routine regardless of what the user may have previously specified.
42 */
43 static PetscErrorCode TaoSolve_NTR(Tao tao)
44 {
45   TAO_NTR            *tr = (TAO_NTR *)tao->data;
46   KSPType            ksp_type;
47   PetscBool          is_nash,is_stcg,is_gltr,is_bfgs,is_jacobi,is_symmetric,sym_set;
48   KSPConvergedReason ksp_reason;
49   PC                 pc;
50   PetscReal          fmin, ftrial, prered, actred, kappa, sigma, beta;
51   PetscReal          tau, tau_1, tau_2, tau_max, tau_min, max_radius;
52   PetscReal          f, gnorm;
53 
54   PetscReal          norm_d;
55   PetscErrorCode     ierr;
56   PetscInt           bfgsUpdates = 0;
57   PetscInt           needH;
58 
59   PetscInt           i_max = 5;
60   PetscInt           j_max = 1;
61   PetscInt           i, j, N, n, its;
62 
63   PetscFunctionBegin;
64   if (tao->XL || tao->XU || tao->ops->computebounds) {
65     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");CHKERRQ(ierr);
66   }
67 
68   ierr = KSPGetType(tao->ksp,&ksp_type);CHKERRQ(ierr);
69   ierr = PetscStrcmp(ksp_type,KSPCGNASH,&is_nash);CHKERRQ(ierr);
70   ierr = PetscStrcmp(ksp_type,KSPCGSTCG,&is_stcg);CHKERRQ(ierr);
71   ierr = PetscStrcmp(ksp_type,KSPCGGLTR,&is_gltr);CHKERRQ(ierr);
72   if (!is_nash && !is_stcg && !is_gltr) {
73     SETERRQ(PETSC_COMM_SELF,1,"TAO_NTR requires nash, stcg, or gltr for the KSP");
74   }
75 
76   /* Initialize the radius and modify if it is too large or small */
77   tao->trust = tao->trust0;
78   tao->trust = PetscMax(tao->trust, tr->min_radius);
79   tao->trust = PetscMin(tao->trust, tr->max_radius);
80 
81 /* Allocate the vectors needed for the BFGS approximation */
82 ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr);
83 ierr = PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs);CHKERRQ(ierr);
84 ierr = PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi);CHKERRQ(ierr);
85 if (is_bfgs) {
86   tr->bfgs_pre = pc;
87   ierr = PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M);CHKERRQ(ierr);
88   ierr = VecGetLocalSize(tao->solution, &n);CHKERRQ(ierr);
89   ierr = VecGetSize(tao->solution, &N);CHKERRQ(ierr);
90   ierr = MatSetSizes(tr->M, n, n, N, N);CHKERRQ(ierr);
91   ierr = MatLMVMAllocate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr);
92   ierr = MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric);CHKERRQ(ierr);
93   if (!sym_set || !is_symmetric) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
94 } else if (is_jacobi) {
95   ierr = PCJacobiSetUseAbs(pc,PETSC_TRUE);CHKERRQ(ierr);
96 }
97 
98   /* Check convergence criteria */
99   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
100   ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
101   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Inf or NaN");
102   needH = 1;
103 
104   tao->reason = TAO_CONTINUE_ITERATING;
105   ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
106   ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);CHKERRQ(ierr);
107   ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
108   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
109 
110   /*  Initialize trust-region radius */
111   switch(tr->init_type) {
112   case NTR_INIT_CONSTANT:
113     /*  Use the initial radius specified */
114     break;
115 
116   case NTR_INIT_INTERPOLATION:
117     /*  Use the initial radius specified */
118     max_radius = 0.0;
119 
120     for (j = 0; j < j_max; ++j) {
121       fmin = f;
122       sigma = 0.0;
123 
124       if (needH) {
125         ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
126         needH = 0;
127       }
128 
129       for (i = 0; i < i_max; ++i) {
130 
131         ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr);
132         ierr = VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);CHKERRQ(ierr);
133         ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
134 
135         if (PetscIsInfOrNanReal(ftrial)) {
136           tau = tr->gamma1_i;
137         }
138         else {
139           if (ftrial < fmin) {
140             fmin = ftrial;
141             sigma = -tao->trust / gnorm;
142           }
143 
144           ierr = MatMult(tao->hessian, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
145           ierr = VecDot(tao->gradient, tao->stepdirection, &prered);CHKERRQ(ierr);
146 
147           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
148           actred = f - ftrial;
149           if ((PetscAbsScalar(actred) <= tr->epsilon) &&
150               (PetscAbsScalar(prered) <= tr->epsilon)) {
151             kappa = 1.0;
152           }
153           else {
154             kappa = actred / prered;
155           }
156 
157           tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
158           tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
159           tau_min = PetscMin(tau_1, tau_2);
160           tau_max = PetscMax(tau_1, tau_2);
161 
162           if (PetscAbsScalar(kappa - 1.0) <= tr->mu1_i) {
163             /*  Great agreement */
164             max_radius = PetscMax(max_radius, tao->trust);
165 
166             if (tau_max < 1.0) {
167               tau = tr->gamma3_i;
168             }
169             else if (tau_max > tr->gamma4_i) {
170               tau = tr->gamma4_i;
171             }
172             else {
173               tau = tau_max;
174             }
175           }
176           else if (PetscAbsScalar(kappa - 1.0) <= tr->mu2_i) {
177             /*  Good agreement */
178             max_radius = PetscMax(max_radius, tao->trust);
179 
180             if (tau_max < tr->gamma2_i) {
181               tau = tr->gamma2_i;
182             }
183             else if (tau_max > tr->gamma3_i) {
184               tau = tr->gamma3_i;
185             }
186             else {
187               tau = tau_max;
188             }
189           }
190           else {
191             /*  Not good agreement */
192             if (tau_min > 1.0) {
193               tau = tr->gamma2_i;
194             }
195             else if (tau_max < tr->gamma1_i) {
196               tau = tr->gamma1_i;
197             }
198             else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
199               tau = tr->gamma1_i;
200             }
201             else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) &&
202                      ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
203               tau = tau_1;
204             }
205             else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) &&
206                      ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
207               tau = tau_2;
208             }
209             else {
210               tau = tau_max;
211             }
212           }
213         }
214         tao->trust = tau * tao->trust;
215       }
216 
217       if (fmin < f) {
218         f = fmin;
219         ierr = VecAXPY(tao->solution, sigma, tao->gradient);CHKERRQ(ierr);
220         ierr = TaoComputeGradient(tao,tao->solution, tao->gradient);CHKERRQ(ierr);
221 
222         ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
223 
224         if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
225         needH = 1;
226 
227         ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
228         ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0);CHKERRQ(ierr);
229         ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
230         if (tao->reason != TAO_CONTINUE_ITERATING) {
231           PetscFunctionReturn(0);
232         }
233       }
234     }
235     tao->trust = PetscMax(tao->trust, max_radius);
236 
237     /*  Modify the radius if it is too large or small */
238     tao->trust = PetscMax(tao->trust, tr->min_radius);
239     tao->trust = PetscMin(tao->trust, tr->max_radius);
240     break;
241 
242   default:
243     /*  Norm of the first direction will initialize radius */
244     tao->trust = 0.0;
245     break;
246   }
247 
248   /* Have not converged; continue with Newton method */
249   while (tao->reason == TAO_CONTINUE_ITERATING) {
250     ++tao->niter;
251     tao->ksp_its=0;
252     /* Compute the Hessian */
253     if (needH) {
254       ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr);
255       needH = 0;
256     }
257 
258     if (tr->bfgs_pre) {
259       /* Update the limited memory preconditioner */
260       ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr);
261       ++bfgsUpdates;
262     }
263 
264     while (tao->reason == TAO_CONTINUE_ITERATING) {
265       ierr = KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);CHKERRQ(ierr);
266 
267       /* Solve the trust region subproblem */
268       ierr = KSPCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
269       ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
270       ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
271       tao->ksp_its+=its;
272       tao->ksp_tot_its+=its;
273       ierr = KSPCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
274 
275       if (0.0 == tao->trust) {
276         /* Radius was uninitialized; use the norm of the direction */
277         if (norm_d > 0.0) {
278           tao->trust = norm_d;
279 
280           /* Modify the radius if it is too large or small */
281           tao->trust = PetscMax(tao->trust, tr->min_radius);
282           tao->trust = PetscMin(tao->trust, tr->max_radius);
283         }
284         else {
285           /* The direction was bad; set radius to default value and re-solve
286              the trust-region subproblem to get a direction */
287           tao->trust = tao->trust0;
288 
289           /* Modify the radius if it is too large or small */
290           tao->trust = PetscMax(tao->trust, tr->min_radius);
291           tao->trust = PetscMin(tao->trust, tr->max_radius);
292 
293           ierr = KSPCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr);
294           ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr);
295           ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr);
296           tao->ksp_its+=its;
297           tao->ksp_tot_its+=its;
298           ierr = KSPCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr);
299 
300           if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
301         }
302       }
303       ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr);
304       ierr = KSPGetConvergedReason(tao->ksp, &ksp_reason);CHKERRQ(ierr);
305       if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
306         /* Preconditioner is numerically indefinite; reset the
307            approximate if using BFGS preconditioning. */
308         ierr = MatLMVMReset(tr->M, PETSC_FALSE);CHKERRQ(ierr);
309         ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr);
310         bfgsUpdates = 1;
311       }
312 
313       if (NTR_UPDATE_REDUCTION == tr->update_type) {
314         /* Get predicted reduction */
315         ierr = KSPCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
316         if (prered >= 0.0) {
317           /* The predicted reduction has the wrong sign.  This cannot
318              happen in infinite precision arithmetic.  Step should
319              be rejected! */
320           tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
321         }
322         else {
323           /* Compute trial step and function value */
324           ierr = VecCopy(tao->solution,tr->W);CHKERRQ(ierr);
325           ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr);
326           ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
327 
328           if (PetscIsInfOrNanReal(ftrial)) {
329             tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
330           } else {
331             /* Compute and actual reduction */
332             actred = f - ftrial;
333             prered = -prered;
334             if ((PetscAbsScalar(actred) <= tr->epsilon) &&
335                 (PetscAbsScalar(prered) <= tr->epsilon)) {
336               kappa = 1.0;
337             }
338             else {
339               kappa = actred / prered;
340             }
341 
342             /* Accept or reject the step and update radius */
343             if (kappa < tr->eta1) {
344               /* Reject the step */
345               tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
346             }
347             else {
348               /* Accept the step */
349               if (kappa < tr->eta2) {
350                 /* Marginal bad step */
351                 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
352               }
353               else if (kappa < tr->eta3) {
354                 /* Reasonable step */
355                 tao->trust = tr->alpha3 * tao->trust;
356               }
357               else if (kappa < tr->eta4) {
358                 /* Good step */
359                 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
360               }
361               else {
362                 /* Very good step */
363                 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
364               }
365               break;
366             }
367           }
368         }
369       }
370       else {
371         /* Get predicted reduction */
372         ierr = KSPCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr);
373         if (prered >= 0.0) {
374           /* The predicted reduction has the wrong sign.  This cannot
375              happen in infinite precision arithmetic.  Step should
376              be rejected! */
377           tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
378         }
379         else {
380           ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr);
381           ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr);
382           ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr);
383           if (PetscIsInfOrNanReal(ftrial)) {
384             tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
385           }
386           else {
387             ierr = VecDot(tao->gradient, tao->stepdirection, &beta);CHKERRQ(ierr);
388             actred = f - ftrial;
389             prered = -prered;
390             if ((PetscAbsScalar(actred) <= tr->epsilon) &&
391                 (PetscAbsScalar(prered) <= tr->epsilon)) {
392               kappa = 1.0;
393             }
394             else {
395               kappa = actred / prered;
396             }
397 
398             tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
399             tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
400             tau_min = PetscMin(tau_1, tau_2);
401             tau_max = PetscMax(tau_1, tau_2);
402 
403             if (kappa >= 1.0 - tr->mu1) {
404               /* Great agreement; accept step and update radius */
405               if (tau_max < 1.0) {
406                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
407               }
408               else if (tau_max > tr->gamma4) {
409                 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
410               }
411               else {
412                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
413               }
414               break;
415             }
416             else if (kappa >= 1.0 - tr->mu2) {
417               /* Good agreement */
418 
419               if (tau_max < tr->gamma2) {
420                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
421               }
422               else if (tau_max > tr->gamma3) {
423                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
424               }
425               else if (tau_max < 1.0) {
426                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
427               }
428               else {
429                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
430               }
431               break;
432             }
433             else {
434               /* Not good agreement */
435               if (tau_min > 1.0) {
436                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
437               }
438               else if (tau_max < tr->gamma1) {
439                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
440               }
441               else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
442                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
443               }
444               else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) &&
445                        ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
446                 tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
447               }
448               else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) &&
449                        ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
450                 tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
451               }
452               else {
453                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
454               }
455             }
456           }
457         }
458       }
459 
460       /* The step computed was not good and the radius was decreased.
461          Monitor the radius to terminate. */
462       ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
463       ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);CHKERRQ(ierr);
464       ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
465     }
466 
467     /* The radius may have been increased; modify if it is too large */
468     tao->trust = PetscMin(tao->trust, tr->max_radius);
469 
470     if (tao->reason == TAO_CONTINUE_ITERATING) {
471       ierr = VecCopy(tr->W, tao->solution);CHKERRQ(ierr);
472       f = ftrial;
473       ierr = TaoComputeGradient(tao, tao->solution, tao->gradient);CHKERRQ(ierr);
474       ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
475       if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
476       needH = 1;
477       ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
478       ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust);CHKERRQ(ierr);
479       ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
480     }
481   }
482   PetscFunctionReturn(0);
483 }
484 
485 /*------------------------------------------------------------*/
486 static PetscErrorCode TaoSetUp_NTR(Tao tao)
487 {
488   TAO_NTR *tr = (TAO_NTR *)tao->data;
489   PetscErrorCode ierr;
490 
491   PetscFunctionBegin;
492   if (!tao->gradient) {ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr);}
493   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);}
494   if (!tr->W) {ierr = VecDuplicate(tao->solution, &tr->W);CHKERRQ(ierr);}
495 
496   tr->bfgs_pre = 0;
497   tr->M = 0;
498   PetscFunctionReturn(0);
499 }
500 
501 /*------------------------------------------------------------*/
502 static PetscErrorCode TaoDestroy_NTR(Tao tao)
503 {
504   TAO_NTR        *tr = (TAO_NTR *)tao->data;
505   PetscErrorCode ierr;
506 
507   PetscFunctionBegin;
508   if (tao->setupcalled) {
509     ierr = VecDestroy(&tr->W);CHKERRQ(ierr);
510   }
511   ierr = PetscFree(tao->data);CHKERRQ(ierr);
512   PetscFunctionReturn(0);
513 }
514 
515 /*------------------------------------------------------------*/
516 static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao)
517 {
518   TAO_NTR        *tr = (TAO_NTR *)tao->data;
519   PetscErrorCode ierr;
520 
521   PetscFunctionBegin;
522   ierr = PetscOptionsHead(PetscOptionsObject,"Newton trust region method for unconstrained optimization");CHKERRQ(ierr);
523   ierr = PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL);CHKERRQ(ierr);
524   ierr = PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL);CHKERRQ(ierr);
525   ierr = PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL);CHKERRQ(ierr);
526   ierr = PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL);CHKERRQ(ierr);
527   ierr = PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL);CHKERRQ(ierr);
528   ierr = PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL);CHKERRQ(ierr);
529   ierr = PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL);CHKERRQ(ierr);
530   ierr = PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL);CHKERRQ(ierr);
531   ierr = PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL);CHKERRQ(ierr);
532   ierr = PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL);CHKERRQ(ierr);
533   ierr = PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL);CHKERRQ(ierr);
534   ierr = PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL);CHKERRQ(ierr);
535   ierr = PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL);CHKERRQ(ierr);
536   ierr = PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL);CHKERRQ(ierr);
537   ierr = PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL);CHKERRQ(ierr);
538   ierr = PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL);CHKERRQ(ierr);
539   ierr = PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL);CHKERRQ(ierr);
540   ierr = PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL);CHKERRQ(ierr);
541   ierr = PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL);CHKERRQ(ierr);
542   ierr = PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL);CHKERRQ(ierr);
543   ierr = PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL);CHKERRQ(ierr);
544   ierr = PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL);CHKERRQ(ierr);
545   ierr = PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL);CHKERRQ(ierr);
546   ierr = PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL);CHKERRQ(ierr);
547   ierr = PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL);CHKERRQ(ierr);
548   ierr = PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL);CHKERRQ(ierr);
549   ierr = PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL);CHKERRQ(ierr);
550   ierr = PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL);CHKERRQ(ierr);
551   ierr = PetscOptionsTail();CHKERRQ(ierr);
552   ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr);
553   PetscFunctionReturn(0);
554 }
555 
556 /*------------------------------------------------------------*/
557 /*MC
558   TAONTR - Newton's method with trust region for unconstrained minimization.
559   At each iteration, the Newton trust region method solves the system.
560   NTR expects a KSP solver with a trust region radius.
561             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k
562 
563   Options Database Keys:
564 + -tao_ntr_init_type - "constant","direction","interpolation"
565 . -tao_ntr_update_type - "reduction","interpolation"
566 . -tao_ntr_min_radius - lower bound on trust region radius
567 . -tao_ntr_max_radius - upper bound on trust region radius
568 . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
569 . -tao_ntr_mu1_i - mu1 interpolation init factor
570 . -tao_ntr_mu2_i - mu2 interpolation init factor
571 . -tao_ntr_gamma1_i - gamma1 interpolation init factor
572 . -tao_ntr_gamma2_i - gamma2 interpolation init factor
573 . -tao_ntr_gamma3_i - gamma3 interpolation init factor
574 . -tao_ntr_gamma4_i - gamma4 interpolation init factor
575 . -tao_ntr_theta_i - thetha1 interpolation init factor
576 . -tao_ntr_eta1 - eta1 reduction update factor
577 . -tao_ntr_eta2 - eta2 reduction update factor
578 . -tao_ntr_eta3 - eta3 reduction update factor
579 . -tao_ntr_eta4 - eta4 reduction update factor
580 . -tao_ntr_alpha1 - alpha1 reduction update factor
581 . -tao_ntr_alpha2 - alpha2 reduction update factor
582 . -tao_ntr_alpha3 - alpha3 reduction update factor
583 . -tao_ntr_alpha4 - alpha4 reduction update factor
584 . -tao_ntr_alpha4 - alpha4 reduction update factor
585 . -tao_ntr_mu1 - mu1 interpolation update
586 . -tao_ntr_mu2 - mu2 interpolation update
587 . -tao_ntr_gamma1 - gamma1 interpolcation update
588 . -tao_ntr_gamma2 - gamma2 interpolcation update
589 . -tao_ntr_gamma3 - gamma3 interpolcation update
590 . -tao_ntr_gamma4 - gamma4 interpolation update
591 - -tao_ntr_theta - theta interpolation update
592 
593   Level: beginner
594 M*/
595 
596 PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
597 {
598   TAO_NTR *tr;
599   PetscErrorCode ierr;
600 
601   PetscFunctionBegin;
602 
603   ierr = PetscNewLog(tao,&tr);CHKERRQ(ierr);
604 
605   tao->ops->setup = TaoSetUp_NTR;
606   tao->ops->solve = TaoSolve_NTR;
607   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
608   tao->ops->destroy = TaoDestroy_NTR;
609 
610   /* Override default settings (unless already changed) */
611   if (!tao->max_it_changed) tao->max_it = 50;
612   if (!tao->trust0_changed) tao->trust0 = 100.0;
613   tao->data = (void*)tr;
614 
615   /*  Standard trust region update parameters */
616   tr->eta1 = 1.0e-4;
617   tr->eta2 = 0.25;
618   tr->eta3 = 0.50;
619   tr->eta4 = 0.90;
620 
621   tr->alpha1 = 0.25;
622   tr->alpha2 = 0.50;
623   tr->alpha3 = 1.00;
624   tr->alpha4 = 2.00;
625   tr->alpha5 = 4.00;
626 
627   /*  Interpolation trust region update parameters */
628   tr->mu1 = 0.10;
629   tr->mu2 = 0.50;
630 
631   tr->gamma1 = 0.25;
632   tr->gamma2 = 0.50;
633   tr->gamma3 = 2.00;
634   tr->gamma4 = 4.00;
635 
636   tr->theta = 0.05;
637 
638   /*  Interpolation parameters for initialization */
639   tr->mu1_i = 0.35;
640   tr->mu2_i = 0.50;
641 
642   tr->gamma1_i = 0.0625;
643   tr->gamma2_i = 0.50;
644   tr->gamma3_i = 2.00;
645   tr->gamma4_i = 5.00;
646 
647   tr->theta_i = 0.25;
648 
649   tr->min_radius = 1.0e-10;
650   tr->max_radius = 1.0e10;
651   tr->epsilon    = 1.0e-6;
652 
653   tr->init_type       = NTR_INIT_INTERPOLATION;
654   tr->update_type     = NTR_UPDATE_REDUCTION;
655 
656   /* Set linear solver to default for trust region */
657   ierr = KSPCreate(((PetscObject)tao)->comm,&tao->ksp);CHKERRQ(ierr);
658   ierr = PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1);CHKERRQ(ierr);
659   ierr = KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);CHKERRQ(ierr);
660   ierr = KSPAppendOptionsPrefix(tao->ksp,"tao_ntr_");CHKERRQ(ierr);
661   ierr = KSPSetType(tao->ksp,KSPCGSTCG);CHKERRQ(ierr);
662   PetscFunctionReturn(0);
663 }
664