Actual source code: ntl.c
petsc-3.9.4 2018-09-11
1: #include <../src/tao/matrix/lmvmmat.h>
2: #include <../src/tao/unconstrained/impls/ntl/ntlimpl.h>
4: #include <petscksp.h>
6: #define NTL_PC_NONE 0
7: #define NTL_PC_AHESS 1
8: #define NTL_PC_BFGS 2
9: #define NTL_PC_PETSC 3
10: #define NTL_PC_TYPES 4
12: #define BFGS_SCALE_AHESS 0
13: #define BFGS_SCALE_BFGS 1
14: #define BFGS_SCALE_TYPES 2
16: #define NTL_INIT_CONSTANT 0
17: #define NTL_INIT_DIRECTION 1
18: #define NTL_INIT_INTERPOLATION 2
19: #define NTL_INIT_TYPES 3
21: #define NTL_UPDATE_REDUCTION 0
22: #define NTL_UPDATE_INTERPOLATION 1
23: #define NTL_UPDATE_TYPES 2
25: static const char *NTL_PC[64] = {"none","ahess","bfgs","petsc"};
27: static const char *BFGS_SCALE[64] = {"ahess","bfgs"};
29: static const char *NTL_INIT[64] = {"constant","direction","interpolation"};
31: static const char *NTL_UPDATE[64] = {"reduction","interpolation"};
33: /* Routine for BFGS preconditioner */
35: static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x)
36: {
38: Mat M;
44: PCShellGetContext(pc,(void**)&M);
45: MatLMVMSolve(M, b, x);
46: return(0);
47: }
49: /* Implements Newton's Method with a trust-region, line-search approach for
50: solving unconstrained minimization problems. A More'-Thuente line search
51: is used to guarantee that the bfgs preconditioner remains positive
52: definite. */
54: #define NTL_NEWTON 0
55: #define NTL_BFGS 1
56: #define NTL_SCALED_GRADIENT 2
57: #define NTL_GRADIENT 3
59: static PetscErrorCode TaoSolve_NTL(Tao tao)
60: {
61: TAO_NTL *tl = (TAO_NTL *)tao->data;
62: KSPType ksp_type;
63: PetscBool is_nash,is_stcg,is_gltr;
64: KSPConvergedReason ksp_reason;
65: PC pc;
66: TaoLineSearchConvergedReason ls_reason;
68: PetscReal fmin, ftrial, prered, actred, kappa, sigma;
69: PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius;
70: PetscReal f, fold, gdx, gnorm;
71: PetscReal step = 1.0;
73: PetscReal delta;
74: PetscReal norm_d = 0.0;
75: PetscErrorCode ierr;
76: PetscInt stepType;
77: PetscInt its;
79: PetscInt bfgsUpdates = 0;
80: PetscInt needH;
82: PetscInt i_max = 5;
83: PetscInt j_max = 1;
84: PetscInt i, j, n, N;
86: PetscInt tr_reject;
89: if (tao->XL || tao->XU || tao->ops->computebounds) {
90: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n");
91: }
93: KSPGetType(tao->ksp,&ksp_type);
94: PetscStrcmp(ksp_type,KSPCGNASH,&is_nash);
95: PetscStrcmp(ksp_type,KSPCGSTCG,&is_stcg);
96: PetscStrcmp(ksp_type,KSPCGGLTR,&is_gltr);
97: if (!is_nash && !is_stcg && !is_gltr) {
98: SETERRQ(PETSC_COMM_SELF,1,"TAO_NTR requires nash, stcg, or gltr for the KSP");
99: }
101: /* Initialize the radius and modify if it is too large or small */
102: tao->trust = tao->trust0;
103: tao->trust = PetscMax(tao->trust, tl->min_radius);
104: tao->trust = PetscMin(tao->trust, tl->max_radius);
106: if (NTL_PC_BFGS == tl->pc_type && !tl->M) {
107: VecGetLocalSize(tao->solution,&n);
108: VecGetSize(tao->solution,&N);
109: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tl->M);
110: MatLMVMAllocateVectors(tl->M,tao->solution);
111: }
113: /* Check convergence criteria */
114: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
115: VecNorm(tao->gradient, NORM_2, &gnorm);
116: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
117: needH = 1;
119: tao->reason = TAO_CONTINUE_ITERATING;
120: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
121: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
122: (*tao->ops->convergencetest)(tao,tao->cnvP);
123: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
125: /* Create vectors for the limited memory preconditioner */
126: if ((NTL_PC_BFGS == tl->pc_type) && (BFGS_SCALE_BFGS != tl->bfgs_scale_type)) {
127: if (!tl->Diag) {
128: VecDuplicate(tao->solution, &tl->Diag);
129: }
130: }
132: /* Modify the preconditioner to use the bfgs approximation */
133: KSPGetPC(tao->ksp, &pc);
134: switch(tl->pc_type) {
135: case NTL_PC_NONE:
136: PCSetType(pc, PCNONE);
137: PCSetFromOptions(pc);
138: break;
140: case NTL_PC_AHESS:
141: PCSetType(pc, PCJACOBI);
142: PCSetFromOptions(pc);
143: PCJacobiSetUseAbs(pc,PETSC_TRUE);
144: break;
146: case NTL_PC_BFGS:
147: PCSetType(pc, PCSHELL);
148: PCSetFromOptions(pc);
149: PCShellSetName(pc, "bfgs");
150: PCShellSetContext(pc, tl->M);
151: PCShellSetApply(pc, MatLMVMSolveShell);
152: break;
154: default:
155: /* Use the pc method set by pc_type */
156: break;
157: }
159: /* Initialize trust-region radius */
160: switch(tl->init_type) {
161: case NTL_INIT_CONSTANT:
162: /* Use the initial radius specified */
163: break;
165: case NTL_INIT_INTERPOLATION:
166: /* Use the initial radius specified */
167: max_radius = 0.0;
169: for (j = 0; j < j_max; ++j) {
170: fmin = f;
171: sigma = 0.0;
173: if (needH) {
174: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
175: needH = 0;
176: }
178: for (i = 0; i < i_max; ++i) {
179: VecCopy(tao->solution, tl->W);
180: VecAXPY(tl->W, -tao->trust/gnorm, tao->gradient);
182: TaoComputeObjective(tao, tl->W, &ftrial);
183: if (PetscIsInfOrNanReal(ftrial)) {
184: tau = tl->gamma1_i;
185: } else {
186: if (ftrial < fmin) {
187: fmin = ftrial;
188: sigma = -tao->trust / gnorm;
189: }
191: MatMult(tao->hessian, tao->gradient, tao->stepdirection);
192: VecDot(tao->gradient, tao->stepdirection, &prered);
194: prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
195: actred = f - ftrial;
196: if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) {
197: kappa = 1.0;
198: } else {
199: kappa = actred / prered;
200: }
202: tau_1 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred);
203: tau_2 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred);
204: tau_min = PetscMin(tau_1, tau_2);
205: tau_max = PetscMax(tau_1, tau_2);
207: if (PetscAbsScalar(kappa - 1.0) <= tl->mu1_i) {
208: /* Great agreement */
209: max_radius = PetscMax(max_radius, tao->trust);
211: if (tau_max < 1.0) {
212: tau = tl->gamma3_i;
213: } else if (tau_max > tl->gamma4_i) {
214: tau = tl->gamma4_i;
215: } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) {
216: tau = tau_1;
217: } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) {
218: tau = tau_2;
219: } else {
220: tau = tau_max;
221: }
222: } else if (PetscAbsScalar(kappa - 1.0) <= tl->mu2_i) {
223: /* Good agreement */
224: max_radius = PetscMax(max_radius, tao->trust);
226: if (tau_max < tl->gamma2_i) {
227: tau = tl->gamma2_i;
228: } else if (tau_max > tl->gamma3_i) {
229: tau = tl->gamma3_i;
230: } else {
231: tau = tau_max;
232: }
233: } else {
234: /* Not good agreement */
235: if (tau_min > 1.0) {
236: tau = tl->gamma2_i;
237: } else if (tau_max < tl->gamma1_i) {
238: tau = tl->gamma1_i;
239: } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) {
240: tau = tl->gamma1_i;
241: } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) {
242: tau = tau_1;
243: } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) {
244: tau = tau_2;
245: } else {
246: tau = tau_max;
247: }
248: }
249: }
250: tao->trust = tau * tao->trust;
251: }
253: if (fmin < f) {
254: f = fmin;
255: VecAXPY(tao->solution, sigma, tao->gradient);
256: TaoComputeGradient(tao, tao->solution, tao->gradient);
258: VecNorm(tao->gradient, NORM_2, &gnorm);
259: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
260: needH = 1;
262: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
263: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
264: (*tao->ops->convergencetest)(tao,tao->cnvP);
265: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
266: }
267: }
268: tao->trust = PetscMax(tao->trust, max_radius);
270: /* Modify the radius if it is too large or small */
271: tao->trust = PetscMax(tao->trust, tl->min_radius);
272: tao->trust = PetscMin(tao->trust, tl->max_radius);
273: break;
275: default:
276: /* Norm of the first direction will initialize radius */
277: tao->trust = 0.0;
278: break;
279: }
281: /* Set initial scaling for the BFGS preconditioner
282: This step is done after computing the initial trust-region radius
283: since the function value may have decreased */
284: if (NTL_PC_BFGS == tl->pc_type) {
285: if (f != 0.0) {
286: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
287: } else {
288: delta = 2.0 / (gnorm*gnorm);
289: }
290: MatLMVMSetDelta(tl->M, delta);
291: }
293: /* Set counter for gradient/reset steps */
294: tl->ntrust = 0;
295: tl->newt = 0;
296: tl->bfgs = 0;
297: tl->sgrad = 0;
298: tl->grad = 0;
300: /* Have not converged; continue with Newton method */
301: while (tao->reason == TAO_CONTINUE_ITERATING) {
302: ++tao->niter;
303: tao->ksp_its=0;
304: /* Compute the Hessian */
305: if (needH) {
306: TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);
307: }
309: if (NTL_PC_BFGS == tl->pc_type) {
310: if (BFGS_SCALE_AHESS == tl->bfgs_scale_type) {
311: /* Obtain diagonal for the bfgs preconditioner */
312: MatGetDiagonal(tao->hessian, tl->Diag);
313: VecAbs(tl->Diag);
314: VecReciprocal(tl->Diag);
315: MatLMVMSetScale(tl->M, tl->Diag);
316: }
318: /* Update the limited memory preconditioner */
319: MatLMVMUpdate(tl->M,tao->solution, tao->gradient);
320: ++bfgsUpdates;
321: }
322: KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);
323: /* Solve the Newton system of equations */
324: KSPCGSetRadius(tao->ksp,tl->max_radius);
325: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
326: KSPGetIterationNumber(tao->ksp,&its);
327: tao->ksp_its+=its;
328: tao->ksp_tot_its+=its;
329: KSPCGGetNormD(tao->ksp, &norm_d);
331: if (0.0 == tao->trust) {
332: /* Radius was uninitialized; use the norm of the direction */
333: if (norm_d > 0.0) {
334: tao->trust = norm_d;
336: /* Modify the radius if it is too large or small */
337: tao->trust = PetscMax(tao->trust, tl->min_radius);
338: tao->trust = PetscMin(tao->trust, tl->max_radius);
339: } else {
340: /* The direction was bad; set radius to default value and re-solve
341: the trust-region subproblem to get a direction */
342: tao->trust = tao->trust0;
344: /* Modify the radius if it is too large or small */
345: tao->trust = PetscMax(tao->trust, tl->min_radius);
346: tao->trust = PetscMin(tao->trust, tl->max_radius);
348: KSPCGSetRadius(tao->ksp,tl->max_radius);
349: KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);
350: KSPGetIterationNumber(tao->ksp,&its);
351: tao->ksp_its+=its;
352: tao->ksp_tot_its+=its;
353: KSPCGGetNormD(tao->ksp, &norm_d);
355: if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero");
356: }
357: }
359: VecScale(tao->stepdirection, -1.0);
360: KSPGetConvergedReason(tao->ksp, &ksp_reason);
361: if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (NTL_PC_BFGS == tl->pc_type) && (bfgsUpdates > 1)) {
362: /* Preconditioner is numerically indefinite; reset the
363: approximate if using BFGS preconditioning. */
365: if (f != 0.0) {
366: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
367: } else {
368: delta = 2.0 / (gnorm*gnorm);
369: }
370: MatLMVMSetDelta(tl->M, delta);
371: MatLMVMReset(tl->M);
372: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
373: bfgsUpdates = 1;
374: }
376: /* Check trust-region reduction conditions */
377: tr_reject = 0;
378: if (NTL_UPDATE_REDUCTION == tl->update_type) {
379: /* Get predicted reduction */
380: KSPCGGetObjFcn(tao->ksp,&prered);
381: if (prered >= 0.0) {
382: /* The predicted reduction has the wrong sign. This cannot
383: happen in infinite precision arithmetic. Step should
384: be rejected! */
385: tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
386: tr_reject = 1;
387: } else {
388: /* Compute trial step and function value */
389: VecCopy(tao->solution, tl->W);
390: VecAXPY(tl->W, 1.0, tao->stepdirection);
391: TaoComputeObjective(tao, tl->W, &ftrial);
393: if (PetscIsInfOrNanReal(ftrial)) {
394: tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
395: tr_reject = 1;
396: } else {
397: /* Compute and actual reduction */
398: actred = f - ftrial;
399: prered = -prered;
400: if ((PetscAbsScalar(actred) <= tl->epsilon) &&
401: (PetscAbsScalar(prered) <= tl->epsilon)) {
402: kappa = 1.0;
403: } else {
404: kappa = actred / prered;
405: }
407: /* Accept of reject the step and update radius */
408: if (kappa < tl->eta1) {
409: /* Reject the step */
410: tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d);
411: tr_reject = 1;
412: } else {
413: /* Accept the step */
414: if (kappa < tl->eta2) {
415: /* Marginal bad step */
416: tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d);
417: } else if (kappa < tl->eta3) {
418: /* Reasonable step */
419: tao->trust = tl->alpha3 * tao->trust;
420: } else if (kappa < tl->eta4) {
421: /* Good step */
422: tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust);
423: } else {
424: /* Very good step */
425: tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust);
426: }
427: }
428: }
429: }
430: } else {
431: /* Get predicted reduction */
432: KSPCGGetObjFcn(tao->ksp,&prered);
433: if (prered >= 0.0) {
434: /* The predicted reduction has the wrong sign. This cannot
435: happen in infinite precision arithmetic. Step should
436: be rejected! */
437: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
438: tr_reject = 1;
439: } else {
440: VecCopy(tao->solution, tl->W);
441: VecAXPY(tl->W, 1.0, tao->stepdirection);
442: TaoComputeObjective(tao, tl->W, &ftrial);
443: if (PetscIsInfOrNanReal(ftrial)) {
444: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
445: tr_reject = 1;
446: } else {
447: VecDot(tao->gradient, tao->stepdirection, &gdx);
449: actred = f - ftrial;
450: prered = -prered;
451: if ((PetscAbsScalar(actred) <= tl->epsilon) &&
452: (PetscAbsScalar(prered) <= tl->epsilon)) {
453: kappa = 1.0;
454: } else {
455: kappa = actred / prered;
456: }
458: tau_1 = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred);
459: tau_2 = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred);
460: tau_min = PetscMin(tau_1, tau_2);
461: tau_max = PetscMax(tau_1, tau_2);
463: if (kappa >= 1.0 - tl->mu1) {
464: /* Great agreement; accept step and update radius */
465: if (tau_max < 1.0) {
466: tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
467: } else if (tau_max > tl->gamma4) {
468: tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d);
469: } else {
470: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
471: }
472: } else if (kappa >= 1.0 - tl->mu2) {
473: /* Good agreement */
475: if (tau_max < tl->gamma2) {
476: tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
477: } else if (tau_max > tl->gamma3) {
478: tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d);
479: } else if (tau_max < 1.0) {
480: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
481: } else {
482: tao->trust = PetscMax(tao->trust, tau_max * norm_d);
483: }
484: } else {
485: /* Not good agreement */
486: if (tau_min > 1.0) {
487: tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d);
488: } else if (tau_max < tl->gamma1) {
489: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
490: } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) {
491: tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d);
492: } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) {
493: tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
494: } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) {
495: tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
496: } else {
497: tao->trust = tau_max * PetscMin(tao->trust, norm_d);
498: }
499: tr_reject = 1;
500: }
501: }
502: }
503: }
505: if (tr_reject) {
506: /* The trust-region constraints rejected the step. Apply a linesearch.
507: Check for descent direction. */
508: VecDot(tao->stepdirection, tao->gradient, &gdx);
509: if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) {
510: /* Newton step is not descent or direction produced Inf or NaN */
512: if (NTL_PC_BFGS != tl->pc_type) {
513: /* We don't have the bfgs matrix around and updated
514: Must use gradient direction in this case */
515: VecCopy(tao->gradient, tao->stepdirection);
516: VecScale(tao->stepdirection, -1.0);
517: ++tl->grad;
518: stepType = NTL_GRADIENT;
519: } else {
520: /* Attempt to use the BFGS direction */
521: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
522: VecScale(tao->stepdirection, -1.0);
524: /* Check for success (descent direction) */
525: VecDot(tao->stepdirection, tao->gradient, &gdx);
526: if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) {
527: /* BFGS direction is not descent or direction produced not a number
528: We can assert bfgsUpdates > 1 in this case because
529: the first solve produces the scaled gradient direction,
530: which is guaranteed to be descent */
532: /* Use steepest descent direction (scaled) */
533: if (f != 0.0) {
534: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
535: } else {
536: delta = 2.0 / (gnorm*gnorm);
537: }
538: MatLMVMSetDelta(tl->M, delta);
539: MatLMVMReset(tl->M);
540: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
541: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
542: VecScale(tao->stepdirection, -1.0);
544: bfgsUpdates = 1;
545: ++tl->sgrad;
546: stepType = NTL_SCALED_GRADIENT;
547: } else {
548: if (1 == bfgsUpdates) {
549: /* The first BFGS direction is always the scaled gradient */
550: ++tl->sgrad;
551: stepType = NTL_SCALED_GRADIENT;
552: } else {
553: ++tl->bfgs;
554: stepType = NTL_BFGS;
555: }
556: }
557: }
558: } else {
559: /* Computed Newton step is descent */
560: ++tl->newt;
561: stepType = NTL_NEWTON;
562: }
564: /* Perform the linesearch */
565: fold = f;
566: VecCopy(tao->solution, tl->Xold);
567: VecCopy(tao->gradient, tl->Gold);
569: step = 1.0;
570: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);
571: TaoAddLineSearchCounts(tao);
573: while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */
574: /* Linesearch failed */
575: f = fold;
576: VecCopy(tl->Xold, tao->solution);
577: VecCopy(tl->Gold, tao->gradient);
579: switch(stepType) {
580: case NTL_NEWTON:
581: /* Failed to obtain acceptable iterate with Newton step */
583: if (NTL_PC_BFGS != tl->pc_type) {
584: /* We don't have the bfgs matrix around and being updated
585: Must use gradient direction in this case */
586: VecCopy(tao->gradient, tao->stepdirection);
587: ++tl->grad;
588: stepType = NTL_GRADIENT;
589: } else {
590: /* Attempt to use the BFGS direction */
591: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
594: /* Check for success (descent direction) */
595: VecDot(tao->stepdirection, tao->gradient, &gdx);
596: if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) {
597: /* BFGS direction is not descent or direction produced
598: not a number. We can assert bfgsUpdates > 1 in this case
599: Use steepest descent direction (scaled) */
601: if (f != 0.0) {
602: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
603: } else {
604: delta = 2.0 / (gnorm*gnorm);
605: }
606: MatLMVMSetDelta(tl->M, delta);
607: MatLMVMReset(tl->M);
608: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
609: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
611: bfgsUpdates = 1;
612: ++tl->sgrad;
613: stepType = NTL_SCALED_GRADIENT;
614: } else {
615: if (1 == bfgsUpdates) {
616: /* The first BFGS direction is always the scaled gradient */
617: ++tl->sgrad;
618: stepType = NTL_SCALED_GRADIENT;
619: } else {
620: ++tl->bfgs;
621: stepType = NTL_BFGS;
622: }
623: }
624: }
625: break;
627: case NTL_BFGS:
628: /* Can only enter if pc_type == NTL_PC_BFGS
629: Failed to obtain acceptable iterate with BFGS step
630: Attempt to use the scaled gradient direction */
632: if (f != 0.0) {
633: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
634: } else {
635: delta = 2.0 / (gnorm*gnorm);
636: }
637: MatLMVMSetDelta(tl->M, delta);
638: MatLMVMReset(tl->M);
639: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
640: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
642: bfgsUpdates = 1;
643: ++tl->sgrad;
644: stepType = NTL_SCALED_GRADIENT;
645: break;
647: case NTL_SCALED_GRADIENT:
648: /* Can only enter if pc_type == NTL_PC_BFGS
649: The scaled gradient step did not produce a new iterate;
650: attemp to use the gradient direction.
651: Need to make sure we are not using a different diagonal scaling */
652: MatLMVMSetScale(tl->M, tl->Diag);
653: MatLMVMSetDelta(tl->M, 1.0);
654: MatLMVMReset(tl->M);
655: MatLMVMUpdate(tl->M, tao->solution, tao->gradient);
656: MatLMVMSolve(tl->M, tao->gradient, tao->stepdirection);
658: bfgsUpdates = 1;
659: ++tl->grad;
660: stepType = NTL_GRADIENT;
661: break;
662: }
663: VecScale(tao->stepdirection, -1.0);
665: /* This may be incorrect; linesearch has values for stepmax and stepmin
666: that should be reset. */
667: step = 1.0;
668: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason);
669: TaoAddLineSearchCounts(tao);
670: }
672: if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) {
673: /* Failed to find an improving point */
674: f = fold;
675: VecCopy(tl->Xold, tao->solution);
676: VecCopy(tl->Gold, tao->gradient);
677: tao->trust = 0.0;
678: step = 0.0;
679: tao->reason = TAO_DIVERGED_LS_FAILURE;
680: break;
681: } else if (stepType == NTL_NEWTON) {
682: if (step < tl->nu1) {
683: /* Very bad step taken; reduce radius */
684: tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
685: } else if (step < tl->nu2) {
686: /* Reasonably bad step taken; reduce radius */
687: tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust);
688: } else if (step < tl->nu3) {
689: /* Reasonable step was taken; leave radius alone */
690: if (tl->omega3 < 1.0) {
691: tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust);
692: } else if (tl->omega3 > 1.0) {
693: tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust);
694: }
695: } else if (step < tl->nu4) {
696: /* Full step taken; increase the radius */
697: tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust);
698: } else {
699: /* More than full step taken; increase the radius */
700: tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust);
701: }
702: } else {
703: /* Newton step was not good; reduce the radius */
704: tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust);
705: }
706: } else {
707: /* Trust-region step is accepted */
708: VecCopy(tl->W, tao->solution);
709: f = ftrial;
710: TaoComputeGradient(tao, tao->solution, tao->gradient);
711: ++tl->ntrust;
712: }
714: /* The radius may have been increased; modify if it is too large */
715: tao->trust = PetscMin(tao->trust, tl->max_radius);
717: /* Check for converged */
718: VecNorm(tao->gradient, NORM_2, &gnorm);
719: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number");
720: needH = 1;
722: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
723: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
724: (*tao->ops->convergencetest)(tao,tao->cnvP);
725: }
726: return(0);
727: }
729: /* ---------------------------------------------------------- */
730: static PetscErrorCode TaoSetUp_NTL(Tao tao)
731: {
732: TAO_NTL *tl = (TAO_NTL *)tao->data;
736: if (!tao->gradient) {VecDuplicate(tao->solution, &tao->gradient); }
737: if (!tao->stepdirection) {VecDuplicate(tao->solution, &tao->stepdirection);}
738: if (!tl->W) { VecDuplicate(tao->solution, &tl->W);}
739: if (!tl->Xold) { VecDuplicate(tao->solution, &tl->Xold);}
740: if (!tl->Gold) { VecDuplicate(tao->solution, &tl->Gold);}
741: tl->Diag = 0;
742: tl->M = 0;
743: return(0);
744: }
746: /*------------------------------------------------------------*/
747: static PetscErrorCode TaoDestroy_NTL(Tao tao)
748: {
749: TAO_NTL *tl = (TAO_NTL *)tao->data;
753: if (tao->setupcalled) {
754: VecDestroy(&tl->W);
755: VecDestroy(&tl->Xold);
756: VecDestroy(&tl->Gold);
757: }
758: VecDestroy(&tl->Diag);
759: MatDestroy(&tl->M);
760: PetscFree(tao->data);
761: return(0);
762: }
764: /*------------------------------------------------------------*/
765: static PetscErrorCode TaoSetFromOptions_NTL(PetscOptionItems *PetscOptionsObject,Tao tao)
766: {
767: TAO_NTL *tl = (TAO_NTL *)tao->data;
771: PetscOptionsHead(PetscOptionsObject,"Newton trust region with line search method for unconstrained optimization");
772: PetscOptionsEList("-tao_ntl_pc_type", "pc type", "", NTL_PC, NTL_PC_TYPES, NTL_PC[tl->pc_type], &tl->pc_type,NULL);
773: PetscOptionsEList("-tao_ntl_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tl->bfgs_scale_type], &tl->bfgs_scale_type,NULL);
774: PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type,NULL);
775: PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type,NULL);
776: PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1,NULL);
777: PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2,NULL);
778: PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3,NULL);
779: PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4,NULL);
780: PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1,NULL);
781: PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2,NULL);
782: PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3,NULL);
783: PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4,NULL);
784: PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5,NULL);
785: PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1,NULL);
786: PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2,NULL);
787: PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3,NULL);
788: PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4,NULL);
789: PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1,NULL);
790: PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2,NULL);
791: PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3,NULL);
792: PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4,NULL);
793: PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5,NULL);
794: PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i,NULL);
795: PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i,NULL);
796: PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i,NULL);
797: PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i,NULL);
798: PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i,NULL);
799: PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i,NULL);
800: PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i,NULL);
801: PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1,NULL);
802: PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2,NULL);
803: PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1,NULL);
804: PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2,NULL);
805: PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3,NULL);
806: PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4,NULL);
807: PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta,NULL);
808: PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius,NULL);
809: PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius,NULL);
810: PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon,NULL);
811: PetscOptionsTail();
812: TaoLineSearchSetFromOptions(tao->linesearch);
813: KSPSetFromOptions(tao->ksp);
814: return(0);
815: }
817: /*------------------------------------------------------------*/
818: static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer)
819: {
820: TAO_NTL *tl = (TAO_NTL *)tao->data;
821: PetscInt nrejects;
822: PetscBool isascii;
826: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
827: if (isascii) {
828: PetscViewerASCIIPushTab(viewer);
829: if (NTL_PC_BFGS == tl->pc_type && tl->M) {
830: MatLMVMGetRejects(tl->M, &nrejects);
831: PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", nrejects);
832: }
833: PetscViewerASCIIPrintf(viewer, "Trust-region steps: %D\n", tl->ntrust);
834: PetscViewerASCIIPrintf(viewer, "Newton search steps: %D\n", tl->newt);
835: PetscViewerASCIIPrintf(viewer, "BFGS search steps: %D\n", tl->bfgs);
836: PetscViewerASCIIPrintf(viewer, "Scaled gradient search steps: %D\n", tl->sgrad);
837: PetscViewerASCIIPrintf(viewer, "Gradient search steps: %D\n", tl->grad);
838: PetscViewerASCIIPopTab(viewer);
839: }
840: return(0);
841: }
843: /* ---------------------------------------------------------- */
844: /*MC
845: TAONTR - Newton's method with trust region and linesearch
846: for unconstrained minimization.
847: At each iteration, the Newton trust region method solves the system for d
848: and performs a line search in the d direction:
850: min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k
852: Options Database Keys:
853: + -tao_ntl_pc_type - "none","ahess","bfgs","petsc"
854: . -tao_ntl_bfgs_scale_type - type of scaling with bfgs pc, "ahess" or "bfgs"
855: . -tao_ntl_init_type - "constant","direction","interpolation"
856: . -tao_ntl_update_type - "reduction","interpolation"
857: . -tao_ntl_min_radius - lower bound on trust region radius
858: . -tao_ntl_max_radius - upper bound on trust region radius
859: . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction
860: . -tao_ntl_mu1_i - mu1 interpolation init factor
861: . -tao_ntl_mu2_i - mu2 interpolation init factor
862: . -tao_ntl_gamma1_i - gamma1 interpolation init factor
863: . -tao_ntl_gamma2_i - gamma2 interpolation init factor
864: . -tao_ntl_gamma3_i - gamma3 interpolation init factor
865: . -tao_ntl_gamma4_i - gamma4 interpolation init factor
866: . -tao_ntl_theta_i - thetha1 interpolation init factor
867: . -tao_ntl_eta1 - eta1 reduction update factor
868: . -tao_ntl_eta2 - eta2 reduction update factor
869: . -tao_ntl_eta3 - eta3 reduction update factor
870: . -tao_ntl_eta4 - eta4 reduction update factor
871: . -tao_ntl_alpha1 - alpha1 reduction update factor
872: . -tao_ntl_alpha2 - alpha2 reduction update factor
873: . -tao_ntl_alpha3 - alpha3 reduction update factor
874: . -tao_ntl_alpha4 - alpha4 reduction update factor
875: . -tao_ntl_alpha4 - alpha4 reduction update factor
876: . -tao_ntl_mu1 - mu1 interpolation update
877: . -tao_ntl_mu2 - mu2 interpolation update
878: . -tao_ntl_gamma1 - gamma1 interpolcation update
879: . -tao_ntl_gamma2 - gamma2 interpolcation update
880: . -tao_ntl_gamma3 - gamma3 interpolcation update
881: . -tao_ntl_gamma4 - gamma4 interpolation update
882: - -tao_ntl_theta - theta1 interpolation update
884: Level: beginner
885: M*/
887: PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao)
888: {
889: TAO_NTL *tl;
891: const char *morethuente_type = TAOLINESEARCHMT;
894: PetscNewLog(tao,&tl);
895: tao->ops->setup = TaoSetUp_NTL;
896: tao->ops->solve = TaoSolve_NTL;
897: tao->ops->view = TaoView_NTL;
898: tao->ops->setfromoptions = TaoSetFromOptions_NTL;
899: tao->ops->destroy = TaoDestroy_NTL;
901: /* Override default settings (unless already changed) */
902: if (!tao->max_it_changed) tao->max_it = 50;
903: if (!tao->trust0_changed) tao->trust0 = 100.0;
905: tao->data = (void*)tl;
907: /* Default values for trust-region radius update based on steplength */
908: tl->nu1 = 0.25;
909: tl->nu2 = 0.50;
910: tl->nu3 = 1.00;
911: tl->nu4 = 1.25;
913: tl->omega1 = 0.25;
914: tl->omega2 = 0.50;
915: tl->omega3 = 1.00;
916: tl->omega4 = 2.00;
917: tl->omega5 = 4.00;
919: /* Default values for trust-region radius update based on reduction */
920: tl->eta1 = 1.0e-4;
921: tl->eta2 = 0.25;
922: tl->eta3 = 0.50;
923: tl->eta4 = 0.90;
925: tl->alpha1 = 0.25;
926: tl->alpha2 = 0.50;
927: tl->alpha3 = 1.00;
928: tl->alpha4 = 2.00;
929: tl->alpha5 = 4.00;
931: /* Default values for trust-region radius update based on interpolation */
932: tl->mu1 = 0.10;
933: tl->mu2 = 0.50;
935: tl->gamma1 = 0.25;
936: tl->gamma2 = 0.50;
937: tl->gamma3 = 2.00;
938: tl->gamma4 = 4.00;
940: tl->theta = 0.05;
942: /* Default values for trust region initialization based on interpolation */
943: tl->mu1_i = 0.35;
944: tl->mu2_i = 0.50;
946: tl->gamma1_i = 0.0625;
947: tl->gamma2_i = 0.5;
948: tl->gamma3_i = 2.0;
949: tl->gamma4_i = 5.0;
951: tl->theta_i = 0.25;
953: /* Remaining parameters */
954: tl->min_radius = 1.0e-10;
955: tl->max_radius = 1.0e10;
956: tl->epsilon = 1.0e-6;
958: tl->pc_type = NTL_PC_BFGS;
959: tl->bfgs_scale_type = BFGS_SCALE_AHESS;
960: tl->init_type = NTL_INIT_INTERPOLATION;
961: tl->update_type = NTL_UPDATE_REDUCTION;
963: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
964: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
965: TaoLineSearchSetType(tao->linesearch, morethuente_type);
966: TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
967: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
968: KSPCreate(((PetscObject)tao)->comm,&tao->ksp);
969: PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
970: KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix);
971: KSPSetType(tao->ksp,KSPCGSTCG);
972: return(0);
973: }