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, &reg);
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, &reg);
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(&reg,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: }