1 #include <petsctaolinesearch.h> 2 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h> 3 4 #define LMVM_STEP_BFGS 0 5 #define LMVM_STEP_GRAD 1 6 7 static PetscErrorCode TaoSolve_LMVM(Tao tao) 8 { 9 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 10 PetscReal f, fold, gdx, gnorm; 11 PetscReal step = 1.0; 12 PetscInt stepType = LMVM_STEP_GRAD, nupdates; 13 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 14 15 PetscFunctionBegin; 16 17 if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n")); 18 19 /* Check convergence criteria */ 20 PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient)); 21 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 22 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, step)); 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 (!lmP->recycle) { 33 lmP->bfgs = 0; 34 lmP->grad = 0; 35 PetscCall(MatLMVMReset(lmP->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 PetscTryTypeMethod(tao, update, tao->niter, tao->user_update); 42 43 /* Compute direction */ 44 if (lmP->H0) { 45 PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0)); 46 stepType = LMVM_STEP_BFGS; 47 } 48 PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient)); 49 PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D)); 50 PetscCall(MatLMVMGetUpdateCount(lmP->M, &nupdates)); 51 if (nupdates > 0) stepType = LMVM_STEP_BFGS; 52 53 /* Check for success (descent direction) */ 54 PetscCall(VecDot(lmP->D, tao->gradient, &gdx)); 55 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 56 /* Step is not descent or direction produced not a number 57 We can assert bfgsUpdates > 1 in this case because 58 the first solve produces the scaled gradient direction, 59 which is guaranteed to be descent 60 61 Use steepest descent direction (scaled) 62 */ 63 64 PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE)); 65 PetscCall(MatLMVMClearJ0(lmP->M)); 66 PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient)); 67 PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D)); 68 69 /* On a reset, the direction cannot be not a number; it is a 70 scaled gradient step. No need to check for this condition. */ 71 stepType = LMVM_STEP_GRAD; 72 } 73 PetscCall(VecScale(lmP->D, -1.0)); 74 75 /* Perform the linesearch */ 76 fold = f; 77 PetscCall(VecCopy(tao->solution, lmP->Xold)); 78 PetscCall(VecCopy(tao->gradient, lmP->Gold)); 79 80 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status)); 81 PetscCall(TaoAddLineSearchCounts(tao)); 82 83 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) { 84 /* Reset factors and use scaled gradient step */ 85 f = fold; 86 PetscCall(VecCopy(lmP->Xold, tao->solution)); 87 PetscCall(VecCopy(lmP->Gold, tao->gradient)); 88 89 /* Failed to obtain acceptable iterate with BFGS step */ 90 /* Attempt to use the scaled gradient direction */ 91 92 PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE)); 93 PetscCall(MatLMVMClearJ0(lmP->M)); 94 PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient)); 95 PetscCall(MatSolve(lmP->M, tao->solution, tao->gradient)); 96 97 /* On a reset, the direction cannot be not a number; it is a 98 scaled gradient step. No need to check for this condition. */ 99 stepType = LMVM_STEP_GRAD; 100 PetscCall(VecScale(lmP->D, -1.0)); 101 102 /* Perform the linesearch */ 103 PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status)); 104 PetscCall(TaoAddLineSearchCounts(tao)); 105 } 106 107 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 108 /* Failed to find an improving point */ 109 f = fold; 110 PetscCall(VecCopy(lmP->Xold, tao->solution)); 111 PetscCall(VecCopy(lmP->Gold, tao->gradient)); 112 step = 0.0; 113 tao->reason = TAO_DIVERGED_LS_FAILURE; 114 } else { 115 /* LS found valid step, so tally up step type */ 116 switch (stepType) { 117 case LMVM_STEP_BFGS: 118 ++lmP->bfgs; 119 break; 120 case LMVM_STEP_GRAD: 121 ++lmP->grad; 122 break; 123 default: 124 break; 125 } 126 /* Compute new gradient norm */ 127 PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm)); 128 } 129 130 /* Check convergence */ 131 tao->niter++; 132 PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its)); 133 PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step)); 134 PetscUseTypeMethod(tao, convergencetest, tao->cnvP); 135 } 136 PetscFunctionReturn(0); 137 } 138 139 static PetscErrorCode TaoSetUp_LMVM(Tao tao) 140 { 141 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 142 PetscInt n, N; 143 PetscBool is_set, is_spd; 144 145 PetscFunctionBegin; 146 /* Existence of tao->solution checked in TaoSetUp() */ 147 if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient)); 148 if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection)); 149 if (!lmP->D) PetscCall(VecDuplicate(tao->solution, &lmP->D)); 150 if (!lmP->Xold) PetscCall(VecDuplicate(tao->solution, &lmP->Xold)); 151 if (!lmP->Gold) PetscCall(VecDuplicate(tao->solution, &lmP->Gold)); 152 153 /* Create matrix for the limited memory approximation */ 154 PetscCall(VecGetLocalSize(tao->solution, &n)); 155 PetscCall(VecGetSize(tao->solution, &N)); 156 PetscCall(MatSetSizes(lmP->M, n, n, N, N)); 157 PetscCall(MatLMVMAllocate(lmP->M, tao->solution, tao->gradient)); 158 PetscCall(MatIsSPDKnown(lmP->M, &is_set, &is_spd)); 159 PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite."); 160 161 /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */ 162 if (lmP->H0) PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0)); 163 PetscFunctionReturn(0); 164 } 165 166 /* ---------------------------------------------------------- */ 167 static PetscErrorCode TaoDestroy_LMVM(Tao tao) 168 { 169 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 170 171 PetscFunctionBegin; 172 if (tao->setupcalled) { 173 PetscCall(VecDestroy(&lmP->Xold)); 174 PetscCall(VecDestroy(&lmP->Gold)); 175 PetscCall(VecDestroy(&lmP->D)); 176 } 177 PetscCall(MatDestroy(&lmP->M)); 178 if (lmP->H0) PetscCall(PetscObjectDereference((PetscObject)lmP->H0)); 179 PetscCall(PetscFree(tao->data)); 180 PetscFunctionReturn(0); 181 } 182 183 /*------------------------------------------------------------*/ 184 static PetscErrorCode TaoSetFromOptions_LMVM(Tao tao, PetscOptionItems *PetscOptionsObject) 185 { 186 TAO_LMVM *lm = (TAO_LMVM *)tao->data; 187 188 PetscFunctionBegin; 189 PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for unconstrained optimization"); 190 PetscCall(PetscOptionsBool("-tao_lmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", lm->recycle, &lm->recycle, NULL)); 191 PetscCall(TaoLineSearchSetFromOptions(tao->linesearch)); 192 PetscCall(MatSetFromOptions(lm->M)); 193 PetscOptionsHeadEnd(); 194 PetscFunctionReturn(0); 195 } 196 197 /*------------------------------------------------------------*/ 198 static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer) 199 { 200 TAO_LMVM *lm = (TAO_LMVM *)tao->data; 201 PetscBool isascii; 202 PetscInt recycled_its; 203 204 PetscFunctionBegin; 205 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii)); 206 if (isascii) { 207 PetscCall(PetscViewerASCIIPrintf(viewer, " Gradient steps: %" PetscInt_FMT "\n", lm->grad)); 208 if (lm->recycle) { 209 PetscCall(PetscViewerASCIIPrintf(viewer, " Recycle: on\n")); 210 recycled_its = lm->bfgs + lm->grad; 211 PetscCall(PetscViewerASCIIPrintf(viewer, " Total recycled iterations: %" PetscInt_FMT "\n", recycled_its)); 212 } 213 } 214 PetscFunctionReturn(0); 215 } 216 217 /* ---------------------------------------------------------- */ 218 219 /*MC 220 TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton 221 optimization solver for unconstrained minimization. It solves 222 the Newton step 223 Hkdk = - gk 224 225 using an approximation Bk in place of Hk, where Bk is composed using 226 the BFGS update formula. A More-Thuente line search is then used 227 to computed the steplength in the dk direction 228 229 Options Database Keys: 230 + -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls 231 - -tao_lmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation 232 233 Level: beginner 234 M*/ 235 236 PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao) 237 { 238 TAO_LMVM *lmP; 239 const char *morethuente_type = TAOLINESEARCHMT; 240 241 PetscFunctionBegin; 242 tao->ops->setup = TaoSetUp_LMVM; 243 tao->ops->solve = TaoSolve_LMVM; 244 tao->ops->view = TaoView_LMVM; 245 tao->ops->setfromoptions = TaoSetFromOptions_LMVM; 246 tao->ops->destroy = TaoDestroy_LMVM; 247 248 PetscCall(PetscNew(&lmP)); 249 lmP->D = NULL; 250 lmP->M = NULL; 251 lmP->Xold = NULL; 252 lmP->Gold = NULL; 253 lmP->H0 = NULL; 254 lmP->recycle = PETSC_FALSE; 255 256 tao->data = (void *)lmP; 257 /* Override default settings (unless already changed) */ 258 if (!tao->max_it_changed) tao->max_it = 2000; 259 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 260 261 PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch)); 262 PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1)); 263 PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type)); 264 PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao)); 265 PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix)); 266 267 PetscCall(KSPInitializePackage()); 268 PetscCall(MatCreate(((PetscObject)tao)->comm, &lmP->M)); 269 PetscCall(PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1)); 270 PetscCall(MatSetType(lmP->M, MATLMVMBFGS)); 271 PetscCall(MatSetOptionsPrefix(lmP->M, "tao_lmvm_")); 272 PetscFunctionReturn(0); 273 } 274