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. This version includes a 7 line search fall-back in the event of a trust region failure. 8 9 The linear system solve has to be done with a conjugate gradient method. 10 */ 11 12 static PetscErrorCode TaoSolve_BNTL(Tao tao) 13 { 14 PetscErrorCode ierr; 15 TAO_BNK *bnk = (TAO_BNK *)tao->data; 16 TaoLineSearchConvergedReason ls_reason; 17 18 Vec active_step; 19 PetscReal oldTrust, prered, actred, stepNorm, gdx, delta, steplen; 20 PetscBool stepAccepted = PETSC_TRUE; 21 PetscInt stepType, bfgsUpdates, updateType; 22 23 PetscFunctionBegin; 24 /* Project the current point onto the feasible set */ 25 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 26 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 27 28 /* Project the initial point onto the feasible region */ 29 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 30 31 /* Check convergence criteria */ 32 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, bnk->unprojected_gradient);CHKERRQ(ierr); 33 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 34 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 35 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 36 37 tao->reason = TAO_CONTINUE_ITERATING; 38 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 39 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,tao->trust);CHKERRQ(ierr); 40 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 41 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 42 43 /* Initialize the preconditioner and trust radius */ 44 ierr = TaoBNKInitialize(tao);CHKERRQ(ierr); 45 46 /* Have not converged; continue with Newton method */ 47 while (tao->reason == TAO_CONTINUE_ITERATING) { 48 tao->niter++; 49 tao->ksp_its=0; 50 /* Compute the Hessian */ 51 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 52 /* Update the BFGS preconditioner */ 53 if (BNK_PC_BFGS == bnk->pc_type) { 54 if (BFGS_SCALE_PHESS == bnk->bfgs_scale_type) { 55 /* Obtain diagonal for the bfgs preconditioner */ 56 ierr = MatGetDiagonal(tao->hessian, bnk->Diag);CHKERRQ(ierr); 57 ierr = VecAbs(bnk->Diag);CHKERRQ(ierr); 58 ierr = VecReciprocal(bnk->Diag);CHKERRQ(ierr); 59 ierr = MatLMVMSetScale(bnk->M,bnk->Diag);CHKERRQ(ierr); 60 } 61 /* Update the limited memory preconditioner and get existing # of updates */ 62 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 63 } 64 65 /* Use the common BNK kernel to compute the Newton step (for inactive variables only) */ 66 ierr = TaoBNKComputeStep(tao, PETSC_FALSE, &stepType);CHKERRQ(ierr); 67 68 /* Store current solution before it changes */ 69 oldTrust = tao->trust; 70 bnk->fold = bnk->f; 71 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 72 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 73 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 74 75 /* Temporarily accept the step and project it into the bounds */ 76 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 77 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 78 79 /* Check if the projection changed the step direction */ 80 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 81 ierr = VecAXPY(tao->stepdirection, -1.0, bnk->Xold);CHKERRQ(ierr); 82 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 83 if (stepNorm != bnk->dnorm) { 84 /* Projection changed the step, so we have to recompute predicted reduction. 85 However, we deliberately do not change the step norm and the trust radius 86 in order for the safeguard to more closely mimic a piece-wise linesearch 87 along the bounds. */ 88 ierr = MatMult(tao->hessian, tao->stepdirection, bnk->Xwork);CHKERRQ(ierr); 89 ierr = VecAYPX(bnk->Xwork, -0.5, tao->gradient);CHKERRQ(ierr); 90 ierr = VecDot(bnk->Xwork, tao->stepdirection, &prered); 91 } else { 92 /* Step did not change, so we can just recover the pre-computed prediction */ 93 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 94 } 95 prered = -prered; 96 97 /* Compute the actual reduction and update the trust radius */ 98 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 99 actred = bnk->fold - bnk->f; 100 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 101 102 if (stepAccepted) { 103 /* Step is good, evaluate the gradient and the hessian */ 104 steplen = 1.0; 105 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 106 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 107 } else { 108 /* Trust-region rejected the step. Revert the solution. */ 109 bnk->f = bnk->fold; 110 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 111 112 /* Now check to make sure the Newton step is a descent direction... */ 113 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 114 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 115 /* Newton step is not descent or direction produced Inf or NaN */ 116 --bnk->newt; 117 if (BNK_PC_BFGS != bnk->pc_type) { 118 /* We don't have the BFGS matrix around and updated 119 Must use gradient direction in this case */ 120 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 121 ++bnk->grad; 122 stepType = BNK_GRADIENT; 123 } else { 124 /* We have the BFGS matrix, so attempt to use the BFGS direction */ 125 ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 126 127 /* Check for success (descent direction) 128 NOTE: Negative gdx here means not a descent direction because 129 the fall-back step is missing a negative sign. */ 130 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 131 if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) { 132 /* BFGS direction is not descent or direction produced not a number 133 We can assert bfgsUpdates > 1 in this case because 134 the first solve produces the scaled gradient direction, 135 which is guaranteed to be descent */ 136 137 /* Use steepest descent direction (scaled) */ 138 if (bnk->f != 0.0) { 139 delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm); 140 } else { 141 delta = 2.0 / (bnk->gnorm*bnk->gnorm); 142 } 143 ierr = MatLMVMSetDelta(bnk->M, delta);CHKERRQ(ierr); 144 ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr); 145 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 146 ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 147 148 ++bnk->sgrad; 149 stepType = BNK_SCALED_GRADIENT; 150 } else { 151 ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr); 152 if (1 == bfgsUpdates) { 153 /* The first BFGS direction is always the scaled gradient */ 154 ++bnk->sgrad; 155 stepType = BNK_SCALED_GRADIENT; 156 } else { 157 ++bnk->bfgs; 158 stepType = BNK_BFGS; 159 } 160 } 161 } 162 } 163 /* Make sure that the step is zero for actively bounded variables */ 164 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 165 ierr = VecGetSubVector(tao->stepdirection, bnk->active_idx, &active_step);CHKERRQ(ierr); 166 ierr = VecSet(active_step, 0.0);CHKERRQ(ierr); 167 ierr = VecRestoreSubVector(tao->stepdirection, bnk->active_idx, &active_step);CHKERRQ(ierr); 168 169 /* Trigger the line search */ 170 ierr = TaoBNKPerformLineSearch(tao, stepType, &steplen, &ls_reason);CHKERRQ(ierr); 171 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 172 /* Line search failed, revert solution and terminate */ 173 bnk->f = bnk->fold; 174 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 175 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 176 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 177 tao->trust = 0.0; 178 tao->reason = TAO_DIVERGED_LS_FAILURE; 179 } else { 180 /* Line search succeeded so we should update the trust radius based on the LS step length */ 181 updateType = bnk->update_type; 182 bnk->update_type = BNK_UPDATE_STEP; 183 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 184 bnk->update_type = updateType; 185 } 186 } 187 188 /* Check for termination */ 189 ierr = TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 190 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 191 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 192 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 193 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 194 } 195 PetscFunctionReturn(0); 196 } 197 198 /*------------------------------------------------------------*/ 199 200 PETSC_EXTERN PetscErrorCode TaoCreate_BNTL(Tao tao) 201 { 202 TAO_BNK *bnk; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 207 tao->ops->solve=TaoSolve_BNTL; 208 209 bnk = (TAO_BNK *)tao->data; 210 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 211 bnk->sval = 0.0; /* disable Hessian shifting */ 212 PetscFunctionReturn(0); 213 }