Actual source code: owlqn.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/owlqn/owlqn.h>
5: #define OWLQN_BFGS 0
6: #define OWLQN_SCALED_GRADIENT 1
7: #define OWLQN_GRADIENT 2
9: static PetscErrorCode ProjDirect_OWLQN(Vec d, Vec g)
10: {
11: PetscErrorCode ierr;
12: const PetscReal *gptr;
13: PetscReal *dptr;
14: PetscInt low,high,low1,high1,i;
17: ierr=VecGetOwnershipRange(d,&low,&high);
18: ierr=VecGetOwnershipRange(g,&low1,&high1);
20: VecGetArrayRead(g,&gptr);
21: VecGetArray(d,&dptr);
22: for (i = 0; i < high-low; i++) {
23: if (dptr[i] * gptr[i] <= 0.0 ) {
24: dptr[i] = 0.0;
25: }
26: }
27: VecRestoreArray(d,&dptr);
28: VecRestoreArrayRead(g,&gptr);
29: return(0);
30: }
32: static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda)
33: {
34: PetscErrorCode ierr;
35: const PetscReal *xptr;
36: PetscReal *gptr;
37: PetscInt low,high,low1,high1,i;
40: ierr=VecGetOwnershipRange(x,&low,&high);
41: ierr=VecGetOwnershipRange(gv,&low1,&high1);
43: VecGetArrayRead(x,&xptr);
44: VecGetArray(gv,&gptr);
45: for (i = 0; i < high-low; i++) {
46: if (xptr[i] < 0.0) gptr[i] = gptr[i] - lambda;
47: else if (xptr[i] > 0.0) gptr[i] = gptr[i] + lambda;
48: else if (gptr[i] + lambda < 0.0) gptr[i] = gptr[i] + lambda;
49: else if (gptr[i] - lambda > 0.0) gptr[i] = gptr[i] - lambda;
50: else gptr[i] = 0.0;
51: }
52: VecRestoreArray(gv,&gptr);
53: VecRestoreArrayRead(x,&xptr);
54: return(0);
55: }
57: static PetscErrorCode TaoSolve_OWLQN(Tao tao)
58: {
59: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
60: PetscReal f, fold, gdx, gnorm;
61: PetscReal step = 1.0;
62: PetscReal delta;
63: PetscErrorCode ierr;
64: PetscInt stepType;
65: PetscInt iter = 0;
66: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
67: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
70: if (tao->XL || tao->XU || tao->ops->computebounds) {
71: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by owlqn algorithm\n");
72: }
74: /* Check convergence criteria */
75: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
77: VecCopy(tao->gradient, lmP->GV);
79: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
81: VecNorm(lmP->GV,NORM_2,&gnorm);
83: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
85: TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);
86: if (reason != TAO_CONTINUE_ITERATING) return(0);
88: /* Set initial scaling for the function */
89: if (f != 0.0) {
90: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
91: } else {
92: delta = 2.0 / (gnorm*gnorm);
93: }
94: MatLMVMSetDelta(lmP->M,delta);
96: /* Set counter for gradient/reset steps */
97: lmP->bfgs = 0;
98: lmP->sgrad = 0;
99: lmP->grad = 0;
101: /* Have not converged; continue with Newton method */
102: while (reason == TAO_CONTINUE_ITERATING) {
103: /* Compute direction */
104: MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
105: MatLMVMSolve(lmP->M, lmP->GV, lmP->D);
107: ProjDirect_OWLQN(lmP->D,lmP->GV);
109: ++lmP->bfgs;
111: /* Check for success (descent direction) */
112: VecDot(lmP->D, lmP->GV , &gdx);
113: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
115: /* Step is not descent or direction produced not a number
116: We can assert bfgsUpdates > 1 in this case because
117: the first solve produces the scaled gradient direction,
118: which is guaranteed to be descent
120: Use steepest descent direction (scaled) */
121: ++lmP->grad;
123: if (f != 0.0) {
124: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
125: } else {
126: delta = 2.0 / (gnorm*gnorm);
127: }
128: MatLMVMSetDelta(lmP->M, delta);
129: MatLMVMReset(lmP->M);
130: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
131: MatLMVMSolve(lmP->M,lmP->GV, lmP->D);
133: ProjDirect_OWLQN(lmP->D,lmP->GV);
135: lmP->bfgs = 1;
136: ++lmP->sgrad;
137: stepType = OWLQN_SCALED_GRADIENT;
138: } else {
139: if (1 == lmP->bfgs) {
140: /* The first BFGS direction is always the scaled gradient */
141: ++lmP->sgrad;
142: stepType = OWLQN_SCALED_GRADIENT;
143: } else {
144: ++lmP->bfgs;
145: stepType = OWLQN_BFGS;
146: }
147: }
149: VecScale(lmP->D, -1.0);
151: /* Perform the linesearch */
152: fold = f;
153: VecCopy(tao->solution, lmP->Xold);
154: VecCopy(tao->gradient, lmP->Gold);
156: TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step,&ls_status);
157: TaoAddLineSearchCounts(tao);
159: while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) {
161: /* Reset factors and use scaled gradient step */
162: f = fold;
163: VecCopy(lmP->Xold, tao->solution);
164: VecCopy(lmP->Gold, tao->gradient);
165: VecCopy(tao->gradient, lmP->GV);
167: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
169: switch(stepType) {
170: case OWLQN_BFGS:
171: /* Failed to obtain acceptable iterate with BFGS step
172: Attempt to use the scaled gradient direction */
174: if (f != 0.0) {
175: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
176: } else {
177: delta = 2.0 / (gnorm*gnorm);
178: }
179: MatLMVMSetDelta(lmP->M, delta);
180: MatLMVMReset(lmP->M);
181: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
182: MatLMVMSolve(lmP->M, lmP->GV, lmP->D);
184: ProjDirect_OWLQN(lmP->D,lmP->GV);
186: lmP->bfgs = 1;
187: ++lmP->sgrad;
188: stepType = OWLQN_SCALED_GRADIENT;
189: break;
191: case OWLQN_SCALED_GRADIENT:
192: /* The scaled gradient step did not produce a new iterate;
193: attempt to use the gradient direction.
194: Need to make sure we are not using a different diagonal scaling */
195: MatLMVMSetDelta(lmP->M, 1.0);
196: MatLMVMReset(lmP->M);
197: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
198: MatLMVMSolve(lmP->M, lmP->GV, lmP->D);
200: ProjDirect_OWLQN(lmP->D,lmP->GV);
202: lmP->bfgs = 1;
203: ++lmP->grad;
204: stepType = OWLQN_GRADIENT;
205: break;
206: }
207: VecScale(lmP->D, -1.0);
210: /* Perform the linesearch */
211: TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
212: TaoAddLineSearchCounts(tao);
213: }
215: if ((int)ls_status < 0) {
216: /* Failed to find an improving point*/
217: f = fold;
218: VecCopy(lmP->Xold, tao->solution);
219: VecCopy(lmP->Gold, tao->gradient);
220: VecCopy(tao->gradient, lmP->GV);
221: step = 0.0;
222: } else {
223: /* a little hack here, because that gv is used to store g */
224: VecCopy(lmP->GV, tao->gradient);
225: }
227: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
229: /* Check for termination */
231: VecNorm(lmP->GV,NORM_2,&gnorm);
233: iter++;
234: TaoMonitor(tao,iter,f,gnorm,0.0,step,&reason);
236: if ((int)ls_status < 0) break;
237: }
238: return(0);
239: }
241: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
242: {
243: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
244: PetscInt n,N;
248: /* Existence of tao->solution checked in TaoSetUp() */
249: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
250: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
251: if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D); }
252: if (!lmP->GV) {VecDuplicate(tao->solution,&lmP->GV); }
253: if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold); }
254: if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold); }
256: /* Create matrix for the limited memory approximation */
257: VecGetLocalSize(tao->solution,&n);
258: VecGetSize(tao->solution,&N);
259: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M);
260: MatLMVMAllocateVectors(lmP->M,tao->solution);
261: return(0);
262: }
264: /* ---------------------------------------------------------- */
265: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
266: {
267: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
271: if (tao->setupcalled) {
272: VecDestroy(&lmP->Xold);
273: VecDestroy(&lmP->Gold);
274: VecDestroy(&lmP->D);
275: MatDestroy(&lmP->M);
276: VecDestroy(&lmP->GV);
277: }
278: PetscFree(tao->data);
279: return(0);
280: }
282: /*------------------------------------------------------------*/
283: static PetscErrorCode TaoSetFromOptions_OWLQN(PetscOptionItems *PetscOptionsObject,Tao tao)
284: {
285: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
289: PetscOptionsHead(PetscOptionsObject,"Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
290: PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,NULL);
291: PetscOptionsTail();
292: TaoLineSearchSetFromOptions(tao->linesearch);
293: return(0);
294: }
296: /*------------------------------------------------------------*/
297: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
298: {
299: TAO_OWLQN *lm = (TAO_OWLQN *)tao->data;
300: PetscBool isascii;
304: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
305: if (isascii) {
306: PetscViewerASCIIPushTab(viewer);
307: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
308: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
309: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
310: PetscViewerASCIIPopTab(viewer);
311: }
312: return(0);
313: }
315: /* ---------------------------------------------------------- */
316: /*MC
317: TAOOWLQN - orthant-wise limited memory quasi-newton algorithm
319: . - tao_owlqn_lambda - regulariser weight
321: Level: beginner
322: M*/
325: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
326: {
327: TAO_OWLQN *lmP;
328: const char *owarmijo_type = TAOLINESEARCHOWARMIJO;
332: tao->ops->setup = TaoSetUp_OWLQN;
333: tao->ops->solve = TaoSolve_OWLQN;
334: tao->ops->view = TaoView_OWLQN;
335: tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
336: tao->ops->destroy = TaoDestroy_OWLQN;
338: PetscNewLog(tao,&lmP);
339: lmP->D = 0;
340: lmP->M = 0;
341: lmP->GV = 0;
342: lmP->Xold = 0;
343: lmP->Gold = 0;
344: lmP->lambda = 1.0;
346: tao->data = (void*)lmP;
347: /* Override default settings (unless already changed) */
348: if (!tao->max_it_changed) tao->max_it = 2000;
349: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
351: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
352: TaoLineSearchSetType(tao->linesearch,owarmijo_type);
353: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
354: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
355: return(0);
356: }