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