xref: /petsc/src/tao/unconstrained/impls/owlqn/owlqn.c (revision 30f8f9828c55eeb8db61ae34754e9b56df1e9b4e) 
1 #include <petsctaolinesearch.h>
2 #include <../src/tao/matrix/lmvmmat.h>
3 #include <../src/tao/unconstrained/impls/owlqn/owlqn.h>
4 
5 #define OWLQN_BFGS                0
6 #define OWLQN_SCALED_GRADIENT     1
7 #define OWLQN_GRADIENT            2
8 
9 #undef __FUNCT__
10 #define __FUNCT__ "ProjDirect_OWLQN"
11 static PetscErrorCode ProjDirect_OWLQN(Vec d, Vec g)
12 {
13   PetscErrorCode ierr;
14   PetscReal      *gptr,*dptr;
15   PetscInt       low,high,low1,high1,i;
16 
17   PetscFunctionBegin;
18   ierr=VecGetOwnershipRange(d,&low,&high);CHKERRQ(ierr);
19   ierr=VecGetOwnershipRange(g,&low1,&high1);CHKERRQ(ierr);
20 
21   ierr = VecGetArray(g,&gptr);CHKERRQ(ierr);
22   ierr = VecGetArray(d,&dptr);CHKERRQ(ierr);
23   for (i = 0; i < high-low; i++) {
24     if (dptr[i] * gptr[i] <= 0.0 ) {
25       dptr[i] = 0.0;
26     }
27   }
28   ierr = VecRestoreArray(d,&dptr);CHKERRQ(ierr);
29   ierr = VecRestoreArray(g,&gptr);CHKERRQ(ierr);
30   PetscFunctionReturn(0);
31 }
32 
33 #undef __FUNCT__
34 #define __FUNCT__ "ComputePseudoGrad_OWLQN"
35 static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda)
36 {
37   PetscErrorCode ierr;
38   PetscReal      *xptr,*gptr;
39   PetscInt       low,high,low1,high1,i;
40 
41   PetscFunctionBegin;
42   ierr=VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr);
43   ierr=VecGetOwnershipRange(gv,&low1,&high1);CHKERRQ(ierr);
44 
45   ierr = VecGetArray(x,&xptr);CHKERRQ(ierr);
46   ierr = VecGetArray(gv,&gptr);CHKERRQ(ierr);
47   for (i = 0; i < high-low; i++) {
48     if (xptr[i] < 0.0)               gptr[i] = gptr[i] - lambda;
49     else if (xptr[i] > 0.0)          gptr[i] = gptr[i] + lambda;
50     else if (gptr[i] + lambda < 0.0) gptr[i] = gptr[i] + lambda;
51     else if (gptr[i] - lambda > 0.0) gptr[i] = gptr[i] - lambda;
52     else                             gptr[i] = 0.0;
53   }
54   ierr = VecRestoreArray(gv,&gptr);CHKERRQ(ierr);
55   ierr = VecRestoreArray(x,&xptr);CHKERRQ(ierr);
56   PetscFunctionReturn(0);
57 }
58 
59 #undef __FUNCT__
60 #define __FUNCT__ "TaoSolve_OWLQN"
61 static PetscErrorCode TaoSolve_OWLQN(Tao tao)
62 {
63   TAO_OWLQN                      *lmP = (TAO_OWLQN *)tao->data;
64   PetscReal                      f, fold, gdx, gnorm;
65   PetscReal                      step = 1.0;
66   PetscReal                      delta;
67   PetscErrorCode                 ierr;
68   PetscInt                       stepType;
69   PetscInt                       iter = 0;
70   TaoTerminationReason     reason = TAO_CONTINUE_ITERATING;
71   TaoLineSearchTerminationReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
72 
73   PetscFunctionBegin;
74 
75   if (tao->XL || tao->XU || tao->ops->computebounds) {
76     ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by owlqn algorithm\n");CHKERRQ(ierr);
77   }
78 
79   /* Check convergence criteria */
80   ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr);
81 
82   ierr = VecCopy(tao->gradient, lmP->GV);CHKERRQ(ierr);
83 
84   ierr = ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);CHKERRQ(ierr);
85 
86   ierr = VecNorm(lmP->GV,NORM_2,&gnorm);CHKERRQ(ierr);
87 
88   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
89 
90   ierr = TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);CHKERRQ(ierr);
91   if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0);
92 
93   /* Set initial scaling for the function */
94   if (f != 0.0) {
95     delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
96   } else {
97     delta = 2.0 / (gnorm*gnorm);
98   }
99   ierr = MatLMVMSetDelta(lmP->M,delta);CHKERRQ(ierr);
100 
101   /* Set counter for gradient/reset steps */
102   lmP->bfgs = 0;
103   lmP->sgrad = 0;
104   lmP->grad = 0;
105 
106   /* Have not converged; continue with Newton method */
107   while (reason == TAO_CONTINUE_ITERATING) {
108     /* Compute direction */
109     ierr = MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr);
110     ierr = MatLMVMSolve(lmP->M, lmP->GV, lmP->D);CHKERRQ(ierr);
111 
112     ierr = ProjDirect_OWLQN(lmP->D,lmP->GV);CHKERRQ(ierr);
113 
114     ++lmP->bfgs;
115 
116     /* Check for success (descent direction) */
117     ierr = VecDot(lmP->D, lmP->GV , &gdx);CHKERRQ(ierr);
118     if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
119 
120       /* Step is not descent or direction produced not a number
121          We can assert bfgsUpdates > 1 in this case because
122          the first solve produces the scaled gradient direction,
123          which is guaranteed to be descent
124 
125          Use steepest descent direction (scaled) */
126       ++lmP->grad;
127 
128       if (f != 0.0) {
129         delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
130       } else {
131         delta = 2.0 / (gnorm*gnorm);
132       }
133       ierr = MatLMVMSetDelta(lmP->M, delta);CHKERRQ(ierr);
134       ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr);
135       ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr);
136       ierr = MatLMVMSolve(lmP->M,lmP->GV, lmP->D);CHKERRQ(ierr);
137 
138       ierr = ProjDirect_OWLQN(lmP->D,lmP->GV);CHKERRQ(ierr);
139 
140       lmP->bfgs = 1;
141       ++lmP->sgrad;
142       stepType = OWLQN_SCALED_GRADIENT;
143     } else {
144       if (1 == lmP->bfgs) {
145         /* The first BFGS direction is always the scaled gradient */
146         ++lmP->sgrad;
147         stepType = OWLQN_SCALED_GRADIENT;
148       } else {
149         ++lmP->bfgs;
150         stepType = OWLQN_BFGS;
151       }
152     }
153 
154     ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr);
155 
156     /* Perform the linesearch */
157     fold = f;
158     ierr = VecCopy(tao->solution, lmP->Xold);CHKERRQ(ierr);
159     ierr = VecCopy(tao->gradient, lmP->Gold);CHKERRQ(ierr);
160 
161     ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step,&ls_status);CHKERRQ(ierr);
162     ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
163 
164     while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) {
165 
166       /* Reset factors and use scaled gradient step */
167       f = fold;
168       ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr);
169       ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr);
170       ierr = VecCopy(tao->gradient, lmP->GV);CHKERRQ(ierr);
171 
172       ierr = ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);CHKERRQ(ierr);
173 
174       switch(stepType) {
175       case OWLQN_BFGS:
176         /* Failed to obtain acceptable iterate with BFGS step
177            Attempt to use the scaled gradient direction */
178 
179         if (f != 0.0) {
180           delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
181         } else {
182           delta = 2.0 / (gnorm*gnorm);
183         }
184         ierr = MatLMVMSetDelta(lmP->M, delta);CHKERRQ(ierr);
185         ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr);
186         ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr);
187         ierr = MatLMVMSolve(lmP->M, lmP->GV, lmP->D);CHKERRQ(ierr);
188 
189         ierr = ProjDirect_OWLQN(lmP->D,lmP->GV);CHKERRQ(ierr);
190 
191         lmP->bfgs = 1;
192         ++lmP->sgrad;
193         stepType = OWLQN_SCALED_GRADIENT;
194         break;
195 
196       case OWLQN_SCALED_GRADIENT:
197         /* The scaled gradient step did not produce a new iterate;
198            attempt to use the gradient direction.
199            Need to make sure we are not using a different diagonal scaling */
200         ierr = MatLMVMSetDelta(lmP->M, 1.0);CHKERRQ(ierr);
201         ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr);
202         ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr);
203         ierr = MatLMVMSolve(lmP->M, lmP->GV, lmP->D);CHKERRQ(ierr);
204 
205         ierr = ProjDirect_OWLQN(lmP->D,lmP->GV);CHKERRQ(ierr);
206 
207         lmP->bfgs = 1;
208         ++lmP->grad;
209         stepType = OWLQN_GRADIENT;
210         break;
211       }
212       ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr);
213 
214 
215       /* Perform the linesearch */
216       ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);CHKERRQ(ierr);
217       ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr);
218     }
219 
220     if ((int)ls_status < 0) {
221       /* Failed to find an improving point*/
222       f = fold;
223       ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr);
224       ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr);
225       ierr = VecCopy(tao->gradient, lmP->GV);CHKERRQ(ierr);
226       step = 0.0;
227     } else {
228       /* a little hack here, because that gv is used to store g */
229       ierr = VecCopy(lmP->GV, tao->gradient);CHKERRQ(ierr);
230     }
231 
232     ierr = ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);CHKERRQ(ierr);
233 
234     /* Check for termination */
235 
236     ierr = VecNorm(lmP->GV,NORM_2,&gnorm);CHKERRQ(ierr);
237 
238     iter++;
239     ierr = TaoMonitor(tao,iter,f,gnorm,0.0,step,&reason);CHKERRQ(ierr);
240 
241     if ((int)ls_status < 0) break;
242   }
243   PetscFunctionReturn(0);
244 }
245 
246 #undef __FUNCT__
247 #define __FUNCT__ "TaoSetUp_OWLQN"
248 static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
249 {
250   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;
251   PetscInt       n,N;
252   PetscErrorCode ierr;
253 
254   PetscFunctionBegin;
255   /* Existence of tao->solution checked in TaoSetUp() */
256   if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr);  }
257   if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr);  }
258   if (!lmP->D) {ierr = VecDuplicate(tao->solution,&lmP->D);CHKERRQ(ierr);  }
259   if (!lmP->GV) {ierr = VecDuplicate(tao->solution,&lmP->GV);CHKERRQ(ierr);  }
260   if (!lmP->Xold) {ierr = VecDuplicate(tao->solution,&lmP->Xold);CHKERRQ(ierr);  }
261   if (!lmP->Gold) {ierr = VecDuplicate(tao->solution,&lmP->Gold);CHKERRQ(ierr);  }
262 
263   /* Create matrix for the limited memory approximation */
264   ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr);
265   ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr);
266   ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M);CHKERRQ(ierr);
267   ierr = MatLMVMAllocateVectors(lmP->M,tao->solution);CHKERRQ(ierr);
268   PetscFunctionReturn(0);
269 }
270 
271 /* ---------------------------------------------------------- */
272 #undef __FUNCT__
273 #define __FUNCT__ "TaoDestroy_OWLQN"
274 static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
275 {
276   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;
277   PetscErrorCode ierr;
278 
279   PetscFunctionBegin;
280   if (tao->setupcalled) {
281     ierr = VecDestroy(&lmP->Xold);CHKERRQ(ierr);
282     ierr = VecDestroy(&lmP->Gold);CHKERRQ(ierr);
283     ierr = VecDestroy(&lmP->D);CHKERRQ(ierr);
284     ierr = MatDestroy(&lmP->M);CHKERRQ(ierr);
285     ierr = VecDestroy(&lmP->GV);CHKERRQ(ierr);
286   }
287   ierr = PetscFree(tao->data);CHKERRQ(ierr);
288   PetscFunctionReturn(0);
289 }
290 
291 /*------------------------------------------------------------*/
292 #undef __FUNCT__
293 #define __FUNCT__ "TaoSetFromOptions_OWLQN"
294 static PetscErrorCode TaoSetFromOptions_OWLQN(Tao tao)
295 {
296   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;
297   PetscBool      flg;
298   PetscErrorCode ierr;
299 
300   PetscFunctionBegin;
301   ierr = PetscOptionsHead("Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");CHKERRQ(ierr);
302   ierr = PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,&flg); CHKERRQ(ierr);
303   ierr = PetscOptionsTail();CHKERRQ(ierr);
304   ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr);
305   PetscFunctionReturn(0);
306 }
307 
308 /*------------------------------------------------------------*/
309 #undef __FUNCT__
310 #define __FUNCT__ "TaoView_OWLQN"
311 static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
312 {
313   TAO_OWLQN      *lm = (TAO_OWLQN *)tao->data;
314   PetscBool      isascii;
315   PetscErrorCode ierr;
316 
317   PetscFunctionBegin;
318   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr);
319   if (isascii) {
320     ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr);
321     ierr = PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);CHKERRQ(ierr);
322     ierr = PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);CHKERRQ(ierr);
323     ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);CHKERRQ(ierr);
324     ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr);
325   }
326   PetscFunctionReturn(0);
327 }
328 
329 /* ---------------------------------------------------------- */
330 
331 EXTERN_C_BEGIN
332 #undef __FUNCT__
333 #define __FUNCT__ "TaoCreate_OWLQN"
334 PetscErrorCode TaoCreate_OWLQN(Tao tao)
335 {
336   TAO_OWLQN      *lmP;
337   const char     *owarmijo_type = TAOLINESEARCH_OWARMIJO;
338   PetscErrorCode ierr;
339 
340   PetscFunctionBegin;
341   tao->ops->setup = TaoSetUp_OWLQN;
342   tao->ops->solve = TaoSolve_OWLQN;
343   tao->ops->view = TaoView_OWLQN;
344   tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
345   tao->ops->destroy = TaoDestroy_OWLQN;
346 
347   ierr = PetscNewLog(tao,&lmP);CHKERRQ(ierr);
348   lmP->D = 0;
349   lmP->M = 0;
350   lmP->GV = 0;
351   lmP->Xold = 0;
352   lmP->Gold = 0;
353   lmP->lambda = 1.0;
354 
355   tao->data = (void*)lmP;
356   tao->max_it = 2000;
357   tao->max_funcs = 4000;
358   tao->fatol = 1e-4;
359   tao->frtol = 1e-4;
360 
361   ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr);
362   ierr = TaoLineSearchSetType(tao->linesearch,owarmijo_type);CHKERRQ(ierr);
363   ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr);
364   PetscFunctionReturn(0);
365 }
366 EXTERN_C_END
367 
368