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