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 tao->niter++; 48 tao->ksp_its=0; 49 /* Compute the Hessian */ 50 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 51 /* Update the BFGS preconditioner */ 52 if (BNK_PC_BFGS == bnk->pc_type) { 53 if (BFGS_SCALE_PHESS == bnk->bfgs_scale_type) { 54 /* Obtain diagonal for the bfgs preconditioner */ 55 ierr = MatGetDiagonal(tao->hessian, bnk->Diag);CHKERRQ(ierr); 56 ierr = VecAbs(bnk->Diag);CHKERRQ(ierr); 57 ierr = VecReciprocal(bnk->Diag);CHKERRQ(ierr); 58 ierr = MatLMVMSetScale(bnk->M,bnk->Diag);CHKERRQ(ierr); 59 } 60 /* Update the limited memory preconditioner and get existing # of updates */ 61 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 62 } 63 64 /* Use the common BNK kernel to compute the Newton step (for inactive variables only) */ 65 ierr = TaoBNKComputeStep(tao, PETSC_FALSE, &stepType);CHKERRQ(ierr); 66 67 /* Store current solution before it changes */ 68 oldTrust = tao->trust; 69 bnk->fold = bnk->f; 70 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 71 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 72 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 73 74 /* Temporarily accept the step and project it into the bounds */ 75 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 76 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 77 78 /* Check if the projection changed the step direction */ 79 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 80 ierr = VecAXPY(tao->stepdirection, -1.0, bnk->Xold);CHKERRQ(ierr); 81 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 82 if (stepNorm != bnk->dnorm) { 83 /* Projection changed the step, so we have to recompute predicted reduction. 84 However, we deliberately do not change the step norm and the trust radius 85 in order for the safeguard to more closely mimic a piece-wise linesearch 86 along the bounds. */ 87 ierr = MatMult(tao->hessian, tao->stepdirection, bnk->Xwork);CHKERRQ(ierr); 88 ierr = VecAYPX(bnk->Xwork, -0.5, tao->gradient);CHKERRQ(ierr); 89 ierr = VecDot(bnk->Xwork, tao->stepdirection, &prered); 90 } else { 91 /* Step did not change, so we can just recover the pre-computed prediction */ 92 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 93 } 94 prered = -prered; 95 96 /* Compute the actual reduction and update the trust radius */ 97 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 98 actred = bnk->fold - bnk->f; 99 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 100 101 if (stepAccepted) { 102 /* Step is good, evaluate the gradient and the hessian */ 103 steplen = 1.0; 104 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 105 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 106 } else { 107 /* Trust-region rejected the step. Revert the solution. */ 108 bnk->f = bnk->fold; 109 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 110 111 /* Now check to make sure the Newton step is a descent direction... */ 112 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 113 if ((gdx >= 0.0) || PetscIsInfOrNanReal(gdx)) { 114 /* Newton step is not descent or direction produced Inf or NaN */ 115 --bnk->newt; 116 if (BNK_PC_BFGS != bnk->pc_type) { 117 /* We don't have the BFGS matrix around and updated 118 Must use gradient direction in this case */ 119 ierr = VecCopy(tao->gradient, tao->stepdirection);CHKERRQ(ierr); 120 ++bnk->grad; 121 stepType = BNK_GRADIENT; 122 } else { 123 /* We have the BFGS matrix, so attempt to use the BFGS direction */ 124 ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 125 126 /* Check for success (descent direction) 127 NOTE: Negative gdx here means not a descent direction because 128 the fall-back step is missing a negative sign. */ 129 ierr = VecDot(tao->stepdirection, tao->gradient, &gdx);CHKERRQ(ierr); 130 if ((gdx <= 0) || PetscIsInfOrNanReal(gdx)) { 131 /* BFGS direction is not descent or direction produced not a number 132 We can assert bfgsUpdates > 1 in this case because 133 the first solve produces the scaled gradient direction, 134 which is guaranteed to be descent */ 135 136 /* Use steepest descent direction (scaled) */ 137 if (bnk->f != 0.0) { 138 delta = 2.0 * PetscAbsScalar(bnk->f) / (bnk->gnorm*bnk->gnorm); 139 } else { 140 delta = 2.0 / (bnk->gnorm*bnk->gnorm); 141 } 142 ierr = MatLMVMSetDelta(bnk->M, delta);CHKERRQ(ierr); 143 ierr = MatLMVMReset(bnk->M);CHKERRQ(ierr); 144 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 145 ierr = MatLMVMSolve(bnk->M, bnk->unprojected_gradient, tao->stepdirection);CHKERRQ(ierr); 146 147 ++bnk->sgrad; 148 stepType = BNK_SCALED_GRADIENT; 149 } else { 150 ierr = MatLMVMGetUpdates(bnk->M, &bfgsUpdates);CHKERRQ(ierr); 151 if (1 == bfgsUpdates) { 152 /* The first BFGS direction is always the scaled gradient */ 153 ++bnk->sgrad; 154 stepType = BNK_SCALED_GRADIENT; 155 } else { 156 ++bnk->bfgs; 157 stepType = BNK_BFGS; 158 } 159 } 160 } 161 } 162 /* Make sure the safeguarded fall-back step is zero for actively bounded variables */ 163 ierr = VecScale(tao->stepdirection, -1.0);CHKERRQ(ierr); 164 ierr = TaoBNKBoundStep(tao, tao->stepdirection);CHKERRQ(ierr); 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 tao->trust = oldTrust; 181 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, stepType, &stepAccepted);CHKERRQ(ierr); 182 bnk->update_type = updateType; 183 } 184 } 185 186 /* Check for termination */ 187 ierr = TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 188 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 189 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 190 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 191 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 192 } 193 PetscFunctionReturn(0); 194 } 195 196 /*------------------------------------------------------------*/ 197 198 PETSC_EXTERN PetscErrorCode TaoCreate_BNTL(Tao tao) 199 { 200 TAO_BNK *bnk; 201 PetscErrorCode ierr; 202 203 PetscFunctionBegin; 204 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 205 tao->ops->solve=TaoSolve_BNTL; 206 207 bnk = (TAO_BNK *)tao->data; 208 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 209 bnk->sval = 0.0; /* disable Hessian shifting */ 210 PetscFunctionReturn(0); 211 }