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 /* 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 *); 17 18 /* Embedded TAOBNCG */ 19 Tao bncg; 20 TAO_BNCG *bncg_ctx; 21 PetscInt max_cg_its, tot_cg_its; 22 Vec bncg_sol; 23 24 /* Allocated vectors */ 25 Vec W, Xwork, Gwork, Xold, Gold; 26 Vec unprojected_gradient, unprojected_gradient_old; 27 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; 32 33 /* Scalar values for the solution and step */ 34 PetscReal fold, f, gnorm, dnorm; 35 36 /* Parameters for active set estimation */ 37 PetscReal as_tol; 38 PetscReal as_step; 39 40 /* BFGS preconditioner data */ 41 PC bfgs_pre; 42 Mat M; 43 Vec Diag_min, Diag_max; 44 45 /* Parameters when updating the perturbation added to the Hessian matrix 46 according to the following scheme: 47 48 pert = sval; 49 50 do until convergence 51 shift Hessian by pert 52 solve Newton system 53 54 if (linear solver failed or did not compute a descent direction) 55 use steepest descent direction and increase perturbation 56 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 66 67 pert = min(psfac * pert, pmsfac * norm(G)) 68 if (pert < pmin) 69 pert = 0 70 fi 71 fi 72 73 perform line search 74 function and gradient evaluation 75 check convergence 76 od 77 */ 78 PetscReal sval; /* Starting perturbation value, default zero */ 79 80 PetscReal imin; /* Minimum perturbation added during initialization */ 81 PetscReal imax; /* Maximum perturbation added during initialization */ 82 PetscReal imfac; /* Merit function factor during initialization */ 83 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 */ 91 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 */ 109 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 */ 115 116 /* Parameters when updating the trust-region radius based on reduction 117 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 */ 135 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 */ 141 142 /* Parameters when updating the trust-region radius based on interpolation 143 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 */ 162 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 */ 167 168 PetscReal theta; /* factor used for interpolation */ 169 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 */ 173 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 */ 178 179 PetscReal theta_i; /* factor used for interpolation */ 180 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 */ 186 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 */ 191 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 resetksp; 206 207 /* Implementation specific context */ 208 void *ctx; 209 } TAO_BNK; 210 211 #define BNK_NEWTON 0 212 #define BNK_BFGS 1 213 #define BNK_SCALED_GRADIENT 2 214 #define BNK_GRADIENT 3 215 216 #define BNK_INIT_CONSTANT 0 217 #define BNK_INIT_DIRECTION 1 218 #define BNK_INIT_INTERPOLATION 2 219 #define BNK_INIT_TYPES 3 220 221 #define BNK_UPDATE_STEP 0 222 #define BNK_UPDATE_REDUCTION 1 223 #define BNK_UPDATE_INTERPOLATION 2 224 #define BNK_UPDATE_TYPES 3 225 226 #define BNK_AS_NONE 0 227 #define BNK_AS_BERTSEKAS 1 228 #define BNK_AS_TYPES 2 229 230 PETSC_INTERN PetscErrorCode TaoCreate_BNK(Tao); 231 PETSC_INTERN PetscErrorCode TaoSetUp_BNK(Tao); 232 PETSC_INTERN PetscErrorCode TaoSetFromOptions_BNK(Tao, PetscOptionItems *); 233 PETSC_INTERN PetscErrorCode TaoDestroy_BNK(Tao); 234 PETSC_INTERN PetscErrorCode TaoView_BNK(Tao, PetscViewer); 235 236 PETSC_INTERN PetscErrorCode TaoSolve_BNLS(Tao); 237 PETSC_INTERN PetscErrorCode TaoSolve_BNTR(Tao); 238 PETSC_INTERN PetscErrorCode TaoSolve_BNTL(Tao); 239 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); 252 253 #endif /* if !defined(__TAO_BNK_H) */ 254