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