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, PetscCtx ctx)
222: {
223: TAO_BRGN *gn = (TAO_BRGN *)ctx;
225: PetscFunctionBegin;
226: /* Update basic tao information from the subsolver */
227: {
228: gn->parent->objective_term.term->nobj = tao->objective_term.term->nobj;
229: gn->parent->objective_term.term->ngrad = tao->objective_term.term->ngrad;
230: gn->parent->objective_term.term->nobjgrad = tao->objective_term.term->nobjgrad;
231: gn->parent->objective_term.term->nhess = tao->objective_term.term->nhess;
232: }
233: gn->parent->nres = tao->nres;
234: gn->parent->niter = tao->niter;
235: gn->parent->ksp_its = tao->ksp_its;
236: gn->parent->ksp_tot_its = tao->ksp_tot_its;
237: gn->parent->fc = tao->fc;
238: PetscCall(TaoGetConvergedReason(tao, &gn->parent->reason));
239: /* Update the solution vectors */
240: if (iter == 0) {
241: PetscCall(VecSet(gn->x_old, 0.0));
242: } else {
243: PetscCall(VecCopy(tao->solution, gn->x_old));
244: PetscCall(VecCopy(tao->solution, gn->parent->solution));
245: }
246: /* Update the gradient */
247: PetscCall(VecCopy(tao->gradient, gn->parent->gradient));
249: /* Update damping parameter for LM */
250: if (gn->reg_type == TAOBRGN_REGULARIZATION_LM) {
251: if (iter > 0) {
252: if (gn->fc_old > tao->fc) {
253: gn->lambda = gn->lambda * gn->downhill_lambda_change;
254: } else {
255: /* uphill step */
256: gn->lambda = gn->lambda * gn->uphill_lambda_change;
257: }
258: }
259: gn->fc_old = tao->fc;
260: }
262: /* Call general purpose update function */
263: PetscTryTypeMethod(gn->parent, update, gn->parent->niter, gn->parent->user_update);
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static PetscErrorCode TaoBRGNGetRegularizationType_BRGN(Tao tao, TaoBRGNRegularizationType *type)
268: {
269: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
271: PetscFunctionBegin;
272: *type = gn->reg_type;
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: /*@
277: TaoBRGNGetRegularizationType - Get the `TaoBRGNRegularizationType` of a `TAOBRGN`
279: Not collective
281: Input Parameter:
282: . tao - a `Tao` of type `TAOBRGN`
284: Output Parameter:
285: . type - the `TaoBRGNRegularizationType`
287: Level: advanced
289: .seealso: [](ch_tao), `Tao`, `TAOBRGN`, `TaoBRGNRegularizationType`, `TaoBRGNSetRegularizationType()`
290: @*/
291: PetscErrorCode TaoBRGNGetRegularizationType(Tao tao, TaoBRGNRegularizationType *type)
292: {
293: PetscFunctionBegin;
295: PetscAssertPointer(type, 2);
296: PetscUseMethod((PetscObject)tao, "TaoBRGNGetRegularizationType_C", (Tao, TaoBRGNRegularizationType *), (tao, type));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode TaoBRGNSetRegularizationType_BRGN(Tao tao, TaoBRGNRegularizationType type)
301: {
302: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
304: PetscFunctionBegin;
305: gn->reg_type = type;
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: /*@
310: TaoBRGNSetRegularizationType - Set the `TaoBRGNRegularizationType` of a `TAOBRGN`
312: Logically collective
314: Input Parameters:
315: + tao - a `Tao` of type `TAOBRGN`
316: - type - the `TaoBRGNRegularizationType`
318: Level: advanced
320: .seealso: [](ch_tao), `Tao`, `TAOBRGN`, `TaoBRGNRegularizationType`, `TaoBRGNGetRegularizationType`
321: @*/
322: PetscErrorCode TaoBRGNSetRegularizationType(Tao tao, TaoBRGNRegularizationType type)
323: {
324: PetscFunctionBegin;
327: PetscTryMethod((PetscObject)tao, "TaoBRGNSetRegularizationType_C", (Tao, TaoBRGNRegularizationType), (tao, type));
328: PetscFunctionReturn(PETSC_SUCCESS);
329: }
331: static PetscErrorCode TaoSolve_BRGN(Tao tao)
332: {
333: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
335: PetscFunctionBegin;
336: PetscCall(TaoSolve(gn->subsolver));
337: /* Update basic tao information from the subsolver */
338: /* NOTE: Due to subsolver nature of BRGN, TaoTerm cannot properly track counts */
339: tao->objective_term.term->nobj = gn->subsolver->objective_term.term->nobj;
340: tao->objective_term.term->ngrad = gn->subsolver->objective_term.term->ngrad;
341: tao->objective_term.term->nobjgrad = gn->subsolver->objective_term.term->nobjgrad;
342: tao->objective_term.term->nhess = gn->subsolver->objective_term.term->nhess;
344: tao->nres = gn->subsolver->nres;
345: tao->niter = gn->subsolver->niter;
346: tao->ksp_its = gn->subsolver->ksp_its;
347: tao->ksp_tot_its = gn->subsolver->ksp_tot_its;
348: PetscCall(TaoGetConvergedReason(gn->subsolver, &tao->reason));
349: /* Update vectors */
350: PetscCall(VecCopy(gn->subsolver->solution, tao->solution));
351: PetscCall(VecCopy(gn->subsolver->gradient, tao->gradient));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: static PetscErrorCode TaoSetFromOptions_BRGN(Tao tao, PetscOptionItems PetscOptionsObject)
356: {
357: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
358: TaoLineSearch ls;
360: PetscFunctionBegin;
361: 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.");
362: 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));
363: PetscCall(PetscOptionsReal("-tao_brgn_regularizer_weight", "regularizer weight (default 1e-4)", "", gn->lambda, &gn->lambda, NULL));
364: 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));
365: 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));
366: 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));
367: PetscCall(PetscOptionsEnum("-tao_brgn_regularization_type", "regularization type", "", TaoBRGNRegularizationTypes, (PetscEnum)gn->reg_type, (PetscEnum *)&gn->reg_type, NULL));
368: PetscOptionsHeadEnd();
369: /* set unit line search direction as the default when using the lm regularizer */
370: if (gn->reg_type == TAOBRGN_REGULARIZATION_LM) {
371: PetscCall(TaoGetLineSearch(gn->subsolver, &ls));
372: PetscCall(TaoLineSearchSetType(ls, TAOLINESEARCHUNIT));
373: }
374: PetscCall(TaoSetFromOptions(gn->subsolver));
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: static PetscErrorCode TaoView_BRGN(Tao tao, PetscViewer viewer)
379: {
380: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
381: PetscBool isascii;
383: PetscFunctionBegin;
384: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
385: if (isascii) {
386: PetscCall(PetscViewerASCIIPushTab(viewer));
387: PetscCall(PetscViewerASCIIPrintf(viewer, "Regularizer weight: %g\n", (double)gn->lambda));
388: PetscCall(PetscViewerASCIIPrintf(viewer, "BRGN Regularization Type: %s\n", TaoBRGNRegularizationTypes[gn->reg_type]));
389: switch (gn->reg_type) {
390: case TAOBRGN_REGULARIZATION_L1DICT:
391: PetscCall(PetscViewerASCIIPrintf(viewer, "L1 smooth epsilon: %g\n", (double)gn->epsilon));
392: break;
393: case TAOBRGN_REGULARIZATION_LM:
394: PetscCall(PetscViewerASCIIPrintf(viewer, "Downhill trust region decrease factor:: %g\n", (double)gn->downhill_lambda_change));
395: PetscCall(PetscViewerASCIIPrintf(viewer, "Uphill trust region increase factor:: %g\n", (double)gn->uphill_lambda_change));
396: break;
397: case TAOBRGN_REGULARIZATION_L2PROX:
398: case TAOBRGN_REGULARIZATION_L2PURE:
399: case TAOBRGN_REGULARIZATION_USER:
400: default:
401: break;
402: }
403: PetscCall(PetscViewerASCIIPopTab(viewer));
404: }
405: PetscCall(PetscViewerASCIIPushTab(viewer));
406: PetscCall(TaoView(gn->subsolver, viewer));
407: PetscCall(PetscViewerASCIIPopTab(viewer));
408: PetscFunctionReturn(PETSC_SUCCESS);
409: }
411: static PetscErrorCode TaoSetUp_BRGN(Tao tao)
412: {
413: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
414: PetscBool is_bnls, is_bntr, is_bntl;
415: PetscInt n, N, K; /* dict has size K*N*/
417: PetscFunctionBegin;
418: PetscCheck(tao->ls_res, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualRoutine() must be called before setup!");
419: PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNLS, &is_bnls));
420: PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTR, &is_bntr));
421: PetscCall(PetscObjectTypeCompare((PetscObject)gn->subsolver, TAOBNTL, &is_bntl));
422: PetscCheck((!is_bnls && !is_bntr && !is_bntl) || tao->ls_jac, PetscObjectComm((PetscObject)tao), PETSC_ERR_ORDER, "TaoSetResidualJacobianRoutine() must be called before setup!");
423: if (!tao->gradient) PetscCall(VecDuplicate(tao->solution, &tao->gradient));
424: if (!gn->x_work) PetscCall(VecDuplicate(tao->solution, &gn->x_work));
425: if (!gn->r_work) PetscCall(VecDuplicate(tao->ls_res, &gn->r_work));
426: if (!gn->x_old) PetscCall(VecDuplicate(tao->solution, &gn->x_old));
428: if (TAOBRGN_REGULARIZATION_L1DICT == gn->reg_type) {
429: if (!gn->y) {
430: if (gn->D) {
431: 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 */
432: PetscCall(MatCreateVecs(gn->D, NULL, &gn->y));
433: } else {
434: PetscCall(VecDuplicate(tao->solution, &gn->y)); /* If user does not setup dict matrix, use identity matrix, K=N */
435: }
436: }
437: if (!gn->y_work) PetscCall(VecDuplicate(gn->y, &gn->y_work));
438: if (!gn->diag) PetscCall(VecDuplicate(gn->y, &gn->diag));
439: }
440: if (TAOBRGN_REGULARIZATION_LM == gn->reg_type) {
441: if (!gn->diag) PetscCall(MatCreateVecs(tao->ls_jac, &gn->diag, NULL));
442: if (!gn->damping) PetscCall(MatCreateVecs(tao->ls_jac, &gn->damping, NULL));
443: }
445: if (!tao->setupcalled) {
446: /* Hessian setup */
447: if (gn->mat_explicit) {
448: PetscCall(TaoComputeResidualJacobian(tao, tao->solution, tao->ls_jac, tao->ls_jac_pre));
449: PetscCall(MatTransposeMatMult(tao->ls_jac, tao->ls_jac, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &gn->H));
450: } else {
451: PetscCall(VecGetLocalSize(tao->solution, &n));
452: PetscCall(VecGetSize(tao->solution, &N));
453: PetscCall(MatCreate(PetscObjectComm((PetscObject)tao), &gn->H));
454: PetscCall(MatSetSizes(gn->H, n, n, N, N));
455: PetscCall(MatSetType(gn->H, MATSHELL));
456: PetscCall(MatSetOption(gn->H, MAT_SYMMETRIC, PETSC_TRUE));
457: PetscCall(MatShellSetOperation(gn->H, MATOP_MULT, (PetscErrorCodeFn *)GNHessianProd));
458: PetscCall(MatShellSetContext(gn->H, gn));
459: }
460: PetscCall(MatSetUp(gn->H));
461: /* Subsolver setup,include initial vector and dictionary D */
462: PetscCall(TaoSetUpdate(gn->subsolver, GNHookFunction, gn));
463: PetscCall(TaoSetSolution(gn->subsolver, tao->solution));
464: if (tao->bounded) PetscCall(TaoSetVariableBounds(gn->subsolver, tao->XL, tao->XU));
465: PetscCall(TaoSetResidualRoutine(gn->subsolver, tao->ls_res, tao->ops->computeresidual, tao->user_lsresP));
466: PetscCall(TaoSetJacobianResidualRoutine(gn->subsolver, tao->ls_jac, tao->ls_jac, tao->ops->computeresidualjacobian, tao->user_lsjacP));
467: PetscCall(TaoSetObjectiveAndGradient(gn->subsolver, NULL, GNObjectiveGradientEval, gn));
468: PetscCall(TaoSetHessian(gn->subsolver, gn->H, gn->H, GNComputeHessian, gn));
469: /* Propagate some options down */
470: PetscCall(TaoSetTolerances(gn->subsolver, tao->gatol, tao->grtol, tao->gttol));
471: PetscCall(TaoSetMaximumIterations(gn->subsolver, tao->max_it));
472: PetscCall(TaoSetMaximumFunctionEvaluations(gn->subsolver, tao->max_funcs));
473: PetscCall(TaoSetUp(gn->subsolver));
474: }
475: PetscFunctionReturn(PETSC_SUCCESS);
476: }
478: static PetscErrorCode TaoDestroy_BRGN(Tao tao)
479: {
480: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
482: PetscFunctionBegin;
483: if (tao->setupcalled) {
484: PetscCall(VecDestroy(&tao->gradient));
485: PetscCall(VecDestroy(&gn->x_work));
486: PetscCall(VecDestroy(&gn->r_work));
487: PetscCall(VecDestroy(&gn->x_old));
488: PetscCall(VecDestroy(&gn->diag));
489: PetscCall(VecDestroy(&gn->y));
490: PetscCall(VecDestroy(&gn->y_work));
491: }
492: PetscCall(VecDestroy(&gn->damping));
493: PetscCall(VecDestroy(&gn->diag));
494: PetscCall(MatDestroy(&gn->H));
495: PetscCall(MatDestroy(&gn->D));
496: PetscCall(MatDestroy(&gn->Hreg));
497: PetscCall(TaoDestroy(&gn->subsolver));
498: gn->parent = NULL;
499: PetscCall(PetscFree(tao->data));
500: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetRegularizationType_C", NULL));
501: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizationType_C", NULL));
502: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetDampingVector_C", NULL));
503: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetDictionaryMatrix_C", NULL));
504: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetSubsolver_C", NULL));
505: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerWeight_C", NULL));
506: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetL1SmoothEpsilon_C", NULL));
507: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerObjectiveAndGradientRoutine_C", NULL));
508: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerHessianRoutine_C", NULL));
509: PetscFunctionReturn(PETSC_SUCCESS);
510: }
512: /*@
513: TaoBRGNGetSubsolver - Get the pointer to the subsolver inside a `TAOBRGN`
515: Collective
517: Input Parameters:
518: + tao - the Tao solver context
519: - subsolver - the `Tao` sub-solver context
521: Level: advanced
523: .seealso: `Tao`, `Mat`, `TAOBRGN`
524: @*/
525: PetscErrorCode TaoBRGNGetSubsolver(Tao tao, Tao *subsolver)
526: {
527: PetscFunctionBegin;
529: PetscUseMethod((PetscObject)tao, "TaoBRGNGetSubsolver_C", (Tao, Tao *), (tao, subsolver));
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: static PetscErrorCode TaoBRGNGetSubsolver_BRGN(Tao tao, Tao *subsolver)
534: {
535: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
537: PetscFunctionBegin;
538: *subsolver = gn->subsolver;
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: /*@
543: TaoBRGNSetRegularizerWeight - Set the regularizer weight for the Gauss-Newton least-squares algorithm
545: Collective
547: Input Parameters:
548: + tao - the `Tao` solver context
549: - lambda - L1-norm regularizer weight
551: Level: beginner
553: .seealso: `Tao`, `Mat`, `TAOBRGN`
554: @*/
555: PetscErrorCode TaoBRGNSetRegularizerWeight(Tao tao, PetscReal lambda)
556: {
557: PetscFunctionBegin;
560: PetscTryMethod((PetscObject)tao, "TaoBRGNSetRegularizerWeight_C", (Tao, PetscReal), (tao, lambda));
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: static PetscErrorCode TaoBRGNSetRegularizerWeight_BRGN(Tao tao, PetscReal lambda)
565: {
566: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
568: PetscFunctionBegin;
569: gn->lambda = lambda;
570: PetscFunctionReturn(PETSC_SUCCESS);
571: }
573: /*@
574: TaoBRGNSetL1SmoothEpsilon - Set the L1-norm smooth approximation parameter for L1-regularized least-squares algorithm
576: Collective
578: Input Parameters:
579: + tao - the `Tao` solver context
580: - epsilon - L1-norm smooth approximation parameter
582: Level: advanced
584: .seealso: `Tao`, `Mat`, `TAOBRGN`
585: @*/
586: PetscErrorCode TaoBRGNSetL1SmoothEpsilon(Tao tao, PetscReal epsilon)
587: {
588: PetscFunctionBegin;
591: PetscTryMethod((PetscObject)tao, "TaoBRGNSetL1SmoothEpsilon_C", (Tao, PetscReal), (tao, epsilon));
592: PetscFunctionReturn(PETSC_SUCCESS);
593: }
595: static PetscErrorCode TaoBRGNSetL1SmoothEpsilon_BRGN(Tao tao, PetscReal epsilon)
596: {
597: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
599: PetscFunctionBegin;
600: gn->epsilon = epsilon;
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /*@
605: TaoBRGNSetDictionaryMatrix - bind the dictionary matrix from user application context to gn->D, for compressed sensing (with least-squares problem)
607: Input Parameters:
608: + tao - the `Tao` context
609: - dict - the user specified dictionary matrix. We allow to set a `NULL` dictionary, which means identity matrix by default
611: Level: advanced
613: .seealso: `Tao`, `Mat`, `TAOBRGN`
614: @*/
615: PetscErrorCode TaoBRGNSetDictionaryMatrix(Tao tao, Mat dict)
616: {
617: PetscFunctionBegin;
619: PetscTryMethod((PetscObject)tao, "TaoBRGNSetDictionaryMatrix_C", (Tao, Mat), (tao, dict));
620: PetscFunctionReturn(PETSC_SUCCESS);
621: }
623: static PetscErrorCode TaoBRGNSetDictionaryMatrix_BRGN(Tao tao, Mat dict)
624: {
625: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
627: PetscFunctionBegin;
628: if (dict) {
630: PetscCheckSameComm(tao, 1, dict, 2);
631: PetscCall(PetscObjectReference((PetscObject)dict));
632: }
633: PetscCall(MatDestroy(&gn->D));
634: gn->D = dict;
635: PetscFunctionReturn(PETSC_SUCCESS);
636: }
638: /*@C
639: TaoBRGNSetRegularizerObjectiveAndGradientRoutine - Sets the user-defined regularizer call-back
640: function into the algorithm.
642: Input Parameters:
643: + tao - the Tao context
644: . func - function pointer for the regularizer value and gradient evaluation
645: - ctx - user context for the regularizer
647: Calling sequence:
648: + tao - the `Tao` context
649: . u - the location at which to compute the objective and gradient
650: . val - location to store objective function value
651: . g - location to store gradient
652: - ctx - user context for the regularizer Hessian
654: Level: advanced
656: .seealso: `Tao`, `Mat`, `TAOBRGN`
657: @*/
658: PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine(Tao tao, PetscErrorCode (*func)(Tao tao, Vec u, PetscReal *val, Vec g, PetscCtx ctx), PetscCtx ctx)
659: {
660: PetscFunctionBegin;
662: PetscTryMethod((PetscObject)tao, "TaoBRGNSetRegularizerObjectiveAndGradientRoutine_C", (Tao, PetscErrorCode (*)(Tao, Vec, PetscReal *, Vec, void *), void *), (tao, func, ctx));
663: PetscFunctionReturn(PETSC_SUCCESS);
664: }
666: static PetscErrorCode TaoBRGNSetRegularizerObjectiveAndGradientRoutine_BRGN(Tao tao, PetscErrorCode (*func)(Tao tao, Vec u, PetscReal *val, Vec g, PetscCtx ctx), PetscCtx ctx)
667: {
668: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
670: PetscFunctionBegin;
671: if (ctx) gn->reg_obj_ctx = ctx;
672: if (func) gn->regularizerobjandgrad = func;
673: PetscFunctionReturn(PETSC_SUCCESS);
674: }
676: /*@C
677: TaoBRGNSetRegularizerHessianRoutine - Sets the user-defined regularizer call-back
678: function into the algorithm.
680: Input Parameters:
681: + tao - the `Tao` context
682: . Hreg - user-created matrix for the Hessian of the regularization term
683: . func - function pointer for the regularizer Hessian evaluation
684: - ctx - user context for the regularizer Hessian
686: Calling sequence:
687: + tao - the `Tao` context
688: . u - the location at which to compute the Hessian
689: . Hreg - user-created matrix for the Hessian of the regularization term
690: - ctx - user context for the regularizer Hessian
692: Level: advanced
694: .seealso: `Tao`, `Mat`, `TAOBRGN`
695: @*/
696: PetscErrorCode TaoBRGNSetRegularizerHessianRoutine(Tao tao, Mat Hreg, PetscErrorCode (*func)(Tao tao, Vec u, Mat Hreg, PetscCtx ctx), PetscCtx ctx)
697: {
698: PetscFunctionBegin;
700: PetscTryMethod((PetscObject)tao, "TaoBRGNSetRegularizerHessianRoutine_C", (Tao, Mat, PetscErrorCode (*)(Tao, Vec, Mat, void *), void *), (tao, Hreg, func, ctx));
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
704: static PetscErrorCode TaoBRGNSetRegularizerHessianRoutine_BRGN(Tao tao, Mat Hreg, PetscErrorCode (*func)(Tao tao, Vec u, Mat Hreg, PetscCtx ctx), PetscCtx ctx)
705: {
706: TAO_BRGN *gn = (TAO_BRGN *)tao->data;
708: PetscFunctionBegin;
709: if (Hreg) {
711: PetscCheckSameComm(tao, 1, Hreg, 2);
712: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONG, "NULL Hessian detected! User must provide valid Hessian for the regularizer.");
713: if (ctx) gn->reg_hess_ctx = ctx;
714: if (func) gn->regularizerhessian = func;
715: if (Hreg) {
716: PetscCall(PetscObjectReference((PetscObject)Hreg));
717: PetscCall(MatDestroy(&gn->Hreg));
718: gn->Hreg = Hreg;
719: }
720: PetscFunctionReturn(PETSC_SUCCESS);
721: }
723: /*MC
724: TAOBRGN - Bounded Regularized Gauss-Newton method for solving nonlinear least-squares
725: problems with bound constraints. This algorithm is a thin wrapper around `TAOBNTL`
726: that constructs the Gauss-Newton problem with the user-provided least-squares
727: residual and Jacobian. The algorithm offers an L2-norm ("l2pure"), L2-norm proximal point ("l2prox")
728: regularizer, and L1-norm dictionary regularizer ("l1dict"), where we approximate the
729: L1-norm ||x||_1 by sum_i(sqrt(x_i^2+epsilon^2)-epsilon) with a small positive number epsilon.
730: Also offered is the "lm" regularizer which uses a scaled diagonal of J^T J.
731: With the "lm" regularizer, `TAOBRGN` is a Levenberg-Marquardt optimizer.
732: The user can also provide own regularization function.
734: Options Database Keys:
735: + -tao_brgn_regularization_type - regularization type ("user", "l2prox", "l2pure", "l1dict", "lm") (default "l2prox")
736: . -tao_brgn_regularizer_weight - regularizer weight (default 1e-4)
737: - -tao_brgn_l1_smooth_epsilon - L1-norm smooth approximation parameter: ||x||_1 = sum(sqrt(x.^2+epsilon^2)-epsilon) (default 1e-6)
739: Level: beginner
741: .seealso: `Tao`, `TaoBRGNGetSubsolver()`, `TaoBRGNSetRegularizerWeight()`, `TaoBRGNSetL1SmoothEpsilon()`, `TaoBRGNSetDictionaryMatrix()`,
742: `TaoBRGNSetRegularizerObjectiveAndGradientRoutine()`, `TaoBRGNSetRegularizerHessianRoutine()`
743: M*/
744: PETSC_EXTERN PetscErrorCode TaoCreate_BRGN(Tao tao)
745: {
746: TAO_BRGN *gn;
747: const char *prefix;
749: PetscFunctionBegin;
750: PetscCall(PetscNew(&gn));
752: tao->ops->destroy = TaoDestroy_BRGN;
753: tao->ops->setup = TaoSetUp_BRGN;
754: tao->ops->setfromoptions = TaoSetFromOptions_BRGN;
755: tao->ops->view = TaoView_BRGN;
756: tao->ops->solve = TaoSolve_BRGN;
757: tao->uses_gradient = PETSC_TRUE;
759: PetscCall(TaoParametersInitialize(tao));
761: tao->data = gn;
762: gn->reg_type = TAOBRGN_REGULARIZATION_L2PROX;
763: gn->lambda = 1e-4;
764: gn->epsilon = 1e-6;
765: gn->downhill_lambda_change = 1. / 5.;
766: gn->uphill_lambda_change = 1.5;
767: gn->parent = tao;
769: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)tao, &prefix));
770: PetscCall(TaoCreate(PetscObjectComm((PetscObject)tao), &gn->subsolver));
771: PetscCall(TaoSetType(gn->subsolver, TAOBNLS));
772: PetscCall(TaoSetOptionsPrefix(gn->subsolver, prefix));
773: PetscCall(TaoAppendOptionsPrefix(gn->subsolver, "tao_brgn_subsolver_"));
774: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetRegularizationType_C", TaoBRGNGetRegularizationType_BRGN));
775: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizationType_C", TaoBRGNSetRegularizationType_BRGN));
776: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetDampingVector_C", TaoBRGNGetDampingVector_BRGN));
777: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetDictionaryMatrix_C", TaoBRGNSetDictionaryMatrix_BRGN));
778: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNGetSubsolver_C", TaoBRGNGetSubsolver_BRGN));
779: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerWeight_C", TaoBRGNSetRegularizerWeight_BRGN));
780: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetL1SmoothEpsilon_C", TaoBRGNSetL1SmoothEpsilon_BRGN));
781: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerObjectiveAndGradientRoutine_C", TaoBRGNSetRegularizerObjectiveAndGradientRoutine_BRGN));
782: PetscCall(PetscObjectComposeFunction((PetscObject)tao, "TaoBRGNSetRegularizerHessianRoutine_C", TaoBRGNSetRegularizerHessianRoutine_BRGN));
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }