Actual source code: owlqn.c

  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:   const PetscReal *gptr;
 11:   PetscReal       *dptr;
 12:   PetscInt        low,high,low1,high1,i;

 14:   VecGetOwnershipRange(d,&low,&high);
 15:   VecGetOwnershipRange(g,&low1,&high1);

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

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

 35:   VecGetOwnershipRange(x,&low,&high);
 36:   VecGetOwnershipRange(gv,&low1,&high1);

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

 52: static PetscErrorCode TaoSolve_OWLQN(Tao tao)
 53: {
 54:   TAO_OWLQN                    *lmP = (TAO_OWLQN *)tao->data;
 55:   PetscReal                    f, fold, gdx, gnorm;
 56:   PetscReal                    step = 1.0;
 57:   PetscReal                    delta;
 58:   PetscInt                     stepType;
 59:   PetscInt                     iter = 0;
 60:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

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

 66:   /* Check convergence criteria */
 67:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 68:   VecCopy(tao->gradient, lmP->GV);
 69:   ComputePseudoGrad_OWLQN(tao->solution,lmP->GV,lmP->lambda);
 70:   VecNorm(lmP->GV,NORM_2,&gnorm);

 73:   tao->reason = TAO_CONTINUE_ITERATING;
 74:   TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
 75:   TaoMonitor(tao,iter,f,gnorm,0.0,step);
 76:   (*tao->ops->convergencetest)(tao,tao->cnvP);
 77:   if (tao->reason != TAO_CONTINUE_ITERATING) return 0;

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

 83:   /* Set counter for gradient/reset steps */
 84:   lmP->bfgs = 0;
 85:   lmP->sgrad = 0;
 86:   lmP->grad = 0;

 88:   /* Have not converged; continue with Newton method */
 89:   while (tao->reason == TAO_CONTINUE_ITERATING) {
 90:     /* Call general purpose update function */
 91:     if (tao->ops->update) {
 92:       (*tao->ops->update)(tao, tao->niter, tao->user_update);
 93:     }

 95:     /* Compute direction */
 96:     MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
 97:     MatSolve(lmP->M, lmP->GV, lmP->D);

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

101:     ++lmP->bfgs;

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

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

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

115:       delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
116:       MatLMVMSetJ0Scale(lmP->M, delta);
117:       MatLMVMReset(lmP->M, PETSC_FALSE);
118:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
119:       MatSolve(lmP->M,lmP->GV, lmP->D);

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

123:       lmP->bfgs = 1;
124:       ++lmP->sgrad;
125:       stepType = OWLQN_SCALED_GRADIENT;
126:     } else {
127:       if (1 == lmP->bfgs) {
128:         /* The first BFGS direction is always the scaled gradient */
129:         ++lmP->sgrad;
130:         stepType = OWLQN_SCALED_GRADIENT;
131:       } else {
132:         ++lmP->bfgs;
133:         stepType = OWLQN_BFGS;
134:       }
135:     }

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

139:     /* Perform the linesearch */
140:     fold = f;
141:     VecCopy(tao->solution, lmP->Xold);
142:     VecCopy(tao->gradient, lmP->Gold);

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

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

149:       /* Reset factors and use scaled gradient step */
150:       f = fold;
151:       VecCopy(lmP->Xold, tao->solution);
152:       VecCopy(lmP->Gold, tao->gradient);
153:       VecCopy(tao->gradient, lmP->GV);

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

157:       switch(stepType) {
158:       case OWLQN_BFGS:
159:         /* Failed to obtain acceptable iterate with BFGS step
160:            Attempt to use the scaled gradient direction */

162:         delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm*gnorm);
163:         MatLMVMSetJ0Scale(lmP->M, delta);
164:         MatLMVMReset(lmP->M, PETSC_FALSE);
165:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
166:         MatSolve(lmP->M, lmP->GV, lmP->D);

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

170:         lmP->bfgs = 1;
171:         ++lmP->sgrad;
172:         stepType = OWLQN_SCALED_GRADIENT;
173:         break;

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

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

186:         lmP->bfgs = 1;
187:         ++lmP->grad;
188:         stepType = OWLQN_GRADIENT;
189:         break;
190:       }
191:       VecScale(lmP->D, -1.0);

193:       /* Perform the linesearch */
194:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
195:       TaoAddLineSearchCounts(tao);
196:     }

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

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

212:     /* Check for termination */

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

216:     iter++;
217:     TaoLogConvergenceHistory(tao,f,gnorm,0.0,tao->ksp_its);
218:     TaoMonitor(tao,iter,f,gnorm,0.0,step);
219:     (*tao->ops->convergencetest)(tao,tao->cnvP);

221:     if ((int)ls_status < 0) break;
222:   }
223:   return 0;
224: }

226: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
227: {
228:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;
229:   PetscInt       n,N;

231:   /* Existence of tao->solution checked in TaoSetUp() */
232:   if (!tao->gradient) VecDuplicate(tao->solution,&tao->gradient);
233:   if (!tao->stepdirection) VecDuplicate(tao->solution,&tao->stepdirection);
234:   if (!lmP->D) VecDuplicate(tao->solution,&lmP->D);
235:   if (!lmP->GV) VecDuplicate(tao->solution,&lmP->GV);
236:   if (!lmP->Xold) VecDuplicate(tao->solution,&lmP->Xold);
237:   if (!lmP->Gold) VecDuplicate(tao->solution,&lmP->Gold);

239:   /* Create matrix for the limited memory approximation */
240:   VecGetLocalSize(tao->solution,&n);
241:   VecGetSize(tao->solution,&N);
242:   MatCreateLMVMBFGS(((PetscObject)tao)->comm,n,N,&lmP->M);
243:   MatLMVMAllocate(lmP->M,tao->solution,tao->gradient);
244:   return 0;
245: }

247: /* ---------------------------------------------------------- */
248: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
249: {
250:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;

252:   if (tao->setupcalled) {
253:     VecDestroy(&lmP->Xold);
254:     VecDestroy(&lmP->Gold);
255:     VecDestroy(&lmP->D);
256:     MatDestroy(&lmP->M);
257:     VecDestroy(&lmP->GV);
258:   }
259:   PetscFree(tao->data);
260:   return 0;
261: }

263: /*------------------------------------------------------------*/
264: static PetscErrorCode TaoSetFromOptions_OWLQN(PetscOptionItems *PetscOptionsObject,Tao tao)
265: {
266:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;

268:   PetscOptionsHead(PetscOptionsObject,"Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
269:   PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,NULL);
270:   PetscOptionsTail();
271:   TaoLineSearchSetFromOptions(tao->linesearch);
272:   return 0;
273: }

275: /*------------------------------------------------------------*/
276: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
277: {
278:   TAO_OWLQN      *lm = (TAO_OWLQN *)tao->data;
279:   PetscBool      isascii;

281:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
282:   if (isascii) {
283:     PetscViewerASCIIPushTab(viewer);
284:     PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
285:     PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
286:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
287:     PetscViewerASCIIPopTab(viewer);
288:   }
289:   return 0;
290: }

292: /* ---------------------------------------------------------- */
293: /*MC
294:   TAOOWLQN - orthant-wise limited memory quasi-newton algorithm

296: . - tao_owlqn_lambda - regulariser weight

298:   Level: beginner
299: M*/

301: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
302: {
303:   TAO_OWLQN      *lmP;
304:   const char     *owarmijo_type = TAOLINESEARCHOWARMIJO;

306:   tao->ops->setup = TaoSetUp_OWLQN;
307:   tao->ops->solve = TaoSolve_OWLQN;
308:   tao->ops->view = TaoView_OWLQN;
309:   tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
310:   tao->ops->destroy = TaoDestroy_OWLQN;

312:   PetscNewLog(tao,&lmP);
313:   lmP->D = NULL;
314:   lmP->M = NULL;
315:   lmP->GV = NULL;
316:   lmP->Xold = NULL;
317:   lmP->Gold = NULL;
318:   lmP->lambda = 1.0;

320:   tao->data = (void*)lmP;
321:   /* Override default settings (unless already changed) */
322:   if (!tao->max_it_changed) tao->max_it = 2000;
323:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

325:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
326:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
327:   TaoLineSearchSetType(tao->linesearch,owarmijo_type);
328:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
329:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
330:   return 0;
331: }