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