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:   KSPDestroy(&tao->ksp);
 22:   PetscFree(tao->data);
 23:   return 0;
 24: }

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

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

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

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

 53: /* ---------------------------------------------------------- */
 54: static PetscErrorCode TaoSetup_TRON(Tao tao)
 55: {
 56:   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:   return 0;
 66: }

 68: static PetscErrorCode TaoSolve_TRON(Tao tao)
 69: {
 70:   TAO_TRON                    *tron = (TAO_TRON *)tao->data;
 71:   PetscInt                     its;
 72:   TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING;
 73:   PetscReal                    prered, actred, delta, f, f_new, rhok, gdx, xdiff, stepsize;

 75:   tron->pgstepsize = 1.0;
 76:   tao->trust       = tao->trust0;
 77:   /*   Project the current point onto the feasible set */
 78:   TaoComputeVariableBounds(tao);
 79:   TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU);

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

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

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

 93:   /* Initialize trust region radius */
 94:   tao->trust = tao->trust0;
 95:   if (tao->trust <= 0) tao->trust = PetscMax(tron->gnorm * tron->gnorm, 1.0);

 97:   /* Initialize step sizes for the line searches */
 98:   tron->pgstepsize = 1.0;
 99:   tron->stepsize   = tao->trust;

101:   tao->reason = TAO_CONTINUE_ITERATING;
102:   TaoLogConvergenceHistory(tao, tron->f, tron->gnorm, 0.0, tao->ksp_its);
103:   TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, tron->stepsize);
104:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
105:   while (tao->reason == TAO_CONTINUE_ITERATING) {
106:     /* Call general purpose update function */
107:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);

109:     /* Perform projected gradient iterations */
110:     TronGradientProjections(tao, tron);

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

115:     tao->ksp_its      = 0;
116:     f                 = tron->f;
117:     delta             = tao->trust;
118:     tron->n_free_last = tron->n_free;
119:     TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre);

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

126:     /* If no free variables */
127:     if (tron->n_free == 0) {
128:       VecNorm(tao->gradient, NORM_2, &tron->gnorm);
129:       TaoLogConvergenceHistory(tao, tron->f, tron->gnorm, 0.0, tao->ksp_its);
130:       TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta);
131:       PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
132:       if (!tao->reason) tao->reason = TAO_CONVERGED_STEPTOL;
133:       break;
134:     }
135:     /* use free_local to mask/submat gradient, hessian, stepdirection */
136:     TaoVecGetSubVec(tao->gradient, tron->Free_Local, tao->subset_type, 0.0, &tron->R);
137:     TaoVecGetSubVec(tao->gradient, tron->Free_Local, tao->subset_type, 0.0, &tron->DXFree);
138:     VecSet(tron->DXFree, 0.0);
139:     VecScale(tron->R, -1.0);
140:     TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub);
141:     if (tao->hessian == tao->hessian_pre) {
142:       MatDestroy(&tron->Hpre_sub);
143:       PetscObjectReference((PetscObject)(tron->H_sub));
144:       tron->Hpre_sub = tron->H_sub;
145:     } else {
146:       TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type, &tron->Hpre_sub);
147:     }
148:     KSPReset(tao->ksp);
149:     KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub);
150:     while (1) {
151:       /* Approximately solve the reduced linear system */
152:       KSPCGSetRadius(tao->ksp, delta);

154:       KSPSolve(tao->ksp, tron->R, tron->DXFree);
155:       KSPGetIterationNumber(tao->ksp, &its);
156:       tao->ksp_its += its;
157:       tao->ksp_tot_its += its;
158:       VecSet(tao->stepdirection, 0.0);

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

163:       VecDot(tao->gradient, tao->stepdirection, &gdx);
164:       VecCopy(tao->solution, tron->X_New);
165:       VecCopy(tao->gradient, tron->G_New);

167:       stepsize = 1.0;
168:       f_new    = f;

170:       TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
171:       TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection, &stepsize, &ls_reason);
172:       TaoAddLineSearchCounts(tao);

174:       MatMult(tao->hessian, tao->stepdirection, tron->Work);
175:       VecAYPX(tron->Work, 0.5, tao->gradient);
176:       VecDot(tao->stepdirection, tron->Work, &prered);
177:       actred = f_new - f;
178:       if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) {
179:         rhok = 1.0;
180:       } else if (actred < 0) {
181:         rhok = PetscAbs(-actred / prered);
182:       } else {
183:         rhok = 0.0;
184:       }

186:       /* Compare actual improvement to the quadratic model */
187:       if (rhok > tron->eta1) { /* Accept the point */
188:         /* d = x_new - x */
189:         VecCopy(tron->X_New, tao->stepdirection);
190:         VecAXPY(tao->stepdirection, -1.0, tao->solution);

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

195:         /* Adjust trust region size */
196:         if (rhok < tron->eta2) {
197:           delta = PetscMin(xdiff, delta) * tron->sigma1;
198:         } else if (rhok > tron->eta4) {
199:           delta = PetscMin(xdiff, delta) * tron->sigma3;
200:         } else if (rhok > tron->eta3) {
201:           delta = PetscMin(xdiff, delta) * tron->sigma2;
202:         }
203:         VecBoundGradientProjection(tron->G_New, tron->X_New, tao->XL, tao->XU, tao->gradient);
204:         ISDestroy(&tron->Free_Local);
205:         VecWhichInactive(tao->XL, tron->X_New, tao->gradient, tao->XU, PETSC_TRUE, &tron->Free_Local);
206:         f = f_new;
207:         VecNorm(tao->gradient, NORM_2, &tron->gnorm);
208:         VecCopy(tron->X_New, tao->solution);
209:         VecCopy(tron->G_New, tao->gradient);
210:         break;
211:       } else if (delta <= 1e-30) {
212:         break;
213:       } else {
214:         delta /= 4.0;
215:       }
216:     } /* end linear solve loop */

218:     tron->f      = f;
219:     tron->actred = actred;
220:     tao->trust   = delta;
221:     tao->niter++;
222:     TaoLogConvergenceHistory(tao, tron->f, tron->gnorm, 0.0, tao->ksp_its);
223:     TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, stepsize);
224:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
225:   } /* END MAIN LOOP  */
226:   return 0;
227: }

229: static PetscErrorCode TronGradientProjections(Tao tao, TAO_TRON *tron)
230: {
231:   PetscInt                     i;
232:   TaoLineSearchConvergedReason ls_reason;
233:   PetscReal                    actred = -1.0, actred_max = 0.0;
234:   PetscReal                    f_new;
235:   /*
236:      The gradient and function value passed into and out of this
237:      routine should be current and correct.

239:      The free, active, and binding variables should be already identified
240:   */

242:   for (i = 0; i < tron->maxgpits; ++i) {
243:     if (-actred <= (tron->pg_ftol) * actred_max) break;

245:     ++tron->gp_iterates;
246:     ++tron->total_gp_its;
247:     f_new = tron->f;

249:     VecCopy(tao->gradient, tao->stepdirection);
250:     VecScale(tao->stepdirection, -1.0);
251:     TaoLineSearchSetInitialStepLength(tao->linesearch, tron->pgstepsize);
252:     TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, &tron->pgstepsize, &ls_reason);
253:     TaoAddLineSearchCounts(tao);

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

258:     /* Update the iterate */
259:     actred     = f_new - tron->f;
260:     actred_max = PetscMax(actred_max, -(f_new - tron->f));
261:     tron->f    = f_new;
262:   }
263:   return 0;
264: }

266: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU)
267: {
268:   TAO_TRON *tron = (TAO_TRON *)tao->data;


275:   VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, tron->Work);
276:   VecCopy(tron->Work, DXL);
277:   VecAXPY(DXL, -1.0, tao->gradient);
278:   VecSet(DXU, 0.0);
279:   VecPointwiseMax(DXL, DXL, DXU);

281:   VecCopy(tao->gradient, DXU);
282:   VecAXPY(DXU, -1.0, tron->Work);
283:   VecSet(tron->Work, 0.0);
284:   VecPointwiseMin(DXU, tron->Work, DXU);
285:   return 0;
286: }

288: /*------------------------------------------------------------*/
289: /*MC
290:   TAOTRON - The TRON algorithm is an active-set Newton trust region method
291:   for bound-constrained minimization.

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

297:   Level: beginner
298: M*/
299: PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
300: {
301:   TAO_TRON   *tron;
302:   const char *morethuente_type = TAOLINESEARCHMT;

304:   tao->ops->setup          = TaoSetup_TRON;
305:   tao->ops->solve          = TaoSolve_TRON;
306:   tao->ops->view           = TaoView_TRON;
307:   tao->ops->setfromoptions = TaoSetFromOptions_TRON;
308:   tao->ops->destroy        = TaoDestroy_TRON;
309:   tao->ops->computedual    = TaoComputeDual_TRON;

311:   PetscNew(&tron);
312:   tao->data = (void *)tron;

314:   /* Override default settings (unless already changed) */
315:   if (!tao->max_it_changed) tao->max_it = 50;
316:   if (!tao->trust0_changed) tao->trust0 = 1.0;
317:   if (!tao->steptol_changed) tao->steptol = 0.0;

319:   /* Initialize pointers and variables */
320:   tron->n        = 0;
321:   tron->maxgpits = 3;
322:   tron->pg_ftol  = 0.001;

324:   tron->eta1 = 1.0e-4;
325:   tron->eta2 = 0.25;
326:   tron->eta3 = 0.50;
327:   tron->eta4 = 0.90;

329:   tron->sigma1 = 0.5;
330:   tron->sigma2 = 2.0;
331:   tron->sigma3 = 4.0;

333:   tron->gp_iterates  = 0; /* Cumulative number */
334:   tron->total_gp_its = 0;
335:   tron->n_free       = 0;

337:   tron->DXFree     = NULL;
338:   tron->R          = NULL;
339:   tron->X_New      = NULL;
340:   tron->G_New      = NULL;
341:   tron->Work       = NULL;
342:   tron->Free_Local = NULL;
343:   tron->H_sub      = NULL;
344:   tron->Hpre_sub   = NULL;
345:   tao->subset_type = TAO_SUBSET_SUBVEC;

347:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
348:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
349:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
350:   TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
351:   TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);

353:   KSPCreate(((PetscObject)tao)->comm, &tao->ksp);
354:   PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1);
355:   KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);
356:   KSPSetType(tao->ksp, KSPSTCG);
357:   return 0;
358: }