Actual source code: lmvm.c
petsc-3.14.6 2021-03-30
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: PetscErrorCode ierr;
13: PetscInt stepType = LMVM_STEP_GRAD, nupdates;
14: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
18: if (tao->XL || tao->XU || tao->ops->computebounds) {
19: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");
20: }
22: /* Check convergence criteria */
23: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
24: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
26: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");
28: tao->reason = TAO_CONTINUE_ITERATING;
29: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
30: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
31: (*tao->ops->convergencetest)(tao,tao->cnvP);
32: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
34: /* Set counter for gradient/reset steps */
35: if (!lmP->recycle) {
36: lmP->bfgs = 0;
37: lmP->grad = 0;
38: MatLMVMReset(lmP->M, PETSC_FALSE);
39: }
41: /* Have not converged; continue with Newton method */
42: while (tao->reason == TAO_CONTINUE_ITERATING) {
43: /* Call general purpose update function */
44: if (tao->ops->update) {
45: (*tao->ops->update)(tao, tao->niter, tao->user_update);
46: }
48: /* Compute direction */
49: if (lmP->H0) {
50: MatLMVMSetJ0(lmP->M, lmP->H0);
51: stepType = LMVM_STEP_BFGS;
52: }
53: MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
54: MatSolve(lmP->M, tao->gradient, lmP->D);
55: MatLMVMGetUpdateCount(lmP->M, &nupdates);
56: if (nupdates > 0) stepType = LMVM_STEP_BFGS;
58: /* Check for success (descent direction) */
59: VecDot(lmP->D, tao->gradient, &gdx);
60: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
61: /* Step is not descent or direction produced not a number
62: We can assert bfgsUpdates > 1 in this case because
63: the first solve produces the scaled gradient direction,
64: which is guaranteed to be descent
66: Use steepest descent direction (scaled)
67: */
69: MatLMVMReset(lmP->M, PETSC_FALSE);
70: MatLMVMClearJ0(lmP->M);
71: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
72: MatSolve(lmP->M,tao->gradient, lmP->D);
74: /* On a reset, the direction cannot be not a number; it is a
75: scaled gradient step. No need to check for this condition. */
76: stepType = LMVM_STEP_GRAD;
77: }
78: VecScale(lmP->D, -1.0);
80: /* Perform the linesearch */
81: fold = f;
82: VecCopy(tao->solution, lmP->Xold);
83: VecCopy(tao->gradient, lmP->Gold);
85: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);
86: TaoAddLineSearchCounts(tao);
88: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) {
89: /* Reset factors and use scaled gradient step */
90: f = fold;
91: VecCopy(lmP->Xold, tao->solution);
92: VecCopy(lmP->Gold, tao->gradient);
94: /* Failed to obtain acceptable iterate with BFGS step */
95: /* Attempt to use the scaled gradient direction */
97: MatLMVMReset(lmP->M, PETSC_FALSE);
98: MatLMVMClearJ0(lmP->M);
99: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
100: MatSolve(lmP->M, tao->solution, tao->gradient);
102: /* On a reset, the direction cannot be not a number; it is a
103: scaled gradient step. No need to check for this condition. */
104: stepType = LMVM_STEP_GRAD;
105: VecScale(lmP->D, -1.0);
107: /* Perform the linesearch */
108: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);
109: TaoAddLineSearchCounts(tao);
110: }
112: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
113: /* Failed to find an improving point */
114: f = fold;
115: VecCopy(lmP->Xold, tao->solution);
116: VecCopy(lmP->Gold, tao->gradient);
117: step = 0.0;
118: tao->reason = TAO_DIVERGED_LS_FAILURE;
119: } else {
120: /* LS found valid step, so tally up step type */
121: switch (stepType) {
122: case LMVM_STEP_BFGS:
123: ++lmP->bfgs;
124: break;
125: case LMVM_STEP_GRAD:
126: ++lmP->grad;
127: break;
128: default:
129: break;
130: }
131: /* Compute new gradient norm */
132: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
133: }
135: /* Check convergence */
136: tao->niter++;
137: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
138: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
139: (*tao->ops->convergencetest)(tao,tao->cnvP);
140: }
141: return(0);
142: }
144: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
145: {
146: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
147: PetscInt n,N;
149: PetscBool is_spd;
152: /* Existence of tao->solution checked in TaoSetUp() */
153: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
154: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
155: if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D); }
156: if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold); }
157: if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold); }
159: /* Create matrix for the limited memory approximation */
160: VecGetLocalSize(tao->solution,&n);
161: VecGetSize(tao->solution,&N);
162: MatSetSizes(lmP->M, n, n, N, N);
163: MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);
164: MatGetOption(lmP->M, MAT_SPD, &is_spd);
165: if (!is_spd) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite.");
167: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
168: if (lmP->H0) {
169: MatLMVMSetJ0(lmP->M, lmP->H0);
170: }
172: return(0);
173: }
175: /* ---------------------------------------------------------- */
176: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
177: {
178: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
182: if (tao->setupcalled) {
183: VecDestroy(&lmP->Xold);
184: VecDestroy(&lmP->Gold);
185: VecDestroy(&lmP->D);
186: }
187: MatDestroy(&lmP->M);
188: if (lmP->H0) {
189: PetscObjectDereference((PetscObject)lmP->H0);
190: }
191: PetscFree(tao->data);
193: return(0);
194: }
196: /*------------------------------------------------------------*/
197: static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject,Tao tao)
198: {
199: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
203: PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");
204: PetscOptionsBool("-tao_lmvm_recycle","enable recycling of the BFGS matrix between subsequent TaoSolve() calls","",lm->recycle,&lm->recycle,NULL);
205: TaoLineSearchSetFromOptions(tao->linesearch);
206: MatSetFromOptions(lm->M);
207: PetscOptionsTail();
208: return(0);
209: }
211: /*------------------------------------------------------------*/
212: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
213: {
214: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
215: PetscBool isascii;
216: PetscInt recycled_its;
220: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
221: if (isascii) {
222: PetscViewerASCIIPrintf(viewer, " Gradient steps: %D\n", lm->grad);
223: if (lm->recycle) {
224: PetscViewerASCIIPrintf(viewer, " Recycle: on\n");
225: recycled_its = lm->bfgs + lm->grad;
226: PetscViewerASCIIPrintf(viewer, " Total recycled iterations: %D\n", recycled_its);
227: }
228: }
229: return(0);
230: }
232: /* ---------------------------------------------------------- */
234: /*MC
235: TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
236: optimization solver for unconstrained minimization. It solves
237: the Newton step
238: Hkdk = - gk
240: using an approximation Bk in place of Hk, where Bk is composed using
241: the BFGS update formula. A More-Thuente line search is then used
242: to computed the steplength in the dk direction
244: Options Database Keys:
245: + -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls
246: - -tao_lmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation
248: Level: beginner
249: M*/
251: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
252: {
253: TAO_LMVM *lmP;
254: const char *morethuente_type = TAOLINESEARCHMT;
258: tao->ops->setup = TaoSetUp_LMVM;
259: tao->ops->solve = TaoSolve_LMVM;
260: tao->ops->view = TaoView_LMVM;
261: tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
262: tao->ops->destroy = TaoDestroy_LMVM;
264: PetscNewLog(tao,&lmP);
265: lmP->D = NULL;
266: lmP->M = NULL;
267: lmP->Xold = NULL;
268: lmP->Gold = NULL;
269: lmP->H0 = NULL;
270: lmP->recycle = PETSC_FALSE;
272: tao->data = (void*)lmP;
273: /* Override default settings (unless already changed) */
274: if (!tao->max_it_changed) tao->max_it = 2000;
275: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
277: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
278: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
279: TaoLineSearchSetType(tao->linesearch,morethuente_type);
280: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
281: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
283: KSPInitializePackage();
284: MatCreate(((PetscObject)tao)->comm, &lmP->M);
285: PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1);
286: MatSetType(lmP->M, MATLMVMBFGS);
287: MatSetOptionsPrefix(lmP->M, "tao_lmvm_");
288: return(0);
289: }