1 #include <../src/tao/bound/impls/tron/tron.h> 2 #include <../src/tao/matrix/submatfree.h> 3 4 /* TRON Routines */ 5 static PetscErrorCode TronGradientProjections(Tao,TAO_TRON*); 6 /*------------------------------------------------------------*/ 7 static PetscErrorCode TaoDestroy_TRON(Tao tao) 8 { 9 TAO_TRON *tron = (TAO_TRON *)tao->data; 10 11 PetscFunctionBegin; 12 PetscCall(VecDestroy(&tron->X_New)); 13 PetscCall(VecDestroy(&tron->G_New)); 14 PetscCall(VecDestroy(&tron->Work)); 15 PetscCall(VecDestroy(&tron->DXFree)); 16 PetscCall(VecDestroy(&tron->R)); 17 PetscCall(VecDestroy(&tron->diag)); 18 PetscCall(VecScatterDestroy(&tron->scatter)); 19 PetscCall(ISDestroy(&tron->Free_Local)); 20 PetscCall(MatDestroy(&tron->H_sub)); 21 PetscCall(MatDestroy(&tron->Hpre_sub)); 22 PetscCall(KSPDestroy(&tao->ksp)); 23 PetscCall(PetscFree(tao->data)); 24 PetscFunctionReturn(0); 25 } 26 27 /*------------------------------------------------------------*/ 28 static PetscErrorCode TaoSetFromOptions_TRON(PetscOptionItems *PetscOptionsObject,Tao tao) 29 { 30 TAO_TRON *tron = (TAO_TRON *)tao->data; 31 PetscBool flg; 32 33 PetscFunctionBegin; 34 PetscOptionsHeadBegin(PetscOptionsObject,"Newton Trust Region Method for bound constrained optimization"); 35 PetscCall(PetscOptionsInt("-tao_tron_maxgpits","maximum number of gradient projections per TRON iterate","TaoSetMaxGPIts",tron->maxgpits,&tron->maxgpits,&flg)); 36 PetscOptionsHeadEnd(); 37 PetscCall(KSPSetFromOptions(tao->ksp)); 38 PetscFunctionReturn(0); 39 } 40 41 /*------------------------------------------------------------*/ 42 static PetscErrorCode TaoView_TRON(Tao tao, PetscViewer viewer) 43 { 44 TAO_TRON *tron = (TAO_TRON *)tao->data; 45 PetscBool isascii; 46 47 PetscFunctionBegin; 48 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii)); 49 if (isascii) { 50 PetscCall(PetscViewerASCIIPrintf(viewer,"Total PG its: %" PetscInt_FMT ",",tron->total_gp_its)); 51 PetscCall(PetscViewerASCIIPrintf(viewer,"PG tolerance: %g \n",(double)tron->pg_ftol)); 52 } 53 PetscFunctionReturn(0); 54 } 55 56 /* ---------------------------------------------------------- */ 57 static PetscErrorCode TaoSetup_TRON(Tao tao) 58 { 59 TAO_TRON *tron = (TAO_TRON *)tao->data; 60 61 PetscFunctionBegin; 62 /* Allocate some arrays */ 63 PetscCall(VecDuplicate(tao->solution, &tron->diag)); 64 PetscCall(VecDuplicate(tao->solution, &tron->X_New)); 65 PetscCall(VecDuplicate(tao->solution, &tron->G_New)); 66 PetscCall(VecDuplicate(tao->solution, &tron->Work)); 67 PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 68 PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 69 PetscFunctionReturn(0); 70 } 71 72 static PetscErrorCode TaoSolve_TRON(Tao tao) 73 { 74 TAO_TRON *tron = (TAO_TRON *)tao->data; 75 PetscInt its; 76 TaoLineSearchConvergedReason ls_reason = TAOLINESEARCH_CONTINUE_ITERATING; 77 PetscReal prered,actred,delta,f,f_new,rhok,gdx,xdiff,stepsize; 78 79 PetscFunctionBegin; 80 tron->pgstepsize = 1.0; 81 tao->trust = tao->trust0; 82 /* Project the current point onto the feasible set */ 83 PetscCall(TaoComputeVariableBounds(tao)); 84 PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU)); 85 86 /* Project the initial point onto the feasible region */ 87 PetscCall(VecMedian(tao->XL,tao->solution,tao->XU,tao->solution)); 88 89 /* Compute the objective function and gradient */ 90 PetscCall(TaoComputeObjectiveAndGradient(tao,tao->solution,&tron->f,tao->gradient)); 91 PetscCall(VecNorm(tao->gradient,NORM_2,&tron->gnorm)); 92 PetscCheck(!PetscIsInfOrNanReal(tron->f) && !PetscIsInfOrNanReal(tron->gnorm),PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 93 94 /* Project the gradient and calculate the norm */ 95 PetscCall(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient)); 96 PetscCall(VecNorm(tao->gradient,NORM_2,&tron->gnorm)); 97 98 /* Initialize trust region radius */ 99 tao->trust=tao->trust0; 100 if (tao->trust <= 0) { 101 tao->trust=PetscMax(tron->gnorm*tron->gnorm,1.0); 102 } 103 104 /* Initialize step sizes for the line searches */ 105 tron->pgstepsize=1.0; 106 tron->stepsize=tao->trust; 107 108 tao->reason = TAO_CONTINUE_ITERATING; 109 PetscCall(TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its)); 110 PetscCall(TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,tron->stepsize)); 111 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 112 while (tao->reason==TAO_CONTINUE_ITERATING) { 113 /* Call general purpose update function */ 114 if (tao->ops->update) PetscCall((*tao->ops->update)(tao, tao->niter, tao->user_update)); 115 116 /* Perform projected gradient iterations */ 117 PetscCall(TronGradientProjections(tao,tron)); 118 119 PetscCall(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient)); 120 PetscCall(VecNorm(tao->gradient,NORM_2,&tron->gnorm)); 121 122 tao->ksp_its=0; 123 f=tron->f; delta=tao->trust; 124 tron->n_free_last = tron->n_free; 125 PetscCall(TaoComputeHessian(tao,tao->solution,tao->hessian,tao->hessian_pre)); 126 127 /* Generate index set (IS) of which bound constraints are active */ 128 PetscCall(ISDestroy(&tron->Free_Local)); 129 PetscCall(VecWhichInactive(tao->XL,tao->solution,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local)); 130 PetscCall(ISGetSize(tron->Free_Local, &tron->n_free)); 131 132 /* If no free variables */ 133 if (tron->n_free == 0) { 134 PetscCall(VecNorm(tao->gradient,NORM_2,&tron->gnorm)); 135 PetscCall(TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its)); 136 PetscCall(TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,delta)); 137 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 138 if (!tao->reason) { 139 tao->reason = TAO_CONVERGED_STEPTOL; 140 } 141 break; 142 } 143 /* use free_local to mask/submat gradient, hessian, stepdirection */ 144 PetscCall(TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->R)); 145 PetscCall(TaoVecGetSubVec(tao->gradient,tron->Free_Local,tao->subset_type,0.0,&tron->DXFree)); 146 PetscCall(VecSet(tron->DXFree,0.0)); 147 PetscCall(VecScale(tron->R, -1.0)); 148 PetscCall(TaoMatGetSubMat(tao->hessian, tron->Free_Local, tron->diag, tao->subset_type, &tron->H_sub)); 149 if (tao->hessian == tao->hessian_pre) { 150 PetscCall(MatDestroy(&tron->Hpre_sub)); 151 PetscCall(PetscObjectReference((PetscObject)(tron->H_sub))); 152 tron->Hpre_sub = tron->H_sub; 153 } else { 154 PetscCall(TaoMatGetSubMat(tao->hessian_pre, tron->Free_Local, tron->diag, tao->subset_type,&tron->Hpre_sub)); 155 } 156 PetscCall(KSPReset(tao->ksp)); 157 PetscCall(KSPSetOperators(tao->ksp, tron->H_sub, tron->Hpre_sub)); 158 while (1) { 159 160 /* Approximately solve the reduced linear system */ 161 PetscCall(KSPCGSetRadius(tao->ksp,delta)); 162 163 PetscCall(KSPSolve(tao->ksp, tron->R, tron->DXFree)); 164 PetscCall(KSPGetIterationNumber(tao->ksp,&its)); 165 tao->ksp_its+=its; 166 tao->ksp_tot_its+=its; 167 PetscCall(VecSet(tao->stepdirection,0.0)); 168 169 /* Add dxfree matrix to compute step direction vector */ 170 PetscCall(VecISAXPY(tao->stepdirection,tron->Free_Local,1.0,tron->DXFree)); 171 172 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 173 PetscCall(VecCopy(tao->solution, tron->X_New)); 174 PetscCall(VecCopy(tao->gradient, tron->G_New)); 175 176 stepsize=1.0;f_new=f; 177 178 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch,1.0)); 179 PetscCall(TaoLineSearchApply(tao->linesearch, tron->X_New, &f_new, tron->G_New, tao->stepdirection,&stepsize,&ls_reason)); 180 PetscCall(TaoAddLineSearchCounts(tao)); 181 182 PetscCall(MatMult(tao->hessian, tao->stepdirection, tron->Work)); 183 PetscCall(VecAYPX(tron->Work, 0.5, tao->gradient)); 184 PetscCall(VecDot(tao->stepdirection, tron->Work, &prered)); 185 actred = f_new - f; 186 if ((PetscAbsScalar(actred) <= 1e-6) && (PetscAbsScalar(prered) <= 1e-6)) { 187 rhok = 1.0; 188 } else if (actred<0) { 189 rhok=PetscAbs(-actred/prered); 190 } else { 191 rhok=0.0; 192 } 193 194 /* Compare actual improvement to the quadratic model */ 195 if (rhok > tron->eta1) { /* Accept the point */ 196 /* d = x_new - x */ 197 PetscCall(VecCopy(tron->X_New, tao->stepdirection)); 198 PetscCall(VecAXPY(tao->stepdirection, -1.0, tao->solution)); 199 200 PetscCall(VecNorm(tao->stepdirection, NORM_2, &xdiff)); 201 xdiff *= stepsize; 202 203 /* Adjust trust region size */ 204 if (rhok < tron->eta2) { 205 delta = PetscMin(xdiff,delta)*tron->sigma1; 206 } else if (rhok > tron->eta4) { 207 delta= PetscMin(xdiff,delta)*tron->sigma3; 208 } else if (rhok > tron->eta3) { 209 delta=PetscMin(xdiff,delta)*tron->sigma2; 210 } 211 PetscCall(VecBoundGradientProjection(tron->G_New,tron->X_New, tao->XL, tao->XU, tao->gradient)); 212 PetscCall(ISDestroy(&tron->Free_Local)); 213 PetscCall(VecWhichInactive(tao->XL,tron->X_New,tao->gradient,tao->XU,PETSC_TRUE,&tron->Free_Local)); 214 f=f_new; 215 PetscCall(VecNorm(tao->gradient,NORM_2,&tron->gnorm)); 216 PetscCall(VecCopy(tron->X_New, tao->solution)); 217 PetscCall(VecCopy(tron->G_New, tao->gradient)); 218 break; 219 } 220 else if (delta <= 1e-30) { 221 break; 222 } 223 else { 224 delta /= 4.0; 225 } 226 } /* end linear solve loop */ 227 228 tron->f=f; tron->actred=actred; tao->trust=delta; 229 tao->niter++; 230 PetscCall(TaoLogConvergenceHistory(tao,tron->f,tron->gnorm,0.0,tao->ksp_its)); 231 PetscCall(TaoMonitor(tao,tao->niter,tron->f,tron->gnorm,0.0,stepsize)); 232 PetscCall((*tao->ops->convergencetest)(tao,tao->cnvP)); 233 } /* END MAIN LOOP */ 234 PetscFunctionReturn(0); 235 } 236 237 static PetscErrorCode TronGradientProjections(Tao tao,TAO_TRON *tron) 238 { 239 PetscInt i; 240 TaoLineSearchConvergedReason ls_reason; 241 PetscReal actred=-1.0,actred_max=0.0; 242 PetscReal f_new; 243 /* 244 The gradient and function value passed into and out of this 245 routine should be current and correct. 246 247 The free, active, and binding variables should be already identified 248 */ 249 PetscFunctionBegin; 250 251 for (i=0;i<tron->maxgpits;++i) { 252 253 if (-actred <= (tron->pg_ftol)*actred_max) break; 254 255 ++tron->gp_iterates; 256 ++tron->total_gp_its; 257 f_new=tron->f; 258 259 PetscCall(VecCopy(tao->gradient,tao->stepdirection)); 260 PetscCall(VecScale(tao->stepdirection,-1.0)); 261 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch,tron->pgstepsize)); 262 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection,&tron->pgstepsize, &ls_reason)); 263 PetscCall(TaoAddLineSearchCounts(tao)); 264 265 PetscCall(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tao->gradient)); 266 PetscCall(VecNorm(tao->gradient,NORM_2,&tron->gnorm)); 267 268 /* Update the iterate */ 269 actred = f_new - tron->f; 270 actred_max = PetscMax(actred_max,-(f_new - tron->f)); 271 tron->f = f_new; 272 } 273 PetscFunctionReturn(0); 274 } 275 276 static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) 277 { 278 279 TAO_TRON *tron = (TAO_TRON *)tao->data; 280 281 PetscFunctionBegin; 282 PetscValidHeaderSpecific(tao,TAO_CLASSID,1); 283 PetscValidHeaderSpecific(DXL,VEC_CLASSID,2); 284 PetscValidHeaderSpecific(DXU,VEC_CLASSID,3); 285 PetscCheck(tron->Work && tao->gradient,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist."); 286 287 PetscCall(VecBoundGradientProjection(tao->gradient,tao->solution,tao->XL,tao->XU,tron->Work)); 288 PetscCall(VecCopy(tron->Work,DXL)); 289 PetscCall(VecAXPY(DXL,-1.0,tao->gradient)); 290 PetscCall(VecSet(DXU,0.0)); 291 PetscCall(VecPointwiseMax(DXL,DXL,DXU)); 292 293 PetscCall(VecCopy(tao->gradient,DXU)); 294 PetscCall(VecAXPY(DXU,-1.0,tron->Work)); 295 PetscCall(VecSet(tron->Work,0.0)); 296 PetscCall(VecPointwiseMin(DXU,tron->Work,DXU)); 297 PetscFunctionReturn(0); 298 } 299 300 /*------------------------------------------------------------*/ 301 /*MC 302 TAOTRON - The TRON algorithm is an active-set Newton trust region method 303 for bound-constrained minimization. 304 305 Options Database Keys: 306 + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate 307 - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets 308 309 Level: beginner 310 M*/ 311 PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao) 312 { 313 TAO_TRON *tron; 314 const char *morethuente_type = TAOLINESEARCHMT; 315 316 PetscFunctionBegin; 317 tao->ops->setup = TaoSetup_TRON; 318 tao->ops->solve = TaoSolve_TRON; 319 tao->ops->view = TaoView_TRON; 320 tao->ops->setfromoptions = TaoSetFromOptions_TRON; 321 tao->ops->destroy = TaoDestroy_TRON; 322 tao->ops->computedual = TaoComputeDual_TRON; 323 324 PetscCall(PetscNewLog(tao,&tron)); 325 tao->data = (void*)tron; 326 327 /* Override default settings (unless already changed) */ 328 if (!tao->max_it_changed) tao->max_it = 50; 329 if (!tao->trust0_changed) tao->trust0 = 1.0; 330 if (!tao->steptol_changed) tao->steptol = 0.0; 331 332 /* Initialize pointers and variables */ 333 tron->n = 0; 334 tron->maxgpits = 3; 335 tron->pg_ftol = 0.001; 336 337 tron->eta1 = 1.0e-4; 338 tron->eta2 = 0.25; 339 tron->eta3 = 0.50; 340 tron->eta4 = 0.90; 341 342 tron->sigma1 = 0.5; 343 tron->sigma2 = 2.0; 344 tron->sigma3 = 4.0; 345 346 tron->gp_iterates = 0; /* Cumulative number */ 347 tron->total_gp_its = 0; 348 tron->n_free = 0; 349 350 tron->DXFree=NULL; 351 tron->R=NULL; 352 tron->X_New=NULL; 353 tron->G_New=NULL; 354 tron->Work=NULL; 355 tron->Free_Local=NULL; 356 tron->H_sub=NULL; 357 tron->Hpre_sub=NULL; 358 tao->subset_type = TAO_SUBSET_SUBVEC; 359 360 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 361 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 362 PetscCall(TaoLineSearchSetType(tao->linesearch,morethuente_type)); 363 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch,tao)); 364 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix)); 365 366 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 367 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 368 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix)); 369 PetscCall(KSPSetType(tao->ksp,KSPSTCG)); 370 PetscFunctionReturn(0); 371 } 372