Actual source code: brgn.c
petsc-3.12.5 2020-03-29
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_L2PURE 2
6: #define BRGN_REGULARIZATION_L1DICT 3
7: #define BRGN_REGULARIZATION_TYPES 4
9: static const char *BRGN_REGULARIZATION_TABLE[64] = {"user","l2prox","l2pure","l1dict"};
11: static PetscErrorCode GNHessianProd(Mat H,Vec in,Vec out)
12: {
13: TAO_BRGN *gn;
14: PetscErrorCode ierr;
15:
17: MatShellGetContext(H,&gn);
18: MatMult(gn->subsolver->ls_jac,in,gn->r_work);
19: MatMultTranspose(gn->subsolver->ls_jac,gn->r_work,out);
20: switch (gn->reg_type) {
21: case BRGN_REGULARIZATION_USER:
22: MatMult(gn->Hreg,in,gn->x_work);
23: VecAXPY(out,gn->lambda,gn->x_work);
24: break;
25: case BRGN_REGULARIZATION_L2PURE:
26: VecAXPY(out,gn->lambda,in);
27: break;
28: case BRGN_REGULARIZATION_L2PROX:
29: VecAXPY(out,gn->lambda,in);
30: break;
31: case BRGN_REGULARIZATION_L1DICT:
32: /* out = out + lambda*D'*(diag.*(D*in)) */
33: if (gn->D) {
34: MatMult(gn->D,in,gn->y);/* y = D*in */
35: } else {
36: VecCopy(in,gn->y);
37: }
38: VecPointwiseMult(gn->y_work,gn->diag,gn->y); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */
39: if (gn->D) {
40: MatMultTranspose(gn->D,gn->y_work,gn->x_work); /* x_work = D'*(diag.*(D*in)) */
41: } else {
42: VecCopy(gn->y_work,gn->x_work);
43: }
44: VecAXPY(out,gn->lambda,gn->x_work);
45: break;
46: }
47: return(0);
48: }
50: static PetscErrorCode GNObjectiveGradientEval(Tao tao,Vec X,PetscReal *fcn,Vec G,void *ptr)
51: {
52: TAO_BRGN *gn = (TAO_BRGN *)ptr;
53: PetscInt K; /* dimension of D*X */
54: PetscScalar yESum;
55: PetscErrorCode ierr;
56: PetscReal f_reg;
57:
59: /* compute objective *fcn*/
60: /* compute first term 0.5*||ls_res||_2^2 */
61: TaoComputeResidual(tao,X,tao->ls_res);
62: VecDot(tao->ls_res,tao->ls_res,fcn);
63: *fcn *= 0.5;
64: /* compute gradient G */
65: TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);
66: MatMultTranspose(tao->ls_jac,tao->ls_res,G);
67: /* add the regularization contribution */
68: switch (gn->reg_type) {
69: case BRGN_REGULARIZATION_USER:
70: (*gn->regularizerobjandgrad)(tao,X,&f_reg,gn->x_work,gn->reg_obj_ctx);
71: *fcn += gn->lambda*f_reg;
72: VecAXPY(G,gn->lambda,gn->x_work);
73: break;
74: case BRGN_REGULARIZATION_L2PURE:
75: /* compute f = f + lambda*0.5*xk'*xk */
76: VecDot(X,X,&f_reg);
77: *fcn += gn->lambda*0.5*f_reg;
78: /* compute G = G + lambda*xk */
79: VecAXPY(G,gn->lambda,X);
80: break;
81: case BRGN_REGULARIZATION_L2PROX:
82: /* compute f = f + lambda*0.5*(xk - xkm1)'*(xk - xkm1) */
83: VecAXPBYPCZ(gn->x_work,1.0,-1.0,0.0,X,gn->x_old);
84: VecDot(gn->x_work,gn->x_work,&f_reg);
85: *fcn += gn->lambda*0.5*f_reg;
86: /* compute G = G + lambda*(xk - xkm1) */
87: VecAXPBYPCZ(G,gn->lambda,-gn->lambda,1.0,X,gn->x_old);
88: break;
89: case BRGN_REGULARIZATION_L1DICT:
90: /* compute f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/
91: if (gn->D) {
92: MatMult(gn->D,X,gn->y);/* y = D*x */
93: } else {
94: VecCopy(X,gn->y);
95: }
96: VecPointwiseMult(gn->y_work,gn->y,gn->y);
97: VecShift(gn->y_work,gn->epsilon*gn->epsilon);
98: VecSqrtAbs(gn->y_work); /* gn->y_work = sqrt(y.^2+epsilon^2) */
99: VecSum(gn->y_work,&yESum);
100: VecGetSize(gn->y,&K);
101: *fcn += gn->lambda*(yESum - K*gn->epsilon);
102: /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */
103: VecPointwiseDivide(gn->y_work,gn->y,gn->y_work); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */
104: if (gn->D) {
105: MatMultTranspose(gn->D,gn->y_work,gn->x_work);
106: } else {
107: VecCopy(gn->y_work,gn->x_work);
108: }
109: VecAXPY(G,gn->lambda,gn->x_work);
110: break;
111: }
112: return(0);
113: }
115: static PetscErrorCode GNComputeHessian(Tao tao,Vec X,Mat H,Mat Hpre,void *ptr)
116: {
117: TAO_BRGN *gn = (TAO_BRGN *)ptr;
119:
121: TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);
123: switch (gn->reg_type) {
124: case BRGN_REGULARIZATION_USER:
125: (*gn->regularizerhessian)(tao,X,gn->Hreg,gn->reg_hess_ctx);
126: break;
127: case BRGN_REGULARIZATION_L2PURE:
128: break;
129: case BRGN_REGULARIZATION_L2PROX:
130: break;
131: case BRGN_REGULARIZATION_L1DICT:
132: /* 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 */
133: if (gn->D) {
134: MatMult(gn->D,X,gn->y);/* y = D*x */
135: } else {
136: VecCopy(X,gn->y);
137: }
138: VecPointwiseMult(gn->y_work,gn->y,gn->y);
139: VecShift(gn->y_work,gn->epsilon*gn->epsilon);
140: VecCopy(gn->y_work,gn->diag); /* gn->diag = y.^2+epsilon^2 */
141: VecSqrtAbs(gn->y_work); /* gn->y_work = sqrt(y.^2+epsilon^2) */
142: VecPointwiseMult(gn->diag,gn->y_work,gn->diag);/* gn->diag = sqrt(y.^2+epsilon^2).^3 */
143: VecReciprocal(gn->diag);
144: VecScale(gn->diag,gn->epsilon*gn->epsilon);
145: break;
146: }
147: return(0);
148: }
150: static PetscErrorCode GNHookFunction(Tao tao,PetscInt iter, void *ctx)
151: {
152: TAO_BRGN *gn = (TAO_BRGN *)ctx;
153: PetscErrorCode ierr;
154:
156: /* Update basic tao information from the subsolver */
157: gn->parent->nfuncs = tao->nfuncs;
158: gn->parent->ngrads = tao->ngrads;
159: gn->parent->nfuncgrads = tao->nfuncgrads;
160: gn->parent->nhess = tao->nhess;
161: gn->parent->niter = tao->niter;
162: gn->parent->ksp_its = tao->ksp_its;
163: gn->parent->ksp_tot_its = tao->ksp_tot_its;
164: TaoGetConvergedReason(tao,&gn->parent->reason);
165: /* Update the solution vectors */
166: if (iter == 0) {
167: VecSet(gn->x_old,0.0);
168: } else {
169: VecCopy(tao->solution,gn->x_old);
170: VecCopy(tao->solution,gn->parent->solution);
171: }
172: /* Update the gradient */
173: VecCopy(tao->gradient,gn->parent->gradient);
174: /* Call general purpose update function */
175: if (gn->parent->ops->update) {
176: (*gn->parent->ops->update)(gn->parent,gn->parent->niter,gn->parent->user_update);
177: }
178: return(0);
179: }
181: static PetscErrorCode TaoSolve_BRGN(Tao tao)
182: {
183: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
184: PetscErrorCode ierr;
187: TaoSolve(gn->subsolver);
188: /* Update basic tao information from the subsolver */
189: tao->nfuncs = gn->subsolver->nfuncs;
190: tao->ngrads = gn->subsolver->ngrads;
191: tao->nfuncgrads = gn->subsolver->nfuncgrads;
192: tao->nhess = gn->subsolver->nhess;
193: tao->niter = gn->subsolver->niter;
194: tao->ksp_its = gn->subsolver->ksp_its;
195: tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
196: TaoGetConvergedReason(gn->subsolver,&tao->reason);
197: /* Update vectors */
198: VecCopy(gn->subsolver->solution,tao->solution);
199: VecCopy(gn->subsolver->gradient,tao->gradient);
200: return(0);
201: }
203: static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao)
204: {
205: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
206: PetscErrorCode ierr;
209: 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.");
210: PetscOptionsReal("-tao_brgn_regularizer_weight","regularizer weight (default 1e-4)","",gn->lambda,&gn->lambda,NULL);
211: 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);
212: PetscOptionsEList("-tao_brgn_regularization_type","regularization type", "",BRGN_REGULARIZATION_TABLE,BRGN_REGULARIZATION_TYPES,BRGN_REGULARIZATION_TABLE[gn->reg_type],&gn->reg_type,NULL);
213: PetscOptionsTail();
214: TaoSetFromOptions(gn->subsolver);
215: return(0);
216: }
218: static PetscErrorCode TaoView_BRGN(Tao tao,PetscViewer viewer)
219: {
220: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
221: PetscErrorCode ierr;
224: PetscViewerASCIIPushTab(viewer);
225: TaoView(gn->subsolver,viewer);
226: PetscViewerASCIIPopTab(viewer);
227: return(0);
228: }
230: static PetscErrorCode TaoSetUp_BRGN(Tao tao)
231: {
232: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
233: PetscErrorCode ierr;
234: PetscBool is_bnls,is_bntr,is_bntl;
235: PetscInt i,n,N,K; /* dict has size K*N*/
238: if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualRoutine() must be called before setup!");
239: PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNLS,&is_bnls);
240: PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTR,&is_bntr);
241: PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTL,&is_bntl);
242: if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualJacobianRoutine() must be called before setup!");
243: if (!tao->gradient) {
244: VecDuplicate(tao->solution,&tao->gradient);
245: }
246: if (!gn->x_work) {
247: VecDuplicate(tao->solution,&gn->x_work);
248: }
249: if (!gn->r_work) {
250: VecDuplicate(tao->ls_res,&gn->r_work);
251: }
252: if (!gn->x_old) {
253: VecDuplicate(tao->solution,&gn->x_old);
254: VecSet(gn->x_old,0.0);
255: }
256:
257: if (BRGN_REGULARIZATION_L1DICT == gn->reg_type) {
258: if (gn->D) {
259: 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 */
260: } else {
261: VecGetSize(tao->solution,&K); /* If user does not setup dict matrix, use identiy matrix, K=N */
262: }
263: if (!gn->y) {
264: VecCreate(PETSC_COMM_SELF,&gn->y);
265: VecSetSizes(gn->y,PETSC_DECIDE,K);
266: VecSetFromOptions(gn->y);
267: VecSet(gn->y,0.0);
269: }
270: if (!gn->y_work) {
271: VecDuplicate(gn->y,&gn->y_work);
272: }
273: if (!gn->diag) {
274: VecDuplicate(gn->y,&gn->diag);
275: VecSet(gn->diag,0.0);
276: }
277: }
279: if (!tao->setupcalled) {
280: /* Hessian setup */
281: VecGetLocalSize(tao->solution,&n);
282: VecGetSize(tao->solution,&N);
283: MatSetSizes(gn->H,n,n,N,N);
284: MatSetType(gn->H,MATSHELL);
285: MatSetUp(gn->H);
286: MatShellSetOperation(gn->H,MATOP_MULT,(void (*)(void))GNHessianProd);
287: MatShellSetContext(gn->H,(void*)gn);
288: /* Subsolver setup,include initial vector and dicttionary D */
289: TaoSetUpdate(gn->subsolver,GNHookFunction,(void*)gn);
290: TaoSetInitialVector(gn->subsolver,tao->solution);
291: if (tao->bounded) {
292: TaoSetVariableBounds(gn->subsolver,tao->XL,tao->XU);
293: }
294: TaoSetResidualRoutine(gn->subsolver,tao->ls_res,tao->ops->computeresidual,tao->user_lsresP);
295: TaoSetJacobianResidualRoutine(gn->subsolver,tao->ls_jac,tao->ls_jac,tao->ops->computeresidualjacobian,tao->user_lsjacP);
296: TaoSetObjectiveAndGradientRoutine(gn->subsolver,GNObjectiveGradientEval,(void*)gn);
297: TaoSetHessianRoutine(gn->subsolver,gn->H,gn->H,GNComputeHessian,(void*)gn);
298: /* Propagate some options down */
299: TaoSetTolerances(gn->subsolver,tao->gatol,tao->grtol,tao->gttol);
300: TaoSetMaximumIterations(gn->subsolver,tao->max_it);
301: TaoSetMaximumFunctionEvaluations(gn->subsolver,tao->max_funcs);
302: for (i=0; i<tao->numbermonitors; ++i) {
303: TaoSetMonitor(gn->subsolver,tao->monitor[i],tao->monitorcontext[i],tao->monitordestroy[i]);
304: PetscObjectReference((PetscObject)(tao->monitorcontext[i]));
305: }
306: TaoSetUp(gn->subsolver);
307: }
308: return(0);
309: }
311: static PetscErrorCode TaoDestroy_BRGN(Tao tao)
312: {
313: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
314: PetscErrorCode ierr;
317: if (tao->setupcalled) {
318: VecDestroy(&tao->gradient);
319: VecDestroy(&gn->x_work);
320: VecDestroy(&gn->r_work);
321: VecDestroy(&gn->x_old);
322: VecDestroy(&gn->diag);
323: VecDestroy(&gn->y);
324: VecDestroy(&gn->y_work);
325: }
326: MatDestroy(&gn->H);
327: MatDestroy(&gn->D);
328: MatDestroy(&gn->Hreg);
329: TaoDestroy(&gn->subsolver);
330: gn->parent = NULL;
331: PetscFree(tao->data);
332: return(0);
333: }
335: /*MC
336: TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares
337: problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL
338: that constructs the Gauss-Newton problem with the user-provided least-squares
339: residual and Jacobian. The algorithm offers an L2-norm ("l2pure"), L2-norm proximal point ("l2prox")
340: regularizer, and L1-norm dictionary regularizer ("l1dict"), where we approximate the
341: L1-norm ||x||_1 by sum_i(sqrt(x_i^2+epsilon^2)-epsilon) with a small positive number epsilon.
342: The user can also provide own regularization function.
344: Options Database Keys:
345: + -tao_brgn_regularization_type - regularization type ("user", "l2prox", "l2pure", "l1dict") (default "l2prox")
346: . -tao_brgn_regularizer_weight - regularizer weight (default 1e-4)
347: - -tao_brgn_l1_smooth_epsilon - L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)
349: Level: beginner
350: M*/
351: PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao)
352: {
353: TAO_BRGN *gn;
355:
357: PetscNewLog(tao,&gn);
358:
359: tao->ops->destroy = TaoDestroy_BRGN;
360: tao->ops->setup = TaoSetUp_BRGN;
361: tao->ops->setfromoptions = TaoSetFromOptions_BRGN;
362: tao->ops->view = TaoView_BRGN;
363: tao->ops->solve = TaoSolve_BRGN;
364:
365: tao->data = (void*)gn;
366: gn->reg_type = BRGN_REGULARIZATION_L2PROX;
367: gn->lambda = 1e-4;
368: gn->epsilon = 1e-6;
369: gn->parent = tao;
370:
371: MatCreate(PetscObjectComm((PetscObject)tao),&gn->H);
372: MatSetOptionsPrefix(gn->H,"tao_brgn_hessian_");
373:
374: TaoCreate(PetscObjectComm((PetscObject)tao),&gn->subsolver);
375: TaoSetType(gn->subsolver,TAOBNLS);
376: TaoSetOptionsPrefix(gn->subsolver,"tao_brgn_subsolver_");
377: return(0);
378: }
380: /*@
381: TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN
383: Collective on Tao
385: Level: advanced
386:
387: Input Parameters:
388: + tao - the Tao solver context
389: - subsolver - the Tao sub-solver context
390: @*/
391: PetscErrorCode TaoBRGNGetSubsolver(Tao tao,Tao *subsolver)
392: {
393: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
394:
396: *subsolver = gn->subsolver;
397: return(0);
398: }
400: /*@
401: TaoBRGNSetRegularizerWeight - Set the regularizer weight for the Gauss-Newton least-squares algorithm
403: Collective on Tao
404:
405: Input Parameters:
406: + tao - the Tao solver context
407: - lambda - L1-norm regularizer weight
409: Level: beginner
410: @*/
411: PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao,PetscReal lambda)
412: {
413: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
414:
415: /* Initialize lambda here */
418: gn->lambda = lambda;
419: return(0);
420: }
422: /*@
423: TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm
425: Collective on Tao
426:
427: Input Parameters:
428: + tao - the Tao solver context
429: - epsilon - L1-norm smooth approximation parameter
431: Level: advanced
432: @*/
433: PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao,PetscReal epsilon)
434: {
435: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
436:
437: /* Initialize epsilon here */
440: gn->epsilon = epsilon;
441: return(0);
442: }
444: /*@
445: TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)
447: Input Parameters:
448: + tao - the Tao context
449: - dict - the user specified dictionary matrix. We allow to set a null dictionary, which means identity matrix by default
451: Level: advanced
452: @*/
453: PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao,Mat dict)
454: {
455: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
459: if (dict) {
462: PetscObjectReference((PetscObject)dict);
463: }
464: MatDestroy(&gn->D);
465: gn->D = dict;
466: return(0);
467: }
469: /*@C
470: TaoBRGNSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back
471: function into the algorithm.
473: Input Parameters:
474: + tao - the Tao context
475: . func - function pointer for the regularizer value and gradient evaluation
476: - ctx - user context for the regularizer
478: Level: advanced
479: @*/
480: PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao,PetscErrorCode (*func)(Tao,Vec,PetscReal *,Vec,void*),void *ctx)
481: {
482: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
486: if (ctx) {
487: gn->reg_obj_ctx = ctx;
488: }
489: if (func) {
490: gn->regularizerobjandgrad = func;
491: }
492: return(0);
493: }
495: /*@C
496: TaoBRGNSetRegularizerHessianRoutine - Sets the user-defined regularizer call-back
497: function into the algorithm.
499: Input Parameters:
500: + tao - the Tao context
501: . Hreg - user-created matrix for the Hessian of the regularization term
502: . func - function pointer for the regularizer Hessian evaluation
503: - ctx - user context for the regularizer Hessian
505: Level: advanced
506: @*/
507: PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao,Mat Hreg,PetscErrorCode (*func)(Tao,Vec,Mat,void*),void *ctx)
508: {
509: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
514: if (Hreg) {
517: } else SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONG,"NULL Hessian detected! User must provide valid Hessian for the regularizer.");
518: if (ctx) {
519: gn->reg_hess_ctx = ctx;
520: }
521: if (func) {
522: gn->regularizerhessian = func;
523: }
524: if (Hreg) {
525: PetscObjectReference((PetscObject)Hreg);
526: MatDestroy(&gn->Hreg);
527: gn->Hreg = Hreg;
528: }
529: return(0);
530: }