xref: /petsc/src/tao/unconstrained/impls/nls/nlsimpl.h (revision 1daac38e32cc923d29756e34a751df471eb1037d)
1*1daac38eSTodd Munson /*
2*1daac38eSTodd Munson Context for a Newton line search method (unconstrained minimization)
3*1daac38eSTodd Munson */
4*1daac38eSTodd Munson 
5*1daac38eSTodd Munson #if !defined(__TAO_NLS_H)
6*1daac38eSTodd Munson #define __TAO_NLS_H
7*1daac38eSTodd Munson #include <petsc/private/taoimpl.h>
8*1daac38eSTodd Munson 
9*1daac38eSTodd Munson typedef struct {
10*1daac38eSTodd Munson   Mat M;
11*1daac38eSTodd Munson 
12*1daac38eSTodd Munson   Vec D;
13*1daac38eSTodd Munson   Vec W;
14*1daac38eSTodd Munson 
15*1daac38eSTodd Munson   Vec Xold;
16*1daac38eSTodd Munson   Vec Gold;
17*1daac38eSTodd Munson   Vec Diag;
18*1daac38eSTodd Munson 
19*1daac38eSTodd Munson   /* Parameters when updating the perturbation added to the Hessian matrix
20*1daac38eSTodd Munson      according to the following scheme:
21*1daac38eSTodd Munson 
22*1daac38eSTodd Munson      pert = sval;
23*1daac38eSTodd Munson 
24*1daac38eSTodd Munson      do until convergence
25*1daac38eSTodd Munson        shift Hessian by pert
26*1daac38eSTodd Munson        solve Newton system
27*1daac38eSTodd Munson 
28*1daac38eSTodd Munson        if (linear solver failed or did not compute a descent direction)
29*1daac38eSTodd Munson          use steepest descent direction and increase perturbation
30*1daac38eSTodd Munson 
31*1daac38eSTodd Munson          if (0 == pert)
32*1daac38eSTodd Munson            initialize perturbation
33*1daac38eSTodd Munson            pert = min(imax, max(imin, imfac * norm(G)))
34*1daac38eSTodd Munson          else
35*1daac38eSTodd Munson            increase perturbation
36*1daac38eSTodd Munson            pert = min(pmax, max(pgfac * pert, pmgfac * norm(G)))
37*1daac38eSTodd Munson          fi
38*1daac38eSTodd Munson        else
39*1daac38eSTodd Munson          use linear solver direction and decrease perturbation
40*1daac38eSTodd Munson 
41*1daac38eSTodd Munson          pert = min(psfac * pert, pmsfac * norm(G))
42*1daac38eSTodd Munson          if (pert < pmin)
43*1daac38eSTodd Munson            pert = 0
44*1daac38eSTodd Munson          fi
45*1daac38eSTodd Munson        fi
46*1daac38eSTodd Munson 
47*1daac38eSTodd Munson        perform line search
48*1daac38eSTodd Munson        function and gradient evaluation
49*1daac38eSTodd Munson        check convergence
50*1daac38eSTodd Munson      od
51*1daac38eSTodd Munson   */
52*1daac38eSTodd Munson   PetscReal sval;               /*  Starting perturbation value, default zero */
53*1daac38eSTodd Munson 
54*1daac38eSTodd Munson   PetscReal imin;               /*  Minimum perturbation added during initialization  */
55*1daac38eSTodd Munson   PetscReal imax;               /*  Maximum perturbation added during initialization */
56*1daac38eSTodd Munson   PetscReal imfac;              /*  Merit function factor during initialization */
57*1daac38eSTodd Munson 
58*1daac38eSTodd Munson   PetscReal pmin;               /*  Minimim perturbation value */
59*1daac38eSTodd Munson   PetscReal pmax;               /*  Maximum perturbation value */
60*1daac38eSTodd Munson   PetscReal pgfac;              /*  Perturbation growth factor */
61*1daac38eSTodd Munson   PetscReal psfac;              /*  Perturbation shrink factor */
62*1daac38eSTodd Munson   PetscReal pmgfac;             /*  Merit function growth factor */
63*1daac38eSTodd Munson   PetscReal pmsfac;             /*  Merit function shrink factor */
64*1daac38eSTodd Munson 
65*1daac38eSTodd Munson   /* Parameters when updating the trust-region radius based on steplength
66*1daac38eSTodd Munson      if   step < nu1            (very bad step)
67*1daac38eSTodd Munson        radius = omega1 * min(norm(d), radius)
68*1daac38eSTodd Munson      elif step < nu2            (bad step)
69*1daac38eSTodd Munson        radius = omega2 * min(norm(d), radius)
70*1daac38eSTodd Munson      elif step < nu3            (okay step)
71*1daac38eSTodd Munson        radius = omega3 * radius;
72*1daac38eSTodd Munson      elif step < nu4            (good step)
73*1daac38eSTodd Munson        radius = max(omega4 * norm(d), radius)
74*1daac38eSTodd Munson      else                       (very good step)
75*1daac38eSTodd Munson        radius = max(omega5 * norm(d), radius)
76*1daac38eSTodd Munson      fi
77*1daac38eSTodd Munson   */
78*1daac38eSTodd Munson   PetscReal nu1;                /*  used to compute trust-region radius */
79*1daac38eSTodd Munson   PetscReal nu2;                /*  used to compute trust-region radius */
80*1daac38eSTodd Munson   PetscReal nu3;                /*  used to compute trust-region radius */
81*1daac38eSTodd Munson   PetscReal nu4;                /*  used to compute trust-region radius */
82*1daac38eSTodd Munson 
83*1daac38eSTodd Munson   PetscReal omega1;             /*  factor used for trust-region update */
84*1daac38eSTodd Munson   PetscReal omega2;             /*  factor used for trust-region update */
85*1daac38eSTodd Munson   PetscReal omega3;             /*  factor used for trust-region update */
86*1daac38eSTodd Munson   PetscReal omega4;             /*  factor used for trust-region update */
87*1daac38eSTodd Munson   PetscReal omega5;             /*  factor used for trust-region update */
88*1daac38eSTodd Munson 
89*1daac38eSTodd Munson   /* Parameters when updating the trust-region radius based on reduction
90*1daac38eSTodd Munson 
91*1daac38eSTodd Munson      kappa = ared / pred
92*1daac38eSTodd Munson      if   kappa < eta1          (very bad step)
93*1daac38eSTodd Munson        radius = alpha1 * min(norm(d), radius)
94*1daac38eSTodd Munson      elif kappa < eta2          (bad step)
95*1daac38eSTodd Munson        radius = alpha2 * min(norm(d), radius)
96*1daac38eSTodd Munson      elif kappa < eta3          (okay step)
97*1daac38eSTodd Munson        radius = alpha3 * radius;
98*1daac38eSTodd Munson      elif kappa < eta4          (good step)
99*1daac38eSTodd Munson        radius = max(alpha4 * norm(d), radius)
100*1daac38eSTodd Munson      else                       (very good step)
101*1daac38eSTodd Munson        radius = max(alpha5 * norm(d), radius)
102*1daac38eSTodd Munson      fi
103*1daac38eSTodd Munson   */
104*1daac38eSTodd Munson   PetscReal eta1;               /*  used to compute trust-region radius */
105*1daac38eSTodd Munson   PetscReal eta2;               /*  used to compute trust-region radius */
106*1daac38eSTodd Munson   PetscReal eta3;               /*  used to compute trust-region radius */
107*1daac38eSTodd Munson   PetscReal eta4;               /*  used to compute trust-region radius */
108*1daac38eSTodd Munson 
109*1daac38eSTodd Munson   PetscReal alpha1;             /*  factor used for trust-region update */
110*1daac38eSTodd Munson   PetscReal alpha2;             /*  factor used for trust-region update */
111*1daac38eSTodd Munson   PetscReal alpha3;             /*  factor used for trust-region update */
112*1daac38eSTodd Munson   PetscReal alpha4;             /*  factor used for trust-region update */
113*1daac38eSTodd Munson   PetscReal alpha5;             /*  factor used for trust-region update */
114*1daac38eSTodd Munson 
115*1daac38eSTodd Munson   /* Parameters when updating the trust-region radius based on interpolation
116*1daac38eSTodd Munson 
117*1daac38eSTodd Munson      kappa = ared / pred
118*1daac38eSTodd Munson      if   kappa >= 1.0 - mu1    (very good step)
119*1daac38eSTodd Munson        choose tau in [gamma3, gamma4]
120*1daac38eSTodd Munson        radius = max(tau * norm(d), radius)
121*1daac38eSTodd Munson      elif kappa >= 1.0 - mu2    (good step)
122*1daac38eSTodd Munson        choose tau in [gamma2, gamma3]
123*1daac38eSTodd Munson        if (tau >= 1.0)
124*1daac38eSTodd Munson          radius = max(tau * norm(d), radius)
125*1daac38eSTodd Munson        else
126*1daac38eSTodd Munson          radius = tau * min(norm(d), radius)
127*1daac38eSTodd Munson        fi
128*1daac38eSTodd Munson      else                       (bad step)
129*1daac38eSTodd Munson        choose tau in [gamma1, 1.0]
130*1daac38eSTodd Munson        radius = tau * min(norm(d), radius)
131*1daac38eSTodd Munson      fi
132*1daac38eSTodd Munson   */
133*1daac38eSTodd Munson   PetscReal mu1;                /*  used for model agreement in interpolation */
134*1daac38eSTodd Munson   PetscReal mu2;                /*  used for model agreement in interpolation */
135*1daac38eSTodd Munson 
136*1daac38eSTodd Munson   PetscReal gamma1;             /*  factor used for interpolation */
137*1daac38eSTodd Munson   PetscReal gamma2;             /*  factor used for interpolation */
138*1daac38eSTodd Munson   PetscReal gamma3;             /*  factor used for interpolation */
139*1daac38eSTodd Munson   PetscReal gamma4;             /*  factor used for interpolation */
140*1daac38eSTodd Munson 
141*1daac38eSTodd Munson   PetscReal theta;              /*  factor used for interpolation */
142*1daac38eSTodd Munson 
143*1daac38eSTodd Munson   /*  Parameters when initializing trust-region radius based on interpolation */
144*1daac38eSTodd Munson   PetscReal mu1_i;              /*  used for model agreement in interpolation */
145*1daac38eSTodd Munson   PetscReal mu2_i;              /*  used for model agreement in interpolation */
146*1daac38eSTodd Munson 
147*1daac38eSTodd Munson   PetscReal gamma1_i;           /*  factor used for interpolation */
148*1daac38eSTodd Munson   PetscReal gamma2_i;           /*  factor used for interpolation */
149*1daac38eSTodd Munson   PetscReal gamma3_i;           /*  factor used for interpolation */
150*1daac38eSTodd Munson   PetscReal gamma4_i;           /*  factor used for interpolation */
151*1daac38eSTodd Munson 
152*1daac38eSTodd Munson   PetscReal theta_i;            /*  factor used for interpolation */
153*1daac38eSTodd Munson 
154*1daac38eSTodd Munson   /*  Other parameters */
155*1daac38eSTodd Munson   PetscReal min_radius;         /*  lower bound on initial radius value */
156*1daac38eSTodd Munson   PetscReal max_radius;         /*  upper bound on trust region radius */
157*1daac38eSTodd Munson   PetscReal epsilon;            /*  tolerance used when computing ared/pred */
158*1daac38eSTodd Munson 
159*1daac38eSTodd Munson   PetscInt newt;                /*  Newton directions attempted */
160*1daac38eSTodd Munson   PetscInt bfgs;                /*  BFGS directions attempted */
161*1daac38eSTodd Munson   PetscInt sgrad;               /*  Scaled gradient directions attempted */
162*1daac38eSTodd Munson   PetscInt grad;                /*  Gradient directions attempted */
163*1daac38eSTodd Munson 
164*1daac38eSTodd Munson 
165*1daac38eSTodd Munson   PetscInt pc_type;             /*  Preconditioner for the code */
166*1daac38eSTodd Munson   PetscInt bfgs_scale_type;     /*  Scaling matrix to used for the bfgs preconditioner */
167*1daac38eSTodd Munson   PetscInt init_type;           /*  Trust-region initialization method */
168*1daac38eSTodd Munson   PetscInt update_type;         /*  Trust-region update method */
169*1daac38eSTodd Munson 
170*1daac38eSTodd Munson   PetscInt ksp_atol;
171*1daac38eSTodd Munson   PetscInt ksp_rtol;
172*1daac38eSTodd Munson   PetscInt ksp_ctol;
173*1daac38eSTodd Munson   PetscInt ksp_negc;
174*1daac38eSTodd Munson   PetscInt ksp_dtol;
175*1daac38eSTodd Munson   PetscInt ksp_iter;
176*1daac38eSTodd Munson   PetscInt ksp_othr;
177*1daac38eSTodd Munson } TAO_NLS;
178*1daac38eSTodd Munson 
179*1daac38eSTodd Munson #endif /* if !defined(__TAO_NLS_H) */
180