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