1 #include <../src/tao/unconstrained/impls/ntl/ntlimpl.h> 2 3 #include <petscksp.h> 4 5 #define NTL_INIT_CONSTANT 0 6 #define NTL_INIT_DIRECTION 1 7 #define NTL_INIT_INTERPOLATION 2 8 #define NTL_INIT_TYPES 3 9 10 #define NTL_UPDATE_REDUCTION 0 11 #define NTL_UPDATE_INTERPOLATION 1 12 #define NTL_UPDATE_TYPES 2 13 14 static const char *NTL_INIT[64] = {"constant","direction","interpolation"}; 15 16 static const char *NTL_UPDATE[64] = {"reduction","interpolation"}; 17 18 /* Implements Newton's Method with a trust-region, line-search approach for 19 solving unconstrained minimization problems. A More'-Thuente line search 20 is used to guarantee that the bfgs preconditioner remains positive 21 definite. */ 22 23 #define NTL_NEWTON 0 24 #define NTL_BFGS 1 25 #define NTL_SCALED_GRADIENT 2 26 #define NTL_GRADIENT 3 27 28 static PetscErrorCode TaoSolve_NTL(Tao tao) 29 { 30 TAO_NTL *tl = (TAO_NTL *)tao->data; 31 KSPType ksp_type; 32 PetscBool is_nash,is_stcg,is_gltr,is_bfgs,is_jacobi,is_symmetric,sym_set; 33 KSPConvergedReason ksp_reason; 34 PC pc; 35 TaoLineSearchConvergedReason ls_reason; 36 37 PetscReal fmin, ftrial, prered, actred, kappa, sigma; 38 PetscReal tau, tau_1, tau_2, tau_max, tau_min, max_radius; 39 PetscReal f, fold, gdx, gnorm; 40 PetscReal step = 1.0; 41 42 PetscReal norm_d = 0.0; 43 PetscInt stepType; 44 PetscInt its; 45 46 PetscInt bfgsUpdates = 0; 47 PetscInt needH; 48 49 PetscInt i_max = 5; 50 PetscInt j_max = 1; 51 PetscInt i, j, n, N; 52 53 PetscInt tr_reject; 54 55 PetscFunctionBegin; 56 if (tao->XL || tao->XU || tao->ops->computebounds) { 57 PetscCall(PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by ntl algorithm\n")); 58 } 59 60 PetscCall(KSPGetType(tao->ksp,&ksp_type)); 61 PetscCall(PetscStrcmp(ksp_type,KSPNASH,&is_nash)); 62 PetscCall(PetscStrcmp(ksp_type,KSPSTCG,&is_stcg)); 63 PetscCall(PetscStrcmp(ksp_type,KSPGLTR,&is_gltr)); 64 PetscCheck(is_nash || is_stcg || is_gltr,PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"TAO_NTR requires nash, stcg, or gltr for the KSP"); 65 66 /* Initialize the radius and modify if it is too large or small */ 67 tao->trust = tao->trust0; 68 tao->trust = PetscMax(tao->trust, tl->min_radius); 69 tao->trust = PetscMin(tao->trust, tl->max_radius); 70 71 /* Allocate the vectors needed for the BFGS approximation */ 72 PetscCall(KSPGetPC(tao->ksp, &pc)); 73 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCLMVM, &is_bfgs)); 74 PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCJACOBI, &is_jacobi)); 75 if (is_bfgs) { 76 tl->bfgs_pre = pc; 77 PetscCall(PCLMVMGetMatLMVM(tl->bfgs_pre, &tl->M)); 78 PetscCall(VecGetLocalSize(tao->solution, &n)); 79 PetscCall(VecGetSize(tao->solution, &N)); 80 PetscCall(MatSetSizes(tl->M, n, n, N, N)); 81 PetscCall(MatLMVMAllocate(tl->M, tao->solution, tao->gradient)); 82 PetscCall(MatIsSymmetricKnown(tl->M, &sym_set, &is_symmetric)); 83 PetscCheck(sym_set && is_symmetric,PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix in the LMVM preconditioner must be symmetric."); 84 } else if (is_jacobi) { 85 PetscCall(PCJacobiSetUseAbs(pc,PETSC_TRUE)); 86 } 87 88 /* Check convergence criteria */ 89 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient)); 90 PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm)); 91 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 92 needH = 1; 93 94 tao->reason = TAO_CONTINUE_ITERATING; 95 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 96 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,step)); 97 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 98 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 99 100 /* Initialize trust-region radius */ 101 switch(tl->init_type) { 102 case NTL_INIT_CONSTANT: 103 /* Use the initial radius specified */ 104 break; 105 106 case NTL_INIT_INTERPOLATION: 107 /* Use the initial radius specified */ 108 max_radius = 0.0; 109 110 for (j = 0; j < j_max; ++j) { 111 fmin = f; 112 sigma = 0.0; 113 114 if (needH) { 115 PetscCall(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 116 needH = 0; 117 } 118 119 for (i = 0; i < i_max; ++i) { 120 PetscCall(VecCopy(tao->solution, tl->W)); 121 PetscCall(VecAXPY(tl->W, -tao->trust/gnorm, tao->gradient)); 122 123 PetscCall(TaoComputeObjective(tao, tl->W, &ftrial)); 124 if (PetscIsInfOrNanReal(ftrial)) { 125 tau = tl->gamma1_i; 126 } else { 127 if (ftrial < fmin) { 128 fmin = ftrial; 129 sigma = -tao->trust / gnorm; 130 } 131 132 PetscCall(MatMult(tao->hessian, tao->gradient, tao->stepdirection)); 133 PetscCall(VecDot(tao->gradient, tao->stepdirection, &prered)); 134 135 prered = tao->trust * (gnorm - 0.5 * tao->trust * prered / (gnorm * gnorm)); 136 actred = f - ftrial; 137 if ((PetscAbsScalar(actred) <= tl->epsilon) && (PetscAbsScalar(prered) <= tl->epsilon)) { 138 kappa = 1.0; 139 } else { 140 kappa = actred / prered; 141 } 142 143 tau_1 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust + (1.0 - tl->theta_i) * prered - actred); 144 tau_2 = tl->theta_i * gnorm * tao->trust / (tl->theta_i * gnorm * tao->trust - (1.0 + tl->theta_i) * prered + actred); 145 tau_min = PetscMin(tau_1, tau_2); 146 tau_max = PetscMax(tau_1, tau_2); 147 148 if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu1_i) { 149 /* Great agreement */ 150 max_radius = PetscMax(max_radius, tao->trust); 151 152 if (tau_max < 1.0) { 153 tau = tl->gamma3_i; 154 } else if (tau_max > tl->gamma4_i) { 155 tau = tl->gamma4_i; 156 } else if (tau_1 >= 1.0 && tau_1 <= tl->gamma4_i && tau_2 < 1.0) { 157 tau = tau_1; 158 } else if (tau_2 >= 1.0 && tau_2 <= tl->gamma4_i && tau_1 < 1.0) { 159 tau = tau_2; 160 } else { 161 tau = tau_max; 162 } 163 } else if (PetscAbsScalar(kappa - (PetscReal)1.0) <= tl->mu2_i) { 164 /* Good agreement */ 165 max_radius = PetscMax(max_radius, tao->trust); 166 167 if (tau_max < tl->gamma2_i) { 168 tau = tl->gamma2_i; 169 } else if (tau_max > tl->gamma3_i) { 170 tau = tl->gamma3_i; 171 } else { 172 tau = tau_max; 173 } 174 } else { 175 /* Not good agreement */ 176 if (tau_min > 1.0) { 177 tau = tl->gamma2_i; 178 } else if (tau_max < tl->gamma1_i) { 179 tau = tl->gamma1_i; 180 } else if ((tau_min < tl->gamma1_i) && (tau_max >= 1.0)) { 181 tau = tl->gamma1_i; 182 } else if ((tau_1 >= tl->gamma1_i) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1_i) || (tau_2 >= 1.0))) { 183 tau = tau_1; 184 } else if ((tau_2 >= tl->gamma1_i) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1_i) || (tau_2 >= 1.0))) { 185 tau = tau_2; 186 } else { 187 tau = tau_max; 188 } 189 } 190 } 191 tao->trust = tau * tao->trust; 192 } 193 194 if (fmin < f) { 195 f = fmin; 196 PetscCall(VecAXPY(tao->solution, sigma, tao->gradient)); 197 PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient)); 198 199 PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm)); 200 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 201 needH = 1; 202 203 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 204 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,step)); 205 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 206 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 207 } 208 } 209 tao->trust = PetscMax(tao->trust, max_radius); 210 211 /* Modify the radius if it is too large or small */ 212 tao->trust = PetscMax(tao->trust, tl->min_radius); 213 tao->trust = PetscMin(tao->trust, tl->max_radius); 214 break; 215 216 default: 217 /* Norm of the first direction will initialize radius */ 218 tao->trust = 0.0; 219 break; 220 } 221 222 /* Set counter for gradient/reset steps */ 223 tl->ntrust = 0; 224 tl->newt = 0; 225 tl->bfgs = 0; 226 tl->grad = 0; 227 228 /* Have not converged; continue with Newton method */ 229 while (tao->reason == TAO_CONTINUE_ITERATING) { 230 /* Call general purpose update function */ 231 if (tao->ops->update) { 232 PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 233 } 234 ++tao->niter; 235 tao->ksp_its=0; 236 /* Compute the Hessian */ 237 if (needH) { 238 PetscCall(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 239 } 240 241 if (tl->bfgs_pre) { 242 /* Update the limited memory preconditioner */ 243 PetscCall(MatLMVMUpdate(tl->M,tao->solution, tao->gradient)); 244 ++bfgsUpdates; 245 } 246 PetscCall(KSPSetOperators(tao->ksp, tao->hessian, tao->hessian_pre)); 247 /* Solve the Newton system of equations */ 248 PetscCall(KSPCGSetRadius(tao->ksp,tl->max_radius)); 249 PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection)); 250 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 251 tao->ksp_its+=its; 252 tao->ksp_tot_its+=its; 253 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 254 255 if (0.0 == tao->trust) { 256 /* Radius was uninitialized; use the norm of the direction */ 257 if (norm_d > 0.0) { 258 tao->trust = norm_d; 259 260 /* Modify the radius if it is too large or small */ 261 tao->trust = PetscMax(tao->trust, tl->min_radius); 262 tao->trust = PetscMin(tao->trust, tl->max_radius); 263 } else { 264 /* The direction was bad; set radius to default value and re-solve 265 the trust-region subproblem to get a direction */ 266 tao->trust = tao->trust0; 267 268 /* Modify the radius if it is too large or small */ 269 tao->trust = PetscMax(tao->trust, tl->min_radius); 270 tao->trust = PetscMin(tao->trust, tl->max_radius); 271 272 PetscCall(KSPCGSetRadius(tao->ksp,tl->max_radius)); 273 PetscCall(KSPSolve(tao->ksp, tao->gradient, tao->stepdirection)); 274 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 275 tao->ksp_its+=its; 276 tao->ksp_tot_its+=its; 277 PetscCall(KSPCGGetNormD(tao->ksp, &norm_d)); 278 279 PetscCheck(norm_d != 0.0,PetscObjectComm((PetscObject)tao),PETSC_ERR_PLIB, "Initial direction zero"); 280 } 281 } 282 283 PetscCall(VecScale(tao->stepdirection, -1.0)); 284 PetscCall(KSPGetConvergedReason(tao->ksp, &ksp_reason)); 285 if ((KSP_DIVERGED_INDEFINITE_PC == ksp_reason) && (tl->bfgs_pre)) { 286 /* Preconditioner is numerically indefinite; reset the 287 approximate if using BFGS preconditioning. */ 288 PetscCall(MatLMVMReset(tl->M, PETSC_FALSE)); 289 PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient)); 290 bfgsUpdates = 1; 291 } 292 293 /* Check trust-region reduction conditions */ 294 tr_reject = 0; 295 if (NTL_UPDATE_REDUCTION == tl->update_type) { 296 /* Get predicted reduction */ 297 PetscCall(KSPCGGetObjFcn(tao->ksp,&prered)); 298 if (prered >= 0.0) { 299 /* The predicted reduction has the wrong sign. This cannot 300 happen in infinite precision arithmetic. Step should 301 be rejected! */ 302 tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d); 303 tr_reject = 1; 304 } else { 305 /* Compute trial step and function value */ 306 PetscCall(VecCopy(tao->solution, tl->W)); 307 PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection)); 308 PetscCall(TaoComputeObjective(tao, tl->W, &ftrial)); 309 310 if (PetscIsInfOrNanReal(ftrial)) { 311 tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d); 312 tr_reject = 1; 313 } else { 314 /* Compute and actual reduction */ 315 actred = f - ftrial; 316 prered = -prered; 317 if ((PetscAbsScalar(actred) <= tl->epsilon) && 318 (PetscAbsScalar(prered) <= tl->epsilon)) { 319 kappa = 1.0; 320 } else { 321 kappa = actred / prered; 322 } 323 324 /* Accept of reject the step and update radius */ 325 if (kappa < tl->eta1) { 326 /* Reject the step */ 327 tao->trust = tl->alpha1 * PetscMin(tao->trust, norm_d); 328 tr_reject = 1; 329 } else { 330 /* Accept the step */ 331 if (kappa < tl->eta2) { 332 /* Marginal bad step */ 333 tao->trust = tl->alpha2 * PetscMin(tao->trust, norm_d); 334 } else if (kappa < tl->eta3) { 335 /* Reasonable step */ 336 tao->trust = tl->alpha3 * tao->trust; 337 } else if (kappa < tl->eta4) { 338 /* Good step */ 339 tao->trust = PetscMax(tl->alpha4 * norm_d, tao->trust); 340 } else { 341 /* Very good step */ 342 tao->trust = PetscMax(tl->alpha5 * norm_d, tao->trust); 343 } 344 } 345 } 346 } 347 } else { 348 /* Get predicted reduction */ 349 PetscCall(KSPCGGetObjFcn(tao->ksp,&prered)); 350 if (prered >= 0.0) { 351 /* The predicted reduction has the wrong sign. This cannot 352 happen in infinite precision arithmetic. Step should 353 be rejected! */ 354 tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d); 355 tr_reject = 1; 356 } else { 357 PetscCall(VecCopy(tao->solution, tl->W)); 358 PetscCall(VecAXPY(tl->W, 1.0, tao->stepdirection)); 359 PetscCall(TaoComputeObjective(tao, tl->W, &ftrial)); 360 if (PetscIsInfOrNanReal(ftrial)) { 361 tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d); 362 tr_reject = 1; 363 } else { 364 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 365 366 actred = f - ftrial; 367 prered = -prered; 368 if ((PetscAbsScalar(actred) <= tl->epsilon) && 369 (PetscAbsScalar(prered) <= tl->epsilon)) { 370 kappa = 1.0; 371 } else { 372 kappa = actred / prered; 373 } 374 375 tau_1 = tl->theta * gdx / (tl->theta * gdx - (1.0 - tl->theta) * prered + actred); 376 tau_2 = tl->theta * gdx / (tl->theta * gdx + (1.0 + tl->theta) * prered - actred); 377 tau_min = PetscMin(tau_1, tau_2); 378 tau_max = PetscMax(tau_1, tau_2); 379 380 if (kappa >= 1.0 - tl->mu1) { 381 /* Great agreement; accept step and update radius */ 382 if (tau_max < 1.0) { 383 tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d); 384 } else if (tau_max > tl->gamma4) { 385 tao->trust = PetscMax(tao->trust, tl->gamma4 * norm_d); 386 } else { 387 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 388 } 389 } else if (kappa >= 1.0 - tl->mu2) { 390 /* Good agreement */ 391 392 if (tau_max < tl->gamma2) { 393 tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d); 394 } else if (tau_max > tl->gamma3) { 395 tao->trust = PetscMax(tao->trust, tl->gamma3 * norm_d); 396 } else if (tau_max < 1.0) { 397 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 398 } else { 399 tao->trust = PetscMax(tao->trust, tau_max * norm_d); 400 } 401 } else { 402 /* Not good agreement */ 403 if (tau_min > 1.0) { 404 tao->trust = tl->gamma2 * PetscMin(tao->trust, norm_d); 405 } else if (tau_max < tl->gamma1) { 406 tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d); 407 } else if ((tau_min < tl->gamma1) && (tau_max >= 1.0)) { 408 tao->trust = tl->gamma1 * PetscMin(tao->trust, norm_d); 409 } else if ((tau_1 >= tl->gamma1) && (tau_1 < 1.0) && ((tau_2 < tl->gamma1) || (tau_2 >= 1.0))) { 410 tao->trust = tau_1 * PetscMin(tao->trust, norm_d); 411 } else if ((tau_2 >= tl->gamma1) && (tau_2 < 1.0) && ((tau_1 < tl->gamma1) || (tau_2 >= 1.0))) { 412 tao->trust = tau_2 * PetscMin(tao->trust, norm_d); 413 } else { 414 tao->trust = tau_max * PetscMin(tao->trust, norm_d); 415 } 416 tr_reject = 1; 417 } 418 } 419 } 420 } 421 422 if (tr_reject) { 423 /* The trust-region constraints rejected the step. Apply a linesearch. 424 Check for descent direction. */ 425 PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx)); 426 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 427 /* Newton step is not descent or direction produced Inf or NaN */ 428 429 if (!tl->bfgs_pre) { 430 /* We don't have the bfgs matrix around and updated 431 Must use gradient direction in this case */ 432 PetscCall(VecCopy(tao->gradient, tao->stepdirection)); 433 PetscCall(VecScale(tao->stepdirection, -1.0)); 434 ++tl->grad; 435 stepType = NTL_GRADIENT; 436 } else { 437 /* Attempt to use the BFGS direction */ 438 PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection)); 439 PetscCall(VecScale(tao->stepdirection, -1.0)); 440 441 /* Check for success (descent direction) */ 442 PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx)); 443 if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) { 444 /* BFGS direction is not descent or direction produced not a number 445 We can assert bfgsUpdates > 1 in this case because 446 the first solve produces the scaled gradient direction, 447 which is guaranteed to be descent */ 448 449 /* Use steepest descent direction (scaled) */ 450 PetscCall(MatLMVMReset(tl->M, PETSC_FALSE)); 451 PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient)); 452 PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection)); 453 PetscCall(VecScale(tao->stepdirection, -1.0)); 454 455 bfgsUpdates = 1; 456 ++tl->grad; 457 stepType = NTL_GRADIENT; 458 } else { 459 PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates)); 460 if (1 == bfgsUpdates) { 461 /* The first BFGS direction is always the scaled gradient */ 462 ++tl->grad; 463 stepType = NTL_GRADIENT; 464 } else { 465 ++tl->bfgs; 466 stepType = NTL_BFGS; 467 } 468 } 469 } 470 } else { 471 /* Computed Newton step is descent */ 472 ++tl->newt; 473 stepType = NTL_NEWTON; 474 } 475 476 /* Perform the linesearch */ 477 fold = f; 478 PetscCall(VecCopy(tao->solution, tl->Xold)); 479 PetscCall(VecCopy(tao->gradient, tl->Gold)); 480 481 step = 1.0; 482 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason)); 483 PetscCall(TaoAddLineSearchCounts(tao)); 484 485 while (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER && stepType != NTL_GRADIENT) { /* Linesearch failed */ 486 /* Linesearch failed */ 487 f = fold; 488 PetscCall(VecCopy(tl->Xold, tao->solution)); 489 PetscCall(VecCopy(tl->Gold, tao->gradient)); 490 491 switch(stepType) { 492 case NTL_NEWTON: 493 /* Failed to obtain acceptable iterate with Newton step */ 494 495 if (tl->bfgs_pre) { 496 /* We don't have the bfgs matrix around and being updated 497 Must use gradient direction in this case */ 498 PetscCall(VecCopy(tao->gradient, tao->stepdirection)); 499 ++tl->grad; 500 stepType = NTL_GRADIENT; 501 } else { 502 /* Attempt to use the BFGS direction */ 503 PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection)); 504 505 /* Check for success (descent direction) */ 506 PetscCall(VecDot(tao->stepdirection, tao->gradient, &gdx)); 507 if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) { 508 /* BFGS direction is not descent or direction produced 509 not a number. We can assert bfgsUpdates > 1 in this case 510 Use steepest descent direction (scaled) */ 511 PetscCall(MatLMVMReset(tl->M, PETSC_FALSE)); 512 PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient)); 513 PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection)); 514 515 bfgsUpdates = 1; 516 ++tl->grad; 517 stepType = NTL_GRADIENT; 518 } else { 519 PetscCall(MatLMVMGetUpdateCount(tl->M, &bfgsUpdates)); 520 if (1 == bfgsUpdates) { 521 /* The first BFGS direction is always the scaled gradient */ 522 ++tl->grad; 523 stepType = NTL_GRADIENT; 524 } else { 525 ++tl->bfgs; 526 stepType = NTL_BFGS; 527 } 528 } 529 } 530 break; 531 532 case NTL_BFGS: 533 /* Can only enter if pc_type == NTL_PC_BFGS 534 Failed to obtain acceptable iterate with BFGS step 535 Attempt to use the scaled gradient direction */ 536 PetscCall(MatLMVMReset(tl->M, PETSC_FALSE)); 537 PetscCall(MatLMVMUpdate(tl->M, tao->solution, tao->gradient)); 538 PetscCall(MatSolve(tl->M, tao->gradient, tao->stepdirection)); 539 540 bfgsUpdates = 1; 541 ++tl->grad; 542 stepType = NTL_GRADIENT; 543 break; 544 } 545 PetscCall(VecScale(tao->stepdirection, -1.0)); 546 547 /* This may be incorrect; linesearch has values for stepmax and stepmin 548 that should be reset. */ 549 step = 1.0; 550 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_reason)); 551 PetscCall(TaoAddLineSearchCounts(tao)); 552 } 553 554 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 555 /* Failed to find an improving point */ 556 f = fold; 557 PetscCall(VecCopy(tl->Xold, tao->solution)); 558 PetscCall(VecCopy(tl->Gold, tao->gradient)); 559 tao->trust = 0.0; 560 step = 0.0; 561 tao->reason = TAO_DIVERGED_LS_FAILURE; 562 break; 563 } else if (stepType == NTL_NEWTON) { 564 if (step < tl->nu1) { 565 /* Very bad step taken; reduce radius */ 566 tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust); 567 } else if (step < tl->nu2) { 568 /* Reasonably bad step taken; reduce radius */ 569 tao->trust = tl->omega2 * PetscMin(norm_d, tao->trust); 570 } else if (step < tl->nu3) { 571 /* Reasonable step was taken; leave radius alone */ 572 if (tl->omega3 < 1.0) { 573 tao->trust = tl->omega3 * PetscMin(norm_d, tao->trust); 574 } else if (tl->omega3 > 1.0) { 575 tao->trust = PetscMax(tl->omega3 * norm_d, tao->trust); 576 } 577 } else if (step < tl->nu4) { 578 /* Full step taken; increase the radius */ 579 tao->trust = PetscMax(tl->omega4 * norm_d, tao->trust); 580 } else { 581 /* More than full step taken; increase the radius */ 582 tao->trust = PetscMax(tl->omega5 * norm_d, tao->trust); 583 } 584 } else { 585 /* Newton step was not good; reduce the radius */ 586 tao->trust = tl->omega1 * PetscMin(norm_d, tao->trust); 587 } 588 } else { 589 /* Trust-region step is accepted */ 590 PetscCall(VecCopy(tl->W, tao->solution)); 591 f = ftrial; 592 PetscCall(TaoComputeGradient(tao, tao->solution, tao->gradient)); 593 ++tl->ntrust; 594 } 595 596 /* The radius may have been increased; modify if it is too large */ 597 tao->trust = PetscMin(tao->trust, tl->max_radius); 598 599 /* Check for converged */ 600 PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm)); 601 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER,"User provided compute function generated Not-a-Number"); 602 needH = 1; 603 604 PetscCall(TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its)); 605 PetscCall(TaoMonitor(tao,tao->niter,f,gnorm,0.0,step)); 606 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 607 } 608 PetscFunctionReturn(0); 609 } 610 611 /* ---------------------------------------------------------- */ 612 static PetscErrorCode TaoSetUp_NTL(Tao tao) 613 { 614 TAO_NTL *tl = (TAO_NTL *)tao->data; 615 616 PetscFunctionBegin; 617 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 618 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 619 if (!tl->W) PetscCall(VecDuplicate(tao->solution, &tl->W)); 620 if (!tl->Xold) PetscCall(VecDuplicate(tao->solution, &tl->Xold)); 621 if (!tl->Gold) PetscCall(VecDuplicate(tao->solution, &tl->Gold)); 622 tl->bfgs_pre = NULL; 623 tl->M = NULL; 624 PetscFunctionReturn(0); 625 } 626 627 /*------------------------------------------------------------*/ 628 static PetscErrorCode TaoDestroy_NTL(Tao tao) 629 { 630 TAO_NTL *tl = (TAO_NTL *)tao->data; 631 632 PetscFunctionBegin; 633 if (tao->setupcalled) { 634 PetscCall(VecDestroy(&tl->W)); 635 PetscCall(VecDestroy(&tl->Xold)); 636 PetscCall(VecDestroy(&tl->Gold)); 637 } 638 PetscCall(PetscFree(tao->data)); 639 PetscFunctionReturn(0); 640 } 641 642 /*------------------------------------------------------------*/ 643 static PetscErrorCode TaoSetFromOptions_NTL(PetscOptionItems *PetscOptionsObject,Tao tao) 644 { 645 TAO_NTL *tl = (TAO_NTL *)tao->data; 646 647 PetscFunctionBegin; 648 PetscOptionsHeadBegin(PetscOptionsObject,"Newton trust region with line search method for unconstrained optimization"); 649 PetscCall(PetscOptionsEList("-tao_ntl_init_type", "radius initialization type", "", NTL_INIT, NTL_INIT_TYPES, NTL_INIT[tl->init_type], &tl->init_type,NULL)); 650 PetscCall(PetscOptionsEList("-tao_ntl_update_type", "radius update type", "", NTL_UPDATE, NTL_UPDATE_TYPES, NTL_UPDATE[tl->update_type], &tl->update_type,NULL)); 651 PetscCall(PetscOptionsReal("-tao_ntl_eta1", "poor steplength; reduce radius", "", tl->eta1, &tl->eta1,NULL)); 652 PetscCall(PetscOptionsReal("-tao_ntl_eta2", "reasonable steplength; leave radius alone", "", tl->eta2, &tl->eta2,NULL)); 653 PetscCall(PetscOptionsReal("-tao_ntl_eta3", "good steplength; increase radius", "", tl->eta3, &tl->eta3,NULL)); 654 PetscCall(PetscOptionsReal("-tao_ntl_eta4", "excellent steplength; greatly increase radius", "", tl->eta4, &tl->eta4,NULL)); 655 PetscCall(PetscOptionsReal("-tao_ntl_alpha1", "", "", tl->alpha1, &tl->alpha1,NULL)); 656 PetscCall(PetscOptionsReal("-tao_ntl_alpha2", "", "", tl->alpha2, &tl->alpha2,NULL)); 657 PetscCall(PetscOptionsReal("-tao_ntl_alpha3", "", "", tl->alpha3, &tl->alpha3,NULL)); 658 PetscCall(PetscOptionsReal("-tao_ntl_alpha4", "", "", tl->alpha4, &tl->alpha4,NULL)); 659 PetscCall(PetscOptionsReal("-tao_ntl_alpha5", "", "", tl->alpha5, &tl->alpha5,NULL)); 660 PetscCall(PetscOptionsReal("-tao_ntl_nu1", "poor steplength; reduce radius", "", tl->nu1, &tl->nu1,NULL)); 661 PetscCall(PetscOptionsReal("-tao_ntl_nu2", "reasonable steplength; leave radius alone", "", tl->nu2, &tl->nu2,NULL)); 662 PetscCall(PetscOptionsReal("-tao_ntl_nu3", "good steplength; increase radius", "", tl->nu3, &tl->nu3,NULL)); 663 PetscCall(PetscOptionsReal("-tao_ntl_nu4", "excellent steplength; greatly increase radius", "", tl->nu4, &tl->nu4,NULL)); 664 PetscCall(PetscOptionsReal("-tao_ntl_omega1", "", "", tl->omega1, &tl->omega1,NULL)); 665 PetscCall(PetscOptionsReal("-tao_ntl_omega2", "", "", tl->omega2, &tl->omega2,NULL)); 666 PetscCall(PetscOptionsReal("-tao_ntl_omega3", "", "", tl->omega3, &tl->omega3,NULL)); 667 PetscCall(PetscOptionsReal("-tao_ntl_omega4", "", "", tl->omega4, &tl->omega4,NULL)); 668 PetscCall(PetscOptionsReal("-tao_ntl_omega5", "", "", tl->omega5, &tl->omega5,NULL)); 669 PetscCall(PetscOptionsReal("-tao_ntl_mu1_i", "", "", tl->mu1_i, &tl->mu1_i,NULL)); 670 PetscCall(PetscOptionsReal("-tao_ntl_mu2_i", "", "", tl->mu2_i, &tl->mu2_i,NULL)); 671 PetscCall(PetscOptionsReal("-tao_ntl_gamma1_i", "", "", tl->gamma1_i, &tl->gamma1_i,NULL)); 672 PetscCall(PetscOptionsReal("-tao_ntl_gamma2_i", "", "", tl->gamma2_i, &tl->gamma2_i,NULL)); 673 PetscCall(PetscOptionsReal("-tao_ntl_gamma3_i", "", "", tl->gamma3_i, &tl->gamma3_i,NULL)); 674 PetscCall(PetscOptionsReal("-tao_ntl_gamma4_i", "", "", tl->gamma4_i, &tl->gamma4_i,NULL)); 675 PetscCall(PetscOptionsReal("-tao_ntl_theta_i", "", "", tl->theta_i, &tl->theta_i,NULL)); 676 PetscCall(PetscOptionsReal("-tao_ntl_mu1", "", "", tl->mu1, &tl->mu1,NULL)); 677 PetscCall(PetscOptionsReal("-tao_ntl_mu2", "", "", tl->mu2, &tl->mu2,NULL)); 678 PetscCall(PetscOptionsReal("-tao_ntl_gamma1", "", "", tl->gamma1, &tl->gamma1,NULL)); 679 PetscCall(PetscOptionsReal("-tao_ntl_gamma2", "", "", tl->gamma2, &tl->gamma2,NULL)); 680 PetscCall(PetscOptionsReal("-tao_ntl_gamma3", "", "", tl->gamma3, &tl->gamma3,NULL)); 681 PetscCall(PetscOptionsReal("-tao_ntl_gamma4", "", "", tl->gamma4, &tl->gamma4,NULL)); 682 PetscCall(PetscOptionsReal("-tao_ntl_theta", "", "", tl->theta, &tl->theta,NULL)); 683 PetscCall(PetscOptionsReal("-tao_ntl_min_radius", "lower bound on initial radius", "", tl->min_radius, &tl->min_radius,NULL)); 684 PetscCall(PetscOptionsReal("-tao_ntl_max_radius", "upper bound on radius", "", tl->max_radius, &tl->max_radius,NULL)); 685 PetscCall(PetscOptionsReal("-tao_ntl_epsilon", "tolerance used when computing actual and predicted reduction", "", tl->epsilon, &tl->epsilon,NULL)); 686 PetscOptionsHeadEnd(); 687 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 688 PetscCall(KSPSetFromOptions(tao->ksp)); 689 PetscFunctionReturn(0); 690 } 691 692 /*------------------------------------------------------------*/ 693 static PetscErrorCode TaoView_NTL(Tao tao, PetscViewer viewer) 694 { 695 TAO_NTL *tl = (TAO_NTL *)tao->data; 696 PetscBool isascii; 697 698 PetscFunctionBegin; 699 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 700 if (isascii) { 701 PetscCall(PetscViewerASCIIPushTab(viewer)); 702 PetscCall(PetscViewerASCIIPrintf(viewer, "Trust-region steps: %" PetscInt_FMT "\n", tl->ntrust)); 703 PetscCall(PetscViewerASCIIPrintf(viewer, "Newton search steps: %" PetscInt_FMT "\n", tl->newt)); 704 PetscCall(PetscViewerASCIIPrintf(viewer, "BFGS search steps: %" PetscInt_FMT "\n", tl->bfgs)); 705 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient search steps: %" PetscInt_FMT "\n", tl->grad)); 706 PetscCall(PetscViewerASCIIPopTab(viewer)); 707 } 708 PetscFunctionReturn(0); 709 } 710 711 /* ---------------------------------------------------------- */ 712 /*MC 713 TAONTL - Newton's method with trust region globalization and line search fallback. 714 At each iteration, the Newton trust region method solves the system for d 715 and performs a line search in the d direction: 716 717 min_d .5 dT Hk d + gkT d, s.t. ||d|| < Delta_k 718 719 Options Database Keys: 720 + -tao_ntl_init_type - "constant","direction","interpolation" 721 . -tao_ntl_update_type - "reduction","interpolation" 722 . -tao_ntl_min_radius - lower bound on trust region radius 723 . -tao_ntl_max_radius - upper bound on trust region radius 724 . -tao_ntl_epsilon - tolerance for accepting actual / predicted reduction 725 . -tao_ntl_mu1_i - mu1 interpolation init factor 726 . -tao_ntl_mu2_i - mu2 interpolation init factor 727 . -tao_ntl_gamma1_i - gamma1 interpolation init factor 728 . -tao_ntl_gamma2_i - gamma2 interpolation init factor 729 . -tao_ntl_gamma3_i - gamma3 interpolation init factor 730 . -tao_ntl_gamma4_i - gamma4 interpolation init factor 731 . -tao_ntl_theta_i - theta1 interpolation init factor 732 . -tao_ntl_eta1 - eta1 reduction update factor 733 . -tao_ntl_eta2 - eta2 reduction update factor 734 . -tao_ntl_eta3 - eta3 reduction update factor 735 . -tao_ntl_eta4 - eta4 reduction update factor 736 . -tao_ntl_alpha1 - alpha1 reduction update factor 737 . -tao_ntl_alpha2 - alpha2 reduction update factor 738 . -tao_ntl_alpha3 - alpha3 reduction update factor 739 . -tao_ntl_alpha4 - alpha4 reduction update factor 740 . -tao_ntl_alpha4 - alpha4 reduction update factor 741 . -tao_ntl_mu1 - mu1 interpolation update 742 . -tao_ntl_mu2 - mu2 interpolation update 743 . -tao_ntl_gamma1 - gamma1 interpolcation update 744 . -tao_ntl_gamma2 - gamma2 interpolcation update 745 . -tao_ntl_gamma3 - gamma3 interpolcation update 746 . -tao_ntl_gamma4 - gamma4 interpolation update 747 - -tao_ntl_theta - theta1 interpolation update 748 749 Level: beginner 750 M*/ 751 PETSC_EXTERN PetscErrorCode TaoCreate_NTL(Tao tao) 752 { 753 TAO_NTL *tl; 754 const char *morethuente_type = TAOLINESEARCHMT; 755 756 PetscFunctionBegin; 757 PetscCall(PetscNewLog(tao,&tl)); 758 tao->ops->setup = TaoSetUp_NTL; 759 tao->ops->solve = TaoSolve_NTL; 760 tao->ops->view = TaoView_NTL; 761 tao->ops->setfromoptions = TaoSetFromOptions_NTL; 762 tao->ops->destroy = TaoDestroy_NTL; 763 764 /* Override default settings (unless already changed) */ 765 if (!tao->max_it_changed) tao->max_it = 50; 766 if (!tao->trust0_changed) tao->trust0 = 100.0; 767 768 tao->data = (void*)tl; 769 770 /* Default values for trust-region radius update based on steplength */ 771 tl->nu1 = 0.25; 772 tl->nu2 = 0.50; 773 tl->nu3 = 1.00; 774 tl->nu4 = 1.25; 775 776 tl->omega1 = 0.25; 777 tl->omega2 = 0.50; 778 tl->omega3 = 1.00; 779 tl->omega4 = 2.00; 780 tl->omega5 = 4.00; 781 782 /* Default values for trust-region radius update based on reduction */ 783 tl->eta1 = 1.0e-4; 784 tl->eta2 = 0.25; 785 tl->eta3 = 0.50; 786 tl->eta4 = 0.90; 787 788 tl->alpha1 = 0.25; 789 tl->alpha2 = 0.50; 790 tl->alpha3 = 1.00; 791 tl->alpha4 = 2.00; 792 tl->alpha5 = 4.00; 793 794 /* Default values for trust-region radius update based on interpolation */ 795 tl->mu1 = 0.10; 796 tl->mu2 = 0.50; 797 798 tl->gamma1 = 0.25; 799 tl->gamma2 = 0.50; 800 tl->gamma3 = 2.00; 801 tl->gamma4 = 4.00; 802 803 tl->theta = 0.05; 804 805 /* Default values for trust region initialization based on interpolation */ 806 tl->mu1_i = 0.35; 807 tl->mu2_i = 0.50; 808 809 tl->gamma1_i = 0.0625; 810 tl->gamma2_i = 0.5; 811 tl->gamma3_i = 2.0; 812 tl->gamma4_i = 5.0; 813 814 tl->theta_i = 0.25; 815 816 /* Remaining parameters */ 817 tl->min_radius = 1.0e-10; 818 tl->max_radius = 1.0e10; 819 tl->epsilon = 1.0e-6; 820 821 tl->init_type = NTL_INIT_INTERPOLATION; 822 tl->update_type = NTL_UPDATE_REDUCTION; 823 824 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 825 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch,(PetscObject)tao,1)); 826 PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type)); 827 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 828 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 829 PetscCall(KSPCreate(((PetscObject)tao)->comm,&tao->ksp)); 830 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp,(PetscObject)tao,1)); 831 PetscCall(KSPSetOptionsPrefix(tao->ksp,tao->hdr.prefix)); 832 PetscCall(KSPAppendOptionsPrefix(tao->ksp,"tao_ntl_")); 833 PetscCall(KSPSetType(tao->ksp,KSPSTCG)); 834 PetscFunctionReturn(0); 835 } 836