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);