1 #include <petsctaolinesearch.h> 2 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h> 3 #include <../src/tao/bound/impls/blmvm/blmvm.h> 4 5 #define LMVM_STEP_BFGS 0 6 #define LMVM_STEP_GRAD 1 7 8 static PetscErrorCode TaoSolve_LMVM(Tao tao) 9 { 10 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 11 PetscReal f, fold, gdx, gnorm; 12 PetscReal step = 1.0; 13 PetscErrorCode ierr; 14 PetscInt stepType = LMVM_STEP_GRAD, nupdates; 15 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 16 17 PetscFunctionBegin; 18 19 if (tao->XL || tao->XU || tao->ops->computebounds) { 20 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");CHKERRQ(ierr); 21 } 22 23 /* Check convergence criteria */ 24 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr); 25 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 26 27 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN"); 28 29 tao->reason = TAO_CONTINUE_ITERATING; 30 ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 31 ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);CHKERRQ(ierr); 32 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 33 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 34 35 /* Set counter for gradient/reset steps */ 36 if (!lmP->recycle) { 37 lmP->bfgs = 0; 38 lmP->grad = 0; 39 ierr = MatLMVMReset(lmP->M, PETSC_FALSE); CHKERRQ(ierr); 40 } 41 42 /* Have not converged; continue with Newton method */ 43 while (tao->reason == TAO_CONTINUE_ITERATING) { 44 /* Compute direction */ 45 if (lmP->H0) { 46 ierr = MatLMVMSetJ0(lmP->M, lmP->H0);CHKERRQ(ierr); 47 stepType = LMVM_STEP_BFGS; 48 } 49 ierr = MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr); 50 ierr = MatSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr); 51 ierr = MatLMVMGetUpdateCount(lmP->M, &nupdates); CHKERRQ(ierr); 52 if (nupdates > 0) stepType = LMVM_STEP_BFGS; 53 54 /* Check for success (descent direction) */ 55 ierr = VecDot(lmP->D, tao->gradient, &gdx);CHKERRQ(ierr); 56 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 57 /* Step is not descent or direction produced not a number 58 We can assert bfgsUpdates > 1 in this case because 59 the first solve produces the scaled gradient direction, 60 which is guaranteed to be descent 61 62 Use steepest descent direction (scaled) 63 */ 64 65 ierr = MatLMVMResetJ0(lmP->M);CHKERRQ(ierr); 66 ierr = MatLMVMReset(lmP->M, PETSC_FALSE);CHKERRQ(ierr); 67 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 68 ierr = MatSolve(lmP->M,tao->gradient, lmP->D);CHKERRQ(ierr); 69 70 /* On a reset, the direction cannot be not a number; it is a 71 scaled gradient step. No need to check for this condition. */ 72 stepType = LMVM_STEP_GRAD; 73 } 74 ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr); 75 76 /* Perform the linesearch */ 77 fold = f; 78 ierr = VecCopy(tao->solution, lmP->Xold);CHKERRQ(ierr); 79 ierr = VecCopy(tao->gradient, lmP->Gold);CHKERRQ(ierr); 80 81 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);CHKERRQ(ierr); 82 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 83 84 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) { 85 /* Reset factors and use scaled gradient step */ 86 f = fold; 87 ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr); 88 ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr); 89 90 /* Failed to obtain acceptable iterate with BFGS step */ 91 /* Attempt to use the scaled gradient direction */ 92 93 ierr = MatLMVMResetJ0(lmP->M);CHKERRQ(ierr); 94 ierr = MatLMVMReset(lmP->M, PETSC_FALSE);CHKERRQ(ierr); 95 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 96 ierr = MatSolve(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 97 98 /* On a reset, the direction cannot be not a number; it is a 99 scaled gradient step. No need to check for this condition. */ 100 stepType = LMVM_STEP_GRAD; 101 ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr); 102 103 /* Perform the linesearch */ 104 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);CHKERRQ(ierr); 105 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 106 } 107 108 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 109 /* Failed to find an improving point */ 110 f = fold; 111 ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr); 112 ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr); 113 step = 0.0; 114 tao->reason = TAO_DIVERGED_LS_FAILURE; 115 } else { 116 /* LS found valid step, so tally up step type */ 117 switch (stepType) { 118 case LMVM_STEP_BFGS: 119 ++lmP->bfgs; 120 break; 121 case LMVM_STEP_GRAD: 122 ++lmP->grad; 123 break; 124 default: 125 break; 126 } 127 /* Compute new gradient norm */ 128 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 129 } 130 131 /* Check convergence */ 132 tao->niter++; 133 ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 134 ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);CHKERRQ(ierr); 135 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 136 } 137 PetscFunctionReturn(0); 138 } 139 140 static PetscErrorCode TaoSetUp_LMVM(Tao tao) 141 { 142 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 143 PetscInt n,N; 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 /* Existence of tao->solution checked in TaoSetUp() */ 148 if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr); } 149 if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr); } 150 if (!lmP->D) {ierr = VecDuplicate(tao->solution,&lmP->D);CHKERRQ(ierr); } 151 if (!lmP->Xold) {ierr = VecDuplicate(tao->solution,&lmP->Xold);CHKERRQ(ierr); } 152 if (!lmP->Gold) {ierr = VecDuplicate(tao->solution,&lmP->Gold);CHKERRQ(ierr); } 153 154 /* Create matrix for the limited memory approximation */ 155 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 156 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 157 ierr = MatSetSizes(lmP->M, n, n, N, N);CHKERRQ(ierr); 158 ierr = MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr); 159 160 /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */ 161 if (lmP->H0) { 162 ierr = MatLMVMSetJ0(lmP->M, lmP->H0);CHKERRQ(ierr); 163 } 164 165 PetscFunctionReturn(0); 166 } 167 168 /* ---------------------------------------------------------- */ 169 static PetscErrorCode TaoDestroy_LMVM(Tao tao) 170 { 171 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 172 PetscErrorCode ierr; 173 174 PetscFunctionBegin; 175 if (tao->setupcalled) { 176 ierr = VecDestroy(&lmP->Xold);CHKERRQ(ierr); 177 ierr = VecDestroy(&lmP->Gold);CHKERRQ(ierr); 178 ierr = VecDestroy(&lmP->D);CHKERRQ(ierr); 179 } 180 ierr = MatDestroy(&lmP->M);CHKERRQ(ierr); 181 if (lmP->H0) { 182 ierr = PetscObjectDereference((PetscObject)lmP->H0);CHKERRQ(ierr); 183 } 184 ierr = PetscFree(tao->data);CHKERRQ(ierr); 185 186 PetscFunctionReturn(0); 187 } 188 189 /*------------------------------------------------------------*/ 190 static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject,Tao tao) 191 { 192 TAO_LMVM *lm = (TAO_LMVM *)tao->data; 193 PetscErrorCode ierr; 194 195 PetscFunctionBegin; 196 ierr = PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");CHKERRQ(ierr); 197 ierr = PetscOptionsBool("-tao_lmvm_recycle","enable recycling of the BFGS matrix between subsequent TaoSolve() calls","",lm->recycle,&lm->recycle,NULL);CHKERRQ(ierr); 198 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 199 ierr = MatSetFromOptions(lm->M);CHKERRQ(ierr); 200 ierr = PetscOptionsTail();CHKERRQ(ierr); 201 PetscFunctionReturn(0); 202 } 203 204 /*------------------------------------------------------------*/ 205 static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer) 206 { 207 TAO_LMVM *lm = (TAO_LMVM *)tao->data; 208 PetscBool isascii; 209 PetscInt recycled_its; 210 PetscErrorCode ierr; 211 212 PetscFunctionBegin; 213 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 214 if (isascii) { 215 ierr = PetscViewerASCIIPrintf(viewer, " Gradient steps: %D\n", lm->grad);CHKERRQ(ierr); 216 if (lm->recycle) { 217 ierr = PetscViewerASCIIPrintf(viewer, " Recycle: on\n");CHKERRQ(ierr); 218 recycled_its = lm->bfgs + lm->grad; 219 ierr = PetscViewerASCIIPrintf(viewer, " Total recycled iterations: %D\n", recycled_its);CHKERRQ(ierr); 220 } 221 } 222 PetscFunctionReturn(0); 223 } 224 225 /* ---------------------------------------------------------- */ 226 227 /*MC 228 TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton 229 optimization solver for unconstrained minimization. It solves 230 the Newton step 231 Hkdk = - gk 232 233 using an approximation Bk in place of Hk, where Bk is composed using 234 the BFGS update formula. A More-Thuente line search is then used 235 to computed the steplength in the dk direction 236 Options Database Keys: 237 . -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls 238 239 Level: beginner 240 M*/ 241 242 PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao) 243 { 244 TAO_LMVM *lmP; 245 const char *morethuente_type = TAOLINESEARCHMT; 246 PetscErrorCode ierr; 247 248 PetscFunctionBegin; 249 tao->ops->setup = TaoSetUp_LMVM; 250 tao->ops->solve = TaoSolve_LMVM; 251 tao->ops->view = TaoView_LMVM; 252 tao->ops->setfromoptions = TaoSetFromOptions_LMVM; 253 tao->ops->destroy = TaoDestroy_LMVM; 254 255 ierr = PetscNewLog(tao,&lmP);CHKERRQ(ierr); 256 lmP->D = 0; 257 lmP->M = 0; 258 lmP->Xold = 0; 259 lmP->Gold = 0; 260 lmP->H0 = NULL; 261 lmP->recycle = PETSC_FALSE; 262 263 tao->data = (void*)lmP; 264 /* Override default settings (unless already changed) */ 265 if (!tao->max_it_changed) tao->max_it = 2000; 266 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 267 268 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr); 269 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr); 270 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr); 271 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 272 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr); 273 274 ierr = KSPInitializePackage();CHKERRQ(ierr); 275 ierr = MatCreate(((PetscObject)tao)->comm, &lmP->M);CHKERRQ(ierr); 276 ierr = PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1);CHKERRQ(ierr); 277 ierr = MatSetType(lmP->M, MATLMVMBFGS);CHKERRQ(ierr); 278 ierr = MatSetOptionsPrefix(lmP->M, "tao_lmvm_");CHKERRQ(ierr); 279 PetscFunctionReturn(0); 280 } 281 282 PETSC_EXTERN PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0) 283 { 284 TAO_LMVM *lmP; 285 TAO_BLMVM *blmP; 286 const TaoType type; 287 PetscBool is_lmvm, is_blmvm; 288 PetscErrorCode ierr; 289 290 ierr = TaoGetType(tao, &type);CHKERRQ(ierr); 291 ierr = PetscStrcmp(type, TAOLMVM, &is_lmvm);CHKERRQ(ierr); 292 ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr); 293 294 if (is_lmvm) { 295 lmP = (TAO_LMVM *)tao->data; 296 ierr = PetscObjectReference((PetscObject)H0);CHKERRQ(ierr); 297 lmP->H0 = H0; 298 } else if (is_blmvm) { 299 blmP = (TAO_BLMVM *)tao->data; 300 ierr = PetscObjectReference((PetscObject)H0);CHKERRQ(ierr); 301 blmP->H0 = H0; 302 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM."); 303 PetscFunctionReturn(0); 304 } 305 306 PETSC_EXTERN PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0) 307 { 308 TAO_LMVM *lmP; 309 TAO_BLMVM *blmP; 310 const TaoType type; 311 PetscBool is_lmvm, is_blmvm; 312 Mat M; 313 314 PetscErrorCode ierr; 315 316 ierr = TaoGetType(tao, &type);CHKERRQ(ierr); 317 ierr = PetscStrcmp(type, TAOLMVM, &is_lmvm);CHKERRQ(ierr); 318 ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr); 319 320 if (is_lmvm) { 321 lmP = (TAO_LMVM *)tao->data; 322 M = lmP->M; 323 } else if (is_blmvm) { 324 blmP = (TAO_BLMVM *)tao->data; 325 M = blmP->M; 326 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM."); 327 ierr = MatLMVMGetJ0(M, H0);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330 331 PETSC_EXTERN PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp) 332 { 333 TAO_LMVM *lmP; 334 TAO_BLMVM *blmP; 335 const TaoType type; 336 PetscBool is_lmvm, is_blmvm; 337 Mat M; 338 PetscErrorCode ierr; 339 340 ierr = TaoGetType(tao, &type);CHKERRQ(ierr); 341 ierr = PetscStrcmp(type, TAOLMVM, &is_lmvm);CHKERRQ(ierr); 342 ierr = PetscStrcmp(type, TAOBLMVM, &is_blmvm);CHKERRQ(ierr); 343 344 if (is_lmvm) { 345 lmP = (TAO_LMVM *)tao->data; 346 M = lmP->M; 347 } else if (is_blmvm) { 348 blmP = (TAO_BLMVM *)tao->data; 349 M = blmP->M; 350 } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM."); 351 ierr = MatLMVMGetJ0KSP(M, ksp);CHKERRQ(ierr); 352 PetscFunctionReturn(0); 353 }