Actual source code: owlqn.c

petsc-3.12.5 2020-03-29
Report Typos and Errors
  1:  #include <petsctaolinesearch.h>
  2:  #include <../src/tao/unconstrained/impls/owlqn/owlqn.h>

  4: #define OWLQN_BFGS                0
  5: #define OWLQN_SCALED_GRADIENT     1
  6: #define OWLQN_GRADIENT            2

  8: static PetscErrorCode ProjDirect_OWLQN(Vec d, Vec g)
  9: {
 10:   PetscErrorCode  ierr;
 11:   const PetscReal *gptr;
 12:   PetscReal       *dptr;
 13:   PetscInt        low,high,low1,high1,i;

 16:   ierr=VecGetOwnershipRange(d,&low,&high);
 17:   ierr=VecGetOwnershipRange(g,&low1,&high1);

 19:   VecGetArrayRead(g,&gptr);
 20:   VecGetArray(d,&dptr);
 21:   for (i = 0; i < high-low; i++) {
 22:     if (dptr[i] * gptr[i] <= 0.0 ) {
 23:       dptr[i] = 0.0;
 24:     }
 25:   }
 26:   VecRestoreArray(d,&dptr);
 27:   VecRestoreArrayRead(g,&gptr);
 28:   return(0);
 29: }

 31: static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda)
 32: {
 33:   PetscErrorCode  ierr;
 34:   const PetscReal *xptr;
 35:   PetscReal       *gptr;
 36:   PetscInt        low,high,low1,high1,i;

 39:   ierr=VecGetOwnershipRange(x,&low,&high);
 40:   ierr=VecGetOwnershipRange(gv,&low1,&high1);

 42:   VecGetArrayRead(x,&xptr);
 43:   VecGetArray(gv,&gptr);
 44:   for (i = 0; i < high-low; i++) {
 45:     if (xptr[i] < 0.0)               gptr[i] = gptr[i] - lambda;
 46:     else if (xptr[i] > 0.0)          gptr[i] = gptr[i] + lambda;
 47:     else if (gptr[i] + lambda < 0.0) gptr[i] = gptr[i] + lambda;
 48:     else if (gptr[i] - lambda > 0.0) gptr[i] = gptr[i] - lambda;
 49:     else                             gptr[i] = 0.0;
 50:   }
 51:   VecRestoreArray(gv,&gptr);
 52:   VecRestoreArrayRead(x,&xptr);
 53:   return(0);
 54: }

 56: static PetscErrorCode TaoSolve_OWLQN(Tao tao)
 57: {
 58:   TAO_OWLQN                    *lmP = (TAO_OWLQN *)tao->data;
 59:   PetscReal                    f, fold, gdx, gnorm;
 60:   PetscReal                    step = 1.0;
 61:   PetscReal                    delta;
 62:   PetscErrorCode               ierr;
 63:   PetscInt                     stepType;
 64:   PetscInt                     iter = 0;
 65:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

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

 72:   /* Check convergence criteria */
 73:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);

 75:   VecCopy(tao->gradient, lmP->GV);

 77:   ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);

 79:   VecNorm(lmP->GV,NORM_2,&gnorm);

 81:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf or NaN");

 83:   tao->reason = TAO_CONTINUE_ITERATING;
 84:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
 85:   TaoMonitor(tao,iter,f,gnorm,0.0,step);
 86:   (*tao->ops->convergencetest)(tao,tao->cnvP);
 87:   if (tao->reason != TAO_CONTINUE_ITERATING) return(0);

 89:   /* Set initial scaling for the function */
 90:   delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
 91:   MatLMVMSetJ0Scale(lmP->M, delta);

 93:   /* Set counter for gradient/reset steps */
 94:   lmP->bfgs = 0;
 95:   lmP->sgrad = 0;
 96:   lmP->grad = 0;

 98:   /* Have not converged; continue with Newton method */
 99:   while (tao->reason == TAO_CONTINUE_ITERATING) {
100:     /* Call general purpose update function */
101:     if (tao->ops->update) {
102:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
103:     }

105:     /* Compute direction */
106:     MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
107:     MatSolve(lmP->M, lmP->GV, lmP->D);

109:     ProjDirect_OWLQN(lmP->D,lmP->GV);

111:     ++lmP->bfgs;

113:     /* Check for success (descent direction) */
114:     VecDot(lmP->D, lmP->GV , &gdx);
115:     if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {

117:       /* Step is not descent or direction produced not a number
118:          We can assert bfgsUpdates > 1 in this case because
119:          the first solve produces the scaled gradient direction,
120:          which is guaranteed to be descent

122:          Use steepest descent direction (scaled) */
123:       ++lmP->grad;

125:       delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
126:       MatLMVMSetJ0Scale(lmP->M, delta);
127:       MatLMVMReset(lmP->M, PETSC_FALSE);
128:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
129:       MatSolve(lmP->M,lmP->GV, lmP->D);

131:       ProjDirect_OWLQN(lmP->D,lmP->GV);

133:       lmP->bfgs = 1;
134:       ++lmP->sgrad;
135:       stepType = OWLQN_SCALED_GRADIENT;
136:     } else {
137:       if (1 == lmP->bfgs) {
138:         /* The first BFGS direction is always the scaled gradient */
139:         ++lmP->sgrad;
140:         stepType = OWLQN_SCALED_GRADIENT;
141:       } else {
142:         ++lmP->bfgs;
143:         stepType = OWLQN_BFGS;
144:       }
145:     }

147:     VecScale(lmP->D, -1.0);

149:     /* Perform the linesearch */
150:     fold = f;
151:     VecCopy(tao->solution, lmP->Xold);
152:     VecCopy(tao->gradient, lmP->Gold);

154:     TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step,&ls_status);
155:     TaoAddLineSearchCounts(tao);

157:     while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) {

159:       /* Reset factors and use scaled gradient step */
160:       f = fold;
161:       VecCopy(lmP->Xold, tao->solution);
162:       VecCopy(lmP->Gold, tao->gradient);
163:       VecCopy(tao->gradient, lmP->GV);

165:       ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);

167:       switch(stepType) {
168:       case OWLQN_BFGS:
169:         /* Failed to obtain acceptable iterate with BFGS step
170:            Attempt to use the scaled gradient direction */

172:         delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
173:         MatLMVMSetJ0Scale(lmP->M, delta);
174:         MatLMVMReset(lmP->M, PETSC_FALSE);
175:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
176:         MatSolve(lmP->M, lmP->GV, lmP->D);

178:         ProjDirect_OWLQN(lmP->D,lmP->GV);

180:         lmP->bfgs = 1;
181:         ++lmP->sgrad;
182:         stepType = OWLQN_SCALED_GRADIENT;
183:         break;

185:       case OWLQN_SCALED_GRADIENT:
186:         /* The scaled gradient step did not produce a new iterate;
187:            attempt to use the gradient direction.
188:            Need to make sure we are not using a different diagonal scaling */
189:         MatLMVMSetJ0Scale(lmP->M, 1.0);
190:         MatLMVMReset(lmP->M, PETSC_FALSE);
191:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
192:         MatSolve(lmP->M, lmP->GV, lmP->D);

194:         ProjDirect_OWLQN(lmP->D,lmP->GV);

196:         lmP->bfgs = 1;
197:         ++lmP->grad;
198:         stepType = OWLQN_GRADIENT;
199:         break;
200:       }
201:       VecScale(lmP->D, -1.0);

203:       /* Perform the linesearch */
204:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
205:       TaoAddLineSearchCounts(tao);
206:     }

208:     if ((int)ls_status < 0) {
209:       /* Failed to find an improving point*/
210:       f = fold;
211:       VecCopy(lmP->Xold, tao->solution);
212:       VecCopy(lmP->Gold, tao->gradient);
213:       VecCopy(tao->gradient, lmP->GV);
214:       step = 0.0;
215:     } else {
216:       /* a little hack here, because that gv is used to store g */
217:       VecCopy(lmP->GV, tao->gradient);
218:     }

220:     ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);

222:     /* Check for termination */

224:     VecNorm(lmP->GV,NORM_2,&gnorm);

226:     iter++;
227:     TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
228:     TaoMonitor(tao,iter,f,gnorm,0.0,step);
229:     (*tao->ops->convergencetest)(tao,tao->cnvP);

231:     if ((int)ls_status < 0) break;
232:   }
233:   return(0);
234: }

236: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
237: {
238:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;
239:   PetscInt       n,N;

243:   /* Existence of tao->solution checked in TaoSetUp() */
244:   if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);  }
245:   if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection);  }
246:   if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D);  }
247:   if (!lmP->GV) {VecDuplicate(tao->solution,&lmP->GV);  }
248:   if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold);  }
249:   if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold);  }

251:   /* Create matrix for the limited memory approximation */
252:   VecGetLocalSize(tao->solution,&n);
253:   VecGetSize(tao->solution,&N);
254:   MatCreateLMVMBFGS(((PetscObject)tao)->comm,n,N,&lmP->M);
255:   MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);
256:   return(0);
257: }

259: /* ---------------------------------------------------------- */
260: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
261: {
262:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;

266:   if (tao->setupcalled) {
267:     VecDestroy(&lmP->Xold);
268:     VecDestroy(&lmP->Gold);
269:     VecDestroy(&lmP->D);
270:     MatDestroy(&lmP->M);
271:     VecDestroy(&lmP->GV);
272:   }
273:   PetscFree(tao->data);
274:   return(0);
275: }

277: /*------------------------------------------------------------*/
278: static PetscErrorCode TaoSetFromOptions_OWLQN(PetscOptionItems *PetscOptionsObject,Tao tao)
279: {
280:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;

284:   PetscOptionsHead(PetscOptionsObject,"Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
285:   PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,NULL);
286:   PetscOptionsTail();
287:   TaoLineSearchSetFromOptions(tao->linesearch);
288:   return(0);
289: }

291: /*------------------------------------------------------------*/
292: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
293: {
294:   TAO_OWLQN      *lm = (TAO_OWLQN *)tao->data;
295:   PetscBool      isascii;

299:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
300:   if (isascii) {
301:     PetscViewerASCIIPushTab(viewer);
302:     PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
303:     PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
304:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
305:     PetscViewerASCIIPopTab(viewer);
306:   }
307:   return(0);
308: }

310: /* ---------------------------------------------------------- */
311: /*MC
312:   TAOOWLQN - orthant-wise limited memory quasi-newton algorithm

314: . - tao_owlqn_lambda - regulariser weight

316:   Level: beginner
317: M*/


320: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
321: {
322:   TAO_OWLQN      *lmP;
323:   const char     *owarmijo_type = TAOLINESEARCHOWARMIJO;

327:   tao->ops->setup = TaoSetUp_OWLQN;
328:   tao->ops->solve = TaoSolve_OWLQN;
329:   tao->ops->view = TaoView_OWLQN;
330:   tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
331:   tao->ops->destroy = TaoDestroy_OWLQN;

333:   PetscNewLog(tao,&lmP);
334:   lmP->D = 0;
335:   lmP->M = 0;
336:   lmP->GV = 0;
337:   lmP->Xold = 0;
338:   lmP->Gold = 0;
339:   lmP->lambda = 1.0;

341:   tao->data = (void*)lmP;
342:   /* Override default settings (unless already changed) */
343:   if (!tao->max_it_changed) tao->max_it = 2000;
344:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

346:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
347:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
348:   TaoLineSearchSetType(tao->linesearch,owarmijo_type);
349:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
350:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
351:   return(0);
352: }