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