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