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;


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

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


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

 30:   /*  Set counter for gradient/reset steps */
 31:   if (!lmP->recycle) {
 32:     lmP->bfgs = 0;
 33:     lmP->grad = 0;
 34:     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:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);

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

 52:     /*  Check for success (descent direction) */
 53:     VecDot(lmP->D, tao->gradient, &gdx);
 54:     if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
 55:       /* Step is not descent or direction produced not a number
 56:          We can assert bfgsUpdates > 1 in this case because
 57:          the first solve produces the scaled gradient direction,
 58:          which is guaranteed to be descent

 60:          Use steepest descent direction (scaled)
 61:       */

 63:       MatLMVMReset(lmP->M, PETSC_FALSE);
 64:       MatLMVMClearJ0(lmP->M);
 65:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
 66:       MatSolve(lmP->M, tao->gradient, lmP->D);

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

 74:     /*  Perform the linesearch */
 75:     fold = f;
 76:     VecCopy(tao->solution, lmP->Xold);
 77:     VecCopy(tao->gradient, lmP->Gold);

 79:     TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);
 80:     TaoAddLineSearchCounts(tao);

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

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

 91:       MatLMVMReset(lmP->M, PETSC_FALSE);
 92:       MatLMVMClearJ0(lmP->M);
 93:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
 94:       MatSolve(lmP->M, tao->solution, tao->gradient);

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

101:       /*  Perform the linesearch */
102:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);
103:       TaoAddLineSearchCounts(tao);
104:     }

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

129:     /* Check convergence */
130:     tao->niter++;
131:     TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
132:     TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step);
133:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
134:   }
135:   return 0;
136: }

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

144:   /* Existence of tao->solution checked in TaoSetUp() */
145:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
146:   if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
147:   if (!lmP->D) VecDuplicate(tao->solution, &lmP->D);
148:   if (!lmP->Xold) VecDuplicate(tao->solution, &lmP->Xold);
149:   if (!lmP->Gold) VecDuplicate(tao->solution, &lmP->Gold);

151:   /*  Create matrix for the limited memory approximation */
152:   VecGetLocalSize(tao->solution, &n);
153:   VecGetSize(tao->solution, &N);
154:   MatSetSizes(lmP->M, n, n, N, N);
155:   MatLMVMAllocate(lmP->M, tao->solution, tao->gradient);
156:   MatIsSPDKnown(lmP->M, &is_set, &is_spd);

159:   /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
160:   if (lmP->H0) MatLMVMSetJ0(lmP->M, lmP->H0);
161:   return 0;
162: }

164: /* ---------------------------------------------------------- */
165: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
166: {
167:   TAO_LMVM *lmP = (TAO_LMVM *)tao->data;

169:   if (tao->setupcalled) {
170:     VecDestroy(&lmP->Xold);
171:     VecDestroy(&lmP->Gold);
172:     VecDestroy(&lmP->D);
173:   }
174:   MatDestroy(&lmP->M);
175:   if (lmP->H0) PetscObjectDereference((PetscObject)lmP->H0);
176:   PetscFree(tao->data);
177:   return 0;
178: }

180: /*------------------------------------------------------------*/
181: static PetscErrorCode TaoSetFromOptions_LMVM(Tao tao, PetscOptionItems *PetscOptionsObject)
182: {
183:   TAO_LMVM *lm = (TAO_LMVM *)tao->data;

185:   PetscOptionsHeadBegin(PetscOptionsObject, "Limited-memory variable-metric method for unconstrained optimization");
186:   PetscOptionsBool("-tao_lmvm_recycle", "enable recycling of the BFGS matrix between subsequent TaoSolve() calls", "", lm->recycle, &lm->recycle, NULL);
187:   TaoLineSearchSetFromOptions(tao->linesearch);
188:   MatSetFromOptions(lm->M);
189:   PetscOptionsHeadEnd();
190:   return 0;
191: }

193: /*------------------------------------------------------------*/
194: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
195: {
196:   TAO_LMVM *lm = (TAO_LMVM *)tao->data;
197:   PetscBool isascii;
198:   PetscInt  recycled_its;

200:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
201:   if (isascii) {
202:     PetscViewerASCIIPrintf(viewer, "  Gradient steps: %" PetscInt_FMT "\n", lm->grad);
203:     if (lm->recycle) {
204:       PetscViewerASCIIPrintf(viewer, "  Recycle: on\n");
205:       recycled_its = lm->bfgs + lm->grad;
206:       PetscViewerASCIIPrintf(viewer, "  Total recycled iterations: %" PetscInt_FMT "\n", recycled_its);
207:     }
208:   }
209:   return 0;
210: }

212: /* ---------------------------------------------------------- */

214: /*MC
215:   TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
216:   optimization solver for unconstrained minimization. It solves
217:   the Newton step
218:           Hkdk = - gk

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

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

228:   Level: beginner
229: M*/

231: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
232: {
233:   TAO_LMVM   *lmP;
234:   const char *morethuente_type = TAOLINESEARCHMT;

236:   tao->ops->setup          = TaoSetUp_LMVM;
237:   tao->ops->solve          = TaoSolve_LMVM;
238:   tao->ops->view           = TaoView_LMVM;
239:   tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
240:   tao->ops->destroy        = TaoDestroy_LMVM;

242:   PetscNew(&lmP);
243:   lmP->D       = NULL;
244:   lmP->M       = NULL;
245:   lmP->Xold    = NULL;
246:   lmP->Gold    = NULL;
247:   lmP->H0      = NULL;
248:   lmP->recycle = PETSC_FALSE;

250:   tao->data = (void *)lmP;
251:   /* Override default settings (unless already changed) */
252:   if (!tao->max_it_changed) tao->max_it = 2000;
253:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

255:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
256:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
257:   TaoLineSearchSetType(tao->linesearch, morethuente_type);
258:   TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
259:   TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);

261:   KSPInitializePackage();
262:   MatCreate(((PetscObject)tao)->comm, &lmP->M);
263:   PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1);
264:   MatSetType(lmP->M, MATLMVMBFGS);
265:   MatSetOptionsPrefix(lmP->M, "tao_lmvm_");
266:   return 0;
267: }