Actual source code: lmvm.c
petsc-3.8.4 2018-03-24
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
9: static PetscErrorCode TaoSolve_LMVM(Tao tao)
10: {
11: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
12: PetscReal f, fold, gdx, gnorm;
13: PetscReal step = 1.0;
14: PetscReal delta;
15: PetscErrorCode ierr;
16: PetscInt stepType;
17: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
18: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
22: if (tao->XL || tao->XU || tao->ops->computebounds) {
23: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");
24: }
26: /* Check convergence criteria */
27: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
28: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
30: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
32: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);
33: if (reason != TAO_CONTINUE_ITERATING) return(0);
35: /* Set initial scaling for the function */
36: if (f != 0.0) {
37: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
38: } else {
39: delta = 2.0 / (gnorm*gnorm);
40: }
41: MatLMVMSetDelta(lmP->M,delta);
43: /* Set counter for gradient/reset steps */
44: lmP->bfgs = 0;
45: lmP->sgrad = 0;
46: lmP->grad = 0;
48: /* Have not converged; continue with Newton method */
49: while (reason == TAO_CONTINUE_ITERATING) {
50: /* Compute direction */
51: MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
52: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
53: ++lmP->bfgs;
55: /* Check for success (descent direction) */
56: VecDot(lmP->D, tao->gradient, &gdx);
57: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
58: /* Step is not descent or direction produced not a number
59: We can assert bfgsUpdates > 1 in this case because
60: the first solve produces the scaled gradient direction,
61: which is guaranteed to be descent
63: Use steepest descent direction (scaled)
64: */
66: ++lmP->grad;
68: if (f != 0.0) {
69: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
70: } else {
71: delta = 2.0 / (gnorm*gnorm);
72: }
73: MatLMVMSetDelta(lmP->M, delta);
74: MatLMVMReset(lmP->M);
75: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
76: MatLMVMSolve(lmP->M,tao->gradient, lmP->D);
78: /* On a reset, the direction cannot be not a number; it is a
79: scaled gradient step. No need to check for this condition. */
81: lmP->bfgs = 1;
82: ++lmP->sgrad;
83: stepType = LMVM_SCALED_GRADIENT;
84: } else {
85: if (1 == lmP->bfgs) {
86: /* The first BFGS direction is always the scaled gradient */
87: ++lmP->sgrad;
88: stepType = LMVM_SCALED_GRADIENT;
89: } else {
90: ++lmP->bfgs;
91: stepType = LMVM_BFGS;
92: }
93: }
94: VecScale(lmP->D, -1.0);
96: /* Perform the linesearch */
97: fold = f;
98: VecCopy(tao->solution, lmP->Xold);
99: VecCopy(tao->gradient, lmP->Gold);
101: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);
102: TaoAddLineSearchCounts(tao);
104: while (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_GRADIENT)) {
105: /* Linesearch failed */
106: /* Reset factors and use scaled gradient step */
107: f = fold;
108: VecCopy(lmP->Xold, tao->solution);
109: VecCopy(lmP->Gold, tao->gradient);
111: switch(stepType) {
112: case LMVM_BFGS:
113: /* Failed to obtain acceptable iterate with BFGS step */
114: /* Attempt to use the scaled gradient direction */
116: if (f != 0.0) {
117: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
118: } else {
119: delta = 2.0 / (gnorm*gnorm);
120: }
121: MatLMVMSetDelta(lmP->M, delta);
122: MatLMVMReset(lmP->M);
123: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
124: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
126: /* On a reset, the direction cannot be not a number; it is a
127: scaled gradient step. No need to check for this condition. */
129: lmP->bfgs = 1;
130: ++lmP->sgrad;
131: stepType = LMVM_SCALED_GRADIENT;
132: break;
134: case LMVM_SCALED_GRADIENT:
135: /* The scaled gradient step did not produce a new iterate;
136: attempt to use the gradient direction.
137: Need to make sure we are not using a different diagonal scaling */
138: MatLMVMSetDelta(lmP->M, 1.0);
139: MatLMVMReset(lmP->M);
140: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
141: MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
143: lmP->bfgs = 1;
144: ++lmP->grad;
145: stepType = LMVM_GRADIENT;
146: break;
147: }
148: VecScale(lmP->D, -1.0);
150: /* Perform the linesearch */
151: TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);
152: TaoAddLineSearchCounts(tao);
153: }
155: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
156: /* Failed to find an improving point */
157: f = fold;
158: VecCopy(lmP->Xold, tao->solution);
159: VecCopy(lmP->Gold, tao->gradient);
160: step = 0.0;
161: reason = TAO_DIVERGED_LS_FAILURE;
162: tao->reason = TAO_DIVERGED_LS_FAILURE;
163: }
165: /* Check for termination */
166: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
168: tao->niter++;
169: TaoMonitor(tao,tao->niter,f,gnorm,0.0,step,&reason);
170: }
171: return(0);
172: }
174: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
175: {
176: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
177: PetscInt n,N;
179: KSP H0ksp;
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);
195: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
196: if (lmP->H0) {
197: const char *prefix;
198: PC H0pc;
200: MatLMVMSetH0(lmP->M, lmP->H0);
201: MatLMVMGetH0KSP(lmP->M, &H0ksp);
203: TaoGetOptionsPrefix(tao, &prefix);
204: KSPSetOptionsPrefix(H0ksp, prefix);
205: PetscObjectAppendOptionsPrefix((PetscObject)H0ksp, "tao_h0_");
206: KSPGetPC(H0ksp, &H0pc);
207: PetscObjectAppendOptionsPrefix((PetscObject)H0pc, "tao_h0_");
209: KSPSetFromOptions(H0ksp);
210: KSPSetUp(H0ksp);
211: }
213: return(0);
214: }
216: /* ---------------------------------------------------------- */
217: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
218: {
219: TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
223: if (tao->setupcalled) {
224: VecDestroy(&lmP->Xold);
225: VecDestroy(&lmP->Gold);
226: VecDestroy(&lmP->D);
227: MatDestroy(&lmP->M);
228: }
230: if (lmP->H0) {
231: PetscObjectDereference((PetscObject)lmP->H0);
232: }
234: PetscFree(tao->data);
236: return(0);
237: }
239: /*------------------------------------------------------------*/
240: static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject,Tao tao)
241: {
245: PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");
246: TaoLineSearchSetFromOptions(tao->linesearch);
247: PetscOptionsTail();
248: return(0);
249: }
251: /*------------------------------------------------------------*/
252: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
253: {
254: TAO_LMVM *lm = (TAO_LMVM *)tao->data;
255: PetscBool isascii;
259: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
260: if (isascii) {
261: PetscViewerASCIIPushTab(viewer);
262: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
263: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
264: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
265: PetscViewerASCIIPopTab(viewer);
266: }
267: return(0);
268: }
270: /* ---------------------------------------------------------- */
272: /*MC
273: TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
274: optimization solver for unconstrained minimization. It solves
275: the Newton step
276: Hkdk = - gk
278: using an approximation Bk in place of Hk, where Bk is composed using
279: the BFGS update formula. A More-Thuente line search is then used
280: to computed the steplength in the dk direction
281: Options Database Keys:
282: + -tao_lmm_vectors - number of vectors to use for approximation
283: . -tao_lmm_scale_type - "none","scalar","broyden"
284: . -tao_lmm_limit_type - "none","average","relative","absolute"
285: . -tao_lmm_rescale_type - "none","scalar","gl"
286: . -tao_lmm_limit_mu - mu limiting factor
287: . -tao_lmm_limit_nu - nu limiting factor
288: . -tao_lmm_delta_min - minimum delta value
289: . -tao_lmm_delta_max - maximum delta value
290: . -tao_lmm_broyden_phi - phi factor for Broyden scaling
291: . -tao_lmm_scalar_alpha - alpha factor for scalar scaling
292: . -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal
293: . -tao_lmm_rescale_beta - beta factor for rescaling diagonal
294: . -tao_lmm_scalar_history - amount of history for scalar scaling
295: . -tao_lmm_rescale_history - amount of history for rescaling diagonal
296: - -tao_lmm_eps - rejection tolerance
298: Level: beginner
299: M*/
301: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
302: {
303: TAO_LMVM *lmP;
304: const char *morethuente_type = TAOLINESEARCHMT;
308: tao->ops->setup = TaoSetUp_LMVM;
309: tao->ops->solve = TaoSolve_LMVM;
310: tao->ops->view = TaoView_LMVM;
311: tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
312: tao->ops->destroy = TaoDestroy_LMVM;
314: PetscNewLog(tao,&lmP);
315: lmP->D = 0;
316: lmP->M = 0;
317: lmP->Xold = 0;
318: lmP->Gold = 0;
319: lmP->H0 = NULL;
321: tao->data = (void*)lmP;
322: /* Override default settings (unless already changed) */
323: if (!tao->max_it_changed) tao->max_it = 2000;
324: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
326: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
327: TaoLineSearchSetType(tao->linesearch,morethuente_type);
328: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
329: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
330: return(0);
331: }