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