1 #include <petsctaolinesearch.h>
2 #include <../src/tao/unconstrained/impls/owlqn/owlqn.h>
3
4 #define OWLQN_BFGS 0
5 #define OWLQN_SCALED_GRADIENT 1
6 #define OWLQN_GRADIENT 2
7
ProjDirect_OWLQN(Vec d,Vec g)8 static PetscErrorCode ProjDirect_OWLQN(Vec d, Vec g)
9 {
10 const PetscReal *gptr;
11 PetscReal *dptr;
12 PetscInt low, high, low1, high1, i;
13
14 PetscFunctionBegin;
15 PetscCall(VecGetOwnershipRange(d, &low, &high));
16 PetscCall(VecGetOwnershipRange(g, &low1, &high1));
17
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 }
27
ComputePseudoGrad_OWLQN(Vec x,Vec gv,PetscReal lambda)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;
33
34 PetscFunctionBegin;
35 PetscCall(VecGetOwnershipRange(x, &low, &high));
36 PetscCall(VecGetOwnershipRange(gv, &low1, &high1));
37
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 }
51
TaoSolve_OWLQN(Tao tao)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;
60
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"));
63
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");
70
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);
76
77 /* Set initial scaling for the function */
78 delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm * gnorm);
79 PetscCall(MatLMVMSetJ0Scale(lmP->M, delta));
80
81 /* Set counter for gradient/reset steps */
82 lmP->bfgs = 0;
83 lmP->sgrad = 0;
84 lmP->grad = 0;
85
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);
90
91 /* Compute direction */
92 PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
93 PetscCall(MatSolve(lmP->M, lmP->GV, lmP->D));
94
95 PetscCall(ProjDirect_OWLQN(lmP->D, lmP->GV));
96
97 ++lmP->bfgs;
98
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
106
107 Use steepest descent direction (scaled) */
108 ++lmP->grad;
109
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));
115
116 PetscCall(ProjDirect_OWLQN(lmP->D, lmP->GV));
117
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 }
131
132 PetscCall(VecScale(lmP->D, -1.0));
133
134 /* Perform the linesearch */
135 fold = f;
136 PetscCall(VecCopy(tao->solution, lmP->Xold));
137 PetscCall(VecCopy(tao->gradient, lmP->Gold));
138
139 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status));
140 PetscCall(TaoAddLineSearchCounts(tao));
141
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));
148
149 PetscCall(ComputePseudoGrad_OWLQN(tao->solution, lmP->GV, lmP->lambda));
150
151 switch (stepType) {
152 case OWLQN_BFGS:
153 /* Failed to obtain acceptable iterate with BFGS step
154 Attempt to use the scaled gradient direction */
155
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));
161
162 PetscCall(ProjDirect_OWLQN(lmP->D, lmP->GV));
163
164 lmP->bfgs = 1;
165 ++lmP->sgrad;
166 stepType = OWLQN_SCALED_GRADIENT;
167 break;
168
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));
177
178 PetscCall(ProjDirect_OWLQN(lmP->D, lmP->GV));
179
180 lmP->bfgs = 1;
181 ++lmP->grad;
182 stepType = OWLQN_GRADIENT;
183 break;
184 }
185 PetscCall(VecScale(lmP->D, -1.0));
186
187 /* Perform the linesearch */
188 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status));
189 PetscCall(TaoAddLineSearchCounts(tao));
190 }
191
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 }
203
204 PetscCall(ComputePseudoGrad_OWLQN(tao->solution, lmP->GV, lmP->lambda));
205
206 /* Check for termination */
207
208 PetscCall(VecNorm(lmP->GV, NORM_2, &gnorm));
209
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);
214
215 if ((int)ls_status < 0) break;
216 }
217 PetscFunctionReturn(PETSC_SUCCESS);
218 }
219
TaoSetUp_OWLQN(Tao tao)220 static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
221 {
222 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
223 PetscInt n, N;
224
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));
233
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 }
241
TaoDestroy_OWLQN(Tao tao)242 static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
243 {
244 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
245
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 }
257
TaoSetFromOptions_OWLQN(Tao tao,PetscOptionItems PetscOptionsObject)258 static PetscErrorCode TaoSetFromOptions_OWLQN(Tao tao, PetscOptionItems PetscOptionsObject)
259 {
260 TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
261
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 }
269
TaoView_OWLQN(Tao tao,PetscViewer viewer)270 static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
271 {
272 TAO_OWLQN *lm = (TAO_OWLQN *)tao->data;
273 PetscBool isascii;
274
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 }
286
287 /*MC
288 TAOOWLQN - orthant-wise limited memory quasi-newton algorithm
289
290 . - tao_owlqn_lambda - regulariser weight
291
292 Level: beginner
293 M*/
294
TaoCreate_OWLQN(Tao tao)295 PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
296 {
297 TAO_OWLQN *lmP;
298 const char *owarmijo_type = TAOLINESEARCHOWARMIJO;
299
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;
306
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;
314
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);
320
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 }
328