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) dptr[i] = 0.0;
 21:   }
 22:   VecRestoreArray(d, &dptr);
 23:   VecRestoreArrayRead(g, &gptr);
 24:   return 0;
 25: }

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

 33:   VecGetOwnershipRange(x, &low, &high);
 34:   VecGetOwnershipRange(gv, &low1, &high1);

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

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

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

 62:   /* Check convergence criteria */
 63:   TaoComputeObjectiveAndGradient(tao, tao->solution, &f, tao->gradient);
 64:   VecCopy(tao->gradient, lmP->GV);
 65:   ComputePseudoGrad_OWLQN(tao->solution, lmP->GV, lmP->lambda);
 66:   VecNorm(lmP->GV, NORM_2, &gnorm);

 69:   tao->reason = TAO_CONTINUE_ITERATING;
 70:   TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
 71:   TaoMonitor(tao, iter, f, gnorm, 0.0, step);
 72:   PetscUseTypeMethod(tao, convergencetest, tao->cnvP);
 73:   if (tao->reason != TAO_CONTINUE_ITERATING) return 0;

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

 79:   /* Set counter for gradient/reset steps */
 80:   lmP->bfgs  = 0;
 81:   lmP->sgrad = 0;
 82:   lmP->grad  = 0;

 84:   /* Have not converged; continue with Newton method */
 85:   while (tao->reason == TAO_CONTINUE_ITERATING) {
 86:     /* Call general purpose update function */
 87:     PetscTryTypeMethod(tao, update, tao->niter, tao->user_update);

 89:     /* Compute direction */
 90:     MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
 91:     MatSolve(lmP->M, lmP->GV, lmP->D);

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

 95:     ++lmP->bfgs;

 97:     /* Check for success (descent direction) */
 98:     VecDot(lmP->D, lmP->GV, &gdx);
 99:     if ((gdx <= 0.0) || PetscIsInfOrNanReal(gdx)) {
100:       /* Step is not descent or direction produced not a number
101:          We can assert bfgsUpdates > 1 in this case because
102:          the first solve produces the scaled gradient direction,
103:          which is guaranteed to be descent

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

108:       delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm * gnorm);
109:       MatLMVMSetJ0Scale(lmP->M, delta);
110:       MatLMVMReset(lmP->M, PETSC_FALSE);
111:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
112:       MatSolve(lmP->M, lmP->GV, lmP->D);

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

116:       lmP->bfgs = 1;
117:       ++lmP->sgrad;
118:       stepType = OWLQN_SCALED_GRADIENT;
119:     } else {
120:       if (1 == lmP->bfgs) {
121:         /* The first BFGS direction is always the scaled gradient */
122:         ++lmP->sgrad;
123:         stepType = OWLQN_SCALED_GRADIENT;
124:       } else {
125:         ++lmP->bfgs;
126:         stepType = OWLQN_BFGS;
127:       }
128:     }

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

132:     /* Perform the linesearch */
133:     fold = f;
134:     VecCopy(tao->solution, lmP->Xold);
135:     VecCopy(tao->gradient, lmP->Gold);

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

140:     while (((int)ls_status < 0) && (stepType != OWLQN_GRADIENT)) {
141:       /* Reset factors and use scaled gradient step */
142:       f = fold;
143:       VecCopy(lmP->Xold, tao->solution);
144:       VecCopy(lmP->Gold, tao->gradient);
145:       VecCopy(tao->gradient, lmP->GV);

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

149:       switch (stepType) {
150:       case OWLQN_BFGS:
151:         /* Failed to obtain acceptable iterate with BFGS step
152:            Attempt to use the scaled gradient direction */

154:         delta = 2.0 * PetscMax(1.0, PetscAbsScalar(f)) / (gnorm * gnorm);
155:         MatLMVMSetJ0Scale(lmP->M, delta);
156:         MatLMVMReset(lmP->M, PETSC_FALSE);
157:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
158:         MatSolve(lmP->M, lmP->GV, lmP->D);

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

162:         lmP->bfgs = 1;
163:         ++lmP->sgrad;
164:         stepType = OWLQN_SCALED_GRADIENT;
165:         break;

167:       case OWLQN_SCALED_GRADIENT:
168:         /* The scaled gradient step did not produce a new iterate;
169:            attempt to use the gradient direction.
170:            Need to make sure we are not using a different diagonal scaling */
171:         MatLMVMSetJ0Scale(lmP->M, 1.0);
172:         MatLMVMReset(lmP->M, PETSC_FALSE);
173:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
174:         MatSolve(lmP->M, lmP->GV, lmP->D);

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

178:         lmP->bfgs = 1;
179:         ++lmP->grad;
180:         stepType = OWLQN_GRADIENT;
181:         break;
182:       }
183:       VecScale(lmP->D, -1.0);

185:       /* Perform the linesearch */
186:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
187:       TaoAddLineSearchCounts(tao);
188:     }

190:     if ((int)ls_status < 0) {
191:       /* Failed to find an improving point*/
192:       f = fold;
193:       VecCopy(lmP->Xold, tao->solution);
194:       VecCopy(lmP->Gold, tao->gradient);
195:       VecCopy(tao->gradient, lmP->GV);
196:       step = 0.0;
197:     } else {
198:       /* a little hack here, because that gv is used to store g */
199:       VecCopy(lmP->GV, tao->gradient);
200:     }

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

204:     /* Check for termination */

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

208:     iter++;
209:     TaoLogConvergenceHistory(tao, f, gnorm, 0.0, tao->ksp_its);
210:     TaoMonitor(tao, iter, f, gnorm, 0.0, step);
211:     PetscUseTypeMethod(tao, convergencetest, tao->cnvP);

213:     if ((int)ls_status < 0) break;
214:   }
215:   return 0;
216: }

218: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
219: {
220:   TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;
221:   PetscInt   n, N;

223:   /* Existence of tao->solution checked in TaoSetUp() */
224:   if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
225:   if (!tao->stepdirection) VecDuplicate(tao->solution, &tao->stepdirection);
226:   if (!lmP->D) VecDuplicate(tao->solution, &lmP->D);
227:   if (!lmP->GV) VecDuplicate(tao->solution, &lmP->GV);
228:   if (!lmP->Xold) VecDuplicate(tao->solution, &lmP->Xold);
229:   if (!lmP->Gold) VecDuplicate(tao->solution, &lmP->Gold);

231:   /* Create matrix for the limited memory approximation */
232:   VecGetLocalSize(tao->solution, &n);
233:   VecGetSize(tao->solution, &N);
234:   MatCreateLMVMBFGS(((PetscObject)tao)->comm, n, N, &lmP->M);
235:   MatLMVMAllocate(lmP->M, tao->solution, tao->gradient);
236:   return 0;
237: }

239: /* ---------------------------------------------------------- */
240: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
241: {
242:   TAO_OWLQN *lmP = (TAO_OWLQN *)tao->data;

244:   if (tao->setupcalled) {
245:     VecDestroy(&lmP->Xold);
246:     VecDestroy(&lmP->Gold);
247:     VecDestroy(&lmP->D);
248:     MatDestroy(&lmP->M);
249:     VecDestroy(&lmP->GV);
250:   }
251:   PetscFree(tao->data);
252:   return 0;
253: }

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

260:   PetscOptionsHeadBegin(PetscOptionsObject, "Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
261:   PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight", "", 100, &lmP->lambda, NULL);
262:   PetscOptionsHeadEnd();
263:   TaoLineSearchSetFromOptions(tao->linesearch);
264:   return 0;
265: }

267: /*------------------------------------------------------------*/
268: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
269: {
270:   TAO_OWLQN *lm = (TAO_OWLQN *)tao->data;
271:   PetscBool  isascii;

273:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
274:   if (isascii) {
275:     PetscViewerASCIIPushTab(viewer);
276:     PetscViewerASCIIPrintf(viewer, "BFGS steps: %" PetscInt_FMT "\n", lm->bfgs);
277:     PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %" PetscInt_FMT "\n", lm->sgrad);
278:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %" PetscInt_FMT "\n", lm->grad);
279:     PetscViewerASCIIPopTab(viewer);
280:   }
281:   return 0;
282: }

284: /* ---------------------------------------------------------- */
285: /*MC
286:   TAOOWLQN - orthant-wise limited memory quasi-newton algorithm

288: . - tao_owlqn_lambda - regulariser weight

290:   Level: beginner
291: M*/

293: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
294: {
295:   TAO_OWLQN  *lmP;
296:   const char *owarmijo_type = TAOLINESEARCHOWARMIJO;

298:   tao->ops->setup          = TaoSetUp_OWLQN;
299:   tao->ops->solve          = TaoSolve_OWLQN;
300:   tao->ops->view           = TaoView_OWLQN;
301:   tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
302:   tao->ops->destroy        = TaoDestroy_OWLQN;

304:   PetscNew(&lmP);
305:   lmP->D      = NULL;
306:   lmP->M      = NULL;
307:   lmP->GV     = NULL;
308:   lmP->Xold   = NULL;
309:   lmP->Gold   = NULL;
310:   lmP->lambda = 1.0;

312:   tao->data = (void *)lmP;
313:   /* Override default settings (unless already changed) */
314:   if (!tao->max_it_changed) tao->max_it = 2000;
315:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;

317:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
318:   PetscObjectIncrementTabLevel((PetscObject)tao->linesearch, (PetscObject)tao, 1);
319:   TaoLineSearchSetType(tao->linesearch, owarmijo_type);
320:   TaoLineSearchUseTaoRoutines(tao->linesearch, tao);
321:   TaoLineSearchSetOptionsPrefix(tao->linesearch, tao->hdr.prefix);
322:   return 0;
323: }