Actual source code: brgn.c
petsc-3.11.4 2019-09-28
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: }