Actual source code: brgn.c
1: #include <../src/tao/leastsquares/impls/brgn/brgn.h>
3: const char *const TaoBRGNRegularizationTypes[] = {"user", "l2prox", "l2pure", "l1dict", "lm", "TaoBRGNRegularizationType", "TAOBRGN_REGULARIZATION_", NULL};
5: static PetscErrorCode GNHessianProd(Mat H, Vec in, Vec out)
6: {
7: TAO_BRGN *gn;
9: PetscFunctionBegin;
10: PetscCall(MatShellGetContext(H, &gn));
11: PetscCall(MatMult(gn->subsolver->ls_jac, in, gn->r_work));
12: PetscCall(MatMultTranspose(gn->subsolver->ls_jac, gn->r_work, out));
13: switch (gn->reg_type) {
14: case TAOBRGN_REGULARIZATION_USER:
15: PetscCall(MatMult(gn->Hreg, in, gn->x_work));
16: PetscCall(VecAXPY(out, gn->lambda, gn->x_work));
17: break;
18: case TAOBRGN_REGULARIZATION_L2PURE:
19: PetscCall(VecAXPY(out, gn->lambda, in));
20: break;
21: case TAOBRGN_REGULARIZATION_L2PROX:
22: PetscCall(VecAXPY(out, gn->lambda, in));
23: break;
24: case TAOBRGN_REGULARIZATION_L1DICT:
25: /* out = out + lambda*D'*(diag.*(D*in)) */
26: if (gn->D) {
27: PetscCall(MatMult(gn->D, in, gn->y)); /* y = D*in */
28: } else {
29: PetscCall(VecCopy(in, gn->y));
30: }
31: PetscCall(VecPointwiseMult(gn->y_work, gn->diag, gn->y)); /* y_work = diag.*(D*in), where diag = epsilon^2 ./ sqrt(x.^2+epsilon^2).^3 */
32: if (gn->D) {
33: PetscCall(MatMultTranspose(gn->D, gn->y_work, gn->x_work)); /* x_work = D'*(diag.*(D*in)) */
34: } else {
35: PetscCall(VecCopy(gn->y_work, gn->x_work));
36: }
37: PetscCall(VecAXPY(out, gn->lambda, gn->x_work));
38: break;
39: case TAOBRGN_REGULARIZATION_LM:
40: PetscCall(VecPointwiseMult(gn->x_work, gn->damping, in));
41: PetscCall(VecAXPY(out, 1, gn->x_work));
42: break;
43: }
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
46: static PetscErrorCode ComputeDamping(TAO_BRGN *gn)
47: {
48: const PetscScalar *diag_ary;
49: PetscScalar *damping_ary;
50: PetscInt i, n;
52: PetscFunctionBegin;
53: /* update damping */
54: PetscCall(VecGetArray(gn->damping, &damping_ary));
55: PetscCall(VecGetArrayRead(gn->diag, &diag_ary));
56: PetscCall(VecGetLocalSize(gn->damping, &n));
57: for (i = 0; i < n; i++) damping_ary[i] = PetscClipInterval(diag_ary[i], PETSC_SQRT_MACHINE_EPSILON, PetscSqrtReal(PETSC_MAX_REAL));
58: PetscCall(VecScale(gn->damping, gn->lambda));
59: PetscCall(VecRestoreArray(gn->damping, &damping_ary));
60: PetscCall(VecRestoreArrayRead(gn->diag, &diag_ary));
61: PetscFunctionReturn(PETSC_SUCCESS);
62: }
63: /*@
64: TaoBRGNGetDampingVector - Get the damping vector $\mathrm{diag}(J^T J)$ from a `TAOBRGN` with `TAOBRGN_REGULARIZATION_LM` regularization
66: Collective
68: Input Parameter:
69: . tao - a `Tao` of type `TAOBRGN` with `TAOBRGN_REGULARIZATION_LM` regularization
71: Output Parameter:
72: . d - the damping vector
74: Level: developer
76: .seealso: [](ch_tao), `Tao`, `TAOBRGN`, `TaoBRGNRegularzationTypes`
77: @*/
78: PetscErrorCode TaoBRGNGetDampingVector(Tao tao, Vec *d)
79: {
80: PetscFunctionBegin;
82: PetscAssertPointer(d, 2);
83: PetscUseMethod((PetscObject)tao, "TaoBRGNGetDampingVector_C", (Tao, Vec *), (tao, d));
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: static PetscErrorCode TaoBRGNGetDampingVector_BRGN(Tao tao, Vec *d)
88: {
89: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
91: PetscFunctionBegin;
92: PetscCheck(gn->reg_type == TAOBRGN_REGULARIZATION_LM, PetscObjectComm((PetscObject)tao), PETSC_ERR_SUP, "Damping vector is only available if regularization type is lm.");
93: *d = gn->damping;
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode GNObjectiveGradientEval(Tao tao, Vec X, PetscReal *fcn, Vec G, void *ptr)
98: {
99: TAO_BRGN *gn = (TAO_BRGN *)ptr;
100: PetscInt K; /* dimension of D*X */
101: PetscScalar yESum;
102: PetscReal f_reg;
104: PetscFunctionBegin;
105: /* compute objective *fcn*/
106: /* compute first term 0.5*||ls_res||_2^2 */
107: PetscCall(TaoComputeResidual(tao, X, tao->ls_res));
108: PetscCall(VecDot(tao->ls_res, tao->ls_res, fcn));
109: *fcn *= 0.5;
110: /* compute gradient G */
111: PetscCall(TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre));
112: PetscCall(MatMultTranspose(tao->ls_jac, tao->ls_res, G));
113: /* add the regularization contribution */
114: switch (gn->reg_type) {
115: case TAOBRGN_REGULARIZATION_USER:
116: PetscCall((*gn->regularizerobjandgrad)(tao, X, &f_reg, gn->x_work, gn->reg_obj_ctx));
117: *fcn += gn->lambda * f_reg;
118: PetscCall(VecAXPY(G, gn->lambda, gn->x_work));
119: break;
120: case TAOBRGN_REGULARIZATION_L2PURE:
121: /* compute f = f + lambda*0.5*xk'*xk */
122: PetscCall(VecDot(X, X, &f_reg));
123: *fcn += gn->lambda * 0.5 * f_reg;
124: /* compute G = G + lambda*xk */
125: PetscCall(VecAXPY(G, gn->lambda, X));
126: break;
127: case TAOBRGN_REGULARIZATION_L2PROX:
128: /* compute f = f + lambda*0.5*(xk - xkm1)'*(xk - xkm1) */
129: PetscCall(VecAXPBYPCZ(gn->x_work, 1.0, -1.0, 0.0, X, gn->x_old));
130: PetscCall(VecDot(gn->x_work, gn->x_work, &f_reg));
131: *fcn += gn->lambda * 0.5 * f_reg;
132: /* compute G = G + lambda*(xk - xkm1) */
133: PetscCall(VecAXPBYPCZ(G, gn->lambda, -gn->lambda, 1.0, X, gn->x_old));
134: break;
135: case TAOBRGN_REGULARIZATION_L1DICT:
136: /* compute f = f + lambda*sum(sqrt(y.^2+epsilon^2) - epsilon), where y = D*x*/
137: if (gn->D) {
138: PetscCall(MatMult(gn->D, X, gn->y)); /* y = D*x */
139: } else {
140: PetscCall(VecCopy(X, gn->y));
141: }
142: PetscCall(VecPointwiseMult(gn->y_work, gn->y, gn->y));
143: PetscCall(VecShift(gn->y_work, gn->epsilon * gn->epsilon));
144: PetscCall(VecSqrtAbs(gn->y_work)); /* gn->y_work = sqrt(y.^2+epsilon^2) */
145: PetscCall(VecSum(gn->y_work, &yESum));
146: PetscCall(VecGetSize(gn->y, &K));
147: *fcn += gn->lambda * (yESum - K * gn->epsilon);
148: /* compute G = G + lambda*D'*(y./sqrt(y.^2+epsilon^2)),where y = D*x */
149: PetscCall(VecPointwiseDivide(gn->y_work, gn->y, gn->y_work)); /* reuse y_work = y./sqrt(y.^2+epsilon^2) */
150: if (gn->D) {
151: PetscCall(MatMultTranspose(gn->D, gn->y_work, gn->x_work));
152: } else {
153: PetscCall(VecCopy(gn->y_work, gn->x_work));
154: }
155: PetscCall(VecAXPY(G, gn->lambda, gn->x_work));
156: break;
157: case TAOBRGN_REGULARIZATION_LM:
158: break;
159: default:
160: break;
161: }
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode GNComputeHessian(Tao tao, Vec X, Mat H, Mat Hpre, void *ptr)
166: {
167: TAO_BRGN *gn = (TAO_BRGN *)ptr;
168: PetscInt i, n, cstart, cend;
169: PetscScalar *cnorms, *diag_ary;
171: PetscFunctionBegin;
172: PetscCall(TaoComputeResidualJacobian(tao, X, tao->ls_jac, tao->ls_jac_pre));
173: if (gn->mat_explicit) PetscCall(MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_REUSE_MATRIX, PETSC_DETERMINE, &gn->H));
175: switch (gn->reg_type) {
176: case TAOBRGN_REGULARIZATION_USER:
177: PetscCall((*gn->regularizerhessian)(tao, X, gn->Hreg, gn->reg_hess_ctx));
178: if (gn->mat_explicit) PetscCall(MatAXPY(gn->H, 1.0, gn->Hreg, DIFFERENT_NONZERO_PATTERN));
179: break;
180: case TAOBRGN_REGULARIZATION_L2PURE:
181: if (gn->mat_explicit) PetscCall(MatShift(gn->H, gn->lambda));
182: break;
183: case TAOBRGN_REGULARIZATION_L2PROX:
184: if (gn->mat_explicit) PetscCall(MatShift(gn->H, gn->lambda));
185: break;
186: case TAOBRGN_REGULARIZATION_L1DICT:
187: /* 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 */
188: if (gn->D) {
189: PetscCall(MatMult(gn->D, X, gn->y)); /* y = D*x */
190: } else {
191: PetscCall(VecCopy(X, gn->y));
192: }
193: PetscCall(VecPointwiseMult(gn->y_work, gn->y, gn->y));
194: PetscCall(VecShift(gn->y_work, gn->epsilon * gn->epsilon));
195: PetscCall(VecCopy(gn->y_work, gn->diag)); /* gn->diag = y.^2+epsilon^2 */
196: PetscCall(VecSqrtAbs(gn->y_work)); /* gn->y_work = sqrt(y.^2+epsilon^2) */
197: PetscCall(VecPointwiseMult(gn->diag, gn->y_work, gn->diag)); /* gn->diag = sqrt(y.^2+epsilon^2).^3 */
198: PetscCall(VecReciprocal(gn->diag));
199: PetscCall(VecScale(gn->diag, gn->epsilon * gn->epsilon));
200: if (gn->mat_explicit) PetscCall(MatDiagonalSet(gn->H, gn->diag, ADD_VALUES));
201: break;
202: case TAOBRGN_REGULARIZATION_LM:
203: /* compute diagonal of J^T J */
204: PetscCall(MatGetSize(gn->parent->ls_jac, NULL, &n));
205: PetscCall(PetscMalloc1(n, &cnorms));
206: PetscCall(MatGetColumnNorms(gn->parent->ls_jac, NORM_2, cnorms));
207: PetscCall(MatGetOwnershipRangeColumn(gn->parent->ls_jac, &cstart, &cend));
208: PetscCall(VecGetArray(gn->diag, &diag_ary));
209: for (i = 0; i < cend - cstart; i++) diag_ary[i] = cnorms[cstart + i] * cnorms[cstart + i];
210: PetscCall(VecRestoreArray(gn->diag, &diag_ary));
211: PetscCall(PetscFree(cnorms));
212: PetscCall(ComputeDamping(gn));
213: if (gn->mat_explicit) PetscCall(MatDiagonalSet(gn->H, gn->damping, ADD_VALUES));
214: break;
215: default:
216: break;
217: }
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: static PetscErrorCode GNHookFunction(Tao tao, PetscInt iter, void *ctx)
222: {
223: TAO_BRGN *gn = (TAO_BRGN *)ctx;
225: PetscFunctionBegin;
226: /* Update basic tao information from the subsolver */
227: gn->parent->nfuncs = tao->nfuncs;
228: gn->parent->ngrads = tao->ngrads;
229: gn->parent->nfuncgrads = tao->nfuncgrads;
230: gn->parent->nhess = tao->nhess;
231: gn->parent->niter = tao->niter;
232: gn->parent->ksp_its = tao->ksp_its;
233: gn->parent->ksp_tot_its = tao->ksp_tot_its;
234: gn->parent->fc = tao->fc;
235: PetscCall(TaoGetConvergedReason(tao, &gn->parent->reason));
236: /* Update the solution vectors */
237: if (iter == 0) {
238: PetscCall(VecSet(gn->x_old, 0.0));
239: } else {
240: PetscCall(VecCopy(tao->solution, gn->x_old));
241: PetscCall(VecCopy(tao->solution, gn->parent->solution));
242: }
243: /* Update the gradient */
244: PetscCall(VecCopy(tao->gradient, gn->parent->gradient));
246: /* Update damping parameter for LM */
247: if (gn->reg_type == TAOBRGN_REGULARIZATION_LM) {
248: if (iter > 0) {
249: if (gn->fc_old > tao->fc) {
250: gn->lambda = gn->lambda * gn->downhill_lambda_change;
251: } else {
252: /* uphill step */
253: gn->lambda = gn->lambda * gn->uphill_lambda_change;
254: }
255: }
256: gn->fc_old = tao->fc;
257: }
259: /* Call general purpose update function */
260: if (gn->parent->ops->update) PetscCall((*gn->parent->ops->update)(gn->parent, gn->parent->niter, gn->parent->user_update));
261: PetscFunctionReturn(PETSC_SUCCESS);
262: }
264: static PetscErrorCode TaoBRGNGetRegularizationType_BRGN(Tao tao, TaoBRGNRegularizationType *type)
265: {
266: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
268: PetscFunctionBegin;
269: *type = gn->reg_type;
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: TaoBRGNGetRegularizationType - Get the `TaoBRGNRegularizationType` of a `TAOBRGN`
276: Not collective
278: Input Parameter:
279: . tao - a `Tao` of type `TAOBRGN`
281: Output Parameter:
282: . type - the `TaoBRGNRegularizationType`
284: Level: advanced
286: .seealso: [](ch_tao), `Tao`, `TAOBRGN`, `TaoBRGNRegularizationType`, `TaoBRGNSetRegularizationType()`
287: @*/
288: PetscErrorCode TaoBRGNGetRegularizationType(Tao tao, TaoBRGNRegularizationType *type)
289: {
290: PetscFunctionBegin;
292: PetscAssertPointer(type, 2);
293: PetscUseMethod((PetscObject)tao, "TaoBRGNGetRegularizationType_C", (Tao, TaoBRGNRegularizationType *), (tao, type));
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: static PetscErrorCode TaoBRGNSetRegularizationType_BRGN(Tao tao, TaoBRGNRegularizationType type)
298: {
299: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
301: PetscFunctionBegin;
302: gn->reg_type = type;
303: PetscFunctionReturn(PETSC_SUCCESS);
304: }
306: /*@
307: TaoBRGNSetRegularizationType - Set the `TaoBRGNRegularizationType` of a `TAOBRGN`
309: Logically collective
311: Input Parameters:
312: + tao - a `Tao` of type `TAOBRGN`
313: - type - the `TaoBRGNRegularizationType`
315: Level: advanced
317: .seealso: [](ch_tao), `Tao`, `TAOBRGN`, `TaoBRGNRegularizationType`, `TaoBRGNGetRegularizationType`
318: @*/
319: PetscErrorCode TaoBRGNSetRegularizationType(Tao tao, TaoBRGNRegularizationType type)
320: {
321: PetscFunctionBegin;
324: PetscTryMethod((PetscObject)tao, "TaoBRGNSetRegularizationType_C", (Tao, TaoBRGNRegularizationType), (tao, type));
325: PetscFunctionReturn(PETSC_SUCCESS);
326: }
328: static PetscErrorCode TaoSolve_BRGN(Tao tao)
329: {
330: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
332: PetscFunctionBegin;
333: PetscCall(TaoSolve(gn->subsolver));
334: /* Update basic tao information from the subsolver */
335: tao->nfuncs = gn->subsolver->nfuncs;
336: tao->ngrads = gn->subsolver->ngrads;
337: tao->nfuncgrads = gn->subsolver->nfuncgrads;
338: tao->nhess = gn->subsolver->nhess;
339: tao->niter = gn->subsolver->niter;
340: tao->ksp_its = gn->subsolver->ksp_its;
341: tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
342: PetscCall(TaoGetConvergedReason(gn->subsolver, &tao->reason));
343: /* Update vectors */
344: PetscCall(VecCopy(gn->subsolver->solution, tao->solution));
345: PetscCall(VecCopy(gn->subsolver->gradient, tao->gradient));
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: static PetscErrorCode TaoSetFromOptions_BRGN(Tao tao, PetscOptionItems PetscOptionsObject)
350: {
351: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
352: TaoLineSearch ls;
354: PetscFunctionBegin;
355: 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.");
356: 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));
357: PetscCall(PetscOptionsReal("-tao_brgn_regularizer_weight", "regularizer weight (default 1e-4)", "", gn->lambda, &gn->lambda, NULL));
358: 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));
359: 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));
360: 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));
361: PetscCall(PetscOptionsEnum("-tao_brgn_regularization_type", "regularization type", "", TaoBRGNRegularizationTypes, (PetscEnum)gn->reg_type, (PetscEnum *)&gn->reg_type, NULL));
362: PetscOptionsHeadEnd();
363: /* set unit line search direction as the default when using the lm regularizer */
364: if (gn->reg_type == TAOBRGN_REGULARIZATION_LM) {
365: PetscCall(TaoGetLineSearch(gn->subsolver, &ls));
366: PetscCall(TaoLineSearchSetType(ls, TAOLINESEARCHUNIT));
367: }
368: PetscCall(TaoSetFromOptions(gn->subsolver));
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: static PetscErrorCode TaoView_BRGN(Tao tao, PetscViewer viewer)
373: {
374: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
375: PetscBool isascii;
377: PetscFunctionBegin;
378: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
379: if (isascii) {
380: PetscCall(PetscViewerASCIIPushTab(viewer));
381: PetscCall(PetscViewerASCIIPrintf(viewer, "Regularizer weight: %g\n", (double)gn->lambda));
382: PetscCall(PetscViewerASCIIPrintf(viewer, "BRGN Regularization Type: %s\n", TaoBRGNRegularizationTypes[gn->reg_type]));
383: switch (gn->reg_type) {
384: case TAOBRGN_REGULARIZATION_L1DICT:
385: PetscCall(PetscViewerASCIIPrintf(viewer, "L1 smooth epsilon: %g\n", (double)gn->epsilon));
386: break;
387: case TAOBRGN_REGULARIZATION_LM:
388: PetscCall(PetscViewerASCIIPrintf(viewer, "Downhill trust region decrease factor:: %g\n", (double)gn->downhill_lambda_change));
389: PetscCall(PetscViewerASCIIPrintf(viewer, "Uphill trust region increase factor:: %g\n", (double)gn->uphill_lambda_change));
390: break;
391: case TAOBRGN_REGULARIZATION_L2PROX:
392: case TAOBRGN_REGULARIZATION_L2PURE:
393: case TAOBRGN_REGULARIZATION_USER:
394: default:
395: break;
396: }
397: PetscCall(PetscViewerASCIIPopTab(viewer));
398: }
399: PetscCall(PetscViewerASCIIPushTab(viewer));
400: PetscCall(TaoView(gn->subsolver, viewer));
401: PetscCall(PetscViewerASCIIPopTab(viewer));
402: PetscFunctionReturn(PETSC_SUCCESS);
403: }
405: static PetscErrorCode TaoSetUp_BRGN(Tao tao)
406: {
407: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
408: PetscBool is_bnls, is_bntr, is_bntl;
409: PetscInt i, n, N, K; /* dict has size K*N*/
411: PetscFunctionBegin;
412: PetscCheck(tao->ls_res, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualRoutine() must be called before setup!");
413: PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNLS, &is_bnls));
414: PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTR, &is_bntr));
415: PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTL, &is_bntl));
416: PetscCheck((!is_bnls && !is_bntr && !is_bntl) || tao->ls_jac, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualJacobianRoutine() must be called before setup!");
417: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
418: if (!gn->x_work) PetscCall(VecDuplicate(tao->solution, &gn->x_work));
419: if (!gn->r_work) PetscCall(VecDuplicate(tao->ls_res, &gn->r_work));
420: if (!gn->x_old) {
421: PetscCall(VecDuplicate(tao->solution, &gn->x_old));
422: PetscCall(VecSet(gn->x_old, 0.0));
423: }
425: if (TAOBRGN_REGULARIZATION_L1DICT == gn->reg_type) {
426: if (!gn->y) {
427: if (gn->D) {
428: 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 */
429: PetscCall(MatCreateVecs(gn->D, NULL, &gn->y));
430: } else {
431: PetscCall(VecDuplicate(tao->solution, &gn->y)); /* If user does not setup dict matrix, use identity matrix, K=N */
432: }
433: PetscCall(VecSet(gn->y, 0.0));
434: }
435: if (!gn->y_work) PetscCall(VecDuplicate(gn->y, &gn->y_work));
436: if (!gn->diag) {
437: PetscCall(VecDuplicate(gn->y, &gn->diag));
438: PetscCall(VecSet(gn->diag, 0.0));
439: }
440: }
441: if (TAOBRGN_REGULARIZATION_LM == gn->reg_type) {
442: if (!gn->diag) PetscCall(MatCreateVecs(tao->ls_jac, &gn->diag, NULL));
443: if (!gn->damping) PetscCall(MatCreateVecs(tao->ls_jac, &gn->damping, NULL));
444: }
446: if (!tao->setupcalled) {
447: /* Hessian setup */
448: if (gn->mat_explicit) {
449: PetscCall(TaoComputeResidualJacobian(tao, tao->solution, tao->ls_jac, tao->ls_jac_pre));
450: PetscCall(MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &gn->H));
451: } else {
452: PetscCall(VecGetLocalSize(tao->solution, &n));
453: PetscCall(VecGetSize(tao->solution, &N));
454: PetscCall(MatCreate(PetscObjectComm((PetscObject)tao), &gn->H));
455: PetscCall(MatSetSizes(gn->H, n, n, N, N));
456: PetscCall(MatSetType(gn->H, MATSHELL));
457: PetscCall(MatSetOption(gn->H, MAT_SYMMETRIC, PETSC_TRUE));
458: PetscCall(MatShellSetOperation(gn->H, MATOP_MULT, (void (*)(void))GNHessianProd));
459: PetscCall(MatShellSetContext(gn->H, gn));
460: }
461: PetscCall(MatSetUp(gn->H));
462: /* Subsolver setup,include initial vector and dictionary D */
463: PetscCall(TaoSetUpdate(gn->subsolver, GNHookFunction, gn));
464: PetscCall(TaoSetSolution(gn->subsolver, tao->solution));
465: if (tao->bounded) PetscCall(TaoSetVariableBounds(gn->subsolver, tao->XL, tao->XU));
466: PetscCall(TaoSetResidualRoutine(gn->subsolver, tao->ls_res, tao->ops->computeresidual, tao->user_lsresP));
467: PetscCall(TaoSetJacobianResidualRoutine(gn->subsolver, tao->ls_jac, tao->ls_jac, tao->ops->computeresidualjacobian, tao->user_lsjacP));
468: PetscCall(TaoSetObjectiveAndGradient(gn->subsolver, NULL, GNObjectiveGradientEval, gn));
469: PetscCall(TaoSetHessian(gn->subsolver, gn->H, gn->H, GNComputeHessian, gn));
470: /* Propagate some options down */
471: PetscCall(TaoSetTolerances(gn->subsolver, tao->gatol, tao->grtol, tao->gttol));
472: PetscCall(TaoSetMaximumIterations(gn->subsolver, tao->max_it));
473: PetscCall(TaoSetMaximumFunctionEvaluations(gn->subsolver, tao->max_funcs));
474: for (i = 0; i < tao->numbermonitors; ++i) {
475: PetscCall(TaoMonitorSet(gn->subsolver, tao->monitor[i], tao->monitorcontext[i], tao->monitordestroy[i]));
476: PetscCall(PetscObjectReference((PetscObject)tao->monitorcontext[i]));
477: }
478: PetscCall(TaoSetUp(gn->subsolver));
479: }
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: static PetscErrorCode TaoDestroy_BRGN(Tao tao)
484: {
485: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
487: PetscFunctionBegin;
488: if (tao->setupcalled) {
489: PetscCall(VecDestroy(&tao->gradient));
490: PetscCall(VecDestroy(&gn->x_work));
491: PetscCall(VecDestroy(&gn->r_work));
492: PetscCall(VecDestroy(&gn->x_old));
493: PetscCall(VecDestroy(&gn->diag));
494: PetscCall(VecDestroy(&gn->y));
495: PetscCall(VecDestroy(&gn->y_work));
496: }
497: PetscCall(VecDestroy(&gn->damping));
498: PetscCall(VecDestroy(&gn->diag));
499: PetscCall(MatDestroy(&gn->H));
500: PetscCall(MatDestroy(&gn->D));
501: PetscCall(MatDestroy(&gn->Hreg));
502: PetscCall(TaoDestroy(&gn->subsolver));
503: gn->parent = NULL;
504: PetscCall(PetscFree(tao->data));
505: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetRegularizationType_C", NULL));
506: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizationType_C", NULL));
507: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetDampingVector_C", NULL));
508: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetDictionaryMatrix_C", NULL));
509: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetSubsolver_C", NULL));
510: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerWeight_C", NULL));
511: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetL1SmoothEpsilon_C", NULL));
512: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerObjectiveAndGradientRoutine_C", NULL));
513: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerHessianRoutine_C", NULL));
514: PetscFunctionReturn(PETSC_SUCCESS);
515: }
517: /*@
518: TaoBRGNGetSubsolver - Get the pointer to the subsolver inside a `TAOBRGN`
520: Collective
522: Input Parameters:
523: + tao - the Tao solver context
524: - subsolver - the `Tao` sub-solver context
526: Level: advanced
528: .seealso: `Tao`, `Mat`, `TAOBRGN`
529: @*/
530: PetscErrorCode TaoBRGNGetSubsolver(Tao tao, Tao *subsolver)
531: {
532: PetscFunctionBegin;
534: PetscUseMethod((PetscObject)tao, "TaoBRGNGetSubsolver_C", (Tao, Tao *), (tao, subsolver));
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: static PetscErrorCode TaoBRGNGetSubsolver_BRGN(Tao tao, Tao *subsolver)
539: {
540: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
542: PetscFunctionBegin;
543: *subsolver = gn->subsolver;
544: PetscFunctionReturn(PETSC_SUCCESS);
545: }
547: /*@
548: TaoBRGNSetRegularizerWeight - Set the regularizer weight for the Gauss-Newton least-squares algorithm
550: Collective
552: Input Parameters:
553: + tao - the `Tao` solver context
554: - lambda - L1-norm regularizer weight
556: Level: beginner
558: .seealso: `Tao`, `Mat`, `TAOBRGN`
559: @*/
560: PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao, PetscReal lambda)
561: {
562: PetscFunctionBegin;
565: PetscTryMethod((PetscObject)tao, "TaoBRGNSetRegularizerWeight_C", (Tao, PetscReal), (tao, lambda));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: static PetscErrorCode TaoBRGNSetRegularizerWeight_BRGN(Tao tao, PetscReal lambda)
570: {
571: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
573: PetscFunctionBegin;
574: gn->lambda = lambda;
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: /*@
579: TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm
581: Collective
583: Input Parameters:
584: + tao - the `Tao` solver context
585: - epsilon - L1-norm smooth approximation parameter
587: Level: advanced
589: .seealso: `Tao`, `Mat`, `TAOBRGN`
590: @*/
591: PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao, PetscReal epsilon)
592: {
593: PetscFunctionBegin;
596: PetscTryMethod((PetscObject)tao, "TaoBRGNSetL1SmoothEpsilon_C", (Tao, PetscReal), (tao, epsilon));
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: static PetscErrorCode TaoBRGNSetL1SmoothEpsilon_BRGN(Tao tao, PetscReal epsilon)
601: {
602: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
604: PetscFunctionBegin;
605: gn->epsilon = epsilon;
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: /*@
610: TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)
612: Input Parameters:
613: + tao - the `Tao` context
614: - dict - the user specified dictionary matrix. We allow to set a `NULL` dictionary, which means identity matrix by default
616: Level: advanced
618: .seealso: `Tao`, `Mat`, `TAOBRGN`
619: @*/
620: PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao, Mat dict)
621: {
622: PetscFunctionBegin;
624: PetscTryMethod((PetscObject)tao, "TaoBRGNSetDictionaryMatrix_C", (Tao, Mat), (tao, dict));
625: PetscFunctionReturn(PETSC_SUCCESS);
626: }
628: static PetscErrorCode TaoBRGNSetDictionaryMatrix_BRGN(Tao tao, Mat dict)
629: {
630: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
632: PetscFunctionBegin;
633: if (dict) {
635: PetscCheckSameComm(tao, 1, dict, 2);
636: PetscCall(PetscObjectReference((PetscObject)dict));
637: }
638: PetscCall(MatDestroy(&gn->D));
639: gn->D = dict;
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: /*@C
644: TaoBRGNSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back
645: function into the algorithm.
647: Input Parameters:
648: + tao - the Tao context
649: . func - function pointer for the regularizer value and gradient evaluation
650: - ctx - user context for the regularizer
652: Calling sequence:
653: + tao - the `Tao` context
654: . u - the location at which to compute the objective and gradient
655: . val - location to store objective function value
656: . g - location to store gradient
657: - ctx - user context for the regularizer Hessian
659: Level: advanced
661: .seealso: `Tao`, `Mat`, `TAOBRGN`
662: @*/
663: PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao, PetscErrorCode (*func)(Tao tao, Vec u, PetscReal *val, Vec g, void *ctx), void *ctx)
664: {
665: PetscFunctionBegin;
667: PetscTryMethod((PetscObject)tao, "TaoBRGNSetRegularizerObjectiveAndGradientRoutine_C", (Tao, PetscErrorCode (*)(Tao, Vec, PetscReal *, Vec, void *), void *), (tao, func, ctx));
668: PetscFunctionReturn(PETSC_SUCCESS);
669: }
671: static PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine_BRGN(Tao tao, PetscErrorCode (*func)(Tao tao, Vec u, PetscReal *val, Vec g, void *ctx), void *ctx)
672: {
673: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
675: PetscFunctionBegin;
676: if (ctx) gn->reg_obj_ctx = ctx;
677: if (func) gn->regularizerobjandgrad = func;
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: /*@C
682: TaoBRGNSetRegularizerHessianRoutine - Sets the user-defined regularizer call-back
683: function into the algorithm.
685: Input Parameters:
686: + tao - the `Tao` context
687: . Hreg - user-created matrix for the Hessian of the regularization term
688: . func - function pointer for the regularizer Hessian evaluation
689: - ctx - user context for the regularizer Hessian
691: Calling sequence:
692: + tao - the `Tao` context
693: . u - the location at which to compute the Hessian
694: . Hreg - user-created matrix for the Hessian of the regularization term
695: - ctx - user context for the regularizer Hessian
697: Level: advanced
699: .seealso: `Tao`, `Mat`, `TAOBRGN`
700: @*/
701: PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao, Mat Hreg, PetscErrorCode (*func)(Tao tao, Vec u, Mat Hreg, void *ctx), void *ctx)
702: {
703: PetscFunctionBegin;
705: PetscTryMethod((PetscObject)tao, "TaoBRGNSetRegularizerHessianRoutine_C", (Tao, Mat, PetscErrorCode (*)(Tao, Vec, Mat, void *), void *), (tao, Hreg, func, ctx));
706: PetscFunctionReturn(PETSC_SUCCESS);
707: }
709: static PetscErrorCode TaoBRGNSetRegularizerHessianRoutine_BRGN(Tao tao, Mat Hreg, PetscErrorCode (*func)(Tao tao, Vec u, Mat Hreg, void *ctx), void *ctx)
710: {
711: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
713: PetscFunctionBegin;
714: if (Hreg) {
716: PetscCheckSameComm(tao, 1, Hreg, 2);
717: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "NULL Hessian detected! User must provide valid Hessian for the regularizer.");
718: if (ctx) gn->reg_hess_ctx = ctx;
719: if (func) gn->regularizerhessian = func;
720: if (Hreg) {
721: PetscCall(PetscObjectReference((PetscObject)Hreg));
722: PetscCall(MatDestroy(&gn->Hreg));
723: gn->Hreg = Hreg;
724: }
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: /*MC
729: TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares
730: problems with bound constraints. This algorithm is a thin wrapper around `TAOBNTL`
731: that constructs the Gauss-Newton problem with the user-provided least-squares
732: residual and Jacobian. The algorithm offers an L2-norm ("l2pure"), L2-norm proximal point ("l2prox")
733: regularizer, and L1-norm dictionary regularizer ("l1dict"), where we approximate the
734: L1-norm ||x||_1 by sum_i(sqrt(x_i^2+epsilon^2)-epsilon) with a small positive number epsilon.
735: Also offered is the "lm" regularizer which uses a scaled diagonal of J^T J.
736: With the "lm" regularizer, `TAOBRGN` is a Levenberg-Marquardt optimizer.
737: The user can also provide own regularization function.
739: Options Database Keys:
740: + -tao_brgn_regularization_type - regularization type ("user", "l2prox", "l2pure", "l1dict", "lm") (default "l2prox")
741: . -tao_brgn_regularizer_weight - regularizer weight (default 1e-4)
742: - -tao_brgn_l1_smooth_epsilon - L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)
744: Level: beginner
746: .seealso: `Tao`, `TaoBRGNGetSubsolver()`, `TaoBRGNSetRegularizerWeight()`, `TaoBRGNSetL1SmoothEpsilon()`, `TaoBRGNSetDictionaryMatrix()`,
747: `TaoBRGNSetRegularizerObjectiveAndGradientRoutine()`, `TaoBRGNSetRegularizerHessianRoutine()`
748: M*/
749: PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao)
750: {
751: TAO_BRGN *gn;
753: PetscFunctionBegin;
754: PetscCall(PetscNew(&gn));
756: tao->ops->destroy = TaoDestroy_BRGN;
757: tao->ops->setup = TaoSetUp_BRGN;
758: tao->ops->setfromoptions = TaoSetFromOptions_BRGN;
759: tao->ops->view = TaoView_BRGN;
760: tao->ops->solve = TaoSolve_BRGN;
762: PetscCall(TaoParametersInitialize(tao));
764: tao->data = gn;
765: gn->reg_type = TAOBRGN_REGULARIZATION_L2PROX;
766: gn->lambda = 1e-4;
767: gn->epsilon = 1e-6;
768: gn->downhill_lambda_change = 1. / 5.;
769: gn->uphill_lambda_change = 1.5;
770: gn->parent = tao;
772: PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao), &gn->subsolver));
773: PetscCall(TaoSetType(gn->subsolver, TAOBNLS));
774: PetscCall(TaoSetOptionsPrefix(gn->subsolver, "tao_brgn_subsolver_"));
775: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetRegularizationType_C", TaoBRGNGetRegularizationType_BRGN));
776: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizationType_C", TaoBRGNSetRegularizationType_BRGN));
777: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetDampingVector_C", TaoBRGNGetDampingVector_BRGN));
778: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetDictionaryMatrix_C", TaoBRGNSetDictionaryMatrix_BRGN));
779: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetSubsolver_C", TaoBRGNGetSubsolver_BRGN));
780: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerWeight_C", TaoBRGNSetRegularizerWeight_BRGN));
781: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetL1SmoothEpsilon_C", TaoBRGNSetL1SmoothEpsilon_BRGN));
782: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerObjectiveAndGradientRoutine_C", TaoBRGNSetRegularizerObjectiveAndGradientRoutine_BRGN));
783: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerHessianRoutine_C", TaoBRGNSetRegularizerHessianRoutine_BRGN));
784: PetscFunctionReturn(PETSC_SUCCESS);
785: }