Actual source code: lmvm.c
petsc-3.7.5 2017-01-01
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: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
32: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
34: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);
35: if (reason != TAO_CONTINUE_ITERATING) return(0);
37: /* Set initial scaling for the function */
38: if (f != 0.0) {
39: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
40: } else {
41: delta = 2.0 / (gnorm*gnorm);
42: }
43: MatLMVMSetDelta(lmP->M,delta);
45: /* Set counter for gradient/reset steps */
46: lmP->bfgs = 0;
47: lmP->sgrad = 0;
48: lmP->grad = 0;
50: /* Have not converged; continue with Newton method */
51: while (reason == TAO_CONTINUE_ITERATING) {
52: /* Compute direction */
53: MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
54: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
55: ++lmP->bfgs;
57: /* Check for success (descent direction) */
58: VecDot(lmP->D, tao->gradient, &gdx);
59: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
60: /* Step is not descent or direction produced not a number
61: We can assert bfgsUpdates > 1 in this case because
62: the first solve produces the scaled gradient direction,
63: which is guaranteed to be descent
65: Use steepest descent direction (scaled)
66: */
68: ++lmP->grad;
70: if (f != 0.0) {
71: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
72: } else {
73: delta = 2.0 / (gnorm*gnorm);
74: }
75: MatLMVMSetDelta(lmP->M, delta);
76: MatLMVMReset(lmP->M);
77: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
78: MatLMVMSolve(lmP->M,tao->gradient, lmP->D);
80: /* On a reset, the direction cannot be not a number; it is a
81: scaled gradient step. No need to check for this condition. */
83: lmP->bfgs = 1;
84: ++lmP->sgrad;
85: stepType = LMVM_SCALED_GRADIENT;
86: } else {
87: if (1 == lmP->bfgs) {
88: /* The first BFGS direction is always the scaled gradient */
89: ++lmP->sgrad;
90: stepType = LMVM_SCALED_GRADIENT;
91: } else {
92: ++lmP->bfgs;
93: stepType = LMVM_BFGS;
94: }
95: }
96: VecScale(lmP->D, -1.0);
98: /* Perform the linesearch */
99: fold = f;
100: VecCopy(tao->solution, lmP->Xold);
101: VecCopy(tao->gradient, lmP->Gold);
103: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);
104: TaoAddLineSearchCounts(tao);
106: while (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_GRADIENT)) {
107: /* Linesearch failed */
108: /* Reset factors and use scaled gradient step */
109: f = fold;
110: VecCopy(lmP->Xold, tao->solution);
111: VecCopy(lmP->Gold, tao->gradient);
113: switch(stepType) {
114: case LMVM_BFGS:
115: /* Failed to obtain acceptable iterate with BFGS step */
116: /* Attempt to use the scaled gradient direction */
118: if (f != 0.0) {
119: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
120: } else {
121: delta = 2.0 / (gnorm*gnorm);
122: }
123: MatLMVMSetDelta(lmP->M, delta);
124: MatLMVMReset(lmP->M);
125: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
126: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
128: /* On a reset, the direction cannot be not a number; it is a
129: scaled gradient step. No need to check for this condition. */
131: lmP->bfgs = 1;
132: ++lmP->sgrad;
133: stepType = LMVM_SCALED_GRADIENT;
134: break;
136: case LMVM_SCALED_GRADIENT:
137: /* The scaled gradient step did not produce a new iterate;
138: attempt to use the gradient direction.
139: Need to make sure we are not using a different diagonal scaling */
140: MatLMVMSetDelta(lmP->M, 1.0);
141: MatLMVMReset(lmP->M);
142: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
143: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
145: lmP->bfgs = 1;
146: ++lmP->grad;
147: stepType = LMVM_GRADIENT;
148: break;
149: }
150: VecScale(lmP->D, -1.0);
152: /* Perform the linesearch */
153: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);
154: TaoAddLineSearchCounts(tao);
155: }
157: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
158: /* Failed to find an improving point */
159: f = fold;
160: VecCopy(lmP->Xold, tao->solution);
161: VecCopy(lmP->Gold, tao->gradient);
162: step = 0.0;
163: reason = TAO_DIVERGED_LS_FAILURE;
164: tao->reason = TAO_DIVERGED_LS_FAILURE;
165: }
167: /* Check for termination */
168: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
170: tao->niter++;
171: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step,&reason);
172: }
173: return(0);
174: }
178: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
179: {
180: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
181: PetscInt n,N;
183: KSP H0ksp;
186: /* Existence of tao->solution checked in TaoSetUp() */
187: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
188: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
189: if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D); }
190: if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold); }
191: if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold); }
193: /* Create matrix for the limited memory approximation */
194: VecGetLocalSize(tao->solution,&n);
195: VecGetSize(tao->solution,&N);
196: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M);
197: MatLMVMAllocateVectors(lmP->M,tao->solution);
199: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
200: if (lmP->H0) {
201: const char *prefix;
202: PC H0pc;
204: MatLMVMSetH0(lmP->M, lmP->H0);
205: MatLMVMGetH0KSP(lmP->M, &H0ksp);
207: TaoGetOptionsPrefix(tao, &prefix);
208: KSPSetOptionsPrefix(H0ksp, prefix);
209: PetscObjectAppendOptionsPrefix((PetscObject)H0ksp, "tao_h0_");
210: KSPGetPC(H0ksp, &H0pc);
211: PetscObjectAppendOptionsPrefix((PetscObject)H0pc, "tao_h0_");
213: KSPSetFromOptions(H0ksp);
214: KSPSetUp(H0ksp);
215: }
217: return(0);
218: }
220: /* ---------------------------------------------------------- */
223: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
224: {
225: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
229: if (tao->setupcalled) {
230: VecDestroy(&lmP->Xold);
231: VecDestroy(&lmP->Gold);
232: VecDestroy(&lmP->D);
233: MatDestroy(&lmP->M);
234: }
236: if (lmP->H0) {
237: PetscObjectDereference((PetscObject)lmP->H0);
238: }
240: PetscFree(tao->data);
242: return(0);
243: }
245: /*------------------------------------------------------------*/
248: static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject,Tao tao)
249: {
253: PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");
254: TaoLineSearchSetFromOptions(tao->linesearch);
255: PetscOptionsTail();
256: return(0);
257: }
259: /*------------------------------------------------------------*/
262: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
263: {
264: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
265: PetscBool isascii;
269: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
270: if (isascii) {
271: PetscViewerASCIIPushTab(viewer);
272: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
273: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
274: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
275: PetscViewerASCIIPopTab(viewer);
276: }
277: return(0);
278: }
280: /* ---------------------------------------------------------- */
282: /*MC
283: TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
284: optimization solver for unconstrained minimization. It solves
285: the Newton step
286: Hkdk = - gk
288: using an approximation Bk in place of Hk, where Bk is composed using
289: the BFGS update formula. A More-Thuente line search is then used
290: to computed the steplength in the dk direction
291: Options Database Keys:
292: + -tao_lmm_vectors - number of vectors to use for approximation
293: . -tao_lmm_scale_type - "none","scalar","broyden"
294: . -tao_lmm_limit_type - "none","average","relative","absolute"
295: . -tao_lmm_rescale_type - "none","scalar","gl"
296: . -tao_lmm_limit_mu - mu limiting factor
297: . -tao_lmm_limit_nu - nu limiting factor
298: . -tao_lmm_delta_min - minimum delta value
299: . -tao_lmm_delta_max - maximum delta value
300: . -tao_lmm_broyden_phi - phi factor for Broyden scaling
301: . -tao_lmm_scalar_alpha - alpha factor for scalar scaling
302: . -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal
303: . -tao_lmm_rescale_beta - beta factor for rescaling diagonal
304: . -tao_lmm_scalar_history - amount of history for scalar scaling
305: . -tao_lmm_rescale_history - amount of history for rescaling diagonal
306: - -tao_lmm_eps - rejection tolerance
308: Level: beginner
309: M*/
313: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
314: {
315: TAO_LMVM *lmP;
316: const char *morethuente_type = TAOLINESEARCHMT;
320: tao->ops->setup = TaoSetUp_LMVM;
321: tao->ops->solve = TaoSolve_LMVM;
322: tao->ops->view = TaoView_LMVM;
323: tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
324: tao->ops->destroy = TaoDestroy_LMVM;
326: PetscNewLog(tao,&lmP);
327: lmP->D = 0;
328: lmP->M = 0;
329: lmP->Xold = 0;
330: lmP->Gold = 0;
331: lmP->H0 = NULL;
333: tao->data = (void*)lmP;
334: /* Override default settings (unless already changed) */
335: if (!tao->max_it_changed) tao->max_it = 2000;
336: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
338: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
339: TaoLineSearchSetType(tao->linesearch,morethuente_type);
340: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
341: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
342: return(0);
343: }