Actual source code: gpcg.c

      1: #include <petscksp.h>
      2: #include <../src/tao/quadratic/impls/gpcg/gpcg.h>

      4: static PetscErrorCode GPCGGradProjections(Tao tao);
      5: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch, Vec, PetscReal *, Vec, void *);

      7: static PetscErrorCode TaoDestroy_GPCG(Tao tao)
      8: {
      9:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;

     11:   /* Free allocated memory in GPCG structure */
     12:   PetscFunctionBegin;
     13:   PetscCall(VecDestroy(&gpcg->B));
     14:   PetscCall(VecDestroy(&gpcg->Work));
     15:   PetscCall(VecDestroy(&gpcg->X_New));
     16:   PetscCall(VecDestroy(&gpcg->G_New));
     17:   PetscCall(VecDestroy(&gpcg->DXFree));
     18:   PetscCall(VecDestroy(&gpcg->R));
     19:   PetscCall(VecDestroy(&gpcg->PG));
     20:   PetscCall(MatDestroy(&gpcg->Hsub));
     21:   PetscCall(MatDestroy(&gpcg->Hsub_pre));
     22:   PetscCall(ISDestroy(&gpcg->Free_Local));
     23:   PetscCall(KSPDestroy(&tao->ksp));
     24:   PetscCall(PetscFree(tao->data));
     25:   PetscFunctionReturn(PETSC_SUCCESS);
     26: }

     28: static PetscErrorCode TaoSetFromOptions_GPCG(Tao tao, PetscOptionItems PetscOptionsObject)
     29: {
     30:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
     31:   PetscBool flg;

     33:   PetscFunctionBegin;
     34:   PetscOptionsHeadBegin(PetscOptionsObject, "Gradient Projection, Conjugate Gradient method for bound constrained optimization");
     35:   PetscCall(PetscOptionsInt("-tao_gpcg_maxpgits", "maximum number of gradient projections per GPCG iterate", NULL, gpcg->maxgpits, &gpcg->maxgpits, &flg));
     36:   PetscOptionsHeadEnd();
     37:   PetscCall(KSPSetFromOptions(tao->ksp));
     38:   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
     39:   PetscFunctionReturn(PETSC_SUCCESS);
     40: }

     42: static PetscErrorCode TaoView_GPCG(Tao tao, PetscViewer viewer)
     43: {
     44:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
     45:   PetscBool isascii;

     47:   PetscFunctionBegin;
     48:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
     49:   if (isascii) {
     50:     PetscCall(PetscViewerASCIIPrintf(viewer, "Total PG its: %" PetscInt_FMT ",", gpcg->total_gp_its));
     51:     PetscCall(PetscViewerASCIIPrintf(viewer, "PG tolerance: %g \n", (double)gpcg->pg_ftol));
     52:   }
     53:   PetscCall(TaoLineSearchView(tao->linesearch, viewer));
     54:   PetscFunctionReturn(PETSC_SUCCESS);
     55: }

     57: /* GPCGObjectiveAndGradient()
     58:    Compute f=0.5 * x'Hx + b'x + c
     59:            g=Hx + b
     60: */
     61: static PetscErrorCode GPCGObjectiveAndGradient(TaoLineSearch ls, Vec X, PetscReal *f, Vec G, void *tptr)
     62: {
     63:   Tao       tao  = (Tao)tptr;
     64:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;
     65:   PetscReal f1, f2;

     67:   PetscFunctionBegin;
     68:   PetscCall(MatMult(tao->hessian, X, G));
     69:   PetscCall(VecDot(G, X, &f1));
     70:   PetscCall(VecDot(gpcg->B, X, &f2));
     71:   PetscCall(VecAXPY(G, 1.0, gpcg->B));
     72:   *f = f1 / 2.0 + f2 + gpcg->c;
     73:   PetscFunctionReturn(PETSC_SUCCESS);
     74: }

     76: static PetscErrorCode TaoSetup_GPCG(Tao tao)
     77: {
     78:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;

     80:   PetscFunctionBegin;
     81:   /* Allocate some arrays */
     82:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
     83:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));

     85:   PetscCall(VecDuplicate(tao->solution, &gpcg->B));
     86:   PetscCall(VecDuplicate(tao->solution, &gpcg->Work));
     87:   PetscCall(VecDuplicate(tao->solution, &gpcg->X_New));
     88:   PetscCall(VecDuplicate(tao->solution, &gpcg->G_New));
     89:   PetscCall(VecDuplicate(tao->solution, &gpcg->DXFree));
     90:   PetscCall(VecDuplicate(tao->solution, &gpcg->R));
     91:   PetscCall(VecDuplicate(tao->solution, &gpcg->PG));
     92:   PetscFunctionReturn(PETSC_SUCCESS);
     93: }

     95: static PetscErrorCode TaoSolve_GPCG(Tao tao)
     96: {
     97:   TAO_GPCG                    *gpcg = (TAO_GPCG *)tao->data;
     98:   PetscInt                     its;
     99:   PetscReal                    actred, f, f_new, gnorm, gdx, stepsize, xtb;
    100:   PetscReal                    xtHx;
    101:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

    103:   PetscFunctionBegin;
    104:   PetscCall(TaoComputeVariableBounds(tao));
    105:   PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
    106:   PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));

    108:   /* Using f = .5*x'Hx + x'b + c and g=Hx + b,  compute b,c */
    109:   PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
    110:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
    111:   PetscCall(VecCopy(tao->gradient, gpcg->B));
    112:   PetscCall(MatMult(tao->hessian, tao->solution, gpcg->Work));
    113:   PetscCall(VecDot(gpcg->Work, tao->solution, &xtHx));
    114:   PetscCall(VecAXPY(gpcg->B, -1.0, gpcg->Work));
    115:   PetscCall(VecDot(gpcg->B, tao->solution, &xtb));
    116:   gpcg->c = f - xtHx / 2.0 - xtb;
    117:   PetscCall(ISDestroy(&gpcg->Free_Local));
    118:   PetscCall(VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &gpcg->Free_Local));

    120:   /* Project the gradient and calculate the norm */
    121:   PetscCall(VecCopy(tao->gradient, gpcg->G_New));
    122:   PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG));
    123:   PetscCall(VecNorm(gpcg->PG, NORM_2, &gpcg->gnorm));
    124:   tao->step = 1.0;
    125:   gpcg->f   = f;

    127:   /* Check Stopping Condition      */
    128:   tao->reason = TAO_CONTINUE_ITERATING;
    129:   PetscCall(TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its));
    130:   PetscCall(TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step));
    131:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

    133:   while (tao->reason == TAO_CONTINUE_ITERATING) {
    134:     /* Call general purpose update function */
    135:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
    136:     tao->ksp_its = 0;

    138:     PetscCall(GPCGGradProjections(tao));
    139:     PetscCall(ISGetSize(gpcg->Free_Local, &gpcg->n_free));

    141:     f     = gpcg->f;
    142:     gnorm = gpcg->gnorm;

    144:     PetscCall(KSPReset(tao->ksp));

    146:     if (gpcg->n_free > 0) {
    147:       /* Create a reduced linear system */
    148:       PetscCall(VecDestroy(&gpcg->R));
    149:       PetscCall(VecDestroy(&gpcg->DXFree));
    150:       PetscCall(TaoVecGetSubVec(tao->gradient, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R));
    151:       PetscCall(VecScale(gpcg->R, -1.0));
    152:       PetscCall(TaoVecGetSubVec(tao->stepdirection, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->DXFree));
    153:       PetscCall(VecSet(gpcg->DXFree, 0.0));

    155:       PetscCall(TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub));

    157:       if (tao->hessian_pre == tao->hessian) {
    158:         PetscCall(MatDestroy(&gpcg->Hsub_pre));
    159:         PetscCall(PetscObjectReference((PetscObject)gpcg->Hsub));
    160:         gpcg->Hsub_pre = gpcg->Hsub;
    161:       } else {
    162:         PetscCall(TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre));
    163:       }

    165:       PetscCall(KSPReset(tao->ksp));
    166:       PetscCall(KSPSetOperators(tao->ksp, gpcg->Hsub, gpcg->Hsub_pre));

    168:       PetscCall(KSPSolve(tao->ksp, gpcg->R, gpcg->DXFree));
    169:       PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    170:       tao->ksp_its += its;
    171:       tao->ksp_tot_its += its;
    172:       PetscCall(VecSet(tao->stepdirection, 0.0));
    173:       PetscCall(VecISAXPY(tao->stepdirection, gpcg->Free_Local, 1.0, gpcg->DXFree));

    175:       PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
    176:       PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
    177:       f_new = f;
    178:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, &stepsize, &ls_status));

    180:       actred = f_new - f;

    182:       /* Evaluate the function and gradient at the new point */
    183:       PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG));
    184:       PetscCall(VecNorm(gpcg->PG, NORM_2, &gnorm));
    185:       f = f_new;
    186:       PetscCall(ISDestroy(&gpcg->Free_Local));
    187:       PetscCall(VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &gpcg->Free_Local));
    188:     } else {
    189:       actred     = 0;
    190:       gpcg->step = 1.0;
    191:       /* if there were no free variables, no cg method */
    192:     }

    194:     tao->niter++;
    195:     gpcg->f      = f;
    196:     gpcg->gnorm  = gnorm;
    197:     gpcg->actred = actred;
    198:     PetscCall(TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its));
    199:     PetscCall(TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step));
    200:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    201:     if (tao->reason != TAO_CONTINUE_ITERATING) break;
    202:   } /* END MAIN LOOP  */
    203:   PetscFunctionReturn(PETSC_SUCCESS);
    204: }

    206: static PetscErrorCode GPCGGradProjections(Tao tao)
    207: {
    208:   TAO_GPCG                    *gpcg = (TAO_GPCG *)tao->data;
    209:   PetscInt                     i;
    210:   PetscReal                    actred = -1.0, actred_max = 0.0, gAg, gtg = gpcg->gnorm, alpha;
    211:   PetscReal                    f_new, gdx, stepsize;
    212:   Vec                          DX = tao->stepdirection, XL = tao->XL, XU = tao->XU, Work = gpcg->Work;
    213:   Vec                          X = tao->solution, G = tao->gradient;
    214:   TaoLineSearchConvergedReason lsflag = TAOLINESEARCH_CONTINUE_ITERATING;

    216:   /*
    217:      The free, active, and binding variables should be already identified
    218:   */
    219:   PetscFunctionBegin;
    220:   for (i = 0; i < gpcg->maxgpits; i++) {
    221:     if (-actred <= (gpcg->pg_ftol) * actred_max) break;
    222:     PetscCall(VecBoundGradientProjection(G, X, XL, XU, DX));
    223:     PetscCall(VecScale(DX, -1.0));
    224:     PetscCall(VecDot(DX, G, &gdx));

    226:     PetscCall(MatMult(tao->hessian, DX, Work));
    227:     PetscCall(VecDot(DX, Work, &gAg));

    229:     gpcg->gp_iterates++;
    230:     gpcg->total_gp_its++;

    232:     gtg = -gdx;
    233:     if (PetscAbsReal(gAg) == 0.0) {
    234:       alpha = 1.0;
    235:     } else {
    236:       alpha = PetscAbsReal(gtg / gAg);
    237:     }
    238:     PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, alpha));
    239:     f_new = gpcg->f;
    240:     PetscCall(TaoLineSearchApply(tao->linesearch, X, &f_new, G, DX, &stepsize, &lsflag));

    242:     /* Update the iterate */
    243:     actred     = f_new - gpcg->f;
    244:     actred_max = PetscMax(actred_max, -(f_new - gpcg->f));
    245:     gpcg->f    = f_new;
    246:     PetscCall(ISDestroy(&gpcg->Free_Local));
    247:     PetscCall(VecWhichInactive(XL, X, tao->gradient, XU, PETSC_TRUE, &gpcg->Free_Local));
    248:   }

    250:   gpcg->gnorm = gtg;
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: } /* End gradient projections */

    254: static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
    255: {
    256:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;

    258:   PetscFunctionBegin;
    259:   PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work));
    260:   PetscCall(VecCopy(gpcg->Work, DXL));
    261:   PetscCall(VecAXPY(DXL, -1.0, tao->gradient));
    262:   PetscCall(VecSet(DXU, 0.0));
    263:   PetscCall(VecPointwiseMax(DXL, DXL, DXU));

    265:   PetscCall(VecCopy(tao->gradient, DXU));
    266:   PetscCall(VecAXPY(DXU, -1.0, gpcg->Work));
    267:   PetscCall(VecSet(gpcg->Work, 0.0));
    268:   PetscCall(VecPointwiseMin(DXU, gpcg->Work, DXU));
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: }

    272: /*MC
    273:   TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
    274:             conjugate-gradient based method for bound-constrained minimization

    276:   Options Database Keys:
    277: + -tao_gpcg_maxpgits - maximum number of gradient projections for GPCG iterate
    278: - -tao_subset_type   - "subvec","mask","matrix-free", strategies for handling active-sets

    280:   Level: beginner

    282: .seealso: `Tao`, `TaoType`, `TAOTRON`, `TAOBQPIP`, `TAOLINESEARCHGPCG`
    283: M*/
    284: PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
    285: {
    286:   TAO_GPCG *gpcg;

    288:   PetscFunctionBegin;
    289:   tao->ops->setup            = TaoSetup_GPCG;
    290:   tao->ops->solve            = TaoSolve_GPCG;
    291:   tao->ops->view             = TaoView_GPCG;
    292:   tao->ops->setfromoptions   = TaoSetFromOptions_GPCG;
    293:   tao->ops->destroy          = TaoDestroy_GPCG;
    294:   tao->ops->computedual      = TaoComputeDual_GPCG;
    295:   tao->uses_gradient         = PETSC_TRUE;
    296:   tao->uses_hessian_matrices = PETSC_TRUE;

    298:   PetscCall(PetscNew(&gpcg));
    299:   tao->data = (void *)gpcg;

    301:   /* Override default settings (unless already changed) */
    302:   PetscCall(TaoParametersInitialize(tao));
    303:   PetscObjectParameterSetDefault(tao, max_it, 500);
    304:   PetscObjectParameterSetDefault(tao, max_funcs, 100000);
    305:   PetscObjectParameterSetDefault(tao, gatol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);
    306:   PetscObjectParameterSetDefault(tao, grtol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);

    308:   /* Initialize pointers and variables */
    309:   gpcg->n        = 0;
    310:   gpcg->maxgpits = 8;
    311:   gpcg->pg_ftol  = 0.1;

    313:   gpcg->gp_iterates  = 0; /* Cumulative number */
    314:   gpcg->total_gp_its = 0;

    316:   /* Initialize pointers and variables */
    317:   gpcg->n_bind      = 0;
    318:   gpcg->n_free      = 0;
    319:   gpcg->n_upper     = 0;
    320:   gpcg->n_lower     = 0;
    321:   gpcg->subset_type = TAO_SUBSET_MASK;
    322:   gpcg->Hsub        = NULL;
    323:   gpcg->Hsub_pre    = NULL;

    325:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
    326:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
    327:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
    328:   PetscCall(KSPSetType(tao->ksp, KSPNASH));

    330:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
    331:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
    332:   PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG));
    333:   PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao));
    334:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
    335:   PetscFunctionReturn(PETSC_SUCCESS);
    336: }