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: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
95: while (tao->reason == TAO_CONTINUE_ITERATING) {
96: /* Call general purpose update function */
97: PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);
99: /* compute bt = Remp(Wt-1) - <Wt-1, At> */
100: VecDot(W, G, &bt);
101: bt = f - bt;
103: /* First gather the gradient to the rank-0 node */
104: VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
105: VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
107: /* Bring up the inner solver */
108: if (rank == 0) {
109: ensure_df_space(tao->niter + 1, &df);
110: make_grad_node(bmrm->local_w, &pgrad);
111: tail_glist->next = pgrad;
112: tail_glist = pgrad;
114: df.a[tao->niter] = 1.0;
115: df.f[tao->niter] = -bt;
116: df.u[tao->niter] = 1.0;
117: df.l[tao->niter] = 0.0;
119: /* set up the Q */
120: pgrad = grad_list.next;
121: for (i = 0; i <= tao->niter; i++) {
123: VecDot(pgrad->V, bmrm->local_w, ®);
124: df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda;
125: pgrad = pgrad->next;
126: }
128: if (tao->niter > 0) {
129: df.x[tao->niter] = 0.0;
130: solve(&df);
131: } else df.x[0] = 1.0;
133: /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */
134: jtwt = 0.0;
135: VecSet(bmrm->local_w, 0.0);
136: pgrad = grad_list.next;
137: for (i = 0; i <= tao->niter; i++) {
138: jtwt -= df.x[i] * df.f[i];
139: VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V);
140: pgrad = pgrad->next;
141: }
143: VecNorm(bmrm->local_w, NORM_2, ®);
144: reg = 0.5 * lambda * reg * reg;
145: jtwt -= reg;
146: } /* end if rank == 0 */
148: /* scatter the new W to all nodes */
149: VecScatterBegin(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE);
150: VecScatterEnd(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE);
152: TaoComputeObjectiveAndGradient(tao, W, &f, G);
154: MPI_Bcast(&jtwt, 1, MPIU_REAL, 0, comm);
155: MPI_Bcast(®, 1, MPIU_REAL, 0, comm);
157: jw = reg + f; /* J(w) = regularizer + Remp(w) */
158: if (jw < min_jw) min_jw = jw;
159: if (jtwt > max_jtwt) max_jtwt = jtwt;
161: pre_epsilon = epsilon;
162: epsilon = min_jw - jtwt;
164: if (rank == 0) {
165: if (innerSolverTol > epsilon) innerSolverTol = epsilon;
166: else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7;
168: /* if the annealing doesn't work well, lower the inner solver tolerance */
169: if (pre_epsilon < epsilon) innerSolverTol *= 0.2;
171: df.tol = innerSolverTol * 0.5;
172: }
174: tao->niter++;
175: TaoLogConvergenceHistory(tao, min_jw, epsilon, 0.0, tao->ksp_its);
176: TaoMonitor(tao, tao->niter, min_jw, epsilon, 0.0, tao->step);
177: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
178: }
180: /* free all the memory */
181: if (rank == 0) {
182: destroy_grad_list(&grad_list);
183: destroy_df_solver(&df);
184: }
186: VecDestroy(&bmrm->local_w);
187: VecScatterDestroy(&bmrm->scatter);
188: return 0;
189: }
191: /* ---------------------------------------------------------- */
193: static PetscErrorCode TaoSetup_BMRM(Tao tao)
194: {
195: /* Allocate some arrays */
196: if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
197: return 0;
198: }
200: /*------------------------------------------------------------*/
201: static PetscErrorCode TaoDestroy_BMRM(Tao tao)
202: {
203: PetscFree(tao->data);
204: return 0;
205: }
207: static PetscErrorCode TaoSetFromOptions_BMRM(Tao tao, PetscOptionItems *PetscOptionsObject)
208: {
209: TAO_BMRM *bmrm = (TAO_BMRM *)tao->data;
211: PetscOptionsHeadBegin(PetscOptionsObject, "BMRM for regularized risk minimization");
212: PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight", "", 100, &bmrm->lambda, NULL);
213: PetscOptionsHeadEnd();
214: return 0;
215: }
217: /*------------------------------------------------------------*/
218: static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer)
219: {
220: PetscBool isascii;
222: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
223: if (isascii) {
224: PetscViewerASCIIPushTab(viewer);
225: PetscViewerASCIIPopTab(viewer);
226: }
227: return 0;
228: }
230: /*------------------------------------------------------------*/
231: /*MC
232: TAOBMRM - bundle method for regularized risk minimization
234: Options Database Keys:
235: . - tao_bmrm_lambda - regulariser weight
237: Level: beginner
238: M*/
240: PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao)
241: {
242: TAO_BMRM *bmrm;
244: tao->ops->setup = TaoSetup_BMRM;
245: tao->ops->solve = TaoSolve_BMRM;
246: tao->ops->view = TaoView_BMRM;
247: tao->ops->setfromoptions = TaoSetFromOptions_BMRM;
248: tao->ops->destroy = TaoDestroy_BMRM;
250: PetscNew(&bmrm);
251: bmrm->lambda = 1.0;
252: tao->data = (void *)bmrm;
254: /* Override default settings (unless already changed) */
255: if (!tao->max_it_changed) tao->max_it = 2000;
256: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
257: if (!tao->gatol_changed) tao->gatol = 1.0e-12;
258: if (!tao->grtol_changed) tao->grtol = 1.0e-12;
260: return 0;
261: }
263: PetscErrorCode init_df_solver(TAO_DF *df)
264: {
265: PetscInt i, n = INCRE_DIM;
267: /* default values */
268: df->maxProjIter = 200;
269: df->maxPGMIter = 300000;
270: df->b = 1.0;
272: /* memory space required by Dai-Fletcher */
273: df->cur_num_cp = n;
274: PetscMalloc1(n, &df->f);
275: PetscMalloc1(n, &df->a);
276: PetscMalloc1(n, &df->l);
277: PetscMalloc1(n, &df->u);
278: PetscMalloc1(n, &df->x);
279: PetscMalloc1(n, &df->Q);
281: for (i = 0; i < n; i++) PetscMalloc1(n, &df->Q[i]);
283: PetscMalloc1(n, &df->g);
284: PetscMalloc1(n, &df->y);
285: PetscMalloc1(n, &df->tempv);
286: PetscMalloc1(n, &df->d);
287: PetscMalloc1(n, &df->Qd);
288: PetscMalloc1(n, &df->t);
289: PetscMalloc1(n, &df->xplus);
290: PetscMalloc1(n, &df->tplus);
291: PetscMalloc1(n, &df->sk);
292: PetscMalloc1(n, &df->yk);
294: PetscMalloc1(n, &df->ipt);
295: PetscMalloc1(n, &df->ipt2);
296: PetscMalloc1(n, &df->uv);
297: return 0;
298: }
300: PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df)
301: {
302: PetscReal *tmp, **tmp_Q;
303: PetscInt i, n, old_n;
305: df->dim = dim;
306: if (dim <= df->cur_num_cp) return 0;
308: old_n = df->cur_num_cp;
309: df->cur_num_cp += INCRE_DIM;
310: n = df->cur_num_cp;
312: /* memory space required by dai-fletcher */
313: PetscMalloc1(n, &tmp);
314: PetscArraycpy(tmp, df->f, old_n);
315: PetscFree(df->f);
316: df->f = tmp;
318: PetscMalloc1(n, &tmp);
319: PetscArraycpy(tmp, df->a, old_n);
320: PetscFree(df->a);
321: df->a = tmp;
323: PetscMalloc1(n, &tmp);
324: PetscArraycpy(tmp, df->l, old_n);
325: PetscFree(df->l);
326: df->l = tmp;
328: PetscMalloc1(n, &tmp);
329: PetscArraycpy(tmp, df->u, old_n);
330: PetscFree(df->u);
331: df->u = tmp;
333: PetscMalloc1(n, &tmp);
334: PetscArraycpy(tmp, df->x, old_n);
335: PetscFree(df->x);
336: df->x = tmp;
338: PetscMalloc1(n, &tmp_Q);
339: for (i = 0; i < n; i++) {
340: PetscMalloc1(n, &tmp_Q[i]);
341: if (i < old_n) {
342: PetscArraycpy(tmp_Q[i], df->Q[i], old_n);
343: PetscFree(df->Q[i]);
344: }
345: }
347: PetscFree(df->Q);
348: df->Q = tmp_Q;
350: PetscFree(df->g);
351: PetscMalloc1(n, &df->g);
353: PetscFree(df->y);
354: PetscMalloc1(n, &df->y);
356: PetscFree(df->tempv);
357: PetscMalloc1(n, &df->tempv);
359: PetscFree(df->d);
360: PetscMalloc1(n, &df->d);
362: PetscFree(df->Qd);
363: PetscMalloc1(n, &df->Qd);
365: PetscFree(df->t);
366: PetscMalloc1(n, &df->t);
368: PetscFree(df->xplus);
369: PetscMalloc1(n, &df->xplus);
371: PetscFree(df->tplus);
372: PetscMalloc1(n, &df->tplus);
374: PetscFree(df->sk);
375: PetscMalloc1(n, &df->sk);
377: PetscFree(df->yk);
378: PetscMalloc1(n, &df->yk);
380: PetscFree(df->ipt);
381: PetscMalloc1(n, &df->ipt);
383: PetscFree(df->ipt2);
384: PetscMalloc1(n, &df->ipt2);
386: PetscFree(df->uv);
387: PetscMalloc1(n, &df->uv);
388: return 0;
389: }
391: PetscErrorCode destroy_df_solver(TAO_DF *df)
392: {
393: PetscInt i;
395: PetscFree(df->f);
396: PetscFree(df->a);
397: PetscFree(df->l);
398: PetscFree(df->u);
399: PetscFree(df->x);
401: for (i = 0; i < df->cur_num_cp; i++) PetscFree(df->Q[i]);
402: PetscFree(df->Q);
403: PetscFree(df->ipt);
404: PetscFree(df->ipt2);
405: PetscFree(df->uv);
406: PetscFree(df->g);
407: PetscFree(df->y);
408: PetscFree(df->tempv);
409: PetscFree(df->d);
410: PetscFree(df->Qd);
411: PetscFree(df->t);
412: PetscFree(df->xplus);
413: PetscFree(df->tplus);
414: PetscFree(df->sk);
415: PetscFree(df->yk);
416: return 0;
417: }
419: /* Piecewise linear monotone target function for the Dai-Fletcher projector */
420: PetscReal phi(PetscReal *x, PetscInt n, PetscReal lambda, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u)
421: {
422: PetscReal r = 0.0;
423: PetscInt i;
425: for (i = 0; i < n; i++) {
426: x[i] = -c[i] + lambda * a[i];
427: if (x[i] > u[i]) x[i] = u[i];
428: else if (x[i] < l[i]) x[i] = l[i];
429: r += a[i] * x[i];
430: }
431: return r - b;
432: }
434: /** Modified Dai-Fletcher QP projector solves the problem:
435: *
436: * minimise 0.5*x'*x - c'*x
437: * subj to a'*x = b
438: * l \leq x \leq u
439: *
440: * \param c The point to be projected onto feasible set
441: */
442: PetscInt project(PetscInt n, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u, PetscReal *x, PetscReal *lam_ext, TAO_DF *df)
443: {
444: PetscReal lambda, lambdal, lambdau, dlambda, lambda_new;
445: PetscReal r, rl, ru, s;
446: PetscInt innerIter;
447: PetscBool nonNegativeSlack = PETSC_FALSE;
449: *lam_ext = 0;
450: lambda = 0;
451: dlambda = 0.5;
452: innerIter = 1;
454: /* \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b)
455: *
456: * Optimality conditions for \phi:
457: *
458: * 1. lambda <= 0
459: * 2. r <= 0
460: * 3. r*lambda == 0
461: */
463: /* Bracketing Phase */
464: r = phi(x, n, lambda, a, b, c, l, u);
466: if (nonNegativeSlack) {
467: /* inequality constraint, i.e., with \xi >= 0 constraint */
468: if (r < TOL_R) return 0;
469: } else {
470: /* equality constraint ,i.e., without \xi >= 0 constraint */
471: if (PetscAbsReal(r) < TOL_R) return 0;
472: }
474: if (r < 0.0) {
475: lambdal = lambda;
476: rl = r;
477: lambda = lambda + dlambda;
478: r = phi(x, n, lambda, a, b, c, l, u);
479: while (r < 0.0 && dlambda < BMRM_INFTY) {
480: lambdal = lambda;
481: s = rl / r - 1.0;
482: if (s < 0.1) s = 0.1;
483: dlambda = dlambda + dlambda / s;
484: lambda = lambda + dlambda;
485: rl = r;
486: r = phi(x, n, lambda, a, b, c, l, u);
487: }
488: lambdau = lambda;
489: ru = r;
490: } else {
491: lambdau = lambda;
492: ru = r;
493: lambda = lambda - dlambda;
494: r = phi(x, n, lambda, a, b, c, l, u);
495: while (r > 0.0 && dlambda > -BMRM_INFTY) {
496: lambdau = lambda;
497: s = ru / r - 1.0;
498: if (s < 0.1) s = 0.1;
499: dlambda = dlambda + dlambda / s;
500: lambda = lambda - dlambda;
501: ru = r;
502: r = phi(x, n, lambda, a, b, c, l, u);
503: }
504: lambdal = lambda;
505: rl = r;
506: }
510: if (ru == 0) return innerIter;
512: /* Secant Phase */
513: s = 1.0 - rl / ru;
514: dlambda = dlambda / s;
515: lambda = lambdau - dlambda;
516: r = phi(x, n, lambda, a, b, c, l, u);
518: while (PetscAbsReal(r) > TOL_R && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda)) && innerIter < df->maxProjIter) {
519: innerIter++;
520: if (r > 0.0) {
521: if (s <= 2.0) {
522: lambdau = lambda;
523: ru = r;
524: s = 1.0 - rl / ru;
525: dlambda = (lambdau - lambdal) / s;
526: lambda = lambdau - dlambda;
527: } else {
528: s = ru / r - 1.0;
529: if (s < 0.1) s = 0.1;
530: dlambda = (lambdau - lambda) / s;
531: lambda_new = 0.75 * lambdal + 0.25 * lambda;
532: if (lambda_new < (lambda - dlambda)) lambda_new = lambda - dlambda;
533: lambdau = lambda;
534: ru = r;
535: lambda = lambda_new;
536: s = (lambdau - lambdal) / (lambdau - lambda);
537: }
538: } else {
539: if (s >= 2.0) {
540: lambdal = lambda;
541: rl = r;
542: s = 1.0 - rl / ru;
543: dlambda = (lambdau - lambdal) / s;
544: lambda = lambdau - dlambda;
545: } else {
546: s = rl / r - 1.0;
547: if (s < 0.1) s = 0.1;
548: dlambda = (lambda - lambdal) / s;
549: lambda_new = 0.75 * lambdau + 0.25 * lambda;
550: if (lambda_new > (lambda + dlambda)) lambda_new = lambda + dlambda;
551: lambdal = lambda;
552: rl = r;
553: lambda = lambda_new;
554: s = (lambdau - lambdal) / (lambdau - lambda);
555: }
556: }
557: r = phi(x, n, lambda, a, b, c, l, u);
558: }
560: *lam_ext = lambda;
561: if (innerIter >= df->maxProjIter) PetscInfo(NULL, "WARNING: DaiFletcher max iterations\n");
562: return innerIter;
563: }
565: PetscErrorCode solve(TAO_DF *df)
566: {
567: PetscInt i, j, innerIter, it, it2, luv, info;
568: PetscReal gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam = 0.0, lam_ext;
569: PetscReal DELTAsv, ProdDELTAsv;
570: PetscReal c, *tempQ;
571: PetscReal *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol;
572: PetscReal *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd;
573: PetscReal *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk;
574: PetscReal **Q = df->Q, *f = df->f, *t = df->t;
575: PetscInt dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv;
577: /* variables for the adaptive nonmonotone linesearch */
578: PetscInt L, llast;
579: PetscReal fr, fbest, fv, fc, fv0;
581: c = BMRM_INFTY;
583: DELTAsv = EPS_SV;
584: if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F;
585: else ProdDELTAsv = EPS_SV;
587: for (i = 0; i < dim; i++) tempv[i] = -x[i];
589: lam_ext = 0.0;
591: /* Project the initial solution */
592: project(dim, a, b, tempv, l, u, x, &lam_ext, df);
594: /* Compute gradient
595: g = Q*x + f; */
597: it = 0;
598: for (i = 0; i < dim; i++) {
599: if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i;
600: }
602: PetscArrayzero(t, dim);
603: for (i = 0; i < it; i++) {
604: tempQ = Q[ipt[i]];
605: for (j = 0; j < dim; j++) t[j] += (tempQ[j] * x[ipt[i]]);
606: }
607: for (i = 0; i < dim; i++) g[i] = t[i] + f[i];
609: /* y = -(x_{k} - g_{k}) */
610: for (i = 0; i < dim; i++) y[i] = g[i] - x[i];
612: /* Project x_{k} - g_{k} */
613: project(dim, a, b, y, l, u, tempv, &lam_ext, df);
615: /* y = P(x_{k} - g_{k}) - x_{k} */
616: max = ALPHA_MIN;
617: for (i = 0; i < dim; i++) {
618: y[i] = tempv[i] - x[i];
619: if (PetscAbsReal(y[i]) > max) max = PetscAbsReal(y[i]);
620: }
622: if (max < tol * 1e-3) return 0;
624: alpha = 1.0 / max;
626: /* fv0 = f(x_{0}). Recall t = Q x_{k} */
627: fv0 = 0.0;
628: for (i = 0; i < dim; i++) fv0 += x[i] * (0.5 * t[i] + f[i]);
630: /* adaptive nonmonotone linesearch */
631: L = 2;
632: fr = ALPHA_MAX;
633: fbest = fv0;
634: fc = fv0;
635: llast = 0;
636: akold = bkold = 0.0;
638: /* Iterator begins */
639: for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) {
640: /* tempv = -(x_{k} - alpha*g_{k}) */
641: for (i = 0; i < dim; i++) tempv[i] = alpha * g[i] - x[i];
643: /* Project x_{k} - alpha*g_{k} */
644: project(dim, a, b, tempv, l, u, y, &lam_ext, df);
646: /* gd = \inner{d_{k}}{g_{k}}
647: d = P(x_{k} - alpha*g_{k}) - x_{k}
648: */
649: gd = 0.0;
650: for (i = 0; i < dim; i++) {
651: d[i] = y[i] - x[i];
652: gd += d[i] * g[i];
653: }
655: /* Gradient computation */
657: /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */
659: it = it2 = 0;
660: for (i = 0; i < dim; i++) {
661: if (PetscAbsReal(d[i]) > (ProdDELTAsv * 1.0e-2)) ipt[it++] = i;
662: }
663: for (i = 0; i < dim; i++) {
664: if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i;
665: }
667: PetscArrayzero(Qd, dim);
668: /* compute Qd = Q*d */
669: if (it < it2) {
670: for (i = 0; i < it; i++) {
671: tempQ = Q[ipt[i]];
672: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]);
673: }
674: } else { /* compute Qd = Q*y-t */
675: for (i = 0; i < it2; i++) {
676: tempQ = Q[ipt2[i]];
677: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]);
678: }
679: for (j = 0; j < dim; j++) Qd[j] -= t[j];
680: }
682: /* ak = inner{d_{k}}{d_{k}} */
683: ak = 0.0;
684: for (i = 0; i < dim; i++) ak += d[i] * d[i];
686: bk = 0.0;
687: for (i = 0; i < dim; i++) bk += d[i] * Qd[i];
689: if (bk > EPS * ak && gd < 0.0) lamnew = -gd / bk;
690: else lamnew = 1.0;
692: /* fv is computing f(x_{k} + d_{k}) */
693: fv = 0.0;
694: for (i = 0; i < dim; i++) {
695: xplus[i] = x[i] + d[i];
696: tplus[i] = t[i] + Qd[i];
697: fv += xplus[i] * (0.5 * tplus[i] + f[i]);
698: }
700: /* fr is fref */
701: if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)) {
702: fv = 0.0;
703: for (i = 0; i < dim; i++) {
704: xplus[i] = x[i] + lamnew * d[i];
705: tplus[i] = t[i] + lamnew * Qd[i];
706: fv += xplus[i] * (0.5 * tplus[i] + f[i]);
707: }
708: }
710: for (i = 0; i < dim; i++) {
711: sk[i] = xplus[i] - x[i];
712: yk[i] = tplus[i] - t[i];
713: x[i] = xplus[i];
714: t[i] = tplus[i];
715: g[i] = t[i] + f[i];
716: }
718: /* update the line search control parameters */
719: if (fv < fbest) {
720: fbest = fv;
721: fc = fv;
722: llast = 0;
723: } else {
724: fc = (fc > fv ? fc : fv);
725: llast++;
726: if (llast == L) {
727: fr = fc;
728: fc = fv;
729: llast = 0;
730: }
731: }
733: ak = bk = 0.0;
734: for (i = 0; i < dim; i++) {
735: ak += sk[i] * sk[i];
736: bk += sk[i] * yk[i];
737: }
739: if (bk <= EPS * ak) alpha = ALPHA_MAX;
740: else {
741: if (bkold < EPS * akold) alpha = ak / bk;
742: else alpha = (akold + ak) / (bkold + bk);
744: if (alpha > ALPHA_MAX) alpha = ALPHA_MAX;
745: else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN;
746: }
748: akold = ak;
749: bkold = bk;
751: /* stopping criterion based on KKT conditions */
752: /* at optimal, gradient of lagrangian w.r.t. x is zero */
754: bk = 0.0;
755: for (i = 0; i < dim; i++) bk += x[i] * x[i];
757: if (PetscSqrtReal(ak) < tol * 10 * PetscSqrtReal(bk)) {
758: it = 0;
759: luv = 0;
760: kktlam = 0.0;
761: for (i = 0; i < dim; i++) {
762: /* x[i] is active hence lagrange multipliers for box constraints
763: are zero. The lagrange multiplier for ineq. const. is then
764: defined as below
765: */
766: if ((x[i] > DELTAsv) && (x[i] < c - DELTAsv)) {
767: ipt[it++] = i;
768: kktlam = kktlam - a[i] * g[i];
769: } else uv[luv++] = i;
770: }
772: if (it == 0 && PetscSqrtReal(ak) < tol * 0.5 * PetscSqrtReal(bk)) return 0;
773: else {
774: kktlam = kktlam / it;
775: info = 1;
776: for (i = 0; i < it; i++) {
777: if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) {
778: info = 0;
779: break;
780: }
781: }
782: if (info == 1) {
783: for (i = 0; i < luv; i++) {
784: if (x[uv[i]] <= DELTAsv) {
785: /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may
786: not be zero. So, the gradient without beta is > 0
787: */
788: if (g[uv[i]] + kktlam * a[uv[i]] < -tol) {
789: info = 0;
790: break;
791: }
792: } else {
793: /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may
794: not be zero. So, the gradient without eta is < 0
795: */
796: if (g[uv[i]] + kktlam * a[uv[i]] > tol) {
797: info = 0;
798: break;
799: }
800: }
801: }
802: }
804: if (info == 1) return 0;
805: }
806: }
807: }
808: return 0;
809: }