Actual source code: blmvm.c

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

      5: /*------------------------------------------------------------*/
      6: static PetscErrorCode TaoSolve_BLMVM(Tao tao)
      7: {
      8:   TAO_BLMVM                   *blmP      = (TAO_BLMVM *)tao->data;
      9:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
     10:   PetscReal                    f, fold, gdx, gnorm, gnorm2;
     11:   PetscReal                    stepsize = 1.0, delta;

     13:   PetscFunctionBegin;
     14:   /*  Project initial point onto bounds */
     15:   PetscCall(TaoComputeVariableBounds(tao));
     16:   PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
     17:   PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));

     19:   /* Check convergence criteria */
     20:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, blmP->unprojected_gradient));
     21:   PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));

     23:   PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
     24:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

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

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

     39:   /* Have not converged; continue with Newton method */
     40:   while (tao->reason == TAO_CONTINUE_ITERATING) {
     41:     /* Call general purpose update function */
     42:     if (tao->ops->update) {
     43:       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
     44:       PetscCall(TaoComputeObjective(tao, tao->solution, &f));
     45:     }
     46:     /* Compute direction */
     47:     gnorm2 = gnorm * gnorm;
     48:     if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
     49:     if (f == 0.0) {
     50:       delta = 2.0 / gnorm2;
     51:     } else {
     52:       delta = 2.0 * PetscAbsScalar(f) / gnorm2;
     53:     }
     54:     PetscCall(MatLMVMSymBroydenSetDelta(blmP->M, delta));
     55:     PetscCall(MatLMVMUpdate(blmP->M, tao->solution, tao->gradient));
     56:     PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
     57:     PetscCall(VecBoundGradientProjection(tao->stepdirection, tao->solution, tao->XL, tao->XU, tao->gradient));

     59:     /* Check for success (descent direction) */
     60:     PetscCall(VecDot(blmP->unprojected_gradient, tao->gradient, &gdx));
     61:     if (gdx <= 0) {
     62:       /* Step is not descent or solve was not successful
     63:          Use steepest descent direction (scaled) */
     64:       ++blmP->grad;

     66:       PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
     67:       PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient));
     68:       PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
     69:     }
     70:     PetscCall(VecScale(tao->stepdirection, -1.0));

     72:     /* Perform the linesearch */
     73:     fold = f;
     74:     PetscCall(VecCopy(tao->solution, blmP->Xold));
     75:     PetscCall(VecCopy(blmP->unprojected_gradient, blmP->Gold));
     76:     PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
     77:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status));
     78:     PetscCall(TaoAddLineSearchCounts(tao));

     80:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
     81:       /* Linesearch failed
     82:          Reset factors and use scaled (projected) gradient step */
     83:       ++blmP->reset;

     85:       f = fold;
     86:       PetscCall(VecCopy(blmP->Xold, tao->solution));
     87:       PetscCall(VecCopy(blmP->Gold, blmP->unprojected_gradient));

     89:       PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
     90:       PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient));
     91:       PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
     92:       PetscCall(VecScale(tao->stepdirection, -1.0));

     94:       /* This may be incorrect; linesearch has values for stepmax and stepmin
     95:          that should be reset. */
     96:       PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
     97:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status));
     98:       PetscCall(TaoAddLineSearchCounts(tao));

    100:       if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
    101:         tao->reason = TAO_DIVERGED_LS_FAILURE;
    102:         break;
    103:       }
    104:     }

    106:     /* Check for converged */
    107:     PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
    108:     PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
    109:     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
    110:     tao->niter++;
    111:     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
    112:     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize));
    113:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    114:   }
    115:   PetscFunctionReturn(PETSC_SUCCESS);
    116: }

    118: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
    119: {
    120:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;

    122:   PetscFunctionBegin;
    123:   /* Existence of tao->solution checked in TaoSetup() */
    124:   PetscCall(VecDuplicate(tao->solution, &blmP->Xold));
    125:   PetscCall(VecDuplicate(tao->solution, &blmP->Gold));
    126:   PetscCall(VecDuplicate(tao->solution, &blmP->unprojected_gradient));
    127:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
    128:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
    129:   /* Allocate matrix for the limited memory approximation */
    130:   PetscCall(MatLMVMAllocate(blmP->M, tao->solution, blmP->unprojected_gradient));

    132:   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
    133:   if (blmP->H0) PetscCall(MatLMVMSetJ0(blmP->M, blmP->H0));
    134:   PetscFunctionReturn(PETSC_SUCCESS);
    135: }

    137: /* ---------------------------------------------------------- */
    138: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
    139: {
    140:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;

    142:   PetscFunctionBegin;
    143:   if (tao->setupcalled) {
    144:     PetscCall(VecDestroy(&blmP->unprojected_gradient));
    145:     PetscCall(VecDestroy(&blmP->Xold));
    146:     PetscCall(VecDestroy(&blmP->Gold));
    147:   }
    148:   PetscCall(MatDestroy(&blmP->M));
    149:   if (blmP->H0) PetscCall(PetscObjectDereference((PetscObject)blmP->H0));
    150:   PetscCall(PetscFree(tao->data));
    151:   PetscFunctionReturn(PETSC_SUCCESS);
    152: }

    154: /*------------------------------------------------------------*/
    155: static PetscErrorCode TaoSetFromOptions_BLMVM(Tao tao, PetscOptionItems *PetscOptionsObject)
    156: {
    157:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
    158:   PetscBool  is_spd, is_set;

    160:   PetscFunctionBegin;
    161:   PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for bound constrained optimization");
    162:   PetscCall(PetscOptionsBool("-tao_blmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", blmP->recycle, &blmP->recycle, NULL));
    163:   PetscOptionsHeadEnd();
    164:   PetscCall(MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix));
    165:   PetscCall(MatAppendOptionsPrefix(blmP->M, "tao_blmvm_"));
    166:   PetscCall(MatSetFromOptions(blmP->M));
    167:   PetscCall(MatIsSPDKnown(blmP->M, &is_set, &is_spd));
    168:   PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite");
    169:   PetscFunctionReturn(PETSC_SUCCESS);
    170: }

    172: /*------------------------------------------------------------*/
    173: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
    174: {
    175:   TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
    176:   PetscBool  isascii;

    178:   PetscFunctionBegin;
    179:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
    180:   if (isascii) {
    181:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad));
    182:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
    183:     PetscCall(MatView(lmP->M, viewer));
    184:     PetscCall(PetscViewerPopFormat(viewer));
    185:   }
    186:   PetscFunctionReturn(PETSC_SUCCESS);
    187: }

    189: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
    190: {
    191:   TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;

    193:   PetscFunctionBegin;
    197:   PetscCheck(tao->gradient && blm->unprojected_gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist.");

    199:   PetscCall(VecCopy(tao->gradient, DXL));
    200:   PetscCall(VecAXPY(DXL, -1.0, blm->unprojected_gradient));
    201:   PetscCall(VecSet(DXU, 0.0));
    202:   PetscCall(VecPointwiseMax(DXL, DXL, DXU));

    204:   PetscCall(VecCopy(blm->unprojected_gradient, DXU));
    205:   PetscCall(VecAXPY(DXU, -1.0, tao->gradient));
    206:   PetscCall(VecAXPY(DXU, 1.0, DXL));
    207:   PetscFunctionReturn(PETSC_SUCCESS);
    208: }

    210: /* ---------------------------------------------------------- */
    211: /*MC
    212:   TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
    213:          for nonlinear minimization with bound constraints. It is an extension
    214:          of `TAOLMVM`

    216:   Options Database Key:
    217: .     -tao_lmm_recycle - enable recycling of LMVM information between subsequent `TaoSolve()` calls

    219:   Level: beginner

    221: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
    222: M*/
    223: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
    224: {
    225:   TAO_BLMVM  *blmP;
    226:   const char *morethuente_type = TAOLINESEARCHMT;

    228:   PetscFunctionBegin;
    229:   tao->ops->setup          = TaoSetup_BLMVM;
    230:   tao->ops->solve          = TaoSolve_BLMVM;
    231:   tao->ops->view           = TaoView_BLMVM;
    232:   tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
    233:   tao->ops->destroy        = TaoDestroy_BLMVM;
    234:   tao->ops->computedual    = TaoComputeDual_BLMVM;

    236:   PetscCall(PetscNew(&blmP));
    237:   blmP->H0      = NULL;
    238:   blmP->recycle = PETSC_FALSE;
    239:   tao->data     = (void *)blmP;

    241:   /* Override default settings (unless already changed) */
    242:   PetscCall(TaoParametersInitialize(tao));
    243:   PetscObjectParameterSetDefault(tao, max_it, 2000);
    244:   PetscObjectParameterSetDefault(tao, max_funcs, 4000);

    246:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
    247:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
    248:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
    249:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));

    251:   PetscCall(KSPInitializePackage());
    252:   PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M));
    253:   PetscCall(MatSetType(blmP->M, MATLMVMBFGS));
    254:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1));
    255:   PetscFunctionReturn(PETSC_SUCCESS);
    256: }

    258: /*@
    259:   TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent `TaoSolve()` calls.

    261:   Input Parameters:
    262: + tao - the `Tao` solver context
    263: - flg - Boolean flag for recycling (`PETSC_TRUE` or `PETSC_FALSE`)

    265:   Level: intermediate

    267: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`
    268: @*/
    269: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
    270: {
    271:   TAO_LMVM  *lmP;
    272:   TAO_BLMVM *blmP;
    273:   PetscBool  is_lmvm, is_blmvm;

    275:   PetscFunctionBegin;
    276:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
    277:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
    278:   if (is_lmvm) {
    279:     lmP          = (TAO_LMVM *)tao->data;
    280:     lmP->recycle = flg;
    281:   } else if (is_blmvm) {
    282:     blmP          = (TAO_BLMVM *)tao->data;
    283:     blmP->recycle = flg;
    284:   }
    285:   PetscFunctionReturn(PETSC_SUCCESS);
    286: }

    288: /*@
    289:   TaoLMVMSetH0 - Set the initial Hessian for the QN approximation

    291:   Input Parameters:
    292: + tao - the `Tao` solver context
    293: - H0  - `Mat` object for the initial Hessian

    295:   Level: advanced

    297: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
    298: @*/
    299: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
    300: {
    301:   TAO_LMVM  *lmP;
    302:   TAO_BLMVM *blmP;
    303:   PetscBool  is_lmvm, is_blmvm;

    305:   PetscFunctionBegin;
    306:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
    307:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
    308:   if (is_lmvm) {
    309:     lmP = (TAO_LMVM *)tao->data;
    310:     PetscCall(PetscObjectReference((PetscObject)H0));
    311:     lmP->H0 = H0;
    312:   } else if (is_blmvm) {
    313:     blmP = (TAO_BLMVM *)tao->data;
    314:     PetscCall(PetscObjectReference((PetscObject)H0));
    315:     blmP->H0 = H0;
    316:   }
    317:   PetscFunctionReturn(PETSC_SUCCESS);
    318: }

    320: /*@
    321:   TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian

    323:   Input Parameter:
    324: . tao - the `Tao` solver context

    326:   Output Parameter:
    327: . H0 - `Mat` object for the initial Hessian

    329:   Level: advanced

    331: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
    332: @*/
    333: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
    334: {
    335:   TAO_LMVM  *lmP;
    336:   TAO_BLMVM *blmP;
    337:   PetscBool  is_lmvm, is_blmvm;
    338:   Mat        M;

    340:   PetscFunctionBegin;
    341:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
    342:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
    343:   if (is_lmvm) {
    344:     lmP = (TAO_LMVM *)tao->data;
    345:     M   = lmP->M;
    346:   } else if (is_blmvm) {
    347:     blmP = (TAO_BLMVM *)tao->data;
    348:     M    = blmP->M;
    349:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
    350:   PetscCall(MatLMVMGetJ0(M, H0));
    351:   PetscFunctionReturn(PETSC_SUCCESS);
    352: }

    354: /*@
    355:   TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian

    357:   Input Parameter:
    358: . tao - the `Tao` solver context

    360:   Output Parameter:
    361: . ksp - `KSP` solver context for the initial Hessian

    363:   Level: advanced

    365: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`
    366: @*/
    367: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
    368: {
    369:   TAO_LMVM  *lmP;
    370:   TAO_BLMVM *blmP;
    371:   PetscBool  is_lmvm, is_blmvm;
    372:   Mat        M;

    374:   PetscFunctionBegin;
    375:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
    376:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
    377:   if (is_lmvm) {
    378:     lmP = (TAO_LMVM *)tao->data;
    379:     M   = lmP->M;
    380:   } else if (is_blmvm) {
    381:     blmP = (TAO_BLMVM *)tao->data;
    382:     M    = blmP->M;
    383:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
    384:   PetscCall(MatLMVMGetJ0KSP(M, ksp));
    385:   PetscFunctionReturn(PETSC_SUCCESS);
    386: }