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