Hamza Cherkaoui
email: hamza.cherkaoui@cea.fr

Jeremias Sulam
email: jsulam1@jhu.edu

Thomas Moreau
email: thomas.moreau@inria.fr

Learning to solve TV regularized problems with unrolled algorithms

Total Variation (TV) is a popular regularization strategy that promotes piece-wise constant signals by constraining the 1 -norm of the first order derivative of the estimated signal. The resulting optimization problem is usually solved using iterative algorithms such as proximal gradient descent, primal-dual algorithms or ADMM. However, such methods can require a very large number of iterations to converge to a suitable solution. In this paper, we accelerate such iterative algorithms by unfolding proximal gradient descent solvers in order to learn their parameters for 1D TV regularized problems. While this could be done using the synthesis formulation, we demonstrate that this leads to slower performances. The main difficulty in applying such methods in the analysis formulation lies in proposing a way to compute the derivatives through the proximal operator. As our main contribution, we develop and characterize two approaches to do so, describe their benefits and limitations, and discuss the regime where they can actually improve over iterative procedures. We validate those findings with experiments on synthetic and real data.

Introduction

Ill-posed inverse problems appear naturally in signal and image processing and machine learning, requiring extra regularization techniques. Total Variation (TV) is a popular regularization strategy with a long history [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], and has found a large number of applications in neuro-imaging [START_REF] Fikret | Total activation: fMRI deconvolution through spatio-temporal regularization[END_REF], medical imaging reconstruction [START_REF] Tian | Low-dose CT reconstruction via edge-preserving total variation regularization[END_REF], among myriad applications [START_REF] Rodríguez | Total Variation Regularization Algorithms for Images Corrupted with Different Noise Models: A Review[END_REF][START_REF] Darbon | Image Restoration with Discrete Constrained Total Variation Part I: Fast and Exact Optimization[END_REF]. TV promotes piece-wise constant estimates by penalizing the 1 -norm of the first order derivative of the estimated signal, and it provides a simple, yet efficient regularization technique.

TV-regularized problems are typically convex, and so a wide variety of algorithms are in principle applicable. Since the 1 norm in the TV term is non-smooth, Proximal Gradient Descent (PGD) is the most popular choice [START_REF] Rockafellar | Monotone Operators and the Proximal Point Algorithm[END_REF]). Yet, the computation for the corresponding proximal operator (denoted prox-TV) represents a major difficulty in this case as it does not have a closed-form analytic solution. For 1D problems, it is possible to rely on dynamic programming to compute prox-TV, such as the taut string algorithm [START_REF] Davies | Local Extremes, Runs, Strings and Multiresolution[END_REF]Condat, 2013a). Another alternative consists in computing the proximal operator with iterative first order algorithm [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF][START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF][START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF][START_REF] Condat | A Primal-Dual Splitting Method for Convex Optimization Involving Lipschitzian, Proximable and Linear Composite Terms[END_REF]. Other algorithms to solve TV-regularized 34th Conference on Neural Information Processing Systems (NeurIPS 2020), Vancouver, Canada.

problems rely on primal dual algorithms [START_REF] Chambolle | A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging[END_REF][START_REF] Condat | A Primal-Dual Splitting Method for Convex Optimization Involving Lipschitzian, Proximable and Linear Composite Terms[END_REF] or Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers[END_REF]. These algorithms typically use one sequence of estimates for each term in the objective and try to make them as close as possible while minimizing the associated term. While these algorithms are efficient for denoising problems -where one is mainly concerned with good reconstruction -they can result in estimate that are not very well regularized if the two sequences are not close enough.

When on fixed computational budget, iterative optimization methods can become impractical as they often require many iterations to give a satisfactory estimate. To accelerate the resolution of these problems with a finite (and small) number of iterations, one can resort to unrolled and learned optimization algorithms (see [START_REF] Monga | Algorithm Unrolling: Interpretable, Efficient Deep Learning for Signal and Image Processing[END_REF] for a review). In their seminal work, [START_REF] Gregor | Learning Fast Approximations of Sparse Coding[END_REF] proposed the Learned ISTA (LISTA), where the parameters of an unfolded Iterative Shrinkage-Thresholding Algorithm (ISTA) are learned with gradient descent and back-propagation. This allows to accelerate the approximate solution of a Lasso problem [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF], with a fixed number of iteration, for signals from a certain distribution. The core principle behind the success of this approach is that the network parameters can adaptively leverage the sensing matrix structure [START_REF] Moreau | Understanding Neural Sparse Coding with Matrix Factorization[END_REF] as well as the input distribution [START_REF] Giryes | Tradeoffs between Convergence Speed and Reconstruction Accuracy in Inverse Problems[END_REF][START_REF] Ablin | Learning step sizes for unfolded sparse coding[END_REF]. Many extensions of this original idea have been proposed to learn different algorithms [START_REF] Sprechmann | Learning Efficient Structured Sparse Models[END_REF][START_REF] Sprechmann | Efficient Supervised Sparse Analysis and Synthesis Operators[END_REF][START_REF] Borgerding | AMP-Inspired Deep Networks for Sparse Linear Inverse Problems[END_REF] or for different classes of problem [START_REF] Xin | Maximal Sparsity with Deep Networks?[END_REF][START_REF] Giryes | Tradeoffs between Convergence Speed and Reconstruction Accuracy in Inverse Problems[END_REF][START_REF] Sulam | On Multi-Layer Basis Pursuit, Efficient Algorithms and Convolutional Neural Networks[END_REF]. The motif in most of these adaptations is that all operations in the learned algorithms are either linear or separable, thus resulting in sub-differentials that are easy to compute and implement via back-propagation. Algorithm unrolling is also used in the context of bi-level optimization problems such as hyper-parameter selection. Here, the unrolled architecture provides a way to compute the derivative of the inner optimization problem solution compared to another variable such as the regularisation parameter using back-propagation [START_REF] Bertrand | Implicit differentiation of Lasso-type models for hyperparameter optimization[END_REF].

The focus of this paper is to apply algorithm unrolling to TV-regularized problems in the 1D case. While one could indeed apply the LISTA approach directly to the synthesis formulation of these problems, we show in this paper that using such formulation leads to slower iterative or learned algorithms compared to their analysis counterparts. The extension of learnable algorithms to the analysis formulation is not trivial, as the inner proximal operator does not have an analytical or separable expression. We propose two architectures that can learn TV-solvers in their analysis form directly based on PGD. The first architecture uses an exact algorithm to compute the prox-TV and we derive the formulation of its weak Jacobian in order to learn the network's parameters. Our second method rely on a nested LISTA network in order to approximate the prox-TV itself in a differentiable way. This latter approach can be linked to inexact proximal gradient methods [START_REF] Schmidt | Convergence Rates of Inexact Proximal-Gradient Methods for Convex Optimization[END_REF][START_REF] Machart | Optimal Computational Trade-Off of Inexact Proximal Methods[END_REF]. These results are backed with numerical experiments on synthetic and real data. Concurrently to our work, [START_REF] Lecouat | Designing and Learning Trainable Priors with Non-Cooperative Games[END_REF] also proposed an approach to differentiate the solution of TV-regularized problems. While their work can be applied in the context of 2D signals, they rely on smoothing the regularization term using Moreau-Yosida regularization, which results in smoother estimates from theirs learned networks. In contrast, our work allows to compute sharper signals but can only be applied to 1D signals.

The rest of the paper is organized as follows. In Section 2, we describe the different formulations for TV-regularized problems and their complexity. We also recall central ideas of algorithm unfolding. Section 3 introduces our two approaches for learnable network architectures based on PGD. Finally, the two proposed methods are evaluated on real and synthetic data in Section 4.

Notations For a vector x ∈ R k , we denote x q its q -norm. For a matrix A ∈ R m×k , we denote A 2 its 2 -norm, which corresponds to its largest singular value and A † denotes its pseudoinverse. For an ordered subset of indices S ⊂ {1, . . . , k}, x S denote the vector in R |S| with element (x S ) t = x it for i t ∈ S. For a matrix A ∈ R m×k , A :,S denotes the sub-matrix [A :,i1 , . . . A :,i |S| ] composed with the columns A :,it of index i t ∈ S of A. For the rest of the paper, we refer to the operators

D ∈ R k-1×k , D ∈ R k×k , L ∈ R k×k and R ∈ R k×k as: D =      -1 1 0 . . . 0 0 -1 1 . . . . . . . . . . . . . . . . . . 0 0 . . . 0 -1 1      D =      1 0 . . . 0 -1 1 . . . . . . . . . . . . 0 0 -1 1      L =      1 0 . . . 0 1 1 . . . . . . . . . . . . . . . 0 1 . . . 1 1      R =      0 0 . . . 0 0 1 . . . . . . . . . . . . . . . 0 0 . . . 0 1      2
2 Solving TV-regularized problems

We begin by detailing the TV-regularized problem that will be the main focus of our work. Consider a latent vector u ∈ R k , a design matrix A ∈ R m×k and the corresponding observation x ∈ R m . The original formulation of the TV-regularized regression problem is referred to as the analysis formulation [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. For a given regularization parameter λ > 0, it reads

min u∈R k P (u) = 1 2 x -Au 2 2 + λ u T V , (1)
where u T V = Du 1 , and D ∈ R k-1×k stands for the first order finite difference operator, as defined above. The problem in (1) can be seen as a special case of a Generalized Lasso problem [START_REF] Tibshirani | The solution path of the generalized lasso[END_REF]; one in which the analysis operator is D. Note that problem P is convex, but the T V -norm is non-smooth. In these cases, a practical alternative is the PGD, which iterates between a gradient descent step and the prox-TV. This algorithm's iterates read

u (t+1) = prox λ ρ • T V u (t) - 1 ρ A (Au (t) -x) , (2)
where ρ = A 2 2 and the prox-TV is defined as

prox µ • T V (y) = arg min u∈R k F y (u) = 1 2 y -u 2 2 + µ u T V . (3)
Problem (3) does not have a closed-form solution, and one needs to resort to iterative techniques to compute it. In our case, as the problem is 1D, the prox-TV problem can be addressed with a dynamic programming approach, such as the taut-string algorithm (Condat, 2013a). This scales as O(k) in all practical situations and is thus much more efficient than other optimization based iterative algorithms [START_REF] Rockafellar | Monotone Operators and the Proximal Point Algorithm[END_REF][START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF][START_REF] Condat | A Primal-Dual Splitting Method for Convex Optimization Involving Lipschitzian, Proximable and Linear Composite Terms[END_REF] for which each iteration is O(k 2 ) at best.

With a generic matrix A ∈ R m×k , the PGD algorithm is known to have a sublinear convergence rate [START_REF] Combettes | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF]. More precisely, for any initialization u (0) and solution u * , the iterates satisfy

P (u (t) ) -P (u * ) ≤ ρ 2t u (0) -u * 2 2 , (4)
where u * is a solution of the problem in (1). Note that the constant ρ can have a significant effect. Indeed, it is clear from (4) that doubling ρ leads to consider doubling the number of iterations.

Synthesis formulation

An alternative formulation for TV-regularized problems relies on removing the analysis operator D from the 1 -norm and translating it into a synthesis expression [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF]. Removing D from the non-smooth term simplifies the expression of the proximal operator by making it separable, as in the Lasso. The operator D is not directly invertible but keeping the first value of the vector u allows for perfect reconstruction. This motivates the definition of the operator D ∈ R k×k , and its inverse L ∈ R k×k , as defined previously. Naturally, L is the discrete integration operator. Considering the change of variable z = Du, and using the operator R ∈ R k×k , the problem in (1) is equivalent to

min z∈R k S(z) = 1 2 x -ALz 2 2 + λ Rz 1 . (5)
Note that for any z ∈ R k , S(z) = P (Lz). There is thus an exact equivalence between solutions from the synthesis and the analysis formulation, and the solution for the analysis can be obtained with u * = Lz * . The benefit of this formulation is that the problem above now reduces to a Lasso problem [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF]. In this case, the PGD algorithm is reduced to the ISTA with a closed-form proximal operator (the soft-thresholding). Note that this simple formulation is only possible in 1D where the first order derivative space is unconstrained. In larger dimensions, the derivative must be constrained to verify the Fubini's formula that enforces the symmetry of integration over dimensions.

While it is also possible to derive synthesis formulation in higher dimension [START_REF] Elad | Analysis versus synthesis in signal priors[END_REF], this does not lead to simplistic proximal operator. Dimension k

10 1 10 3 10 5 AL 2 2 / A 2 2 Mean E AL 2 2 A 2 2 Proposition 2.1 Conjecture 2.2 Figure 1: Evolution of E AL 2 2 A 2 2
w.r.t the dimension k for random matrices A with i.i.d normally distributed entries. In light blue is the confidence interval [0.1, 0.9] computed with the quantiles. We observe that it scales as O(k 2 ) and that our conjectured bound seems tight.

For this synthesis formulation, with a generic matrix A ∈ R m×k , the PGD algorithm has also a sublinear convergence rate [START_REF] Beck | A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems[END_REF] such that

P (u (t) ) -P (u * ) ≤ 2 ρ t u (0) -u * 2 2 , (6)
with ρ = AL 2 2 (see Subsection F.1 for full derivation). While the rate of this algorithm is the same as in the analysis formulation -in O( 1 t ) -the constant ρ related to the operator norm differs. We now present two results that will characterize the value of ρ.

E AL 2 2 A 2 2 ≥ 2k + 1 4π 2 + o(1)
The full proof can be found in Subsection F.3. The lower bound is constructed by using A T A A 2 2 u 1 u 1 for a unit vector u 1 and computing explicitely the expectation for rank one matrices. To assess the tightness of this bound, we evaluated numerically E AL 2 2 A 2 2 on a set of 1000 matrices sampled with i.i.d normally distributed entries. The results are displayed w.r.t the dimension k in Figure 1. It is clear that the lower bound from Proposition 2.1 is not tight. This is expected as we consider only the leading eigenvector of A to derive it in the proof. The following conjecture gives a tighter bound.

Conjecture 2.2 (Expectation for the ratio

AL 2 2 A 2 2
). Under the same conditions as in Proposition 2.1, the expectation of AL 2 2 / A 2 2 is given by

E AL 2 2 A 2 2 = (2k + 1) 2 16π 2 + o(1) .
We believe this conjecture can potentially be proven with analogous developments as those in Proposition 2.1, but integrating over all dimensions. However, a main difficulty lies in the fact that integration over all eigenvectors have to be carried out jointly as they are not independent. This is subject of current ongoing work.

Finally, we can expect that ρ/ρ scales as Θ(k 2 ). This leads to the observation that ρ 2 ρ in large enough dimension. As a result, the analysis formulation should be much more efficient in terms of iterations than the synthesis formulation -as long as the prox-TVcan be dealt with efficiently.

Unrolled iterative algorithms

As shown by [START_REF] Gregor | Learning Fast Approximations of Sparse Coding[END_REF], ISTA is equivalent to a recurrent neural network (RNN) with a particular structure. This observation can be generalized to PGD algorithms for any penalized least squares problem of the form

u * (x) = arg min u L(x, u) = 1 2 x -Bu 2 2 + λg(u) , (7)
W x x prox µg u * W u (a) PGD -Recurrent Neural Network x W (0) x prox µ (1) g W (1) u W (1) x prox µ (2) g W (2) u W (2) x prox µ (3) g u (3) (b) LPGD -Unfolded network for Learned PGD with T = 3
Θ (t) = {W (t)
x , W

u , µ (t) } can be learned by minimizing the loss (10) to approximate good solution of (7) on average.

where g is proper and convex, as depicted in Figure 2a. By unrolling this architecture with T layers, we obtain a network φ Θ (T ) (x) = u (T ) -illustrated in Figure 2b -with parameters

Θ (T ) = {W (t) x , W (t) u , µ (t) } T
t=1 , defined by the following recursion

u (0) = B † x ; u (t) = prox µ (t) g (W (t) x x + W (t) u u (t-1) ) . (8)
As underlined by ( 4), a good estimate u (0) is crucial in order to have a fast convergence toward u * (x). However, this chosen initialization is mitigated by the first layer of the network which learns to set a good initial guess for u (1) . For a network with T layers, one recovers exactly the T -th iteration of PGD if the weights are chosen constant equal to

W (t) x = 1 ρ B , W (t) u = (Id - 1 ρ B B) , µ (t) = λ ρ , with ρ = B 2 2 . (9)
In practice, this choice of parameters are used as initialization for a posterior training stage. In many practical applications, one is interested in minimizing the loss (7) for a fixed B and a particular distribution over the space of x, P. As a result, the goal of this training stage is to find parameters Θ (T ) that minimize the risk, or expected loss, E[L(x, φ Θ (T ) (x))] over P. Since one does not have access to this distribution, and following an empirical risk minimization approach with a given training set {x 1 , . . . x N } (assumed sampled i.i.d from P), the network is trained by minimizing min

Θ (T ) 1 N N i=1 L(x i , φ Θ (T ) (x i )) . ( 10
)
Note that when T → +∞, the presented initialization in (9) gives a global minimizer of the loss for all x i , as the network converges to exact PGD. When T is fixed, however, the output of the network is not a minimizer of (7) in general. Minimizing this empirical risk can therefore find a weight configuration that reduces the sub-optimality of the network relative to (7) over the input distribution used to train the network. In such a way, the network learns an algorithm to approximate the solution of (7) for a particular class or distributions of signals. It is important to note here that while this procedure can accelerate the resolution the problem, the learned algorithm will only be valid for inputs x i coming from the same input distribution P as the training samples. The algorithm might not converge for samples which are too different from the training set, unlike the iterative algorithm which is guaranteed to converge for any sample.

This network architecture design can be directly applied to TV regularized problems if the synthesis formulation ( 5) is used. Indeed, in this case PGD reduces to the ISTA algorithm, with B = AL and prox µg = ST(•, µ) becomes simply a soft-thresholding operator (which is only applied on the coordinates {2, . . . k}, following the definition of R). However, as discussed in Proposition 2.1, the conditioning of the synthesis problem makes the estimation of the solution slow, increasing the number of network layers needed to get a good estimate of the solution. In the next section, we will extend these learning-based ideas directly to the analysis formulation by deriving a way to obtain exact and approximate expressions for the sub-differential of the non-separable prox-TV.

3 Back-propagating through TV proximal operator

Our two approaches to define learnable networks based on PGD for TV-regularized problems in the analysis formulation differ on the computation of the prox-TV and its derivatives. Our first approach consists in directly computing the weak derivatives of the exact proximal operator while the second one uses a differentiable approximation.

Derivative of prox-TV

While there is no analytic solution to the prox-TV, it can be computed exactly (numerically) for 1D problems using the taut-string algorithm (Condat, 2013a). This operator can thus be applied at each layer of the network, reproducing the architecture described in Figure 2b. We define the LPGD-Taut network φ Θ (T ) (x) with the following recursion formula

φ Θ (T ) (x) = prox µ (T ) • T V W (T ) x x + W (T ) u φ Θ (T -1) (x) (11)
To be able to learn the parameters through gradient descent, one needs to compute the derivatives of (10) w.r.t the parameters Θ (T ) . Denoting h = W (t)

x x + W (t) u φ Θ (t-1) (x) and u = prox µ (t) • T V (h)
, the application of the chain rule (as implemented efficiently by automatic differentiation) results in

∂L ∂h = J x (h, µ (t) ) ∂L ∂u , and
∂L ∂µ (t) = J µ (h, µ (t) ) ∂L ∂u , (12)
where J x (h, µ) ∈ R k×k and J µ (h, µ) ∈ R k×1 denotes the weak Jacobian of the output of the proximal operator u with respect to the first and second input respectively. We now give the analytic formulation of these weak Jacobians in the following proposition.

Proposition 3.1. [Weak Jacobian of prox-TV] Let x ∈ R k and u = prox µ • T V (x)
, and denote by S the support of z = Du. Then, the weak Jacobian J x and J µ of the prox-TV relative to x and µ can be computed as

J x (x, µ) = L :,S (L :,S L :,S ) -1 L :,S and J µ (x, µ) = -L :,S (L :,S L :,S ) -1 sign(Du) S
The proof of this proposition can be found in Subsection G.1. Note that the dependency in the inputs is only through S and sign(Du), where u is a short-hand for prox µ • T V (x). As a result, computing these weak Jacobians can be done efficiently by simply storing sign(Du) as a mask, as it would be done for a RELU or the soft-thresholding activations, and requiring just 2(k -1) bits. With these expressions, it is thus possible to compute gradient relatively to all parameters in the network, and employ them via back-propagation.

Unrolled prox-TV

As an alternative to the previous approach, we propose to use the LISTA network to approximate the prox-TV (3). The prox-TV can be reformulated with a synthesis approach resulting in a Lasso i.e.

z * = arg min z 1 2 h -Lz 2 2 + µ Rz 1 (13)
The proximal operator solution can then be retrieved with prox µ • T V (h) = Lz * . This problem can be solved using ISTA, and approximated efficiently with a LISTA network [START_REF] Gregor | Learning Fast Approximations of Sparse Coding[END_REF].

For the resulting architecture -dubbed LPGD-LISTA -prox µ • T V (h) is replaced by a nested LISTA network with a fixed number of layers T in defined recursively with z (0) = Dh and

z ( +1) = ST W ( ,t) z z ( ) + W ( ,t) h Φ Θ (t) , µ ( ,t) ρ . (14)
Here, W

( ,t) z , W ( ,t) h
, µ ( ,t) are the weights of the nested LISTA network for layer . They are initialized with weights chosen as in (9) to ensure that the initial state approximates the prox-TV. Note that the weigths of each of these inner layers are also learned through back-propagation during training.

The choice of this architecture provides a differentiable (approximate) proximal operator. Indeed, the LISTA network is composed only of linear and soft-thresholding layers -standard tools for deep-learning libraries. The gradient of the network's parameters can thus be computed using classic automatic differentiation. Moreover, if the inner network is not trained, the gradient computed with this method will converge toward the gradient computed using Proposition 3.1 as T in goes to ∞ (see Proposition G.2). Thus, in this untrained setting with infinitely many inner layers, the network is equivalent to LPGD-Taut as the output of the layer also converges toward the exact proximal operator. We see that synthesis formulations are outperformed by the analysis counter part. Both our methods are able to accelerate the resolution of (20), at least in the first iterations.

Connections to inexact PGD A drawback of approximating the prox-TV via an iterative procedure is, precisely, that it is not exact. This optimization error results from a trade-off between computational cost and convergence rate. Using results from [START_REF] Machart | Optimal Computational Trade-Off of Inexact Proximal Methods[END_REF], one can compute the scaling of T and T in to reach an error level of δ with an untrained network. Proposition G.3 shows that without learning, T should scale as O( 1t ) and T in should be larger than O(ln(1 δ )). This scaling gives potential guidelines to set these parameters, as one can expect that learning the parameters of the network would reduce these requirement.

Experiments

All experiments are performed in Python using PyTorch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF]. We used the implementation 1 of [START_REF] Barbero | Modular proximal optimization for multidimensional total-variation regularization[END_REF] to compute TV proximal operator using taut-string algorithm. The code to reproduce the figures is available online2 .

In all experiments, we initialize u 0 = A † x. Moreover, we employed a normalized λ reg as a penalty parameter: we first compute the value of λ max (which is the minimal value for which z = 0 is solution of ( 5)) and we refer to λ as the ratio so that λ reg = λλ max , with λ ∈ [0, 1] (see Appendix D). As the computational complexity of all compared algorithms is the same except for the proximal operator, we compare them in term of iterations.

Simulation

We generate n = 2000 times series and used half for training and other half for testing and comparing the different algorithms. We train all the network's parameters jointly -those to approximate the gradient for each iteration along with those to define the inner proximal operator. The full training process is described in Appendix A. We set the length of the source signals (u i ) n i=1 ∈ R n×k to k = 8 with a support of |S| = 2 non-zero coefficients (larger dimensions will be showcased in the real data application). We generate A ∈ R m×k as a Gaussian matrix with m = 5, obtaining then (u i ) n i=1 ∈ R n×p . Moreover, we add Gaussian noise to measurements x i = Au i with a signal to noise ratio (SNR) of 1.0.

We compare our proposed methods, LPGD-Taut network and the LPGD-LISTA with T in = 50 inner layers to PGD and Accelerated PGD with the analysis formulation. For completeness, we also add the FISTA algorithm for the synthesis formulation in order to illustrate Proposition 2.1 along with its learned version.

Figure 3 presents the risk (or expected function value, P ) of each algorithm as a function of the number of layers or, equivalently, iterations. For the learned algorithms, the curves in t display the performances of a network with t layer trained specifically. We observe that all the synthesis formulation algorithms are slower than their analysis counterparts, empirically validating Proposition 2.1. More interestingly, in a high sparsity context, after a certain point, the error sharply increase.

Moreover, both of the proposed methods accelerate the resolution of (20) in a low iteration regime. However, when the regularization parameter is high (λ = 0.8), we observe that the performance of the LPGD-LISTA tends to plateau. It is possible that such a high level of sparsity require more than 50 layers for the inner network (which computes the prox-TV). According to Section 3.2, the error associated with this proximity step hinders the global convergence, making the loss function decrease slowly. Increasing the number of inner layers would alleviate this issue, though at the expense of increased computational burden for both training and runtime. For LPGD-Taut, while the Taut-string algorithm ensures that the recovered support is exact for the proximal step, the overall support can be badly estimated in the first iterations. This can lead to un-informative gradients as they greatly depend on the support of the solution in this case, and explain the reduced performances of the network in the high sparsity setting.

Inexact prox-TV With the same data (x i ) n i=1 ∈ R n×m , we empirically investigate the error of the prox-TV

(t) k = F u (t) (z (t) ) -F u (t) (z *
) and evaluate it for c with different number of layers (T ∈ [20, 50]). We also investigate the case where the parameter of the nested LISTA in LPGD-LISTA are trained compared to their initialization in untrained version.

Figure 4 depicts the error k for each layer. We see that learning the parameters of the unrolled prox-TV in LPGD-LISTA barely improves the performance. More interestingly, we observe that in a high sparsity setting the error sharply increases after a certain number of layers. This is likely cause by the high sparsity of the estimates, the small numbers of iterations of the inner network (between 20 and 50) are insufficient to obtain an accurate solution to the proximal operator. This is in accordance with inexact PGD theory which predict that such algorithm has no exact convergence guarantees [START_REF] Schmidt | Convergence Rates of Inexact Proximal-Gradient Methods for Convex Optimization[END_REF].

fMRI data deconvolution

Functional magnetic resonance imaging (fMRI) is a non-invasive method for recording the brain activity by dynamically measuring blood oxygenation level-dependent (BOLD) contrast, denoted here x. The latter reflects the local changes in the deoxyhemoglobin concentration in the brain [START_REF] Ogawa | Intrinsic signal changes accompanying sensory stimulation: Functional brain mapping with magnetic resonance imaging[END_REF] and thus indirectly measures neural activity through the neurovascular coupling. This coupling is usually modelled as a linear and time-invariant system and characterized by its impulse response, the so-called haemodynamic response function (HRF), denoted here h. Recent developments propose to estimate either the neural activity signal independently [START_REF] Fikret | Total activation: fMRI deconvolution through spatio-temporal regularization[END_REF]Cherkaoui et al., 2019b) or jointly with the HRF (Cherkaoui et al., 2019a;[START_REF] Farouj | Bold Signal Deconvolution Under Uncertain HAEModynamics: A Semi-Blind Approach[END_REF]. Estimating the neural activity signal with a fixed HRF is akin to a deconvolution problem regularized with TV-norm,

min u∈R k P (u) = 1 2 h * u -x 2 2 + λ u T V (15)
To demonstrate the usefulness of our approach with real data, where the training set has not the exact same distribution than the testing set, we compare the LPGD-Taut to Accelerated PGD for the analysis formulation on this deconvolution problem. We choose two subjects from the UK Bio Bank (UKBB) dataset [START_REF] Sudlow | UK Biobank: An Open Access Resource for Identifying the Causes of a Wide Range of Complex Diseases of Middle and Old Age[END_REF], perform the usual fMRI processing and reduce the dimension of the problem to retain only 8000 time-series of 250 time-frames, corresponding to a record of 3 minute 03 seconds. The full preprocessing pipeline is described in Appendix B. We train the LPGD taut-string network solver on the first subject and Figure 5 reports the performance of the two algorithms on the second subject for λ = 0.1. The performance is reported relatively to the number of iteration as the computational complexity of each iteration or layer for both methods is equivalent. It is clear that LPGD-Taut converges faster than the Accelerated PGD even on real data. In particular, acceleration is higher when the regularization parameter λ is smaller. As mentioned previously, this acceleration is likely to be caused by the better learning capacity of the network in a low sparsity context. The same experiment is repeated for λ = 0.8 in

E P x (u (t) ) -P x (u * )
Accelerated PGD -analysis LPGD-Taut Figure 5: Performance comparison (λ = 0.1) between our analytic prox-TV derivative method and the PGD in the analysis formulation for the HRF deconvolution problem with fMRI data. Our proposed method outperform the FISTA algorithm in the analysis formulation.

Conclusion

This paper studies the optimization of TV-regularized problems via learned PGD. We demonstrated, both analytically and numerically, that it is better to address these problems in their original analysis formulation rather than resort to the simpler (alas slower) synthesis version. We then proposed two different algorithms that allow for the efficient computation and derivation of the required prox-TV, exactly or approximately. Our experiments on synthetic and real data demonstrate that our learned networks for prox-TV provide a significant advantage in convergence speed.

Finally, we believe that the principles presented in this paper could be generalized and deployed in other optimization problems, involving not just the TV-norm but more general analysis-type priors. In particular, this paper only apply for 1D TV problems because the equivalence between Lasso and TV is not exact in higher dimension. In this case, we believe exploiting a dual formulation [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF] for the problem could allow us to derive similar learnable algorithms.

Broader Impact

This work attempts to shed some understanding into empirical phenomena in signal processing -in our case, piecewise constant approximations. As such, it is our hope that this work encourages fellow researchers to invest in the study and development of principled machine learning tools. Besides these, we do not foresee any other immediate societal consequences.

A Network training process strategy

Here, we give a more detailed description of the training procedure used in Section 4.

Optimization algorithm for training

In our experiments, all networks are trained using Gradient Descent (GD) with back-tracking line search. The gradients are computed using automatic differentiation in Pytorch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] for most layers and the weak Jacobian proposed in Subsection 3.1 for the back-propagation through the prox-TV. The learning is stopped once a step-size of η limit = 10 -20 is reached in the back-tracking step. For LPGD-LISTA, the weights of the inner LISTA computing the prox-TV are trained jointly with the parameters of the outer unrolled PGD.

Weight initialization All layers for an unrolled algorithm are initialized using the values of weights in ( 9) that ensure the output of the layer with T layers corresponds to the output of T iterations of the original algorithm. To further stabilize the training, we use a layer-wise approach. When training a network with T 1 + T 2 layers after having trained a network with T 1 layers, the first T 1 layers in the new network are initialized with the weights of the one trained previously, and the remaining layers are initialized using weights value from ( 9). This ensures that the initial value of the loss for the new network is smaller than the one from the shallower one if the unrolled algorithm is monotonous (as it is the case for PGD).

B Real fMRI data acquisition parameters and preprocessing strategy

In this section, we complete the description of the resting-state fMRI (rs-fMRI) data used for the experiment of Fig. 5. For this experiment, we investigate the 6 min long rs-fMRI acquisition (TR=0.735 s) from the UK Bio Bank dataset [START_REF] Sudlow | UK Biobank: An Open Access Resource for Identifying the Causes of a Wide Range of Complex Diseases of Middle and Old Age[END_REF]. The following pre-processing steps were applied on the images: motion correction, grand-mean intensity normalisation, high-pass temporal filtering, Echo planar imaging unwarping, Gradient Distortion Correction unwarping and structured artefacts removal by Independant Components Analysis. More details on the processing pipeline can found in [START_REF] Alfaro-Almagro | Image Processing and Quality Control for the first 10,000 Brain Imaging Datasets from UK Biobank[END_REF].

On top of this preprocessing, we perform a standard fMRI preprocessing proposed in the python package Nilearn3 . This standard pipeline includes to detrend the data, standardize it and filter high and low frequencies to reduce the presence of noise.

C Real fMRI data experiment addition results

Here, we provide an extra experiment for Subsection 4.2 with λ = 0.8λ max and recall the previous one with λ = 0.1λ max to help performance comparison in different regularization regime.

E P x (u (t) ) -P x (u * ) Accelerated PGD -analysis LPGD-Taut (a) λ = 0.1λmax 0 5 10 15 Layers t 10 -8 E P x (u (t) ) -P x (u * ) Accelerated PGD -analysis LPGD-Taut (b) λ = 0.8λmax Figure C
.1: Performance comparison between LPGD-Taut and iterative PGD for the analysis formulation for the HRF deconvolution problem with fMRI data. Our proposed method outperform the FISTA algorithm in the analysis formulation. We notice a slight degradation of the acceleration in this high sparsity context.

We can see that the performance drop in the higher sparsity setting compared to the performances with λ = 0.1λ max but LPGD-Taut still outperforms iterative algorithms for this task on real data.

D Computing λ max for the TV regularized problem

The definition of λ max is the smallest value for the regularisation parameter λ such that the solution of the T V -regularized problem is constant. This corresponds to the definition of λ max in the Lasso, which is the smallest regularisation parameter such that 0 is solution. We here derive its analytic value which is used to rescale all experiments. This is important to define an equiregularisation set for the training and testing samples, to have a coherent and generalizable training. Proposition D.1. The value of λ max for the TV-regularized problem is

λ max = A (Ac1 -x) ∞ where c = p i=1 S i x i p i=1 S 2 i and S i = k j=1 A i,j .
Proof. We first derive the constant solution of the 2 -regression problem associated to (1). For c ∈ R, we consider a constant vector c1. The best constant solution for the 2 -regression problem is obtained by solving

min c∈R f x (c) = 1 2 x -cA1 2 2 . ( 16
)
The first order optimality condition in c reads

∇f x (c) = n i=1 ( k j=1 A i,j )(c k j=1 A i,j -x i ) = n i=1 S i (cS i -x i ) = 0 , (17)
and thus c =

p i=1 S i x i p i=1 S 2 i .
Then, we look at the conditions on λ to ensure that the constant solution c1 is solution of the regularized problem. The first order conditions on the regularized problem reads

0 ∈ ∂P x (c1) = A (Ac1 -x) + λ∂ Dc1 1 (18)
Next, we develop the previous equality:

∀j ∈ {2, . . . k}, A j (Ac1 -x) ∈ λ∂( Dc1 1 ) j = [-λ, λ] since Dc1 = 0 (19)
Thus, the constrains are all satisfied for λ ≥ λ max , with λ max = A (Ac1x) ∞ and as c is solution for the unregularized problem reduced to a constant, c1 is solution of the TV-regularized problem for all λ ≥ λ max .

E Dual formulation

In this work, we devote our effort in the analysis formulation depicted in (1). In this section, we propose to investigate the dual formulation corresponding to (1) in order to rationalize our choice to focus on approaches that solve the prox-TV with an iterative method.

Dual derivation First, we derive the dual of the analysis formulation for the prox-TV. Proposition E.1. [Dual re-parametrization for the analysis formulation TV problem (1)]

Considering the primal analysis problem (with operator and variables defined as previously)

P x (u) = 1 2 x -Au 2 2 + λ Du 1 (20)
Then, the dual formulation reads:

p = -min v 1 2 A † D v 2 2 -v DA † x (21) s.t. v ∞ ≤ λ (22)
Proof. Defining, f and g, such as f (u) = 1 2 x -Au 2 2 and g(u) = λ u 1 and by denoting p the minimum of (20) w.r.t u, the problem reads:

p = min u f (u) + g(Du) (23)
With the Fenchel-Rockafellar duality theorem, we derive the dual re-parametrization:

p = -min v f * (-D v) + g * (v) (24)
Note, that in this case we have the equality with p since the problem (1) is µ-strongly convex, with µ = 1 2 . We have g * (v) =min u g(u)v u. With a component-wise minimization, we obtain g * (v) i = δ |vi|≤λ with δ being the convex indicator. Thus, we deduce that g * (v) = δ v ∞≤λ .

Then, we have f * (v) =min u f (u)v u. By cancelling the gradient we obtain:

f * (v) = 1 2 A † v 2 2 + v A †
x This allows use to conclude the demonstration. Note that, if we set A = Id, we obtain the same problem as [START_REF] Barbero | Modular proximal optimization for multidimensional total-variation regularization[END_REF][START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF].

Performance comparison

We propose to compare the performance of different iterative solvers to assess their performance.

We generate n = 1000 times series to compare the performance between the different algorithms. We set the length of the source signals (u i ) n i=1 ∈ R n×k to k = 40 with a support of |S| = 4 non-zero coefficients. We generate A ∈ R m×k as a Gaussian matrix with m = 40, obtaining then (u i ) n i=1 ∈ R n×p . Moreover, we add Gaussian noise to measurements x i = Au i with a signal to noise ratio (SNR) of 1.0.

We select the PGD and its accelerated version with the synthesis primal formulation (5) ("Synthesis primal A/PGD"), the PGD and its accelerated version with the analysis primal formulation ("Analysis primal A/PGD"). We consider also the PGD and its accelerated version [START_REF] Chambolle | An Algorithm for Total Variation Minimization and Applications[END_REF], for the analysis dual formulation ("Analysis dual A/PGD") and finally we add the primal/dual algorithm [START_REF] Condat | A Primal-Dual Splitting Method for Convex Optimization Involving Lipschitzian, Proximable and Linear Composite Terms[END_REF] for the analysis primal formulation ("Analysis primal dual GD").

Figure E.1 a proposes performance comparison for an exhaustive selection of the algorithm used to solve (1). We see that the analysis primal formulation proposes the best performance for each regularization parameter. We notice that the [START_REF] Condat | A Primal-Dual Splitting Method for Convex Optimization Involving Lipschitzian, Proximable and Linear Composite Terms[END_REF] provides good performance too. Finally, the synthesis primal formulation along with the analysis dual formulation provides the slowest performance. Those results reinforces our choice to focus on the PGD of the analysis primal formulation.

F Proof for Section 2 F.1 Convergence rate of PGD for the synthesis formulation ( 6)

Proof. The convergence rate of ISTA for the synthesis formulation reads Figure E.1: Performance comparison between the iterative solver for the synthesis and analysis formulation with the corresponding primal, dual or primal-dual re-parametrization. We notice that the primal analysis formulation provides the best performance in term of iterations. We also observe that the higher the regularization parameter, the faster the performance for each algorithm.

S(z (t) ) -S(z * ) ≤ ρ 2t z (0) -z * 2 2 . ( 25
We use S(z (t) ) = P (Lz (t) ) = P (u (t) ) to get the correct left-hand side term. For the right hand side, we use z (0) = Du (0) , and z * = Du * , which gives z (0) -z * 2 2 = D(u (0) -u * ) 2 2 ≤ 4 u (0) -u * 2 2 . The last majoration comes from the fact that D 2 2 ≤ 4, as shown per Lemma F.1. This yields Proof. As L is invertible, so is L L. To compute the singular values σ l of L, we will compute the eigenvalues µ l of (L L) -1 and use the relation

P (u (t) ) -P (u * ) ≤ 2 ρ t u (0) -u * 2 2 . (26)
σ l = 1 √ µ l (27)
With simple computations, we obtain

M k = (L L) -1 = L -1 (L ) -1 = D D =       1 -1 0 . . . -1 2 -1 0 . . . . . . . . . . . . 0 -1 2 -1 -1 2       (28)
This matrix is tri-diagonal with a quasi-toepliz structure. Its characteristic polynomial P k (µ) is given by:

P k (µ) = |µ Id -M k | = µ -1 1 0 . . . 1 µ -2 1 0 . . . . . . . . . . . . 0 1 µ -2 1 0 1 µ -2 (29) = (µ -1)Q k-1 (µ) -Q k-2 (µ) (30)
where ( 30) is obtained by developing the determinant relatively to the first line and Q k (µ) is the characteristic polynomial of matrix M k equal to M k except for the the top left coefficient which is replaced by 2

M k =       2 -1 0 . . . -1 2 -1 0 . . . . . . . . . . . . 0 -1 2 -1 0 -1 2       (31)
Using the same development as (30), one can show that Q k verifies the recurrence relationship

Q k (µ) = (µ -2)Q k-1 (µ) -Q k-2 (µ); Q 1 (µ) = 2 -µ, Q 0 (µ) = 1 . ( 32
)
Using this with (30) yields

P k (µ) = Q k (µ) + Q k-1 (µ) (33)
With the change of variable ν = µ-2 2 and denoting

U k (ν) = Q k (2 + 2ν), the recursion becomes U k (ν) = 2νU k-1 (ν) -U k-2 (ν); U 1 (ν) = 2ν, U 0 (µ) = 1 . ( 34
)
This recursion defines the Chebyshev polynomials of the second kind (Chebyshev, 1853) which verifies the following relation

∀θ ∈ [0, 2π], U k (cos(θ)) sin(θ) = sin((k + 1)θ) . (35)
Translating this relationship to

Q k gives ∀θ ∈ [0, 2π], Q k (2 + 2 cos(θ)) sin(θ) = sin((k + 1)θ) . (36)
Using this in (33) shows that for θ ∈ [0, 2π[ P k verify

P k (2 + 2 cos(θ) sin(θ) = sin((k + 1)θ) + sin(kθ) . (37)
The equation sin((k + 1)θ) + sin(kθ) = 0 ,

has l solution in [0, 2π[ that are given by θ l = 2πl 2k+1 for l ∈ {1, . . . n}. As for all l, sin(θ l ) = 0, the values µ l = 2 + 2 cos(θ l ) = 4 cos 2 ( πl 2k+1 ) are the roots of P k and therefor the eigenvalues of M k . Using (27) yields the expected value for σ l .

The singular value of L is thus obtain for l = k and we get

L 2 = σ k = 1 2 cos( πk 2k+1 ) = 1 2 cos( π 2 (1 -1 2k+1 )) , (39)
= 1 2 sin( π 2 1 2k+1 ) = 2k + 1 π + o(1) . (40)
Where the last approximation comes from 1 sin(x) = 1/x + o(1) when x is close to 0.

≥ 2k + 1 4π 2 + o(1)
Proof. Finding the norm of AL can be written as

AL 2 2 = max x∈R k xL A ALx; s.t. x 2 = 1 (41)
From Lemma F.1, we can write L = W ΣV with V , W two unitary matrices and Σ a diagonal matrix with Σ l,l = σ l for all l ∈ {1, .., k}.

First, we consider the case where A A is a rank one matrix with A A = A 2 2 u 1 u 1 , with vector u 1 uniformly sampled from the 2 -ball and fixed A 2 . In this case, as W is unitary, w 1 = W u 1 is also a vector uniformly sampled from the sphere. Also as V is unitary, it is possible to re-parametrize (41) by y = V x such that max

y∈R k A 2 2 y Σw 1 w 1 Σy; s.t. y 2 = 1 (42)
This problem can be maximized by taking y = Σu1 Σu1 2 , which gives

AL 2 2 = A 2 2 Σw 1 2 2
(43) Then, we compute the expectation of Σw 1 2 2 with respect with w 1 , a random vector sampled in the 2 unit ball,

E w1 [ Σw 1 2 2 ] = k l=1 σ 2 l E[u 2 1,i ] = k l=1 1 4 cos 2 πl 2k+1 1 k = 1 2π k l=1 π 2k 1 cos 2 πl 2k+1 . (44)
Here, we made use of the fact that for a random vector u 1 on the sphere in dimension k

, E[u 1,i ] = 1 k
In the last part of the equation, we recognize a Riemann sum for the interval [0, π 2 [. However, x → 1 cos 2 (x) is not integrable on this interval. As the function is positive and monotone, we can still use the integral to highlight the asymptotic behavior of the series. For k large enough, we consider the integral

π 2 -π 2k+1 0 1 cos 2 (x) dx = sin(x) cos(x) π 2 -π 2k+1 0 = cos π 2k+1 sin π 2k+1 = 2k + 1 π + o(1) (45)
Thus, for k large enough, we obtain

E w1 Σw 1 2 2 = 1 2π 2k + 1 π + o(1) (46)
Thus, we get

E AL 2 2 A 2 2 = k + 1 2 π 2 + o(1) (47)
This concludes the case where A A is of rank-1 with uniformly distributed eigenvector.

In the case where A A is larger rank, it is lower bounded by A 2 2 u 1 u 1 where u 1 is its eigenvector associated to its largest eigenvalue, since it is psd. Since A A is a Whishart matrix, its eigenvectors are uniformly distributed on the sphere [START_REF] Silverstein | On the eigenvectors of large dimensional sample covariance matrices[END_REF]. We can thus use the same lower bound as previously for the whole matrix.

G Proof for Section 3 G.1 Proof for Proposition 3.1 Proposition 3.1. [Weak Jacobian of prox-TV] Let x ∈ R k and u = prox µ • T V (x), and denote by S the support of z = Du. Then, the weak Jacobian J x and J µ of the prox-TV relative to x and µ can be computed as J x (x, µ) = L :,S (L :,S L :,S ) -1 L :,S and J µ (x, µ) = -L :,S (L :,S L :,S ) Indeed, when |t| > τ , |z| = |t|τ > 0 and conversely, |z| = 0 when |t| < τ . Moreover, when |t| > τ , we have sign(t) = sign(z) and thus the two expressions for ∂ ST ∂τ match. Using this Lemma G.1, we now give the proof of Proposition 3.1.

Proof. The proof is inspired from the proof from Bertrand et al. (2020, Proposition 1). We denote u(x, µ) = prox µ • T V (x), hence u(x, µ) is defined by u(x, µ) = arg min û 1 2

xû 2 2 + µ û T V (49)

Equivalently, as we have seen previously in (5), using the change of variable û = Lẑ and minimizing over ẑ gives

min ẑ 1 2 x -Lẑ 2 + µ Rẑ 1 . ( 50
)
We denote by z(x, µ) the minimizer of the previous equation. Thus, the solution u(x, µ) of the original problem (49) can be recovered using u(, µ) = Lz(x, µ). Iterative PGD can be used to solve (50) and z(x, mu) is a fixed point of the iterative procedure. That is to say the solution z verifies

z 1 (x, µ) = z 1 (x, µ) -1 ρ (L (Lz(x, µ) -x)) 1 , z i (x, µ) = ST z i (x, µ) -1 ρ (L (Lz(x, µ) -x)) i , µ ρ for i = 2 . . . k . (51)
Using the result from Lemma G.1, we can differentiate (51) and obtain the following equation for the weak Jacobian J x (x, µ) = ∂z ∂x (x, µ) of z(x, µ) relative to x

J x (x, µ) =     1 1 {|z2(x,µ)|>0}
. . .

1 {|z k (x,µ)|>0}     (Id - 1 ρ L L) J x (x, µ) + 1 ρ L Id . ( 52
)
Identifying the non-zero coefficient in the indicator vectors yields J x,S c (x, µ) = 0 J x,S (x, µ) = (Id -1 ρ L :,S L :,S ) J x,S (x, µ) + 1 ρ L :,S .

As, L is invertible, so is L :,S L :,S for any support S and solving the second equations yields the following J x,S = (L :,S L :,S ) -1 L :,S

Using u = Lz and the chain rules yields the expecting result for the weak Jacobian relative to x, noticing that as Ĵx,S c = 0, L Ĵx = L :,S Ĵx,S .

Similarly, concerning, J µ (x, µ), we use the result from Lemma G.1 an differentiale (51) and obtain J µ (x, µ) = ∂z ∂µ (x, µ) of z(x, µ) relative to µ J µ (x, µ) = (

) 55
where γ is the condition number for L i.e.

cos( π 2k+1 ) sin( π 2k+1 ) . We are looking for minimal parameters T and T in such that the error bound in ( 59) is bellow a certain error level δ.

We consider the case where there exists some constants C 0 ≥ u (0)u * 2 and C 1 ≥ max u ( )prox µ ρ (u ( ) ) 2 upper bounding how far the initialization can be compared to the result of the global problem and the sub-problems respectively. We denote α 1 = 3 2 ρ C 1 . The right hand side of ( 59) can be upper bounded by as

ρ 2T   u (0) -u * 2 + 3 T =1 2(1 -γ) Tin u ( -1) -prox µ ρ (u ( -1) ) 2 2 ρ   2 ≤ ρ 2T C 0 + α 1 T (1 -γ) Tin/2 2 (60)
Then, we are looking for T, T in such that this upper bound is lower than δ, i.e.

ρ 2T C 0 + α 1 T (1 -γ) Tin/2 2 ≤ δ (61) ⇔ C 0 + α 1 T (1 -γ) Tin/2 2 - 2δ ρ T ≤ 0 (62) ⇔ C 0 + α 1 T (1 -γ) Tin/2 - 2δ ρ √ T B + α 1 T (1 -γ) Tin/2 + 2δ ρ √ T ≥0 ≤ 0 (63) ⇔C 0 + α 1 T (1 -γ) Tin/2 - 2δ ρ √ T ≤ 0 (64) (65)
Denoting α 2 = 2δ ρ and X = √ T , we get the following function of X and T in f (X, T in ) = α 1 (1γ) Tin/2 X 2α 2 X + C 0 (66)

The inequality f (X, T in ) ≤ 0 has a solution if and only if α 2 2 -4C 0 α 1 (1γ) Tin/2 ≥ 0 i.e.

T in ≥ 2 ln

√ 2ρC1 ln 1 1-γ yields f (X, T in ) = α 2 2 4C 0 X 2 -α 2 X + C 0 = α 2 2 4C 0 (X - 2C 0 α 2 ) 2 for X = 2C0 α2 = √ 2ρC0 √ δ i.e. T = 2ρC 2 0 δ .

Figure 2 :

 2 Figure 2: Algorithm Unrolling -Neural network representation of iterative algorithms. The parameters Θ (t) = {W

Figure 3 :

 3 Figure 3: Performance comparison for different regularisation levels (left) λ = 0.1, (right) λ = 0.8.We see that synthesis formulations are outperformed by the analysis counter part. Both our methods are able to accelerate the resolution of (20), at least in the first iterations.

Figure 4 :

 4 Figure4: Proximal operator error comparison for different regularisation levels (left) λ = 0.1, (right) λ = 0.8. We see that learn the trained unrolled prox-TV barely improve the performance. More interestingly, in a high sparsity context, after a certain point, the error sharply increase.

Figure

F. 2

 2 Computing the spectrum of L Lemma F.1. [Singular values of L] The singular values of L ∈ R k×k are given by

  sign(z 2 (x, µ))1 {|z2(x,µ)|>0} . . . sign(z k (x, µ))1 {|z k (x,µ)|>0}

  Taking the minimal value for T in i.e. T in = 2

  -1 sign(Du) S (Weak derivative of the soft-thresholding; Deledalle et al. 2014). The soft-thresholding operator ST : R × R + → R defined by ST(t, τ ) = sign(t)(|t|τ ) + is weakly differentiable with weak derivatives A very important remark here is to notice that if one denote z = ST(t, τ ), one can rewrite these weak derivatives as ∂ ST ∂t (t, τ ) = 1 {|z|>0} , and ∂ ST ∂τ (t, τ ) =sign(z) • 1 {|z|>0} .

	First, we recall Lemma G.1 to weakly derive the soft-thresholding. Lemma G.1 (48)

∂ ST ∂t (t, τ ) = 1 {|t|>τ } , and ∂ ST ∂τ (t, τ ) =sign(t) • 1 {|t|>τ } ,

where

1 {|t|>τ } = 1, if |t| > τ, 0, otherwise.

Available at https://github.com/albarji/proxTV

Available at https://github.com/hcherkaoui/carpet.

https://nilearn.github.io

Acknowledgement

We gratefully acknowledge discussions with Pierre Ablin, whose suggestions helped us completing some parts of the proofs. H. Cherkaoui is supported by a CEA PhD scholarship. J. Sulam is partially supported by NSF Grant 2007649.

Identifying the non-zero coefficient in the indicator vectors yields J µ,S c (x, µ) = 0 J µ,S (x, µ) = J µ,S c (x, µ) -1 ρ L :,S L :,S J µ,S c (x, µ) -1 ρ sign(z S (x, µ)) .

(56)

As previous, solving the second equation yields the following

Using u = Lz and the chain rules yields the expecting result for the weak Jacobian relative to µ, noticing that as Ĵµ,S c = 0.

G.2 Convergence of the weak Jacobian

Proposition G.2. Linear convergence of the weak Jacobian We consider the mapping z (Tin) : , µR k × R + → R k defined where z (Tin) (x) is defined by recursion

Then the weak

∂µ of this mapping relative to the inputs x and µ converges linearly toward the weak Jacobian J x and J µ of prox µ • T V (x) defined in Proposition 3.1. This mapping defined in (58) corresponds to the inner network in LPGD-LISTA when the weights of the network have not been learned.

Proof. As L is invertible, problem (50) is strongly convex and have a unique solution. We can thus apply the result from Bertrand et al. (2020, Proposition 2) which shows the linear convergence of the weak Jacobian Ĵx = ∂z (T i n) ∂x and Ĵµ = ∂z (T i n) ∂µ for ISTA toward Ĵx and Ĵµ of the synthesis formulation of the prox. Using the linear relationship between the analysis and the synthesis formulations yields the expected result.

G.3 Estimating T in and T to achieve δ error

Using inexact proximal gradient descent results from [START_REF] Schmidt | Convergence Rates of Inexact Proximal-Gradient Methods for Convex Optimization[END_REF] and [START_REF] Machart | Optimal Computational Trade-Off of Inexact Proximal Methods[END_REF], we compute the scaling of T in and T to achieve a given error level δ > 0. Proposition G.3. [Scaling of T and T in w.r.t the error level δ] Let δ the error defined such as P x (u (T ) ) -P x (u * ) ≤ δ. We suppose there exists some constants C 0 ≥ u (0)u * 2 and C 1 ≥ max u ( )prox µ ρ (u ( ) ) 2 Then, T the number of layers for the global network and T in the inner number of layers for the prox-TV scale are given by

Proof. As discussed by [START_REF] Machart | Optimal Computational Trade-Off of Inexact Proximal Methods[END_REF], the global convergence rate of inexact PGD with T in inner iteration is given by