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: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
228: ++tao->niter;
229: tao->ksp_its = 0;
230: /* Compute the Hessian */
231: if (needH) PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
233: if (tl->bfgs_pre) {
234: /* Update the limited memory preconditioner */
235: PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
236: ++bfgsUpdates;
237: }
238: PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));
239: /* Solve the Newton system of equations */
240: PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
241: PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
242: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
243: tao->ksp_its += its;
244: tao->ksp_tot_its += its;
245: PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));
247: if (0.0 == tao->trust) {
248: /* Radius was uninitialized; use the norm of the direction */
249: if (norm_d > 0.0) {
250: tao->trust = norm_d;
252: /* Modify the radius if it is too large or small */
253: tao->trust = PetscMax(tao->trust, tl->min_radius);
254: tao->trust = PetscMin(tao->trust, tl->max_radius);
255: } else {
256: /* The direction was bad; set radius to default value and re-solve
257: the trust-region subproblem to get a direction */
258: tao->trust = tao->trust0;
260: /* Modify the radius if it is too large or small */
261: tao->trust = PetscMax(tao->trust, tl->min_radius);
262: tao->trust = PetscMin(tao->trust, tl->max_radius);
264: PetscCall(KSPCGSetRadius(tao->ksp, tl->max_radius));
265: PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
266: PetscCall(KSPGetIterationNumber(tao->ksp, &its));
267: tao->ksp_its += its;
268: tao->ksp_tot_its += its;
269: PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));
271: PetscCheck(norm_d != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero");
272: }
273: }
275: PetscCall(VecScale(tao->stepdirection, -1.0));
276: PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason));
277: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tl->bfgs_pre)) {
278: /* Preconditioner is numerically indefinite; reset the
279: approximate if using BFGS preconditioning. */
280: PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
281: PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
282: bfgsUpdates = 1;
283: }
285: /* Check trust-region reduction conditions */
286: tr_reject = 0;
287: if (NTL_UPDATE_REDUCTION == tl->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 = tl->alpha1 * PetscMin(tao->trust, norm_d);
295: tr_reject = 1;
296: } else {
297: /* Compute trial step and function value */
298: PetscCall(VecCopy(tao->solution, tl->W));
299: PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
300: PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
302: if (PetscIsInfOrNanReal(ftrial)) {
303: tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
304: tr_reject = 1;
305: } else {
306: /* Compute and actual reduction */
307: actred = f - ftrial;
308: prered = -prered;
309: if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
310: kappa = 1.0;
311: } else {
312: kappa = actred / prered;
313: }
315: /* Accept of reject the step and update radius */
316: if (kappa < tl->eta1) {
317: /* Reject the step */
318: tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
319: tr_reject = 1;
320: } else {
321: /* Accept the step */
322: if (kappa < tl->eta2) {
323: /* Marginal bad step */
324: tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
325: } else if (kappa < tl->eta3) {
326: /* Reasonable step */
327: tao->trust = tl->alpha3 * tao->trust;
328: } else if (kappa < tl->eta4) {
329: /* Good step */
330: tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
331: } else {
332: /* Very good step */
333: tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
334: }
335: }
336: }
337: }
338: } else {
339: /* Get predicted reduction */
340: PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
341: if (prered >= 0.0) {
342: /* The predicted reduction has the wrong sign. This cannot
343: happen in infinite precision arithmetic. Step should
344: be rejected! */
345: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
346: tr_reject = 1;
347: } else {
348: PetscCall(VecCopy(tao->solution, tl->W));
349: PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection));
350: PetscCall(TaoComputeObjective(tao, tl->W, &ftrial));
351: if (PetscIsInfOrNanReal(ftrial)) {
352: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
353: tr_reject = 1;
354: } else {
355: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx));
357: actred = f - ftrial;
358: prered = -prered;
359: if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
360: kappa = 1.0;
361: } else {
362: kappa = actred / prered;
363: }
365: tau_1 = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
366: tau_2 = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
367: tau_min = PetscMin(tau_1, tau_2);
368: tau_max = PetscMax(tau_1, tau_2);
370: if (kappa >= 1.0 - tl->mu1) {
371: /* Great agreement; accept step and update radius */
372: if (tau_max < 1.0) {
373: tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
374: } else if (tau_max > tl->gamma4) {
375: tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
376: } else {
377: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
378: }
379: } else if (kappa >= 1.0 - tl->mu2) {
380: /* Good agreement */
382: if (tau_max < tl->gamma2) {
383: tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
384: } else if (tau_max > tl->gamma3) {
385: tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
386: } else if (tau_max < 1.0) {
387: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
388: } else {
389: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
390: }
391: } else {
392: /* Not good agreement */
393: if (tau_min > 1.0) {
394: tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
395: } else if (tau_max < tl->gamma1) {
396: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
397: } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
398: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
399: } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
400: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
401: } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
402: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
403: } else {
404: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
405: }
406: tr_reject = 1;
407: }
408: }
409: }
410: }
412: if (tr_reject) {
413: /* The trust-region constraints rejected the step. Apply a linesearch.
414: Check for descent direction. */
415: PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
416: if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
417: /* Newton step is not descent or direction produced Inf or NaN */
419: if (!tl->bfgs_pre) {
420: /* We don't have the bfgs matrix around and updated
421: Must use gradient direction in this case */
422: PetscCall(VecCopy(tao->gradient, tao->stepdirection));
423: PetscCall(VecScale(tao->stepdirection, -1.0));
424: ++tl->grad;
425: stepType = NTL_GRADIENT;
426: } else {
427: /* Attempt to use the BFGS direction */
428: PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
429: PetscCall(VecScale(tao->stepdirection, -1.0));
431: /* Check for success (descent direction) */
432: PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
433: if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
434: /* BFGS direction is not descent or direction produced not a number
435: We can assert bfgsUpdates > 1 in this case because
436: the first solve produces the scaled gradient direction,
437: which is guaranteed to be descent */
439: /* Use steepest descent direction (scaled) */
440: PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
441: PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
442: PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
443: PetscCall(VecScale(tao->stepdirection, -1.0));
445: bfgsUpdates = 1;
446: ++tl->grad;
447: stepType = NTL_GRADIENT;
448: } else {
449: PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
450: if (1 == bfgsUpdates) {
451: /* The first BFGS direction is always the scaled gradient */
452: ++tl->grad;
453: stepType = NTL_GRADIENT;
454: } else {
455: ++tl->bfgs;
456: stepType = NTL_BFGS;
457: }
458: }
459: }
460: } else {
461: /* Computed Newton step is descent */
462: ++tl->newt;
463: stepType = NTL_NEWTON;
464: }
466: /* Perform the linesearch */
467: fold = f;
468: PetscCall(VecCopy(tao->solution, tl->Xold));
469: PetscCall(VecCopy(tao->gradient, tl->Gold));
471: step = 1.0;
472: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
473: PetscCall(TaoAddLineSearchCounts(tao));
475: while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */
476: /* Linesearch failed */
477: f = fold;
478: PetscCall(VecCopy(tl->Xold, tao->solution));
479: PetscCall(VecCopy(tl->Gold, tao->gradient));
481: switch (stepType) {
482: case NTL_NEWTON:
483: /* Failed to obtain acceptable iterate with Newton step */
485: if (tl->bfgs_pre) {
486: /* We don't have the bfgs matrix around and being updated
487: Must use gradient direction in this case */
488: PetscCall(VecCopy(tao->gradient, tao->stepdirection));
489: ++tl->grad;
490: stepType = NTL_GRADIENT;
491: } else {
492: /* Attempt to use the BFGS direction */
493: PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
495: /* Check for success (descent direction) */
496: PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
497: if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
498: /* BFGS direction is not descent or direction produced
499: not a number. We can assert bfgsUpdates > 1 in this case
500: Use steepest descent direction (scaled) */
501: PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
502: PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
503: PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
505: bfgsUpdates = 1;
506: ++tl->grad;
507: stepType = NTL_GRADIENT;
508: } else {
509: PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates));
510: if (1 == bfgsUpdates) {
511: /* The first BFGS direction is always the scaled gradient */
512: ++tl->grad;
513: stepType = NTL_GRADIENT;
514: } else {
515: ++tl->bfgs;
516: stepType = NTL_BFGS;
517: }
518: }
519: }
520: break;
522: case NTL_BFGS:
523: /* Can only enter if pc_type == NTL_PC_BFGS
524: Failed to obtain acceptable iterate with BFGS step
525: Attempt to use the scaled gradient direction */
526: PetscCall(MatLMVMReset(tl->M, PETSC_FALSE));
527: PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient));
528: PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection));
530: bfgsUpdates = 1;
531: ++tl->grad;
532: stepType = NTL_GRADIENT;
533: break;
534: }
535: PetscCall(VecScale(tao->stepdirection, -1.0));
537: /* This may be incorrect; linesearch has values for stepmax and stepmin
538: that should be reset. */
539: step = 1.0;
540: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason));
541: PetscCall(TaoAddLineSearchCounts(tao));
542: }
544: if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
545: /* Failed to find an improving point */
546: f = fold;
547: PetscCall(VecCopy(tl->Xold, tao->solution));
548: PetscCall(VecCopy(tl->Gold, tao->gradient));
549: tao->trust = 0.0;
550: step = 0.0;
551: tao->reason = TAO_DIVERGED_LS_FAILURE;
552: break;
553: } else if (stepType == NTL_NEWTON) {
554: if (step < tl->nu1) {
555: /* Very bad step taken; reduce radius */
556: tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
557: } else if (step < tl->nu2) {
558: /* Reasonably bad step taken; reduce radius */
559: tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
560: } else if (step < tl->nu3) {
561: /* Reasonable step was taken; leave radius alone */
562: if (tl->omega3 < 1.0) {
563: tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
564: } else if (tl->omega3 > 1.0) {
565: tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
566: }
567: } else if (step < tl->nu4) {
568: /* Full step taken; increase the radius */
569: tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
570: } else {
571: /* More than full step taken; increase the radius */
572: tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
573: }
574: } else {
575: /* Newton step was not good; reduce the radius */
576: tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
577: }
578: } else {
579: /* Trust-region step is accepted */
580: PetscCall(VecCopy(tl->W, tao->solution));
581: f = ftrial;
582: PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
583: ++tl->ntrust;
584: }
586: /* The radius may have been increased; modify if it is too large */
587: tao->trust = PetscMin(tao->trust, tl->max_radius);
589: /* Check for converged */
590: PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
591: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
592: needH = 1;
594: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
595: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
596: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
597: }
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: /* ---------------------------------------------------------- */
602: static PetscErrorCode TaoSetUp_NTL(Tao tao)
603: {
604: TAO_NTL *tl = (TAO_NTL *)tao->data;
606: PetscFunctionBegin;
607: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
608: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
609: if (!tl->W) PetscCall(VecDuplicate(tao->solution, &tl->W));
610: if (!tl->Xold) PetscCall(VecDuplicate(tao->solution, &tl->Xold));
611: if (!tl->Gold) PetscCall(VecDuplicate(tao->solution, &tl->Gold));
612: tl->bfgs_pre = NULL;
613: tl->M = NULL;
614: PetscFunctionReturn(PETSC_SUCCESS);
615: }
617: /*------------------------------------------------------------*/
618: static PetscErrorCode TaoDestroy_NTL(Tao tao)
619: {
620: TAO_NTL *tl = (TAO_NTL *)tao->data;
622: PetscFunctionBegin;
623: if (tao->setupcalled) {
624: PetscCall(VecDestroy(&tl->W));
625: PetscCall(VecDestroy(&tl->Xold));
626: PetscCall(VecDestroy(&tl->Gold));
627: }
628: PetscCall(KSPDestroy(&tao->ksp));
629: PetscCall(PetscFree(tao->data));
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: /*------------------------------------------------------------*/
634: static PetscErrorCode TaoSetFromOptions_NTL(Tao tao, PetscOptionItems *PetscOptionsObject)
635: {
636: TAO_NTL *tl = (TAO_NTL *)tao->data;
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: }
683: /*------------------------------------------------------------*/
684: static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer)
685: {
686: TAO_NTL *tl = (TAO_NTL *)tao->data;
687: PetscBool isascii;
689: PetscFunctionBegin;
690: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
691: if (isascii) {
692: PetscCall(PetscViewerASCIIPushTab(viewer));
693: PetscCall(PetscViewerASCIIPrintf(viewer, "Trust-region steps: %" PetscInt_FMT "\n", tl->ntrust));
694: PetscCall(PetscViewerASCIIPrintf(viewer, "Newton search steps: %" PetscInt_FMT "\n", tl->newt));
695: PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS search steps: %" PetscInt_FMT "\n", tl->bfgs));
696: PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient search steps: %" PetscInt_FMT "\n", tl->grad));
697: PetscCall(PetscViewerASCIIPopTab(viewer));
698: }
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: /* ---------------------------------------------------------- */
703: /*MC
704: TAONTL - Newton's method with trust region globalization and line search fallback.
705: At each iteration, the Newton trust region method solves the system for d
706: and performs a line search in the d direction:
708: min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k
710: Options Database Keys:
711: + -tao_ntl_init_type - "constant","direction","interpolation"
712: . -tao_ntl_update_type - "reduction","interpolation"
713: . -tao_ntl_min_radius - lower bound on trust region radius
714: . -tao_ntl_max_radius - upper bound on trust region radius
715: . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
716: . -tao_ntl_mu1_i - mu1 interpolation init factor
717: . -tao_ntl_mu2_i - mu2 interpolation init factor
718: . -tao_ntl_gamma1_i - gamma1 interpolation init factor
719: . -tao_ntl_gamma2_i - gamma2 interpolation init factor
720: . -tao_ntl_gamma3_i - gamma3 interpolation init factor
721: . -tao_ntl_gamma4_i - gamma4 interpolation init factor
722: . -tao_ntl_theta_i - theta1 interpolation init factor
723: . -tao_ntl_eta1 - eta1 reduction update factor
724: . -tao_ntl_eta2 - eta2 reduction update factor
725: . -tao_ntl_eta3 - eta3 reduction update factor
726: . -tao_ntl_eta4 - eta4 reduction update factor
727: . -tao_ntl_alpha1 - alpha1 reduction update factor
728: . -tao_ntl_alpha2 - alpha2 reduction update factor
729: . -tao_ntl_alpha3 - alpha3 reduction update factor
730: . -tao_ntl_alpha4 - alpha4 reduction update factor
731: . -tao_ntl_alpha4 - alpha4 reduction update factor
732: . -tao_ntl_mu1 - mu1 interpolation update
733: . -tao_ntl_mu2 - mu2 interpolation update
734: . -tao_ntl_gamma1 - gamma1 interpolcation update
735: . -tao_ntl_gamma2 - gamma2 interpolcation update
736: . -tao_ntl_gamma3 - gamma3 interpolcation update
737: . -tao_ntl_gamma4 - gamma4 interpolation update
738: - -tao_ntl_theta - theta1 interpolation update
740: Level: beginner
741: M*/
742: PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao)
743: {
744: TAO_NTL *tl;
745: const char *morethuente_type = TAOLINESEARCHMT;
747: PetscFunctionBegin;
748: PetscCall(PetscNew(&tl));
749: tao->ops->setup = TaoSetUp_NTL;
750: tao->ops->solve = TaoSolve_NTL;
751: tao->ops->view = TaoView_NTL;
752: tao->ops->setfromoptions = TaoSetFromOptions_NTL;
753: tao->ops->destroy = TaoDestroy_NTL;
755: /* Override default settings (unless already changed) */
756: if (!tao->max_it_changed) tao->max_it = 50;
757: if (!tao->trust0_changed) tao->trust0 = 100.0;
759: tao->data = (void *)tl;
761: /* Default values for trust-region radius update based on steplength */
762: tl->nu1 = 0.25;
763: tl->nu2 = 0.50;
764: tl->nu3 = 1.00;
765: tl->nu4 = 1.25;
767: tl->omega1 = 0.25;
768: tl->omega2 = 0.50;
769: tl->omega3 = 1.00;
770: tl->omega4 = 2.00;
771: tl->omega5 = 4.00;
773: /* Default values for trust-region radius update based on reduction */
774: tl->eta1 = 1.0e-4;
775: tl->eta2 = 0.25;
776: tl->eta3 = 0.50;
777: tl->eta4 = 0.90;
779: tl->alpha1 = 0.25;
780: tl->alpha2 = 0.50;
781: tl->alpha3 = 1.00;
782: tl->alpha4 = 2.00;
783: tl->alpha5 = 4.00;
785: /* Default values for trust-region radius update based on interpolation */
786: tl->mu1 = 0.10;
787: tl->mu2 = 0.50;
789: tl->gamma1 = 0.25;
790: tl->gamma2 = 0.50;
791: tl->gamma3 = 2.00;
792: tl->gamma4 = 4.00;
794: tl->theta = 0.05;
796: /* Default values for trust region initialization based on interpolation */
797: tl->mu1_i = 0.35;
798: tl->mu2_i = 0.50;
800: tl->gamma1_i = 0.0625;
801: tl->gamma2_i = 0.5;
802: tl->gamma3_i = 2.0;
803: tl->gamma4_i = 5.0;
805: tl->theta_i = 0.25;
807: /* Remaining parameters */
808: tl->min_radius = 1.0e-10;
809: tl->max_radius = 1.0e10;
810: tl->epsilon = 1.0e-6;
812: tl->init_type = NTL_INIT_INTERPOLATION;
813: tl->update_type = NTL_UPDATE_REDUCTION;
815: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
816: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
817: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
818: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
819: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
820: PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
821: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
822: PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
823: PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_ntl_"));
824: PetscCall(KSPSetType(tao->ksp, KSPSTCG));
825: PetscFunctionReturn(PETSC_SUCCESS);
826: }