Bruno Lecouat Inria

Jean Ponce
email: jean.ponce@inria.fr

Julien Mairal
email: julien.mairal@inria.fr

A Flexible Framework for Designing Trainable Priors with Adaptive Smoothing and Game Encoding

We introduce a general framework for designing and training neural network layers whose forward passes can be interpreted as solving non-smooth convex optimization problems, and whose architectures are derived from an optimization algorithm. We focus on convex games, solved by local agents represented by the nodes of a graph and interacting through regularization functions. This approach is appealing for solving imaging problems, as it allows the use of classical image priors within deep models that are trainable end to end. The priors used in this presentation include variants of total variation, Laplacian regularization, bilateral filtering, sparse coding on learned dictionaries, and non-local self similarities. Our models are fully interpretable as well as parameter and data efficient. Our experiments demonstrate their effectiveness on a large diversity of tasks ranging from image denoising and compressed sensing for fMRI to dense stereo matching.

Introduction

Despite the undeniable successes of deep learning in domains as varied as image processing [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF] and recognition [START_REF] He | Deep residual learning for image recognition[END_REF], natural language processing [START_REF] Collobert | A unified architecture for natural language processing: Deep neural networks with multitask learning[END_REF], speech [START_REF] Oord | Wavenet: A generative model for raw audio[END_REF] or bioinformatics [START_REF] Alipanahi | Predicting the sequence specificities of dna-and rna-binding proteins by deep learning[END_REF], feedforward neural networks are often maligned as being "black boxes" that, except perhaps for their top classification or regression layers, are difficult or even impossible to interpret. In imaging applications, for example, the elementary operations typically consist of convolutions and pointwise nonlinearities, with many parameters adjusted by backpropagation, and no obvious functional interpretation.

In this paper, we consider instead network architectures explicitly derived from an optimization algorithm, and thus interpretable from a functional point of view. The first instance of this approach we are aware of is LISTA [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF], which provides a fast approximation of sparse coding. Yet, we are not content to design an architecture that provides a fast approximation to a given optimization problem, but we also want to learn a data representation pertinent for the corresponding task. This yields an unusual machine learning paradigm, where one learns the parameters of a parametric objective function used to represent data, while designing an optimization algorithm to minimize it efficiently.

Even though interpretability is not always necessary to achieve good prediction, this point of view, sometimes called algorithm unrolling [START_REF] Fletcher | Plug-in estimation in high-dimensional linear inverse problems: A rigorous analysis[END_REF][START_REF] Monga | Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing[END_REF], has proven successful for solving inverse imaging problems, providing effective and parameter-efficient models. This approach allows the use of domain-specific priors within trainable deep models, leading to a large number of applications such as compressive imaging [START_REF] Sun | Deep admm-net for compressive sensing mri[END_REF][START_REF] Zhang | Ista-net: Interpretable optimization-inspired deep network for image compressive sensing[END_REF], demosaicking [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF], denoising [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF][START_REF] Scetbon | Deep k-svd denoising[END_REF][START_REF] Simon | Rethinking the CSC model for natural images[END_REF], and super-resolution [START_REF] Wang | Deep networks for image super-resolution with sparse prior[END_REF] .

However, existing approaches are often limited to simple image priors such as sparsity induced by the 1 -norm [START_REF] Simon | Rethinking the CSC model for natural images[END_REF], or differentiable regularization functions [START_REF] Lefkimmiatis | Universal denoising networks: a novel cnn architecture for image denoising[END_REF], and a general algorithmic framework for combining complex, possibly non-smooth, regularization functions is still missing. Our paper addresses this issue and is able to leverage a large class of image priors such as total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], the 1 -norm, structured sparse coding [START_REF] Mairal | Non-local sparse models for image restoration[END_REF], or Laplacian regularization, where local optimization problems interact with each others. The interaction can be local among direct neighbors on an image grid, or non-local, capturing for instance similarities between spatially distant image patches [START_REF] Buades | A non-local algorithm for image denoising[END_REF][START_REF] Dabov | BM3D image denoising with shape-adaptive principal component analysis[END_REF].

In this context, we adopt a more general and flexible point of view than the standard convex optimization paradigm, and consider formulations to represent data based on non-cooperative games [START_REF] Nikaidô | Note on non-cooperative convex games[END_REF] potentially involving non-smooth terms, which are tackled by using the Moreau-Yosida regularization technique [START_REF] Hiriart-Urruty | Convex analysis and minimization algorithms I: Fundamentals[END_REF][START_REF] Yu | Better approximation and faster algorithm using the proximal average[END_REF]. Unrolling the resulting optimization algorithm results in a network architecture that can be trained end-to-end and capture any combination of the domain-specific priors mentioned above. This approach includes and improves upon specific trainable sparse coding models based on the 1 -norm for example [START_REF] Simon | Rethinking the CSC model for natural images[END_REF][START_REF] Wang | Deep networks for image super-resolution with sparse prior[END_REF]. More importantly perhaps, it can be used to construct several interesting new image priors: In particular, we show that a trainable variant of total variation and its non-local variant based on self similarities is competitive with the state of the art in imaging tasks, despite using up to 50 times fewer parameters, with corresponding gains in speed. We demonstrate the effectivness and the flexibility of our approach on several imaging tasks, namely denoising, compressed fMRI reconstruction, and stereo matching.

Summary of our contributions. First, we provide a new framework for building trainable variants of a large class of domain-specific image priors . Second, we show that several of these priors match or even outperform existing techniques that use a much larger number of parameters and training data. Finally, we present a set of practical tricks to make optimization-driven layers easy to train.

Background and Related Work

Classical image priors. Inverse imaging problems are often solved by minimizing a data fitting term with respect to model parameters, regularized with a penalty that encourages solutions with a particular structure. In image processing, the community long focused on designing handcrafted priors such as sparse coding on learned dictionaries [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF][START_REF] Mairal | Sparse modeling for image and vision processing[END_REF], diffusion operators [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], total variation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], and non-local self similarities [START_REF] Buades | A non-local algorithm for image denoising[END_REF], which is a key ingredient of successful restoration algorithms such as BM3D [START_REF] Dabov | BM3D image denoising with shape-adaptive principal component analysis[END_REF]. However these methods are now often outperformed by deep learning models [START_REF] Liu | Non-local recurrent network for image restoration[END_REF][START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF][START_REF] Zhang | Ffdnet: Toward a fast and flexible solution for cnn-based image denoising[END_REF], which leverage pairs of corrupted/clean training images in a supervised fashion.

Bilevel optimization.

A simple method for mixing data representation learning with optimization is to use a bi-level formulation [START_REF] Mairal | Task-driven dictionary learning[END_REF]. For instance, assuming that one is given pairs (x i , y i ) i=1...n of corrupted/clean signals with x i and y i in R m , one may consider the following bi-level objective

min θ∈Θ,W∈R m×p 1 n n i=1 L(y i , Wz θ (x i )) where z θ (x i ) ∈ arg min z∈R p h θ (x i , z), (1)
where θ is a set of model parameters, Wz θ (x i ) is a prediction which is compared to y i through a loss function L : R m × R m → R + , and the data representation z θ (x i ) in R p is obtained by minimizing some function h θ . Note that for simplicity, we have considered here a multivariate regression problem, where given a signal x in R m , we want to predict another signal y in R m , but this formulation also applies to classification problems. It was first introduced for sparse coding in [START_REF] Mairal | Task-driven dictionary learning[END_REF][START_REF] Yang | Coupled dictionary training for image super-resolution[END_REF] and it has recently been extended to the case when z θ (x i ) is replaced by an approximate minimizer of h θ .

Unrolled algorithms. A common approach to solving (1) consists in choosing an iterative method for minimizing h θ and then define z θ (x i ) as the output of the optimization method after K iterations.

The sequence of operations performed by the optimization method can be seen as a computational graph and ∇ θ z θ can be computed by automatic differentiation. This often yields neural-network-like computational graphs, which we call optimization-driven layers. Such architectures have found multiple applications such as training of conditional random fields [START_REF] Zheng | Conditional random fields as recurrent neural networks[END_REF], stabilization of generative adversarial networks [START_REF] Metz | Unrolled generative adversarial networks[END_REF], structured prediction [START_REF] Belanger | End-to-end learning for structured prediction energy networks[END_REF], or hyper-parameters tuning [START_REF] Maclaurin | Gradient-based hyperparameter optimization through reversible learning[END_REF]. For image restoration, various optimization problems have been explored including for example sparse coding [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF][START_REF] Simon | Rethinking the CSC model for natural images[END_REF][START_REF] Zhang | Ista-net: Interpretable optimization-inspired deep network for image compressive sensing[END_REF], non linear diffusion [START_REF] Chen | Trainable nonlinear reaction diffusion: A flexible framework for fast and effective image restoration[END_REF] and differential operator regularization [START_REF] Lefkimmiatis | Universal denoising networks: a novel cnn architecture for image denoising[END_REF]. Many inference algorithms have been investigated including proximal gradient descent [START_REF] Lefkimmiatis | Universal denoising networks: a novel cnn architecture for image denoising[END_REF][START_REF] Simon | Rethinking the CSC model for natural images[END_REF], ADMM [START_REF] Fletcher | Plug-in estimation in high-dimensional linear inverse problems: A rigorous analysis[END_REF], half quadratic spliting [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF], or augmented Lagrangian [START_REF] Romano | The little engine that could: Regularization by denoising (red)[END_REF].

3 A General Framework for Learning Optimization-Driven Layers

Proposed Approach

We adopt a more general point of view than [START_REF] Alipanahi | Predicting the sequence specificities of dna-and rna-binding proteins by deep learning[END_REF], where we assume that input signals admit a local "patch" structure (e.g., rectangular image regions) and the data representation encodes individual patches. Assuming that there are m patches in x, we denote by Z (x) = [z 1 (x), . . . , z m (x)] in R p×m the representation of x and by z j (x) the representation of patch j (we omit the dependency on the model parameters θ for simplicity). In imaging applications and as in previous models [START_REF] Simon | Rethinking the CSC model for natural images[END_REF], Z (x) can be seen as a feature map akin to that of a convolutional neural network with p channels.

P i x z i Z x
Figure 1: Our models encode locally an input feature vector. The local optimal solutions z j interact trough the regularization function ψ j θ (Z).

Encoding with non-cooperative convex games. Concretely, given a signal x, we denote by P j x the patch of x centered at position j, where P j is a linear patch extraction operator, and we define the optimal encoding Z (x) of x as a Nash equilibrium of the set of problems min

zj ∈Z h θ (P j x, z j ) + ψ j θ (Z) for j = 1, . . . , m, (2)
where h θ is a a convex reconstruction objective for each patch, parametrized by θ, ψ j θ is a convex regularization function encoding interactions between the variable z j and the remaining ones z l for l = j, and Z is a convex subset of R p . When Z is compact, the problem is a specific instance of a non-cooperative convex game [START_REF] Nikaidô | Note on non-cooperative convex games[END_REF], which is known to admit at least one Nash equilibrium-that is, a solution such that one of the objectives in ( 2) is optimal with respect to its variable z j when the other variables z l for l = j are fixed. The conditions under which an optimization algorithm is guaranteed to return such an equilibrium point are well studied, see Section 3.3, and in many situations the compactness of Z is not required, as also observed in our experiments where we choose Z = R p . For instance, in several practical cases, (2) can be solved by minimizing the sum of m convex terms, a setting called a potential game, which boils down to a classical convex optimization problem.

Application of our Framework to Inverse Imaging Problems

In this section, we show how to leverage our optimization-driven layers for imaging. For the sake of clarity we choose to narrow down the scope of this presentation to imaging, even though our method is not limited to this single application: different modalities including for example genomic/graph data could benefit from our methodology.

Examples of models h θ . We consider two cases in the rest of this presentation:

• Pixel reconstruction: h θ (P j x, z j ) = (x j -z j ) 2 , where x j is the pixel j of x and z j is a scalar, corresponding to patches of size q = 1 × 1 and p = 1. • Patch encoding on a dictionary: h θ (P j x, z j ) = P j x j -Dz j 2 , where D in R q×p is a dictionary, q is the patch size, and p is the number of dictionary elements. This is a classical model where patch j is approximated by a linear, often sparse, combination of dictionary elements [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF].

Only the second choice involves model parameters D (represented by θ). These two loss functions are common in image processing [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], but other losses may be used for other modalities.

Linear reconstruction with a dictionary. Assuming that y and x have the same size m for simplicity, predicting y from a feature map Z (x) is typically achieved by using a learned dictionary matrix W in R q×p where q is the patch size. Then, Wz j (x)3 can be interpreted as a reconstruction of the j-patch of y . Since the patches overlap, we obtain q estimators for every pixel, which can be combined by averaging (neglecting border effects below for simplicity), yielding the prediction where P j is the linear operator that places a patch of size q at position j in a signal of dimension m. Patch averaging is a classical operation in patch-based image restoration algorithms, see [START_REF] Elad | Image denoising via sparse and redundant representations over learned dictionaries[END_REF], which can be interpreted in terms of transposed convolution 4 and admits fast implementations on GPUs.

ŷ(x, θ, W) = 1 q m j=1 P j Wz j (x), (3)
Learning problem. For image restoration, given training pairs of corrupted/clean images {x i , y i } i=1,...,n , we consider the regression problem min θ∈Θ,W∈R q×p y iŷ(x i , θ, W) 2 where ŷ(x i ) is defined in [START_REF] Belanger | End-to-end learning for structured prediction energy networks[END_REF].

Examples of regularization functions ψ j

θ . Our framework allows the use of several regularization functions, which are presented in the table 1 below. We assume that the patches are nodes in a graph, and denote by N j the set of neighbors of the patch j. For natural images, the graph may be a two-dimensional grid with edge weights a j,k that depend on the relative position of the patches j and k, which we denote by a j-k , but it may also be a non-local graph based on some similarity function as in [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF][START_REF] Liu | Non-local recurrent network for image restoration[END_REF]. Concretely, we can consider:

• the distance d j,k NL = diag(κ)(P j x -P k x) 2 between patches j and k of the image x, where κ in R q is a set of parameters to learn, and q is the patch size, and we define normalized weights a j,k NL = e -d j,k NL / l∈Nj e -d j,k NL . • or a distance inspired from the bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF]. In that case we define the distance

d j,k BL = xi-xj 2 2σ 2 d + i-j 2 2σ 2 r
between pixels on a local window N j centered around pixel j, and we define again normalized weights a j,k BL = e -d j,k BL / l∈Nj e -d j,k BL .

Table 1: A non-exhaustive list of regularization functions ψ θ covered by our framework.

ψ j θ (Z) Model parameters Laplacian k∈Nj a j-k z j -z k 2 weights in R |N | Non-local Laplacian k∈Nj a j,k NL z j -z k 2 κ in R q Bilateral filter (BF) k∈Nj a j,k BL z j -z k 2 σ d ∈ R and σ r ∈ R Total variation (TV) k∈Nj a j-k z j -z k 1 weights in R |N | Non-local total variation (NLTV) k∈Nj a j,k NL z j -z k 1 κ in R q Bilateral TV (BLTV) k∈Nj a j,k BL z j -z k 1 σ d ∈ R and σ r ∈ R Weighted 1 -norm (sparse coding) p l=1 λ l |z j [l]| λ in R p Non-local group regularization p l=1 λ l k∈Nj a j,k z k [l] 2 λ in R p and κ in R q Variance reduction Wz j -P j ŷ 2 with ŷ from (3) W from (3)
s Novelty of the proposed formulation and relation to previous work.

• Total variation: to the best of our knowledge, the basic anisotropic TV penalty [START_REF] Chambolle | An introduction to total variation for image analysis[END_REF] does not seem to appear in the literature on unrolled algorithms with end-to-end training. Note also that our TV variant allows learning non-symmetric weights a j,k = a k,j , leading to a non-cooperative game that goes beyond the classical convex optimization framework typically used with the TV penalty. • Non-local TV: the non-local TV penalty presented above is based on a classical formulation [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF],

but can be incorporated within a trainable deep network with non-symmetric weights. • Bilateral filtering: the bilateral filter and its TV variant implemented in this paper are based on classical formulations [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF][START_REF] Farsiu | Advances and challenges in super-resolution[END_REF]. But they have not, to the best of our knowledge, been implemented as trainable priors.

• Sparse coding and variance reduction: the weighted 1 -norm combined with the patch encoding loss h θ yields a sparse coding formulation (SC) that has been well studied within optimizationdriven layers [START_REF] Scetbon | Deep k-svd denoising[END_REF][START_REF] Simon | Rethinking the CSC model for natural images[END_REF]. Yet, the codes z j in the SC setup are obtained by solving independent optimization problems, which has motivated by Simon and Elad [START_REF] Simon | Rethinking the CSC model for natural images[END_REF] to propose instead a Convolutional Sparse Coding model (CSC), where the full image is approximated by a linear combination of small dictionary elements. Unfortunately, as noted in [START_REF] Simon | Rethinking the CSC model for natural images[END_REF], CSC leads to ill-conditioned optimization problems, making a hybrid approach between SC and CSC more effective. Our paper proposes an alternative solution combining the weighted 1 -norm regularization with a variance reduction penalty, which forces the codes z j to reach a consensus when reconstructing the image ŷ.

Our experiments show that this approach outperforms [START_REF] Simon | Rethinking the CSC model for natural images[END_REF] for image denoising. • Non-local group regularization: This regularization function corresponds to a soft variant of the Group Lasso penalty [START_REF] Turlach | Simultaneous variable selection[END_REF], which encourages similar patches to share similar sparsity patterns (set of non-zero elements of the codes z j ). It was originally used in [START_REF] Mairal | Non-local sparse models for image restoration[END_REF] and was recently revisited within optimization-driven layers with an heuristic algorithm [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF]. Our paper provides a better justified algorithmic framework as well as the ability to combine this penalty with other ones. Z t+1 ← Z t -η t H θ (Z t , X); 5: end for 6: Approximate clean images by linear reconstruction ŷ = 1 q m j=1 P j Wz j (X, θ); 7: Compute the 2 reconstruction loss yŷ(x, θ, W) 2 2 on the minibatch; 8: Compute an estimate of the gradients wrt. (θ, W) with auto-diff; 9: Update trainable parameters (θ, W) with Adam;

Differentiability and End-to-end Training

In this section we adress end-to-end training of the optimization-driven layers. Given pairs of training data {x i , y i } i=1,...,n , we consider the learning problem

min θ∈Θ 1 n n i=1 L (y i , g θ (z θ (x i ))) , (4)
where g θ is a differentiable function. We consider the approximation where the codes z j (x) are obtained as the K-th step of an optimization algorithm for solving the problem (2). To obtain these codes, we leverage (i) iterative gradient and extra-gradient methods, which are classical for solving game problems [START_REF] Bertsekas | Parallel and distributed computation: numerical methods[END_REF][START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF], and (ii) a smoothing technique for dealing with the regularization functions ψ j θ above when they are non-smooth. Refer to Algorithm 1 for an overview of the training procedure. We start with the first point when dealing with smooth objectives. Unrolled optimization for convex games. Consider a set of m objective functions of the form

min zj ∈R p h j (Z) with Z = [z 1 , . . . , z m ], (5)
where the functions h j are convex and differentiable and may depend on other parameters than z j . Our objective is to find a zero of the simultaneous gradient

H(Z) = [∇ z1 h 1 (Z), • • • , ∇ zm h m (Z)], (6)
which corresponds to a Nash equilibrium of the game (5). In the rest of this presentation, we consider both the general setting and the simpler case of so-called potential games, for which the equilibrium can be found as the optimum of a single convex objective. This is the case for several of our regularizers, for example the TV penalty with symmetric weights. More details are provided in Appendix A on the nature of the non-cooperative games corresponding to our penalties. Two standard methods studied in the variational inequality litterature [START_REF] Bertsekas | Parallel and distributed computation: numerical methods[END_REF][START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF][START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Mertikopoulos | Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile[END_REF] are the gradient and the extra-gradient [START_REF] Korpelevich | The extragradient method for finding saddle points and other problems[END_REF] methods. The iterates of the basic gradient method are given by where η t > 0 is a step-size. These iterates are known to converge under a condition called strong monotonicity of the operator H, which is related to the concept of strong convexity in optimization, see [START_REF] Bertsekas | Parallel and distributed computation: numerical methods[END_REF]. Because this condition is relatively stringent, the extra-gradient method is often prefered [START_REF] Korpelevich | The extragradient method for finding saddle points and other problems[END_REF], as it is known to converge under weaker conditions, see [START_REF] Juditsky | Solving variational inequalities with stochastic mirror-prox algorithm[END_REF][START_REF] Mertikopoulos | Optimistic mirror descent in saddle-point problems: Going the extra (gradient) mile[END_REF]. The intuition of the method is to compute a look ahead step in order to compute more stable directions of descent:

Z t+1 = Z t -η t H(Z t ), (7)
Extrapolation step Z t+1/2 = Z t -η t H(Z t ) Update step Z t+1 = Z t -η t H(Z t+1/2 ). (8)
In this paper, our strategy is to unroll iterates of one of these two algorithms, and then to use auto differentiation for learning the model parameters θ. Furthermore, parameters that control the optimization process (e.g., step size η t ) can also be learnt with this approach. It should be noted that optimization-driven layers have never been used before in the context of non-cooperative games, to the best of our knowledge, and therefore an empirical study is needed to choose between the strategies [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF] or [START_REF] Chang | Pyramid stereo matching network[END_REF]. In our experiments, extra-gradient descent has always performed at least as well, and sometimes significantly better, than plain gradient descent for comparable computational budgets. See for example Table 2 for a smoothed variant of the TV penalty.

Moreau-Yosida smoothing. The non-smooth regularization functions we consider can be written as a sum of simple terms. Omitting the dependency on θ for simplicity, we may indeed write

ψ j (Z) = r k=1 φ k (L k,j (Z)
) for some r ≥ 1, where L k,j is a linear mapping and φ k is either the 1 -or 2 -norm. For instance, φ k is the 1 -norm with L k,j (Z) = a k,j (z j -z k ) in R p for the TV penalty, and

L k,j (Z) = [ √ a 1,j z 1 (k), . . . , √ a 1,j z q (k)]
in R q with φ k being the 2 -norm for the non-local group regularization. Handling such non-smooth convex terms may be achieved by leveraging the so-called Moreau-Yosida regularization [START_REF] Lemaréchal | Practical aspects of the moreau-yosida regularization: Theoretical preliminaries[END_REF][START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace Hilbertien[END_REF][START_REF] Yosida | Functional analysis[END_REF]]

Φ k (u) = min v φ k (v) + α 2 v -u 2 ,
which defines an optimization problem whose solution is called the proximal operator Prox

φ k /α [u].
As shown in [START_REF] Lemaréchal | Practical aspects of the moreau-yosida regularization: Theoretical preliminaries[END_REF], Φ k is always differentiable and ∇Φ k (u) = α(u -Prox φ k /α [u]), which can be computed in closed form when φ k = 1 or 2 . The positive parameter α controls the trade-off between smoothness (the gradient of Φ k is α-Lipschitz) and the quality of approximation. It is thus natural to define a smoothed approximation Ψ j of ψ j as Ψ j (Z) = r k=1 Φ k (L k,j (Z)). Note that when the proximal operator of ψ j can be computed efficiently, as is the case for the 1 -norm, gradient descent algorithms can typically be adapted to handle the non-smooth penalty without extra computational cost [START_REF] Mairal | Sparse modeling for image and vision processing[END_REF], and there is no need for Moreau-Yosida smoothing. However, the proximal operator of the TV penalty and the non-local group regularization do not admit fast implementations. For the first one, computing the proximal operator requires solving a network flow problem [START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF], whereas the second one is essentially easy to solve when the weights a j,k form non-overlapping groups of variables, leading to a penalty called group Lasso [START_REF] Turlach | Simultaneous variable selection[END_REF]. We are now ready to present our unrolled algorithm as we have previously discussed gradient-based algorithms for solving convex smooth games and a smoothing technique for handling non-smooth terms. Generally, at iteration t, the gradient algorithm (7) performs the following simultaneous updates for all problems j

u (t) k,j ← Prox φ k /α k,t [L k,j (Z (t) )] for k = 1, . . . , r z (t+1) j ← z (t) j -η t ∇ zj h θ P j x, z (t) j + r k=1 α k,t L * k,j L k,j (Z (t) ) -u (t) k,j j
, where L * k,j is the adjoint of the linear mapping L k,j . The computation of the gradients can be implemented with simple operations allowing auto-differentiation in deep learning frameworks. Interestingly, the smoothing parameter α can be made iteration-dependent, and learned along with other model parameters such that the amount of smoothing is chosen automatically.

Tricks of the Trade for Unrolled Optimization

Our strategy is to unroll iterates of our algorithms, and then compute ∇ θ z θ by automatic differentiation. We present here a set of practical rules, some old and some new, facilitating training when h θ is a patch encoding function on a dictionary D.

Initialization. To help the algorithm converge, we choose an initial stepsize η t ≤ 1 L , where L is the Lipschitz constant of ∇ z h θ , which is the classical step-size used by ISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. To do so, inspired by [START_REF] Simon | Rethinking the CSC model for natural images[END_REF] we normalize the initial dictionary by its largest singular value and take η 0 = 1. Note that we can go one step further and normalize the dictionary throughout the training phase. This is in fact equivalent to the spectral normalization that has received some attention recently, notably for generative adversarial networks [START_REF] Miyato | Spectral normalization for generative adversarial networks[END_REF].

Untied parameters. In our framework, ∇ zj h θ (P j x, z j ) = D (Dz j -P j x). It has been suggested in previous work [START_REF] Gregor | Learning fast approximations of sparse coding[END_REF][START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF][START_REF] Simon | Rethinking the CSC model for natural images[END_REF] to introduce an additional parameter C of the same size as D, and consider instead the parametrization C (Dz j -P j x), C acting as a learned preconditioner. Even though the theoretical effect of this modification is not fully understood, it has been observed to accelerate convergence and boost performance for denoising tasks [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF]. In our experiments, we will indicate in which cases we use this heuristic.

Backtracking. A simple way for handling the potential instability of the unrolled algorithm is to use a backtracking scheme which automatically decreases the stepsize when the training loss diverges. This heuristic was used for instance in [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF]. More details are provided in Appendix B.

Barzilai-Borwein method for choosing the stepsize. A different, perhaps more principled, approach to improved stability consists in adaptively choosing an adaptive stepsize η t . The literature on convex optimization proposes a set of effective rules, known as Barzilai-Borwein (BB) step size rules [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF]. Even though these rules were not designed for convex games, they appear to be very effective in practice in the context of our optimization-driven layers. Concretely, they lead to step sizes η t,j = D Ds j 2 / Ds j 2 with s j = z

(t) j -z (t-1) j
for problem j at iteration t. In our experiments, we observed that spectral normalization, backtracking, and Barzilai-Borwein step size were all effective to stabilize training. We have noticed that the spectral normalization impacts negatively the reconstruction accuracy, while the BB method tend to improve it by using larger stepsizes, at the expense of a larger computational cost. This is illustrated in Table 3 for a smoothed variant of sparse coding (we indicate with a crossmark when the algorithm diverges). In addition, we observe that the untied models brings a small boost in reconstruction accuracy.

Experiments

We consider three different tasks, illustrated with various combinations of regularization functions in order to demonstrate the wide applicability of our approach and its flexibility. A software package and additional details are provided in the supplementary material for reproducibility purposes.

Image denoising. For image denoising experiments, we use the standard setting of [START_REF] Zhang | Beyond a gaussian denoiser: Residual learning of deep cnn for image denoising[END_REF] with BSD400 [START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF] as a training set and on BSD68 as a test set. We optimize the parameters of our models using Adam [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] and also use the backtracking strategy described in Section 3.4 that automatically decreases the learning rate by a factor 0.5 when the training loss diverges. For the non-local models, we follow [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF] and update the similarity matrices three times during the inference step. We use the parametrization with the C matrix for our patch-based experiments. We also combine our variance regularization with [START_REF] Lecouat | Fully trainable and interpretable non-local sparse models for image restoration[END_REF]. Additional training details and hyperparameters choices can be found in Appendix B. We report performance in terms of averaged PSNR in Table 4, and more detailed tables with additional results are available in Appendix C for pixel-level models, and for the patch-based models involving a dictionary D. Our models based on non-local sparse approximations perform better than the competing deep learning models with the exception of [START_REF] Liu | Non-local recurrent network for image restoration[END_REF] for σ ≥ 15 with much fewer parameters. In addition, we also observed that our assymetric TV models are almost on par with BM3D while being significantly faster (see Appendix C for more details) with only a very small amount of parameters.

where the degradation matrix is A = PF, P is a diagonal binary sampling matrix for a given sub-sampling pattern, F is the discrete Fourier transform such that the observed corrupted signal x is in the Fourier domain, and Ψ is a regularization function. This problem highlights the ability of our framework to handle both localized and non localized constraints. In our paper, we implemented two models revisiting some well studied priors for compressed sensing in an end-to-end fashion:

• Pixel reconstruction with total variation: we aim at solving the optimization for each node min yi∈R Ayx 2 2 + TV i (y). In the past, total variation has been widely used for MRI [START_REF] Lustig | Sparse mri: The application of compressed sensing for rapid mr imaging[END_REF], often in combination with sparse regularization in the wavelet domain.

• Patch encoding on a dictionary with sparse coding: we solve a collection of optimization problems of the form min zi∈R n Aŷ(z) -x 2 2 + λ z i 1 , with y = 1 n j R j Dz j the average of the overlapping patches. Some previous methods have explored dictionary-based reconstruction [START_REF] Ravishankar | Mr image reconstruction from highly undersampled k-space data by dictionary learning[END_REF], but they were not investigated from a task-driven manner with end-to-end training.

In our experiments, we use the same setting as [START_REF] Sun | Deep admm-net for compressive sensing mri[END_REF] for fair comparison: we train and test our models on the brain MRI dataset studied in that paper. Our models are trained separately for each sampling rate. We used the pseudo radial sampling for the matrix P similarly to the other methods. The reconstruction accuracy are reported in term of PSNR over the test set in Table 5. Our trainable model relying on a trainable TV prior performs surprisingly well given the conceptual simplicity of the prior. Also importantly, it runs significantly faster than all competing methods with a very small number of parameters. Furthermore, our trainable sparse coding method for fMRI gives strong performance and exceeds the state of the art for sampling rates larger than 30%. Note that architecture choices (patch and dictionary size) of our models are the same as for the denoising task, and we did not try to optimize them for the considered task, thus demonstrating the robustness of our approach. Dense Stereo Matching. Our approach can be used to provide a generic regularization module that can easily be integrated into various neural architectures. We showcase its versatility by using it for deep stereo matching [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. Given aligned image pairs, the goal is to compute disparity d for each pixel. Traditionally stereo matching is formulated as minimization of an energy function E data (d) + λE smooth (d) where the data term, E data measures how well d agrees with the input image pairs, E smooth enforces consistency among neighboring pixels' disparities: TV is a commonly chosen regularizer. Recent deep learning methods tackle the problem as a supervised regression to estimate continuous disparity map given pairs of stereo views and ground truth disparity maps [START_REF] Chang | Pyramid stereo matching network[END_REF]. We propose to combine our smoothing TV block with a state-of-the-art deep learning model [START_REF] Chang | Pyramid stereo matching network[END_REF]. In practice, we combine our block with a pretrained model on the SceneFlow [START_REF] Mayer | A large dataset to train convolutional networks for disparity, optical flow, and scene flow estimation[END_REF] dataset, and fine-tune the pretrained model on the kitti2015 [START_REF] Geiger | Are we ready for autonomous driving? the kitti vision benchmark suite[END_REF] train set, following the training procedure described in [START_REF] Chang | Pyramid stereo matching network[END_REF]. We used the original implementation of [START_REF] Chang | Pyramid stereo matching network[END_REF] available online and did not change any hyperparameters. We report in Table 7 the performance on the validation set in term of 3 pixels error which counts predicted pixel as correct if the disparity deviates from the ground truth from 3 pixels or less. We ran the experiment 10 times for each model (with and without the TV regularization). We observed that our TV block introduces very few additional parameters and consistently boosts performances.

Training with few examples. We conducted denoising experiments with less training data and report corresponding results in Table 6. We use the code released by the authors for training DnCNN with less data. Very interestingly the gap between our best model and CNN-based models increases when decreasing the size of the training set. We believe that this is an appealing feature, particularly relevant for applications in medical imaging or microscopy where the amount of training data can be very limited.

Discussion

We have presented a general framework based on non-cooperative games to train end-to-end imaging priors. Our experiments demonstrate the flexibility and the effectiveness of our approach on diverse tasks ranging from image denoising to fMRI reconstruction and dense stereo matching. Beyond image processing, we believe that the issue of interpretability is important. We consider models with a clear mathematical description of the decision function they produce. As a by-product, our models are also more parameter efficient than classical deep learning models. We believe that these are important steps to build systems that should not be seen as black boxes anymore, that produce explanable decisions, and that do not require training a system for days on a huge corpus of annotated data. These are important questions, which we are planning to address explicitly in the future.

Appendix

This supplementary material is organized as follows: In Section A, we discuss additional priors that were not presented in the main paper, but which are in principle compatible with our framework, and we provide more details about potential games. In Section B, we provide implementation details that are useful to reproduce the results of our paper (note that the code is also provided). In Section C, we present additional quantitative results and additional results regarding inference speed of our models that were not included in the main paper for space limitation reasons. Finally, in Section D, we present additional qualitative results (which require zooming on a computer screen).

A Discussion on Models and Priors

A. [START_REF] Alipanahi | Predicting the sequence specificities of dna-and rna-binding proteins by deep learning[END_REF]

Additional Priors

Our framework makes it possible to handle models of the form:

h j (Z) = h θ (P j x j , z j ) + λ r k=1 φ k (L k,j (Z)), (10)
where φ k is a simple convex function that admits a proximal operator in closed form, and L k,j is a linear operator. In the main paper, several regularization functions have been considered, including the total variation, variance reduction, or non-local group regularization penalties. Here, we would like to mention a few additional ones, which are in principle compatible with our framework, but which we did not investigate experimentally. In particular, two of them may be of particular interest, and may be the topic of future work:

• the regularization λ H z j 1 , where H is a matrix, may correspond to several settings. The matrix H may be for instance a wavelet basis, or may by learned, corresponding then to the penalty used in the analysis dictionary learning model from the paper "The cosparse analysis model and algorithms" of Nam et al., 2013. • the regularization λφ(H z j ) where φ is a smooth function is closely related to the model introduced in [START_REF] Lefkimmiatis | Universal denoising networks: a novel cnn architecture for image denoising[END_REF], and to the Field of experts model of Roth and Black from the 2005 paper "Fields of Experts: A Framework for Learning Image Priors", even though the functions used in these other works are not convex.

A.2 Potential Games

A potential game is a non-cooperative convex game whose Nash equilibria correspond to the solutions of a convex optimization problem. We will now consider problems of the form [START_REF] Collobert | A unified architecture for natural language processing: Deep neural networks with multitask learning[END_REF], and show that all penalties that admit some symmetry are in fact potential games. Assuming the functions φ k to be smooth for simplicity, optimality conditions for the convex problems [START_REF] Collobert | A unified architecture for natural language processing: Deep neural networks with multitask learning[END_REF] are, for all j = 1, . . . , m:

∇ zj h θ (P j x j , z j ) + λ r k=1 ∇ zj φk,j (Z) = 0, with φk,j (Z) = φ k (L k,j (Z)). (11)
Let us now assume the following symmetry condition such that if problem l involves a variable z j through a function φk,l (Z), then problem j also involves the same term. Based on this assumption, we may define the potential function

V (Z) := m j=1 h θ (P j x j , z j ) + λ 2 r k=1 φk,j (Z) .
The partial derivative of this potential function with respect to z j is then

∇ zj h θ (P j x j , z j ) + λ 2 m l=1 r k=1 ∇ zj φk,l (Z) = ∇ zj h θ (P j x j , z j ) + λ 2 m l=1 k∈N j,l ∇ zj φk,l (Z),
where N j,l is the set of functions φk,l involving variable z j . The previous gradient can then be simplified into

∇ zj h θ (P j x j , z j ) + λ 2 r j=1 ∇ zj φk,l (Z) + λ 2 l =j k∈N j,l ∇ zj φk,l (Z).
Since the symmetry condition can be expressed as r j=1 φk,l (Z) = l =j k∈N j,l φk,l (Z), the condition ∇V (Z) = 0 is then equivalent to [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF]. Note that we have assumed the functions φ k to be smooth for simplicity, but a similar reasoning can be conducted for non-smooth functions, by using the concept of subgradients.

Examples of potential games.

• the 1 -norm: with r = 1 and φ1,j = z j 1 , since problem j does not involve any variable z l for l = j; • Symmetric TV / Laplacian: problem j may involve a variable z l through a term a j,l z j -z l 1 .

Then, problem l involves the same term a l,j z j -z l 1 under the condition a j,l = a l,j . • Symmetric non local group with r = p and φk,j = λ k [ √ a j,1 z 1 (k), . . . , √ a j,m z m (k)] 2 . Under the condition of symmetric weights a j,l = a l,j , we obtain again a potential game. Potential games are appealing as they provide guarantees about the existence of Nash equilibria without requiring optimizing over a compact set. Yet, we have found that allowing non-symmetric weights often performs better. This is illustrated in Table A1 for a simple denoising experiment. whose size equals the window size used for computing non-local self-similarities; whereas we train pixel level models on the full size images. For fMRI experiments we also trained the models on the full sized images. We apply a mild data augmentation (random rotation by 90 • and horizontal flips).

We optimize the parameters of our models using ADAM [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF]. The learning rate is set to 6 × 10 -4 at initialization and is sequentially lowered during training by a factor of 0.35 every 80 training steps, in the same way for all experiments. Similar to [START_REF] Simon | Rethinking the CSC model for natural images[END_REF], we normalize the initial dictionary D 0 by its largest singular value as explained in the main paper in Section 3.4. We initialize the dictionary C,D and W with the same dictionary obtained with an unsupervised dictionary learning algorithm (using SPAMS library).

We have implemented the backtracking strategy described in Section 3.4 of the main paper for all our algorithms, which automatically decreases the learning rate by a factor 0.8 when the loss function increases too much on the training set, and restore a previous snapshot of the model. Divergence is monitored by computing the loss on the training set every 10 epochs. Training the non-local models for denoising are the longer models to train and takes about 2 days on a Titan RTX GPU. We summarize the chosen hyperparameters for the experiments in Table A2.

C Additional Quantitative Results

C.1 Inference speed

In Table A3 we provide a comparison of our TV models in terms of speed with BM3D for grayscale denoising on the BSD68 dataset. For fair comparison, we reported computation time both on gpu and cpu.

C.2 Image denoising

We provide additional results for grayscale denoising with different variations of the prior introduced in the main paper, as well as combination of different priors. We reported performances for gray denoising in Table A4 for the pixel based models, and in Table A5 for the patch based models. In Table A4 untied κ denotes when we used a different set of learned parameters κ at each stage of the refinement step of the similarity matrix for the non-local models.

Figure 2 :

 2 Figure 2: Architecture of our trainable models for image restoration.

Figure A3 :Figure A4 :

 A3A4 Figure A3: Grayscale denoising for 4 images from the BSD68 dataset. Best seen by zooming on a computer screen.

  Algorithm 1 Pseudocode of the general training procedure for image restoration

1: Sample a minibatch of pairs of corrupted/clean images {(x 0 , y 0 ), • • • , (x K , y K )}; 2: Extract overlapping patches of corrupted images to form tensors

X i = [P 1 x i , • • • , P n x i ]; 3: for t = 1, 2, . .

. , K do

Compute an approximate Nash equilibrium Z of the convex games 4:

Table 2 :

 2 Gradient descent (GD) vs. Extra-gradient. Denoising results in avg. PSNR with σ = 25 on BSD68[START_REF] Martin | A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics[END_REF].

	Method	GD (24 iters) GD (48 iters) Extra-gradient (24 iters)
	Trainable TV symmetric	27.58	27.50	27.82
	Trainable TV assymetric	27.99	27.89	28.24

Table 3

 3

	: Study of stabilization techniques for
	learnt sparse coding. Denoising results in aver-
	age PSNR with σ = 25 on BSD68.
	Method	Psnr (dB) D C,D
	BM3D [11]	28.57
	Sparse Coding (SC)
	SC + Backtracking	28.71 28.83
	SC + Spectral norm	28.69 28.82
	SC + Barzilai-Borwein 28.82 28.86

Table 4 :

 4 Grayscale denoising on BSD68, training on BSD400 for all methods. Performance in terms of average PSNR. Tiny CNN is a CNN baseline with few parameters. See Appendix C for qualitative results.

	Method	Params	5	Noise Level (σ) 15 25	50
	Tiny CNN (ours)	326	35.17 29.42 26.90 24.06
	Tiny CNN (ours)	1200	36.47 30.36 27.70 24.60
	BM3D [11]	-	37.57 31.07 28.57 25.62
	LSCC [34]	-	37.70 31.28 28.71 25.72
	CSCnet [52]	62k	37.69 31.40	28.93 26.04
	GroupSC [26]	68k	37.95 31.71 29.20 26.17
	FFDNet [64]	486k	N/A	31.63 29.19 26.29
	DnCNN [63]	556k	37.68 31.73 29.22 26.23
	NLRN [29]	330k	37.92 31.88 29.41 26.47
	Pixel-reconstruction
	TV symmetric	288	36.91 30.27 27.66 24.51
	TV assymetric -extra-grad	480	37.30 30.76 28.24 25.32
	Laplacian symmetric	288	35.17 28.42 26.14 23.70
	Laplacian assymetric -extra-grad	480	35.20 28.46 26.39 23.77
	Bilateral -extra-grad	146	36.75 29.89 27.20 23.72
	Bilateral TV -extra-grad	146	36.94 30.46 27.78 24.52
	Non-local TV -extra-grad	307	37.53 31.03 28.50 25.26
	Non-local Laplacian -extra-grad	307	37.54 31.00 28.47 25.46
	Patch-reconstruction
	Sparse Coding (SC)	68k	37.84 31.46 28.90 25.84
	Sparse Coding + Variance	68k	37.83 31.49 29.00 26.08
	Sparse Coding + TV	68k	37.84 31.50 29.02 26.10
	Sparse Coding + TV + Variance	68k	37.84 31.51 29.03 26.09
	Non-local group	68k	37.95 31.69 29.19 26.19
	Non-local group + Variance	68k	37.96 31.70 29.22 26.28
	GroupSC + Variance	68k	37.96 31.75 29.24 26.34
	Compressed Sensing for fMRI. Compressed Sensing for functional magnetic resonance imaging
	(fMRI) aims at reconstructing functional MR images from a small number of samples in the Fourier
	space. The corresponding inverse problem is
	min y∈R n Ay -x 2 2 + λΨ(y),

Table 5 :

 5 Compressed sensing for fMRI on the MR brain dataset using a pseudo radial sampling pattern. Performance comparisons in terms of PSNR (dB).

	Method	Params	20 %	30 %	40%	50%	Test time
	TV [30]	-	35.20	37.99	40.00	41.69 0.731s (cpu)
	RecPF [59]	-	35.32	38.06	40.03	41.71 0.315s (cpu)
	SIDWT	-	35.66	38.72	40.88	42.67 7.867s (cpu)
	PANO [46]	-	36.52	39.13	40.31	41.81 35.33s (cpu)
	BM3D-MRI [13]	-	37.98	40.33	41.99	43.47 40.91s (cpu)
	ADMM-net [53]	-	37.17	39.84	41.56	43.00 0.791s (cpu)
	ISTA-net [62]	337k	38.73 40.89	42.52 44.09	0.143s (gpu)
	CS-TV (ours)	140	36.80	39.63	41.58	43.46 0.015s (gpu)
	CS-Sparse coding (ours)	68k	37.80	40.50	42.46	44.16 0.213s (gpu)
	CS-Sparse coding + Variance (ours)	68k	37.79 40.67	42.54	44.17 0.213s (gpu)

Table 6 :

 6 Denoising with less data. Results in terms of average PSNR(dB) on BSD68 with σ = 15. All the models are trained on a similar subset of BSD400 for fair comparaison.

	Method	Params	400	Training images 200 100	50
	DnCNN [63]	556k 31.73 31.65 31.47 31.23
	TV extra-grad	480 30.75 30.72 30.67 30.66
	SC+Var	68k 31.49 31.49 31.47 31.40
	GroupSC+Var	68k 31.75 31.66 31.62 31.54

Table 7 :

 7 Dense stereo matching fine-tuning on kitti2015 train set, performance reported on the kitti2015 validation set.

	Model	3-px error (%)
	PSMNet [8]	2.14 ± 0.04
	PSMNet+TV 12	2.11 ± 0.03
	PSMNet+TV 24	2.11 ± 0.04
	PSMNet+TV extra	2.10 ± 0.03

Table A1 :

 A1 Symmetric vs assymmetric grayscale denoising on BSD68, training on BSD400 for all methods. Performance is measured in terms of average PSNR.

	Method	Params	5	Noise Level (σ) 15 25	50
	TV symmetric	72	36.08 30.21 27.58 24.74
	TV assymetric -extra-grad	480	37.30 30.76 28.24 25.32
	Laplacian symmetric	72	34.88 28.14 25.90 23.45
	Laplacian assymetric -extra-grad	480	35.20 28.46 26.39 23.77
	Non-local group -symmetric	68k	37.94 31.67 29.17 26.16
	Non-local group -assymetric	68k	37.95 31.69 29.20 26.19
	B Implementation Details and Reproducibility

B.1 Training Details

For the training of patch-based models for denoising, we randomly extract patches of size 56 × 56

Table A2 :

 A2 Hyper-parameters chosen for every task.

	Experiment	Gray denoising (patch) Gray denoising (pixel) fMRI
	Patch size	9	-	9
	Dictionary size	256	-	256
	Nr epochs	300	300	150
	Batch size	32	32	1
	K iterations	24	24	24
	Middle averaging			-
	Correlation update frequency f	1/6	1/12	-

Table A3 :

 A3 Inference speed for image denoising.

		Params	Psnr	Speed
	BM3D [11]	-	25.62	7.28s (cpu)
	TV assymetric	240	24.93 0.014s (gpu) / 0.18s (cpu)
	TV assymetric (extra)	480	25.32 0.021s (gpu) / 0.28s (cpu)

Table A4 :

 A4 Pixel level grayscale denoising on BSD68, training on BSD400 for all models. Performance is measured in terms of average PSNR.

	Method	Params	5	Noise Level (σ) 15 25	50
	BM3D [11]	-	37.57 31.07 28.57 25.62
	Tiny CNN	326	35.17 29.42 26.90 24.06
	Tiny CNN	1200	36.47 30.36 27.70 24.60
	TV symmetric	288	36.08 30.21 27.58 24.74
	TV symmetric -extra-grad	144	37.02 30.33 27.82 24.81
	TV assymetric-	240	36.83 30.49 27.99 24.93
	TV assymetric -extra-grad	480

We employ a debiasing dictionary W = D to improve the quality of the reconstructions. Debiaising is commonly used when dealing with 1 penalty which is known to shrink the coefficients Z too much.

torch.nn.functionnal.conv2D_transpose on PyTorch[START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF].

Acknowledgments and Disclosure of Funding

This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). JM and BL were supported by the ERC grant number 714381 (SOLARIS project) and by ANR 3IA MIAI@Grenoble Alpes (ANR-19-P3IA-0003). JP was supported in part by the Louis Vuitton/ENS chair in artificial intelligence and the Inria/NYU collaboration. Finally, this work was granted access to the HPC resources of IDRIS under the allocation 2020-AD011011252 made by GENCI.

Broader Impact

Our main field of application is in image processing, with a focus on image restoration and reconstruction, whose benefits for society are clear and well established, even though a misuse of such a technology is of course possible; for instance, the same total variation penalty may be used in medical imaging, for personal photography, in astronomical imaging, or for restoring images produced by military devices. More specifically, our paper is addressing the issue of interpretability of neural networks, by considering models admitting a functional (mathematical) description of their decision functions, and with less parameters than classical deep learning models. As we mentioned in the discussion section, we believe that these are first steps to build systems producing explanable decisions, and that are more data and energy efficient. These are important issues going beyond image processing, which we would like to address in future work.