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(PETSC_SUCCESS); 25 } 26 27 /*------------------------------------------------------------*/ 28 static PetscErrorCode TaoSetFromOptions_TRON(Tao tao, PetscOptionItems PetscOptionsObject) 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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 infinity 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) tao->trust = PetscMax(tron->gnorm * tron->gnorm, 1.0); 101 102 /* Initialize step sizes for the line searches */ 103 tron->pgstepsize = 1.0; 104 tron->stepsize = tao->trust; 105 106 tao->reason = TAO_CONTINUE_ITERATING; 107 PetscCall(TaoLogConvergenceHistory(tao, tron->f, tron->gnorm, 0.0, tao->ksp_its)); 108 PetscCall(TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, tron->stepsize)); 109 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 110 while (tao->reason == TAO_CONTINUE_ITERATING) { 111 /* Call general purpose update function */ 112 if (tao->ops->update) { 113 PetscUseTypeMethod(tao, update, tao->niter, tao->user_update); 114 PetscCall(TaoComputeObjective(tao, tao->solution, &tron->f)); 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; 125 delta = tao->trust; 126 tron->n_free_last = tron->n_free; 127 PetscCall(TaoComputeHessian(tao, tao->solution, tao->hessian, tao->hessian_pre)); 128 129 /* Generate index set (IS) of which bound constraints are active */ 130 PetscCall(ISDestroy(&tron->Free_Local)); 131 PetscCall(VecWhichInactive(tao->XL, tao->solution, tao->gradient, tao->XU, PETSC_TRUE, &tron->Free_Local)); 132 PetscCall(ISGetSize(tron->Free_Local, &tron->n_free)); 133 134 /* If no free variables */ 135 if (tron->n_free == 0) { 136 PetscCall(VecNorm(tao->gradient, NORM_2, &tron->gnorm)); 137 PetscCall(TaoLogConvergenceHistory(tao, tron->f, tron->gnorm, 0.0, tao->ksp_its)); 138 PetscCall(TaoMonitor(tao, tao->niter, tron->f, tron->gnorm, 0.0, delta)); 139 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 140 if (!tao->reason) tao->reason = TAO_CONVERGED_STEPTOL; 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 /* Approximately solve the reduced linear system */ 160 PetscCall(KSPCGSetRadius(tao->ksp, delta)); 161 162 PetscCall(KSPSolve(tao->ksp, tron->R, tron->DXFree)); 163 PetscCall(KSPGetIterationNumber(tao->ksp, &its)); 164 tao->ksp_its += its; 165 tao->ksp_tot_its += its; 166 PetscCall(VecSet(tao->stepdirection, 0.0)); 167 168 /* Add dxfree matrix to compute step direction vector */ 169 PetscCall(VecISAXPY(tao->stepdirection, tron->Free_Local, 1.0, tron->DXFree)); 170 171 PetscCall(VecDot(tao->gradient, tao->stepdirection, &gdx)); 172 PetscCall(VecCopy(tao->solution, tron->X_New)); 173 PetscCall(VecCopy(tao->gradient, tron->G_New)); 174 175 stepsize = 1.0; 176 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 } else if (delta <= 1e-30) { 220 break; 221 } else { 222 delta /= 4.0; 223 } 224 } /* end linear solve loop */ 225 226 tron->f = f; 227 tron->actred = actred; 228 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 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 233 } /* END MAIN LOOP */ 234 PetscFunctionReturn(PETSC_SUCCESS); 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 250 PetscFunctionBegin; 251 for (i = 0; i < tron->maxgpits; ++i) { 252 if (-actred <= (tron->pg_ftol) * actred_max) break; 253 254 ++tron->gp_iterates; 255 ++tron->total_gp_its; 256 f_new = tron->f; 257 258 PetscCall(VecCopy(tao->gradient, tao->stepdirection)); 259 PetscCall(VecScale(tao->stepdirection, -1.0)); 260 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, tron->pgstepsize)); 261 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f_new, tao->gradient, tao->stepdirection, &tron->pgstepsize, &ls_reason)); 262 PetscCall(TaoAddLineSearchCounts(tao)); 263 264 PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, tao->gradient)); 265 PetscCall(VecNorm(tao->gradient, NORM_2, &tron->gnorm)); 266 267 /* Update the iterate */ 268 actred = f_new - tron->f; 269 actred_max = PetscMax(actred_max, -(f_new - tron->f)); 270 tron->f = f_new; 271 } 272 PetscFunctionReturn(PETSC_SUCCESS); 273 } 274 275 static PetscErrorCode TaoComputeDual_TRON(Tao tao, Vec DXL, Vec DXU) 276 { 277 TAO_TRON *tron = (TAO_TRON *)tao->data; 278 279 PetscFunctionBegin; 280 PetscValidHeaderSpecific(tao, TAO_CLASSID, 1); 281 PetscValidHeaderSpecific(DXL, VEC_CLASSID, 2); 282 PetscValidHeaderSpecific(DXU, VEC_CLASSID, 3); 283 PetscCheck(tron->Work && tao->gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist."); 284 285 PetscCall(VecBoundGradientProjection(tao->gradient, tao->solution, tao->XL, tao->XU, tron->Work)); 286 PetscCall(VecCopy(tron->Work, DXL)); 287 PetscCall(VecAXPY(DXL, -1.0, tao->gradient)); 288 PetscCall(VecSet(DXU, 0.0)); 289 PetscCall(VecPointwiseMax(DXL, DXL, DXU)); 290 291 PetscCall(VecCopy(tao->gradient, DXU)); 292 PetscCall(VecAXPY(DXU, -1.0, tron->Work)); 293 PetscCall(VecSet(tron->Work, 0.0)); 294 PetscCall(VecPointwiseMin(DXU, tron->Work, DXU)); 295 PetscFunctionReturn(PETSC_SUCCESS); 296 } 297 298 /*------------------------------------------------------------*/ 299 /*MC 300 TAOTRON - The TRON algorithm is an active-set Newton trust region method 301 for bound-constrained minimization. 302 303 Options Database Keys: 304 + -tao_tron_maxgpits - maximum number of gradient projections per TRON iterate 305 - -tao_subset_type - "subvec","mask","matrix-free", strategies for handling active-sets 306 307 Level: beginner 308 M*/ 309 PETSC_EXTERN PetscErrorCode TaoCreate_TRON(Tao tao) 310 { 311 TAO_TRON *tron; 312 const char *morethuente_type = TAOLINESEARCHMT; 313 314 PetscFunctionBegin; 315 tao->ops->setup = TaoSetup_TRON; 316 tao->ops->solve = TaoSolve_TRON; 317 tao->ops->view = TaoView_TRON; 318 tao->ops->setfromoptions = TaoSetFromOptions_TRON; 319 tao->ops->destroy = TaoDestroy_TRON; 320 tao->ops->computedual = TaoComputeDual_TRON; 321 322 PetscCall(PetscNew(&tron)); 323 tao->data = (void *)tron; 324 325 /* Override default settings (unless already changed) */ 326 PetscCall(TaoParametersInitialize(tao)); 327 PetscObjectParameterSetDefault(tao, max_it, 50); 328 PetscObjectParameterSetDefault(tao, trust0, 1.0); 329 PetscObjectParameterSetDefault(tao, steptol, 0.0); 330 331 /* Initialize pointers and variables */ 332 tron->n = 0; 333 tron->maxgpits = 3; 334 tron->pg_ftol = 0.001; 335 336 tron->eta1 = 1.0e-4; 337 tron->eta2 = 0.25; 338 tron->eta3 = 0.50; 339 tron->eta4 = 0.90; 340 341 tron->sigma1 = 0.5; 342 tron->sigma2 = 2.0; 343 tron->sigma3 = 4.0; 344 345 tron->gp_iterates = 0; /* Cumulative number */ 346 tron->total_gp_its = 0; 347 tron->n_free = 0; 348 349 tron->DXFree = NULL; 350 tron->R = NULL; 351 tron->X_New = NULL; 352 tron->G_New = NULL; 353 tron->Work = NULL; 354 tron->Free_Local = NULL; 355 tron->H_sub = NULL; 356 tron->Hpre_sub = NULL; 357 tao->subset_type = TAO_SUBSET_SUBVEC; 358 359 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 360 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 361 PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type)); 362 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 363 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix)); 364 365 PetscCall(KSPCreate(((PetscObject)tao)->comm, &tao->ksp)); 366 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->ksp, (PetscObject)tao, 1)); 367 PetscCall(KSPSetOptionsPrefix(tao->ksp, tao->hdr.prefix)); 368 PetscCall(KSPSetType(tao->ksp, KSPSTCG)); 369 PetscFunctionReturn(PETSC_SUCCESS); 370 } 371