Actual source code: ntl.c
1: #include <../src/tao/unconstrained/impls/ntl/ntlimpl.h>
3: #include <petscksp.h>
5: #define NTL_INIT_CONSTANT 0
6: #define NTL_INIT_DIRECTION 1
7: #define NTL_INIT_INTERPOLATION 2
8: #define NTL_INIT_TYPES 3
10: #define NTL_UPDATE_REDUCTION 0
11: #define NTL_UPDATE_INTERPOLATION 1
12: #define NTL_UPDATE_TYPES 2
14: static const char *NTL_INIT[64] = {"constant", "direction", "interpolation"};
16: static const char *NTL_UPDATE[64] = {"reduction", "interpolation"};
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. */
23: #define NTL_NEWTON 0
24: #define NTL_BFGS 1
25: #define NTL_SCALED_GRADIENT 2
26: #define NTL_GRADIENT 3
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;
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;
42: PetscReal norm_d = 0.0;
43: PetscInt stepType;
44: PetscInt its;
46: PetscInt bfgsUpdates = 0;
47: PetscInt needH;
49: PetscInt i_max = 5;
50: PetscInt j_max = 1;
51: PetscInt i, j, n, N;
53: PetscInt tr_reject;
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"));
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");
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);
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));
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 Inf or NaN");
88: needH = 1;
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);
96: /* Initialize trust-region radius */
97: switch (tl->init_type) {
98: case NTL_INIT_CONSTANT:
99: /* Use the initial radius specified */
100: break;
102: case NTL_INIT_INTERPOLATION:
103: /* Use the initial radius specified */
104: max_radius = 0.0;
106: for (j = 0; j < j_max; ++j) {
107: fmin = f;
108: sigma = 0.0;
110: if (needH) {
111: PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
112: needH = 0;
113: }
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));
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: }
128: PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection));
129: PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered));
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: }
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);
144: if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu1_i) {
145: /* Great agreement */
146: max_radius = PetscMax(max_radius, tao->trust);
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);
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: }
190: if (fmin < f) {
191: f = fmin;
192: PetscCall(VecAXPY(tao->solution, sigma, tao->gradient));
193: PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
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 Inf or NaN");
197: needH = 1;
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);
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;
212: default:
213: /* Norm of the first direction will initialize radius */
214: tao->trust = 0.0;
215: break;
216: }
218: /* Set counter for gradient/reset steps */
219: tl->ntrust = 0;
220: tl->newt = 0;
221: tl->bfgs = 0;
222: tl->grad = 0;
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));
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));
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;
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;
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);
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));
274: PetscCheck(norm_d != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero");
275: }
276: }
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: }
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));
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: }
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));
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: }
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);
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 */
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: }
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 Inf or NaN */
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));
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 */
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));
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: }
469: /* Perform the linesearch */
470: fold = f;
471: PetscCall(VecCopy(tao->solution, tl->Xold));
472: PetscCall(VecCopy(tao->gradient, tl->Gold));
474: step = 1.0;
475: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
476: PetscCall(TaoAddLineSearchCounts(tao));
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));
484: switch (stepType) {
485: case NTL_NEWTON:
486: /* Failed to obtain acceptable iterate with Newton step */
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));
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));
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;
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));
533: bfgsUpdates = 1;
534: ++tl->grad;
535: stepType = NTL_GRADIENT;
536: break;
537: }
538: PetscCall(VecScale(tao->stepdirection, -1.0));
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: }
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: }
589: /* The radius may have been increased; modify if it is too large */
590: tao->trust = PetscMin(tao->trust, tl->max_radius);
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;
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: }
604: /* ---------------------------------------------------------- */
605: static PetscErrorCode TaoSetUp_NTL(Tao tao)
606: {
607: TAO_NTL *tl = (TAO_NTL *)tao->data;
609: PetscFunctionBegin;
610: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
611: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
612: if (!tl->W) PetscCall(VecDuplicate(tao->solution, &tl->W));
613: if (!tl->Xold) PetscCall(VecDuplicate(tao->solution, &tl->Xold));
614: if (!tl->Gold) PetscCall(VecDuplicate(tao->solution, &tl->Gold));
615: tl->bfgs_pre = NULL;
616: tl->M = NULL;
617: PetscFunctionReturn(PETSC_SUCCESS);
618: }
620: /*------------------------------------------------------------*/
621: static PetscErrorCode TaoDestroy_NTL(Tao tao)
622: {
623: TAO_NTL *tl = (TAO_NTL *)tao->data;
625: PetscFunctionBegin;
626: if (tao->setupcalled) {
627: PetscCall(VecDestroy(&tl->W));
628: PetscCall(VecDestroy(&tl->Xold));
629: PetscCall(VecDestroy(&tl->Gold));
630: }
631: PetscCall(KSPDestroy(&tao->ksp));
632: PetscCall(PetscFree(tao->data));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: /*------------------------------------------------------------*/
637: static PetscErrorCode TaoSetFromOptions_NTL(Tao tao, PetscOptionItems *PetscOptionsObject)
638: {
639: TAO_NTL *tl = (TAO_NTL *)tao->data;
641: PetscFunctionBegin;
642: PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region with line search method for unconstrained optimization");
643: PetscCall(PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type, NULL));
644: PetscCall(PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type, NULL));
645: PetscCall(PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1, NULL));
646: PetscCall(PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2, NULL));
647: PetscCall(PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3, NULL));
648: PetscCall(PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4, NULL));
649: PetscCall(PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1, NULL));
650: PetscCall(PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2, NULL));
651: PetscCall(PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3, NULL));
652: PetscCall(PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4, NULL));
653: PetscCall(PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5, NULL));
654: PetscCall(PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1, NULL));
655: PetscCall(PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2, NULL));
656: PetscCall(PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3, NULL));
657: PetscCall(PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4, NULL));
658: PetscCall(PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1, NULL));
659: PetscCall(PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2, NULL));
660: PetscCall(PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3, NULL));
661: PetscCall(PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4, NULL));
662: PetscCall(PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5, NULL));
663: PetscCall(PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i, NULL));
664: PetscCall(PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i, NULL));
665: PetscCall(PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i, NULL));
666: PetscCall(PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i, NULL));
667: PetscCall(PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i, NULL));
668: PetscCall(PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i, NULL));
669: PetscCall(PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i, NULL));
670: PetscCall(PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1, NULL));
671: PetscCall(PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2, NULL));
672: PetscCall(PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1, NULL));
673: PetscCall(PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2, NULL));
674: PetscCall(PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3, NULL));
675: PetscCall(PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4, NULL));
676: PetscCall(PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta, NULL));
677: PetscCall(PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius, NULL));
678: PetscCall(PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius, NULL));
679: PetscCall(PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon, NULL));
680: PetscOptionsHeadEnd();
681: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
682: PetscCall(KSPSetFromOptions(tao->ksp));
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: /*------------------------------------------------------------*/
687: static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer)
688: {
689: TAO_NTL *tl = (TAO_NTL *)tao->data;
690: PetscBool isascii;
692: PetscFunctionBegin;
693: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
694: if (isascii) {
695: PetscCall(PetscViewerASCIIPushTab(viewer));
696: PetscCall(PetscViewerASCIIPrintf(viewer, "Trust-region steps: %" PetscInt_FMT "\n", tl->ntrust));
697: PetscCall(PetscViewerASCIIPrintf(viewer, "Newton search steps: %" PetscInt_FMT "\n", tl->newt));
698: PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS search steps: %" PetscInt_FMT "\n", tl->bfgs));
699: PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient search steps: %" PetscInt_FMT "\n", tl->grad));
700: PetscCall(PetscViewerASCIIPopTab(viewer));
701: }
702: PetscFunctionReturn(PETSC_SUCCESS);
703: }
705: /* ---------------------------------------------------------- */
706: /*MC
707: TAONTL - Newton's method with trust region globalization and line search fallback.
708: At each iteration, the Newton trust region method solves the system for d
709: and performs a line search in the d direction:
711: min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k
713: Options Database Keys:
714: + -tao_ntl_init_type - "constant","direction","interpolation"
715: . -tao_ntl_update_type - "reduction","interpolation"
716: . -tao_ntl_min_radius - lower bound on trust region radius
717: . -tao_ntl_max_radius - upper bound on trust region radius
718: . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
719: . -tao_ntl_mu1_i - mu1 interpolation init factor
720: . -tao_ntl_mu2_i - mu2 interpolation init factor
721: . -tao_ntl_gamma1_i - gamma1 interpolation init factor
722: . -tao_ntl_gamma2_i - gamma2 interpolation init factor
723: . -tao_ntl_gamma3_i - gamma3 interpolation init factor
724: . -tao_ntl_gamma4_i - gamma4 interpolation init factor
725: . -tao_ntl_theta_i - theta1 interpolation init factor
726: . -tao_ntl_eta1 - eta1 reduction update factor
727: . -tao_ntl_eta2 - eta2 reduction update factor
728: . -tao_ntl_eta3 - eta3 reduction update factor
729: . -tao_ntl_eta4 - eta4 reduction update factor
730: . -tao_ntl_alpha1 - alpha1 reduction update factor
731: . -tao_ntl_alpha2 - alpha2 reduction update factor
732: . -tao_ntl_alpha3 - alpha3 reduction update factor
733: . -tao_ntl_alpha4 - alpha4 reduction update factor
734: . -tao_ntl_alpha4 - alpha4 reduction update factor
735: . -tao_ntl_mu1 - mu1 interpolation update
736: . -tao_ntl_mu2 - mu2 interpolation update
737: . -tao_ntl_gamma1 - gamma1 interpolcation update
738: . -tao_ntl_gamma2 - gamma2 interpolcation update
739: . -tao_ntl_gamma3 - gamma3 interpolcation update
740: . -tao_ntl_gamma4 - gamma4 interpolation update
741: - -tao_ntl_theta - theta1 interpolation update
743: Level: beginner
744: M*/
745: PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao)
746: {
747: TAO_NTL *tl;
748: const char *morethuente_type = TAOLINESEARCHMT;
750: PetscFunctionBegin;
751: PetscCall(PetscNew(&tl));
752: tao->ops->setup = TaoSetUp_NTL;
753: tao->ops->solve = TaoSolve_NTL;
754: tao->ops->view = TaoView_NTL;
755: tao->ops->setfromoptions = TaoSetFromOptions_NTL;
756: tao->ops->destroy = TaoDestroy_NTL;
758: /* Override default settings (unless already changed) */
759: PetscCall(TaoParametersInitialize(tao));
760: PetscObjectParameterSetDefault(tao, max_it, 50);
761: PetscObjectParameterSetDefault(tao, trust0, 100.0);
763: tao->data = (void *)tl;
765: /* Default values for trust-region radius update based on steplength */
766: tl->nu1 = 0.25;
767: tl->nu2 = 0.50;
768: tl->nu3 = 1.00;
769: tl->nu4 = 1.25;
771: tl->omega1 = 0.25;
772: tl->omega2 = 0.50;
773: tl->omega3 = 1.00;
774: tl->omega4 = 2.00;
775: tl->omega5 = 4.00;
777: /* Default values for trust-region radius update based on reduction */
778: tl->eta1 = 1.0e-4;
779: tl->eta2 = 0.25;
780: tl->eta3 = 0.50;
781: tl->eta4 = 0.90;
783: tl->alpha1 = 0.25;
784: tl->alpha2 = 0.50;
785: tl->alpha3 = 1.00;
786: tl->alpha4 = 2.00;
787: tl->alpha5 = 4.00;
789: /* Default values for trust-region radius update based on interpolation */
790: tl->mu1 = 0.10;
791: tl->mu2 = 0.50;
793: tl->gamma1 = 0.25;
794: tl->gamma2 = 0.50;
795: tl->gamma3 = 2.00;
796: tl->gamma4 = 4.00;
798: tl->theta = 0.05;
800: /* Default values for trust region initialization based on interpolation */
801: tl->mu1_i = 0.35;
802: tl->mu2_i = 0.50;
804: tl->gamma1_i = 0.0625;
805: tl->gamma2_i = 0.5;
806: tl->gamma3_i = 2.0;
807: tl->gamma4_i = 5.0;
809: tl->theta_i = 0.25;
811: /* Remaining parameters */
812: tl->min_radius = 1.0e-10;
813: tl->max_radius = 1.0e10;
814: tl->epsilon = 1.0e-6;
816: tl->init_type = NTL_INIT_INTERPOLATION;
817: tl->update_type = NTL_UPDATE_REDUCTION;
819: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
820: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
821: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
822: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
823: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
824: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
825: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
826: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
827: PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_ntl_"));
828: PetscCall(KSPSetType(tao->ksp, KSPSTCG));
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }