Actual source code: taoterml1.c
1: #include <petsc/private/taoimpl.h>
3: typedef struct _n_TaoTerm_L1 TaoTerm_L1;
5: struct _n_TaoTerm_L1 {
6: PetscReal epsilon;
7: PetscBool epsilon_warning;
8: Vec diff; // caches $x - p$
9: Vec d; // caches $d_i = sqrt((diff_i)^2 + epsilon^2)
10: PetscReal d_epsilon;
11: PetscReal d_sum;
12: PetscObjectId d_x_id, d_p_id;
13: PetscObjectState d_x_state, d_p_state;
14: Vec diag;
15: PetscObjectId diag_x_id, diag_p_id;
16: PetscObjectState diag_x_state, diag_p_state;
17: PetscReal diag_epsilon;
18: };
20: static PetscErrorCode TaoTermDestroy_L1(TaoTerm term)
21: {
22: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
24: PetscFunctionBegin;
25: PetscCall(VecDestroy(&l1->diff));
26: PetscCall(VecDestroy(&l1->d));
27: PetscCall(VecDestroy(&l1->diag));
28: PetscCall(PetscFree(l1));
29: term->data = NULL;
30: PetscCall(TaoTermDestroy_ElementwiseDivergence_Internal(term));
31: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermL1SetEpsilon_C", NULL));
32: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermL1GetEpsilon_C", NULL));
33: PetscFunctionReturn(PETSC_SUCCESS);
34: }
36: static PetscErrorCode TaoTermL1ComputeData(TaoTerm term, Vec x, Vec params, Vec *_diff, Vec *d)
37: {
38: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
39: PetscObjectId x_id, p_id = 0;
40: PetscObjectState x_state, p_state = 0;
41: Vec diff = x;
43: PetscFunctionBegin;
44: PetscCall(PetscObjectGetId((PetscObject)x, &x_id));
45: PetscCall(PetscObjectStateGet((PetscObject)x, &x_state));
46: if (params) {
47: PetscCall(PetscObjectGetId((PetscObject)params, &p_id));
48: PetscCall(PetscObjectStateGet((PetscObject)params, &p_state));
49: }
50: if (l1->d_x_id != x_id || l1->d_x_state != x_state || l1->d_p_id != p_id || l1->d_p_state != p_state || l1->d_epsilon != l1->epsilon) {
51: l1->d_x_id = x_id;
52: l1->d_x_state = x_state;
53: l1->d_p_id = p_id;
54: l1->d_p_state = p_state;
55: l1->d_epsilon = l1->epsilon;
56: diff = x;
57: if (params) {
58: PetscCall(VecIfNotCongruentGetSameLayoutVec(x, &l1->diff));
59: PetscCall(VecWAXPY(l1->diff, -1.0, params, x));
60: diff = l1->diff;
61: }
62: if (l1->epsilon != 0.0) {
63: PetscCall(VecIfNotCongruentGetSameLayoutVec(x, &l1->d));
64: PetscCall(VecPointwiseMult(l1->d, diff, diff));
65: PetscCall(VecShift(l1->d, l1->epsilon * l1->epsilon));
66: PetscCall(VecSqrtAbs(l1->d));
67: }
68: }
69: if (params) diff = l1->diff;
70: *_diff = diff;
71: *d = l1->d;
72: PetscFunctionReturn(PETSC_SUCCESS);
73: }
75: static PetscErrorCode TaoTermL1ComputeDiag(TaoTerm term, Vec x, Vec params, Vec *diag)
76: {
77: Vec diff, d;
78: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
79: PetscObjectId x_id, p_id = 0;
80: PetscObjectState x_state, p_state = 0;
82: PetscFunctionBegin;
83: PetscCall(TaoTermL1ComputeData(term, x, params, &diff, &d));
84: PetscCall(PetscObjectGetId((PetscObject)x, &x_id));
85: PetscCall(PetscObjectStateGet((PetscObject)x, &x_state));
86: if (params) {
87: PetscCall(PetscObjectGetId((PetscObject)params, &p_id));
88: PetscCall(PetscObjectStateGet((PetscObject)params, &p_state));
89: }
90: if (l1->diag_x_id != x_id || l1->diag_x_state != x_state || l1->diag_p_id != p_id || l1->diag_p_state != p_state || l1->diag_epsilon != l1->epsilon) {
91: l1->diag_x_id = x_id;
92: l1->diag_x_state = x_state;
93: l1->diag_p_id = p_id;
94: l1->diag_p_state = p_state;
95: l1->diag_epsilon = l1->epsilon;
96: if (l1->epsilon != 0.0) {
97: PetscCall(VecIfNotCongruentGetSameLayoutVec(x, &l1->diag));
98: PetscCall(VecCopy(d, l1->diag));
99: PetscCall(VecPointwiseMult(l1->diag, l1->diag, d));
100: PetscCall(VecPointwiseMult(l1->diag, l1->diag, d));
101: PetscCall(VecReciprocal(l1->diag));
102: PetscCall(VecScale(l1->diag, l1->epsilon * l1->epsilon));
103: }
104: }
105: *diag = l1->diag;
106: PetscFunctionReturn(PETSC_SUCCESS);
107: }
109: static PetscErrorCode TaoTermComputeObjective_L1_Internal(TaoTerm term, Vec diff, Vec d, PetscReal *value)
110: {
111: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
113: PetscFunctionBegin;
114: if (l1->epsilon == 0.0) {
115: PetscCall(VecNorm(diff, NORM_1, value));
116: } else {
117: PetscScalar sum;
118: PetscInt n;
120: PetscCall(VecGetSize(d, &n));
121: PetscCall(VecSum(d, &sum));
122: *value = PetscRealPart(sum) - n * l1->epsilon;
123: }
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode TaoTermComputeObjective_L1(TaoTerm term, Vec x, Vec params, PetscReal *value)
128: {
129: Vec diff, d;
131: PetscFunctionBegin;
132: PetscCall(TaoTermL1ComputeData(term, x, params, &diff, &d));
133: PetscCall(TaoTermComputeObjective_L1_Internal(term, diff, d, value));
134: PetscFunctionReturn(PETSC_SUCCESS);
135: }
137: static PetscErrorCode TaoTermL1DerivativeCheck(TaoTerm term)
138: {
139: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
141: PetscFunctionBegin;
142: if (l1->epsilon_warning == PETSC_FALSE) {
143: l1->epsilon_warning = PETSC_TRUE;
144: PetscCall(PetscInfo(term, "%s: Asking for derivatives of l1 norm, which is not smooth. Consider smoothing the TaoTerm with TaoTermL1SetEpsilon() or using a derivative-free Tao solver\n", ((PetscObject)term)->prefix));
145: }
146: PetscFunctionReturn(PETSC_SUCCESS);
147: }
149: static PetscErrorCode TaoTermComputeGradient_L1_Internal(TaoTerm term, Vec diff, Vec d, Vec g)
150: {
151: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
153: PetscFunctionBegin;
154: if (l1->epsilon == 0.0) {
155: PetscCall(TaoTermL1DerivativeCheck(term));
156: PetscCall(VecPointwiseSign(g, diff, VEC_SIGN_ZERO_TO_ZERO));
157: } else {
158: PetscCall(VecPointwiseDivide(g, diff, d));
159: }
160: PetscFunctionReturn(PETSC_SUCCESS);
161: }
163: static PetscErrorCode TaoTermComputeGradient_L1(TaoTerm term, Vec x, Vec params, Vec g)
164: {
165: Vec diff, d;
167: PetscFunctionBegin;
168: PetscCall(TaoTermL1ComputeData(term, x, params, &diff, &d));
169: PetscCall(TaoTermComputeGradient_L1_Internal(term, diff, d, g));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode TaoTermComputeObjectiveAndGradient_L1(TaoTerm term, Vec x, Vec params, PetscReal *value, Vec g)
174: {
175: Vec diff, d;
177: PetscFunctionBegin;
178: PetscCall(TaoTermL1ComputeData(term, x, params, &diff, &d));
179: PetscCall(TaoTermComputeObjective_L1_Internal(term, diff, d, value));
180: PetscCall(TaoTermComputeGradient_L1_Internal(term, diff, d, g));
181: PetscFunctionReturn(PETSC_SUCCESS);
182: }
184: static PetscErrorCode TaoTermComputeHessian_L1_Internal(TaoTerm term, Vec diag, Mat H)
185: {
186: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
188: PetscFunctionBegin;
189: if (l1->epsilon == 0.0) {
190: PetscCall(TaoTermL1DerivativeCheck(term));
191: PetscCall(MatZeroEntries(H));
192: PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
193: PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
194: } else {
195: PetscCall(MatDiagonalSet(H, diag, INSERT_VALUES));
196: }
197: PetscFunctionReturn(PETSC_SUCCESS);
198: }
200: static PetscErrorCode TaoTermComputeHessian_L1(TaoTerm term, Vec x, Vec params, Mat H, Mat Hpre)
201: {
202: Vec diag = NULL; /* Appease -Wmaybe-uninitialized */
204: PetscFunctionBegin;
205: if (H == NULL && Hpre == NULL) PetscFunctionReturn(PETSC_SUCCESS);
206: PetscCall(TaoTermL1ComputeDiag(term, x, params, &diag));
207: if (H) PetscCall(TaoTermComputeHessian_L1_Internal(term, diag, H));
208: if (Hpre && Hpre != H) PetscCall(TaoTermComputeHessian_L1_Internal(term, diag, Hpre));
209: PetscFunctionReturn(PETSC_SUCCESS);
210: }
212: static PetscErrorCode TaoTermCreateHessianMatrices_L1(TaoTerm term, Mat *H, Mat *Hpre)
213: {
214: PetscBool is_hdiag, is_hprediag;
216: PetscFunctionBegin;
217: PetscCall(PetscInfo(term, "%s: Creating TAOTERML1 Hessian Matrices. TAOTERML1 only accepts MATDIAGONAL for MatType, overriding any user-set MatType.\n", ((PetscObject)term)->prefix));
218: PetscCall(PetscStrcmp(term->H_mattype, MATDIAGONAL, &is_hdiag));
219: PetscCall(PetscStrcmp(term->Hpre_mattype, MATDIAGONAL, &is_hprediag));
220: if (!is_hdiag) {
221: PetscCall(PetscFree(term->H_mattype));
222: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
223: }
224: if (!is_hprediag) {
225: PetscCall(PetscFree(term->Hpre_mattype));
226: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));
227: }
228: PetscCall(TaoTermCreateHessianMatricesDefault(term, H, Hpre));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: /*@
233: TaoTermL1SetEpsilon - Set an $\epsilon$ smoothing parameter.
235: Logically collective
237: Input Parameters:
238: + term - a `TaoTerm` of type `TAOTERML1`
239: - epsilon - a real number $\geq 0$
241: Options Database Keys:
242: . -tao_term_l1_epsilon <real> - $\epsilon$
244: Level: advanced
246: If $\epsilon = 0$ (the default), then `term` computes $\|x - p\|_1$, but if $\epsilon > 0$, then it computes
247: $\sum_{i=0}^{n-1} \left(\sqrt{(x_i-p_i)^2 + \epsilon^2} - \epsilon\right)$.
249: .seealso: [](sec_tao_term),
250: `TaoTerm`,
251: `TAOTERML1`,
252: `TaoTermL1GetEpsilon()`
253: @*/
254: PetscErrorCode TaoTermL1SetEpsilon(TaoTerm term, PetscReal epsilon)
255: {
256: PetscFunctionBegin;
259: PetscTryMethod(term, "TaoTermL1SetEpsilon_C", (TaoTerm, PetscReal), (term, epsilon));
260: PetscFunctionReturn(PETSC_SUCCESS);
261: }
263: static PetscErrorCode TaoTermL1SetEpsilon_L1(TaoTerm term, PetscReal epsilon)
264: {
265: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
267: PetscFunctionBegin;
268: PetscCheck(epsilon >= 0, PetscObjectComm((PetscObject)term), PETSC_ERR_ARG_OUTOFRANGE, "L1 epsilon (%g) cannot be < 0.0", (double)epsilon);
269: l1->epsilon = epsilon;
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: /*@
274: TaoTermL1GetEpsilon - Get the $\epsilon$ smoothing parameter set by `TaoTermL1SetEpsilon()`.
276: Not collective
278: Input Parameter:
279: . term - a `TaoTerm` of type `TAOTERML1`
281: Output Parameter:
282: . epsilon - the smoothing parameter
284: Level: advanced
286: .seealso: [](sec_tao_term),
287: `TaoTerm`,
288: `TAOTERML1`,
289: `TaoTermL1SetEpsilon()`
290: @*/
291: PetscErrorCode TaoTermL1GetEpsilon(TaoTerm term, PetscReal *epsilon)
292: {
293: PetscFunctionBegin;
295: PetscAssertPointer(epsilon, 2);
296: PetscUseMethod(term, "TaoTermL1GetEpsilon_C", (TaoTerm, PetscReal *), (term, epsilon));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode TaoTermL1GetEpsilon_L1(TaoTerm term, PetscReal *epsilon)
301: {
302: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
304: PetscFunctionBegin;
305: *epsilon = l1->epsilon;
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: static PetscErrorCode TaoTermView_L1(TaoTerm term, PetscViewer viewer)
310: {
311: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
312: PetscBool is_ascii;
314: PetscFunctionBegin;
315: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &is_ascii));
316: if (is_ascii) PetscCall(PetscViewerASCIIPrintf(viewer, "epsilon (tao_term_l1_epsilon): %g\n", (double)l1->epsilon));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }
320: static PetscErrorCode TaoTermSetFromOptions_L1(TaoTerm term, PetscOptionItems PetscOptionsObject)
321: {
322: TaoTerm_L1 *l1 = (TaoTerm_L1 *)term->data;
323: PetscBool is_hdiag, is_hprediag;
325: PetscFunctionBegin;
326: PetscOptionsHeadBegin(PetscOptionsObject, "TaoTerm l1 options");
327: PetscCall(PetscOptionsBoundedReal("-tao_term_l1_epsilon", "smoothing parameter", "TaoTermL1SetEpsilon", l1->epsilon, &l1->epsilon, NULL, 0.0));
328: PetscOptionsHeadEnd();
329: PetscCall(PetscStrcmp(term->H_mattype, MATDIAGONAL, &is_hdiag));
330: PetscCall(PetscStrcmp(term->Hpre_mattype, MATDIAGONAL, &is_hprediag));
331: if (!is_hdiag) {
332: PetscCall(PetscFree(term->H_mattype));
333: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
334: }
335: if (!is_hprediag) {
336: PetscCall(PetscFree(term->Hpre_mattype));
337: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));
338: }
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode TaoTermIsComputeHessianFDPossible_L1(TaoTerm term, PetscBool3 *ispossible)
343: {
344: PetscFunctionBegin;
345: *ispossible = PETSC_BOOL3_FALSE;
346: PetscFunctionReturn(PETSC_SUCCESS);
347: }
349: /*MC
350: TAOTERML1 - A `TaoTerm` that computes $\|x - p\|_1$, for solution $x$ and parameters $p$.
352: Level: intermediate
354: Options Database Keys:
355: . -tao_term_l1_epsilon <real> - (default 0.0) a smoothing parameter (see `TaoTermL1SetEpsilon()`)
357: Notes:
358: This term is `TAOTERM_PARAMETERS_OPTIONAL`. If the parameters argument is `NULL` for
359: evaluation routines the term computes $\|x\|_1$.
361: This term has a smoothing parameter $\epsilon$ that defaults to 0: if $\epsilon > 0$,
362: the term computes a smooth approximation of $\|x - p\|_1$, see `TaoTermL1SetEpsilon()`.
364: The default Hessian creation mode (see `TaoTermGetCreateHessianMode()`) is `H == Hpre` and `TaoTermCreateHessianMatrices()`
365: will create a `MATDIAGONAL` for the Hessian.
367: .seealso: [](sec_tao_term),
368: `TaoTerm`,
369: `TaoTermType`,
370: `TaoTermCreateL1()`,
371: `TaoTermL1GetEpsilon()`,
372: `TaoTermL1SetEpsilon()`,
373: `TAOTERMHALFL2SQUARED`,
374: `TAOTERMQUADRATIC`
375: M*/
376: PETSC_INTERN PetscErrorCode TaoTermCreate_L1(TaoTerm term)
377: {
378: TaoTerm_L1 *l1;
380: PetscFunctionBegin;
381: PetscCall(TaoTermCreate_ElementwiseDivergence_Internal(term));
382: PetscCall(PetscNew(&l1));
383: term->data = (void *)l1;
385: PetscCall(PetscFree(term->H_mattype));
386: PetscCall(PetscFree(term->Hpre_mattype));
388: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->H_mattype));
389: PetscCall(PetscStrallocpy(MATDIAGONAL, (char **)&term->Hpre_mattype));
391: term->ops->destroy = TaoTermDestroy_L1;
392: term->ops->view = TaoTermView_L1;
393: term->ops->setfromoptions = TaoTermSetFromOptions_L1;
394: term->ops->objective = TaoTermComputeObjective_L1;
395: term->ops->gradient = TaoTermComputeGradient_L1;
396: term->ops->objectiveandgradient = TaoTermComputeObjectiveAndGradient_L1;
397: term->ops->hessian = TaoTermComputeHessian_L1;
398: term->ops->createhessianmatrices = TaoTermCreateHessianMatrices_L1;
399: term->ops->iscomputehessianfdpossible = TaoTermIsComputeHessianFDPossible_L1;
401: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermL1SetEpsilon_C", TaoTermL1SetEpsilon_L1));
402: PetscCall(PetscObjectComposeFunction((PetscObject)term, "TaoTermL1GetEpsilon_C", TaoTermL1GetEpsilon_L1));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: /*@
407: TaoTermCreateL1 - Create a `TaoTerm` for the objective function term $\|x - p\|_1$.
409: Collective
411: Input Parameters:
412: + comm - the MPI communicator where the term will be computed
413: . n - the local size of the $x$ and $p$ vectors (or `PETSC_DECIDE`)
414: . N - the global size of the $x$ and $p$ vectors (or `PETSC_DECIDE`)
415: - epsilon - a non-negative smoothing parameter (see `TaoTermL1SetEpsilon()`)
417: Output Parameter:
418: . term - the `TaoTerm`
420: Level: beginner
422: Note:
423: If you would like to add an L1 regularization term $\alpha \|x\|_1$ to the objective function of a `Tao`, do the following\:
424: .vb
425: VecGetLocalSize(x, &n);
426: VecGetSize(x, &N);
427: TaoTermCreateL1(PetscObjectComm((PetscObject)x), n, N, 0.0, &term);
428: TaoAddTerm(tao, "reg_", alpha, term, NULL, NULL);
429: TaoTermDestroy(&term);
430: .ve
431: If you would like to have a dictionary matrix term $\alpha \|D x\|_1$, do the same but pass `D` as the map of the term\:
432: .vb
433: MatGetLocalSize(D, &m, NULL);
434: MatGetSize(D, &M, NULL);
435: TaoTermCreateL1(PetscObjectComm((PetscObject)D), m, M, 0.0, &term);
436: TaoAddTerm(tao, "reg_", alpha, term, NULL, D);
437: TaoTermDestroy(&term);
438: .ve
440: .seealso: [](sec_tao_term),
441: `TaoTerm`,
442: `TAOTERML1`,
443: `TaoTermL1GetEpsilon()`,
444: `TaoTermL1SetEpsilon()`,
445: `TaoTermCreateHalfL2Squared()`,
446: `TaoTermCreateQuadratic()`
447: @*/
448: PetscErrorCode TaoTermCreateL1(MPI_Comm comm, PetscInt n, PetscInt N, PetscReal epsilon, TaoTerm *term)
449: {
450: TaoTerm _term;
452: PetscFunctionBegin;
453: PetscAssertPointer(term, 5);
454: PetscCheck(epsilon >= 0, comm, PETSC_ERR_ARG_OUTOFRANGE, "L1 epsilon (%g) cannot be < 0.0", (double)epsilon);
455: PetscCall(TaoTermCreate(comm, &_term));
456: PetscCall(TaoTermSetType(_term, TAOTERML1));
457: PetscCall(TaoTermSetSolutionSizes(_term, n, N, 1));
458: PetscCall(TaoTermL1SetEpsilon(_term, epsilon));
459: *term = _term;
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }