1 #include <../src/tao/bound/impls/bnk/bnk.h> 2 #include <petscksp.h> 3 4 /* 5 Implements Newton's Method with a line search approach for 6 solving bound constrained minimization problems. 7 8 ------------------------------------------------------------ 9 10 x_0 = VecMedian(x_0) 11 f_0, g_0 = TaoComputeObjectiveAndGradient(x_0) 12 pg_0 = project(g_0) 13 check convergence at pg_0 14 needH = TaoBNKInitialize(default:BNK_INIT_DIRECTION) 15 niter = 0 16 step_accepted = true 17 18 while niter < max_it 19 if needH 20 If max_cg_steps > 0 21 x_k, g_k, pg_k = TaoSolve(BNCG) 22 end 23 24 H_k = TaoComputeHessian(x_k) 25 if pc_type == BNK_PC_BFGS 26 add correction to BFGS approx 27 if scale_type == BNK_SCALE_AHESS 28 D = VecMedian(1e-6, abs(diag(H_k)), 1e6) 29 scale BFGS with VecReciprocal(D) 30 end 31 end 32 needH = False 33 end 34 35 if pc_type = BNK_PC_BFGS 36 B_k = BFGS 37 else 38 B_k = VecMedian(1e-6, abs(diag(H_k)), 1e6) 39 B_k = VecReciprocal(B_k) 40 end 41 w = x_k - VecMedian(x_k - 0.001*B_k*g_k) 42 eps = min(eps, norm2(w)) 43 determine the active and inactive index sets such that 44 L = {i : (x_k)_i <= l_i + eps && (g_k)_i > 0} 45 U = {i : (x_k)_i >= u_i - eps && (g_k)_i < 0} 46 F = {i : l_i = (x_k)_i = u_i} 47 A = {L + U + F} 48 IA = {i : i not in A} 49 50 generate the reduced system Hr_k dr_k = -gr_k for variables in IA 51 if p > 0 52 Hr_k += p* 53 end 54 if pc_type == BNK_PC_BFGS && scale_type == BNK_SCALE_PHESS 55 D = VecMedian(1e-6, abs(diag(Hr_k)), 1e6) 56 scale BFGS with VecReciprocal(D) 57 end 58 solve Hr_k dr_k = -gr_k 59 set d_k to (l - x) for variables in L, (u - x) for variables in U, and 0 for variables in F 60 61 if dot(d_k, pg_k)) >= 0 || norm(d_k) == NaN || norm(d_k) == Inf 62 dr_k = -BFGS*gr_k for variables in I 63 if dot(d_k, pg_k)) >= 0 || norm(d_k) == NaN || norm(d_k) == Inf 64 reset the BFGS preconditioner 65 calculate scale delta and apply it to BFGS 66 dr_k = -BFGS*gr_k for variables in I 67 if dot(d_k, pg_k)) >= 0 || norm(d_k) == NaN || norm(d_k) == Inf 68 dr_k = -gr_k for variables in I 69 end 70 end 71 end 72 73 x_{k+1}, f_{k+1}, g_{k+1}, ls_failed = TaoBNKPerformLineSearch() 74 if ls_failed 75 f_{k+1} = f_k 76 x_{k+1} = x_k 77 g_{k+1} = g_k 78 pg_{k+1} = pg_k 79 terminate 80 else 81 pg_{k+1} = project(g_{k+1}) 82 count the accepted step type (Newton, BFGS, scaled grad or grad) 83 end 84 85 niter += 1 86 check convergence at pg_{k+1} 87 end 88 */ 89 90 PetscErrorCode TaoSolve_BNLS(Tao tao) { 91 TAO_BNK *bnk = (TAO_BNK *)tao->data; 92 KSPConvergedReason ksp_reason; 93 TaoLineSearchConvergedReason ls_reason; 94 PetscReal steplen = 1.0, resnorm; 95 PetscBool cgTerminate, needH = PETSC_TRUE, stepAccepted, shift = PETSC_TRUE; 96 PetscInt stepType; 97 98 PetscFunctionBegin; 99 /* Initialize the preconditioner, KSP solver and trust radius/line search */ 100 tao->reason = TAO_CONTINUE_ITERATING; 101 PetscCall(TaoBNKInitialize(tao, bnk->init_type, &needH)); 102 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 103 104 /* Have not converged; continue with Newton method */ 105 while (tao->reason == TAO_CONTINUE_ITERATING) { 106 /* Call general purpose update function */ 107 if (tao->ops->update) { 108 PetscUseTypeMethod(tao, update, tao->niter, tao->user_update); 109 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, bnk->unprojected_gradient)); 110 } 111 112 if (needH && bnk->inactive_idx) { 113 /* Take BNCG steps (if enabled) to trade-off Hessian evaluations for more gradient evaluations */ 114 PetscCall(TaoBNKTakeCGSteps(tao, &cgTerminate)); 115 if (cgTerminate) { 116 tao->reason = bnk->bncg->reason; 117 PetscFunctionReturn(0); 118 } 119 /* Compute the hessian and update the BFGS preconditioner at the new iterate */ 120 PetscCall((*bnk->computehessian)(tao)); 121 needH = PETSC_FALSE; 122 } 123 124 /* Use the common BNK kernel to compute the safeguarded Newton step (for inactive variables only) */ 125 PetscCall((*bnk->computestep)(tao, shift, &ksp_reason, &stepType)); 126 PetscCall(TaoBNKSafeguardStep(tao, ksp_reason, &stepType)); 127 128 /* Store current solution before it changes */ 129 bnk->fold = bnk->f; 130 PetscCall(VecCopy(tao->solution, bnk->Xold)); 131 PetscCall(VecCopy(tao->gradient, bnk->Gold)); 132 PetscCall(VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old)); 133 134 /* Trigger the line search */ 135 PetscCall(TaoBNKPerformLineSearch(tao, &stepType, &steplen, &ls_reason)); 136 137 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 138 /* Failed to find an improving point */ 139 needH = PETSC_FALSE; 140 bnk->f = bnk->fold; 141 PetscCall(VecCopy(bnk->Xold, tao->solution)); 142 PetscCall(VecCopy(bnk->Gold, tao->gradient)); 143 PetscCall(VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient)); 144 steplen = 0.0; 145 tao->reason = TAO_DIVERGED_LS_FAILURE; 146 } else { 147 /* new iterate so we need to recompute the Hessian */ 148 needH = PETSC_TRUE; 149 /* compute the projected gradient */ 150 PetscCall(TaoBNKEstimateActiveSet(tao, bnk->as_type)); 151 PetscCall(VecCopy(bnk->unprojected_gradient, tao->gradient)); 152 PetscCall(VecISSet(tao->gradient, bnk->active_idx, 0.0)); 153 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm)); 154 /* update the trust radius based on the step length */ 155 PetscCall(TaoBNKUpdateTrustRadius(tao, 0.0, 0.0, BNK_UPDATE_STEP, stepType, &stepAccepted)); 156 /* count the accepted step type */ 157 PetscCall(TaoBNKAddStepCounts(tao, stepType)); 158 /* active BNCG recycling for next iteration */ 159 PetscCall(TaoSetRecycleHistory(bnk->bncg, PETSC_TRUE)); 160 } 161 162 /* Check for termination */ 163 PetscCall(VecFischer(tao->solution, bnk->unprojected_gradient, tao->XL, tao->XU, bnk->W)); 164 PetscCall(VecNorm(bnk->W, NORM_2, &resnorm)); 165 PetscCheck(!PetscIsInfOrNanReal(resnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 166 ++tao->niter; 167 PetscCall(TaoLogConvergenceHistory(tao, bnk->f, resnorm, 0.0, tao->ksp_its)); 168 PetscCall(TaoMonitor(tao, tao->niter, bnk->f, resnorm, 0.0, steplen)); 169 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 170 } 171 PetscFunctionReturn(0); 172 } 173 174 /*------------------------------------------------------------*/ 175 /*MC 176 TAOBNLS - Bounded Newton Line Search for nonlinear minimization with bound constraints. 177 178 Options Database Keys: 179 + -tao_bnk_max_cg_its - maximum number of bounded conjugate-gradient iterations taken in each Newton loop 180 . -tao_bnk_init_type - trust radius initialization method ("constant", "direction", "interpolation") 181 . -tao_bnk_update_type - trust radius update method ("step", "direction", "interpolation") 182 - -tao_bnk_as_type - active-set estimation method ("none", "bertsekas") 183 184 Level: beginner 185 M*/ 186 PETSC_EXTERN PetscErrorCode TaoCreate_BNLS(Tao tao) { 187 TAO_BNK *bnk; 188 189 PetscFunctionBegin; 190 PetscCall(TaoCreate_BNK(tao)); 191 tao->ops->solve = TaoSolve_BNLS; 192 193 bnk = (TAO_BNK *)tao->data; 194 bnk->init_type = BNK_INIT_DIRECTION; 195 bnk->update_type = BNK_UPDATE_STEP; /* trust region updates based on line search step length */ 196 PetscFunctionReturn(0); 197 } 198