Actual source code: bnk.h
1: /*
2: Context for bounded Newton-Krylov type optimization algorithms
3: */
5: #pragma once
6: #include <petsc/private/taoimpl.h>
7: #include <../src/tao/bound/impls/bncg/bncg.h>
9: typedef struct {
10: /* Function pointer for hessian evaluation
11: NOTE: This is necessary so that quasi-Newton-Krylov methods can "evaluate"
12: a quasi-Newton approximation while full Newton-Krylov methods call-back to
13: the application's Hessian */
14: PetscErrorCode (*computehessian)(Tao);
15: PetscErrorCode (*computestep)(Tao, PetscBool, KSPConvergedReason *, PetscInt *);
17: /* Embedded TAOBNCG */
18: Tao bncg;
19: TAO_BNCG *bncg_ctx;
20: PetscInt max_cg_its, tot_cg_its;
21: Vec bncg_sol;
23: /* Allocated vectors */
24: Vec W, Xwork, Gwork, Xold, Gold;
25: Vec unprojected_gradient, unprojected_gradient_old;
27: /* Unallocated matrices and vectors */
28: Mat H_inactive, Hpre_inactive;
29: Vec X_inactive, G_inactive, inactive_work, active_work;
30: IS inactive_idx, active_idx, active_lower, active_upper, active_fixed;
32: /* Scalar values for the solution and step */
33: PetscReal fold, f, gnorm, dnorm;
35: /* Parameters for active set estimation */
36: PetscReal as_tol;
37: PetscReal as_step;
39: /* BFGS preconditioner data */
40: PC bfgs_pre;
41: Mat M;
42: Vec Diag_min, Diag_max;
44: /* Parameters when updating the perturbation added to the Hessian matrix
45: according to the following scheme:
47: pert = sval;
49: do until convergence
50: shift Hessian by pert
51: solve Newton system
53: if (linear solver failed or did not compute a descent direction)
54: use steepest descent direction and increase perturbation
56: if (0 == pert)
57: initialize perturbation
58: pert = min(imax, max(imin, imfac * norm(G)))
59: else
60: increase perturbation
61: pert = min(pmax, max(pgfac * pert, pmgfac * norm(G)))
62: fi
63: else
64: use linear solver direction and decrease perturbation
66: pert = min(psfac * pert, pmsfac * norm(G))
67: if (pert < pmin)
68: pert = 0
69: fi
70: fi
72: perform line search
73: function and gradient evaluation
74: check convergence
75: od
76: */
77: PetscReal sval; /* Starting perturbation value, default zero */
79: PetscReal imin; /* Minimum perturbation added during initialization */
80: PetscReal imax; /* Maximum perturbation added during initialization */
81: PetscReal imfac; /* Merit function factor during initialization */
83: PetscReal pert; /* Current perturbation value */
84: PetscReal pmin; /* Minimum perturbation value */
85: PetscReal pmax; /* Maximum perturbation value */
86: PetscReal pgfac; /* Perturbation growth factor */
87: PetscReal psfac; /* Perturbation shrink factor */
88: PetscReal pmgfac; /* Merit function growth factor */
89: PetscReal pmsfac; /* Merit function shrink factor */
91: /* Parameters when updating the trust-region radius based on steplength
92: if step < nu1 (very bad step)
93: radius = omega1 * min(norm(d), radius)
94: elif step < nu2 (bad step)
95: radius = omega2 * min(norm(d), radius)
96: elif step < nu3 (okay step)
97: radius = omega3 * radius;
98: elif step < nu4 (good step)
99: radius = max(omega4 * norm(d), radius)
100: else (very good step)
101: radius = max(omega5 * norm(d), radius)
102: fi
103: */
104: PetscReal nu1; /* used to compute trust-region radius */
105: PetscReal nu2; /* used to compute trust-region radius */
106: PetscReal nu3; /* used to compute trust-region radius */
107: PetscReal nu4; /* used to compute trust-region radius */
109: PetscReal omega1; /* factor used for trust-region update */
110: PetscReal omega2; /* factor used for trust-region update */
111: PetscReal omega3; /* factor used for trust-region update */
112: PetscReal omega4; /* factor used for trust-region update */
113: PetscReal omega5; /* factor used for trust-region update */
115: /* Parameters when updating the trust-region radius based on reduction
117: kappa = ared / pred
118: if kappa < eta1 (very bad step)
119: radius = alpha1 * min(norm(d), radius)
120: elif kappa < eta2 (bad step)
121: radius = alpha2 * min(norm(d), radius)
122: elif kappa < eta3 (okay step)
123: radius = alpha3 * radius;
124: elif kappa < eta4 (good step)
125: radius = max(alpha4 * norm(d), radius)
126: else (very good step)
127: radius = max(alpha5 * norm(d), radius)
128: fi
129: */
130: PetscReal eta1; /* used to compute trust-region radius */
131: PetscReal eta2; /* used to compute trust-region radius */
132: PetscReal eta3; /* used to compute trust-region radius */
133: PetscReal eta4; /* used to compute trust-region radius */
135: PetscReal alpha1; /* factor used for trust-region update */
136: PetscReal alpha2; /* factor used for trust-region update */
137: PetscReal alpha3; /* factor used for trust-region update */
138: PetscReal alpha4; /* factor used for trust-region update */
139: PetscReal alpha5; /* factor used for trust-region update */
141: /* Parameters when updating the trust-region radius based on interpolation
143: kappa = ared / pred
144: if kappa >= 1.0 - mu1 (very good step)
145: choose tau in [gamma3, gamma4]
146: radius = max(tau * norm(d), radius)
147: elif kappa >= 1.0 - mu2 (good step)
148: choose tau in [gamma2, gamma3]
149: if (tau >= 1.0)
150: radius = max(tau * norm(d), radius)
151: else
152: radius = tau * min(norm(d), radius)
153: fi
154: else (bad step)
155: choose tau in [gamma1, 1.0]
156: radius = tau * min(norm(d), radius)
157: fi
158: */
159: PetscReal mu1; /* used for model agreement in interpolation */
160: PetscReal mu2; /* used for model agreement in interpolation */
162: PetscReal gamma1; /* factor used for interpolation */
163: PetscReal gamma2; /* factor used for interpolation */
164: PetscReal gamma3; /* factor used for interpolation */
165: PetscReal gamma4; /* factor used for interpolation */
167: PetscReal theta; /* factor used for interpolation */
169: /* Parameters when initializing trust-region radius based on interpolation */
170: PetscReal mu1_i; /* used for model agreement in interpolation */
171: PetscReal mu2_i; /* used for model agreement in interpolation */
173: PetscReal gamma1_i; /* factor used for interpolation */
174: PetscReal gamma2_i; /* factor used for interpolation */
175: PetscReal gamma3_i; /* factor used for interpolation */
176: PetscReal gamma4_i; /* factor used for interpolation */
178: PetscReal theta_i; /* factor used for interpolation */
180: /* Other parameters */
181: PetscReal min_radius; /* lower bound on initial radius value */
182: PetscReal max_radius; /* upper bound on trust region radius */
183: PetscReal epsilon; /* tolerance used when computing ared/pred */
184: PetscReal dmin, dmax; /* upper and lower bounds for the Hessian diagonal vector */
186: PetscInt newt; /* Newton directions attempted */
187: PetscInt bfgs; /* BFGS directions attempted */
188: PetscInt sgrad; /* Scaled gradient directions attempted */
189: PetscInt grad; /* Gradient directions attempted */
191: PetscInt as_type; /* Active set estimation method */
192: PetscInt bfgs_scale_type; /* Scaling matrix to used for the bfgs preconditioner */
193: PetscInt init_type; /* Trust-region initialization method */
194: PetscInt update_type; /* Trust-region update method */
196: /* Trackers for KSP solution type and convergence reasons */
197: PetscInt ksp_atol;
198: PetscInt ksp_rtol;
199: PetscInt ksp_ctol;
200: PetscInt ksp_negc;
201: PetscInt ksp_dtol;
202: PetscInt ksp_iter;
203: PetscInt ksp_othr;
204: PetscBool resetksp;
206: /* Implementation specific context */
207: void *ctx;
208: } TAO_BNK;
210: #define BNK_NEWTON 0
211: #define BNK_BFGS 1
212: #define BNK_SCALED_GRADIENT 2
213: #define BNK_GRADIENT 3
215: #define BNK_INIT_CONSTANT 0
216: #define BNK_INIT_DIRECTION 1
217: #define BNK_INIT_INTERPOLATION 2
218: #define BNK_INIT_TYPES 3
220: #define BNK_UPDATE_STEP 0
221: #define BNK_UPDATE_REDUCTION 1
222: #define BNK_UPDATE_INTERPOLATION 2
223: #define BNK_UPDATE_TYPES 3
225: #define BNK_AS_NONE 0
226: #define BNK_AS_BERTSEKAS 1
227: #define BNK_AS_TYPES 2
229: PETSC_INTERN PetscErrorCode TaoCreate_BNK(Tao);
230: PETSC_INTERN PetscErrorCode TaoSetUp_BNK(Tao);
231: PETSC_INTERN PetscErrorCode TaoSetFromOptions_BNK(Tao, PetscOptionItems *);
232: PETSC_INTERN PetscErrorCode TaoDestroy_BNK(Tao);
233: PETSC_INTERN PetscErrorCode TaoView_BNK(Tao, PetscViewer);
235: PETSC_INTERN PetscErrorCode TaoSolve_BNLS(Tao);
236: PETSC_INTERN PetscErrorCode TaoSolve_BNTR(Tao);
237: PETSC_INTERN PetscErrorCode TaoSolve_BNTL(Tao);
239: PETSC_INTERN PetscErrorCode TaoBNKPreconBFGS(PC, Vec, Vec);
240: PETSC_INTERN PetscErrorCode TaoBNKInitialize(Tao, PetscInt, PetscBool *);
241: PETSC_INTERN PetscErrorCode TaoBNKEstimateActiveSet(Tao, PetscInt);
242: PETSC_INTERN PetscErrorCode TaoBNKComputeHessian(Tao);
243: PETSC_INTERN PetscErrorCode TaoBNKBoundStep(Tao, PetscInt, Vec);
244: PETSC_INTERN PetscErrorCode TaoBNKTakeCGSteps(Tao, PetscBool *);
245: PETSC_INTERN PetscErrorCode TaoBNKComputeStep(Tao, PetscBool, KSPConvergedReason *, PetscInt *);
246: PETSC_INTERN PetscErrorCode TaoBNKRecomputePred(Tao, Vec, PetscReal *);
247: PETSC_INTERN PetscErrorCode TaoBNKSafeguardStep(Tao, KSPConvergedReason, PetscInt *);
248: PETSC_INTERN PetscErrorCode TaoBNKPerformLineSearch(Tao, PetscInt *, PetscReal *, TaoLineSearchConvergedReason *);
249: PETSC_INTERN PetscErrorCode TaoBNKUpdateTrustRadius(Tao, PetscReal, PetscReal, PetscInt, PetscInt, PetscBool *);
250: PETSC_INTERN PetscErrorCode TaoBNKAddStepCounts(Tao, PetscInt);