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