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 #undef __FUNCT__ 10 #define __FUNCT__ "TaoSolve_LMVM" 11 static PetscErrorCode TaoSolve_LMVM(Tao tao) 12 { 13 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 14 PetscReal f, fold, gdx, gnorm; 15 PetscReal step = 1.0; 16 PetscReal delta; 17 PetscErrorCode ierr; 18 PetscInt stepType; 19 TaoConvergedReason reason = TAO_CONTINUE_ITERATING; 20 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 21 22 PetscFunctionBegin; 23 24 if (tao->XL || tao->XU || tao->ops->computebounds) { 25 ierr = PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");CHKERRQ(ierr); 26 } 27 28 /* Check convergence criteria */ 29 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr); 30 ierr = VecNorm(tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 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 lmP->bfgs = 0; 46 lmP->sgrad = 0; 47 lmP->grad = 0; 48 49 /* Have not converged; continue with Newton method */ 50 while (reason == TAO_CONTINUE_ITERATING) { 51 /* Compute direction */ 52 ierr = MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr); 53 ierr = MatLMVMSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr); 54 ++lmP->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 ++lmP->grad; 68 69 if (f != 0.0) { 70 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 71 } else { 72 delta = 2.0 / (gnorm*gnorm); 73 } 74 ierr = MatLMVMSetDelta(lmP->M, delta);CHKERRQ(ierr); 75 ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr); 76 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 77 ierr = MatLMVMSolve(lmP->M,tao->gradient, lmP->D);CHKERRQ(ierr); 78 79 /* On a reset, the direction cannot be not a number; it is a 80 scaled gradient step. No need to check for this condition. */ 81 82 lmP->bfgs = 1; 83 ++lmP->sgrad; 84 stepType = LMVM_SCALED_GRADIENT; 85 } else { 86 if (1 == lmP->bfgs) { 87 /* The first BFGS direction is always the scaled gradient */ 88 ++lmP->sgrad; 89 stepType = LMVM_SCALED_GRADIENT; 90 } else { 91 ++lmP->bfgs; 92 stepType = LMVM_BFGS; 93 } 94 } 95 ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr); 96 97 /* Perform the linesearch */ 98 fold = f; 99 ierr = VecCopy(tao->solution, lmP->Xold);CHKERRQ(ierr); 100 ierr = VecCopy(tao->gradient, lmP->Gold);CHKERRQ(ierr); 101 102 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);CHKERRQ(ierr); 103 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 104 105 while (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_GRADIENT)) { 106 /* Linesearch failed */ 107 /* Reset factors and use scaled gradient step */ 108 f = fold; 109 ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr); 110 ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr); 111 112 switch(stepType) { 113 case LMVM_BFGS: 114 /* Failed to obtain acceptable iterate with BFGS step */ 115 /* Attempt to use the scaled gradient direction */ 116 117 if (f != 0.0) { 118 delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm); 119 } else { 120 delta = 2.0 / (gnorm*gnorm); 121 } 122 ierr = MatLMVMSetDelta(lmP->M, delta);CHKERRQ(ierr); 123 ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr); 124 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 125 ierr = MatLMVMSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr); 126 127 /* On a reset, the direction cannot be not a number; it is a 128 scaled gradient step. No need to check for this condition. */ 129 130 lmP->bfgs = 1; 131 ++lmP->sgrad; 132 stepType = LMVM_SCALED_GRADIENT; 133 break; 134 135 case LMVM_SCALED_GRADIENT: 136 /* The scaled gradient step did not produce a new iterate; 137 attempt to use the gradient direction. 138 Need to make sure we are not using a different diagonal scaling */ 139 ierr = MatLMVMSetDelta(lmP->M, 1.0);CHKERRQ(ierr); 140 ierr = MatLMVMReset(lmP->M);CHKERRQ(ierr); 141 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 142 ierr = MatLMVMSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr); 143 144 lmP->bfgs = 1; 145 ++lmP->grad; 146 stepType = LMVM_GRADIENT; 147 break; 148 } 149 ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr); 150 151 /* Perform the linesearch */ 152 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);CHKERRQ(ierr); 153 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 154 } 155 156 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 157 /* Failed to find an improving point */ 158 f = fold; 159 ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr); 160 ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr); 161 step = 0.0; 162 reason = TAO_DIVERGED_LS_FAILURE; 163 tao->reason = TAO_DIVERGED_LS_FAILURE; 164 } 165 /* Check for termination */ 166 ierr = VecNorm(tao->gradient, NORM_2, &gnorm);CHKERRQ(ierr); 167 tao->niter++; 168 ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,step,&reason);CHKERRQ(ierr); 169 } 170 PetscFunctionReturn(0); 171 } 172 173 #undef __FUNCT__ 174 #define __FUNCT__ "TaoSetUp_LMVM" 175 static PetscErrorCode TaoSetUp_LMVM(Tao tao) 176 { 177 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 178 PetscInt n,N; 179 PetscErrorCode ierr; 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 PetscFunctionReturn(0); 195 } 196 197 /* ---------------------------------------------------------- */ 198 #undef __FUNCT__ 199 #define __FUNCT__ "TaoDestroy_LMVM" 200 static PetscErrorCode TaoDestroy_LMVM(Tao tao) 201 { 202 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 203 PetscErrorCode ierr; 204 205 PetscFunctionBegin; 206 if (tao->setupcalled) { 207 ierr = VecDestroy(&lmP->Xold);CHKERRQ(ierr); 208 ierr = VecDestroy(&lmP->Gold);CHKERRQ(ierr); 209 ierr = VecDestroy(&lmP->D);CHKERRQ(ierr); 210 ierr = MatDestroy(&lmP->M);CHKERRQ(ierr); 211 } 212 ierr = PetscFree(tao->data);CHKERRQ(ierr); 213 PetscFunctionReturn(0); 214 } 215 216 /*------------------------------------------------------------*/ 217 #undef __FUNCT__ 218 #define __FUNCT__ "TaoSetFromOptions_LMVM" 219 static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptions *PetscOptionsObject,Tao tao) 220 { 221 PetscErrorCode ierr; 222 223 PetscFunctionBegin; 224 ierr = PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");CHKERRQ(ierr); 225 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 226 ierr = PetscOptionsTail();CHKERRQ(ierr); 227 PetscFunctionReturn(0); 228 } 229 230 /*------------------------------------------------------------*/ 231 #undef __FUNCT__ 232 #define __FUNCT__ "TaoView_LMVM" 233 static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer) 234 { 235 TAO_LMVM *lm = (TAO_LMVM *)tao->data; 236 PetscBool isascii; 237 PetscErrorCode ierr; 238 239 PetscFunctionBegin; 240 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 241 if (isascii) { 242 ierr = PetscViewerASCIIPushTab(viewer);CHKERRQ(ierr); 243 ierr = PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);CHKERRQ(ierr); 244 ierr = PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);CHKERRQ(ierr); 245 ierr = PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);CHKERRQ(ierr); 246 ierr = PetscViewerASCIIPopTab(viewer);CHKERRQ(ierr); 247 } 248 PetscFunctionReturn(0); 249 } 250 251 /* ---------------------------------------------------------- */ 252 253 /*MC 254 TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton 255 optimization solver for unconstrained minimization. It solves 256 the Newton step 257 Hkdk = - gk 258 259 using an approximation Bk in place of Hk, where Bk is composed using 260 the BFGS update formula. A More-Thuente line search is then used 261 to computed the steplength in the dk direction 262 Options Database Keys: 263 + -tao_lmm_vectors - number of vectors to use for approximation 264 . -tao_lmm_scale_type - "none","scalar","broyden" 265 . -tao_lmm_limit_type - "none","average","relative","absolute" 266 . -tao_lmm_rescale_type - "none","scalar","gl" 267 . -tao_lmm_limit_mu - mu limiting factor 268 . -tao_lmm_limit_nu - nu limiting factor 269 . -tao_lmm_delta_min - minimum delta value 270 . -tao_lmm_delta_max - maximum delta value 271 . -tao_lmm_broyden_phi - phi factor for Broyden scaling 272 . -tao_lmm_scalar_alpha - alpha factor for scalar scaling 273 . -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal 274 . -tao_lmm_rescale_beta - beta factor for rescaling diagonal 275 . -tao_lmm_scalar_history - amount of history for scalar scaling 276 . -tao_lmm_rescale_history - amount of history for rescaling diagonal 277 - -tao_lmm_eps - rejection tolerance 278 279 Level: beginner 280 M*/ 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "TaoCreate_LMVM" 284 PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao) 285 { 286 TAO_LMVM *lmP; 287 const char *morethuente_type = TAOLINESEARCHMT; 288 PetscErrorCode ierr; 289 290 PetscFunctionBegin; 291 tao->ops->setup = TaoSetUp_LMVM; 292 tao->ops->solve = TaoSolve_LMVM; 293 tao->ops->view = TaoView_LMVM; 294 tao->ops->setfromoptions = TaoSetFromOptions_LMVM; 295 tao->ops->destroy = TaoDestroy_LMVM; 296 297 ierr = PetscNewLog(tao,&lmP);CHKERRQ(ierr); 298 lmP->D = 0; 299 lmP->M = 0; 300 lmP->Xold = 0; 301 lmP->Gold = 0; 302 303 tao->data = (void*)lmP; 304 /* Override default settings (unless already changed) */ 305 if (!tao->max_it_changed) tao->max_it = 2000; 306 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 307 if (!tao->fatol_changed) tao->fatol = 1.0e-4; 308 if (!tao->frtol_changed) tao->frtol = 1.0e-4; 309 310 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr); 311 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr); 312 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 313 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317