Actual source code: owlqn.c

petsc-3.14.6 2021-03-30
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);
 74:   VecCopy(tao->gradient, lmP->GV);
 75:   ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
 76:   VecNorm(lmP->GV,NORM_2,&gnorm);
 77:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_USER, "User provided compute function generated Inf or NaN");

 79:   tao->reason = TAO_CONTINUE_ITERATING;
 80:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
 81:   TaoMonitor(tao,iter,f,gnorm,0.0,step);
 82:   (*tao->ops->convergencetest)(tao,tao->cnvP);
 83:   if (tao->reason != TAO_CONTINUE_ITERATING) return(0);

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

 89:   /* Set counter for gradient/reset steps */
 90:   lmP->bfgs = 0;
 91:   lmP->sgrad = 0;
 92:   lmP->grad = 0;

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

101:     /* Compute direction */
102:     MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
103:     MatSolve(lmP->M, lmP->GV, lmP->D);

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

107:     ++lmP->bfgs;

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

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

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

121:       delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
122:       MatLMVMSetJ0Scale(lmP->M, delta);
123:       MatLMVMReset(lmP->M, PETSC_FALSE);
124:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
125:       MatSolve(lmP->M,lmP->GV, lmP->D);

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

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

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

145:     /* Perform the linesearch */
146:     fold = f;
147:     VecCopy(tao->solution, lmP->Xold);
148:     VecCopy(tao->gradient, lmP->Gold);

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

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

155:       /* Reset factors and use scaled gradient step */
156:       f = fold;
157:       VecCopy(lmP->Xold, tao->solution);
158:       VecCopy(lmP->Gold, tao->gradient);
159:       VecCopy(tao->gradient, lmP->GV);

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

163:       switch(stepType) {
164:       case OWLQN_BFGS:
165:         /* Failed to obtain acceptable iterate with BFGS step
166:            Attempt to use the scaled gradient direction */

168:         delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
169:         MatLMVMSetJ0Scale(lmP->M, delta);
170:         MatLMVMReset(lmP->M, PETSC_FALSE);
171:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
172:         MatSolve(lmP->M, lmP->GV, lmP->D);

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

176:         lmP->bfgs = 1;
177:         ++lmP->sgrad;
178:         stepType = OWLQN_SCALED_GRADIENT;
179:         break;

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

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

192:         lmP->bfgs = 1;
193:         ++lmP->grad;
194:         stepType = OWLQN_GRADIENT;
195:         break;
196:       }
197:       VecScale(lmP->D, -1.0);

199:       /* Perform the linesearch */
200:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
201:       TaoAddLineSearchCounts(tao);
202:     }

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

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

218:     /* Check for termination */

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

222:     iter++;
223:     TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
224:     TaoMonitor(tao,iter,f,gnorm,0.0,step);
225:     (*tao->ops->convergencetest)(tao,tao->cnvP);

227:     if ((int)ls_status < 0) break;
228:   }
229:   return(0);
230: }

232: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
233: {
234:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;
235:   PetscInt       n,N;

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

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

255: /* ---------------------------------------------------------- */
256: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
257: {
258:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;

262:   if (tao->setupcalled) {
263:     VecDestroy(&lmP->Xold);
264:     VecDestroy(&lmP->Gold);
265:     VecDestroy(&lmP->D);
266:     MatDestroy(&lmP->M);
267:     VecDestroy(&lmP->GV);
268:   }
269:   PetscFree(tao->data);
270:   return(0);
271: }

273: /*------------------------------------------------------------*/
274: static PetscErrorCode TaoSetFromOptions_OWLQN(PetscOptionItems *PetscOptionsObject,Tao tao)
275: {
276:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;

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

287: /*------------------------------------------------------------*/
288: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
289: {
290:   TAO_OWLQN      *lm = (TAO_OWLQN *)tao->data;
291:   PetscBool      isascii;

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

306: /* ---------------------------------------------------------- */
307: /*MC
308:   TAOOWLQN - orthant-wise limited memory quasi-newton algorithm

310: . - tao_owlqn_lambda - regulariser weight

312:   Level: beginner
313: M*/


316: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
317: {
318:   TAO_OWLQN      *lmP;
319:   const char     *owarmijo_type = TAOLINESEARCHOWARMIJO;

323:   tao->ops->setup = TaoSetUp_OWLQN;
324:   tao->ops->solve = TaoSolve_OWLQN;
325:   tao->ops->view = TaoView_OWLQN;
326:   tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
327:   tao->ops->destroy = TaoDestroy_OWLQN;

329:   PetscNewLog(tao,&lmP);
330:   lmP->D = NULL;
331:   lmP->M = NULL;
332:   lmP->GV = NULL;
333:   lmP->Xold = NULL;
334:   lmP->Gold = NULL;
335:   lmP->lambda = 1.0;

337:   tao->data = (void*)lmP;
338:   /* Override default settings (unless already changed) */
339:   if (!tao->max_it_changed) tao->max_it = 2000;
340:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

342:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
343:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
344:   TaoLineSearchSetType(tao->linesearch,owarmijo_type);
345:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
346:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
347:   return(0);
348: }