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