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