1 #include <../src/tao/unconstrained/impls/ntr/ntrimpl.h> 2 3 #include <petscksp.h> 4 5 #define NTR_INIT_CONSTANT 0 6 #define NTR_INIT_DIRECTION 1 7 #define NTR_INIT_INTERPOLATION 2 8 #define NTR_INIT_TYPES 3 9 10 #define NTR_UPDATE_REDUCTION 0 11 #define NTR_UPDATE_INTERPOLATION 1 12 #define NTR_UPDATE_TYPES 2 13 14 static const char *NTR_INIT[64] = {"constant","direction","interpolation"}; 15 16 static const char *NTR_UPDATE[64] = {"reduction","interpolation"}; 17 18 /* 19 TaoSolve_NTR - Implements Newton's Method with a trust region approach 20 for solving unconstrained minimization problems. 21 22 The basic algorithm is taken from MINPACK-2 (dstrn). 23 24 TaoSolve_NTR computes a local minimizer of a twice differentiable function 25 f by applying a trust region variant of Newton's method. At each stage 26 of the algorithm, we use the prconditioned conjugate gradient method to 27 determine an approximate minimizer of the quadratic equation 28 29 q(s) = <s, Hs + g> 30 31 subject to the trust region constraint 32 33 || s ||_M <= radius, 34 35 where radius is the trust region radius and M is a symmetric positive 36 definite matrix (the preconditioner). Here g is the gradient and H 37 is the Hessian matrix. 38 39 Note: TaoSolve_NTR MUST use the iterative solver KSPNASH, KSPSTCG, 40 or KSPGLTR. Thus, we set KSPNASH, KSPSTCG, or KSPGLTR in this 41 routine regardless of what the user may have previously specified. 42 */ 43 static PetscErrorCode TaoSolve_NTR(Tao tao) 44 { 45 TAO_NTR *tr = (TAO_NTR *)tao->data; 46 KSPType ksp_type; 47 PetscBool is_nash,is_stcg,is_gltr,is_bfgs,is_jacobi,is_symmetric,sym_set; 48 KSPConvergedReason ksp_reason; 49 PC pc; 50 PetscReal fmin, ftrial, prered, actred, kappa, sigma, beta; 51 PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius; 52 PetscReal f, gnorm; 53 54 PetscReal norm_d; 55 PetscInt bfgsUpdates = 0; 56 PetscInt needH; 57 58 PetscInt i_max = 5; 59 PetscInt j_max = 1; 60 PetscInt i, j, N, n, its; 61 62 PetscFunctionBegin; 63 if (tao->XL || tao->XU || tao->ops->computebounds) { 64 PetscCall(PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by ntr algorithm\n")); 65 } 66 67 PetscCall(KSPGetType(tao->ksp,&ksp_type)); 68 PetscCall(PetscStrcmp(ksp_type,KSPNASH,&is_nash)); 69 PetscCall(PetscStrcmp(ksp_type,KSPSTCG,&is_stcg)); 70 PetscCall(PetscStrcmp(ksp_type,KSPGLTR,&is_gltr)); 71 PetscCheck(is_nash || is_stcg || is_gltr,PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"TAO_NTR requires nash, stcg, or gltr for the KSP"); 72 73 /* Initialize the radius and modify if it is too large or small */ 74 tao->trust = tao->trust0; 75 tao->trust = PetscMax(tao->trust, tr->min_radius); 76 tao->trust = PetscMin(tao->trust, tr->max_radius); 77 78 /* Allocate the vectors needed for the BFGS approximation */ 79 PetscCall(KSPGetPC(tao->ksp, &pc)); 80 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs)); 81 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi)); 82 if (is_bfgs) { 83 tr->bfgs_pre = pc; 84 PetscCall(PCLMVMGetMatLMVM(tr->bfgs_pre, &tr->M)); 85 PetscCall(VecGetLocalSize(tao->solution, &n)); 86 PetscCall(VecGetSize(tao->solution, &N)); 87 PetscCall(MatSetSizes(tr->M, n, n, N, N)); 88 PetscCall(MatLMVMAllocate(tr->M, tao->solution, tao->gradient)); 89 PetscCall(MatIsSymmetricKnown(tr->M, &sym_set, &is_symmetric)); 90 PetscCheck(sym_set && is_symmetric,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric."); 91 } else if (is_jacobi) PetscCall(PCJacobiSetUseAbs(pc,PETSC_TRUE)); 92 93 /* Check convergence criteria */ 94 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient)); 95 PetscCall(TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm)); 96 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER,"User provided compute function generated Inf or NaN"); 97 needH = 1; 98 99 tao->reason = TAO_CONTINUE_ITERATING; 100 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 101 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0)); 102 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 103 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 104 105 /* Initialize trust-region radius */ 106 switch (tr->init_type) { 107 case NTR_INIT_CONSTANT: 108 /* Use the initial radius specified */ 109 break; 110 111 case NTR_INIT_INTERPOLATION: 112 /* Use the initial radius specified */ 113 max_radius = 0.0; 114 115 for (j = 0; j < j_max; ++j) { 116 fmin = f; 117 sigma = 0.0; 118 119 if (needH) { 120 PetscCall(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 121 needH = 0; 122 } 123 124 for (i = 0; i < i_max; ++i) { 125 126 PetscCall(VecCopy(tao->solution, tr->W)); 127 PetscCall(VecAXPY(tr->W, -tao->trust/gnorm, tao->gradient)); 128 PetscCall(TaoComputeObjective(tao, tr->W, &ftrial)); 129 130 if (PetscIsInfOrNanReal(ftrial)) { 131 tau = tr->gamma1_i; 132 } 133 else { 134 if (ftrial < fmin) { 135 fmin = ftrial; 136 sigma = -tao->trust / gnorm; 137 } 138 139 PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection)); 140 PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered)); 141 142 prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm)); 143 actred = f - ftrial; 144 if ((PetscAbsScalar(actred) <= tr->epsilon) && 145 (PetscAbsScalar(prered) <= tr->epsilon)) { 146 kappa = 1.0; 147 } 148 else { 149 kappa = actred / prered; 150 } 151 152 tau_1 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust + (1.0 - tr->theta_i) * prered - actred); 153 tau_2 = tr->theta_i * gnorm * tao->trust / (tr->theta_i * gnorm * tao->trust - (1.0 + tr->theta_i) * prered + actred); 154 tau_min = PetscMin(tau_1, tau_2); 155 tau_max = PetscMax(tau_1, tau_2); 156 157 if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu1_i) { 158 /* Great agreement */ 159 max_radius = PetscMax(max_radius, tao->trust); 160 161 if (tau_max < 1.0) { 162 tau = tr->gamma3_i; 163 } 164 else if (tau_max > tr->gamma4_i) { 165 tau = tr->gamma4_i; 166 } 167 else { 168 tau = tau_max; 169 } 170 } 171 else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tr->mu2_i) { 172 /* Good agreement */ 173 max_radius = PetscMax(max_radius, tao->trust); 174 175 if (tau_max < tr->gamma2_i) { 176 tau = tr->gamma2_i; 177 } 178 else if (tau_max > tr->gamma3_i) { 179 tau = tr->gamma3_i; 180 } 181 else { 182 tau = tau_max; 183 } 184 } 185 else { 186 /* Not good agreement */ 187 if (tau_min > 1.0) { 188 tau = tr->gamma2_i; 189 } 190 else if (tau_max < tr->gamma1_i) { 191 tau = tr->gamma1_i; 192 } 193 else if ((tau_min < tr->gamma1_i) && (tau_max >= 1.0)) { 194 tau = tr->gamma1_i; 195 } 196 else if ((tau_1 >= tr->gamma1_i) && (tau_1 < 1.0) && 197 ((tau_2 < tr->gamma1_i) || (tau_2 >= 1.0))) { 198 tau = tau_1; 199 } 200 else if ((tau_2 >= tr->gamma1_i) && (tau_2 < 1.0) && 201 ((tau_1 < tr->gamma1_i) || (tau_2 >= 1.0))) { 202 tau = tau_2; 203 } 204 else { 205 tau = tau_max; 206 } 207 } 208 } 209 tao->trust = tau * tao->trust; 210 } 211 212 if (fmin < f) { 213 f = fmin; 214 PetscCall(VecAXPY(tao->solution, sigma, tao->gradient)); 215 PetscCall(TaoComputeGradient(tao,tao->solution, tao->gradient)); 216 217 PetscCall(TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm)); 218 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 219 needH = 1; 220 221 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 222 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,1.0)); 223 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 224 if (tao->reason != TAO_CONTINUE_ITERATING) { 225 PetscFunctionReturn(0); 226 } 227 } 228 } 229 tao->trust = PetscMax(tao->trust, max_radius); 230 231 /* Modify the radius if it is too large or small */ 232 tao->trust = PetscMax(tao->trust, tr->min_radius); 233 tao->trust = PetscMin(tao->trust, tr->max_radius); 234 break; 235 236 default: 237 /* Norm of the first direction will initialize radius */ 238 tao->trust = 0.0; 239 break; 240 } 241 242 /* Have not converged; continue with Newton method */ 243 while (tao->reason == TAO_CONTINUE_ITERATING) { 244 /* Call general purpose update function */ 245 if (tao->ops->update) PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 246 ++tao->niter; 247 tao->ksp_its=0; 248 /* Compute the Hessian */ 249 if (needH) { 250 PetscCall(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 251 needH = 0; 252 } 253 254 if (tr->bfgs_pre) { 255 /* Update the limited memory preconditioner */ 256 PetscCall(MatLMVMUpdate(tr->M, tao->solution, tao->gradient)); 257 ++bfgsUpdates; 258 } 259 260 while (tao->reason == TAO_CONTINUE_ITERATING) { 261 PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre)); 262 263 /* Solve the trust region subproblem */ 264 PetscCall(KSPCGSetRadius(tao->ksp,tao->trust)); 265 PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection)); 266 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 267 tao->ksp_its+=its; 268 tao->ksp_tot_its+=its; 269 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 270 271 if (0.0 == tao->trust) { 272 /* Radius was uninitialized; use the norm of the direction */ 273 if (norm_d > 0.0) { 274 tao->trust = norm_d; 275 276 /* Modify the radius if it is too large or small */ 277 tao->trust = PetscMax(tao->trust, tr->min_radius); 278 tao->trust = PetscMin(tao->trust, tr->max_radius); 279 } 280 else { 281 /* The direction was bad; set radius to default value and re-solve 282 the trust-region subproblem to get a direction */ 283 tao->trust = tao->trust0; 284 285 /* Modify the radius if it is too large or small */ 286 tao->trust = PetscMax(tao->trust, tr->min_radius); 287 tao->trust = PetscMin(tao->trust, tr->max_radius); 288 289 PetscCall(KSPCGSetRadius(tao->ksp,tao->trust)); 290 PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection)); 291 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 292 tao->ksp_its+=its; 293 tao->ksp_tot_its+=its; 294 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 295 296 PetscCheck(norm_d != 0.0,PetscObjectComm((PetscObject)tao),PETSC_ERR_PLIB, "Initial direction zero"); 297 } 298 } 299 PetscCall(VecScale(tao->stepdirection, -1.0)); 300 PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason)); 301 if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tr->bfgs_pre)) { 302 /* Preconditioner is numerically indefinite; reset the 303 approximate if using BFGS preconditioning. */ 304 PetscCall(MatLMVMReset(tr->M, PETSC_FALSE)); 305 PetscCall(MatLMVMUpdate(tr->M, tao->solution, tao->gradient)); 306 bfgsUpdates = 1; 307 } 308 309 if (NTR_UPDATE_REDUCTION == tr->update_type) { 310 /* Get predicted reduction */ 311 PetscCall(KSPCGGetObjFcn(tao->ksp,&prered)); 312 if (prered >= 0.0) { 313 /* The predicted reduction has the wrong sign. This cannot 314 happen in infinite precision arithmetic. Step should 315 be rejected! */ 316 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 317 } 318 else { 319 /* Compute trial step and function value */ 320 PetscCall(VecCopy(tao->solution,tr->W)); 321 PetscCall(VecAXPY(tr->W, 1.0, tao->stepdirection)); 322 PetscCall(TaoComputeObjective(tao, tr->W, &ftrial)); 323 324 if (PetscIsInfOrNanReal(ftrial)) { 325 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 326 } else { 327 /* Compute and actual reduction */ 328 actred = f - ftrial; 329 prered = -prered; 330 if ((PetscAbsScalar(actred) <= tr->epsilon) && 331 (PetscAbsScalar(prered) <= tr->epsilon)) { 332 kappa = 1.0; 333 } 334 else { 335 kappa = actred / prered; 336 } 337 338 /* Accept or reject the step and update radius */ 339 if (kappa < tr->eta1) { 340 /* Reject the step */ 341 tao->trust = tr->alpha1 * PetscMin(tao->trust, norm_d); 342 } 343 else { 344 /* Accept the step */ 345 if (kappa < tr->eta2) { 346 /* Marginal bad step */ 347 tao->trust = tr->alpha2 * PetscMin(tao->trust, norm_d); 348 } 349 else if (kappa < tr->eta3) { 350 /* Reasonable step */ 351 tao->trust = tr->alpha3 * tao->trust; 352 } 353 else if (kappa < tr->eta4) { 354 /* Good step */ 355 tao->trust = PetscMax(tr->alpha4 * norm_d, tao->trust); 356 } 357 else { 358 /* Very good step */ 359 tao->trust = PetscMax(tr->alpha5 * norm_d, tao->trust); 360 } 361 break; 362 } 363 } 364 } 365 } 366 else { 367 /* Get predicted reduction */ 368 PetscCall(KSPCGGetObjFcn(tao->ksp,&prered)); 369 if (prered >= 0.0) { 370 /* The predicted reduction has the wrong sign. This cannot 371 happen in infinite precision arithmetic. Step should 372 be rejected! */ 373 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 374 } 375 else { 376 PetscCall(VecCopy(tao->solution, tr->W)); 377 PetscCall(VecAXPY(tr->W, 1.0, tao->stepdirection)); 378 PetscCall(TaoComputeObjective(tao, tr->W, &ftrial)); 379 if (PetscIsInfOrNanReal(ftrial)) { 380 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 381 } 382 else { 383 PetscCall(VecDot(tao->gradient, tao->stepdirection, &beta)); 384 actred = f - ftrial; 385 prered = -prered; 386 if ((PetscAbsScalar(actred) <= tr->epsilon) && 387 (PetscAbsScalar(prered) <= tr->epsilon)) { 388 kappa = 1.0; 389 } 390 else { 391 kappa = actred / prered; 392 } 393 394 tau_1 = tr->theta * beta / (tr->theta * beta - (1.0 - tr->theta) * prered + actred); 395 tau_2 = tr->theta * beta / (tr->theta * beta + (1.0 + tr->theta) * prered - actred); 396 tau_min = PetscMin(tau_1, tau_2); 397 tau_max = PetscMax(tau_1, tau_2); 398 399 if (kappa >= 1.0 - tr->mu1) { 400 /* Great agreement; accept step and update radius */ 401 if (tau_max < 1.0) { 402 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d); 403 } 404 else if (tau_max > tr->gamma4) { 405 tao->trust = PetscMax(tao->trust, tr->gamma4 * norm_d); 406 } 407 else { 408 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 409 } 410 break; 411 } 412 else if (kappa >= 1.0 - tr->mu2) { 413 /* Good agreement */ 414 415 if (tau_max < tr->gamma2) { 416 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d); 417 } 418 else if (tau_max > tr->gamma3) { 419 tao->trust = PetscMax(tao->trust, tr->gamma3 * norm_d); 420 } 421 else if (tau_max < 1.0) { 422 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 423 } 424 else { 425 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 426 } 427 break; 428 } 429 else { 430 /* Not good agreement */ 431 if (tau_min > 1.0) { 432 tao->trust = tr->gamma2 * PetscMin(tao->trust, norm_d); 433 } 434 else if (tau_max < tr->gamma1) { 435 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 436 } 437 else if ((tau_min < tr->gamma1) && (tau_max >= 1.0)) { 438 tao->trust = tr->gamma1 * PetscMin(tao->trust, norm_d); 439 } 440 else if ((tau_1 >= tr->gamma1) && (tau_1 < 1.0) && 441 ((tau_2 < tr->gamma1) || (tau_2 >= 1.0))) { 442 tao->trust = tau_1 * PetscMin(tao->trust, norm_d); 443 } 444 else if ((tau_2 >= tr->gamma1) && (tau_2 < 1.0) && 445 ((tau_1 < tr->gamma1) || (tau_2 >= 1.0))) { 446 tao->trust = tau_2 * PetscMin(tao->trust, norm_d); 447 } 448 else { 449 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 450 } 451 } 452 } 453 } 454 } 455 456 /* The step computed was not good and the radius was decreased. 457 Monitor the radius to terminate. */ 458 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 459 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust)); 460 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 461 } 462 463 /* The radius may have been increased; modify if it is too large */ 464 tao->trust = PetscMin(tao->trust, tr->max_radius); 465 466 if (tao->reason == TAO_CONTINUE_ITERATING) { 467 PetscCall(VecCopy(tr->W, tao->solution)); 468 f = ftrial; 469 PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient)); 470 PetscCall(TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm)); 471 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 472 needH = 1; 473 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 474 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,tao->trust)); 475 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 476 } 477 } 478 PetscFunctionReturn(0); 479 } 480 481 /*------------------------------------------------------------*/ 482 static PetscErrorCode TaoSetUp_NTR(Tao tao) 483 { 484 TAO_NTR *tr = (TAO_NTR *)tao->data; 485 486 PetscFunctionBegin; 487 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 488 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 489 if (!tr->W) PetscCall(VecDuplicate(tao->solution, &tr->W)); 490 491 tr->bfgs_pre = NULL; 492 tr->M = NULL; 493 PetscFunctionReturn(0); 494 } 495 496 /*------------------------------------------------------------*/ 497 static PetscErrorCode TaoDestroy_NTR(Tao tao) 498 { 499 TAO_NTR *tr = (TAO_NTR *)tao->data; 500 501 PetscFunctionBegin; 502 if (tao->setupcalled) { 503 PetscCall(VecDestroy(&tr->W)); 504 } 505 PetscCall(KSPDestroy(&tao->ksp)); 506 PetscCall(PetscFree(tao->data)); 507 PetscFunctionReturn(0); 508 } 509 510 /*------------------------------------------------------------*/ 511 static PetscErrorCode TaoSetFromOptions_NTR(PetscOptionItems *PetscOptionsObject,Tao tao) 512 { 513 TAO_NTR *tr = (TAO_NTR *)tao->data; 514 515 PetscFunctionBegin; 516 PetscOptionsHeadBegin(PetscOptionsObject,"Newton trust region method for unconstrained optimization"); 517 PetscCall(PetscOptionsEList("-tao_ntr_init_type", "tao->trust initialization type", "", NTR_INIT, NTR_INIT_TYPES, NTR_INIT[tr->init_type], &tr->init_type,NULL)); 518 PetscCall(PetscOptionsEList("-tao_ntr_update_type", "radius update type", "", NTR_UPDATE, NTR_UPDATE_TYPES, NTR_UPDATE[tr->update_type], &tr->update_type,NULL)); 519 PetscCall(PetscOptionsReal("-tao_ntr_eta1", "step is unsuccessful if actual reduction < eta1 * predicted reduction", "", tr->eta1, &tr->eta1,NULL)); 520 PetscCall(PetscOptionsReal("-tao_ntr_eta2", "", "", tr->eta2, &tr->eta2,NULL)); 521 PetscCall(PetscOptionsReal("-tao_ntr_eta3", "", "", tr->eta3, &tr->eta3,NULL)); 522 PetscCall(PetscOptionsReal("-tao_ntr_eta4", "", "", tr->eta4, &tr->eta4,NULL)); 523 PetscCall(PetscOptionsReal("-tao_ntr_alpha1", "", "", tr->alpha1, &tr->alpha1,NULL)); 524 PetscCall(PetscOptionsReal("-tao_ntr_alpha2", "", "", tr->alpha2, &tr->alpha2,NULL)); 525 PetscCall(PetscOptionsReal("-tao_ntr_alpha3", "", "", tr->alpha3, &tr->alpha3,NULL)); 526 PetscCall(PetscOptionsReal("-tao_ntr_alpha4", "", "", tr->alpha4, &tr->alpha4,NULL)); 527 PetscCall(PetscOptionsReal("-tao_ntr_alpha5", "", "", tr->alpha5, &tr->alpha5,NULL)); 528 PetscCall(PetscOptionsReal("-tao_ntr_mu1", "", "", tr->mu1, &tr->mu1,NULL)); 529 PetscCall(PetscOptionsReal("-tao_ntr_mu2", "", "", tr->mu2, &tr->mu2,NULL)); 530 PetscCall(PetscOptionsReal("-tao_ntr_gamma1", "", "", tr->gamma1, &tr->gamma1,NULL)); 531 PetscCall(PetscOptionsReal("-tao_ntr_gamma2", "", "", tr->gamma2, &tr->gamma2,NULL)); 532 PetscCall(PetscOptionsReal("-tao_ntr_gamma3", "", "", tr->gamma3, &tr->gamma3,NULL)); 533 PetscCall(PetscOptionsReal("-tao_ntr_gamma4", "", "", tr->gamma4, &tr->gamma4,NULL)); 534 PetscCall(PetscOptionsReal("-tao_ntr_theta", "", "", tr->theta, &tr->theta,NULL)); 535 PetscCall(PetscOptionsReal("-tao_ntr_mu1_i", "", "", tr->mu1_i, &tr->mu1_i,NULL)); 536 PetscCall(PetscOptionsReal("-tao_ntr_mu2_i", "", "", tr->mu2_i, &tr->mu2_i,NULL)); 537 PetscCall(PetscOptionsReal("-tao_ntr_gamma1_i", "", "", tr->gamma1_i, &tr->gamma1_i,NULL)); 538 PetscCall(PetscOptionsReal("-tao_ntr_gamma2_i", "", "", tr->gamma2_i, &tr->gamma2_i,NULL)); 539 PetscCall(PetscOptionsReal("-tao_ntr_gamma3_i", "", "", tr->gamma3_i, &tr->gamma3_i,NULL)); 540 PetscCall(PetscOptionsReal("-tao_ntr_gamma4_i", "", "", tr->gamma4_i, &tr->gamma4_i,NULL)); 541 PetscCall(PetscOptionsReal("-tao_ntr_theta_i", "", "", tr->theta_i, &tr->theta_i,NULL)); 542 PetscCall(PetscOptionsReal("-tao_ntr_min_radius", "lower bound on initial trust-region radius", "", tr->min_radius, &tr->min_radius,NULL)); 543 PetscCall(PetscOptionsReal("-tao_ntr_max_radius", "upper bound on trust-region radius", "", tr->max_radius, &tr->max_radius,NULL)); 544 PetscCall(PetscOptionsReal("-tao_ntr_epsilon", "tolerance used when computing actual and predicted reduction", "", tr->epsilon, &tr->epsilon,NULL)); 545 PetscOptionsHeadEnd(); 546 PetscCall(KSPSetFromOptions(tao->ksp)); 547 PetscFunctionReturn(0); 548 } 549 550 /*------------------------------------------------------------*/ 551 /*MC 552 TAONTR - Newton's method with trust region for unconstrained minimization. 553 At each iteration, the Newton trust region method solves the system. 554 NTR expects a KSP solver with a trust region radius. 555 min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k 556 557 Options Database Keys: 558 + -tao_ntr_init_type - "constant","direction","interpolation" 559 . -tao_ntr_update_type - "reduction","interpolation" 560 . -tao_ntr_min_radius - lower bound on trust region radius 561 . -tao_ntr_max_radius - upper bound on trust region radius 562 . -tao_ntr_epsilon - tolerance for accepting actual / predicted reduction 563 . -tao_ntr_mu1_i - mu1 interpolation init factor 564 . -tao_ntr_mu2_i - mu2 interpolation init factor 565 . -tao_ntr_gamma1_i - gamma1 interpolation init factor 566 . -tao_ntr_gamma2_i - gamma2 interpolation init factor 567 . -tao_ntr_gamma3_i - gamma3 interpolation init factor 568 . -tao_ntr_gamma4_i - gamma4 interpolation init factor 569 . -tao_ntr_theta_i - theta1 interpolation init factor 570 . -tao_ntr_eta1 - eta1 reduction update factor 571 . -tao_ntr_eta2 - eta2 reduction update factor 572 . -tao_ntr_eta3 - eta3 reduction update factor 573 . -tao_ntr_eta4 - eta4 reduction update factor 574 . -tao_ntr_alpha1 - alpha1 reduction update factor 575 . -tao_ntr_alpha2 - alpha2 reduction update factor 576 . -tao_ntr_alpha3 - alpha3 reduction update factor 577 . -tao_ntr_alpha4 - alpha4 reduction update factor 578 . -tao_ntr_alpha4 - alpha4 reduction update factor 579 . -tao_ntr_mu1 - mu1 interpolation update 580 . -tao_ntr_mu2 - mu2 interpolation update 581 . -tao_ntr_gamma1 - gamma1 interpolcation update 582 . -tao_ntr_gamma2 - gamma2 interpolcation update 583 . -tao_ntr_gamma3 - gamma3 interpolcation update 584 . -tao_ntr_gamma4 - gamma4 interpolation update 585 - -tao_ntr_theta - theta interpolation update 586 587 Level: beginner 588 M*/ 589 590 PETSC_EXTERN PetscErrorCode TaoCreate_NTR(Tao tao) 591 { 592 TAO_NTR *tr; 593 594 PetscFunctionBegin; 595 596 PetscCall(PetscNewLog(tao,&tr)); 597 598 tao->ops->setup = TaoSetUp_NTR; 599 tao->ops->solve = TaoSolve_NTR; 600 tao->ops->setfromoptions = TaoSetFromOptions_NTR; 601 tao->ops->destroy = TaoDestroy_NTR; 602 603 /* Override default settings (unless already changed) */ 604 if (!tao->max_it_changed) tao->max_it = 50; 605 if (!tao->trust0_changed) tao->trust0 = 100.0; 606 tao->data = (void*)tr; 607 608 /* Standard trust region update parameters */ 609 tr->eta1 = 1.0e-4; 610 tr->eta2 = 0.25; 611 tr->eta3 = 0.50; 612 tr->eta4 = 0.90; 613 614 tr->alpha1 = 0.25; 615 tr->alpha2 = 0.50; 616 tr->alpha3 = 1.00; 617 tr->alpha4 = 2.00; 618 tr->alpha5 = 4.00; 619 620 /* Interpolation trust region update parameters */ 621 tr->mu1 = 0.10; 622 tr->mu2 = 0.50; 623 624 tr->gamma1 = 0.25; 625 tr->gamma2 = 0.50; 626 tr->gamma3 = 2.00; 627 tr->gamma4 = 4.00; 628 629 tr->theta = 0.05; 630 631 /* Interpolation parameters for initialization */ 632 tr->mu1_i = 0.35; 633 tr->mu2_i = 0.50; 634 635 tr->gamma1_i = 0.0625; 636 tr->gamma2_i = 0.50; 637 tr->gamma3_i = 2.00; 638 tr->gamma4_i = 5.00; 639 640 tr->theta_i = 0.25; 641 642 tr->min_radius = 1.0e-10; 643 tr->max_radius = 1.0e10; 644 tr->epsilon = 1.0e-6; 645 646 tr->init_type = NTR_INIT_INTERPOLATION; 647 tr->update_type = NTR_UPDATE_REDUCTION; 648 649 /* Set linear solver to default for trust region */ 650 PetscCall(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 651 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1)); 652 PetscCall(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 653 PetscCall(KSPAppendOptionsPrefix(tao->ksp,"tao_ntr_")); 654 PetscCall(KSPSetType(tao->ksp,KSPSTCG)); 655 PetscFunctionReturn(0); 656 } 657