1 #include <petsctaolinesearch.h> /*I "petsctaolinesearch.h" I*/ 2 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h> 3 #include <../src/tao/bound/impls/blmvm/blmvm.h> 4 5 /*------------------------------------------------------------*/ 6 static PetscErrorCode TaoSolve_BLMVM(Tao tao) { 7 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 8 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 9 PetscReal f, fold, gdx, gnorm, gnorm2; 10 PetscReal stepsize = 1.0, delta; 11 12 PetscFunctionBegin; 13 /* Project initial point onto bounds */ 14 PetscCall(TaoComputeVariableBounds(tao)); 15 PetscCall(VecMedian(tao->XL, tao->solution, tao->XU, tao->solution)); 16 PetscCall(TaoLineSearchSetVariableBounds(tao->linesearch, tao->XL, tao->XU)); 17 18 /* Check convergence criteria */ 19 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, blmP->unprojected_gradient)); 20 PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient)); 21 22 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 23 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 24 25 tao->reason = TAO_CONTINUE_ITERATING; 26 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 27 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize)); 28 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 29 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 30 31 /* Set counter for gradient/reset steps */ 32 if (!blmP->recycle) { 33 blmP->grad = 0; 34 blmP->reset = 0; 35 PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE)); 36 } 37 38 /* Have not converged; continue with Newton method */ 39 while (tao->reason == TAO_CONTINUE_ITERATING) { 40 /* Call general purpose update function */ 41 if (tao->ops->update) { 42 PetscUseTypeMethod(tao, update, tao->niter, tao->user_update); 43 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient)); 44 } 45 /* Compute direction */ 46 gnorm2 = gnorm * gnorm; 47 if (gnorm2 == 0.0) gnorm2 = PETSC_MACHINE_EPSILON; 48 if (f == 0.0) { 49 delta = 2.0 / gnorm2; 50 } else { 51 delta = 2.0 * PetscAbsScalar(f) / gnorm2; 52 } 53 PetscCall(MatLMVMSymBroydenSetDelta(blmP->M, delta)); 54 PetscCall(MatLMVMUpdate(blmP->M, tao->solution, tao->gradient)); 55 PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection)); 56 PetscCall(VecBoundGradientProjection(tao->stepdirection, tao->solution, tao->XL, tao->XU, tao->gradient)); 57 58 /* Check for success (descent direction) */ 59 PetscCall(VecDot(blmP->unprojected_gradient, tao->gradient, &gdx)); 60 if (gdx <= 0) { 61 /* Step is not descent or solve was not successful 62 Use steepest descent direction (scaled) */ 63 ++blmP->grad; 64 65 PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE)); 66 PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient)); 67 PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection)); 68 } 69 PetscCall(VecScale(tao->stepdirection, -1.0)); 70 71 /* Perform the linesearch */ 72 fold = f; 73 PetscCall(VecCopy(tao->solution, blmP->Xold)); 74 PetscCall(VecCopy(blmP->unprojected_gradient, blmP->Gold)); 75 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0)); 76 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status)); 77 PetscCall(TaoAddLineSearchCounts(tao)); 78 79 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 80 /* Linesearch failed 81 Reset factors and use scaled (projected) gradient step */ 82 ++blmP->reset; 83 84 f = fold; 85 PetscCall(VecCopy(blmP->Xold, tao->solution)); 86 PetscCall(VecCopy(blmP->Gold, blmP->unprojected_gradient)); 87 88 PetscCall(MatLMVMReset(blmP->M, PETSC_FALSE)); 89 PetscCall(MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient)); 90 PetscCall(MatSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection)); 91 PetscCall(VecScale(tao->stepdirection, -1.0)); 92 93 /* This may be incorrect; linesearch has values for stepmax and stepmin 94 that should be reset. */ 95 PetscCall(TaoLineSearchSetInitialStepLength(tao->linesearch, 1.0)); 96 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status)); 97 PetscCall(TaoAddLineSearchCounts(tao)); 98 99 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 100 tao->reason = TAO_DIVERGED_LS_FAILURE; 101 break; 102 } 103 } 104 105 /* Check for converged */ 106 PetscCall(VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient)); 107 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 108 PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Not-a-Number"); 109 tao->niter++; 110 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 111 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize)); 112 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 113 } 114 PetscFunctionReturn(0); 115 } 116 117 static PetscErrorCode TaoSetup_BLMVM(Tao tao) { 118 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 119 120 PetscFunctionBegin; 121 /* Existence of tao->solution checked in TaoSetup() */ 122 PetscCall(VecDuplicate(tao->solution, &blmP->Xold)); 123 PetscCall(VecDuplicate(tao->solution, &blmP->Gold)); 124 PetscCall(VecDuplicate(tao->solution, &blmP->unprojected_gradient)); 125 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 126 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 127 /* Allocate matrix for the limited memory approximation */ 128 PetscCall(MatLMVMAllocate(blmP->M, tao->solution, blmP->unprojected_gradient)); 129 130 /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */ 131 if (blmP->H0) PetscCall(MatLMVMSetJ0(blmP->M, blmP->H0)); 132 PetscFunctionReturn(0); 133 } 134 135 /* ---------------------------------------------------------- */ 136 static PetscErrorCode TaoDestroy_BLMVM(Tao tao) { 137 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 138 139 PetscFunctionBegin; 140 if (tao->setupcalled) { 141 PetscCall(VecDestroy(&blmP->unprojected_gradient)); 142 PetscCall(VecDestroy(&blmP->Xold)); 143 PetscCall(VecDestroy(&blmP->Gold)); 144 } 145 PetscCall(MatDestroy(&blmP->M)); 146 if (blmP->H0) PetscObjectDereference((PetscObject)blmP->H0); 147 PetscCall(PetscFree(tao->data)); 148 PetscFunctionReturn(0); 149 } 150 151 /*------------------------------------------------------------*/ 152 static PetscErrorCode TaoSetFromOptions_BLMVM(Tao tao, PetscOptionItems *PetscOptionsObject) { 153 TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data; 154 PetscBool is_spd, is_set; 155 156 PetscFunctionBegin; 157 PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for bound constrained optimization"); 158 PetscCall(PetscOptionsBool("-tao_blmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", blmP->recycle, &blmP->recycle, NULL)); 159 PetscOptionsHeadEnd(); 160 PetscCall(MatSetOptionsPrefix(blmP->M, ((PetscObject)tao)->prefix)); 161 PetscCall(MatAppendOptionsPrefix(blmP->M, "tao_blmvm_")); 162 PetscCall(MatSetFromOptions(blmP->M)); 163 PetscCall(MatIsSPDKnown(blmP->M, &is_set, &is_spd)); 164 PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix must be symmetric positive-definite"); 165 PetscFunctionReturn(0); 166 } 167 168 /*------------------------------------------------------------*/ 169 static PetscErrorCode TaoView_BLMVM(Tao tao, PetscViewer viewer) { 170 TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data; 171 PetscBool isascii; 172 173 PetscFunctionBegin; 174 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 175 if (isascii) { 176 PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lmP->grad)); 177 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_INFO)); 178 PetscCall(MatView(lmP->M, viewer)); 179 PetscCall(PetscViewerPopFormat(viewer)); 180 } 181 PetscFunctionReturn(0); 182 } 183 184 static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU) { 185 TAO_BLMVM *blm = (TAO_BLMVM *)tao->data; 186 187 PetscFunctionBegin; 188 PetscValidHeaderSpecific(tao, TAO_CLASSID, 1); 189 PetscValidHeaderSpecific(DXL, VEC_CLASSID, 2); 190 PetscValidHeaderSpecific(DXU, VEC_CLASSID, 3); 191 PetscCheck(tao->gradient && blm->unprojected_gradient, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Dual variables don't exist yet or no longer exist."); 192 193 PetscCall(VecCopy(tao->gradient, DXL)); 194 PetscCall(VecAXPY(DXL, -1.0, blm->unprojected_gradient)); 195 PetscCall(VecSet(DXU, 0.0)); 196 PetscCall(VecPointwiseMax(DXL, DXL, DXU)); 197 198 PetscCall(VecCopy(blm->unprojected_gradient, DXU)); 199 PetscCall(VecAXPY(DXU, -1.0, tao->gradient)); 200 PetscCall(VecAXPY(DXU, 1.0, DXL)); 201 PetscFunctionReturn(0); 202 } 203 204 /* ---------------------------------------------------------- */ 205 /*MC 206 TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method 207 for nonlinear minimization with bound constraints. It is an extension 208 of TAOLMVM 209 210 Options Database Keys: 211 . -tao_lmm_recycle - enable recycling of LMVM information between subsequent TaoSolve calls 212 213 Level: beginner 214 M*/ 215 PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao) { 216 TAO_BLMVM *blmP; 217 const char *morethuente_type = TAOLINESEARCHMT; 218 219 PetscFunctionBegin; 220 tao->ops->setup = TaoSetup_BLMVM; 221 tao->ops->solve = TaoSolve_BLMVM; 222 tao->ops->view = TaoView_BLMVM; 223 tao->ops->setfromoptions = TaoSetFromOptions_BLMVM; 224 tao->ops->destroy = TaoDestroy_BLMVM; 225 tao->ops->computedual = TaoComputeDual_BLMVM; 226 227 PetscCall(PetscNew(&blmP)); 228 blmP->H0 = NULL; 229 blmP->recycle = PETSC_FALSE; 230 tao->data = (void *)blmP; 231 232 /* Override default settings (unless already changed) */ 233 if (!tao->max_it_changed) tao->max_it = 2000; 234 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 235 236 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 237 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 238 PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type)); 239 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 240 241 PetscCall(KSPInitializePackage()); 242 PetscCall(MatCreate(((PetscObject)tao)->comm, &blmP->M)); 243 PetscCall(MatSetType(blmP->M, MATLMVMBFGS)); 244 PetscCall(PetscObjectIncrementTabLevel((PetscObject)blmP->M, (PetscObject)tao, 1)); 245 PetscFunctionReturn(0); 246 } 247 248 /*@ 249 TaoLMVMRecycle - Enable/disable recycling of the QN history between subsequent TaoSolve calls. 250 251 Input Parameters: 252 + tao - the Tao solver context 253 - flg - Boolean flag for recycling (PETSC_TRUE or PETSC_FALSE) 254 255 Level: intermediate 256 @*/ 257 PetscErrorCode TaoLMVMRecycle(Tao tao, PetscBool flg) { 258 TAO_LMVM *lmP; 259 TAO_BLMVM *blmP; 260 PetscBool is_lmvm, is_blmvm; 261 262 PetscFunctionBegin; 263 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm)); 264 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm)); 265 if (is_lmvm) { 266 lmP = (TAO_LMVM *)tao->data; 267 lmP->recycle = flg; 268 } else if (is_blmvm) { 269 blmP = (TAO_BLMVM *)tao->data; 270 blmP->recycle = flg; 271 } 272 PetscFunctionReturn(0); 273 } 274 275 /*@ 276 TaoLMVMSetH0 - Set the initial Hessian for the QN approximation 277 278 Input Parameters: 279 + tao - the Tao solver context 280 - H0 - Mat object for the initial Hessian 281 282 Level: advanced 283 284 .seealso: `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()` 285 @*/ 286 PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0) { 287 TAO_LMVM *lmP; 288 TAO_BLMVM *blmP; 289 PetscBool is_lmvm, is_blmvm; 290 291 PetscFunctionBegin; 292 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm)); 293 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm)); 294 if (is_lmvm) { 295 lmP = (TAO_LMVM *)tao->data; 296 PetscCall(PetscObjectReference((PetscObject)H0)); 297 lmP->H0 = H0; 298 } else if (is_blmvm) { 299 blmP = (TAO_BLMVM *)tao->data; 300 PetscCall(PetscObjectReference((PetscObject)H0)); 301 blmP->H0 = H0; 302 } 303 PetscFunctionReturn(0); 304 } 305 306 /*@ 307 TaoLMVMGetH0 - Get the matrix object for the QN initial Hessian 308 309 Input Parameters: 310 . tao - the Tao solver context 311 312 Output Parameters: 313 . H0 - Mat object for the initial Hessian 314 315 Level: advanced 316 317 .seealso: `TaoLMVMSetH0()`, `TaoLMVMGetH0KSP()` 318 @*/ 319 PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0) { 320 TAO_LMVM *lmP; 321 TAO_BLMVM *blmP; 322 PetscBool is_lmvm, is_blmvm; 323 Mat M; 324 325 PetscFunctionBegin; 326 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm)); 327 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm)); 328 if (is_lmvm) { 329 lmP = (TAO_LMVM *)tao->data; 330 M = lmP->M; 331 } else if (is_blmvm) { 332 blmP = (TAO_BLMVM *)tao->data; 333 M = blmP->M; 334 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM."); 335 PetscCall(MatLMVMGetJ0(M, H0)); 336 PetscFunctionReturn(0); 337 } 338 339 /*@ 340 TaoLMVMGetH0KSP - Get the iterative solver for applying the inverse of the QN initial Hessian 341 342 Input Parameters: 343 . tao - the Tao solver context 344 345 Output Parameters: 346 . ksp - KSP solver context for the initial Hessian 347 348 Level: advanced 349 350 .seealso: `TaoLMVMGetH0()`, `TaoLMVMGetH0KSP()` 351 @*/ 352 PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp) { 353 TAO_LMVM *lmP; 354 TAO_BLMVM *blmP; 355 PetscBool is_lmvm, is_blmvm; 356 Mat M; 357 358 PetscFunctionBegin; 359 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOLMVM, &is_lmvm)); 360 PetscCall(PetscObjectTypeCompare((PetscObject)tao, TAOBLMVM, &is_blmvm)); 361 if (is_lmvm) { 362 lmP = (TAO_LMVM *)tao->data; 363 M = lmP->M; 364 } else if (is_blmvm) { 365 blmP = (TAO_BLMVM *)tao->data; 366 M = blmP->M; 367 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "This routine applies to TAO_LMVM and TAO_BLMVM."); 368 PetscCall(MatLMVMGetJ0KSP(M, ksp)); 369 PetscFunctionReturn(0); 370 } 371