Actual source code: owlqn.c
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: const PetscReal *gptr;
11: PetscReal *dptr;
12: PetscInt low,high,low1,high1,i;
14: VecGetOwnershipRange(d,&low,&high);
15: VecGetOwnershipRange(g,&low1,&high1);
17: VecGetArrayRead(g,&gptr);
18: VecGetArray(d,&dptr);
19: for (i = 0; i < high-low; i++) {
20: if (dptr[i] * gptr[i] <= 0.0) {
21: dptr[i] = 0.0;
22: }
23: }
24: VecRestoreArray(d,&dptr);
25: VecRestoreArrayRead(g,&gptr);
26: return 0;
27: }
29: static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda)
30: {
31: const PetscReal *xptr;
32: PetscReal *gptr;
33: PetscInt low,high,low1,high1,i;
35: VecGetOwnershipRange(x,&low,&high);
36: VecGetOwnershipRange(gv,&low1,&high1);
38: VecGetArrayRead(x,&xptr);
39: VecGetArray(gv,&gptr);
40: for (i = 0; i < high-low; i++) {
41: if (xptr[i] < 0.0) gptr[i] = gptr[i] - lambda;
42: else if (xptr[i] > 0.0) gptr[i] = gptr[i] + lambda;
43: else if (gptr[i] + lambda < 0.0) gptr[i] = gptr[i] + lambda;
44: else if (gptr[i] - lambda > 0.0) gptr[i] = gptr[i] - lambda;
45: else gptr[i] = 0.0;
46: }
47: VecRestoreArray(gv,&gptr);
48: VecRestoreArrayRead(x,&xptr);
49: return 0;
50: }
52: static PetscErrorCode TaoSolve_OWLQN(Tao tao)
53: {
54: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
55: PetscReal f, fold, gdx, gnorm;
56: PetscReal step = 1.0;
57: PetscReal delta;
58: PetscInt stepType;
59: PetscInt iter = 0;
60: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
62: if (tao->XL || tao->XU || tao->ops->computebounds) {
63: PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by owlqn algorithm\n");
64: }
66: /* Check convergence criteria */
67: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
68: VecCopy(tao->gradient, lmP->GV);
69: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
70: VecNorm(lmP->GV,NORM_2,&gnorm);
73: tao->reason = TAO_CONTINUE_ITERATING;
74: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
75: TaoMonitor(tao,iter,f,gnorm,0.0,step);
76: (*tao->ops->convergencetest)(tao,tao->cnvP);
77: if (tao->reason != TAO_CONTINUE_ITERATING) return 0;
79: /* Set initial scaling for the function */
80: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
81: MatLMVMSetJ0Scale(lmP->M, delta);
83: /* Set counter for gradient/reset steps */
84: lmP->bfgs = 0;
85: lmP->sgrad = 0;
86: lmP->grad = 0;
88: /* Have not converged; continue with Newton method */
89: while (tao->reason == TAO_CONTINUE_ITERATING) {
90: /* Call general purpose update function */
91: if (tao->ops->update) {
92: (*tao->ops->update)(tao, tao->niter, tao->user_update);
93: }
95: /* Compute direction */
96: MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
97: MatSolve(lmP->M, lmP->GV, lmP->D);
99: ProjDirect_OWLQN(lmP->D,lmP->GV);
101: ++lmP->bfgs;
103: /* Check for success (descent direction) */
104: VecDot(lmP->D, lmP->GV , &gdx);
105: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
107: /* Step is not descent or direction produced not a number
108: We can assert bfgsUpdates > 1 in this case because
109: the first solve produces the scaled gradient direction,
110: which is guaranteed to be descent
112: Use steepest descent direction (scaled) */
113: ++lmP->grad;
115: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
116: MatLMVMSetJ0Scale(lmP->M, delta);
117: MatLMVMReset(lmP->M, PETSC_FALSE);
118: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
119: MatSolve(lmP->M,lmP->GV, lmP->D);
121: ProjDirect_OWLQN(lmP->D,lmP->GV);
123: lmP->bfgs = 1;
124: ++lmP->sgrad;
125: stepType = OWLQN_SCALED_GRADIENT;
126: } else {
127: if (1 == lmP->bfgs) {
128: /* The first BFGS direction is always the scaled gradient */
129: ++lmP->sgrad;
130: stepType = OWLQN_SCALED_GRADIENT;
131: } else {
132: ++lmP->bfgs;
133: stepType = OWLQN_BFGS;
134: }
135: }
137: VecScale(lmP->D, -1.0);
139: /* Perform the linesearch */
140: fold = f;
141: VecCopy(tao->solution, lmP->Xold);
142: VecCopy(tao->gradient, lmP->Gold);
144: TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step,&ls_status);
145: TaoAddLineSearchCounts(tao);
147: while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) {
149: /* Reset factors and use scaled gradient step */
150: f = fold;
151: VecCopy(lmP->Xold, tao->solution);
152: VecCopy(lmP->Gold, tao->gradient);
153: VecCopy(tao->gradient, lmP->GV);
155: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
157: switch(stepType) {
158: case OWLQN_BFGS:
159: /* Failed to obtain acceptable iterate with BFGS step
160: Attempt to use the scaled gradient direction */
162: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
163: MatLMVMSetJ0Scale(lmP->M, delta);
164: MatLMVMReset(lmP->M, PETSC_FALSE);
165: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
166: MatSolve(lmP->M, lmP->GV, lmP->D);
168: ProjDirect_OWLQN(lmP->D,lmP->GV);
170: lmP->bfgs = 1;
171: ++lmP->sgrad;
172: stepType = OWLQN_SCALED_GRADIENT;
173: break;
175: case OWLQN_SCALED_GRADIENT:
176: /* The scaled gradient step did not produce a new iterate;
177: attempt to use the gradient direction.
178: Need to make sure we are not using a different diagonal scaling */
179: MatLMVMSetJ0Scale(lmP->M, 1.0);
180: MatLMVMReset(lmP->M, PETSC_FALSE);
181: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
182: MatSolve(lmP->M, lmP->GV, lmP->D);
184: ProjDirect_OWLQN(lmP->D,lmP->GV);
186: lmP->bfgs = 1;
187: ++lmP->grad;
188: stepType = OWLQN_GRADIENT;
189: break;
190: }
191: VecScale(lmP->D, -1.0);
193: /* Perform the linesearch */
194: TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
195: TaoAddLineSearchCounts(tao);
196: }
198: if ((int)ls_status < 0) {
199: /* Failed to find an improving point*/
200: f = fold;
201: VecCopy(lmP->Xold, tao->solution);
202: VecCopy(lmP->Gold, tao->gradient);
203: VecCopy(tao->gradient, lmP->GV);
204: step = 0.0;
205: } else {
206: /* a little hack here, because that gv is used to store g */
207: VecCopy(lmP->GV, tao->gradient);
208: }
210: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
212: /* Check for termination */
214: VecNorm(lmP->GV,NORM_2,&gnorm);
216: iter++;
217: TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
218: TaoMonitor(tao,iter,f,gnorm,0.0,step);
219: (*tao->ops->convergencetest)(tao,tao->cnvP);
221: if ((int)ls_status < 0) break;
222: }
223: return 0;
224: }
226: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
227: {
228: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
229: PetscInt n,N;
231: /* Existence of tao->solution checked in TaoSetUp() */
232: if (!tao->gradient) VecDuplicate(tao->solution,&tao->gradient);
233: if (!tao->stepdirection) VecDuplicate(tao->solution,&tao->stepdirection);
234: if (!lmP->D) VecDuplicate(tao->solution,&lmP->D);
235: if (!lmP->GV) VecDuplicate(tao->solution,&lmP->GV);
236: if (!lmP->Xold) VecDuplicate(tao->solution,&lmP->Xold);
237: if (!lmP->Gold) VecDuplicate(tao->solution,&lmP->Gold);
239: /* Create matrix for the limited memory approximation */
240: VecGetLocalSize(tao->solution,&n);
241: VecGetSize(tao->solution,&N);
242: MatCreateLMVMBFGS(((PetscObject)tao)->comm,n,N,&lmP->M);
243: MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);
244: return 0;
245: }
247: /* ---------------------------------------------------------- */
248: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
249: {
250: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
252: if (tao->setupcalled) {
253: VecDestroy(&lmP->Xold);
254: VecDestroy(&lmP->Gold);
255: VecDestroy(&lmP->D);
256: MatDestroy(&lmP->M);
257: VecDestroy(&lmP->GV);
258: }
259: PetscFree(tao->data);
260: return 0;
261: }
263: /*------------------------------------------------------------*/
264: static PetscErrorCode TaoSetFromOptions_OWLQN(PetscOptionItems *PetscOptionsObject,Tao tao)
265: {
266: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
268: PetscOptionsHead(PetscOptionsObject,"Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
269: PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,NULL);
270: PetscOptionsTail();
271: TaoLineSearchSetFromOptions(tao->linesearch);
272: return 0;
273: }
275: /*------------------------------------------------------------*/
276: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
277: {
278: TAO_OWLQN *lm = (TAO_OWLQN *)tao->data;
279: PetscBool isascii;
281: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
282: if (isascii) {
283: PetscViewerASCIIPushTab(viewer);
284: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
285: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
286: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
287: PetscViewerASCIIPopTab(viewer);
288: }
289: return 0;
290: }
292: /* ---------------------------------------------------------- */
293: /*MC
294: TAOOWLQN - orthant-wise limited memory quasi-newton algorithm
296: . - tao_owlqn_lambda - regulariser weight
298: Level: beginner
299: M*/
301: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
302: {
303: TAO_OWLQN *lmP;
304: const char *owarmijo_type = TAOLINESEARCHOWARMIJO;
306: tao->ops->setup = TaoSetUp_OWLQN;
307: tao->ops->solve = TaoSolve_OWLQN;
308: tao->ops->view = TaoView_OWLQN;
309: tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
310: tao->ops->destroy = TaoDestroy_OWLQN;
312: PetscNewLog(tao,&lmP);
313: lmP->D = NULL;
314: lmP->M = NULL;
315: lmP->GV = NULL;
316: lmP->Xold = NULL;
317: lmP->Gold = NULL;
318: lmP->lambda = 1.0;
320: tao->data = (void*)lmP;
321: /* Override default settings (unless already changed) */
322: if (!tao->max_it_changed) tao->max_it = 2000;
323: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
325: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
326: PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
327: TaoLineSearchSetType(tao->linesearch,owarmijo_type);
328: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
329: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
330: return 0;
331: }