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