1 #include <../src/tao/unconstrained/impls/bmrm/bmrm.h> 2 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*); 9 10 11 /*------------------------------------------------------------*/ 12 /* The main solver function 13 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) 16 17 Regularizer assumed to be L2 norm = lambda*0.5*W'W () 18 */ 19 20 static PetscErrorCode make_grad_node(Vec X, Vec_Chain **p) 21 { 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ierr = PetscNew(p);CHKERRQ(ierr); 26 ierr = VecDuplicate(X, &(*p)->V);CHKERRQ(ierr); 27 ierr = VecCopy(X, (*p)->V);CHKERRQ(ierr); 28 (*p)->next = NULL; 29 PetscFunctionReturn(0); 30 } 31 32 static PetscErrorCode destroy_grad_list(Vec_Chain *head) 33 { 34 PetscErrorCode ierr; 35 Vec_Chain *p = head->next, *q; 36 37 PetscFunctionBegin; 38 while(p) { 39 q = p->next; 40 ierr = VecDestroy(&p->V);CHKERRQ(ierr); 41 ierr = PetscFree(p);CHKERRQ(ierr); 42 p = q; 43 } 44 head->next = NULL; 45 PetscFunctionReturn(0); 46 } 47 48 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; 55 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; 65 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; 71 72 PetscFunctionBegin; 73 ierr = PetscObjectGetComm((PetscObject)tao,&comm);CHKERRQ(ierr); 74 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 75 lambda = bmrm->lambda; 76 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; 83 84 if (!rank) { 85 ierr = init_df_solver(&df);CHKERRQ(ierr); 86 grad_list.next = NULL; 87 tail_glist = &grad_list; 88 } 89 90 df.tol = 1e-6; 91 reason = TAO_CONTINUE_ITERATING; 92 93 /*-----------------Algorithm Begins------------------------*/ 94 /* make the scatter */ 95 ierr = VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w);CHKERRQ(ierr); 96 ierr = VecAssemblyBegin(bmrm->local_w);CHKERRQ(ierr); 97 ierr = VecAssemblyEnd(bmrm->local_w);CHKERRQ(ierr); 98 99 /* NOTE: In application pass the sub-gradient of Remp(W) */ 100 ierr = TaoComputeObjectiveAndGradient(tao, W, &f, G);CHKERRQ(ierr); 101 ierr = TaoMonitor(tao,tao->niter,f,1.0,0.0,tao->step,&reason);CHKERRQ(ierr); 102 while (reason == TAO_CONTINUE_ITERATING) { 103 /* compute bt = Remp(Wt-1) - <Wt-1, At> */ 104 ierr = VecDot(W, G, &bt);CHKERRQ(ierr); 105 bt = f - bt; 106 107 /* First gather the gradient to the master node */ 108 ierr = VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 109 ierr = VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);CHKERRQ(ierr); 110 111 /* Bring up the inner solver */ 112 if (!rank) { 113 ierr = ensure_df_space(tao->niter+1, &df);CHKERRQ(ierr); 114 ierr = make_grad_node(bmrm->local_w, &pgrad);CHKERRQ(ierr); 115 tail_glist->next = pgrad; 116 tail_glist = pgrad; 117 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; 122 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 ierr = VecDot(pgrad->V, bmrm->local_w, ®);CHKERRQ(ierr); 128 df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda; 129 pgrad = pgrad->next; 130 } 131 132 if (tao->niter > 0) { 133 df.x[tao->niter] = 0.0; 134 ierr = solve(&df);CHKERRQ(ierr); 135 } else 136 df.x[0] = 1.0; 137 138 /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */ 139 jtwt = 0.0; 140 ierr = VecSet(bmrm->local_w, 0.0);CHKERRQ(ierr); 141 pgrad = grad_list.next; 142 for (i=0; i<=tao->niter; i++) { 143 jtwt -= df.x[i] * df.f[i]; 144 ierr = VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V);CHKERRQ(ierr); 145 pgrad = pgrad->next; 146 } 147 148 ierr = VecNorm(bmrm->local_w, NORM_2, ®);CHKERRQ(ierr); 149 reg = 0.5*lambda*reg*reg; 150 jtwt -= reg; 151 } /* end if rank == 0 */ 152 153 /* scatter the new W to all nodes */ 154 ierr = VecScatterBegin(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 155 ierr = VecScatterEnd(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 156 157 ierr = TaoComputeObjectiveAndGradient(tao, W, &f, G);CHKERRQ(ierr); 158 159 ierr = MPI_Bcast(&jtwt,1,MPIU_REAL,0,comm);CHKERRQ(ierr); 160 ierr = MPI_Bcast(®,1,MPIU_REAL,0,comm);CHKERRQ(ierr); 161 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; 165 166 pre_epsilon = epsilon; 167 epsilon = min_jw - jtwt; 168 169 if (!rank) { 170 if (innerSolverTol > epsilon) innerSolverTol = epsilon; 171 else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7; 172 173 /* if the annealing doesn't work well, lower the inner solver tolerance */ 174 if(pre_epsilon < epsilon) innerSolverTol *= 0.2; 175 176 df.tol = innerSolverTol*0.5; 177 } 178 179 tao->niter++; 180 ierr = TaoMonitor(tao,tao->niter,min_jw,epsilon,0.0,tao->step,&reason);CHKERRQ(ierr); 181 } 182 183 /* free all the memory */ 184 if (!rank) { 185 ierr = destroy_grad_list(&grad_list);CHKERRQ(ierr); 186 ierr = destroy_df_solver(&df);CHKERRQ(ierr); 187 } 188 189 ierr = VecDestroy(&bmrm->local_w);CHKERRQ(ierr); 190 ierr = VecScatterDestroy(&bmrm->scatter);CHKERRQ(ierr); 191 PetscFunctionReturn(0); 192 } 193 194 195 /* ---------------------------------------------------------- */ 196 197 static PetscErrorCode TaoSetup_BMRM(Tao tao) 198 { 199 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 /* Allocate some arrays */ 204 if (!tao->gradient) { 205 ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr); 206 } 207 PetscFunctionReturn(0); 208 } 209 210 /*------------------------------------------------------------*/ 211 static PetscErrorCode TaoDestroy_BMRM(Tao tao) 212 { 213 PetscErrorCode ierr; 214 215 PetscFunctionBegin; 216 ierr = PetscFree(tao->data);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 220 static PetscErrorCode TaoSetFromOptions_BMRM(PetscOptionItems *PetscOptionsObject,Tao tao) 221 { 222 PetscErrorCode ierr; 223 TAO_BMRM* bmrm = (TAO_BMRM*)tao->data; 224 225 PetscFunctionBegin; 226 ierr = PetscOptionsHead(PetscOptionsObject,"BMRM for regularized risk minimization");CHKERRQ(ierr); 227 ierr = PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight","", 100,&bmrm->lambda,NULL);CHKERRQ(ierr); 228 ierr = PetscOptionsTail();CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231 232 /*------------------------------------------------------------*/ 233 static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer) 234 { 235 PetscBool isascii; 236 PetscErrorCode ierr; 237 238 PetscFunctionBegin; 239 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 240 if (isascii) { 241 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 242 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 243 } 244 PetscFunctionReturn(0); 245 } 246 247 /*------------------------------------------------------------*/ 248 /*MC 249 TAOBMRM - bundle method for regularized risk minimization 250 251 Options Database Keys: 252 . - tao_bmrm_lambda - regulariser weight 253 254 Level: beginner 255 M*/ 256 257 PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao) 258 { 259 TAO_BMRM *bmrm; 260 PetscErrorCode ierr; 261 262 PetscFunctionBegin; 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; 268 269 ierr = PetscNewLog(tao,&bmrm);CHKERRQ(ierr); 270 bmrm->lambda = 1.0; 271 tao->data = (void*)bmrm; 272 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; 278 279 PetscFunctionReturn(0); 280 } 281 282 PetscErrorCode init_df_solver(TAO_DF *df) 283 { 284 PetscInt i, n = INCRE_DIM; 285 PetscErrorCode ierr; 286 287 PetscFunctionBegin; 288 /* default values */ 289 df->maxProjIter = 200; 290 df->maxPGMIter = 300000; 291 df->b = 1.0; 292 293 /* memory space required by Dai-Fletcher */ 294 df->cur_num_cp = n; 295 ierr = PetscMalloc1(n, &df->f);CHKERRQ(ierr); 296 ierr = PetscMalloc1(n, &df->a);CHKERRQ(ierr); 297 ierr = PetscMalloc1(n, &df->l);CHKERRQ(ierr); 298 ierr = PetscMalloc1(n, &df->u);CHKERRQ(ierr); 299 ierr = PetscMalloc1(n, &df->x);CHKERRQ(ierr); 300 ierr = PetscMalloc1(n, &df->Q);CHKERRQ(ierr); 301 302 for (i = 0; i < n; i ++) { 303 ierr = PetscMalloc1(n, &df->Q[i]);CHKERRQ(ierr); 304 } 305 306 ierr = PetscMalloc1(n, &df->g);CHKERRQ(ierr); 307 ierr = PetscMalloc1(n, &df->y);CHKERRQ(ierr); 308 ierr = PetscMalloc1(n, &df->tempv);CHKERRQ(ierr); 309 ierr = PetscMalloc1(n, &df->d);CHKERRQ(ierr); 310 ierr = PetscMalloc1(n, &df->Qd);CHKERRQ(ierr); 311 ierr = PetscMalloc1(n, &df->t);CHKERRQ(ierr); 312 ierr = PetscMalloc1(n, &df->xplus);CHKERRQ(ierr); 313 ierr = PetscMalloc1(n, &df->tplus);CHKERRQ(ierr); 314 ierr = PetscMalloc1(n, &df->sk);CHKERRQ(ierr); 315 ierr = PetscMalloc1(n, &df->yk);CHKERRQ(ierr); 316 317 ierr = PetscMalloc1(n, &df->ipt);CHKERRQ(ierr); 318 ierr = PetscMalloc1(n, &df->ipt2);CHKERRQ(ierr); 319 ierr = PetscMalloc1(n, &df->uv);CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 323 PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df) 324 { 325 PetscErrorCode ierr; 326 PetscReal *tmp, **tmp_Q; 327 PetscInt i, n, old_n; 328 329 PetscFunctionBegin; 330 df->dim = dim; 331 if (dim <= df->cur_num_cp) PetscFunctionReturn(0); 332 333 old_n = df->cur_num_cp; 334 df->cur_num_cp += INCRE_DIM; 335 n = df->cur_num_cp; 336 337 /* memory space required by dai-fletcher */ 338 ierr = PetscMalloc1(n, &tmp);CHKERRQ(ierr); 339 ierr = PetscMemcpy(tmp, df->f, sizeof(PetscReal)*old_n);CHKERRQ(ierr); 340 ierr = PetscFree(df->f);CHKERRQ(ierr); 341 df->f = tmp; 342 343 ierr = PetscMalloc1(n, &tmp);CHKERRQ(ierr); 344 ierr = PetscMemcpy(tmp, df->a, sizeof(PetscReal)*old_n);CHKERRQ(ierr); 345 ierr = PetscFree(df->a);CHKERRQ(ierr); 346 df->a = tmp; 347 348 ierr = PetscMalloc1(n, &tmp);CHKERRQ(ierr); 349 ierr = PetscMemcpy(tmp, df->l, sizeof(PetscReal)*old_n);CHKERRQ(ierr); 350 ierr = PetscFree(df->l);CHKERRQ(ierr); 351 df->l = tmp; 352 353 ierr = PetscMalloc1(n, &tmp);CHKERRQ(ierr); 354 ierr = PetscMemcpy(tmp, df->u, sizeof(PetscReal)*old_n);CHKERRQ(ierr); 355 ierr = PetscFree(df->u);CHKERRQ(ierr); 356 df->u = tmp; 357 358 ierr = PetscMalloc1(n, &tmp);CHKERRQ(ierr); 359 ierr = PetscMemcpy(tmp, df->x, sizeof(PetscReal)*old_n);CHKERRQ(ierr); 360 ierr = PetscFree(df->x);CHKERRQ(ierr); 361 df->x = tmp; 362 363 ierr = PetscMalloc1(n, &tmp_Q);CHKERRQ(ierr); 364 for (i = 0; i < n; i ++) { 365 ierr = PetscMalloc1(n, &tmp_Q[i]);CHKERRQ(ierr); 366 if (i < old_n) { 367 ierr = PetscMemcpy(tmp_Q[i], df->Q[i], sizeof(PetscReal)*old_n);CHKERRQ(ierr); 368 ierr = PetscFree(df->Q[i]);CHKERRQ(ierr); 369 } 370 } 371 372 ierr = PetscFree(df->Q);CHKERRQ(ierr); 373 df->Q = tmp_Q; 374 375 ierr = PetscFree(df->g);CHKERRQ(ierr); 376 ierr = PetscMalloc1(n, &df->g);CHKERRQ(ierr); 377 378 ierr = PetscFree(df->y);CHKERRQ(ierr); 379 ierr = PetscMalloc1(n, &df->y);CHKERRQ(ierr); 380 381 ierr = PetscFree(df->tempv);CHKERRQ(ierr); 382 ierr = PetscMalloc1(n, &df->tempv);CHKERRQ(ierr); 383 384 ierr = PetscFree(df->d);CHKERRQ(ierr); 385 ierr = PetscMalloc1(n, &df->d);CHKERRQ(ierr); 386 387 ierr = PetscFree(df->Qd);CHKERRQ(ierr); 388 ierr = PetscMalloc1(n, &df->Qd);CHKERRQ(ierr); 389 390 ierr = PetscFree(df->t);CHKERRQ(ierr); 391 ierr = PetscMalloc1(n, &df->t);CHKERRQ(ierr); 392 393 ierr = PetscFree(df->xplus);CHKERRQ(ierr); 394 ierr = PetscMalloc1(n, &df->xplus);CHKERRQ(ierr); 395 396 ierr = PetscFree(df->tplus);CHKERRQ(ierr); 397 ierr = PetscMalloc1(n, &df->tplus);CHKERRQ(ierr); 398 399 ierr = PetscFree(df->sk);CHKERRQ(ierr); 400 ierr = PetscMalloc1(n, &df->sk);CHKERRQ(ierr); 401 402 ierr = PetscFree(df->yk);CHKERRQ(ierr); 403 ierr = PetscMalloc1(n, &df->yk);CHKERRQ(ierr); 404 405 ierr = PetscFree(df->ipt);CHKERRQ(ierr); 406 ierr = PetscMalloc1(n, &df->ipt);CHKERRQ(ierr); 407 408 ierr = PetscFree(df->ipt2);CHKERRQ(ierr); 409 ierr = PetscMalloc1(n, &df->ipt2);CHKERRQ(ierr); 410 411 ierr = PetscFree(df->uv);CHKERRQ(ierr); 412 ierr = PetscMalloc1(n, &df->uv);CHKERRQ(ierr); 413 PetscFunctionReturn(0); 414 } 415 416 PetscErrorCode destroy_df_solver(TAO_DF *df) 417 { 418 PetscErrorCode ierr; 419 PetscInt i; 420 421 PetscFunctionBegin; 422 ierr = PetscFree(df->f);CHKERRQ(ierr); 423 ierr = PetscFree(df->a);CHKERRQ(ierr); 424 ierr = PetscFree(df->l);CHKERRQ(ierr); 425 ierr = PetscFree(df->u);CHKERRQ(ierr); 426 ierr = PetscFree(df->x);CHKERRQ(ierr); 427 428 for (i = 0; i < df->cur_num_cp; i ++) { 429 ierr = PetscFree(df->Q[i]);CHKERRQ(ierr); 430 } 431 ierr = PetscFree(df->Q);CHKERRQ(ierr); 432 ierr = PetscFree(df->ipt);CHKERRQ(ierr); 433 ierr = PetscFree(df->ipt2);CHKERRQ(ierr); 434 ierr = PetscFree(df->uv);CHKERRQ(ierr); 435 ierr = PetscFree(df->g);CHKERRQ(ierr); 436 ierr = PetscFree(df->y);CHKERRQ(ierr); 437 ierr = PetscFree(df->tempv);CHKERRQ(ierr); 438 ierr = PetscFree(df->d);CHKERRQ(ierr); 439 ierr = PetscFree(df->Qd);CHKERRQ(ierr); 440 ierr = PetscFree(df->t);CHKERRQ(ierr); 441 ierr = PetscFree(df->xplus);CHKERRQ(ierr); 442 ierr = PetscFree(df->tplus);CHKERRQ(ierr); 443 ierr = PetscFree(df->sk);CHKERRQ(ierr); 444 ierr = PetscFree(df->yk);CHKERRQ(ierr); 445 PetscFunctionReturn(0); 446 } 447 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; 453 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 } 462 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; 477 PetscErrorCode ierr; 478 479 *lam_ext = 0; 480 lambda = 0; 481 dlambda = 0.5; 482 innerIter = 1; 483 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 */ 492 493 /* Bracketing Phase */ 494 r = phi(x, n, lambda, a, b, c, l, u); 495 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 } 503 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 } 537 538 if(PetscAbsReal(dlambda) > BMRM_INFTY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"L2N2_DaiFletcherPGM detected Infeasible QP problem!"); 539 540 if(ru == 0){ 541 return innerIter; 542 } 543 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); 549 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 } 595 596 *lam_ext = lambda; 597 if(innerIter >= df->maxProjIter) { 598 ierr = PetscPrintf(PETSC_COMM_SELF, "WARNING: DaiFletcher max iterations\n");CHKERRQ(ierr); 599 } 600 return innerIter; 601 } 602 603 604 PetscErrorCode solve(TAO_DF *df) 605 { 606 PetscErrorCode ierr; 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; 616 617 /*** variables for the adaptive nonmonotone linesearch ***/ 618 PetscInt L, llast; 619 PetscReal fr, fbest, fv, fc, fv0; 620 621 c = BMRM_INFTY; 622 623 DELTAsv = EPS_SV; 624 if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F; 625 else ProdDELTAsv = EPS_SV; 626 627 for (i = 0; i < dim; i++) tempv[i] = -x[i]; 628 629 lam_ext = 0.0; 630 631 /* Project the initial solution */ 632 projcount += project(dim, a, b, tempv, l, u, x, &lam_ext, df); 633 634 /* Compute gradient 635 g = Q*x + f; */ 636 637 it = 0; 638 for (i = 0; i < dim; i++) { 639 if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i; 640 } 641 642 ierr = PetscMemzero(t, dim*sizeof(PetscReal));CHKERRQ(ierr); 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 } 650 651 652 /* y = -(x_{k} - g_{k}) */ 653 for (i = 0; i < dim; i++){ 654 y[i] = g[i] - x[i]; 655 } 656 657 /* Project x_{k} - g_{k} */ 658 projcount += project(dim, a, b, y, l, u, tempv, &lam_ext, df); 659 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 } 666 667 if (max < tol*1e-3){ 668 return 0; 669 } 670 671 alpha = 1.0 / max; 672 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]); 676 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; 684 685 /*** Iterator begins ***/ 686 for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) { 687 688 /* tempv = -(x_{k} - alpha*g_{k}) */ 689 for (i = 0; i < dim; i++) tempv[i] = alpha*g[i] - x[i]; 690 691 /* Project x_{k} - alpha*g_{k} */ 692 projcount += project(dim, a, b, tempv, l, u, y, &lam_ext, df); 693 694 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 } 703 704 /* Gradient computation */ 705 706 /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */ 707 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 } 715 716 ierr = PetscMemzero(Qd, dim*sizeof(PetscReal));CHKERRQ(ierr); 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 } 730 731 /* ak = inner{d_{k}}{d_{k}} */ 732 ak = 0.0; 733 for (i = 0; i < dim; i++) ak += d[i] * d[i]; 734 735 bk = 0.0; 736 for (i = 0; i < dim; i++) bk += d[i]*Qd[i]; 737 738 if (bk > EPS*ak && gd < 0.0) lamnew = -gd/bk; 739 else lamnew = 1.0; 740 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 } 748 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 } 759 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 } 767 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 } 782 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 } 788 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); 793 794 if (alpha > ALPHA_MAX) alpha = ALPHA_MAX; 795 else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN; 796 } 797 798 akold = ak; 799 bkold = bk; 800 801 /*** stopping criterion based on KKT conditions ***/ 802 /* at optimal, gradient of lagrangian w.r.t. x is zero */ 803 804 bk = 0.0; 805 for (i = 0; i < dim; i++) bk += x[i] * x[i]; 806 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 } 821 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 } 853 854 if (info == 1) return 0; 855 } 856 } 857 } 858 return 0; 859 } 860 861 862