Actual source code: lmvm.c
petsc-3.6.4 2016-04-12
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/matrix/lmvmmat.h>
3: #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>
5: #define LMVM_BFGS 0
6: #define LMVM_SCALED_GRADIENT 1
7: #define LMVM_GRADIENT 2
11: static PetscErrorCode TaoSolve_LMVM(Tao tao)
12: {
13: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
14: PetscReal f, fold, gdx, gnorm;
15: PetscReal step = 1.0;
16: PetscReal delta;
17: PetscErrorCode ierr;
18: PetscInt stepType;
19: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
20: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
24: if (tao->XL || tao->XU || tao->ops->computebounds) {
25: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");
26: }
28: /* Check convergence criteria */
29: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
30: VecNorm(tao->gradient,NORM_2,&gnorm);
31: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
33: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);
34: if (reason != TAO_CONTINUE_ITERATING) return(0);
36: /* Set initial scaling for the function */
37: if (f != 0.0) {
38: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
39: } else {
40: delta = 2.0 / (gnorm*gnorm);
41: }
42: MatLMVMSetDelta(lmP->M,delta);
44: /* Set counter for gradient/reset steps */
45: lmP->bfgs = 0;
46: lmP->sgrad = 0;
47: lmP->grad = 0;
49: /* Have not converged; continue with Newton method */
50: while (reason == TAO_CONTINUE_ITERATING) {
51: /* Compute direction */
52: MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
53: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
54: ++lmP->bfgs;
56: /* Check for success (descent direction) */
57: VecDot(lmP->D, tao->gradient, &gdx);
58: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
59: /* Step is not descent or direction produced not a number
60: We can assert bfgsUpdates > 1 in this case because
61: the first solve produces the scaled gradient direction,
62: which is guaranteed to be descent
64: Use steepest descent direction (scaled)
65: */
67: ++lmP->grad;
69: if (f != 0.0) {
70: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
71: } else {
72: delta = 2.0 / (gnorm*gnorm);
73: }
74: MatLMVMSetDelta(lmP->M, delta);
75: MatLMVMReset(lmP->M);
76: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
77: MatLMVMSolve(lmP->M,tao->gradient, lmP->D);
79: /* On a reset, the direction cannot be not a number; it is a
80: scaled gradient step. No need to check for this condition. */
82: lmP->bfgs = 1;
83: ++lmP->sgrad;
84: stepType = LMVM_SCALED_GRADIENT;
85: } else {
86: if (1 == lmP->bfgs) {
87: /* The first BFGS direction is always the scaled gradient */
88: ++lmP->sgrad;
89: stepType = LMVM_SCALED_GRADIENT;
90: } else {
91: ++lmP->bfgs;
92: stepType = LMVM_BFGS;
93: }
94: }
95: VecScale(lmP->D, -1.0);
97: /* Perform the linesearch */
98: fold = f;
99: VecCopy(tao->solution, lmP->Xold);
100: VecCopy(tao->gradient, lmP->Gold);
102: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);
103: TaoAddLineSearchCounts(tao);
105: while (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_GRADIENT)) {
106: /* Linesearch failed */
107: /* Reset factors and use scaled gradient step */
108: f = fold;
109: VecCopy(lmP->Xold, tao->solution);
110: VecCopy(lmP->Gold, tao->gradient);
112: switch(stepType) {
113: case LMVM_BFGS:
114: /* Failed to obtain acceptable iterate with BFGS step */
115: /* Attempt to use the scaled gradient direction */
117: if (f != 0.0) {
118: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
119: } else {
120: delta = 2.0 / (gnorm*gnorm);
121: }
122: MatLMVMSetDelta(lmP->M, delta);
123: MatLMVMReset(lmP->M);
124: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
125: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
127: /* On a reset, the direction cannot be not a number; it is a
128: scaled gradient step. No need to check for this condition. */
130: lmP->bfgs = 1;
131: ++lmP->sgrad;
132: stepType = LMVM_SCALED_GRADIENT;
133: break;
135: case LMVM_SCALED_GRADIENT:
136: /* The scaled gradient step did not produce a new iterate;
137: attempt to use the gradient direction.
138: Need to make sure we are not using a different diagonal scaling */
139: MatLMVMSetDelta(lmP->M, 1.0);
140: MatLMVMReset(lmP->M);
141: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
142: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
144: lmP->bfgs = 1;
145: ++lmP->grad;
146: stepType = LMVM_GRADIENT;
147: break;
148: }
149: VecScale(lmP->D, -1.0);
151: /* Perform the linesearch */
152: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);
153: TaoAddLineSearchCounts(tao);
154: }
156: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
157: /* Failed to find an improving point */
158: f = fold;
159: VecCopy(lmP->Xold, tao->solution);
160: VecCopy(lmP->Gold, tao->gradient);
161: step = 0.0;
162: reason = TAO_DIVERGED_LS_FAILURE;
163: tao->reason = TAO_DIVERGED_LS_FAILURE;
164: }
165: /* Check for termination */
166: VecNorm(tao->gradient, NORM_2, &gnorm);
167: tao->niter++;
168: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step,&reason);
169: }
170: return(0);
171: }
175: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
176: {
177: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
178: PetscInt n,N;
182: /* Existence of tao->solution checked in TaoSetUp() */
183: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
184: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
185: if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D); }
186: if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold); }
187: if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold); }
189: /* Create matrix for the limited memory approximation */
190: VecGetLocalSize(tao->solution,&n);
191: VecGetSize(tao->solution,&N);
192: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M);
193: MatLMVMAllocateVectors(lmP->M,tao->solution);
194: return(0);
195: }
197: /* ---------------------------------------------------------- */
200: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
201: {
202: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
206: if (tao->setupcalled) {
207: VecDestroy(&lmP->Xold);
208: VecDestroy(&lmP->Gold);
209: VecDestroy(&lmP->D);
210: MatDestroy(&lmP->M);
211: }
212: PetscFree(tao->data);
213: return(0);
214: }
216: /*------------------------------------------------------------*/
219: static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptions *PetscOptionsObject,Tao tao)
220: {
224: PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");
225: TaoLineSearchSetFromOptions(tao->linesearch);
226: PetscOptionsTail();
227: return(0);
228: }
230: /*------------------------------------------------------------*/
233: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
234: {
235: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
236: PetscBool isascii;
240: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
241: if (isascii) {
242: PetscViewerASCIIPushTab(viewer);
243: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
244: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
245: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
246: PetscViewerASCIIPopTab(viewer);
247: }
248: return(0);
249: }
251: /* ---------------------------------------------------------- */
253: /*MC
254: TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
255: optimization solver for unconstrained minimization. It solves
256: the Newton step
257: Hkdk = - gk
259: using an approximation Bk in place of Hk, where Bk is composed using
260: the BFGS update formula. A More-Thuente line search is then used
261: to computed the steplength in the dk direction
262: Options Database Keys:
263: + -tao_lmm_vectors - number of vectors to use for approximation
264: . -tao_lmm_scale_type - "none","scalar","broyden"
265: . -tao_lmm_limit_type - "none","average","relative","absolute"
266: . -tao_lmm_rescale_type - "none","scalar","gl"
267: . -tao_lmm_limit_mu - mu limiting factor
268: . -tao_lmm_limit_nu - nu limiting factor
269: . -tao_lmm_delta_min - minimum delta value
270: . -tao_lmm_delta_max - maximum delta value
271: . -tao_lmm_broyden_phi - phi factor for Broyden scaling
272: . -tao_lmm_scalar_alpha - alpha factor for scalar scaling
273: . -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal
274: . -tao_lmm_rescale_beta - beta factor for rescaling diagonal
275: . -tao_lmm_scalar_history - amount of history for scalar scaling
276: . -tao_lmm_rescale_history - amount of history for rescaling diagonal
277: - -tao_lmm_eps - rejection tolerance
279: Level: beginner
280: M*/
284: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
285: {
286: TAO_LMVM *lmP;
287: const char *morethuente_type = TAOLINESEARCHMT;
291: tao->ops->setup = TaoSetUp_LMVM;
292: tao->ops->solve = TaoSolve_LMVM;
293: tao->ops->view = TaoView_LMVM;
294: tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
295: tao->ops->destroy = TaoDestroy_LMVM;
297: PetscNewLog(tao,&lmP);
298: lmP->D = 0;
299: lmP->M = 0;
300: lmP->Xold = 0;
301: lmP->Gold = 0;
303: tao->data = (void*)lmP;
304: /* Override default settings (unless already changed) */
305: if (!tao->max_it_changed) tao->max_it = 2000;
306: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
307: if (!tao->fatol_changed) tao->fatol = 1.0e-4;
308: if (!tao->frtol_changed) tao->frtol = 1.0e-4;
310: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
311: TaoLineSearchSetType(tao->linesearch,morethuente_type);
312: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
313: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
314: return(0);
315: }