1 /* 2 Context for bounded Newton-Krylov type optimization algorithms 3 */ 4 5 #pragma once 6 #include <petsc/private/taoimpl.h> 7 #include <../src/tao/bound/impls/bncg/bncg.h> 8 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 *); 16 17 /* Embedded TAOBNCG */ 18 Tao bncg; 19 TAO_BNCG *bncg_ctx; 20 PetscInt max_cg_its, tot_cg_its; 21 Vec bncg_sol; 22 23 /* Allocated vectors */ 24 Vec W, Xwork, Gwork, Xold, Gold; 25 Vec unprojected_gradient, unprojected_gradient_old; 26 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; 31 32 /* Scalar values for the solution and step */ 33 PetscReal fold, f, gnorm, dnorm; 34 35 /* Parameters for active set estimation */ 36 PetscReal as_tol; 37 PetscReal as_step; 38 39 /* BFGS preconditioner data */ 40 PC bfgs_pre; 41 Mat M; 42 Vec Diag_min, Diag_max; 43 44 /* Parameters when updating the perturbation added to the Hessian matrix 45 according to the following scheme: 46 47 pert = sval; 48 49 do until convergence 50 shift Hessian by pert 51 solve Newton system 52 53 if (linear solver failed or did not compute a descent direction) 54 use steepest descent direction and increase perturbation 55 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 65 66 pert = min(psfac * pert, pmsfac * norm(G)) 67 if (pert < pmin) 68 pert = 0 69 fi 70 fi 71 72 perform line search 73 function and gradient evaluation 74 check convergence 75 od 76 */ 77 PetscReal sval; /* Starting perturbation value, default zero */ 78 79 PetscReal imin; /* Minimum perturbation added during initialization */ 80 PetscReal imax; /* Maximum perturbation added during initialization */ 81 PetscReal imfac; /* Merit function factor during initialization */ 82 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 */ 90 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 */ 108 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 */ 114 115 /* Parameters when updating the trust-region radius based on reduction 116 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 */ 134 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 */ 140 141 /* Parameters when updating the trust-region radius based on interpolation 142 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 */ 161 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 */ 166 167 PetscReal theta; /* factor used for interpolation */ 168 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 */ 172 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 */ 177 178 PetscReal theta_i; /* factor used for interpolation */ 179 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 */ 185 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 */ 190 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 */ 195 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; 205 206 /* Implementation specific context */ 207 PetscCtx ctx; 208 } TAO_BNK; 209 210 #define BNK_NEWTON 0 211 #define BNK_BFGS 1 212 #define BNK_SCALED_GRADIENT 2 213 #define BNK_GRADIENT 3 214 215 #define BNK_INIT_CONSTANT 0 216 #define BNK_INIT_DIRECTION 1 217 #define BNK_INIT_INTERPOLATION 2 218 #define BNK_INIT_TYPES 3 219 220 #define BNK_UPDATE_STEP 0 221 #define BNK_UPDATE_REDUCTION 1 222 #define BNK_UPDATE_INTERPOLATION 2 223 #define BNK_UPDATE_TYPES 3 224 225 #define BNK_AS_NONE 0 226 #define BNK_AS_BERTSEKAS 1 227 #define BNK_AS_TYPES 2 228 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); 234 235 PETSC_INTERN PetscErrorCode TaoSolve_BNLS(Tao); 236 PETSC_INTERN PetscErrorCode TaoSolve_BNTR(Tao); 237 PETSC_INTERN PetscErrorCode TaoSolve_BNTL(Tao); 238 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); 251