Actual source code: brgn.c

petsc-3.11.2 2019-05-18
Report Typos and Errors
  1:  #include <../src/tao/leastsquares/impls/brgn/brgn.h>

  3: #define BRGN_REGULARIZATION_USER    0
  4: #define BRGN_REGULARIZATION_L2PROX  1
  5: #define BRGN_REGULARIZATION_L1DICT  2
  6: #define BRGN_REGULARIZATION_TYPES   3

  8: static const char *BRGN_REGULARIZATION_TABLE[64] = {"user","l2prox","l1dict"};

 10: static PetscErrorCode GNHessianProd(Mat H,Vec in,Vec out)
 11: {
 12:   TAO_BRGN              *gn;
 13:   PetscErrorCode        ierr;
 14: 
 16:   MatShellGetContext(H,&gn);
 17:   MatMult(gn->subsolver->ls_jac,in,gn->r_work);
 18:   MatMultTranspose(gn->subsolver->ls_jac,gn->r_work,out);
 19:   switch (gn->reg_type) {
 20:   case BRGN_REGULARIZATION_USER:
 21:     MatMult(gn->Hreg,in,gn->x_work);
 22:     VecAXPY(out,gn->lambda,gn->x_work);
 23:     break;
 24:   case BRGN_REGULARIZATION_L2PROX:
 25:     VecAXPY(out,gn->lambda,in);
 26:     break;
 27:   case BRGN_REGULARIZATION_L1DICT:
 28:     /* out = out + lambda*D'*(diag.*(D*in)) */
 29:     if (gn->D) {
 30:       MatMult(gn->D,in,gn->y);/* y = D*in */
 31:     } else {
 32:       VecCopy(in,gn->y);
 33:     }
 34:     VecPointwiseMult(gn->y_work,gn->diag,gn->y);   /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */
 35:     if (gn->D) {
 36:       MatMultTranspose(gn->D,gn->y_work,gn->x_work); /* x_work = D'*(diag.*(D*in)) */
 37:     } else {
 38:       VecCopy(gn->y_work,gn->x_work);
 39:     }
 40:     VecAXPY(out,gn->lambda,gn->x_work);
 41:     break;
 42:   }
 43:   return(0);
 44: }

 46: static PetscErrorCode GNObjectiveGradientEval(Tao tao,Vec X,PetscReal *fcn,Vec G,void *ptr)
 47: {
 48:   TAO_BRGN              *gn = (TAO_BRGN *)ptr;
 49:   PetscInt              K;                    /* dimension of D*X */
 50:   PetscScalar           yESum;
 51:   PetscErrorCode        ierr;
 52:   PetscReal             f_reg;
 53: 
 55:   /* compute objective *fcn*/
 56:   /* compute first term 0.5*||ls_res||_2^2 */
 57:   TaoComputeResidual(tao,X,tao->ls_res);
 58:   VecDot(tao->ls_res,tao->ls_res,fcn);
 59:   *fcn *= 0.5;
 60:   /* compute gradient G */
 61:   TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);
 62:   MatMultTranspose(tao->ls_jac,tao->ls_res,G);
 63:   /* add the regularization contribution */
 64:   switch (gn->reg_type) {
 65:   case BRGN_REGULARIZATION_USER:
 66:     (*gn->regularizerobjandgrad)(tao,X,&f_reg,gn->x_work,gn->reg_obj_ctx);
 67:     *fcn += gn->lambda*f_reg;
 68:     VecAXPY(G,gn->lambda,gn->x_work);
 69:     break;
 70:   case BRGN_REGULARIZATION_L2PROX:
 71:     /* compute f = f + lambda*0.5*(xk - xkm1)'*(xk - xkm1) */
 72:     VecAXPBYPCZ(gn->x_work,1.0,-1.0,0.0,X,gn->x_old);
 73:     VecDot(gn->x_work,gn->x_work,&f_reg);
 74:     *fcn += gn->lambda*0.5*f_reg;
 75:     /* compute G = G + lambda*(xk - xkm1) */
 76:     VecAXPBYPCZ(G,gn->lambda,-gn->lambda,1.0,X,gn->x_old);
 77:     break;
 78:   case BRGN_REGULARIZATION_L1DICT:
 79:     /* compute f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/
 80:     if (gn->D) {
 81:       MatMult(gn->D,X,gn->y);/* y = D*x */
 82:     } else {
 83:       VecCopy(X,gn->y);
 84:     }
 85:     VecPointwiseMult(gn->y_work,gn->y,gn->y);
 86:     VecShift(gn->y_work,gn->epsilon*gn->epsilon);
 87:     VecSqrtAbs(gn->y_work);  /* gn->y_work = sqrt(y.^2+epsilon^2) */
 88:     VecSum(gn->y_work,&yESum);
 89:     VecGetSize(gn->y,&K);
 90:     *fcn += gn->lambda*(yESum - K*gn->epsilon);
 91:     /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */
 92:     VecPointwiseDivide(gn->y_work,gn->y,gn->y_work); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */
 93:     if (gn->D) {
 94:       MatMultTranspose(gn->D,gn->y_work,gn->x_work);
 95:     } else {
 96:       VecCopy(gn->y_work,gn->x_work);
 97:     }
 98:     VecAXPY(G,gn->lambda,gn->x_work);
 99:     break;
100:   }
101:   return(0);
102: }

104: static PetscErrorCode GNComputeHessian(Tao tao,Vec X,Mat H,Mat Hpre,void *ptr)
105: {
106:   TAO_BRGN              *gn = (TAO_BRGN *)ptr;
108: 
110:   TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);

112:   switch (gn->reg_type) {
113:   case BRGN_REGULARIZATION_USER:
114:     (*gn->regularizerhessian)(tao,X,gn->Hreg,gn->reg_hess_ctx);
115:     break;
116:   case BRGN_REGULARIZATION_L2PROX:
117:     break;
118:   case BRGN_REGULARIZATION_L1DICT:
119:     /* calculate and store diagonal matrix as a vector: diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3* --> diag = epsilon^2 ./ sqrt(y.^2+epsilon^2).^3,where y = D*x */
120:     if (gn->D) {
121:       MatMult(gn->D,X,gn->y);/* y = D*x */
122:     } else {
123:       VecCopy(X,gn->y);
124:     }
125:     VecPointwiseMult(gn->y_work,gn->y,gn->y);
126:     VecShift(gn->y_work,gn->epsilon*gn->epsilon);
127:     VecCopy(gn->y_work,gn->diag);                  /* gn->diag = y.^2+epsilon^2 */
128:     VecSqrtAbs(gn->y_work);                        /* gn->y_work = sqrt(y.^2+epsilon^2) */
129:     VecPointwiseMult(gn->diag,gn->y_work,gn->diag);/* gn->diag = sqrt(y.^2+epsilon^2).^3 */
130:     VecReciprocal(gn->diag);
131:     VecScale(gn->diag,gn->epsilon*gn->epsilon);
132:     break;
133:   }
134:   return(0);
135: }

137: static PetscErrorCode GNHookFunction(Tao tao,PetscInt iter, void *ctx)
138: {
139:   TAO_BRGN              *gn = (TAO_BRGN *)ctx;
140:   PetscErrorCode        ierr;
141: 
143:   /* Update basic tao information from the subsolver */
144:   gn->parent->nfuncs = tao->nfuncs;
145:   gn->parent->ngrads = tao->ngrads;
146:   gn->parent->nfuncgrads = tao->nfuncgrads;
147:   gn->parent->nhess = tao->nhess;
148:   gn->parent->niter = tao->niter;
149:   gn->parent->ksp_its = tao->ksp_its;
150:   gn->parent->ksp_tot_its = tao->ksp_tot_its;
151:   TaoGetConvergedReason(tao,&gn->parent->reason);
152:   /* Update the solution vectors */
153:   if (iter == 0) {
154:     VecSet(gn->x_old,0.0);
155:   } else {
156:     VecCopy(tao->solution,gn->x_old);
157:     VecCopy(tao->solution,gn->parent->solution);
158:   }
159:   /* Update the gradient */
160:   VecCopy(tao->gradient,gn->parent->gradient);
161:   /* Call general purpose update function */
162:   if (gn->parent->ops->update) {
163:     (*gn->parent->ops->update)(gn->parent,gn->parent->niter,gn->parent->user_update);
164:   }
165:   return(0);
166: }

168: static PetscErrorCode TaoSolve_BRGN(Tao tao)
169: {
170:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
171:   PetscErrorCode        ierr;

174:   TaoSolve(gn->subsolver);
175:   /* Update basic tao information from the subsolver */
176:   tao->nfuncs = gn->subsolver->nfuncs;
177:   tao->ngrads = gn->subsolver->ngrads;
178:   tao->nfuncgrads = gn->subsolver->nfuncgrads;
179:   tao->nhess = gn->subsolver->nhess;
180:   tao->niter = gn->subsolver->niter;
181:   tao->ksp_its = gn->subsolver->ksp_its;
182:   tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
183:   TaoGetConvergedReason(gn->subsolver,&tao->reason);
184:   /* Update vectors */
185:   VecCopy(gn->subsolver->solution,tao->solution);
186:   VecCopy(gn->subsolver->gradient,tao->gradient);
187:   return(0);
188: }

190: static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao)
191: {
192:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
193:   PetscErrorCode        ierr;

196:   PetscOptionsHead(PetscOptionsObject,"least-squares problems with regularizer: ||f(x)||^2 + lambda*g(x), g(x) = ||xk-xkm1||^2 or ||Dx||_1 or user defined function.");
197:   PetscOptionsReal("-tao_brgn_regularizer_weight","regularizer weight (default 1e-4)","",gn->lambda,&gn->lambda,NULL);
198:   PetscOptionsReal("-tao_brgn_l1_smooth_epsilon","L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)","",gn->epsilon,&gn->epsilon,NULL);
199:   PetscOptionsEList("-tao_brgn_regularization_type","regularization type", "",BRGN_REGULARIZATION_TABLE,BRGN_REGULARIZATION_TYPES,BRGN_REGULARIZATION_TABLE[gn->reg_type],&gn->reg_type,NULL);
200:   PetscOptionsTail();
201:   TaoSetFromOptions(gn->subsolver);
202:   return(0);
203: }

205: static PetscErrorCode TaoView_BRGN(Tao tao,PetscViewer viewer)
206: {
207:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
208:   PetscErrorCode        ierr;

211:   PetscViewerASCIIPushTab(viewer);
212:   TaoView(gn->subsolver,viewer);
213:   PetscViewerASCIIPopTab(viewer);
214:   return(0);
215: }

217: static PetscErrorCode TaoSetUp_BRGN(Tao tao)
218: {
219:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
220:   PetscErrorCode        ierr;
221:   PetscBool             is_bnls,is_bntr,is_bntl;
222:   PetscInt              i,n,N,K; /* dict has size K*N*/

225:   if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualRoutine() must be called before setup!");
226:   PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNLS,&is_bnls);
227:   PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTR,&is_bntr);
228:   PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTL,&is_bntl);
229:   if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualJacobianRoutine() must be called before setup!");
230:   if (!tao->gradient) {
231:     VecDuplicate(tao->solution,&tao->gradient);
232:   }
233:   if (!gn->x_work) {
234:     VecDuplicate(tao->solution,&gn->x_work);
235:   }
236:   if (!gn->r_work) {
237:     VecDuplicate(tao->ls_res,&gn->r_work);
238:   }
239:   if (!gn->x_old) {
240:     VecDuplicate(tao->solution,&gn->x_old);
241:     VecSet(gn->x_old,0.0);
242:   }
243: 
244:   if (BRGN_REGULARIZATION_L1DICT == gn->reg_type) {
245:     if (gn->D) {
246:       MatGetSize(gn->D,&K,&N); /* Shell matrices still must have sizes defined. K = N for identity matrix, K=N-1 or N for gradient matrix */
247:     } else {
248:       VecGetSize(tao->solution,&K); /* If user does not setup dict matrix, use identiy matrix, K=N */
249:     }
250:     if (!gn->y) {
251:       VecCreate(PETSC_COMM_SELF,&gn->y);
252:       VecSetSizes(gn->y,PETSC_DECIDE,K);
253:       VecSetFromOptions(gn->y);
254:       VecSet(gn->y,0.0);

256:     }
257:     if (!gn->y_work) {
258:       VecDuplicate(gn->y,&gn->y_work);
259:     }
260:     if (!gn->diag) {
261:       VecDuplicate(gn->y,&gn->diag);
262:       VecSet(gn->diag,0.0);
263:     }
264:   }

266:   if (!tao->setupcalled) {
267:     /* Hessian setup */
268:     VecGetLocalSize(tao->solution,&n);
269:     VecGetSize(tao->solution,&N);
270:     MatSetSizes(gn->H,n,n,N,N);
271:     MatSetType(gn->H,MATSHELL);
272:     MatSetUp(gn->H);
273:     MatShellSetOperation(gn->H,MATOP_MULT,(void (*)(void))GNHessianProd);
274:     MatShellSetContext(gn->H,(void*)gn);
275:     /* Subsolver setup,include initial vector and dicttionary D */
276:     TaoSetUpdate(gn->subsolver,GNHookFunction,(void*)gn);
277:     TaoSetInitialVector(gn->subsolver,tao->solution);
278:     if (tao->bounded) {
279:       TaoSetVariableBounds(gn->subsolver,tao->XL,tao->XU);
280:     }
281:     TaoSetResidualRoutine(gn->subsolver,tao->ls_res,tao->ops->computeresidual,tao->user_lsresP);
282:     TaoSetJacobianResidualRoutine(gn->subsolver,tao->ls_jac,tao->ls_jac,tao->ops->computeresidualjacobian,tao->user_lsjacP);
283:     TaoSetObjectiveAndGradientRoutine(gn->subsolver,GNObjectiveGradientEval,(void*)gn);
284:     TaoSetHessianRoutine(gn->subsolver,gn->H,gn->H,GNComputeHessian,(void*)gn);
285:     /* Propagate some options down */
286:     TaoSetTolerances(gn->subsolver,tao->gatol,tao->grtol,tao->gttol);
287:     TaoSetMaximumIterations(gn->subsolver,tao->max_it);
288:     TaoSetMaximumFunctionEvaluations(gn->subsolver,tao->max_funcs);
289:     for (i=0; i<tao->numbermonitors; ++i) {
290:       TaoSetMonitor(gn->subsolver,tao->monitor[i],tao->monitorcontext[i],tao->monitordestroy[i]);
291:       PetscObjectReference((PetscObject)(tao->monitorcontext[i]));
292:     }
293:     TaoSetUp(gn->subsolver);
294:   }
295:   return(0);
296: }

298: static PetscErrorCode TaoDestroy_BRGN(Tao tao)
299: {
300:   TAO_BRGN              *gn = (TAO_BRGN *)tao->data;
301:   PetscErrorCode        ierr;

304:   if (tao->setupcalled) {
305:     VecDestroy(&tao->gradient);
306:     VecDestroy(&gn->x_work);
307:     VecDestroy(&gn->r_work);
308:     VecDestroy(&gn->x_old);
309:     VecDestroy(&gn->diag);
310:     VecDestroy(&gn->y);
311:     VecDestroy(&gn->y_work);
312:   }
313:   MatDestroy(&gn->H);
314:   MatDestroy(&gn->D);
315:   MatDestroy(&gn->Hreg);
316:   TaoDestroy(&gn->subsolver);
317:   gn->parent = NULL;
318:   PetscFree(tao->data);
319:   return(0);
320: }

322: /*MC
323:   TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares 
324:             problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL 
325:             that constructs the Gauss-Newton problem with the user-provided least-squares 
326:             residual and Jacobian. The algorithm offers both an L2-norm proximal point ("l2prox") 
327:             regularizer, and a L1-norm dictionary regularizer ("l1dict"), where we approximate the 
328:             L1-norm ||x||_1 by sum_i(sqrt(x_i^2+epsilon^2)-epsilon) with a small positive number epsilon.
329:             The user can also provide own regularization function.

331:   Options Database Keys:
332:   + -tao_brgn_regularization_type - regularization type ("user", "l2prox", "l1dict") (default "l2prox")
333:   . -tao_brgn_regularizer_weight  - regularizer weight (default 1e-4)
334:   - -tao_brgn_l1_smooth_epsilon   - L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)

336:   Level: beginner
337: M*/
338: PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao)
339: {
340:   TAO_BRGN       *gn;
342: 
344:   PetscNewLog(tao,&gn);
345: 
346:   tao->ops->destroy = TaoDestroy_BRGN;
347:   tao->ops->setup = TaoSetUp_BRGN;
348:   tao->ops->setfromoptions = TaoSetFromOptions_BRGN;
349:   tao->ops->view = TaoView_BRGN;
350:   tao->ops->solve = TaoSolve_BRGN;
351: 
352:   tao->data = (void*)gn;
353:   gn->reg_type = BRGN_REGULARIZATION_L2PROX;
354:   gn->lambda = 1e-4;
355:   gn->epsilon = 1e-6;
356:   gn->parent = tao;
357: 
358:   MatCreate(PetscObjectComm((PetscObject)tao),&gn->H);
359:   MatSetOptionsPrefix(gn->H,"tao_brgn_hessian_");
360: 
361:   TaoCreate(PetscObjectComm((PetscObject)tao),&gn->subsolver);
362:   TaoSetType(gn->subsolver,TAOBNLS);
363:   TaoSetOptionsPrefix(gn->subsolver,"tao_brgn_subsolver_");
364:   return(0);
365: }

367: /*@
368:   TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN

370:   Collective on Tao

372:   Level: advanced
373:   
374:   Input Parameters:
375: +  tao - the Tao solver context
376: -  subsolver - the Tao sub-solver context
377: @*/
378: PetscErrorCode TaoBRGNGetSubsolver(Tao tao,Tao *subsolver)
379: {
380:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;
381: 
383:   *subsolver = gn->subsolver;
384:   return(0);
385: }

387: /*@
388:   TaoBRGNSetRegularizerWeight - Set the regularizer weight for the Gauss-Newton least-squares algorithm

390:   Collective on Tao
391:   
392:   Input Parameters:
393: +  tao - the Tao solver context
394: -  lambda - L1-norm regularizer weight

396:   Level: beginner
397: @*/
398: PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao,PetscReal lambda)
399: {
400:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;
401: 
402:   /* Initialize lambda here */

405:   gn->lambda = lambda;
406:   return(0);
407: }

409: /*@
410:   TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm

412:   Collective on Tao
413:   
414:   Input Parameters:
415: +  tao - the Tao solver context
416: -  epsilon - L1-norm smooth approximation parameter

418:   Level: advanced
419: @*/
420: PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao,PetscReal epsilon)
421: {
422:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;
423: 
424:   /* Initialize epsilon here */

427:   gn->epsilon = epsilon;
428:   return(0);
429: }

431: /*@
432:    TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)

434:    Input Parameters:
435: +  tao  - the Tao context
436: .  dict - the user specified dictionary matrix.  We allow to set a null dictionary, which means identity matrix by default

438:     Level: advanced
439: @*/
440: PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao,Mat dict)
441: {
442:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;
446:   if (dict) {
449:     PetscObjectReference((PetscObject)dict);
450:   }
451:   MatDestroy(&gn->D);
452:   gn->D = dict;
453:   return(0);
454: }

456: /*@C
457:    TaoBRGNSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back 
458:    function into the algorithm.

460:    Input Parameters:
461:    + tao - the Tao context
462:    . func - function pointer for the regularizer value and gradient evaluation
463:    - ctx - user context for the regularizer

465:    Level: advanced
466: @*/
467: PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao,PetscErrorCode (*func)(Tao,Vec,PetscReal *,Vec,void*),void *ctx)
468: {
469:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;

473:   if (ctx) {
474:     gn->reg_obj_ctx = ctx;
475:   }
476:   if (func) {
477:     gn->regularizerobjandgrad = func;
478:   }
479:   return(0);
480: }

482: /*@C
483:    TaoBRGNSetRegularizerHessianRoutine - Sets the user-defined regularizer call-back 
484:    function into the algorithm.

486:    Input Parameters:
487:    + tao - the Tao context
488:    . Hreg - user-created matrix for the Hessian of the regularization term
489:    . func - function pointer for the regularizer Hessian evaluation
490:    - ctx - user context for the regularizer Hessian

492:    Level: advanced
493: @*/
494: PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao,Mat Hreg,PetscErrorCode (*func)(Tao,Vec,Mat,void*),void *ctx)
495: {
496:   TAO_BRGN       *gn = (TAO_BRGN *)tao->data;

501:   if (Hreg) {
504:   } else SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONG,"NULL Hessian detected! User must provide valid Hessian for the regularizer.");
505:   if (ctx) {
506:     gn->reg_hess_ctx = ctx;
507:   }
508:   if (func) {
509:     gn->regularizerhessian = func;
510:   }
511:   if (Hreg) {
512:     PetscObjectReference((PetscObject)Hreg);
513:     MatDestroy(&gn->Hreg);
514:     gn->Hreg = Hreg;
515:   }
516:   return(0);
517: }