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