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