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