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