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