1 #include <../src/tao/bound/impls/bnk/bnk.h> 2 #include <petscksp.h> 3 4 /* 5 Implements Newton's Method with a trust region approach for solving 6 bound constrained minimization problems. 7 8 ------------------------------------------------------------ 9 10 initialize trust radius (default: BNK_INIT_INTERPOLATION) 11 x_0 = VecMedian(x_0) 12 f_0, g_0 = TaoComputeObjectiveAndGradient(x_0) 13 pg_0 = VecBoundGradientProjection(g_0) 14 check convergence at pg_0 15 niter = 0 16 step_accepted = true 17 18 while niter <= max_it 19 if step_accepted 20 niter += 1 21 H_k = TaoComputeHessian(x_k) 22 if pc_type == BNK_PC_BFGS 23 add correction to BFGS approx 24 if scale_type == BNK_SCALE_AHESS 25 D = VecMedian(1e-6, abs(diag(H_k)), 1e6) 26 scale BFGS with VecReciprocal(D) 27 end 28 end 29 end 30 31 if pc_type = BNK_PC_BFGS 32 B_k = BFGS 33 else 34 B_k = VecMedian(1e-6, abs(diag(H_k)), 1e6) 35 B_k = VecReciprocal(B_k) 36 end 37 w = x_k - VecMedian(x_k - 0.001*B_k*g_k) 38 eps = min(eps, norm2(w)) 39 determine the active and inactive index sets such that 40 L = {i : (x_k)_i <= l_i + eps && (g_k)_i > 0} 41 U = {i : (x_k)_i >= u_i - eps && (g_k)_i < 0} 42 F = {i : l_i = (x_k)_i = u_i} 43 A = {L + U + F} 44 I = {i : i not in A} 45 46 generate the reduced system Hr_k dr_k = -gr_k for variables in I 47 if pc_type == BNK_PC_BFGS && scale_type == BNK_SCALE_PHESS 48 D = VecMedian(1e-6, abs(diag(Hr_k)), 1e6) 49 scale BFGS with VecReciprocal(D) 50 end 51 solve Hr_k dr_k = -gr_k 52 set d_k to (l - x) for variables in L, (u - x) for variables in U, and 0 for variables in F 53 54 x_{k+1} = VecMedian(x_k + d_k) 55 s = x_{k+1} - x_k 56 prered = dot(s, 0.5*gr_k - Hr_k*s) 57 f_{k+1} = TaoComputeObjective(x_{k+1}) 58 actred = f_k - f_{k+1} 59 60 oldTrust = trust 61 step_accepted, trust = TaoBNKUpdateTrustRadius(default: BNK_UPDATE_REDUCTION) 62 if step_accepted 63 g_{k+1} = TaoComputeGradient(x_{k+1}) 64 pg_{k+1} = VecBoundGradientProjection(g_{k+1}) 65 count the accepted Newton step 66 else 67 f_{k+1} = f_k 68 x_{k+1} = x_k 69 g_{k+1} = g_k 70 pg_{k+1} = pg_k 71 if trust == oldTrust 72 terminate because we cannot shrink the radius any further 73 end 74 end 75 76 check convergence at pg_{k+1} 77 end 78 */ 79 80 static PetscErrorCode TaoSolve_BNTR(Tao tao) 81 { 82 PetscErrorCode ierr; 83 TAO_BNK *bnk = (TAO_BNK *)tao->data; 84 KSPConvergedReason ksp_reason; 85 86 PetscReal resnorm, oldTrust, prered, actred, stepNorm, steplen; 87 PetscBool cgTerminate, needH = PETSC_TRUE, stepAccepted, shift = PETSC_FALSE; 88 PetscInt stepType = BNK_NEWTON; 89 90 PetscFunctionBegin; 91 /* Initialize the preconditioner, KSP solver and trust radius/line search */ 92 tao->reason = TAO_CONTINUE_ITERATING; 93 ierr = TaoBNKInitialize(tao, bnk->init_type, &needH);CHKERRQ(ierr); 94 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 95 96 /* Have not converged; continue with Newton method */ 97 while (tao->reason == TAO_CONTINUE_ITERATING) { 98 tao->niter++; 99 100 if (needH) { 101 /* Take BNCG steps (if enabled) to trade-off Hessian evaluations for more gradient evaluations */ 102 ierr = TaoBNKTakeCGSteps(tao, &cgTerminate);CHKERRQ(ierr); 103 if (cgTerminate) { 104 tao->reason = bnk->bncg->reason; 105 PetscFunctionReturn(0); 106 } 107 /* Compute the hessian and update the BFGS preconditioner at the new iterate */ 108 ierr = TaoBNKComputeHessian(tao);CHKERRQ(ierr); 109 needH = PETSC_FALSE; 110 } 111 112 /* Store current solution before it changes */ 113 oldTrust = tao->trust; 114 bnk->fold = bnk->f; 115 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 116 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 117 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 118 119 /* Enter into trust region loops */ 120 stepAccepted = PETSC_FALSE; 121 while (!stepAccepted && tao->reason == TAO_CONTINUE_ITERATING) { 122 tao->ksp_its=0; 123 124 /* Use the common BNK kernel to compute the Newton step (for inactive variables only) */ 125 ierr = TaoBNKComputeStep(tao, shift, &ksp_reason);CHKERRQ(ierr); 126 127 /* Temporarily accept the step and project it into the bounds */ 128 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 129 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 130 131 /* Check if the projection changed the step direction */ 132 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 133 ierr = VecAXPY(tao->stepdirection, -1.0, bnk->Xold);CHKERRQ(ierr); 134 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 135 if (stepNorm != bnk->dnorm) { 136 /* Projection changed the step, so we have to recompute predicted reduction. 137 However, we deliberately do not change the step norm and the trust radius 138 in order for the safeguard to more closely mimic a piece-wise linesearch 139 along the bounds. */ 140 ierr = TaoBNKRecomputePred(tao, tao->stepdirection, &prered);CHKERRQ(ierr); 141 } else { 142 /* Step did not change, so we can just recover the pre-computed prediction */ 143 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 144 } 145 prered = -prered; 146 147 /* Compute the actual reduction and update the trust radius */ 148 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 149 actred = bnk->fold - bnk->f; 150 oldTrust = tao->trust; 151 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, bnk->update_type, stepType, &stepAccepted);CHKERRQ(ierr); 152 153 if (stepAccepted) { 154 /* Step is good, evaluate the gradient and flip the need-Hessian switch */ 155 steplen = 1.0; 156 needH = PETSC_TRUE; 157 ++bnk->newt; 158 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 159 ierr = TaoBNKEstimateActiveSet(tao, bnk->as_type);CHKERRQ(ierr); 160 ierr = VecCopy(bnk->unprojected_gradient, tao->gradient);CHKERRQ(ierr); 161 ierr = VecISSet(tao->gradient, bnk->active_idx, 0.0);CHKERRQ(ierr); 162 ierr = VecNorm(tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 163 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 164 } else { 165 /* Step is bad, revert old solution and re-solve with new radius*/ 166 steplen = 0.0; 167 needH = PETSC_FALSE; 168 bnk->f = bnk->fold; 169 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 170 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 171 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 172 if (oldTrust == tao->trust) { 173 /* Can't change the radius anymore so just terminate */ 174 tao->reason = TAO_DIVERGED_TR_REDUCTION; 175 } 176 } 177 178 /* Check for termination */ 179 ierr = VecFischer(tao->solution, bnk->unprojected_gradient, tao->XL, tao->XU, bnk->Gwork);CHKERRQ(ierr); 180 ierr = VecNorm(bnk->Gwork, NORM_2, &resnorm);CHKERRQ(ierr); 181 ierr = TaoLogConvergenceHistory(tao, bnk->f, resnorm, 0.0, tao->ksp_its);CHKERRQ(ierr); 182 ierr = TaoMonitor(tao, tao->niter, bnk->f, resnorm, 0.0, steplen);CHKERRQ(ierr); 183 ierr = (*tao->ops->convergencetest)(tao, tao->cnvP);CHKERRQ(ierr); 184 } 185 } 186 PetscFunctionReturn(0); 187 } 188 189 /*------------------------------------------------------------*/ 190 191 PETSC_INTERN PetscErrorCode TaoSetUp_BNTR(Tao tao) 192 { 193 TAO_BNK *bnk = (TAO_BNK *)tao->data; 194 PetscErrorCode ierr; 195 196 PetscFunctionBegin; 197 ierr = TaoSetUp_BNK(tao);CHKERRQ(ierr); 198 if (!bnk->is_nash && !bnk->is_stcg && !bnk->is_gltr) SETERRQ(PETSC_COMM_SELF,1,"Must use a trust-region CG method for KSP (KSPNASH, KSPSTCG, KSPGLTR)"); 199 PetscFunctionReturn(0); 200 } 201 202 /*------------------------------------------------------------*/ 203 204 PETSC_INTERN PetscErrorCode TaoCreate_BNTR(Tao tao) 205 { 206 TAO_BNK *bnk; 207 PetscErrorCode ierr; 208 209 PetscFunctionBegin; 210 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 211 tao->ops->solve=TaoSolve_BNTR; 212 tao->ops->setup=TaoSetUp_BNTR; 213 214 bnk = (TAO_BNK *)tao->data; 215 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 216 PetscFunctionReturn(0); 217 }