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