Jérôme Bolte

Edouard Pauwels

Samuel Vaiter

Automatic differentiation of nonsmooth iterative algorithms

Differentiation along algorithms, i.e., piggyback propagation of derivatives, is now routinely used to differentiate iterative solvers in differentiable programming. Asymptotics is well understood for many smooth problems but the nondifferentiable case is hardly considered. Is there a limiting object for nonsmooth piggyback automatic differentiation (AD)? Does it have any variational meaning and can it be used effectively in machine learning? Is there a connection with classical derivative? All these questions are addressed under appropriate nonexpansivity conditions in the framework of conservative derivatives which has proved useful in understanding nonsmooth AD. For nonsmooth piggyback iterations, we characterize the attractor set of nonsmooth piggyback iterations as a set-valued fixed point which remains in the conservative framework. This has various consequences and in particular almost everywhere convergence of classical derivatives. Our results are illustrated on parametric convex optimization problems with forward-backward, Douglas-Rachford and Alternating Direction of Multiplier algorithms as well as the Heavy-Ball method.

Dedicated to the memory of Andreas Griewank -a pioneer in automatic differentiation and optimization -who passed away on September 2021.

1 Introduction

x k (θ) J x k (θ)
x(θ) J pb x (θ) nonsmooth autodiff k → +∞ derivative? limit?

Figure 1: We study existence and meaning of J pb x as a derivative of x, compatible with automatic differentiation of the iterates (x k (θ)) k∈N .

Differentiable programming. We consider a Lipschitz function F : R p × R m → R p , representing an iterative algorithm, parameterized by θ ∈ R m , with Lipschitz initialization x 0 : θ → x 0 (θ) and

x k+1 (θ) = F (x k (θ), θ) = F θ (x k (θ)), (1)
where F θ := F (•, θ), under the assumption that x k (θ) converges to the unique fixed point of F θ : x(θ) = fix(F θ ). Such recursion represent for instance algorithms to solve an optimization problem min x h(x) (e.g. empirical risk minimization), such as gradient descent: F (x, θ) = x -θ∇h(x). But (1) could also be a fixed-point equation such as a deep equilibrium network [START_REF] Bai | Deep equilibrium models[END_REF].

In the last years, a paradigm shift occurred: such algorithms are now implemented in algorithmic differentiation (AD)-friendly frameworks such as Tensorflow [START_REF] Abadi | TensorFlow: Large-scale machine learning on heterogeneous systems[END_REF], PyTorch [START_REF] Paszke | Pytorch: An imperative style, high-performance deep learning library[END_REF] or JAX [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF] to name a few. Assuming that F is differentiable, it is possible to compute iteratively the derivatives of x k+1 with respect to θ using the differential calculus chain rule resulting in so called "piggyback" recursion:

∂ ∂θ x k+1 (θ) = ∂ 1 F (x k (θ), θ) • ∂ ∂θ x k (θ) + ∂ 2 F (x k (θ), θ), (2)
where ∂ ∂θ x k is the Jacobian of x k with respect to θ. In practice, automatic differentiation frameworks do not compute the full Jacobian, but compute either vector-Jacobian products in reverse-mode (or backpropagation) [START_REF] Rumelhart | Learning representations by back-propagating errors[END_REF] or Jacobian-vector products in forward mode [START_REF] Wengert | A simple automatic derivative evaluation program[END_REF]. We rather consider the full Jacobian, and therefore, our findings apply to both modes. We focus on two issues arising with nonsmooth recursions, illustrated in Figure 1. (i) what can be said about the chain rule [START_REF] Ablin | Super-efficiency of automatic differentiation for functions defined as a minimum[END_REF] and its asymptotics when the function F is not smooth (for example a projected gradient step)? (ii) how to interpret its asymptotics as a notion of derivative for x, the fixed point of F θ ? We propose a joint answer to both questions, providing a solid theoretical ground to the idea of algorithmic differentiation of numerical solvers involving nonsmooth components in a differentiable programming context. Related works. Algorithmic use of the chain rule [START_REF] Ablin | Super-efficiency of automatic differentiation for functions defined as a minimum[END_REF] to differentiate programs takes its root in [START_REF] Wengert | A simple automatic derivative evaluation program[END_REF], where forward differentiation was first proposed, and later in reverse mode [START_REF] Linnainmaa | The representation of the cumulative rounding error of an algorithm as a Taylor expansion of the local rounding errors[END_REF]. Along with the practical development of automatic differentiation, the question on how to prove the convergence of the iterative sequence (2) was investigated, notably in the optimization community as reviewed in [START_REF] Griewank | Piggyback differentiation and optimization[END_REF]. This is an important paper containing several ideas in differentiable programming rediscovered/reused later: implicit differentiation [START_REF] Pedregosa | Hyperparameter optimization with approximate gradient[END_REF][START_REF] Rajeswaran | Meta-learning with implicit gradients[END_REF] and its stability [START_REF] Blondel | Efficient and modular implicit differentiation[END_REF], adjoint fixed point iteration [START_REF] Bai | Deep equilibrium models[END_REF] that is a key aspect of the deep equilibrium network and linear convergence of (2) as discussed below. Notably, the linear convergence of the Jacobians was investigated in [START_REF] Charles | Automatic differentiation and iterative processes[END_REF][START_REF] Griewank | Derivative convergence for iterative equation solvers[END_REF] for the forward mode and in [START_REF] Christianson | Reverse accumulation and attractive fixed points[END_REF]Theorem 2.3] for the reverse mode. This was more recently investigated -for C 2 functions -in the imaging community for primal-dual algorithms [START_REF] Chambolle | Learning consistent discretizations of the total variation[END_REF][START_REF] Bogensperger | Convergence of a Piggybackstyle method for the differentiation of solutions of standard saddle-point problems[END_REF] and in the machine learning community for gradient descent [START_REF] Mehmood | Automatic differentiation of some first-order methods in parametric optimization[END_REF][START_REF] Lorraine | Optimizing millions of hyperparameters by implicit differentiation[END_REF] and the Heavy-ball [START_REF] Mehmood | Automatic differentiation of some first-order methods in parametric optimization[END_REF] method. Note that in the specific context where F solves a min-min problem, the authors in [START_REF] Ablin | Super-efficiency of automatic differentiation for functions defined as a minimum[END_REF] proved the linear convergence of the Jacobians. The use of automatic differentiaton for nonsmooth functions was justified by the development of the notion of conservative Jacobians [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF][START_REF] Bolte | A mathematical model for automatic differentiation in machine learning[END_REF] with a nonsmooth version of the chain rule for compositional models. The correctness of automatic differentiation was also investigated in [START_REF] Lee | On correctness of automatic differentiation for non-differentiable functions[END_REF] for a large class of functions that are piecewise analytic, and also in [START_REF] Sham | Provably correct automatic sub-differentiation for qualified programs[END_REF] where a qualification condition is used to compute a Clarke Jacobian. Along with automatic differentiation, a natural way to differentiate a model such as [START_REF] Abadi | TensorFlow: Large-scale machine learning on heterogeneous systems[END_REF] is by implicit differentiation, recently applied in several works [START_REF] Bai | Deep equilibrium models[END_REF][START_REF] Agrawal | Differentiable convex optimization layers[END_REF][START_REF] Ghaoui | Implicit deep learning[END_REF]. To study these models with nonsmooth functions, an implicit function theorem [START_REF] Bolte | Nonsmooth Implicit Differentiation for Machine Learning and Optimization[END_REF] was proved for path-differentiable functions.

Contributions: Under suitable nonexpansivity assumptions, our contributions are as follows.

• We address both questions illustrated in Figure 1 for nonsmooth recursions: set-valued nonsmooth extensions of the piggyback recursion (2) have a well defined limit, described as the fixed point of subset map (Theorem 1), it is conservative for the fixed point map x. This is a nonsmooth "infinite" chain rule for AD (Theorem 2).

• For almost all θ, despite nonsmoothness, recursion (2) is well defined using the classical Jacobian and converges to the classical Jacobian of the fixed point x (Corollary 2). This has implications for both forward and reverse modes of AD.

• For a large class of functions (Lipschitz-gradient selection), it is possible to give a quantitative rate estimate (Corollary 3), namely to prove linear convergence of the derivatives.

• We show that these results can be applied to proximal splitting algorithms in nonsmooth convex optimization. We include forward-backward (Proposition 2), as well Douglas-Rachford (Proposition 3) and ADMM, a numerical illustration of the convergence of derivatives is given in Figure 2.

• We also illustrate that, contrarily to the smooth case, nonsmooth piggy back derivatives of momentum methods such as the Heavy-ball algorithm, may diverge even if the iterates converge linearly (Proposition 4).

Notations. A function

f : R p → R m is locally Lipschtiz if, for each x ∈ R n ,
there exists a neighborhood of x on which f is Lipschitz. Denoting by R ⊆ R p , the full measure set where f is differentiable, the Clarke Jacobian [START_REF] Frank | Optimization and nonsmooth analysis[END_REF] at x ∈ R p is defined as

Jac c f (x) = conv M ∈ R m×p , ∃(y k ) k≥0 s.t. lim k→∞ y k = x, y k ∈ R, lim k→∞ ∂f ∂y (y k ) = M .
(3) The Clarke subdifferential, ∂ c f is defined similarly. Given two matrices A, B with compatible dimension, [A, B] is their concatenation. For a set X , convX is its convex hull (the smallest convex set containing X ). The symbol B denotes a unit ball, the corresponding norm should be clear from the context.

Nonsmooth piggyback differentiation

We first show how the use of the notion of conservative Jacobians allow us to justify rigorously the existence of a nonsmooth equivalent of piggyback iterations in (2) that is compatible with AD.

Conservative Jacobians. Conservative Jacobians were introduced in [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF] as a generalization of derivatives to study automatic differentiation of nonsmooth functions. Given a locally Lipschitz continuous function f : R p → R m , we say that the set-valued J : R p ⇒ R p×m is a conservative Jacobian for the path differentiable f if J has a closed graph, is locally bounded and nowhere empty with

d dt f (γ(t)) = J(γ(t)) γ(t) a.e. ( 4
)
for any γ : [0, 1] → R p absolutely continuous with respect to the Lebesgue measure. Conservative gradients are defined similarly. We refer to [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF] for extensive examples and properties of this class of function. Let us mention that the classes convex functions, definable functions, or semialgebraic functions are contained in the set of path differentiable functions. Given D f : R p ⇒ R p , a conservative gradient for f : R p → R, we have:

• (Clarke subgradient), for all x ∈ R p , ∂ c f (x) ⊂ conv(D f (x)). • (Gradient almost everywhere) D f (x) = {∇f (x)} for almost all x ∈ R p .
• (Calculus) differential calculus rules preserve conservativity, e.g. sum and compositions of conservative Jacobians are conservative Jacobians. An important point is that D f can be used as a first order optimization oracle for methods of gradient type, while preserving usual convergence guaranties [START_REF] Bolte | A mathematical model for automatic differentiation in machine learning[END_REF].

Piggyback differentiation of recursive algorithms. The following is standing throughout the text Assumption 1 (The conservative Jacobian of the iteration mapping is a contraction) F is locally Lipschitz, path differentiable, jointly in (x, θ), and J F is a conservative Jacobian for F . There exists 0 ≤ ρ < 1, such that for any (x, θ) ∈ R p × R m and any pair [A, B] ∈ J F (x, θ), with A ∈ R p×p and B ∈ R p×m , the operator norm of A is at most ρ.

Under Assumption 1, F θ is a strict contraction so that (x k (θ)) k∈N converges linearly to x(θ) = fix(F θ ) the unique fixed point of the iteration mapping F θ . More precisely, for all k ∈ N, we have

x k -x(θ) ≤ ρ k x 0 -F θ (x 0 ) 1 -ρ .
Furthermore, for every k ∈ N, let us define the following set-valued piggyback recursion:

J x k+1 (θ) = {AJ + B, [A, B] ∈ J F (x k (θ), θ), J ∈ J x k (θ)} . (PB)
We will show in Section 3 that (PB) plays the same role as (2) in the nonsmooth setting. Note that one can recursively evaluates a sequence J k ∈ J x k , k ∈ N as follows

J k+1 = A k J k + B k where [A k , B k ] ∈ J F (x k (θ), θ), (5)
which corresponds to the operations actually implemented in nonsmooth AD frameworks.

Remark 1 (Local contractions) Assumption 1 may be relaxed locally as follows: for all θ, the fixed point set fix(F θ ) of the iteration mapping F θ is a singleton xθ such that x k (θ) → x(θ) as k → ∞, and the operator norm condition on J F in Assumption 1 holds at the point (x(θ), θ). By graph closedness of J F , in a neighborhood of (x(θ), θ), F θ is a strict contraction and the operator norm condition on J F holds, possibly with a larger contraction factor ρ. After finitely many steps, the iterates (x k ) k∈N remain in this neighborhood and all our convergence results hold, due to their asymptotic nature.

3 Asymptotics of nonsmooth piggyback differentiation Affine iterations by packets of matrices. Let J ⊂ R p×(p+m) be a compact subset of matrices such that any matrix of the form [A, B] ∈ J with A ∈ R p×p is such that A has operator norm at most ρ < 1. We let J act naturally on the matrices of size p × m as follows J : R p×m ⇒ R p×m the function from R p×m to subsets of R p×m which is defined for each X ∈ R p×n as follows: J (X) = {AX + B, [A, B] ∈ J }. This defines a set-valued map through, for any X ⊂ R p×m ,

J (X ) = {AX + B, [A, B] ∈ J, X ∈ X }. (6)
On the model of recursions of the form (PB), we consider sequences (X k ) k∈N of subsets of R p×m satisfying the recursion X k+1 = J (X k ).

We have the following instance of the Banach-Picard theorem (proved in Appendix A).

Theorem 1 (Set-valued affine contractions) Let J ⊂ R p×(p+m) be a compact subset of matrices as above with ρ < 1. Then there is a unique nonempty compact set fix(J ) ⊂ R p×m satisfying fix(J ) = J (fix(J )), where the action of J is given in [START_REF] Heinz H Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF].

Let (X k ) k∈N be a sequence of compact subsets of R p×m , such that X 0 = ∅, and satisfying the recursion [START_REF] Bertrand | Implicit differentiation of lasso-type models for hyperparameter optimization[END_REF]. We have for all k

∈ N dist(X k , fix(J )) ≤ ρ k dist(X 0 , J (X 0 )) 1 -ρ ,
where dist is the Hausdorff distance related to the Euclidean norm on p × m matrices.

An infinite chain rule and its consequences

Define the following set-valued map based on the fix operator from Theorem 1,

J pb x : θ ⇒ fix [J F (x(θ), θ)]
. where x(θ) is the unique fixed point of the algorithmic recursion. Note that since x(θ) = fix(F θ ), we also have equivalently that J pb

x is the fixed-point of the Jacobian of the fixed-point of F θ :

J pb x : θ ⇒ fix [J F (fix(F θ ), θ)] .
We have the following (proved in Appendix B).

Theorem 2 (A conservative mapping for the fixed point map) Under Assumption 1, J pb x is well-defined, and is a conservative Jacobian for the fixed point map x.

Combining with Theorem 1 ensures the convergence of the set-valued piggyback iterations (PB).

Corollary 1 (Convergence of the piggyback derivatives) Under Assumption 1, for all θ, the recursion (PB) satisfies

lim k→∞ gap(J x k (θ), J pb x (θ)) = 0. (8)
Unrolling the expression of J x k , we can rewrite (8) as a set-valued product such that

lim K→+∞ gap K k=0 J F (x k (θ), θ), J pb x (θ) = 0.
In plain words, this a limit-derivative exchange result: Asymptotically, the gap between the automatic differentiation of x k and the derivative of the limit is zero. This implies in particular that the recursion (5) produces bounded sequences and all its accumulation points are in J pb x . Using the fact that conservative Jacobians equal classical Jacobians almost everywhere [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF], this implies convergence of derivatives in a classical sense.

Corollary 2 (Convergence of the classical piggyback derivatives) Under Assumption 1, for almost all θ, the classical Jacobian ∂ ∂θ x k (θ), is well defined for all k and converges towards the classical Jacobian of x:

lim k→∞ ∂ ∂θ x k (θ) = ∂ ∂θ x(θ).
Remark 2 (Connection to implicit differentiation) The authors in [START_REF] Bolte | Nonsmooth Implicit Differentiation for Machine Learning and Optimization[END_REF] proved a qualification-free version of the implicit function theorem. Assuming that for every [A, B] ∈ J(x(θ), θ), the matrix I -A is invertible, we have that

J imp x : θ ⇒ (I -A) -1 B, [A, B] ∈ J F (x(θ), θ) (9)
is a conservative Jacobian for x. Under Assumption 1, one has J imp x (θ) ⊂ J pb x (θ) for all θ. Unfortunately, as soon as F is not differentiable, the inclusion may be strict, see details in Appendix C.

Consequence for algorithmic differentiation

Given k ∈ N, θ ∈ R m , wk ∈ R p , the following algorithms allow us to compute ẋk = J k θ using the forward mode of automatic differentation (Jacobian Vector Products, JVP), and θT k = wT k J k using the backward mode of automatic differentiation (Vector Jacobian Products, VJP).

Algorithm 1: Algorithmic differentiation of recursion (1), forward and reverse modes

Input:

k ∈ N, θ ∈ R m , θ ∈ R m , wk ∈ R p , initialization function x 0 (θ), recursion function F (x, θ), conservative Jacobians J F (x, θ) and J x0 (θ). Initialize: x 0 = x 0 (θ) ∈ R p .
Forward mode (JVP):

ẋ0 = J θ, J ∈ J x0 (θ). for i = 1, . . . , k do x i = F (x i-1 , θ) ẋi = A i-1 ẋi-1 + B i-1 θ [A i-1 , B i-1 ] ∈ J F (x i-1 , θ) Return: ẋk Reverse mode (VJP): θk = 0. for i = 1, . . . , k do x i = F (x i-1 , θ) for i = k, . . . , 1 do θk = θk + B T i-1 wi wi-1 = A T i-1 wi [A i-1 , B i-1 ] ∈ J F (x i-1 , θ) θk = θk + J T w0 , J ∈ J x0 (θ)

Return: θk

The following result is a consequence of Corollary 2 combined with algorithmic differentiation arguments, its proof is given in Appendix C.

Proposition 1 (Convergence of VJP and JVP

) Let k ∈ N, θ ∈ R m , wk ∈ R p , x k ∈ R p , ẋk ∈ R p , θT k ∈ R m be as in Agorithm 1 under Assumption 1. Then for almost all θ ∈ R m , ẋk → ∂ x ∂θ θ. Assume furthermore that, as k → ∞, wk → w (for example, wk = ∇ (x k ) for a C 1 loss ), then for almost all θ ∈ R m , θT k → wT ∂ x ∂θ . Remark 3
In addition to Proposition 1, in both cases, for all θ, all accumulation points of both ẋk and θT k are elements of J pb x θ and wT J pb x respectively.

Linear convergence rate for semialgebraic piecewise smooth selection function

Semialgebraic functions are ubiquitous in machine learning (piecewise polynomials, 1 , 2 norms, determinant matrix rank . . . ). We refer the reader to [START_REF] Bolte | A mathematical model for automatic differentiation in machine learning[END_REF] for a thorough discussion of their extensions, and use in machine learning. For more technical details, see [START_REF] Coste | Istituti editoriali e poligrafici internazionali[END_REF][START_REF] Coste | An introduction to semialgebraic geometry[END_REF] for introductory material on semialgebraic and o-minimal geometry.

Lipschitz gradient selection functions. Let F : R p → R q be semialgebraic and continuous. We say that F has a Lipschitz gradient selection (s, F 1 , . . . , F m ) if s : R p → (1, . . . , m) is semialgebraic and there exists L ≥ 0 such that for i = 1 . . . , m, F i : R p → R p is semialgebraic with L-Lipschitz Jacobian, and for all x ∈ R p , F (x) = F s(x) (x).

For any x ∈ R p , set I(x) = {i ∈ {1, . . . , m} , F (x) = F i (x)}. The set-valued map J s F : R p ⇒ R p×q given by

J s F : x ⇒ conv ∂F i ∂x (x), i ∈ I(x) ,
is a conservative Jacobian for F as shown in [START_REF] Bolte | A mathematical model for automatic differentiation in machine learning[END_REF]. Here ∂Fi ∂x denotes the classical Jacobian of F i . Let us stress that such a structure is ubiquitous in applications [START_REF] Bolte | A mathematical model for automatic differentiation in machine learning[END_REF][START_REF] Lee | On correctness of automatic differentiation for non-differentiable functions[END_REF].

Rate of convergence. We may now strengthen Corollary 1 by proving the linear convergence of piggyback derivatives towards the fixed point. The following is a consequence of the fact that the proposed selection conservative Jacobians of Lipschitz gradient selection functions are Lipschitz-like (Lemma 3 in Appendix D.1). Note that semialgebraicity is only used as a sufficient condition to ensure conservativity of the selection Jacobian together with this Lipschitz like property. It could be relaxed if it can be guaranteed by other means, in particular one could consider the broader class of definable functions in order to handle log-likelihood data fitting terms.

Corollary 3 (Linear convergence of piggyback derivatives) In addition to Assumption 1, assume that F has a Lipschitz gradient selection structure as above. Then, for any θ and > 0, there exists C > 0 such that the recursion (PB) with

J F = J s F satisfies gap(J x k (θ), J pb x (θ)) ≤ C( √ ρ + ) k , ∀k ∈ N.
Moreover, classical Jacobians in Corollary 2 converge at a linear rate for almost all θ.

Application to proximal splitting methods in convex optimization

Consider the composite parametric convex optimization problem, where θ ∈ R m represents parameters and x ∈ R p is the decision variable

x(θ) = arg min x f (x, θ) + g(x, θ).
The purpose of this section is to construct examples of function F used in recursion (1) based on known algorithms. The following assumption will be standing throughout the section.

Assumption 2 f is semialgebraic, convex, its gradient with respecto to x for fixed θ, ∇ x f , is locally Lipschitz jointly in (x, θ) and L-Lipschitz in x for fixed θ. Semialgebraicity implies that ∇ x f is path-differentiable jointly in (x, θ), we denote by J 2 f its Clarke Jacobian. g is semialgebraic, convex in x for fixed θ, and lower semicontinuous. For all α > 0, we assume that G α (x, θ) → prox αg(•,θ) (x) is locally Lipschitz jointly in (x, θ). semialgebraicity implies that it is also path differentiable jointly in (x, θ), we denote by J Gα its Clarke Jacobian.

This assumption covers a very large diversity of problems in convex optimization as most gradient and prox operations used in practice are semialgebraic. Under Assumption 2, we will provide sufficient conditions on f and g for Assumption 1 to hold for different algorithmic recursions. These are therefore sufficient for the validity of the convergence results in Corollary 1 and Corollary 2, Proposition 1, as well Corollary 3 in the piecewise selection case. The proofs are postponed to Appendix E.

Splitting algorithms

In this section we provide sufficient condition for Assumption 1 to hold. The underlying conservative Jacobian is obtained by combining Clarke Jacobians of elementary algorithmic operations (gradient, proximal operator in Assumption 2), using the compositional rules of differential calculus [START_REF] Bolte | A mathematical model for automatic differentiation in machine learning[END_REF] and implicit differentiation [START_REF] Bolte | Nonsmooth Implicit Differentiation for Machine Learning and Optimization[END_REF]. Using [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF], such Jacobians are conservative by semialgebraicity and their combination provide conservative Jacobians for the corresponding algorithmic recursion F . These objects are explicitly constructed in Appendix E.

Forward-backward algorithm. The forward-backward iterations are given for α > 0 by

x k+1 = prox αg(•,θ) (x k -α∇ x f (x k , θ)) . ( 10
)
Proposition 2 Under Assumption 2 with 0 < α < 2 L , denote by F α : R p×m → R p the forward-backward recursion in [START_REF] Bolte | Nonsmooth Implicit Differentiation for Machine Learning and Optimization[END_REF]. For µ > 0, if either f or g is µ-strongly convex in x for all θ, then F α is a strict contraction and Assumption 1 holds. Douglas-Rachford. Given α > 0, the algorithm goes as follows

y k+1 = 1 2 (I + R αf (•,θ) R αg(•,θ) )y k , (11)
where R αf (•,θ) = 2prox αf (•,θ) -I is the reflected proximal operator, which is 1-Lipschitz (and similarly for g). Following [START_REF] Heinz H Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]Theorem 26.11], if the problem has a minimizer, then (y k ) k∈N converges to a fixed point of ( 11), ȳ such that x = prox αg (ȳ) is a solution to the optimization problem. Following [25, Theorem 1], if f is strongly convex, then R αf (•,θ) is ρ-Lipschitz for some ρ < 1 and our differentiation result applies to Douglas-Rachford splitting in this setting.

Proposition 3 Under Assumption 2 with α > 0, denote by F α : R p×m → R p the Douglas-Rachford recursion in [START_REF] Bolte | A mathematical model for automatic differentiation in machine learning[END_REF]. If f is µ-strongly convex in x for all θ, then F α is a strict contraction and Assumption 1 holds.

Alternating Direction Method of Multipliers algorithms. Consider the separable convex problem

min u,v φ θ (u) + ψ θ (v) subject to A θ u + B θ v = c θ . (12)
The alternating direction method of multipliers (ADMM) algorithm combines two partial minimization of an augmented Lagrangian, and a dual update:

u k+1 = arg min u φ θ (u) + x A θ u + α 2 A θ u + B θ v k -c θ 2 2 v k+1 = arg min v ψ θ (v) + x B θ v + α 2 A θ u k+1 + B θ v k -c θ 2 2 x k+1 = x k + α(A θ u k+1 + B θ v k+1 -c θ ). (13)
As observed in [START_REF] Gabay | Applications of the method of multipliers to variational inequalities[END_REF], the ADMM algorithm can be seen as the Douglas-Rachford splitting method applied to the Fenchel dual of problem [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF] (see Appendix E.3 for more details). More precisely, ADMM updates are equivalent to Douglas-Rachford iterations applied to the following problem min

x c θ x + φ * θ (-A θ x) f (x,θ) + ψ * θ (-B θ x) g(x,θ) . (14)
Therefore, if φ θ is strongly convex with Lipschitz gradient and A θ is injective, then ADMM converges linearly and one is able to combine derivatives of proximal operators to differentiate ADMM.

Numerical illustration.

We now detail how Figure 2 discussed in the introduction is obtained, and how it illustrates our theoretical results. We consider four scenarios (Ridge, Lasso, Sparse inverse covariance selection and trend filtering) corresponding to the four columns. For each of them, the first line shows the empirical linear rate of the iterates x k and the second line shows the empirical linear rate of the derivative ∂ ∂θ x k . All experiments are repeated 100 times and we report the median along with first and last deciles.

Forward-Backward for the Ridge. The Ridge estimator is defined for θ > 0 as x(θ) = arg min x∈R p 1 2 Axb 2 2 + θ x 2 2 Among several possibilities to solve it, one can use the Forward-Backward algorithm applied to f : (x, θ) → 1 2 Axb 2 2 and g : θ x 2 2 . Since g is strongly convex, the operator F α is strongly convex, and thus Proposition 2 may be applied.

Forward-Backward algorithm for the Lasso. Consider the Forward-Backward algorithm applied to the Lasso problem [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], with parameter θ > 0, x(θ) ∈ arg min

x∈R p 1 2 Ax - b 2 2 + θ x 1 = arg min x 1 2L Ax -b 2 2 + θ L x 1
, where L is any upper bound on the operator norm of A T A. The gradient of the quadratic part is 1 Lipschitz so we may consider the forward backward algorithm [START_REF] Bolte | Nonsmooth Implicit Differentiation for Machine Learning and Optimization[END_REF], with unit step size with f : (x, θ) → 1 2L Axb 2 2 and g : (x, θ) → θ L x 1 . A well known qualification condition involving a generalized support at optimality ensures uniqueness of the Lasso solution [START_REF] Efron | Least angle regression[END_REF][START_REF] Mairal | Complexity analysis of the lasso regularization path[END_REF]. This conditions holds for generic problem data [START_REF] Ryan | The lasso problem and uniqueness[END_REF]. Following [START_REF] Bolte | Nonsmooth Implicit Differentiation for Machine Learning and Optimization[END_REF]Proposition 5], under this qualification condition, the implicit conservative Jacobian J F is such that, at the solution x * , J F (x * ) has an operator norm of at most 1, and the matrix set I -J F only contains invertible matrices. This means that there exists ρ < 1, such that any M ∈ J F (x * ) has operator norm at most ρ. Following Remark 1, all our convergence results apply qualitatively. Note that we recover the results of [7, Proposition 2] for the Lasso. Douglas-Rachford for the Sparse Inverse Covariance Selection. The Sparse Inverse Covariance Selection [START_REF] Martin | High-dimensional graphical model selection using \ell_1-regularized logistic regression[END_REF][START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] reads x(θ) ∈ arg min x∈R n×n tr(Cx)-log det x+θ i,j |x i,j |, where C is a symmetric positive matrix and θ > 0. It is possible to apply Douglas-Rachford to f : (x, θ) → tr(Cx)log det x and g : (x, θ) → θ x 1,1 . It is known that f is locally strongly convex, indeed x →log det x is a standard self-concordant barrier in semidefinite programming [START_REF] Nesterov | Interior-point polynomial algorithms in convex programming[END_REF]. Following Remark 1, all our convergence results apply qualitatively.

ADMM for Trend Filtering. Introduced in [START_REF] Tibshirani | Adaptive piecewise polynomial estimation via trend filtering[END_REF] in statistics as a generalization of the Total Variation, the trend filtering estimator with observation θ ∈ R p reads x(θ) = arg min x∈R p 1 2 xθ 2 2 + λ D (k) x 1 , where D (k) is a forward finite-difference approximation of a differential operator of order k (here k = 2). Using ψ θ : u → λ u 1 (strongly convex),

φ θ : v → v -θ 2 2 , A θ = -I (injective), B θ = D (k)
, and c θ = 0, we can apply the ADMM to solve trend filtering.

Failure of automatic differentiation for inertial methods

In this section we consider the Heavy-Ball method for strongly convex objectives, in its global linear convergence regime. When applied to a C 2 objective, accumulation of derivatives converges to the derivative of the solution map [START_REF] Griewank | Piggyback differentiation and optimization[END_REF][START_REF] Mehmood | Automatic differentiation of some first-order methods in parametric optimization[END_REF][START_REF] Lorraine | Optimizing millions of hyperparameters by implicit differentiation[END_REF]. However, we provide a C 1,1 strongly convex parametric objective with path differentiable derivative, such that forward propagation of derivatives along the Heavy Ball algorithm contains diverging vectors for a given parameter value. In this example, one may obtain a conservative Jacobian by other means, such as implicit differentiation or algorithmic differentiation of the gradient descent algorithm, both avoiding this divergent behaviors.

Heavy-ball algorithm and global convergence

Consider a function f : R p × R m → R, and β > 0, for simplicity, when the second argument is fixed we write f θ : x → f (x, θ). Set for all x, y, θ, F (x, y, θ) = (x -∇f θ (x) + β(xy), x), consider the Heavy-Ball algorithm (x k+1 , y k+1 ) = F (x k , y k , θ) for k ∈ N.

If f θ is µ-strongly convex with L-Lipschitz gradient, then, choosing α = 1/L and β <

1 2 µ 2L + µ 2
4L 2 + 2 , the algorithm will converge globally at a linear rate to the unique solution, x(θ) [23, Theorem 4], local convergence is due to Polyak [START_REF] Polyak | Introduction to optimization[END_REF]. Furthermore, if in addition f is C 2 forward propagation of derivatives converge to the derivative of the solution [START_REF] Griewank | Piggyback differentiation and optimization[END_REF][START_REF] Griewank | Evaluating derivatives: principles and techniques of algorithmic differentiation[END_REF][START_REF] Mehmood | Automatic differentiation of some first-order methods in parametric optimization[END_REF].

A diverging Jacobian accumulation

Details and proof of the following result are given in Section F. Proposition 4 (Piggyback differentiation fails for the Heavy Ball method) Consider f : R 2 → R, such that for all θ ∈ R, f (x, θ) = x 2 /2 if x ≥ 0 and f (x, θ) = x 2 /8 if x < 0. Assume that α = 1 and β = 3/4. Then the heavy ball algorithm converges globally to 0 and ∇f is path differentiable. The Clarke Jacobian of F with respect to (x, y) at (0, 0, 0)

is J F (0, 0, 0) = conv {M 1 , M 2 }, where the product M 1 M 1 M 2 M 2 has eigenvalue -9/8.
The presence of an eigenvalue with modulus greater than 1 may produce divergence in (PB). Set

f 1 : (x, θ) → x 2 /2 if x ≥ 0 x 2 /8 if x < 0. f 2 : (x, θ) → x 2 /2 if x > 0 x 2 /8 if x ≤ 0.
Note that f 1 and f 2 are both equivalent to f as they implement the same function. With initializations x(θ) = y(θ) = θ, we run a few iterations of the Heavy Ball algorithm for θ = 0, and implement (PB) alternating between two steps on f 1 and two steps on f 2 and differentiate the resulting sequence (x k ) k∈N with respect to θ using algorithmic differentiation. The divergence phenomenon predicted by Proposition 4 is illustrated in Figure 3, while the true derivative is 0 (the sequence is constant).

Conclusion

We have developed a flexible theoretical framework to describe convergence of piggyback differentiation applied to nonsmooth recursions -providing, in particular, a rigorous meaning to the differentiation of nonsmooth solvers. The relevance of our approach is illustrated on some major composite convex optimization problems through widely used methods as forward-backward, Douglas-Rachford or ADMM algorithms. Our framework allows however to consider many other abstract algorithmic recursions and provides thus theoretical ground for more general problems such as variational inequalities or saddle point problems as in [START_REF] Chambolle | Learning consistent discretizations of the total variation[END_REF][START_REF] Bogensperger | Convergence of a Piggybackstyle method for the differentiation of solutions of standard saddle-point problems[END_REF]. As a matter for future work, we shall consider relaxing Assumption 1 to study a wider class of methods, e.g., when F is not a strict contraction. This is the appendix for "Convergence of piggyback differentiation of nonsmooth iterative solvers".

Appendices

A Properties of affine iterations on compact subsets

A.1 Banach-Picard theorem: Proof of Theorem 1

For a compact set, Z we denote by Z sup the maximal norm of elements in Z:

Z sup = sup z∈Z z .
In order to prove our fixed point result, we need first the following lemma.

Lemma 1 (Bounding Hausdorff distances) Let X , Y, Z ⊂ R p be nonempty compact sets, such that X ⊂ Y + Z and Y ⊂ X + Z then dist(X , Y) ≤ Z sup Proof : The first inclusion says that for any x ∈ X , there is y(x) ∈ Y, z(x) ∈ Z such that x = y(x) + z(x). We deduce that for any x ∈ X

min y∈Y x -y = min y∈Y y(x) -z(x) -y ≤ z(x) ≤ max z∈Z z
Therefore, max x∈X min y∈Y x-y ≤ max z∈Z z . Symmetrically, max y∈Y min x∈X x-y ≤ max z∈Z z and the result follows. We now prove Theorem 1.

Proof of Theorem 1: Recall that the action of J on matrices is defined in [START_REF] Heinz H Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF] and by A and B the projections of J on the first p and last l columns respectively, that is A = {A, ∃B, [A, B] ∈ J } and similarly for B. Note that A is a compact set and that all matrices in A have an operator norm of at most ρ. We claim that the restriction of J to compact subsets is a strict contraction in Hausdorff metric. Indeed, for any X , Y compact subsets of R p×m , we have by using Lemma 1 and noting that J preserves the inclusion,

J (X ) ⊂ J (Y + dist(X , Y)B) ⊂ J (Y) + dist(X , Y)AB ⊂ J (Y) + ρdist(X , Y)B J (Y) ⊂ J (X + dist(X , Y)B) ⊂ J (X ) + dist(X , Y)AB ⊂ J (X ) + ρdist(X , Y)B
where the last inclusion follows because AB ⊂ ρB, where B is the unit ball (for the Euclidean norm) of p × m matrices, since by assumption all square matrices in A have operator norm at most ρ. We deduce that dist(J (X ), J (Y)) ≤ ρdist(X , Y) using Lemma 1, that is the action of J on subsets of p × m matrices is ρ Lipschitz with respect to Hausdorff metric. Let us show that (X k ) k∈N remains in a bounded set, we have for all k

X k+1 sup ≤ AX k + B sup ≤ AX k sup + B sup ≤ ρ X k sup + B sup ,
which entails

X k+1 sup - B sup 1 -ρ ≤ ρ X k sup - B sup 1 -ρ .
We distinguish two cases

• if X k sup > B sup
1-ρ , then X k+1 sup gets either closer to

B sup 1-ρ or below it, in particular it decreases. • if X k sup ≤ B sup 1-ρ then X k+1 sup ≤ B sup
1-ρ and we remain below this threshold for all k.

All in all, for all k ∈ N,

X k+1 sup ≤ max X k sup , B sup 1 -ρ ≤ . . . ≤ max X 0 sup , B sup 1 -ρ , and lim sup k X k sup ≤ B sup
1-ρ . We have shown that the sequence remains in a bounded set so that the recursion actually takes place in a compact set C ⊂ R p×m which contains all the iterates in its interior, we consider the restriction of the topology to this subset. By [START_REF] Charalambos | Infinite Dimensional Analysis: A Hitchhiker's Guide[END_REF]Theorem 3.85], the closed subsets form a complete metric space. The result is an application of Banach-Picard theorem (for example [START_REF] Halsey | Real analysis[END_REF]Section 10.3]). In particular (see [START_REF] Charalambos | Infinite Dimensional Analysis: A Hitchhiker's Guide[END_REF]Theorem 3.88]), L is the unique fixed point of J and it is closed and bounded, hence compact. Note that we can consider larger compact sets to take into account larger initializations, the fixed point remains the same. Indeed for a larger compact C containing C, L is in the interior of C and is still a fixed point of J when the topology is restricted to C and this fixed point must be unique.

A.2 Properties of the fixed-set mapping

We now equip the set of matrices R p×(p×m) with the norm [A, B] p,m = max{ A op , B } where A ∈ R p×p and B ∈ R p×m . The set of compact subsets of R p×(p+m) is endowed with the corresponding Hausdorff distance.

Definition 1 (Affine contraction sets) For ρ ∈ [0, 1), we denote by C ρ , the set of compact subsets of matrices in R p× (p+m) such that for all S ⊂ R p×(p+m) , S ⊂ C ρ and all M ∈ S, we have A op ≤ ρ where A ∈ R p×p is the matrix made of the first p columns of M .

Given J ∈ C ρ , we denote by fix(J ) the unique fixed point of the corresponding iteration mapping as defined in Theorem 1. We have the following Proposition 5 (Monotonicity of the fixed set) Given J ∈ C ρ and J ∈ C ρ (as in Definition 1), such that J ⊂ J , we have fix(J ) ⊂ fix( J ).

Proof : Setting X 0 = fix(J ), we have

X 0 = J (X 0 ) ⊂ J (X 0 ),
and the result follows by the same argument as in the last paragraph of the proof of Theorem 1.

Proposition 6 (The fixed-set mapping is locally Lipschitz continuous) The function fix is locally Lipschitz continuous on C ρ (as in Definition 1). More precisely, for any

J 0 ∈ C ρ and J ∈ C ρ , dist (fix(J 0 ), fix(J )) ≤ 1 1 -ρ + sup [A0,B0]∈J0 B 0 (1 -ρ) 2 dist(J 0 , J )
Proof : Given J 0 ∈ C ρ and J ∈ C ρ , we remark that J ⊂ J 0 + dist(J 0 , J )B pm , where dist and B pm are considered with respect to the norm • pm . This means

J ⊂ {[A 0 , B 0 ] + [C, D], [A 0 , B 0 ] ∈ J 0 , [C, D] p,m ≤ dist(J 0 , J )}
We have

J (fix(J 0 )) = {AX + B, [A, B] ∈ J , X ∈ fix(J 0 )} ⊂ {A 0 X + B 0 , [A 0 , B 0 ] ∈ J 0 , X ∈ fix(J 0 )} + {CX + D, [C, D] mp ≤ dist(J 0 , J ), X ∈ fix(J 0 )} = J 0 (fix(J 0 )) + {CX + D, [C, D] mp ≤ dist(J 0 , J ), X ∈ fix(J 0 )} = fix(J 0 ) + {CX + D, [C, D] mp ≤ dist(J 0 , J ), X ∈ fix(J 0 )} .
This sets one inclusion. Similarly, we have

fix(J 0 ) = J 0 (fix(J 0 )) ⊂ J (fix(J 0 )) + {CX + D, [C, D] mp ≤ dist(J 0 , J ), X ∈ fix(J 0 )} . Recall that [C, D] mp = max{ C op , D }, we have for any [C, D] with [C, D] mp ≤ dist(J 0 , J ) and X ∈ fix(J 0 ), CX + D ≤ C op fix(J 0 ) sup + D ≤ dist(J 0 , J )(1 + fix(J 0 ) sup ).
We deduce using Lemma 1 that dist(fix(J 0 ), J (fix

(J 0 ))) ≤ dist(J 0 , J )(1 + fix(J 0 ) sup ).
Setting X 0 = fix(J 0 ), invoking Theorem 1 with J and k = 0, we have

dist(fix(J 0 ), fix(J )) ≤ dist(J 0 , J )(1 + fix(J 0 ) sup ) 1 -ρ ≤ dist(J 0 , J ) (1 -ρ + sup [A0,B0]∈J0 B 0 ) (1 -ρ) 2 .

A.3 Perturbed iterations

The following proposition shows that the linear convergence property is actually stable to perturbations. It will be useful to show that all potential limits of unrolling algorithmic differentiation recursions are contained in the corresponding fixed point set.

Proposition 7 (Perturbed set sequences) Let ρ < 1 and > 0 such that ρ + < 1. Let (J k ) k∈N be a sequence in C ρ+ and J ∈ C ρ (as in Definition 1). Assume that for all k ∈ N gap pm (J k , J ) ≤ or in other words J k ⊂ J + B pm where B pm is the unit ball of the norm • pm . Then the recursion on compact sets

X k+1 = J k (X k ) satisfies for all k ∈ N gap(X k , fix(J )) ≤ (ρ + ) k (1 + ρ + ) X 0 sup + sup [A,B]∈ J B + 1 -ρ - + (1 -ρ + sup [A,B]∈ J B ) (1 -ρ) 2 .
In other words, X k ⊂ fix( J ) + C(ρ, , k)B where C(ρ, , k) is the constant above.

Proof :

Set J := {J + [C, D], J ∈ J , [C, D] mp ≤ }.
Denote by ( Xk ) k∈N the sequence satisfying the recursion, Xk+1 = J ( Xk ) with X 0 = X0 . We have

X 1 = J (X 0 ) ⊂ J (X 0 ) = X1
and by recursion X k ⊂ Xk for all k ∈ N. By Theorem 1, we have

dist( Xk , fix(J )) ≤ (ρ + ) k dist(X 0 , J (X 0 )) 1 -ρ - .
We deduce from Proposition 6 that for all k ∈ N, dist( Xk , fix( J ))

≤ dist( Xk , fix(J )) + dist(fix(J ), fix( J ))

≤ (ρ + ) k dist(X 0 , J (X 0 )) 1 -ρ - + (1 -ρ + sup [A,B]∈ J B ) (1 -ρ) 2 dist(J , J ) ≤ (ρ + ) k (1 + ρ + ) X 0 sup + sup [A,B]∈ J B + 1 -ρ - + (1 -ρ + sup [A,B]∈ J B ) (1 -ρ) 2 .
And the result follows because

max X∈X k min L∈fix( J ) X -L ≤ max X∈ Xk min L∈fix( J ) X -L ≤ dist( Xk , fix( J )).

This allows to obtain explicit convergence results as follows

Corollary 4 (Limit of iterations with vanishing perturbations) Let ρ < 1 and J ∈ C ρ (as in Definition 1). Let (J k ) k∈N be a sequence of matrices such that for all k ∈ N gap pm (J k , J ) ≤ k where ( k ) k∈N is a positive sequence such that there exists a constant a > 0 such that k ≤ aρ k for all k ∈ N. Then for the recursion on compact sets of p × m matrices

X k+1 = J k (X k )
There are constants C, c > 0 such that for all k ∈ N gap(X k , fix(J )) ≤ Ce -ck .

Furthermore, one can take c = log

1 √ ρ+
for arbitrary > 0.

Proof : We consider K ∈ N such that k ≤ for all k ∈ N where + ρ < 1. Without loss of generality, we may assume that K = 0. Using the same notations as in the proof of Proposition 7, we have X k ⊂ Xk for all k ∈ N. Furthermore, it follows from the same arguments as in the proof of Theorem 1 that

X k sup ≤ Xk sup ≤ M, (15)
for a constant M > 0. Now choose k ∈ N, applying Proposition 7 shifting the initialization 0 to k, we have for all m ∈ N max

X∈X k+m min L∈fix(J ) X -L ≤ (ρ + k ) m (1 + ρ + k ) X k sup + sup [A,B]∈ J B + k 1 -ρ -k + k (1 -ρ + sup [A,B]∈ J B ) (1 -ρ) 2 ≤ (ρ + ) m (1 + ρ + )M + sup [A,B]∈ J B + 1 -ρ - + aρ k (1 -ρ + sup [A,B]∈ J B ) (1 -ρ) 2 ,
where we have used the bound [START_REF] Christianson | Reverse accumulation and attractive fixed points[END_REF] and the fact that k ≤ and k ≤ aρ k . Setting

u = (1+ρ+ )M +sup [A,B]∈ J B + 1-ρ- and v = a (1-ρ+sup [A,B]∈ J B ) (1-ρ) 2
we have

max X∈X 2k min L∈fix(J ) X -L ≤ u(ρ + ) k + vρ k ≤ (u + v)(ρ + ) 2k/2 ≤ u + v (ρ + ) 1/2 (ρ + ) 2k/2 , max X∈X 2k+1 min L∈fix(J ) X -L ≤ u(ρ + ) k+1 + vρ k ≤ u + v (ρ + ) 1/2 (ρ + ) (2k+1)/2 .
Since k was arbitrary, this proves the desired result.

B Existence of a conservative Jacobian for autodiff

B.1 Regularity of J pb x
We recall the main notations and elements of Assumption 1. We assume that F is locally Lipschitz, path differentiable, and denote by J F : R p+m ⇒ R p×(p+m) a conservative Jacobian for F . Now assume that any pair [A, B] ∈ J F (x, θ) is such that the operator norm of A is at most ρ < 1, that is for all x and θ, J F (x, θ) ∈ C ρ (as in Definition 1). Define the following set-valued map

J pb x : θ ⇒ fix [J F (x(θ), θ)] .
Here, x(θ) = fix(F θ ) is the unique fixed point of the algorithmic recursion so that we actually have

J pb x : θ ⇒ fix [J F (fix(F θ ), θ)] .
We have the following Lemma 2 (Regularity of J pb

x ) The mapping J pb x is nonempty valued, locally bounded and has a closed graph.

Proof : The fact that J pb

x is locally bounded and non empty valued comes from the fact that J F is locally bounded with nonempty values and x is locally Lipschitz combined with Theorem 1. By local Lipschitz continuity of x and the fact that J F has a closed graph, the set-valued map θ ⇒ J F (x(θ), θ) also has a closed graph. By continuity of fix(J ) with respect to the Hausdorff distance, see Proposition 6, J pb

x has a closed graph.

B.2 Proof of Theorem 2

Proof : Following Remark 2, we set

J imp x : θ ⇒ (I -A) -1 B, [A, B] ∈ J F (x(θ), θ) ,
a conservative Jacobian for x and L 0 = J imp x . Now set by recursion for all k ∈ N

L k+1 : θ ⇒ J F (x(θ), θ)(L k (θ)).
Recall that this means for all θ ∈ R m and k ∈ N

L k+1 (θ) = {AL + B, [A, B] ∈ J F (x(θ), θ), L ∈ L k (θ)}.
Since F (x(θ), θ) = x(θ) for all θ, J F is conservative for F and L 0 is conservative for x, we have by induction that for all k ∈ N, L k is conservative for x. Fix l : R m → R m an arbitrary Borel measurable selection in J pb x , that is l(θ) ∈ J pb x (θ) for all θ ∈ R m . Such a selection exist by [START_REF] Charalambos | Infinite Dimensional Analysis: A Hitchhiker's Guide[END_REF]Theorem 18.20] because J pb

x has a closed graph by Lemma 2. Set for all k ∈ N a measurable selection

l k : θ → arg min z∈L k (θ) z -l(θ) .
The function (z, θ) → zl(θ) is Caratheodory (continuous in z, measurable in θ), so such a selection exists (Aliprantis Theorem 18.19). By Theorem 1, we have that dist(L k (θ), J pb x (θ)) tends to 0 as k grows, for all θ ∈ R m , where the convergence is in Hausdorff distance. Actually since all set-valued objects are locally bounded, the convergence occurs uniformly on every compact. This implies in particular that l k converges pointwise to l. Fix an absolutely continuous path γ : [0, 1] → R m . We have for all k ∈ N, by conservativity,

x(γ(1)) -x(γ(0)) = 1 0 l k (γ(t)) γ(t)dt.
Furthermore, l k • γ is measurable, converges pointwise to l • γ and l k • γ can be uniformly bounded, let K be such a bound. The integrable function g : t → K γ(t) dominates the integrand and l k • γ × γ converges pointwise to l • γ × γ. By the dominated convergence theorem (see [START_REF] Halsey | Real analysis[END_REF]Section 4.4] ), we have

x(γ(1)) -x(γ(0)) = 1 0 l(γ(t)) γ(t)dt.
L has a Castaing representation with a dense sequence of measurable selection [START_REF] Charalambos | Infinite Dimensional Analysis: A Hitchhiker's Guide[END_REF]Theorem 18.14]. Since l was an arbitrary measurable selection in L, conservativity of L follows by [START_REF] Marx | Path differentiability of ode flows[END_REF]Lemma 8].

B.3 Proof of Corollary 1

Proof : Fix θ. We have x k (θ) → x(θ), so that for any > 0, there exists K ∈ N such that J F (x k (θ), θ) ⊂ J F (x(θ), θ) + B for all k ≥ K. The result is then a consequence of Proposition 7, letting → 0. The last part is due to the conservativity of J pb

x which must be a singleton almost everywhere, equal to the classical Jacobian.

B.4 Proof of Corollary 3

Proof : Define (L k ) k∈N , a sequence of conservative Jacobians for x as in the begining of the proof of Theorem 2 in Appendix B.2. By [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF]Theorem 1], for each k ∈ N, there is a full measure set S k ⊂ R m such that L k (θ) = ∂ x ∂θ (θ) for all θ ∈ S k . Similarly, there exists a full measure set S -1 ⊂ R m such that J pb x (θ) = ∂ x ∂θ (θ) for all θ ∈ S -1 . Setting S = ∩ +∞ i=-1 S i , S has full measure and for all θ ∈ S and for all k ∈ N,

J pb x (θ) = ∂ x ∂θ (θ) L k (θ) = ∂ x ∂θ (θ) .
Following the proof of Theorem 2 in Appendix B.2, L k converges to J pb x in Hausdorff distance, which means that convergence occurs in the classical sense since all sets in the sequence are singletons.

B.5 Proof of Proposition 1

Proof : Under the setting of Corollary 2, for almost all θ ∈ R m , recursion (PB) or ( 5) reduce to the following, and all k ∈ N

J k+1 = A k J k + B k (16)
where

J k = ∂x k ∂θ , A k = ∂F ∂x (x k , θ
) and B k = ∂F ∂θ (x k , θ) are classical Jacobians and J k converges to the classical Jacobian of ∂ x ∂θ (θ). Fix such a θ ∈ R m and k ∈ N, k ≥ 1. With the notation of Algorithm 1, for the forward mode, multiplying [START_REF] Frank | Optimization and nonsmooth analysis[END_REF] on the right by θ, we have for all i ∈ 1, . . . k

J i θ = A i-1 J i-1 θ + B i-1 θ.
Setting ẋi = J i θ, this is exactly the recursion implemented by Algorithm 1 in forward mode. Corollary 2 and the result follows from convergence of J k . As for the backward mode a simple recursion shows that

J k = A k-1 A k-2 . . . A 0 J 0 + A k-1 A k-2 . . . A 1 B 0 + . . . + A k-1 A k-2 . . . A i B i-1 + . . . + A k-1 B k-2 + B k-1 . (17)
Setting B -1 = J 0 , we may rewrite equivalently,

J k = B k-1 + k-1 i=0   i j=k-1 A j   B i-1 . (18)
Transposing and multiplying on the right by wk , we have

J T k wk = B T k-1 wk + k-1 i=0 B T i-1   k-1 j=i A T j   wk . (19)
We set for all i = 0, . . . , k -1,

wi = k-1 j=i A T j wk . (20)
We have the backward recursion relation, for i = k, . . . , 1 wi-1 = A T i-1 wi , which is the recursion implemented by Algorithm 1 in reverse mode. Combining [START_REF] Efron | Least angle regression[END_REF] and [START_REF] Ghaoui | Implicit deep learning[END_REF], we obtain

J T k wk = B T k-1 wk + k-1 i=0 B i-1 wT = k i=1 B T i-1 wi + J T 0 w0 ,
which is the quantity accumulated in θk in Algorithm 1. This proves that θT k returned by the backward mode is indeed equal to wT k J k and the convergence follows from convergence of both wk and J k as k → ∞.

C Connection with implicit differentiation

Recall that for all θ

J imp x (θ) = (I -A) -1 B, [A, B] ∈ J F (x(θ), θ) = {M, ∃[A, B] ∈ J F (x(θ), θ) M = AM + B} .
Setting J = J F (x(θ), θ), we have therefore that J imp x (θ) ⊂ J (J imp x (θ)). By recursion, for all k ∈ N, J imp x (θ) ⊂ J k (J imp x (θ)) and passing to the limit using Theorem 1, J imp x (θ) ⊂ fix(J ) = J pb x (θ). In particular, if F is continuously differentiable, then (PB) with classical Jacobians converges towards a classical implicit derivative. However, the inclusion J imp x (θ) ⊂ J pb x (θ) may be strict as the following example shows.

Example 1 Set J = {[A, B], A ∈ A, B ∈ B}, where A = λ+1 4 0 0 2-λ 4 , λ ∈ [0, 1] B = 1 1 .
We set

T = (I -A) -1 B = 4 3-λ 4 2+λ , λ ∈ [0, 1] .
As already observed, we have T ⊂ AT + B, but the inclusion is strict. Therefore T is not a fixed point of the affine iteration and it is only contained in it. Indeed, we have

1+1 4 0 0 2-1 4 4 3-0 4 2+0 + 1 1 = 5 3 3 2 ∈ AT + B.
However solving for λ

(J s F (x), J s F (x 0 )) ≤ L x -x 0 , ∀x, x -x 0 ≤ R,
where L is the Lipschitz constant given by the selection structure of F .

Proof : Fix x 0 ∈ R p and consider the function g which associates to r > 0 a subset of {1, . . . , m} defined as

g(r) = ∪ x-x0 ≤r I(x).
The function g is semialgebraic and therefore it admits a limit as r → 0. The function g is actually piecewise constant so that the limit is reached for some R > 0 by semialgebraicity. This means that there is R > 0 and an index set I ⊂ {1, . . . , m} such that I(x) ⊂ I for all x such that xx 0 ≤ R. Furthermore, for each i ∈ I and all 0 < r ≤ R, there exists x such that xx 0 ≤ r and F i (x) = F (x). By continuity of each component F i , we have for each i ∈ I, F i (x 0 ) = F (x 0 ), that is I ⊂ I(x 0 ). We deduce that for each x such that xx 0 ≤ R and i ∈ I(x), we have

min V ∈J s F (x0) V - ∂F i ∂x (x) ≤ ∂F i ∂x (x 0 ) - ∂F i ∂x (x) ≤ L x -x 0 .
Fix any Z ∈ J s F (x), it is a convex combination of ∂Fi ∂x (x) for i ∈ I(x) so by convexity of the distance, we have min

V ∈J s F (x0) V -Z ≤ L x -x 0 ,
which proves the result since this allows to bound the supremum over Z ∈ J s F (x) by the desired quantity.

D.2 Proof of Corollary 3

Proof : This is a consequence of linear convergence of the recursion x k+1 = F (x k , θ) combined with Lemma 3 and Corollary 4.

E Proximal splitting algortihms in convex optimization

E.1 Proof of Proposition 2

Proof : We consider the gradient step operation H α : (x, θ) → x -α∇ x f (x, θ). We have for all (x, θ),

F α (x, θ) = G α (H α (x, θ), θ).
By Assumption 2, both G α and H α are 1-Lipschitz in x for fixed θ and we are going to show that if either f or g satisfy the strong convexity condition, the corresponding map is a strict contraction in x for fixed θ. Furthermore, the mapping Jac c Hα : (

x, θ) ⇒ [I -αA, -αB], [A, B] ∈ J 2 f (x, θ)
is the Clarke Jacobian of H α . By Assumption 2, all the functions are path-differentiable [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF] and one may obtain a conservative jacobian for F by applying differential calculus rules [START_REF] Bolte | Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning[END_REF]. We set for all (x, θ) a conservative Jacobian for F α ,

J Fα (x, θ) = [C(I -αA), -αCB + D], [A, B] ∈ J 2 f (x, θ), [C, D] ∈ J Gα (x -α∇ x f (x, θ), θ) (21)
Whenever ∇ x f is differentiable at (x, θ), the first p columns of its Jacobian form a symmetric positive definite square matrix with eigenvalues at most L. This implies that the matrix (I -αA) in ( 21) is symmetirc with eignevalues in [-1, 1] and strictly greater than -1.

Similarly, whenever G α is differentiable, since it is 1-Lipschitz in x for fixed θ and the gradient of a C 1 function, the first p columns of its Jacobian form a symmetric positive definite square matrix with eigenvalues at most 1. This implies that the matrix C in ( 21) is symmetric with eignevalues in [0, 1]. In addition, we have the following;

• Assume that for all θ, f is µ-strongly convex. In this case, similarly as above the matrix (I -αA) in ( 21) has eigenvalue in (-1, 1) for all (x, θ).

• Assume that for all θ, g is µ-strongly convex. In this case, similarly as above the matrix C in [START_REF] Friedman | Sparse inverse covariance estimation with the graphical lasso[END_REF] has eigenvalue in [0, 1/(1 + αµ)] for all (x, θ) [START_REF] Heinz H Bauschke | Convex analysis and monotone operator theory in Hilbert spaces[END_REF]Proposition 23.13].

In both cases, the product C(I -αA) in ( 21) has operator norm strictly smaller than 1 and Assumption 1 holds.

E.2 Proof of Proposition 3

Proof : From [6, Proposition 23.11], both R αf and R αg are 1-Lipschitz. We are going to show that R αf is a strict contraction and the result will follow. Since f is C 1,1 in x, we have for all θ ∈ R m , z = prox αf (•,θ) (x) ⇔ z + α∇ x f (z, θ)x = 0

Set H α (z, x, θ) = z + α∇ x f (z, θ)x, we have that Jac c Hα (z, x, θ) ⇒ {[I + αA, -I, αB]}

is the Clarke Jacobian of H α . Similarly as in Appendix E.1, by strong convexity of f , the matrix I + αA in ( 22) is symmetric with eigenvalues strictly greater than 0 and smaller than 1. By implicit differential calculus rule in [10, Theorem 2], the mapping

J prox αf (•,θ) (x, θ) ⇒ [(I + αA) -1 , -α(I + αA) -1 B], [A, B] ∈ J 2 f (prox αf (•,θ) , θ) (23)
is conservative for (x, θ) → prox αf (•,θ) . Furthermore, the matrix (I + αA) -1 in ( 23) is symmetric eigenvalues in (0, 1). This entails that the mapping

is conservative for R αf (•,θ) and the matrix 2(I + αA) -1 -I is symmetric with eigenvalues in (-1, 1).

Similarly, the mapping

J R αg(•,θ) (x, θ) ⇒ [2C -I, 2D -I], [C, D] ∈ J prox αg(x,θ) (25)
is the Clarke Jacobian of R αg(•,θ) and the matrix 2C -I in ( 25) is symmetric with eigenvalues in [-1, 1]. One may combine J R αf (•,θ) and J R αg•,θ) , using differential calculus rule to obtain a conservative Jacobian J Fα for F α , such that for all (x, θ) and [E, F ] ∈ J Fα (x, θ), the square matrix E is of the form I 2 + ((I + αA) -1 -I)(2C -I) where A is from [START_REF] Charles | Automatic differentiation and iterative processes[END_REF] and C is from [START_REF] Giselsson | Linear convergence and metric selection for douglasrachford splitting and admm[END_REF]. Such a matrix E has operator norm strictly smaller than 1 which is Assumption 1.

E.3 Equivalence between ADMM and dual Douglas-Rachford

We need the following lemma. Then, h is convex, lower semicontinuous, closed, and

prox αh (x) = x + α(Aû -v) (26)
where (û, v) ∈ arg min u,v

F (u) + G(v) + x (Au -v) + α 2 Au -v 2 2 .
The material contained in this section is already known in the litterature accross several papers and lecture notes, but for the sake of completeness, we include a full derivation of the equivalence.

In this appendix, we drop the dependency to the variable θ since we are only concerned on the behaviour with respect to x. We recall that the iteration of Douglas-Rachford are defined by an initialization y 0 and the recursion x k+1 = prox f (y k )

y k+1 = y k + prox g (2x k+1 -y k ) -x k+1 . (27)
By denoting xk = x k+1 and ỹk = y k , we can rewrite the updates of Douglas-Rachford (given x0 and ỹ0 ) as

ỹk+1 = ỹk + prox g (2x k -ỹk ) -xk . xk+1 = prox f (ỹ k+1 ) (28)
Introducing the variable r = prox g (2xŷ), this is also equivalent to

This formulation will be convenient to show how to retrieve the equations of ADMM [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF].

The dual problem of ( 12) is given by ( 14)

max x -f (x) -g(x). ( 31
)
where f (x) = φ (-Ax) + c x and g(x) = ψ(-Bx)

We consider the update rules given by (30), i.e., r = prox αg (x + w)

x = prox αf (rw) Finaly, combining the expression of r and x, we obtain ŵ = αBv.

Figure 2 :

 2 Figure 2: Illustration of the linear convergence of proximal splitting methods. (First line) Distance of the iterates to the fixed point. (Second line) Distance of the piggyback Jacobians to the Jacobian of the fixed point. The acronyms are FB for Forward-Backward, DR for Douglas-Rachford and ADMM for Alternating Direction Method of Multipliers. In all cases, despite nonsmoothness, piggyback Jacobians converge, illustrating Corollary 2. Blue lines represent the median of 100 repetitions with random data, and the blue shaded area represents the first and last deciles.

Figure 3 :

 3 Figure 3: Behaviour of automatic differentiation for first-order methods on a quadratic function. (Left) Stability of the propagation of derivatives for the fixed step-size gradient descent. (Right) Instability of the propagation of Heavy-Ball initialized. Both methods are initialized at optimum.

,

  the first equation requires λ =3 5 while the second requires λ = 2 3 which shows that the given vector does not belong to T .D Semialgebraic Lipschitz gradient selection functionsD.1 Lipschitz property of conservative Jacobians of selectionsLemma 3 (Conservative Jacobians of selections are Lipschitz-like) Let F be continuous, semialgebraic with Lipschitz gradient selection. Then for each x 0 ∈ R p , there exists R > 0 such that gap

J

  R αf (•,θ) (x, θ) ⇒ [2(I + αA) -1 -I, -2α(I + αA) -1 B -I], [A, B] ∈ J 2 f (prox αf (•,θ) , θ)

Lemma 4

 4 Let F, G two convex, lower semicontinuous and closed functions and h defined byh(x) = F * (-A x) + G * (x).

  rk+1 = prox g (2x kŷk ) xk+1 = prox f (ŷ k + rk+1xk ) ŷk+1 = ŷk + rk+1xk (29)Using the change of variable ŵk = xkŷk , we haverk+1 = prox g (x k + ŵk ) xk+1 = prox f (r k+1ŵk ) ŵk+1 = ŵk + xk+1 -rk+1 .

2 Au -c + w/α 2 2 .

 22 )ŵ = w + xr.(34)Applying Lemma 4 to F = φ and G = ι c , we rewrite (32) byr = x + w + α(Aûc) where û = arg min u φ(x) + x (Auv) + αUsing the same lemma to F = ψ and G = 0, we rewrite (33) byx = x + α(Aû + Bvc) where v = arg min v ψ(v) + x Bv + α 2 Aû + Bvc .

Acknowledgments J. B. and E. P. acknowledge the financial support of the AI Interdisciplinary Institute ANITI funding under the grant agreement ANR-19-PI3A-0004, Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant numbers FA9550-19-1-7026, and ANR MaSDOL 19-CE23-0017-01. J. B. also acknowledges the support of ANR Chess, grant ANR-17-EURE-0010, TSE-P and the Centre Lagrange. S. V. acknowledges the support ANR GraVa, grant ANR-18-CE40-0005.

F Inertial methods

Let us first recall notations from Section 5. Consider a function f : R p × R m → R, and β > 0, for simplicity, when the second argument is fixed we write f θ : x → f (x, θ). Set for all x, y, θ, F (x, y, θ) = (x -∇f θ (x) + β(xy), x), consider the Heavy-Ball algorithm (x k+1 , y k+1 ) = F (x k , y k , θ) for k ∈ N. If f θ is µ-strongly convex with L-Lipschitz gradient, then, choosing α = 1/L and β < 1 The Jacobian of F for the Heavy-Ball agorithm (in x, y) is of the form

, then the Hessian can be replaced by a set-valued conservative Jacobian of the gradient: J ∇f θ .

Proof of Proposition 4:

Recall that the function f : R 2 → R is given by

We have f (x) = x for t ≥ 0 and f (x) = x 4 for t < 0, therefore, f is 1-Lipschitz. The Clarke subdifferential of f is { 1 4 } for t < 0, {1} for t > 0 and the segment 1 4 , 1 at t = 0. Finally, f is µ = 1 4 strongly convex and has L = 1 Lipschitz gradient and the unique fixed point of the Heavy-Ball algorithm applied to f (•, θ) is x = y = θ. Choosing α = 1, β = 0.75, we have

Therefore, the heavy ball algorithm with this choice of parameter converges linearly to the unique solution which is 0, a fixed point of the iteration mapping. Set

At (0, 0, 0), the last column of the Jacobian of F is (0, 0) and the first two columns are given by

where

Therefore, the Clarke Jacobian of F (with respect to x, y) at (0, 0, 0) is given by

We have

which has two eigenvalues -9 8 < -1 and -9 32 . Setting for any θ ∈ R x 0 (θ) = θ, y 0 (θ) = θ, we have for all k ∈ N x k (θ) = y k (θ) = θ, in other words, this is the unique fixed point of the Heavy-Ball algorithm.

Given l ∈ N, the forward propagation recursion in (PB) presented in Figure 3 satisfies for k = 8l

This products will diverge diverge due to the eigenvalue of (M 1 M 1 M 2 M 2 ) 2 strictly above 1.

In other words, for all k, J x 8k given by (PB) contains elements which magnitude diverges at a geometric rate. We conclude that, for all k ∈ N, J x k contains elements which magnitude diverge at a geometric rate. This illustrates the failure of forward derivative propagation on f (•, θ): the Heavy Ball algorithm is stable and globally linearly convergent, its fixed point is differentiable (it is actually constant in θ), yet there is a parametric initialization x(θ), y(θ) such that forward propagation of derivatives produces diverging elements for θ = 0. Note that implicit differentiation provides the correct derivative, which is 0, since x(θ) = 0 is the unique fixed point of the gradient iterations. Forward derivative propagation on the gradient descent algorithms also results in the limit in 0 derivative since it only contains element which converge to 0 at a geometric rate.

Le us emphasize again that such pathology would not happen if f was C 2 . Indeed, in this case, J 2 f would be single valued and the divergence phenomenon would not appear. This illustrate a fundamental difference between C 1,1 and C 2 objectives in terms of forward derivative propagation for second order inertial methods.

G Experiments details

All the experiments where run on a MacBook M1 Pro (arm64), on Python 3.9 and numpy 1.21 for a compute time inferior to one hour. They are repeated 100 times, and we report the median as a blue line and the first and last deciles as a blue shaded area. The solutions are computed with 2000 iterations, and the curves are reported for the 1000 first iterations. The differentiation of all methods is performed in forward-mode with jacfwd of the module jax.

Forward-Backward for the Ridge. The dimensions of the problem are n = 500, p = 300. The design matrix is Gaussian, i.e., X i,j

H Assets used

Our numerical experiments rely on:

• numpy [START_REF] Harris | Array programming with NumPy[END_REF], released under BSD-3 license.

• matplotlib [START_REF] Hunter | Matplotlib: A 2d graphics environment[END_REF], released under PSF license.

• jax [START_REF] Bradbury | JAX: composable transformations of Python+NumPy programs[END_REF], released under Apache-2.0 license.


Imane Lboukili

Xavier Descombes

Georgios Stamatas

Automatic cell identification and analysis on in vivo reflectance confocal microscopy images of the human epidermis

Keywords: epidermis, image segmentation, reflectance confocal microscopy, object detection

Reflectance confocal microscopy (RCM) allows real-time in vivo visualization of the skin at cellular level. The study of RCM images provides information on the topological and geometrical properties of the epidermis. These may change in each layer of the epidermis, depending on the subject's age and the presence of certain dermatological conditions. Studying RCM images requires manual identification of cells to derive these properties which is timeconsuming and subject to human error, highlighting the need for an automated cell identification method. We propose an automated pipeline to analyze the structure of the skin in RCM images. The first step is to identify the region of interest (ROI) containing the epidermal cells. The second step is to identify individual cells in the segmented tissue area using an image filter. We then use prior biological knowledge to process the resulting detected cells, removing cells that are too small and reapplying the used filter locally on detected regions that are too big to be considered as a single cell. The results are evaluated both on simulated data and on manually annotated real RCM data. This study shows that automatic cell identification can be achieved, with an accuracy (precision and recall) that matches the interexpert variability.

INTRODUCTION

Reflectance Confocal Microscopy 1,2 (RCM) is a real-time non-invasive in vivo technology that allows the visualization of the skin epidermis and upper layers of the dermis at cellular level. It is non-invasive thus making it a technique of choice for repeated sampling on a skin site without damage, when studying the changes in skin structure over-time, or when an invasive biopsy cannot be considered, e.g., cosmetology, study of baby skin, etc. It is made of the collection of signals arising from light reflections at the interface of microstructures with different indices of refraction. In skin, such microstructures are keratin fibers, melanosomes, collagen fibers and cell membranes. Therefore, it provides information on the geometrical and topological properties of the skin, which play important roles in the makeup of the skin barrier.

In most cases, analysis of RCM stacks is done manually, providing qualitative observations. However, manual analysis is time-consuming, intensive, and subject to human interpretation and inter-expert differences. Thus, we could benefit from automated methods to quantitatively analyze RCM images. An important first step in any quantitative study of skin is cell detection. Unfortunately, it is challenging and requires a robust generic algorithm to alleviate images nonuniformity and noise.

The epidermis is made of four distinct layers. From deepest to superficial, they are stratum basale (SB), stratum spinosum (SS), stratum granulosum (SG), and stratum corneum (SC).

In RCM images of light-pigmented skin, the SC appears as large bright islands surrounded by dark grooves. It is made of dead but biochemically active cells [START_REF] Murphrey | Histology, Stratum Corneum[END_REF] . As we cannot observe individual cells on RCM images of the SC, our method will not be applied to these images.

The SG and SS appear as islands of viable cells, called keratinocytes, arranged in a honeycomb pattern. Granular cells are typically larger than spinous cells and as such they have a higher density [START_REF] Bensaci | Geometrical and topological analysis of in vivo confocal microscopy images reveals dynamic maturation of epidermal structures during the first years of life[END_REF] .

Finally, the SB is made of the smallest keratinocytes, as their differentiation starts in the SB and continues as cells migrate towards the skin surface getting bigger and flatter. The SB is attached to the dermis on the dermal-epidermal junction, and thus we can sometimes observe the top of dermal papillae on RCM images of the SB. Additionally, melanin-producing melanocytes are scattered through the basal layer. Organelles filled with melanin, called melanosomes, are transferred from melanocytes to keratinocytes. Melanin has a high reflectance [START_REF] Rajadhyaksha | In vivo confocal scanning laser microscopy of human skin: melanin provides strong contrast[END_REF][START_REF] Sánchez-Mateos | Reflectance-Mode Confocal Microscopy in Dermatological Oncology[END_REF] thus making RCM images of the SB more noisy and of lower quality than images of the SG and SS. For these reasons, we focus on the automated detection of keratinocytes in the granular and spinous layers.

On RCM images of the SG and SS, keratinocytes are characterized by a dark center and a grainy cytoplasm due to micro-structures, surrounded by bright grainy membranes (see Figure 1).

Attempts at the automated identification of epidermal cells on RCM images have been made [START_REF] Harris | A Pulse Coupled Neural Network Segmentation Algorithm for Reflectance Confocal Images of Epithelial Tissue[END_REF][START_REF] Gareau | Automated identification of epidermal keratinocytes in reflectance confocal microscopy[END_REF] . Unfortunately, the amount of noise and heterogeneity of RCM images hinders the development of accurate segmentation methods. We propose a method to automatically detect keratinocytes centers on RCM images of the SG and SS. We compare our results to a manually obtained ground truth and obtain an accuracy on par with cross-experts identification.

Figure 1. RCM image of the stratum spinosum of light-pigmented skin. In blue, the border between tissue and background. In pink, non-informative areas. In orange, bright spots. And in red, epidermal cells.

DATA

In vivo RCM images were captured on the volar forearm of 80 participants: 60 children aged 3 months to 10 years old and 20 adults aged 25 to 40 years old. All participants have light-pigmented skin, with Fitzpatrick types between I and IV. Inclusion criteria required that the participants were in good health, with no history of skin disease and had not applied any products on the observed area the day of the study.

Images were captured using a Vivascope 1500 reflectance confocal microscope with a z-resolution of 5 micrometers and xy-resolution of 1 micrometer. Images started at the SC and progressed down towards the SB. The size of each image is 1000*1000 pixels.

Image classification in one of the four epidermal layers was obtained using a hybrid deep learning algorthim [START_REF] Kaur | Hybrid deep learning for Reflectance Confocal Microscopy skin images[END_REF] , allowing us to focus only images of the stratum granulosum and stratum spinosum.

IMAGE ANALYSIS WORKFLOW 3.1. Identification of the region of interest

RCM images tend to be noisy and non-uniform, which hinders the development of automated segmentation methods. To guide our cell detection, we started by identifying the region of interest (ROI), i.e., the regions containing epidermal cells. To do so, the black background was first identified (see Figure 2). Indeed, islands of cells surrounded by dark empty areas are observed on RCM images. They are due to grooves on the skin surface called micro-relief lines [START_REF] Rajadhyaksha | In vivo confocal scanning laser microscopy of human skin II: advances in instrumentation and comparison with histology[END_REF] (see Figure 1). In order to identify these furrows and begin building a binary mask of the ROI, a Morphological Geodesic Active Contour 11 algorithm (so called Snake) was applied to each image. This method employs morphological operators to detect visible contours based on their intrinsic geometric measures, even if they are noisy or partially unclear.

After identifying the micro-relief lines on each RCM image, non-informative areas were detected. These a reas can be seen in the tissue on RCM images and are due to low contrast and a drop in signal-to-noise ratio (see Figure 1). A texture classification was applied to the images, by training a Support Vector Machine (see Figure 2) on four features of the Grey Level Co-occurrence Matrix [START_REF] Haralick | Textural Features for Image Classification[END_REF] (GLCM) which successfully discriminates between the informative and non-informative areas. These features are, (a) homogeneity which measures the closeness of the GLCM distribution to its diagonal (reflecting correlation), (b) contrast which measures the local variations in the GLCM; (c) dissimilarity which measures the similarity between pixels, and (d) energy which measures the uniformity of the area.

Our third step in ROI identification was to remove the bright spots sometimes observable in RCM images (Figure 2). Indeed, RCM images of the SG and SS may contain bright areas due to the presence of keratin in hair shafts, from cornified cells at the periphery of the islands (see Figure 1). This was accomplished by applying a succession of dilations and erosions on the RCM image where the background and the non-informative areas were removed, and which had been blurred with a Gaussian filter and binarized.

Figure 2. Steps of the identification of the region of interest. A morphological snake was used to identify the borders with the background, followed by a Support Vector Machine algorithm to detect the non-informative areas, and a succession of morphological operations to remove bright spots to an RCM image at stratum granulosum level. In blue the region of interest mask.

Identification of individual cells

After identifying the ROI on the RCM image, a median filter was used to remove noise, followed by a local normalization which renders the variance and mean of the denoised image stationary (see Figure 3). Then, the resulting image was filtered with Sato tubeness filter [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF] to detect white continuous ridges, here the bright cell membranes (see Figure 3). The filter parameters were chosen to approximate the width and length of a cell membrane in the SG and SS. To the filter output, a median filter and a local normalization were applied, while making sure that the ROI binary mask is respected (see Figure 3). A Gabor filter was then applied to the previous image to refine membrane detection by convoluting the image by a sinusoidal signal of varying frequencies and orientations (see Figure 3). The output of the Gabor filter was equalized with a histogram equalization to adjust the image contrast, followed by a Gaussian adaptive thresholding (see Figure 3) which dynamically and locally changes the binarization threshold over the entire image to account for changes in contrast and brightness, assuming that smaller regions of an RCM image are more likely to be similar. A connected-components analysis was used on the obtained binary image to remove any small blobs in the detected membranes, followed by a second connected-components analysis on the inverse of the image to close any holes in the membranes due to the graininess of the image and of the cell membranes (see Figure 3). Finally, the clean binary image was skeletonized, and any spurious branches were removed from the skeleton (see Figure 3). A median filter and a local normalization were applied to the image with the ROI mask, followed by a Sato filter. Its output was filtered with a median filtered and locally normalized, and a Gabor filter was applied to it. A threshold was applied on the output after histogram equalization and small blobs were removed with a connected components analysis. The result was then skeletonized, and spurious branches were removed. ROI: region of interest.

Post-processing

After obtaining the skeleton, the ROI mask was cleaned to remove any remaining spurs. This was accomplished by applying a morphological closing to the skeleton. Individual contours, i.e., detected keratinocytes, were detected on the skeleton. To improve the detection, very small contours were removed (𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 < 100 for SG and 𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 < 50 for SS) as well as long contours at the border with the background, i.e., 𝑎𝑎𝑒𝑒𝑒𝑒𝑎𝑎𝑒𝑒𝑒𝑒𝑎𝑎𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 > 0.85 (see Figure 4 A). These thresholds were determined empirically. The remaining contours were divided into two groups, (1) large contours with an 𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 > 1000 for SG and 𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 > 120 for SS, and (2) small contours with an 𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 ≤ 1000 for SG and 𝑎𝑎𝑎𝑎𝑎𝑎𝑎𝑎 ≤ 120 for SS. On each area of the original image determined by a large contour, a Gabor filter followed by a Sato filter [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF] were applied with different parameters than previously used. The output was then binarized with Otsu thresholding [START_REF] Otsu | A Threshold Selection Method from Gray-Level Histograms[END_REF] for SG images and with Gaussian adaptive thresholding for SS images, and small blobs were removed with a connected-components analysis. The subsequent binary image was skeletonized, and its contours detected. On images of the SG, obtained contours with an area smaller than 110 were merged with their neighbors (see Figure 4 B). On images of the SS, where cells, and therefore detected contours, are much smaller, if the second filter iteration still failed to detect more than one contour, as many ellipses as possible were fitted within the detected contour, with these ellipses having the median minor and major axis length of the small contours (see Figure 4 C). These new contours were then combined with the previously found small ones and their cell centers were detected.

Accuracy evaluation

The obtained cell centers were used to initiate a marker-controlled watershed [START_REF] Meyer | Morphological segmentation[END_REF] on the ROI. This method considers the input image as a topographic surface which is flooded starting from set seeds or markers, i.e., the detected cell centers, and returns a labelled gray-scale image, where each label is a catching basin, i.e., a detected cell. This la beled image was then compared against manually detected cell centers using the software d-accuracy [START_REF] Debreuve | DAccuracy[END_REF] which evaluates several indexes of the detection quality (see Figure 5). Two accuracy metrics were evaluated, (a) precision, (fraction of correctly detected cells among all the detected cells) and (b) recall (the fraction of accurately detected keratinocytes among all cells defined in the ground truth).

Figure 5. Detection accuracy evaluation steps. A marker-controlled watershed was applied to the detected cell centers and the resulting labels were compared to the manually detected ground truth, in pink. The returned metrics were precision and accuracy.

RESULTS

We evaluate the performance of our approach with respect to two experts. When compared to the first expert, our cell detection approach has a precision of 77% and of 85% on SG and SS images respectively and a recall of 85% and 67% on SG and SS confocal images respectively. And when compared to the second expert, our cell detection approach has a precision of 64% and of 75% on SG and SS images respectively and a recall of 82% and 30% on SG and SS confocal images respectively (see Table 1).

When looking into the differences between the experts, we notice than expert 2 is less sensitive in his detection especially in SS images. The obtained results are more consistent with expert 1 and prove to be accurate compared to inter-expert variability.

DISCUSSION

Our approach gives reasonable results on RCM images of the granular and spinous layers. Performance can be hindered by the presence of multiple layers per image, which makes parametrization of the different steps complicated, although, parameters for a given layer are the same. Our method is a multi-step approach, with multiple parameters each one influencing cell detection and its accuracy. Additionally, the noise and non-uniformity of the images has a great impact on the method performance. Steps like median filtering, local normalization and ROI determination decrease the impact of noise on the results but do not remove it completely.

Computational time is about 10mn depending on the size of the ROI, which is a major advantage when compared to the time required to identify keratinocytes manually on RCM images.

We have shown that the automated detection of keratinocytes on RCM images of the SG and SS is achievable, which is an important step towards the quantitative study of these images and of skin. This prevents from tedious manual work, which is the current practice in RCM images analysis.

The presented approach on confocal images will be useful in uncovering new insights in the study of skin physiology, skin maturation and skin ageing [START_REF] Bensaci | Geometrical and topological analysis of in vivo confocal microscopy images reveals dynamic maturation of epidermal structures during the first years of life[END_REF][START_REF] Stamatas | Infant skin physiology and development during the first years of life: a review of recent findings based on in vivo studies[END_REF][START_REF] Stamatas | Infant skin maturation: Structural changes revealed by in vivo reflect[END_REF][START_REF] Cinotti | Structural skin changes in elderly people investigated by reflectance confocal microscopy[END_REF][START_REF] Kawasaki | Age-related morphometric changes of inner structures of the skin assessed by in vivo reflectance confocal microscopy[END_REF] and skin diseases observable with RCM, e.g., melanomas [START_REF] Gareau | Automated detection of malignant features in confocal microscopy on superficial spreading melanoma versus nevi[END_REF][START_REF] Bozkurt | A Multiresolution Convolutional Neural Network with Partial Label Training for Annotating Reflectance Confocal Microscopy Images of Skin[END_REF][START_REF] Zorgui | A Convolutional Neural Network for Lentigo Diagnosis[END_REF][START_REF] Halimi | An unsupervised Bayesian approach for the joint reconstruction and classification of cutaneous reflectance confocal microscopy images[END_REF][START_REF] Kose | A machine learning method for identifying morphological patterns in reflectance confocal microscopy mosaics of melanocytic skin lesions in-vivo[END_REF][START_REF] Koller | In vivo reflectance confocal microscopy: automated diagnostic image analysis of melanocytic skin tumours[END_REF][START_REF] Kollias | Optical non-invasive approaches to diagnosis of skin diseases[END_REF] .

AKNOWLEDGMENTS

This research was fully funded by Johnson & Johnson Santé Beauté France.

Figure 3 .

 3 Figure 3. Steps of the identification individual cells.A median filter and a local normalization were applied to the image with the ROI mask, followed by a Sato filter. Its output was filtered with a median filtered and locally normalized, and a Gabor filter was applied to it. A threshold was applied on the output after histogram equalization and small blobs were removed with a connected components analysis. The result was then skeletonized, and spurious branches were removed. ROI: region of interest.

Figure 4 .

 4 Figure 4. Post-processing steps. (A) The skeleton obtained after the previous step was cleaned, and contours were detected. Small contours were removed, as well as long contours close to the border with the background. The remaining contours were divided into two groups, small and big contours. Big contours were filtered again to

Table 1 .

 1 Table of the cell detection accuracy on stratum granulosum and stratum spinosum RCM images for two different experts

			𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷 =	𝑹𝑹𝑷𝑷𝑷𝑷𝑭𝑭𝑭𝑭𝑭𝑭 =
			𝑻𝑻𝑷𝑷𝑻𝑻𝑷𝑷 𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝒔𝒔𝑷𝑷𝒔𝒔𝑷𝑷𝑷𝑷	𝑻𝑻𝑷𝑷𝑻𝑻𝑷𝑷 𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝒔𝒔𝑷𝑷𝒔𝒔𝑷𝑷𝑷𝑷
			𝑻𝑻𝑷𝑷𝑻𝑻𝑷𝑷 𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝒔𝒔𝑷𝑷𝒔𝒔𝑷𝑷𝑷𝑷 + 𝑭𝑭𝑭𝑭𝑭𝑭𝑷𝑷𝑷𝑷 𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝒔𝒔𝑷𝑷𝒔𝒔𝑷𝑷𝑷𝑷	𝑻𝑻𝑷𝑷𝑻𝑻𝑷𝑷 𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝑷𝒔𝒔𝑷𝑷𝒔𝒔𝑷𝑷𝑷𝑷 + 𝑭𝑭𝑭𝑭𝑭𝑭𝑷𝑷𝑷𝑷 𝑵𝑵𝑷𝑷𝑵𝑵𝑭𝑭𝒔𝒔𝑷𝑷𝒔𝒔𝑷𝑷𝑷𝑷
	Stratum	Detection vs.	77%	85%
	granulosum	Expert 1
	RCM image	Detection vs.	64%	82%
		Expert 2
		Expert 1 vs.	52%	48%
		Expert 2
	Stratum	Detection vs.	85%	67%
	spinosum	Expert 1
	RCM image	Detection vs.	75%	30%
		Expert 2
		Expert 1 vs.	20%	28%
		Expert 2