Actual source code: lmvm.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <petsctaolinesearch.h>
  2: #include <../src/tao/matrix/lmvmmat.h>
  3: #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>

  5: #define LMVM_BFGS                0
  6: #define LMVM_SCALED_GRADIENT     1
  7: #define LMVM_GRADIENT            2

 11: static PetscErrorCode TaoSolve_LMVM(Tao tao)
 12: {
 13:   TAO_LMVM                     *lmP = (TAO_LMVM *)tao->data;
 14:   PetscReal                    f, fold, gdx, gnorm;
 15:   PetscReal                    step = 1.0;
 16:   PetscReal                    delta;
 17:   PetscErrorCode               ierr;
 18:   PetscInt                     stepType;
 19:   TaoConvergedReason           reason = TAO_CONTINUE_ITERATING;
 20:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;


 24:   if (tao->XL || tao->XU || tao->ops->computebounds) {
 25:     PetscPrintf(((PetscObject)tao)->comm,"WARNING: Variable bounds have been set but will be ignored by lmvm algorithm\n");
 26:   }

 28:   /*  Check convergence criteria */
 29:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 30:   VecNorm(tao->gradient,NORM_2,&gnorm);
 31:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");

 33:   TaoMonitor(tao, tao->niter, f, gnorm, 0.0, step, &reason);
 34:   if (reason != TAO_CONTINUE_ITERATING) return(0);

 36:   /*  Set initial scaling for the function */
 37:   if (f != 0.0) {
 38:     delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
 39:   } else {
 40:     delta = 2.0 / (gnorm*gnorm);
 41:   }
 42:   MatLMVMSetDelta(lmP->M,delta);

 44:   /*  Set counter for gradient/reset steps */
 45:   lmP->bfgs = 0;
 46:   lmP->sgrad = 0;
 47:   lmP->grad = 0;

 49:   /*  Have not converged; continue with Newton method */
 50:   while (reason == TAO_CONTINUE_ITERATING) {
 51:     /*  Compute direction */
 52:     MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
 53:     MatLMVMSolve(lmP->M, tao->gradient, lmP->D);
 54:     ++lmP->bfgs;

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

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

 67:       ++lmP->grad;

 69:       if (f != 0.0) {
 70:         delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
 71:       } else {
 72:         delta = 2.0 / (gnorm*gnorm);
 73:       }
 74:       MatLMVMSetDelta(lmP->M, delta);
 75:       MatLMVMReset(lmP->M);
 76:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
 77:       MatLMVMSolve(lmP->M,tao->gradient, lmP->D);

 79:       /* On a reset, the direction cannot be not a number; it is a
 80:          scaled gradient step.  No need to check for this condition. */

 82:       lmP->bfgs = 1;
 83:       ++lmP->sgrad;
 84:       stepType = LMVM_SCALED_GRADIENT;
 85:     } else {
 86:       if (1 == lmP->bfgs) {
 87:         /*  The first BFGS direction is always the scaled gradient */
 88:         ++lmP->sgrad;
 89:         stepType = LMVM_SCALED_GRADIENT;
 90:       } else {
 91:         ++lmP->bfgs;
 92:         stepType = LMVM_BFGS;
 93:       }
 94:     }
 95:     VecScale(lmP->D, -1.0);

 97:     /*  Perform the linesearch */
 98:     fold = f;
 99:     VecCopy(tao->solution, lmP->Xold);
100:     VecCopy(tao->gradient, lmP->Gold);

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

105:     while (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER && (stepType != LMVM_GRADIENT)) {
106:       /*  Linesearch failed */
107:       /*  Reset factors and use scaled gradient step */
108:       f = fold;
109:       VecCopy(lmP->Xold, tao->solution);
110:       VecCopy(lmP->Gold, tao->gradient);

112:       switch(stepType) {
113:       case LMVM_BFGS:
114:         /*  Failed to obtain acceptable iterate with BFGS step */
115:         /*  Attempt to use the scaled gradient direction */

117:         if (f != 0.0) {
118:           delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
119:         } else {
120:           delta = 2.0 / (gnorm*gnorm);
121:         }
122:         MatLMVMSetDelta(lmP->M, delta);
123:         MatLMVMReset(lmP->M);
124:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
125:         MatLMVMSolve(lmP->M, tao->gradient, lmP->D);

127:         /* On a reset, the direction cannot be not a number; it is a
128:            scaled gradient step.  No need to check for this condition. */

130:         lmP->bfgs = 1;
131:         ++lmP->sgrad;
132:         stepType = LMVM_SCALED_GRADIENT;
133:         break;

135:       case LMVM_SCALED_GRADIENT:
136:         /* The scaled gradient step did not produce a new iterate;
137:            attempt to use the gradient direction.
138:            Need to make sure we are not using a different diagonal scaling */
139:         MatLMVMSetDelta(lmP->M, 1.0);
140:         MatLMVMReset(lmP->M);
141:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
142:         MatLMVMSolve(lmP->M, tao->gradient, lmP->D);

144:         lmP->bfgs = 1;
145:         ++lmP->grad;
146:         stepType = LMVM_GRADIENT;
147:         break;
148:       }
149:       VecScale(lmP->D, -1.0);

151:       /*  Perform the linesearch */
152:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, tao->gradient, lmP->D, &step, &ls_status);
153:       TaoAddLineSearchCounts(tao);
154:     }

156:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
157:       /*  Failed to find an improving point */
158:       f = fold;
159:       VecCopy(lmP->Xold, tao->solution);
160:       VecCopy(lmP->Gold, tao->gradient);
161:       step = 0.0;
162:       reason = TAO_DIVERGED_LS_FAILURE;
163:       tao->reason = TAO_DIVERGED_LS_FAILURE;
164:     }
165:     /*  Check for termination */
166:     VecNorm(tao->gradient, NORM_2, &gnorm);
167:     tao->niter++;
168:     TaoMonitor(tao,tao->niter,f,gnorm,0.0,step,&reason);
169:   }
170:   return(0);
171: }

175: static PetscErrorCode TaoSetUp_LMVM(Tao tao)
176: {
177:   TAO_LMVM       *lmP = (TAO_LMVM *)tao->data;
178:   PetscInt       n,N;

182:   /* Existence of tao->solution checked in TaoSetUp() */
183:   if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);  }
184:   if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection);  }
185:   if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D);  }
186:   if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold);  }
187:   if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold);  }

189:   /*  Create matrix for the limited memory approximation */
190:   VecGetLocalSize(tao->solution,&n);
191:   VecGetSize(tao->solution,&N);
192:   MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M);
193:   MatLMVMAllocateVectors(lmP->M,tao->solution);
194:   return(0);
195: }

197: /* ---------------------------------------------------------- */
200: static PetscErrorCode TaoDestroy_LMVM(Tao tao)
201: {
202:   TAO_LMVM       *lmP = (TAO_LMVM *)tao->data;

206:   if (tao->setupcalled) {
207:     VecDestroy(&lmP->Xold);
208:     VecDestroy(&lmP->Gold);
209:     VecDestroy(&lmP->D);
210:     MatDestroy(&lmP->M);
211:   }
212:   PetscFree(tao->data);
213:   return(0);
214: }

216: /*------------------------------------------------------------*/
219: static PetscErrorCode TaoSetFromOptions_LMVM(PetscOptions *PetscOptionsObject,Tao tao)
220: {

224:   PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for unconstrained optimization");
225:   TaoLineSearchSetFromOptions(tao->linesearch);
226:   PetscOptionsTail();
227:   return(0);
228: }

230: /*------------------------------------------------------------*/
233: static PetscErrorCode TaoView_LMVM(Tao tao, PetscViewer viewer)
234: {
235:   TAO_LMVM       *lm = (TAO_LMVM *)tao->data;
236:   PetscBool      isascii;

240:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
241:   if (isascii) {
242:     PetscViewerASCIIPushTab(viewer);
243:     PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
244:     PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
245:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
246:     PetscViewerASCIIPopTab(viewer);
247:   }
248:   return(0);
249: }

251: /* ---------------------------------------------------------- */

253: /*MC
254:      TAOLMVM - Limited Memory Variable Metric method is a quasi-Newton
255:      optimization solver for unconstrained minimization. It solves
256:      the Newton step
257:               Hkdk = - gk

259:      using an approximation Bk in place of Hk, where Bk is composed using
260:      the BFGS update formula. A More-Thuente line search is then used
261:      to computed the steplength in the dk direction
262:   Options Database Keys:
263: +     -tao_lmm_vectors - number of vectors to use for approximation
264: .     -tao_lmm_scale_type - "none","scalar","broyden"
265: .     -tao_lmm_limit_type - "none","average","relative","absolute"
266: .     -tao_lmm_rescale_type - "none","scalar","gl"
267: .     -tao_lmm_limit_mu - mu limiting factor
268: .     -tao_lmm_limit_nu - nu limiting factor
269: .     -tao_lmm_delta_min - minimum delta value
270: .     -tao_lmm_delta_max - maximum delta value
271: .     -tao_lmm_broyden_phi - phi factor for Broyden scaling
272: .     -tao_lmm_scalar_alpha - alpha factor for scalar scaling
273: .     -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal
274: .     -tao_lmm_rescale_beta - beta factor for rescaling diagonal
275: .     -tao_lmm_scalar_history - amount of history for scalar scaling
276: .     -tao_lmm_rescale_history - amount of history for rescaling diagonal
277: -     -tao_lmm_eps - rejection tolerance

279:   Level: beginner
280: M*/

284: PETSC_EXTERN PetscErrorCode TaoCreate_LMVM(Tao tao)
285: {
286:   TAO_LMVM       *lmP;
287:   const char     *morethuente_type = TAOLINESEARCHMT;

291:   tao->ops->setup = TaoSetUp_LMVM;
292:   tao->ops->solve = TaoSolve_LMVM;
293:   tao->ops->view = TaoView_LMVM;
294:   tao->ops->setfromoptions = TaoSetFromOptions_LMVM;
295:   tao->ops->destroy = TaoDestroy_LMVM;

297:   PetscNewLog(tao,&lmP);
298:   lmP->D = 0;
299:   lmP->M = 0;
300:   lmP->Xold = 0;
301:   lmP->Gold = 0;

303:   tao->data = (void*)lmP;
304:   /* Override default settings (unless already changed) */
305:   if (!tao->max_it_changed) tao->max_it = 2000;
306:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
307:   if (!tao->fatol_changed) tao->fatol = 1.0e-4;
308:   if (!tao->frtol_changed) tao->frtol = 1.0e-4;

310:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
311:   TaoLineSearchSetType(tao->linesearch,morethuente_type);
312:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
313:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
314:   return(0);
315: }