1 #include <../src/tao/matrix/lmvmmat.h> 2 #include <../src/tao/unconstrained/impls/ntr/ntr.h> 3 4 #include <petscksp.h> 5 6 #define NTR_KSP_NASH 0 7 #define NTR_KSP_STCG 1 8 #define NTR_KSP_GLTR 2 9 #define NTR_KSP_TYPES 3 10 11 #define NTR_PC_NONE 0 12 #define NTR_PC_AHESS 1 13 #define NTR_PC_BFGS 2 14 #define NTR_PC_PETSC 3 15 #define NTR_PC_TYPES 4 16 17 #define BFGS_SCALE_AHESS 0 18 #define BFGS_SCALE_BFGS 1 19 #define BFGS_SCALE_TYPES 2 20 21 #define NTR_INIT_CONSTANT 0 22 #define NTR_INIT_DIRECTION 1 23 #define NTR_INIT_INTERPOLATION 2 24 #define NTR_INIT_TYPES 3 25 26 #define NTR_UPDATE_REDUCTION 0 27 #define NTR_UPDATE_INTERPOLATION 1 28 #define NTR_UPDATE_TYPES 2 29 30 static const char *NTR_KSP[64] = { "nash", "stcg", "gltr"}; 31 32 static const char *NTR_PC[64] = { "none", "ahess", "bfgs", "petsc"}; 33 34 static const char *BFGS_SCALE[64] = { "ahess", "bfgs"}; 35 36 static const char *NTR_INIT[64] = { "constant", "direction", "interpolation"}; 37 38 static const char *NTR_UPDATE[64] = { "reduction", "interpolation"}; 39 40 /* Routine for BFGS preconditioner */ 41 static PetscErrorCode MatLMVMSolveShell(PC pc, Vec xin, Vec xout); 42 43 /* 44 TaoSolve_NTR - Implements Newton's Method with a trust region approach 45 for solving unconstrained minimization problems. 46 47 The basic algorithm is taken from MINPACK-2 (dstrn). 48 49 TaoSolve_NTR computes a local minimizer of a twice differentiable function 50 f by applying a trust region variant of Newton's method. At each stage 51 of the algorithm, we use the prconditioned conjugate gradient method to 52 determine an approximate minimizer of the quadratic equation 53 54 q(s) = <s, Hs + g> 55 56 subject to the trust region constraint 57 58 || s ||_M <= radius, 59 60 where radius is the trust region radius and M is a symmetric positive 61 definite matrix (the preconditioner). Here g is the gradient and H 62 is the Hessian matrix. 63 64 Note: TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG, 65 or KSPGLTR. Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this 66 routine regardless of what the user may have previously specified. 67 */ 68 #undef __FUNCT__ 69 #define __FUNCT__ "TaoSolve_NTR" 70 static PetscErrorCode TaoSolve_NTR(Tao tao) 71 { 72 TAO_NTR *tr = (TAO_NTR *)tao->data; 73 PC pc; 74 KSPConvergedReason ksp_reason; 75 TaoConvergedReason reason; 76 PetscReal fmin, ftrial, prered, actred, kappa, sigma, beta; 77 PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius; 78 PetscReal f, gnorm; 79 80 PetscReal delta; 81 PetscReal norm_d; 82 PetscErrorCode ierr; 83 PetscInt bfgsUpdates = 0; 84 PetscInt needH; 85 86 PetscInt i_max = 5; 87 PetscInt j_max = 1; 88 PetscInt i, j, N, n, its; 89 90 PetscFunctionBegin; 91 if (tao->XL || tao->XU || tao->ops->computebounds) { 92 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n");CHKERRQ(ierr); 93 } 94 95 tao->trust = tao->trust0; 96 97 /* Modify the radius if it is too large or small */ 98 tao->trust = PetscMax(tao->trust, tr->min_radius); 99 tao->trust = PetscMin(tao->trust, tr->max_radius); 100 101 102 if (NTR_PC_BFGS == tr->pc_type && !tr->M) { 103 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 104 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 105 ierr = MatCreateLMVM(((PetscObject)tao)->comm,n,N,&tr->M);CHKERRQ(ierr); 106 ierr = MatLMVMAllocateVectors(tr->M,tao->solution);CHKERRQ(ierr); 107 } 108 109 /* Check convergence criteria */ 110 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr); 111 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 112 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 113 needH = 1; 114 115 ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr); 116 if (reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 117 118 /* Create vectors for the limited memory preconditioner */ 119 if ((NTR_PC_BFGS == tr->pc_type) && 120 (BFGS_SCALE_BFGS != tr->bfgs_scale_type)) { 121 if (!tr->Diag) { 122 ierr = VecDuplicate(tao->solution, &tr->Diag);CHKERRQ(ierr); 123 } 124 } 125 126 switch(tr->ksp_type) { 127 case NTR_KSP_NASH: 128 ierr = KSPSetType(tao->ksp, KSPNASH);CHKERRQ(ierr); 129 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 130 break; 131 132 case NTR_KSP_STCG: 133 ierr = KSPSetType(tao->ksp, KSPSTCG);CHKERRQ(ierr); 134 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 135 break; 136 137 default: 138 ierr = KSPSetType(tao->ksp, KSPGLTR);CHKERRQ(ierr); 139 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 140 break; 141 } 142 143 /* Modify the preconditioner to use the bfgs approximation */ 144 ierr = KSPGetPC(tao->ksp, &pc);CHKERRQ(ierr); 145 switch(tr->pc_type) { 146 case NTR_PC_NONE: 147 ierr = PCSetType(pc, PCNONE);CHKERRQ(ierr); 148 ierr = PCSetFromOptions(pc);CHKERRQ(ierr); 149 break; 150 151 case NTR_PC_AHESS: 152 ierr = PCSetType(pc, PCJACOBI);CHKERRQ(ierr); 153 ierr = PCSetFromOptions(pc);CHKERRQ(ierr); 154 ierr = PCJacobiSetUseAbs(pc,PETSC_TRUE);CHKERRQ(ierr); 155 break; 156 157 case NTR_PC_BFGS: 158 ierr = PCSetType(pc, PCSHELL);CHKERRQ(ierr); 159 ierr = PCSetFromOptions(pc);CHKERRQ(ierr); 160 ierr = PCShellSetName(pc, "bfgs");CHKERRQ(ierr); 161 ierr = PCShellSetContext(pc, tr->M);CHKERRQ(ierr); 162 ierr = PCShellSetApply(pc, MatLMVMSolveShell);CHKERRQ(ierr); 163 break; 164 165 default: 166 /* Use the pc method set by pc_type */ 167 break; 168 } 169 170 /* Initialize trust-region radius */ 171 switch(tr->init_type) { 172 case NTR_INIT_CONSTANT: 173 /* Use the initial radius specified */ 174 break; 175 176 case NTR_INIT_INTERPOLATION: 177 /* Use the initial radius specified */ 178 max_radius = 0.0; 179 180 for (j = 0; j < j_max; ++j) { 181 fmin = f; 182 sigma = 0.0; 183 184 if (needH) { 185 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 186 needH = 0; 187 } 188 189 for (i = 0; i < i_max; ++i) { 190 191 ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr); 192 ierr = VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient);CHKERRQ(ierr); 193 ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr); 194 195 if (PetscIsInfOrNanReal(ftrial)) { 196 tau = tr->gamma1_i; 197 } 198 else { 199 if (ftrial < fmin) { 200 fmin = ftrial; 201 sigma = -tao->trust / gnorm; 202 } 203 204 ierr = MatMult(tao->hessian, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 205 ierr = VecDot(tao->gradient, tao->stepdirection, &prered);CHKERRQ(ierr); 206 207 prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm)); 208 actred = f - ftrial; 209 if ((PetscAbsScalar(actred) <= tr->epsilon) && 210 (PetscAbsScalar(prered) <= tr->epsilon)) { 211 kappa = 1.0; 212 } 213 else { 214 kappa = actred / prered; 215 } 216 217 tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred); 218 tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred); 219 tau_min = PetscMin(tau_1, tau_2); 220 tau_max = PetscMax(tau_1, tau_2); 221 222 if (PetscAbsScalar(kappa - 1.0) <= tr->mu1_i) { 223 /* Great agreement */ 224 max_radius = PetscMax(max_radius, tao->trust); 225 226 if (tau_max < 1.0) { 227 tau = tr->gamma3_i; 228 } 229 else if (tau_max > tr->gamma4_i) { 230 tau = tr->gamma4_i; 231 } 232 else { 233 tau = tau_max; 234 } 235 } 236 else if (PetscAbsScalar(kappa - 1.0) <= tr->mu2_i) { 237 /* Good agreement */ 238 max_radius = PetscMax(max_radius, tao->trust); 239 240 if (tau_max < tr->gamma2_i) { 241 tau = tr->gamma2_i; 242 } 243 else if (tau_max > tr->gamma3_i) { 244 tau = tr->gamma3_i; 245 } 246 else { 247 tau = tau_max; 248 } 249 } 250 else { 251 /* Not good agreement */ 252 if (tau_min > 1.0) { 253 tau = tr->gamma2_i; 254 } 255 else if (tau_max < tr->gamma1_i) { 256 tau = tr->gamma1_i; 257 } 258 else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) { 259 tau = tr->gamma1_i; 260 } 261 else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) && 262 ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) { 263 tau = tau_1; 264 } 265 else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) && 266 ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) { 267 tau = tau_2; 268 } 269 else { 270 tau = tau_max; 271 } 272 } 273 } 274 tao->trust = tau * tao->trust; 275 } 276 277 if (fmin < f) { 278 f = fmin; 279 ierr = VecAXPY(tao->solution, sigma, tao->gradient);CHKERRQ(ierr); 280 ierr = TaoComputeGradient(tao,tao->solution, tao->gradient);CHKERRQ(ierr); 281 282 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 283 284 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 285 needH = 1; 286 287 ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, 1.0, &reason);CHKERRQ(ierr); 288 if (reason != TAO_CONTINUE_ITERATING) { 289 PetscFunctionReturn(0); 290 } 291 } 292 } 293 tao->trust = PetscMax(tao->trust, max_radius); 294 295 /* Modify the radius if it is too large or small */ 296 tao->trust = PetscMax(tao->trust, tr->min_radius); 297 tao->trust = PetscMin(tao->trust, tr->max_radius); 298 break; 299 300 default: 301 /* Norm of the first direction will initialize radius */ 302 tao->trust = 0.0; 303 break; 304 } 305 306 /* Set initial scaling for the BFGS preconditioner 307 This step is done after computing the initial trust-region radius 308 since the function value may have decreased */ 309 if (NTR_PC_BFGS == tr->pc_type) { 310 if (f != 0.0) { 311 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 312 } 313 else { 314 delta = 2.0 / (gnorm*gnorm); 315 } 316 ierr = MatLMVMSetDelta(tr->M,delta);CHKERRQ(ierr); 317 } 318 319 /* Have not converged; continue with Newton method */ 320 while (reason == TAO_CONTINUE_ITERATING) { 321 ++tao->niter; 322 tao->ksp_its=0; 323 /* Compute the Hessian */ 324 if (needH) { 325 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 326 needH = 0; 327 } 328 329 if (NTR_PC_BFGS == tr->pc_type) { 330 if (BFGS_SCALE_AHESS == tr->bfgs_scale_type) { 331 /* Obtain diagonal for the bfgs preconditioner */ 332 ierr = MatGetDiagonal(tao->hessian, tr->Diag);CHKERRQ(ierr); 333 ierr = VecAbs(tr->Diag);CHKERRQ(ierr); 334 ierr = VecReciprocal(tr->Diag);CHKERRQ(ierr); 335 ierr = MatLMVMSetScale(tr->M,tr->Diag);CHKERRQ(ierr); 336 } 337 338 /* Update the limited memory preconditioner */ 339 ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr); 340 ++bfgsUpdates; 341 } 342 343 while (reason == TAO_CONTINUE_ITERATING) { 344 ierr = KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre);CHKERRQ(ierr); 345 346 /* Solve the trust region subproblem */ 347 if (NTR_KSP_NASH == tr->ksp_type) { 348 ierr = KSPNASHSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 349 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 350 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 351 tao->ksp_its+=its; 352 tao->ksp_tot_its+=its; 353 ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 354 } else if (NTR_KSP_STCG == tr->ksp_type) { 355 ierr = KSPSTCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 356 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 357 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 358 tao->ksp_its+=its; 359 tao->ksp_tot_its+=its; 360 ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 361 } else { /* NTR_KSP_GLTR */ 362 ierr = KSPGLTRSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 363 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 364 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 365 tao->ksp_its+=its; 366 tao->ksp_tot_its+=its; 367 ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 368 } 369 370 if (0.0 == tao->trust) { 371 /* Radius was uninitialized; use the norm of the direction */ 372 if (norm_d > 0.0) { 373 tao->trust = norm_d; 374 375 /* Modify the radius if it is too large or small */ 376 tao->trust = PetscMax(tao->trust, tr->min_radius); 377 tao->trust = PetscMin(tao->trust, tr->max_radius); 378 } 379 else { 380 /* The direction was bad; set radius to default value and re-solve 381 the trust-region subproblem to get a direction */ 382 tao->trust = tao->trust0; 383 384 /* Modify the radius if it is too large or small */ 385 tao->trust = PetscMax(tao->trust, tr->min_radius); 386 tao->trust = PetscMin(tao->trust, tr->max_radius); 387 388 if (NTR_KSP_NASH == tr->ksp_type) { 389 ierr = KSPNASHSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 390 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 391 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 392 tao->ksp_its+=its; 393 tao->ksp_tot_its+=its; 394 ierr = KSPNASHGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 395 } else if (NTR_KSP_STCG == tr->ksp_type) { 396 ierr = KSPSTCGSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 397 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 398 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 399 tao->ksp_its+=its; 400 tao->ksp_tot_its+=its; 401 ierr = KSPSTCGGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 402 } else { /* NTR_KSP_GLTR */ 403 ierr = KSPGLTRSetRadius(tao->ksp,tao->trust);CHKERRQ(ierr); 404 ierr = KSPSolve(tao->ksp, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 405 ierr = KSPGetIterationNumber(tao->ksp,&its);CHKERRQ(ierr); 406 tao->ksp_its+=its; 407 tao->ksp_tot_its+=its; 408 ierr = KSPGLTRGetNormD(tao->ksp, &norm_d);CHKERRQ(ierr); 409 } 410 411 if (norm_d == 0.0) SETERRQ(PETSC_COMM_SELF,1, "Initial direction zero"); 412 } 413 } 414 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 415 ierr = KSPGetConvergedReason(tao->ksp, &ksp_reason);CHKERRQ(ierr); 416 if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && 417 (NTR_PC_BFGS == tr->pc_type) && (bfgsUpdates > 1)) { 418 /* Preconditioner is numerically indefinite; reset the 419 approximate if using BFGS preconditioning. */ 420 421 if (f != 0.0) { 422 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 423 } 424 else { 425 delta = 2.0 / (gnorm*gnorm); 426 } 427 ierr = MatLMVMSetDelta(tr->M, delta);CHKERRQ(ierr); 428 ierr = MatLMVMReset(tr->M);CHKERRQ(ierr); 429 ierr = MatLMVMUpdate(tr->M, tao->solution, tao->gradient);CHKERRQ(ierr); 430 bfgsUpdates = 1; 431 } 432 433 if (NTR_UPDATE_REDUCTION == tr->update_type) { 434 /* Get predicted reduction */ 435 if (NTR_KSP_NASH == tr->ksp_type) { 436 ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 437 } else if (NTR_KSP_STCG == tr->ksp_type) { 438 ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 439 } else { /* gltr */ 440 ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 441 } 442 443 if (prered >= 0.0) { 444 /* The predicted reduction has the wrong sign. This cannot 445 happen in infinite precision arithmetic. Step should 446 be rejected! */ 447 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 448 } 449 else { 450 /* Compute trial step and function value */ 451 ierr = VecCopy(tao->solution,tr->W);CHKERRQ(ierr); 452 ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr); 453 ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr); 454 455 if (PetscIsInfOrNanReal(ftrial)) { 456 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 457 } else { 458 /* Compute and actual reduction */ 459 actred = f - ftrial; 460 prered = -prered; 461 if ((PetscAbsScalar(actred) <= tr->epsilon) && 462 (PetscAbsScalar(prered) <= tr->epsilon)) { 463 kappa = 1.0; 464 } 465 else { 466 kappa = actred / prered; 467 } 468 469 /* Accept or reject the step and update radius */ 470 if (kappa < tr->eta1) { 471 /* Reject the step */ 472 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 473 } 474 else { 475 /* Accept the step */ 476 if (kappa < tr->eta2) { 477 /* Marginal bad step */ 478 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d); 479 } 480 else if (kappa < tr->eta3) { 481 /* Reasonable step */ 482 tao->trust = tr->alpha3 * tao->trust; 483 } 484 else if (kappa < tr->eta4) { 485 /* Good step */ 486 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust); 487 } 488 else { 489 /* Very good step */ 490 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust); 491 } 492 break; 493 } 494 } 495 } 496 } 497 else { 498 /* Get predicted reduction */ 499 if (NTR_KSP_NASH == tr->ksp_type) { 500 ierr = KSPNASHGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 501 } else if (NTR_KSP_STCG == tr->ksp_type) { 502 ierr = KSPSTCGGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 503 } else { /* gltr */ 504 ierr = KSPGLTRGetObjFcn(tao->ksp,&prered);CHKERRQ(ierr); 505 } 506 507 if (prered >= 0.0) { 508 /* The predicted reduction has the wrong sign. This cannot 509 happen in infinite precision arithmetic. Step should 510 be rejected! */ 511 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 512 } 513 else { 514 ierr = VecCopy(tao->solution, tr->W);CHKERRQ(ierr); 515 ierr = VecAXPY(tr->W, 1.0, tao->stepdirection);CHKERRQ(ierr); 516 ierr = TaoComputeObjective(tao, tr->W, &ftrial);CHKERRQ(ierr); 517 if (PetscIsInfOrNanReal(ftrial)) { 518 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 519 } 520 else { 521 ierr = VecDot(tao->gradient, tao->stepdirection, &beta);CHKERRQ(ierr); 522 actred = f - ftrial; 523 prered = -prered; 524 if ((PetscAbsScalar(actred) <= tr->epsilon) && 525 (PetscAbsScalar(prered) <= tr->epsilon)) { 526 kappa = 1.0; 527 } 528 else { 529 kappa = actred / prered; 530 } 531 532 tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred); 533 tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred); 534 tau_min = PetscMin(tau_1, tau_2); 535 tau_max = PetscMax(tau_1, tau_2); 536 537 if (kappa >= 1.0 - tr->mu1) { 538 /* Great agreement; accept step and update radius */ 539 if (tau_max < 1.0) { 540 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d); 541 } 542 else if (tau_max > tr->gamma4) { 543 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d); 544 } 545 else { 546 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 547 } 548 break; 549 } 550 else if (kappa >= 1.0 - tr->mu2) { 551 /* Good agreement */ 552 553 if (tau_max < tr->gamma2) { 554 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d); 555 } 556 else if (tau_max > tr->gamma3) { 557 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d); 558 } 559 else if (tau_max < 1.0) { 560 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 561 } 562 else { 563 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 564 } 565 break; 566 } 567 else { 568 /* Not good agreement */ 569 if (tau_min > 1.0) { 570 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d); 571 } 572 else if (tau_max < tr->gamma1) { 573 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 574 } 575 else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) { 576 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 577 } 578 else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && 579 ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) { 580 tao->trust = tau_1 * PetscMin(tao->trust, norm_d); 581 } 582 else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && 583 ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) { 584 tao->trust = tau_2 * PetscMin(tao->trust, norm_d); 585 } 586 else { 587 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 588 } 589 } 590 } 591 } 592 } 593 594 /* The step computed was not good and the radius was decreased. 595 Monitor the radius to terminate. */ 596 ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust, &reason);CHKERRQ(ierr); 597 } 598 599 /* The radius may have been increased; modify if it is too large */ 600 tao->trust = PetscMin(tao->trust, tr->max_radius); 601 602 if (reason == TAO_CONTINUE_ITERATING) { 603 ierr = VecCopy(tr->W, tao->solution);CHKERRQ(ierr); 604 f = ftrial; 605 ierr = TaoComputeGradient(tao, tao->solution, tao->gradient);CHKERRQ(ierr); 606 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 607 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 608 needH = 1; 609 ierr = TaoMonitor(tao, tao->niter, f, gnorm, 0.0, tao->trust, &reason);CHKERRQ(ierr); 610 } 611 } 612 PetscFunctionReturn(0); 613 } 614 615 /*------------------------------------------------------------*/ 616 #undef __FUNCT__ 617 #define __FUNCT__ "TaoSetUp_NTR" 618 static PetscErrorCode TaoSetUp_NTR(Tao tao) 619 { 620 TAO_NTR *tr = (TAO_NTR *)tao->data; 621 PetscErrorCode ierr; 622 623 PetscFunctionBegin; 624 625 if (!tao->gradient) {ierr = VecDuplicate(tao->solution, &tao->gradient);CHKERRQ(ierr);} 626 if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution, &tao->stepdirection);CHKERRQ(ierr);} 627 if (!tr->W) {ierr = VecDuplicate(tao->solution, &tr->W);CHKERRQ(ierr);} 628 629 tr->Diag = 0; 630 tr->M = 0; 631 632 633 PetscFunctionReturn(0); 634 } 635 636 /*------------------------------------------------------------*/ 637 #undef __FUNCT__ 638 #define __FUNCT__ "TaoDestroy_NTR" 639 static PetscErrorCode TaoDestroy_NTR(Tao tao) 640 { 641 TAO_NTR *tr = (TAO_NTR *)tao->data; 642 PetscErrorCode ierr; 643 644 PetscFunctionBegin; 645 if (tao->setupcalled) { 646 ierr = VecDestroy(&tr->W);CHKERRQ(ierr); 647 } 648 ierr = MatDestroy(&tr->M);CHKERRQ(ierr); 649 ierr = VecDestroy(&tr->Diag);CHKERRQ(ierr); 650 ierr = PetscFree(tao->data);CHKERRQ(ierr); 651 PetscFunctionReturn(0); 652 } 653 654 /*------------------------------------------------------------*/ 655 #undef __FUNCT__ 656 #define __FUNCT__ "TaoSetFromOptions_NTR" 657 static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao) 658 { 659 TAO_NTR *tr = (TAO_NTR *)tao->data; 660 PetscErrorCode ierr; 661 662 PetscFunctionBegin; 663 ierr = PetscOptionsHead(PetscOptionsObject,"Newton trust region method for unconstrained optimization");CHKERRQ(ierr); 664 ierr = PetscOptionsEList("-tao_ntr_ksp_type", "ksp type", "", NTR_KSP, NTR_KSP_TYPES, NTR_KSP[tr->ksp_type], &tr->ksp_type,NULL);CHKERRQ(ierr); 665 ierr = PetscOptionsEList("-tao_ntr_pc_type", "pc type", "", NTR_PC, NTR_PC_TYPES, NTR_PC[tr->pc_type], &tr->pc_type,NULL);CHKERRQ(ierr); 666 ierr = PetscOptionsEList("-tao_ntr_bfgs_scale_type", "bfgs scale type", "", BFGS_SCALE, BFGS_SCALE_TYPES, BFGS_SCALE[tr->bfgs_scale_type], &tr->bfgs_scale_type,NULL);CHKERRQ(ierr); 667 ierr = PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL);CHKERRQ(ierr); 668 ierr = PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL);CHKERRQ(ierr); 669 ierr = PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL);CHKERRQ(ierr); 670 ierr = PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL);CHKERRQ(ierr); 671 ierr = PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL);CHKERRQ(ierr); 672 ierr = PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL);CHKERRQ(ierr); 673 ierr = PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL);CHKERRQ(ierr); 674 ierr = PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL);CHKERRQ(ierr); 675 ierr = PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL);CHKERRQ(ierr); 676 ierr = PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL);CHKERRQ(ierr); 677 ierr = PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL);CHKERRQ(ierr); 678 ierr = PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL);CHKERRQ(ierr); 679 ierr = PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL);CHKERRQ(ierr); 680 ierr = PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL);CHKERRQ(ierr); 681 ierr = PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL);CHKERRQ(ierr); 682 ierr = PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL);CHKERRQ(ierr); 683 ierr = PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL);CHKERRQ(ierr); 684 ierr = PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL);CHKERRQ(ierr); 685 ierr = PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL);CHKERRQ(ierr); 686 ierr = PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL);CHKERRQ(ierr); 687 ierr = PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL);CHKERRQ(ierr); 688 ierr = PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL);CHKERRQ(ierr); 689 ierr = PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL);CHKERRQ(ierr); 690 ierr = PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL);CHKERRQ(ierr); 691 ierr = PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL);CHKERRQ(ierr); 692 ierr = PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL);CHKERRQ(ierr); 693 ierr = PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL);CHKERRQ(ierr); 694 ierr = PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL);CHKERRQ(ierr); 695 ierr = PetscOptionsTail();CHKERRQ(ierr); 696 ierr = KSPSetFromOptions(tao->ksp);CHKERRQ(ierr); 697 PetscFunctionReturn(0); 698 } 699 700 /*------------------------------------------------------------*/ 701 #undef __FUNCT__ 702 #define __FUNCT__ "TaoView_NTR" 703 static PetscErrorCode TaoView_NTR(Tao tao, PetscViewer viewer) 704 { 705 TAO_NTR *tr = (TAO_NTR *)tao->data; 706 PetscErrorCode ierr; 707 PetscInt nrejects; 708 PetscBool isascii; 709 710 PetscFunctionBegin; 711 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);CHKERRQ(ierr); 712 if (isascii) { 713 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 714 if (NTR_PC_BFGS == tr->pc_type && tr->M) { 715 ierr = MatLMVMGetRejects(tr->M, &nrejects);CHKERRQ(ierr); 716 ierr = PetscViewerASCIIPrintf(viewer, "Rejected matrix updates: %D\n", nrejects);CHKERRQ(ierr); 717 } 718 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 719 } 720 PetscFunctionReturn(0); 721 } 722 723 /*------------------------------------------------------------*/ 724 /*MC 725 TAONTR - Newton's method with trust region for unconstrained minimization. 726 At each iteration, the Newton trust region method solves the system. 727 NTR expects a KSP solver with a trust region radius. 728 min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k 729 730 Options Database Keys: 731 + -tao_ntr_ksp_type - "nash","stcg","gltr" 732 . -tao_ntr_pc_type - "none","ahess","bfgs","petsc" 733 . -tao_ntr_bfgs_scale_type - type of scaling with bfgs pc, "ahess" or "bfgs" 734 . -tao_ntr_init_type - "constant","direction","interpolation" 735 . -tao_ntr_update_type - "reduction","interpolation" 736 . -tao_ntr_min_radius - lower bound on trust region radius 737 . -tao_ntr_max_radius - upper bound on trust region radius 738 . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction 739 . -tao_ntr_mu1_i - mu1 interpolation init factor 740 . -tao_ntr_mu2_i - mu2 interpolation init factor 741 . -tao_ntr_gamma1_i - gamma1 interpolation init factor 742 . -tao_ntr_gamma2_i - gamma2 interpolation init factor 743 . -tao_ntr_gamma3_i - gamma3 interpolation init factor 744 . -tao_ntr_gamma4_i - gamma4 interpolation init factor 745 . -tao_ntr_theta_i - thetha1 interpolation init factor 746 . -tao_ntr_eta1 - eta1 reduction update factor 747 . -tao_ntr_eta2 - eta2 reduction update factor 748 . -tao_ntr_eta3 - eta3 reduction update factor 749 . -tao_ntr_eta4 - eta4 reduction update factor 750 . -tao_ntr_alpha1 - alpha1 reduction update factor 751 . -tao_ntr_alpha2 - alpha2 reduction update factor 752 . -tao_ntr_alpha3 - alpha3 reduction update factor 753 . -tao_ntr_alpha4 - alpha4 reduction update factor 754 . -tao_ntr_alpha4 - alpha4 reduction update factor 755 . -tao_ntr_mu1 - mu1 interpolation update 756 . -tao_ntr_mu2 - mu2 interpolation update 757 . -tao_ntr_gamma1 - gamma1 interpolcation update 758 . -tao_ntr_gamma2 - gamma2 interpolcation update 759 . -tao_ntr_gamma3 - gamma3 interpolcation update 760 . -tao_ntr_gamma4 - gamma4 interpolation update 761 - -tao_ntr_theta - theta interpolation update 762 763 Level: beginner 764 M*/ 765 766 #undef __FUNCT__ 767 #define __FUNCT__ "TaoCreate_NTR" 768 PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao) 769 { 770 TAO_NTR *tr; 771 PetscErrorCode ierr; 772 773 PetscFunctionBegin; 774 775 ierr = PetscNewLog(tao,&tr);CHKERRQ(ierr); 776 777 tao->ops->setup = TaoSetUp_NTR; 778 tao->ops->solve = TaoSolve_NTR; 779 tao->ops->view = TaoView_NTR; 780 tao->ops->setfromoptions = TaoSetFromOptions_NTR; 781 tao->ops->destroy = TaoDestroy_NTR; 782 783 /* Override default settings (unless already changed) */ 784 if (!tao->max_it_changed) tao->max_it = 50; 785 if (!tao->trust0_changed) tao->trust0 = 100.0; 786 tao->data = (void*)tr; 787 788 /* Standard trust region update parameters */ 789 tr->eta1 = 1.0e-4; 790 tr->eta2 = 0.25; 791 tr->eta3 = 0.50; 792 tr->eta4 = 0.90; 793 794 tr->alpha1 = 0.25; 795 tr->alpha2 = 0.50; 796 tr->alpha3 = 1.00; 797 tr->alpha4 = 2.00; 798 tr->alpha5 = 4.00; 799 800 /* Interpolation parameters */ 801 tr->mu1_i = 0.35; 802 tr->mu2_i = 0.50; 803 804 tr->gamma1_i = 0.0625; 805 tr->gamma2_i = 0.50; 806 tr->gamma3_i = 2.00; 807 tr->gamma4_i = 5.00; 808 809 tr->theta_i = 0.25; 810 811 /* Interpolation trust region update parameters */ 812 tr->mu1 = 0.10; 813 tr->mu2 = 0.50; 814 815 tr->gamma1 = 0.25; 816 tr->gamma2 = 0.50; 817 tr->gamma3 = 2.00; 818 tr->gamma4 = 4.00; 819 820 tr->theta = 0.05; 821 822 tr->min_radius = 1.0e-10; 823 tr->max_radius = 1.0e10; 824 tr->epsilon = 1.0e-6; 825 826 tr->ksp_type = NTR_KSP_STCG; 827 tr->pc_type = NTR_PC_BFGS; 828 tr->bfgs_scale_type = BFGS_SCALE_AHESS; 829 tr->init_type = NTR_INIT_INTERPOLATION; 830 tr->update_type = NTR_UPDATE_REDUCTION; 831 832 833 /* Set linear solver to default for trust region */ 834 ierr = KSPCreate(((PetscObject)tao)->comm, &tao->ksp);CHKERRQ(ierr); 835 ierr = KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix);CHKERRQ(ierr); 836 PetscFunctionReturn(0); 837 } 838 839 840 #undef __FUNCT__ 841 #define __FUNCT__ "MatLMVMSolveShell" 842 static PetscErrorCode MatLMVMSolveShell(PC pc, Vec b, Vec x) 843 { 844 PetscErrorCode ierr; 845 Mat M; 846 PetscFunctionBegin; 847 PetscValidHeaderSpecific(pc,PC_CLASSID,1); 848 PetscValidHeaderSpecific(b,VEC_CLASSID,2); 849 PetscValidHeaderSpecific(x,VEC_CLASSID,3); 850 ierr = PCShellGetContext(pc,(void**)&M);CHKERRQ(ierr); 851 ierr = MatLMVMSolve(M, b, x);CHKERRQ(ierr); 852 PetscFunctionReturn(0); 853 } 854