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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    169: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
    170: {
    171:   TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
    172:   PetscBool  isascii;

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

    185: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
    186: {
    187:   TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;

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

    195:   PetscCall(VecCopy(tao->gradient, DXL));
    196:   PetscCall(VecAXPY(DXL, -1.0, blm->unprojected_gradient));
    197:   PetscCall(VecSet(DXU, 0.0));
    198:   PetscCall(VecPointwiseMax(DXL, DXL, DXU));

    200:   PetscCall(VecCopy(blm->unprojected_gradient, DXU));
    201:   PetscCall(VecAXPY(DXU, -1.0, tao->gradient));
    202:   PetscCall(VecAXPY(DXU, 1.0, DXL));
    203:   PetscFunctionReturn(PETSC_SUCCESS);
    204: }

    206: /*MC
    207:   TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
    208:          for nonlinear minimization with bound constraints. It is an extension
    209:          of `TAOLMVM`

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

    214:   Level: beginner

    216: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
    217: M*/
    218: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
    219: {
    220:   TAO_BLMVM  *blmP;
    221:   const char *morethuente_type = TAOLINESEARCHMT;

    223:   PetscFunctionBegin;
    224:   tao->ops->setup          = TaoSetup_BLMVM;
    225:   tao->ops->solve          = TaoSolve_BLMVM;
    226:   tao->ops->view           = TaoView_BLMVM;
    227:   tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
    228:   tao->ops->destroy        = TaoDestroy_BLMVM;
    229:   tao->ops->computedual    = TaoComputeDual_BLMVM;
    230:   tao->uses_gradient       = PETSC_TRUE;

    232:   PetscCall(PetscNew(&blmP));
    233:   blmP->H0      = NULL;
    234:   blmP->recycle = PETSC_FALSE;
    235:   tao->data     = (void *)blmP;

    237:   /* Override default settings (unless already changed) */
    238:   PetscCall(TaoParametersInitialize(tao));
    239:   PetscObjectParameterSetDefault(tao, max_it, 2000);
    240:   PetscObjectParameterSetDefault(tao, max_funcs, 4000);

    242:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
    243:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
    244:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
    245:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));

    247:   PetscCall(KSPInitializePackage());
    248:   PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M));
    249:   PetscCall(MatSetType(blmP->M, MATLMVMBFGS));
    250:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1));
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

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

    257:   Input Parameters:
    258: + tao - the `Tao` solver context
    259: - flg - Boolean flag for recycling (`PETSC_TRUE` or `PETSC_FALSE`)

    261:   Level: intermediate

    263: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`
    264: @*/
    265: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
    266: {
    267:   TAO_LMVM  *lmP;
    268:   TAO_BLMVM *blmP;
    269:   PetscBool  is_lmvm, is_blmvm;

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

    284: /*@
    285:   TaoLMVMSetH0 - Set the initial Hessian for the QN approximation

    287:   Input Parameters:
    288: + tao - the `Tao` solver context
    289: - H0  - `Mat` object for the initial Hessian

    291:   Level: advanced

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

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

    316: /*@
    317:   TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian

    319:   Input Parameter:
    320: . tao - the `Tao` solver context

    322:   Output Parameter:
    323: . H0 - `Mat` object for the initial Hessian

    325:   Level: advanced

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

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

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

    353:   Input Parameter:
    354: . tao - the `Tao` solver context

    356:   Output Parameter:
    357: . ksp - `KSP` solver context for the initial Hessian

    359:   Level: advanced

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

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