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: PetscNew(p);
22: VecDuplicate(X, &(*p)->V);
23: VecCopy(X, (*p)->V);
24: (*p)->next = NULL;
25: return 0;
26: }
28: static PetscErrorCode destroy_grad_list(Vec_Chain *head)
29: {
30: Vec_Chain *p = head->next, *q;
32: while (p) {
33: q = p->next;
34: VecDestroy(&p->V);
35: PetscFree(p);
36: p = q;
37: }
38: head->next = NULL;
39: return 0;
40: }
42: static PetscErrorCode TaoSolve_BMRM(Tao tao)
43: {
44: TAO_DF df;
45: TAO_BMRM *bmrm = (TAO_BMRM*)tao->data;
47: /* Values and pointers to parts of the optimization problem */
48: PetscReal f = 0.0;
49: Vec W = tao->solution;
50: Vec G = tao->gradient;
51: PetscReal lambda;
52: PetscReal bt;
53: Vec_Chain grad_list, *tail_glist, *pgrad;
54: PetscInt i;
55: PetscMPIInt rank;
57: /* Used in converged criteria check */
58: PetscReal reg;
59: PetscReal jtwt = 0.0, max_jtwt, pre_epsilon, epsilon, jw, min_jw;
60: PetscReal innerSolverTol;
61: MPI_Comm comm;
63: PetscObjectGetComm((PetscObject)tao,&comm);
64: MPI_Comm_rank(comm, &rank);
65: lambda = bmrm->lambda;
67: /* Check Stopping Condition */
68: tao->step = 1.0;
69: max_jtwt = -BMRM_INFTY;
70: min_jw = BMRM_INFTY;
71: innerSolverTol = 1.0;
72: epsilon = 0.0;
74: if (rank == 0) {
75: init_df_solver(&df);
76: grad_list.next = NULL;
77: tail_glist = &grad_list;
78: }
80: df.tol = 1e-6;
81: tao->reason = TAO_CONTINUE_ITERATING;
83: /*-----------------Algorithm Begins------------------------*/
84: /* make the scatter */
85: VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w);
86: VecAssemblyBegin(bmrm->local_w);
87: VecAssemblyEnd(bmrm->local_w);
89: /* NOTE: In application pass the sub-gradient of Remp(W) */
90: TaoComputeObjectiveAndGradient(tao, W, &f, G);
91: TaoLogConvergenceHistory(tao,f,1.0,0.0,tao->ksp_its);
92: TaoMonitor(tao,tao->niter,f,1.0,0.0,tao->step);
93: (*tao->ops->convergencetest)(tao,tao->cnvP);
95: while (tao->reason == TAO_CONTINUE_ITERATING) {
96: /* Call general purpose update function */
97: if (tao->ops->update) {
98: (*tao->ops->update)(tao, tao->niter, tao->user_update);
99: }
101: /* compute bt = Remp(Wt-1) - <Wt-1, At> */
102: VecDot(W, G, &bt);
103: bt = f - bt;
105: /* First gather the gradient to the rank-0 node */
106: VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
107: VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
109: /* Bring up the inner solver */
110: if (rank == 0) {
111: ensure_df_space(tao->niter+1, &df);
112: make_grad_node(bmrm->local_w, &pgrad);
113: tail_glist->next = pgrad;
114: tail_glist = pgrad;
116: df.a[tao->niter] = 1.0;
117: df.f[tao->niter] = -bt;
118: df.u[tao->niter] = 1.0;
119: df.l[tao->niter] = 0.0;
121: /* set up the Q */
122: pgrad = grad_list.next;
123: for (i=0; i<=tao->niter; i++) {
125: VecDot(pgrad->V, bmrm->local_w, ®);
126: df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda;
127: pgrad = pgrad->next;
128: }
130: if (tao->niter > 0) {
131: df.x[tao->niter] = 0.0;
132: solve(&df);
133: } else
134: df.x[0] = 1.0;
136: /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */
137: jtwt = 0.0;
138: 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: VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V);
143: pgrad = pgrad->next;
144: }
146: 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: VecScatterBegin(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);
153: VecScatterEnd(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);
155: TaoComputeObjectiveAndGradient(tao, W, &f, G);
157: MPI_Bcast(&jtwt,1,MPIU_REAL,0,comm);
158: 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: TaoLogConvergenceHistory(tao,min_jw,epsilon,0.0,tao->ksp_its);
179: TaoMonitor(tao,tao->niter,min_jw,epsilon,0.0,tao->step);
180: (*tao->ops->convergencetest)(tao,tao->cnvP);
181: }
183: /* free all the memory */
184: if (rank == 0) {
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: }
194: /* ---------------------------------------------------------- */
196: static PetscErrorCode TaoSetup_BMRM(Tao tao)
197: {
198: /* Allocate some arrays */
199: if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
200: return 0;
201: }
203: /*------------------------------------------------------------*/
204: static PetscErrorCode TaoDestroy_BMRM(Tao tao)
205: {
206: PetscFree(tao->data);
207: return 0;
208: }
210: static PetscErrorCode TaoSetFromOptions_BMRM(PetscOptionItems *PetscOptionsObject,Tao tao)
211: {
212: TAO_BMRM* bmrm = (TAO_BMRM*)tao->data;
214: PetscOptionsHead(PetscOptionsObject,"BMRM for regularized risk minimization");
215: PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight","", 100,&bmrm->lambda,NULL);
216: PetscOptionsTail();
217: return 0;
218: }
220: /*------------------------------------------------------------*/
221: static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer)
222: {
223: PetscBool isascii;
225: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
226: if (isascii) {
227: PetscViewerASCIIPushTab(viewer);
228: PetscViewerASCIIPopTab(viewer);
229: }
230: return 0;
231: }
233: /*------------------------------------------------------------*/
234: /*MC
235: TAOBMRM - bundle method for regularized risk minimization
237: Options Database Keys:
238: . - tao_bmrm_lambda - regulariser weight
240: Level: beginner
241: M*/
243: PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao)
244: {
245: TAO_BMRM *bmrm;
247: tao->ops->setup = TaoSetup_BMRM;
248: tao->ops->solve = TaoSolve_BMRM;
249: tao->ops->view = TaoView_BMRM;
250: tao->ops->setfromoptions = TaoSetFromOptions_BMRM;
251: tao->ops->destroy = TaoDestroy_BMRM;
253: PetscNewLog(tao,&bmrm);
254: bmrm->lambda = 1.0;
255: tao->data = (void*)bmrm;
257: /* Override default settings (unless already changed) */
258: if (!tao->max_it_changed) tao->max_it = 2000;
259: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
260: if (!tao->gatol_changed) tao->gatol = 1.0e-12;
261: if (!tao->grtol_changed) tao->grtol = 1.0e-12;
263: return 0;
264: }
266: PetscErrorCode init_df_solver(TAO_DF *df)
267: {
268: PetscInt i, n = INCRE_DIM;
270: /* default values */
271: df->maxProjIter = 200;
272: df->maxPGMIter = 300000;
273: df->b = 1.0;
275: /* memory space required by Dai-Fletcher */
276: df->cur_num_cp = n;
277: PetscMalloc1(n, &df->f);
278: PetscMalloc1(n, &df->a);
279: PetscMalloc1(n, &df->l);
280: PetscMalloc1(n, &df->u);
281: PetscMalloc1(n, &df->x);
282: PetscMalloc1(n, &df->Q);
284: for (i = 0; i < n; i ++) {
285: PetscMalloc1(n, &df->Q[i]);
286: }
288: PetscMalloc1(n, &df->g);
289: PetscMalloc1(n, &df->y);
290: PetscMalloc1(n, &df->tempv);
291: PetscMalloc1(n, &df->d);
292: PetscMalloc1(n, &df->Qd);
293: PetscMalloc1(n, &df->t);
294: PetscMalloc1(n, &df->xplus);
295: PetscMalloc1(n, &df->tplus);
296: PetscMalloc1(n, &df->sk);
297: PetscMalloc1(n, &df->yk);
299: PetscMalloc1(n, &df->ipt);
300: PetscMalloc1(n, &df->ipt2);
301: PetscMalloc1(n, &df->uv);
302: return 0;
303: }
305: PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df)
306: {
307: PetscReal *tmp, **tmp_Q;
308: PetscInt i, n, old_n;
310: df->dim = dim;
311: if (dim <= df->cur_num_cp) return 0;
313: old_n = df->cur_num_cp;
314: df->cur_num_cp += INCRE_DIM;
315: n = df->cur_num_cp;
317: /* memory space required by dai-fletcher */
318: PetscMalloc1(n, &tmp);
319: PetscArraycpy(tmp, df->f, old_n);
320: PetscFree(df->f);
321: df->f = tmp;
323: PetscMalloc1(n, &tmp);
324: PetscArraycpy(tmp, df->a, old_n);
325: PetscFree(df->a);
326: df->a = tmp;
328: PetscMalloc1(n, &tmp);
329: PetscArraycpy(tmp, df->l, old_n);
330: PetscFree(df->l);
331: df->l = tmp;
333: PetscMalloc1(n, &tmp);
334: PetscArraycpy(tmp, df->u, old_n);
335: PetscFree(df->u);
336: df->u = tmp;
338: PetscMalloc1(n, &tmp);
339: PetscArraycpy(tmp, df->x, old_n);
340: PetscFree(df->x);
341: df->x = tmp;
343: PetscMalloc1(n, &tmp_Q);
344: for (i = 0; i < n; i ++) {
345: PetscMalloc1(n, &tmp_Q[i]);
346: if (i < old_n) {
347: PetscArraycpy(tmp_Q[i], df->Q[i], old_n);
348: PetscFree(df->Q[i]);
349: }
350: }
352: PetscFree(df->Q);
353: df->Q = tmp_Q;
355: PetscFree(df->g);
356: PetscMalloc1(n, &df->g);
358: PetscFree(df->y);
359: PetscMalloc1(n, &df->y);
361: PetscFree(df->tempv);
362: PetscMalloc1(n, &df->tempv);
364: PetscFree(df->d);
365: PetscMalloc1(n, &df->d);
367: PetscFree(df->Qd);
368: PetscMalloc1(n, &df->Qd);
370: PetscFree(df->t);
371: PetscMalloc1(n, &df->t);
373: PetscFree(df->xplus);
374: PetscMalloc1(n, &df->xplus);
376: PetscFree(df->tplus);
377: PetscMalloc1(n, &df->tplus);
379: PetscFree(df->sk);
380: PetscMalloc1(n, &df->sk);
382: PetscFree(df->yk);
383: PetscMalloc1(n, &df->yk);
385: PetscFree(df->ipt);
386: PetscMalloc1(n, &df->ipt);
388: PetscFree(df->ipt2);
389: PetscMalloc1(n, &df->ipt2);
391: PetscFree(df->uv);
392: PetscMalloc1(n, &df->uv);
393: return 0;
394: }
396: PetscErrorCode destroy_df_solver(TAO_DF *df)
397: {
398: PetscInt i;
400: PetscFree(df->f);
401: PetscFree(df->a);
402: PetscFree(df->l);
403: PetscFree(df->u);
404: PetscFree(df->x);
406: for (i = 0; i < df->cur_num_cp; i ++) {
407: PetscFree(df->Q[i]);
408: }
409: PetscFree(df->Q);
410: PetscFree(df->ipt);
411: PetscFree(df->ipt2);
412: PetscFree(df->uv);
413: PetscFree(df->g);
414: PetscFree(df->y);
415: PetscFree(df->tempv);
416: PetscFree(df->d);
417: PetscFree(df->Qd);
418: PetscFree(df->t);
419: PetscFree(df->xplus);
420: PetscFree(df->tplus);
421: PetscFree(df->sk);
422: PetscFree(df->yk);
423: return 0;
424: }
426: /* Piecewise linear monotone target function for the Dai-Fletcher projector */
427: PetscReal phi(PetscReal *x,PetscInt n,PetscReal lambda,PetscReal *a,PetscReal b,PetscReal *c,PetscReal *l,PetscReal *u)
428: {
429: PetscReal r = 0.0;
430: PetscInt i;
432: for (i = 0; i < n; i++) {
433: x[i] = -c[i] + lambda*a[i];
434: if (x[i] > u[i]) x[i] = u[i];
435: else if (x[i] < l[i]) x[i] = l[i];
436: r += a[i]*x[i];
437: }
438: return r - b;
439: }
441: /** Modified Dai-Fletcher QP projector solves the problem:
442: *
443: * minimise 0.5*x'*x - c'*x
444: * subj to a'*x = b
445: * l \leq x \leq u
446: *
447: * \param c The point to be projected onto feasible set
448: */
449: PetscInt project(PetscInt n,PetscReal *a,PetscReal b,PetscReal *c,PetscReal *l,PetscReal *u,PetscReal *x,PetscReal *lam_ext,TAO_DF *df)
450: {
451: PetscReal lambda, lambdal, lambdau, dlambda, lambda_new;
452: PetscReal r, rl, ru, s;
453: PetscInt innerIter;
454: PetscBool nonNegativeSlack = PETSC_FALSE;
456: *lam_ext = 0;
457: lambda = 0;
458: dlambda = 0.5;
459: innerIter = 1;
461: /* \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b)
462: *
463: * Optimality conditions for \phi:
464: *
465: * 1. lambda <= 0
466: * 2. r <= 0
467: * 3. r*lambda == 0
468: */
470: /* Bracketing Phase */
471: r = phi(x, n, lambda, a, b, c, l, u);
473: if (nonNegativeSlack) {
474: /* inequality constraint, i.e., with \xi >= 0 constraint */
475: if (r < TOL_R) return 0;
476: } else {
477: /* equality constraint ,i.e., without \xi >= 0 constraint */
478: if (PetscAbsReal(r) < TOL_R) return 0;
479: }
481: if (r < 0.0) {
482: lambdal = lambda;
483: rl = r;
484: lambda = lambda + dlambda;
485: r = phi(x, n, lambda, a, b, c, l, u);
486: while (r < 0.0 && dlambda < BMRM_INFTY) {
487: lambdal = lambda;
488: s = rl/r - 1.0;
489: if (s < 0.1) s = 0.1;
490: dlambda = dlambda + dlambda/s;
491: lambda = lambda + dlambda;
492: rl = r;
493: r = phi(x, n, lambda, a, b, c, l, u);
494: }
495: lambdau = lambda;
496: ru = r;
497: } else {
498: lambdau = lambda;
499: ru = r;
500: lambda = lambda - dlambda;
501: r = phi(x, n, lambda, a, b, c, l, u);
502: while (r > 0.0 && dlambda > -BMRM_INFTY) {
503: lambdau = lambda;
504: s = ru/r - 1.0;
505: if (s < 0.1) s = 0.1;
506: dlambda = dlambda + dlambda/s;
507: lambda = lambda - dlambda;
508: ru = r;
509: r = phi(x, n, lambda, a, b, c, l, u);
510: }
511: lambdal = lambda;
512: rl = r;
513: }
517: if (ru == 0) {
518: return innerIter;
519: }
521: /* Secant Phase */
522: s = 1.0 - rl/ru;
523: dlambda = dlambda/s;
524: lambda = lambdau - dlambda;
525: r = phi(x, n, lambda, a, b, c, l, u);
527: while (PetscAbsReal(r) > TOL_R
528: && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda))
529: && 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))
544: lambda_new = lambda - dlambda;
545: lambdau = lambda;
546: ru = r;
547: lambda = lambda_new;
548: s = (lambdau - lambdal) / (lambdau - lambda);
549: }
550: } else {
551: if (s >= 2.0) {
552: lambdal = lambda;
553: rl = r;
554: s = 1.0 - rl/ru;
555: dlambda = (lambdau - lambdal) / s;
556: lambda = lambdau - dlambda;
557: } else {
558: s = rl/r - 1.0;
559: if (s < 0.1) s = 0.1;
560: dlambda = (lambda-lambdal) / s;
561: lambda_new = 0.75*lambdau + 0.25*lambda;
562: if (lambda_new > (lambda + dlambda))
563: lambda_new = lambda + dlambda;
564: lambdal = lambda;
565: rl = r;
566: lambda = lambda_new;
567: s = (lambdau - lambdal) / (lambdau-lambda);
568: }
569: }
570: r = phi(x, n, lambda, a, b, c, l, u);
571: }
573: *lam_ext = lambda;
574: if (innerIter >= df->maxProjIter) {
575: PetscInfo(NULL,"WARNING: DaiFletcher max iterations\n");
576: }
577: return innerIter;
578: }
580: PetscErrorCode solve(TAO_DF *df)
581: {
582: PetscInt i, j, innerIter, it, it2, luv, info, lscount = 0;
583: PetscReal gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam=0.0, lam_ext;
584: PetscReal DELTAsv, ProdDELTAsv;
585: PetscReal c, *tempQ;
586: PetscReal *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol;
587: PetscReal *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd;
588: PetscReal *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk;
589: PetscReal **Q = df->Q, *f = df->f, *t = df->t;
590: PetscInt dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv;
592: /* variables for the adaptive nonmonotone linesearch */
593: PetscInt L, llast;
594: PetscReal fr, fbest, fv, fc, fv0;
596: c = BMRM_INFTY;
598: DELTAsv = EPS_SV;
599: if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F;
600: else ProdDELTAsv = EPS_SV;
602: for (i = 0; i < dim; i++) tempv[i] = -x[i];
604: lam_ext = 0.0;
606: /* Project the initial solution */
607: project(dim, a, b, tempv, l, u, x, &lam_ext, df);
609: /* Compute gradient
610: g = Q*x + f; */
612: it = 0;
613: for (i = 0; i < dim; i++) {
614: if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i;
615: }
617: PetscArrayzero(t, dim);
618: for (i = 0; i < it; i++) {
619: tempQ = Q[ipt[i]];
620: for (j = 0; j < dim; j++) t[j] += (tempQ[j]*x[ipt[i]]);
621: }
622: for (i = 0; i < dim; i++) {
623: g[i] = t[i] + f[i];
624: }
626: /* y = -(x_{k} - g_{k}) */
627: for (i = 0; i < dim; i++) {
628: y[i] = g[i] - x[i];
629: }
631: /* Project x_{k} - g_{k} */
632: project(dim, a, b, y, l, u, tempv, &lam_ext, df);
634: /* y = P(x_{k} - g_{k}) - x_{k} */
635: max = ALPHA_MIN;
636: for (i = 0; i < dim; i++) {
637: y[i] = tempv[i] - x[i];
638: if (PetscAbsReal(y[i]) > max) max = PetscAbsReal(y[i]);
639: }
641: if (max < tol*1e-3) {
642: return 0;
643: }
645: alpha = 1.0 / max;
647: /* fv0 = f(x_{0}). Recall t = Q x_{k} */
648: fv0 = 0.0;
649: for (i = 0; i < dim; i++) fv0 += x[i] * (0.5*t[i] + f[i]);
651: /* adaptive nonmonotone linesearch */
652: L = 2;
653: fr = ALPHA_MAX;
654: fbest = fv0;
655: fc = fv0;
656: llast = 0;
657: akold = bkold = 0.0;
659: /* Iterator begins */
660: for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) {
662: /* tempv = -(x_{k} - alpha*g_{k}) */
663: for (i = 0; i < dim; i++) tempv[i] = alpha*g[i] - x[i];
665: /* Project x_{k} - alpha*g_{k} */
666: project(dim, a, b, tempv, l, u, y, &lam_ext, df);
668: /* gd = \inner{d_{k}}{g_{k}}
669: d = P(x_{k} - alpha*g_{k}) - x_{k}
670: */
671: gd = 0.0;
672: for (i = 0; i < dim; i++) {
673: d[i] = y[i] - x[i];
674: gd += d[i] * g[i];
675: }
677: /* Gradient computation */
679: /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */
681: it = it2 = 0;
682: for (i = 0; i < dim; i++) {
683: if (PetscAbsReal(d[i]) > (ProdDELTAsv*1.0e-2)) ipt[it++] = i;
684: }
685: for (i = 0; i < dim; i++) {
686: if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i;
687: }
689: PetscArrayzero(Qd, dim);
690: /* compute Qd = Q*d */
691: if (it < it2) {
692: for (i = 0; i < it; i++) {
693: tempQ = Q[ipt[i]];
694: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]);
695: }
696: } else { /* compute Qd = Q*y-t */
697: for (i = 0; i < it2; i++) {
698: tempQ = Q[ipt2[i]];
699: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]);
700: }
701: for (j = 0; j < dim; j++) Qd[j] -= t[j];
702: }
704: /* ak = inner{d_{k}}{d_{k}} */
705: ak = 0.0;
706: for (i = 0; i < dim; i++) ak += d[i] * d[i];
708: bk = 0.0;
709: for (i = 0; i < dim; i++) bk += d[i]*Qd[i];
711: if (bk > EPS*ak && gd < 0.0) lamnew = -gd/bk;
712: else lamnew = 1.0;
714: /* fv is computing f(x_{k} + d_{k}) */
715: fv = 0.0;
716: for (i = 0; i < dim; i++) {
717: xplus[i] = x[i] + d[i];
718: tplus[i] = t[i] + Qd[i];
719: fv += xplus[i] * (0.5*tplus[i] + f[i]);
720: }
722: /* fr is fref */
723: if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)) {
724: lscount++;
725: fv = 0.0;
726: for (i = 0; i < dim; i++) {
727: xplus[i] = x[i] + lamnew*d[i];
728: tplus[i] = t[i] + lamnew*Qd[i];
729: fv += xplus[i] * (0.5*tplus[i] + f[i]);
730: }
731: }
733: for (i = 0; i < dim; i++) {
734: sk[i] = xplus[i] - x[i];
735: yk[i] = tplus[i] - t[i];
736: x[i] = xplus[i];
737: t[i] = tplus[i];
738: g[i] = t[i] + f[i];
739: }
741: /* update the line search control parameters */
742: if (fv < fbest) {
743: fbest = fv;
744: fc = fv;
745: llast = 0;
746: } else {
747: fc = (fc > fv ? fc : fv);
748: llast++;
749: if (llast == L) {
750: fr = fc;
751: fc = fv;
752: llast = 0;
753: }
754: }
756: ak = bk = 0.0;
757: for (i = 0; i < dim; i++) {
758: ak += sk[i] * sk[i];
759: bk += sk[i] * yk[i];
760: }
762: if (bk <= EPS*ak) alpha = ALPHA_MAX;
763: else {
764: if (bkold < EPS*akold) alpha = ak/bk;
765: else alpha = (akold+ak)/(bkold+bk);
767: if (alpha > ALPHA_MAX) alpha = ALPHA_MAX;
768: else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN;
769: }
771: akold = ak;
772: bkold = bk;
774: /* stopping criterion based on KKT conditions */
775: /* at optimal, gradient of lagrangian w.r.t. x is zero */
777: bk = 0.0;
778: for (i = 0; i < dim; i++) bk += x[i] * x[i];
780: if (PetscSqrtReal(ak) < tol*10 * PetscSqrtReal(bk)) {
781: it = 0;
782: luv = 0;
783: kktlam = 0.0;
784: for (i = 0; i < dim; i++) {
785: /* x[i] is active hence lagrange multipliers for box constraints
786: are zero. The lagrange multiplier for ineq. const. is then
787: defined as below
788: */
789: if ((x[i] > DELTAsv) && (x[i] < c-DELTAsv)) {
790: ipt[it++] = i;
791: kktlam = kktlam - a[i]*g[i];
792: } else uv[luv++] = i;
793: }
795: if (it == 0 && PetscSqrtReal(ak) < tol*0.5 * PetscSqrtReal(bk)) return 0;
796: else {
797: kktlam = kktlam/it;
798: info = 1;
799: for (i = 0; i < it; i++) {
800: if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) {
801: info = 0;
802: break;
803: }
804: }
805: if (info == 1) {
806: for (i = 0; i < luv; i++) {
807: if (x[uv[i]] <= DELTAsv) {
808: /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may
809: not be zero. So, the gradient without beta is > 0
810: */
811: if (g[uv[i]] + kktlam*a[uv[i]] < -tol) {
812: info = 0;
813: break;
814: }
815: } else {
816: /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may
817: not be zero. So, the gradient without eta is < 0
818: */
819: if (g[uv[i]] + kktlam*a[uv[i]] > tol) {
820: info = 0;
821: break;
822: }
823: }
824: }
825: }
827: if (info == 1) return 0;
828: }
829: }
830: }
831: return 0;
832: }