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