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: /*------------------------------------------------------------*/
      8: static PetscErrorCode TaoDestroy_GPCG(Tao tao)
      9: {
     10:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;

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

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

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

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

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

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

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

     79: /* ---------------------------------------------------------- */
     80: static PetscErrorCode TaoSetup_GPCG(Tao tao)
     81: {
     82:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;

     84:   PetscFunctionBegin;
     85:   /* Allocate some arrays */
     86:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
     87:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));

     89:   PetscCall(VecDuplicate(tao->solution, &gpcg->B));
     90:   PetscCall(VecDuplicate(tao->solution, &gpcg->Work));
     91:   PetscCall(VecDuplicate(tao->solution, &gpcg->X_New));
     92:   PetscCall(VecDuplicate(tao->solution, &gpcg->G_New));
     93:   PetscCall(VecDuplicate(tao->solution, &gpcg->DXFree));
     94:   PetscCall(VecDuplicate(tao->solution, &gpcg->R));
     95:   PetscCall(VecDuplicate(tao->solution, &gpcg->PG));
     96:   /*
     97:     if (gpcg->ksp_type == GPCG_KSP_NASH) {
     98:         PetscCall(KSPSetType(tao->ksp,KSPNASH));
     99:       } else if (gpcg->ksp_type == GPCG_KSP_STCG) {
    100:         PetscCall(KSPSetType(tao->ksp,KSPSTCG));
    101:       } else {
    102:         PetscCall(KSPSetType(tao->ksp,KSPGLTR));
    103:       }
    104:       if (tao->ksp->ops->setfromoptions) {
    105:         (*tao->ksp->ops->setfromoptions)(tao->ksp);
    106:       }

    108:     }
    109:   */
    110:   PetscFunctionReturn(PETSC_SUCCESS);
    111: }

    113: static PetscErrorCode TaoSolve_GPCG(Tao tao)
    114: {
    115:   TAO_GPCG                    *gpcg = (TAO_GPCG *)tao->data;
    116:   PetscInt                     its;
    117:   PetscReal                    actred, f, f_new, gnorm, gdx, stepsize, xtb;
    118:   PetscReal                    xtHx;
    119:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

    121:   PetscFunctionBegin;
    122:   PetscCall(TaoComputeVariableBounds(tao));
    123:   PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
    124:   PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));

    126:   /* Using f = .5*x'Hx + x'b + c and g=Hx + b,  compute b,c */
    127:   PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));
    128:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
    129:   PetscCall(VecCopy(tao->gradient, gpcg->B));
    130:   PetscCall(MatMult(tao->hessian, tao->solution, gpcg->Work));
    131:   PetscCall(VecDot(gpcg->Work, tao->solution, &xtHx));
    132:   PetscCall(VecAXPY(gpcg->B, -1.0, gpcg->Work));
    133:   PetscCall(VecDot(gpcg->B, tao->solution, &xtb));
    134:   gpcg->c = f - xtHx / 2.0 - xtb;
    135:   if (gpcg->Free_Local) PetscCall(ISDestroy(&gpcg->Free_Local));
    136:   PetscCall(VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &gpcg->Free_Local));

    138:   /* Project the gradient and calculate the norm */
    139:   PetscCall(VecCopy(tao->gradient, gpcg->G_New));
    140:   PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG));
    141:   PetscCall(VecNorm(gpcg->PG, NORM_2, &gpcg->gnorm));
    142:   tao->step = 1.0;
    143:   gpcg->f   = f;

    145:   /* Check Stopping Condition      */
    146:   tao->reason = TAO_CONTINUE_ITERATING;
    147:   PetscCall(TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its));
    148:   PetscCall(TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step));
    149:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

    151:   while (tao->reason == TAO_CONTINUE_ITERATING) {
    152:     /* Call general purpose update function */
    153:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
    154:     tao->ksp_its = 0;

    156:     PetscCall(GPCGGradProjections(tao));
    157:     PetscCall(ISGetSize(gpcg->Free_Local, &gpcg->n_free));

    159:     f     = gpcg->f;
    160:     gnorm = gpcg->gnorm;

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

    164:     if (gpcg->n_free > 0) {
    165:       /* Create a reduced linear system */
    166:       PetscCall(VecDestroy(&gpcg->R));
    167:       PetscCall(VecDestroy(&gpcg->DXFree));
    168:       PetscCall(TaoVecGetSubVec(tao->gradient, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R));
    169:       PetscCall(VecScale(gpcg->R, -1.0));
    170:       PetscCall(TaoVecGetSubVec(tao->stepdirection, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->DXFree));
    171:       PetscCall(VecSet(gpcg->DXFree, 0.0));

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

    175:       if (tao->hessian_pre == tao->hessian) {
    176:         PetscCall(MatDestroy(&gpcg->Hsub_pre));
    177:         PetscCall(PetscObjectReference((PetscObject)gpcg->Hsub));
    178:         gpcg->Hsub_pre = gpcg->Hsub;
    179:       } else {
    180:         PetscCall(TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre));
    181:       }

    183:       PetscCall(KSPReset(tao->ksp));
    184:       PetscCall(KSPSetOperators(tao->ksp, gpcg->Hsub, gpcg->Hsub_pre));

    186:       PetscCall(KSPSolve(tao->ksp, gpcg->R, gpcg->DXFree));
    187:       PetscCall(KSPGetIterationNumber(tao->ksp, &its));
    188:       tao->ksp_its += its;
    189:       tao->ksp_tot_its += its;
    190:       PetscCall(VecSet(tao->stepdirection, 0.0));
    191:       PetscCall(VecISAXPY(tao->stepdirection, gpcg->Free_Local, 1.0, gpcg->DXFree));

    193:       PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx));
    194:       PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
    195:       f_new = f;
    196:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, &stepsize, &ls_status));

    198:       actred = f_new - f;

    200:       /* Evaluate the function and gradient at the new point */
    201:       PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG));
    202:       PetscCall(VecNorm(gpcg->PG, NORM_2, &gnorm));
    203:       f = f_new;
    204:       PetscCall(ISDestroy(&gpcg->Free_Local));
    205:       PetscCall(VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &gpcg->Free_Local));
    206:     } else {
    207:       actred     = 0;
    208:       gpcg->step = 1.0;
    209:       /* if there were no free variables, no cg method */
    210:     }

    212:     tao->niter++;
    213:     gpcg->f      = f;
    214:     gpcg->gnorm  = gnorm;
    215:     gpcg->actred = actred;
    216:     PetscCall(TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its));
    217:     PetscCall(TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step));
    218:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
    219:     if (tao->reason != TAO_CONTINUE_ITERATING) break;
    220:   } /* END MAIN LOOP  */
    221:   PetscFunctionReturn(PETSC_SUCCESS);
    222: }

    224: static PetscErrorCode GPCGGradProjections(Tao tao)
    225: {
    226:   TAO_GPCG                    *gpcg = (TAO_GPCG *)tao->data;
    227:   PetscInt                     i;
    228:   PetscReal                    actred = -1.0, actred_max = 0.0, gAg, gtg = gpcg->gnorm, alpha;
    229:   PetscReal                    f_new, gdx, stepsize;
    230:   Vec                          DX = tao->stepdirection, XL = tao->XL, XU = tao->XU, Work = gpcg->Work;
    231:   Vec                          X = tao->solution, G = tao->gradient;
    232:   TaoLineSearchConvergedReason lsflag = TAOLINESEARCH_CONTINUE_ITERATING;

    234:   /*
    235:      The free, active, and binding variables should be already identified
    236:   */
    237:   PetscFunctionBegin;
    238:   for (i = 0; i < gpcg->maxgpits; i++) {
    239:     if (-actred <= (gpcg->pg_ftol) * actred_max) break;
    240:     PetscCall(VecBoundGradientProjection(G, X, XL, XU, DX));
    241:     PetscCall(VecScale(DX, -1.0));
    242:     PetscCall(VecDot(DX, G, &gdx));

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

    247:     gpcg->gp_iterates++;
    248:     gpcg->total_gp_its++;

    250:     gtg = -gdx;
    251:     if (PetscAbsReal(gAg) == 0.0) {
    252:       alpha = 1.0;
    253:     } else {
    254:       alpha = PetscAbsReal(gtg / gAg);
    255:     }
    256:     PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, alpha));
    257:     f_new = gpcg->f;
    258:     PetscCall(TaoLineSearchApply(tao->linesearch, X, &f_new, G, DX, &stepsize, &lsflag));

    260:     /* Update the iterate */
    261:     actred     = f_new - gpcg->f;
    262:     actred_max = PetscMax(actred_max, -(f_new - gpcg->f));
    263:     gpcg->f    = f_new;
    264:     PetscCall(ISDestroy(&gpcg->Free_Local));
    265:     PetscCall(VecWhichInactive(XL, X, tao->gradient, XU, PETSC_TRUE, &gpcg->Free_Local));
    266:   }

    268:   gpcg->gnorm = gtg;
    269:   PetscFunctionReturn(PETSC_SUCCESS);
    270: } /* End gradient projections */

    272: static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
    273: {
    274:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;

    276:   PetscFunctionBegin;
    277:   PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work));
    278:   PetscCall(VecCopy(gpcg->Work, DXL));
    279:   PetscCall(VecAXPY(DXL, -1.0, tao->gradient));
    280:   PetscCall(VecSet(DXU, 0.0));
    281:   PetscCall(VecPointwiseMax(DXL, DXL, DXU));

    283:   PetscCall(VecCopy(tao->gradient, DXU));
    284:   PetscCall(VecAXPY(DXU, -1.0, gpcg->Work));
    285:   PetscCall(VecSet(gpcg->Work, 0.0));
    286:   PetscCall(VecPointwiseMin(DXU, gpcg->Work, DXU));
    287:   PetscFunctionReturn(PETSC_SUCCESS);
    288: }

    290: /*------------------------------------------------------------*/
    291: /*MC
    292:   TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
    293:         conjugate-gradient based method for bound-constrained minimization

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

    299:   Level: beginner
    300: M*/
    301: PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
    302: {
    303:   TAO_GPCG *gpcg;

    305:   PetscFunctionBegin;
    306:   tao->ops->setup          = TaoSetup_GPCG;
    307:   tao->ops->solve          = TaoSolve_GPCG;
    308:   tao->ops->view           = TaoView_GPCG;
    309:   tao->ops->setfromoptions = TaoSetFromOptions_GPCG;
    310:   tao->ops->destroy        = TaoDestroy_GPCG;
    311:   tao->ops->computedual    = TaoComputeDual_GPCG;

    313:   PetscCall(PetscNew(&gpcg));
    314:   tao->data = (void *)gpcg;

    316:   /* Override default settings (unless already changed) */
    317:   PetscCall(TaoParametersInitialize(tao));
    318:   PetscObjectParameterSetDefault(tao, max_it, 500);
    319:   PetscObjectParameterSetDefault(tao, max_funcs, 100000);
    320:   PetscObjectParameterSetDefault(tao, gatol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);
    321:   PetscObjectParameterSetDefault(tao, grtol, PetscDefined(USE_REAL_SINGLE) ? 1e-6 : 1e-12);

    323:   /* Initialize pointers and variables */
    324:   gpcg->n        = 0;
    325:   gpcg->maxgpits = 8;
    326:   gpcg->pg_ftol  = 0.1;

    328:   gpcg->gp_iterates  = 0; /* Cumulative number */
    329:   gpcg->total_gp_its = 0;

    331:   /* Initialize pointers and variables */
    332:   gpcg->n_bind      = 0;
    333:   gpcg->n_free      = 0;
    334:   gpcg->n_upper     = 0;
    335:   gpcg->n_lower     = 0;
    336:   gpcg->subset_type = TAO_SUBSET_MASK;
    337:   gpcg->Hsub        = NULL;
    338:   gpcg->Hsub_pre    = NULL;

    340:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
    341:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
    342:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
    343:   PetscCall(KSPSetType(tao->ksp, KSPNASH));

    345:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
    346:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
    347:   PetscCall(TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG));
    348:   PetscCall(TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao));
    349:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
    350:   PetscFunctionReturn(PETSC_SUCCESS);
    351: }