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;
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_L2PURE:
25: VecAXPY(out, gn->lambda, in);
26: break;
27: case BRGN_REGULARIZATION_L2PROX:
28: VecAXPY(out, gn->lambda, in);
29: break;
30: case BRGN_REGULARIZATION_L1DICT:
31: /* out = out + lambda*D'*(diag.*(D*in)) */
32: if (gn->D) {
33: MatMult(gn->D, in, gn->y); /* y = D*in */
34: } else {
35: VecCopy(in, gn->y);
36: }
37: VecPointwiseMult(gn->y_work, gn->diag, gn->y); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */
38: if (gn->D) {
39: MatMultTranspose(gn->D, gn->y_work, gn->x_work); /* x_work = D'*(diag.*(D*in)) */
40: } else {
41: VecCopy(gn->y_work, gn->x_work);
42: }
43: VecAXPY(out, gn->lambda, gn->x_work);
44: break;
45: case BRGN_REGULARIZATION_LM:
46: VecPointwiseMult(gn->x_work, gn->damping, in);
47: VecAXPY(out, 1, gn->x_work);
48: break;
49: }
50: return 0;
51: }
52: static PetscErrorCode ComputeDamping(TAO_BRGN *gn)
53: {
54: const PetscScalar *diag_ary;
55: PetscScalar *damping_ary;
56: PetscInt i, n;
58: /* update damping */
59: VecGetArray(gn->damping, &damping_ary);
60: VecGetArrayRead(gn->diag, &diag_ary);
61: VecGetLocalSize(gn->damping, &n);
62: for (i = 0; i < n; i++) damping_ary[i] = PetscClipInterval(diag_ary[i], PETSC_SQRT_MACHINE_EPSILON, PetscSqrtReal(PETSC_MAX_REAL));
63: VecScale(gn->damping, gn->lambda);
64: VecRestoreArray(gn->damping, &damping_ary);
65: VecRestoreArrayRead(gn->diag, &diag_ary);
66: return 0;
67: }
69: PetscErrorCode TaoBRGNGetDampingVector(Tao tao, Vec *d)
70: {
71: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
74: *d = gn->damping;
75: return 0;
76: }
78: static PetscErrorCode GNObjectiveGradientEval(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr)
79: {
80: TAO_BRGN *gn = (TAO_BRGN *)ptr;
81: PetscInt K; /* dimension of D*X */
82: PetscScalar yESum;
83: PetscReal f_reg;
85: /* compute objective *fcn*/
86: /* compute first term 0.5*||ls_res||_2^2 */
87: TaoComputeResidual(tao, X, tao->ls_res);
88: VecDot(tao->ls_res, tao->ls_res, fcn);
89: *fcn *= 0.5;
90: /* compute gradient G */
91: TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);
92: MatMultTranspose(tao->ls_jac, tao->ls_res, G);
93: /* add the regularization contribution */
94: switch (gn->reg_type) {
95: case BRGN_REGULARIZATION_USER:
96: (*gn->regularizerobjandgrad)(tao, X, &f_reg, gn->x_work, gn->reg_obj_ctx);
97: *fcn += gn->lambda * f_reg;
98: VecAXPY(G, gn->lambda, gn->x_work);
99: break;
100: case BRGN_REGULARIZATION_L2PURE:
101: /* compute f = f + lambda*0.5*xk'*xk */
102: VecDot(X, X, &f_reg);
103: *fcn += gn->lambda * 0.5 * f_reg;
104: /* compute G = G + lambda*xk */
105: VecAXPY(G, gn->lambda, X);
106: break;
107: case BRGN_REGULARIZATION_L2PROX:
108: /* compute f = f + lambda*0.5*(xk - xkm1)'*(xk - xkm1) */
109: VecAXPBYPCZ(gn->x_work, 1.0, -1.0, 0.0, X, gn->x_old);
110: VecDot(gn->x_work, gn->x_work, &f_reg);
111: *fcn += gn->lambda * 0.5 * f_reg;
112: /* compute G = G + lambda*(xk - xkm1) */
113: VecAXPBYPCZ(G, gn->lambda, -gn->lambda, 1.0, X, gn->x_old);
114: break;
115: case BRGN_REGULARIZATION_L1DICT:
116: /* compute f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/
117: if (gn->D) {
118: MatMult(gn->D, X, gn->y); /* y = D*x */
119: } else {
120: VecCopy(X, gn->y);
121: }
122: VecPointwiseMult(gn->y_work, gn->y, gn->y);
123: VecShift(gn->y_work, gn->epsilon * gn->epsilon);
124: VecSqrtAbs(gn->y_work); /* gn->y_work = sqrt(y.^2+epsilon^2) */
125: VecSum(gn->y_work, &yESum);
126: VecGetSize(gn->y, &K);
127: *fcn += gn->lambda * (yESum - K * gn->epsilon);
128: /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */
129: VecPointwiseDivide(gn->y_work, gn->y, gn->y_work); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */
130: if (gn->D) {
131: MatMultTranspose(gn->D, gn->y_work, gn->x_work);
132: } else {
133: VecCopy(gn->y_work, gn->x_work);
134: }
135: VecAXPY(G, gn->lambda, gn->x_work);
136: break;
137: }
138: return 0;
139: }
141: static PetscErrorCode GNComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
142: {
143: TAO_BRGN *gn = (TAO_BRGN *)ptr;
144: PetscInt i, n, cstart, cend;
145: PetscScalar *cnorms, *diag_ary;
147: TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre);
148: if (gn->mat_explicit) MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_REUSE_MATRIX, PETSC_DEFAULT, &gn->H);
150: switch (gn->reg_type) {
151: case BRGN_REGULARIZATION_USER:
152: (*gn->regularizerhessian)(tao, X, gn->Hreg, gn->reg_hess_ctx);
153: if (gn->mat_explicit) MatAXPY(gn->H, 1.0, gn->Hreg, DIFFERENT_NONZERO_PATTERN);
154: break;
155: case BRGN_REGULARIZATION_L2PURE:
156: if (gn->mat_explicit) MatShift(gn->H, gn->lambda);
157: break;
158: case BRGN_REGULARIZATION_L2PROX:
159: if (gn->mat_explicit) MatShift(gn->H, gn->lambda);
160: break;
161: case BRGN_REGULARIZATION_L1DICT:
162: /* 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 */
163: if (gn->D) {
164: MatMult(gn->D, X, gn->y); /* y = D*x */
165: } else {
166: VecCopy(X, gn->y);
167: }
168: VecPointwiseMult(gn->y_work, gn->y, gn->y);
169: VecShift(gn->y_work, gn->epsilon * gn->epsilon);
170: VecCopy(gn->y_work, gn->diag); /* gn->diag = y.^2+epsilon^2 */
171: VecSqrtAbs(gn->y_work); /* gn->y_work = sqrt(y.^2+epsilon^2) */
172: VecPointwiseMult(gn->diag, gn->y_work, gn->diag); /* gn->diag = sqrt(y.^2+epsilon^2).^3 */
173: VecReciprocal(gn->diag);
174: VecScale(gn->diag, gn->epsilon * gn->epsilon);
175: if (gn->mat_explicit) MatDiagonalSet(gn->H, gn->diag, ADD_VALUES);
176: break;
177: case BRGN_REGULARIZATION_LM:
178: /* compute diagonal of J^T J */
179: MatGetSize(gn->parent->ls_jac, NULL, &n);
180: PetscMalloc1(n, &cnorms);
181: MatGetColumnNorms(gn->parent->ls_jac, NORM_2, cnorms);
182: MatGetOwnershipRangeColumn(gn->parent->ls_jac, &cstart, &cend);
183: VecGetArray(gn->diag, &diag_ary);
184: for (i = 0; i < cend - cstart; i++) diag_ary[i] = cnorms[cstart + i] * cnorms[cstart + i];
185: VecRestoreArray(gn->diag, &diag_ary);
186: PetscFree(cnorms);
187: ComputeDamping(gn);
188: if (gn->mat_explicit) MatDiagonalSet(gn->H, gn->damping, ADD_VALUES);
189: break;
190: }
191: return 0;
192: }
194: static PetscErrorCode GNHookFunction(Tao tao, PetscInt iter, void *ctx)
195: {
196: TAO_BRGN *gn = (TAO_BRGN *)ctx;
198: /* Update basic tao information from the subsolver */
199: gn->parent->nfuncs = tao->nfuncs;
200: gn->parent->ngrads = tao->ngrads;
201: gn->parent->nfuncgrads = tao->nfuncgrads;
202: gn->parent->nhess = tao->nhess;
203: gn->parent->niter = tao->niter;
204: gn->parent->ksp_its = tao->ksp_its;
205: gn->parent->ksp_tot_its = tao->ksp_tot_its;
206: gn->parent->fc = tao->fc;
207: TaoGetConvergedReason(tao, &gn->parent->reason);
208: /* Update the solution vectors */
209: if (iter == 0) {
210: VecSet(gn->x_old, 0.0);
211: } else {
212: VecCopy(tao->solution, gn->x_old);
213: VecCopy(tao->solution, gn->parent->solution);
214: }
215: /* Update the gradient */
216: VecCopy(tao->gradient, gn->parent->gradient);
218: /* Update damping parameter for LM */
219: if (gn->reg_type == BRGN_REGULARIZATION_LM) {
220: if (iter > 0) {
221: if (gn->fc_old > tao->fc) {
222: gn->lambda = gn->lambda * gn->downhill_lambda_change;
223: } else {
224: /* uphill step */
225: gn->lambda = gn->lambda * gn->uphill_lambda_change;
226: }
227: }
228: gn->fc_old = tao->fc;
229: }
231: /* Call general purpose update function */
232: if (gn->parent->ops->update) (*gn->parent->ops->update)(gn->parent, gn->parent->niter, gn->parent->user_update);
233: return 0;
234: }
236: static PetscErrorCode TaoSolve_BRGN(Tao tao)
237: {
238: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
240: TaoSolve(gn->subsolver);
241: /* Update basic tao information from the subsolver */
242: tao->nfuncs = gn->subsolver->nfuncs;
243: tao->ngrads = gn->subsolver->ngrads;
244: tao->nfuncgrads = gn->subsolver->nfuncgrads;
245: tao->nhess = gn->subsolver->nhess;
246: tao->niter = gn->subsolver->niter;
247: tao->ksp_its = gn->subsolver->ksp_its;
248: tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
249: TaoGetConvergedReason(gn->subsolver, &tao->reason);
250: /* Update vectors */
251: VecCopy(gn->subsolver->solution, tao->solution);
252: VecCopy(gn->subsolver->gradient, tao->gradient);
253: return 0;
254: }
256: static PetscErrorCode TaoSetFromOptions_BRGN(Tao tao, PetscOptionItems *PetscOptionsObject)
257: {
258: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
259: TaoLineSearch ls;
261: PetscOptionsHeadBegin(PetscOptionsObject, "least-squares problems with regularizer: ||f(x)||^2 + lambda*g(x), g(x) = ||xk-xkm1||^2 or ||Dx||_1 or user defined function.");
262: PetscOptionsBool("-tao_brgn_mat_explicit", "switches the Hessian construction to be an explicit matrix rather than MATSHELL", "", gn->mat_explicit, &gn->mat_explicit, NULL);
263: PetscOptionsReal("-tao_brgn_regularizer_weight", "regularizer weight (default 1e-4)", "", gn->lambda, &gn->lambda, NULL);
264: 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);
265: 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);
266: 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);
267: PetscOptionsEList("-tao_brgn_regularization_type", "regularization type", "", BRGN_REGULARIZATION_TABLE, BRGN_REGULARIZATION_TYPES, BRGN_REGULARIZATION_TABLE[gn->reg_type], &gn->reg_type, NULL);
268: PetscOptionsHeadEnd();
269: /* set unit line search direction as the default when using the lm regularizer */
270: if (gn->reg_type == BRGN_REGULARIZATION_LM) {
271: TaoGetLineSearch(gn->subsolver, &ls);
272: TaoLineSearchSetType(ls, TAOLINESEARCHUNIT);
273: }
274: TaoSetFromOptions(gn->subsolver);
275: return 0;
276: }
278: static PetscErrorCode TaoView_BRGN(Tao tao, PetscViewer viewer)
279: {
280: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
282: PetscViewerASCIIPushTab(viewer);
283: TaoView(gn->subsolver, viewer);
284: PetscViewerASCIIPopTab(viewer);
285: return 0;
286: }
288: static PetscErrorCode TaoSetUp_BRGN(Tao tao)
289: {
290: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
291: PetscBool is_bnls, is_bntr, is_bntl;
292: PetscInt i, n, N, K; /* dict has size K*N*/
295: PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNLS, &is_bnls);
296: PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTR, &is_bntr);
297: PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTL, &is_bntl);
299: if (!tao->gradient) VecDuplicate(tao->solution, &tao->gradient);
300: if (!gn->x_work) VecDuplicate(tao->solution, &gn->x_work);
301: if (!gn->r_work) VecDuplicate(tao->ls_res, &gn->r_work);
302: if (!gn->x_old) {
303: VecDuplicate(tao->solution, &gn->x_old);
304: VecSet(gn->x_old, 0.0);
305: }
307: if (BRGN_REGULARIZATION_L1DICT == gn->reg_type) {
308: if (!gn->y) {
309: if (gn->D) {
310: 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 */
311: MatCreateVecs(gn->D, NULL, &gn->y);
312: } else {
313: VecDuplicate(tao->solution, &gn->y); /* If user does not setup dict matrix, use identiy matrix, K=N */
314: }
315: VecSet(gn->y, 0.0);
316: }
317: if (!gn->y_work) VecDuplicate(gn->y, &gn->y_work);
318: if (!gn->diag) {
319: VecDuplicate(gn->y, &gn->diag);
320: VecSet(gn->diag, 0.0);
321: }
322: }
323: if (BRGN_REGULARIZATION_LM == gn->reg_type) {
324: if (!gn->diag) MatCreateVecs(tao->ls_jac, &gn->diag, NULL);
325: if (!gn->damping) MatCreateVecs(tao->ls_jac, &gn->damping, NULL);
326: }
328: if (!tao->setupcalled) {
329: /* Hessian setup */
330: if (gn->mat_explicit) {
331: TaoComputeResidualJacobian(tao, tao->solution, tao->ls_jac, tao->ls_jac_pre);
332: MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &gn->H);
333: } else {
334: VecGetLocalSize(tao->solution, &n);
335: VecGetSize(tao->solution, &N);
336: MatCreate(PetscObjectComm((PetscObject)tao), &gn->H);
337: MatSetSizes(gn->H, n, n, N, N);
338: MatSetType(gn->H, MATSHELL);
339: MatSetOption(gn->H, MAT_SYMMETRIC, PETSC_TRUE);
340: MatShellSetOperation(gn->H, MATOP_MULT, (void (*)(void))GNHessianProd);
341: MatShellSetContext(gn->H, gn);
342: }
343: MatSetUp(gn->H);
344: /* Subsolver setup,include initial vector and dictionary D */
345: TaoSetUpdate(gn->subsolver, GNHookFunction, gn);
346: TaoSetSolution(gn->subsolver, tao->solution);
347: if (tao->bounded) TaoSetVariableBounds(gn->subsolver, tao->XL, tao->XU);
348: TaoSetResidualRoutine(gn->subsolver, tao->ls_res, tao->ops->computeresidual, tao->user_lsresP);
349: TaoSetJacobianResidualRoutine(gn->subsolver, tao->ls_jac, tao->ls_jac, tao->ops->computeresidualjacobian, tao->user_lsjacP);
350: TaoSetObjectiveAndGradient(gn->subsolver, NULL, GNObjectiveGradientEval, gn);
351: TaoSetHessian(gn->subsolver, gn->H, gn->H, GNComputeHessian, gn);
352: /* Propagate some options down */
353: TaoSetTolerances(gn->subsolver, tao->gatol, tao->grtol, tao->gttol);
354: TaoSetMaximumIterations(gn->subsolver, tao->max_it);
355: TaoSetMaximumFunctionEvaluations(gn->subsolver, tao->max_funcs);
356: for (i = 0; i < tao->numbermonitors; ++i) {
357: TaoSetMonitor(gn->subsolver, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i]);
358: PetscObjectReference((PetscObject)(tao->monitorcontext[i]));
359: }
360: TaoSetUp(gn->subsolver);
361: }
362: return 0;
363: }
365: static PetscErrorCode TaoDestroy_BRGN(Tao tao)
366: {
367: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
369: if (tao->setupcalled) {
370: VecDestroy(&tao->gradient);
371: VecDestroy(&gn->x_work);
372: VecDestroy(&gn->r_work);
373: VecDestroy(&gn->x_old);
374: VecDestroy(&gn->diag);
375: VecDestroy(&gn->y);
376: VecDestroy(&gn->y_work);
377: }
378: VecDestroy(&gn->damping);
379: VecDestroy(&gn->diag);
380: MatDestroy(&gn->H);
381: MatDestroy(&gn->D);
382: MatDestroy(&gn->Hreg);
383: TaoDestroy(&gn->subsolver);
384: gn->parent = NULL;
385: PetscFree(tao->data);
386: return 0;
387: }
389: /*MC
390: TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares
391: problems with bound constraints. This algorithm is a thin wrapper around TAOBNTL
392: that constructs the Gauss-Newton problem with the user-provided least-squares
393: residual and Jacobian. The algorithm offers an L2-norm ("l2pure"), L2-norm proximal point ("l2prox")
394: regularizer, and L1-norm dictionary regularizer ("l1dict"), where we approximate the
395: L1-norm ||x||_1 by sum_i(sqrt(x_i^2+epsilon^2)-epsilon) with a small positive number epsilon.
396: Also offered is the "lm" regularizer which uses a scaled diagonal of J^T J.
397: With the "lm" regularizer, BRGN is a Levenberg-Marquardt optimizer.
398: The user can also provide own regularization function.
400: Options Database Keys:
401: + -tao_brgn_regularization_type - regularization type ("user", "l2prox", "l2pure", "l1dict", "lm") (default "l2prox")
402: . -tao_brgn_regularizer_weight - regularizer weight (default 1e-4)
403: - -tao_brgn_l1_smooth_epsilon - L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)
405: Level: beginner
406: M*/
407: PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao)
408: {
409: TAO_BRGN *gn;
411: PetscNew(&gn);
413: tao->ops->destroy = TaoDestroy_BRGN;
414: tao->ops->setup = TaoSetUp_BRGN;
415: tao->ops->setfromoptions = TaoSetFromOptions_BRGN;
416: tao->ops->view = TaoView_BRGN;
417: tao->ops->solve = TaoSolve_BRGN;
419: tao->data = gn;
420: gn->reg_type = BRGN_REGULARIZATION_L2PROX;
421: gn->lambda = 1e-4;
422: gn->epsilon = 1e-6;
423: gn->downhill_lambda_change = 1. / 5.;
424: gn->uphill_lambda_change = 1.5;
425: gn->parent = tao;
427: TaoCreate(PetscObjectComm((PetscObject)tao), &gn->subsolver);
428: TaoSetType(gn->subsolver, TAOBNLS);
429: TaoSetOptionsPrefix(gn->subsolver, "tao_brgn_subsolver_");
430: return 0;
431: }
433: /*@
434: TaoBRGNGetSubsolver - Get the pointer to the subsolver inside BRGN
436: Collective
438: Level: advanced
440: Input Parameters:
441: + tao - the Tao solver context
442: - subsolver - the Tao sub-solver context
443: @*/
444: PetscErrorCode TaoBRGNGetSubsolver(Tao tao, Tao *subsolver)
445: {
446: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
448: *subsolver = gn->subsolver;
449: return 0;
450: }
452: /*@
453: TaoBRGNSetRegularizerWeight - Set the regularizer weight for the Gauss-Newton least-squares algorithm
455: Collective
457: Input Parameters:
458: + tao - the Tao solver context
459: - lambda - L1-norm regularizer weight
461: Level: beginner
462: @*/
463: PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao, PetscReal lambda)
464: {
465: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
467: /* Initialize lambda here */
469: gn->lambda = lambda;
470: return 0;
471: }
473: /*@
474: TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm
476: Collective
478: Input Parameters:
479: + tao - the Tao solver context
480: - epsilon - L1-norm smooth approximation parameter
482: Level: advanced
483: @*/
484: PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao, PetscReal epsilon)
485: {
486: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
488: /* Initialize epsilon here */
490: gn->epsilon = epsilon;
491: return 0;
492: }
494: /*@
495: TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)
497: Input Parameters:
498: + tao - the Tao context
499: - dict - the user specified dictionary matrix. We allow to set a null dictionary, which means identity matrix by default
501: Level: advanced
502: @*/
503: PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao, Mat dict)
504: {
505: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
507: if (dict) {
510: PetscObjectReference((PetscObject)dict);
511: }
512: MatDestroy(&gn->D);
513: gn->D = dict;
514: return 0;
515: }
517: /*@C
518: TaoBRGNSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back
519: function into the algorithm.
521: Input Parameters:
522: + tao - the Tao context
523: . func - function pointer for the regularizer value and gradient evaluation
524: - ctx - user context for the regularizer
526: Level: advanced
527: @*/
528: PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao, PetscErrorCode (*func)(Tao, Vec, PetscReal *, Vec, void *), void *ctx)
529: {
530: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
533: if (ctx) gn->reg_obj_ctx = ctx;
534: if (func) gn->regularizerobjandgrad = func;
535: return 0;
536: }
538: /*@C
539: TaoBRGNSetRegularizerHessianRoutine - Sets the user-defined regularizer call-back
540: function into the algorithm.
542: Input Parameters:
543: + tao - the Tao context
544: . Hreg - user-created matrix for the Hessian of the regularization term
545: . func - function pointer for the regularizer Hessian evaluation
546: - ctx - user context for the regularizer Hessian
548: Level: advanced
549: @*/
550: PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao, Mat Hreg, PetscErrorCode (*func)(Tao, Vec, Mat, void *), void *ctx)
551: {
552: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
555: if (Hreg) {
558: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "NULL Hessian detected! User must provide valid Hessian for the regularizer.");
559: if (ctx) gn->reg_hess_ctx = ctx;
560: if (func) gn->regularizerhessian = func;
561: if (Hreg) {
562: PetscObjectReference((PetscObject)Hreg);
563: MatDestroy(&gn->Hreg);
564: gn->Hreg = Hreg;
565: }
566: return 0;
567: }