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