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 should be done with a conjugate gradient 10 method, although any method can be used. 11 */ 12 13 static PetscErrorCode TaoSolve_BNTL(Tao tao) 14 { 15 PetscErrorCode ierr; 16 TAO_BNK *bnk = (TAO_BNK *)tao->data; 17 TaoLineSearchConvergedReason ls_reason; 18 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 49 if (stepAccepted) { 50 tao->niter++; 51 tao->ksp_its=0; 52 /* Compute the Hessian */ 53 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 54 /* Update the BFGS preconditioner */ 55 if (BNK_PC_BFGS == bnk->pc_type) { 56 if (BFGS_SCALE_PHESS == bnk->bfgs_scale_type) { 57 /* Obtain diagonal for the bfgs preconditioner */ 58 ierr = MatGetDiagonal(tao->hessian, bnk->Diag);CHKERRQ(ierr); 59 ierr = VecAbs(bnk->Diag);CHKERRQ(ierr); 60 ierr = VecReciprocal(bnk->Diag);CHKERRQ(ierr); 61 ierr = MatLMVMSetScale(bnk->M,bnk->Diag);CHKERRQ(ierr); 62 } 63 /* Update the limited memory preconditioner and get existing # of updates */ 64 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 65 } 66 } 67 68 /* Use the common BNK kernel to compute the raw Newton step */ 69 ierr = TaoBNKComputeStep(tao, PETSC_FALSE, &stepType);CHKERRQ(ierr); 70 71 /* Store current solution before it changes */ 72 oldTrust = tao->trust; 73 bnk->fold = bnk->f; 74 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 75 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 76 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 77 78 /* Temporarily accept the step and project it into the bounds */ 79 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 80 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 81 82 /* Check if the projection changed the step direction */ 83 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 84 ierr = VecAXPBY(tao->stepdirection, -1.0, 1.0, bnk->Xold);CHKERRQ(ierr); 85 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 86 if (stepNorm != bnk->dnorm) { 87 /* Projection changed the step, so we have to adjust trust radius and recompute predicted reduction */ 88 bnk->dnorm = stepNorm; 89 tao->trust = bnk->dnorm; 90 ierr = MatMult(tao->hessian, tao->stepdirection, bnk->Xwork);CHKERRQ(ierr); 91 ierr = VecAYPX(bnk->Xwork, -0.5, tao->gradient);CHKERRQ(ierr); 92 ierr = VecDot(bnk->Xwork, tao->stepdirection, &prered); 93 } else { 94 /* Step did not change, so we can just recover the pre-computed prediction */ 95 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 96 } 97 prered = -prered; 98 99 /* Compute the actual reduction and update the trust radius */ 100 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 101 actred = bnk->fold - bnk->f; 102 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 103 104 if (stepAccepted) { 105 /* Step is good, evaluate the gradient and the hessian */ 106 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 107 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 108 } else { 109 /* Trust-region rejected the step. Revert the solution. */ 110 bnk->f = bnk->fold; 111 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 112 113 /* Now check to make sure the Newton step is a descent direction... */ 114 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 115 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 116 /* Newton step is not descent or direction produced Inf or NaN */ 117 --bnk->newt; 118 if (BNK_PC_BFGS != bnk->pc_type) { 119 /* We don't have the BFGS matrix around and updated 120 Must use gradient direction in this case */ 121 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 122 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 123 ++bnk->grad; 124 stepType = BNK_GRADIENT; 125 } else { 126 /* We have the BFGS matrix, so attempt to use the BFGS direction */ 127 ierr = MatLMVMSolve(bnk->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 128 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 129 130 /* Check for success (descent direction) */ 131 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 132 if ((gdx >= 0) || PetscIsInfOrNanReal(gdx)) { 133 /* BFGS direction is not descent or direction produced not a number 134 We can assert bfgsUpdates > 1 in this case because 135 the first solve produces the scaled gradient direction, 136 which is guaranteed to be descent */ 137 138 /* Use steepest descent direction (scaled) */ 139 if (bnk->f != 0.0) { 140 delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm); 141 } else { 142 delta = 2.0 / (bnk->gnorm*bnk->gnorm); 143 } 144 ierr = MatLMVMSetDelta(bnk->M, delta);CHKERRQ(ierr); 145 ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr); 146 ierr = MatLMVMUpdate(bnk->M, tao->solution, tao->gradient);CHKERRQ(ierr); 147 ierr = MatLMVMSolve(bnk->M, tao->gradient, tao->stepdirection);CHKERRQ(ierr); 148 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 149 150 ++bnk->sgrad; 151 stepType = BNK_SCALED_GRADIENT; 152 } else { 153 ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr); 154 if (1 == bfgsUpdates) { 155 /* The first BFGS direction is always the scaled gradient */ 156 ++bnk->sgrad; 157 stepType = BNK_SCALED_GRADIENT; 158 } else { 159 ++bnk->bfgs; 160 stepType = BNK_BFGS; 161 } 162 } 163 } 164 } 165 166 /* Trigger the line search */ 167 ierr = TaoBNKPerformLineSearch(tao, stepType, &steplen, &ls_reason);CHKERRQ(ierr); 168 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 169 /* Line search failed, revert solution and terminate */ 170 bnk->f = bnk->fold; 171 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 172 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 173 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 174 tao->trust = 0.0; 175 tao->reason = TAO_DIVERGED_LS_FAILURE; 176 } else { 177 /* Line search succeeded so we should update the trust radius based on the LS step length */ 178 updateType = bnk->update_type; 179 bnk->update_type = BNK_UPDATE_STEP; 180 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 181 bnk->update_type = updateType; 182 } 183 } 184 185 /* Check for termination */ 186 ierr = TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 187 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 188 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 189 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,tao->trust);CHKERRQ(ierr); 190 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 191 } 192 PetscFunctionReturn(0); 193 } 194 195 PETSC_EXTERN PetscErrorCode TaoCreate_BNTL(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_BNTL; 203 204 bnk = (TAO_BNK *)tao->data; 205 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 206 bnk->sval = 0.0; /* disable Hessian shifting */ 207 PetscFunctionReturn(0); 208 }