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:   PetscFunctionBegin;
 12:   PetscCall(VecDestroy(&tron->X_New));
 13:   PetscCall(VecDestroy(&tron->G_New));
 14:   PetscCall(VecDestroy(&tron->Work));
 15:   PetscCall(VecDestroy(&tron->DXFree));
 16:   PetscCall(VecDestroy(&tron->R));
 17:   PetscCall(VecDestroy(&tron->diag));
 18:   PetscCall(VecScatterDestroy(&tron->scatter));
 19:   PetscCall(ISDestroy(&tron->Free_Local));
 20:   PetscCall(MatDestroy(&tron->H_sub));
 21:   PetscCall(MatDestroy(&tron->Hpre_sub));
 22:   PetscCall(KSPDestroy(&tao->ksp));
 23:   PetscCall(PetscFree(tao->data));
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

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

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

 41: /*------------------------------------------------------------*/
 42: static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer)
 43: {
 44:   TAO_TRON *tron = (TAO_TRON *)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 ",", tron->total_gp_its));
 51:     PetscCall(PetscViewerASCIIPrintf(viewer, "PG tolerance: %g \n", (double)tron->pg_ftol));
 52:   }
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: /* ---------------------------------------------------------- */
 57: static PetscErrorCode TaoSetup_TRON(Tao tao)
 58: {
 59:   TAO_TRON *tron = (TAO_TRON *)tao->data;

 61:   PetscFunctionBegin;
 62:   /* Allocate some arrays */
 63:   PetscCall(VecDuplicate(tao->solution, &tron->diag));
 64:   PetscCall(VecDuplicate(tao->solution, &tron->X_New));
 65:   PetscCall(VecDuplicate(tao->solution, &tron->G_New));
 66:   PetscCall(VecDuplicate(tao->solution, &tron->Work));
 67:   PetscCall(VecDuplicate(tao->solution, &tao->gradient));
 68:   PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
 69:   PetscFunctionReturn(PETSC_SUCCESS);
 70: }

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

 79:   PetscFunctionBegin;
 80:   tron->pgstepsize = 1.0;
 81:   tao->trust       = tao->trust0;
 82:   /*   Project the current point onto the feasible set */
 83:   PetscCall(TaoComputeVariableBounds(tao));
 84:   PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));

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

 89:   /* Compute the objective function and gradient */
 90:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &tron->f, tao->gradient));
 91:   PetscCall(VecNorm(tao->gradient, NORM_2, &tron->gnorm));
 92:   PetscCheck(!PetscIsInfOrNanReal(tron->f) && !PetscIsInfOrNanReal(tron->gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

 94:   /* Project the gradient and calculate the norm */
 95:   PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
 96:   PetscCall(VecNorm(tao->gradient, NORM_2, &tron->gnorm));

 98:   /* Initialize trust region radius */
 99:   tao->trust = tao->trust0;
100:   if (tao->trust <= 0) tao->trust = PetscMax(tron->gnorm * tron->gnorm, 1.0);

102:   /* Initialize step sizes for the line searches */
103:   tron->pgstepsize = 1.0;
104:   tron->stepsize   = tao->trust;

106:   tao->reason = TAO_CONTINUE_ITERATING;
107:   PetscCall(TaoLogConvergenceHistory(tao, tron->f, tron->gnorm, 0.0, tao->ksp_its));
108:   PetscCall(TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, tron->stepsize));
109:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
110:   while (tao->reason == TAO_CONTINUE_ITERATING) {
111:     /* Call general purpose update function */
112:     if (tao->ops->update) {
113:       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
114:       PetscCall(TaoComputeObjective(tao, tao->solution, &tron->f));
115:     }

117:     /* Perform projected gradient iterations */
118:     PetscCall(TronGradientProjections(tao, tron));

120:     PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
121:     PetscCall(VecNorm(tao->gradient, NORM_2, &tron->gnorm));

123:     tao->ksp_its      = 0;
124:     f                 = tron->f;
125:     delta             = tao->trust;
126:     tron->n_free_last = tron->n_free;
127:     PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre));

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

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

162:       PetscCall(KSPSolve(tao->ksp, tron->R, tron->DXFree));
163:       PetscCall(KSPGetIterationNumber(tao->ksp, &its));
164:       tao->ksp_its += its;
165:       tao->ksp_tot_its += its;
166:       PetscCall(VecSet(tao->stepdirection, 0.0));

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

171:       PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx));
172:       PetscCall(VecCopy(tao->solution, tron->X_New));
173:       PetscCall(VecCopy(tao->gradient, tron->G_New));

175:       stepsize = 1.0;
176:       f_new    = f;

178:       PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
179:       PetscCall(TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection, &stepsize, &ls_reason));
180:       PetscCall(TaoAddLineSearchCounts(tao));

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

194:       /* Compare actual improvement to the quadratic model */
195:       if (rhok > tron->eta1) { /* Accept the point */
196:         /* d = x_new - x */
197:         PetscCall(VecCopy(tron->X_New, tao->stepdirection));
198:         PetscCall(VecAXPY(tao->stepdirection, -1.0, tao->solution));

200:         PetscCall(VecNorm(tao->stepdirection, NORM_2, &xdiff));
201:         xdiff *= stepsize;

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

226:     tron->f      = f;
227:     tron->actred = actred;
228:     tao->trust   = delta;
229:     tao->niter++;
230:     PetscCall(TaoLogConvergenceHistory(tao, tron->f, tron->gnorm, 0.0, tao->ksp_its));
231:     PetscCall(TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, stepsize));
232:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
233:   } /* END MAIN LOOP  */
234:   PetscFunctionReturn(PETSC_SUCCESS);
235: }

237: static PetscErrorCode TronGradientProjections(Tao tao, TAO_TRON *tron)
238: {
239:   PetscInt                     i;
240:   TaoLineSearchConvergedReason ls_reason;
241:   PetscReal                    actred = -1.0, actred_max = 0.0;
242:   PetscReal                    f_new;
243:   /*
244:      The gradient and function value passed into and out of this
245:      routine should be current and correct.

247:      The free, active, and binding variables should be already identified
248:   */

250:   PetscFunctionBegin;
251:   for (i = 0; i < tron->maxgpits; ++i) {
252:     if (-actred <= (tron->pg_ftol) * actred_max) break;

254:     ++tron->gp_iterates;
255:     ++tron->total_gp_its;
256:     f_new = tron->f;

258:     PetscCall(VecCopy(tao->gradient, tao->stepdirection));
259:     PetscCall(VecScale(tao->stepdirection, -1.0));
260:     PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, tron->pgstepsize));
261:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, &tron->pgstepsize, &ls_reason));
262:     PetscCall(TaoAddLineSearchCounts(tao));

264:     PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
265:     PetscCall(VecNorm(tao->gradient, NORM_2, &tron->gnorm));

267:     /* Update the iterate */
268:     actred     = f_new - tron->f;
269:     actred_max = PetscMax(actred_max, -(f_new - tron->f));
270:     tron->f    = f_new;
271:   }
272:   PetscFunctionReturn(PETSC_SUCCESS);
273: }

275: static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU)
276: {
277:   TAO_TRON *tron = (TAO_TRON *)tao->data;

279:   PetscFunctionBegin;
283:   PetscCheck(tron->Work && tao->gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist.");

285:   PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, tron->Work));
286:   PetscCall(VecCopy(tron->Work, DXL));
287:   PetscCall(VecAXPY(DXL, -1.0, tao->gradient));
288:   PetscCall(VecSet(DXU, 0.0));
289:   PetscCall(VecPointwiseMax(DXL, DXL, DXU));

291:   PetscCall(VecCopy(tao->gradient, DXU));
292:   PetscCall(VecAXPY(DXU, -1.0, tron->Work));
293:   PetscCall(VecSet(tron->Work, 0.0));
294:   PetscCall(VecPointwiseMin(DXU, tron->Work, DXU));
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: /*------------------------------------------------------------*/
299: /*MC
300:   TAOTRON - The TRON algorithm is an active-set Newton trust region method
301:   for bound-constrained minimization.

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

307:   Level: beginner
308: M*/
309: PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao)
310: {
311:   TAO_TRON   *tron;
312:   const char *morethuente_type = TAOLINESEARCHMT;

314:   PetscFunctionBegin;
315:   tao->ops->setup          = TaoSetup_TRON;
316:   tao->ops->solve          = TaoSolve_TRON;
317:   tao->ops->view           = TaoView_TRON;
318:   tao->ops->setfromoptions = TaoSetFromOptions_TRON;
319:   tao->ops->destroy        = TaoDestroy_TRON;
320:   tao->ops->computedual    = TaoComputeDual_TRON;

322:   PetscCall(PetscNew(&tron));
323:   tao->data = (void *)tron;

325:   /* Override default settings (unless already changed) */
326:   if (!tao->max_it_changed) tao->max_it = 50;
327:   if (!tao->trust0_changed) tao->trust0 = 1.0;
328:   if (!tao->steptol_changed) tao->steptol = 0.0;

330:   /* Initialize pointers and variables */
331:   tron->n        = 0;
332:   tron->maxgpits = 3;
333:   tron->pg_ftol  = 0.001;

335:   tron->eta1 = 1.0e-4;
336:   tron->eta2 = 0.25;
337:   tron->eta3 = 0.50;
338:   tron->eta4 = 0.90;

340:   tron->sigma1 = 0.5;
341:   tron->sigma2 = 2.0;
342:   tron->sigma3 = 4.0;

344:   tron->gp_iterates  = 0; /* Cumulative number */
345:   tron->total_gp_its = 0;
346:   tron->n_free       = 0;

348:   tron->DXFree     = NULL;
349:   tron->R          = NULL;
350:   tron->X_New      = NULL;
351:   tron->G_New      = NULL;
352:   tron->Work       = NULL;
353:   tron->Free_Local = NULL;
354:   tron->H_sub      = NULL;
355:   tron->Hpre_sub   = NULL;
356:   tao->subset_type = TAO_SUBSET_SUBVEC;

358:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
359:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
360:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
361:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
362:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));

364:   PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp));
365:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1));
366:   PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix));
367:   PetscCall(KSPSetType(tao->ksp, KSPSTCG));
368:   PetscFunctionReturn(PETSC_SUCCESS);
369: }