Actual source code: tron.c

  1: #include <../src/tao/bound/impls/tron/tron.h>
  2: #include <../src/tao/matrix/submatfree.h>

  4: /* TRON Routines */
  5: static PetscErrorCode TronGradientProjections(Tao,TAO_TRON*);
  6: /*------------------------------------------------------------*/
  7: static PetscErrorCode TaoDestroy_TRON(Tao tao)
  8: {
  9:   TAO_TRON       *tron = (TAO_TRON *)tao->data;

 11:   VecDestroy(&tron->X_New);
 12:   VecDestroy(&tron->G_New);
 13:   VecDestroy(&tron->Work);
 14:   VecDestroy(&tron->DXFree);
 15:   VecDestroy(&tron->R);
 16:   VecDestroy(&tron->diag);
 17:   VecScatterDestroy(&tron->scatter);
 18:   ISDestroy(&tron->Free_Local);
 19:   MatDestroy(&tron->H_sub);
 20:   MatDestroy(&tron->Hpre_sub);
 21:   PetscFree(tao->data);
 22:   return 0;
 23: }

 25: /*------------------------------------------------------------*/
 26: static PetscErrorCode TaoSetFromOptions_TRON(PetscOptionItems *PetscOptionsObject,Tao tao)
 27: {
 28:   TAO_TRON       *tron = (TAO_TRON *)tao->data;
 29:   PetscBool      flg;

 31:   PetscOptionsHead(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization");
 32:   PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg);
 33:   PetscOptionsTail();
 34:   KSPSetFromOptions(tao->ksp);
 35:   return 0;
 36: }

 38: /*------------------------------------------------------------*/
 39: static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer)
 40: {
 41:   TAO_TRON         *tron = (TAO_TRON *)tao->data;
 42:   PetscBool        isascii;

 44:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
 45:   if (isascii) {
 46:     PetscViewerASCIIPrintf(viewer,"Total PG its: %D,",tron->total_gp_its);
 47:     PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol);
 48:   }
 49:   return 0;
 50: }

 52: /* ---------------------------------------------------------- */
 53: static PetscErrorCode TaoSetup_TRON(Tao tao)
 54: {
 55:   TAO_TRON       *tron = (TAO_TRON *)tao->data;


 58:   /* Allocate some arrays */
 59:   VecDuplicate(tao->solution, &tron->diag);
 60:   VecDuplicate(tao->solution, &tron->X_New);
 61:   VecDuplicate(tao->solution, &tron->G_New);
 62:   VecDuplicate(tao->solution, &tron->Work);
 63:   VecDuplicate(tao->solution, &tao->gradient);
 64:   VecDuplicate(tao->solution, &tao->stepdirection);
 65:   if (!tao->XL) {
 66:     VecDuplicate(tao->solution, &tao->XL);
 67:     VecSet(tao->XL, PETSC_NINFINITY);
 68:   }
 69:   if (!tao->XU) {
 70:     VecDuplicate(tao->solution, &tao->XU);
 71:     VecSet(tao->XU, PETSC_INFINITY);
 72:   }
 73:   return 0;
 74: }

 76: static PetscErrorCode TaoSolve_TRON(Tao tao)
 77: {
 78:   TAO_TRON                     *tron = (TAO_TRON *)tao->data;
 79:   PetscInt                     its;
 80:   TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
 81:   PetscReal                    prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize;

 83:   tron->pgstepsize = 1.0;
 84:   tao->trust = tao->trust0;
 85:   /*   Project the current point onto the feasible set */
 86:   TaoComputeVariableBounds(tao);
 87:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);

 89:   /* Project the initial point onto the feasible region */
 90:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);

 92:   /* Compute the objective function and gradient */
 93:   TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient);
 94:   VecNorm(tao->gradient,NORM_2,&tron->gnorm);

 97:   /* Project the gradient and calculate the norm */
 98:   VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
 99:   VecNorm(tao->gradient,NORM_2,&tron->gnorm);

101:   /* Initialize trust region radius */
102:   tao->trust=tao->trust0;
103:   if (tao->trust <= 0) {
104:     tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0);
105:   }

107:   /* Initialize step sizes for the line searches */
108:   tron->pgstepsize=1.0;
109:   tron->stepsize=tao->trust;

111:   tao->reason = TAO_CONTINUE_ITERATING;
112:   TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
113:   TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize);
114:   (*tao->ops->convergencetest)(tao,tao->cnvP);
115:   while (tao->reason==TAO_CONTINUE_ITERATING) {
116:     /* Call general purpose update function */
117:     if (tao->ops->update) {
118:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
119:     }

121:     /* Perform projected gradient iterations */
122:     TronGradientProjections(tao,tron);

124:     VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
125:     VecNorm(tao->gradient,NORM_2,&tron->gnorm);

127:     tao->ksp_its=0;
128:     f=tron->f; delta=tao->trust;
129:     tron->n_free_last = tron->n_free;
130:     TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);

132:     /* Generate index set (IS) of which bound constraints are active */
133:     ISDestroy(&tron->Free_Local);
134:     VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
135:     ISGetSize(tron->Free_Local, &tron->n_free);

137:     /* If no free variables */
138:     if (tron->n_free == 0) {
139:       VecNorm(tao->gradient,NORM_2,&tron->gnorm);
140:       TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
141:       TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta);
142:       (*tao->ops->convergencetest)(tao,tao->cnvP);
143:       if (!tao->reason) {
144:         tao->reason = TAO_CONVERGED_STEPTOL;
145:       }
146:       break;
147:     }
148:     /* use free_local to mask/submat gradient, hessian, stepdirection */
149:     TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R);
150:     TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree);
151:     VecSet(tron->DXFree,0.0);
152:     VecScale(tron->R, -1.0);
153:     TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);
154:     if (tao->hessian == tao->hessian_pre) {
155:       MatDestroy(&tron->Hpre_sub);
156:       PetscObjectReference((PetscObject)(tron->H_sub));
157:       tron->Hpre_sub = tron->H_sub;
158:     } else {
159:       TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub);
160:     }
161:     KSPReset(tao->ksp);
162:     KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);
163:     while (1) {

165:       /* Approximately solve the reduced linear system */
166:       KSPCGSetRadius(tao->ksp,delta);

168:       KSPSolve(tao->ksp, tron->R, tron->DXFree);
169:       KSPGetIterationNumber(tao->ksp,&its);
170:       tao->ksp_its+=its;
171:       tao->ksp_tot_its+=its;
172:       VecSet(tao->stepdirection,0.0);

174:       /* Add dxfree matrix to compute step direction vector */
175:       VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree);

177:       VecDot(tao->gradient, tao->stepdirection, &gdx);
178:       VecCopy(tao->solution, tron->X_New);
179:       VecCopy(tao->gradient, tron->G_New);

181:       stepsize=1.0;f_new=f;

183:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
184:       TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason);
185:       TaoAddLineSearchCounts(tao);

187:       MatMult(tao->hessian, tao->stepdirection, tron->Work);
188:       VecAYPX(tron->Work, 0.5, tao->gradient);
189:       VecDot(tao->stepdirection, tron->Work, &prered);
190:       actred = f_new - f;
191:       if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) {
192:         rhok = 1.0;
193:       } else if (actred<0) {
194:         rhok=PetscAbs(-actred/prered);
195:       } else {
196:         rhok=0.0;
197:       }

199:       /* Compare actual improvement to the quadratic model */
200:       if (rhok > tron->eta1) { /* Accept the point */
201:         /* d = x_new - x */
202:         VecCopy(tron->X_New, tao->stepdirection);
203:         VecAXPY(tao->stepdirection, -1.0, tao->solution);

205:         VecNorm(tao->stepdirection, NORM_2, &xdiff);
206:         xdiff *= stepsize;

208:         /* Adjust trust region size */
209:         if (rhok < tron->eta2) {
210:           delta = PetscMin(xdiff,delta)*tron->sigma1;
211:         } else if (rhok > tron->eta4) {
212:           delta= PetscMin(xdiff,delta)*tron->sigma3;
213:         } else if (rhok > tron->eta3) {
214:           delta=PetscMin(xdiff,delta)*tron->sigma2;
215:         }
216:         VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient);
217:         ISDestroy(&tron->Free_Local);
218:         VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local);
219:         f=f_new;
220:         VecNorm(tao->gradient,NORM_2,&tron->gnorm);
221:         VecCopy(tron->X_New, tao->solution);
222:         VecCopy(tron->G_New, tao->gradient);
223:         break;
224:       }
225:       else if (delta <= 1e-30) {
226:         break;
227:       }
228:       else {
229:         delta /= 4.0;
230:       }
231:     } /* end linear solve loop */

233:     tron->f=f; tron->actred=actred; tao->trust=delta;
234:     tao->niter++;
235:     TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its);
236:     TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,stepsize);
237:     (*tao->ops->convergencetest)(tao,tao->cnvP);
238:   }  /* END MAIN LOOP  */
239:   return 0;
240: }

242: static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron)
243: {
244:   PetscErrorCode               ierr;
245:   PetscInt                     i;
246:   TaoLineSearchConvergedReason ls_reason;
247:   PetscReal                    actred=-1.0,actred_max=0.0;
248:   PetscReal                    f_new;
249:   /*
250:      The gradient and function value passed into and out of this
251:      routine should be current and correct.

253:      The free, active, and binding variables should be already identified
254:   */

256:   for (i=0;i<tron->maxgpits;++i) {

258:     if (-actred <= (tron->pg_ftol)*actred_max) break;

260:     ++tron->gp_iterates;
261:     ++tron->total_gp_its;
262:     f_new=tron->f;

264:     VecCopy(tao->gradient,tao->stepdirection);
265:     VecScale(tao->stepdirection,-1.0);
266:     TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize);
267:     TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,
268:                               &tron->pgstepsize, &ls_reason);
269:     TaoAddLineSearchCounts(tao);

271:     VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient);
272:     VecNorm(tao->gradient,NORM_2,&tron->gnorm);

274:     /* Update the iterate */
275:     actred = f_new - tron->f;
276:     actred_max = PetscMax(actred_max,-(f_new - tron->f));
277:     tron->f = f_new;
278:   }
279:   return 0;
280: }

282: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU)
283: {

285:   TAO_TRON       *tron = (TAO_TRON *)tao->data;


292:   VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work);
293:   VecCopy(tron->Work,DXL);
294:   VecAXPY(DXL,-1.0,tao->gradient);
295:   VecSet(DXU,0.0);
296:   VecPointwiseMax(DXL,DXL,DXU);

298:   VecCopy(tao->gradient,DXU);
299:   VecAXPY(DXU,-1.0,tron->Work);
300:   VecSet(tron->Work,0.0);
301:   VecPointwiseMin(DXU,tron->Work,DXU);
302:   return 0;
303: }

305: /*------------------------------------------------------------*/
306: /*MC
307:   TAOTRON - The TRON algorithm is an active-set Newton trust region method
308:   for bound-constrained minimization.

310:   Options Database Keys:
311: + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate
312: - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets

314:   Level: beginner
315: M*/
316: PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
317: {
318:   TAO_TRON       *tron;
319:   const char     *morethuente_type = TAOLINESEARCHMT;

321:   tao->ops->setup          = TaoSetup_TRON;
322:   tao->ops->solve          = TaoSolve_TRON;
323:   tao->ops->view           = TaoView_TRON;
324:   tao->ops->setfromoptions = TaoSetFromOptions_TRON;
325:   tao->ops->destroy        = TaoDestroy_TRON;
326:   tao->ops->computedual    = TaoComputeDual_TRON;

328:   PetscNewLog(tao,&tron);
329:   tao->data = (void*)tron;

331:   /* Override default settings (unless already changed) */
332:   if (!tao->max_it_changed) tao->max_it = 50;
333:   if (!tao->trust0_changed) tao->trust0 = 1.0;
334:   if (!tao->steptol_changed) tao->steptol = 0.0;

336:   /* Initialize pointers and variables */
337:   tron->n            = 0;
338:   tron->maxgpits     = 3;
339:   tron->pg_ftol      = 0.001;

341:   tron->eta1         = 1.0e-4;
342:   tron->eta2         = 0.25;
343:   tron->eta3         = 0.50;
344:   tron->eta4         = 0.90;

346:   tron->sigma1       = 0.5;
347:   tron->sigma2       = 2.0;
348:   tron->sigma3       = 4.0;

350:   tron->gp_iterates  = 0; /* Cumulative number */
351:   tron->total_gp_its = 0;
352:   tron->n_free       = 0;

354:   tron->DXFree=NULL;
355:   tron->R=NULL;
356:   tron->X_New=NULL;
357:   tron->G_New=NULL;
358:   tron->Work=NULL;
359:   tron->Free_Local=NULL;
360:   tron->H_sub=NULL;
361:   tron->Hpre_sub=NULL;
362:   tao->subset_type = TAO_SUBSET_SUBVEC;

364:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
365:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
366:   TaoLineSearchSetType(tao->linesearch,morethuente_type);
367:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
368:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

370:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
371:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
372:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
373:   KSPSetType(tao->ksp,KSPSTCG);
374:   return 0;
375: }