Actual source code: bnk.h

petsc-3.12.5 2020-03-29
Report Typos and Errors
  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 */
196: 
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;
205:   PetscBool is_nash, is_stcg, is_gltr;

207:   /* Implementation specific context */
208:   void* ctx;
209: } TAO_BNK;

211: #define BNK_NEWTON              0
212: #define BNK_BFGS                1
213: #define BNK_SCALED_GRADIENT     2
214: #define BNK_GRADIENT            3

216: #define BNK_INIT_CONSTANT         0
217: #define BNK_INIT_DIRECTION        1
218: #define BNK_INIT_INTERPOLATION    2
219: #define BNK_INIT_TYPES            3

221: #define BNK_UPDATE_STEP           0
222: #define BNK_UPDATE_REDUCTION      1
223: #define BNK_UPDATE_INTERPOLATION  2
224: #define BNK_UPDATE_TYPES          3

226: #define BNK_AS_NONE        0
227: #define BNK_AS_BERTSEKAS   1
228: #define BNK_AS_TYPES       2

230: PETSC_INTERN PetscErrorCode TaoCreate_BNK(Tao);
231: PETSC_INTERN PetscErrorCode TaoSetUp_BNK(Tao);
232: PETSC_INTERN PetscErrorCode TaoSetFromOptions_BNK(PetscOptionItems*, Tao);
233: PETSC_INTERN PetscErrorCode TaoDestroy_BNK(Tao);
234: PETSC_INTERN PetscErrorCode TaoView_BNK(Tao, PetscViewer);

236: PETSC_INTERN PetscErrorCode TaoSolve_BNLS(Tao);
237: PETSC_INTERN PetscErrorCode TaoSolve_BNTR(Tao);
238: PETSC_INTERN PetscErrorCode TaoSolve_BNTL(Tao);

240: PETSC_INTERN PetscErrorCode TaoBNKPreconBFGS(PC, Vec, Vec);
241: PETSC_INTERN PetscErrorCode TaoBNKInitialize(Tao, PetscInt, PetscBool*);
242: PETSC_INTERN PetscErrorCode TaoBNKEstimateActiveSet(Tao, PetscInt);
243: PETSC_INTERN PetscErrorCode TaoBNKComputeHessian(Tao);
244: PETSC_INTERN PetscErrorCode TaoBNKBoundStep(Tao, PetscInt, Vec);
245: PETSC_INTERN PetscErrorCode TaoBNKTakeCGSteps(Tao, PetscBool*);
246: PETSC_INTERN PetscErrorCode TaoBNKComputeStep(Tao, PetscBool, KSPConvergedReason*, PetscInt*);
247: PETSC_INTERN PetscErrorCode TaoBNKRecomputePred(Tao, Vec, PetscReal*);
248: PETSC_INTERN PetscErrorCode TaoBNKSafeguardStep(Tao, KSPConvergedReason, PetscInt*);
249: PETSC_INTERN PetscErrorCode TaoBNKPerformLineSearch(Tao, PetscInt*, PetscReal*, TaoLineSearchConvergedReason*);
250: PETSC_INTERN PetscErrorCode TaoBNKUpdateTrustRadius(Tao, PetscReal, PetscReal, PetscInt, PetscInt, PetscBool*);
251: PETSC_INTERN PetscErrorCode TaoBNKAddStepCounts(Tao, PetscInt);

253: #endif /* if !defined(__TAO_BNK_H) */