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