xref: /petsc/src/tao/unconstrained/impls/lmvm/lmvm.c (revision d5ae23803a94ca2b9c66785d81449705f0cc1e9e) 
1 #include <petsctaolinesearch.h>
2 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>
3 #include <../src/tao/bound/impls/blmvm/blmvm.h>
4 
5 #define LMVM_STEP_BFGS     0
6 #define LMVM_STEP_GRAD     1
7 
8 static PetscErrorCode TaoSolve_LMVM(Tao tao)
9 {
10   TAO_LMVM                     *lmP = (TAO_LMVM *)tao->data;
11   PetscReal                    f, fold, gdx, gnorm;
12   PetscReal                    step = 1.0;
13   PetscErrorCode               ierr;
14   PetscInt                     stepType = LMVM_STEP_GRAD, nupdates;
15   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
16 
17   PetscFunctionBegin;
18 
19   if (tao->XL || tao->XU || tao->ops->computebounds) {
20     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");CHKERRQ(ierr);
21   }
22 
23   /*  Check convergence criteria */
24   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
25   ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
26 
27   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
28 
29   tao->reason = TAO_CONTINUE_ITERATING;
30   ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
31   ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);CHKERRQ(ierr);
32   ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
33   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
34 
35   /*  Set counter for gradient/reset steps */
36   if (!lmP->recycle) {
37     lmP->bfgs = 0;
38     lmP->grad = 0;
39     ierr = MatLMVMReset(lmP->M, PETSC_FALSE); CHKERRQ(ierr);
40   }
41 
42   /*  Have not converged; continue with Newton method */
43   while (tao->reason == TAO_CONTINUE_ITERATING) {
44     /*  Compute direction */
45     if (lmP->H0) {
46       ierr = MatLMVMSetJ0(lmP->M, lmP->H0);CHKERRQ(ierr);
47       stepType = LMVM_STEP_BFGS;
48     }
49     ierr = MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr);
50     ierr = MatSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr);
51     ierr = MatLMVMGetUpdateCount(lmP->M, &nupdates); CHKERRQ(ierr);
52     if (nupdates > 0) stepType = LMVM_STEP_BFGS;
53 
54     /*  Check for success (descent direction) */
55     ierr = VecDot(lmP->D, tao->gradient, &gdx);CHKERRQ(ierr);
56     if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
57       /* Step is not descent or direction produced not a number
58          We can assert bfgsUpdates > 1 in this case because
59          the first solve produces the scaled gradient direction,
60          which is guaranteed to be descent
61 
62          Use steepest descent direction (scaled)
63       */
64 
65       ierr = MatLMVMReset(lmP->M, PETSC_FALSE);CHKERRQ(ierr);
66       ierr = MatLMVMClearJ0(lmP->M);CHKERRQ(ierr);
67       ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr);
68       ierr = MatSolve(lmP->M,tao->gradient, lmP->D);CHKERRQ(ierr);
69 
70       /* On a reset, the direction cannot be not a number; it is a
71          scaled gradient step.  No need to check for this condition. */
72       stepType = LMVM_STEP_GRAD;
73     }
74     ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr);
75 
76     /*  Perform the linesearch */
77     fold = f;
78     ierr = VecCopy(tao->solution, lmP->Xold);CHKERRQ(ierr);
79     ierr = VecCopy(tao->gradient, lmP->Gold);CHKERRQ(ierr);
80 
81     ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);CHKERRQ(ierr);
82     ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
83 
84     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) {
85       /*  Reset factors and use scaled gradient step */
86       f = fold;
87       ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr);
88       ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr);
89 
90       /*  Failed to obtain acceptable iterate with BFGS step */
91       /*  Attempt to use the scaled gradient direction */
92 
93       ierr = MatLMVMReset(lmP->M, PETSC_FALSE);CHKERRQ(ierr);
94       ierr = MatLMVMClearJ0(lmP->M);CHKERRQ(ierr);
95       ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr);
96       ierr = MatSolve(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr);
97 
98       /* On a reset, the direction cannot be not a number; it is a
99           scaled gradient step.  No need to check for this condition. */
100       stepType = LMVM_STEP_GRAD;
101       ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr);
102 
103       /*  Perform the linesearch */
104       ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);CHKERRQ(ierr);
105       ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
106     }
107 
108     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
109       /*  Failed to find an improving point */
110       f = fold;
111       ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr);
112       ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr);
113       step = 0.0;
114       tao->reason = TAO_DIVERGED_LS_FAILURE;
115     } else {
116       /* LS found valid step, so tally up step type */
117       switch (stepType) {
118       case LMVM_STEP_BFGS:
119         ++lmP->bfgs;
120         break;
121       case LMVM_STEP_GRAD:
122         ++lmP->grad;
123         break;
124       default:
125         break;
126       }
127       /*  Compute new gradient norm */
128       ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr);
129     }
130 
131     /* Check convergence */
132     tao->niter++;
133     ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr);
134     ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);CHKERRQ(ierr);
135     ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr);
136   }
137   PetscFunctionReturn(0);
138 }
139 
140 static PetscErrorCode TaoSetUp_LMVM(Tao tao)
141 {
142   TAO_LMVM       *lmP = (TAO_LMVM *)tao->data;
143   PetscInt       n,N;
144   PetscErrorCode ierr;
145   PetscBool      is_spd;
146 
147   PetscFunctionBegin;
148   /* Existence of tao->solution checked in TaoSetUp() */
149   if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);  }
150   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr);  }
151   if (!lmP->D) {ierr = VecDuplicate(tao->solution,&lmP->D);CHKERRQ(ierr);  }
152   if (!lmP->Xold) {ierr = VecDuplicate(tao->solution,&lmP->Xold);CHKERRQ(ierr);  }
153   if (!lmP->Gold) {ierr = VecDuplicate(tao->solution,&lmP->Gold);CHKERRQ(ierr);  }
154 
155   /*  Create matrix for the limited memory approximation */
156   ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
157   ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
158   ierr = MatSetSizes(lmP->M, n, n, N, N);CHKERRQ(ierr);
159   ierr = MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr);
160   ierr = MatGetOption(lmP->M, MAT_SPD, &is_spd);CHKERRQ(ierr);
161   if (!is_spd) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite.");
162 
163   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
164   if (lmP->H0) {
165     ierr = MatLMVMSetJ0(lmP->M, lmP->H0);CHKERRQ(ierr);
166   }
167 
168   PetscFunctionReturn(0);
169 }
170 
171 /* ---------------------------------------------------------- */
172 static PetscErrorCode TaoDestroy_LMVM(Tao tao)
173 {
174   TAO_LMVM       *lmP = (TAO_LMVM *)tao->data;
175   PetscErrorCode ierr;
176 
177   PetscFunctionBegin;
178   if (tao->setupcalled) {
179     ierr = VecDestroy(&lmP->Xold);CHKERRQ(ierr);
180     ierr = VecDestroy(&lmP->Gold);CHKERRQ(ierr);
181     ierr = VecDestroy(&lmP->D);CHKERRQ(ierr);
182   }
183   ierr = MatDestroy(&lmP->M);CHKERRQ(ierr);
184   if (lmP->H0) {
185     ierr = PetscObjectDereference((PetscObject)lmP->H0);CHKERRQ(ierr);
186   }
187   ierr = PetscFree(tao->data);CHKERRQ(ierr);
188 
189   PetscFunctionReturn(0);
190 }
191 
192 /*------------------------------------------------------------*/
193 static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject,Tao tao)
194 {
195   TAO_LMVM       *lm = (TAO_LMVM *)tao->data;
196   PetscErrorCode ierr;
197 
198   PetscFunctionBegin;
199   ierr = PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");CHKERRQ(ierr);
200   ierr = PetscOptionsBool("-tao_lmvm_recycle","enable recycling of the BFGS matrix between subsequent TaoSolve() calls","",lm->recycle,&lm->recycle,NULL);CHKERRQ(ierr);
201   ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
202   ierr = MatSetFromOptions(lm->M);CHKERRQ(ierr);
203   ierr = PetscOptionsTail();CHKERRQ(ierr);
204   PetscFunctionReturn(0);
205 }
206 
207 /*------------------------------------------------------------*/
208 static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
209 {
210   TAO_LMVM       *lm = (TAO_LMVM *)tao->data;
211   PetscBool      isascii;
212   PetscInt       recycled_its;
213   PetscErrorCode ierr;
214 
215   PetscFunctionBegin;
216   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
217   if (isascii) {
218     ierr = PetscViewerASCIIPrintf(viewer, "  Gradient steps: %D\n", lm->grad);CHKERRQ(ierr);
219     if (lm->recycle) {
220       ierr = PetscViewerASCIIPrintf(viewer, "  Recycle: on\n");CHKERRQ(ierr);
221       recycled_its = lm->bfgs + lm->grad;
222       ierr = PetscViewerASCIIPrintf(viewer, "  Total recycled iterations: %D\n", recycled_its);CHKERRQ(ierr);
223     }
224   }
225   PetscFunctionReturn(0);
226 }
227 
228 /* ---------------------------------------------------------- */
229 
230 /*MC
231   TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
232   optimization solver for unconstrained minimization. It solves
233   the Newton step
234           Hkdk = - gk
235 
236   using an approximation Bk in place of Hk, where Bk is composed using
237   the BFGS update formula. A More-Thuente line search is then used
238   to computed the steplength in the dk direction
239 
240   Options Database Keys:
241 .   -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls
242 .   -tao_lmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation
243 
244   Level: beginner
245 M*/
246 
247 PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
248 {
249   TAO_LMVM       *lmP;
250   const char     *morethuente_type = TAOLINESEARCHMT;
251   PetscErrorCode ierr;
252 
253   PetscFunctionBegin;
254   tao->ops->setup = TaoSetUp_LMVM;
255   tao->ops->solve = TaoSolve_LMVM;
256   tao->ops->view = TaoView_LMVM;
257   tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
258   tao->ops->destroy = TaoDestroy_LMVM;
259 
260   ierr = PetscNewLog(tao,&lmP);CHKERRQ(ierr);
261   lmP->D = 0;
262   lmP->M = 0;
263   lmP->Xold = 0;
264   lmP->Gold = 0;
265   lmP->H0   = NULL;
266   lmP->recycle = PETSC_FALSE;
267 
268   tao->data = (void*)lmP;
269   /* Override default settings (unless already changed) */
270   if (!tao->max_it_changed) tao->max_it = 2000;
271   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
272 
273   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr);
274   ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr);
275   ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr);
276   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr);
277   ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr);
278 
279   ierr = KSPInitializePackage();CHKERRQ(ierr);
280   ierr = MatCreate(((PetscObject)tao)->comm, &lmP->M);CHKERRQ(ierr);
281   ierr = PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1);CHKERRQ(ierr);
282   ierr = MatSetType(lmP->M, MATLMVMBFGS);CHKERRQ(ierr);
283   ierr = MatSetOptionsPrefix(lmP->M, "tao_lmvm_");CHKERRQ(ierr);
284   PetscFunctionReturn(0);
285 }
286 
287 PETSC_EXTERN PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
288 {
289   TAO_LMVM       *lmP;
290   TAO_BLMVM      *blmP;
291   TaoType        type;
292   PetscBool      is_lmvm, is_blmvm;
293   PetscErrorCode ierr;
294 
295   ierr = TaoGetType(tao, &type);CHKERRQ(ierr);
296   ierr = PetscStrcmp(type, TAOLMVM,  &is_lmvm);CHKERRQ(ierr);
297   ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr);
298 
299   if (is_lmvm) {
300     lmP = (TAO_LMVM *)tao->data;
301     ierr = PetscObjectReference((PetscObject)H0);CHKERRQ(ierr);
302     lmP->H0 = H0;
303   } else if (is_blmvm) {
304     blmP = (TAO_BLMVM *)tao->data;
305     ierr = PetscObjectReference((PetscObject)H0);CHKERRQ(ierr);
306     blmP->H0 = H0;
307   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");
308   PetscFunctionReturn(0);
309 }
310 
311 PETSC_EXTERN PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
312 {
313   TAO_LMVM       *lmP;
314   TAO_BLMVM      *blmP;
315   TaoType        type;
316   PetscBool      is_lmvm, is_blmvm;
317   Mat            M;
318 
319   PetscErrorCode ierr;
320 
321   ierr = TaoGetType(tao, &type);CHKERRQ(ierr);
322   ierr = PetscStrcmp(type, TAOLMVM,  &is_lmvm);CHKERRQ(ierr);
323   ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr);
324 
325   if (is_lmvm) {
326     lmP = (TAO_LMVM *)tao->data;
327     M = lmP->M;
328   } else if (is_blmvm) {
329     blmP = (TAO_BLMVM *)tao->data;
330     M = blmP->M;
331   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");
332   ierr = MatLMVMGetJ0(M, H0);CHKERRQ(ierr);
333   PetscFunctionReturn(0);
334 }
335 
336 PETSC_EXTERN PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
337 {
338   TAO_LMVM       *lmP;
339   TAO_BLMVM      *blmP;
340   TaoType        type;
341   PetscBool      is_lmvm, is_blmvm;
342   Mat            M;
343   PetscErrorCode ierr;
344 
345   ierr = TaoGetType(tao, &type);CHKERRQ(ierr);
346   ierr = PetscStrcmp(type, TAOLMVM,  &is_lmvm);CHKERRQ(ierr);
347   ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr);
348 
349   if (is_lmvm) {
350     lmP = (TAO_LMVM *)tao->data;
351     M = lmP->M;
352   } else if (is_blmvm) {
353     blmP = (TAO_BLMVM *)tao->data;
354     M = blmP->M;
355   } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");
356   ierr = MatLMVMGetJ0KSP(M, ksp);CHKERRQ(ierr);
357   PetscFunctionReturn(0);
358 }