Actual source code: owlqn.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/owlqn/owlqn.h>

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

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

 19:   ierr=VecGetOwnershipRange(d,&low,&high);
 20:   ierr=VecGetOwnershipRange(g,&low1,&high1);

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

 36: static PetscErrorCode ComputePseudoGrad_OWLQN(Vec x, Vec gv, PetscReal lambda)
 37: {
 38:   PetscErrorCode  ierr;
 39:   const PetscReal *xptr;
 40:   PetscReal       *gptr;
 41:   PetscInt        low,high,low1,high1,i;

 44:   ierr=VecGetOwnershipRange(x,&low,&high);
 45:   ierr=VecGetOwnershipRange(gv,&low1,&high1);

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

 63: static PetscErrorCode TaoSolve_OWLQN(Tao tao)
 64: {
 65:   TAO_OWLQN                    *lmP = (TAO_OWLQN *)tao->data;
 66:   PetscReal                    f, fold, gdx, gnorm;
 67:   PetscReal                    step = 1.0;
 68:   PetscReal                    delta;
 69:   PetscErrorCode               ierr;
 70:   PetscInt                     stepType;
 71:   PetscInt                     iter = 0;
 72:   TaoConvergedReason           reason = TAO_CONTINUE_ITERATING;
 73:   TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

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

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

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

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

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

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

 91:   TaoMonitor(tao, iter, f, gnorm, 0.0, step, &reason);
 92:   if (reason != TAO_CONTINUE_ITERATING) return(0);

 94:   /* Set initial scaling for the function */
 95:   if (f != 0.0) {
 96:     delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
 97:   } else {
 98:     delta = 2.0 / (gnorm*gnorm);
 99:   }
100:   MatLMVMSetDelta(lmP->M,delta);

102:   /* Set counter for gradient/reset steps */
103:   lmP->bfgs = 0;
104:   lmP->sgrad = 0;
105:   lmP->grad = 0;

107:   /* Have not converged; continue with Newton method */
108:   while (reason == TAO_CONTINUE_ITERATING) {
109:     /* Compute direction */
110:     MatLMVMUpdate(lmP->M,tao->solution,tao->gradient);
111:     MatLMVMSolve(lmP->M, lmP->GV, lmP->D);

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

115:     ++lmP->bfgs;

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

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

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

129:       if (f != 0.0) {
130:         delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
131:       } else {
132:         delta = 2.0 / (gnorm*gnorm);
133:       }
134:       MatLMVMSetDelta(lmP->M, delta);
135:       MatLMVMReset(lmP->M);
136:       MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
137:       MatLMVMSolve(lmP->M,lmP->GV, lmP->D);

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

141:       lmP->bfgs = 1;
142:       ++lmP->sgrad;
143:       stepType = OWLQN_SCALED_GRADIENT;
144:     } else {
145:       if (1 == lmP->bfgs) {
146:         /* The first BFGS direction is always the scaled gradient */
147:         ++lmP->sgrad;
148:         stepType = OWLQN_SCALED_GRADIENT;
149:       } else {
150:         ++lmP->bfgs;
151:         stepType = OWLQN_BFGS;
152:       }
153:     }

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

157:     /* Perform the linesearch */
158:     fold = f;
159:     VecCopy(tao->solution, lmP->Xold);
160:     VecCopy(tao->gradient, lmP->Gold);

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

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

167:       /* Reset factors and use scaled gradient step */
168:       f = fold;
169:       VecCopy(lmP->Xold, tao->solution);
170:       VecCopy(lmP->Gold, tao->gradient);
171:       VecCopy(tao->gradient, lmP->GV);

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

175:       switch(stepType) {
176:       case OWLQN_BFGS:
177:         /* Failed to obtain acceptable iterate with BFGS step
178:            Attempt to use the scaled gradient direction */

180:         if (f != 0.0) {
181:           delta = 2.0 * PetscAbsScalar(f) / (gnorm*gnorm);
182:         } else {
183:           delta = 2.0 / (gnorm*gnorm);
184:         }
185:         MatLMVMSetDelta(lmP->M, delta);
186:         MatLMVMReset(lmP->M);
187:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
188:         MatLMVMSolve(lmP->M, lmP->GV, lmP->D);

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

192:         lmP->bfgs = 1;
193:         ++lmP->sgrad;
194:         stepType = OWLQN_SCALED_GRADIENT;
195:         break;

197:       case OWLQN_SCALED_GRADIENT:
198:         /* The scaled gradient step did not produce a new iterate;
199:            attempt to use the gradient direction.
200:            Need to make sure we are not using a different diagonal scaling */
201:         MatLMVMSetDelta(lmP->M, 1.0);
202:         MatLMVMReset(lmP->M);
203:         MatLMVMUpdate(lmP->M, tao->solution, tao->gradient);
204:         MatLMVMSolve(lmP->M, lmP->GV, lmP->D);

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

208:         lmP->bfgs = 1;
209:         ++lmP->grad;
210:         stepType = OWLQN_GRADIENT;
211:         break;
212:       }
213:       VecScale(lmP->D, -1.0);


216:       /* Perform the linesearch */
217:       TaoLineSearchApply(tao->linesearch, tao->solution, &f, lmP->GV, lmP->D, &step, &ls_status);
218:       TaoAddLineSearchCounts(tao);
219:     }

221:     if ((int)ls_status < 0) {
222:       /* Failed to find an improving point*/
223:       f = fold;
224:       VecCopy(lmP->Xold, tao->solution);
225:       VecCopy(lmP->Gold, tao->gradient);
226:       VecCopy(tao->gradient, lmP->GV);
227:       step = 0.0;
228:     } else {
229:       /* a little hack here, because that gv is used to store g */
230:       VecCopy(lmP->GV, tao->gradient);
231:     }

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

235:     /* Check for termination */

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

239:     iter++;
240:     TaoMonitor(tao,iter,f,gnorm,0.0,step,&reason);

242:     if ((int)ls_status < 0) break;
243:   }
244:   return(0);
245: }

249: static PetscErrorCode TaoSetUp_OWLQN(Tao tao)
250: {
251:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;
252:   PetscInt       n,N;

256:   /* Existence of tao->solution checked in TaoSetUp() */
257:   if (!tao->gradient) {VecDuplicate(tao->solution,&tao->gradient);  }
258:   if (!tao->stepdirection) {VecDuplicate(tao->solution,&tao->stepdirection);  }
259:   if (!lmP->D) {VecDuplicate(tao->solution,&lmP->D);  }
260:   if (!lmP->GV) {VecDuplicate(tao->solution,&lmP->GV);  }
261:   if (!lmP->Xold) {VecDuplicate(tao->solution,&lmP->Xold);  }
262:   if (!lmP->Gold) {VecDuplicate(tao->solution,&lmP->Gold);  }

264:   /* Create matrix for the limited memory approximation */
265:   VecGetLocalSize(tao->solution,&n);
266:   VecGetSize(tao->solution,&N);
267:   MatCreateLMVM(((PetscObject)tao)->comm,n,N,&lmP->M);
268:   MatLMVMAllocateVectors(lmP->M,tao->solution);
269:   return(0);
270: }

272: /* ---------------------------------------------------------- */
275: static PetscErrorCode TaoDestroy_OWLQN(Tao tao)
276: {
277:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;

281:   if (tao->setupcalled) {
282:     VecDestroy(&lmP->Xold);
283:     VecDestroy(&lmP->Gold);
284:     VecDestroy(&lmP->D);
285:     MatDestroy(&lmP->M);
286:     VecDestroy(&lmP->GV);
287:   }
288:   PetscFree(tao->data);
289:   return(0);
290: }

292: /*------------------------------------------------------------*/
295: static PetscErrorCode TaoSetFromOptions_OWLQN(PetscOptions *PetscOptionsObject,Tao tao)
296: {
297:   TAO_OWLQN      *lmP = (TAO_OWLQN *)tao->data;

301:   PetscOptionsHead(PetscOptionsObject,"Orthant-Wise Limited-memory method for Quasi-Newton unconstrained optimization");
302:   PetscOptionsReal("-tao_owlqn_lambda", "regulariser weight","", 100,&lmP->lambda,NULL);
303:   PetscOptionsTail();
304:   TaoLineSearchSetFromOptions(tao->linesearch);
305:   return(0);
306: }

308: /*------------------------------------------------------------*/
311: static PetscErrorCode TaoView_OWLQN(Tao tao, PetscViewer viewer)
312: {
313:   TAO_OWLQN      *lm = (TAO_OWLQN *)tao->data;
314:   PetscBool      isascii;

318:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
319:   if (isascii) {
320:     PetscViewerASCIIPushTab(viewer);
321:     PetscViewerASCIIPrintf(viewer, "BFGS steps: %D\n", lm->bfgs);
322:     PetscViewerASCIIPrintf(viewer, "Scaled gradient steps: %D\n", lm->sgrad);
323:     PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lm->grad);
324:     PetscViewerASCIIPopTab(viewer);
325:   }
326:   return(0);
327: }

329: /* ---------------------------------------------------------- */
330: /*MC
331:   TAOOWLQN - orthant-wise limited memory quasi-newton algorithm

333: . - tao_owlqn_lambda - regulariser weight

335:   Level: beginner
336: M*/


341: PETSC_EXTERN PetscErrorCode TaoCreate_OWLQN(Tao tao)
342: {
343:   TAO_OWLQN      *lmP;
344:   const char     *owarmijo_type = TAOLINESEARCHOWARMIJO;

348:   tao->ops->setup = TaoSetUp_OWLQN;
349:   tao->ops->solve = TaoSolve_OWLQN;
350:   tao->ops->view = TaoView_OWLQN;
351:   tao->ops->setfromoptions = TaoSetFromOptions_OWLQN;
352:   tao->ops->destroy = TaoDestroy_OWLQN;

354:   PetscNewLog(tao,&lmP);
355:   lmP->D = 0;
356:   lmP->M = 0;
357:   lmP->GV = 0;
358:   lmP->Xold = 0;
359:   lmP->Gold = 0;
360:   lmP->lambda = 1.0;

362:   tao->data = (void*)lmP;
363:   /* Override default settings (unless already changed) */
364:   if (!tao->max_it_changed) tao->max_it = 2000;
365:   if (!tao->max_funcs_changed) tao->max_funcs = 4000;
366:   if (!tao->fatol_changed) tao->fatol = 1.0e-4;
367:   if (!tao->frtol_changed) tao->frtol = 1.0e-4;

369:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch);
370:   TaoLineSearchSetType(tao->linesearch,owarmijo_type);
371:   TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
372:   TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
373:   return(0);
374: }