Actual source code: bncg.c
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/bound/impls/bncg/bncg.h>
3: #include <petscksp.h>
5: const char *const TaoBNCGTypes[] = {"gd", "pcgd", "hs", "fr", "prp", "prp_plus", "dy", "hz", "dk", "kd", "ssml_bfgs", "ssml_dfp", "ssml_brdn", "TAOBNCGType", "TAO_BNCG_", NULL};
7: #define CG_AS_NONE 0
8: #define CG_AS_BERTSEKAS 1
9: #define CG_AS_SIZE 2
11: static const char *CG_AS_TYPE[64] = {"none", "bertsekas"};
13: PetscErrorCode TaoBNCGEstimateActiveSet(Tao tao, PetscInt asType)
14: {
15: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
17: PetscFunctionBegin;
18: PetscCall(ISDestroy(&cg->inactive_old));
19: if (cg->inactive_idx) {
20: PetscCall(ISDuplicate(cg->inactive_idx, &cg->inactive_old));
21: PetscCall(ISCopy(cg->inactive_idx, cg->inactive_old));
22: }
23: switch (asType) {
24: case CG_AS_NONE:
25: PetscCall(ISDestroy(&cg->inactive_idx));
26: PetscCall(VecWhichInactive(tao->XL, tao->solution, cg->unprojected_gradient, tao->XU, PETSC_TRUE, &cg->inactive_idx));
27: PetscCall(ISDestroy(&cg->active_idx));
28: PetscCall(ISComplementVec(cg->inactive_idx, tao->solution, &cg->active_idx));
29: break;
30: case CG_AS_BERTSEKAS:
31: /* Use gradient descent to estimate the active set */
32: PetscCall(VecCopy(cg->unprojected_gradient, cg->W));
33: PetscCall(VecScale(cg->W, -1.0));
34: PetscCall(TaoEstimateActiveBounds(tao->solution, tao->XL, tao->XU, cg->unprojected_gradient, cg->W, cg->work, cg->as_step, &cg->as_tol, &cg->active_lower, &cg->active_upper, &cg->active_fixed, &cg->active_idx, &cg->inactive_idx));
35: break;
36: default:
37: break;
38: }
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: PetscErrorCode TaoBNCGBoundStep(Tao tao, PetscInt asType, Vec step)
43: {
44: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
46: PetscFunctionBegin;
47: switch (asType) {
48: case CG_AS_NONE:
49: if (cg->active_idx) PetscCall(VecISSet(step, cg->active_idx, 0.0));
50: break;
51: case CG_AS_BERTSEKAS:
52: PetscCall(TaoBoundStep(tao->solution, tao->XL, tao->XU, cg->active_lower, cg->active_upper, cg->active_fixed, 1.0, step));
53: break;
54: default:
55: break;
56: }
57: PetscFunctionReturn(PETSC_SUCCESS);
58: }
60: static PetscErrorCode TaoSolve_BNCG(Tao tao)
61: {
62: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
63: PetscReal step = 1.0, gnorm, gnorm2, resnorm;
64: PetscInt nDiff;
66: PetscFunctionBegin;
67: /* Project the current point onto the feasible set */
68: PetscCall(TaoComputeVariableBounds(tao));
69: PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU));
71: /* Project the initial point onto the feasible region */
72: PetscCall(TaoBoundSolution(tao->solution, tao->XL, tao->XU, 0.0, &nDiff, tao->solution));
74: if (nDiff > 0 || !tao->recycle) PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &cg->f, cg->unprojected_gradient));
75: PetscCall(VecNorm(cg->unprojected_gradient, NORM_2, &gnorm));
76: PetscCheck(!PetscIsInfOrNanReal(cg->f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
78: /* Estimate the active set and compute the projected gradient */
79: PetscCall(TaoBNCGEstimateActiveSet(tao, cg->as_type));
81: /* Project the gradient and calculate the norm */
82: PetscCall(VecCopy(cg->unprojected_gradient, tao->gradient));
83: if (cg->active_idx) PetscCall(VecISSet(tao->gradient, cg->active_idx, 0.0));
84: PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
85: gnorm2 = gnorm * gnorm;
87: /* Initialize counters */
88: tao->niter = 0;
89: cg->ls_fails = cg->descent_error = 0;
90: cg->resets = -1;
91: cg->skipped_updates = cg->pure_gd_steps = 0;
92: cg->iter_quad = 0;
94: /* Convergence test at the starting point. */
95: tao->reason = TAO_CONTINUE_ITERATING;
97: PetscCall(VecFischer(tao->solution, cg->unprojected_gradient, tao->XL, tao->XU, cg->W));
98: PetscCall(VecNorm(cg->W, NORM_2, &resnorm));
99: PetscCheck(!PetscIsInfOrNanReal(resnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
100: PetscCall(TaoLogConvergenceHistory(tao, cg->f, resnorm, 0.0, tao->ksp_its));
101: PetscCall(TaoMonitor(tao, tao->niter, cg->f, resnorm, 0.0, step));
102: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
103: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
104: /* Calculate initial direction. */
105: if (!tao->recycle) {
106: /* We are not recycling a solution/history from a past TaoSolve */
107: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
108: }
109: /* Initial gradient descent step. Scaling by 1.0 also does a decent job for some problems. */
110: while (1) {
111: /* Call general purpose update function */
112: if (tao->ops->update) {
113: PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
114: PetscCall(TaoComputeObjective(tao, tao->solution, &cg->f));
115: }
116: PetscCall(TaoBNCGConductIteration(tao, gnorm));
117: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
118: }
119: }
121: static PetscErrorCode TaoSetUp_BNCG(Tao tao)
122: {
123: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
125: PetscFunctionBegin;
126: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
127: if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
128: if (!cg->W) PetscCall(VecDuplicate(tao->solution, &cg->W));
129: if (!cg->work) PetscCall(VecDuplicate(tao->solution, &cg->work));
130: if (!cg->sk) PetscCall(VecDuplicate(tao->solution, &cg->sk));
131: if (!cg->yk) PetscCall(VecDuplicate(tao->gradient, &cg->yk));
132: if (!cg->X_old) PetscCall(VecDuplicate(tao->solution, &cg->X_old));
133: if (!cg->G_old) PetscCall(VecDuplicate(tao->gradient, &cg->G_old));
134: if (cg->diag_scaling) {
135: PetscCall(VecDuplicate(tao->solution, &cg->d_work));
136: PetscCall(VecDuplicate(tao->solution, &cg->y_work));
137: PetscCall(VecDuplicate(tao->solution, &cg->g_work));
138: }
139: if (!cg->unprojected_gradient) PetscCall(VecDuplicate(tao->gradient, &cg->unprojected_gradient));
140: if (!cg->unprojected_gradient_old) PetscCall(VecDuplicate(tao->gradient, &cg->unprojected_gradient_old));
141: PetscCall(MatLMVMAllocate(cg->B, cg->sk, cg->yk));
142: if (cg->pc) PetscCall(MatLMVMSetJ0(cg->B, cg->pc));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: static PetscErrorCode TaoDestroy_BNCG(Tao tao)
147: {
148: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
150: PetscFunctionBegin;
151: if (tao->setupcalled) {
152: PetscCall(VecDestroy(&cg->W));
153: PetscCall(VecDestroy(&cg->work));
154: PetscCall(VecDestroy(&cg->X_old));
155: PetscCall(VecDestroy(&cg->G_old));
156: PetscCall(VecDestroy(&cg->unprojected_gradient));
157: PetscCall(VecDestroy(&cg->unprojected_gradient_old));
158: PetscCall(VecDestroy(&cg->g_work));
159: PetscCall(VecDestroy(&cg->d_work));
160: PetscCall(VecDestroy(&cg->y_work));
161: PetscCall(VecDestroy(&cg->sk));
162: PetscCall(VecDestroy(&cg->yk));
163: }
164: PetscCall(ISDestroy(&cg->active_lower));
165: PetscCall(ISDestroy(&cg->active_upper));
166: PetscCall(ISDestroy(&cg->active_fixed));
167: PetscCall(ISDestroy(&cg->active_idx));
168: PetscCall(ISDestroy(&cg->inactive_idx));
169: PetscCall(ISDestroy(&cg->inactive_old));
170: PetscCall(ISDestroy(&cg->new_inactives));
171: PetscCall(MatDestroy(&cg->B));
172: PetscCall(MatDestroy(&cg->pc));
173: PetscCall(PetscFree(tao->data));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode TaoSetFromOptions_BNCG(Tao tao, PetscOptionItems PetscOptionsObject)
178: {
179: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
181: PetscFunctionBegin;
182: PetscOptionsHeadBegin(PetscOptionsObject, "Nonlinear Conjugate Gradient method for unconstrained optimization");
183: PetscCall(PetscOptionsEnum("-tao_bncg_type", "CG update formula", "TaoBNCGTypes", TaoBNCGTypes, (PetscEnum)cg->cg_type, (PetscEnum *)&cg->cg_type, NULL));
184: if (cg->cg_type != TAO_BNCG_SSML_BFGS) cg->alpha = -1.0; /* Setting defaults for non-BFGS methods. User can change it below. */
185: if (TAO_BNCG_GD == cg->cg_type) {
186: cg->cg_type = TAO_BNCG_PCGD;
187: /* Set scaling equal to none or, at best, scalar scaling. */
188: cg->unscaled_restart = PETSC_TRUE;
189: cg->diag_scaling = PETSC_FALSE;
190: }
191: PetscCall(PetscOptionsReal("-tao_bncg_hz_eta", "(developer) cutoff tolerance for HZ", "", cg->hz_eta, &cg->hz_eta, NULL));
192: PetscCall(PetscOptionsReal("-tao_bncg_eps", "(developer) cutoff value for restarts", "", cg->epsilon, &cg->epsilon, NULL));
193: PetscCall(PetscOptionsReal("-tao_bncg_dk_eta", "(developer) cutoff tolerance for DK", "", cg->dk_eta, &cg->dk_eta, NULL));
194: PetscCall(PetscOptionsReal("-tao_bncg_xi", "(developer) Parameter in the KD method", "", cg->xi, &cg->xi, NULL));
195: PetscCall(PetscOptionsReal("-tao_bncg_theta", "(developer) update parameter for the Broyden method", "", cg->theta, &cg->theta, NULL));
196: PetscCall(PetscOptionsReal("-tao_bncg_hz_theta", "(developer) parameter for the HZ (2006) method", "", cg->hz_theta, &cg->hz_theta, NULL));
197: PetscCall(PetscOptionsReal("-tao_bncg_alpha", "(developer) parameter for the scalar scaling", "", cg->alpha, &cg->alpha, NULL));
198: PetscCall(PetscOptionsReal("-tao_bncg_bfgs_scale", "(developer) update parameter for bfgs/brdn CG methods", "", cg->bfgs_scale, &cg->bfgs_scale, NULL));
199: PetscCall(PetscOptionsReal("-tao_bncg_dfp_scale", "(developer) update parameter for bfgs/brdn CG methods", "", cg->dfp_scale, &cg->dfp_scale, NULL));
200: PetscCall(PetscOptionsBool("-tao_bncg_diag_scaling", "Enable diagonal Broyden-like preconditioning", "", cg->diag_scaling, &cg->diag_scaling, NULL));
201: PetscCall(PetscOptionsBool("-tao_bncg_dynamic_restart", "(developer) use dynamic restarts as in HZ, DK, KD", "", cg->use_dynamic_restart, &cg->use_dynamic_restart, NULL));
202: PetscCall(PetscOptionsBool("-tao_bncg_unscaled_restart", "(developer) use unscaled gradient restarts", "", cg->unscaled_restart, &cg->unscaled_restart, NULL));
203: PetscCall(PetscOptionsReal("-tao_bncg_zeta", "(developer) Free parameter for the Kou-Dai method", "", cg->zeta, &cg->zeta, NULL));
204: PetscCall(PetscOptionsInt("-tao_bncg_min_quad", "(developer) Number of iterations with approximate quadratic behavior needed for restart", "", cg->min_quad, &cg->min_quad, NULL));
205: PetscCall(PetscOptionsInt("-tao_bncg_min_restart_num", "(developer) Number of iterations between restarts (times dimension)", "", cg->min_restart_num, &cg->min_restart_num, NULL));
206: PetscCall(PetscOptionsBool("-tao_bncg_spaced_restart", "(developer) Enable regular steepest descent restarting every fixed number of iterations", "", cg->spaced_restart, &cg->spaced_restart, NULL));
207: PetscCall(PetscOptionsBool("-tao_bncg_no_scaling", "Disable all scaling except in restarts", "", cg->no_scaling, &cg->no_scaling, NULL));
208: if (cg->no_scaling) {
209: cg->diag_scaling = PETSC_FALSE;
210: cg->alpha = -1.0;
211: }
212: if (cg->alpha == -1.0 && cg->cg_type == TAO_BNCG_KD && !cg->diag_scaling) { /* Some more default options that appear to be good. */
213: cg->neg_xi = PETSC_TRUE;
214: }
215: PetscCall(PetscOptionsBool("-tao_bncg_neg_xi", "(developer) Use negative xi when it might be a smaller descent direction than necessary", "", cg->neg_xi, &cg->neg_xi, NULL));
216: PetscCall(PetscOptionsEList("-tao_bncg_as_type", "active set estimation method", "", CG_AS_TYPE, CG_AS_SIZE, CG_AS_TYPE[cg->as_type], &cg->as_type, NULL));
217: PetscCall(PetscOptionsReal("-tao_bncg_as_tol", "(developer) initial tolerance used when estimating actively bounded variables", "", cg->as_tol, &cg->as_tol, NULL));
218: PetscCall(PetscOptionsReal("-tao_bncg_as_step", "(developer) step length used when estimating actively bounded variables", "", cg->as_step, &cg->as_step, NULL));
219: PetscCall(PetscOptionsReal("-tao_bncg_delta_min", "(developer) minimum scaling factor used for scaled gradient restarts", "", cg->delta_min, &cg->delta_min, NULL));
220: PetscCall(PetscOptionsReal("-tao_bncg_delta_max", "(developer) maximum scaling factor used for scaled gradient restarts", "", cg->delta_max, &cg->delta_max, NULL));
222: PetscOptionsHeadEnd();
223: PetscCall(MatSetOptionsPrefix(cg->B, ((PetscObject)tao)->prefix));
224: PetscCall(MatAppendOptionsPrefix(cg->B, "tao_bncg_"));
225: PetscCall(MatSetFromOptions(cg->B));
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: static PetscErrorCode TaoView_BNCG(Tao tao, PetscViewer viewer)
230: {
231: PetscBool isascii;
232: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
234: PetscFunctionBegin;
235: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
236: if (isascii) {
237: PetscCall(PetscViewerASCIIPushTab(viewer));
238: PetscCall(PetscViewerASCIIPrintf(viewer, "CG Type: %s\n", TaoBNCGTypes[cg->cg_type]));
239: PetscCall(PetscViewerASCIIPrintf(viewer, "Skipped Stepdirection Updates: %" PetscInt_FMT "\n", cg->skipped_updates));
240: PetscCall(PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %" PetscInt_FMT "\n", cg->resets));
241: PetscCall(PetscViewerASCIIPrintf(viewer, "Pure gradient steps: %" PetscInt_FMT "\n", cg->pure_gd_steps));
242: PetscCall(PetscViewerASCIIPrintf(viewer, "Not a descent direction: %" PetscInt_FMT "\n", cg->descent_error));
243: PetscCall(PetscViewerASCIIPrintf(viewer, "Line search fails: %" PetscInt_FMT "\n", cg->ls_fails));
244: if (cg->diag_scaling) {
245: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
246: if (isascii) {
247: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO));
248: PetscCall(MatView(cg->B, viewer));
249: PetscCall(PetscViewerPopFormat(viewer));
250: }
251: }
252: PetscCall(PetscViewerASCIIPopTab(viewer));
253: }
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: PetscErrorCode TaoBNCGComputeScalarScaling(PetscReal yty, PetscReal yts, PetscReal sts, PetscReal *scale, PetscReal alpha)
258: {
259: PetscReal a, b, c, sig1, sig2;
261: PetscFunctionBegin;
262: *scale = 0.0;
263: if (1.0 == alpha) *scale = yts / yty;
264: else if (0.0 == alpha) *scale = sts / yts;
265: else if (-1.0 == alpha) *scale = 1.0;
266: else {
267: a = yty;
268: b = yts;
269: c = sts;
270: a *= alpha;
271: b *= -(2.0 * alpha - 1.0);
272: c *= alpha - 1.0;
273: sig1 = (-b + PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
274: sig2 = (-b - PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
275: /* accept the positive root as the scalar */
276: if (sig1 > 0.0) *scale = sig1;
277: else if (sig2 > 0.0) *scale = sig2;
278: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
279: }
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: /*MC
284: TAOBNCG - Bound-constrained Nonlinear Conjugate Gradient method.
286: Options Database Keys:
287: + -tao_bncg_recycle - enable recycling the latest calculated gradient vector in subsequent `TaoSolve()` calls (currently disabled)
288: . -tao_bncg_eta r - restart tolerance
289: . -tao_bncg_type taocg_type - cg formula
290: . -tao_bncg_as_type (none|bertsekas) - active set estimation method
291: . -tao_bncg_as_tol r - tolerance used in Bertsekas active-set estimation
292: . -tao_bncg_as_step r - trial step length used in Bertsekas active-set estimation
293: . -tao_bncg_eps r - cutoff used for determining whether or not we restart based on steplength each iteration,
294: as well as determining whether or not we continue using the last stepdirection. Defaults to machine precision.
295: . -tao_bncg_theta r - convex combination parameter for the Broyden method
296: . -tao_bncg_hz_eta r - cutoff tolerance for the beta term in the `hz`, `dk` methods
297: . -tao_bncg_dk_eta r - cutoff tolerance for the beta term in the `hz`, `dk` methods
298: . -tao_bncg_xi r - Multiplicative constant of the gamma term in the `kd` method
299: . -tao_bncg_hz_theta r - Multiplicative constant of the theta term for the `hz` method
300: . -tao_bncg_bfgs_scale r - Scaling parameter of the BFGS contribution to the scalar Broyden method
301: . -tao_bncg_dfp_scale r - Scaling parameter of the dfp contribution to the scalar Broyden method
302: . -tao_bncg_diag_scaling b - Whether or not to use diagonal initialization/preconditioning for the CG methods. Default True.
303: . -tao_bncg_dynamic_restart b - use dynamic restart strategy in the `hz`, `dk`, `kd` methods
304: . -tao_bncg_unscaled_restart b - whether or not to scale the gradient when doing gradient descent restarts
305: . -tao_bncg_zeta r - Scaling parameter in the `kd` method
306: . -tao_bncg_delta_min r - Minimum bound for rescaling during restarted gradient descent steps
307: . -tao_bncg_delta_max r - Maximum bound for rescaling during restarted gradient descent steps
308: . -tao_bncg_min_quad i - Number of quadratic-like steps in a row necessary to do a dynamic restart
309: . -tao_bncg_min_restart_num i - This number, x, makes sure there is a gradient descent step every $x*n$ iterations, where `n` is the dimension of the problem
310: . -tao_bncg_spaced_restart (true|false) - whether or not to do gradient descent steps every x*n iterations
311: . -tao_bncg_no_scaling b - If true, eliminates all scaling, including defaults.
312: - -tao_bncg_neg_xi b - Whether or not to use negative xi in the `kd` method under certain conditions
314: Note:
315: CG formulas are:
316: + `gd` - Gradient Descent
317: . `fr` - Fletcher-Reeves
318: . `pr` - Polak-Ribiere-Polyak
319: . `prp` - Polak-Ribiere-Plus
320: . `hs` - Hestenes-Steifel
321: . `dy` - Dai-Yuan
322: . `ssml_bfgs` - Self-Scaling Memoryless BFGS
323: . `ssml_dfp` - Self-Scaling Memoryless DFP
324: . `ssml_brdn` - Self-Scaling Memoryless Broyden
325: . `hz` - Hager-Zhang (CG_DESCENT 5.3)
326: . `dk` - Dai-Kou (2013)
327: - `kd` - Kou-Dai (2015)
329: Level: beginner
331: The various algorithmic factors can only be supplied via the options database
333: .seealso: `Tao`, `TAONTR`, `TAONTL`, `TAONM`, `TAOCG`, `TaoType`, `TaoCreate()`
334: M*/
336: PETSC_EXTERN PetscErrorCode TaoCreate_BNCG(Tao tao)
337: {
338: TAO_BNCG *cg;
339: const char *morethuente_type = TAOLINESEARCHMT;
341: PetscFunctionBegin;
342: tao->ops->setup = TaoSetUp_BNCG;
343: tao->ops->solve = TaoSolve_BNCG;
344: tao->ops->view = TaoView_BNCG;
345: tao->ops->setfromoptions = TaoSetFromOptions_BNCG;
346: tao->ops->destroy = TaoDestroy_BNCG;
347: tao->uses_gradient = PETSC_TRUE;
349: /* Override default settings (unless already changed) */
350: PetscObjectParameterSetDefault(tao, max_it, 2000);
351: PetscObjectParameterSetDefault(tao, max_funcs, 4000);
353: /* Note: nondefault values should be used for nonlinear conjugate gradient */
354: /* method. In particular, gtol should be less than 0.5; the value used in */
355: /* Nocedal and Wright is 0.10. We use the default values for the */
356: /* linesearch because it seems to work better. */
357: PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
358: PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
359: PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
360: PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
362: PetscCall(PetscNew(&cg));
363: tao->data = (void *)cg;
364: PetscCall(KSPInitializePackage());
365: PetscCall(MatCreate(PetscObjectComm((PetscObject)tao), &cg->B));
366: PetscCall(PetscObjectIncrementTabLevel((PetscObject)cg->B, (PetscObject)tao, 1));
367: PetscCall(MatSetType(cg->B, MATLMVMDIAGBROYDEN));
369: cg->dk_eta = 0.5;
370: cg->hz_eta = 0.4;
371: cg->delta_min = 1e-7;
372: cg->delta_max = 100;
373: cg->theta = 1.0;
374: cg->hz_theta = 1.0;
375: cg->dfp_scale = 1.0;
376: cg->bfgs_scale = 1.0;
377: cg->zeta = 0.1;
378: cg->min_quad = 6;
379: cg->min_restart_num = 6; /* As in CG_DESCENT and KD2015*/
380: cg->xi = 1.0;
381: cg->neg_xi = PETSC_TRUE;
382: cg->tol_quad = 1e-8;
383: cg->as_step = 0.001;
384: cg->as_tol = 0.001;
385: cg->eps_23 = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0); /* Just a little tighter*/
386: cg->as_type = CG_AS_BERTSEKAS;
387: cg->cg_type = TAO_BNCG_SSML_BFGS;
388: cg->alpha = 1.0;
389: cg->diag_scaling = PETSC_TRUE;
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: PetscErrorCode TaoBNCGResetUpdate(Tao tao, PetscReal gnormsq)
394: {
395: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
396: PetscReal scaling;
398: PetscFunctionBegin;
399: ++cg->resets;
400: scaling = 2.0 * PetscMax(1.0, PetscAbsScalar(cg->f)) / PetscMax(gnormsq, cg->eps_23);
401: scaling = PetscMin(cg->delta_max, PetscMax(cg->delta_min, scaling));
402: if (cg->unscaled_restart) {
403: scaling = 1.0;
404: ++cg->pure_gd_steps;
405: }
406: PetscCall(VecAXPBY(tao->stepdirection, -scaling, 0.0, tao->gradient));
407: /* Also want to reset our diagonal scaling with each restart */
408: if (cg->diag_scaling) PetscCall(MatLMVMReset(cg->B, PETSC_FALSE));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: PetscErrorCode TaoBNCGCheckDynamicRestart(Tao tao, PetscReal stepsize, PetscReal gd, PetscReal gd_old, PetscBool *dynrestart, PetscReal fold)
413: {
414: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
415: PetscReal quadinterp;
417: PetscFunctionBegin;
418: if (cg->f < cg->min_quad / 10) {
419: *dynrestart = PETSC_FALSE;
420: PetscFunctionReturn(PETSC_SUCCESS);
421: } /* just skip this since this strategy doesn't work well for functions near zero */
422: quadinterp = 2.0 * (cg->f - fold) / (stepsize * (gd + gd_old));
423: if (PetscAbs(quadinterp - 1.0) < cg->tol_quad) ++cg->iter_quad;
424: else {
425: cg->iter_quad = 0;
426: *dynrestart = PETSC_FALSE;
427: }
428: if (cg->iter_quad >= cg->min_quad) {
429: cg->iter_quad = 0;
430: *dynrestart = PETSC_TRUE;
431: }
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: PETSC_INTERN PetscErrorCode TaoBNCGStepDirectionUpdate(Tao tao, PetscReal gnorm2, PetscReal step, PetscReal fold, PetscReal gnorm2_old, PetscReal dnorm, PetscBool pcgd_fallback)
436: {
437: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
438: PetscReal gamma = 1.0, tau_k, beta;
439: PetscReal tmp = 1.0, ynorm, ynorm2 = 1.0, snorm = 1.0, dk_yk = 1.0, gd;
440: PetscReal gkp1_yk, gd_old, tau_bfgs, tau_dfp, gkp1D_yk, gtDg;
441: PetscInt dim;
442: PetscBool cg_restart = PETSC_FALSE;
444: PetscFunctionBegin;
445: /* Local curvature check to see if we need to restart */
446: if (tao->niter >= 1 || tao->recycle) {
447: PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
448: PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
449: ynorm2 = ynorm * ynorm;
450: PetscCall(VecDot(cg->yk, tao->stepdirection, &dk_yk));
451: if (step * dnorm < PETSC_MACHINE_EPSILON || step * dk_yk < PETSC_MACHINE_EPSILON) {
452: cg_restart = PETSC_TRUE;
453: ++cg->skipped_updates;
454: }
455: }
456: /* If the user wants regular restarts, do it every min_restart_num * n iterations, where n = dimension */
457: if (cg->spaced_restart) {
458: PetscCall(VecGetSize(tao->gradient, &dim));
459: if (0 == tao->niter % (cg->min_restart_num * dim)) cg_restart = PETSC_TRUE;
460: }
461: /* Compute the diagonal scaling vector if applicable */
462: if (cg->diag_scaling) PetscCall(MatLMVMUpdate(cg->B, tao->solution, tao->gradient));
464: /* A note on diagonal scaling (to be added to paper):
465: For the FR, PR, PRP, and DY methods, the diagonally scaled versions
466: must be derived as a preconditioned CG method rather than as
467: a Hessian initialization like in the Broyden methods. */
469: /* In that case, one writes the objective function as
470: f(x) \equiv f(Ay). Gradient evaluations yield g(x_k) = A g(Ay_k) = A g(x_k).
471: Furthermore, the direction d_k \equiv (x_k - x_{k-1})/step according to
472: HZ (2006) becomes A^{-1} d_k, such that d_k^T g_k remains the
473: same under preconditioning. Note that A is diagonal, such that A^T = A. */
475: /* This yields questions like what the dot product d_k^T y_k
476: should look like. HZ mistakenly treats that as the same under
477: preconditioning, but that is not necessarily true. */
479: /* Observe y_k \equiv g_k - g_{k-1}, and under the P.C. transformation,
480: we get d_k^T y_k = (d_k^T A_k^{-T} A_k g_k - d_k^T A_k^{-T} A_{k-1} g_{k-1}),
481: yielding d_k^T y_k = d_k^T g_k - d_k^T (A_k^{-T} A_{k-1} g_{k-1}), which is
482: NOT the same if our matrix used to construct the preconditioner is updated between iterations.
483: This same issue is found when considering dot products of the form g_{k+1}^T y_k. */
485: /* Compute CG step direction */
486: if (cg_restart) {
487: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
488: } else if (pcgd_fallback) {
489: /* Just like preconditioned CG */
490: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
491: PetscCall(VecAXPBY(tao->stepdirection, -1.0, 0.0, cg->g_work));
492: } else if (ynorm2 > PETSC_MACHINE_EPSILON) {
493: switch (cg->cg_type) {
494: case TAO_BNCG_PCGD:
495: if (!cg->diag_scaling) {
496: if (!cg->no_scaling) {
497: cg->sts = step * step * dnorm * dnorm;
498: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->sts, &tau_k, cg->alpha));
499: } else {
500: tau_k = 1.0;
501: ++cg->pure_gd_steps;
502: }
503: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, 0.0, tao->gradient));
504: } else {
505: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
506: PetscCall(VecAXPBY(tao->stepdirection, -1.0, 0.0, cg->g_work));
507: }
508: break;
510: case TAO_BNCG_HS:
511: /* Classic Hestenes-Stiefel method, modified with scalar and diagonal preconditioning. */
512: if (!cg->diag_scaling) {
513: cg->sts = step * step * dnorm * dnorm;
514: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
515: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->sts, &tau_k, cg->alpha));
516: beta = tau_k * gkp1_yk / dk_yk;
517: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
518: } else {
519: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
520: PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
521: beta = gkp1_yk / dk_yk;
522: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
523: }
524: break;
526: case TAO_BNCG_FR:
527: PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
528: PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
529: PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
530: ynorm2 = ynorm * ynorm;
531: PetscCall(VecDot(cg->yk, tao->stepdirection, &dk_yk));
532: if (!cg->diag_scaling) {
533: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, step * step * dnorm * dnorm, &tau_k, cg->alpha));
534: beta = tau_k * gnorm2 / gnorm2_old;
535: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
536: } else {
537: PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old)); /* Before it's updated */
538: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
539: PetscCall(VecDot(tao->gradient, cg->g_work, &tmp));
540: beta = tmp / gnorm2_old;
541: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
542: }
543: break;
545: case TAO_BNCG_PRP:
546: snorm = step * dnorm;
547: if (!cg->diag_scaling) {
548: PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
549: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
550: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, snorm * snorm, &tau_k, cg->alpha));
551: beta = tau_k * gkp1_yk / gnorm2_old;
552: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
553: } else {
554: PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old));
555: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
556: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
557: beta = gkp1_yk / gnorm2_old;
558: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
559: }
560: break;
562: case TAO_BNCG_PRP_PLUS:
563: PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
564: PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
565: ynorm2 = ynorm * ynorm;
566: if (!cg->diag_scaling) {
567: PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
568: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
569: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, snorm * snorm, &tau_k, cg->alpha));
570: beta = tau_k * gkp1_yk / gnorm2_old;
571: beta = PetscMax(beta, 0.0);
572: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
573: } else {
574: PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old)); /* Old gtDg */
575: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
576: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
577: beta = gkp1_yk / gnorm2_old;
578: beta = PetscMax(beta, 0.0);
579: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
580: }
581: break;
583: case TAO_BNCG_DY:
584: /* Dai, Yu-Hong, and Yaxiang Yuan. "A nonlinear conjugate gradient method with a strong global convergence property."
585: SIAM Journal on optimization 10, no. 1 (1999): 177-182. */
586: if (!cg->diag_scaling) {
587: PetscCall(VecDot(tao->stepdirection, tao->gradient, &gd));
588: PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
589: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->yts, &tau_k, cg->alpha));
590: beta = tau_k * gnorm2 / (gd - gd_old);
591: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
592: } else {
593: PetscCall(MatMult(cg->B, tao->stepdirection, cg->d_work));
594: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
595: PetscCall(VecDot(cg->g_work, tao->gradient, >Dg));
596: PetscCall(VecDot(tao->stepdirection, cg->G_old, &gd_old));
597: PetscCall(VecDot(cg->d_work, cg->g_work, &dk_yk));
598: dk_yk = dk_yk - gd_old;
599: beta = gtDg / dk_yk;
600: PetscCall(VecScale(cg->d_work, beta));
601: PetscCall(VecWAXPY(tao->stepdirection, -1.0, cg->g_work, cg->d_work));
602: }
603: break;
605: case TAO_BNCG_HZ:
606: /* Hager, William W., and Hongchao Zhang. "Algorithm 851: CG_DESCENT, a conjugate gradient method with guaranteed descent."
607: ACM Transactions on Mathematical Software (TOMS) 32, no. 1 (2006): 113-137. */
608: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
609: PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
610: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
611: snorm = dnorm * step;
612: cg->yts = step * dk_yk;
613: if (cg->use_dynamic_restart) PetscCall(TaoBNCGCheckDynamicRestart(tao, step, gd, gd_old, &cg->dynamic_restart, fold));
614: if (cg->dynamic_restart) {
615: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
616: } else {
617: if (!cg->diag_scaling) {
618: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
619: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
620: /* Supplying cg->alpha = -1.0 will give the CG_DESCENT 5.3 special case of tau_k = 1.0 */
621: tmp = gd / dk_yk;
622: beta = tau_k * (gkp1_yk / dk_yk - ynorm2 * gd / (dk_yk * dk_yk));
623: /* Bound beta as in CG_DESCENT 5.3, as implemented, with the third comparison from DK 2013 */
624: beta = PetscMax(PetscMax(beta, cg->hz_eta * tau_k * gd_old / (dnorm * dnorm)), cg->dk_eta * tau_k * gd / (dnorm * dnorm));
625: /* d <- -t*g + beta*t*d */
626: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
627: } else {
628: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless BFGS step */
629: cg->yty = ynorm2;
630: cg->sts = snorm * snorm;
631: /* Apply the diagonal scaling to all my vectors */
632: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
633: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
634: PetscCall(MatSolve(cg->B, tao->stepdirection, cg->d_work));
635: /* Construct the constant ytDgkp1 */
636: PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
637: /* Construct the constant for scaling Dkyk in the update */
638: tmp = gd / dk_yk;
639: PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
640: tau_k = -tau_k * gd / (dk_yk * dk_yk);
641: /* beta is the constant which adds the dk contribution */
642: beta = gkp1_yk / dk_yk + cg->hz_theta * tau_k; /* HZ; (1.15) from DK 2013 */
643: /* From HZ2013, modified to account for diagonal scaling*/
644: PetscCall(VecDot(cg->G_old, cg->d_work, &gd_old));
645: PetscCall(VecDot(tao->stepdirection, cg->g_work, &gd));
646: beta = PetscMax(PetscMax(beta, cg->hz_eta * gd_old / (dnorm * dnorm)), cg->dk_eta * gd / (dnorm * dnorm));
647: /* Do the update */
648: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
649: }
650: }
651: break;
653: case TAO_BNCG_DK:
654: /* Dai, Yu-Hong, and Cai-Xia Kou. "A nonlinear conjugate gradient algorithm with an optimal property and an improved Wolfe line search."
655: SIAM Journal on Optimization 23, no. 1 (2013): 296-320. */
656: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
657: PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
658: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
659: snorm = step * dnorm;
660: cg->yts = dk_yk * step;
661: if (!cg->diag_scaling) {
662: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
663: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
664: /* Use cg->alpha = -1.0 to get tau_k = 1.0 as in CG_DESCENT 5.3 */
665: tmp = gd / dk_yk;
666: beta = tau_k * (gkp1_yk / dk_yk - ynorm2 * gd / (dk_yk * dk_yk) + gd / (dnorm * dnorm)) - step * gd / dk_yk;
667: beta = PetscMax(PetscMax(beta, cg->hz_eta * tau_k * gd_old / (dnorm * dnorm)), cg->dk_eta * tau_k * gd / (dnorm * dnorm));
668: /* d <- -t*g + beta*t*d */
669: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, 0.0, beta, tao->gradient, cg->yk));
670: } else {
671: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless BFGS step */
672: cg->yty = ynorm2;
673: cg->sts = snorm * snorm;
674: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
675: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
676: PetscCall(MatSolve(cg->B, tao->stepdirection, cg->d_work));
677: /* Construct the constant ytDgkp1 */
678: PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
679: PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
680: tau_k = tau_k * gd / (dk_yk * dk_yk);
681: tmp = gd / dk_yk;
682: /* beta is the constant which adds the dk contribution */
683: beta = gkp1_yk / dk_yk - step * tmp - tau_k;
684: /* Update this for the last term in beta */
685: PetscCall(VecDot(cg->y_work, tao->stepdirection, &dk_yk));
686: beta += tmp * dk_yk / (dnorm * dnorm); /* projection of y_work onto dk */
687: PetscCall(VecDot(tao->stepdirection, cg->g_work, &gd));
688: PetscCall(VecDot(cg->G_old, cg->d_work, &gd_old));
689: beta = PetscMax(PetscMax(beta, cg->hz_eta * gd_old / (dnorm * dnorm)), cg->dk_eta * gd / (dnorm * dnorm));
690: /* Do the update */
691: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
692: }
693: break;
695: case TAO_BNCG_KD:
696: /* Kou, Cai-Xia, and Yu-Hong Dai. "A modified self-scaling memoryless Broyden-Fletcher-Goldfarb-Shanno method for unconstrained optimization."
697: Journal of Optimization Theory and Applications 165, no. 1 (2015): 209-224. */
698: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
699: PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
700: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
701: snorm = step * dnorm;
702: cg->yts = dk_yk * step;
703: if (cg->use_dynamic_restart) PetscCall(TaoBNCGCheckDynamicRestart(tao, step, gd, gd_old, &cg->dynamic_restart, fold));
704: if (cg->dynamic_restart) {
705: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
706: } else {
707: if (!cg->diag_scaling) {
708: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
709: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
710: beta = tau_k * (gkp1_yk / dk_yk - ynorm2 * gd / (dk_yk * dk_yk)) - step * gd / dk_yk;
711: if (beta < cg->zeta * tau_k * gd / (dnorm * dnorm)) /* 0.1 is KD's zeta parameter */
712: {
713: beta = cg->zeta * tau_k * gd / (dnorm * dnorm);
714: gamma = 0.0;
715: } else {
716: if (gkp1_yk < 0 && cg->neg_xi) gamma = -1.0 * gd / dk_yk;
717: /* This seems to be very effective when there's no tau_k scaling.
718: This guarantees a large descent step every iteration, going through DK 2015 Lemma 3.1's proof but allowing for negative xi */
719: else gamma = cg->xi * gd / dk_yk;
720: }
721: /* d <- -t*g + beta*t*d + t*tmp*yk */
722: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, gamma * tau_k, beta, tao->gradient, cg->yk));
723: } else {
724: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless BFGS step */
725: cg->yty = ynorm2;
726: cg->sts = snorm * snorm;
727: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
728: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
729: /* Construct the constant ytDgkp1 */
730: PetscCall(VecDot(cg->yk, cg->g_work, &gkp1D_yk));
731: /* Construct the constant for scaling Dkyk in the update */
732: gamma = gd / dk_yk;
733: /* tau_k = -ytDy/(ytd)^2 * gd */
734: PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
735: tau_k = tau_k * gd / (dk_yk * dk_yk);
736: /* beta is the constant which adds the d_k contribution */
737: beta = gkp1D_yk / dk_yk - step * gamma - tau_k;
738: /* Here is the requisite check */
739: PetscCall(VecDot(tao->stepdirection, cg->g_work, &tmp));
740: if (cg->neg_xi) {
741: /* modified KD implementation */
742: if (gkp1D_yk / dk_yk < 0) gamma = -1.0 * gd / dk_yk;
743: else gamma = cg->xi * gd / dk_yk;
744: if (beta < cg->zeta * tmp / (dnorm * dnorm)) {
745: beta = cg->zeta * tmp / (dnorm * dnorm);
746: gamma = 0.0;
747: }
748: } else { /* original KD 2015 implementation */
749: if (beta < cg->zeta * tmp / (dnorm * dnorm)) {
750: beta = cg->zeta * tmp / (dnorm * dnorm);
751: gamma = 0.0;
752: } else gamma = cg->xi * gd / dk_yk;
753: }
754: /* Do the update in two steps */
755: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
756: PetscCall(VecAXPY(tao->stepdirection, gamma, cg->y_work));
757: }
758: }
759: break;
761: case TAO_BNCG_SSML_BFGS:
762: /* Perry, J. M. "A class of conjugate gradient algorithms with a two-step variable-metric memory."
763: Discussion Papers 269 (1977). */
764: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
765: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
766: snorm = step * dnorm;
767: cg->yts = dk_yk * step;
768: cg->yty = ynorm2;
769: cg->sts = snorm * snorm;
770: if (!cg->diag_scaling) {
771: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
772: PetscCall(TaoBNCGComputeScalarScaling(cg->yty, cg->yts, cg->sts, &tau_k, cg->alpha));
773: tmp = gd / dk_yk;
774: beta = tau_k * (gkp1_yk / dk_yk - cg->yty * gd / (dk_yk * dk_yk)) - step * tmp;
775: /* d <- -t*g + beta*t*d + t*tmp*yk */
776: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp * tau_k, beta, tao->gradient, cg->yk));
777: } else {
778: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless BFGS step */
779: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
780: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
781: /* compute scalar gamma */
782: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
783: PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
784: gamma = gd / dk_yk;
785: /* Compute scalar beta */
786: beta = (gkp1_yk / dk_yk - gd * tmp / (dk_yk * dk_yk)) - step * gd / dk_yk;
787: /* Compute stepdirection d_kp1 = gamma*Dkyk + beta*dk - Dkgkp1 */
788: PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
789: }
790: break;
792: case TAO_BNCG_SSML_DFP:
793: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
794: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
795: snorm = step * dnorm;
796: cg->yts = dk_yk * step;
797: cg->yty = ynorm2;
798: cg->sts = snorm * snorm;
799: if (!cg->diag_scaling) {
800: /* Instead of a regular convex combination, we will solve a quadratic formula. */
801: PetscCall(TaoBNCGComputeScalarScaling(cg->yty, cg->yts, cg->sts, &tau_k, cg->alpha));
802: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
803: tau_k = cg->dfp_scale * tau_k;
804: tmp = tau_k * gkp1_yk / cg->yty;
805: beta = -step * gd / dk_yk;
806: /* d <- -t*g + beta*d + tmp*yk */
807: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp, beta, tao->gradient, cg->yk));
808: } else {
809: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless DFP step */
810: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
811: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
812: /* compute scalar gamma */
813: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
814: PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
815: gamma = (gkp1_yk / tmp);
816: /* Compute scalar beta */
817: beta = -step * gd / dk_yk;
818: /* Compute stepdirection d_kp1 = gamma*Dkyk + beta*dk - Dkgkp1 */
819: PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
820: }
821: break;
823: case TAO_BNCG_SSML_BRDN:
824: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
825: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
826: snorm = step * dnorm;
827: cg->yts = step * dk_yk;
828: cg->yty = ynorm2;
829: cg->sts = snorm * snorm;
830: if (!cg->diag_scaling) {
831: /* Instead of a regular convex combination, we will solve a quadratic formula. */
832: PetscCall(TaoBNCGComputeScalarScaling(cg->yty, step * dk_yk, snorm * snorm, &tau_bfgs, cg->bfgs_scale));
833: PetscCall(TaoBNCGComputeScalarScaling(cg->yty, step * dk_yk, snorm * snorm, &tau_dfp, cg->dfp_scale));
834: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
835: tau_k = cg->theta * tau_bfgs + (1.0 - cg->theta) * tau_dfp;
836: /* If bfgs_scale = 1.0, it should reproduce the bfgs tau_bfgs. If bfgs_scale = 0.0,
837: it should reproduce the tau_dfp scaling. Same with dfp_scale. */
838: tmp = cg->theta * tau_bfgs * gd / dk_yk + (1 - cg->theta) * tau_dfp * gkp1_yk / cg->yty;
839: beta = cg->theta * tau_bfgs * (gkp1_yk / dk_yk - cg->yty * gd / (dk_yk * dk_yk)) - step * gd / dk_yk;
840: /* d <- -t*g + beta*d + tmp*yk */
841: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp, beta, tao->gradient, cg->yk));
842: } else {
843: /* We have diagonal scaling enabled */
844: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
845: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
846: /* compute scalar gamma */
847: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
848: PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
849: gamma = cg->theta * gd / dk_yk + (1 - cg->theta) * (gkp1_yk / tmp);
850: /* Compute scalar beta */
851: beta = cg->theta * (gkp1_yk / dk_yk - gd * tmp / (dk_yk * dk_yk)) - step * gd / dk_yk;
852: /* Compute stepdirection dkp1 = gamma*Dkyk + beta*dk - Dkgkp1 */
853: PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
854: }
855: break;
857: default:
858: break;
859: }
860: }
861: PetscFunctionReturn(PETSC_SUCCESS);
862: }
864: PETSC_INTERN PetscErrorCode TaoBNCGConductIteration(Tao tao, PetscReal gnorm)
865: {
866: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
867: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
868: PetscReal step = 1.0, gnorm2, gd, dnorm = 0.0;
869: PetscReal gnorm2_old, f_old, resnorm, gnorm_old;
870: PetscBool pcgd_fallback = PETSC_FALSE;
872: PetscFunctionBegin;
873: /* We are now going to perform a line search along the direction. */
874: /* Store solution and gradient info before it changes */
875: PetscCall(VecCopy(tao->solution, cg->X_old));
876: PetscCall(VecCopy(tao->gradient, cg->G_old));
877: PetscCall(VecCopy(cg->unprojected_gradient, cg->unprojected_gradient_old));
879: gnorm_old = gnorm;
880: gnorm2_old = gnorm_old * gnorm_old;
881: f_old = cg->f;
882: /* Perform bounded line search. If we are recycling a solution from a previous */
883: /* TaoSolve, then we want to immediately skip to calculating a new direction rather than performing a linesearch */
884: if (!(tao->recycle && 0 == tao->niter)) {
885: /* Above logic: the below code happens every iteration, except for the first iteration of a recycled TaoSolve */
886: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
887: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
888: PetscCall(TaoAddLineSearchCounts(tao));
890: /* Check linesearch failure */
891: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
892: ++cg->ls_fails;
893: if (cg->cg_type == TAO_BNCG_GD) {
894: /* Nothing left to do but fail out of the optimization */
895: step = 0.0;
896: tao->reason = TAO_DIVERGED_LS_FAILURE;
897: } else {
898: /* Restore previous point, perform preconditioned GD and regular GD steps at the last good point */
899: PetscCall(VecCopy(cg->X_old, tao->solution));
900: PetscCall(VecCopy(cg->G_old, tao->gradient));
901: PetscCall(VecCopy(cg->unprojected_gradient_old, cg->unprojected_gradient));
902: gnorm = gnorm_old;
903: gnorm2 = gnorm2_old;
904: cg->f = f_old;
906: /* Fall back on preconditioned CG (so long as you're not already using it) */
907: if (cg->cg_type != TAO_BNCG_PCGD && cg->diag_scaling) {
908: pcgd_fallback = PETSC_TRUE;
909: PetscCall(TaoBNCGStepDirectionUpdate(tao, gnorm2, step, f_old, gnorm2_old, dnorm, pcgd_fallback));
911: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
912: PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
914: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
915: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
916: PetscCall(TaoAddLineSearchCounts(tao));
918: pcgd_fallback = PETSC_FALSE;
919: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
920: /* Going to perform a regular gradient descent step. */
921: ++cg->ls_fails;
922: step = 0.0;
923: }
924: }
925: /* Fall back on the scaled gradient step */
926: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
927: ++cg->ls_fails;
928: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
929: PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
930: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
931: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
932: PetscCall(TaoAddLineSearchCounts(tao));
933: }
935: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
936: /* Nothing left to do but fail out of the optimization */
937: ++cg->ls_fails;
938: step = 0.0;
939: tao->reason = TAO_DIVERGED_LS_FAILURE;
940: } else {
941: /* One of the fallbacks worked. Set them both back equal to false. */
942: pcgd_fallback = PETSC_FALSE;
943: }
944: }
945: }
946: /* Convergence test for line search failure */
947: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
949: /* Standard convergence test */
950: PetscCall(VecFischer(tao->solution, cg->unprojected_gradient, tao->XL, tao->XU, cg->W));
951: PetscCall(VecNorm(cg->W, NORM_2, &resnorm));
952: PetscCheck(!PetscIsInfOrNanReal(resnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
953: PetscCall(TaoLogConvergenceHistory(tao, cg->f, resnorm, 0.0, tao->ksp_its));
954: PetscCall(TaoMonitor(tao, tao->niter, cg->f, resnorm, 0.0, step));
955: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
956: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
957: }
958: /* Assert we have an updated step and we need at least one more iteration. */
959: /* Calculate the next direction */
960: /* Estimate the active set at the new solution */
961: PetscCall(TaoBNCGEstimateActiveSet(tao, cg->as_type));
962: /* Compute the projected gradient and its norm */
963: PetscCall(VecCopy(cg->unprojected_gradient, tao->gradient));
964: if (cg->active_idx) PetscCall(VecISSet(tao->gradient, cg->active_idx, 0.0));
965: PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
966: gnorm2 = gnorm * gnorm;
968: /* Calculate some quantities used in the StepDirectionUpdate. */
969: PetscCall(VecNorm(tao->stepdirection, NORM_2, &dnorm));
970: /* Update the step direction. */
971: PetscCall(TaoBNCGStepDirectionUpdate(tao, gnorm2, step, f_old, gnorm2_old, dnorm, pcgd_fallback));
972: ++tao->niter;
973: PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
975: if (cg->cg_type != TAO_BNCG_GD) {
976: /* Figure out which previously active variables became inactive this iteration */
977: PetscCall(ISDestroy(&cg->new_inactives));
978: if (cg->inactive_idx && cg->inactive_old) PetscCall(ISDifference(cg->inactive_idx, cg->inactive_old, &cg->new_inactives));
979: /* Selectively reset the CG step those freshly inactive variables */
980: if (cg->new_inactives) {
981: PetscCall(VecGetSubVector(tao->stepdirection, cg->new_inactives, &cg->inactive_step));
982: PetscCall(VecGetSubVector(cg->unprojected_gradient, cg->new_inactives, &cg->inactive_grad));
983: PetscCall(VecCopy(cg->inactive_grad, cg->inactive_step));
984: PetscCall(VecScale(cg->inactive_step, -1.0));
985: PetscCall(VecRestoreSubVector(tao->stepdirection, cg->new_inactives, &cg->inactive_step));
986: PetscCall(VecRestoreSubVector(cg->unprojected_gradient, cg->new_inactives, &cg->inactive_grad));
987: }
988: /* Verify that this is a descent direction */
989: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
990: PetscCall(VecNorm(tao->stepdirection, NORM_2, &dnorm));
991: if (PetscIsInfOrNanReal(gd) || (gd / (dnorm * dnorm) <= -1e10 || gd / (dnorm * dnorm) >= -1e-10)) {
992: /* Not a descent direction, so we reset back to projected gradient descent */
993: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
994: PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
995: ++cg->descent_error;
996: } else {
997: }
998: }
999: PetscFunctionReturn(PETSC_SUCCESS);
1000: }
1002: PETSC_INTERN PetscErrorCode TaoBNCGSetH0(Tao tao, Mat H0)
1003: {
1004: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1005: PetscBool same;
1007: PetscFunctionBegin;
1008: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBNCG, &same));
1009: if (same) {
1010: PetscCall(PetscObjectReference((PetscObject)H0));
1011: cg->pc = H0;
1012: }
1013: PetscFunctionReturn(PETSC_SUCCESS);
1014: }
1016: /*@
1017: TaoBNCGGetType - Return the type for the `TAOBNCG` solver
1019: Input Parameter:
1020: . tao - the `Tao` solver context
1022: Output Parameter:
1023: . type - `TAOBNCG` type
1025: Level: advanced
1027: .seealso: `Tao`, `TAOBNCG`, `TaoBNCGSetType()`, `TaoBNCGType`
1028: @*/
1029: PetscErrorCode TaoBNCGGetType(Tao tao, TaoBNCGType *type)
1030: {
1031: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1032: PetscBool same;
1034: PetscFunctionBegin;
1035: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBNCG, &same));
1036: PetscCheck(same, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "TAO solver is not BNCG type");
1037: *type = cg->cg_type;
1038: PetscFunctionReturn(PETSC_SUCCESS);
1039: }
1041: /*@
1042: TaoBNCGSetType - Set the type for the `TAOBNCG` solver
1044: Input Parameters:
1045: + tao - the `Tao` solver context
1046: - type - `TAOBNCG` type
1048: Level: advanced
1050: .seealso: `Tao`, `TAOBNCG`, `TaoBNCGGetType()`, `TaoBNCGType`
1051: @*/
1052: PetscErrorCode TaoBNCGSetType(Tao tao, TaoBNCGType type)
1053: {
1054: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1055: PetscBool same;
1057: PetscFunctionBegin;
1058: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBNCG, &same));
1059: if (same) cg->cg_type = type;
1060: PetscFunctionReturn(PETSC_SUCCESS);
1061: }