Actual source code: taocg.c

      1: #include <petsctaolinesearch.h>
      2: #include <../src/tao/unconstrained/impls/cg/taocg.h>

      4: #define CG_FletcherReeves   0
      5: #define CG_PolakRibiere     1
      6: #define CG_PolakRibierePlus 2
      7: #define CG_HestenesStiefel  3
      8: #define CG_DaiYuan          4
      9: #define CG_Types            5

     11: static const char *CG_Table[64] = {"fr", "pr", "prp", "hs", "dy"};

     13: static PetscErrorCode TaoSolve_CG(Tao tao)
     14: {
     15:   TAO_CG                      *cgP       = (TAO_CG *)tao->data;
     16:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
     17:   PetscReal                    step      = 1.0, f, gnorm, gnorm2, delta, gd, ginner, beta;
     18:   PetscReal                    gd_old, gnorm2_old, f_old;

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

     23:   /*  Check convergence criteria */
     24:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
     25:   PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
     26:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");

     28:   tao->reason = TAO_CONTINUE_ITERATING;
     29:   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
     30:   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
     31:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
     32:   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);

     34:   /*  Set initial direction to -gradient */
     35:   PetscCall(VecCopy(tao->gradient, tao->stepdirection));
     36:   PetscCall(VecScale(tao->stepdirection, -1.0));
     37:   gnorm2 = gnorm * gnorm;

     39:   /*  Set initial scaling for the function */
     40:   if (f != 0.0) {
     41:     delta = 2.0 * PetscAbsScalar(f) / gnorm2;
     42:     delta = PetscMax(delta, cgP->delta_min);
     43:     delta = PetscMin(delta, cgP->delta_max);
     44:   } else {
     45:     delta = 2.0 / gnorm2;
     46:     delta = PetscMax(delta, cgP->delta_min);
     47:     delta = PetscMin(delta, cgP->delta_max);
     48:   }
     49:   /*  Set counter for gradient and reset steps */
     50:   cgP->ngradsteps  = 0;
     51:   cgP->nresetsteps = 0;

     53:   while (1) {
     54:     /* Call general purpose update function */
     55:     if (tao->ops->update) {
     56:       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
     57:       PetscCall(TaoComputeObjective(tao, tao->solution, &f));
     58:     }
     59:     /*  Save the current gradient information */
     60:     f_old      = f;
     61:     gnorm2_old = gnorm2;
     62:     PetscCall(VecCopy(tao->solution, cgP->X_old));
     63:     PetscCall(VecCopy(tao->gradient, cgP->G_old));
     64:     PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
     65:     if ((gd >= 0) || PetscIsInfOrNanReal(gd)) {
     66:       ++cgP->ngradsteps;
     67:       if (f != 0.0) {
     68:         delta = 2.0 * PetscAbsScalar(f) / gnorm2;
     69:         delta = PetscMax(delta, cgP->delta_min);
     70:         delta = PetscMin(delta, cgP->delta_max);
     71:       } else {
     72:         delta = 2.0 / gnorm2;
     73:         delta = PetscMax(delta, cgP->delta_min);
     74:         delta = PetscMin(delta, cgP->delta_max);
     75:       }

     77:       PetscCall(VecCopy(tao->gradient, tao->stepdirection));
     78:       PetscCall(VecScale(tao->stepdirection, -1.0));
     79:     }

     81:     /*  Search direction for improving point */
     82:     PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta));
     83:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status));
     84:     PetscCall(TaoAddLineSearchCounts(tao));
     85:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
     86:       /*  Linesearch failed */
     87:       /*  Reset factors and use scaled gradient step */
     88:       ++cgP->nresetsteps;
     89:       f      = f_old;
     90:       gnorm2 = gnorm2_old;
     91:       PetscCall(VecCopy(cgP->X_old, tao->solution));
     92:       PetscCall(VecCopy(cgP->G_old, tao->gradient));

     94:       if (f != 0.0) {
     95:         delta = 2.0 * PetscAbsScalar(f) / gnorm2;
     96:         delta = PetscMax(delta, cgP->delta_min);
     97:         delta = PetscMin(delta, cgP->delta_max);
     98:       } else {
     99:         delta = 2.0 / gnorm2;
    100:         delta = PetscMax(delta, cgP->delta_min);
    101:         delta = PetscMin(delta, cgP->delta_max);
    102:       }

    104:       PetscCall(VecCopy(tao->gradient, tao->stepdirection));
    105:       PetscCall(VecScale(tao->stepdirection, -1.0));

    107:       PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta));
    108:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status));
    109:       PetscCall(TaoAddLineSearchCounts(tao));

    111:       if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
    112:         /*  Linesearch failed again */
    113:         /*  switch to unscaled gradient */
    114:         f = f_old;
    115:         PetscCall(VecCopy(cgP->X_old, tao->solution));
    116:         PetscCall(VecCopy(cgP->G_old, tao->gradient));
    117:         delta = 1.0;
    118:         PetscCall(VecCopy(tao->gradient, tao->stepdirection));
    119:         PetscCall(VecScale(tao->stepdirection, -1.0));

    121:         PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta));
    122:         PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status));
    123:         PetscCall(TaoAddLineSearchCounts(tao));
    124:         if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
    125:           /*  Line search failed for last time -- give up */
    126:           f = f_old;
    127:           PetscCall(VecCopy(cgP->X_old, tao->solution));
    128:           PetscCall(VecCopy(cgP->G_old, tao->gradient));
    129:           step        = 0.0;
    130:           tao->reason = TAO_DIVERGED_LS_FAILURE;
    131:         }
    132:       }
    133:     }

    135:     /*  Check for bad value */
    136:     PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
    137:     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User-provided compute function generated infinity or NaN");

    139:     /*  Check for termination */
    140:     gnorm2 = gnorm * gnorm;
    141:     tao->niter++;
    142:     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
    143:     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
    144:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    145:     if (tao->reason != TAO_CONTINUE_ITERATING) break;

    147:     /*  Check for restart condition */
    148:     PetscCall(VecDot(tao->gradient, cgP->G_old, &ginner));
    149:     if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) {
    150:       /*  Gradients far from orthogonal; use steepest descent direction */
    151:       beta = 0.0;
    152:     } else {
    153:       /*  Gradients close to orthogonal; use conjugate gradient formula */
    154:       switch (cgP->cg_type) {
    155:       case CG_FletcherReeves:
    156:         beta = gnorm2 / gnorm2_old;
    157:         break;

    159:       case CG_PolakRibiere:
    160:         beta = (gnorm2 - ginner) / gnorm2_old;
    161:         break;

    163:       case CG_PolakRibierePlus:
    164:         beta = PetscMax((gnorm2 - ginner) / gnorm2_old, 0.0);
    165:         break;

    167:       case CG_HestenesStiefel:
    168:         PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
    169:         PetscCall(VecDot(cgP->G_old, tao->stepdirection, &gd_old));
    170:         beta = (gnorm2 - ginner) / (gd - gd_old);
    171:         break;

    173:       case CG_DaiYuan:
    174:         PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
    175:         PetscCall(VecDot(cgP->G_old, tao->stepdirection, &gd_old));
    176:         beta = gnorm2 / (gd - gd_old);
    177:         break;

    179:       default:
    180:         beta = 0.0;
    181:         break;
    182:       }
    183:     }

    185:     /*  Compute the direction d=-g + beta*d */
    186:     PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient));

    188:     /*  update initial steplength choice */
    189:     delta = 1.0;
    190:     delta = PetscMax(delta, cgP->delta_min);
    191:     delta = PetscMin(delta, cgP->delta_max);
    192:   }
    193:   PetscFunctionReturn(PETSC_SUCCESS);
    194: }

    196: static PetscErrorCode TaoSetUp_CG(Tao tao)
    197: {
    198:   TAO_CG *cgP = (TAO_CG *)tao->data;

    200:   PetscFunctionBegin;
    201:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
    202:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
    203:   if (!cgP->X_old) PetscCall(VecDuplicate(tao->solution, &cgP->X_old));
    204:   if (!cgP->G_old) PetscCall(VecDuplicate(tao->gradient, &cgP->G_old));
    205:   PetscFunctionReturn(PETSC_SUCCESS);
    206: }

    208: static PetscErrorCode TaoDestroy_CG(Tao tao)
    209: {
    210:   TAO_CG *cgP = (TAO_CG *)tao->data;

    212:   PetscFunctionBegin;
    213:   if (tao->setupcalled) {
    214:     PetscCall(VecDestroy(&cgP->X_old));
    215:     PetscCall(VecDestroy(&cgP->G_old));
    216:   }
    217:   PetscCall(TaoLineSearchDestroy(&tao->linesearch));
    218:   PetscCall(PetscFree(tao->data));
    219:   PetscFunctionReturn(PETSC_SUCCESS);
    220: }

    222: static PetscErrorCode TaoSetFromOptions_CG(Tao tao, PetscOptionItems PetscOptionsObject)
    223: {
    224:   TAO_CG *cgP = (TAO_CG *)tao->data;

    226:   PetscFunctionBegin;
    227:   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
    228:   PetscOptionsHeadBegin(PetscOptionsObject, "Nonlinear Conjugate Gradient method for unconstrained optimization");
    229:   PetscCall(PetscOptionsReal("-tao_cg_eta", "restart tolerance", "", cgP->eta, &cgP->eta, NULL));
    230:   PetscCall(PetscOptionsEList("-tao_cg_type", "cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, NULL));
    231:   PetscCall(PetscOptionsReal("-tao_cg_delta_min", "minimum delta value", "", cgP->delta_min, &cgP->delta_min, NULL));
    232:   PetscCall(PetscOptionsReal("-tao_cg_delta_max", "maximum delta value", "", cgP->delta_max, &cgP->delta_max, NULL));
    233:   PetscOptionsHeadEnd();
    234:   PetscFunctionReturn(PETSC_SUCCESS);
    235: }

    237: static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer)
    238: {
    239:   PetscBool isascii;
    240:   TAO_CG   *cgP = (TAO_CG *)tao->data;

    242:   PetscFunctionBegin;
    243:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    244:   if (isascii) {
    245:     PetscCall(PetscViewerASCIIPushTab(viewer));
    246:     PetscCall(PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type]));
    247:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", cgP->ngradsteps));
    248:     PetscCall(PetscViewerASCIIPrintf(viewer, "Reset steps: %" PetscInt_FMT "\n", cgP->nresetsteps));
    249:     PetscCall(PetscViewerASCIIPopTab(viewer));
    250:   }
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

    254: /*MC
    255:    TAOCG -   Nonlinear conjugate gradient method is an extension of the
    256:                nonlinear conjugate gradient solver for nonlinear optimization.

    258:    Options Database Keys:
    259: +      -tao_cg_eta r           - restart tolerance
    260: .      -tao_cg_type taocg_type - cg formula, one of `fr`, `pr`, `prp`, `hs`, `dy`
    261: .      -tao_cg_delta_min r     - minimum delta value
    262: -      -tao_cg_delta_max r     - maximum delta value

    264:   Level: beginner

    266:    Note:
    267:    The cg formulas are:
    268: +  `fr`  - Fletcher-Reeves
    269: .  `pr`  - Polak-Ribiere
    270: .  `prp` - Polak-Ribiere-Plus
    271: .  `hs`  - Hestenes-Steifel
    272: -  `dy`  - Dai-Yuan

    274: .seealso: `Tao`, `TAONTR`, `TAONLS`, `TaoType`, `TaoCreate()`
    275: M*/

    277: PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao)
    278: {
    279:   TAO_CG     *cgP;
    280:   const char *morethuente_type = TAOLINESEARCHMT;

    282:   PetscFunctionBegin;
    283:   tao->ops->setup          = TaoSetUp_CG;
    284:   tao->ops->solve          = TaoSolve_CG;
    285:   tao->ops->view           = TaoView_CG;
    286:   tao->ops->setfromoptions = TaoSetFromOptions_CG;
    287:   tao->ops->destroy        = TaoDestroy_CG;
    288:   tao->uses_gradient       = PETSC_TRUE;

    290:   /* Override default settings (unless already changed) */
    291:   PetscCall(TaoParametersInitialize(tao));
    292:   PetscObjectParameterSetDefault(tao, max_it, 2000);
    293:   PetscObjectParameterSetDefault(tao, max_funcs, 4000);

    295:   /*  Note: nondefault values should be used for nonlinear conjugate gradient  */
    296:   /*  method.  In particular, gtol should be less that 0.5; the value used in  */
    297:   /*  Nocedal and Wright is 0.10.  We use the default values for the  */
    298:   /*  linesearch because it seems to work better. */
    299:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
    300:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
    301:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
    302:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
    303:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));

    305:   PetscCall(PetscNew(&cgP));
    306:   tao->data      = (void *)cgP;
    307:   cgP->eta       = 0.1;
    308:   cgP->delta_min = 1e-7;
    309:   cgP->delta_max = 100;
    310:   cgP->cg_type   = CG_PolakRibierePlus;
    311:   PetscFunctionReturn(PETSC_SUCCESS);
    312: }