1 #include <petsctaolinesearch.h> 2 #include <../src/tao/unconstrained/impls/lmvm/lmvm.h> 3 4 #define LMVM_STEP_BFGS 0 5 #define LMVM_STEP_GRAD 1 6 7 static PetscErrorCode TaoSolve_LMVM(Tao tao) 8 { 9 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 10 PetscReal f, fold, gdx, gnorm; 11 PetscReal step = 1.0; 12 PetscErrorCode ierr; 13 PetscInt stepType = LMVM_STEP_GRAD, nupdates; 14 TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING; 15 16 PetscFunctionBegin; 17 18 if (tao->XL || tao->XU || tao->ops->computebounds) { 19 ierr = PetscInfo(tao,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");CHKERRQ(ierr); 20 } 21 22 /* Check convergence criteria */ 23 ierr = TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);CHKERRQ(ierr); 24 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 25 26 if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN"); 27 28 tao->reason = TAO_CONTINUE_ITERATING; 29 ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 30 ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);CHKERRQ(ierr); 31 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 32 if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(0); 33 34 /* Set counter for gradient/reset steps */ 35 if (!lmP->recycle) { 36 lmP->bfgs = 0; 37 lmP->grad = 0; 38 ierr = MatLMVMReset(lmP->M, PETSC_FALSE);CHKERRQ(ierr); 39 } 40 41 /* Have not converged; continue with Newton method */ 42 while (tao->reason == TAO_CONTINUE_ITERATING) { 43 /* Call general purpose update function */ 44 if (tao->ops->update) { 45 ierr = (*tao->ops->update)(tao, tao->niter, tao->user_update);CHKERRQ(ierr); 46 } 47 48 /* Compute direction */ 49 if (lmP->H0) { 50 ierr = MatLMVMSetJ0(lmP->M, lmP->H0);CHKERRQ(ierr); 51 stepType = LMVM_STEP_BFGS; 52 } 53 ierr = MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr); 54 ierr = MatSolve(lmP->M, tao->gradient, lmP->D);CHKERRQ(ierr); 55 ierr = MatLMVMGetUpdateCount(lmP->M, &nupdates);CHKERRQ(ierr); 56 if (nupdates > 0) stepType = LMVM_STEP_BFGS; 57 58 /* Check for success (descent direction) */ 59 ierr = VecDot(lmP->D, tao->gradient, &gdx);CHKERRQ(ierr); 60 if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) { 61 /* Step is not descent or direction produced not a number 62 We can assert bfgsUpdates > 1 in this case because 63 the first solve produces the scaled gradient direction, 64 which is guaranteed to be descent 65 66 Use steepest descent direction (scaled) 67 */ 68 69 ierr = MatLMVMReset(lmP->M, PETSC_FALSE);CHKERRQ(ierr); 70 ierr = MatLMVMClearJ0(lmP->M);CHKERRQ(ierr); 71 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 72 ierr = MatSolve(lmP->M,tao->gradient, lmP->D);CHKERRQ(ierr); 73 74 /* On a reset, the direction cannot be not a number; it is a 75 scaled gradient step. No need to check for this condition. */ 76 stepType = LMVM_STEP_GRAD; 77 } 78 ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr); 79 80 /* Perform the linesearch */ 81 fold = f; 82 ierr = VecCopy(tao->solution, lmP->Xold);CHKERRQ(ierr); 83 ierr = VecCopy(tao->gradient, lmP->Gold);CHKERRQ(ierr); 84 85 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step,&ls_status);CHKERRQ(ierr); 86 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 87 88 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) { 89 /* Reset factors and use scaled gradient step */ 90 f = fold; 91 ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr); 92 ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr); 93 94 /* Failed to obtain acceptable iterate with BFGS step */ 95 /* Attempt to use the scaled gradient direction */ 96 97 ierr = MatLMVMReset(lmP->M, PETSC_FALSE);CHKERRQ(ierr); 98 ierr = MatLMVMClearJ0(lmP->M);CHKERRQ(ierr); 99 ierr = MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 100 ierr = MatSolve(lmP->M, tao->solution, tao->gradient);CHKERRQ(ierr); 101 102 /* On a reset, the direction cannot be not a number; it is a 103 scaled gradient step. No need to check for this condition. */ 104 stepType = LMVM_STEP_GRAD; 105 ierr = VecScale(lmP->D, -1.0);CHKERRQ(ierr); 106 107 /* Perform the linesearch */ 108 ierr = TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);CHKERRQ(ierr); 109 ierr = TaoAddLineSearchCounts(tao);CHKERRQ(ierr); 110 } 111 112 if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) { 113 /* Failed to find an improving point */ 114 f = fold; 115 ierr = VecCopy(lmP->Xold, tao->solution);CHKERRQ(ierr); 116 ierr = VecCopy(lmP->Gold, tao->gradient);CHKERRQ(ierr); 117 step = 0.0; 118 tao->reason = TAO_DIVERGED_LS_FAILURE; 119 } else { 120 /* LS found valid step, so tally up step type */ 121 switch (stepType) { 122 case LMVM_STEP_BFGS: 123 ++lmP->bfgs; 124 break; 125 case LMVM_STEP_GRAD: 126 ++lmP->grad; 127 break; 128 default: 129 break; 130 } 131 /* Compute new gradient norm */ 132 ierr = TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);CHKERRQ(ierr); 133 } 134 135 /* Check convergence */ 136 tao->niter++; 137 ierr = TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);CHKERRQ(ierr); 138 ierr = TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);CHKERRQ(ierr); 139 ierr = (*tao->ops->convergencetest)(tao,tao->cnvP);CHKERRQ(ierr); 140 } 141 PetscFunctionReturn(0); 142 } 143 144 static PetscErrorCode TaoSetUp_LMVM(Tao tao) 145 { 146 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 147 PetscInt n,N; 148 PetscErrorCode ierr; 149 PetscBool is_spd; 150 151 PetscFunctionBegin; 152 /* Existence of tao->solution checked in TaoSetUp() */ 153 if (!tao->gradient) {ierr = VecDuplicate(tao->solution,&tao->gradient);CHKERRQ(ierr); } 154 if (!tao->stepdirection) {ierr = VecDuplicate(tao->solution,&tao->stepdirection);CHKERRQ(ierr); } 155 if (!lmP->D) {ierr = VecDuplicate(tao->solution,&lmP->D);CHKERRQ(ierr); } 156 if (!lmP->Xold) {ierr = VecDuplicate(tao->solution,&lmP->Xold);CHKERRQ(ierr); } 157 if (!lmP->Gold) {ierr = VecDuplicate(tao->solution,&lmP->Gold);CHKERRQ(ierr); } 158 159 /* Create matrix for the limited memory approximation */ 160 ierr = VecGetLocalSize(tao->solution,&n);CHKERRQ(ierr); 161 ierr = VecGetSize(tao->solution,&N);CHKERRQ(ierr); 162 ierr = MatSetSizes(lmP->M, n, n, N, N);CHKERRQ(ierr); 163 ierr = MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);CHKERRQ(ierr); 164 ierr = MatGetOption(lmP->M, MAT_SPD, &is_spd);CHKERRQ(ierr); 165 if (!is_spd) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite."); 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 } 171 172 PetscFunctionReturn(0); 173 } 174 175 /* ---------------------------------------------------------- */ 176 static PetscErrorCode TaoDestroy_LMVM(Tao tao) 177 { 178 TAO_LMVM *lmP = (TAO_LMVM *)tao->data; 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 if (tao->setupcalled) { 183 ierr = VecDestroy(&lmP->Xold);CHKERRQ(ierr); 184 ierr = VecDestroy(&lmP->Gold);CHKERRQ(ierr); 185 ierr = VecDestroy(&lmP->D);CHKERRQ(ierr); 186 } 187 ierr = MatDestroy(&lmP->M);CHKERRQ(ierr); 188 if (lmP->H0) { 189 ierr = PetscObjectDereference((PetscObject)lmP->H0);CHKERRQ(ierr); 190 } 191 ierr = PetscFree(tao->data);CHKERRQ(ierr); 192 193 PetscFunctionReturn(0); 194 } 195 196 /*------------------------------------------------------------*/ 197 static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject,Tao tao) 198 { 199 TAO_LMVM *lm = (TAO_LMVM *)tao->data; 200 PetscErrorCode ierr; 201 202 PetscFunctionBegin; 203 ierr = PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");CHKERRQ(ierr); 204 ierr = PetscOptionsBool("-tao_lmvm_recycle","enable recycling of the BFGS matrix between subsequent TaoSolve() calls","",lm->recycle,&lm->recycle,NULL);CHKERRQ(ierr); 205 ierr = TaoLineSearchSetFromOptions(tao->linesearch);CHKERRQ(ierr); 206 ierr = MatSetFromOptions(lm->M);CHKERRQ(ierr); 207 ierr = PetscOptionsTail();CHKERRQ(ierr); 208 PetscFunctionReturn(0); 209 } 210 211 /*------------------------------------------------------------*/ 212 static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer) 213 { 214 TAO_LMVM *lm = (TAO_LMVM *)tao->data; 215 PetscBool isascii; 216 PetscInt recycled_its; 217 PetscErrorCode ierr; 218 219 PetscFunctionBegin; 220 ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);CHKERRQ(ierr); 221 if (isascii) { 222 ierr = PetscViewerASCIIPrintf(viewer, " Gradient steps: %D\n", lm->grad);CHKERRQ(ierr); 223 if (lm->recycle) { 224 ierr = PetscViewerASCIIPrintf(viewer, " Recycle: on\n");CHKERRQ(ierr); 225 recycled_its = lm->bfgs + lm->grad; 226 ierr = PetscViewerASCIIPrintf(viewer, " Total recycled iterations: %D\n", recycled_its);CHKERRQ(ierr); 227 } 228 } 229 PetscFunctionReturn(0); 230 } 231 232 /* ---------------------------------------------------------- */ 233 234 /*MC 235 TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton 236 optimization solver for unconstrained minimization. It solves 237 the Newton step 238 Hkdk = - gk 239 240 using an approximation Bk in place of Hk, where Bk is composed using 241 the BFGS update formula. A More-Thuente line search is then used 242 to computed the steplength in the dk direction 243 244 Options Database Keys: 245 + -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls 246 - -tao_lmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation 247 248 Level: beginner 249 M*/ 250 251 PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao) 252 { 253 TAO_LMVM *lmP; 254 const char *morethuente_type = TAOLINESEARCHMT; 255 PetscErrorCode ierr; 256 257 PetscFunctionBegin; 258 tao->ops->setup = TaoSetUp_LMVM; 259 tao->ops->solve = TaoSolve_LMVM; 260 tao->ops->view = TaoView_LMVM; 261 tao->ops->setfromoptions = TaoSetFromOptions_LMVM; 262 tao->ops->destroy = TaoDestroy_LMVM; 263 264 ierr = PetscNewLog(tao,&lmP);CHKERRQ(ierr); 265 lmP->D = NULL; 266 lmP->M = NULL; 267 lmP->Xold = NULL; 268 lmP->Gold = NULL; 269 lmP->H0 = NULL; 270 lmP->recycle = PETSC_FALSE; 271 272 tao->data = (void*)lmP; 273 /* Override default settings (unless already changed) */ 274 if (!tao->max_it_changed) tao->max_it = 2000; 275 if (!tao->max_funcs_changed) tao->max_funcs = 4000; 276 277 ierr = TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);CHKERRQ(ierr); 278 ierr = PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);CHKERRQ(ierr); 279 ierr = TaoLineSearchSetType(tao->linesearch,morethuente_type);CHKERRQ(ierr); 280 ierr = TaoLineSearchUseTaoRoutines(tao->linesearch,tao);CHKERRQ(ierr); 281 ierr = TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);CHKERRQ(ierr); 282 283 ierr = KSPInitializePackage();CHKERRQ(ierr); 284 ierr = MatCreate(((PetscObject)tao)->comm, &lmP->M);CHKERRQ(ierr); 285 ierr = PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1);CHKERRQ(ierr); 286 ierr = MatSetType(lmP->M, MATLMVMBFGS);CHKERRQ(ierr); 287 ierr = MatSetOptionsPrefix(lmP->M, "tao_lmvm_");CHKERRQ(ierr); 288 PetscFunctionReturn(0); 289 } 290