Cesare Molinari

Mathurin Massias

Lorenzo Rosasco

Silvia Villa

Iterative regularization for convex regularizers

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Machine learning often reduces to estimating some model parameters. This approach raises at least two orders of questions: first, multiple solutions may exist, amongst which a specific one must be selected; second, potential instabilities with respect to noise and sampling must be controlled. A classical way to achieve both goals is to consider explicitly penalized or constrained objective functions. In machine learning, this leads to regularized empirical risk minimization [START_REF] Shalev-Shwartz | Understanding machine learning: From theory to algorithms[END_REF]. A more recent approach is based on directly exploiting an iterative optimization procedure for an unconstrained/unpenalized problem. This approach is shared by several related ideas. One is implicit regularization [START_REF] Mahoney | Approximate computation and implicit regularization for very large-scale data analysis[END_REF][START_REF] Gunasekar | Implicit regularization in matrix factorization[END_REF], stemming from the observation that the bias is controlled increasing the number of iterations, just like in penalized methods it is controlled decreasing the penalty parameter. Another one is early stopping [START_REF] Yao | On early stopping in gradient descent learning[END_REF][START_REF] Raskutti | Early stopping and non-parametric regression: An optimal data-dependent stopping rule[END_REF], putting emphasis on the fact that running the iterates to convergence might lead to instabilities in the presence of noise. Yet another, and more classical, idea is iterative regularization, where both aspects (convergence and stability) are considered to be relevant [START_REF] Engl | Regularization of inverse problems[END_REF][START_REF] Kaltenbacher | Iterative regularization methods for nonlinear ill-posed problems[END_REF]. This approach naturally blends modeling and numerical aspects, often improving computational efficiency, while retaining good prediction accuracy [START_REF] Yao | On early stopping in gradient descent learning[END_REF]. Another reason of interest is that iterative regularization may be one of the mechanisms explaining generalization in deep learning [START_REF] Neyshabur | Geometry of optimization and implicit regularization in deep learning[END_REF][START_REF] Gunasekar | Implicit regularization in matrix factorization[END_REF][START_REF] Arora | Implicit regularization in deep matrix factorization[END_REF][START_REF] Vaškevičius | The statistical complexity of early stopped mirror descent[END_REF].

A classic illustrative example is gradient descent for linear least squares. The latter, if suitably initialized, converges (is biased) to the minimal Euclidean norm solution. Moreover, its stability is controlled along the iterative process, allowing to derive early stopping criterions depending on the noise [START_REF] Engl | Regularization of inverse problems[END_REF][START_REF] Raskutti | Early stopping and non-parametric regression: An optimal data-dependent stopping rule[END_REF]. There are a number of developments of these basic results. For example, one line of work has considered extensions to other gradient-based methods, such as stochastic and accelerated gradient descent [START_REF] Zhang | Boosting with early stopping: Convergence and consistency[END_REF][START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF][START_REF] Rosasco | Learning with incremental iterative regularization[END_REF][START_REF] Pagliana | Implicit regularization of accelerated methods in Hilbert spaces[END_REF]. Another line of work has considered classification problems [START_REF] Gunasekar | Implicit regularization in matrix factorization[END_REF][START_REF] Soudry | The implicit bias of gradient descent on separable data[END_REF] and also nonlinear models, such as deep networks [START_REF] Neyshabur | Geometry of optimization and implicit regularization in deep learning[END_REF], see also [START_REF] Kaltenbacher | Iterative regularization methods for nonlinear ill-posed problems[END_REF] for results for non linear inverse problems.

In this work, we are interested in iterative regularization procedures where the considered bias is not the Euclidean norm but rather a general convex functional. The question is to determine whether or not there exists an iteration analogous to gradient descent for such general bias. This question has been stud-ied when the bias is strongly convex. In this case, linearized Bregman iterations (a.k.a. mirror descent) can be used [START_REF] Burger | Error estimation for Bregman iterations and inverse scale space methods in image restoration[END_REF][START_REF] Gunasekar | Characterizing implicit bias in terms of optimization geometry[END_REF]. For this approach accelerated algorithms [START_REF] Matet | Don't relax: early stopping for convex regularization[END_REF] have also been considered and studied. Finally, approaches have also been studied based on diagonal methods [START_REF] Garrigos | Iterative regularization via dual diagonal descent[END_REF] and their acceleration [START_REF] Calatroni | Accelerated iterative regularization via dual diagonal descent[END_REF]. The general convex case, even for linear models, is much less understood. There have been studies for ADMM/Bregman iteration, but the procedure requires solving a nontrivial optimization problem at each iteration [START_REF] Burger | Error estimation for Bregman iterations and inverse scale space methods in image restoration[END_REF]. Further, stability and convergence results are proved only in terms of Bregman divergence, which in general is a weak result. Interestingly, various recent results study iterative regularization for sparse recovery, where the bias is defined by an 1 norm [START_REF] Agarwal | Stochastic optimization and sparse statistical recovery: Optimal algorithms for high dimensions[END_REF][START_REF] Osher | Sparse recovery via differential inclusions[END_REF][START_REF] Vaškevičius | Implicit regularization for optimal sparse recovery[END_REF].

In this paper we propose and study an efficient algorithm for general convex bias beyond 1 norm. Indeed, we adapt the Chambolle and Pock (CP) algorithm, popular in imaging [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], and study its iterative regularization properties. The CP algorithm is a first order primal-dual method, thus easy to implement and requiring only matrix-vector multiplications and proximity operators. In the setting of linear models with worst case errors, our analysis provides dimension free convergence and stability results in terms of both Bregman divergence and approximate feasibility. A combination of these two results allows to derive strong convergence results in the 1 norm case. The proof relies on results from the analysis of primal-dual methods with errors [START_REF] Rasch | Inexact first-order primal-dual algorithms[END_REF]. From our general results, several special cases can be derived and we discuss as an example sparse recovery, proving dimension free estimates in norm. In the experimental section, we investigate the proposed method and show state-of-theart performances with significant computation savings compared to the Tikhonov approach.

Notation, The set of integers from 1 to n is rns. Let f : R n Ñ RYt`8u and J : R p Ñ RYt`8u be proper, convex, and lower semicontinuous. The subdifferential of J at w P R p is BJpwq. The Bregman divergence associated to J is denoted D θ J pw, w 1 q :" Jpwq Jpw 1 q ´xθ, w ´w1 y, where θ P BJpw 1 q. The Fenchel-Legendre conjugate of f is f ‹ pθq :" sup w xw, θy´f pwq. The indicator function ι tyu is equal to zero if the argument equals y and `8 otherwise. For a proper convex lower semicontinuous function J, prox J pxq " arg min y Jpyq ` x ´y 2 {2.

2 Over-parametrization, implicit and explicit regularization

The basic problem of supervised learning is to find a relationship to predict outputs y from inputs x,

x Þ Ñ f pxq « y ,
given a limited number of pairs px i , y i q n i"1 with, e.g. x i P R d and y i P R. The search for a solution is typically restricted to a set of parametrized functions f w , with w P R p . A prototype example are linear models where p " d and f w pxq " xw, xy, or more generally f w pxq " ř p j"1 w j φ j pxq, for some dictionary φ j : R d Ñ R, j " 1, . . . , p [START_REF] Hastie | The Elements of Statistical Learning[END_REF][START_REF] Shalev-Shwartz | Understanding machine learning: From theory to algorithms[END_REF]. In modern applications, it is often the case that the number of parameters p is vastly larger than the number of available data points n, a regime called over-parametrized. Excluding degenerate cases, one can then expect to find a solution w capable of interpolating the data, that is satisfying,

f w px i q " y i , @i P rns . (1)
In the sequel we consider the case of a linear f w pxq " xw, xy. A popular method to find a solution to (1) is gradient descent on least squares, also called Landweber iteration:

w k " w k´1 ´γX J pXw k´1 ´yq , (2)
where X and y are the data matrix and the outputs vector, respectively (see Section 3 for more details). It is well known [START_REF] Engl | Regularization of inverse problems[END_REF] that, if initialized at w 0 " 0 and with γ ă 2{ X 2 op , the iterations of gradient descent converge to a specific solution, namely arg min

wPR p w s.t. Xw " y . (3)
This means that amongst all solutions, the algorithm is implicitly biased towards that with small norm. The bias is implicit in the sense that there is no explicit penalization or constraint in the iterations (2). This approach can be contrasted to explicit penalization (Tikhonov regularization), w pλq " arg min

wPR p λ w 2 ` y ´Xw 2 , ( 4
)
where the minimal norm solution (3) is obtained for λ going to zero. It is well known that for Tikhonov regularization larger values of λ improve stability. Interestingly, the same effect can also be achieved with gradient descent (2), by not running the iterations until convergence, a technique often referred to as early stopping [START_REF] Engl | Regularization of inverse problems[END_REF][START_REF] Yao | On early stopping in gradient descent learning[END_REF]. In this view the number of iterations k plays the role of a regularization parameter just like λ in Tikhonov regularization (or rather 1{λ). Iterative regularization is particularly appealing in the large scale setting, where substantial computational savings are expected: early stopping needs a finite number of iterations (2), while Tikhonov regularization requires solving exactly Problem (4) for multiple values of λ.

It is natural to ask whether the above iterative regularization scheme applies to biases beyond the Euclidean norm. For a strongly convex J, an answer is given by considering the mirror descent algorithm [START_REF] Nemirovsky | Problem complexity and method efficiency in optimization[END_REF][START_REF] Beck | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF] with respect to the Bregman divergence induced by J. The bias J is not used to define an explicit penalization of an empirical risk, but it appears in the mirror descent algorithm, and in this sense is "less implicit". The results in [START_REF] Benning | Gradient descent in a generalised Bregman distance framework[END_REF] and [START_REF] Gunasekar | Characterizing implicit bias in terms of optimization geometry[END_REF] show that mirror descent is implicitly biased towards the solution of the following problem

min wPR p Jpwq s.t. Xw " y , (5)
and exhibit similar regularization and stability properties to the one of the gradient descent algorithm. In both [START_REF] Benning | Gradient descent in a generalised Bregman distance framework[END_REF] and [START_REF] Gunasekar | Characterizing implicit bias in terms of optimization geometry[END_REF], the key technical assumption is strong convexity of J leaving open the question of dealing with biases that are only convex. In this paper, we take steps to fill in this gap studying an efficient approach for which we characterize the iterative regularization properties.

Problem setting and proposed algorithm

We begin describing the algorithm we consider and its derivation. We first set some notation and the main assumption. In the following, X is an n by p matrix and y an n-dimensional vector. We consider the case where y is unknown, and a vector y δ is available such that y ´yδ ď δ, where δ ě 0 can be interpreted as the noise level.

Assumption 1. We assume that the bias J : R p Ñ R Y t`8u is proper, convex, and lower semicontinuous. We also assume that Problem (5) has at least one solution (in particular the linear equation has at least one solution for the exact data y).

The latter is satisfied, for instance, if J is coercive, n ď p and X has rank n. Then Problem (5) is feasible and has a solution for every y. Our assumption is weaker: we do not even require feasability of the noisy constraint Xw " y δ .

Note that we use a vectorial notation for simplicity but our results are dimension free and sharp for an infinite dimensional setting where X is a linear bounded operator between separable Hilbert spaces.

Proposed algorithm

Consider the following iterations, with initialization w 0 P R p , θ 0 " θ ´1 P R n , and parameters τ , σ such that στ X 2 op ă 1:

# w k`1 " prox τ J pw k ´τ X J p2θ k ´θk´1 qq , θ k`1 " θ k `σpXw k`1 ´yq . ( 6
)
If θ 0 " 0, since θ k " σ ř k i"1 pXw i ´yq, this algorithm can be rewritten without θ k :

w k`1 " prox τ J ´wk ´τ σX J `řk i"1 pXw i ´yq`Xw k ´y˘¯.
In terms of computations the algorithm ( 6) is very similar to the forward-backward/proximal gradient algorithm [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF]. The difference is that the gradient term is here replaced by the sum of past gradients. We instantiate algorithm (6) for two popular choices of J. For J " ¨ 2 , the updates read:

# w k`1 " 1 1`τ pw k ´τ X J p2θ k ´θk´1 qq , θ k`1 " θ k `σpXw k`1 ´yq .
Notice that, though involving very similar computations, the algorithm does not reduce to gradient descent iterations (2).

For J " ¨ 1 , denoting by STp¨, τ q the softthresholding operator of parameter τ , the iterations (6) read:

# w k`1 " STpw k ´τ X J p2θ k ´θk´1 q, τ q , θ k`1 " θ k `σpXw k`1 ´yq .
Also in this case, it is similar -yet not equivalent -to a popular algorithm to solve the Tikhonov problem: the Iterative Soft-Thresholding Algorithm [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF].

Proposition 1. Under Assumption 1, the iterations (6) converge to a point pw ‹ , θ ‹ q such that Xw ‹ " y.

Additionally, w ‹ is a minimizer of J amongst all interpolating solutions, meaning that it solves Problem (5).

Notice that for Landweber iteration, converging to the minimal norm solution (3) requires w 0 " 0. Our procedure has no such requirement: for any initialization, the iterates w k converge to a solution of Problem ( 5).

We now provide some examples of the above setting.

Example 2 (Sparse recovery). Choosing J " ¨ 1 corresponds to finding the minimal 1 -norm solution to a linear system, and in this case Problem (5) is known as Basis Pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. The relaxed approach of Problem (4) in this case yields the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]. 1 -based approaches have had a tremendous impact in imaging, signal processing and machine learning in the last decades [START_REF] Hastie | Statistical Learning with Sparsity: The Lasso and Generalizations[END_REF].

Example 3 (Low rank matrix completion). In several applications, such as recommendation systems, it is useful to recover a low rank matrix, starting from the observation of a subset of its entries [START_REF] Candès | Exact matrix completion via convex optimization[END_REF]. A convex formulation is:

min W PR p 1 ˆp2 W ˚s.t. W ij " Y ij @pi, jq P D , ( 7
)
where ¨ ˚is the nuclear norm and D Ă rp 1 s ˆrp 2 s is the set of observed entries of the matrix Y . In that case, X is not a design matrix: it is a (self-adjoint and linear) masking operator from R p1ˆp2 to R p1ˆp2 , such that pXW q ij has value W ij if pi, jq P D and 0 otherwise; the constraints write XW " XY .

Example 4 (Total variation). In many imaging tasks such as deblurring and denoising, regularization through total variation allows to simultaneously preserve edges whilst removing noise in flat regions [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. The problem of Total Variation is min W PR p 1 ˆp2 ∇W 1 s.t. XW " Y , where X is usually a blurring operator. The problem can be reformulated as:

min W Ωp W q s.t. X W " Ỹ , with W " ˆW U ˙, Ωp W q " }U } 1 , X " ˆX 0 ∇ ´Id ˙and Ỹ " ˆY 0 ˙.
Other examples include 8 norm for antiparsimony, 1 { 2 mixed-norm for group sparsity, or Kullback-Leibler divergence. Finally, for any linear operator K and convex function J such that prox J is available, we can handle problems of the form min w JpKwq s.t.

Xw " y .

Indeed, we can introduce an auxiliary variable as in Example 4, and avoid the computation of prox J˝K .

Chambolle-Pock algorithm

In this section we prove Proposition 1 by casting ( 6) as an instance of the Chambolle-Pock algorithm (Chambolle and Pock, 2011) which solves:

min w f pXwq `gpwq .
Hence, for f " ι tyu and g " J, it can minimize a convex function on a set defined by linear equalities, as in Problem (5). For this choice of f and g, it instantiates as (6) (see Appendix B.2).

Amongst other assets, algorithm (6) only involves matrix-vector multiplications, and the computation of prox J , available in closed-form in many cases. The only tunable parameters are two step-sizes, τ and σ, which are easy to set. As usual for this class of methods, called primal-dual, the Lagrangian is a useful tool to establish convergence results. The Lagrangian of Problem ( 5) is

Lpw, θq " Jpwq `xθ, Xw ´yy , (8)
where θ P R n is the dual variable. Under a technical condition (Appendix B.1), w ‹ is a solution of Problem (5) if and only if there exists a dual variable θ ‹ such that pw ‹ , θ ‹ q is a saddle-point for the Lagrangian, namely, iff for every pw, θq P R p ˆRn ,

Lpw ‹ , θq ď Lpw ‹ , θ ‹ q ď Lpw, θ ‹ q . ( 9
)
The variable θ is in our setting just an auxiliary variable, and we will be interested in convergence properties of w k towards w ‹ .

Other algorithms As mentioned in the introduction, other algorithms could be considered, e.g. ADMM/Bregman iteration. However, we are not aware of methods that can be efficiently implemented in our general setting. In Appendix A. we provide an extensive review discussing the connection with a number of different approaches and related works.

Theoretical analysis

In this section, we analyze the convergence properties of Algorithm ( 6). First, we need to choose a suitable criterion to estimate the approximation properties of the iterates. In general, it is not reasonable to expect a rate of convergence for the distance between the iterates and the solution. Indeed, since the problem is only convex, it is well known that the convergence in distance can be arbitrarily slow. In Section 4.1, we explain why a reasonable choice is given by the duality gap together with the residual norm (respectively, Lpw k , θ ‹ q ´Lpw ‹ , θ k q and Xw k ´y ). For these two quantities, we derive:

• convergence rates in the exact case, i.e. when the data y is available (Proposition 6);

• early-stopping bounds in the inexact case, i.e. when the accessible data is only y δ with y δ ´y ď δ (Proposition 7 and Corollary 8).

In Section 4.4, we apply our analysis to the specific choice of J equal to the 1 -norm. In this particular case, combining the previous results, we even obtain bounds directly on w k ´w‹ .

Measure of optimality

To discuss which criterion is significant to study the algorithm convergence, we recall from (9) that, if

Lpw 1 , θq ´Lpw, θ 1 q ď 0 (10) for every pw, θq P R p ˆRn , then pw 1 , θ 1 q is a primaldual solution. In general, it is difficult to prove that Equation ( 10) holds for every pw, θq P R p ˆRn . Then, given a saddle-point pw ‹ , θ ‹ q and a generic pw 1 , θ 1 q P R p ˆRn , it is popular to consider the quantity

Lpw 1 , θ ‹ q ´Lpw ‹ , θ 1 q ě 0 . ( 11
)
To establish the optimality of pw 1 , θ 1 q, it is not enough to ensure Lpw 1 , θ ‹ q ´Lpw ‹ , θ 1 q " 0. Lemma 5, proved in Appendix C.1, shows that this condition, when coupled with Xw 1 " y, implies that w 1 is a solution of Pb (5).

Lemma 5. Let pw ‹ , θ ‹ q be a primal-dual solution and pw 1 , θ 1 q a point in R p ˆRn such that Lpw 1 , θ ‹ q Ĺpw ‹ , θ 1 q " 0 and Xw 1 " y. Then pw 1 , θ ‹ q is a primaldual solution.

Thus, the quantities Lpw k , θ ‹ q ´Lpw ‹ , θ k q and Xw k ´y , studied together, are a reasonable measure of optimality for the iterate w k .

Note that Lpw k , θ ‹ q ´Lpw ‹ , θ k q is the error measure used in a series of papers dealing with regularization of inverse problems with general convex regularizers, see e.g. [START_REF] Burger | Error estimation for Bregman iterations and inverse scale space methods in image restoration[END_REF]. It is well known that if J is strongly convex then this quantity controls the distance in norm (Remark 9) and therefore is a proper measure of convergence. If J is only convex, this measure of error can be quite weak. In Section 4.4 we point out the limitations of this quantity when dealing with J " ¨ 1 . For this choice of J, Lp0, θ ‹ q ´Lpw ‹ , θ k q is 0 for any θ k ; as shown on Figure 1, this quantity is also 0 when w k and w ‹ have the same support and sign.

Exact case

First consider the iterates pw k , θ k q obtained by applying iterations (6) to the exact problem, namely where the data y is available. Let pw ‹ , θ ‹ q be a saddle-point for the Lagrangian. Denoting the primal-dual variables by z " pw, θq, we have z k " pw k , θ k q for the iterates of the algorithm and z ‹ " pw ‹ , θ ‹ q for the saddle-point. For τ and σ ą 0, define V as the following square weighted norm on R p ˆRn :

V pzq :" 1 2τ w 2 `1 2σ θ 2 . ( 12
)
For the averaged iterates w k :" 1 k ř k t"1 w t and θ k :"

1 k ř k t"1 θ t ,
we have the following rates.

Proposition 6 (Convergence rates). Under Assumption 1, let ε P p0, 1q and assume that the step-sizes are such that στ ď ε{ X 2 op . Then

Lpw k , θ ‹ q ´Lpw ‹ , θ k q ď V pz 0 ´z‹ q k and Xw k ´y 2 ď 2p1 `εqV pz 0 ´z‹ q σεp1 ´εqk .
The first result is classical (see [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]). Alternatively, it can be obtained by setting δ " 0 in Proposition 7, where we study the more general inexact case. To the best of our knowledge, the second bound is new and can also be derived by setting δ " 0 in Proposition 7. A similar result, in the more specific case of primal-dual coordinate descent, can be found in [START_REF] Fercoq | A coordinate descent primal-dual algorithm with large step size and possibly non separable functions[END_REF]. Note that both results of Proposition 6 are true for every primal-dual solution. On the other hand, the left-hand-side in the second equation does not depend on the selection of z ‹ and so the bound can be improved by taking the inf over all primal-dual solutions.

Inexact case

We now consider the iterates pw k , θ k q, and their averaged versions pw k , θ k q, obtained by applying iterations (6) to the noisy problem, where y is replaced by y δ with y δ ´y ď δ. In Proposition 7, we derive early-stopping bounds for the iterates, in terms of duality gap Lpw k , θ ‹ q ´Lpw ‹ , θ k q and residual norm Xw k ´y . We highlight that, despite the error in the data y δ , both quantities are defined in terms of y and hence related to the noiseless problem. In particular, pw ‹ , θ ‹ q is a saddle-point for the noiseless Lagrangian. We have the following estimates, whose proofs are given in Appendix C.3.

Proposition 7 (Stability). Under Assumption 1, let ε P p0, 1q and assume that the step-sizes are such that στ ď ε{ X 2 op . Then, and

Lpw k , θ ‹ q ´Lpw ‹ , θ k q ď 1 k ´aV pz 0 ´z‹ q `?2σδk ¯2 (
Xw k ´y 2 ď 2p1 `εq σεp1 ´εq « a 2σV pz 0 ´z‹ qδ `σε 1 ´ε δ 2 `2σδ 2 k `1 k V pz 0 ´z‹ q ff . ( 14
)
Note that, in the exact case δ " 0, we recover the convergence results stated in Proposition 6. Moreover, we have the following corollary.

Corollary 8 (Early-stopping). Under the assumptions of Proposition 7, choose k " c{δ for some c ą 0.

Then there exist constants C, C 1 and C 2 such that Lpw k , θ ‹ q ´Lpw ‹ , θ k q ď Cδ and

Xw k ´y 2 ď C 1 δ `C2 δ 2 .
The constants appearing in the Corollary are the ones from Proposition 7. They only depend on the saddlepoint z ‹ , the initialization z 0 and the step-sizes τ, σ.

We next add some remarks.

Remark 9. When J is γ-strongly convex, in particular when Jp¨q " 1 2 ¨ 2 , both the residual norm and the distance between the averaged iterate and the solution can be controlled by Lpw k , θ ‹ q´Lpw ‹ , θ k q. Indeed, recalling Equation (25),

Xw k ´y 2 ď X 2 w k ´w‹ 2 ď 2 X 2 γ D ´XJ θ ‹ J `wk , w ‹ " 2 X 2 γ " Lpw k , θ ‹ q ´Lpw ‹ , θ k q ı .
In particular, the previous early-stopping bounds are of the same order of the ones obtained by dual gradient descent in [START_REF] Matet | Don't relax: early stopping for convex regularization[END_REF].

Remark 10. Similar estimates have been obtained in [START_REF] Burger | Error estimation for Bregman iterations and inverse scale space methods in image restoration[END_REF], both for the Tikhonov variational scheme and for the Bregman iteration (also called inverse scale space method) with stopping-criteria given by the discrepancy principle. In the first case (see Theorem 3.1), for a suitable choice of the regularization parameter, the authors get similar estimates for the Tikhonov regularized solution w λ : D s J pw λ , w ‹ q ď Cδ and Xw λ ´y 2 ď C 1 δ 2 , where D s J is the symmetric Bregman divergence. For the Bregman iteration (see Theorem 4.2), they get an early-stopping bound on D p k J pw ‹ , w k q, where p k P BJpw k q. Note that they do not get any estimate for the quantity D ´XJ θ ‹ J `wk , w ‹ neither for the residual norm. Moreover, the method requires to solve, at each iteration, an optimization problem with the same complexity of the original one.

Proof Sketch The proof of Proposition 7 is inspired by [START_REF] Rasch | Inexact first-order primal-dual algorithms[END_REF]. In this paper, the kind of errors allowed in the prox of the nonextrapolated step (θ k update) are more general than the ones allowed for the extrapolated step (w k update).

Here, we study stability properties of algorithm (6) when y is replaced by y δ . This change can be read as an inexact proximity operator in the update of θ computation; in order to have this error in the nonextrapolated step, we study algorithm (6), that is CP algorithm applied to the dual problem. We summarize here the main steps. In Lemma 15, we derive a "descent property" for every step t, which we then cumulate summing from t " 1 to t " k and using two different approximations (Lemmas 16 and 17). The two bounds that we get are similar, but independent. The first one has the following form,

1 2σ θ k ´θ‹ 2 `k ÿ t"1 rLpw t , θ ‹ q ´Lpw ‹ , θ t qs ď V pz 0 ´z‹ q `δ k ÿ t"1 θ t ´θ‹ . ( 15
)
We use the latter twice. First we combine it with Lemma 14, a discrete version of Bihari's Lemma. This allows to estimate, for every 1 ď t ď k, the quantity θ t ´θ‹ ď 2σδk `a2σV pz 0 ´z‹ q .

(16)

Then we use again Equation ( 15), joint with the previous information, to find a bound on ř k t"1 rLpw t , θ ‹ q ´Lpw ‹ , θ t qs. The second inequality (see Lemma 17) has the following form,

σα 2η k ÿ t"1 Xw t ´y 2 ď V pz 0 ´z‹ q `δ k ÿ t"1 θ t ´θ `1 2 σ pη ´1q δ 2 k . ( 17
)
Using again the bound on θ t ´θ‹ and choosing η " p1 `εq { p1 ´εq, we find an estimate for ř k t"1 Xw t ´y 2 . In both cases, we get the claim on the averaged iterates by Jensen's inequality.

Remark 11. Following Rasch and Chambolle (2020), our analysis can straightforwardly be extended to the case in which prox J in iterations (6) is computed approximately, with type-2 errors. This leads to bounds similar in spirit to those of Proposition 7, but involving the cumulative prox error up to iteration k and thus yielding slightly different stopping criteria.

An example: sparse recovery

In the case of sparse recovery (J " ¨ 1 ), controlling the duality gap and the feasability yields a bound on the distance to the minimizer, thanks to the following result (Grasmair et al., 2011, Lemma 3.10).

Lemma 12. Let pw ‹ , θ ‹ q be such that Xw ‹ " y and ´XJ θ ‹ P B ¨ 1 pw ‹ q. With Γ :" tj P rps : |X J :j θ ‹ | " 1u, assume that X Γ (X restricted to columns whose indices lie in Γ) is injective. Let m :" max jRΓ |X J :j θ ‹ | ă 1. Then, for all w P R p ,

w ´w‹ ď X ´1 Γ op Xw ´y `1` X ´1 Γ op X op 1´m D ´XJ θ ‹ ¨ 1 pw, w ‹ q . ( 18
)
Note that, under the assumptions of Lemma 12, the primal solution to Problem ( 19) is unique (see (Grasmair et al., 2011, Thm 4.7)). Combining the latter with Corollary 8 yields a strong early-stopping result.

Corollary 13 (Early-stopping for J " ¨ 1 ).

Under the assumptions of Proposition 7 and Lemma 12, choose k " c{δ for c ą 0. Then there exist constants C 1 and C 2 such that

w k ´w‹ ď C 1 ? δ `C2 δ .
The constants C 1 , C 2 depend on the saddle-point z ‹ , the initialization z 0 , the step-sizes τ, σ and the norms of X and X ´1 Γ . A completely different approach has There exists a stopping time: these distances reach a minimum before converging to their limit. been considered, for the same problem, in [START_REF] Vaškevičius | Implicit regularization for optimal sparse recovery[END_REF]. A related approach, based on dynamical systems, has been proposed in [START_REF] Osher | Sparse recovery via differential inclusions[END_REF]. Similar results for the Tikhonov regularization approach can be found in Schöpfer and Lorenz (2019b).

Empirical analysis

We stress that there is no implicit regularization result dealing with any non strongly convex J to compare to. The only competitor is therefore the Tikhonov approach. The code with scripts to reproduce the experiments (relying heavily on numpy [START_REF] Harris | Array programming with numpy[END_REF] and numba [START_REF] Lam | Numba: A llvm-based python jit compiler[END_REF]) is available at lcsl.github.io/iterreg.

Sparse recovery

Random data for this experiment are generated as follows: pn, pq " p200, 500q, columns of X are Gaussian with covpX :i , X :j q " 0.2 |i´j| , y " Xw 0 where w 0 has 75 equal non zero entries, scaled such that y " 20 (in order to have a meaningful range of values for δ). Note that the linear system Xw " y δ has solutions for any y δ , since X is full-rank. The noiseless solution w ‹ is determined by running algorithm (6) up to convergence, on y. For the considered values of δ, y δ is created by adding i.i.d. Gaussian noise to y, so that y ´yδ " δ. We denote by w δ k the iterates of algorithm (6) ran on y δ .

Existence of stopping time. In the first experiment, we highlight the existence of an optimal iterate in terms of distance to w ‹ . Figure 2 shows semiconvergence: before converging to their limit, the iterates get close to w ‹ . Note that this is stronger than the results of Corollary 13, since the optimal iteration here is the minimizer of the distance, and not some iterate for which there exists an upper bound on the distance to w ‹ . As expected, as δ decreases, the optimal iteration k increases and the optimal iterate w δ k is closer to w ‹ .

Dependency of empirical stopping time on δ.

In the same setting as above, for 20 values of δ between 0.1 and 6, we generate 100 values of y δ . We run algorithm (6) for 5000 iterations on y δ and determine the empirical best stopping time as k ‹ pδq " arg min k w δ k ´w‹ ă `8. Figure 3 shows the mean of the inverse empirical stopping time as a function of δ, where a clear linear trend (k " c{δ) appears as suggested by Proposition 7 and Corollary 13.

In real settings, w ‹ and δ are unknown, and so is the stopping time. It can still be by cross validation or similar procedure, as is usually done for the optimal λ in explicit regularization.

Comparison with Tikhonov approach on real data. The most popular regularization approach is to solve Problem (4) (here, the Lasso) for typically 1 100 values of λ geometrically chosen as λ t " 10 ´3t{99 X J y 8 for t " 0, . . . , 99. In Figure 4 we compare the Lasso regularization path to the Basis Pursuit optimization path of the Chambolle-Pock algorithm. The dataset for this experiment is rcv1train from libsvm [START_REF] Fan | Liblinear: A library for large linear classification[END_REF], with pn, pq " p20 242, 26 683q. The figure of merit is the prediction mean squared error on left out data, using 4-fold cross validation (dashed color lines), with the average across the folds in black. The horizontal line marks the λ (resp. the iteration k) for which the Lasso path (resp. the optimization path of iterations ( 6)) reaches its minimum MSE on the test fold.

The first observation is that the Basis Pursuit solution (both the end of the optimization (k " `8) and regularization paths (λ " 0)) performs very poorly, having a MSE greater than the one obtained by the 0 solution For the bottom plot, this would also be visible if 1 default grid in scikit-learn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF] the number of iterations of Algorithm ( 6) was picked greater than 500, which we do not do for readability of the figure. It is therefore necessary to early stop.

The second observation is that the minimal MSEs on both paths are similar: 0.19 for Lasso path, 0.21 for optimization path of Algorithm (6). The main point is however that it takes 20 iterations of algorithm (6) to reach its best iterate, while the optimal λ for the Lasso is around λ max {100. If the default grid of 100 values between λ max and λ max {1000 was used, this means that 66 Lassos must be solved, each one needing hundreds or thousands of iterations to converge. This is reflected in the timings: 0.5 s for Algorithm (6) vs 50 s for Tikhonov, eventhough we use a state-of-the-art coordinate descent + working set approach to solve the Lasso, with warm-start (using the solution for λ t´1 as initialization for problem with λ t ).

Low rank matrix completion

Random data for this experiment is generated as follows: the matrices are d ˆd with d " 20. Y is equal to U V J with U and V of size d ˆ5, whose entries are i.i.d Gaussian (Y is rank 5). We scale Y such that Y " 20. Recall that in low rank matrix completion (Example 3), X corresponds to a masking operator (the observed entries); to determine which entries are observed, we uniformly draw d 2 {5 observed couples pi, jq P rds ˆrds. Figure 5 shows the same type of results as Figure 2: iterates first approach the noiseless solution, then get further away, justifying early stopping of the iterates. For this experiment, we use higher values for δ to better highlight the semiconver- gence, as curves get flatter for e.g. δ " 5.7. Note that in that case, the algorithm can still be early stopped to save computations.

Conclusion

We have studied implicit regularization for convex bias, not necessarily strongly convex nor smooth. We proposed to use the Chambolle-Pock algorithm and we analyzed both convergence and stability to deterministic worst case noise. Our general analysis was specialized, as an example, to the problem of sparse recovery. The approach was investigated empirically both for sparse recovery and matrix completion, showing great timing improvements over relaxation approaches. A future development is to consider more specific noise models than the worst-case, such as stochastic noise. We emphasize again that our results hold in infinite dimension. It would be interesting to specialize our analysis in the finite dimensional setting, when the noisy solution always exists (in the least-square sense) and so the iterates produced by the algorithm are bounded. Moreover, it would be interesting to consider additional assumptions such as sparsity. Considering the role of initialization or nonlinear models would also be of interest. Finally, it would complete the analysis to obtain lower bounds for this class of problems, to confirm the sharpness of our results.

A Detailed discussion of related works

The idea of exploiting the implicit regularizing properties of optimization algorithms is not new, and has been studied in three different related areas, often under the name of iterative regularization: inverse problems, image restoration, and machine learning. The related results can roughly be divided in those assuming strong convexity of J and those assuming only convexity of J. Related approaches to implicit regularization include diagonal strategies and exact regularization approaches. Extensions to general data fits and non-convex/nonlinear problems have been considered. In the following we briefly review existing results.

' Gradient and stochastic descent. The study of implicit regularization properties of gradient descent, known in the inverse problem community as Landweber method, goes back to the 50s [START_REF] Engl | Regularization of inverse problems[END_REF]. The classical result shows that gradient descent applied to least squares and initialized at 0 converges to the minimal norm solution of the linear equation (1). Accelerated versions have been also studied under the name of ν-method [START_REF] Engl | Regularization of inverse problems[END_REF]. Generalization towards more general regularizers, apart from p norms with p ą 1, has not been considered much by this community, while there is a rich literature in the non-convex setting for nonlinear inverse problems [START_REF] Kaltenbacher | Iterative regularization methods for nonlinear ill-posed problems[END_REF]. These ideas have been extended to machine learning considering regularizing properties of gradient descent [START_REF] Yao | On early stopping in gradient descent learning[END_REF] and its stochastic versions [START_REF] Moulines | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF][START_REF] Rosasco | Learning with incremental iterative regularization[END_REF].

' Linearized Bregman and Mirror descent. The interest in more general regularizers has been mainly motivated by imaging applications and total variation regularization. Starting from [START_REF] Osher | An iterative regularization method for total variation-based image restoration[END_REF] there is an entire line of work devoted to iterative regularization for general convex regularizers (see e.g. [START_REF] Burger | Error estimation for Bregman iterations and inverse scale space methods in image restoration[END_REF] and references therein). We briefly review the available algorithms and their advantages and limitations. If strong convexity of J is assumed, the algorithm of choice is mirror descent [START_REF] Nemirovski | Problem complexity and method efficiency in optimization[END_REF][START_REF] Teboulle | Mirror descent and nonlinear projected subgradient methods for convex optimization[END_REF]. It has been popularized in the inverse/imaging problems community under the name of linearized Bregman iteration [START_REF] Yin | Bregman iterative algorithms for l1-minimization with applications to compressed sensing[END_REF]. It has been shown that this algorithm in combination with a discrepancy type stopping rule regularizes ill posed problems. The stability and regularization properties of the accelerated variant of the algorithm have been studied also in [START_REF] Matet | Don't relax: early stopping for convex regularization[END_REF], using a different approach, based on the interpretation of the method as a gradient descent applied to the dual problem (21). Similar ideas can be found in (Schöpfer and Lorenz, 2019a).

' Bregman iteration and ADMM. If the regularizer J is not strongly convex, but only convex, as in our case, the algorithm above cannot be applied. The algorithm of choice is in this context ADMM [START_REF] Boyd | Distributed optimization and statistical learning via the alternating direction method of multipliers[END_REF], which has been studied in the imaging community under the name of Bregman iteration. Its regularization properties can be found in [START_REF] Burger | Error estimation for Bregman iterations and inverse scale space methods in image restoration[END_REF]. However, this method has a main drawback: at each iteration the solution of a nontrivial optimization problem of the form min }Xw ´y} 2 `Jpwq `xw, ηy, for η P R p is required, and in general la subroutine is needed at each iteration. In the setting where n is big, this can have a high computational cost. The extension of this approach to nonlinear inverse problems has been considered in [START_REF] Bachmayr | Iterative total variation schemes for nonlinear inverse problems[END_REF].

' Bregmanized Operator Splitting and linearized/preconditioned ADMM. These are variants of Bregman iteration and ADMM very similar to the CP algorithm: they rely on preconditioning to avoid the solution of a difficult optimization problem at each iteration. These have been used empirically as regularizing procedures in the context of inverse and imaging problems [START_REF] Zhang | A unified primaldual algorithm framework based on Bregman iteration[END_REF][START_REF] Zhang | Bregmanized nonlocal regularization for deconvolution and sparse reconstruction[END_REF]. We are not aware of any theoretical quantitative stability result.

' Diagonal approaches The implicit regularization techniques described above are well-suited for problems where the quadratic data fit is appropriate. If other losses are used, this approach completely neglect them. A way to circumvent this problem is to use a diagonal strategy. The idea is to combine an optimization algorithm with a sequence of approximations of the original problem (19) which change at each iteration [START_REF] Bahraoui | Convergence of diagonally stationary sequences in convex optimization[END_REF]. Convergence rates and stability of diagonal approaches for inverse problems have been considered in [START_REF] Garrigos | Iterative regularization via dual diagonal descent[END_REF]; [START_REF] Calatroni | Accelerated iterative regularization via dual diagonal descent[END_REF].

' Sparse recovery and compressed sensing In the context of sparse recovery the implicit regularization approach has been considered in [START_REF] Osher | Sparse recovery via differential inclusions[END_REF], and also in [START_REF] Vaškevičius | Implicit regularization for optimal sparse recovery[END_REF]. Matching pursuit [START_REF] Mallat | Matching pursuit with timefrequency dictionaries[END_REF] is a computational procedure which can be used to select relevant components, but it is not clear from the theoretical point of view how to early stop the iterations .

' Exact regularization Another possible approach is to use the notion of exact regularization [START_REF] Friedlander | Exact regularization of convex programs[END_REF][START_REF] Schopfer | Exact regularization of polyhedral norms[END_REF]. The latter refers to solving min w Jpwq `αQpwq, s.t. Xw " y , where Q is strongly convex and to showing that there exists a value of α such that this new problem and ?? have the same minimizer. Then, known iterative regularization results of the strongly convex case [START_REF] Matet | Don't relax: early stopping for convex regularization[END_REF]

Assume that Problem (19) admits a solution w ‹ satisfying the following qualification condition,

pDθ ‹ P R n q ´XJ θ ‹ P BJpw ‹ q . (QC)
Reasoning as in the proof of the Fenchel-Rockafellar duality theorem (Peypouquet, 2015, Thm. 3.51) it follows that strong duality holds, and θ ‹ is a solution of Problem (21). Primal-dual solutions are thus characterized by the first order conditions, ´XJ θ ‹ P BJpw ‹ q and Xw ‹ " y .

We stress the fact that we assume the existence of a solution w ‹ satisfying (QC), but the primal solution is not necessarily unique. On the other hand, from strong duality we get also that, for every primal solution, there exists a dual one such that (QC) (and so Equations ( 9) and ( 22)) is verified.

In this general case, the Chambolle-Pock's algorithm (with interpolation parameter equal to 1) is given by y k`1 " prox τ g ‹ py k `τ Kp2x k ´xk´1 qq,

x k`1 " prox σf px k ´σK J y k`1 q.

Notice that the CP algorithm, except for the interpolation, treats the primal and the dual problem in a symmetric way. In particular, we can cast the method both for Problems (23) and ( 24). In order to apply the latter to our dual problem, we set f " xy, ¨y, g " J ‹ and K " ´XJ . Then g ‹ " J, prox σf pθq " θ ´σy and we recover Equation ( 6):

w k`1 " prox τ J `wk ´τ X J p2θ k ´θk´1 q ˘, θ k`1 " θ k `σ `Xw k`1 ´yδ ˘.

The latter uses, in the update of the variable w, an interpolation of θ with the value at the previous step.

As we already remarked, we could also apply the CP algorithm directly to the primal problem, setting f " J, g " ι tyu and K " X. Then g ‹ " xy, ¨y and prox τ g ‹ pθq " θ ´τ y, leading to the following method:

θ k`1 " θ k `τ `X p2w k ´wk´1 q ´yδ ˘, w k`1 " prox σJ `wk ´σX J θ k`1 ˘.

In this case, in the update of the variable θ, we use an interpolation of w. In general, the two versions should not differ in a significant manner. Nevertheless, the error we consider affects the data y δ and so its nature is not symmetric. Then, a different choice for the interpolation can play a role. In this work, we put emphasis in Algorithm (6) because it is the one for which we have proximal errors in the non-extrapolated step.

C Proofs

C.1 Lemmas Lemma 5. Let pw ‹ , θ ‹ q be a primal-dual solution and pw 1 , θ 1 q a point in R p ˆRn such that Lpw 1 , θ ‹ q´Lpw ‹ , θ 1 q " 0 and Xw 1 " y. Then pw 1 , θ ‹ q is a primal-dual solution.

Proof.

Step 1: the duality gap is the Bregman divergence. Indeed, using ´XJ θ ‹ P BJpw ‹ q and Xw ‹ " y:

Lpw 1 , θ ‹ q ´Lpw ‹ , θ 1 q " Jpw 1 q ´Jpw ‹ q `xθ ‹ , Xw 1 ´yy ´xθ 1 , Xw ‹ ´yy " Jpw 1 q ´Jpw ‹ q `xX J θ ‹ , w 1 ´w‹ y "

D ´XJ θ ‹ J pw 1 , w ‹ q , (25)
Step 2: Zero duality gap plus feasibility implies primal optimality We show that if v P BJpw ‹ q and D v J pw 1 , w ‹ q " 0, then v P BJpw 1 q. Indeed, Jpw 1 q ´Jpw ‹ q ´xv, w 1 ´w‹ y " 0 and so, for all z P R p , Jpzq ě Jpw ‹ q `xv, z ´w‹ y " Jpw 1 q ´xv, w 1 ´w‹ y `xv, z ´w‹ y " Jpw 1 q `xv, z ´w1 y.

(26)

The statement follows by applying step 2 with v " ´XJ θ ‹ .

Next, we recall the result that allows us to control the non-vanishing error. It is a discrete version of Bihari's Lemma and a particular case of Lemma 1 in [START_REF] Schmidt | Convergence rates of inexact proximal-gradient methods for convex optimization[END_REF], where the proof can be found.

Lemma 14. Assume that pu j q is a non-negative sequence and that λ ě 0, S ě 0 with S ě u 2 0 . If u 2 t ď S `λ ř t j"1 u j , then

u t ď λt 2 `«S `ˆλt 2 ˙2ff 1 2 .
So, in particular, u t ď λt `?S.

C.2 Preliminary estimates

Lemma 15 (One step estimate). Defining θk :" 2θ k ´θk´1 , the updates of (6) for the noisy problem read as:

w k`1 " prox τ J ´wk ´τ X J θk ¯, (27)
θ k`1 " θ k `σ `Xw k`1 ´yδ ˘. (28)
Then, for any pw, θq P R p ˆRn , we have the following estimate:

V pz k`1 ´zq ´V pz k ´zq `V pz k`1 ´zk q `rLpw k`1 , θq ´Lpw, θ k`1 qs `xθ k`1 ´θ, y δ ´yy `xθ k`1 ´θ k , X pw ´wk`1 qy ď 0.

(29)

Figure 1 :

 1 Figure1: For J " ¨ 1 , convergence of primal-dual iterates w k towards w ‹ , measured in norm, feasability and Bregman divergence. The Bregman divergence quickly vanishes up to numerical errors (the iterates have the same sign as the solution); yet the iterates are still far from the solution.

Figure 2 :

 2 Figure 2: Distance between noisy Basis Pursuit iterates w δk and noiseless solution w ‹ , for various values of δ. There exists a stopping time: these distances reach a minimum before converging to their limit.

Figure 3 :

 3 Figure 3: Influence of δ on the stopping time. In agreement with theory, the empirical stopping time roughly scales as 1{δ.

Figure 4 :

 4 Figure 4: Comparison of Tikhonov regularization path (top) and optimization path of Algorithm (6) (bottom) on rcv1 with 4-fold cross-validation. Minimal value reached: 0.19 (top), 0.21 (bottom). Computation time up to optimal parameter: 50 s (top); 0.5 s (bottom).

Figure 5 :

 5 Figure 5: Distance between low rank matrix completion iterates w δk and noiseless solution w ‹ , for various values of δ. There exists a stopping time: a minimum before distance is reached before the limit.

f

  ‹ p´K J yq `g‹ pyq ( .

  can be applied.

	B Duality and Chambolle-Pock algorithm
	B.1 Duality
	The Chambolle-Pock algorithm belongs to the class of primal-dual methods, designed to jointly solve ?? (the
	primal problem), and its dual. ?? rewrites as
		arg min	Jpwq `ιtyu pXwq ,	(19)
		wPR p
	with Lagrangian
		Lpw, θq " Jpwq `xθ, Xw ´yy .	(20)
	Fenchel-Rockafellar duality (Peypouquet, 2015, Sec. 3.6.2) can be applied to compute the dual problem; observing
	that ι ‹ tyu pθq " xy, θy, this dual reads:
		!	)
	arg min θPR n	J ‹ p´X J θq `ι‹ tyu pθq	" arg min θPR n

J ‹ p´X J θq `xy, θy ( .

Proceedings of the 24 th International Conference on Artificial Intelligence and Statistics (AISTATS) 2021, San Diego, California, USA. PMLR: Volume 130. Copyright 2021 by the author(s). Contact: cecio.molinari@gmail.com

Acknowledgments

This material is based upon work supported by the Center for Brains, Minds and Machines (CBMM), funded by NSF STC award CCF-1231216, and the Italian Institute of Technology. Part of this work has been carried out at the Machine Learning Genoa (MaLGa) center, Università di Genova (IT). L. R. acknowledges the financial support of the European Research Council (grant SLING 819789), the AFOSR projects FA9550-17-1-0390 and BAA-AFRL-AFOSR-2016-0007 (European Office of Aerospace Research and Development), and the EU H2020-MSCA-RISE project NoMADS -DLV-777826. S. V. acknowledges the support of INDAM-GNAMPA, Project 2019: "Equazioni integro-differenziali: aspetti teorici e applicazioni".

Proof. Consider first Equation ( 27) and the firm non-expasiveness of the proximal-point. Then we get that, for any w P R p , 0 ě w k`1 ´w

2 ´ ´wk ´τ X J θk ¯´w 2 ` w k`1 ´´w k ´τ X J θk ¯ 2 `2τ rJpw k`1 q ´Jpwqs " w k`1 ´w 2 ´ w k ´w 2 ` w k`1 ´wk 2 `2τ rJpw k`1 q ´Jpwqs `2τ xX J θk , w k ´wy `2τ xX J θk , w k`1 ´wk y " w k`1 ´w 2 ´ w k ´w 2 ` w k`1 ´wk 2 `2τ rJpw k`1 q ´Jpwqs `2τ x θk , X pw k`1 ´wqy.

Now consider Equation ( 27) and notice that the dual update can be re-written as θ k`1 " prox σxy δ ,¨y pθ k `σXw k`1 q. Similarly as before, for any θ P R n ,

Recall that z :" pw, θq and the definition of V in Equation ( 12). Divide the first inequality by 2τ , the second one by 2σ and sum-up, to get

To conclude, compute

Proof. We start from Equation ( 29), switching the index from k to t and evaluating pw, θq at the saddle-point pw ‹ , θ ‹ q. Recall that θt :" 2θ t ´θt´1 , to get

where in the last estimate we used Cauchy-Schwartz and Young inequalities, the latter with parameter σ. Then, using the definition of ω :" 1 ´τ σ X 2 op , we have

Imposing θ ´1 " θ 0 , summing-up the latter from t " 0 to t " k ´1 and using the telescopic property, we get

where in the last inequality we used again Cauchy-Schwartz and Young inequalities with parameter σ. Reordering, we obtain the claim.

Lemma 17 (Second cumulative estimate). For ε ą 0 and η " 1`ε 1´ε ě 1, define ω :" ε ´στ X 2 op . Then we have

Proof. In a similar fashion as in the previous proof, we start again from Equation ( 29), switching the index from k to t and evaluating pw, θq at the saddle-point pw ‹ , θ ‹ q. Since θt " θ t `pθ t ´θt´1 q " θ t `σpXw t ´yq and θ t`1 ´θt " σpXw t`1 ´yδ q, we get

So,

V pz t`1 ´z‹ q ´V pz t ´z‹ q `1 2τ w t`1 ´wt 2 `σ 2 Xw t`1 ´y 2 `rLpw t`1 , θ ‹ q ´L pw ‹ , θ t`1 qs ď σxX pw t`1 ´wt q , Xw t`1 ´yy `xθ t`1 ´θ‹ , y ´yδ y `σxXw t`1 ´y, y δ ´yy ´σ 2 y δ ´y

In the last inequality we used three times Cauchy-Schwartz inequality, the bound on the error given by y δ ´y ď δ and two times Young inequality with parameters ε ą 0 and η " 1`ε 1´ε ą 0. Then, re-ordering and recalling the definitions of ω :" ε ´στ X 2 op , we obtain

Summing-up the latter from t " 0 to t " k ´1, by telescopic property, we get

By trivial manipulations, we get the claim.

C.3 Proof of Proposition 7

Proposition 7 (Stability). Under Assumption 1, let ε P p0, 1q and assume that the step-sizes are such that στ ď ε{ X 2 op . Then,

and

Proof. Inequality in Equation (30) holds true for every k ě 1. Then, recalling that Lpw, θ ‹ q ´Lpw ‹ , θq ě 0 for every pw, θq P R p ˆRn and that ω ě 0 by assumption, for every t ě 1 we have that

Apply Lemma 14 to Equation (32) with u j " θ j ´θ‹ , S " 2σV pz 0 ´zq and λ " 2σδ, to get θ t ´θ‹ ď 2σδt `a2σV pz 0 ´zq.

In particular, for 1 ď t ď k, we have θ t ´θ‹ ď 2σδk `a2σV pz 0 ´zq.

Insert the latter in Equation ( 30), to obtain

´2σδk `a2σV pz 0 ´zq "

V pz 0 ´zq `δk a 2σV pz 0 ´zq `2σδ 2 k 2 ď ´aV pz 0 ´zq `?2σδk ¯2 .

Iterative regularization with convex regularizers By Jensen's inequality, we get the claim.

For the second result, recall that, from Equation ( 33 θ t ´θ‹ `σ pη ´1q δ 2 k 2 `V pz 0 ´z‹ q ď a 2σV pz 0 ´z‹ qδk `2σδ 2 k 2 `σ pη ´1q δ 2 k 2 `V pz 0 ´z‹ q.

By Jensen's inequality, rearranging the terms, and taking η " 1`ε 1´ε , we get the claim.


Edouard Januel

Jérôme De Seze

Patrick Vermersch

Elisabeth Maillart

Bertrand Bourre

Julie Pique

Xavier Moisset

Caroline Bensa

Adil Maarouf

Jean Pelletier

Sandra Vukusic

Bertrand Audoin

Céline Louapre

Stéphane Beltran

Eric Berger

Kevin Bigaut

Nathalie Derache

Salimata Gassama

Olivier Heinzlef

Laurent Kremer

Philippe Kerschen

Aude Maurousset

Vito A G Ricigliano

Pierre-Yves Garcia

Post-vaccine COVID-19 in patients with multiple sclerosis or neuromyelitis optica

Introduction

Anti-CD20 (ocrelizumab; rituximab) are immunosuppressive therapies that profoundly deplete B lymphocytes. They are widely used to treat multiple sclerosis (MS) and Neuromyelitis Optica Spectrum Disorders (NMOSD) and have proven to be highly effective in reducing relapse rate and disability accumulation [START_REF] Hauser | Ocrelizumab versus Interferon Beta-1a in Relapsing Multiple Sclerosis[END_REF], with an increased risk of severe infections over time. During the recent Sars-CoV-2 pandemic, ocrelizumab and rituximab have been associated with a higher risk of severe COVID-19 (Odds Ratio 2.05[95% Confidence Interval 1.39;3.02]) [START_REF] Sormani | DMTs and Covid-19 severity in MS: a pooled analysis from Italy and France[END_REF]. Moreover, ocrelizumab exposure is also known to attenuate humoral responses to pneumococcal, influenzae, and tetanus toxoid vaccines [START_REF] Bar-Or | Effect of ocrelizumab on vaccine responses in patients with multiple sclerosis: The VELOCE study[END_REF], raising suspicion that B cell depletion may also impair vaccine efficacy. In a recent Israelian cohort study, only 22.7% (10 of 44) and 3.8% (1 of 26) of MS patients treated by anti-CD20 and Fingolimod respectively developed anti-spike protein-based positive serology following 2 doses of BNT162b2 vaccine [START_REF] Achiron | Humoral immune response to COVID-19 mRNA vaccine in patients with multiple sclerosis treated with highefficacy disease-modifying therapies[END_REF].

The first COVID-19 vaccine (BNT162b2) was approved in France and Europe on December 21, 2020. On As anti-CD20 therapies exposure is associated with COVID-19 severity and with lower anti-spike humoral immunity after vaccination, we aimed to assess if COVID-19 had occurred among previously vaccinated patients, and determine the characteristics of these patients.

Case series:

We used the COVISEP registry, the French database of COVID-19 in patients with MS or NMOSD, comprising 1650 MS and 67 NMOSD COVID-19 infected patients, on July 30, 2021. We extracted all cases of COVID-19 occurring since December 21, 2020, corresponding to 460 cases.

Sixteen MS patients were infected with Covid-19 more than 7 days after receiving their first dose of COVID-19 vaccine, and before receiving the second dose. Of them, 6 patients were treated with anti-CD20, 3 with natalizumab, 1 with fingolimod, teriflunomide, glatiramer, cladribine and dimethylfumarate and 2 had no treatment.

Eighteen patients were diagnosed with COVID-19 while they had received two doses of BNT162b2 vaccine, 17 had MS and 1 had aquaporine-4 positive NMOSD. The characteristics of these patients are summarized in Table 1. Median age (P25;P75) of post-vaccine COVID-19 patients was 36.8 (31.0;44.4) years, 14 (77.8%) were women, and only one had a comorbidity (Patient 13 had diabetes). All patients were diagnosed due to typical COVID-19 symptoms occurrence. COVID-19 first symptoms occurred at a median (P25;P75) of 21 (9;60) days after the second dose. Thirteen of the 18 patients (72.2%) were treated with anti-CD20 therapies, with a median (P25;P75) of 5 ( 4 In the COVISEP registry, 18 patients have been infected with COVID-19 after 2 doses of BNT162b2 vaccine, 13 of them were treated with anti-CD20 therapies, and four with fingolimod, indicating that COVID-19 occurrence after a double vaccine dose was not uncommon among these patients. By comparison, in the BNT162b2 randomized trial, including exclusively non-immunocompromised patients, with a median follow up of two months, 9 patients on 18556 were infected after two doses. In contrast, only one case of post-vaccine COVID-19 occurred among patients treated with other DMTs and none among untreated patients. In France in June 2021, among the 59301 MS treated patients, 6185 MS patients were with Ocrelizumab (10.4% of MS treated patient) [START_REF] Gers | Groupement pour l'Élaboration et la Réalisation de Statistiques[END_REF]. The number of patients treated with Rituximab in France was not available, but these patients represented one third of MS and NMOSD anti-CD20 treated patients at the Pitié-Salpêtrière Hospital. Thus, we can assume that approximately 9000 MS or NMOSD patients were treated with anti-CD20 in France. Considering fingolimod, 10352 patients were treated with this drug in June 2021 (17.5% of MS treated patients) [START_REF] Gers | Groupement pour l'Élaboration et la Réalisation de Statistiques[END_REF].

Sixty percent of French MS patients had received a double dose of COVID-19 vaccine on July 18, 2021( 6), but the exact proportion of vaccinated patients according to their treatment is not available. Furthermore, the variation of COVID-19 incidence rate and the appearance of COVID-19 new variants since the end of BNT162b2 randomized trial did not allow us to calculate and compare the expected number of MS and NMOSD post-vaccine COVID-19 cases per disease modifying therapies. However, in the COVISEP registry, the strikingly high proportion of anti-CD20 treated patients among post-vaccine cases (74.0%) while they were 20.0% among non-vaccinated COVID-19 is a strong argument to suspect that anti-CD20 exposure may be associated with a lower COVID-19 vaccination efficacy.

Due to the concerns regarding low immunogenicity of COVID-19 vaccines among anti-CD20 treated patients, French health authorities recommended on April 6, 2021 to perform a third dose four weeks after the second dose in anti-CD20 treated population. Our results reinforce this recommendation, even if we cannot be sure that this third dose would guarantee the same protection as vaccinated nonimmunocompromised patients. In our series, one 44-year-old anti-CD20 treated patient was infected with COVID-19 despite triple dose vaccination and another 43-year-old patient treated with ocrelizumab developed COVID-19 while she had been already infected more than one year before, and in the interim had received two doses of BNT162b2 vaccine. In contrast, among non-immunocompromised patients with previous evidence of COVID-19 infection, a single dose of mRNA vaccine is usually sufficient to induce a strong immunological response [START_REF] Krammer | Antibody Responses in Seropositive Persons after a Single Dose of SARS-CoV-2 mRNA Vaccine[END_REF].

Among non-anti-CD20 treated patients, four COVID-19 post-vaccine patients were treated with fingolimod (23.5% of MS post-vaccine COVID-19, compared with 11.5% among non-vaccinated COVID-19). Even if this ratio difference is not as high as for anti-CD20, fingolimod represent a substantial part of post-vaccine COVID-19, while only one patient on interferon, and none on other DMT or untreated patient were recorded in the COVISEP registry.

To our knowledge, another case of COVID-19 after a double dose vaccination was previously described in a 52-year-old patient who had been vaccinated two weeks after his last ocrelizumab infusion [START_REF] Chilimuri | COVID-19 Vaccine Failure in a Patient with Multiple Sclerosis on Ocrelizumab[END_REF]. In our series, 3 patients had been vaccinated less than 3 months after the last anti-CD20 infusion. This short delay between infusion and vaccination may enhance the risk of vaccine inefficacy. Three patients without prior COVID-19 infection developed symptoms within 7 days after the second vaccine dose, thus before obtaining vaccine full efficacy based on the BNT162b2 randomized trial [START_REF] Polack | Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine[END_REF]. However, in the BNT162b2 randomized trial, efficacy within 7 days after second dose was estimated to be 91% [START_REF] Polack | Safety and Efficacy of the BNT162b2 mRNA Covid-19 Vaccine[END_REF].

In this preliminary study, we observed that 72% of the patients with COVID-19 after double dose vaccine were treated with anti-CD20 therapies, while they represented only 20% of non-vaccinated COVID-19. This could prompt us to suspect that ocrelizumab and rituximab may not only impair COVID-19 vaccine humoral response, but more importantly may impair vaccination efficacy. However, we must keep in mind that due to public health decisions, MS patients treated with anti-CD20 were prioritized for COVID-19 vaccination since February 11, 2021. Thus, COVID-19 vaccine inefficacy in anti-CD20 treated patients may have been overstated in our study, as they may have been more quickly vaccinated than other DMT. However, this difference is decreasing over time, indeed, on 18 July 2021, 60% of all MS patients were vaccinated. Vaccination effectiveness assessment for other MS or NMOSD immunosuppressive treatments (particularly fingolimod) warrants further longitudinal follow-up. As a limitation, we did not have information about COVID-19 vaccine serological response for post-vaccine COVID-19 cases.

Importantly, despite occurring in immunocompromised patients, COVID-19 severity was mild in our series, except for one 67-year-old patient, and none of the patients died.

Conclusion

We describe 18 cases of COVID-19 infection occurring despite double dose BNT162b2 vaccine among patients with MS or NMOSD in the COVISEP registry, 13 of them were treated with anti-CD20 therapies. COVID-19 severity was mild. These preliminary findings stress the need for a prospective clinical and biological follow up on COVID-19 vaccine efficacy among this population exposed to immunosuppressive therapies, and further studies are needed to assess the impact of prior vaccination on COVID-19 severity among immunocompromised patients.

.

February 11 ,

 11 2021, MS patients treated with anti-CD20 therapies became a priority population for vaccination. Later on, vaccination was recommended for all MS patients older than 50 years on March 2, 2021, and for all adult MS patients on May 12, 2021, although most patients could have access to vaccination in early April 2021.

  ;6) cycles. Only one had low concentration of IgG (Patient 16: IgG 5.4 g/L). In contrast, among non-vaccinated COVID-19 patients in the COVISEP registry, 20.0% were treated with anti-CD20. Anti-CD20 exposure was associated with postvaccine COVID-19 occurrence among MS and NMOSD patients (OR 10.4; CI [3.7;29.3], p<0.001). One patient (patient 18), treated with ocrelizumab, had previously been infected with Covid-19 thirteen months before, with laboratory evidence of anti-S COVID-19 immunity. Four of the 17 MS post-vaccine patients (23.5%) were treated with fingolimod, while the were 11.5% among MS non-vaccinated COVID-19 (OR for post-vaccine COVID-19 occurrence in MS patients treated with fingolimod: 2.3; CI [0.75;7.1], p=0.146). Only one patient was treated with another disease modifying therapy (DMT) (patient 17: interferon beta), and none among non-treated patients. In three patients without prior infection, COVID-19 first symptoms began within 7 days after the second dose. Covid-19 severity was generally mild, but one patient (patient 1) required intensive care hospitalization, without the need for mechanical ventilation. Discussion:

Table 1 :

 1 Characteristics of patients with COVID-19 occurrence despite double dose vaccination Expanded Disability Status Scale. MS: Multiple sclerosis. NMO: Neuromyelitis optica. Comorbidities: cardiovascular, pulmonary, diabetes, obesity (BMI>30), smoking. *Patient N°18 had previously been infected with COVID-19 13 months before second infection, with laboratory evidence of anti-COVID-19 immunity.

	Case	Diag	Age	EDSS	Comorbi	Gender Diseas	Numbe	Immun	Delay	Vaccin	Delay 2 nd or 3 rd	COVID-19
		nosis	(year		dity		e	r	of	e		Last anti-	e type	Dose vaccine -	severity
			s)				modify	anti-		globuli	CD20	and	COVID-19 first
							ing	CD20		n	G	infusion	numbe	symptoms
							therap	cycles	levels		-Vaccine	r	of	(days)
							y	before	before	1 st dose	doses
								vaccina	COVID-	(weeks)
								tion		19
										infecti
										on
										(g/L)
		COVID-19 symptoms > 7 days after BNT162b2 2 nd dose
	1	MS	67	5.5	None	F	Rituxi	6		9.1		17	BNT16	14	Hospitalized
							mab						2b2 2	and receiving
													doses	noninvasive
														ventilation or
														high-flow
														oxygen
	2	MS	24	1	None	F	Ocreliz	1		8.9		30	BNT16	13	Not
							umab						2b2 2	hospitalized,
													doses	limitation	on
														activities
	3	MS	35	0	None	F	Rituxi	2		10.2		98	BNT16	20	Not
							mab						2b2 2	hospitalized,
													doses	limitation	on
														activities
	4	MS	31	3	None	M	Rituxi	8		12.1		3	BNT16	55	Not
							mab						2b2 2	hospitalized,
													doses	limitation	on
														activities
	5	NMO	56	2.5	None	F	Rituxi	5		8.26		13	BNT16	9	Not
		AQP					mab						2b2 2	hospitalized,
		4											doses	limitation	on
														activities
	6	MS	38	0	None	F	Rituxi	5		8.79		19	BNT16	12	Not
							mab						2b2 2	hospitalized, no
													doses	limitations on
														activities
	7	MS	47	6	None	F	Ocreliz	4		8.05		18	BNT16	57	Not
							umab						2b2 2	hospitalized,
													doses	limitation	on
														activities
	8	MS	30	0	None	M	Ocreliz	4		10.96		19	BNT16	86	Not
							umab						2b2 2	hospitalized,
													doses	limitation	on
														activities
	9	MS	48	2.5	None	F	Fingoli	-		-		-	BNT16	23	Hospitalized
							mod						2b2 2	and did not
													doses	required
														supplemental
														oxygen
	10	MS	24	0	None	F	Fingoli	-		-		-	BNT16	45	Not
							mod						2b2 2	hospitalized,
													doses	limitation	on
														activities
	11	MS	32	0	None	F	Ocreliz	6		NA		13	BNT16	61	Not
							umab						2b2 2	hospitalized,
													doses	limitation	on
														activities
	12	MS	28	0	None	F	Fingoli	-		-		-	BNT16	89	Not
							mod						2b2 2	hospitalized,
													doses	limitation	on
														activities
	13	MS	44	6	Diabetes	M	Rituxi	4		11.9		13	BNT16	29 (after 3 rd	Not
							mab						2b2 3	dose)	hospitalized,
													doses	limitation	on
														activities