1 #include <petsctaolinesearch.h> 2 #include <../src/tao/unconstrained/impls/cg/taocg.h> 3 4 #define CG_FletcherReeves 0 5 #define CG_PolakRibiere 1 6 #define CG_PolakRibierePlus 2 7 #define CG_HestenesStiefel 3 8 #define CG_DaiYuan 4 9 #define CG_Types 5 10 11 static const char *CG_Table[64] = {"fr", "pr", "prp", "hs", "dy"}; 12 13 static PetscErrorCode TaoSolve_CG(Tao tao) 14 { 15 TAO_CG *cgP = (TAO_CG *)tao->data; 16 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 17 PetscReal step = 1.0, f, gnorm, gnorm2, delta, gd, ginner, beta; 18 PetscReal gd_old, gnorm2_old, f_old; 19 20 PetscFunctionBegin; 21 if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by cg algorithm\n")); 22 23 /* Check convergence criteria */ 24 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient)); 25 PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm)); 26 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN"); 27 28 tao->reason = TAO_CONTINUE_ITERATING; 29 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 30 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step)); 31 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 32 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS); 33 34 /* Set initial direction to -gradient */ 35 PetscCall(VecCopy(tao->gradient, tao->stepdirection)); 36 PetscCall(VecScale(tao->stepdirection, -1.0)); 37 gnorm2 = gnorm * gnorm; 38 39 /* Set initial scaling for the function */ 40 if (f != 0.0) { 41 delta = 2.0 * PetscAbsScalar(f) / gnorm2; 42 delta = PetscMax(delta, cgP->delta_min); 43 delta = PetscMin(delta, cgP->delta_max); 44 } else { 45 delta = 2.0 / gnorm2; 46 delta = PetscMax(delta, cgP->delta_min); 47 delta = PetscMin(delta, cgP->delta_max); 48 } 49 /* Set counter for gradient and reset steps */ 50 cgP->ngradsteps = 0; 51 cgP->nresetsteps = 0; 52 53 while (1) { 54 /* Call general purpose update function */ 55 if (tao->ops->update) { 56 PetscUseTypeMethod(tao, update, tao->niter, tao->user_update); 57 PetscCall(TaoComputeObjective(tao, tao->solution, &f)); 58 } 59 /* Save the current gradient information */ 60 f_old = f; 61 gnorm2_old = gnorm2; 62 PetscCall(VecCopy(tao->solution, cgP->X_old)); 63 PetscCall(VecCopy(tao->gradient, cgP->G_old)); 64 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd)); 65 if ((gd >= 0) || PetscIsInfOrNanReal(gd)) { 66 ++cgP->ngradsteps; 67 if (f != 0.0) { 68 delta = 2.0 * PetscAbsScalar(f) / gnorm2; 69 delta = PetscMax(delta, cgP->delta_min); 70 delta = PetscMin(delta, cgP->delta_max); 71 } else { 72 delta = 2.0 / gnorm2; 73 delta = PetscMax(delta, cgP->delta_min); 74 delta = PetscMin(delta, cgP->delta_max); 75 } 76 77 PetscCall(VecCopy(tao->gradient, tao->stepdirection)); 78 PetscCall(VecScale(tao->stepdirection, -1.0)); 79 } 80 81 /* Search direction for improving point */ 82 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta)); 83 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status)); 84 PetscCall(TaoAddLineSearchCounts(tao)); 85 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 86 /* Linesearch failed */ 87 /* Reset factors and use scaled gradient step */ 88 ++cgP->nresetsteps; 89 f = f_old; 90 gnorm2 = gnorm2_old; 91 PetscCall(VecCopy(cgP->X_old, tao->solution)); 92 PetscCall(VecCopy(cgP->G_old, tao->gradient)); 93 94 if (f != 0.0) { 95 delta = 2.0 * PetscAbsScalar(f) / gnorm2; 96 delta = PetscMax(delta, cgP->delta_min); 97 delta = PetscMin(delta, cgP->delta_max); 98 } else { 99 delta = 2.0 / gnorm2; 100 delta = PetscMax(delta, cgP->delta_min); 101 delta = PetscMin(delta, cgP->delta_max); 102 } 103 104 PetscCall(VecCopy(tao->gradient, tao->stepdirection)); 105 PetscCall(VecScale(tao->stepdirection, -1.0)); 106 107 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta)); 108 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status)); 109 PetscCall(TaoAddLineSearchCounts(tao)); 110 111 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 112 /* Linesearch failed again */ 113 /* switch to unscaled gradient */ 114 f = f_old; 115 PetscCall(VecCopy(cgP->X_old, tao->solution)); 116 PetscCall(VecCopy(cgP->G_old, tao->gradient)); 117 delta = 1.0; 118 PetscCall(VecCopy(tao->solution, tao->stepdirection)); 119 PetscCall(VecScale(tao->stepdirection, -1.0)); 120 121 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, delta)); 122 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, tao->stepdirection, &step, &ls_status)); 123 PetscCall(TaoAddLineSearchCounts(tao)); 124 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 125 /* Line search failed for last time -- give up */ 126 f = f_old; 127 PetscCall(VecCopy(cgP->X_old, tao->solution)); 128 PetscCall(VecCopy(cgP->G_old, tao->gradient)); 129 step = 0.0; 130 tao->reason = TAO_DIVERGED_LS_FAILURE; 131 } 132 } 133 } 134 135 /* Check for bad value */ 136 PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm)); 137 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User-provided compute function generated infinity or NaN"); 138 139 /* Check for termination */ 140 gnorm2 = gnorm * gnorm; 141 tao->niter++; 142 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 143 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step)); 144 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 145 if (tao->reason != TAO_CONTINUE_ITERATING) break; 146 147 /* Check for restart condition */ 148 PetscCall(VecDot(tao->gradient, cgP->G_old, &ginner)); 149 if (PetscAbsScalar(ginner) >= cgP->eta * gnorm2) { 150 /* Gradients far from orthogonal; use steepest descent direction */ 151 beta = 0.0; 152 } else { 153 /* Gradients close to orthogonal; use conjugate gradient formula */ 154 switch (cgP->cg_type) { 155 case CG_FletcherReeves: 156 beta = gnorm2 / gnorm2_old; 157 break; 158 159 case CG_PolakRibiere: 160 beta = (gnorm2 - ginner) / gnorm2_old; 161 break; 162 163 case CG_PolakRibierePlus: 164 beta = PetscMax((gnorm2 - ginner) / gnorm2_old, 0.0); 165 break; 166 167 case CG_HestenesStiefel: 168 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd)); 169 PetscCall(VecDot(cgP->G_old, tao->stepdirection, &gd_old)); 170 beta = (gnorm2 - ginner) / (gd - gd_old); 171 break; 172 173 case CG_DaiYuan: 174 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd)); 175 PetscCall(VecDot(cgP->G_old, tao->stepdirection, &gd_old)); 176 beta = gnorm2 / (gd - gd_old); 177 break; 178 179 default: 180 beta = 0.0; 181 break; 182 } 183 } 184 185 /* Compute the direction d=-g + beta*d */ 186 PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, tao->gradient)); 187 188 /* update initial steplength choice */ 189 delta = 1.0; 190 delta = PetscMax(delta, cgP->delta_min); 191 delta = PetscMin(delta, cgP->delta_max); 192 } 193 PetscFunctionReturn(PETSC_SUCCESS); 194 } 195 196 static PetscErrorCode TaoSetUp_CG(Tao tao) 197 { 198 TAO_CG *cgP = (TAO_CG *)tao->data; 199 200 PetscFunctionBegin; 201 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 202 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 203 if (!cgP->X_old) PetscCall(VecDuplicate(tao->solution, &cgP->X_old)); 204 if (!cgP->G_old) PetscCall(VecDuplicate(tao->gradient, &cgP->G_old)); 205 PetscFunctionReturn(PETSC_SUCCESS); 206 } 207 208 static PetscErrorCode TaoDestroy_CG(Tao tao) 209 { 210 TAO_CG *cgP = (TAO_CG *)tao->data; 211 212 PetscFunctionBegin; 213 if (tao->setupcalled) { 214 PetscCall(VecDestroy(&cgP->X_old)); 215 PetscCall(VecDestroy(&cgP->G_old)); 216 } 217 PetscCall(TaoLineSearchDestroy(&tao->linesearch)); 218 PetscCall(PetscFree(tao->data)); 219 PetscFunctionReturn(PETSC_SUCCESS); 220 } 221 222 static PetscErrorCode TaoSetFromOptions_CG(Tao tao, PetscOptionItems PetscOptionsObject) 223 { 224 TAO_CG *cgP = (TAO_CG *)tao->data; 225 226 PetscFunctionBegin; 227 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 228 PetscOptionsHeadBegin(PetscOptionsObject, "Nonlinear Conjugate Gradient method for unconstrained optimization"); 229 PetscCall(PetscOptionsReal("-tao_cg_eta", "restart tolerance", "", cgP->eta, &cgP->eta, NULL)); 230 PetscCall(PetscOptionsEList("-tao_cg_type", "cg formula", "", CG_Table, CG_Types, CG_Table[cgP->cg_type], &cgP->cg_type, NULL)); 231 PetscCall(PetscOptionsReal("-tao_cg_delta_min", "minimum delta value", "", cgP->delta_min, &cgP->delta_min, NULL)); 232 PetscCall(PetscOptionsReal("-tao_cg_delta_max", "maximum delta value", "", cgP->delta_max, &cgP->delta_max, NULL)); 233 PetscOptionsHeadEnd(); 234 PetscFunctionReturn(PETSC_SUCCESS); 235 } 236 237 static PetscErrorCode TaoView_CG(Tao tao, PetscViewer viewer) 238 { 239 PetscBool isascii; 240 TAO_CG *cgP = (TAO_CG *)tao->data; 241 242 PetscFunctionBegin; 243 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 244 if (isascii) { 245 PetscCall(PetscViewerASCIIPushTab(viewer)); 246 PetscCall(PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", CG_Table[cgP->cg_type])); 247 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", cgP->ngradsteps)); 248 PetscCall(PetscViewerASCIIPrintf(viewer, "Reset steps: %" PetscInt_FMT "\n", cgP->nresetsteps)); 249 PetscCall(PetscViewerASCIIPopTab(viewer)); 250 } 251 PetscFunctionReturn(PETSC_SUCCESS); 252 } 253 254 /*MC 255 TAOCG - Nonlinear conjugate gradient method is an extension of the 256 nonlinear conjugate gradient solver for nonlinear optimization. 257 258 Options Database Keys: 259 + -tao_cg_eta <r> - restart tolerance 260 . -tao_cg_type <taocg_type> - cg formula 261 . -tao_cg_delta_min <r> - minimum delta value 262 - -tao_cg_delta_max <r> - maximum delta value 263 264 Notes: 265 CG formulas are: 266 "fr" - Fletcher-Reeves 267 "pr" - Polak-Ribiere 268 "prp" - Polak-Ribiere-Plus 269 "hs" - Hestenes-Steifel 270 "dy" - Dai-Yuan 271 Level: beginner 272 M*/ 273 274 PETSC_EXTERN PetscErrorCode TaoCreate_CG(Tao tao) 275 { 276 TAO_CG *cgP; 277 const char *morethuente_type = TAOLINESEARCHMT; 278 279 PetscFunctionBegin; 280 tao->ops->setup = TaoSetUp_CG; 281 tao->ops->solve = TaoSolve_CG; 282 tao->ops->view = TaoView_CG; 283 tao->ops->setfromoptions = TaoSetFromOptions_CG; 284 tao->ops->destroy = TaoDestroy_CG; 285 286 /* Override default settings (unless already changed) */ 287 PetscCall(TaoParametersInitialize(tao)); 288 PetscObjectParameterSetDefault(tao, max_it, 2000); 289 PetscObjectParameterSetDefault(tao, max_funcs, 4000); 290 291 /* Note: nondefault values should be used for nonlinear conjugate gradient */ 292 /* method. In particular, gtol should be less that 0.5; the value used in */ 293 /* Nocedal and Wright is 0.10. We use the default values for the */ 294 /* linesearch because it seems to work better. */ 295 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 296 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 297 PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type)); 298 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 299 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix)); 300 301 PetscCall(PetscNew(&cgP)); 302 tao->data = (void *)cgP; 303 cgP->eta = 0.1; 304 cgP->delta_min = 1e-7; 305 cgP->delta_max = 100; 306 cgP->cg_type = CG_PolakRibierePlus; 307 PetscFunctionReturn(PETSC_SUCCESS); 308 } 309