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:   VecDestroy(&gpcg->B);
 14:   VecDestroy(&gpcg->Work);
 15:   VecDestroy(&gpcg->X_New);
 16:   VecDestroy(&gpcg->G_New);
 17:   VecDestroy(&gpcg->DXFree);
 18:   VecDestroy(&gpcg->R);
 19:   VecDestroy(&gpcg->PG);
 20:   MatDestroy(&gpcg->Hsub);
 21:   MatDestroy(&gpcg->Hsub_pre);
 22:   ISDestroy(&gpcg->Free_Local);
 23:   KSPDestroy(&tao->ksp);
 24:   PetscFree(tao->data);
 25:   return 0;
 26: }

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

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

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

 48:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
 49:   if (isascii) {
 50:     PetscViewerASCIIPrintf(viewer, "Total PG its: %" PetscInt_FMT ",", gpcg->total_gp_its);
 51:     PetscViewerASCIIPrintf(viewer, "PG tolerance: %g \n", (double)gpcg->pg_ftol);
 52:   }
 53:   TaoLineSearchView(tao->linesearch, viewer);
 54:   return 0;
 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:   MatMult(tao->hessian, X, G);
 68:   VecDot(G, X, &f1);
 69:   VecDot(gpcg->B, X, &f2);
 70:   VecAXPY(G, 1.0, gpcg->B);
 71:   *f = f1 / 2.0 + f2 + gpcg->c;
 72:   return 0;
 73: }

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

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

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

103:     }
104:   */
105:   return 0;
106: }

108: static PetscErrorCode TaoSolve_GPCG(Tao tao)
109: {
110:   TAO_GPCG                    *gpcg = (TAO_GPCG *)tao->data;
111:   PetscInt                     its;
112:   PetscReal                    actred, f, f_new, gnorm, gdx, stepsize, xtb;
113:   PetscReal                    xtHx;
114:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;


117:   TaoComputeVariableBounds(tao);
118:   VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
119:   TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU);

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

133:   /* Project the gradient and calculate the norm */
134:   VecCopy(tao->gradient, gpcg->G_New);
135:   VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->PG);
136:   VecNorm(gpcg->PG, NORM_2, &gpcg->gnorm);
137:   tao->step = 1.0;
138:   gpcg->f   = f;

140:   /* Check Stopping Condition      */
141:   tao->reason = TAO_CONTINUE_ITERATING;
142:   TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its);
143:   TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step);
144:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

146:   while (tao->reason == TAO_CONTINUE_ITERATING) {
147:     /* Call general purpose update function */
148:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
149:     tao->ksp_its = 0;

151:     GPCGGradProjections(tao);
152:     ISGetSize(gpcg->Free_Local, &gpcg->n_free);

154:     f     = gpcg->f;
155:     gnorm = gpcg->gnorm;

157:     KSPReset(tao->ksp);

159:     if (gpcg->n_free > 0) {
160:       /* Create a reduced linear system */
161:       VecDestroy(&gpcg->R);
162:       VecDestroy(&gpcg->DXFree);
163:       TaoVecGetSubVec(tao->gradient, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->R);
164:       VecScale(gpcg->R, -1.0);
165:       TaoVecGetSubVec(tao->stepdirection, gpcg->Free_Local, tao->subset_type, 0.0, &gpcg->DXFree);
166:       VecSet(gpcg->DXFree, 0.0);

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

170:       if (tao->hessian_pre == tao->hessian) {
171:         MatDestroy(&gpcg->Hsub_pre);
172:         PetscObjectReference((PetscObject)gpcg->Hsub);
173:         gpcg->Hsub_pre = gpcg->Hsub;
174:       } else {
175:         TaoMatGetSubMat(tao->hessian, gpcg->Free_Local, gpcg->Work, tao->subset_type, &gpcg->Hsub_pre);
176:       }

178:       KSPReset(tao->ksp);
179:       KSPSetOperators(tao->ksp, gpcg->Hsub, gpcg->Hsub_pre);

181:       KSPSolve(tao->ksp, gpcg->R, gpcg->DXFree);
182:       KSPGetIterationNumber(tao->ksp, &its);
183:       tao->ksp_its += its;
184:       tao->ksp_tot_its += its;
185:       VecSet(tao->stepdirection, 0.0);
186:       VecISAXPY(tao->stepdirection, gpcg->Free_Local, 1.0, gpcg->DXFree);

188:       VecDot(tao->stepdirection, tao->gradient, &gdx);
189:       TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
190:       f_new = f;
191:       TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, &stepsize, &ls_status);

193:       actred = f_new - f;

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

207:     tao->niter++;
208:     gpcg->f      = f;
209:     gpcg->gnorm  = gnorm;
210:     gpcg->actred = actred;
211:     TaoLogConvergenceHistory(tao, f, gpcg->gnorm, 0.0, tao->ksp_its);
212:     TaoMonitor(tao, tao->niter, f, gpcg->gnorm, 0.0, tao->step);
213:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
214:     if (tao->reason != TAO_CONTINUE_ITERATING) break;
215:   } /* END MAIN LOOP  */

217:   return 0;
218: }

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

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

239:     MatMult(tao->hessian, DX, Work);
240:     VecDot(DX, Work, &gAg);

242:     gpcg->gp_iterates++;
243:     gpcg->total_gp_its++;

245:     gtg = -gdx;
246:     if (PetscAbsReal(gAg) == 0.0) {
247:       alpha = 1.0;
248:     } else {
249:       alpha = PetscAbsReal(gtg / gAg);
250:     }
251:     TaoLineSearchSetInitialStepLength(tao->linesearch, alpha);
252:     f_new = gpcg->f;
253:     TaoLineSearchApply(tao->linesearch, X, &f_new, G, DX, &stepsize, &lsflag);

255:     /* Update the iterate */
256:     actred     = f_new - gpcg->f;
257:     actred_max = PetscMax(actred_max, -(f_new - gpcg->f));
258:     gpcg->f    = f_new;
259:     ISDestroy(&gpcg->Free_Local);
260:     VecWhichInactive(XL, X, tao->gradient, XU, PETSC_TRUE, &gpcg->Free_Local);
261:   }

263:   gpcg->gnorm = gtg;
264:   return 0;
265: } /* End gradient projections */

267: static PetscErrorCode TaoComputeDual_GPCG(Tao tao, Vec DXL, Vec DXU)
268: {
269:   TAO_GPCG *gpcg = (TAO_GPCG *)tao->data;

271:   VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, gpcg->Work);
272:   VecCopy(gpcg->Work, DXL);
273:   VecAXPY(DXL, -1.0, tao->gradient);
274:   VecSet(DXU, 0.0);
275:   VecPointwiseMax(DXL, DXL, DXU);

277:   VecCopy(tao->gradient, DXU);
278:   VecAXPY(DXU, -1.0, gpcg->Work);
279:   VecSet(gpcg->Work, 0.0);
280:   VecPointwiseMin(DXU, gpcg->Work, DXU);
281:   return 0;
282: }

284: /*------------------------------------------------------------*/
285: /*MC
286:   TAOGPCG - gradient projected conjugate gradient algorithm is an active-set
287:         conjugate-gradient based method for bound-constrained minimization

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

293:   Level: beginner
294: M*/
295: PETSC_EXTERN PetscErrorCode TaoCreate_GPCG(Tao tao)
296: {
297:   TAO_GPCG *gpcg;

299:   tao->ops->setup          = TaoSetup_GPCG;
300:   tao->ops->solve          = TaoSolve_GPCG;
301:   tao->ops->view           = TaoView_GPCG;
302:   tao->ops->setfromoptions = TaoSetFromOptions_GPCG;
303:   tao->ops->destroy        = TaoDestroy_GPCG;
304:   tao->ops->computedual    = TaoComputeDual_GPCG;

306:   PetscNew(&gpcg);
307:   tao->data = (void *)gpcg;

309:   /* Override default settings (unless already changed) */
310:   if (!tao->max_it_changed) tao->max_it = 500;
311:   if (!tao->max_funcs_changed) tao->max_funcs = 100000;
312: #if defined(PETSC_USE_REAL_SINGLE)
313:   if (!tao->gatol_changed) tao->gatol = 1e-6;
314:   if (!tao->grtol_changed) tao->grtol = 1e-6;
315: #else
316:   if (!tao->gatol_changed) tao->gatol = 1e-12;
317:   if (!tao->grtol_changed) tao->grtol = 1e-12;
318: #endif

320:   /* Initialize pointers and variables */
321:   gpcg->n        = 0;
322:   gpcg->maxgpits = 8;
323:   gpcg->pg_ftol  = 0.1;

325:   gpcg->gp_iterates  = 0; /* Cumulative number */
326:   gpcg->total_gp_its = 0;

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

337:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
338:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
339:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
340:   KSPSetType(tao->ksp, KSPNASH);

342:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
343:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
344:   TaoLineSearchSetType(tao->linesearch, TAOLINESEARCHGPCG);
345:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch, GPCGObjectiveAndGradient, tao);
346:   TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);
347:   return 0;
348: }