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: static PetscErrorCode TaoSolve_BLMVM(Tao tao)
6: {
7: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
8: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
9: PetscReal f, fold, gdx, gnorm, gnorm2;
10: PetscReal stepsize = 1.0, delta;
12: PetscFunctionBegin;
13: /* Project initial point onto bounds */
14: PetscCall(TaoComputeVariableBounds(tao));
15: PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution));
16: PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));
18: /* Check convergence criteria */
19: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, blmP->unprojected_gradient));
20: PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
22: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
23: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
25: tao->reason = TAO_CONTINUE_ITERATING;
26: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
27: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize));
28: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
29: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
31: /* Set counter for gradient/reset steps */
32: if (!blmP->recycle) {
33: blmP->grad = 0;
34: blmP->reset = 0;
35: PetscCall(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: PetscCall(TaoComputeObjective(tao, tao->solution, &f));
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: PetscCall(MatLMVMSymBroydenSetDelta(blmP->M, delta));
54: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, tao->gradient));
55: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
56: PetscCall(VecBoundGradientProjection(tao->stepdirection, tao->solution, tao->XL, tao->XU, tao->gradient));
58: /* Check for success (descent direction) */
59: PetscCall(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: PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
66: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient));
67: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
68: }
69: PetscCall(VecScale(tao->stepdirection, -1.0));
71: /* Perform the linesearch */
72: fold = f;
73: PetscCall(VecCopy(tao->solution, blmP->Xold));
74: PetscCall(VecCopy(blmP->unprojected_gradient, blmP->Gold));
75: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
76: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status));
77: PetscCall(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: PetscCall(VecCopy(blmP->Xold, tao->solution));
86: PetscCall(VecCopy(blmP->Gold, blmP->unprojected_gradient));
88: PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE));
89: PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient));
90: PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection));
91: PetscCall(VecScale(tao->stepdirection, -1.0));
93: /* This may be incorrect; linesearch has values for stepmax and stepmin
94: that should be reset. */
95: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
96: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status));
97: PetscCall(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: PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient));
107: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
108: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number");
109: tao->niter++;
110: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
111: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize));
112: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
113: }
114: PetscFunctionReturn(PETSC_SUCCESS);
115: }
117: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
118: {
119: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
121: PetscFunctionBegin;
122: /* Existence of tao->solution checked in TaoSetup() */
123: PetscCall(VecDuplicate(tao->solution, &blmP->Xold));
124: PetscCall(VecDuplicate(tao->solution, &blmP->Gold));
125: PetscCall(VecDuplicate(tao->solution, &blmP->unprojected_gradient));
126: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
127: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
128: /* Allocate matrix for the limited memory approximation */
129: PetscCall(MatLMVMAllocate(blmP->M, tao->solution, blmP->unprojected_gradient));
131: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
132: if (blmP->H0) PetscCall(MatLMVMSetJ0(blmP->M, blmP->H0));
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
137: {
138: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
140: PetscFunctionBegin;
141: if (tao->setupcalled) {
142: PetscCall(VecDestroy(&blmP->unprojected_gradient));
143: PetscCall(VecDestroy(&blmP->Xold));
144: PetscCall(VecDestroy(&blmP->Gold));
145: }
146: PetscCall(MatDestroy(&blmP->M));
147: PetscCall(PetscObjectDereference((PetscObject)blmP->H0));
148: PetscCall(PetscFree(tao->data));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
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: PetscFunctionBegin;
158: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for bound constrained optimization");
159: PetscCall(PetscOptionsBool("-tao_blmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", blmP->recycle, &blmP->recycle, NULL));
160: PetscOptionsHeadEnd();
161: PetscCall(MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix));
162: PetscCall(MatAppendOptionsPrefix(blmP->M, "tao_blmvm_"));
163: PetscCall(MatSetFromOptions(blmP->M));
164: PetscCall(MatIsSPDKnown(blmP->M, &is_set, &is_spd));
165: PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite");
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer)
170: {
171: TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
172: PetscBool isascii;
174: PetscFunctionBegin;
175: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
176: if (isascii) {
177: PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad));
178: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
179: PetscCall(MatView(lmP->M, viewer));
180: PetscCall(PetscViewerPopFormat(viewer));
181: }
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
186: {
187: TAO_BLMVM *blm = (TAO_BLMVM *)tao->data;
189: PetscFunctionBegin;
193: PetscCheck(tao->gradient && blm->unprojected_gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist.");
195: PetscCall(VecCopy(tao->gradient, DXL));
196: PetscCall(VecAXPY(DXL, -1.0, blm->unprojected_gradient));
197: PetscCall(VecSet(DXU, 0.0));
198: PetscCall(VecPointwiseMax(DXL, DXL, DXU));
200: PetscCall(VecCopy(blm->unprojected_gradient, DXU));
201: PetscCall(VecAXPY(DXU, -1.0, tao->gradient));
202: PetscCall(VecAXPY(DXU, 1.0, DXL));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*MC
207: TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
208: for nonlinear minimization with bound constraints. It is an extension
209: of `TAOLMVM`
211: Options Database Key:
212: . -tao_lmm_recycle - enable recycling of LMVM information between subsequent `TaoSolve()` calls
214: Level: beginner
216: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
217: M*/
218: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
219: {
220: TAO_BLMVM *blmP;
221: const char *morethuente_type = TAOLINESEARCHMT;
223: PetscFunctionBegin;
224: tao->ops->setup = TaoSetup_BLMVM;
225: tao->ops->solve = TaoSolve_BLMVM;
226: tao->ops->view = TaoView_BLMVM;
227: tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
228: tao->ops->destroy = TaoDestroy_BLMVM;
229: tao->ops->computedual = TaoComputeDual_BLMVM;
230: tao->uses_gradient = PETSC_TRUE;
232: PetscCall(PetscNew(&blmP));
233: tao->data = (void *)blmP;
235: /* Override default settings (unless already changed) */
236: PetscObjectParameterSetDefault(tao, max_it, 2000);
237: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
239: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
240: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
241: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
242: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
244: PetscCall(KSPInitializePackage());
245: PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M));
246: PetscCall(MatSetType(blmP->M, MATLMVMBFGS));
247: PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1));
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: /*@
252: TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent `TaoSolve()` calls.
254: Input Parameters:
255: + tao - the `Tao` solver context
256: - flg - Boolean flag for recycling (`PETSC_TRUE` or `PETSC_FALSE`)
258: Level: intermediate
260: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`
261: @*/
262: PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg)
263: {
264: TAO_LMVM *lmP;
265: TAO_BLMVM *blmP;
266: PetscBool is_lmvm, is_blmvm;
268: PetscFunctionBegin;
269: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
270: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
271: if (is_lmvm) {
272: lmP = (TAO_LMVM *)tao->data;
273: lmP->recycle = flg;
274: } else if (is_blmvm) {
275: blmP = (TAO_BLMVM *)tao->data;
276: blmP->recycle = flg;
277: }
278: PetscFunctionReturn(PETSC_SUCCESS);
279: }
281: /*@
282: TaoLMVMSetH0 - Set the initial Hessian for the QN approximation
284: Input Parameters:
285: + tao - the `Tao` solver context
286: - H0 - `Mat` object for the initial Hessian
288: Level: advanced
290: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()`
291: @*/
292: PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
293: {
294: TAO_LMVM *lmP;
295: TAO_BLMVM *blmP;
296: PetscBool is_lmvm, is_blmvm;
298: PetscFunctionBegin;
299: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
300: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
301: if (is_lmvm) {
302: lmP = (TAO_LMVM *)tao->data;
303: PetscCall(PetscObjectReference((PetscObject)H0));
304: lmP->H0 = H0;
305: } else if (is_blmvm) {
306: blmP = (TAO_BLMVM *)tao->data;
307: PetscCall(PetscObjectReference((PetscObject)H0));
308: blmP->H0 = H0;
309: }
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: /*@
314: TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian
316: Input Parameter:
317: . tao - the `Tao` solver context
319: Output Parameter:
320: . H0 - `Mat` object for the initial Hessian
322: Level: advanced
324: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()`
325: @*/
326: PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
327: {
328: TAO_LMVM *lmP;
329: TAO_BLMVM *blmP;
330: PetscBool is_lmvm, is_blmvm;
331: Mat M;
333: PetscFunctionBegin;
334: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
335: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
336: if (is_lmvm) {
337: lmP = (TAO_LMVM *)tao->data;
338: M = lmP->M;
339: } else if (is_blmvm) {
340: blmP = (TAO_BLMVM *)tao->data;
341: M = blmP->M;
342: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
343: PetscCall(MatLMVMGetJ0(M, H0));
344: PetscFunctionReturn(PETSC_SUCCESS);
345: }
347: /*@
348: TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian
350: Input Parameter:
351: . tao - the `Tao` solver context
353: Output Parameter:
354: . ksp - `KSP` solver context for the initial Hessian
356: Level: advanced
358: .seealso: `Tao`, `TAOLMVM`, `TAOBLMVM`, `TaoLMVMGetH0()`
359: @*/
360: PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
361: {
362: TAO_LMVM *lmP;
363: TAO_BLMVM *blmP;
364: PetscBool is_lmvm, is_blmvm;
365: Mat M;
367: PetscFunctionBegin;
368: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm));
369: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm));
370: if (is_lmvm) {
371: lmP = (TAO_LMVM *)tao->data;
372: M = lmP->M;
373: } else if (is_blmvm) {
374: blmP = (TAO_BLMVM *)tao->data;
375: M = blmP->M;
376: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM.");
377: PetscCall(MatLMVMGetJ0KSP(M, ksp));
378: PetscFunctionReturn(PETSC_SUCCESS);
379: }