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: PetscCall(TaoParametersInitialize(tao));
243: PetscObjectParameterSetDefault(tao, max_it, 2000);
244: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
246: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
247: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
248: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
249: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
251: PetscCall(KSPInitializePackage());
252: PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M));
253: PetscCall(MatSetType(blmP->M, MATLMVMBFGS));
254: PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: /*@
259: TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent `TaoSolve()` calls.
261: Input Parameters:
262: + tao - the `Tao` solver context
263: - flg - Boolean flag for recycling (`PETSC_TRUE` or `PETSC_FALSE`)
265: Level: intermediate
267: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`
268: @*/
269: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
270: {
271: TAO_LMVM *lmP;
272: TAO_BLMVM *blmP;
273: PetscBool is_lmvm, is_blmvm;
275: PetscFunctionBegin;
276: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
277: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
278: if (is_lmvm) {
279: lmP = (TAO_LMVM *)tao->data;
280: lmP->recycle = flg;
281: } else if (is_blmvm) {
282: blmP = (TAO_BLMVM *)tao->data;
283: blmP->recycle = flg;
284: }
285: PetscFunctionReturn(PETSC_SUCCESS);
286: }
288: /*@
289: TaoLMVMSetH0 - Set the initial Hessian for the QN approximation
291: Input Parameters:
292: + tao - the `Tao` solver context
293: - H0 - `Mat` object for the initial Hessian
295: Level: advanced
297: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
298: @*/
299: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
300: {
301: TAO_LMVM *lmP;
302: TAO_BLMVM *blmP;
303: PetscBool is_lmvm, is_blmvm;
305: PetscFunctionBegin;
306: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
307: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
308: if (is_lmvm) {
309: lmP = (TAO_LMVM *)tao->data;
310: PetscCall(PetscObjectReference((PetscObject)H0));
311: lmP->H0 = H0;
312: } else if (is_blmvm) {
313: blmP = (TAO_BLMVM *)tao->data;
314: PetscCall(PetscObjectReference((PetscObject)H0));
315: blmP->H0 = H0;
316: }
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: /*@
321: TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian
323: Input Parameter:
324: . tao - the `Tao` solver context
326: Output Parameter:
327: . H0 - `Mat` object for the initial Hessian
329: Level: advanced
331: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
332: @*/
333: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
334: {
335: TAO_LMVM *lmP;
336: TAO_BLMVM *blmP;
337: PetscBool is_lmvm, is_blmvm;
338: Mat M;
340: PetscFunctionBegin;
341: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
342: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
343: if (is_lmvm) {
344: lmP = (TAO_LMVM *)tao->data;
345: M = lmP->M;
346: } else if (is_blmvm) {
347: blmP = (TAO_BLMVM *)tao->data;
348: M = blmP->M;
349: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
350: PetscCall(MatLMVMGetJ0(M, H0));
351: PetscFunctionReturn(PETSC_SUCCESS);
352: }
354: /*@
355: TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian
357: Input Parameter:
358: . tao - the `Tao` solver context
360: Output Parameter:
361: . ksp - `KSP` solver context for the initial Hessian
363: Level: advanced
365: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`
366: @*/
367: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
368: {
369: TAO_LMVM *lmP;
370: TAO_BLMVM *blmP;
371: PetscBool is_lmvm, is_blmvm;
372: Mat M;
374: PetscFunctionBegin;
375: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
376: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
377: if (is_lmvm) {
378: lmP = (TAO_LMVM *)tao->data;
379: M = lmP->M;
380: } else if (is_blmvm) {
381: blmP = (TAO_BLMVM *)tao->data;
382: M = blmP->M;
383: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
384: PetscCall(MatLMVMGetJ0KSP(M, ksp));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }