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