1 #include <../src/tao/unconstrained/impls/bmrm/bmrm.h>
2
3 static PetscErrorCode init_df_solver(TAO_DF *);
4 static PetscErrorCode ensure_df_space(PetscInt, TAO_DF *);
5 static PetscErrorCode destroy_df_solver(TAO_DF *);
6 static PetscReal phi(PetscReal *, PetscInt, PetscReal, PetscReal *, PetscReal, PetscReal *, PetscReal *, PetscReal *);
7 static PetscInt project(PetscInt, PetscReal *, PetscReal, PetscReal *, PetscReal *, PetscReal *, PetscReal *, PetscReal *, TAO_DF *);
8
solve(TAO_DF * df)9 static PetscErrorCode solve(TAO_DF *df)
10 {
11 PetscInt i, j, innerIter, it, it2, luv, info;
12 PetscReal gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam = 0.0, lam_ext;
13 PetscReal DELTAsv, ProdDELTAsv;
14 PetscReal c, *tempQ;
15 PetscReal *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol;
16 PetscReal *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd;
17 PetscReal *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk;
18 PetscReal **Q = df->Q, *f = df->f, *t = df->t;
19 PetscInt dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv;
20
21 PetscCheck(dim >= 0, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Expected dim %" PetscInt_FMT " >= 0", dim);
22 /* variables for the adaptive nonmonotone linesearch */
23 PetscInt L, llast;
24 PetscReal fr, fbest, fv, fc, fv0;
25
26 c = BMRM_INFTY;
27
28 DELTAsv = EPS_SV;
29 if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F;
30 else ProdDELTAsv = EPS_SV;
31
32 for (i = 0; i < dim; i++) tempv[i] = -x[i];
33
34 lam_ext = 0.0;
35
36 /* Project the initial solution */
37 project(dim, a, b, tempv, l, u, x, &lam_ext, df);
38
39 /* Compute gradient
40 g = Q*x + f; */
41
42 it = 0;
43 for (i = 0; i < dim; i++) {
44 if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i;
45 }
46
47 PetscCall(PetscArrayzero(t, dim));
48 for (i = 0; i < it; i++) {
49 tempQ = Q[ipt[i]];
50 for (j = 0; j < dim; j++) t[j] += (tempQ[j] * x[ipt[i]]);
51 }
52 for (i = 0; i < dim; i++) g[i] = t[i] + f[i];
53
54 /* y = -(x_{k} - g_{k}) */
55 for (i = 0; i < dim; i++) y[i] = g[i] - x[i];
56
57 /* Project x_{k} - g_{k} */
58 project(dim, a, b, y, l, u, tempv, &lam_ext, df);
59
60 /* y = P(x_{k} - g_{k}) - x_{k} */
61 max = ALPHA_MIN;
62 for (i = 0; i < dim; i++) {
63 y[i] = tempv[i] - x[i];
64 if (PetscAbsReal(y[i]) > max) max = PetscAbsReal(y[i]);
65 }
66
67 if (max < tol * 1e-3) return PETSC_SUCCESS;
68
69 alpha = 1.0 / max;
70
71 /* fv0 = f(x_{0}). Recall t = Q x_{k} */
72 fv0 = 0.0;
73 for (i = 0; i < dim; i++) fv0 += x[i] * (0.5 * t[i] + f[i]);
74
75 /* adaptive nonmonotone linesearch */
76 L = 2;
77 fr = ALPHA_MAX;
78 fbest = fv0;
79 fc = fv0;
80 llast = 0;
81 akold = bkold = 0.0;
82
83 /* Iterator begins */
84 for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) {
85 /* tempv = -(x_{k} - alpha*g_{k}) */
86 for (i = 0; i < dim; i++) tempv[i] = alpha * g[i] - x[i];
87
88 /* Project x_{k} - alpha*g_{k} */
89 project(dim, a, b, tempv, l, u, y, &lam_ext, df);
90
91 /* gd = \inner{d_{k}}{g_{k}}
92 d = P(x_{k} - alpha*g_{k}) - x_{k}
93 */
94 gd = 0.0;
95 for (i = 0; i < dim; i++) {
96 d[i] = y[i] - x[i];
97 gd += d[i] * g[i];
98 }
99
100 /* Gradient computation */
101
102 /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */
103
104 it = it2 = 0;
105 for (i = 0; i < dim; i++) {
106 if (PetscAbsReal(d[i]) > (ProdDELTAsv * 1.0e-2)) ipt[it++] = i;
107 }
108 for (i = 0; i < dim; i++) {
109 if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i;
110 }
111
112 PetscCall(PetscArrayzero(Qd, dim));
113 /* compute Qd = Q*d */
114 if (it < it2) {
115 for (i = 0; i < it; i++) {
116 tempQ = Q[ipt[i]];
117 for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]);
118 }
119 } else { /* compute Qd = Q*y-t */
120 for (i = 0; i < it2; i++) {
121 tempQ = Q[ipt2[i]];
122 for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]);
123 }
124 for (j = 0; j < dim; j++) Qd[j] -= t[j];
125 }
126
127 /* ak = inner{d_{k}}{d_{k}} */
128 ak = 0.0;
129 for (i = 0; i < dim; i++) ak += d[i] * d[i];
130
131 bk = 0.0;
132 for (i = 0; i < dim; i++) bk += d[i] * Qd[i];
133
134 if (bk > EPS * ak && gd < 0.0) lamnew = -gd / bk;
135 else lamnew = 1.0;
136
137 /* fv is computing f(x_{k} + d_{k}) */
138 fv = 0.0;
139 for (i = 0; i < dim; i++) {
140 xplus[i] = x[i] + d[i];
141 tplus[i] = t[i] + Qd[i];
142 fv += xplus[i] * (0.5 * tplus[i] + f[i]);
143 }
144
145 /* fr is fref */
146 if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)) {
147 fv = 0.0;
148 for (i = 0; i < dim; i++) {
149 xplus[i] = x[i] + lamnew * d[i];
150 tplus[i] = t[i] + lamnew * Qd[i];
151 fv += xplus[i] * (0.5 * tplus[i] + f[i]);
152 }
153 }
154
155 for (i = 0; i < dim; i++) {
156 sk[i] = xplus[i] - x[i];
157 yk[i] = tplus[i] - t[i];
158 x[i] = xplus[i];
159 t[i] = tplus[i];
160 g[i] = t[i] + f[i];
161 }
162
163 /* update the line search control parameters */
164 if (fv < fbest) {
165 fbest = fv;
166 fc = fv;
167 llast = 0;
168 } else {
169 fc = (fc > fv ? fc : fv);
170 llast++;
171 if (llast == L) {
172 fr = fc;
173 fc = fv;
174 llast = 0;
175 }
176 }
177
178 ak = bk = 0.0;
179 for (i = 0; i < dim; i++) {
180 ak += sk[i] * sk[i];
181 bk += sk[i] * yk[i];
182 }
183
184 if (bk <= EPS * ak) alpha = ALPHA_MAX;
185 else {
186 if (bkold < EPS * akold) alpha = ak / bk;
187 else alpha = (akold + ak) / (bkold + bk);
188
189 if (alpha > ALPHA_MAX) alpha = ALPHA_MAX;
190 else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN;
191 }
192
193 akold = ak;
194 bkold = bk;
195
196 /* stopping criterion based on KKT conditions */
197 /* at optimal, gradient of lagrangian w.r.t. x is zero */
198
199 bk = 0.0;
200 for (i = 0; i < dim; i++) bk += x[i] * x[i];
201
202 if (PetscSqrtReal(ak) < tol * 10 * PetscSqrtReal(bk)) {
203 it = 0;
204 luv = 0;
205 kktlam = 0.0;
206 for (i = 0; i < dim; i++) {
207 /* x[i] is active hence lagrange multipliers for box constraints
208 are zero. The lagrange multiplier for ineq. const. is then
209 defined as below
210 */
211 if ((x[i] > DELTAsv) && (x[i] < c - DELTAsv)) {
212 ipt[it++] = i;
213 kktlam = kktlam - a[i] * g[i];
214 } else uv[luv++] = i;
215 }
216
217 if (it == 0 && PetscSqrtReal(ak) < tol * 0.5 * PetscSqrtReal(bk)) return PETSC_SUCCESS;
218 else {
219 kktlam = kktlam / it;
220 info = 1;
221 for (i = 0; i < it; i++) {
222 if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) {
223 info = 0;
224 break;
225 }
226 }
227 if (info == 1) {
228 for (i = 0; i < luv; i++) {
229 if (x[uv[i]] <= DELTAsv) {
230 /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may
231 not be zero. So, the gradient without beta is > 0
232 */
233 if (g[uv[i]] + kktlam * a[uv[i]] < -tol) {
234 info = 0;
235 break;
236 }
237 } else {
238 /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may
239 not be zero. So, the gradient without eta is < 0
240 */
241 if (g[uv[i]] + kktlam * a[uv[i]] > tol) {
242 info = 0;
243 break;
244 }
245 }
246 }
247 }
248
249 if (info == 1) return PETSC_SUCCESS;
250 }
251 }
252 }
253 return PETSC_SUCCESS;
254 }
255
256 /* The main solver function
257
258 f = Remp(W) This is what the user provides us from the application layer
259 So the ComputeGradient function for instance should get us back the subgradient of Remp(W)
260
261 Regularizer assumed to be L2 norm = lambda*0.5*W'W ()
262 */
263
make_grad_node(Vec X,Vec_Chain ** p)264 static PetscErrorCode make_grad_node(Vec X, Vec_Chain **p)
265 {
266 PetscFunctionBegin;
267 PetscCall(PetscNew(p));
268 PetscCall(VecDuplicate(X, &(*p)->V));
269 PetscCall(VecCopy(X, (*p)->V));
270 (*p)->next = NULL;
271 PetscFunctionReturn(PETSC_SUCCESS);
272 }
273
destroy_grad_list(Vec_Chain * head)274 static PetscErrorCode destroy_grad_list(Vec_Chain *head)
275 {
276 Vec_Chain *p = head->next, *q;
277
278 PetscFunctionBegin;
279 while (p) {
280 q = p->next;
281 PetscCall(VecDestroy(&p->V));
282 PetscCall(PetscFree(p));
283 p = q;
284 }
285 head->next = NULL;
286 PetscFunctionReturn(PETSC_SUCCESS);
287 }
288
TaoSolve_BMRM(Tao tao)289 static PetscErrorCode TaoSolve_BMRM(Tao tao)
290 {
291 TAO_DF df;
292 TAO_BMRM *bmrm = (TAO_BMRM *)tao->data;
293
294 /* Values and pointers to parts of the optimization problem */
295 PetscReal f = 0.0;
296 Vec W = tao->solution;
297 Vec G = tao->gradient;
298 PetscReal lambda;
299 PetscReal bt;
300 Vec_Chain grad_list, *tail_glist, *pgrad;
301 PetscInt i;
302 PetscMPIInt rank;
303
304 /* Used in converged criteria check */
305 PetscReal reg;
306 PetscReal jtwt = 0.0, max_jtwt, pre_epsilon, epsilon, jw, min_jw;
307 PetscReal innerSolverTol;
308 MPI_Comm comm;
309
310 PetscFunctionBegin;
311 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm));
312 PetscCallMPI(MPI_Comm_rank(comm, &rank));
313 lambda = bmrm->lambda;
314
315 /* Check Stopping Condition */
316 tao->step = 1.0;
317 max_jtwt = -BMRM_INFTY;
318 min_jw = BMRM_INFTY;
319 innerSolverTol = 1.0;
320 epsilon = 0.0;
321
322 if (rank == 0) {
323 PetscCall(init_df_solver(&df));
324 grad_list.next = NULL;
325 tail_glist = &grad_list;
326 }
327
328 df.tol = 1e-6;
329 tao->reason = TAO_CONTINUE_ITERATING;
330
331 /*-----------------Algorithm Begins------------------------*/
332 /* make the scatter */
333 PetscCall(VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w));
334 PetscCall(VecAssemblyBegin(bmrm->local_w));
335 PetscCall(VecAssemblyEnd(bmrm->local_w));
336
337 /* NOTE: In application pass the sub-gradient of Remp(W) */
338 PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G));
339 PetscCall(TaoLogConvergenceHistory(tao, f, 1.0, 0.0, tao->ksp_its));
340 PetscCall(TaoMonitor(tao, tao->niter, f, 1.0, 0.0, tao->step));
341 PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
342
343 while (tao->reason == TAO_CONTINUE_ITERATING) {
344 /* Call general purpose update function */
345 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
346
347 /* compute bt = Remp(Wt-1) - <Wt-1, At> */
348 PetscCall(VecDot(W, G, &bt));
349 bt = f - bt;
350
351 /* First gather the gradient to the rank-0 node */
352 PetscCall(VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD));
353 PetscCall(VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD));
354
355 /* Bring up the inner solver */
356 if (rank == 0) {
357 PetscCall(ensure_df_space(tao->niter + 1, &df));
358 PetscCall(make_grad_node(bmrm->local_w, &pgrad));
359 tail_glist->next = pgrad;
360 tail_glist = pgrad;
361
362 df.a[tao->niter] = 1.0;
363 df.f[tao->niter] = -bt;
364 df.u[tao->niter] = 1.0;
365 df.l[tao->niter] = 0.0;
366
367 /* set up the Q */
368 pgrad = grad_list.next;
369 for (i = 0; i <= tao->niter; i++) {
370 PetscCheck(pgrad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assert that there are at least tao->niter+1 pgrad available");
371 PetscCall(VecDot(pgrad->V, bmrm->local_w, ®));
372 df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda;
373 pgrad = pgrad->next;
374 }
375
376 if (tao->niter > 0) {
377 df.x[tao->niter] = 0.0;
378 PetscCall(solve(&df));
379 } else df.x[0] = 1.0;
380
381 /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */
382 jtwt = 0.0;
383 PetscCall(VecSet(bmrm->local_w, 0.0));
384 pgrad = grad_list.next;
385 for (i = 0; i <= tao->niter; i++) {
386 jtwt -= df.x[i] * df.f[i];
387 PetscCall(VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V));
388 pgrad = pgrad->next;
389 }
390
391 PetscCall(VecNorm(bmrm->local_w, NORM_2, ®));
392 reg = 0.5 * lambda * reg * reg;
393 jtwt -= reg;
394 } /* end if rank == 0 */
395
396 /* scatter the new W to all nodes */
397 PetscCall(VecScatterBegin(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE));
398 PetscCall(VecScatterEnd(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE));
399
400 PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G));
401
402 PetscCallMPI(MPI_Bcast(&jtwt, 1, MPIU_REAL, 0, comm));
403 PetscCallMPI(MPI_Bcast(®, 1, MPIU_REAL, 0, comm));
404
405 jw = reg + f; /* J(w) = regularizer + Remp(w) */
406 if (jw < min_jw) min_jw = jw;
407 if (jtwt > max_jtwt) max_jtwt = jtwt;
408
409 pre_epsilon = epsilon;
410 epsilon = min_jw - jtwt;
411
412 if (rank == 0) {
413 if (innerSolverTol > epsilon) innerSolverTol = epsilon;
414 else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7;
415
416 /* if the annealing doesn't work well, lower the inner solver tolerance */
417 if (pre_epsilon < epsilon) innerSolverTol *= 0.2;
418
419 df.tol = innerSolverTol * 0.5;
420 }
421
422 tao->niter++;
423 PetscCall(TaoLogConvergenceHistory(tao, min_jw, epsilon, 0.0, tao->ksp_its));
424 PetscCall(TaoMonitor(tao, tao->niter, min_jw, epsilon, 0.0, tao->step));
425 PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
426 }
427
428 /* free all the memory */
429 if (rank == 0) {
430 PetscCall(destroy_grad_list(&grad_list));
431 PetscCall(destroy_df_solver(&df));
432 }
433
434 PetscCall(VecDestroy(&bmrm->local_w));
435 PetscCall(VecScatterDestroy(&bmrm->scatter));
436 PetscFunctionReturn(PETSC_SUCCESS);
437 }
438
TaoSetup_BMRM(Tao tao)439 static PetscErrorCode TaoSetup_BMRM(Tao tao)
440 {
441 PetscFunctionBegin;
442 /* Allocate some arrays */
443 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
444 PetscFunctionReturn(PETSC_SUCCESS);
445 }
446
TaoDestroy_BMRM(Tao tao)447 static PetscErrorCode TaoDestroy_BMRM(Tao tao)
448 {
449 PetscFunctionBegin;
450 PetscCall(PetscFree(tao->data));
451 PetscFunctionReturn(PETSC_SUCCESS);
452 }
453
TaoSetFromOptions_BMRM(Tao tao,PetscOptionItems PetscOptionsObject)454 static PetscErrorCode TaoSetFromOptions_BMRM(Tao tao, PetscOptionItems PetscOptionsObject)
455 {
456 TAO_BMRM *bmrm = (TAO_BMRM *)tao->data;
457
458 PetscFunctionBegin;
459 PetscOptionsHeadBegin(PetscOptionsObject, "BMRM for regularized risk minimization");
460 PetscCall(PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight", "", 100, &bmrm->lambda, NULL));
461 PetscOptionsHeadEnd();
462 PetscFunctionReturn(PETSC_SUCCESS);
463 }
464
TaoView_BMRM(Tao tao,PetscViewer viewer)465 static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer)
466 {
467 PetscBool isascii;
468
469 PetscFunctionBegin;
470 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
471 if (isascii) {
472 PetscCall(PetscViewerASCIIPushTab(viewer));
473 PetscCall(PetscViewerASCIIPopTab(viewer));
474 }
475 PetscFunctionReturn(PETSC_SUCCESS);
476 }
477
478 /*MC
479 TAOBMRM - bundle method for regularized risk minimization
480
481 Options Database Keys:
482 . - tao_bmrm_lambda - regulariser weight
483
484 Level: beginner
485 M*/
486
TaoCreate_BMRM(Tao tao)487 PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao)
488 {
489 TAO_BMRM *bmrm;
490
491 PetscFunctionBegin;
492 tao->ops->setup = TaoSetup_BMRM;
493 tao->ops->solve = TaoSolve_BMRM;
494 tao->ops->view = TaoView_BMRM;
495 tao->ops->setfromoptions = TaoSetFromOptions_BMRM;
496 tao->ops->destroy = TaoDestroy_BMRM;
497
498 PetscCall(PetscNew(&bmrm));
499 bmrm->lambda = 1.0;
500 tao->data = (void *)bmrm;
501
502 /* Override default settings (unless already changed) */
503 PetscCall(TaoParametersInitialize(tao));
504 PetscObjectParameterSetDefault(tao, max_it, 2000);
505 PetscObjectParameterSetDefault(tao, max_funcs, 4000);
506 PetscObjectParameterSetDefault(tao, gatol, 1.0e-12);
507 PetscObjectParameterSetDefault(tao, grtol, 1.0e-12);
508 PetscFunctionReturn(PETSC_SUCCESS);
509 }
510
init_df_solver(TAO_DF * df)511 static PetscErrorCode init_df_solver(TAO_DF *df)
512 {
513 PetscInt i, n = INCRE_DIM;
514
515 PetscFunctionBegin;
516 /* default values */
517 df->maxProjIter = 200;
518 df->maxPGMIter = 300000;
519 df->b = 1.0;
520
521 /* memory space required by Dai-Fletcher */
522 df->cur_num_cp = n;
523 PetscCall(PetscMalloc1(n, &df->f));
524 PetscCall(PetscMalloc1(n, &df->a));
525 PetscCall(PetscMalloc1(n, &df->l));
526 PetscCall(PetscMalloc1(n, &df->u));
527 PetscCall(PetscMalloc1(n, &df->x));
528 PetscCall(PetscMalloc1(n, &df->Q));
529
530 for (i = 0; i < n; i++) PetscCall(PetscMalloc1(n, &df->Q[i]));
531
532 PetscCall(PetscMalloc1(n, &df->g));
533 PetscCall(PetscMalloc1(n, &df->y));
534 PetscCall(PetscMalloc1(n, &df->tempv));
535 PetscCall(PetscMalloc1(n, &df->d));
536 PetscCall(PetscMalloc1(n, &df->Qd));
537 PetscCall(PetscMalloc1(n, &df->t));
538 PetscCall(PetscMalloc1(n, &df->xplus));
539 PetscCall(PetscMalloc1(n, &df->tplus));
540 PetscCall(PetscMalloc1(n, &df->sk));
541 PetscCall(PetscMalloc1(n, &df->yk));
542
543 PetscCall(PetscMalloc1(n, &df->ipt));
544 PetscCall(PetscMalloc1(n, &df->ipt2));
545 PetscCall(PetscMalloc1(n, &df->uv));
546 PetscFunctionReturn(PETSC_SUCCESS);
547 }
548
ensure_df_space(PetscInt dim,TAO_DF * df)549 static PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df)
550 {
551 PetscReal *tmp, **tmp_Q;
552 PetscInt i, n, old_n;
553
554 PetscFunctionBegin;
555 df->dim = dim;
556 if (dim <= df->cur_num_cp) PetscFunctionReturn(PETSC_SUCCESS);
557
558 old_n = df->cur_num_cp;
559 df->cur_num_cp += INCRE_DIM;
560 n = df->cur_num_cp;
561
562 /* memory space required by dai-fletcher */
563 PetscCall(PetscMalloc1(n, &tmp));
564 PetscCall(PetscArraycpy(tmp, df->f, old_n));
565 PetscCall(PetscFree(df->f));
566 df->f = tmp;
567
568 PetscCall(PetscMalloc1(n, &tmp));
569 PetscCall(PetscArraycpy(tmp, df->a, old_n));
570 PetscCall(PetscFree(df->a));
571 df->a = tmp;
572
573 PetscCall(PetscMalloc1(n, &tmp));
574 PetscCall(PetscArraycpy(tmp, df->l, old_n));
575 PetscCall(PetscFree(df->l));
576 df->l = tmp;
577
578 PetscCall(PetscMalloc1(n, &tmp));
579 PetscCall(PetscArraycpy(tmp, df->u, old_n));
580 PetscCall(PetscFree(df->u));
581 df->u = tmp;
582
583 PetscCall(PetscMalloc1(n, &tmp));
584 PetscCall(PetscArraycpy(tmp, df->x, old_n));
585 PetscCall(PetscFree(df->x));
586 df->x = tmp;
587
588 PetscCall(PetscMalloc1(n, &tmp_Q));
589 for (i = 0; i < n; i++) {
590 PetscCall(PetscMalloc1(n, &tmp_Q[i]));
591 if (i < old_n) {
592 PetscCall(PetscArraycpy(tmp_Q[i], df->Q[i], old_n));
593 PetscCall(PetscFree(df->Q[i]));
594 }
595 }
596
597 PetscCall(PetscFree(df->Q));
598 df->Q = tmp_Q;
599
600 PetscCall(PetscFree(df->g));
601 PetscCall(PetscMalloc1(n, &df->g));
602
603 PetscCall(PetscFree(df->y));
604 PetscCall(PetscMalloc1(n, &df->y));
605
606 PetscCall(PetscFree(df->tempv));
607 PetscCall(PetscMalloc1(n, &df->tempv));
608
609 PetscCall(PetscFree(df->d));
610 PetscCall(PetscMalloc1(n, &df->d));
611
612 PetscCall(PetscFree(df->Qd));
613 PetscCall(PetscMalloc1(n, &df->Qd));
614
615 PetscCall(PetscFree(df->t));
616 PetscCall(PetscMalloc1(n, &df->t));
617
618 PetscCall(PetscFree(df->xplus));
619 PetscCall(PetscMalloc1(n, &df->xplus));
620
621 PetscCall(PetscFree(df->tplus));
622 PetscCall(PetscMalloc1(n, &df->tplus));
623
624 PetscCall(PetscFree(df->sk));
625 PetscCall(PetscMalloc1(n, &df->sk));
626
627 PetscCall(PetscFree(df->yk));
628 PetscCall(PetscMalloc1(n, &df->yk));
629
630 PetscCall(PetscFree(df->ipt));
631 PetscCall(PetscMalloc1(n, &df->ipt));
632
633 PetscCall(PetscFree(df->ipt2));
634 PetscCall(PetscMalloc1(n, &df->ipt2));
635
636 PetscCall(PetscFree(df->uv));
637 PetscCall(PetscMalloc1(n, &df->uv));
638 PetscFunctionReturn(PETSC_SUCCESS);
639 }
640
destroy_df_solver(TAO_DF * df)641 static PetscErrorCode destroy_df_solver(TAO_DF *df)
642 {
643 PetscInt i;
644
645 PetscFunctionBegin;
646 PetscCall(PetscFree(df->f));
647 PetscCall(PetscFree(df->a));
648 PetscCall(PetscFree(df->l));
649 PetscCall(PetscFree(df->u));
650 PetscCall(PetscFree(df->x));
651
652 for (i = 0; i < df->cur_num_cp; i++) PetscCall(PetscFree(df->Q[i]));
653 PetscCall(PetscFree(df->Q));
654 PetscCall(PetscFree(df->ipt));
655 PetscCall(PetscFree(df->ipt2));
656 PetscCall(PetscFree(df->uv));
657 PetscCall(PetscFree(df->g));
658 PetscCall(PetscFree(df->y));
659 PetscCall(PetscFree(df->tempv));
660 PetscCall(PetscFree(df->d));
661 PetscCall(PetscFree(df->Qd));
662 PetscCall(PetscFree(df->t));
663 PetscCall(PetscFree(df->xplus));
664 PetscCall(PetscFree(df->tplus));
665 PetscCall(PetscFree(df->sk));
666 PetscCall(PetscFree(df->yk));
667 PetscFunctionReturn(PETSC_SUCCESS);
668 }
669
670 /* Piecewise linear monotone target function for the Dai-Fletcher projector */
phi(PetscReal * x,PetscInt n,PetscReal lambda,PetscReal * a,PetscReal b,PetscReal * c,PetscReal * l,PetscReal * u)671 static PetscReal phi(PetscReal *x, PetscInt n, PetscReal lambda, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u)
672 {
673 PetscReal r = 0.0;
674 PetscInt i;
675
676 for (i = 0; i < n; i++) {
677 x[i] = -c[i] + lambda * a[i];
678 if (x[i] > u[i]) x[i] = u[i];
679 else if (x[i] < l[i]) x[i] = l[i];
680 r += a[i] * x[i];
681 }
682 return r - b;
683 }
684
685 /** Modified Dai-Fletcher QP projector solves the problem:
686 *
687 * minimise 0.5*x'*x - c'*x
688 * subj to a'*x = b
689 * l \leq x \leq u
690 *
691 * \param c The point to be projected onto feasible set
692 */
project(PetscInt n,PetscReal * a,PetscReal b,PetscReal * c,PetscReal * l,PetscReal * u,PetscReal * x,PetscReal * lam_ext,TAO_DF * df)693 static PetscInt project(PetscInt n, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u, PetscReal *x, PetscReal *lam_ext, TAO_DF *df)
694 {
695 PetscReal lambda, lambdal, lambdau, dlambda, lambda_new;
696 PetscReal r, rl, ru, s;
697 PetscInt innerIter;
698 PetscBool nonNegativeSlack = PETSC_FALSE;
699
700 *lam_ext = 0;
701 lambda = 0;
702 dlambda = 0.5;
703 innerIter = 1;
704
705 /* \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b)
706 *
707 * Optimality conditions for \phi:
708 *
709 * 1. lambda <= 0
710 * 2. r <= 0
711 * 3. r*lambda == 0
712 */
713
714 /* Bracketing Phase */
715 r = phi(x, n, lambda, a, b, c, l, u);
716
717 if (nonNegativeSlack) {
718 /* inequality constraint, i.e., with \xi >= 0 constraint */
719 if (r < TOL_R) return PETSC_SUCCESS;
720 } else {
721 /* equality constraint ,i.e., without \xi >= 0 constraint */
722 if (PetscAbsReal(r) < TOL_R) return PETSC_SUCCESS;
723 }
724
725 if (r < 0.0) {
726 lambdal = lambda;
727 rl = r;
728 lambda = lambda + dlambda;
729 r = phi(x, n, lambda, a, b, c, l, u);
730 while (r < 0.0 && dlambda < BMRM_INFTY) {
731 lambdal = lambda;
732 s = rl / r - 1.0;
733 if (s < 0.1) s = 0.1;
734 dlambda = dlambda + dlambda / s;
735 lambda = lambda + dlambda;
736 rl = r;
737 r = phi(x, n, lambda, a, b, c, l, u);
738 }
739 lambdau = lambda;
740 ru = r;
741 } else {
742 lambdau = lambda;
743 ru = r;
744 lambda = lambda - dlambda;
745 r = phi(x, n, lambda, a, b, c, l, u);
746 while (r > 0.0 && dlambda > -BMRM_INFTY) {
747 lambdau = lambda;
748 s = ru / r - 1.0;
749 if (s < 0.1) s = 0.1;
750 dlambda = dlambda + dlambda / s;
751 lambda = lambda - dlambda;
752 ru = r;
753 r = phi(x, n, lambda, a, b, c, l, u);
754 }
755 lambdal = lambda;
756 rl = r;
757 }
758
759 PetscCheckAbort(PetscAbsReal(dlambda) <= BMRM_INFTY, PETSC_COMM_SELF, PETSC_ERR_PLIB, "L2N2_DaiFletcherPGM detected Infeasible QP problem!");
760
761 if (ru == 0) return innerIter;
762
763 /* Secant Phase */
764 s = 1.0 - rl / ru;
765 dlambda = dlambda / s;
766 lambda = lambdau - dlambda;
767 r = phi(x, n, lambda, a, b, c, l, u);
768
769 while (PetscAbsReal(r) > TOL_R && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda)) && innerIter < df->maxProjIter) {
770 innerIter++;
771 if (r > 0.0) {
772 if (s <= 2.0) {
773 lambdau = lambda;
774 ru = r;
775 s = 1.0 - rl / ru;
776 dlambda = (lambdau - lambdal) / s;
777 lambda = lambdau - dlambda;
778 } else {
779 s = ru / r - 1.0;
780 if (s < 0.1) s = 0.1;
781 dlambda = (lambdau - lambda) / s;
782 lambda_new = 0.75 * lambdal + 0.25 * lambda;
783 if (lambda_new < (lambda - dlambda)) lambda_new = lambda - dlambda;
784 lambdau = lambda;
785 ru = r;
786 lambda = lambda_new;
787 s = (lambdau - lambdal) / (lambdau - lambda);
788 }
789 } else {
790 if (s >= 2.0) {
791 lambdal = lambda;
792 rl = r;
793 s = 1.0 - rl / ru;
794 dlambda = (lambdau - lambdal) / s;
795 lambda = lambdau - dlambda;
796 } else {
797 s = rl / r - 1.0;
798 if (s < 0.1) s = 0.1;
799 dlambda = (lambda - lambdal) / s;
800 lambda_new = 0.75 * lambdau + 0.25 * lambda;
801 if (lambda_new > (lambda + dlambda)) lambda_new = lambda + dlambda;
802 lambdal = lambda;
803 rl = r;
804 lambda = lambda_new;
805 s = (lambdau - lambdal) / (lambdau - lambda);
806 }
807 }
808 r = phi(x, n, lambda, a, b, c, l, u);
809 }
810
811 *lam_ext = lambda;
812 if (innerIter >= df->maxProjIter) PetscCallAbort(PETSC_COMM_SELF, PetscInfo(NULL, "WARNING: DaiFletcher max iterations\n"));
813 return innerIter;
814 }
815