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 /* The main solver function 12 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) 15 16 Regularizer assumed to be L2 norm = lambda*0.5*W'W () 17 */ 18 19 static PetscErrorCode make_grad_node(Vec X, Vec_Chain **p) { 20 PetscFunctionBegin; 21 PetscCall(PetscNew(p)); 22 PetscCall(VecDuplicate(X, &(*p)->V)); 23 PetscCall(VecCopy(X, (*p)->V)); 24 (*p)->next = NULL; 25 PetscFunctionReturn(0); 26 } 27 28 static PetscErrorCode destroy_grad_list(Vec_Chain *head) { 29 Vec_Chain *p = head->next, *q; 30 31 PetscFunctionBegin; 32 while (p) { 33 q = p->next; 34 PetscCall(VecDestroy(&p->V)); 35 PetscCall(PetscFree(p)); 36 p = q; 37 } 38 head->next = NULL; 39 PetscFunctionReturn(0); 40 } 41 42 static PetscErrorCode TaoSolve_BMRM(Tao tao) { 43 TAO_DF df; 44 TAO_BMRM *bmrm = (TAO_BMRM *)tao->data; 45 46 /* Values and pointers to parts of the optimization problem */ 47 PetscReal f = 0.0; 48 Vec W = tao->solution; 49 Vec G = tao->gradient; 50 PetscReal lambda; 51 PetscReal bt; 52 Vec_Chain grad_list, *tail_glist, *pgrad; 53 PetscInt i; 54 PetscMPIInt rank; 55 56 /* Used in converged criteria check */ 57 PetscReal reg; 58 PetscReal jtwt = 0.0, max_jtwt, pre_epsilon, epsilon, jw, min_jw; 59 PetscReal innerSolverTol; 60 MPI_Comm comm; 61 62 PetscFunctionBegin; 63 PetscCall(PetscObjectGetComm((PetscObject)tao, &comm)); 64 PetscCallMPI(MPI_Comm_rank(comm, &rank)); 65 lambda = bmrm->lambda; 66 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; 73 74 if (rank == 0) { 75 PetscCall(init_df_solver(&df)); 76 grad_list.next = NULL; 77 tail_glist = &grad_list; 78 } 79 80 df.tol = 1e-6; 81 tao->reason = TAO_CONTINUE_ITERATING; 82 83 /*-----------------Algorithm Begins------------------------*/ 84 /* make the scatter */ 85 PetscCall(VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w)); 86 PetscCall(VecAssemblyBegin(bmrm->local_w)); 87 PetscCall(VecAssemblyEnd(bmrm->local_w)); 88 89 /* NOTE: In application pass the sub-gradient of Remp(W) */ 90 PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G)); 91 PetscCall(TaoLogConvergenceHistory(tao, f, 1.0, 0.0, tao->ksp_its)); 92 PetscCall(TaoMonitor(tao, tao->niter, f, 1.0, 0.0, tao->step)); 93 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 94 95 while (tao->reason == TAO_CONTINUE_ITERATING) { 96 /* Call general purpose update function */ 97 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update); 98 99 /* compute bt = Remp(Wt-1) - <Wt-1, At> */ 100 PetscCall(VecDot(W, G, &bt)); 101 bt = f - bt; 102 103 /* First gather the gradient to the rank-0 node */ 104 PetscCall(VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD)); 105 PetscCall(VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD)); 106 107 /* Bring up the inner solver */ 108 if (rank == 0) { 109 PetscCall(ensure_df_space(tao->niter + 1, &df)); 110 PetscCall(make_grad_node(bmrm->local_w, &pgrad)); 111 tail_glist->next = pgrad; 112 tail_glist = pgrad; 113 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; 118 119 /* set up the Q */ 120 pgrad = grad_list.next; 121 for (i = 0; i <= tao->niter; i++) { 122 PetscCheck(pgrad, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Assert that there are at least tao->niter+1 pgrad available"); 123 PetscCall(VecDot(pgrad->V, bmrm->local_w, ®)); 124 df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda; 125 pgrad = pgrad->next; 126 } 127 128 if (tao->niter > 0) { 129 df.x[tao->niter] = 0.0; 130 PetscCall(solve(&df)); 131 } else df.x[0] = 1.0; 132 133 /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */ 134 jtwt = 0.0; 135 PetscCall(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 PetscCall(VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V)); 140 pgrad = pgrad->next; 141 } 142 143 PetscCall(VecNorm(bmrm->local_w, NORM_2, ®)); 144 reg = 0.5 * lambda * reg * reg; 145 jtwt -= reg; 146 } /* end if rank == 0 */ 147 148 /* scatter the new W to all nodes */ 149 PetscCall(VecScatterBegin(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE)); 150 PetscCall(VecScatterEnd(bmrm->scatter, bmrm->local_w, W, INSERT_VALUES, SCATTER_REVERSE)); 151 152 PetscCall(TaoComputeObjectiveAndGradient(tao, W, &f, G)); 153 154 PetscCallMPI(MPI_Bcast(&jtwt, 1, MPIU_REAL, 0, comm)); 155 PetscCallMPI(MPI_Bcast(®, 1, MPIU_REAL, 0, comm)); 156 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; 160 161 pre_epsilon = epsilon; 162 epsilon = min_jw - jtwt; 163 164 if (rank == 0) { 165 if (innerSolverTol > epsilon) innerSolverTol = epsilon; 166 else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7; 167 168 /* if the annealing doesn't work well, lower the inner solver tolerance */ 169 if (pre_epsilon < epsilon) innerSolverTol *= 0.2; 170 171 df.tol = innerSolverTol * 0.5; 172 } 173 174 tao->niter++; 175 PetscCall(TaoLogConvergenceHistory(tao, min_jw, epsilon, 0.0, tao->ksp_its)); 176 PetscCall(TaoMonitor(tao, tao->niter, min_jw, epsilon, 0.0, tao->step)); 177 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 178 } 179 180 /* free all the memory */ 181 if (rank == 0) { 182 PetscCall(destroy_grad_list(&grad_list)); 183 PetscCall(destroy_df_solver(&df)); 184 } 185 186 PetscCall(VecDestroy(&bmrm->local_w)); 187 PetscCall(VecScatterDestroy(&bmrm->scatter)); 188 PetscFunctionReturn(0); 189 } 190 191 /* ---------------------------------------------------------- */ 192 193 static PetscErrorCode TaoSetup_BMRM(Tao tao) { 194 PetscFunctionBegin; 195 /* Allocate some arrays */ 196 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 197 PetscFunctionReturn(0); 198 } 199 200 /*------------------------------------------------------------*/ 201 static PetscErrorCode TaoDestroy_BMRM(Tao tao) { 202 PetscFunctionBegin; 203 PetscCall(PetscFree(tao->data)); 204 PetscFunctionReturn(0); 205 } 206 207 static PetscErrorCode TaoSetFromOptions_BMRM(Tao tao, PetscOptionItems *PetscOptionsObject) { 208 TAO_BMRM *bmrm = (TAO_BMRM *)tao->data; 209 210 PetscFunctionBegin; 211 PetscOptionsHeadBegin(PetscOptionsObject, "BMRM for regularized risk minimization"); 212 PetscCall(PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight", "", 100, &bmrm->lambda, NULL)); 213 PetscOptionsHeadEnd(); 214 PetscFunctionReturn(0); 215 } 216 217 /*------------------------------------------------------------*/ 218 static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer) { 219 PetscBool isascii; 220 221 PetscFunctionBegin; 222 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 223 if (isascii) { 224 PetscCall(PetscViewerASCIIPushTab(viewer)); 225 PetscCall(PetscViewerASCIIPopTab(viewer)); 226 } 227 PetscFunctionReturn(0); 228 } 229 230 /*------------------------------------------------------------*/ 231 /*MC 232 TAOBMRM - bundle method for regularized risk minimization 233 234 Options Database Keys: 235 . - tao_bmrm_lambda - regulariser weight 236 237 Level: beginner 238 M*/ 239 240 PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao) { 241 TAO_BMRM *bmrm; 242 243 PetscFunctionBegin; 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; 249 250 PetscCall(PetscNewLog(tao, &bmrm)); 251 bmrm->lambda = 1.0; 252 tao->data = (void *)bmrm; 253 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; 259 260 PetscFunctionReturn(0); 261 } 262 263 PetscErrorCode init_df_solver(TAO_DF *df) { 264 PetscInt i, n = INCRE_DIM; 265 266 PetscFunctionBegin; 267 /* default values */ 268 df->maxProjIter = 200; 269 df->maxPGMIter = 300000; 270 df->b = 1.0; 271 272 /* memory space required by Dai-Fletcher */ 273 df->cur_num_cp = n; 274 PetscCall(PetscMalloc1(n, &df->f)); 275 PetscCall(PetscMalloc1(n, &df->a)); 276 PetscCall(PetscMalloc1(n, &df->l)); 277 PetscCall(PetscMalloc1(n, &df->u)); 278 PetscCall(PetscMalloc1(n, &df->x)); 279 PetscCall(PetscMalloc1(n, &df->Q)); 280 281 for (i = 0; i < n; i++) { PetscCall(PetscMalloc1(n, &df->Q[i])); } 282 283 PetscCall(PetscMalloc1(n, &df->g)); 284 PetscCall(PetscMalloc1(n, &df->y)); 285 PetscCall(PetscMalloc1(n, &df->tempv)); 286 PetscCall(PetscMalloc1(n, &df->d)); 287 PetscCall(PetscMalloc1(n, &df->Qd)); 288 PetscCall(PetscMalloc1(n, &df->t)); 289 PetscCall(PetscMalloc1(n, &df->xplus)); 290 PetscCall(PetscMalloc1(n, &df->tplus)); 291 PetscCall(PetscMalloc1(n, &df->sk)); 292 PetscCall(PetscMalloc1(n, &df->yk)); 293 294 PetscCall(PetscMalloc1(n, &df->ipt)); 295 PetscCall(PetscMalloc1(n, &df->ipt2)); 296 PetscCall(PetscMalloc1(n, &df->uv)); 297 PetscFunctionReturn(0); 298 } 299 300 PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df) { 301 PetscReal *tmp, **tmp_Q; 302 PetscInt i, n, old_n; 303 304 PetscFunctionBegin; 305 df->dim = dim; 306 if (dim <= df->cur_num_cp) PetscFunctionReturn(0); 307 308 old_n = df->cur_num_cp; 309 df->cur_num_cp += INCRE_DIM; 310 n = df->cur_num_cp; 311 312 /* memory space required by dai-fletcher */ 313 PetscCall(PetscMalloc1(n, &tmp)); 314 PetscCall(PetscArraycpy(tmp, df->f, old_n)); 315 PetscCall(PetscFree(df->f)); 316 df->f = tmp; 317 318 PetscCall(PetscMalloc1(n, &tmp)); 319 PetscCall(PetscArraycpy(tmp, df->a, old_n)); 320 PetscCall(PetscFree(df->a)); 321 df->a = tmp; 322 323 PetscCall(PetscMalloc1(n, &tmp)); 324 PetscCall(PetscArraycpy(tmp, df->l, old_n)); 325 PetscCall(PetscFree(df->l)); 326 df->l = tmp; 327 328 PetscCall(PetscMalloc1(n, &tmp)); 329 PetscCall(PetscArraycpy(tmp, df->u, old_n)); 330 PetscCall(PetscFree(df->u)); 331 df->u = tmp; 332 333 PetscCall(PetscMalloc1(n, &tmp)); 334 PetscCall(PetscArraycpy(tmp, df->x, old_n)); 335 PetscCall(PetscFree(df->x)); 336 df->x = tmp; 337 338 PetscCall(PetscMalloc1(n, &tmp_Q)); 339 for (i = 0; i < n; i++) { 340 PetscCall(PetscMalloc1(n, &tmp_Q[i])); 341 if (i < old_n) { 342 PetscCall(PetscArraycpy(tmp_Q[i], df->Q[i], old_n)); 343 PetscCall(PetscFree(df->Q[i])); 344 } 345 } 346 347 PetscCall(PetscFree(df->Q)); 348 df->Q = tmp_Q; 349 350 PetscCall(PetscFree(df->g)); 351 PetscCall(PetscMalloc1(n, &df->g)); 352 353 PetscCall(PetscFree(df->y)); 354 PetscCall(PetscMalloc1(n, &df->y)); 355 356 PetscCall(PetscFree(df->tempv)); 357 PetscCall(PetscMalloc1(n, &df->tempv)); 358 359 PetscCall(PetscFree(df->d)); 360 PetscCall(PetscMalloc1(n, &df->d)); 361 362 PetscCall(PetscFree(df->Qd)); 363 PetscCall(PetscMalloc1(n, &df->Qd)); 364 365 PetscCall(PetscFree(df->t)); 366 PetscCall(PetscMalloc1(n, &df->t)); 367 368 PetscCall(PetscFree(df->xplus)); 369 PetscCall(PetscMalloc1(n, &df->xplus)); 370 371 PetscCall(PetscFree(df->tplus)); 372 PetscCall(PetscMalloc1(n, &df->tplus)); 373 374 PetscCall(PetscFree(df->sk)); 375 PetscCall(PetscMalloc1(n, &df->sk)); 376 377 PetscCall(PetscFree(df->yk)); 378 PetscCall(PetscMalloc1(n, &df->yk)); 379 380 PetscCall(PetscFree(df->ipt)); 381 PetscCall(PetscMalloc1(n, &df->ipt)); 382 383 PetscCall(PetscFree(df->ipt2)); 384 PetscCall(PetscMalloc1(n, &df->ipt2)); 385 386 PetscCall(PetscFree(df->uv)); 387 PetscCall(PetscMalloc1(n, &df->uv)); 388 PetscFunctionReturn(0); 389 } 390 391 PetscErrorCode destroy_df_solver(TAO_DF *df) { 392 PetscInt i; 393 394 PetscFunctionBegin; 395 PetscCall(PetscFree(df->f)); 396 PetscCall(PetscFree(df->a)); 397 PetscCall(PetscFree(df->l)); 398 PetscCall(PetscFree(df->u)); 399 PetscCall(PetscFree(df->x)); 400 401 for (i = 0; i < df->cur_num_cp; i++) { PetscCall(PetscFree(df->Q[i])); } 402 PetscCall(PetscFree(df->Q)); 403 PetscCall(PetscFree(df->ipt)); 404 PetscCall(PetscFree(df->ipt2)); 405 PetscCall(PetscFree(df->uv)); 406 PetscCall(PetscFree(df->g)); 407 PetscCall(PetscFree(df->y)); 408 PetscCall(PetscFree(df->tempv)); 409 PetscCall(PetscFree(df->d)); 410 PetscCall(PetscFree(df->Qd)); 411 PetscCall(PetscFree(df->t)); 412 PetscCall(PetscFree(df->xplus)); 413 PetscCall(PetscFree(df->tplus)); 414 PetscCall(PetscFree(df->sk)); 415 PetscCall(PetscFree(df->yk)); 416 PetscFunctionReturn(0); 417 } 418 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 PetscReal r = 0.0; 422 PetscInt i; 423 424 for (i = 0; i < n; i++) { 425 x[i] = -c[i] + lambda * a[i]; 426 if (x[i] > u[i]) x[i] = u[i]; 427 else if (x[i] < l[i]) x[i] = l[i]; 428 r += a[i] * x[i]; 429 } 430 return r - b; 431 } 432 433 /** Modified Dai-Fletcher QP projector solves the problem: 434 * 435 * minimise 0.5*x'*x - c'*x 436 * subj to a'*x = b 437 * l \leq x \leq u 438 * 439 * \param c The point to be projected onto feasible set 440 */ 441 PetscInt project(PetscInt n, PetscReal *a, PetscReal b, PetscReal *c, PetscReal *l, PetscReal *u, PetscReal *x, PetscReal *lam_ext, TAO_DF *df) { 442 PetscReal lambda, lambdal, lambdau, dlambda, lambda_new; 443 PetscReal r, rl, ru, s; 444 PetscInt innerIter; 445 PetscBool nonNegativeSlack = PETSC_FALSE; 446 447 *lam_ext = 0; 448 lambda = 0; 449 dlambda = 0.5; 450 innerIter = 1; 451 452 /* \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b) 453 * 454 * Optimality conditions for \phi: 455 * 456 * 1. lambda <= 0 457 * 2. r <= 0 458 * 3. r*lambda == 0 459 */ 460 461 /* Bracketing Phase */ 462 r = phi(x, n, lambda, a, b, c, l, u); 463 464 if (nonNegativeSlack) { 465 /* inequality constraint, i.e., with \xi >= 0 constraint */ 466 if (r < TOL_R) return 0; 467 } else { 468 /* equality constraint ,i.e., without \xi >= 0 constraint */ 469 if (PetscAbsReal(r) < TOL_R) return 0; 470 } 471 472 if (r < 0.0) { 473 lambdal = lambda; 474 rl = r; 475 lambda = lambda + dlambda; 476 r = phi(x, n, lambda, a, b, c, l, u); 477 while (r < 0.0 && dlambda < BMRM_INFTY) { 478 lambdal = lambda; 479 s = rl / r - 1.0; 480 if (s < 0.1) s = 0.1; 481 dlambda = dlambda + dlambda / s; 482 lambda = lambda + dlambda; 483 rl = r; 484 r = phi(x, n, lambda, a, b, c, l, u); 485 } 486 lambdau = lambda; 487 ru = r; 488 } else { 489 lambdau = lambda; 490 ru = r; 491 lambda = lambda - dlambda; 492 r = phi(x, n, lambda, a, b, c, l, u); 493 while (r > 0.0 && dlambda > -BMRM_INFTY) { 494 lambdau = lambda; 495 s = ru / r - 1.0; 496 if (s < 0.1) s = 0.1; 497 dlambda = dlambda + dlambda / s; 498 lambda = lambda - dlambda; 499 ru = r; 500 r = phi(x, n, lambda, a, b, c, l, u); 501 } 502 lambdal = lambda; 503 rl = r; 504 } 505 506 PetscCheck(PetscAbsReal(dlambda) <= BMRM_INFTY, PETSC_COMM_SELF, PETSC_ERR_PLIB, "L2N2_DaiFletcherPGM detected Infeasible QP problem!"); 507 508 if (ru == 0) { return innerIter; } 509 510 /* Secant Phase */ 511 s = 1.0 - rl / ru; 512 dlambda = dlambda / s; 513 lambda = lambdau - dlambda; 514 r = phi(x, n, lambda, a, b, c, l, u); 515 516 while (PetscAbsReal(r) > TOL_R && dlambda > TOL_LAM * (1.0 + PetscAbsReal(lambda)) && innerIter < df->maxProjIter) { 517 innerIter++; 518 if (r > 0.0) { 519 if (s <= 2.0) { 520 lambdau = lambda; 521 ru = r; 522 s = 1.0 - rl / ru; 523 dlambda = (lambdau - lambdal) / s; 524 lambda = lambdau - dlambda; 525 } else { 526 s = ru / r - 1.0; 527 if (s < 0.1) s = 0.1; 528 dlambda = (lambdau - lambda) / s; 529 lambda_new = 0.75 * lambdal + 0.25 * lambda; 530 if (lambda_new < (lambda - dlambda)) lambda_new = lambda - dlambda; 531 lambdau = lambda; 532 ru = r; 533 lambda = lambda_new; 534 s = (lambdau - lambdal) / (lambdau - lambda); 535 } 536 } else { 537 if (s >= 2.0) { 538 lambdal = lambda; 539 rl = r; 540 s = 1.0 - rl / ru; 541 dlambda = (lambdau - lambdal) / s; 542 lambda = lambdau - dlambda; 543 } else { 544 s = rl / r - 1.0; 545 if (s < 0.1) s = 0.1; 546 dlambda = (lambda - lambdal) / s; 547 lambda_new = 0.75 * lambdau + 0.25 * lambda; 548 if (lambda_new > (lambda + dlambda)) lambda_new = lambda + dlambda; 549 lambdal = lambda; 550 rl = r; 551 lambda = lambda_new; 552 s = (lambdau - lambdal) / (lambdau - lambda); 553 } 554 } 555 r = phi(x, n, lambda, a, b, c, l, u); 556 } 557 558 *lam_ext = lambda; 559 if (innerIter >= df->maxProjIter) { PetscCall(PetscInfo(NULL, "WARNING: DaiFletcher max iterations\n")); } 560 return innerIter; 561 } 562 563 PetscErrorCode solve(TAO_DF *df) { 564 PetscInt i, j, innerIter, it, it2, luv, info; 565 PetscReal gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam = 0.0, lam_ext; 566 PetscReal DELTAsv, ProdDELTAsv; 567 PetscReal c, *tempQ; 568 PetscReal *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol; 569 PetscReal *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd; 570 PetscReal *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk; 571 PetscReal **Q = df->Q, *f = df->f, *t = df->t; 572 PetscInt dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv; 573 574 /* variables for the adaptive nonmonotone linesearch */ 575 PetscInt L, llast; 576 PetscReal fr, fbest, fv, fc, fv0; 577 578 c = BMRM_INFTY; 579 580 DELTAsv = EPS_SV; 581 if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F; 582 else ProdDELTAsv = EPS_SV; 583 584 for (i = 0; i < dim; i++) tempv[i] = -x[i]; 585 586 lam_ext = 0.0; 587 588 /* Project the initial solution */ 589 project(dim, a, b, tempv, l, u, x, &lam_ext, df); 590 591 /* Compute gradient 592 g = Q*x + f; */ 593 594 it = 0; 595 for (i = 0; i < dim; i++) { 596 if (PetscAbsReal(x[i]) > ProdDELTAsv) ipt[it++] = i; 597 } 598 599 PetscCall(PetscArrayzero(t, dim)); 600 for (i = 0; i < it; i++) { 601 tempQ = Q[ipt[i]]; 602 for (j = 0; j < dim; j++) t[j] += (tempQ[j] * x[ipt[i]]); 603 } 604 for (i = 0; i < dim; i++) { g[i] = t[i] + f[i]; } 605 606 /* y = -(x_{k} - g_{k}) */ 607 for (i = 0; i < dim; i++) { y[i] = g[i] - x[i]; } 608 609 /* Project x_{k} - g_{k} */ 610 project(dim, a, b, y, l, u, tempv, &lam_ext, df); 611 612 /* y = P(x_{k} - g_{k}) - x_{k} */ 613 max = ALPHA_MIN; 614 for (i = 0; i < dim; i++) { 615 y[i] = tempv[i] - x[i]; 616 if (PetscAbsReal(y[i]) > max) max = PetscAbsReal(y[i]); 617 } 618 619 if (max < tol * 1e-3) { return 0; } 620 621 alpha = 1.0 / max; 622 623 /* fv0 = f(x_{0}). Recall t = Q x_{k} */ 624 fv0 = 0.0; 625 for (i = 0; i < dim; i++) fv0 += x[i] * (0.5 * t[i] + f[i]); 626 627 /* adaptive nonmonotone linesearch */ 628 L = 2; 629 fr = ALPHA_MAX; 630 fbest = fv0; 631 fc = fv0; 632 llast = 0; 633 akold = bkold = 0.0; 634 635 /* Iterator begins */ 636 for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) { 637 /* tempv = -(x_{k} - alpha*g_{k}) */ 638 for (i = 0; i < dim; i++) tempv[i] = alpha * g[i] - x[i]; 639 640 /* Project x_{k} - alpha*g_{k} */ 641 project(dim, a, b, tempv, l, u, y, &lam_ext, df); 642 643 /* gd = \inner{d_{k}}{g_{k}} 644 d = P(x_{k} - alpha*g_{k}) - x_{k} 645 */ 646 gd = 0.0; 647 for (i = 0; i < dim; i++) { 648 d[i] = y[i] - x[i]; 649 gd += d[i] * g[i]; 650 } 651 652 /* Gradient computation */ 653 654 /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */ 655 656 it = it2 = 0; 657 for (i = 0; i < dim; i++) { 658 if (PetscAbsReal(d[i]) > (ProdDELTAsv * 1.0e-2)) ipt[it++] = i; 659 } 660 for (i = 0; i < dim; i++) { 661 if (PetscAbsReal(y[i]) > ProdDELTAsv) ipt2[it2++] = i; 662 } 663 664 PetscCall(PetscArrayzero(Qd, dim)); 665 /* compute Qd = Q*d */ 666 if (it < it2) { 667 for (i = 0; i < it; i++) { 668 tempQ = Q[ipt[i]]; 669 for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]); 670 } 671 } else { /* compute Qd = Q*y-t */ 672 for (i = 0; i < it2; i++) { 673 tempQ = Q[ipt2[i]]; 674 for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]); 675 } 676 for (j = 0; j < dim; j++) Qd[j] -= t[j]; 677 } 678 679 /* ak = inner{d_{k}}{d_{k}} */ 680 ak = 0.0; 681 for (i = 0; i < dim; i++) ak += d[i] * d[i]; 682 683 bk = 0.0; 684 for (i = 0; i < dim; i++) bk += d[i] * Qd[i]; 685 686 if (bk > EPS * ak && gd < 0.0) lamnew = -gd / bk; 687 else lamnew = 1.0; 688 689 /* fv is computing f(x_{k} + d_{k}) */ 690 fv = 0.0; 691 for (i = 0; i < dim; i++) { 692 xplus[i] = x[i] + d[i]; 693 tplus[i] = t[i] + Qd[i]; 694 fv += xplus[i] * (0.5 * tplus[i] + f[i]); 695 } 696 697 /* fr is fref */ 698 if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)) { 699 fv = 0.0; 700 for (i = 0; i < dim; i++) { 701 xplus[i] = x[i] + lamnew * d[i]; 702 tplus[i] = t[i] + lamnew * Qd[i]; 703 fv += xplus[i] * (0.5 * tplus[i] + f[i]); 704 } 705 } 706 707 for (i = 0; i < dim; i++) { 708 sk[i] = xplus[i] - x[i]; 709 yk[i] = tplus[i] - t[i]; 710 x[i] = xplus[i]; 711 t[i] = tplus[i]; 712 g[i] = t[i] + f[i]; 713 } 714 715 /* update the line search control parameters */ 716 if (fv < fbest) { 717 fbest = fv; 718 fc = fv; 719 llast = 0; 720 } else { 721 fc = (fc > fv ? fc : fv); 722 llast++; 723 if (llast == L) { 724 fr = fc; 725 fc = fv; 726 llast = 0; 727 } 728 } 729 730 ak = bk = 0.0; 731 for (i = 0; i < dim; i++) { 732 ak += sk[i] * sk[i]; 733 bk += sk[i] * yk[i]; 734 } 735 736 if (bk <= EPS * ak) alpha = ALPHA_MAX; 737 else { 738 if (bkold < EPS * akold) alpha = ak / bk; 739 else alpha = (akold + ak) / (bkold + bk); 740 741 if (alpha > ALPHA_MAX) alpha = ALPHA_MAX; 742 else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN; 743 } 744 745 akold = ak; 746 bkold = bk; 747 748 /* stopping criterion based on KKT conditions */ 749 /* at optimal, gradient of lagrangian w.r.t. x is zero */ 750 751 bk = 0.0; 752 for (i = 0; i < dim; i++) bk += x[i] * x[i]; 753 754 if (PetscSqrtReal(ak) < tol * 10 * PetscSqrtReal(bk)) { 755 it = 0; 756 luv = 0; 757 kktlam = 0.0; 758 for (i = 0; i < dim; i++) { 759 /* x[i] is active hence lagrange multipliers for box constraints 760 are zero. The lagrange multiplier for ineq. const. is then 761 defined as below 762 */ 763 if ((x[i] > DELTAsv) && (x[i] < c - DELTAsv)) { 764 ipt[it++] = i; 765 kktlam = kktlam - a[i] * g[i]; 766 } else uv[luv++] = i; 767 } 768 769 if (it == 0 && PetscSqrtReal(ak) < tol * 0.5 * PetscSqrtReal(bk)) return 0; 770 else { 771 kktlam = kktlam / it; 772 info = 1; 773 for (i = 0; i < it; i++) { 774 if (PetscAbsReal(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) { 775 info = 0; 776 break; 777 } 778 } 779 if (info == 1) { 780 for (i = 0; i < luv; i++) { 781 if (x[uv[i]] <= DELTAsv) { 782 /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may 783 not be zero. So, the gradient without beta is > 0 784 */ 785 if (g[uv[i]] + kktlam * a[uv[i]] < -tol) { 786 info = 0; 787 break; 788 } 789 } else { 790 /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may 791 not be zero. So, the gradient without eta is < 0 792 */ 793 if (g[uv[i]] + kktlam * a[uv[i]] > tol) { 794 info = 0; 795 break; 796 } 797 } 798 } 799 } 800 801 if (info == 1) return 0; 802 } 803 } 804 } 805 return 0; 806 } 807