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:   PetscFunctionBegin;
 14:   /*  Project initial point onto bounds */
 15:   PetscCall(TaoComputeVariableBounds(tao));
 16:   PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
 17:   PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));

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

 23:   PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
 24:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

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

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

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

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

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

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

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

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

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

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

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

106:     /* Check for converged */
107:     PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
108:     PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
109:     PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
110:     tao->niter++;
111:     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
112:     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize));
113:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
114:   }
115:   PetscFunctionReturn(PETSC_SUCCESS);
116: }

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

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

132:   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
133:   if (blmP->H0) PetscCall(MatLMVMSetJ0(blmP->M, blmP->H0));
134:   PetscFunctionReturn(PETSC_SUCCESS);
135: }

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

142:   PetscFunctionBegin;
143:   if (tao->setupcalled) {
144:     PetscCall(VecDestroy(&blmP->unprojected_gradient));
145:     PetscCall(VecDestroy(&blmP->Xold));
146:     PetscCall(VecDestroy(&blmP->Gold));
147:   }
148:   PetscCall(MatDestroy(&blmP->M));
149:   if (blmP->H0) PetscCall(PetscObjectDereference((PetscObject)blmP->H0));
150:   PetscCall(PetscFree(tao->data));
151:   PetscFunctionReturn(PETSC_SUCCESS);
152: }

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

160:   PetscFunctionBegin;
161:   PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for bound constrained optimization");
162:   PetscCall(PetscOptionsBool("-tao_blmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", blmP->recycle, &blmP->recycle, NULL));
163:   PetscOptionsHeadEnd();
164:   PetscCall(MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix));
165:   PetscCall(MatAppendOptionsPrefix(blmP->M, "tao_blmvm_"));
166:   PetscCall(MatSetFromOptions(blmP->M));
167:   PetscCall(MatIsSPDKnown(blmP->M, &is_set, &is_spd));
168:   PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite");
169:   PetscFunctionReturn(PETSC_SUCCESS);
170: }

172: /*------------------------------------------------------------*/
173: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
174: {
175:   TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
176:   PetscBool  isascii;

178:   PetscFunctionBegin;
179:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
180:   if (isascii) {
181:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad));
182:     PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
183:     PetscCall(MatView(lmP->M, viewer));
184:     PetscCall(PetscViewerPopFormat(viewer));
185:   }
186:   PetscFunctionReturn(PETSC_SUCCESS);
187: }

189: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
190: {
191:   TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;

193:   PetscFunctionBegin;
197:   PetscCheck(tao->gradient && blm->unprojected_gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist.");

199:   PetscCall(VecCopy(tao->gradient, DXL));
200:   PetscCall(VecAXPY(DXL, -1.0, blm->unprojected_gradient));
201:   PetscCall(VecSet(DXU, 0.0));
202:   PetscCall(VecPointwiseMax(DXL, DXL, DXU));

204:   PetscCall(VecCopy(blm->unprojected_gradient, DXU));
205:   PetscCall(VecAXPY(DXU, -1.0, tao->gradient));
206:   PetscCall(VecAXPY(DXU, 1.0, DXL));
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }

210: /* ---------------------------------------------------------- */
211: /*MC
212:   TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
213:          for nonlinear minimization with bound constraints. It is an extension
214:          of `TAOLMVM`

216:   Options Database Key:
217: .     -tao_lmm_recycle - enable recycling of LMVM information between subsequent `TaoSolve()` calls

219:   Level: beginner

221: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
222: M*/
223: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
224: {
225:   TAO_BLMVM  *blmP;
226:   const char *morethuente_type = TAOLINESEARCHMT;

228:   PetscFunctionBegin;
229:   tao->ops->setup          = TaoSetup_BLMVM;
230:   tao->ops->solve          = TaoSolve_BLMVM;
231:   tao->ops->view           = TaoView_BLMVM;
232:   tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
233:   tao->ops->destroy        = TaoDestroy_BLMVM;
234:   tao->ops->computedual    = TaoComputeDual_BLMVM;

236:   PetscCall(PetscNew(&blmP));
237:   blmP->H0      = NULL;
238:   blmP->recycle = PETSC_FALSE;
239:   tao->data     = (void *)blmP;

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

245:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
246:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
247:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
248:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));

250:   PetscCall(KSPInitializePackage());
251:   PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M));
252:   PetscCall(MatSetType(blmP->M, MATLMVMBFGS));
253:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

257: /*@
258:   TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent `TaoSolve()` calls.

260:   Input Parameters:
261: + tao - the `Tao` solver context
262: - flg - Boolean flag for recycling (`PETSC_TRUE` or `PETSC_FALSE`)

264:   Level: intermediate

266: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`
267: @*/
268: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
269: {
270:   TAO_LMVM  *lmP;
271:   TAO_BLMVM *blmP;
272:   PetscBool  is_lmvm, is_blmvm;

274:   PetscFunctionBegin;
275:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
276:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
277:   if (is_lmvm) {
278:     lmP          = (TAO_LMVM *)tao->data;
279:     lmP->recycle = flg;
280:   } else if (is_blmvm) {
281:     blmP          = (TAO_BLMVM *)tao->data;
282:     blmP->recycle = flg;
283:   }
284:   PetscFunctionReturn(PETSC_SUCCESS);
285: }

287: /*@
288:   TaoLMVMSetH0 - Set the initial Hessian for the QN approximation

290:   Input Parameters:
291: + tao - the `Tao` solver context
292: - H0  - `Mat` object for the initial Hessian

294:   Level: advanced

296: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
297: @*/
298: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
299: {
300:   TAO_LMVM  *lmP;
301:   TAO_BLMVM *blmP;
302:   PetscBool  is_lmvm, is_blmvm;

304:   PetscFunctionBegin;
305:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
306:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
307:   if (is_lmvm) {
308:     lmP = (TAO_LMVM *)tao->data;
309:     PetscCall(PetscObjectReference((PetscObject)H0));
310:     lmP->H0 = H0;
311:   } else if (is_blmvm) {
312:     blmP = (TAO_BLMVM *)tao->data;
313:     PetscCall(PetscObjectReference((PetscObject)H0));
314:     blmP->H0 = H0;
315:   }
316:   PetscFunctionReturn(PETSC_SUCCESS);
317: }

319: /*@
320:   TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian

322:   Input Parameter:
323: . tao - the `Tao` solver context

325:   Output Parameter:
326: . H0 - `Mat` object for the initial Hessian

328:   Level: advanced

330: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
331: @*/
332: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
333: {
334:   TAO_LMVM  *lmP;
335:   TAO_BLMVM *blmP;
336:   PetscBool  is_lmvm, is_blmvm;
337:   Mat        M;

339:   PetscFunctionBegin;
340:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
341:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
342:   if (is_lmvm) {
343:     lmP = (TAO_LMVM *)tao->data;
344:     M   = lmP->M;
345:   } else if (is_blmvm) {
346:     blmP = (TAO_BLMVM *)tao->data;
347:     M    = blmP->M;
348:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
349:   PetscCall(MatLMVMGetJ0(M, H0));
350:   PetscFunctionReturn(PETSC_SUCCESS);
351: }

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

356:   Input Parameter:
357: . tao - the `Tao` solver context

359:   Output Parameter:
360: . ksp - `KSP` solver context for the initial Hessian

362:   Level: advanced

364: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`
365: @*/
366: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
367: {
368:   TAO_LMVM  *lmP;
369:   TAO_BLMVM *blmP;
370:   PetscBool  is_lmvm, is_blmvm;
371:   Mat        M;

373:   PetscFunctionBegin;
374:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
375:   PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
376:   if (is_lmvm) {
377:     lmP = (TAO_LMVM *)tao->data;
378:     M   = lmP->M;
379:   } else if (is_blmvm) {
380:     blmP = (TAO_BLMVM *)tao->data;
381:     M    = blmP->M;
382:   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
383:   PetscCall(MatLMVMGetJ0KSP(M, ksp));
384:   PetscFunctionReturn(PETSC_SUCCESS);
385: }