Actual source code: ntr.c

      1: #include <../src/tao/unconstrained/impls/ntr/ntrimpl.h>

      3: #include <petscksp.h>

      5: #define NTR_INIT_CONSTANT      0
      6: #define NTR_INIT_DIRECTION     1
      7: #define NTR_INIT_INTERPOLATION 2
      8: #define NTR_INIT_TYPES         3

     10: #define NTR_UPDATE_REDUCTION     0
     11: #define NTR_UPDATE_INTERPOLATION 1
     12: #define NTR_UPDATE_TYPES         2

     14: static const char *NTR_INIT[64] = {"constant", "direction", "interpolation"};

     16: static const char *NTR_UPDATE[64] = {"reduction", "interpolation"};

     18: /*
     19:    TaoSolve_NTR - Implements Newton's Method with a trust region approach
     20:    for solving unconstrained minimization problems.

     22:    The basic algorithm is taken from MINPACK-2 (dstrn).

     24:    TaoSolve_NTR computes a local minimizer of a twice differentiable function
     25:    f by applying a trust region variant of Newton's method.  At each stage
     26:    of the algorithm, we use the prconditioned conjugate gradient method to
     27:    determine an approximate minimizer of the quadratic equation

     29:         q(s) = <s, Hs + g>

     31:    subject to the trust region constraint

     33:         || s ||_M <= radius,

     35:    where radius is the trust region radius and M is a symmetric positive
     36:    definite matrix (the preconditioner).  Here g is the gradient and H
     37:    is the Hessian matrix.

     39:    Note:  TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG,
     40:           or KSPGLTR.  Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this
     41:           routine regardless of what the user may have previously specified.
     42: */
     43: static PetscErrorCode TaoSolve_NTR(Tao tao)
     44: {
     45:   TAO_NTR           *tr = (TAO_NTR *)tao->data;
     46:   KSPType            ksp_type;
     47:   PetscBool          is_nash, is_stcg, is_gltr, is_bfgs, is_jacobi, is_symmetric, sym_set;
     48:   KSPConvergedReason ksp_reason;
     49:   PC                 pc;
     50:   PetscReal          fmin, ftrial, prered, actred, kappa, sigma, beta;
     51:   PetscReal          tau, tau_1, tau_2, tau_max, tau_min, max_radius;
     52:   PetscReal          f, gnorm;

     54:   PetscReal norm_d;
     55:   PetscInt  needH;

     57:   PetscInt i_max = 5;
     58:   PetscInt j_max = 1;
     59:   PetscInt i, j, N, n, its;

     61:   PetscFunctionBegin;
     62:   if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n"));

     64:   PetscCall(KSPGetType(tao->ksp, &ksp_type));
     65:   PetscCall(PetscStrcmp(ksp_type, KSPNASH, &is_nash));
     66:   PetscCall(PetscStrcmp(ksp_type, KSPSTCG, &is_stcg));
     67:   PetscCall(PetscStrcmp(ksp_type, KSPGLTR, &is_gltr));
     68:   PetscCheck(is_nash || is_stcg || is_gltr, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "TAO_NTR requires nash, stcg, or gltr for the KSP");

     70:   /* Initialize the radius and modify if it is too large or small */
     71:   tao->trust = tao->trust0;
     72:   tao->trust = PetscMax(tao->trust, tr->min_radius);
     73:   tao->trust = PetscMin(tao->trust, tr->max_radius);

     75:   /* Allocate the vectors needed for the BFGS approximation */
     76:   PetscCall(KSPGetPC(tao->ksp, &pc));
     77:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs));
     78:   PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi));
     79:   if (is_bfgs) {
     80:     tr->bfgs_pre = pc;
     81:     PetscCall(PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M));
     82:     PetscCall(VecGetLocalSize(tao->solution, &n));
     83:     PetscCall(VecGetSize(tao->solution, &N));
     84:     PetscCall(MatSetSizes(tr->M, n, n, N, N));
     85:     PetscCall(MatLMVMAllocate(tr->M, tao->solution, tao->gradient));
     86:     PetscCall(MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric));
     87:     PetscCheck(sym_set && is_symmetric, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric.");
     88:   } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc, PETSC_TRUE));

     90:   /* Check convergence criteria */
     91:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
     92:   PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
     93:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
     94:   needH = 1;

     96:   tao->reason = TAO_CONTINUE_ITERATING;
     97:   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
     98:   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0));
     99:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    100:   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);

    102:   /*  Initialize trust-region radius */
    103:   switch (tr->init_type) {
    104:   case NTR_INIT_CONSTANT:
    105:     /*  Use the initial radius specified */
    106:     break;

    108:   case NTR_INIT_INTERPOLATION:
    109:     /*  Use the initial radius specified */
    110:     max_radius = 0.0;

    112:     for (j = 0; j < j_max; ++j) {
    113:       fmin  = f;
    114:       sigma = 0.0;

    116:       if (needH) {
    117:         PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
    118:         needH = 0;
    119:       }

    121:       for (i = 0; i < i_max; ++i) {
    122:         PetscCall(VecCopy(tao->solution, tr->W));
    123:         PetscCall(VecAXPY(tr->W, -tao->trust / gnorm, tao->gradient));
    124:         PetscCall(TaoComputeObjective(tao, tr->W, &ftrial));

    126:         if (PetscIsInfOrNanReal(ftrial)) {
    127:           tau = tr->gamma1_i;
    128:         } else {
    129:           if (ftrial < fmin) {
    130:             fmin  = ftrial;
    131:             sigma = -tao->trust / gnorm;
    132:           }

    134:           PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection));
    135:           PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered));

    137:           prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm));
    138:           actred = f - ftrial;
    139:           if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
    140:             kappa = 1.0;
    141:           } else {
    142:             kappa = actred / prered;
    143:           }

    145:           tau_1   = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred);
    146:           tau_2   = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred);
    147:           tau_min = PetscMin(tau_1, tau_2);
    148:           tau_max = PetscMax(tau_1, tau_2);

    150:           if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu1_i) {
    151:             /*  Great agreement */
    152:             max_radius = PetscMax(max_radius, tao->trust);

    154:             if (tau_max < 1.0) {
    155:               tau = tr->gamma3_i;
    156:             } else if (tau_max > tr->gamma4_i) {
    157:               tau = tr->gamma4_i;
    158:             } else {
    159:               tau = tau_max;
    160:             }
    161:           } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu2_i) {
    162:             /*  Good agreement */
    163:             max_radius = PetscMax(max_radius, tao->trust);

    165:             if (tau_max < tr->gamma2_i) {
    166:               tau = tr->gamma2_i;
    167:             } else if (tau_max > tr->gamma3_i) {
    168:               tau = tr->gamma3_i;
    169:             } else {
    170:               tau = tau_max;
    171:             }
    172:           } else {
    173:             /*  Not good agreement */
    174:             if (tau_min > 1.0) {
    175:               tau = tr->gamma2_i;
    176:             } else if (tau_max < tr->gamma1_i) {
    177:               tau = tr->gamma1_i;
    178:             } else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) {
    179:               tau = tr->gamma1_i;
    180:             } else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) {
    181:               tau = tau_1;
    182:             } else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) {
    183:               tau = tau_2;
    184:             } else {
    185:               tau = tau_max;
    186:             }
    187:           }
    188:         }
    189:         tao->trust = tau * tao->trust;
    190:       }

    192:       if (fmin < f) {
    193:         f = fmin;
    194:         PetscCall(VecAXPY(tao->solution, sigma, tao->gradient));
    195:         PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));

    197:         PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
    198:         PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
    199:         needH = 1;

    201:         PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
    202:         PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0));
    203:         PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    204:         if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
    205:       }
    206:     }
    207:     tao->trust = PetscMax(tao->trust, max_radius);

    209:     /*  Modify the radius if it is too large or small */
    210:     tao->trust = PetscMax(tao->trust, tr->min_radius);
    211:     tao->trust = PetscMin(tao->trust, tr->max_radius);
    212:     break;

    214:   default:
    215:     /*  Norm of the first direction will initialize radius */
    216:     tao->trust = 0.0;
    217:     break;
    218:   }

    220:   /* Have not converged; continue with Newton method */
    221:   while (tao->reason == TAO_CONTINUE_ITERATING) {
    222:     /* Call general purpose update function */
    223:     if (tao->ops->update) {
    224:       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
    225:       PetscCall(TaoComputeObjective(tao, tao->solution, &f));
    226:     }
    227:     ++tao->niter;
    228:     tao->ksp_its = 0;
    229:     /* Compute the Hessian */
    230:     if (needH) {
    231:       PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
    232:       needH = 0;
    233:     }

    235:     if (tr->bfgs_pre) {
    236:       /* Update the limited memory preconditioner */
    237:       PetscCall(MatLMVMUpdate(tr->M, tao->solution, tao->gradient));
    238:     }

    240:     while (tao->reason == TAO_CONTINUE_ITERATING) {
    241:       PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre));

    243:       /* Solve the trust region subproblem */
    244:       PetscCall(KSPCGSetRadius(tao->ksp, tao->trust));
    245:       PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
    246:       PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    247:       tao->ksp_its += its;
    248:       tao->ksp_tot_its += its;
    249:       PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));

    251:       if (0.0 == tao->trust) {
    252:         /* Radius was uninitialized; use the norm of the direction */
    253:         if (norm_d > 0.0) {
    254:           tao->trust = norm_d;

    256:           /* Modify the radius if it is too large or small */
    257:           tao->trust = PetscMax(tao->trust, tr->min_radius);
    258:           tao->trust = PetscMin(tao->trust, tr->max_radius);
    259:         } else {
    260:           /* The direction was bad; set radius to default value and re-solve
    261:              the trust-region subproblem to get a direction */
    262:           tao->trust = tao->trust0;

    264:           /* Modify the radius if it is too large or small */
    265:           tao->trust = PetscMax(tao->trust, tr->min_radius);
    266:           tao->trust = PetscMin(tao->trust, tr->max_radius);

    268:           PetscCall(KSPCGSetRadius(tao->ksp, tao->trust));
    269:           PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection));
    270:           PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    271:           tao->ksp_its += its;
    272:           tao->ksp_tot_its += its;
    273:           PetscCall(KSPCGGetNormD(tao->ksp, &norm_d));

    275:           PetscCheck(norm_d != 0.0, PetscObjectComm((PetscObject)tao), PETSC_ERR_PLIB, "Initial direction zero");
    276:         }
    277:       }
    278:       PetscCall(VecScale(tao->stepdirection, -1.0));
    279:       PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason));
    280:       if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) {
    281:         /* Preconditioner is numerically indefinite; reset the
    282:            approximate if using BFGS preconditioning. */
    283:         PetscCall(MatLMVMReset(tr->M, PETSC_FALSE));
    284:         PetscCall(MatLMVMUpdate(tr->M, tao->solution, tao->gradient));
    285:       }

    287:       if (NTR_UPDATE_REDUCTION == tr->update_type) {
    288:         /* Get predicted reduction */
    289:         PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
    290:         if (prered >= 0.0) {
    291:           /* The predicted reduction has the wrong sign.  This cannot
    292:              happen in infinite precision arithmetic.  Step should
    293:              be rejected! */
    294:           tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
    295:         } else {
    296:           /* Compute trial step and function value */
    297:           PetscCall(VecCopy(tao->solution, tr->W));
    298:           PetscCall(VecAXPY(tr->W, 1.0, tao->stepdirection));
    299:           PetscCall(TaoComputeObjective(tao, tr->W, &ftrial));

    301:           if (PetscIsInfOrNanReal(ftrial)) {
    302:             tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
    303:           } else {
    304:             /* Compute and actual reduction */
    305:             actred = f - ftrial;
    306:             prered = -prered;
    307:             if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
    308:               kappa = 1.0;
    309:             } else {
    310:               kappa = actred / prered;
    311:             }

    313:             /* Accept or reject the step and update radius */
    314:             if (kappa < tr->eta1) {
    315:               /* Reject the step */
    316:               tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d);
    317:             } else {
    318:               /* Accept the step */
    319:               if (kappa < tr->eta2) {
    320:                 /* Marginal bad step */
    321:                 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d);
    322:               } else if (kappa < tr->eta3) {
    323:                 /* Reasonable step */
    324:                 tao->trust = tr->alpha3 * tao->trust;
    325:               } else if (kappa < tr->eta4) {
    326:                 /* Good step */
    327:                 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust);
    328:               } else {
    329:                 /* Very good step */
    330:                 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust);
    331:               }
    332:               break;
    333:             }
    334:           }
    335:         }
    336:       } else {
    337:         /* Get predicted reduction */
    338:         PetscCall(KSPCGGetObjFcn(tao->ksp, &prered));
    339:         if (prered >= 0.0) {
    340:           /* The predicted reduction has the wrong sign.  This cannot
    341:              happen in infinite precision arithmetic.  Step should
    342:              be rejected! */
    343:           tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
    344:         } else {
    345:           PetscCall(VecCopy(tao->solution, tr->W));
    346:           PetscCall(VecAXPY(tr->W, 1.0, tao->stepdirection));
    347:           PetscCall(TaoComputeObjective(tao, tr->W, &ftrial));
    348:           if (PetscIsInfOrNanReal(ftrial)) {
    349:             tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
    350:           } else {
    351:             PetscCall(VecDot(tao->gradient, tao->stepdirection, &beta));
    352:             actred = f - ftrial;
    353:             prered = -prered;
    354:             if ((PetscAbsScalar(actred) <= tr->epsilon) && (PetscAbsScalar(prered) <= tr->epsilon)) {
    355:               kappa = 1.0;
    356:             } else {
    357:               kappa = actred / prered;
    358:             }

    360:             tau_1   = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred);
    361:             tau_2   = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred);
    362:             tau_min = PetscMin(tau_1, tau_2);
    363:             tau_max = PetscMax(tau_1, tau_2);

    365:             if (kappa >= 1.0 - tr->mu1) {
    366:               /* Great agreement; accept step and update radius */
    367:               if (tau_max < 1.0) {
    368:                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
    369:               } else if (tau_max > tr->gamma4) {
    370:                 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d);
    371:               } else {
    372:                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
    373:               }
    374:               break;
    375:             } else if (kappa >= 1.0 - tr->mu2) {
    376:               /* Good agreement */

    378:               if (tau_max < tr->gamma2) {
    379:                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
    380:               } else if (tau_max > tr->gamma3) {
    381:                 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d);
    382:               } else if (tau_max < 1.0) {
    383:                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
    384:               } else {
    385:                 tao->trust = PetscMax(tao->trust, tau_max * norm_d);
    386:               }
    387:               break;
    388:             } else {
    389:               /* Not good agreement */
    390:               if (tau_min > 1.0) {
    391:                 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d);
    392:               } else if (tau_max < tr->gamma1) {
    393:                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
    394:               } else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) {
    395:                 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d);
    396:               } else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) {
    397:                 tao->trust = tau_1 * PetscMin(tao->trust, norm_d);
    398:               } else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) {
    399:                 tao->trust = tau_2 * PetscMin(tao->trust, norm_d);
    400:               } else {
    401:                 tao->trust = tau_max * PetscMin(tao->trust, norm_d);
    402:               }
    403:             }
    404:           }
    405:         }
    406:       }

    408:       /* The step computed was not good and the radius was decreased.
    409:          Monitor the radius to terminate. */
    410:       PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
    411:       PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust));
    412:       PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    413:     }

    415:     /* The radius may have been increased; modify if it is too large */
    416:     tao->trust = PetscMin(tao->trust, tr->max_radius);

    418:     if (tao->reason == TAO_CONTINUE_ITERATING) {
    419:       PetscCall(VecCopy(tr->W, tao->solution));
    420:       f = ftrial;
    421:       PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient));
    422:       PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
    423:       PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
    424:       needH = 1;
    425:       PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
    426:       PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust));
    427:       PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    428:     }
    429:   }
    430:   PetscFunctionReturn(PETSC_SUCCESS);
    431: }

    433: /*------------------------------------------------------------*/
    434: static PetscErrorCode TaoSetUp_NTR(Tao tao)
    435: {
    436:   TAO_NTR *tr = (TAO_NTR *)tao->data;

    438:   PetscFunctionBegin;
    439:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
    440:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
    441:   if (!tr->W) PetscCall(VecDuplicate(tao->solution, &tr->W));

    443:   tr->bfgs_pre = NULL;
    444:   tr->M        = NULL;
    445:   PetscFunctionReturn(PETSC_SUCCESS);
    446: }

    448: /*------------------------------------------------------------*/
    449: static PetscErrorCode TaoDestroy_NTR(Tao tao)
    450: {
    451:   TAO_NTR *tr = (TAO_NTR *)tao->data;

    453:   PetscFunctionBegin;
    454:   if (tao->setupcalled) PetscCall(VecDestroy(&tr->W));
    455:   PetscCall(KSPDestroy(&tao->ksp));
    456:   PetscCall(PetscFree(tao->data));
    457:   PetscFunctionReturn(PETSC_SUCCESS);
    458: }

    460: /*------------------------------------------------------------*/
    461: static PetscErrorCode TaoSetFromOptions_NTR(Tao tao, PetscOptionItems PetscOptionsObject)
    462: {
    463:   TAO_NTR *tr = (TAO_NTR *)tao->data;

    465:   PetscFunctionBegin;
    466:   PetscOptionsHeadBegin(PetscOptionsObject, "Newton trust region method for unconstrained optimization");
    467:   PetscCall(PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type, NULL));
    468:   PetscCall(PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type, NULL));
    469:   PetscCall(PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1, NULL));
    470:   PetscCall(PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2, NULL));
    471:   PetscCall(PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3, NULL));
    472:   PetscCall(PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4, NULL));
    473:   PetscCall(PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1, NULL));
    474:   PetscCall(PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2, NULL));
    475:   PetscCall(PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3, NULL));
    476:   PetscCall(PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4, NULL));
    477:   PetscCall(PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5, NULL));
    478:   PetscCall(PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1, NULL));
    479:   PetscCall(PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2, NULL));
    480:   PetscCall(PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1, NULL));
    481:   PetscCall(PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2, NULL));
    482:   PetscCall(PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3, NULL));
    483:   PetscCall(PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4, NULL));
    484:   PetscCall(PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta, NULL));
    485:   PetscCall(PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i, NULL));
    486:   PetscCall(PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i, NULL));
    487:   PetscCall(PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i, NULL));
    488:   PetscCall(PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i, NULL));
    489:   PetscCall(PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i, NULL));
    490:   PetscCall(PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i, NULL));
    491:   PetscCall(PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i, NULL));
    492:   PetscCall(PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius, NULL));
    493:   PetscCall(PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius, NULL));
    494:   PetscCall(PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon, NULL));
    495:   PetscOptionsHeadEnd();
    496:   PetscCall(KSPSetFromOptions(tao->ksp));
    497:   PetscFunctionReturn(PETSC_SUCCESS);
    498: }

    500: /*------------------------------------------------------------*/
    501: /*MC
    502:   TAONTR - Newton's method with trust region for unconstrained minimization.
    503:   At each iteration, the Newton trust region method solves the system.
    504:   NTR expects a KSP solver with a trust region radius.
    505:             min_d  .5 dT Hk d + gkT d,  s.t.   ||d|| < Delta_k

    507:   Options Database Keys:
    508: + -tao_ntr_init_type - "constant","direction","interpolation"
    509: . -tao_ntr_update_type - "reduction","interpolation"
    510: . -tao_ntr_min_radius - lower bound on trust region radius
    511: . -tao_ntr_max_radius - upper bound on trust region radius
    512: . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction
    513: . -tao_ntr_mu1_i - mu1 interpolation init factor
    514: . -tao_ntr_mu2_i - mu2 interpolation init factor
    515: . -tao_ntr_gamma1_i - gamma1 interpolation init factor
    516: . -tao_ntr_gamma2_i - gamma2 interpolation init factor
    517: . -tao_ntr_gamma3_i - gamma3 interpolation init factor
    518: . -tao_ntr_gamma4_i - gamma4 interpolation init factor
    519: . -tao_ntr_theta_i - theta1 interpolation init factor
    520: . -tao_ntr_eta1 - eta1 reduction update factor
    521: . -tao_ntr_eta2 - eta2 reduction update factor
    522: . -tao_ntr_eta3 - eta3 reduction update factor
    523: . -tao_ntr_eta4 - eta4 reduction update factor
    524: . -tao_ntr_alpha1 - alpha1 reduction update factor
    525: . -tao_ntr_alpha2 - alpha2 reduction update factor
    526: . -tao_ntr_alpha3 - alpha3 reduction update factor
    527: . -tao_ntr_alpha4 - alpha4 reduction update factor
    528: . -tao_ntr_alpha4 - alpha4 reduction update factor
    529: . -tao_ntr_mu1 - mu1 interpolation update
    530: . -tao_ntr_mu2 - mu2 interpolation update
    531: . -tao_ntr_gamma1 - gamma1 interpolcation update
    532: . -tao_ntr_gamma2 - gamma2 interpolcation update
    533: . -tao_ntr_gamma3 - gamma3 interpolcation update
    534: . -tao_ntr_gamma4 - gamma4 interpolation update
    535: - -tao_ntr_theta - theta interpolation update

    537:   Level: beginner
    538: M*/

    540: PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao)
    541: {
    542:   TAO_NTR *tr;

    544:   PetscFunctionBegin;
    545:   PetscCall(PetscNew(&tr));

    547:   tao->ops->setup          = TaoSetUp_NTR;
    548:   tao->ops->solve          = TaoSolve_NTR;
    549:   tao->ops->setfromoptions = TaoSetFromOptions_NTR;
    550:   tao->ops->destroy        = TaoDestroy_NTR;

    552:   /* Override default settings (unless already changed) */
    553:   PetscCall(TaoParametersInitialize(tao));
    554:   PetscObjectParameterSetDefault(tao, max_it, 50);
    555:   PetscObjectParameterSetDefault(tao, trust0, 100.0);
    556:   tao->data = (void *)tr;

    558:   /*  Standard trust region update parameters */
    559:   tr->eta1 = 1.0e-4;
    560:   tr->eta2 = 0.25;
    561:   tr->eta3 = 0.50;
    562:   tr->eta4 = 0.90;

    564:   tr->alpha1 = 0.25;
    565:   tr->alpha2 = 0.50;
    566:   tr->alpha3 = 1.00;
    567:   tr->alpha4 = 2.00;
    568:   tr->alpha5 = 4.00;

    570:   /*  Interpolation trust region update parameters */
    571:   tr->mu1 = 0.10;
    572:   tr->mu2 = 0.50;

    574:   tr->gamma1 = 0.25;
    575:   tr->gamma2 = 0.50;
    576:   tr->gamma3 = 2.00;
    577:   tr->gamma4 = 4.00;

    579:   tr->theta = 0.05;

    581:   /*  Interpolation parameters for initialization */
    582:   tr->mu1_i = 0.35;
    583:   tr->mu2_i = 0.50;

    585:   tr->gamma1_i = 0.0625;
    586:   tr->gamma2_i = 0.50;
    587:   tr->gamma3_i = 2.00;
    588:   tr->gamma4_i = 5.00;

    590:   tr->theta_i = 0.25;

    592:   tr->min_radius = 1.0e-10;
    593:   tr->max_radius = 1.0e10;
    594:   tr->epsilon    = 1.0e-6;

    596:   tr->init_type   = NTR_INIT_INTERPOLATION;
    597:   tr->update_type = NTR_UPDATE_REDUCTION;

    599:   /* Set linear solver to default for trust region */
    600:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
    601:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
    602:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
    603:   PetscCall(KSPAppendOptionsPrefix(tao->ksp, "tao_ntr_"));
    604:   PetscCall(KSPSetType(tao->ksp, KSPSTCG));
    605:   PetscFunctionReturn(PETSC_SUCCESS);
    606: }