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