Actual source code: lmvm.c
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>
4: #define LMVM_STEP_BFGS 0
5: #define LMVM_STEP_GRAD 1
7: static PetscErrorCode TaoSolve_LMVM(Tao tao)
8: {
9: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
10: PetscReal f, fold, gdx, gnorm;
11: PetscReal step = 1.0;
12: PetscInt stepType = LMVM_STEP_GRAD, nupdates;
13: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
15: PetscFunctionBegin;
16: if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n"));
18: /* Check convergence criteria */
19: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
20: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
22: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
24: tao->reason = TAO_CONTINUE_ITERATING;
25: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
26: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
27: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
28: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
30: /* Set counter for gradient/reset steps */
31: if (!lmP->recycle) {
32: lmP->bfgs = 0;
33: lmP->grad = 0;
34: PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
35: }
37: /* Have not converged; continue with Newton method */
38: while (tao->reason == TAO_CONTINUE_ITERATING) {
39: /* Call general purpose update function */
40: if (tao->ops->update) {
41: PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
42: PetscCall(TaoComputeObjective(tao, tao->solution, &f));
43: }
45: /* Compute direction */
46: if (lmP->H0) {
47: PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0));
48: stepType = LMVM_STEP_BFGS;
49: }
50: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
51: PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D));
52: PetscCall(MatLMVMGetUpdateCount(lmP->M, &nupdates));
53: if (nupdates > 0) stepType = LMVM_STEP_BFGS;
55: /* Check for success (descent direction) */
56: PetscCall(VecDotRealPart(lmP->D, tao->gradient, &gdx));
57: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
58: /* Step is not descent or direction produced not a number
59: We can assert bfgsUpdates > 1 in this case because
60: the first solve produces the scaled gradient direction,
61: which is guaranteed to be descent
63: Use steepest descent direction (scaled)
64: */
66: PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
67: PetscCall(MatLMVMClearJ0(lmP->M));
68: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
69: PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D));
71: /* On a reset, the direction cannot be not a number; it is a
72: scaled gradient step. No need to check for this condition. */
73: stepType = LMVM_STEP_GRAD;
74: }
75: PetscCall(VecScale(lmP->D, -1.0));
77: /* Perform the linesearch */
78: fold = f;
79: PetscCall(VecCopy(tao->solution, lmP->Xold));
80: PetscCall(VecCopy(tao->gradient, lmP->Gold));
82: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status));
83: PetscCall(TaoAddLineSearchCounts(tao));
85: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) {
86: /* Reset factors and use scaled gradient step */
87: f = fold;
88: PetscCall(VecCopy(lmP->Xold, tao->solution));
89: PetscCall(VecCopy(lmP->Gold, tao->gradient));
91: /* Failed to obtain acceptable iterate with BFGS step */
92: /* Attempt to use the scaled gradient direction */
94: PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
95: PetscCall(MatLMVMClearJ0(lmP->M));
96: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
97: PetscCall(MatSolve(lmP->M, tao->solution, tao->gradient));
99: /* On a reset, the direction cannot be not a number; it is a
100: scaled gradient step. No need to check for this condition. */
101: stepType = LMVM_STEP_GRAD;
102: PetscCall(VecScale(lmP->D, -1.0));
104: /* Perform the linesearch */
105: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status));
106: PetscCall(TaoAddLineSearchCounts(tao));
107: }
109: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
110: /* Failed to find an improving point */
111: f = fold;
112: PetscCall(VecCopy(lmP->Xold, tao->solution));
113: PetscCall(VecCopy(lmP->Gold, tao->gradient));
114: step = 0.0;
115: tao->reason = TAO_DIVERGED_LS_FAILURE;
116: } else {
117: /* LS found valid step, so tally up step type */
118: switch (stepType) {
119: case LMVM_STEP_BFGS:
120: ++lmP->bfgs;
121: break;
122: case LMVM_STEP_GRAD:
123: ++lmP->grad;
124: break;
125: default:
126: break;
127: }
128: /* Compute new gradient norm */
129: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
130: }
132: /* Check convergence */
133: tao->niter++;
134: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
135: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
136: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
137: }
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
142: {
143: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
144: PetscInt n, N;
145: PetscBool is_set, is_spd;
147: PetscFunctionBegin;
148: /* Existence of tao->solution checked in TaoSetUp() */
149: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
150: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
151: if (!lmP->D) PetscCall(VecDuplicate(tao->solution, &lmP->D));
152: if (!lmP->Xold) PetscCall(VecDuplicate(tao->solution, &lmP->Xold));
153: if (!lmP->Gold) PetscCall(VecDuplicate(tao->solution, &lmP->Gold));
155: /* Create matrix for the limited memory approximation */
156: PetscCall(VecGetLocalSize(tao->solution, &n));
157: PetscCall(VecGetSize(tao->solution, &N));
158: PetscCall(MatSetSizes(lmP->M, n, n, N, N));
159: PetscCall(MatLMVMAllocate(lmP->M, tao->solution, tao->gradient));
160: PetscCall(MatIsSPDKnown(lmP->M, &is_set, &is_spd));
161: PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite.");
163: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
164: if (lmP->H0) PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0));
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: /* ---------------------------------------------------------- */
169: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
170: {
171: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
173: PetscFunctionBegin;
174: if (tao->setupcalled) {
175: PetscCall(VecDestroy(&lmP->Xold));
176: PetscCall(VecDestroy(&lmP->Gold));
177: PetscCall(VecDestroy(&lmP->D));
178: }
179: PetscCall(MatDestroy(&lmP->M));
180: if (lmP->H0) PetscCall(PetscObjectDereference((PetscObject)lmP->H0));
181: PetscCall(PetscFree(tao->data));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /*------------------------------------------------------------*/
186: static PetscErrorCode TaoSetFromOptions_LMVM(Tao tao, PetscOptionItems *PetscOptionsObject)
187: {
188: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
190: PetscFunctionBegin;
191: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for unconstrained optimization");
192: PetscCall(PetscOptionsBool("-tao_lmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", lm->recycle, &lm->recycle, NULL));
193: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
194: PetscCall(MatSetFromOptions(lm->M));
195: PetscOptionsHeadEnd();
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: /*------------------------------------------------------------*/
200: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
201: {
202: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
203: PetscBool isascii;
204: PetscInt recycled_its;
206: PetscFunctionBegin;
207: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
208: if (isascii) {
209: PetscCall(PetscViewerASCIIPushTab(viewer));
210: PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lm->grad));
211: if (lm->recycle) {
212: PetscCall(PetscViewerASCIIPrintf(viewer, "Recycle: on\n"));
213: recycled_its = lm->bfgs + lm->grad;
214: PetscCall(PetscViewerASCIIPrintf(viewer, "Total recycled iterations: %" PetscInt_FMT "\n", recycled_its));
215: }
216: PetscCall(PetscViewerASCIIPrintf(viewer, "LMVM Matrix:\n"));
217: PetscCall(PetscViewerASCIIPushTab(viewer));
218: PetscCall(MatView(lm->M, viewer));
219: PetscCall(PetscViewerASCIIPopTab(viewer));
220: PetscCall(PetscViewerASCIIPopTab(viewer));
221: }
222: PetscFunctionReturn(PETSC_SUCCESS);
223: }
225: /* ---------------------------------------------------------- */
227: /*MC
228: TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
229: optimization solver for unconstrained minimization. It solves
230: the Newton step
231: Hkdk = - gk
233: using an approximation Bk in place of Hk, where Bk is composed using
234: the BFGS update formula. A More-Thuente line search is then used
235: to computed the steplength in the dk direction
237: Options Database Keys:
238: + -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls
239: - -tao_lmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation
241: Level: beginner
242: M*/
244: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
245: {
246: TAO_LMVM *lmP;
247: const char *morethuente_type = TAOLINESEARCHMT;
249: PetscFunctionBegin;
250: tao->ops->setup = TaoSetUp_LMVM;
251: tao->ops->solve = TaoSolve_LMVM;
252: tao->ops->view = TaoView_LMVM;
253: tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
254: tao->ops->destroy = TaoDestroy_LMVM;
256: PetscCall(PetscNew(&lmP));
257: lmP->D = NULL;
258: lmP->M = NULL;
259: lmP->Xold = NULL;
260: lmP->Gold = NULL;
261: lmP->H0 = NULL;
262: lmP->recycle = PETSC_FALSE;
264: tao->data = (void *)lmP;
265: /* Override default settings (unless already changed) */
266: PetscCall(TaoParametersInitialize(tao));
267: PetscObjectParameterSetDefault(tao, max_it, 2000);
268: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
270: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
271: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
272: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
273: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
274: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
276: PetscCall(KSPInitializePackage());
277: PetscCall(MatCreate(((PetscObject)tao)->comm, &lmP->M));
278: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1));
279: PetscCall(MatSetType(lmP->M, MATLMVMBFGS));
280: PetscCall(MatSetOptionsPrefix(lmP->M, "tao_lmvm_"));
281: PetscFunctionReturn(PETSC_SUCCESS);
282: }