Actual source code: bnk.h
1: /*
2: Context for bounded Newton-Krylov type optimization algorithms
3: */
7: #include <petsc/private/taoimpl.h>
8: #include <../src/tao/bound/impls/bncg/bncg.h>
10: typedef struct {
11: /* Function pointer for hessian evaluation
12: NOTE: This is necessary so that quasi-Newton-Krylov methods can "evaluate"
13: a quasi-Newton approximation while full Newton-Krylov methods call-back to
14: the application's Hessian */
15: PetscErrorCode (*computehessian)(Tao);
16: PetscErrorCode (*computestep)(Tao, PetscBool, KSPConvergedReason*, PetscInt*);
18: /* Embedded TAOBNCG */
19: Tao bncg;
20: TAO_BNCG *bncg_ctx;
21: PetscInt max_cg_its, tot_cg_its;
22: Vec bncg_sol;
24: /* Allocated vectors */
25: Vec W, Xwork, Gwork, Xold, Gold;
26: Vec unprojected_gradient, unprojected_gradient_old;
28: /* Unallocated matrices and vectors */
29: Mat H_inactive, Hpre_inactive;
30: Vec X_inactive, G_inactive, inactive_work, active_work;
31: IS inactive_idx, active_idx, active_lower, active_upper, active_fixed;
33: /* Scalar values for the solution and step */
34: PetscReal fold, f, gnorm, dnorm;
36: /* Parameters for active set estimation */
37: PetscReal as_tol;
38: PetscReal as_step;
40: /* BFGS preconditioner data */
41: PC bfgs_pre;
42: Mat M;
43: Vec Diag_min, Diag_max;
45: /* Parameters when updating the perturbation added to the Hessian matrix
46: according to the following scheme:
48: pert = sval;
50: do until convergence
51: shift Hessian by pert
52: solve Newton system
54: if (linear solver failed or did not compute a descent direction)
55: use steepest descent direction and increase perturbation
57: if (0 == pert)
58: initialize perturbation
59: pert = min(imax, max(imin, imfac * norm(G)))
60: else
61: increase perturbation
62: pert = min(pmax, max(pgfac * pert, pmgfac * norm(G)))
63: fi
64: else
65: use linear solver direction and decrease perturbation
67: pert = min(psfac * pert, pmsfac * norm(G))
68: if (pert < pmin)
69: pert = 0
70: fi
71: fi
73: perform line search
74: function and gradient evaluation
75: check convergence
76: od
77: */
78: PetscReal sval; /* Starting perturbation value, default zero */
80: PetscReal imin; /* Minimum perturbation added during initialization */
81: PetscReal imax; /* Maximum perturbation added during initialization */
82: PetscReal imfac; /* Merit function factor during initialization */
84: PetscReal pert; /* Current perturbation value */
85: PetscReal pmin; /* Minimim perturbation value */
86: PetscReal pmax; /* Maximum perturbation value */
87: PetscReal pgfac; /* Perturbation growth factor */
88: PetscReal psfac; /* Perturbation shrink factor */
89: PetscReal pmgfac; /* Merit function growth factor */
90: PetscReal pmsfac; /* Merit function shrink factor */
92: /* Parameters when updating the trust-region radius based on steplength
93: if step < nu1 (very bad step)
94: radius = omega1 * min(norm(d), radius)
95: elif step < nu2 (bad step)
96: radius = omega2 * min(norm(d), radius)
97: elif step < nu3 (okay step)
98: radius = omega3 * radius;
99: elif step < nu4 (good step)
100: radius = max(omega4 * norm(d), radius)
101: else (very good step)
102: radius = max(omega5 * norm(d), radius)
103: fi
104: */
105: PetscReal nu1; /* used to compute trust-region radius */
106: PetscReal nu2; /* used to compute trust-region radius */
107: PetscReal nu3; /* used to compute trust-region radius */
108: PetscReal nu4; /* used to compute trust-region radius */
110: PetscReal omega1; /* factor used for trust-region update */
111: PetscReal omega2; /* factor used for trust-region update */
112: PetscReal omega3; /* factor used for trust-region update */
113: PetscReal omega4; /* factor used for trust-region update */
114: PetscReal omega5; /* factor used for trust-region update */
116: /* Parameters when updating the trust-region radius based on reduction
118: kappa = ared / pred
119: if kappa < eta1 (very bad step)
120: radius = alpha1 * min(norm(d), radius)
121: elif kappa < eta2 (bad step)
122: radius = alpha2 * min(norm(d), radius)
123: elif kappa < eta3 (okay step)
124: radius = alpha3 * radius;
125: elif kappa < eta4 (good step)
126: radius = max(alpha4 * norm(d), radius)
127: else (very good step)
128: radius = max(alpha5 * norm(d), radius)
129: fi
130: */
131: PetscReal eta1; /* used to compute trust-region radius */
132: PetscReal eta2; /* used to compute trust-region radius */
133: PetscReal eta3; /* used to compute trust-region radius */
134: PetscReal eta4; /* used to compute trust-region radius */
136: PetscReal alpha1; /* factor used for trust-region update */
137: PetscReal alpha2; /* factor used for trust-region update */
138: PetscReal alpha3; /* factor used for trust-region update */
139: PetscReal alpha4; /* factor used for trust-region update */
140: PetscReal alpha5; /* factor used for trust-region update */
142: /* Parameters when updating the trust-region radius based on interpolation
144: kappa = ared / pred
145: if kappa >= 1.0 - mu1 (very good step)
146: choose tau in [gamma3, gamma4]
147: radius = max(tau * norm(d), radius)
148: elif kappa >= 1.0 - mu2 (good step)
149: choose tau in [gamma2, gamma3]
150: if (tau >= 1.0)
151: radius = max(tau * norm(d), radius)
152: else
153: radius = tau * min(norm(d), radius)
154: fi
155: else (bad step)
156: choose tau in [gamma1, 1.0]
157: radius = tau * min(norm(d), radius)
158: fi
159: */
160: PetscReal mu1; /* used for model agreement in interpolation */
161: PetscReal mu2; /* used for model agreement in interpolation */
163: PetscReal gamma1; /* factor used for interpolation */
164: PetscReal gamma2; /* factor used for interpolation */
165: PetscReal gamma3; /* factor used for interpolation */
166: PetscReal gamma4; /* factor used for interpolation */
168: PetscReal theta; /* factor used for interpolation */
170: /* Parameters when initializing trust-region radius based on interpolation */
171: PetscReal mu1_i; /* used for model agreement in interpolation */
172: PetscReal mu2_i; /* used for model agreement in interpolation */
174: PetscReal gamma1_i; /* factor used for interpolation */
175: PetscReal gamma2_i; /* factor used for interpolation */
176: PetscReal gamma3_i; /* factor used for interpolation */
177: PetscReal gamma4_i; /* factor used for interpolation */
179: PetscReal theta_i; /* factor used for interpolation */
181: /* Other parameters */
182: PetscReal min_radius; /* lower bound on initial radius value */
183: PetscReal max_radius; /* upper bound on trust region radius */
184: PetscReal epsilon; /* tolerance used when computing ared/pred */
185: PetscReal dmin, dmax; /* upper and lower bounds for the Hessian diagonal vector */
187: PetscInt newt; /* Newton directions attempted */
188: PetscInt bfgs; /* BFGS directions attempted */
189: PetscInt sgrad; /* Scaled gradient directions attempted */
190: PetscInt grad; /* Gradient directions attempted */
192: PetscInt as_type; /* Active set estimation method */
193: PetscInt bfgs_scale_type; /* Scaling matrix to used for the bfgs preconditioner */
194: PetscInt init_type; /* Trust-region initialization method */
195: PetscInt update_type; /* Trust-region update method */
197: /* Trackers for KSP solution type and convergence reasons */
198: PetscInt ksp_atol;
199: PetscInt ksp_rtol;
200: PetscInt ksp_ctol;
201: PetscInt ksp_negc;
202: PetscInt ksp_dtol;
203: PetscInt ksp_iter;
204: PetscInt ksp_othr;
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(PetscOptionItems*, Tao);
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);
252: #endif /* if !defined(__TAO_BNK_H) */