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