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