Actual source code: lmvm.c

petsc-3.14.6 2021-03-30
Report Typos and Errors
  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:   PetscErrorCode               ierr;
 13:   PetscInt                     stepType = LMVM_STEP_GRAD, nupdates;
 14:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;


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

 22:   /*  Check convergence criteria */
 23:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 24:   TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);

 26:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

 28:   tao->reason = TAO_CONTINUE_ITERATING;
 29:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
 30:   TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
 31:   (*tao->ops->convergencetest)(tao,tao->cnvP);
 32:   if (tao->reason != TAO_CONTINUE_ITERATING) return(0);

 34:   /*  Set counter for gradient/reset steps */
 35:   if (!lmP->recycle) {
 36:     lmP->bfgs = 0;
 37:     lmP->grad = 0;
 38:     MatLMVMReset(lmP->M, PETSC_FALSE);
 39:   }

 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:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
 46:     }

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

 58:     /*  Check for success (descent direction) */
 59:     VecDot(lmP->D, tao->gradient, &gdx);
 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

 66:          Use steepest descent direction (scaled)
 67:       */

 69:       MatLMVMReset(lmP->M, PETSC_FALSE);
 70:       MatLMVMClearJ0(lmP->M);
 71:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
 72:       MatSolve(lmP->M,tao->gradient, lmP->D);

 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:     VecScale(lmP->D, -1.0);

 80:     /*  Perform the linesearch */
 81:     fold = f;
 82:     VecCopy(tao->solution, lmP->Xold);
 83:     VecCopy(tao->gradient, lmP->Gold);

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

 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:       VecCopy(lmP->Xold, tao->solution);
 92:       VecCopy(lmP->Gold, tao->gradient);

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

 97:       MatLMVMReset(lmP->M, PETSC_FALSE);
 98:       MatLMVMClearJ0(lmP->M);
 99:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
100:       MatSolve(lmP->M, tao->solution, tao->gradient);

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:       VecScale(lmP->D, -1.0);

107:       /*  Perform the linesearch */
108:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);
109:       TaoAddLineSearchCounts(tao);
110:     }

112:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
113:       /*  Failed to find an improving point */
114:       f = fold;
115:       VecCopy(lmP->Xold, tao->solution);
116:       VecCopy(lmP->Gold, tao->gradient);
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:       TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
133:     }

135:     /* Check convergence */
136:     tao->niter++;
137:     TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
138:     TaoMonitor(tao,tao->niter,f,gnorm,0.0,step);
139:     (*tao->ops->convergencetest)(tao,tao->cnvP);
140:   }
141:   return(0);
142: }

144: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
145: {
146:   TAO_LMVM       *lmP = (TAO_LMVM *)tao->data;
147:   PetscInt       n,N;
149:   PetscBool      is_spd;

152:   /* Existence of tao->solution checked in TaoSetUp() */
153:   if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);  }
154:   if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection);  }
155:   if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D);  }
156:   if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold);  }
157:   if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold);  }

159:   /*  Create matrix for the limited memory approximation */
160:   VecGetLocalSize(tao->solution,&n);
161:   VecGetSize(tao->solution,&N);
162:   MatSetSizes(lmP->M, n, n, N, N);
163:   MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);
164:   MatGetOption(lmP->M, MAT_SPD, &is_spd);
165:   if (!is_spd) SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_INCOMP, "LMVM matrix is not symmetric positive-definite.");

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:     MatLMVMSetJ0(lmP->M, lmP->H0);
170:   }

172:   return(0);
173: }

175: /* ---------------------------------------------------------- */
176: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
177: {
178:   TAO_LMVM       *lmP = (TAO_LMVM *)tao->data;

182:   if (tao->setupcalled) {
183:     VecDestroy(&lmP->Xold);
184:     VecDestroy(&lmP->Gold);
185:     VecDestroy(&lmP->D);
186:   }
187:   MatDestroy(&lmP->M);
188:   if (lmP->H0) {
189:     PetscObjectDereference((PetscObject)lmP->H0);
190:   }
191:   PetscFree(tao->data);

193:   return(0);
194: }

196: /*------------------------------------------------------------*/
197: static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptionItems *PetscOptionsObject,Tao tao)
198: {
199:   TAO_LMVM       *lm = (TAO_LMVM *)tao->data;

203:   PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");
204:   PetscOptionsBool("-tao_lmvm_recycle","enable recycling of the BFGS matrix between subsequent TaoSolve() calls","",lm->recycle,&lm->recycle,NULL);
205:   TaoLineSearchSetFromOptions(tao->linesearch);
206:   MatSetFromOptions(lm->M);
207:   PetscOptionsTail();
208:   return(0);
209: }

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;

220:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
221:   if (isascii) {
222:     PetscViewerASCIIPrintf(viewer, "  Gradient steps: %D\n", lm->grad);
223:     if (lm->recycle) {
224:       PetscViewerASCIIPrintf(viewer, "  Recycle: on\n");
225:       recycled_its = lm->bfgs + lm->grad;
226:       PetscViewerASCIIPrintf(viewer, "  Total recycled iterations: %D\n", recycled_its);
227:     }
228:   }
229:   return(0);
230: }

232: /* ---------------------------------------------------------- */

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

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

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

248:   Level: beginner
249: M*/

251: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
252: {
253:   TAO_LMVM       *lmP;
254:   const char     *morethuente_type = TAOLINESEARCHMT;

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;

264:   PetscNewLog(tao,&lmP);
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;

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;

277:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
278:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
279:   TaoLineSearchSetType(tao->linesearch,morethuente_type);
280:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
281:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);

283:   KSPInitializePackage();
284:   MatCreate(((PetscObject)tao)->comm, &lmP->M);
285:   PetscObjectIncrementTabLevel((PetscObject)lmP->M, (PetscObject)tao, 1);
286:   MatSetType(lmP->M, MATLMVMBFGS);
287:   MatSetOptionsPrefix(lmP->M, "tao_lmvm_");
288:   return(0);
289: }