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, shift = PETSC_FALSE; 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, bnk->init_type);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 34 /* Compute the hessian and update the BFGS preconditioner at the new iterate*/ 35 ierr = TaoBNKComputeHessian(tao);CHKERRQ(ierr); 36 37 /* Use the common BNK kernel to compute the Newton step (for inactive variables only) */ 38 ierr = TaoBNKComputeStep(tao, shift, &ksp_reason);CHKERRQ(ierr); 39 40 /* Store current solution before it changes */ 41 oldTrust = tao->trust; 42 bnk->fold = bnk->f; 43 ierr = VecCopy(tao->solution, bnk->Xold);CHKERRQ(ierr); 44 ierr = VecCopy(tao->gradient, bnk->Gold);CHKERRQ(ierr); 45 ierr = VecCopy(bnk->unprojected_gradient, bnk->unprojected_gradient_old);CHKERRQ(ierr); 46 47 /* Temporarily accept the step and project it into the bounds */ 48 ierr = VecAXPY(tao->solution, 1.0, tao->stepdirection);CHKERRQ(ierr); 49 ierr = VecMedian(tao->XL, tao->solution, tao->XU, tao->solution);CHKERRQ(ierr); 50 51 /* Check if the projection changed the step direction */ 52 ierr = VecCopy(tao->solution, tao->stepdirection);CHKERRQ(ierr); 53 ierr = VecAXPY(tao->stepdirection, -1.0, bnk->Xold);CHKERRQ(ierr); 54 ierr = VecNorm(tao->stepdirection, NORM_2, &stepNorm);CHKERRQ(ierr); 55 if (stepNorm != bnk->dnorm) { 56 /* Projection changed the step, so we have to recompute predicted reduction. 57 However, we deliberately do not change the step norm and the trust radius 58 in order for the safeguard to more closely mimic a piece-wise linesearch 59 along the bounds. */ 60 ierr = MatMult(bnk->H_inactive, tao->stepdirection, bnk->Xwork);CHKERRQ(ierr); 61 ierr = VecAYPX(bnk->Xwork, -0.5, tao->gradient);CHKERRQ(ierr); 62 ierr = VecDot(bnk->Xwork, tao->stepdirection, &prered); 63 } else { 64 /* Step did not change, so we can just recover the pre-computed prediction */ 65 ierr = KSPCGGetObjFcn(tao->ksp, &prered);CHKERRQ(ierr); 66 } 67 prered = -prered; 68 69 /* Compute the actual reduction and update the trust radius */ 70 ierr = TaoComputeObjective(tao, tao->solution, &bnk->f);CHKERRQ(ierr); 71 actred = bnk->fold - bnk->f; 72 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, bnk->update_type, stepType, &stepAccepted);CHKERRQ(ierr); 73 74 if (stepAccepted) { 75 /* Step is good, evaluate the gradient and the hessian */ 76 steplen = 1.0; 77 ++bnk->newt; 78 ierr = TaoComputeGradient(tao, tao->solution, bnk->unprojected_gradient);CHKERRQ(ierr); 79 ierr = VecBoundGradientProjection(bnk->unprojected_gradient,tao->solution,tao->XL,tao->XU,tao->gradient);CHKERRQ(ierr); 80 } else { 81 /* Trust-region rejected the step. Revert the solution. */ 82 bnk->f = bnk->fold; 83 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 84 85 /* Trigger the line search */ 86 ierr = TaoBNKSafeguardStep(tao, ksp_reason, &stepType);CHKERRQ(ierr); 87 ierr = TaoBNKPerformLineSearch(tao, stepType, &steplen, &ls_reason);CHKERRQ(ierr); 88 if (ls_reason != TAOLINESEARCH_SUCCESS && ls_reason != TAOLINESEARCH_SUCCESS_USER) { 89 /* Line search failed, revert solution and terminate */ 90 bnk->f = bnk->fold; 91 ierr = VecCopy(bnk->Xold, tao->solution);CHKERRQ(ierr); 92 ierr = VecCopy(bnk->Gold, tao->gradient);CHKERRQ(ierr); 93 ierr = VecCopy(bnk->unprojected_gradient_old, bnk->unprojected_gradient);CHKERRQ(ierr); 94 tao->trust = 0.0; 95 tao->reason = TAO_DIVERGED_LS_FAILURE; 96 } else { 97 /* Line search succeeded so we should update the trust radius based on the LS step length */ 98 tao->trust = oldTrust; 99 ierr = TaoBNKUpdateTrustRadius(tao, prered, actred, BNK_UPDATE_STEP, stepType, &stepAccepted);CHKERRQ(ierr); 100 /* count the accepted step type */ 101 ierr = TaoBNKAddStepCounts(tao, stepType);CHKERRQ(ierr); 102 } 103 } 104 105 /* Check for termination */ 106 ierr = TaoGradientNorm(tao, tao->gradient, NORM_2, &bnk->gnorm);CHKERRQ(ierr); 107 if (PetscIsInfOrNanReal(bnk->gnorm)) SETERRQ(PETSC_COMM_SELF,1,"User provided compute function generated Not-a-Number"); 108 ierr = TaoLogConvergenceHistory(tao,bnk->f,bnk->gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 109 ierr = TaoMonitor(tao,tao->niter,bnk->f,bnk->gnorm,0.0,steplen);CHKERRQ(ierr); 110 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 111 } 112 PetscFunctionReturn(0); 113 } 114 115 /*------------------------------------------------------------*/ 116 117 PETSC_EXTERN PetscErrorCode TaoCreate_BNTL(Tao tao) 118 { 119 TAO_BNK *bnk; 120 PetscErrorCode ierr; 121 122 PetscFunctionBegin; 123 ierr = TaoCreate_BNK(tao);CHKERRQ(ierr); 124 tao->ops->solve=TaoSolve_BNTL; 125 126 bnk = (TAO_BNK *)tao->data; 127 bnk->update_type = BNK_UPDATE_REDUCTION; /* trust region updates based on predicted/actual reduction */ 128 bnk->sval = 0.0; /* disable Hessian shifting */ 129 PetscFunctionReturn(0); 130 }