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
TaoSolve_CG(Tao tao)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->gradient, 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
TaoSetUp_CG(Tao tao)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
TaoDestroy_CG(Tao tao)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
TaoSetFromOptions_CG(Tao tao,PetscOptionItems PetscOptionsObject)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
TaoView_CG(Tao tao,PetscViewer viewer)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
TaoCreate_CG(Tao tao)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