Actual source code: lmvm.c

  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>

  4: #define LMVM_STEP_BFGS 0
  5: #define LMVM_STEP_GRAD 1

  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:   PetscInt                     stepType  = LMVM_STEP_GRAD, nupdates;
 13:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

 15:   PetscFunctionBegin;
 16:   if (tao->XL || tao->XU || tao->ops->computebounds) PetscCall(PetscInfo(tao, "WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n"));

 18:   /*  Check convergence criteria */
 19:   PetscCall(TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient));
 20:   PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));

 22:   PetscCheck(!PetscIsInfOrNanReal(f) && !PetscIsInfOrNanReal(gnorm), PetscObjectComm((PetscObject)tao), PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

 24:   tao->reason = TAO_CONTINUE_ITERATING;
 25:   PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
 26:   PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
 27:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 28:   if (tao->reason != TAO_CONTINUE_ITERATING) PetscFunctionReturn(PETSC_SUCCESS);

 30:   /*  Set counter for gradient/reset steps */
 31:   if (!lmP->recycle) {
 32:     lmP->bfgs = 0;
 33:     lmP->grad = 0;
 34:     PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
 35:   }

 37:   /*  Have not converged; continue with Newton method */
 38:   while (tao->reason == TAO_CONTINUE_ITERATING) {
 39:     /* Call general purpose update function */
 40:     if (tao->ops->update) {
 41:       PetscUseTypeMethod(tao, update, tao->niter, tao->user_update);
 42:       PetscCall(TaoComputeObjective(tao, tao->solution, &f));
 43:     }

 45:     /*  Compute direction */
 46:     if (lmP->H0) {
 47:       PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0));
 48:       stepType = LMVM_STEP_BFGS;
 49:     }
 50:     PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
 51:     PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D));
 52:     PetscCall(MatLMVMGetUpdateCount(lmP->M, &nupdates));
 53:     if (nupdates > 0) stepType = LMVM_STEP_BFGS;

 55:     /*  Check for success (descent direction) */
 56:     PetscCall(VecDotRealPart(lmP->D, tao->gradient, &gdx));
 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

 63:          Use steepest descent direction (scaled)
 64:       */

 66:       PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
 67:       PetscCall(MatLMVMClearJ0(lmP->M));
 68:       PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
 69:       PetscCall(MatSolve(lmP->M, tao->gradient, lmP->D));

 71:       /* On a reset, the direction cannot be not a number; it is a
 72:          scaled gradient step.  No need to check for this condition. */
 73:       stepType = LMVM_STEP_GRAD;
 74:     }
 75:     PetscCall(VecScale(lmP->D, -1.0));

 77:     /*  Perform the linesearch */
 78:     fold = f;
 79:     PetscCall(VecCopy(tao->solution, lmP->Xold));
 80:     PetscCall(VecCopy(tao->gradient, lmP->Gold));

 82:     PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status));
 83:     PetscCall(TaoAddLineSearchCounts(tao));

 85:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_STEP_GRAD)) {
 86:       /*  Reset factors and use scaled gradient step */
 87:       f = fold;
 88:       PetscCall(VecCopy(lmP->Xold, tao->solution));
 89:       PetscCall(VecCopy(lmP->Gold, tao->gradient));

 91:       /*  Failed to obtain acceptable iterate with BFGS step */
 92:       /*  Attempt to use the scaled gradient direction */

 94:       PetscCall(MatLMVMReset(lmP->M, PETSC_FALSE));
 95:       PetscCall(MatLMVMClearJ0(lmP->M));
 96:       PetscCall(MatLMVMUpdate(lmP->M, tao->solution, tao->gradient));
 97:       PetscCall(MatSolve(lmP->M, tao->solution, tao->gradient));

 99:       /* On a reset, the direction cannot be not a number; it is a
100:           scaled gradient step.  No need to check for this condition. */
101:       stepType = LMVM_STEP_GRAD;
102:       PetscCall(VecScale(lmP->D, -1.0));

104:       /*  Perform the linesearch */
105:       PetscCall(TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status));
106:       PetscCall(TaoAddLineSearchCounts(tao));
107:     }

109:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
110:       /*  Failed to find an improving point */
111:       f = fold;
112:       PetscCall(VecCopy(lmP->Xold, tao->solution));
113:       PetscCall(VecCopy(lmP->Gold, tao->gradient));
114:       step        = 0.0;
115:       tao->reason = TAO_DIVERGED_LS_FAILURE;
116:     } else {
117:       /* LS found valid step, so tally up step type */
118:       switch (stepType) {
119:       case LMVM_STEP_BFGS:
120:         ++lmP->bfgs;
121:         break;
122:       case LMVM_STEP_GRAD:
123:         ++lmP->grad;
124:         break;
125:       default:
126:         break;
127:       }
128:       /*  Compute new gradient norm */
129:       PetscCall(TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm));
130:     }

132:     /* Check convergence */
133:     tao->niter++;
134:     PetscCall(TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its));
135:     PetscCall(TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step));
136:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
137:   }
138:   PetscFunctionReturn(PETSC_SUCCESS);
139: }

141: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
142: {
143:   TAO_LMVM *lmP = (TAO_LMVM *)tao->data;
144:   PetscInt  n, N;
145:   PetscBool is_set, is_spd;

147:   PetscFunctionBegin;
148:   /* Existence of tao->solution checked in TaoSetUp() */
149:   if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
150:   if (!tao->stepdirection) PetscCall(VecDuplicate(tao->solution, &tao->stepdirection));
151:   if (!lmP->D) PetscCall(VecDuplicate(tao->solution, &lmP->D));
152:   if (!lmP->Xold) PetscCall(VecDuplicate(tao->solution, &lmP->Xold));
153:   if (!lmP->Gold) PetscCall(VecDuplicate(tao->solution, &lmP->Gold));

155:   /*  Create matrix for the limited memory approximation */
156:   PetscCall(VecGetLocalSize(tao->solution, &n));
157:   PetscCall(VecGetSize(tao->solution, &N));
158:   PetscCall(MatSetSizes(lmP->M, n, n, N, N));
159:   PetscCall(MatLMVMAllocate(lmP->M, tao->solution, tao->gradient));
160:   PetscCall(MatIsSPDKnown(lmP->M, &is_set, &is_spd));
161:   PetscCheck(is_set && is_spd, PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite.");

163:   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
164:   if (lmP->H0) PetscCall(MatLMVMSetJ0(lmP->M, lmP->H0));
165:   PetscFunctionReturn(PETSC_SUCCESS);
166: }

168: /* ---------------------------------------------------------- */
169: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
170: {
171:   TAO_LMVM *lmP = (TAO_LMVM *)tao->data;

173:   PetscFunctionBegin;
174:   if (tao->setupcalled) {
175:     PetscCall(VecDestroy(&lmP->Xold));
176:     PetscCall(VecDestroy(&lmP->Gold));
177:     PetscCall(VecDestroy(&lmP->D));
178:   }
179:   PetscCall(MatDestroy(&lmP->M));
180:   if (lmP->H0) PetscCall(PetscObjectDereference((PetscObject)lmP->H0));
181:   PetscCall(PetscFree(tao->data));
182:   PetscFunctionReturn(PETSC_SUCCESS);
183: }

185: /*------------------------------------------------------------*/
186: static PetscErrorCode TaoSetFromOptions_LMVM(Tao tao, PetscOptionItems *PetscOptionsObject)
187: {
188:   TAO_LMVM *lm = (TAO_LMVM *)tao->data;

190:   PetscFunctionBegin;
191:   PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for unconstrained optimization");
192:   PetscCall(PetscOptionsBool("-tao_lmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", lm->recycle, &lm->recycle, NULL));
193:   PetscCall(TaoLineSearchSetFromOptions(tao->linesearch));
194:   PetscCall(MatSetFromOptions(lm->M));
195:   PetscOptionsHeadEnd();
196:   PetscFunctionReturn(PETSC_SUCCESS);
197: }

199: /*------------------------------------------------------------*/
200: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
201: {
202:   TAO_LMVM *lm = (TAO_LMVM *)tao->data;
203:   PetscBool isascii;
204:   PetscInt  recycled_its;

206:   PetscFunctionBegin;
207:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
208:   if (isascii) {
209:     PetscCall(PetscViewerASCIIPushTab(viewer));
210:     PetscCall(PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lm->grad));
211:     if (lm->recycle) {
212:       PetscCall(PetscViewerASCIIPrintf(viewer, "Recycle: on\n"));
213:       recycled_its = lm->bfgs + lm->grad;
214:       PetscCall(PetscViewerASCIIPrintf(viewer, "Total recycled iterations: %" PetscInt_FMT "\n", recycled_its));
215:     }
216:     PetscCall(PetscViewerASCIIPrintf(viewer, "LMVM Matrix:\n"));
217:     PetscCall(PetscViewerASCIIPushTab(viewer));
218:     PetscCall(MatView(lm->M, viewer));
219:     PetscCall(PetscViewerASCIIPopTab(viewer));
220:     PetscCall(PetscViewerASCIIPopTab(viewer));
221:   }
222:   PetscFunctionReturn(PETSC_SUCCESS);
223: }

225: /* ---------------------------------------------------------- */

227: /*MC
228:   TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
229:   optimization solver for unconstrained minimization. It solves
230:   the Newton step
231:           Hkdk = - gk

233:   using an approximation Bk in place of Hk, where Bk is composed using
234:   the BFGS update formula. A More-Thuente line search is then used
235:   to computed the steplength in the dk direction

237:   Options Database Keys:
238: +   -tao_lmvm_recycle - enable recycling LMVM updates between TaoSolve() calls
239: -   -tao_lmvm_no_scale - (developer) disables diagonal Broyden scaling on the LMVM approximation

241:   Level: beginner
242: M*/

244: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
245: {
246:   TAO_LMVM   *lmP;
247:   const char *morethuente_type = TAOLINESEARCHMT;

249:   PetscFunctionBegin;
250:   tao->ops->setup          = TaoSetUp_LMVM;
251:   tao->ops->solve          = TaoSolve_LMVM;
252:   tao->ops->view           = TaoView_LMVM;
253:   tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
254:   tao->ops->destroy        = TaoDestroy_LMVM;

256:   PetscCall(PetscNew(&lmP));
257:   lmP->D       = NULL;
258:   lmP->M       = NULL;
259:   lmP->Xold    = NULL;
260:   lmP->Gold    = NULL;
261:   lmP->H0      = NULL;
262:   lmP->recycle = PETSC_FALSE;

264:   tao->data = (void *)lmP;
265:   /* Override default settings (unless already changed) */
266:   if (!tao->max_it_changed) tao->max_it = 2000;
267:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

269:   PetscCall(TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch));
270:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1));
271:   PetscCall(TaoLineSearchSetType(tao->linesearch, morethuente_type));
272:   PetscCall(TaoLineSearchUseTaoRoutines(tao->linesearch, tao));
273:   PetscCall(TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix));

275:   PetscCall(KSPInitializePackage());
276:   PetscCall(MatCreate(((PetscObject)tao)->comm, &lmP->M));
277:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1));
278:   PetscCall(MatSetType(lmP->M, MATLMVMBFGS));
279:   PetscCall(MatSetOptionsPrefix(lmP->M, "tao_lmvm_"));
280:   PetscFunctionReturn(PETSC_SUCCESS);
281: }