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: PetscFunctionBegin;
15: PetscCall(VecGetOwnershipRange(d, &low, &high));
16: PetscCall(VecGetOwnershipRange(g, &low1, &high1));
18: PetscCall(VecGetArrayRead(g, &gptr));
19: PetscCall(VecGetArray(d, &dptr));
20: for (i = 0; i < high - low; i++) {
21: if (dptr[i] * gptr[i] <= 0.0) dptr[i] = 0.0;
22: }
23: PetscCall(VecRestoreArray(d, &dptr));
24: PetscCall(VecRestoreArrayRead(g, &gptr));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda)
29: {
30: const PetscReal *xptr;
31: PetscReal *gptr;
32: PetscInt low, high, low1, high1, i;
34: PetscFunctionBegin;
35: PetscCall(VecGetOwnershipRange(x, &low, &high));
36: PetscCall(VecGetOwnershipRange(gv, &low1, &high1));
38: PetscCall(VecGetArrayRead(x, &xptr));
39: PetscCall(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: PetscCall(VecRestoreArray(gv, &gptr));
48: PetscCall(VecRestoreArrayRead(x, &xptr));
49: PetscFunctionReturn(PETSC_SUCCESS);
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: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
61: PetscFunctionBegin;
62: if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by owlqn algorithm\n"));
64: /* Check convergence criteria */
65: PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
66: PetscCall(VecCopy(tao->gradient, lmP->GV));
67: PetscCall(ComputePseudoGrad_OWLQN(tao->solution, lmP->GV, lmP->lambda));
68: PetscCall(VecNorm(lmP->GV, NORM_2, &gnorm));
69: PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
71: tao->reason = TAO_CONTINUE_ITERATING;
72: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
73: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
74: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
75: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
77: /* Set initial scaling for the function */
78: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm * gnorm);
79: PetscCall(MatLMVMSetJ0Scale(lmP->M, delta));
81: /* Set counter for gradient/reset steps */
82: lmP->bfgs = 0;
83: lmP->sgrad = 0;
84: lmP->grad = 0;
86: /* Have not converged; continue with Newton method */
87: while (tao->reason == TAO_CONTINUE_ITERATING) {
88: /* Call general purpose update function */
89: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
91: /* Compute direction */
92: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
93: PetscCall(MatSolve(lmP->M, lmP->GV, lmP->D));
95: PetscCall(ProjDirect_OWLQN(lmP->D, lmP->GV));
97: ++lmP->bfgs;
99: /* Check for success (descent direction) */
100: PetscCall(VecDot(lmP->D, lmP->GV, &gdx));
101: if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
102: /* Step is not descent or direction produced not a number
103: We can assert bfgsUpdates > 1 in this case because
104: the first solve produces the scaled gradient direction,
105: which is guaranteed to be descent
107: Use steepest descent direction (scaled) */
108: ++lmP->grad;
110: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm * gnorm);
111: PetscCall(MatLMVMSetJ0Scale(lmP->M, delta));
112: PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
113: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
114: PetscCall(MatSolve(lmP->M, lmP->GV, lmP->D));
116: PetscCall(ProjDirect_OWLQN(lmP->D, lmP->GV));
118: lmP->bfgs = 1;
119: ++lmP->sgrad;
120: stepType = OWLQN_SCALED_GRADIENT;
121: } else {
122: if (1 == lmP->bfgs) {
123: /* The first BFGS direction is always the scaled gradient */
124: ++lmP->sgrad;
125: stepType = OWLQN_SCALED_GRADIENT;
126: } else {
127: ++lmP->bfgs;
128: stepType = OWLQN_BFGS;
129: }
130: }
132: PetscCall(VecScale(lmP->D, -1.0));
134: /* Perform the linesearch */
135: fold = f;
136: PetscCall(VecCopy(tao->solution, lmP->Xold));
137: PetscCall(VecCopy(tao->gradient, lmP->Gold));
139: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status));
140: PetscCall(TaoAddLineSearchCounts(tao));
142: while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) {
143: /* Reset factors and use scaled gradient step */
144: f = fold;
145: PetscCall(VecCopy(lmP->Xold, tao->solution));
146: PetscCall(VecCopy(lmP->Gold, tao->gradient));
147: PetscCall(VecCopy(tao->gradient, lmP->GV));
149: PetscCall(ComputePseudoGrad_OWLQN(tao->solution, lmP->GV, lmP->lambda));
151: switch (stepType) {
152: case OWLQN_BFGS:
153: /* Failed to obtain acceptable iterate with BFGS step
154: Attempt to use the scaled gradient direction */
156: delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm * gnorm);
157: PetscCall(MatLMVMSetJ0Scale(lmP->M, delta));
158: PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
159: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
160: PetscCall(MatSolve(lmP->M, lmP->GV, lmP->D));
162: PetscCall(ProjDirect_OWLQN(lmP->D, lmP->GV));
164: lmP->bfgs = 1;
165: ++lmP->sgrad;
166: stepType = OWLQN_SCALED_GRADIENT;
167: break;
169: case OWLQN_SCALED_GRADIENT:
170: /* The scaled gradient step did not produce a new iterate;
171: attempt to use the gradient direction.
172: Need to make sure we are not using a different diagonal scaling */
173: PetscCall(MatLMVMSetJ0Scale(lmP->M, 1.0));
174: PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
175: PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
176: PetscCall(MatSolve(lmP->M, lmP->GV, lmP->D));
178: PetscCall(ProjDirect_OWLQN(lmP->D, lmP->GV));
180: lmP->bfgs = 1;
181: ++lmP->grad;
182: stepType = OWLQN_GRADIENT;
183: break;
184: }
185: PetscCall(VecScale(lmP->D, -1.0));
187: /* Perform the linesearch */
188: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status));
189: PetscCall(TaoAddLineSearchCounts(tao));
190: }
192: if ((int)ls_status < 0) {
193: /* Failed to find an improving point*/
194: f = fold;
195: PetscCall(VecCopy(lmP->Xold, tao->solution));
196: PetscCall(VecCopy(lmP->Gold, tao->gradient));
197: PetscCall(VecCopy(tao->gradient, lmP->GV));
198: step = 0.0;
199: } else {
200: /* a little hack here, because that gv is used to store g */
201: PetscCall(VecCopy(lmP->GV, tao->gradient));
202: }
204: PetscCall(ComputePseudoGrad_OWLQN(tao->solution, lmP->GV, lmP->lambda));
206: /* Check for termination */
208: PetscCall(VecNorm(lmP->GV, NORM_2, &gnorm));
210: ++tao->niter;
211: PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
212: PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
213: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
215: if ((int)ls_status < 0) break;
216: }
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
221: {
222: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
223: PetscInt n, N;
225: PetscFunctionBegin;
226: /* Existence of tao->solution checked in TaoSetUp() */
227: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
228: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
229: if (!lmP->D) PetscCall(VecDuplicate(tao->solution, &lmP->D));
230: if (!lmP->GV) PetscCall(VecDuplicate(tao->solution, &lmP->GV));
231: if (!lmP->Xold) PetscCall(VecDuplicate(tao->solution, &lmP->Xold));
232: if (!lmP->Gold) PetscCall(VecDuplicate(tao->solution, &lmP->Gold));
234: /* Create matrix for the limited memory approximation */
235: PetscCall(VecGetLocalSize(tao->solution, &n));
236: PetscCall(VecGetSize(tao->solution, &N));
237: PetscCall(MatCreateLMVMBFGS(((PetscObject)tao)->comm, n, N, &lmP->M));
238: PetscCall(MatLMVMAllocate(lmP->M, tao->solution, tao->gradient));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
243: {
244: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
246: PetscFunctionBegin;
247: if (tao->setupcalled) {
248: PetscCall(VecDestroy(&lmP->Xold));
249: PetscCall(VecDestroy(&lmP->Gold));
250: PetscCall(VecDestroy(&lmP->D));
251: PetscCall(MatDestroy(&lmP->M));
252: PetscCall(VecDestroy(&lmP->GV));
253: }
254: PetscCall(PetscFree(tao->data));
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: static PetscErrorCode TaoSetFromOptions_OWLQN(Tao tao, PetscOptionItems PetscOptionsObject)
259: {
260: TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
262: PetscFunctionBegin;
263: PetscOptionsHeadBegin(PetscOptionsObject, "Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
264: PetscCall(PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight", "", 100, &lmP->lambda, NULL));
265: PetscOptionsHeadEnd();
266: PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
271: {
272: TAO_OWLQN *lm = (TAO_OWLQN *)tao->data;
273: PetscBool isascii;
275: PetscFunctionBegin;
276: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
277: if (isascii) {
278: PetscCall(PetscViewerASCIIPushTab(viewer));
279: PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", lm->bfgs));
280: PetscCall(PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %" PetscInt_FMT "\n", lm->sgrad));
281: PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lm->grad));
282: PetscCall(PetscViewerASCIIPopTab(viewer));
283: }
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*MC
288: TAOOWLQN - orthant-wise limited memory quasi-newton algorithm
290: . - tao_owlqn_lambda - regulariser weight
292: Level: beginner
293: M*/
295: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
296: {
297: TAO_OWLQN *lmP;
298: const char *owarmijo_type = TAOLINESEARCHOWARMIJO;
300: PetscFunctionBegin;
301: tao->ops->setup = TaoSetUp_OWLQN;
302: tao->ops->solve = TaoSolve_OWLQN;
303: tao->ops->view = TaoView_OWLQN;
304: tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
305: tao->ops->destroy = TaoDestroy_OWLQN;
307: PetscCall(PetscNew(&lmP));
308: lmP->D = NULL;
309: lmP->M = NULL;
310: lmP->GV = NULL;
311: lmP->Xold = NULL;
312: lmP->Gold = NULL;
313: lmP->lambda = 1.0;
315: tao->data = (void *)lmP;
316: /* Override default settings (unless already changed) */
317: PetscCall(TaoParametersInitialize(tao));
318: PetscObjectParameterSetDefault(tao, max_it, 2000);
319: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
321: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
322: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
323: PetscCall(TaoLineSearchSetType(tao->linesearch, owarmijo_type));
324: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
325: PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));
326: PetscFunctionReturn(PETSC_SUCCESS);
327: }