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 PetscReal oldTrust, prered, actred, stepNorm, gdx, delta, steplen; 19 PetscBool stepAccepted = PETSC_TRUE; 20 PetscInt stepType, bfgsUpdates, updateType; 21 22 PetscFunctionBegin; 23 /* Project the current point onto the feasible set */ 24 ierr = TaoComputeVariableBounds(tao);CHKERRQ(ierr); 25 ierr = TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);CHKERRQ(ierr); 26 27 /* Project the initial point onto the feasible region */ 28 ierr = VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);CHKERRQ(ierr); 29 30 /* Check convergence criteria */ 31 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &bnk->f, bnk->unprojected_gradient);CHKERRQ(ierr); 32 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 33 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&bnk->gnorm);CHKERRQ(ierr); 34 if (PetscIsInfOrNanReal(bnk->f) || PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 35 36 tao->reason = TAO_CONTINUE_ITERATING; 37 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 38 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,tao->trust);CHKERRQ(ierr); 39 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 40 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 41 42 /* Initialize the preconditioner and trust radius */ 43 ierr = TaoBNKInitialize(tao);CHKERRQ(ierr); 44 45 /* Have not converged; continue with Newton method */ 46 while (tao->reason == TAO_CONTINUE_ITERATING) { 47 48 if (stepAccepted) { 49 tao->niter++; 50 tao->ksp_its=0; 51 /* Compute the Hessian */ 52 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 53 /* Update the BFGS preconditioner */ 54 if (BNK_PC_BFGS == bnk->pc_type) { 55 if (BFGS_SCALE_PHESS == bnk->bfgs_scale_type) { 56 /* Obtain diagonal for the bfgs preconditioner */ 57 ierr = MatGetDiagonal(tao->hessian, bnk->Diag);CHKERRQ(ierr); 58 ierr = VecAbs(bnk->Diag);CHKERRQ(ierr); 59 ierr = VecReciprocal(bnk->Diag);CHKERRQ(ierr); 60 ierr = MatLMVMSetScale(bnk->M,bnk->Diag);CHKERRQ(ierr); 61 } 62 /* Update the limited memory preconditioner and get existing # of updates */ 63 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 64 } 65 } 66 67 /* Use the common BNK kernel to compute the raw Newton step */ 68 ierr = TaoBNKComputeStep(tao, PETSC_FALSE, &stepType);CHKERRQ(ierr); 69 70 /* Store current solution before it changes */ 71 oldTrust = tao->trust; 72 bnk->fold = bnk->f; 73 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 74 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 75 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 76 77 /* Temporarily accept the step and project it into the bounds */ 78 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 79 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 80 81 /* Check if the projection changed the step direction */ 82 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 83 ierr = VecAXPBY(tao->stepdirection, -1.0, 1.0, bnk->Xold);CHKERRQ(ierr); 84 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 85 if (stepNorm != bnk->dnorm) { 86 /* Projection changed the step, so we have to adjust trust radius and recompute predicted reduction */ 87 bnk->dnorm = stepNorm; 88 tao->trust = bnk->dnorm; 89 ierr = MatMult(tao->hessian, tao->stepdirection, bnk->Xwork);CHKERRQ(ierr); 90 ierr = VecAYPX(bnk->Xwork, -0.5, tao->gradient);CHKERRQ(ierr); 91 ierr = VecDot(bnk->Xwork, tao->stepdirection, &prered); 92 } else { 93 /* Step did not change, so we can just recover the pre-computed prediction */ 94 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 95 } 96 prered = -prered; 97 98 /* Compute the actual reduction and update the trust radius */ 99 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 100 actred = bnk->fold - bnk->f; 101 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 102 103 if (stepAccepted) { 104 /* Step is good, evaluate the gradient and the hessian */ 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 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 122 ++bnk->grad; 123 stepType = BNK_GRADIENT; 124 } else { 125 /* We have the BFGS matrix, so attempt to use the BFGS direction */ 126 ierr = MatLMVMSolve(bnk->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 127 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 128 129 /* Check for success (descent direction) */ 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, tao->gradient);CHKERRQ(ierr); 146 ierr = MatLMVMSolve(bnk->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 147 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 148 149 ++bnk->sgrad; 150 stepType = BNK_SCALED_GRADIENT; 151 } else { 152 ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr); 153 if (1 == bfgsUpdates) { 154 /* The first BFGS direction is always the scaled gradient */ 155 ++bnk->sgrad; 156 stepType = BNK_SCALED_GRADIENT; 157 } else { 158 ++bnk->bfgs; 159 stepType = BNK_BFGS; 160 } 161 } 162 } 163 } 164 165 /* Trigger the line search */ 166 ierr = TaoBNKPerformLineSearch(tao, stepType, &steplen, &ls_reason);CHKERRQ(ierr); 167 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 168 /* Line search failed, revert solution and terminate */ 169 bnk->f = bnk->fold; 170 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 171 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 172 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 173 tao->trust = 0.0; 174 tao->reason = TAO_DIVERGED_LS_FAILURE; 175 } else { 176 /* Line search succeeded so we should update the trust radius based on the LS step length */ 177 updateType = bnk->update_type; 178 bnk->update_type = BNK_UPDATE_STEP; 179 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 180 bnk->update_type = updateType; 181 } 182 } 183 184 /* Check for termination */ 185 ierr = TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 186 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 187 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 188 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,tao->trust);CHKERRQ(ierr); 189 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 190 } 191 PetscFunctionReturn(0); 192 } 193 194 /*------------------------------------------------------------*/ 195 196 PETSC_EXTERN PetscErrorCode TaoCreate_BNTL(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_BNTL; 204 205 bnk = (TAO_BNK *)tao->data; 206 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 207 bnk->sval = 0.0; /* disable Hessian shifting */ 208 PetscFunctionReturn(0); 209 }