Actual source code: lmvm.c

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

      4: #define LMVM_STEP_BFGS 0
      5: #define LMVM_STEP_GRAD 1

      7: static PetscErrorCode TaoSolve_LMVM(Tao tao)
      8: {
      9:   TAO_LMVM                    *lmP = (TAO_LMVM *)tao->data;
     10:   PetscReal                    f, fold, gdx, gnorm;
     11:   PetscReal                    step      = 1.0;
     12:   PetscInt                     stepType  = LMVM_STEP_GRAD, nupdates;
     13:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

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

     18:   /*  Check convergence criteria */
     19:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
     20:   PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));

     22:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

     24:   tao->reason = TAO_CONTINUE_ITERATING;
     25:   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
     26:   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
     27:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
     28:   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);

     30:   /*  Set counter for gradient/reset steps */
     31:   if (!lmP->recycle) {
     32:     lmP->bfgs = 0;
     33:     lmP->grad = 0;
     34:     PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
     35:   }

     37:   /*  Have not converged; continue with Newton method */
     38:   while (tao->reason == TAO_CONTINUE_ITERATING) {
     39:     /* Call general purpose update function */
     40:     if (tao->ops->update) {
     41:       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
     42:       PetscCall(TaoComputeObjective(tao, tao->solution, &f));
     43:     }

     45:     /*  Compute direction */
     46:     if (lmP->H0) {
     47:       PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0));
     48:       stepType = LMVM_STEP_BFGS;
     49:     }
     50:     PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
     51:     PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D));
     52:     PetscCall(MatLMVMGetUpdateCount(lmP->M, &nupdates));
     53:     if (nupdates > 0) stepType = LMVM_STEP_BFGS;

     55:     /*  Check for success (descent direction) */
     56:     PetscCall(VecDotRealPart(lmP->D, tao->gradient, &gdx));
     57:     if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
     58:       /* Step is not descent or direction produced not a number
     59:          We can assert bfgsUpdates > 1 in this case because
     60:          the first solve produces the scaled gradient direction,
     61:          which is guaranteed to be descent

     63:          Use steepest descent direction (scaled)
     64:       */

     66:       PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
     67:       PetscCall(MatLMVMClearJ0(lmP->M));
     68:       PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
     69:       PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D));

     71:       /* On a reset, the direction cannot be not a number; it is a
     72:          scaled gradient step.  No need to check for this condition. */
     73:       stepType = LMVM_STEP_GRAD;
     74:     }
     75:     PetscCall(VecScale(lmP->D, -1.0));

     77:     /*  Perform the linesearch */
     78:     fold = f;
     79:     PetscCall(VecCopy(tao->solution, lmP->Xold));
     80:     PetscCall(VecCopy(tao->gradient, lmP->Gold));

     82:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status));
     83:     PetscCall(TaoAddLineSearchCounts(tao));

     85:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) {
     86:       /*  Reset factors and use scaled gradient step */
     87:       f = fold;
     88:       PetscCall(VecCopy(lmP->Xold, tao->solution));
     89:       PetscCall(VecCopy(lmP->Gold, tao->gradient));

     91:       /*  Failed to obtain acceptable iterate with BFGS step */
     92:       /*  Attempt to use the scaled gradient direction */

     94:       PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
     95:       PetscCall(MatLMVMClearJ0(lmP->M));
     96:       PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
     97:       PetscCall(MatSolve(lmP->M, tao->solution, tao->gradient));

     99:       /* On a reset, the direction cannot be not a number; it is a
    100:           scaled gradient step.  No need to check for this condition. */
    101:       stepType = LMVM_STEP_GRAD;
    102:       PetscCall(VecScale(lmP->D, -1.0));

    104:       /*  Perform the linesearch */
    105:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status));
    106:       PetscCall(TaoAddLineSearchCounts(tao));
    107:     }

    109:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
    110:       /*  Failed to find an improving point */
    111:       f = fold;
    112:       PetscCall(VecCopy(lmP->Xold, tao->solution));
    113:       PetscCall(VecCopy(lmP->Gold, tao->gradient));
    114:       step        = 0.0;
    115:       tao->reason = TAO_DIVERGED_LS_FAILURE;
    116:     } else {
    117:       /* LS found valid step, so tally up step type */
    118:       switch (stepType) {
    119:       case LMVM_STEP_BFGS:
    120:         ++lmP->bfgs;
    121:         break;
    122:       case LMVM_STEP_GRAD:
    123:         ++lmP->grad;
    124:         break;
    125:       default:
    126:         break;
    127:       }
    128:       /*  Compute new gradient norm */
    129:       PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
    130:     }

    132:     /* Check convergence */
    133:     tao->niter++;
    134:     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
    135:     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
    136:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    137:   }
    138:   PetscFunctionReturn(PETSC_SUCCESS);
    139: }

    141: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
    142: {
    143:   TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
    144:   PetscInt  n, N;
    145:   PetscBool is_set, is_spd;

    147:   PetscFunctionBegin;
    148:   /* Existence of tao->solution checked in TaoSetUp() */
    149:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
    150:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
    151:   if (!lmP->D) PetscCall(VecDuplicate(tao->solution, &lmP->D));
    152:   if (!lmP->Xold) PetscCall(VecDuplicate(tao->solution, &lmP->Xold));
    153:   if (!lmP->Gold) PetscCall(VecDuplicate(tao->solution, &lmP->Gold));

    155:   /*  Create matrix for the limited memory approximation */
    156:   PetscCall(VecGetLocalSize(tao->solution, &n));
    157:   PetscCall(VecGetSize(tao->solution, &N));
    158:   PetscCall(MatSetSizes(lmP->M, n, n, N, N));
    159:   PetscCall(MatLMVMAllocate(lmP->M, tao->solution, tao->gradient));
    160:   PetscCall(MatIsSPDKnown(lmP->M, &is_set, &is_spd));
    161:   PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite.");

    163:   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
    164:   if (lmP->H0) PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0));
    165:   PetscFunctionReturn(PETSC_SUCCESS);
    166: }

    168: /* ---------------------------------------------------------- */
    169: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
    170: {
    171:   TAO_LMVM *lmP = (TAO_LMVM *)tao->data;

    173:   PetscFunctionBegin;
    174:   if (tao->setupcalled) {
    175:     PetscCall(VecDestroy(&lmP->Xold));
    176:     PetscCall(VecDestroy(&lmP->Gold));
    177:     PetscCall(VecDestroy(&lmP->D));
    178:   }
    179:   PetscCall(MatDestroy(&lmP->M));
    180:   if (lmP->H0) PetscCall(PetscObjectDereference((PetscObject)lmP->H0));
    181:   PetscCall(PetscFree(tao->data));
    182:   PetscFunctionReturn(PETSC_SUCCESS);
    183: }

    185: /*------------------------------------------------------------*/
    186: static PetscErrorCode TaoSetFromOptions_LMVM(Tao tao, PetscOptionItems PetscOptionsObject)
    187: {
    188:   TAO_LMVM *lm = (TAO_LMVM *)tao->data;

    190:   PetscFunctionBegin;
    191:   PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for unconstrained optimization");
    192:   PetscCall(PetscOptionsBool("-tao_lmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", lm->recycle, &lm->recycle, NULL));
    193:   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
    194:   PetscCall(MatSetFromOptions(lm->M));
    195:   PetscOptionsHeadEnd();
    196:   PetscFunctionReturn(PETSC_SUCCESS);
    197: }

    199: /*------------------------------------------------------------*/
    200: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
    201: {
    202:   TAO_LMVM *lm = (TAO_LMVM *)tao->data;
    203:   PetscBool isascii;
    204:   PetscInt  recycled_its;

    206:   PetscFunctionBegin;
    207:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    208:   if (isascii) {
    209:     PetscCall(PetscViewerASCIIPushTab(viewer));
    210:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lm->grad));
    211:     if (lm->recycle) {
    212:       PetscCall(PetscViewerASCIIPrintf(viewer, "Recycle: on\n"));
    213:       recycled_its = lm->bfgs + lm->grad;
    214:       PetscCall(PetscViewerASCIIPrintf(viewer, "Total recycled iterations: %" PetscInt_FMT "\n", recycled_its));
    215:     }
    216:     PetscCall(PetscViewerASCIIPrintf(viewer, "LMVM Matrix:\n"));
    217:     PetscCall(PetscViewerASCIIPushTab(viewer));
    218:     PetscCall(MatView(lm->M, viewer));
    219:     PetscCall(PetscViewerASCIIPopTab(viewer));
    220:     PetscCall(PetscViewerASCIIPopTab(viewer));
    221:   }
    222:   PetscFunctionReturn(PETSC_SUCCESS);
    223: }

    225: /* ---------------------------------------------------------- */

    227: /*MC
    228:   TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
    229:   optimization solver for unconstrained minimization. It solves
    230:   the Newton step
    231:           Hkdk = - gk

    233:   using an approximation Bk in place of Hk, where Bk is composed using
    234:   the BFGS update formula. A More-Thuente line search is then used
    235:   to computed the steplength in the dk direction

    237:   Options Database Keys:
    238: +   -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls
    239: -   -tao_lmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation

    241:   Level: beginner
    242: M*/

    244: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
    245: {
    246:   TAO_LMVM   *lmP;
    247:   const char *morethuente_type = TAOLINESEARCHMT;

    249:   PetscFunctionBegin;
    250:   tao->ops->setup          = TaoSetUp_LMVM;
    251:   tao->ops->solve          = TaoSolve_LMVM;
    252:   tao->ops->view           = TaoView_LMVM;
    253:   tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
    254:   tao->ops->destroy        = TaoDestroy_LMVM;

    256:   PetscCall(PetscNew(&lmP));
    257:   lmP->D       = NULL;
    258:   lmP->M       = NULL;
    259:   lmP->Xold    = NULL;
    260:   lmP->Gold    = NULL;
    261:   lmP->H0      = NULL;
    262:   lmP->recycle = PETSC_FALSE;

    264:   tao->data = (void *)lmP;
    265:   /* Override default settings (unless already changed) */
    266:   PetscCall(TaoParametersInitialize(tao));
    267:   PetscObjectParameterSetDefault(tao, max_it, 2000);
    268:   PetscObjectParameterSetDefault(tao, max_funcs, 4000);

    270:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
    271:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
    272:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
    273:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
    274:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));

    276:   PetscCall(KSPInitializePackage());
    277:   PetscCall(MatCreate(((PetscObject)tao)->comm, &lmP->M));
    278:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1));
    279:   PetscCall(MatSetType(lmP->M, MATLMVMBFGS));
    280:   PetscCall(MatSetOptionsPrefix(lmP->M, "tao_lmvm_"));
    281:   PetscFunctionReturn(PETSC_SUCCESS);
    282: }