Actual source code: blmvm.c

  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>
  3: #include <../src/tao/bound/impls/blmvm/blmvm.h>

  5: /*------------------------------------------------------------*/
  6: static PetscErrorCode TaoSolve_BLMVM(Tao tao)
  7: {
  8:   TAO_BLMVM                   *blmP      = (TAO_BLMVM *)tao->data;
  9:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
 10:   PetscReal                    f, fold, gdx, gnorm, gnorm2;
 11:   PetscReal                    stepsize = 1.0, delta;

 13:   /*  Project initial point onto bounds */
 14:   TaoComputeVariableBounds(tao);
 15:   VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);
 16:   TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU);

 18:   /* Check convergence criteria */
 19:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, blmP->unprojected_gradient);
 20:   VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);

 22:   TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);

 25:   tao->reason = TAO_CONTINUE_ITERATING;
 26:   TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
 27:   TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize);
 28:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 29:   if (tao->reason != TAO_CONTINUE_ITERATING) return 0;

 31:   /* Set counter for gradient/reset steps */
 32:   if (!blmP->recycle) {
 33:     blmP->grad  = 0;
 34:     blmP->reset = 0;
 35:     MatLMVMReset(blmP->M, PETSC_FALSE);
 36:   }

 38:   /* Have not converged; continue with Newton method */
 39:   while (tao->reason == TAO_CONTINUE_ITERATING) {
 40:     /* Call general purpose update function */
 41:     if (tao->ops->update) {
 42:       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
 43:       TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 44:     }
 45:     /* Compute direction */
 46:     gnorm2 = gnorm * gnorm;
 47:     if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON;
 48:     if (f == 0.0) {
 49:       delta = 2.0 / gnorm2;
 50:     } else {
 51:       delta = 2.0 * PetscAbsScalar(f) / gnorm2;
 52:     }
 53:     MatLMVMSymBroydenSetDelta(blmP->M, delta);
 54:     MatLMVMUpdate(blmP->M, tao->solution, tao->gradient);
 55:     MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
 56:     VecBoundGradientProjection(tao->stepdirection, tao->solution, tao->XL, tao->XU, tao->gradient);

 58:     /* Check for success (descent direction) */
 59:     VecDot(blmP->unprojected_gradient, tao->gradient, &gdx);
 60:     if (gdx <= 0) {
 61:       /* Step is not descent or solve was not successful
 62:          Use steepest descent direction (scaled) */
 63:       ++blmP->grad;

 65:       MatLMVMReset(blmP->M, PETSC_FALSE);
 66:       MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
 67:       MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
 68:     }
 69:     VecScale(tao->stepdirection, -1.0);

 71:     /* Perform the linesearch */
 72:     fold = f;
 73:     VecCopy(tao->solution, blmP->Xold);
 74:     VecCopy(blmP->unprojected_gradient, blmP->Gold);
 75:     TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
 76:     TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);
 77:     TaoAddLineSearchCounts(tao);

 79:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
 80:       /* Linesearch failed
 81:          Reset factors and use scaled (projected) gradient step */
 82:       ++blmP->reset;

 84:       f = fold;
 85:       VecCopy(blmP->Xold, tao->solution);
 86:       VecCopy(blmP->Gold, blmP->unprojected_gradient);

 88:       MatLMVMReset(blmP->M, PETSC_FALSE);
 89:       MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
 90:       MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
 91:       VecScale(tao->stepdirection, -1.0);

 93:       /* This may be incorrect; linesearch has values for stepmax and stepmin
 94:          that should be reset. */
 95:       TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0);
 96:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);
 97:       TaoAddLineSearchCounts(tao);

 99:       if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
100:         tao->reason = TAO_DIVERGED_LS_FAILURE;
101:         break;
102:       }
103:     }

105:     /* Check for converged */
106:     VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);
107:     TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
109:     tao->niter++;
110:     TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
111:     TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize);
112:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
113:   }
114:   return 0;
115: }

117: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
118: {
119:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;

121:   /* Existence of tao->solution checked in TaoSetup() */
122:   VecDuplicate(tao->solution, &blmP->Xold);
123:   VecDuplicate(tao->solution, &blmP->Gold);
124:   VecDuplicate(tao->solution, &blmP->unprojected_gradient);
125:   if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
126:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
127:   /* Allocate matrix for the limited memory approximation */
128:   MatLMVMAllocate(blmP->M, tao->solution, blmP->unprojected_gradient);

130:   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
131:   if (blmP->H0) MatLMVMSetJ0(blmP->M, blmP->H0);
132:   return 0;
133: }

135: /* ---------------------------------------------------------- */
136: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
137: {
138:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;

140:   if (tao->setupcalled) {
141:     VecDestroy(&blmP->unprojected_gradient);
142:     VecDestroy(&blmP->Xold);
143:     VecDestroy(&blmP->Gold);
144:   }
145:   MatDestroy(&blmP->M);
146:   if (blmP->H0) PetscObjectDereference((PetscObject)blmP->H0);
147:   PetscFree(tao->data);
148:   return 0;
149: }

151: /*------------------------------------------------------------*/
152: static PetscErrorCode TaoSetFromOptions_BLMVM(Tao tao, PetscOptionItems *PetscOptionsObject)
153: {
154:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
155:   PetscBool  is_spd, is_set;

157:   PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for bound constrained optimization");
158:   PetscOptionsBool("-tao_blmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", blmP->recycle, &blmP->recycle, NULL);
159:   PetscOptionsHeadEnd();
160:   MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix);
161:   MatAppendOptionsPrefix(blmP->M, "tao_blmvm_");
162:   MatSetFromOptions(blmP->M);
163:   MatIsSPDKnown(blmP->M, &is_set, &is_spd);
165:   return 0;
166: }

168: /*------------------------------------------------------------*/
169: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
170: {
171:   TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
172:   PetscBool  isascii;

174:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
175:   if (isascii) {
176:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad);
177:     PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO);
178:     MatView(lmP->M, viewer);
179:     PetscViewerPopFormat(viewer);
180:   }
181:   return 0;
182: }

184: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
185: {
186:   TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;


193:   VecCopy(tao->gradient, DXL);
194:   VecAXPY(DXL, -1.0, blm->unprojected_gradient);
195:   VecSet(DXU, 0.0);
196:   VecPointwiseMax(DXL, DXL, DXU);

198:   VecCopy(blm->unprojected_gradient, DXU);
199:   VecAXPY(DXU, -1.0, tao->gradient);
200:   VecAXPY(DXU, 1.0, DXL);
201:   return 0;
202: }

204: /* ---------------------------------------------------------- */
205: /*MC
206:   TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
207:          for nonlinear minimization with bound constraints. It is an extension
208:          of TAOLMVM

210:   Options Database Keys:
211: .     -tao_lmm_recycle - enable recycling of LMVM information between subsequent TaoSolve calls

213:   Level: beginner
214: M*/
215: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
216: {
217:   TAO_BLMVM  *blmP;
218:   const char *morethuente_type = TAOLINESEARCHMT;

220:   tao->ops->setup          = TaoSetup_BLMVM;
221:   tao->ops->solve          = TaoSolve_BLMVM;
222:   tao->ops->view           = TaoView_BLMVM;
223:   tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
224:   tao->ops->destroy        = TaoDestroy_BLMVM;
225:   tao->ops->computedual    = TaoComputeDual_BLMVM;

227:   PetscNew(&blmP);
228:   blmP->H0      = NULL;
229:   blmP->recycle = PETSC_FALSE;
230:   tao->data     = (void *)blmP;

232:   /* Override default settings (unless already changed) */
233:   if (!tao->max_it_changed) tao->max_it = 2000;
234:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

236:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
237:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
238:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
239:   TaoLineSearchUseTaoRoutines(tao->linesearch, tao);

241:   KSPInitializePackage();
242:   MatCreate(((PetscObject)tao)->comm, &blmP->M);
243:   MatSetType(blmP->M, MATLMVMBFGS);
244:   PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1);
245:   return 0;
246: }

248: /*@
249:   TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent TaoSolve calls.

251:   Input Parameters:
252: +  tao  - the Tao solver context
253: -  flg - Boolean flag for recycling (PETSC_TRUE or PETSC_FALSE)

255:   Level: intermediate
256: @*/
257: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
258: {
259:   TAO_LMVM  *lmP;
260:   TAO_BLMVM *blmP;
261:   PetscBool  is_lmvm, is_blmvm;

263:   PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm);
264:   PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm);
265:   if (is_lmvm) {
266:     lmP          = (TAO_LMVM *)tao->data;
267:     lmP->recycle = flg;
268:   } else if (is_blmvm) {
269:     blmP          = (TAO_BLMVM *)tao->data;
270:     blmP->recycle = flg;
271:   }
272:   return 0;
273: }

275: /*@
276:   TaoLMVMSetH0 - Set the initial Hessian for the QN approximation

278:   Input Parameters:
279: +  tao  - the Tao solver context
280: -  H0 - Mat object for the initial Hessian

282:   Level: advanced

284: .seealso: `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
285: @*/
286: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
287: {
288:   TAO_LMVM  *lmP;
289:   TAO_BLMVM *blmP;
290:   PetscBool  is_lmvm, is_blmvm;

292:   PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm);
293:   PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm);
294:   if (is_lmvm) {
295:     lmP = (TAO_LMVM *)tao->data;
296:     PetscObjectReference((PetscObject)H0);
297:     lmP->H0 = H0;
298:   } else if (is_blmvm) {
299:     blmP = (TAO_BLMVM *)tao->data;
300:     PetscObjectReference((PetscObject)H0);
301:     blmP->H0 = H0;
302:   }
303:   return 0;
304: }

306: /*@
307:   TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian

309:   Input Parameters:
310: .  tao  - the Tao solver context

312:   Output Parameters:
313: .  H0 - Mat object for the initial Hessian

315:   Level: advanced

317: .seealso: `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
318: @*/
319: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
320: {
321:   TAO_LMVM  *lmP;
322:   TAO_BLMVM *blmP;
323:   PetscBool  is_lmvm, is_blmvm;
324:   Mat        M;

326:   PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm);
327:   PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm);
328:   if (is_lmvm) {
329:     lmP = (TAO_LMVM *)tao->data;
330:     M   = lmP->M;
331:   } else if (is_blmvm) {
332:     blmP = (TAO_BLMVM *)tao->data;
333:     M    = blmP->M;
334:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
335:   MatLMVMGetJ0(M, H0);
336:   return 0;
337: }

339: /*@
340:   TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian

342:   Input Parameters:
343: .  tao  - the Tao solver context

345:   Output Parameters:
346: .  ksp - KSP solver context for the initial Hessian

348:   Level: advanced

350: .seealso: `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
351: @*/
352: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
353: {
354:   TAO_LMVM  *lmP;
355:   TAO_BLMVM *blmP;
356:   PetscBool  is_lmvm, is_blmvm;
357:   Mat        M;

359:   PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm);
360:   PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm);
361:   if (is_lmvm) {
362:     lmP = (TAO_LMVM *)tao->data;
363:     M   = lmP->M;
364:   } else if (is_blmvm) {
365:     blmP = (TAO_BLMVM *)tao->data;
366:     M    = blmP->M;
367:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
368:   MatLMVMGetJ0KSP(M, ksp);
369:   return 0;
370: }