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;
17: 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"));
19: /* Check convergence criteria */
20: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
21: 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 Inf 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, step));
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 (!lmP->recycle) {
33: lmP->bfgs = 0;
34: lmP->grad = 0;
35: PetscCall(MatLMVMReset(lmP->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: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
43: /* Compute direction */
44: if (lmP->H0) {
45: PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0));
46: stepType = LMVM_STEP_BFGS;
47: }
48: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
49: PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D));
50: PetscCall(MatLMVMGetUpdateCount(lmP->M, &nupdates));
51: if (nupdates > 0) stepType = LMVM_STEP_BFGS;
53: /* Check for success (descent direction) */
54: PetscCall(VecDot(lmP->D, tao->gradient, &gdx));
55: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
56: /* Step is not descent or direction produced not a number
57: We can assert bfgsUpdates > 1 in this case because
58: the first solve produces the scaled gradient direction,
59: which is guaranteed to be descent
61: Use steepest descent direction (scaled)
62: */
64: PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
65: PetscCall(MatLMVMClearJ0(lmP->M));
66: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
67: PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D));
69: /* On a reset, the direction cannot be not a number; it is a
70: scaled gradient step. No need to check for this condition. */
71: stepType = LMVM_STEP_GRAD;
72: }
73: PetscCall(VecScale(lmP->D, -1.0));
75: /* Perform the linesearch */
76: fold = f;
77: PetscCall(VecCopy(tao->solution, lmP->Xold));
78: PetscCall(VecCopy(tao->gradient, lmP->Gold));
80: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status));
81: PetscCall(TaoAddLineSearchCounts(tao));
83: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) {
84: /* Reset factors and use scaled gradient step */
85: f = fold;
86: PetscCall(VecCopy(lmP->Xold, tao->solution));
87: PetscCall(VecCopy(lmP->Gold, tao->gradient));
89: /* Failed to obtain acceptable iterate with BFGS step */
90: /* Attempt to use the scaled gradient direction */
92: PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
93: PetscCall(MatLMVMClearJ0(lmP->M));
94: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
95: PetscCall(MatSolve(lmP->M, tao->solution, tao->gradient));
97: /* On a reset, the direction cannot be not a number; it is a
98: scaled gradient step. No need to check for this condition. */
99: stepType = LMVM_STEP_GRAD;
100: PetscCall(VecScale(lmP->D, -1.0));
102: /* Perform the linesearch */
103: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status));
104: PetscCall(TaoAddLineSearchCounts(tao));
105: }
107: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
108: /* Failed to find an improving point */
109: f = fold;
110: PetscCall(VecCopy(lmP->Xold, tao->solution));
111: PetscCall(VecCopy(lmP->Gold, tao->gradient));
112: step = 0.0;
113: tao->reason = TAO_DIVERGED_LS_FAILURE;
114: } else {
115: /* LS found valid step, so tally up step type */
116: switch (stepType) {
117: case LMVM_STEP_BFGS:
118: ++lmP->bfgs;
119: break;
120: case LMVM_STEP_GRAD:
121: ++lmP->grad;
122: break;
123: default:
124: break;
125: }
126: /* Compute new gradient norm */
127: PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
128: }
130: /* Check convergence */
131: tao->niter++;
132: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
133: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
134: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
135: }
136: PetscFunctionReturn(PETSC_SUCCESS);
137: }
139: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
140: {
141: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
142: PetscInt n, N;
143: PetscBool is_set, is_spd;
145: PetscFunctionBegin;
146: /* Existence of tao->solution checked in TaoSetUp() */
147: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
148: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
149: if (!lmP->D) PetscCall(VecDuplicate(tao->solution, &lmP->D));
150: if (!lmP->Xold) PetscCall(VecDuplicate(tao->solution, &lmP->Xold));
151: if (!lmP->Gold) PetscCall(VecDuplicate(tao->solution, &lmP->Gold));
153: /* Create matrix for the limited memory approximation */
154: PetscCall(VecGetLocalSize(tao->solution, &n));
155: PetscCall(VecGetSize(tao->solution, &N));
156: PetscCall(MatSetSizes(lmP->M, n, n, N, N));
157: PetscCall(MatLMVMAllocate(lmP->M, tao->solution, tao->gradient));
158: PetscCall(MatIsSPDKnown(lmP->M, &is_set, &is_spd));
159: PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite.");
161: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
162: if (lmP->H0) PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
166: /* ---------------------------------------------------------- */
167: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
168: {
169: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
171: PetscFunctionBegin;
172: if (tao->setupcalled) {
173: PetscCall(VecDestroy(&lmP->Xold));
174: PetscCall(VecDestroy(&lmP->Gold));
175: PetscCall(VecDestroy(&lmP->D));
176: }
177: PetscCall(MatDestroy(&lmP->M));
178: if (lmP->H0) PetscCall(PetscObjectDereference((PetscObject)lmP->H0));
179: PetscCall(PetscFree(tao->data));
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*------------------------------------------------------------*/
184: static PetscErrorCode TaoSetFromOptions_LMVM(Tao tao, PetscOptionItems *PetscOptionsObject)
185: {
186: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
188: PetscFunctionBegin;
189: PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for unconstrained optimization");
190: PetscCall(PetscOptionsBool("-tao_lmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", lm->recycle, &lm->recycle, NULL));
191: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
192: PetscCall(MatSetFromOptions(lm->M));
193: PetscOptionsHeadEnd();
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*------------------------------------------------------------*/
198: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
199: {
200: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
201: PetscBool isascii;
202: PetscInt recycled_its;
204: PetscFunctionBegin;
205: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
206: if (isascii) {
207: PetscCall(PetscViewerASCIIPrintf(viewer, " Gradient steps: %" PetscInt_FMT "\n", lm->grad));
208: if (lm->recycle) {
209: PetscCall(PetscViewerASCIIPrintf(viewer, " Recycle: on\n"));
210: recycled_its = lm->bfgs + lm->grad;
211: PetscCall(PetscViewerASCIIPrintf(viewer, " Total recycled iterations: %" PetscInt_FMT "\n", recycled_its));
212: }
213: }
214: PetscFunctionReturn(PETSC_SUCCESS);
215: }
217: /* ---------------------------------------------------------- */
219: /*MC
220: TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
221: optimization solver for unconstrained minimization. It solves
222: the Newton step
223: Hkdk = - gk
225: using an approximation Bk in place of Hk, where Bk is composed using
226: the BFGS update formula. A More-Thuente line search is then used
227: to computed the steplength in the dk direction
229: Options Database Keys:
230: + -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls
231: - -tao_lmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation
233: Level: beginner
234: M*/
236: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
237: {
238: TAO_LMVM *lmP;
239: const char *morethuente_type = TAOLINESEARCHMT;
241: PetscFunctionBegin;
242: tao->ops->setup = TaoSetUp_LMVM;
243: tao->ops->solve = TaoSolve_LMVM;
244: tao->ops->view = TaoView_LMVM;
245: tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
246: tao->ops->destroy = TaoDestroy_LMVM;
248: PetscCall(PetscNew(&lmP));
249: lmP->D = NULL;
250: lmP->M = NULL;
251: lmP->Xold = NULL;
252: lmP->Gold = NULL;
253: lmP->H0 = NULL;
254: lmP->recycle = PETSC_FALSE;
256: tao->data = (void *)lmP;
257: /* Override default settings (unless already changed) */
258: if (!tao->max_it_changed) tao->max_it = 2000;
259: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
261: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
262: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
263: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
264: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
265: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
267: PetscCall(KSPInitializePackage());
268: PetscCall(MatCreate(((PetscObject)tao)->comm, &lmP->M));
269: PetscCall(PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1));
270: PetscCall(MatSetType(lmP->M, MATLMVMBFGS));
271: PetscCall(MatSetOptionsPrefix(lmP->M, "tao_lmvm_"));
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }