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