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