Actual source code: owlqn.c
petsc-3.5.4 2015-05-23
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
11: static PetscErrorCode ProjDirect_OWLQN(Vec d, Vec g)
12: {
13: PetscErrorCode ierr;
14: const PetscReal *gptr;
15: PetscReal *dptr;
16: PetscInt low,high,low1,high1,i;
19: ierr=VecGetOwnershipRange(d,&low,&high);
20: ierr=VecGetOwnershipRange(g,&low1,&high1);
22: VecGetArrayRead(g,&gptr);
23: VecGetArray(d,&dptr);
24: for (i = 0; i < high-low; i++) {
25: if (dptr[i] * gptr[i] <= 0.0 ) {
26: dptr[i] = 0.0;
27: }
28: }
29: VecRestoreArray(d,&dptr);
30: VecRestoreArrayRead(g,&gptr);
31: return(0);
32: }
36: static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda)
37: {
38: PetscErrorCode ierr;
39: const PetscReal *xptr;
40: PetscReal *gptr;
41: PetscInt low,high,low1,high1,i;
44: ierr=VecGetOwnershipRange(x,&low,&high);
45: ierr=VecGetOwnershipRange(gv,&low1,&high1);
47: VecGetArrayRead(x,&xptr);
48: VecGetArray(gv,&gptr);
49: for (i = 0; i < high-low; i++) {
50: if (xptr[i] < 0.0) gptr[i] = gptr[i] - lambda;
51: else if (xptr[i] > 0.0) gptr[i] = gptr[i] + lambda;
52: else if (gptr[i] + lambda < 0.0) gptr[i] = gptr[i] + lambda;
53: else if (gptr[i] - lambda > 0.0) gptr[i] = gptr[i] - lambda;
54: else gptr[i] = 0.0;
55: }
56: VecRestoreArray(gv,&gptr);
57: VecRestoreArrayRead(x,&xptr);
58: return(0);
59: }
63: static PetscErrorCode TaoSolve_OWLQN(Tao tao)
64: {
65: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
66: PetscReal f, fold, gdx, gnorm;
67: PetscReal step = 1.0;
68: PetscReal delta;
69: PetscErrorCode ierr;
70: PetscInt stepType;
71: PetscInt iter = 0;
72: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
73: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
76: if (tao->XL || tao->XU || tao->ops->computebounds) {
77: PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by owlqn algorithm\n");
78: }
80: /* Check convergence criteria */
81: TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
83: VecCopy(tao->gradient, lmP->GV);
85: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
87: VecNorm(lmP->GV,NORM_2,&gnorm);
89: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");
91: TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);
92: if (reason != TAO_CONTINUE_ITERATING) return(0);
94: /* Set initial scaling for the function */
95: if (f != 0.0) {
96: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
97: } else {
98: delta = 2.0 / (gnorm*gnorm);
99: }
100: MatLMVMSetDelta(lmP->M,delta);
102: /* Set counter for gradient/reset steps */
103: lmP->bfgs = 0;
104: lmP->sgrad = 0;
105: lmP->grad = 0;
107: /* Have not converged; continue with Newton method */
108: while (reason == TAO_CONTINUE_ITERATING) {
109: /* Compute direction */
110: MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
111: MatLMVMSolve(lmP->M, lmP->GV, lmP->D);
113: ProjDirect_OWLQN(lmP->D,lmP->GV);
115: ++lmP->bfgs;
117: /* Check for success (descent direction) */
118: VecDot(lmP->D, lmP->GV , &gdx);
119: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
121: /* Step is not descent or direction produced not a number
122: We can assert bfgsUpdates > 1 in this case because
123: the first solve produces the scaled gradient direction,
124: which is guaranteed to be descent
126: Use steepest descent direction (scaled) */
127: ++lmP->grad;
129: if (f != 0.0) {
130: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
131: } else {
132: delta = 2.0 / (gnorm*gnorm);
133: }
134: MatLMVMSetDelta(lmP->M, delta);
135: MatLMVMReset(lmP->M);
136: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
137: MatLMVMSolve(lmP->M,lmP->GV, lmP->D);
139: ProjDirect_OWLQN(lmP->D,lmP->GV);
141: lmP->bfgs = 1;
142: ++lmP->sgrad;
143: stepType = OWLQN_SCALED_GRADIENT;
144: } else {
145: if (1 == lmP->bfgs) {
146: /* The first BFGS direction is always the scaled gradient */
147: ++lmP->sgrad;
148: stepType = OWLQN_SCALED_GRADIENT;
149: } else {
150: ++lmP->bfgs;
151: stepType = OWLQN_BFGS;
152: }
153: }
155: VecScale(lmP->D, -1.0);
157: /* Perform the linesearch */
158: fold = f;
159: VecCopy(tao->solution, lmP->Xold);
160: VecCopy(tao->gradient, lmP->Gold);
162: TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step,&ls_status);
163: TaoAddLineSearchCounts(tao);
165: while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) {
167: /* Reset factors and use scaled gradient step */
168: f = fold;
169: VecCopy(lmP->Xold, tao->solution);
170: VecCopy(lmP->Gold, tao->gradient);
171: VecCopy(tao->gradient, lmP->GV);
173: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
175: switch(stepType) {
176: case OWLQN_BFGS:
177: /* Failed to obtain acceptable iterate with BFGS step
178: Attempt to use the scaled gradient direction */
180: if (f != 0.0) {
181: delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
182: } else {
183: delta = 2.0 / (gnorm*gnorm);
184: }
185: MatLMVMSetDelta(lmP->M, delta);
186: MatLMVMReset(lmP->M);
187: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
188: MatLMVMSolve(lmP->M, lmP->GV, lmP->D);
190: ProjDirect_OWLQN(lmP->D,lmP->GV);
192: lmP->bfgs = 1;
193: ++lmP->sgrad;
194: stepType = OWLQN_SCALED_GRADIENT;
195: break;
197: case OWLQN_SCALED_GRADIENT:
198: /* The scaled gradient step did not produce a new iterate;
199: attempt to use the gradient direction.
200: Need to make sure we are not using a different diagonal scaling */
201: MatLMVMSetDelta(lmP->M, 1.0);
202: MatLMVMReset(lmP->M);
203: MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
204: MatLMVMSolve(lmP->M, lmP->GV, lmP->D);
206: ProjDirect_OWLQN(lmP->D,lmP->GV);
208: lmP->bfgs = 1;
209: ++lmP->grad;
210: stepType = OWLQN_GRADIENT;
211: break;
212: }
213: VecScale(lmP->D, -1.0);
216: /* Perform the linesearch */
217: TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
218: TaoAddLineSearchCounts(tao);
219: }
221: if ((int)ls_status < 0) {
222: /* Failed to find an improving point*/
223: f = fold;
224: VecCopy(lmP->Xold, tao->solution);
225: VecCopy(lmP->Gold, tao->gradient);
226: VecCopy(tao->gradient, lmP->GV);
227: step = 0.0;
228: } else {
229: /* a little hack here, because that gv is used to store g */
230: VecCopy(lmP->GV, tao->gradient);
231: }
233: ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
235: /* Check for termination */
237: VecNorm(lmP->GV,NORM_2,&gnorm);
239: iter++;
240: TaoMonitor(tao,iter,f,gnorm,0.0,step,&reason);
242: if ((int)ls_status < 0) break;
243: }
244: return(0);
245: }
249: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
250: {
251: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
252: PetscInt n,N;
256: /* Existence of tao->solution checked in TaoSetUp() */
257: if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient); }
258: if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection); }
259: if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D); }
260: if (!lmP->GV) {VecDuplicate(tao->solution,&lmP->GV); }
261: if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold); }
262: if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold); }
264: /* Create matrix for the limited memory approximation */
265: VecGetLocalSize(tao->solution,&n);
266: VecGetSize(tao->solution,&N);
267: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M);
268: MatLMVMAllocateVectors(lmP->M,tao->solution);
269: return(0);
270: }
272: /* ---------------------------------------------------------- */
275: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
276: {
277: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
281: if (tao->setupcalled) {
282: VecDestroy(&lmP->Xold);
283: VecDestroy(&lmP->Gold);
284: VecDestroy(&lmP->D);
285: MatDestroy(&lmP->M);
286: VecDestroy(&lmP->GV);
287: }
288: PetscFree(tao->data);
289: return(0);
290: }
292: /*------------------------------------------------------------*/
295: static PetscErrorCode TaoSetFromOptions_OWLQN(Tao tao)
296: {
297: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
298: PetscBool flg;
302: PetscOptionsHead("Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
303: PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,&flg);
304: PetscOptionsTail();
305: TaoLineSearchSetFromOptions(tao->linesearch);
306: return(0);
307: }
309: /*------------------------------------------------------------*/
312: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
313: {
314: TAO_OWLQN *lm = (TAO_OWLQN *)tao->data;
315: PetscBool isascii;
319: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
320: if (isascii) {
321: PetscViewerASCIIPushTab(viewer);
322: PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
323: PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
324: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
325: PetscViewerASCIIPopTab(viewer);
326: }
327: return(0);
328: }
330: /* ---------------------------------------------------------- */
331: /*MC
332: TAOOWLQN - orthant-wise limited memory quasi-newton algorithm
334: . - tao_owlqn_lambda - regulariser weight
336: Level: beginner
337: M*/
340: EXTERN_C_BEGIN
343: PetscErrorCode TaoCreate_OWLQN(Tao tao)
344: {
345: TAO_OWLQN *lmP;
346: const char *owarmijo_type = TAOLINESEARCHOWARMIJO;
350: tao->ops->setup = TaoSetUp_OWLQN;
351: tao->ops->solve = TaoSolve_OWLQN;
352: tao->ops->view = TaoView_OWLQN;
353: tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
354: tao->ops->destroy = TaoDestroy_OWLQN;
356: PetscNewLog(tao,&lmP);
357: lmP->D = 0;
358: lmP->M = 0;
359: lmP->GV = 0;
360: lmP->Xold = 0;
361: lmP->Gold = 0;
362: lmP->lambda = 1.0;
364: tao->data = (void*)lmP;
365: tao->max_it = 2000;
366: tao->max_funcs = 4000;
367: tao->fatol = 1e-4;
368: tao->frtol = 1e-4;
370: TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
371: TaoLineSearchSetType(tao->linesearch,owarmijo_type);
372: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
373: return(0);
374: }
375: EXTERN_C_END