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: {
24: PetscNew(p);
25: VecDuplicate(X, &(*p)->V);
26: VecCopy(X, (*p)->V);
27: (*p)->next = NULL;
28: return(0);
29: }
31: static PetscErrorCode destroy_grad_list(Vec_Chain *head)
32: {
34: Vec_Chain *p = head->next, *q;
37: while (p) {
38: q = p->next;
39: VecDestroy(&p->V);
40: PetscFree(p);
41: p = q;
42: }
43: head->next = NULL;
44: return(0);
45: }
47: static PetscErrorCode TaoSolve_BMRM(Tao tao)
48: {
49: PetscErrorCode ierr;
50: TAO_DF df;
51: TAO_BMRM *bmrm = (TAO_BMRM*)tao->data;
53: /* Values and pointers to parts of the optimization problem */
54: PetscReal f = 0.0;
55: Vec W = tao->solution;
56: Vec G = tao->gradient;
57: PetscReal lambda;
58: PetscReal bt;
59: Vec_Chain grad_list, *tail_glist, *pgrad;
60: PetscInt i;
61: PetscMPIInt rank;
63: /* Used in converged criteria check */
64: PetscReal reg;
65: PetscReal jtwt = 0.0, max_jtwt, pre_epsilon, epsilon, jw, min_jw;
66: PetscReal innerSolverTol;
67: MPI_Comm comm;
70: PetscObjectGetComm((PetscObject)tao,&comm);
71: MPI_Comm_rank(comm, &rank);
72: lambda = bmrm->lambda;
74: /* Check Stopping Condition */
75: tao->step = 1.0;
76: max_jtwt = -BMRM_INFTY;
77: min_jw = BMRM_INFTY;
78: innerSolverTol = 1.0;
79: epsilon = 0.0;
81: if (rank == 0) {
82: init_df_solver(&df);
83: grad_list.next = NULL;
84: tail_glist = &grad_list;
85: }
87: df.tol = 1e-6;
88: tao->reason = TAO_CONTINUE_ITERATING;
90: /*-----------------Algorithm Begins------------------------*/
91: /* make the scatter */
92: VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w);
93: VecAssemblyBegin(bmrm->local_w);
94: VecAssemblyEnd(bmrm->local_w);
96: /* NOTE: In application pass the sub-gradient of Remp(W) */
97: TaoComputeObjectiveAndGradient(tao, W, &f, G);
98: TaoLogConvergenceHistory(tao,f,1.0,0.0,tao->ksp_its);
99: TaoMonitor(tao,tao->niter,f,1.0,0.0,tao->step);
100: (*tao->ops->convergencetest)(tao,tao->cnvP);
102: while (tao->reason == TAO_CONTINUE_ITERATING) {
103: /* Call general purpose update function */
104: if (tao->ops->update) {
105: (*tao->ops->update)(tao, tao->niter, tao->user_update);
106: }
108: /* compute bt = Remp(Wt-1) - <Wt-1, At> */
109: VecDot(W, G, &bt);
110: bt = f - bt;
112: /* First gather the gradient to the rank-0 node */
113: VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
114: VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
116: /* Bring up the inner solver */
117: if (rank == 0) {
118: ensure_df_space(tao->niter+1, &df);
119: make_grad_node(bmrm->local_w, &pgrad);
120: tail_glist->next = pgrad;
121: tail_glist = pgrad;
123: df.a[tao->niter] = 1.0;
124: df.f[tao->niter] = -bt;
125: df.u[tao->niter] = 1.0;
126: df.l[tao->niter] = 0.0;
128: /* set up the Q */
129: pgrad = grad_list.next;
130: for (i=0; i<=tao->niter; i++) {
131: if (!pgrad) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Assert that there are at least tao->niter+1 pgrad available");
132: VecDot(pgrad->V, bmrm->local_w, ®);
133: df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda;
134: pgrad = pgrad->next;
135: }
137: if (tao->niter > 0) {
138: df.x[tao->niter] = 0.0;
139: solve(&df);
140: } else
141: df.x[0] = 1.0;
143: /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */
144: jtwt = 0.0;
145: VecSet(bmrm->local_w, 0.0);
146: pgrad = grad_list.next;
147: for (i=0; i<=tao->niter; i++) {
148: jtwt -= df.x[i] * df.f[i];
149: VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V);
150: pgrad = pgrad->next;
151: }
153: VecNorm(bmrm->local_w, NORM_2, ®);
154: reg = 0.5*lambda*reg*reg;
155: jtwt -= reg;
156: } /* end if rank == 0 */
158: /* scatter the new W to all nodes */
159: VecScatterBegin(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);
160: VecScatterEnd(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);
162: TaoComputeObjectiveAndGradient(tao, W, &f, G);
164: MPI_Bcast(&jtwt,1,MPIU_REAL,0,comm);
165: MPI_Bcast(®,1,MPIU_REAL,0,comm);
167: jw = reg + f; /* J(w) = regularizer + Remp(w) */
168: if (jw < min_jw) min_jw = jw;
169: if (jtwt > max_jtwt) max_jtwt = jtwt;
171: pre_epsilon = epsilon;
172: epsilon = min_jw - jtwt;
174: if (rank == 0) {
175: if (innerSolverTol > epsilon) innerSolverTol = epsilon;
176: else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7;
178: /* if the annealing doesn't work well, lower the inner solver tolerance */
179: if (pre_epsilon < epsilon) innerSolverTol *= 0.2;
181: df.tol = innerSolverTol*0.5;
182: }
184: tao->niter++;
185: TaoLogConvergenceHistory(tao,min_jw,epsilon,0.0,tao->ksp_its);
186: TaoMonitor(tao,tao->niter,min_jw,epsilon,0.0,tao->step);
187: (*tao->ops->convergencetest)(tao,tao->cnvP);
188: }
190: /* free all the memory */
191: if (rank == 0) {
192: destroy_grad_list(&grad_list);
193: destroy_df_solver(&df);
194: }
196: VecDestroy(&bmrm->local_w);
197: VecScatterDestroy(&bmrm->scatter);
198: return(0);
199: }
201: /* ---------------------------------------------------------- */
203: static PetscErrorCode TaoSetup_BMRM(Tao tao)
204: {
209: /* Allocate some arrays */
210: if (!tao->gradient) {
211: VecDuplicate(tao->solution, &tao->gradient);
212: }
213: return(0);
214: }
216: /*------------------------------------------------------------*/
217: static PetscErrorCode TaoDestroy_BMRM(Tao tao)
218: {
222: PetscFree(tao->data);
223: return(0);
224: }
226: static PetscErrorCode TaoSetFromOptions_BMRM(PetscOptionItems *PetscOptionsObject,Tao tao)
227: {
229: TAO_BMRM* bmrm = (TAO_BMRM*)tao->data;
232: PetscOptionsHead(PetscOptionsObject,"BMRM for regularized risk minimization");
233: PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight","", 100,&bmrm->lambda,NULL);
234: PetscOptionsTail();
235: return(0);
236: }
238: /*------------------------------------------------------------*/
239: static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer)
240: {
241: PetscBool isascii;
245: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
246: if (isascii) {
247: PetscViewerASCIIPushTab(viewer);
248: PetscViewerASCIIPopTab(viewer);
249: }
250: return(0);
251: }
253: /*------------------------------------------------------------*/
254: /*MC
255: TAOBMRM - bundle method for regularized risk minimization
257: Options Database Keys:
258: . - tao_bmrm_lambda - regulariser weight
260: Level: beginner
261: M*/
263: PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao)
264: {
265: TAO_BMRM *bmrm;
269: tao->ops->setup = TaoSetup_BMRM;
270: tao->ops->solve = TaoSolve_BMRM;
271: tao->ops->view = TaoView_BMRM;
272: tao->ops->setfromoptions = TaoSetFromOptions_BMRM;
273: tao->ops->destroy = TaoDestroy_BMRM;
275: PetscNewLog(tao,&bmrm);
276: bmrm->lambda = 1.0;
277: tao->data = (void*)bmrm;
279: /* Override default settings (unless already changed) */
280: if (!tao->max_it_changed) tao->max_it = 2000;
281: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
282: if (!tao->gatol_changed) tao->gatol = 1.0e-12;
283: if (!tao->grtol_changed) tao->grtol = 1.0e-12;
285: return(0);
286: }
288: PetscErrorCode init_df_solver(TAO_DF *df)
289: {
290: PetscInt i, n = INCRE_DIM;
294: /* default values */
295: df->maxProjIter = 200;
296: df->maxPGMIter = 300000;
297: df->b = 1.0;
299: /* memory space required by Dai-Fletcher */
300: df->cur_num_cp = n;
301: PetscMalloc1(n, &df->f);
302: PetscMalloc1(n, &df->a);
303: PetscMalloc1(n, &df->l);
304: PetscMalloc1(n, &df->u);
305: PetscMalloc1(n, &df->x);
306: PetscMalloc1(n, &df->Q);
308: for (i = 0; i < n; i ++) {
309: PetscMalloc1(n, &df->Q[i]);
310: }
312: PetscMalloc1(n, &df->g);
313: PetscMalloc1(n, &df->y);
314: PetscMalloc1(n, &df->tempv);
315: PetscMalloc1(n, &df->d);
316: PetscMalloc1(n, &df->Qd);
317: PetscMalloc1(n, &df->t);
318: PetscMalloc1(n, &df->xplus);
319: PetscMalloc1(n, &df->tplus);
320: PetscMalloc1(n, &df->sk);
321: PetscMalloc1(n, &df->yk);
323: PetscMalloc1(n, &df->ipt);
324: PetscMalloc1(n, &df->ipt2);
325: PetscMalloc1(n, &df->uv);
326: return(0);
327: }
329: PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df)
330: {
332: PetscReal *tmp, **tmp_Q;
333: PetscInt i, n, old_n;
336: df->dim = dim;
337: if (dim <= df->cur_num_cp) return(0);
339: old_n = df->cur_num_cp;
340: df->cur_num_cp += INCRE_DIM;
341: n = df->cur_num_cp;
343: /* memory space required by dai-fletcher */
344: PetscMalloc1(n, &tmp);
345: PetscArraycpy(tmp, df->f, old_n);
346: PetscFree(df->f);
347: df->f = tmp;
349: PetscMalloc1(n, &tmp);
350: PetscArraycpy(tmp, df->a, old_n);
351: PetscFree(df->a);
352: df->a = tmp;
354: PetscMalloc1(n, &tmp);
355: PetscArraycpy(tmp, df->l, old_n);
356: PetscFree(df->l);
357: df->l = tmp;
359: PetscMalloc1(n, &tmp);
360: PetscArraycpy(tmp, df->u, old_n);
361: PetscFree(df->u);
362: df->u = tmp;
364: PetscMalloc1(n, &tmp);
365: PetscArraycpy(tmp, df->x, old_n);
366: PetscFree(df->x);
367: df->x = tmp;
369: PetscMalloc1(n, &tmp_Q);
370: for (i = 0; i < n; i ++) {
371: PetscMalloc1(n, &tmp_Q[i]);
372: if (i < old_n) {
373: PetscArraycpy(tmp_Q[i], df->Q[i], old_n);
374: PetscFree(df->Q[i]);
375: }
376: }
378: PetscFree(df->Q);
379: df->Q = tmp_Q;
381: PetscFree(df->g);
382: PetscMalloc1(n, &df->g);
384: PetscFree(df->y);
385: PetscMalloc1(n, &df->y);
387: PetscFree(df->tempv);
388: PetscMalloc1(n, &df->tempv);
390: PetscFree(df->d);
391: PetscMalloc1(n, &df->d);
393: PetscFree(df->Qd);
394: PetscMalloc1(n, &df->Qd);
396: PetscFree(df->t);
397: PetscMalloc1(n, &df->t);
399: PetscFree(df->xplus);
400: PetscMalloc1(n, &df->xplus);
402: PetscFree(df->tplus);
403: PetscMalloc1(n, &df->tplus);
405: PetscFree(df->sk);
406: PetscMalloc1(n, &df->sk);
408: PetscFree(df->yk);
409: PetscMalloc1(n, &df->yk);
411: PetscFree(df->ipt);
412: PetscMalloc1(n, &df->ipt);
414: PetscFree(df->ipt2);
415: PetscMalloc1(n, &df->ipt2);
417: PetscFree(df->uv);
418: PetscMalloc1(n, &df->uv);
419: return(0);
420: }
422: PetscErrorCode destroy_df_solver(TAO_DF *df)
423: {
425: PetscInt i;
428: PetscFree(df->f);
429: PetscFree(df->a);
430: PetscFree(df->l);
431: PetscFree(df->u);
432: PetscFree(df->x);
434: for (i = 0; i < df->cur_num_cp; i ++) {
435: PetscFree(df->Q[i]);
436: }
437: PetscFree(df->Q);
438: PetscFree(df->ipt);
439: PetscFree(df->ipt2);
440: PetscFree(df->uv);
441: PetscFree(df->g);
442: PetscFree(df->y);
443: PetscFree(df->tempv);
444: PetscFree(df->d);
445: PetscFree(df->Qd);
446: PetscFree(df->t);
447: PetscFree(df->xplus);
448: PetscFree(df->tplus);
449: PetscFree(df->sk);
450: PetscFree(df->yk);
451: return(0);
452: }
454: /* Piecewise linear monotone target function for the Dai-Fletcher projector */
455: PetscReal phi(PetscReal *x,PetscInt n,PetscReal lambda,PetscReal *a,PetscReal b,PetscReal *c,PetscReal *l,PetscReal *u)
456: {
457: PetscReal r = 0.0;
458: PetscInt i;
460: for (i = 0; i < n; i++) {
461: x[i] = -c[i] + lambda*a[i];
462: if (x[i] > u[i]) x[i] = u[i];
463: else if (x[i] < l[i]) x[i] = l[i];
464: r += a[i]*x[i];
465: }
466: return r - b;
467: }
469: /** Modified Dai-Fletcher QP projector solves the problem:
470: *
471: * minimise 0.5*x'*x - c'*x
472: * subj to a'*x = b
473: * l \leq x \leq u
474: *
475: * \param c The point to be projected onto feasible set
476: */
477: PetscInt project(PetscInt n,PetscReal *a,PetscReal b,PetscReal *c,PetscReal *l,PetscReal *u,PetscReal *x,PetscReal *lam_ext,TAO_DF *df)
478: {
479: PetscReal lambda, lambdal, lambdau, dlambda, lambda_new;
480: PetscReal r, rl, ru, s;
481: PetscInt innerIter;
482: PetscBool nonNegativeSlack = PETSC_FALSE;
485: *lam_ext = 0;
486: lambda = 0;
487: dlambda = 0.5;
488: innerIter = 1;
490: /* \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b)
491: *
492: * Optimality conditions for \phi:
493: *
494: * 1. lambda <= 0
495: * 2. r <= 0
496: * 3. r*lambda == 0
497: */
499: /* Bracketing Phase */
500: r = phi(x, n, lambda, a, b, c, l, u);
502: if (nonNegativeSlack) {
503: /* inequality constraint, i.e., with \xi >= 0 constraint */
504: if (r < TOL_R) return 0;
505: } else {
506: /* equality constraint ,i.e., without \xi >= 0 constraint */
507: if (PetscAbsReal(r) < TOL_R) return 0;
508: }
510: if (r < 0.0) {
511: lambdal = lambda;
512: rl = r;
513: lambda = lambda + dlambda;
514: r = phi(x, n, lambda, a, b, c, l, u);
515: while (r < 0.0 && dlambda < BMRM_INFTY) {
516: lambdal = lambda;
517: s = rl/r - 1.0;
518: if (s < 0.1) s = 0.1;
519: dlambda = dlambda + dlambda/s;
520: lambda = lambda + dlambda;
521: rl = r;
522: r = phi(x, n, lambda, a, b, c, l, u);
523: }
524: lambdau = lambda;
525: ru = r;
526: } else {
527: lambdau = lambda;
528: ru = r;
529: lambda = lambda - dlambda;
530: r = phi(x, n, lambda, a, b, c, l, u);
531: while (r > 0.0 && dlambda > -BMRM_INFTY) {
532: lambdau = lambda;
533: s = ru/r - 1.0;
534: if (s < 0.1) s = 0.1;
535: dlambda = dlambda + dlambda/s;
536: lambda = lambda - dlambda;
537: ru = r;
538: r = phi(x, n, lambda, a, b, c, l, u);
539: }
540: lambdal = lambda;
541: rl = r;
542: }
544: if (PetscAbsReal(dlambda) > BMRM_INFTY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"L2N2_DaiFletcherPGM detected Infeasible QP problem!");
546: if (ru == 0) {
547: return innerIter;
548: }
550: /* Secant Phase */
551: s = 1.0 - rl/ru;
552: dlambda = dlambda/s;
553: lambda = lambdau - dlambda;
554: r = phi(x, n, lambda, a, b, c, l, u);
556: while (PetscAbsReal(r) > TOL_R
557: && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda))
558: && innerIter < df->maxProjIter) {
559: innerIter++;
560: if (r > 0.0) {
561: if (s <= 2.0) {
562: lambdau = lambda;
563: ru = r;
564: s = 1.0 - rl/ru;
565: dlambda = (lambdau - lambdal) / s;
566: lambda = lambdau - dlambda;
567: } else {
568: s = ru/r-1.0;
569: if (s < 0.1) s = 0.1;
570: dlambda = (lambdau - lambda) / s;
571: lambda_new = 0.75*lambdal + 0.25*lambda;
572: if (lambda_new < (lambda - dlambda))
573: lambda_new = lambda - dlambda;
574: lambdau = lambda;
575: ru = r;
576: lambda = lambda_new;
577: s = (lambdau - lambdal) / (lambdau - lambda);
578: }
579: } else {
580: if (s >= 2.0) {
581: lambdal = lambda;
582: rl = r;
583: s = 1.0 - rl/ru;
584: dlambda = (lambdau - lambdal) / s;
585: lambda = lambdau - dlambda;
586: } else {
587: s = rl/r - 1.0;
588: if (s < 0.1) s = 0.1;
589: dlambda = (lambda-lambdal) / s;
590: lambda_new = 0.75*lambdau + 0.25*lambda;
591: if (lambda_new > (lambda + dlambda))
592: lambda_new = lambda + dlambda;
593: lambdal = lambda;
594: rl = r;
595: lambda = lambda_new;
596: s = (lambdau - lambdal) / (lambdau-lambda);
597: }
598: }
599: r = phi(x, n, lambda, a, b, c, l, u);
600: }
602: *lam_ext = lambda;
603: if (innerIter >= df->maxProjIter) {
604: PetscInfo(NULL,"WARNING: DaiFletcher max iterations\n");
605: }
606: return innerIter;
607: }
609: PetscErrorCode solve(TAO_DF *df)
610: {
612: PetscInt i, j, innerIter, it, it2, luv, info, lscount = 0;
613: PetscReal gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam=0.0, lam_ext;
614: PetscReal DELTAsv, ProdDELTAsv;
615: PetscReal c, *tempQ;
616: PetscReal *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol;
617: PetscReal *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd;
618: PetscReal *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk;
619: PetscReal **Q = df->Q, *f = df->f, *t = df->t;
620: PetscInt dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv;
622: /* variables for the adaptive nonmonotone linesearch */
623: PetscInt L, llast;
624: PetscReal fr, fbest, fv, fc, fv0;
626: c = BMRM_INFTY;
628: DELTAsv = EPS_SV;
629: if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F;
630: else ProdDELTAsv = EPS_SV;
632: for (i = 0; i < dim; i++) tempv[i] = -x[i];
634: lam_ext = 0.0;
636: /* Project the initial solution */
637: project(dim, a, b, tempv, l, u, x, &lam_ext, df);
639: /* Compute gradient
640: g = Q*x + f; */
642: it = 0;
643: for (i = 0; i < dim; i++) {
644: if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i;
645: }
647: PetscArrayzero(t, dim);
648: for (i = 0; i < it; i++) {
649: tempQ = Q[ipt[i]];
650: for (j = 0; j < dim; j++) t[j] += (tempQ[j]*x[ipt[i]]);
651: }
652: for (i = 0; i < dim; i++) {
653: g[i] = t[i] + f[i];
654: }
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: 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: project(dim, a, b, tempv, l, u, y, &lam_ext, df);
698: /* gd = \inner{d_{k}}{g_{k}}
699: d = P(x_{k} - alpha*g_{k}) - x_{k}
700: */
701: gd = 0.0;
702: for (i = 0; i < dim; i++) {
703: d[i] = y[i] - x[i];
704: gd += d[i] * g[i];
705: }
707: /* Gradient computation */
709: /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */
711: it = it2 = 0;
712: for (i = 0; i < dim; i++) {
713: if (PetscAbsReal(d[i]) > (ProdDELTAsv*1.0e-2)) ipt[it++] = i;
714: }
715: for (i = 0; i < dim; i++) {
716: if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i;
717: }
719: PetscArrayzero(Qd, dim);
720: /* compute Qd = Q*d */
721: if (it < it2) {
722: for (i = 0; i < it; i++) {
723: tempQ = Q[ipt[i]];
724: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]);
725: }
726: } else { /* compute Qd = Q*y-t */
727: for (i = 0; i < it2; i++) {
728: tempQ = Q[ipt2[i]];
729: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]);
730: }
731: for (j = 0; j < dim; j++) Qd[j] -= t[j];
732: }
734: /* ak = inner{d_{k}}{d_{k}} */
735: ak = 0.0;
736: for (i = 0; i < dim; i++) ak += d[i] * d[i];
738: bk = 0.0;
739: for (i = 0; i < dim; i++) bk += d[i]*Qd[i];
741: if (bk > EPS*ak && gd < 0.0) lamnew = -gd/bk;
742: else lamnew = 1.0;
744: /* fv is computing f(x_{k} + d_{k}) */
745: fv = 0.0;
746: for (i = 0; i < dim; i++) {
747: xplus[i] = x[i] + d[i];
748: tplus[i] = t[i] + Qd[i];
749: fv += xplus[i] * (0.5*tplus[i] + f[i]);
750: }
752: /* fr is fref */
753: if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)) {
754: lscount++;
755: fv = 0.0;
756: for (i = 0; i < dim; i++) {
757: xplus[i] = x[i] + lamnew*d[i];
758: tplus[i] = t[i] + lamnew*Qd[i];
759: fv += xplus[i] * (0.5*tplus[i] + f[i]);
760: }
761: }
763: for (i = 0; i < dim; i++) {
764: sk[i] = xplus[i] - x[i];
765: yk[i] = tplus[i] - t[i];
766: x[i] = xplus[i];
767: t[i] = tplus[i];
768: g[i] = t[i] + f[i];
769: }
771: /* update the line search control parameters */
772: if (fv < fbest) {
773: fbest = fv;
774: fc = fv;
775: llast = 0;
776: } else {
777: fc = (fc > fv ? fc : fv);
778: llast++;
779: if (llast == L) {
780: fr = fc;
781: fc = fv;
782: llast = 0;
783: }
784: }
786: ak = bk = 0.0;
787: for (i = 0; i < dim; i++) {
788: ak += sk[i] * sk[i];
789: bk += sk[i] * yk[i];
790: }
792: if (bk <= EPS*ak) alpha = ALPHA_MAX;
793: else {
794: if (bkold < EPS*akold) alpha = ak/bk;
795: else alpha = (akold+ak)/(bkold+bk);
797: if (alpha > ALPHA_MAX) alpha = ALPHA_MAX;
798: else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN;
799: }
801: akold = ak;
802: bkold = bk;
804: /* stopping criterion based on KKT conditions */
805: /* at optimal, gradient of lagrangian w.r.t. x is zero */
807: bk = 0.0;
808: for (i = 0; i < dim; i++) bk += x[i] * x[i];
810: if (PetscSqrtReal(ak) < tol*10 * PetscSqrtReal(bk)) {
811: it = 0;
812: luv = 0;
813: kktlam = 0.0;
814: for (i = 0; i < dim; i++) {
815: /* x[i] is active hence lagrange multipliers for box constraints
816: are zero. The lagrange multiplier for ineq. const. is then
817: defined as below
818: */
819: if ((x[i] > DELTAsv) && (x[i] < c-DELTAsv)) {
820: ipt[it++] = i;
821: kktlam = kktlam - a[i]*g[i];
822: } else uv[luv++] = i;
823: }
825: if (it == 0 && PetscSqrtReal(ak) < tol*0.5 * PetscSqrtReal(bk)) return 0;
826: else {
827: kktlam = kktlam/it;
828: info = 1;
829: for (i = 0; i < it; i++) {
830: if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) {
831: info = 0;
832: break;
833: }
834: }
835: if (info == 1) {
836: for (i = 0; i < luv; i++) {
837: if (x[uv[i]] <= DELTAsv) {
838: /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may
839: not be zero. So, the gradient without beta is > 0
840: */
841: if (g[uv[i]] + kktlam*a[uv[i]] < -tol) {
842: info = 0;
843: break;
844: }
845: } else {
846: /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may
847: not be zero. So, the gradient without eta is < 0
848: */
849: if (g[uv[i]] + kktlam*a[uv[i]] > tol) {
850: info = 0;
851: break;
852: }
853: }
854: }
855: }
857: if (info == 1) return 0;
858: }
859: }
860: }
861: return 0;
862: }