Actual source code: brgn.c
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_LM 4
8: #define BRGN_REGULARIZATION_TYPES 5
10: static const char *BRGN_REGULARIZATION_TABLE[64] = {"user","l2prox","l2pure","l1dict","lm"};
12: static PetscErrorCode GNHessianProd(Mat H,Vec in,Vec out)
13: {
14: TAO_BRGN *gn;
15: PetscErrorCode ierr;
18: MatShellGetContext(H,&gn);
19: MatMult(gn->subsolver->ls_jac,in,gn->r_work);
20: MatMultTranspose(gn->subsolver->ls_jac,gn->r_work,out);
21: switch (gn->reg_type) {
22: case BRGN_REGULARIZATION_USER:
23: MatMult(gn->Hreg,in,gn->x_work);
24: VecAXPY(out,gn->lambda,gn->x_work);
25: break;
26: case BRGN_REGULARIZATION_L2PURE:
27: VecAXPY(out,gn->lambda,in);
28: break;
29: case BRGN_REGULARIZATION_L2PROX:
30: VecAXPY(out,gn->lambda,in);
31: break;
32: case BRGN_REGULARIZATION_L1DICT:
33: /* out = out + lambda*D'*(diag.*(D*in)) */
34: if (gn->D) {
35: MatMult(gn->D,in,gn->y);/* y = D*in */
36: } else {
37: VecCopy(in,gn->y);
38: }
39: VecPointwiseMult(gn->y_work,gn->diag,gn->y); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */
40: if (gn->D) {
41: MatMultTranspose(gn->D,gn->y_work,gn->x_work); /* x_work = D'*(diag.*(D*in)) */
42: } else {
43: VecCopy(gn->y_work,gn->x_work);
44: }
45: VecAXPY(out,gn->lambda,gn->x_work);
46: break;
47: case BRGN_REGULARIZATION_LM:
48: VecPointwiseMult(gn->x_work,gn->damping,in);
49: VecAXPY(out,1,gn->x_work);
50: break;
51: }
52: return(0);
53: }
54: static PetscErrorCode ComputeDamping(TAO_BRGN *gn)
55: {
56: const PetscScalar *diag_ary;
57: PetscScalar *damping_ary;
58: PetscInt i,n;
59: PetscErrorCode ierr;
62: /* update damping */
63: VecGetArray(gn->damping,&damping_ary);
64: VecGetArrayRead(gn->diag,&diag_ary);
65: VecGetLocalSize(gn->damping,&n);
66: for (i=0; i<n; i++) {
67: damping_ary[i] = PetscClipInterval(diag_ary[i],PETSC_SQRT_MACHINE_EPSILON,PetscSqrtReal(PETSC_MAX_REAL));
68: }
69: VecScale(gn->damping,gn->lambda);
70: VecRestoreArray(gn->damping,&damping_ary);
71: VecRestoreArrayRead(gn->diag,&diag_ary);
72: return(0);
73: }
75: PetscErrorCode TaoBRGNGetDampingVector(Tao tao,Vec *d)
76: {
77: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
80: if (gn->reg_type != BRGN_REGULARIZATION_LM) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_SUP,"Damping vector is only available if regularization type is lm.");
81: *d = gn->damping;
82: return(0);
83: }
85: static PetscErrorCode GNObjectiveGradientEval(Tao tao,Vec X,PetscReal *fcn,Vec G,void *ptr)
86: {
87: TAO_BRGN *gn = (TAO_BRGN *)ptr;
88: PetscInt K; /* dimension of D*X */
89: PetscScalar yESum;
90: PetscErrorCode ierr;
91: PetscReal f_reg;
94: /* compute objective *fcn*/
95: /* compute first term 0.5*||ls_res||_2^2 */
96: TaoComputeResidual(tao,X,tao->ls_res);
97: VecDot(tao->ls_res,tao->ls_res,fcn);
98: *fcn *= 0.5;
99: /* compute gradient G */
100: TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);
101: MatMultTranspose(tao->ls_jac,tao->ls_res,G);
102: /* add the regularization contribution */
103: switch (gn->reg_type) {
104: case BRGN_REGULARIZATION_USER:
105: (*gn->regularizerobjandgrad)(tao,X,&f_reg,gn->x_work,gn->reg_obj_ctx);
106: *fcn += gn->lambda*f_reg;
107: VecAXPY(G,gn->lambda,gn->x_work);
108: break;
109: case BRGN_REGULARIZATION_L2PURE:
110: /* compute f = f + lambda*0.5*xk'*xk */
111: VecDot(X,X,&f_reg);
112: *fcn += gn->lambda*0.5*f_reg;
113: /* compute G = G + lambda*xk */
114: VecAXPY(G,gn->lambda,X);
115: break;
116: case BRGN_REGULARIZATION_L2PROX:
117: /* compute f = f + lambda*0.5*(xk - xkm1)'*(xk - xkm1) */
118: VecAXPBYPCZ(gn->x_work,1.0,-1.0,0.0,X,gn->x_old);
119: VecDot(gn->x_work,gn->x_work,&f_reg);
120: *fcn += gn->lambda*0.5*f_reg;
121: /* compute G = G + lambda*(xk - xkm1) */
122: VecAXPBYPCZ(G,gn->lambda,-gn->lambda,1.0,X,gn->x_old);
123: break;
124: case BRGN_REGULARIZATION_L1DICT:
125: /* compute f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/
126: if (gn->D) {
127: MatMult(gn->D,X,gn->y);/* y = D*x */
128: } else {
129: VecCopy(X,gn->y);
130: }
131: VecPointwiseMult(gn->y_work,gn->y,gn->y);
132: VecShift(gn->y_work,gn->epsilon*gn->epsilon);
133: VecSqrtAbs(gn->y_work); /* gn->y_work = sqrt(y.^2+epsilon^2) */
134: VecSum(gn->y_work,&yESum);
135: VecGetSize(gn->y,&K);
136: *fcn += gn->lambda*(yESum - K*gn->epsilon);
137: /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */
138: VecPointwiseDivide(gn->y_work,gn->y,gn->y_work); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */
139: if (gn->D) {
140: MatMultTranspose(gn->D,gn->y_work,gn->x_work);
141: } else {
142: VecCopy(gn->y_work,gn->x_work);
143: }
144: VecAXPY(G,gn->lambda,gn->x_work);
145: break;
146: }
147: return(0);
148: }
150: static PetscErrorCode GNComputeHessian(Tao tao,Vec X,Mat H,Mat Hpre,void *ptr)
151: {
152: TAO_BRGN *gn = (TAO_BRGN *)ptr;
153: PetscInt i,n,cstart,cend;
154: PetscScalar *cnorms,*diag_ary;
158: TaoComputeResidualJacobian(tao,X,tao->ls_jac,tao->ls_jac_pre);
159: if (gn->mat_explicit) {
160: MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_REUSE_MATRIX, PETSC_DEFAULT, &gn->H);
161: }
163: switch (gn->reg_type) {
164: case BRGN_REGULARIZATION_USER:
165: (*gn->regularizerhessian)(tao,X,gn->Hreg,gn->reg_hess_ctx);
166: if (gn->mat_explicit) {
167: MatAXPY(gn->H, 1.0, gn->Hreg, DIFFERENT_NONZERO_PATTERN);
168: }
169: break;
170: case BRGN_REGULARIZATION_L2PURE:
171: if (gn->mat_explicit) {
172: MatShift(gn->H, gn->lambda);
173: }
174: break;
175: case BRGN_REGULARIZATION_L2PROX:
176: if (gn->mat_explicit) {
177: MatShift(gn->H, gn->lambda);
178: }
179: break;
180: case BRGN_REGULARIZATION_L1DICT:
181: /* 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 */
182: if (gn->D) {
183: MatMult(gn->D,X,gn->y);/* y = D*x */
184: } else {
185: VecCopy(X,gn->y);
186: }
187: VecPointwiseMult(gn->y_work,gn->y,gn->y);
188: VecShift(gn->y_work,gn->epsilon*gn->epsilon);
189: VecCopy(gn->y_work,gn->diag); /* gn->diag = y.^2+epsilon^2 */
190: VecSqrtAbs(gn->y_work); /* gn->y_work = sqrt(y.^2+epsilon^2) */
191: VecPointwiseMult(gn->diag,gn->y_work,gn->diag);/* gn->diag = sqrt(y.^2+epsilon^2).^3 */
192: VecReciprocal(gn->diag);
193: VecScale(gn->diag,gn->epsilon*gn->epsilon);
194: if (gn->mat_explicit) {
195: MatDiagonalSet(gn->H, gn->diag, ADD_VALUES);
196: }
197: break;
198: case BRGN_REGULARIZATION_LM:
199: /* compute diagonal of J^T J */
200: MatGetSize(gn->parent->ls_jac,NULL,&n);
201: PetscMalloc1(n,&cnorms);
202: MatGetColumnNorms(gn->parent->ls_jac,NORM_2,cnorms);
203: MatGetOwnershipRangeColumn(gn->parent->ls_jac,&cstart,&cend);
204: VecGetArray(gn->diag,&diag_ary);
205: for (i = 0; i < cend-cstart; i++) {
206: diag_ary[i] = cnorms[cstart+i] * cnorms[cstart+i];
207: }
208: VecRestoreArray(gn->diag,&diag_ary);
209: PetscFree(cnorms);
210: ComputeDamping(gn);
211: if (gn->mat_explicit) {
212: MatDiagonalSet(gn->H, gn->damping, ADD_VALUES);
213: }
214: break;
215: }
216: return(0);
217: }
219: static PetscErrorCode GNHookFunction(Tao tao,PetscInt iter, void *ctx)
220: {
221: TAO_BRGN *gn = (TAO_BRGN *)ctx;
222: PetscErrorCode ierr;
225: /* Update basic tao information from the subsolver */
226: gn->parent->nfuncs = tao->nfuncs;
227: gn->parent->ngrads = tao->ngrads;
228: gn->parent->nfuncgrads = tao->nfuncgrads;
229: gn->parent->nhess = tao->nhess;
230: gn->parent->niter = tao->niter;
231: gn->parent->ksp_its = tao->ksp_its;
232: gn->parent->ksp_tot_its = tao->ksp_tot_its;
233: gn->parent->fc = tao->fc;
234: TaoGetConvergedReason(tao,&gn->parent->reason);
235: /* Update the solution vectors */
236: if (iter == 0) {
237: VecSet(gn->x_old,0.0);
238: } else {
239: VecCopy(tao->solution,gn->x_old);
240: VecCopy(tao->solution,gn->parent->solution);
241: }
242: /* Update the gradient */
243: VecCopy(tao->gradient,gn->parent->gradient);
245: /* Update damping parameter for LM */
246: if (gn->reg_type == BRGN_REGULARIZATION_LM) {
247: if (iter > 0) {
248: if (gn->fc_old > tao->fc) {
249: gn->lambda = gn->lambda * gn->downhill_lambda_change;
250: } else {
251: /* uphill step */
252: gn->lambda = gn->lambda * gn->uphill_lambda_change;
253: }
254: }
255: gn->fc_old = tao->fc;
256: }
258: /* Call general purpose update function */
259: if (gn->parent->ops->update) {
260: (*gn->parent->ops->update)(gn->parent,gn->parent->niter,gn->parent->user_update);
261: }
262: return(0);
263: }
265: static PetscErrorCode TaoSolve_BRGN(Tao tao)
266: {
267: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
268: PetscErrorCode ierr;
271: TaoSolve(gn->subsolver);
272: /* Update basic tao information from the subsolver */
273: tao->nfuncs = gn->subsolver->nfuncs;
274: tao->ngrads = gn->subsolver->ngrads;
275: tao->nfuncgrads = gn->subsolver->nfuncgrads;
276: tao->nhess = gn->subsolver->nhess;
277: tao->niter = gn->subsolver->niter;
278: tao->ksp_its = gn->subsolver->ksp_its;
279: tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
280: TaoGetConvergedReason(gn->subsolver,&tao->reason);
281: /* Update vectors */
282: VecCopy(gn->subsolver->solution,tao->solution);
283: VecCopy(gn->subsolver->gradient,tao->gradient);
284: return(0);
285: }
287: static PetscErrorCode TaoSetFromOptions_BRGN(PetscOptionItems *PetscOptionsObject,Tao tao)
288: {
289: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
290: TaoLineSearch ls;
291: PetscErrorCode ierr;
294: 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.");
295: PetscOptionsBool("-tao_brgn_mat_explicit","switches the Hessian construction to be an explicit matrix rather than MATSHELL","",gn->mat_explicit,&gn->mat_explicit,NULL);
296: PetscOptionsReal("-tao_brgn_regularizer_weight","regularizer weight (default 1e-4)","",gn->lambda,&gn->lambda,NULL);
297: 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);
298: PetscOptionsReal("-tao_brgn_lm_downhill_lambda_change","Factor to decrease trust region by on downhill steps","",gn->downhill_lambda_change,&gn->downhill_lambda_change,NULL);
299: PetscOptionsReal("-tao_brgn_lm_uphill_lambda_change","Factor to increase trust region by on uphill steps","",gn->uphill_lambda_change,&gn->uphill_lambda_change,NULL);
300: PetscOptionsEList("-tao_brgn_regularization_type","regularization type", "",BRGN_REGULARIZATION_TABLE,BRGN_REGULARIZATION_TYPES,BRGN_REGULARIZATION_TABLE[gn->reg_type],&gn->reg_type,NULL);
301: PetscOptionsTail();
302: /* set unit line search direction as the default when using the lm regularizer */
303: if (gn->reg_type == BRGN_REGULARIZATION_LM) {
304: TaoGetLineSearch(gn->subsolver,&ls);
305: TaoLineSearchSetType(ls,TAOLINESEARCHUNIT);
306: }
307: TaoSetFromOptions(gn->subsolver);
308: return(0);
309: }
311: static PetscErrorCode TaoView_BRGN(Tao tao,PetscViewer viewer)
312: {
313: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
314: PetscErrorCode ierr;
317: PetscViewerASCIIPushTab(viewer);
318: TaoView(gn->subsolver,viewer);
319: PetscViewerASCIIPopTab(viewer);
320: return(0);
321: }
323: static PetscErrorCode TaoSetUp_BRGN(Tao tao)
324: {
325: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
326: PetscErrorCode ierr;
327: PetscBool is_bnls,is_bntr,is_bntl;
328: PetscInt i,n,N,K; /* dict has size K*N*/
331: if (!tao->ls_res) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualRoutine() must be called before setup!");
332: PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNLS,&is_bnls);
333: PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTR,&is_bntr);
334: PetscObjectTypeCompare((PetscObject)gn->subsolver,TAOBNTL,&is_bntl);
335: if ((is_bnls || is_bntr || is_bntl) && !tao->ls_jac) SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ORDER,"TaoSetResidualJacobianRoutine() must be called before setup!");
336: if (!tao->gradient) {
337: VecDuplicate(tao->solution,&tao->gradient);
338: }
339: if (!gn->x_work) {
340: VecDuplicate(tao->solution,&gn->x_work);
341: }
342: if (!gn->r_work) {
343: VecDuplicate(tao->ls_res,&gn->r_work);
344: }
345: if (!gn->x_old) {
346: VecDuplicate(tao->solution,&gn->x_old);
347: VecSet(gn->x_old,0.0);
348: }
350: if (BRGN_REGULARIZATION_L1DICT == gn->reg_type) {
351: if (!gn->y) {
352: if (gn->D) {
353: 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 */
354: MatCreateVecs(gn->D,NULL,&gn->y);
355: } else {
356: VecDuplicate(tao->solution,&gn->y); /* If user does not setup dict matrix, use identiy matrix, K=N */
357: }
358: VecSet(gn->y,0.0);
359: }
360: if (!gn->y_work) {
361: VecDuplicate(gn->y,&gn->y_work);
362: }
363: if (!gn->diag) {
364: VecDuplicate(gn->y,&gn->diag);
365: VecSet(gn->diag,0.0);
366: }
367: }
368: if (BRGN_REGULARIZATION_LM == gn->reg_type) {
369: if (!gn->diag) {
370: MatCreateVecs(tao->ls_jac,&gn->diag,NULL);
371: }
372: if (!gn->damping) {
373: MatCreateVecs(tao->ls_jac,&gn->damping,NULL);
374: }
375: }
377: if (!tao->setupcalled) {
378: /* Hessian setup */
379: if (gn->mat_explicit) {
380: TaoComputeResidualJacobian(tao,tao->solution,tao->ls_jac,tao->ls_jac_pre);
381: MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &gn->H);
382: } else {
383: VecGetLocalSize(tao->solution,&n);
384: VecGetSize(tao->solution,&N);
385: MatCreate(PetscObjectComm((PetscObject)tao),&gn->H);
386: MatSetSizes(gn->H,n,n,N,N);
387: MatSetType(gn->H,MATSHELL);
388: MatSetOption(gn->H, MAT_SYMMETRIC, PETSC_TRUE);
389: MatShellSetOperation(gn->H,MATOP_MULT,(void (*)(void))GNHessianProd);
390: MatShellSetContext(gn->H,gn);
391: }
392: MatSetUp(gn->H);
393: /* Subsolver setup,include initial vector and dictionary D */
394: TaoSetUpdate(gn->subsolver,GNHookFunction,gn);
395: TaoSetInitialVector(gn->subsolver,tao->solution);
396: if (tao->bounded) {
397: TaoSetVariableBounds(gn->subsolver,tao->XL,tao->XU);
398: }
399: TaoSetResidualRoutine(gn->subsolver,tao->ls_res,tao->ops->computeresidual,tao->user_lsresP);
400: TaoSetJacobianResidualRoutine(gn->subsolver,tao->ls_jac,tao->ls_jac,tao->ops->computeresidualjacobian,tao->user_lsjacP);
401: TaoSetObjectiveAndGradientRoutine(gn->subsolver,GNObjectiveGradientEval,gn);
402: TaoSetHessianRoutine(gn->subsolver,gn->H,gn->H,GNComputeHessian,gn);
403: /* Propagate some options down */
404: TaoSetTolerances(gn->subsolver,tao->gatol,tao->grtol,tao->gttol);
405: TaoSetMaximumIterations(gn->subsolver,tao->max_it);
406: TaoSetMaximumFunctionEvaluations(gn->subsolver,tao->max_funcs);
407: for (i=0; i<tao->numbermonitors; ++i) {
408: TaoSetMonitor(gn->subsolver,tao->monitor[i],tao->monitorcontext[i],tao->monitordestroy[i]);
409: PetscObjectReference((PetscObject)(tao->monitorcontext[i]));
410: }
411: TaoSetUp(gn->subsolver);
412: }
413: return(0);
414: }
416: static PetscErrorCode TaoDestroy_BRGN(Tao tao)
417: {
418: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
419: PetscErrorCode ierr;
422: if (tao->setupcalled) {
423: VecDestroy(&tao->gradient);
424: VecDestroy(&gn->x_work);
425: VecDestroy(&gn->r_work);
426: VecDestroy(&gn->x_old);
427: VecDestroy(&gn->diag);
428: VecDestroy(&gn->y);
429: VecDestroy(&gn->y_work);
430: }
431: VecDestroy(&gn->damping);
432: VecDestroy(&gn->diag);
433: MatDestroy(&gn->H);
434: MatDestroy(&gn->D);
435: MatDestroy(&gn->Hreg);
436: TaoDestroy(&gn->subsolver);
437: gn->parent = NULL;
438: PetscFree(tao->data);
439: return(0);
440: }
442: /*MC
443: TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares
444: problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL
445: that constructs the Gauss-Newton problem with the user-provided least-squares
446: residual and Jacobian. The algorithm offers an L2-norm ("l2pure"), L2-norm proximal point ("l2prox")
447: regularizer, and L1-norm dictionary regularizer ("l1dict"), where we approximate the
448: L1-norm ||x||_1 by sum_i(sqrt(x_i^2+epsilon^2)-epsilon) with a small positive number epsilon.
449: Also offered is the "lm" regularizer which uses a scaled diagonal of J^T J.
450: With the "lm" regularizer, BRGN is a Levenberg-Marquardt optimizer.
451: The user can also provide own regularization function.
453: Options Database Keys:
454: + -tao_brgn_regularization_type - regularization type ("user", "l2prox", "l2pure", "l1dict", "lm") (default "l2prox")
455: . -tao_brgn_regularizer_weight - regularizer weight (default 1e-4)
456: - -tao_brgn_l1_smooth_epsilon - L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)
458: Level: beginner
459: M*/
460: PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao)
461: {
462: TAO_BRGN *gn;
466: PetscNewLog(tao,&gn);
468: tao->ops->destroy = TaoDestroy_BRGN;
469: tao->ops->setup = TaoSetUp_BRGN;
470: tao->ops->setfromoptions = TaoSetFromOptions_BRGN;
471: tao->ops->view = TaoView_BRGN;
472: tao->ops->solve = TaoSolve_BRGN;
474: tao->data = gn;
475: gn->reg_type = BRGN_REGULARIZATION_L2PROX;
476: gn->lambda = 1e-4;
477: gn->epsilon = 1e-6;
478: gn->downhill_lambda_change = 1./5.;
479: gn->uphill_lambda_change = 1.5;
480: gn->parent = tao;
482: TaoCreate(PetscObjectComm((PetscObject)tao),&gn->subsolver);
483: TaoSetType(gn->subsolver,TAOBNLS);
484: TaoSetOptionsPrefix(gn->subsolver,"tao_brgn_subsolver_");
485: return(0);
486: }
488: /*@
489: TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN
491: Collective on Tao
493: Level: advanced
495: Input Parameters:
496: + tao - the Tao solver context
497: - subsolver - the Tao sub-solver context
498: @*/
499: PetscErrorCode TaoBRGNGetSubsolver(Tao tao,Tao *subsolver)
500: {
501: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
504: *subsolver = gn->subsolver;
505: return(0);
506: }
508: /*@
509: TaoBRGNSetRegularizerWeight - Set the regularizer weight for the Gauss-Newton least-squares algorithm
511: Collective on Tao
513: Input Parameters:
514: + tao - the Tao solver context
515: - lambda - L1-norm regularizer weight
517: Level: beginner
518: @*/
519: PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao,PetscReal lambda)
520: {
521: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
523: /* Initialize lambda here */
526: gn->lambda = lambda;
527: return(0);
528: }
530: /*@
531: TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm
533: Collective on Tao
535: Input Parameters:
536: + tao - the Tao solver context
537: - epsilon - L1-norm smooth approximation parameter
539: Level: advanced
540: @*/
541: PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao,PetscReal epsilon)
542: {
543: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
545: /* Initialize epsilon here */
548: gn->epsilon = epsilon;
549: return(0);
550: }
552: /*@
553: TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)
555: Input Parameters:
556: + tao - the Tao context
557: - dict - the user specified dictionary matrix. We allow to set a null dictionary, which means identity matrix by default
559: Level: advanced
560: @*/
561: PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao,Mat dict)
562: {
563: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
567: if (dict) {
570: PetscObjectReference((PetscObject)dict);
571: }
572: MatDestroy(&gn->D);
573: gn->D = dict;
574: return(0);
575: }
577: /*@C
578: TaoBRGNSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back
579: function into the algorithm.
581: Input Parameters:
582: + tao - the Tao context
583: . func - function pointer for the regularizer value and gradient evaluation
584: - ctx - user context for the regularizer
586: Level: advanced
587: @*/
588: PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao,PetscErrorCode (*func)(Tao,Vec,PetscReal *,Vec,void*),void *ctx)
589: {
590: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
594: if (ctx) {
595: gn->reg_obj_ctx = ctx;
596: }
597: if (func) {
598: gn->regularizerobjandgrad = func;
599: }
600: return(0);
601: }
603: /*@C
604: TaoBRGNSetRegularizerHessianRoutine - Sets the user-defined regularizer call-back
605: function into the algorithm.
607: Input Parameters:
608: + tao - the Tao context
609: . Hreg - user-created matrix for the Hessian of the regularization term
610: . func - function pointer for the regularizer Hessian evaluation
611: - ctx - user context for the regularizer Hessian
613: Level: advanced
614: @*/
615: PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao,Mat Hreg,PetscErrorCode (*func)(Tao,Vec,Mat,void*),void *ctx)
616: {
617: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
622: if (Hreg) {
625: } else SETERRQ(PetscObjectComm((PetscObject)tao),PETSC_ERR_ARG_WRONG,"NULL Hessian detected! User must provide valid Hessian for the regularizer.");
626: if (ctx) {
627: gn->reg_hess_ctx = ctx;
628: }
629: if (func) {
630: gn->regularizerhessian = func;
631: }
632: if (Hreg) {
633: PetscObjectReference((PetscObject)Hreg);
634: MatDestroy(&gn->Hreg);
635: gn->Hreg = Hreg;
636: }
637: return(0);
638: }