Antonin Chambolle
email: chambolle@ceremade.dauphine.fr

Thomas Pock
email: pock@icg.tugraz.at

Learning Consistent Discretizations of the Total Variation *

Keywords: Total variation, image denoising, image inpainting, discretization, finite differences, finite elements, learning, bi-level optimization, primal-dual algorithms. AMS subject classifications. 49Q20 35J87 35R35 49M29 94A08 65D18 65N06 65N30 65N50

In this work, we study a general framework of discrete approximations of the total variation for image reconstruction problems. The framework, for which we can show consistency in the sense of Γ-convergence, unifies and extends several existing discretization schemes. In addition, we propose algorithms for learning discretizations of the total variation in order to achieve the best possible reconstruction quality for particular image reconstruction tasks. Interestingly, the learned discretizations significantly differ between the tasks, illustrating that there is no universal best discretization of the total variation.

Introduction.

This paper tries to explore some questions raised in a previous work [START_REF] Chambolle | Crouzeix-Raviart approximation of the total variation on simplicial meshes[END_REF] (see also [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF]), were the authors considered alternative discretizations on the total variation, based on finite elements, for image denoising or (basic) inpainting. In the latter experimental setting, the issue is to recover a missing discontinuity from its boundary datum, and it was observed that even when a discrete approximation is enjoying the property that the discretization of a sharp edge is, whatever the direction, a minimizer of the discrete inpainting problem (as in the case for the line labeled "RT" in Fig. 2.2, see [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF] and [START_REF] Chambolle | Finite differences and finite elements discretizations of the total variation[END_REF]Prop. 4.1]), results could substantially differ in terms of sharpness and precision. In the series of experiments carried on in [START_REF] Chambolle | Crouzeix-Raviart approximation of the total variation on simplicial meshes[END_REF][START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF], the best results for inpainting discontinuities were obtained with a discretization introduced in [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF] and properly implemented in [START_REF] Condat | Discrete total variation: new definition and minimization[END_REF], called the "Condat" discretization. Mimicking the standard definition (by duality) of the total variation (cf [16, §1.1], here to simplify, on a bidimensional bounded open domain Ω ⊂ R 2 ):

(1.1) T V (u) := sup - Ω u div p dx : p ∈ C ∞ c (Ω; R 2 ), |p(x)| ≤ 1 ∀x ∈ Ω
for a function u ∈ L 1 (Ω) (we address only the scalar problem, valid for instance for gray-level images), the idea of [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF][START_REF] Condat | Discrete total variation: new definition and minimization[END_REF] is to define a discrete total variation of a matrix of pixel values (u i,j ) as sup i,j p 1 i+ 1 2 ,j (u i+1,j -u i,j )+p 2 i,j+ 1 2 (u i,j+1 -u i,j ), with a particular set of constraints on the discrete dual variables p = (p 1 i+ 1 2 ,j , p 2 i,j+ 1 2 ). While the usual forward-differences definition, such as considered in [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF], simply requires the constraints (p 1

i+ 1 2 ,j ) 2 + (p 2 i,j+ 1 2
) 2 ≤ 1 for all i, j, the idea of [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF][START_REF] Condat | Discrete total variation: new definition and minimization[END_REF] (see also the variants [START_REF] Destuynder | A dual algorithm for denoising and preserving edges in image processing[END_REF][START_REF] Herrmann | Discrete total variation with finite elements and applications to imaging[END_REF]) is to linearly interpolate the discrete dual variables p = (p 1 i+ 1 2 ,j , p 2 i,j+ 1 2

) into a continuous field p(x) = (p 1 (x), p 2 (x)) and enforce norm constraints for this continuous field in a few well chosen discrete locations.

A first question which is raised by the new definition in [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF][START_REF] Condat | Discrete total variation: new definition and minimization[END_REF] is whether it is consistent, that is, whether this it can be truly regarded as discretization of the total variation functional 1.1. Strangely, the proof is not as obvious as for more standard discretizations. The most rudimentary framework to assess this consistency is the so-called "Γ-convergence" [START_REF] Maso | An introduction to Γ-convergence[END_REF][START_REF] Braides | Γ-convergence for beginners[END_REF], a property which ensures that minimizers of the discrete total variation (plus various types of lower order terms) will converge, as the pixels' size goes to zero, to a minimizer of the continuous one (1.1). This does not give any convergence rate or error bounds, but at least guarantees the consistency of the approach. Error bounds for variants have been discussed in various papers [START_REF] Wang | Error bounds for finite-difference methods for Rudin-Osher-Fatemi image smoothing[END_REF][START_REF] Lai | The Convergence of a Central-Difference Discretization of Rudin-Osher-Fatemi Model for Image Denoising[END_REF][START_REF] Chambolle | Crouzeix-Raviart approximation of the total variation on simplicial meshes[END_REF] (see also the review [START_REF] Chambolle | Finite differences and finite elements discretizations of the total variation[END_REF]) while in [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF] the second author, together with Corentin Caillaud, report an essentially optimal error bound for the "Rudin-Osher-Fatemi" (that is, denoising) problem based on a variant of [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF][START_REF] Condat | Discrete total variation: new definition and minimization[END_REF] where the constraints on the dual fields p(x) are enforced at each "corner" of a pixel, and which corresponds to considering in (1.1) a subset of the admissible fields, namely "Raviart-Thomas" dual fields (which are conforming finite element approximations of fields with divergence [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF]) subject to the square pixels' grid. Up to now, similar error bounds seem quite hard to establish for [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF][START_REF] Condat | Discrete total variation: new definition and minimization[END_REF], we observe however that both the latter and the generalizations we will consider here, obtained by minimizing some empirical error criterion, should in practice behave better than the variant studied in [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF].

In this work, we start by addressing the consistency issue, in a more general framework. Our idea is that the previous approaches are easily generalized, by enforcing constraints on particular averages of the dual fields p(x) (which might correspond to interpolants in particular points as in [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF][START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF][START_REF] Condat | Discrete total variation: new definition and minimization[END_REF], but more generally are obtained by convolutions with kernels of small support). It follows a quite general definition, which also includes as a particular case the classical definition [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]. We show that this definition is consistent in the sense of Γ-convergence (Theorem 2.4).

In a next step, we address the issue of learning the parameters of the model, for different types of tasks. It means that we can try to find the "best" convolution kernels so that the denoising or inpainting problems are "best" solved for a given set of data. The idea is that (heuristically) the error of these kernels should always be better than the error bounds for the standard approaches which fit the framework. We show that such an approach can lead to surprisingly good results for the given tasks, at the expense though of quite more complicated discrete total variations than the ones usually considered (as we replace the standard dual constraint by several -up to 8 in most of our examples and 40 for the denoising of natural images -constraints involving a whole neighborhood of pixels).

In a near future, we believe this framework can be adapted to address more complex energies of singular objects in 2D or 3D, and we plan to extend it for solving harder problems, such as the "total roto-translational variation" described in [START_REF] Chambolle | Total roto-translational variation[END_REF].

The plan of the paper is as follows: in the next Sections 2-2.3 we describe more precisely our general setting. Then in Section 2.4 we state and prove the consistency result. In the second part of the paper, starting with Section 3, we address the learning of the parameters. We describe our strategy and, in Section 4, we apply it to several tasks such as the "denois- ing" and "inpainting problems", which we describe precisely in the corresponding sections. Eventually in the appendix, we (partially) justify the differentiation technique which we have used for optimizing the parameters of our models.

Acknowledgements. The authors would like to thank the referees, for their very accurate reading of the paper and their comments, which have helped improving the final version. Most of this work was done while A.C. was still with CMAP, CNRS and Ecole Polytechnique, Institut Polytechnique de Paris, Palaiseau, France.

Problem and setting.

In this section, we define the general discretization framework and show that several popular discretization schemes are particular instances of this class of discretizations.

An image u ∈ X R M ×N is a matrix u = (u i,j ) 1≤i≤M,1≤j≤N of pixel values (here to simplify we consider scalar images). We introduce the operator D : X → Y R (M -1)×N × R M ×(N -1) defined as usual as Du = (D 1 u, D 2 u) with (2.1)

   (D 1 u) i+ 1 2 ,j = u i+1,j -u i,j i = 1, . . . , M -1, j = 1, . . . , N ; (D 2 u) i,j+ 1 2 = u i,j+1 -u i,j i = 1, . . . , M, j = 1, . . . , N -1 ;
A natural way to discretize the total variation consists in introducing an appropriate norm • Y on Y , typically some sort of "2, 1" norm, and let Du Y be the approximate total variation.

The norm Du Y can of course be defined by duality as

Du Y = sup { p, Du Y : p * Y ≤ 1}
(where • * Y is the dual (or polar) norm of • Y , defined in theory on the dual Y * which here is identified with Y through a Euclidean scalar product •, • Y ). As illustrated in Figure 2.1, the discrete image pixel values u i,j can be thought as pixel-averaged values of an underlying continuous function u, and are localized in the center of the pixels. The discrete dual variables p = (p 1 , p 2 ) with values p 1 i+ 1 2 ,j on vertical edges and p 2 i,j+ 1 2 on horizontal edges can be interpreted as edge-averaged fluxes of an underlying continuous field p(x), localized in the center of the edges.

The most standard choice is, for p ∈ Y , to let:

(2.2) p Y := M i=1 N j=1 (p 1 i+ 1 2 ,j ) 2 + (p 2 i,j+ 1 2 ) 2
where by convention one sets in the norm p 1

M + 1 2 ,j = p 2 i,N + 1 2
= 0 for all i, j. This is known not to give very precise results and suffer anisotropy issues, see [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF] for a recent study. Then, an improvement over (2.2), suggested in [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF][START_REF] Condat | Discrete total variation: new definition and minimization[END_REF], consists in replacing Du Y with another norm defined by means of an averaging operator F : Y → Z, where Z is another Euclidean space, and constraining in the dual formulation above F p rather than p:

(2.3) T V (u) := sup { p, Du Y : F p * Z ≤ 1}
In practice, in this paper, F is of convolution type and Z has more or less the same structure as Y (and in general is a copy of L folds of Y , L ≥ 1). As considered in [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF]Sec. 5], this point of view is particularly natural in the setting where (to simplify, as later in Section 2.4) we have Ω a rectangle, ε > 0 a small parameter with N ∼ M ∼ 1/ε, and u is identified to the piecewise constant "Q0" function

(2.4) i,j u i,j χ (iε-ε,iε)×(jε-ε,jε) (x).
Then, extending the discrete variable p as a Raviart-Thomas field [START_REF] Raviart | A mixed finite element method for 2nd order elliptic problems[END_REF] p(x) subject to the square mesh i,j (iε -ε, iε) × (jε -ε, jε), and with vanishing flux across ∂Ω, the scaled expression ε p, Du Y is precisely equal to -Ω div p(x)u(x)dx, and one expects to recover a sound approximation of (1.1) if the constraint

F p * Z ≤ 1 in (2.
3) is a good approximation of the dual constraint in (1.1) for the field p(x). (We observe also that the expression -Ω div p(x)u(x)dx is well defined for any u ∈ L 1 (Ω), and depends only on the (L 2 ) projection of u on piecewise constant functions of the form (2.4), consisting in replacing u on each small pixel (iε -ε, iε) × (jε -ε, jε) by its average on the pixel as already suggested.)

An obvious interesting case, described in [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF] and which enters our general framework, is when F is chosen in such a way that the discrete constraint in (2.3) is precisely equivalent to |p(x)| ≤ 1 a.e. in Ω for the Raviart-Thomas field p(x), yielding a functional which clearly, for any u ∈ L 1 (Ω), is below (1.1) and can be shown to be a good approximation of the total variation in terms of error bound. It corresponds to the setting described in (2.6)-(2.7) in the next Section 2.1; we refer to [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF]Sec. 5] (cf also [START_REF] Destuynder | A dual algorithm for denoising and preserving edges in image processing[END_REF][START_REF] Chambolle | Finite differences and finite elements discretizations of the total variation[END_REF]), for more details and precise definitions. One objective of this paper is to generalize both [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF] and [START_REF] Condat | Discrete total variation: new definition and minimization[END_REF].

Observe that for r ∈ Y one has:

sup

F p * Z ≤1 p, r Y = sup p inf q p, r Y -q, F p Z + q Z ≤ inf q sup p q Z + p, r -F * q Y = inf q:F * q=r q Z
It is standard that this is an equality (using that the RHS is lower semicontinuous and its convex conjugate is precisely the characteristic of {p : F p * Z ≤ 1}), provided of course the value of the last inf is +∞ when the feasible set is empty, which is when r ∈ imF * = (ker F ) ⊥ .

Hence the discrete total variation (2.3) also has the primal form

(2.5) T V (u) = min q: F * q=Du q Z .
Observe that this formulation of the total variation can also be interpreted as a particular form of a "group-lasso" problem, where one seeks for a sparse representation F * q of the discrete gradient Du based on a convolutional dictionary F * [START_REF] Chambolle | An introduction to continuous optimization for imaging[END_REF]. We believe, that in the spirit of compressed sensing, this formulation enables the recovery of sharp discontinuities in discrete images and breaks the intrinsic band limitations of more classical, digital signal processing inspired discretizations, for example the Shannon total variation proposed in [START_REF] Abergel | The Shannon Total Variation[END_REF].

Some examples.

In the following we mention a few popular discretizations which can all be seen as particular instances of our general discretization framework.

Example 2.1 (Forward differences (FD)). Certainly, the most widely used discretization of the total variation is based on simple forward differences (see for example [START_REF] Chambolle | An algorithm for total variation minimization and applications[END_REF]). In our framework, this corresponds to defining F p ∈ Z Y as

(F p) i,j =   p 1 i+ 1 2 ,j p 2 i,j+ 1 2   .
The corresponding Z-norm is then equivalent to (2.2). As the operator F combines dual variables defined at different spatial locations, it follows that the constraint F p * Z ≤ 1 leads to a bias towards certain directions of the image gradient. This can be clearly observed in Figure 2.2, where some directions appear sharp, while others are overly blurred.

Example 2.2 (Raviart-Thomas discretization of the dual (RT) [START_REF] Destuynder | A dual algorithm for denoising and preserving edges in image processing[END_REF][START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF][START_REF] Chambolle | Finite differences and finite elements discretizations of the total variation[END_REF]). In order to gain more rotational invariance, one can consider four constraints on the dual variables at the four corners of each pixel. In our framework, this corresponds to letting p ∈ Y , F p = (F 1 p, F 2 p, F 3 p, F 4 p) ∈ Z and defining

(2.6) (F 1 p) i-1 2 ,j-1 2 =   p 1 i-1 2 ,j p 2 i,j-1 2   , (F 2 p) i-1 2 ,j+ 1 2 =   p 1 i-1 2 ,j p 2 i,j+ 1 2   , (F 3 p) i+ 1 2 ,j-1 2 =   p 1 i+ 1 2 ,j p 2 i,j-1 2   , (F 4 p) i+ 1 2 ,j+ 1 2 =   p 1 i+ 1 2 ,j p 2 i,j+ 1 2   .
The corresponding Z norm is given by

(2.7) (z 1 , z 2 , z 3 , z 4 ) Z := i,j |z 1 i-1 2 ,j-1 2 | 2 + |z 2 i-1 2 ,j+ 1 2 | 2 + |z 3 i+ 1 2 ,j-1 2 | 2 + |z 4 i+ 1 2 ,j+ 1 2 | 2 ,
where | • | 2 is the standard Euclidean length. The constraint F p * Z ≤ 1 then combines all possible pairs of edges that meet in one of the four corners of each pixel. With this, one ensures feasibility of the Raviart-Thomas interpretation of the dual field as a dual variable in (1.1), as mentioned in the previous section. However, as seen on Figure 2.2, these constraints seem not appropriate for the resolution of edge inpainting problems, as they lead to a strong blur of the inpainted edges with non grid-aligned orientations. This is a bit puzzling, as one can show (in a more precise setting though than the experiments carried on in this paper) that straight line inpainting problems can be solved exactly by the discrete projection of the continuous solution with this particular discretization, see [START_REF] Chambolle | Finite differences and finite elements discretizations of the total variation[END_REF]Prop. 4.1]. It is likely that, as is the case for the variant studied in [START_REF] Chambolle | Crouzeix-Raviart approximation of the total variation on simplicial meshes[END_REF], such problem have non-unique solutions and the optimization, or variants in the setting, favor blurry results over sharper ones. We must observe that the Raviart-Thomas field being admissible in the dual definition of the total variation (more precisely, after some smoothing, in (1.1)), the value of this discrete total variation will always be below the value of the TV of any L 1 function which has average u i,j on each pixel (i, j). In particular, this can lead to an important under-estimation (and therefore unsufficient removal) of a noise oscillating at the same scale as the pixels.

Example 2.3 (Condat's discretization (CD) [START_REF] Hintermüller | Functional-analytic and numerical issues in splitting methods for total variation-based image reconstruction[END_REF][START_REF] Condat | Discrete total variation: new definition and minimization[END_REF]). A more flexible discretization is obtained by linear interpolation of the dual variables at the three grid positions (i, j), (i + 1 2 , j), and (i, j + 1

2 ). This corresponds to letting p ∈ Y , F p = (F 1 p, F 2 p, F 3 p) ∈ Z with, for i = 1, . . . , M , j = 1, . . . , N ,

(F 1 p) i,j =    p 1 i-1 2 ,j +p 1 i+ 1 2 ,j 2 p 2 i,j-1 2 +p 2 i,j+ 1 2 2    ,
(again the convention is that p = 0 when the indices are out of range), and, for i ≤ M -1,

(F 2 p) i+ 1 2 ,j =   p 1 i+ 1 2 ,j p 2 i,j-1 2 +p 2 i,j+ 1 2 +p 2 i+1,j-1 2 +p 2 i+1,j+ 1 2 4   ,
and for j ≤ N -1:

(F 3 p) i,j+ 1 2 =   p 1 i-1 2 ,j +p 1 i+ 1 2 ,j +p 1 i-1 2 ,j+1 +p 1 i+ 1 2 ,j+1 4 p 2 i,j+ 1 2   . The Z norm is then (z 1 , z 2 , z 3 ) Z := i,j |z 1 i,j | 2 + |z 2 i+ 1 2 ,j | 2 + |z 3 i,j+ 1 2 | 2 .
In [START_REF] Condat | Discrete total variation: new definition and minimization[END_REF], another variant (which we call CD4) is proposed which adds another constraint on the grid positions (i + 1 2 , j + 1 2 ). In our framework this corresponds to adding a fourth filter of the form

(F 4 p) i,j =    p 1 i+ 1 2 ,j +p 1 i+ 1 2 ,j+1 2 p 2 i,j+ 1 2 +p 2 i+1,j+ 1 2 2    .
In our experiments (cf. Figure 2.2), the performances of CD4 are comparable to CD, and for simplicity we decided to stick to CD in further experiments. In contrast to RT, the constraint F p * Z ≤ 1 in CD does not guarantee the feasibility of the Raviart-Thomas extension of the discrete dual field, yet it seems this is necessary to better discriminate sharper edges from smoother ones. And indeed, Figure 2.2 confirms the excellent result for inpainting straight edges with various orientations.

Table 2.1 visualizes the filter coefficients of the aforementioned handcrafted methods. Observe that Condat's discretization (CD, CD4) is the only one which relies on some averaging of the dual variables. Figure 2.2 shows the performance of these methods for recovering discontinuities of various orientations.

FD CD CD4 RT

Table 2.1 Visualization of the filters corresponding to the different handcrafted discretization methods. The filter sizes are 2 × 3 for the horizontal direction and 3 × 2 for the vertical direction, which corresponds to a support of 2 × 2 pixels. Black corresponds to filters weights of 0.0, white corresponds to filter weights of 1.0, gray corresponds to filter weights of 0.5 and dark gray corresponds to filter weights of 0.25.

General interpolation operators.

In practice, as in the examples above, we consider general convolution-type operators F . The basic form for such

F is F = (F l ) L l=1 with (2.8) (F l p) i,j = (F l,1 p 1 ) i,j (F l,2 p 2 ) i,j =   ν m,n=-ν ξ l m,n p 1 i+ 1 2 -m,j-n ν m,n=-ν η l m,n p 2 i-m,j+ 1 2 -n   ,
where (ξ l m,n , η l m,n ) are weights with a small support of maximal size (2ν + 1) × (2ν + 1) for some integer ν ≥ 01 . On order to ensure that F l,1 , F l,2 represent valid interpolation kernels, we shall assume that the entries of the filters sum up to one, that is, (2.9)

m,n ξ l m,n = m,n η l m,n = 1 ⇐⇒ F l,1 , F l,2 ∈ C Σ=1 ,
where here C Σ=1 denotes the set of filters whose coefficients sum up to one. It is not necessary to assume that the weights are non-negative, but if not we must assume they are globally bounded (in absolute value) by some fixed constant. A generic point in Z has the form q = (q 1 i,j , . . . , q l i,j ), where each (q l i,j ) is a 2-dimensional vector, for indices i, j ∈ Z 2 ; relevant indices are in {1, ..., M } × {1, ..., N } plus the support of the convolution kernels (ξ l , η l ). Then, the norms we consider are

q Z = i,j,l |q l i,j | 2 , q * Z = max i,j,l |q l i,j | 2 .
Input FD The constraint in the primal problem F * q = Du reads (denoting q l i,j = (q l,1 i,j , q l,2 i,j ) T ):

(D 1 u) i+ 1 2 ,j = m,n,l ξ l m,n q l,1 i+m,j+n , (D 2 u) i,j+ 1 2 ,j = m,n,l
η l m,n q l,2 i+m,j+n .

Boundary conditions.

In this section we briefly discuss how the most popular boundary conditions can be realized in our framework.

Homogeneous Neumann conditions. Homogeneous Neumann conditions are treated as follows: we assume that p ∈ Y , that is, the values of p are defined exactly at the same points as the values of D 1 u, D 2 u. For computing the convolutions, we implicitly assume p = 0 when the indices fall out of range. The variable q ∈ Z, on the other hand, should be nonzero for a larger set of indices, namely i ∈ {-ν, . . . , M + ν}, j ∈ {-ν, . . . , N + ν}.

Homogeneous Dirichlet conditions. For Dirichlet conditions, it is roughly the opposite. Indeed, one can assume that u = 0 and D α u = 0, α = 1, 2, out of the domain, that is, u i,j = 0 if i ≤ 0 or j ≤ 0 or i ≥ M + 1 or j ≥ N + 1. In particular, the variable p multiplying a difference Du which is zero should be unconstrained, i.e., for instance, p 1 i+ 1 2 ,0 for any i. Then, one should consider only convolutions F * p which do not involve unconstrained values of p. This amounts to let q l,α i,j = 0, when the convolution (F l,α p α ) i,j involves values of p out of the support (cf. (2.8)), in particular this implies that (F l,α ) * q l,α = 0 at points where D α u has been forced to zero, in coherence with the constraint F * q = Du.

Heterogeneous Dirichlet conditions. In practice, for solving problems with non homogeneous boundary conditions, we will consider the Neumann discretization introduce above, and will solve our problems with an additional equality constraint on u in a strip of width the size of the convolution operators.

General form: consistency.

To simplify in this section we let N = M , and consider the domain Ω = (0, 1) 2 , divided in N × N pixels of size ε × ε, ε = 1/N . We add the subscript ε to all the notions introduced in the previous section, which now depend on N = 1/ε, which will be sent to ∞. We let

D ε = (1/ε)D : X ε → Y ε , D defined as in (2.1). We define (2.10) T V ε (u) = min ε 2 q Zε : F * ε q = D ε u = sup ε 2 p, D ε u Yε : F ε p * Z ≤ 1 when u ∈ L 1 (Ω) is of the form u = N i,j=1 u i,j χ (iε-ε,iε)×(jε-ε,jε)
, and +∞ else. In (2.10), D ε has the scaling of a discrete gradient, and the ε 2 factor before the sum p, D ε u Yε accounts for the size of the pixels. We only consider, for simplification, the case of Neumann homogeneous boundary conditions as described in the previous section. Then the following holds: Theorem 2.4. Assume the supports and the weights of the convolutions defining

F ε are uniformly bounded. Then T V ε Γ-converges to T V (u) := |Du|(Ω) if u ∈ BV (Ω) , +∞ else.
The convergence both holds in very weak topologies (such as distributional), and in stronger ones such as L p (Ω) for any p < +∞.

Proof. To prove the "lim inf", we assume u ε → u in the distributional sense and that lim inf

ε T V ε (u ε ) is finite. Consider p(x) = (p 1 (x), p 2 (x)) ∈ C ∞ c (Ω; R 2 ) a smooth vector field with compact support and |p(x)| ≤ 1 for every x ∈ Ω. Then Ω u div p dx = lim ε→0 Ω u ε div p dx = lim ε→0 ε 2 D ε u ε , p ε Yε
where the discretized p ε is given by the horizontal and vertical fluxes2

p ε,1 i+ 1 2 ,j = (1/ε) jε (j-1)ε p 1 (iε, y)dy , p ε,2 i,j+ 1 2 = (1/ε) iε (i-1)ε p 2 (x, jε)dx .
We used here that with this definition, which corresponds to the Raviart-Thomas projection of the vector field p, the average of div p on a pixel is precisely the constant div p ε .

Since p is smooth, there exists C depending only on ∇p ∞ and F such that |(F l p ε ) i,j | ≤ 1 + Cε for all i, j, l. Indeed, using the form (2.8), one sees that, denoting

x ε i,j := ((i + 1 2 )ε, (j + 1
2 )ε) the center point of the pixel (iε, (i + 1)ε) × (jε, (j + 1)ε),

(F l,1 p ε,1 ) i,j = ν m,n=-ν ξ l m,n p ε,1 i+ 1 2 -m,j-n = p 1 (x ε i,j ) + ν m,n=-ν ξ l m,n (p ε,1 i+ 1 2 -m,j-n -p 1 (x ε i,j ))
Using that p is smooth and in particular, Lipschitz-continuous, the last term can be bounded by a quantity Cε/ √ 2, for some constant C ∼ ν ∇p ∞ m,n |ξ l m,n | which, by assumption, is uniformly bounded. In the same way,

|(F l,2 p ε,2 ) i,j -p 2 (x ε i,j )| ≤ Cε/ √ 2. Since |p(x ε i,j )| 2 ≤ 1, our claim follows.
We deduce that

F p ε * Z ≤ 1 + Cε, so that p ε /(1 + Cε) is an admissible dual variable in (2.10) so that ε 2 D ε u ε , p ε Yε ≤ (1 + Cε)T V ε (u ε ).
Sending ε to zero and then taking the supremum with respect to all p we obtain that:

(2.11) T V (u) ≤ lim inf ε→0 T V ε (u ε ).
We remark that near the boundaries, p ε vanishes, so this construction remains valid in the framework of the Neumann boundary conditions described in Section 2.3.

Let us now prove the "lim sup". Let us consider u ∈ BV (Ω). We need to build a "recovery" sequence, that is a sequence (u ε ) of discrete functions of the form (2.12)

u ε (x) = N i,j=1 u ε i,j χ (iε-ε,iε)×(jε-ε,jε) (x), for ε = 1/N , N ≥ 1, such that u ε → u in L 1 (Ω) (or actually in some L p (Ω), 1 ≤ p < ∞, if u ∈ L p (Ω) 3 ) and: (2.13) lim sup ε→0 T V ε (u ε ) ≤ T V (u).
As is classical, a first observation is that one may assume that u is bounded, indeed the sequence of truncated functions defined by

u k (x) := max{-k, min{k, u(x)}} , k ≥ 1, satisfies u k → u in L 1 (Ω) and Ω |Du k | → Ω |Du|.
Then, from a recovery sequence for each u k , one easily builds using a diagonal argument a recovery sequence for u. For convenience, one can assume also that u ∈ C ∞ (Ω). Indeed, one can also approximate u with such functions, with again convergence of the total variations, see for instance [16, §1.15-1.16].

In a first step, as Ω = (0, 1) 2 is a square (a rectangle would work as well, for a general Lipschitz domain one should use a more complicated approach), it is clear that one can consider u as an even, 2-periodic functions (first let u(-x, y) = u(x, -y) = u(-x, -y) := u(x, y) for (x, y) ∈ (0, 1) 2 , then extend u 2-periodically in both directions). Observe that the Fourier transform of u, in that case, is real, even, and discrete. We can let, for

∈ Z 2 , û( ) := 1 4 [-1,1] 2 u(x)e -2iπ •x 2 dx = [0,1] 2 u(x) cos(π • x)dx so that u(x) = ∈Z 2 û( )e iπ •x . Then, we consider ρ ∈ C ∞ c (B(0, 1); [0, 1]
) with ρ ≡ 1 on B(0, 1/2) a smooth, symmetric cutoff function (in the Fourier domain), and let for ∈ Z 2 , ûδ ( ) := û( )|ρ(δ /2)| 2 . It turns out that this is equivalent, in the spatial domain, to convolving u with ρ δ * ρ δ , where ρ δ (x) := (1/δ 2 )ρ(x/δ) (is a smoothing kernel which is however not non-negative), and one has u δ → u in L 1 (Ω) with T V (u δ ) → T V (u). The latter holds for instance because since we assumed that u was initially smooth we had in particular [0,1] 2 |∇u| 2 dx < ∞, and passing to the Fourier domain it gets obvious that ∇u δ → ∇u strongly in L 2 (Ω), therefore also in L 1 (Ω).

Hence, again, without loss of generality, one needs just to build a recovery sequence for a function u (real, even, 2-periodic) with vanishing spectrum out of a disk: û( ) = 0 if | | ≥ 2/δ =: R, for some fixed value of δ.

Given N > R an integer and ε = 1/N , let, for (i, j) ∈ Z 2

(2.14)

u ε i,j = u(iε, jε) = N n,m=-N û(n, m)e iπ in+jm N .
This is of course real, even and periodic (u ε i+2N,j = u ε i,j+2N = u ε i,j for all i, j), and it is natural to view the restriction (u ε i,j ) 1≤i,j≤N and its piecewise constant extension (2.12) as a discrete approximation of the original image in the domain Ω = (0, 1) 2 .

In addition, being u (very) smooth, it is obvious that u ε (defined by (2.12)-(2.14)) converges in L 1 (Ω) (and in fact uniformly) to u, in the same way, as

ε → 0 (N → ∞), N i,j=1 u ε i,j -u ε i-1,j ε χ (iε-ε,iε)×(jε-ε,jε) (x) → ∂ 1 u(x)
(also uniformly) and in particular lim sup

ε→0 ε N i,j=1 (u ε i,j -u ε i-1,j ) 2 + (u ε i,j -u ε i,j-1 ) 2 ≤ Ω |∇u|dx = T V (u).
Let us now estimate T V ε (u ε ). We use the definition:

T V ε (u ε ) = ε 2 min { q Zε : F * ε q = D ε u ε }
with F ε a sum of L convolution operators as in (2.8). We need to find a q such that 4

(2.15)

ε 2 q Zε ε N i,j=1 (u ε i,j -u ε i-1,j ) 2 + (u ε i,j -u ε i,j-1 ) 2 .
In particular, it is enough to find q of the form (q 1 , 0, . . . , 0), so that without loss of generality, one can assume L = 1 and drop the superscript l. Therefore, we look for q = (q i,j ) (i, j ∈ Z, q i,j = (q 1 i,j , q 2 i,j ) T ∈ R 2 for all i, j) which satisfies

(2.16) (D 1 ε u ε ) i+ 1 2 ,j = ν m,n=-ν ξ m,n q 1 i+m,j+n , (D 2 ε u ε ) i,j+ 1 2 ,j = ν m,n=-ν η m,n q 2 i+m,j+n at least for 1 ≤ i ≤ N -1, 1 ≤ j ≤ N for the first equality and 1 ≤ i ≤ N , 1 ≤ j ≤ N -1 for
the second, and hoping to have a control on the norm q Z . We compute the discrete Fourier transform of the first equation: for (r, s) ∈ Z 2 , one has

1 2N N i,j=-N u ε i+1,j -u ε i,j ε e -iπ ir+js N = 1 2N N i,j=-N ν m,n=-ν ξ m,n q 1 i+m,j+n e -iπ ir+js N = 1 2N N i,j=-N ν m,n=-ν ξ m,n e iπ mr+ns N q 1 i+m,j+n e -iπ (i+m)r+(j+n)s N = q1 (r, s)   ν m,n=-ν ξ m,n e iπ mr+ns N   .
On the other hand, (2.17)

1 2N N i,j=-N u ε i+1,j -u ε i,j ε e -iπ ir+js N = e iπrε -1 ε û(r, s) vanishes for (r, s) > R. Using that m,n ξ m,n = 1 one sees that ν m,n=-ν ξ m,n e iπ mr+ns N = 1 + ν m,n=-ν ξ m,n (e iπ mr+ns N -1) → 1
as N → ∞ for any r, s with (r, s) ≤ R. In particular, for N large enough and (r,

s) ≤ R, one can just let q1 (r, s) = e iπrε -1 ε û(r, s) 1 + ν m,n=-ν ξ m,n (e iπ mr+ns N - 1)
(and one lets q1 (r, s) = 0 for (r, s) > R). This choices guarantees that the first equation in (2.16) is satisfied, and one does analogously for the second. One could fear that the resulting 4 to simplify the notation we do not denote the dependency of q on ε.

(q i,j ) is a complex vector field, yet since also its real part then should satisfy (2.16), and have therefore the same Fourier transform, this cannot be the case.

It remains to show that with this choice, q Zε satisfies (2.15). One has

q 1 i,j = N r,s=-N q1 (r, s)e iπ ir+js N = (r,s) ≤R e iπrε -1 ε û(r, s) 1 + ν m,n=-ν ξ m,n (e iπ mr+ns N -1)
e iπ ir+js N So that (using (2.17)),

(D 1 ε u ε ) i+ 1 2 ,j -q 1 i,j = (r,s) ≤R ν m,n=-ν ξ m,n (e iπ mr+ns N -1) 1 + ν m,n=-ν ξ m,n (e iπ mr+ns N -1) e iπrε -1 ε û(r, s)e iπ ir+js N and |(D 1 ε u ε ) i+ 1 2 ,j -q 1 i,j | ≤ (r,s) ≤R ν m,n=-ν |ξ m,n || sin(π mr+ns 2 ε)| 1 -ν m,n=-ν |ξ m,n || sin(π mr+ns 2 ε)| | sin πrε 2 | ε |û(r, s)| ≤ Cε
when ε > 0 is small enough. In the same way, one builds q 2 i,j which satisfies (2.16) and such that

|(D 2 ε u ε ) i,j+ 1 2 -q 2 i,j | ≤ Cε.

It follows that (2.15) holds with an error of order O(ε).

Observe to conclude that this construction is compatible with the Neumann boundary conditions described in Section 2.3. One can consider (for each N ) only the values of q i,j for i, j = -ν, . . . , N + ν, the values of p restricted to Y , and restrict F , F * to these values; in this setting the construction builds an admissible recovery sequence.

Compactness. Theorem 2.4 remains incomplete without the following associated compactness result, whose proof is relatively easy and standard. Proposition 2.5. Under the same assumptions as in Theorem 2.4, let (u ε ) ε>0 with:

sup ε>0 T V ε (u ε ) < +∞.
Assume moreover u ε remains bounded (on average, or in some L p (Ω)). Then there is a subsequence (u ε j ) and u ∈ L 1 (Ω) such that u ε j → u in L 1 (Ω).

Proof. The point is that by assumption q ε is uniformly bounded, as well as the discrete derivative D ε u ε = F * q ε . It easily follows that (u ε ) ε have uniformly bounded total variation, and the claim can be deduced.

Learning a better discretization.

In the previous section, we have shown consistency of the general form (2.3) of the discrete total variation, in the spirit of Γ-convergence, for any interpolation (averaging) kernel F . The aim of this section is to investigate numerical methods to learn the "best" interpolation operators F given a set of input images and corresponding ground truth solutions. To do so, we will consider several total variation minimization tasks, where explicit solutions are known (e.g. denoising of a disk, image inpainting) but we will also consider the case of natural image denoising. It will turn out that the learned interpolation operators highly depend on the specific imaging task. That is, an interpolation operator that is found to work well for image denoising will not necessarily work well for image inpainting and vice versa. The reason is that there is a complex interaction between the total variation term and the data fidelity term. For example, in image denoising, the total variation should better smooth nicely the solution (whereas the recovery of sharp edges is mostly driven by the data term), while for image inpainting we expect it to discriminate sharp edges from smoother ones to produce sharper discontinuities.

In what follows, we consider the general class of total variation minimization problem:

(3.1) min

Du=F * q λ q Z + G(u, g),
where G(u, g) is a convex data fidelity term depending on a given input image g, which defines the type of imaging application, e.g. denoising, inpainting, segmentation etc. The parameter λ > 0 is used to control the tradeoff between regularization and data fidelity, when needed.

We also assume that G has a proximal map prox τ G (•), τ > 0 which is simple to compute (in closed form). We denote by u * a minimizer, unique in case G is strongly convex. Since (3.1) is a difficult non-smooth optimization problem, especially in its primal form, we consider the equivalent saddle-point formulation

(3.2) min u,q max p Du -F * q, p + λ q Z + G(u, g),
which can be solved by the (block)-preconditioned primal-dual algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Pock | Diagonal preconditioning for first order primal-dual algorithms[END_REF][START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF], whose main iterations are given in Algorithm 1.

Algorithm 1: Preconditioned primal-dual algorithm for solving (3.2)

• Initialization: u 0 ∈ X, q 0 ∈ Z, p 0 ∈ Y .

• Step sizes: Choose the block-wise step sizes τ u , τ q , σ p such that diag(σ p )

1 2 (D, F * ) diag(τ u , τ q ) 1 2 ≤ 1,
and set θ = 1. • Iterations: For k = 0, . . . , K -1 let

(3.3)              p k+1 = p k + σ p (Du k -F * q k ) pk+1 = p k+1 + θ(p k+1 -p k ) u k+1 = prox τuG (u k -τ u D * pk+1 ) q k+1 = shrink τqλ (q k -τ q F pk+1 ) • Output: Approximate saddle point (u K , q K , p K )
The advantage of block-wise step sizes is that one can easily adapt to different scalings in the linear operators D and F * . The proximal map prox τuG is assumed to be commutable in closed form, see [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] for several examples which will also be used here. The function shrink τ denotes the classical 2,1 norm shrinkage function which is given for each of the pairs q l = (q l,1 , q l,2 ) in q = (q 1 , . . . , q L ) by ql = shrink τ (q l ) ⇐⇒ ql i,j = 1 -

1 max{1, |q l i,j |/τ } ql i,j , l = 1, ..., L
3.1. The learning problem. Clearly, we are here in the classical setting of supervised learning. That is we assume we have given a set of input images G = {g 1 , . . . , g S } and corresponding target images (explicit solutions) T = (t 1 , . . . , t S ). The learning problem is now to find interpolation kernels F such that solutions u * s of (3.2) for input images g s minimize a certain distance measure to the target images t s .

In order to measure the distance (or loss) between u * s and t s we consider a convex and continuously differentiable loss function (u, t). Here, we adopt here the classical quadratic loss function (u, t) = 1 2 u -t 2 , but any other suitable loss function could be considered as well. The overall loss function with respect to the entire training data set is then given by

L(F ) = 1 M N S S s=1 (u * s (F ), t s ),
where we have made explicit the dependence of the solutions u * s (F ) on the interpolation kernels F . We also include a pre-factor 1 M N S in order to normalize the loss function with respect to the total number of image pixels.

The learning problem is now given by the following bilevel optimization problem [START_REF] Kunisch | A bilevel optimization approach for parameter learning in variational models[END_REF] min

F L(F ) + R(F ), (3.4) u * s ∈ arg min u,q max p Du -F * q, p + λ q Z + G(u, g s ), s = 1, . . . , S
Intuitively, this bilevel optimization problem accounts for finding interpolation kernels F such that the solutions of the saddle-point problems (here we are only interested in the primal solution u * s ) minimize the loss function on the entire training set. The function R(F ) is an additional regularization functional that can be used to impose certain constraints on F . As a general requirement for valid interpolation kernels, we must enforce that the filter coefficients in each individual filter sum up to one, that is F l,1 , F l,2 ∈ C Σ=1 , l = 1, . . . , L, where C Σ=1 was defined in (2.9). Hence we set

R(F ) = δ (C Σ=1 ) L,2 (F ) = L l=1 δ C Σ=1 (F l,1 ) + δ C Σ=1 (F l,2 ), here δ C (x) = 0 if x ∈ C, and δ C (x) = ∞ else, is the "indicator function" of a set C.
During learning we will need to compute the orthogonal projection

F = proj (C Σ=1 ) L,2 ( F )
which because of the independence of the constraints can be done for each filter separately onto the set C Σ=1 . The projection of an arbitrary vector x = (x 1 , . . . , xn ) ∈ R n onto the set C Σ=1 is given by

x = proj C Σ=1 (x) ⇐⇒ xi = xi + 1 -n i=1 xi n .
One can also envisage additional constraints in order to make the learned discretization invariant to certain discrete symmetries. We will come back to such constraints when they are needed in the applications. Note that in our framework all constraints are affine or linear, so that, in fact, the optimization can be considered as unconstrained in a proper linear subspace and the reprojection step is equivalent to considering as a descent direction the gradient of the objective restricted to the subspace.

The learning problem causes several difficulties. First, it is highly non-convex since the higher-level problem (loss function) depends on the solution of a non-smooth saddle-point problem. Second, computing gradients of the higher-level problem with respect to F is challenging since it involves the sensitivity of the saddle-point of the lower-level problem with respect to the interpolation kernel F .

Computing derivatives.

The most straightforward approach to compute the derivatives would be to unravel a certain number of iterations of the primal-dual algorithm (Algorithm 1) and to make use of the machinery of automatic differentiation and backpropagation. However, for hard problems such as inpainting, where one would need to unravel several thousands of iterations one easily runs into problems because of memory limitations. In our setting and implementation of the learning, unraveling only 100 iterations takes already more than 5 GB of memory.

In order to allow learning for large problems or a large number of iterations, we rely on standard sensitivity analysis to derive in Appendix A a method to compute derivatives of a loss function with respect to the linear operator appearing in the saddle-point problem. Although the theory is valid for smooth and strongly convex objectives, we use in practice the same strategy to compute our gradients. In order to simplify the presentation in this section, we only show the computation of the gradient of the loss function with respect to one datum (g s , t s ) and we also drop the index s to simplify the notation. The gradient with respect to the entire data set can be easily computed by just summing up the gradients from the individual data examples.

Assume we have computed a saddle-point u * , q * , p * of (3.2) for a certain training example (g, t). Then, according to (A.4), the next step is to solve for the adjoint state U * , Q * , P * which is also the solution of the bi-quadratic saddle point problem

(3.5) min U,Q max P DU -F * Q, P + 1 2 ∇ 2 G(u * )U , U - 1 2 ∇ 2 q * Z Q, Q + ∇ (u * , t), U .
Instead of doing this sequentially, we propose a "Piggyback"-style algorithm [START_REF] Griewank | Piggyback differentiation and optimization[END_REF] that jointly computes u * , q * , p * and its adjoint state U * , Q * , P * . Algorithm 2 gives the main iterations of the Piggyback primal-dual algorithm. It is the adaption to our setting of the Algorithm described Section A.2, whose convergence is proved in a smoother setting.

Algorithm 2: Piggyback primal-dual algorithm for solving (3.2) and its adjoint.

• Initialization: u 0 , U 0 ∈ X, q 0 , Q 0 ∈ Z, p 0 , P 0 ∈ Y .

• Step sizes: Choose the block-wise step sizes τ u , τ q , σ p such that diag(σ p )

1 2 (D, F * ) diag(τ u , τ q ) 1 2 ≤ 1,
and set θ = 1. • Iterations: For each k = 0, . . . , K -1 let

(3.6)                        p k+1 = p k + σ p (Du k -F * q k ), P k+1 = P k + σ p (DU k -F * Q k ) pk+1 = p k+1 + θ(p k+1 -p k ), P k+1 = P k+1 + θ(P k+1 -P k ) ũk+1 = u k -τ u D * pk+1 , Ũ k+1 = U k -τ u (D * P k+1 + ∇ (u k , t)) u k+1 = prox τuG (ũ k+1 ), U k+1 = ∇ prox τuG (ũ k+1 ) • Ũ k+1 qk+1 = q k -τ q F pk+1 , Qk+1 = Q k -τ q F P k+1 q k+1 = shrink τqλ (q k+1 ), Q k+1 = ∇shrink τqλ (q k+1 ) • Qk+1
• Output: Approximate saddle point (u K , q K , p K ) and corresponding adjoint state

(U K , Q K , P K )
Finally, as shown in (A.6), the gradient of (u * (F ), t) with respect to F is approximated by the linear form:

(3.7) ∇ F (u * (F ), t), L ≈ -Q K , Lp K -q K , LP K ⇐⇒ ∇ F (u * (F ), t) ≈ -(Q K ⊗ p K + q K ⊗ P K ).
We refer to the Appendix and [START_REF] Chambolle | Convergence of a Piggyback-style method for the differentiation of solutions of standard saddle-point problems[END_REF] for an estimate on the error of this approximation. Then, the gradient of the loss function on the whole data set is computed as:

(3.8) ∇L(F ) = 1 M N S S s=1 ∇ F (u * s (F ), t s )
Remark 3.1. The expression (3.7) of the gradient in this form is a bit abstract, but can be easily computed by means of automatic differentiation techniques since expressions such as F p or F * q are implemented as convolutions with particular boundary conditions. For this, one first implements the linear form in the first line of (3.7) (which is nothing but a linear first order approximation of the bilevel problem (3.4)) and then computes the gradient of this linear form with respect to the filter coefficients of the linear operators F via automatic differentiation. In all our experiments we use the software package Pytorch to compute the gradient which automatically takes care of the specific boundary conditions in the convolutions.

In a similar fashion, we compute the Jacobian of the proximal maps times vector products ∇ prox τuG (ũ k+1 ) • Ũ k+1 and ∇shrink τqλ (q k+1 ) • Qk+1 appearing in Algorithm 2 by means of automatic differentiation. Note that the respective proximal maps might not be continuously differentiable. This can be accounted for by either smoothing the proximal maps or by selecting an arbitrary sub(super)-gradient which is the standard procedure in the deep learning community for non-smooth activation functions such as the widely used rectified linear unit ReLU(t) = max{0, t}. We found the latter option to work well in practice and hence did not consider smoothing.

Learning algorithm.

Having detailed the computation of the gradient of the loss function with respect to the interpolation kernel F , we now consider the actual learning algorithm. Here, we use an inertial proximal gradient method whose main iterations are presented in Algorithm 3.

Algorithm 3: Proximal gradient method for solving (3.4)

• Initialization:

F -1 = F 0 ∈ (C Σ=1 ) L,2 . • Step sizes: Choose α k > 0, β k ∈ [0, 1). • Iterations: For k = 0, . . . , K -1 let (3.9)    F k = F k + β k F k -F k-1 F k+1 = proj (C Σ=1 ) L,2 F k -α k ∇L( F k )
• Output: Learned interpolation kernels F K Depending on the particular learning problem we will make use of different settings of the step size α k and the inertial parameter β k .

Numerical results.

In this section we make use of the algorithms of the previous section and show how one can learn the "best" interpolation kernels for different imaging applications. We start with the most simple problem which is the recovery (inpainting) of straight discontinuities but will also consider the regularization (smoothing) of the characteristic functions of disks and natural image denoising.

Data generation. In the rest of this section we assume a domain Ω = [-1, 1] 2 which is discretized into N × N square pixels of size h × h, hence h = 2/N . According to (2.4), the values u i,j , with 1 ≤ i, j ≤ N of the discrete images correspond to pixel averages of the continuous functions u(x).

In order to generate the discrete input images and corresponding ground truth solutions, we precompute pseudo-continuous images by first point-sampling the continuous functions on a regular high-resolution grid of BN × BN points and then generate the discrete input images g s and corresponding target images t s , s = 1..., at the target resolution N × N via block averaging with blocks of size B × B. In all experiments we used a block size of B = 64.

We also include in all generated synthetic data random shifts of the grid of the size h/2 in order to simulate partial volume effects which is present in any real-world application. This sort of "data augmentation" is a well known technique in deep learning in order to make the learning more robust and generalize better to new data. In fact in our experiments it turned out that it is a simple way to remove inherent data symmetries which could lead to persistent local minima.

As it is standard in machine learning, we always make use of a training set to learn the filters and then evaluate the learned filters on test data which was not part of the training set. For the synthetic images, it will usually mean that the training set differs from the test set by the random shifts.

Filter sizes. In all synthetic experiments, where the structures to be reconstructed are basically lines or circles, we consider only a small neighborhood of 2×2 pixels, which corresponds to convolution kernels of size 2 × 3 for the horizontal component and 3 × 2 for the vertical component of the dual variable p. See also Figure 2.1 and note that the duals are "living" on the faces of the pixels. For natural image denoising, where locally the image structures are much more complex, we will also consider learning with larger neighborhoods.

Number of filters. We experimented with different numbers L of filter pairs. In order to investigate the influence of L to the reconstruction quality we performed experiments using L ∈ {2, 3, 4, 8} filter pairs. For natural images, we will also perform experiments with a larger number of filter pairs.

Filter symmetries. Depending on the number of filters it is natural to consider different symmetries which should ideally also reflect the symmetries of the data. For example images are known to by translation invariant, which is naturally captured by the convolution operation but images are also largely rotational invariant, which can be enforced by imposing a rotational invariance on the filters to be learned. Enforcing symmetries also reduces the effective parameters to be learned and one can expect that such learned filters generalize better to new data.

In the case of L = 2 or L = 3 we cannot impose a full rotational invariance and hence it is most natural to consider "only" a transpose symmetry as this symmetry group has in fact two elements. Interestingly, it will turn out that such a transpose symmetry is automatically learned (approximately) without enforcing it explicitly.

For L = 4 and L = 8 it is more natural to consider a rotational symmetry by 4 × π 2 , as this symmetry group has four elements.

In order to project the filters onto those symmetry groups, we first perform the projections onto the set C Σ=1 followed by a projection onto the respective symmetry group. As an example, for the transpose symmetry for L = 2, we used the following constraints:

F 1,1 = (F 2,2 ) , F 1,2 = (F 2,1 ) ,
onto which the projection is done by averaging the corresponding filter entries. Overall, this can be shown to give the orthogonal projection onto the intersection of both spaces.

Filter initialization. As the learning problem represents a nonconvex optimization problem, the initialization has an important impact on the outcome. One can consider different random initializations (uniform, Gaussian, ...) quite common in the deep learning literature, but for our problem it turns out to be better to initialize the filters weights by the linear interpolation weights obtained from interpolating the dual variables on certain regular or random positions. Figure 4.1 below shows the initial filters in case L = 8. Note that such initial filters correspond to the choice of the filter weights in the CD discretization.

Implementation and hardware. All algorithms are implemented in Python with the support of GPU processing and automatic differentiation from the Pytorch package. All code will be made available by the authors. The experiments were performed on an Nvidia Titan RTX eGPU equipped with 24GB memory connected to a Dell XPS Laptop with a Quad-Core i9 CPU, 32GB of main memory, and running Ubuntu 20.04. 4.1. Inpainting. We first consider the problem of inpainting straight discontinuities of various orientations from a given boundary datum. As seen on Fig. 2.2, most handcrafted discretizations lead to more or less suboptimal solutions and hence the aim of this section is to learn the best set of interpolation kernels which ensure a sharp discontinuity over a large range of orientations. The general problem (3.1) can be specialized to inpainting by choosing

G(u, g) = (i,j)∈I δ {g i,j } (u i,j ),
where I is a given set which contains the pixel indices corresponding to the given boundary datum and δ C (•) is the indicator function of the convex set C. We also denote the inpainting domain Ī = {1, ..., N } 2 \ I. The proximal map of τ G is given by

û = prox τ G (ū, g) ⇔ ûi,j = g i,j if (i, j) ∈ I ūi,j if (i, j) ∈ Ī.
Training and test data. The training data consists of S = 64 images, each of size N × N pixels, with N = 64 and uniformly sampled orientations θ s = 2sπ/S, s = 0, ..., S -1 of the discontinuities. Figure 4.2 shows a few examples of the training data consisting of the input images g s and the corresponding target images t s . The inpainting domains Ī are shown in gray. As already mentioned above, we also include a small random shift of size h/2 of the discontinuity from the center of the image in order to simulate partial volume effects. The test data was generated in exactly the same way.

Setting of the learning algorithm. First of all, solving the inpainting problem itself is already a very hard convex optimization problem (in fact similar problems are considered to derive worst-cases rates in convex optimization) and hence one typically needs thousands of iterations. Therefore we must expect that learning is even harder!

The gradients of the lower-level problem were computed by performing K = 2000 iterations of the Piggyback primal-dual algorithm (Algorithm 2). We also used a warm-starting strategy for both the variables and their adjoint states in order to improve the accuracy of the gradients. Note that due to memory limitations, it would be impossible to unravel such a high number of iterations of the primal-dual algorithm and compute the gradient via standard backpropagation.

In the learning algorithm (Algorithm 3) we use a diminishing step size rule of the form α k = α 0 √ k+1 for some α 0 > 0, as this choice is known to give an optimal convergence rate for the iterates in subgradient methods. The inertial parameter was set to β k = 0, since for this problem we have to expect errors in the gradient computation, and inertial methods are known to accumulate errors of the gradient.

We stopped the learning algorithm when the loss function did not show any significant improvement, which was usually the case after K = 1000 -2000 iterations.

Results. Table 4.1 plots the learned filters for different settings of the parameter L. First, all learned filters appear significantly different from the known handcrafted ones. Second, the filters L = 2 and L = 3 (where we did not enforce transpose symmetry) show an almost transpose symmetry, very similar to the transpose symmetric filters L = 2 (s) and L = 3 (s). This shows that it is absolutely reasonable to enforce symmetry constraints. Moreover, it reduces the effective number of learnable parameters, indeed in case of L = 4 (s) we effectively learned only 2 filters. Third, it is interesting that the learned filters also contain negative values. For example the filter coefficients of the first filter pair F 1 = (F 1,1 , F 1,2 ) for the setting L = 2 (s) are given by (rounded to three decimals):

F 1,1 = 0.002 0.205 0.334 -0.021 0.161 0.319 , F 1,2 =    -0.029 -0.033 -0.095 1.301 -0.007 -0.137   
We evaluated the learned filters on the test data, which consists of the same covering of orientations but with different partial volume effects. We used a number of K = 105 iterations of the primal-dual algorithm (Algorithm 1) to compute the solutions. Such a high number of iterations is absolutely necessary to obtain high-accuracy solutions, as the convergence of the primal-dual algorithm is very slow (in fact O(1/K)) for the inpainting problem.

L = 2 L = 2 (s) L = 3 L = 3 (s) L = 4 (s) L = 8 (s)
Table 4.1 summarizes and compares the (half) mean squared error (MSE = L(F )) 5 of stateof-the-art handcrafted discretizations and the learned discretizations. For a better readability, we multiply all MSE values by a factor of 10 5 .

First of all, one can see that the MSE is very similar on both the training data and the test data which indicates that overfitting is minor, particularly in cases where we explicitly enforced a certain symmetry constraint. As already discussed, the best handcrafted filters is given by CD, which outperforms the other handcrafted variants by a large margin. One can also see that having a larger number of filters leads to consistently better results, but it is more the symmetry of the filters which is important than the effective number of parameters. ) of the single test images over the range of orientations. Note that a PSNR larger than 50 can already be considered as almost perfect. From the plot one can see that the learned filters are consistently better compared to the handcrafted filters and are much more rotationally invariant. Note that some handcrafted filters are very good in recovering grid-aligned orientations but at the cost of poor results for non-grid-aligned orientations. proximal map, which is needed in the primal dual algorithm is given by û = prox τ G (ū, g) ⇐⇒ ûi,j = ūi,j + τ g i,j 1 + τ .

Data FD RT CD L = 2 L = 2 (s) L = 3 L = 3 (s) L = 4 (s) L = 8 (
It is a well known fact that for the disk denoising problem, the (unique) solution of the ROF model, is given by a disk of the same radius but with a reduced intensity value 1 -2λ rs , depending on the radius r s and the setting of the regularization parameter λ.

Training and test data. The training data comprises S = 64 input images g s of the characteristic functions of disks with radius r s , where we sample r s regularly in the range [0.25, 0.75]. The corresponding target images t s are obtained by reducing the intensity value of the disk to 1 -2λ rs . We used a setting of λ that ensures the existence of a solution for all values of r s . Additionally we incorporated a small random shift of the center of the disk of size h/2 to increase the variability of the data. Figure 4.5 exemplifies a few instances of the training data consisting of the input images and corresponding ground truth solutions. The test data was generated in exactly the same way but with different random shifts.

Setting of the learning algorithm. Due the the strong convexity of the ROF model (in the variable u), the computation of a high-accuracy solution is much easier than for the inpainting problem. Also the learning appears significantly easier in practice. For computing the derivatives, we used K = 200 iterations of the Piggyback algorithm (Algorithm 3), again with a warm-starting strategy for both the variables and their adjoint states. In the learning algorithm we used a constant step size α k = α 0 with α 0 > 0. The inertial parameter β k was set to β k = 1 √ 2 , as this setting is known to lead to a convergent algorithm for optimization problems with 1 α 0 -Lipschitz-continuous gradients [START_REF] Mukkamala | Convex-Concave Backtracking for Inertial Bregman Proximal Gradient Algorithms in Nonconvex Optimization[END_REF]. The learning was stopped when no significant improvement of the objective function was observed, in practice after approximately K = 1000 iterations.

Results. Table 4.3 summarizes the learned filters for the disk denoising problem. At first glance, the filters obtained for disk denoising look very different from the filters learned for inpainting (c.f. Table 4.1). Indeed the coefficients of the first filter pair for the setting L = 2 (s) obtained are given by (rounded to 3 decimals) Moreover, it is interesting that the filters obtained for the setting L = 4 (s) are very similar (up to a translation) to the handcrafted filters of RT. This is also confirmed by the MSE values in Table 4.4 where RT appears to be a very well performing handcrafted filter. From the table one can also see that learning consistently improves on the entire data set (training set and testing set). The best learned variant is given by the setting L = 8 (s). Figure 4.6 additionally plots for the testing data the PSNR values for the different handcrafted and learned filters over the radius of the disk. One can observe that the learned filters consistently perform better compared to the handcrafted filters over the whole range of radii. The simple (de facto standard) discretization FD gives worst results, both in terms of the sharpness of the boundary of the disk and the value of the disk itself. Indeed, one weakness of all handcrafted filters seem to be to identify correctly the value of the disk which results in a large error. This is exactly, where learning improves and finds a discretization which gives a better value of the disk. See Figure 4.6 for a color coding based visualization of this issue. In Figure 4.8, we show the dependence of the MSE on the discretization width h. For this we set the parameter λ such that the intensity value of the regularized disk is equal to 1/2, and we performed the disk denoising experiment for different values of h = 2/N with N ∈ {64, 128, 256, 512, 1024}. The MSE corresponds, in the continuous setting, to the squared L 2 norm of the solution with respect to the expected solution. The plot shows that RT and both learned discretizations L = 4 (s) and L = 8 (s) satisfy an O(h) error bound, although this is rigorously proved, for this experiment, only in the RT case [START_REF] Caillaud | Error estimates for finite differences approximations of the total variation[END_REF].

F 1,1 = -0.

Natural image denoising.

We now propose to learn the "best" filters for denoising natural images. Of course we have to admit that we do not expect the total variation to be competitive against state-of-the art image denoising methods, it is just interesting to see which kind of filters are learned for natural image denoising and how they compare to the filters learned on the synthetic images. We will consider the same setting as in the previous experiments, but will also present results of a larger experiment where we learn 40 filter pairs with an effective neighborhood of 6 × 6 pixels. Despite the fact we are restricted to the setting of consistent discretizations of the total variation, at the lowest scale this introduces a flavor of non-locality which is known to be the key to the success of state-of-the art approaches. We also present results in this setting enforcing an additional 90 • rotational invariance, so that the effective number of filter pairs is reduced to 10. For natural image denoising, we still consider the ROF model, now with Neumann boundary conditions. Since we are dealing with natural images, we do not know (even if we know the noise level) the optimal regularization parameter λ > 0. To overcome this issue, we fix in our optimization a value of λ but replace the sum constraint (2.9) by the looser constraints F ∈ C Σ=µ where:

(4.1) C Σ=µ = {F : ∃µ ∈ R, m,n ξ l m,n = m,n η l m,n = µ, l = 1, . . . L}.
In words, the sum of the coefficients of all the filters of F should be equal to the same value.

Exploiting the 1-homogeneity of the norm, the "optimal" value of the regularization parameter in the original problem (3.1) is then given by λ|µ|, and the normalized filter coefficients, which satisfy (2.9), are simply obtained by dividing all the coefficients of F by µ.

During learning, we project the coefficients onto the constraint (4.1). This of course has a simple closed form. Given m vectors xi = (x i 1 , . . . , xi n ) ∈ R n , i = 1, . . . , m, the projection onto (4.1) reduces to the problem:

min µ,(x i ) m i=1 1 2 m i=1 x i -xi 2 , s.t. n j=1 x i j = µ, i = 1, . . . m.
Assuming that the optimal solution is attained for some μ ∈ R, then the entries of the Training and test data. For the training data, we used S = 64 target patches t s , s = 1, . . . S, each of size 64×64 pixels, which are randomly cropped from images of the well-known BSDS500 database [START_REF] Arbelaez | Contour detection and hierarchical image segmentation[END_REF]. The input patches g s are generated by adding i.i.d. zero-mean Gaussian noise with a standard deviation of σ = 0.05. The test data was generated in exactly the same way, but of course contain a different set of randomly cropped patches. For the handcrafted filters the optimal regularization parameter was determined to be λ = 0.025 and kept constant for both the training and test set. For the learned filters we used λ = 0.02 but the effective regularization parameter was adapted during learning, by the additional freedom to adjust the sum of the filter coefficients.

Setting of the learning algorithm. We used exactly the same setting of the learning algorithm as in the case of the disk denoising experiment.

Results. In Table 4.5 we plot the learned filters for the setting L = 8 (s) using a neighborhood of 2 × 2 pixels, L = 40 and L = 40 (s) using a neighborhood of 6 × 6 pixels. Comparing the setting L = 8 (s) obtained for natural images with the respective filters obtained for disk denoising, we see that the filters for natural images are significantly different. The main reason is that natural photographs do not to contain sharp discontinuities because of the band limitation of optical devices. The larger filters show quite nicely that the learned interpolation kernels correspond to oriented edge-and line-like structures and hence better adapt to the structures of natural images. As mentioned before, another interpretation is that the learned filters correspond to a convolutional dictionary for the local image gradient Du. which, with our assumptions, has for any given linear operator K : X → Y * Y a unique saddle point (x, ŷ). Then, we consider the convex and continuously differentiable loss function L(K) = (x(K), ŷ(K)), where we have made explicit the dependence of the saddle point (x, ŷ) with respect to the linear operator K. Our aim is to find a derivative of L with respect to the linear operator K.

Following a classical approach in sensitivity analysis, we start from the optimality condition of the saddle-point problem (A. Thanks to the strong convexity of g, f * (and the continuity of the Hessians) we can pass to the limit as s → 0 and see that ξ s , η s go to ξ, η satisfying K * η + L * ŷ + D 2 g(x)ξ = 0, -Kξ -Lx + D 2 f * (ŷ)η = 0, where we have exchanged the order of the equations and changed the sign of the second line. The unique solution is given by

ξ η = D 2 g(x) K * -K D 2 f * (ŷ) -1 -L * ŷ Lx
We can now compute the directional derivative L (K; L) = ∇L(K), L , as follows:

(A.3) L (K; L) = ∇ (x, ŷ) T ξ η = ∇ (x, ŷ) T D 2 g(x) K * -K D 2 f * (ŷ) -1 -L * ŷ Lx .
Next, we introduce adjoint variables X, Y , such that

(-X T , Y T ) = ∇ (x, ŷ) T D 2 g(x) K * -K D 2 f * (ŷ) -1 ⇔ ∇ (x, ŷ) = D 2 g(x) -K * K D 2 f * (ŷ)
-X Y , ∇ x (x, ŷ)): however also in the latter we have to substitute (x, ŷ) with known values, which at this point would be (x k , y k ). This would transform (A.7) into:

(A.9)

       X k+1 = ∇ prox τ g (x k -τ K * y k ) • (X k -τ (K * Y k + ∇ x (x k , y k ))) Xk+1 = X k+1 + θ(X k+1 -X k ) Y k+1 = ∇ prox σf * (y k + σK xk+1 ) • (Y k + σ(K Xk+1 + ∇ y (x k , y k ))).
where x k , y k , xk+1 are computed using (A.8). Iterations (A.8) and (A.9) translate to Alg. 2 in the particular setting of this paper.

A.3. Convergence analysis. We denote γ the modulus of convexity of g and δ the modulus of f * . Thanks to [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]Thm. 3], choosing µ = 2 √ γδ/ K , 1/(1 + µ) ≤ θ ≤ 1, τ = µ/(2γ) and σ = µ/(2δ) we have:

γ x k -x 2 + (1 -ω)δ y k -ŷ 2 ≤ Cω k
for ω = (1 + θ)/(2 + µ) < 1 (precisely, ω = θ in case θ = 1/(1 + µ), ω = 1/(1 + µ/2) in case θ = 1). Hence, one should expect X k , Y k to converge to the solution of (A.4) if g, f * and are smooth enough, following the analysis for inexact primal-dual algorithms (see for instance [START_REF] Rasch | Inexact first-order primal-dual algorithms[END_REF]). We can sketch a convergence analysis, following [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF][START_REF] Rasch | Inexact first-order primal-dual algorithms[END_REF], as follows. We write (A.7) as

D 2 g(x k+1 ) • X k+1 = X k -X k+1 τ -(K * Y k + ∇ x (x k , y k ))) D 2 f * (y k+1 ) • Y k+1 = Y k -Y k+1 σ + (K Xk+1 + ∇ y (x k , y k ))).
Subtracting (A.4), we get that where C depends on the Lipschitz constants mentioned above. We observe that we would still obtain a (slower) linear rate if we assumed only Hölder continuity of the same functions. We deduce that

D 2 g(x k+1 ) • (X k+1 -X) = X k -X k+1
γ X k+1 -X 2 ≤ 1 τ X k -X k+1 , X k+1 -X -Y k -Ŷ , K(X k+1 -X) + e k+1 X , X k+1 -X , that is 1 + 2τ γ 2τ X k+1 -X 2 + 1 2τ X k+1 -X k 2 ≤ 1 2τ X k -X 2
-Y k -Ŷ , K(X k+1 -X) + e k+1 X , X k+1 -X and in the same way

1 + 2σδ 2σ Y k+1 -Ŷ 2 + 1 2σ Y k+1 -Y k 2 ≤ 1 2σ Y k -Ŷ 2 + Y k+1 -Ŷ , K( Xk+1 -X) + e k+1 Y , Y k+1 -Ŷ .
As in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] we introduce ∆ k := 1 2τ X k -X 2 + 1 2σ Y k-1 -Ŷ 2 . We sum the previous inequalities to find that

(A.11) (1 + µ)∆ k+1 + 1 2τ X k+1 -X k 2 + 1 2σ Y k -Y k-1 2 ≤ ∆ k -Y k -Ŷ , K(X k+1 -X k ) + θ (Y k-1 -Ŷ k , K(X k -X k-1 ) + θ Y k -Y k-1 , K(X k -X k-1 ) + e k+1 X , X k+1 -X + e k Y , Y k -Ŷ .
The simplest way to get rid of the two error terms (which might not be optimal) is to bound these by e k+1

)∆ k+1 + 1 2τ X k+1 -X k 2 + 1 2σ Y k -Y k-1 2 ≤ ∆ k -Y k -Ŷ , K(X k+1 -X k ) + θ (Y k-1 -Ŷ k , K(X k -X k-1 ) + θ Y k -Y k-1 , K(X k -X k-1 ) + 1 2γ e k+1 X 2 + 1 2δ e k Y 2 .
Letting now ω := (1 + θ)/(2 + µ/2) ∈ (ω, 1), using (A.10) and adapting the proofs of [6, Thm. 3] [24, Thm. 7], we obtain that (X k , Y k ) → ( X, Ŷ ) at rate O(ω k/2 ) (which is strictly slower than the rate for (x k , y k )). As mentioned already, if ∇ , D 2 g, D 2 f * are only Hölder continuous, we may still adapt the above proof to still obtain a (slower) linear convergence rate.

2 i,j+ 1 2 Figure 2 . 1 . 1 i+ 1 2 2 i,j+ 1 2

 22211121 Figure 2.1. General setting of the discretization. Filled circles correspond to discrete pixel values ui,j localized at the centers of the pixels. Unfilled circles correspond to the discrete dual variables, with the values p 1 i+ 1 2 ,j localized on vertical edges and p 2 i,j+ 1 2

Figure 2 . 2 .

 22 Figure 2.2. Reconstruction quality of various "handcrafted" discretization schemes for the problem of inpainting a straight discontinuity. The numbers above the images refer to the reconstruction errors, measured in PSNR.

Figure 4 . 1 .

 41 Figure 4.1. Initialization for L = 8 filter pairs using the weights of linear interpolation, where white corresponds to a weight of 1.0, gray corresponds to a weight of 0.5 and dark gray corresponds to a weight of 0.25.

Figure 4 . 2 .

 42 Figure 4.2. Samples of the training data consisting of the input images gs and corresponding target images ts for learning to inpaint a straight discontinuity. The disk-like inpainting masks Ī are indicated by the gray areas.

Figure 4 .

 4 Figure 4.3 plots the peak signal to noise ratio (PSNR = -10 log 10 2MSE) of the single test images over the range of orientations. Note that a PSNR larger than 50 can already be considered as almost perfect. From the plot one can see that the learned filters are consistently better compared to the handcrafted filters and are much more rotationally invariant. Note

  n = 3 (s) L = 3, n = 3 (s) L = 4, n = 3 (s) L = 8, n = 3 (s)

Figure 4 . 3 .

 43 Figure 4.3. Reconstruction quality (in PSNR) plotted over the angles of the discontinuities. One can see that the learned filters offer a significantly improved isotropy.

Finally, Figure 4

 4 .4 provides false-color error plots of one training example with respect to the corresponding target image. One can see that FD and RT largely fail but CD gives very good results. The learned filters, in particular the setting L = 4 (s), give still better results (although hard to see with the eye).

Figure 4 . 4 .

 44 Figure 4.4. Color coded differences to the ground truth image ts for s = 56, θs ≈ 7π 4 . CD gives the best result among the hand crafted filters but overall, the learned filters using L = 4 (s) gives the best results.

4. 2 .

 2 Disk denoising. In our second experiment, we consider regularizing (denoising) a white disk over a black background using the ROF model and 0-Dirichlet boundary conditions. The ROF model, is obtained from (3.1) by setting G(u, g) = 1 2 u -g 2 2 . The corresponding (a) Input images gs (b) Target images ts

Figure 4 . 5 .

 45 Figure 4.5. Samples of the training data for learning to denoise disks of different sizes.

  Data FD RT CD L = 2 L = 2 (s) L = 3 L = 3 (s) L = 4 (s) L = 8 (s) Train 22.28 1.

Figure 4 . 6 .

 46 Figure 4.6. Reconstruction errors measured in PSNR for different handcrafted and learned filters for the disk denoising problem.

Figure 4 . 7 .

 47 Figure 4.7. Color coded differences to the ground truth image ts with s = 64. RT gives the best result among the handcrafted filters but the value of the disk is well reconstructed only for the learned filter k = 8, s.

= 1 ,Figure 4 . 8 .

 148 Figure 4.8. Plot of the MSE error with respect to the inverse discretization width h -1 of the image domain, for the disk experiment. One sees that the handcrafted discretization RT as well as the two learned discretizations L = 4 (s) and L = 8 (s) satisfy a O(h) error bound.

Figure 4 . 9 .

 49 Figure 4.9. Samples of the training data for natural image denoising. The input images are given by adding i.i.d zero-mean Gaussian noise with a standard deviation of σ = 0.05.

L = 8 2 L 6 L

 826 (s), 2 × = 40 (s), 6 × = 40, 6 × 6 Table 4.6 presents quantitative results and comparisons with hand-crafted discretizations. One can see that the learned discretizations lead to significantly better values than any of the handcrafted discretizations. Interestingly, while the setting L = 40 yields the smallest training MSE, the setting L = 40 (s), enforcing the rotational invariance, gives a slightly better test MSE and seems to generalize better to new data. In Figure 4.10 we finally provide a per-instance comparison between the PSNR values of the learned discretization L = 40 (s) and the handcrafted discretization CD. One can see

  Input g CD, PSNR=28.98 L = 40 (s), PSNR=29.77 (c) Example from the test set

Figure 4 . 10 .

 410 Figure 4.10. Comparison (in PSNR) between the best handcrafted filters (CD) and the learned filters L = 40 (s) for both the training set (a) and the test set (b). Subfigure (c) exemplifies a denoised image patch from the test set where the learned discretization of the total variation gives a significant improvement.



  K x -∇f * (ŷ) = 0 K * ŷ + ∇g(x) = 0 and assume K is perturbed by a small variation sL, |s| 1, in a given direction L. Denoting xs = x + sξ s and ŷs = ŷ + sη s the solution for the linear operator K + sL it follows thatK x + s(Kξ s + Lx s ) -∇f * (ŷ) -s 0 D 2 f * (ŷ + tη s )dt η s = 0, K * ŷ + s(K * η s + L * ŷs ) + ∇g(x) + s 0 D 2 g(x + tξ s )dt ξ s = 0.Dividing by s and making use of the optimality condition (A.2) yields  Kξ s + Lx s -1 s s 0 D 2 f * (ŷ + tη s )dt η s = 0 K * η s + L * ŷs + 1 s s 0 D 2 g(x + tξ s )dt ξ s = 0.

τ- 2 ≤

 2 K * (Y k -Ŷ ) + e k+1 X D 2 f * (y k+1 ) • (Y k+1 -Ŷ ) = Y k -Y k+1 σ + K( Xk+1 -X) + e k+1 Ywhere the error terms are given bye k+1 X = -(∇ x (x k , y k ) -∇ x (x, ŷ) + (D 2 g(x k+1 ) -D 2 g(x)) • X), e k+1 Y = ∇ y (x k , y k ) -∇ y (x, ŷ) -(D 2 f * (y k+1 ) -D 2 f * (ŷ)) • Ŷ .In particular, if we assume additionally that ∇ and D 2 g, D 2 f * are Lipschitz, we have the estimate, for k ≥ 1, Cω k

X 2 /

 2 (2γ) + e k Y 2 /(2δ) + (µ/2)∆ k+1 , thanks to Young's inequality. It yields (A.12) (1 + µ 2

Table 4 . 1

 41 Learned filters for the inpainting problem where symmetric filters are marked by (s).

Table 4 .2 10

 4 5 × the mean squared error (MSE) of handcrafted and learned filters evaluated on both the training and test data.

Table 4 . 3

 43 Learned filters for the disk denoising problem. Symmetric filters are indicated by (s).

Table 4 .

 4 4 10 5 times the mean squared error (MSE) of handcrafted and learned filters for the disk denoising problem.

Table 4 . 5

 45 Learned filters for natural image denoising. Symmetric filters are indicated by (s).

Table 4 .

 4 6 10 4 × the mean squared error (MSE) of handcrafted and learned filters for natural image denoising.

	Data FD RT CD L = 8 (s) L = 40 (s) L = 40
	Train 5.05 5.33 4.87	4.58	4.31	4.22
	Test 4.72 5.05 4.51	4.28	4.10	4.13
	that as expected, a clear majority of image instances are better denoised by the learned
	discretization. We also exemplify provide one typical denoising results, where the learned
	discretization yields an improvement of almost 0.8 dB.

Table 4 . 7

 47 Confusion matrix for the crossover experiments. For a better comparison, all values are here on the scale of 10 5 × the mean squared error (MSE). For comparison the last column depicts the results for the handcrafted filters of CD.

In practice, we will consider rectangular filters, but they can be always implemented by means of larger square filters and setting the respective filter weights to zero.

The notation might seem a bit inconsistent, as the centers of the pixels correspond in fact to the points x ε i,j := ((i + 1 2 )ε, (j + 1 2 )ε) while we choose to denote, as is usual, by (i, j) this location.

To simplify we consider in the sequel only the case p = 1.

Please note the additional factor of 1 2 in L(F ) which we keep in order to make the MSE compatible to our loss function.

* This work was funded by ERC grant HOMOVIS, No. 640156.

We conclude this section by carrying out crossover experiments, where we compare the filters learned on one specific task to the other tasks. The aim of such crossover experiments is to investigate the extent to which the learned discretization schemes can be generalized to other tasks. In order to make a relatively fair comparison, we only compare the filter of the setting L = 8 (s), for which we have obtained good results on all the tasks. In Table 4.7, we provide a "confusion matrix" of the results and we put the respective results of the handcrafted scheme CD as a reference. One can see that the filters obtained for line inpainting seem to generalize best to the other tasks. Surprisingly, the filters obtained for natural image denoising task seem to generalize better than the filters learned for the synthetic disk denoising task. It seems that for disk denoising the role of the quadratic data fidelity term is too strong such that almost "any" consistent discretization of the total variation may give good results. At this point one could of course try to learn optimal filters for several tasks simultaneously, and one can expect that this gives even better generalization capabilities, but we leave this experiments for future work.

Appendix A. Derivatives of saddle-points.

In this appendix, we give a partial justification of the "Piggyback"-style Algorithm 2 the derivative of our loss function. We consider to simplify a smoother objective. This is obtained by assuming that both • Z and G(•, g) in (3.2) are replaced with smooth and strongly convex functions: in such situation our sensitivity analysis shows that Algorithm 2 actually converges to states and adjoint states from which one can recover a gradient of the objective. In general, we can just hope that it produces a limit of gradients which would be computed by smoothing the appropriate functions and then letting the smoothing go to zero, but this remains an open question. On the other hand, the good quality of the gradients we compute and the results of the optimization give a hint that it seems to be the case in our application. In a forthcoming work [START_REF] Chambolle | Convergence of a Piggyback-style method for the differentiation of solutions of standard saddle-point problems[END_REF], we will further analyse this algorithm in a slightly more general setting.

A.1. Theoretical grounds. We consider here a generic problem of type (3.2), yet in a smoother case. Following the notation in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF], we are given two convex, lower semicontinuous functions g, f such that g, f * are smooth (at least C 2 , we will require a bit more in Sec. A.3) and strongly convex (defined on some finite dimensional Euclidean spaces X, Y ). We consider then the standard saddle point problem which we can write as:

Equation (A.4) are the optimality conditions of the bi-quadratic adjoint saddle-point problem:

Denoting by ( X, Ŷ ) the unique solution of this adjoint saddle-point problem, the directional derivative (A.3) is given by

for any L. Using L (K; L) = ∇L(K), L we find that the gradient is given by

In practice, our variable will be only a subset of the operator K, in which case it is convenient to assemble the linear form (A.5) and use automatic differentiation to extract its gradient with respect to the relevant parameters.

A.2. Algorithms. Algorithm 3 in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] (or Alg. 5 in [START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF]) can be implemented to solve both (A.2) and (A.4). For the latter it yields the following iterations: one chooses τ, σ > 0, θ ≤ 1 (particular choices will be discussed later on), choose (X 0 , Y 0 ) and let for each k:

Here (x, ŷ) have to be computed first by solving (A.2), which can be done using the same algorithm (and actually, with the same parameters τ, σ, θ): we pick (x 0 , y 0 ) (possibly as the result of a former optimization) and let for each k:

Now, a natural algorithm consists in replacing in (A.7) the values of (x, ŷ), which can be estimated only up to some precision, with the iterate (x k , y k ) or (x k+1 , y k+1 ) of the iteration (A.8). Denoting prox τ g = (I + τ ∇g) -1 , etc, we observe that for any x, ∇ prox τ g (x) = (I + τ D 2 g(prox τ g (x))) -1 .

Hence ∇ prox τ g (x k -τ K * y k ) = (I + D 2 g(x k+1 )) -1 , etc. Computationally, the idea would be to evaluate such a gradient using automatic differentiation, against the vector X k -τ (K * Y k +