Actual source code: owlqn.c
petsc-3.11.4 2019-09-28
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/unconstrained/impls/owlqn/owlqn.h>
4: #define OWLQN_BFGS 0
5: #define OWLQN_SCALED_GRADIENT 1
6: #define OWLQN_GRADIENT 2
8: static PetscErrorCode ProjDirect_OWLQN(Vec d, Vec g)
9: {
10: PetscErrorCode ierr;
11: const PetscReal *gptr;
12: PetscReal *dptr;
13: PetscInt low,high,low1,high1,i;
16: ierr=VecGetOwnershipRange(d,&low,&high);
17: ierr=VecGetOwnershipRange(g,&low1,&high1);
19: VecGetArrayRead(g,&gptr);
20: VecGetArray(d,&dptr);
21: for (i = 0; i < high-low; i++) {
22: if (dptr[i] * gptr[i] <= 0.0 ) {
23: dptr[i] = 0.0;
24: }
25: }
26: VecRestoreArray(d,&dptr);
27: VecRestoreArrayRead(g,&gptr);
28: return(0);
29: }
31: static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda)
32: {
33: PetscErrorCode ierr;
34: const PetscReal *xptr;
35: PetscReal *gptr;
36: PetscInt low,high,low1,high1,i;
39: ierr=VecGetOwnershipRange(x,&low,&high);
40: ierr=VecGetOwnershipRange(gv,&low1,&high1);
42: VecGetArrayRead(x,&xptr);
43: VecGetArray(gv,&gptr);
44: for (i = 0; i < high-low; i++) {
45: if (xptr[i] < 0.0) gptr[i] = gptr[i] - lambda;
46: else if (xptr[i] > 0.0) gptr[i] = gptr[i] + lambda;
47: else if (gptr[i] + lambda < 0.0) gptr[i] = gptr[i] + lambda;
48: else if (gptr[i] - lambda > 0.0) gptr[i] = gptr[i] - lambda;
49: else gptr[i] = 0.0;
50: }
51: VecRestoreArray(gv,&gptr);
52: VecRestoreArrayRead(x,&xptr);
53: return(0);
54: }
56: static PetscErrorCode TaoSolve_OWLQN(Tao tao)
57: {
58: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
59: PetscReal f, fold, gdx, gnorm;
60: PetscReal step = 1.0;
61: PetscReal delta;
62: PetscErrorCode ierr;
63: PetscInt stepType;
64: PetscInt iter = 0;
65: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
68: if (tao->XL || tao->XU || tao->ops->computebounds) {
69: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by owlqn algorithm\n");
70: }
72: /* Check convergence criteria */
73: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
75: VecCopy(tao->gradient, lmP->GV);
77: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
79: VecNorm(lmP->GV,NORM_2,&gnorm);
81: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
83: tao->reason = TAO_CONTINUE_ITERATING;
84: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
85: TaoMonitor(tao,iter,f,gnorm,0.0,step);
86: (*tao->ops->convergencetest)(tao,tao->cnvP);
87: if (tao->reason != TAO_CONTINUE_ITERATING) return(0);
89: /* Set initial scaling for the function */
90: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
91: MatLMVMSetJ0Scale(lmP->M, delta);
93: /* Set counter for gradient/reset steps */
94: lmP->bfgs = 0;
95: lmP->sgrad = 0;
96: lmP->grad = 0;
98: /* Have not converged; continue with Newton method */
99: while (tao->reason == TAO_CONTINUE_ITERATING) {
100: /* Call general purpose update function */
101: if (tao->ops->update) {
102: (*tao->ops->update)(tao, tao->niter, tao->user_update);
103: }
105: /* Compute direction */
106: MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
107: MatSolve(lmP->M, lmP->GV, lmP->D);
109: ProjDirect_OWLQN(lmP->D,lmP->GV);
111: ++lmP->bfgs;
113: /* Check for success (descent direction) */
114: VecDot(lmP->D, lmP->GV , &gdx);
115: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
117: /* Step is not descent or direction produced not a number
118: We can assert bfgsUpdates > 1 in this case because
119: the first solve produces the scaled gradient direction,
120: which is guaranteed to be descent
122: Use steepest descent direction (scaled) */
123: ++lmP->grad;
125: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
126: MatLMVMSetJ0Scale(lmP->M, delta);
127: MatLMVMReset(lmP->M, PETSC_FALSE);
128: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
129: MatSolve(lmP->M,lmP->GV, lmP->D);
131: ProjDirect_OWLQN(lmP->D,lmP->GV);
133: lmP->bfgs = 1;
134: ++lmP->sgrad;
135: stepType = OWLQN_SCALED_GRADIENT;
136: } else {
137: if (1 == lmP->bfgs) {
138: /* The first BFGS direction is always the scaled gradient */
139: ++lmP->sgrad;
140: stepType = OWLQN_SCALED_GRADIENT;
141: } else {
142: ++lmP->bfgs;
143: stepType = OWLQN_BFGS;
144: }
145: }
147: VecScale(lmP->D, -1.0);
149: /* Perform the linesearch */
150: fold = f;
151: VecCopy(tao->solution, lmP->Xold);
152: VecCopy(tao->gradient, lmP->Gold);
154: TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step,&ls_status);
155: TaoAddLineSearchCounts(tao);
157: while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) {
159: /* Reset factors and use scaled gradient step */
160: f = fold;
161: VecCopy(lmP->Xold, tao->solution);
162: VecCopy(lmP->Gold, tao->gradient);
163: VecCopy(tao->gradient, lmP->GV);
165: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
167: switch(stepType) {
168: case OWLQN_BFGS:
169: /* Failed to obtain acceptable iterate with BFGS step
170: Attempt to use the scaled gradient direction */
172: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
173: MatLMVMSetJ0Scale(lmP->M, delta);
174: MatLMVMReset(lmP->M, PETSC_FALSE);
175: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
176: MatSolve(lmP->M, lmP->GV, lmP->D);
178: ProjDirect_OWLQN(lmP->D,lmP->GV);
180: lmP->bfgs = 1;
181: ++lmP->sgrad;
182: stepType = OWLQN_SCALED_GRADIENT;
183: break;
185: case OWLQN_SCALED_GRADIENT:
186: /* The scaled gradient step did not produce a new iterate;
187: attempt to use the gradient direction.
188: Need to make sure we are not using a different diagonal scaling */
189: MatLMVMSetJ0Scale(lmP->M, 1.0);
190: MatLMVMReset(lmP->M, PETSC_FALSE);
191: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
192: MatSolve(lmP->M, lmP->GV, lmP->D);
194: ProjDirect_OWLQN(lmP->D,lmP->GV);
196: lmP->bfgs = 1;
197: ++lmP->grad;
198: stepType = OWLQN_GRADIENT;
199: break;
200: }
201: VecScale(lmP->D, -1.0);
203: /* Perform the linesearch */
204: TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
205: TaoAddLineSearchCounts(tao);
206: }
208: if ((int)ls_status < 0) {
209: /* Failed to find an improving point*/
210: f = fold;
211: VecCopy(lmP->Xold, tao->solution);
212: VecCopy(lmP->Gold, tao->gradient);
213: VecCopy(tao->gradient, lmP->GV);
214: step = 0.0;
215: } else {
216: /* a little hack here, because that gv is used to store g */
217: VecCopy(lmP->GV, tao->gradient);
218: }
220: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
222: /* Check for termination */
224: VecNorm(lmP->GV,NORM_2,&gnorm);
226: iter++;
227: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
228: TaoMonitor(tao,iter,f,gnorm,0.0,step);
229: (*tao->ops->convergencetest)(tao,tao->cnvP);
231: if ((int)ls_status < 0) break;
232: }
233: return(0);
234: }
236: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
237: {
238: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
239: PetscInt n,N;
243: /* Existence of tao->solution checked in TaoSetUp() */
244: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
245: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
246: if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D); }
247: if (!lmP->GV) {VecDuplicate(tao->solution,&lmP->GV); }
248: if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold); }
249: if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold); }
251: /* Create matrix for the limited memory approximation */
252: VecGetLocalSize(tao->solution,&n);
253: VecGetSize(tao->solution,&N);
254: MatCreateLMVMBFGS(((PetscObject)tao)->comm,n,N,&lmP->M);
255: MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);
256: return(0);
257: }
259: /* ---------------------------------------------------------- */
260: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
261: {
262: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
266: if (tao->setupcalled) {
267: VecDestroy(&lmP->Xold);
268: VecDestroy(&lmP->Gold);
269: VecDestroy(&lmP->D);
270: MatDestroy(&lmP->M);
271: VecDestroy(&lmP->GV);
272: }
273: PetscFree(tao->data);
274: return(0);
275: }
277: /*------------------------------------------------------------*/
278: static PetscErrorCode TaoSetFromOptions_OWLQN(PetscOptionItems *PetscOptionsObject,Tao tao)
279: {
280: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
284: PetscOptionsHead(PetscOptionsObject,"Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
285: PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,NULL);
286: PetscOptionsTail();
287: TaoLineSearchSetFromOptions(tao->linesearch);
288: return(0);
289: }
291: /*------------------------------------------------------------*/
292: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
293: {
294: TAO_OWLQN *lm = (TAO_OWLQN *)tao->data;
295: PetscBool isascii;
299: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
300: if (isascii) {
301: PetscViewerASCIIPushTab(viewer);
302: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
303: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
304: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
305: PetscViewerASCIIPopTab(viewer);
306: }
307: return(0);
308: }
310: /* ---------------------------------------------------------- */
311: /*MC
312: TAOOWLQN - orthant-wise limited memory quasi-newton algorithm
314: . - tao_owlqn_lambda - regulariser weight
316: Level: beginner
317: M*/
320: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
321: {
322: TAO_OWLQN *lmP;
323: const char *owarmijo_type = TAOLINESEARCHOWARMIJO;
327: tao->ops->setup = TaoSetUp_OWLQN;
328: tao->ops->solve = TaoSolve_OWLQN;
329: tao->ops->view = TaoView_OWLQN;
330: tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
331: tao->ops->destroy = TaoDestroy_OWLQN;
333: PetscNewLog(tao,&lmP);
334: lmP->D = 0;
335: lmP->M = 0;
336: lmP->GV = 0;
337: lmP->Xold = 0;
338: lmP->Gold = 0;
339: lmP->lambda = 1.0;
341: tao->data = (void*)lmP;
342: /* Override default settings (unless already changed) */
343: if (!tao->max_it_changed) tao->max_it = 2000;
344: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
346: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
347: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
348: TaoLineSearchSetType(tao->linesearch,owarmijo_type);
349: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
350: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
351: return(0);
352: }