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: if (cg->pc) 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;
348: /* Override default settings (unless already changed) */
349: PetscCall(TaoParametersInitialize(tao));
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->pc = NULL;
371: cg->dk_eta = 0.5;
372: cg->hz_eta = 0.4;
373: cg->dynamic_restart = PETSC_FALSE;
374: cg->unscaled_restart = PETSC_FALSE;
375: cg->no_scaling = PETSC_FALSE;
376: cg->delta_min = 1e-7;
377: cg->delta_max = 100;
378: cg->theta = 1.0;
379: cg->hz_theta = 1.0;
380: cg->dfp_scale = 1.0;
381: cg->bfgs_scale = 1.0;
382: cg->zeta = 0.1;
383: cg->min_quad = 6;
384: cg->min_restart_num = 6; /* As in CG_DESCENT and KD2015*/
385: cg->xi = 1.0;
386: cg->neg_xi = PETSC_TRUE;
387: cg->spaced_restart = PETSC_FALSE;
388: cg->tol_quad = 1e-8;
389: cg->as_step = 0.001;
390: cg->as_tol = 0.001;
391: cg->eps_23 = PetscPowReal(PETSC_MACHINE_EPSILON, 2.0 / 3.0); /* Just a little tighter*/
392: cg->as_type = CG_AS_BERTSEKAS;
393: cg->cg_type = TAO_BNCG_SSML_BFGS;
394: cg->alpha = 1.0;
395: cg->diag_scaling = PETSC_TRUE;
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: PetscErrorCode TaoBNCGResetUpdate(Tao tao, PetscReal gnormsq)
400: {
401: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
402: PetscReal scaling;
404: PetscFunctionBegin;
405: ++cg->resets;
406: scaling = 2.0 * PetscMax(1.0, PetscAbsScalar(cg->f)) / PetscMax(gnormsq, cg->eps_23);
407: scaling = PetscMin(cg->delta_max, PetscMax(cg->delta_min, scaling));
408: if (cg->unscaled_restart) {
409: scaling = 1.0;
410: ++cg->pure_gd_steps;
411: }
412: PetscCall(VecAXPBY(tao->stepdirection, -scaling, 0.0, tao->gradient));
413: /* Also want to reset our diagonal scaling with each restart */
414: if (cg->diag_scaling) PetscCall(MatLMVMReset(cg->B, PETSC_FALSE));
415: PetscFunctionReturn(PETSC_SUCCESS);
416: }
418: PetscErrorCode TaoBNCGCheckDynamicRestart(Tao tao, PetscReal stepsize, PetscReal gd, PetscReal gd_old, PetscBool *dynrestart, PetscReal fold)
419: {
420: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
421: PetscReal quadinterp;
423: PetscFunctionBegin;
424: if (cg->f < cg->min_quad / 10) {
425: *dynrestart = PETSC_FALSE;
426: PetscFunctionReturn(PETSC_SUCCESS);
427: } /* just skip this since this strategy doesn't work well for functions near zero */
428: quadinterp = 2.0 * (cg->f - fold) / (stepsize * (gd + gd_old));
429: if (PetscAbs(quadinterp - 1.0) < cg->tol_quad) ++cg->iter_quad;
430: else {
431: cg->iter_quad = 0;
432: *dynrestart = PETSC_FALSE;
433: }
434: if (cg->iter_quad >= cg->min_quad) {
435: cg->iter_quad = 0;
436: *dynrestart = PETSC_TRUE;
437: }
438: PetscFunctionReturn(PETSC_SUCCESS);
439: }
441: PETSC_INTERN PetscErrorCode TaoBNCGStepDirectionUpdate(Tao tao, PetscReal gnorm2, PetscReal step, PetscReal fold, PetscReal gnorm2_old, PetscReal dnorm, PetscBool pcgd_fallback)
442: {
443: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
444: PetscReal gamma = 1.0, tau_k, beta;
445: PetscReal tmp = 1.0, ynorm, ynorm2 = 1.0, snorm = 1.0, dk_yk = 1.0, gd;
446: PetscReal gkp1_yk, gd_old, tau_bfgs, tau_dfp, gkp1D_yk, gtDg;
447: PetscInt dim;
448: PetscBool cg_restart = PETSC_FALSE;
450: PetscFunctionBegin;
451: /* Local curvature check to see if we need to restart */
452: if (tao->niter >= 1 || tao->recycle) {
453: PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
454: PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
455: ynorm2 = ynorm * ynorm;
456: PetscCall(VecDot(cg->yk, tao->stepdirection, &dk_yk));
457: if (step * dnorm < PETSC_MACHINE_EPSILON || step * dk_yk < PETSC_MACHINE_EPSILON) {
458: cg_restart = PETSC_TRUE;
459: ++cg->skipped_updates;
460: }
461: if (cg->spaced_restart) {
462: PetscCall(VecGetSize(tao->gradient, &dim));
463: if (tao->niter % (dim * cg->min_restart_num)) cg_restart = PETSC_TRUE;
464: }
465: }
466: /* If the user wants regular restarts, do it every 6n iterations, where n=dimension */
467: if (cg->spaced_restart) {
468: PetscCall(VecGetSize(tao->gradient, &dim));
469: if (0 == tao->niter % (6 * dim)) cg_restart = PETSC_TRUE;
470: }
471: /* Compute the diagonal scaling vector if applicable */
472: if (cg->diag_scaling) PetscCall(MatLMVMUpdate(cg->B, tao->solution, tao->gradient));
474: /* A note on diagonal scaling (to be added to paper):
475: For the FR, PR, PRP, and DY methods, the diagonally scaled versions
476: must be derived as a preconditioned CG method rather than as
477: a Hessian initialization like in the Broyden methods. */
479: /* In that case, one writes the objective function as
480: f(x) \equiv f(Ay). Gradient evaluations yield g(x_k) = A g(Ay_k) = A g(x_k).
481: Furthermore, the direction d_k \equiv (x_k - x_{k-1})/step according to
482: HZ (2006) becomes A^{-1} d_k, such that d_k^T g_k remains the
483: same under preconditioning. Note that A is diagonal, such that A^T = A. */
485: /* This yields questions like what the dot product d_k^T y_k
486: should look like. HZ mistakenly treats that as the same under
487: preconditioning, but that is not necessarily true. */
489: /* Observe y_k \equiv g_k - g_{k-1}, and under the P.C. transformation,
490: 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}),
491: 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
492: NOT the same if our matrix used to construct the preconditioner is updated between iterations.
493: This same issue is found when considering dot products of the form g_{k+1}^T y_k. */
495: /* Compute CG step direction */
496: if (cg_restart) {
497: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
498: } else if (pcgd_fallback) {
499: /* Just like preconditioned CG */
500: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
501: PetscCall(VecAXPBY(tao->stepdirection, -1.0, 0.0, cg->g_work));
502: } else if (ynorm2 > PETSC_MACHINE_EPSILON) {
503: switch (cg->cg_type) {
504: case TAO_BNCG_PCGD:
505: if (!cg->diag_scaling) {
506: if (!cg->no_scaling) {
507: cg->sts = step * step * dnorm * dnorm;
508: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->sts, &tau_k, cg->alpha));
509: } else {
510: tau_k = 1.0;
511: ++cg->pure_gd_steps;
512: }
513: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, 0.0, tao->gradient));
514: } else {
515: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
516: PetscCall(VecAXPBY(tao->stepdirection, -1.0, 0.0, cg->g_work));
517: }
518: break;
520: case TAO_BNCG_HS:
521: /* Classic Hestenes-Stiefel method, modified with scalar and diagonal preconditioning. */
522: if (!cg->diag_scaling) {
523: cg->sts = step * step * dnorm * dnorm;
524: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
525: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->sts, &tau_k, cg->alpha));
526: beta = tau_k * gkp1_yk / dk_yk;
527: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
528: } else {
529: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
530: PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
531: beta = gkp1_yk / dk_yk;
532: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
533: }
534: break;
536: case TAO_BNCG_FR:
537: PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
538: PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
539: PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
540: ynorm2 = ynorm * ynorm;
541: PetscCall(VecDot(cg->yk, tao->stepdirection, &dk_yk));
542: if (!cg->diag_scaling) {
543: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, step * step * dnorm * dnorm, &tau_k, cg->alpha));
544: beta = tau_k * gnorm2 / gnorm2_old;
545: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
546: } else {
547: PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old)); /* Before it's updated */
548: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
549: PetscCall(VecDot(tao->gradient, cg->g_work, &tmp));
550: beta = tmp / gnorm2_old;
551: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
552: }
553: break;
555: case TAO_BNCG_PRP:
556: snorm = step * dnorm;
557: if (!cg->diag_scaling) {
558: PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
559: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
560: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, snorm * snorm, &tau_k, cg->alpha));
561: beta = tau_k * gkp1_yk / gnorm2_old;
562: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
563: } else {
564: PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old));
565: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
566: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
567: beta = gkp1_yk / gnorm2_old;
568: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
569: }
570: break;
572: case TAO_BNCG_PRP_PLUS:
573: PetscCall(VecWAXPY(cg->yk, -1.0, cg->G_old, tao->gradient));
574: PetscCall(VecNorm(cg->yk, NORM_2, &ynorm));
575: ynorm2 = ynorm * ynorm;
576: if (!cg->diag_scaling) {
577: PetscCall(VecDot(cg->G_old, cg->G_old, &gnorm2_old));
578: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
579: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, snorm * snorm, &tau_k, cg->alpha));
580: beta = tau_k * gkp1_yk / gnorm2_old;
581: beta = PetscMax(beta, 0.0);
582: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
583: } else {
584: PetscCall(VecDot(cg->G_old, cg->g_work, &gnorm2_old)); /* Old gtDg */
585: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
586: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
587: beta = gkp1_yk / gnorm2_old;
588: beta = PetscMax(beta, 0.0);
589: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
590: }
591: break;
593: case TAO_BNCG_DY:
594: /* Dai, Yu-Hong, and Yaxiang Yuan. "A nonlinear conjugate gradient method with a strong global convergence property."
595: SIAM Journal on optimization 10, no. 1 (1999): 177-182. */
596: if (!cg->diag_scaling) {
597: PetscCall(VecDot(tao->stepdirection, tao->gradient, &gd));
598: PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
599: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, step * dk_yk, cg->yts, &tau_k, cg->alpha));
600: beta = tau_k * gnorm2 / (gd - gd_old);
601: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
602: } else {
603: PetscCall(MatMult(cg->B, tao->stepdirection, cg->d_work));
604: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
605: PetscCall(VecDot(cg->g_work, tao->gradient, >Dg));
606: PetscCall(VecDot(tao->stepdirection, cg->G_old, &gd_old));
607: PetscCall(VecDot(cg->d_work, cg->g_work, &dk_yk));
608: dk_yk = dk_yk - gd_old;
609: beta = gtDg / dk_yk;
610: PetscCall(VecScale(cg->d_work, beta));
611: PetscCall(VecWAXPY(tao->stepdirection, -1.0, cg->g_work, cg->d_work));
612: }
613: break;
615: case TAO_BNCG_HZ:
616: /* Hager, William W., and Hongchao Zhang. "Algorithm 851: CG_DESCENT, a conjugate gradient method with guaranteed descent."
617: ACM Transactions on Mathematical Software (TOMS) 32, no. 1 (2006): 113-137. */
618: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
619: PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
620: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
621: snorm = dnorm * step;
622: cg->yts = step * dk_yk;
623: if (cg->use_dynamic_restart) PetscCall(TaoBNCGCheckDynamicRestart(tao, step, gd, gd_old, &cg->dynamic_restart, fold));
624: if (cg->dynamic_restart) {
625: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
626: } else {
627: if (!cg->diag_scaling) {
628: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
629: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
630: /* Supplying cg->alpha = -1.0 will give the CG_DESCENT 5.3 special case of tau_k = 1.0 */
631: tmp = gd / dk_yk;
632: beta = tau_k * (gkp1_yk / dk_yk - ynorm2 * gd / (dk_yk * dk_yk));
633: /* Bound beta as in CG_DESCENT 5.3, as implemented, with the third comparison from DK 2013 */
634: beta = PetscMax(PetscMax(beta, cg->hz_eta * tau_k * gd_old / (dnorm * dnorm)), cg->dk_eta * tau_k * gd / (dnorm * dnorm));
635: /* d <- -t*g + beta*t*d */
636: PetscCall(VecAXPBY(tao->stepdirection, -tau_k, beta, tao->gradient));
637: } else {
638: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless BFGS step */
639: cg->yty = ynorm2;
640: cg->sts = snorm * snorm;
641: /* Apply the diagonal scaling to all my vectors */
642: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
643: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
644: PetscCall(MatSolve(cg->B, tao->stepdirection, cg->d_work));
645: /* Construct the constant ytDgkp1 */
646: PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
647: /* Construct the constant for scaling Dkyk in the update */
648: tmp = gd / dk_yk;
649: PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
650: tau_k = -tau_k * gd / (dk_yk * dk_yk);
651: /* beta is the constant which adds the dk contribution */
652: beta = gkp1_yk / dk_yk + cg->hz_theta * tau_k; /* HZ; (1.15) from DK 2013 */
653: /* From HZ2013, modified to account for diagonal scaling*/
654: PetscCall(VecDot(cg->G_old, cg->d_work, &gd_old));
655: PetscCall(VecDot(tao->stepdirection, cg->g_work, &gd));
656: beta = PetscMax(PetscMax(beta, cg->hz_eta * gd_old / (dnorm * dnorm)), cg->dk_eta * gd / (dnorm * dnorm));
657: /* Do the update */
658: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
659: }
660: }
661: break;
663: case TAO_BNCG_DK:
664: /* Dai, Yu-Hong, and Cai-Xia Kou. "A nonlinear conjugate gradient algorithm with an optimal property and an improved Wolfe line search."
665: SIAM Journal on Optimization 23, no. 1 (2013): 296-320. */
666: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
667: PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
668: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
669: snorm = step * dnorm;
670: cg->yts = dk_yk * step;
671: if (!cg->diag_scaling) {
672: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
673: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
674: /* Use cg->alpha = -1.0 to get tau_k = 1.0 as in CG_DESCENT 5.3 */
675: tmp = gd / dk_yk;
676: beta = tau_k * (gkp1_yk / dk_yk - ynorm2 * gd / (dk_yk * dk_yk) + gd / (dnorm * dnorm)) - step * gd / dk_yk;
677: beta = PetscMax(PetscMax(beta, cg->hz_eta * tau_k * gd_old / (dnorm * dnorm)), cg->dk_eta * tau_k * gd / (dnorm * dnorm));
678: /* d <- -t*g + beta*t*d */
679: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, 0.0, beta, tao->gradient, cg->yk));
680: } else {
681: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless BFGS step */
682: cg->yty = ynorm2;
683: cg->sts = snorm * snorm;
684: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
685: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
686: PetscCall(MatSolve(cg->B, tao->stepdirection, cg->d_work));
687: /* Construct the constant ytDgkp1 */
688: PetscCall(VecDot(cg->yk, cg->g_work, &gkp1_yk));
689: PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
690: tau_k = tau_k * gd / (dk_yk * dk_yk);
691: tmp = gd / dk_yk;
692: /* beta is the constant which adds the dk contribution */
693: beta = gkp1_yk / dk_yk - step * tmp - tau_k;
694: /* Update this for the last term in beta */
695: PetscCall(VecDot(cg->y_work, tao->stepdirection, &dk_yk));
696: beta += tmp * dk_yk / (dnorm * dnorm); /* projection of y_work onto dk */
697: PetscCall(VecDot(tao->stepdirection, cg->g_work, &gd));
698: PetscCall(VecDot(cg->G_old, cg->d_work, &gd_old));
699: beta = PetscMax(PetscMax(beta, cg->hz_eta * gd_old / (dnorm * dnorm)), cg->dk_eta * gd / (dnorm * dnorm));
700: /* Do the update */
701: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
702: }
703: break;
705: case TAO_BNCG_KD:
706: /* Kou, Cai-Xia, and Yu-Hong Dai. "A modified self-scaling memoryless Broyden-Fletcher-Goldfarb-Shanno method for unconstrained optimization."
707: Journal of Optimization Theory and Applications 165, no. 1 (2015): 209-224. */
708: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
709: PetscCall(VecDot(cg->G_old, tao->stepdirection, &gd_old));
710: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
711: snorm = step * dnorm;
712: cg->yts = dk_yk * step;
713: if (cg->use_dynamic_restart) PetscCall(TaoBNCGCheckDynamicRestart(tao, step, gd, gd_old, &cg->dynamic_restart, fold));
714: if (cg->dynamic_restart) {
715: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
716: } else {
717: if (!cg->diag_scaling) {
718: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
719: PetscCall(TaoBNCGComputeScalarScaling(ynorm2, cg->yts, snorm * snorm, &tau_k, cg->alpha));
720: beta = tau_k * (gkp1_yk / dk_yk - ynorm2 * gd / (dk_yk * dk_yk)) - step * gd / dk_yk;
721: if (beta < cg->zeta * tau_k * gd / (dnorm * dnorm)) /* 0.1 is KD's zeta parameter */
722: {
723: beta = cg->zeta * tau_k * gd / (dnorm * dnorm);
724: gamma = 0.0;
725: } else {
726: if (gkp1_yk < 0 && cg->neg_xi) gamma = -1.0 * gd / dk_yk;
727: /* This seems to be very effective when there's no tau_k scaling.
728: This guarantees a large descent step every iteration, going through DK 2015 Lemma 3.1's proof but allowing for negative xi */
729: else gamma = cg->xi * gd / dk_yk;
730: }
731: /* d <- -t*g + beta*t*d + t*tmp*yk */
732: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, gamma * tau_k, beta, tao->gradient, cg->yk));
733: } else {
734: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless BFGS step */
735: cg->yty = ynorm2;
736: cg->sts = snorm * snorm;
737: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
738: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
739: /* Construct the constant ytDgkp1 */
740: PetscCall(VecDot(cg->yk, cg->g_work, &gkp1D_yk));
741: /* Construct the constant for scaling Dkyk in the update */
742: gamma = gd / dk_yk;
743: /* tau_k = -ytDy/(ytd)^2 * gd */
744: PetscCall(VecDot(cg->yk, cg->y_work, &tau_k));
745: tau_k = tau_k * gd / (dk_yk * dk_yk);
746: /* beta is the constant which adds the d_k contribution */
747: beta = gkp1D_yk / dk_yk - step * gamma - tau_k;
748: /* Here is the requisite check */
749: PetscCall(VecDot(tao->stepdirection, cg->g_work, &tmp));
750: if (cg->neg_xi) {
751: /* modified KD implementation */
752: if (gkp1D_yk / dk_yk < 0) gamma = -1.0 * gd / dk_yk;
753: else gamma = cg->xi * gd / dk_yk;
754: if (beta < cg->zeta * tmp / (dnorm * dnorm)) {
755: beta = cg->zeta * tmp / (dnorm * dnorm);
756: gamma = 0.0;
757: }
758: } else { /* original KD 2015 implementation */
759: if (beta < cg->zeta * tmp / (dnorm * dnorm)) {
760: beta = cg->zeta * tmp / (dnorm * dnorm);
761: gamma = 0.0;
762: } else gamma = cg->xi * gd / dk_yk;
763: }
764: /* Do the update in two steps */
765: PetscCall(VecAXPBY(tao->stepdirection, -1.0, beta, cg->g_work));
766: PetscCall(VecAXPY(tao->stepdirection, gamma, cg->y_work));
767: }
768: }
769: break;
771: case TAO_BNCG_SSML_BFGS:
772: /* Perry, J. M. "A class of conjugate gradient algorithms with a two-step variable-metric memory."
773: Discussion Papers 269 (1977). */
774: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
775: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
776: snorm = step * dnorm;
777: cg->yts = dk_yk * step;
778: cg->yty = ynorm2;
779: cg->sts = snorm * snorm;
780: if (!cg->diag_scaling) {
781: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
782: PetscCall(TaoBNCGComputeScalarScaling(cg->yty, cg->yts, cg->sts, &tau_k, cg->alpha));
783: tmp = gd / dk_yk;
784: beta = tau_k * (gkp1_yk / dk_yk - cg->yty * gd / (dk_yk * dk_yk)) - step * tmp;
785: /* d <- -t*g + beta*t*d + t*tmp*yk */
786: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp * tau_k, beta, tao->gradient, cg->yk));
787: } else {
788: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless BFGS step */
789: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
790: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
791: /* compute scalar gamma */
792: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
793: PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
794: gamma = gd / dk_yk;
795: /* Compute scalar beta */
796: beta = (gkp1_yk / dk_yk - gd * tmp / (dk_yk * dk_yk)) - step * gd / dk_yk;
797: /* Compute stepdirection d_kp1 = gamma*Dkyk + beta*dk - Dkgkp1 */
798: PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
799: }
800: break;
802: case TAO_BNCG_SSML_DFP:
803: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
804: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
805: snorm = step * dnorm;
806: cg->yts = dk_yk * step;
807: cg->yty = ynorm2;
808: cg->sts = snorm * snorm;
809: if (!cg->diag_scaling) {
810: /* Instead of a regular convex combination, we will solve a quadratic formula. */
811: PetscCall(TaoBNCGComputeScalarScaling(cg->yty, cg->yts, cg->sts, &tau_k, cg->alpha));
812: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
813: tau_k = cg->dfp_scale * tau_k;
814: tmp = tau_k * gkp1_yk / cg->yty;
815: beta = -step * gd / dk_yk;
816: /* d <- -t*g + beta*d + tmp*yk */
817: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp, beta, tao->gradient, cg->yk));
818: } else {
819: /* We have diagonal scaling enabled and are taking a diagonally-scaled memoryless DFP step */
820: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
821: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
822: /* compute scalar gamma */
823: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
824: PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
825: gamma = (gkp1_yk / tmp);
826: /* Compute scalar beta */
827: beta = -step * gd / dk_yk;
828: /* Compute stepdirection d_kp1 = gamma*Dkyk + beta*dk - Dkgkp1 */
829: PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
830: }
831: break;
833: case TAO_BNCG_SSML_BRDN:
834: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
835: PetscCall(VecWAXPY(cg->sk, -1.0, cg->X_old, tao->solution));
836: snorm = step * dnorm;
837: cg->yts = step * dk_yk;
838: cg->yty = ynorm2;
839: cg->sts = snorm * snorm;
840: if (!cg->diag_scaling) {
841: /* Instead of a regular convex combination, we will solve a quadratic formula. */
842: PetscCall(TaoBNCGComputeScalarScaling(cg->yty, step * dk_yk, snorm * snorm, &tau_bfgs, cg->bfgs_scale));
843: PetscCall(TaoBNCGComputeScalarScaling(cg->yty, step * dk_yk, snorm * snorm, &tau_dfp, cg->dfp_scale));
844: PetscCall(VecDot(cg->yk, tao->gradient, &gkp1_yk));
845: tau_k = cg->theta * tau_bfgs + (1.0 - cg->theta) * tau_dfp;
846: /* If bfgs_scale = 1.0, it should reproduce the bfgs tau_bfgs. If bfgs_scale = 0.0,
847: it should reproduce the tau_dfp scaling. Same with dfp_scale. */
848: tmp = cg->theta * tau_bfgs * gd / dk_yk + (1 - cg->theta) * tau_dfp * gkp1_yk / cg->yty;
849: beta = cg->theta * tau_bfgs * (gkp1_yk / dk_yk - cg->yty * gd / (dk_yk * dk_yk)) - step * gd / dk_yk;
850: /* d <- -t*g + beta*d + tmp*yk */
851: PetscCall(VecAXPBYPCZ(tao->stepdirection, -tau_k, tmp, beta, tao->gradient, cg->yk));
852: } else {
853: /* We have diagonal scaling enabled */
854: PetscCall(MatSolve(cg->B, tao->gradient, cg->g_work));
855: PetscCall(MatSolve(cg->B, cg->yk, cg->y_work));
856: /* compute scalar gamma */
857: PetscCall(VecDot(cg->g_work, cg->yk, &gkp1_yk));
858: PetscCall(VecDot(cg->y_work, cg->yk, &tmp));
859: gamma = cg->theta * gd / dk_yk + (1 - cg->theta) * (gkp1_yk / tmp);
860: /* Compute scalar beta */
861: beta = cg->theta * (gkp1_yk / dk_yk - gd * tmp / (dk_yk * dk_yk)) - step * gd / dk_yk;
862: /* Compute stepdirection dkp1 = gamma*Dkyk + beta*dk - Dkgkp1 */
863: PetscCall(VecAXPBYPCZ(tao->stepdirection, -1.0, gamma, beta, cg->g_work, cg->y_work));
864: }
865: break;
867: default:
868: break;
869: }
870: }
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: PETSC_INTERN PetscErrorCode TaoBNCGConductIteration(Tao tao, PetscReal gnorm)
875: {
876: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
877: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
878: PetscReal step = 1.0, gnorm2, gd, dnorm = 0.0;
879: PetscReal gnorm2_old, f_old, resnorm, gnorm_old;
880: PetscBool pcgd_fallback = PETSC_FALSE;
882: PetscFunctionBegin;
883: /* We are now going to perform a line search along the direction. */
884: /* Store solution and gradient info before it changes */
885: PetscCall(VecCopy(tao->solution, cg->X_old));
886: PetscCall(VecCopy(tao->gradient, cg->G_old));
887: PetscCall(VecCopy(cg->unprojected_gradient, cg->unprojected_gradient_old));
889: gnorm_old = gnorm;
890: gnorm2_old = gnorm_old * gnorm_old;
891: f_old = cg->f;
892: /* Perform bounded line search. If we are recycling a solution from a previous */
893: /* TaoSolve, then we want to immediately skip to calculating a new direction rather than performing a linesearch */
894: if (!(tao->recycle && 0 == tao->niter)) {
895: /* Above logic: the below code happens every iteration, except for the first iteration of a recycled TaoSolve */
896: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
897: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
898: PetscCall(TaoAddLineSearchCounts(tao));
900: /* Check linesearch failure */
901: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
902: ++cg->ls_fails;
903: if (cg->cg_type == TAO_BNCG_GD) {
904: /* Nothing left to do but fail out of the optimization */
905: step = 0.0;
906: tao->reason = TAO_DIVERGED_LS_FAILURE;
907: } else {
908: /* Restore previous point, perform preconditioned GD and regular GD steps at the last good point */
909: PetscCall(VecCopy(cg->X_old, tao->solution));
910: PetscCall(VecCopy(cg->G_old, tao->gradient));
911: PetscCall(VecCopy(cg->unprojected_gradient_old, cg->unprojected_gradient));
912: gnorm = gnorm_old;
913: gnorm2 = gnorm2_old;
914: cg->f = f_old;
916: /* Fall back on preconditioned CG (so long as you're not already using it) */
917: if (cg->cg_type != TAO_BNCG_PCGD && cg->diag_scaling) {
918: pcgd_fallback = PETSC_TRUE;
919: PetscCall(TaoBNCGStepDirectionUpdate(tao, gnorm2, step, f_old, gnorm2_old, dnorm, pcgd_fallback));
921: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
922: PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
924: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
925: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
926: PetscCall(TaoAddLineSearchCounts(tao));
928: pcgd_fallback = PETSC_FALSE;
929: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
930: /* Going to perform a regular gradient descent step. */
931: ++cg->ls_fails;
932: step = 0.0;
933: }
934: }
935: /* Fall back on the scaled gradient step */
936: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
937: ++cg->ls_fails;
938: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
939: PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
940: PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0));
941: PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &cg->f, cg->unprojected_gradient, tao->stepdirection, &step, &ls_status));
942: PetscCall(TaoAddLineSearchCounts(tao));
943: }
945: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
946: /* Nothing left to do but fail out of the optimization */
947: ++cg->ls_fails;
948: step = 0.0;
949: tao->reason = TAO_DIVERGED_LS_FAILURE;
950: } else {
951: /* One of the fallbacks worked. Set them both back equal to false. */
952: pcgd_fallback = PETSC_FALSE;
953: }
954: }
955: }
956: /* Convergence test for line search failure */
957: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
959: /* Standard convergence test */
960: PetscCall(VecFischer(tao->solution, cg->unprojected_gradient, tao->XL, tao->XU, cg->W));
961: PetscCall(VecNorm(cg->W, NORM_2, &resnorm));
962: PetscCheck(!PetscIsInfOrNanReal(resnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated infinity or NaN");
963: PetscCall(TaoLogConvergenceHistory(tao, cg->f, resnorm, 0.0, tao->ksp_its));
964: PetscCall(TaoMonitor(tao, tao->niter, cg->f, resnorm, 0.0, step));
965: PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
966: if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);
967: }
968: /* Assert we have an updated step and we need at least one more iteration. */
969: /* Calculate the next direction */
970: /* Estimate the active set at the new solution */
971: PetscCall(TaoBNCGEstimateActiveSet(tao, cg->as_type));
972: /* Compute the projected gradient and its norm */
973: PetscCall(VecCopy(cg->unprojected_gradient, tao->gradient));
974: if (cg->active_idx) PetscCall(VecISSet(tao->gradient, cg->active_idx, 0.0));
975: PetscCall(VecNorm(tao->gradient, NORM_2, &gnorm));
976: gnorm2 = gnorm * gnorm;
978: /* Calculate some quantities used in the StepDirectionUpdate. */
979: PetscCall(VecNorm(tao->stepdirection, NORM_2, &dnorm));
980: /* Update the step direction. */
981: PetscCall(TaoBNCGStepDirectionUpdate(tao, gnorm2, step, f_old, gnorm2_old, dnorm, pcgd_fallback));
982: ++tao->niter;
983: PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
985: if (cg->cg_type != TAO_BNCG_GD) {
986: /* Figure out which previously active variables became inactive this iteration */
987: PetscCall(ISDestroy(&cg->new_inactives));
988: if (cg->inactive_idx && cg->inactive_old) PetscCall(ISDifference(cg->inactive_idx, cg->inactive_old, &cg->new_inactives));
989: /* Selectively reset the CG step those freshly inactive variables */
990: if (cg->new_inactives) {
991: PetscCall(VecGetSubVector(tao->stepdirection, cg->new_inactives, &cg->inactive_step));
992: PetscCall(VecGetSubVector(cg->unprojected_gradient, cg->new_inactives, &cg->inactive_grad));
993: PetscCall(VecCopy(cg->inactive_grad, cg->inactive_step));
994: PetscCall(VecScale(cg->inactive_step, -1.0));
995: PetscCall(VecRestoreSubVector(tao->stepdirection, cg->new_inactives, &cg->inactive_step));
996: PetscCall(VecRestoreSubVector(cg->unprojected_gradient, cg->new_inactives, &cg->inactive_grad));
997: }
998: /* Verify that this is a descent direction */
999: PetscCall(VecDot(tao->gradient, tao->stepdirection, &gd));
1000: PetscCall(VecNorm(tao->stepdirection, NORM_2, &dnorm));
1001: if (PetscIsInfOrNanReal(gd) || (gd / (dnorm * dnorm) <= -1e10 || gd / (dnorm * dnorm) >= -1e-10)) {
1002: /* Not a descent direction, so we reset back to projected gradient descent */
1003: PetscCall(TaoBNCGResetUpdate(tao, gnorm2));
1004: PetscCall(TaoBNCGBoundStep(tao, cg->as_type, tao->stepdirection));
1005: ++cg->descent_error;
1006: } else {
1007: }
1008: }
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }
1012: PETSC_INTERN PetscErrorCode TaoBNCGSetH0(Tao tao, Mat H0)
1013: {
1014: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1015: PetscBool same;
1017: PetscFunctionBegin;
1018: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBNCG, &same));
1019: if (same) {
1020: PetscCall(PetscObjectReference((PetscObject)H0));
1021: cg->pc = H0;
1022: }
1023: PetscFunctionReturn(PETSC_SUCCESS);
1024: }
1026: /*@
1027: TaoBNCGGetType - Return the type for the `TAOBNCG` solver
1029: Input Parameter:
1030: . tao - the `Tao` solver context
1032: Output Parameter:
1033: . type - `TAOBNCG` type
1035: Level: advanced
1037: .seealso: `Tao`, `TAOBNCG`, `TaoBNCGSetType()`, `TaoBNCGType`
1038: @*/
1039: PetscErrorCode TaoBNCGGetType(Tao tao, TaoBNCGType *type)
1040: {
1041: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1042: PetscBool same;
1044: PetscFunctionBegin;
1045: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBNCG, &same));
1046: PetscCheck(same, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "TAO solver is not BNCG type");
1047: *type = cg->cg_type;
1048: PetscFunctionReturn(PETSC_SUCCESS);
1049: }
1051: /*@
1052: TaoBNCGSetType - Set the type for the `TAOBNCG` solver
1054: Input Parameters:
1055: + tao - the `Tao` solver context
1056: - type - `TAOBNCG` type
1058: Level: advanced
1060: .seealso: `Tao`, `TAOBNCG`, `TaoBNCGGetType()`, `TaoBNCGType`
1061: @*/
1062: PetscErrorCode TaoBNCGSetType(Tao tao, TaoBNCGType type)
1063: {
1064: TAO_BNCG *cg = (TAO_BNCG *)tao->data;
1065: PetscBool same;
1067: PetscFunctionBegin;
1068: PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBNCG, &same));
1069: if (same) cg->cg_type = type;
1070: PetscFunctionReturn(PETSC_SUCCESS);
1071: }