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 KSPConvergedReason ksp_reason; 17 TaoLineSearchConvergedReason ls_reason; 18 19 PetscReal oldTrust, prered, actred, stepNorm, steplen; 20 PetscBool stepAccepted = PETSC_TRUE; 21 PetscInt stepType; 22 23 PetscFunctionBegin; 24 /* Initialize the preconditioner, KSP solver and trust radius/line search */ 25 tao->reason = TAO_CONTINUE_ITERATING; 26 ierr = TaoBNKInitialize(tao);CHKERRQ(ierr); 27 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 28 29 /* Have not converged; continue with Newton method */ 30 while (tao->reason == TAO_CONTINUE_ITERATING) { 31 tao->niter++; 32 tao->ksp_its=0; 33 /* Compute the Hessian */ 34 ierr = TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre);CHKERRQ(ierr); 35 /* Update the BFGS preconditioner */ 36 if (BNK_PC_BFGS == bnk->pc_type) { 37 if (BFGS_SCALE_PHESS == bnk->bfgs_scale_type) { 38 /* Obtain diagonal for the bfgs preconditioner */ 39 ierr = MatGetDiagonal(tao->hessian, bnk->Diag);CHKERRQ(ierr); 40 ierr = VecAbs(bnk->Diag);CHKERRQ(ierr); 41 ierr = VecReciprocal(bnk->Diag);CHKERRQ(ierr); 42 ierr = MatLMVMSetScale(bnk->M,bnk->Diag);CHKERRQ(ierr); 43 } 44 /* Update the limited memory preconditioner and get existing # of updates */ 45 ierr = MatLMVMUpdate(bnk->M, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 46 } 47 48 /* Use the common BNK kernel to compute the Newton step (for inactive variables only) */ 49 ierr = TaoBNKComputeStep(tao, &ksp_reason);CHKERRQ(ierr); 50 51 /* Store current solution before it changes */ 52 oldTrust = tao->trust; 53 bnk->fold = bnk->f; 54 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 55 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 56 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 57 58 /* Temporarily accept the step and project it into the bounds */ 59 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 60 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 61 62 /* Check if the projection changed the step direction */ 63 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 64 ierr = VecAXPY(tao->stepdirection, -1.0, bnk->Xold);CHKERRQ(ierr); 65 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 66 if (stepNorm != bnk->dnorm) { 67 /* Projection changed the step, so we have to recompute predicted reduction. 68 However, we deliberately do not change the step norm and the trust radius 69 in order for the safeguard to more closely mimic a piece-wise linesearch 70 along the bounds. */ 71 ierr = MatMult(bnk->H_inactive, tao->stepdirection, bnk->Xwork);CHKERRQ(ierr); 72 ierr = VecAYPX(bnk->Xwork, -0.5, tao->gradient);CHKERRQ(ierr); 73 ierr = VecDot(bnk->Xwork, tao->stepdirection, &prered); 74 } else { 75 /* Step did not change, so we can just recover the pre-computed prediction */ 76 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 77 } 78 prered = -prered; 79 80 /* Compute the actual reduction and update the trust radius */ 81 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 82 actred = bnk->fold - bnk->f; 83 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, bnk->update_type, stepType, &stepAccepted);CHKERRQ(ierr); 84 85 if (stepAccepted) { 86 /* Step is good, evaluate the gradient and the hessian */ 87 steplen = 1.0; 88 ++bnk->newt; 89 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 90 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 91 } else { 92 /* Trust-region rejected the step. Revert the solution. */ 93 bnk->f = bnk->fold; 94 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 95 96 /* Trigger the line search */ 97 ierr = TaoBNKSafeguardStep(tao, ksp_reason, &stepType);CHKERRQ(ierr); 98 ierr = TaoBNKPerformLineSearch(tao, stepType, &steplen, &ls_reason);CHKERRQ(ierr); 99 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 100 /* Line search failed, revert solution and terminate */ 101 bnk->f = bnk->fold; 102 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 103 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 104 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 105 tao->trust = 0.0; 106 tao->reason = TAO_DIVERGED_LS_FAILURE; 107 } else { 108 /* Line search succeeded so we should update the trust radius based on the LS step length */ 109 tao->trust = oldTrust; 110 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, BNK_UPDATE_STEP, stepType, &stepAccepted);CHKERRQ(ierr); 111 /* count the accepted step types */ 112 switch (stepType) { 113 case BNK_NEWTON: 114 ++bnk->newt; 115 break; 116 case BNK_BFGS: 117 ++bnk->bfgs; 118 break; 119 case BNK_SCALED_GRADIENT: 120 ++bnk->sgrad; 121 break; 122 case BNK_GRADIENT: 123 ++bnk->grad; 124 break; 125 default: 126 break; 127 } 128 } 129 } 130 131 /* Check for termination */ 132 ierr = TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 133 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 134 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 135 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 136 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 137 } 138 PetscFunctionReturn(0); 139 } 140 141 /*------------------------------------------------------------*/ 142 143 PETSC_EXTERN PetscErrorCode TaoCreate_BNTL(Tao tao) 144 { 145 TAO_BNK *bnk; 146 PetscErrorCode ierr; 147 148 PetscFunctionBegin; 149 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 150 tao->ops->solve=TaoSolve_BNTL; 151 152 bnk = (TAO_BNK *)tao->data; 153 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 154 bnk->sval = 0.0; /* disable Hessian shifting */ 155 PetscFunctionReturn(0); 156 }