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, stepAccepted = PETSC_TRUE, 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, &stepAccepted);CHKERRQ(ierr); 94 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 95 96 /* Have not converged; continue with Newton method */ 97 if (!stepAccepted) {tao->niter = 1;} 98 while (tao->reason == TAO_CONTINUE_ITERATING) { 99 100 if (stepAccepted) { 101 tao->niter++; 102 tao->ksp_its=0; 103 /* Compute the hessian, update the BFGS preconditioner and estimate the active-set at the new iterate */ 104 ierr = TaoBNKComputeHessian(tao);CHKERRQ(ierr); 105 ierr = TaoBNKEstimateActiveSet(tao);CHKERRQ(ierr); 106 } 107 108 /* Take BNCG steps (if enabled) to trade-off Hessian evaluations for more gradient evaluations */ 109 ierr = TaoBNKTakeCGSteps(tao, &cgTerminate);CHKERRQ(ierr); 110 if (cgTerminate) { 111 tao->reason = bnk->bncg->reason; 112 PetscFunctionReturn(0); 113 } 114 115 /* Use the common BNK kernel to compute the Newton step (for inactive variables only) */ 116 ierr = TaoBNKComputeStep(tao, shift, &ksp_reason);CHKERRQ(ierr); 117 118 /* Store current solution before it changes */ 119 oldTrust = tao->trust; 120 bnk->fold = bnk->f; 121 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 122 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 123 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 124 125 /* Temporarily accept the step and project it into the bounds */ 126 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 127 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 128 129 /* Check if the projection changed the step direction */ 130 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 131 ierr = VecAXPY(tao->stepdirection, -1.0, bnk->Xold);CHKERRQ(ierr); 132 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 133 if (stepNorm != bnk->dnorm) { 134 /* Projection changed the step, so we have to recompute predicted reduction. 135 However, we deliberately do not change the step norm and the trust radius 136 in order for the safeguard to more closely mimic a piece-wise linesearch 137 along the bounds. */ 138 ierr = TaoBNKRecomputePred(tao, tao->stepdirection, &prered);CHKERRQ(ierr); 139 } else { 140 /* Step did not change, so we can just recover the pre-computed prediction */ 141 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 142 } 143 prered = -prered; 144 145 /* Compute the actual reduction and update the trust radius */ 146 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 147 actred = bnk->fold - bnk->f; 148 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, bnk->update_type, stepType, &stepAccepted);CHKERRQ(ierr); 149 150 if (stepAccepted) { 151 /* Step is good, evaluate the gradient and the hessian */ 152 steplen = 1.0; 153 ++bnk->newt; 154 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 155 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 156 ierr = VecNorm(tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 157 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 158 } else { 159 /* Step is bad, revert old solution and re-solve with new radius*/ 160 steplen = 0.0; 161 bnk->f = bnk->fold; 162 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 163 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 164 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 165 if (oldTrust == tao->trust) { 166 /* Can't change the radius anymore so just terminate */ 167 tao->reason = TAO_DIVERGED_TR_REDUCTION; 168 } 169 } 170 171 /* Check for termination */ 172 ierr = VecFischer(tao->solution, bnk->unprojected_gradient, tao->XL, tao->XU, bnk->Gwork);CHKERRQ(ierr); 173 ierr = VecNorm(bnk->Gwork, NORM_2, &resnorm);CHKERRQ(ierr); 174 ierr = TaoLogConvergenceHistory(tao, bnk->f, resnorm, 0.0, tao->ksp_its);CHKERRQ(ierr); 175 ierr = TaoMonitor(tao, tao->niter, bnk->f, resnorm, 0.0, steplen);CHKERRQ(ierr); 176 ierr = (*tao->ops->convergencetest)(tao, tao->cnvP);CHKERRQ(ierr); 177 } 178 PetscFunctionReturn(0); 179 } 180 181 /*------------------------------------------------------------*/ 182 183 PETSC_INTERN PetscErrorCode TaoSetUp_BNTR(Tao tao) 184 { 185 TAO_BNK *bnk = (TAO_BNK *)tao->data; 186 PetscErrorCode ierr; 187 188 PetscFunctionBegin; 189 ierr = TaoSetUp_BNK(tao);CHKERRQ(ierr); 190 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)"); 191 PetscFunctionReturn(0); 192 } 193 194 /*------------------------------------------------------------*/ 195 196 PETSC_INTERN PetscErrorCode TaoCreate_BNTR(Tao tao) 197 { 198 TAO_BNK *bnk; 199 PetscErrorCode ierr; 200 201 PetscFunctionBegin; 202 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 203 tao->ops->solve=TaoSolve_BNTR; 204 tao->ops->setup=TaoSetUp_BNTR; 205 206 bnk = (TAO_BNK *)tao->data; 207 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 208 PetscFunctionReturn(0); 209 }