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