Actual source code: bmrm.c
petsc-3.8.4 2018-03-24
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: TaoConvergedReason reason;
53: TAO_DF df;
54: TAO_BMRM *bmrm = (TAO_BMRM*)tao->data;
56: /* Values and pointers to parts of the optimization problem */
57: PetscReal f = 0.0;
58: Vec W = tao->solution;
59: Vec G = tao->gradient;
60: PetscReal lambda;
61: PetscReal bt;
62: Vec_Chain grad_list, *tail_glist, *pgrad;
63: PetscInt i;
64: PetscMPIInt rank;
66: /* Used in converged criteria check */
67: PetscReal reg;
68: PetscReal jtwt, max_jtwt, pre_epsilon, epsilon, jw, min_jw;
69: PetscReal innerSolverTol;
70: MPI_Comm comm;
73: PetscObjectGetComm((PetscObject)tao,&comm);
74: MPI_Comm_rank(comm, &rank);
75: lambda = bmrm->lambda;
77: /* Check Stopping Condition */
78: tao->step = 1.0;
79: max_jtwt = -BMRM_INFTY;
80: min_jw = BMRM_INFTY;
81: innerSolverTol = 1.0;
82: epsilon = 0.0;
84: if (!rank) {
85: init_df_solver(&df);
86: grad_list.next = NULL;
87: tail_glist = &grad_list;
88: }
90: df.tol = 1e-6;
91: reason = TAO_CONTINUE_ITERATING;
93: /*-----------------Algorithm Begins------------------------*/
94: /* make the scatter */
95: VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w);
96: VecAssemblyBegin(bmrm->local_w);
97: VecAssemblyEnd(bmrm->local_w);
99: /* NOTE: In application pass the sub-gradient of Remp(W) */
100: TaoComputeObjectiveAndGradient(tao, W, &f, G);
101: TaoMonitor(tao,tao->niter,f,1.0,0.0,tao->step,&reason);
102: while (reason == TAO_CONTINUE_ITERATING) {
103: /* compute bt = Remp(Wt-1) - <Wt-1, At> */
104: VecDot(W, G, &bt);
105: bt = f - bt;
107: /* First gather the gradient to the master node */
108: VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
109: VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
111: /* Bring up the inner solver */
112: if (!rank) {
113: ensure_df_space(tao->niter+1, &df);
114: make_grad_node(bmrm->local_w, &pgrad);
115: tail_glist->next = pgrad;
116: tail_glist = pgrad;
118: df.a[tao->niter] = 1.0;
119: df.f[tao->niter] = -bt;
120: df.u[tao->niter] = 1.0;
121: df.l[tao->niter] = 0.0;
123: /* set up the Q */
124: pgrad = grad_list.next;
125: for (i=0; i<=tao->niter; i++) {
126: if (!pgrad) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Assert that there are at least tao->niter+1 pgrad available");
127: VecDot(pgrad->V, bmrm->local_w, ®);
128: df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda;
129: pgrad = pgrad->next;
130: }
132: if (tao->niter > 0) {
133: df.x[tao->niter] = 0.0;
134: solve(&df);
135: } else
136: df.x[0] = 1.0;
138: /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */
139: jtwt = 0.0;
140: VecSet(bmrm->local_w, 0.0);
141: pgrad = grad_list.next;
142: for (i=0; i<=tao->niter; i++) {
143: jtwt -= df.x[i] * df.f[i];
144: VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V);
145: pgrad = pgrad->next;
146: }
148: VecNorm(bmrm->local_w, NORM_2, ®);
149: reg = 0.5*lambda*reg*reg;
150: jtwt -= reg;
151: } /* end if rank == 0 */
153: /* scatter the new W to all nodes */
154: VecScatterBegin(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);
155: VecScatterEnd(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);
157: TaoComputeObjectiveAndGradient(tao, W, &f, G);
159: MPI_Bcast(&jtwt,1,MPIU_REAL,0,comm);
160: MPI_Bcast(®,1,MPIU_REAL,0,comm);
162: jw = reg + f; /* J(w) = regularizer + Remp(w) */
163: if (jw < min_jw) min_jw = jw;
164: if (jtwt > max_jtwt) max_jtwt = jtwt;
166: pre_epsilon = epsilon;
167: epsilon = min_jw - jtwt;
169: if (!rank) {
170: if (innerSolverTol > epsilon) innerSolverTol = epsilon;
171: else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7;
173: /* if the annealing doesn't work well, lower the inner solver tolerance */
174: if(pre_epsilon < epsilon) innerSolverTol *= 0.2;
176: df.tol = innerSolverTol*0.5;
177: }
179: tao->niter++;
180: TaoMonitor(tao,tao->niter,min_jw,epsilon,0.0,tao->step,&reason);
181: }
183: /* free all the memory */
184: if (!rank) {
185: destroy_grad_list(&grad_list);
186: destroy_df_solver(&df);
187: }
189: VecDestroy(&bmrm->local_w);
190: VecScatterDestroy(&bmrm->scatter);
191: return(0);
192: }
195: /* ---------------------------------------------------------- */
197: static PetscErrorCode TaoSetup_BMRM(Tao tao)
198: {
203: /* Allocate some arrays */
204: if (!tao->gradient) {
205: VecDuplicate(tao->solution, &tao->gradient);
206: }
207: return(0);
208: }
210: /*------------------------------------------------------------*/
211: static PetscErrorCode TaoDestroy_BMRM(Tao tao)
212: {
216: PetscFree(tao->data);
217: return(0);
218: }
220: static PetscErrorCode TaoSetFromOptions_BMRM(PetscOptionItems *PetscOptionsObject,Tao tao)
221: {
223: TAO_BMRM* bmrm = (TAO_BMRM*)tao->data;
226: PetscOptionsHead(PetscOptionsObject,"BMRM for regularized risk minimization");
227: PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight","", 100,&bmrm->lambda,NULL);
228: PetscOptionsTail();
229: return(0);
230: }
232: /*------------------------------------------------------------*/
233: static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer)
234: {
235: PetscBool isascii;
239: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
240: if (isascii) {
241: PetscViewerASCIIPushTab(viewer);
242: PetscViewerASCIIPopTab(viewer);
243: }
244: return(0);
245: }
247: /*------------------------------------------------------------*/
248: /*MC
249: TAOBMRM - bundle method for regularized risk minimization
251: Options Database Keys:
252: . - tao_bmrm_lambda - regulariser weight
254: Level: beginner
255: M*/
257: PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao)
258: {
259: TAO_BMRM *bmrm;
263: tao->ops->setup = TaoSetup_BMRM;
264: tao->ops->solve = TaoSolve_BMRM;
265: tao->ops->view = TaoView_BMRM;
266: tao->ops->setfromoptions = TaoSetFromOptions_BMRM;
267: tao->ops->destroy = TaoDestroy_BMRM;
269: PetscNewLog(tao,&bmrm);
270: bmrm->lambda = 1.0;
271: tao->data = (void*)bmrm;
273: /* Override default settings (unless already changed) */
274: if (!tao->max_it_changed) tao->max_it = 2000;
275: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
276: if (!tao->gatol_changed) tao->gatol = 1.0e-12;
277: if (!tao->grtol_changed) tao->grtol = 1.0e-12;
279: return(0);
280: }
282: PetscErrorCode init_df_solver(TAO_DF *df)
283: {
284: PetscInt i, n = INCRE_DIM;
288: /* default values */
289: df->maxProjIter = 200;
290: df->maxPGMIter = 300000;
291: df->b = 1.0;
293: /* memory space required by Dai-Fletcher */
294: df->cur_num_cp = n;
295: PetscMalloc1(n, &df->f);
296: PetscMalloc1(n, &df->a);
297: PetscMalloc1(n, &df->l);
298: PetscMalloc1(n, &df->u);
299: PetscMalloc1(n, &df->x);
300: PetscMalloc1(n, &df->Q);
302: for (i = 0; i < n; i ++) {
303: PetscMalloc1(n, &df->Q[i]);
304: }
306: PetscMalloc1(n, &df->g);
307: PetscMalloc1(n, &df->y);
308: PetscMalloc1(n, &df->tempv);
309: PetscMalloc1(n, &df->d);
310: PetscMalloc1(n, &df->Qd);
311: PetscMalloc1(n, &df->t);
312: PetscMalloc1(n, &df->xplus);
313: PetscMalloc1(n, &df->tplus);
314: PetscMalloc1(n, &df->sk);
315: PetscMalloc1(n, &df->yk);
317: PetscMalloc1(n, &df->ipt);
318: PetscMalloc1(n, &df->ipt2);
319: PetscMalloc1(n, &df->uv);
320: return(0);
321: }
323: PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df)
324: {
326: PetscReal *tmp, **tmp_Q;
327: PetscInt i, n, old_n;
330: df->dim = dim;
331: if (dim <= df->cur_num_cp) return(0);
333: old_n = df->cur_num_cp;
334: df->cur_num_cp += INCRE_DIM;
335: n = df->cur_num_cp;
337: /* memory space required by dai-fletcher */
338: PetscMalloc1(n, &tmp);
339: PetscMemcpy(tmp, df->f, sizeof(PetscReal)*old_n);
340: PetscFree(df->f);
341: df->f = tmp;
343: PetscMalloc1(n, &tmp);
344: PetscMemcpy(tmp, df->a, sizeof(PetscReal)*old_n);
345: PetscFree(df->a);
346: df->a = tmp;
348: PetscMalloc1(n, &tmp);
349: PetscMemcpy(tmp, df->l, sizeof(PetscReal)*old_n);
350: PetscFree(df->l);
351: df->l = tmp;
353: PetscMalloc1(n, &tmp);
354: PetscMemcpy(tmp, df->u, sizeof(PetscReal)*old_n);
355: PetscFree(df->u);
356: df->u = tmp;
358: PetscMalloc1(n, &tmp);
359: PetscMemcpy(tmp, df->x, sizeof(PetscReal)*old_n);
360: PetscFree(df->x);
361: df->x = tmp;
363: PetscMalloc1(n, &tmp_Q);
364: for (i = 0; i < n; i ++) {
365: PetscMalloc1(n, &tmp_Q[i]);
366: if (i < old_n) {
367: PetscMemcpy(tmp_Q[i], df->Q[i], sizeof(PetscReal)*old_n);
368: PetscFree(df->Q[i]);
369: }
370: }
372: PetscFree(df->Q);
373: df->Q = tmp_Q;
375: PetscFree(df->g);
376: PetscMalloc1(n, &df->g);
378: PetscFree(df->y);
379: PetscMalloc1(n, &df->y);
381: PetscFree(df->tempv);
382: PetscMalloc1(n, &df->tempv);
384: PetscFree(df->d);
385: PetscMalloc1(n, &df->d);
387: PetscFree(df->Qd);
388: PetscMalloc1(n, &df->Qd);
390: PetscFree(df->t);
391: PetscMalloc1(n, &df->t);
393: PetscFree(df->xplus);
394: PetscMalloc1(n, &df->xplus);
396: PetscFree(df->tplus);
397: PetscMalloc1(n, &df->tplus);
399: PetscFree(df->sk);
400: PetscMalloc1(n, &df->sk);
402: PetscFree(df->yk);
403: PetscMalloc1(n, &df->yk);
405: PetscFree(df->ipt);
406: PetscMalloc1(n, &df->ipt);
408: PetscFree(df->ipt2);
409: PetscMalloc1(n, &df->ipt2);
411: PetscFree(df->uv);
412: PetscMalloc1(n, &df->uv);
413: return(0);
414: }
416: PetscErrorCode destroy_df_solver(TAO_DF *df)
417: {
419: PetscInt i;
422: PetscFree(df->f);
423: PetscFree(df->a);
424: PetscFree(df->l);
425: PetscFree(df->u);
426: PetscFree(df->x);
428: for (i = 0; i < df->cur_num_cp; i ++) {
429: PetscFree(df->Q[i]);
430: }
431: PetscFree(df->Q);
432: PetscFree(df->ipt);
433: PetscFree(df->ipt2);
434: PetscFree(df->uv);
435: PetscFree(df->g);
436: PetscFree(df->y);
437: PetscFree(df->tempv);
438: PetscFree(df->d);
439: PetscFree(df->Qd);
440: PetscFree(df->t);
441: PetscFree(df->xplus);
442: PetscFree(df->tplus);
443: PetscFree(df->sk);
444: PetscFree(df->yk);
445: return(0);
446: }
448: /* Piecewise linear monotone target function for the Dai-Fletcher projector */
449: PetscReal phi(PetscReal *x,PetscInt n,PetscReal lambda,PetscReal *a,PetscReal b,PetscReal *c,PetscReal *l,PetscReal *u)
450: {
451: PetscReal r = 0.0;
452: PetscInt i;
454: for (i = 0; i < n; i++){
455: x[i] = -c[i] + lambda*a[i];
456: if (x[i] > u[i]) x[i] = u[i];
457: else if(x[i] < l[i]) x[i] = l[i];
458: r += a[i]*x[i];
459: }
460: return r - b;
461: }
463: /** Modified Dai-Fletcher QP projector solves the problem:
464: *
465: * minimise 0.5*x'*x - c'*x
466: * subj to a'*x = b
467: * l \leq x \leq u
468: *
469: * \param c The point to be projected onto feasible set
470: */
471: PetscInt project(PetscInt n,PetscReal *a,PetscReal b,PetscReal *c,PetscReal *l,PetscReal *u,PetscReal *x,PetscReal *lam_ext,TAO_DF *df)
472: {
473: PetscReal lambda, lambdal, lambdau, dlambda, lambda_new;
474: PetscReal r, rl, ru, s;
475: PetscInt innerIter;
476: PetscBool nonNegativeSlack = PETSC_FALSE;
479: *lam_ext = 0;
480: lambda = 0;
481: dlambda = 0.5;
482: innerIter = 1;
484: /* \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b)
485: *
486: * Optimality conditions for \phi:
487: *
488: * 1. lambda <= 0
489: * 2. r <= 0
490: * 3. r*lambda == 0
491: */
493: /* Bracketing Phase */
494: r = phi(x, n, lambda, a, b, c, l, u);
496: if(nonNegativeSlack) {
497: /* inequality constraint, i.e., with \xi >= 0 constraint */
498: if (r < TOL_R) return 0;
499: } else {
500: /* equality constraint ,i.e., without \xi >= 0 constraint */
501: if (PetscAbsReal(r) < TOL_R) return 0;
502: }
504: if (r < 0.0){
505: lambdal = lambda;
506: rl = r;
507: lambda = lambda + dlambda;
508: r = phi(x, n, lambda, a, b, c, l, u);
509: while (r < 0.0 && dlambda < BMRM_INFTY) {
510: lambdal = lambda;
511: s = rl/r - 1.0;
512: if (s < 0.1) s = 0.1;
513: dlambda = dlambda + dlambda/s;
514: lambda = lambda + dlambda;
515: rl = r;
516: r = phi(x, n, lambda, a, b, c, l, u);
517: }
518: lambdau = lambda;
519: ru = r;
520: } else {
521: lambdau = lambda;
522: ru = r;
523: lambda = lambda - dlambda;
524: r = phi(x, n, lambda, a, b, c, l, u);
525: while (r > 0.0 && dlambda > -BMRM_INFTY) {
526: lambdau = lambda;
527: s = ru/r - 1.0;
528: if (s < 0.1) s = 0.1;
529: dlambda = dlambda + dlambda/s;
530: lambda = lambda - dlambda;
531: ru = r;
532: r = phi(x, n, lambda, a, b, c, l, u);
533: }
534: lambdal = lambda;
535: rl = r;
536: }
538: if(PetscAbsReal(dlambda) > BMRM_INFTY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"L2N2_DaiFletcherPGM detected Infeasible QP problem!");
540: if(ru == 0){
541: return innerIter;
542: }
544: /* Secant Phase */
545: s = 1.0 - rl/ru;
546: dlambda = dlambda/s;
547: lambda = lambdau - dlambda;
548: r = phi(x, n, lambda, a, b, c, l, u);
550: while (PetscAbsReal(r) > TOL_R
551: && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda))
552: && innerIter < df->maxProjIter){
553: innerIter++;
554: if (r > 0.0){
555: if (s <= 2.0){
556: lambdau = lambda;
557: ru = r;
558: s = 1.0 - rl/ru;
559: dlambda = (lambdau - lambdal) / s;
560: lambda = lambdau - dlambda;
561: } else {
562: s = ru/r-1.0;
563: if (s < 0.1) s = 0.1;
564: dlambda = (lambdau - lambda) / s;
565: lambda_new = 0.75*lambdal + 0.25*lambda;
566: if (lambda_new < (lambda - dlambda))
567: lambda_new = lambda - dlambda;
568: lambdau = lambda;
569: ru = r;
570: lambda = lambda_new;
571: s = (lambdau - lambdal) / (lambdau - lambda);
572: }
573: } else {
574: if (s >= 2.0){
575: lambdal = lambda;
576: rl = r;
577: s = 1.0 - rl/ru;
578: dlambda = (lambdau - lambdal) / s;
579: lambda = lambdau - dlambda;
580: } else {
581: s = rl/r - 1.0;
582: if (s < 0.1) s = 0.1;
583: dlambda = (lambda-lambdal) / s;
584: lambda_new = 0.75*lambdau + 0.25*lambda;
585: if (lambda_new > (lambda + dlambda))
586: lambda_new = lambda + dlambda;
587: lambdal = lambda;
588: rl = r;
589: lambda = lambda_new;
590: s = (lambdau - lambdal) / (lambdau-lambda);
591: }
592: }
593: r = phi(x, n, lambda, a, b, c, l, u);
594: }
596: *lam_ext = lambda;
597: if(innerIter >= df->maxProjIter) {
598: PetscPrintf(PETSC_COMM_SELF, "WARNING: DaiFletcher max iterations\n");
599: }
600: return innerIter;
601: }
604: PetscErrorCode solve(TAO_DF *df)
605: {
607: PetscInt i, j, innerIter, it, it2, luv, info, lscount = 0, projcount = 0;
608: PetscReal gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam=0.0, lam_ext;
609: PetscReal DELTAsv, ProdDELTAsv;
610: PetscReal c, *tempQ;
611: PetscReal *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol;
612: PetscReal *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd;
613: PetscReal *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk;
614: PetscReal **Q = df->Q, *f = df->f, *t = df->t;
615: PetscInt dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv;
617: /*** variables for the adaptive nonmonotone linesearch ***/
618: PetscInt L, llast;
619: PetscReal fr, fbest, fv, fc, fv0;
621: c = BMRM_INFTY;
623: DELTAsv = EPS_SV;
624: if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F;
625: else ProdDELTAsv = EPS_SV;
627: for (i = 0; i < dim; i++) tempv[i] = -x[i];
629: lam_ext = 0.0;
631: /* Project the initial solution */
632: projcount += project(dim, a, b, tempv, l, u, x, &lam_ext, df);
634: /* Compute gradient
635: g = Q*x + f; */
637: it = 0;
638: for (i = 0; i < dim; i++) {
639: if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i;
640: }
642: PetscMemzero(t, dim*sizeof(PetscReal));
643: for (i = 0; i < it; i++){
644: tempQ = Q[ipt[i]];
645: for (j = 0; j < dim; j++) t[j] += (tempQ[j]*x[ipt[i]]);
646: }
647: for (i = 0; i < dim; i++){
648: g[i] = t[i] + f[i];
649: }
652: /* y = -(x_{k} - g_{k}) */
653: for (i = 0; i < dim; i++){
654: y[i] = g[i] - x[i];
655: }
657: /* Project x_{k} - g_{k} */
658: projcount += project(dim, a, b, y, l, u, tempv, &lam_ext, df);
660: /* y = P(x_{k} - g_{k}) - x_{k} */
661: max = ALPHA_MIN;
662: for (i = 0; i < dim; i++){
663: y[i] = tempv[i] - x[i];
664: if (PetscAbsReal(y[i]) > max) max = PetscAbsReal(y[i]);
665: }
667: if (max < tol*1e-3){
668: return 0;
669: }
671: alpha = 1.0 / max;
673: /* fv0 = f(x_{0}). Recall t = Q x_{k} */
674: fv0 = 0.0;
675: for (i = 0; i < dim; i++) fv0 += x[i] * (0.5*t[i] + f[i]);
677: /*** adaptive nonmonotone linesearch ***/
678: L = 2;
679: fr = ALPHA_MAX;
680: fbest = fv0;
681: fc = fv0;
682: llast = 0;
683: akold = bkold = 0.0;
685: /*** Iterator begins ***/
686: for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) {
688: /* tempv = -(x_{k} - alpha*g_{k}) */
689: for (i = 0; i < dim; i++) tempv[i] = alpha*g[i] - x[i];
691: /* Project x_{k} - alpha*g_{k} */
692: projcount += project(dim, a, b, tempv, l, u, y, &lam_ext, df);
695: /* gd = \inner{d_{k}}{g_{k}}
696: d = P(x_{k} - alpha*g_{k}) - x_{k}
697: */
698: gd = 0.0;
699: for (i = 0; i < dim; i++){
700: d[i] = y[i] - x[i];
701: gd += d[i] * g[i];
702: }
704: /* Gradient computation */
706: /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */
708: it = it2 = 0;
709: for (i = 0; i < dim; i++){
710: if (PetscAbsReal(d[i]) > (ProdDELTAsv*1.0e-2)) ipt[it++] = i;
711: }
712: for (i = 0; i < dim; i++) {
713: if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i;
714: }
716: PetscMemzero(Qd, dim*sizeof(PetscReal));
717: /* compute Qd = Q*d */
718: if (it < it2){
719: for (i = 0; i < it; i++){
720: tempQ = Q[ipt[i]];
721: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]);
722: }
723: } else { /* compute Qd = Q*y-t */
724: for (i = 0; i < it2; i++){
725: tempQ = Q[ipt2[i]];
726: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]);
727: }
728: for (j = 0; j < dim; j++) Qd[j] -= t[j];
729: }
731: /* ak = inner{d_{k}}{d_{k}} */
732: ak = 0.0;
733: for (i = 0; i < dim; i++) ak += d[i] * d[i];
735: bk = 0.0;
736: for (i = 0; i < dim; i++) bk += d[i]*Qd[i];
738: if (bk > EPS*ak && gd < 0.0) lamnew = -gd/bk;
739: else lamnew = 1.0;
741: /* fv is computing f(x_{k} + d_{k}) */
742: fv = 0.0;
743: for (i = 0; i < dim; i++){
744: xplus[i] = x[i] + d[i];
745: tplus[i] = t[i] + Qd[i];
746: fv += xplus[i] * (0.5*tplus[i] + f[i]);
747: }
749: /* fr is fref */
750: if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)){
751: lscount++;
752: fv = 0.0;
753: for (i = 0; i < dim; i++){
754: xplus[i] = x[i] + lamnew*d[i];
755: tplus[i] = t[i] + lamnew*Qd[i];
756: fv += xplus[i] * (0.5*tplus[i] + f[i]);
757: }
758: }
760: for (i = 0; i < dim; i++){
761: sk[i] = xplus[i] - x[i];
762: yk[i] = tplus[i] - t[i];
763: x[i] = xplus[i];
764: t[i] = tplus[i];
765: g[i] = t[i] + f[i];
766: }
768: /* update the line search control parameters */
769: if (fv < fbest){
770: fbest = fv;
771: fc = fv;
772: llast = 0;
773: } else {
774: fc = (fc > fv ? fc : fv);
775: llast++;
776: if (llast == L){
777: fr = fc;
778: fc = fv;
779: llast = 0;
780: }
781: }
783: ak = bk = 0.0;
784: for (i = 0; i < dim; i++){
785: ak += sk[i] * sk[i];
786: bk += sk[i] * yk[i];
787: }
789: if (bk <= EPS*ak) alpha = ALPHA_MAX;
790: else {
791: if (bkold < EPS*akold) alpha = ak/bk;
792: else alpha = (akold+ak)/(bkold+bk);
794: if (alpha > ALPHA_MAX) alpha = ALPHA_MAX;
795: else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN;
796: }
798: akold = ak;
799: bkold = bk;
801: /*** stopping criterion based on KKT conditions ***/
802: /* at optimal, gradient of lagrangian w.r.t. x is zero */
804: bk = 0.0;
805: for (i = 0; i < dim; i++) bk += x[i] * x[i];
807: if (PetscSqrtReal(ak) < tol*10 * PetscSqrtReal(bk)){
808: it = 0;
809: luv = 0;
810: kktlam = 0.0;
811: for (i = 0; i < dim; i++){
812: /* x[i] is active hence lagrange multipliers for box constraints
813: are zero. The lagrange multiplier for ineq. const. is then
814: defined as below
815: */
816: if ((x[i] > DELTAsv) && (x[i] < c-DELTAsv)){
817: ipt[it++] = i;
818: kktlam = kktlam - a[i]*g[i];
819: } else uv[luv++] = i;
820: }
822: if (it == 0 && PetscSqrtReal(ak) < tol*0.5 * PetscSqrtReal(bk)) return 0;
823: else {
824: kktlam = kktlam/it;
825: info = 1;
826: for (i = 0; i < it; i++) {
827: if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) {
828: info = 0;
829: break;
830: }
831: }
832: if (info == 1) {
833: for (i = 0; i < luv; i++) {
834: if (x[uv[i]] <= DELTAsv){
835: /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may
836: not be zero. So, the gradient without beta is > 0
837: */
838: if (g[uv[i]] + kktlam*a[uv[i]] < -tol){
839: info = 0;
840: break;
841: }
842: } else {
843: /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may
844: not be zero. So, the gradient without eta is < 0
845: */
846: if (g[uv[i]] + kktlam*a[uv[i]] > tol) {
847: info = 0;
848: break;
849: }
850: }
851: }
852: }
854: if (info == 1) return 0;
855: }
856: }
857: }
858: return 0;
859: }