Quentin Klopfenstein

Quentin Bertrand

Alexandre Gramfort

Joseph Salmon

Samuel Vaiter

Model identification and local linear convergence of coordinate descent

For composite nonsmooth optimization problems, Forward-Backward algorithm achieves model identification (e.g., support identification for the Lasso) after a finite number of iterations, provided the objective function is regular enough. Results concerning coordinate descent are scarcer and model identification has only been shown for specific estimators, the support-vector machine for instance. In this work, we show that cyclic coordinate descent achieves model identification in finite time for a wide class of functions. In addition, we prove explicit local linear convergence rates for coordinate descent. Extensive experiments on various estimators and on real datasets demonstrate that these rates match well empirical results.

Introduction 1.Coordinate descent

Over the last two decades, coordinate descent (CD) algorithms have become a powerful tool to solve large scale optimization problems [START_REF] Friedman | Pathwise coordinate optimization[END_REF][START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]. Many applications coming from machine learning or compressed sensing have lead to optimization problems that can be solved efficiently via CD algorithms: the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], the elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] or support-vector machine [START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF]. All the previously cited estimators are based on an optimization problem which can be 1 written:

x ∈ arg min

x∈R p {Φ(x) f (x) + p j=1 g j (x j ) g(x) } , (1)
with f a convex smooth (i.e., with a Lipschitz gradient) function and g j proper closed and convex functions. In the past twenty years, the popularity of CD algorithms has greatly increased due to the well suited structure of the new optimization problems mentioned above (i.e., separability of the nonsmooth term), as well as the possible parallelization of the algorithms [START_REF] Fercoq | Accelerated, parallel and proximal coordinate descent[END_REF].

The key idea behind CD (Algorithm 1) is to solve small and simple subproblems iteratively until convergence. More formally, for a function Φ : R p → R, the idea is to minimize successively one dimensional functions Φ |xj : R → R, updating only one coordinate at a time, while the others remain unchanged. There exists many variants of CD algorithms, the main branching being:

• The index selection. There are different ways to choose the index of the updated coordinate at each iteration. The main variants can be divided in three categories, cyclic CD [START_REF] Tseng | Block-coordinate gradient descent method for linearly constrained nonsmooth separable optimization[END_REF] when the indices are chosen in the set [p] {1, . . . , p} cyclically. Random CD [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF], where the indices are chosen following a given random distribution. Finally, greedy CD [START_REF] Nutini | Coordinate descent converges faster with the Gauss-Southwell rule than random selection[END_REF] picks an index, optimizing a given criterion: largest decrease of the objective function, or largest gradient norm (Gauss-Southwell rule), for instance.

• The update rule. There also exists several possible schemes for the coordinate update: exact minimization, coordinate gradient descent or prox-linear update (see Shi et al. 2016, Sec. 2.2 for details).

In this work, we will focus on cyclic CD with prox-linear update rule (Algorithm 1): a popular instance, e.g., the one coded in popular packages such as glmnet [START_REF] Friedman | Pathwise coordinate optimization[END_REF] or sklearn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. Among the methods of coordinate selection, random CD has been extensively studied, especially by [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF] for the minimization of a smooth function f . It was the first paper proving global non-asymptotic 1/k convergence rate in the case of a smooth and convex f . This work was later extended to composite optimization f + j g j for nonsmooth separable functions [START_REF] Richtárik | Iteration complexity of randomized block-coordinate descent methods for minimizing a composite function[END_REF][START_REF] Fercoq | Accelerated, parallel and proximal coordinate descent[END_REF]. Refined convergence rates were also shown by [START_REF] Shalev-Shwartz | Stochastic methods for l1-regularized loss minimization[END_REF][START_REF] Shalev-Shwartz | Stochastic dual coordinate ascent methods for regularized loss minimization[END_REF]. These convergence results have then been extended to coordinate descent with equality constraints [START_REF] Necoara | A random coordinate descent algorithm for optimization problems with composite objective function and linear coupled constraints[END_REF] that induce non-separability as found in the SVM dual problem in the presence of the bias term. Different distributions have been considered for the index selection such as uniform distribution [START_REF] Fercoq | Accelerated, parallel and proximal coordinate descent[END_REF][START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF][START_REF] Shalev-Shwartz | Stochastic methods for l1-regularized loss minimization[END_REF][START_REF] Shalev-Shwartz | Stochastic dual coordinate ascent methods for regularized loss minimization[END_REF], importance sampling (Leventhal and [START_REF] Lewis | Randomized methods for linear constraints: convergence rates and conditioning[END_REF][START_REF] Zhang | Solving large scale linear prediction problems using stochastic gradient descent algorithms[END_REF] and arbitrary sampling [START_REF] Necoara | A random coordinate descent algorithm for optimization problems with composite objective function and linear coupled constraints[END_REF]Qu and Richtárik, 2016a,b).

On the opposite, theory on cyclic coordinate descent is more fuzzy, the analysis in the cyclic case being more difficult. First, [START_REF] Luo | On the convergence of the coordinate descent method for convex differentiable minimization[END_REF]; [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]; [START_REF] Tseng | Block-coordinate gradient descent method for linearly constrained nonsmooth separable optimization[END_REF]; [START_REF] Razaviyayn | A unified convergence analysis of block successive minimization methods for nonsmooth optimization[END_REF] have shown convergence results for (block) CD algorithms for nonsmooth optimization problems (without rates1 ). Then, [START_REF] Beck | On the convergence of block coordinate type methods[END_REF] showed 1/k convergence rates for Lipschitz convex functions and linear convergence rates in the strongly convex case. [START_REF] Saha | On the nonasymptotic convergence of cyclic coordinate descent methods[END_REF] proved 1/k convergence rates for composite optimization f + . 1 under "isotonicity" condition. [START_REF] Sun | Improved iteration complexity bounds of cyclic block coordinate descent for convex problems[END_REF]; [START_REF] Hong | Iteration complexity analysis of block coordinate descent methods[END_REF] have extended the latter results and showed 1/k convergence rates with improved constants for composite optimization f + j g j . [START_REF] Li | On faster convergence of cyclic block coordinate descent-type methods for strongly convex minimization[END_REF] have extended the work of [START_REF] Beck | On the convergence of block coordinate type methods[END_REF] to the nonsmooth case and refined their convergence rates in the smooth case. Finally, as far as we know, the work by [START_REF] Xu | A globally convergent algorithm for nonconvex optimization based on block coordinate update[END_REF] is the first one tackling the problem of local linear convergence. They have proved local linear convergence under the very general Kurdyka-Lojasiewicz hypothesis, relaxing convexity assumptions. Following the line of work by [START_REF] Liang | Local linear convergence of forward-backward under partial smoothness[END_REF], we use a more restrictive framework (see Section 1.4) that allows to achieve finer results: model identification as well as improved local convergence results.

Model identification

Nonsmooth optimization problems coming from machine learning such as the Lasso or the support-vector machine (SVM) generally generate solutions lying onto a low-complexity model (see Definition 1 for details). For the Lasso, for example, a solution x has typically only a few non-zeros coefficients: it lies on the model set T x = {u ∈ R p : supp(u) ⊆ supp(x )}, where supp(x) is the support of x, i.e., the set of indices corresponding to the non-zero coefficients. A question of interest in the literature is: does the algorithm achieve model identification after a finite number of iterations? Formally, does it exist K > 0 such that for all k > K, x (k) ∈ T x ? For the Lasso the question boils down to "does it exist K > 0 such that for all k > K, supp(x (k) ) ⊆ supp(x )"? This finite time identification property is paramount for features selection [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF]), but also for potential acceleration methods [START_REF] Massias | Celer: a fast solver for the lasso with dual extrapolation[END_REF] of the CD algorithm, as well as model calibration [START_REF] Bertrand | Implicit differentiation of lasso-type models for hyperparameter optimization[END_REF].

Finite model identification was first proved in [START_REF] Bertsekas | On the Goldstein-Levitin-Polyak gradient projection method[END_REF] for the projected gradient method with non-negative constraints. In this case, after a finite number of steps the sparsity pattern of the iterates is the same as the sparsity pattern of the solution. It means that for k large enough, x (k) i = 0 for all i such that x i = 0. Then, many other results of finite model identification have been shown in different settings and for various algorithms. For the projected gradient descent algorithm, identification was proved for polyhedral constraints [START_REF] Burke | On the identification of active constraints[END_REF], for general convex constraints [START_REF] Wright | Identifiable surfaces in constrained optimization[END_REF], and even non-convex constraints [START_REF] Hare | Identifying active constraints via partial smoothness and prox-regularity[END_REF]. More recently, identification was proved for proximal gradient algorithm [START_REF] Mercier | Lectures on topics in finite element solution of elliptic problems[END_REF][START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], for the 1 regularized problem [START_REF] Hare | Identifying active manifolds in regularization problems[END_REF]. [START_REF] Liang | Local linear convergence of forward-backward under partial smoothness[END_REF][START_REF] Liang | Activity Identification and Local Linear Convergence of Forward-Backward-type Methods[END_REF]; [START_REF] Vaiter | Model consistency of partly smooth regularizers[END_REF] have shown model identification and local linear convergence for proximal gradient descent. These results have then been extended to other popular machine learning algorithms such as SAGA, SVRG [START_REF] Poon | Local convergence properties of SAGA/Prox-SVRG and acceleration[END_REF] and ADMM [START_REF] Poon | Trajectory of alternating direction method of multipliers and adaptive acceleration[END_REF]. To our knowledge, CD has not been extensively studied with a similar generality. Some identification results have been shown for CD, but only on specific models [START_REF] She | Linear convergence and support vector identification of sequential minimal optimization[END_REF][START_REF] Massias | Dual extrapolation for sparse generalized linear models[END_REF] or variants of CD [START_REF] Wright | Accelerated block-coordinate relaxation for regularized optimization[END_REF], in general, under restrictive hypothesis. The authors are not aware of generic model identification results for CD Algorithm 1.

Notation

General notation. We write • the Euclidean norm on vectors. For, x, γ ∈ R p , the weighted norm is denoted x γ p j=1 γ j x 2 j . For a differentiable function ψ : R p → R p , at x ∈ R p , we write J ψ(x) ∈ R p×p the Jacobian of ψ at x. For a set S, we denote by S c its complement. We denote [p] = {1, . . . , p}. Let (e j ) p j=1 be the vectors of the canonical base of R p . We denote the coordinatewise multiplication of two vectors u and v by u v and by u M the row wise multiplication between a vector and a matrix. We denote by B(x, ) the ball of center x and radius . The spectral radius of a matrix M is denoted ρ(M ).

Convex analysis. We recall the definition of the proximity operator of a convex function g, for any γ > 0:

prox γg (x) = arg min y∈R p 1 2γ ||x -y|| 2 + g(y) .
Let C ⊂ R p be a convex set, aff(C) denotes its affine hull, the smallest affine set containing C, and ri(C) denotes its relative interior (the interior of its affine hull). The indicator function of C is the function defined for any x ∈ R p by

δ C (x) = 0 if x ∈ C +∞ otherwise . ( 2
)
The domain of a function f is defined as dom

(f ) = {x ∈ R p : f (x) < +∞}.
For a convex function f , ∂f (x) denotes its subdifferential at x and is given by ∂f

(x) = {s ∈ R p : f (y) ≥ f (x) + s, y -x , ∀y ∈ dom(f )}.
We denote by L j the coordinatewise Lipschitz constants of ∇ j f , i.e., , for all x ∈ R p , h j ∈ R:

||∇ j f (x + e j h j ) -∇ j f (x)|| ≤ L j |h j | . (3)
Coordinate descent. We denote 0 < γ j ≤ 1/L j the local step size and γ = (γ 1 , . . . , γ p ) . To prove model identification we need to "keep track" of the iterates: following the notation from [START_REF] Beck | On the convergence of block coordinate type methods[END_REF] coordinate descent can be written:

Algorithm 1 Proximal coordinate descent input : γ 1 , . . . , γ p ∈ R + , n iter ∈ N, x (0) ∈ R p for k = 0, . . . , n iter do // index selection x (0,k) ← x (k) for j = 1, . . . , p do x (j,k) ← x (j-1,k) x (j,k) j ← prox γj gj x (j-1,k) j -γ j ∇ j f (x (j-1,k) ) x (k+1) ← x (p,k) return x niter+1

Assumptions on composite problem

We consider the optimization problem defined in Equation ( 1) with the following assumptions:

Assumption 1 (Smoothness). f is a convex and differentiable function, with a Lipschitz gradient.

Assumption 2 (Proper, closed, convex). For any j ∈ [p], g j is proper, closed and convex.

Assumption 3 (Existence). The problem admits at least one solution: arg min

x∈R p Φ(x) = ∅ . ( 4
)
Assumption 4 (Non degeneracy). The problem is non-degenerate: for any x ∈ arg min x∈R p Φ(x)

-∇f (x ) ∈ ri (∂g(x )) . (5)
Assumption 4 can be seen as a generalization of qualification constraints (Hare and Lewis, 2007, Sec. 1).

Contributions

With mild assumptions on the g j functions, for the cyclic proximal coordinate descent algorithm:

• We prove finite time model identification (Th. 1).

• We provide local linear convergence rates (Th. 2).

• We illustrate our results on multiple real datasets and estimators (Section 4) showing that our theoritical rates match the empirical ones.

Model identification for CD

As stated before, the solutions of the Lasso are structured. Using an iterative algorithm like Algorithm 1 to find an approximate solution (since we stop after a finite number of iterations) brings the question of structure recovery. For the Lasso, the underlying structure, also called model [START_REF] Candès | Simple bounds for recovering low-complexity models[END_REF], is identified by the Forward-Backward algorithm. It means that after a finite number of iterations, the iterative algorithm leads to an approximate solution that shares a similar structure than the true solution of the optimization problem [START_REF] Liang | Local linear convergence of forward-backward under partial smoothness[END_REF][START_REF] Vaiter | Model consistency of partly smooth regularizers[END_REF][START_REF] Fadili | Sensitivity analysis for mirror-stratifiable convex functions[END_REF]. For the Lasso, the underlying model is related to the notion of support: i.e., the non-zero coefficients for the Lasso, and it can be generalized for the case of completely separable functions as follows:

Definition 1 (Generalized support, [START_REF] Sun | Are we there yet? manifold identification of gradient-related proximal methods[END_REF]. We call generalized support S x ⊆ [p] the set of indices j ∈ [p] where g j is differentiable at x j :

S x {j ∈ [p] : ∂g j (x j ) is a singleton} . ( 6
)
This notion can be unified with the definition of model subspace from Vaiter et al. (2015, Sec. 3.1):

Definition 2 (Model subspace, [START_REF] Vaiter | Model selection with low complexity priors[END_REF]. We denote the model subspace at x:

T x = {u ∈ R p : ∀j ∈ S c
x , u j = 0} .

See Lemma 1 in Appendix A for details.

Examples in machine learning. The 1 norm. The function g(x) = p i=1 |x i | is certainly the most popular nonsmooth convex regularizer promoting sparsity. Indeed, the 1 norm generates structured solution with model subspace [START_REF] Vaiter | Model consistency of partly smooth regularizers[END_REF]. We have that S x = {j ∈ [p] : x j = 0} since | • | is differentiable everywhere but not at 0, and the model subspace reads:

T x = {u ∈ R p : supp(u) ⊆ supp(x)} . ( 8
)
The box constraints indicator function δ [0,C] . This indicator function appears for instance in box constrained optimization problems such as the dual problem of the SVM. Let

I 0 x = {j ∈ [p] : x j = 0} and I C x = {j ∈ [p] : x j = C}, then T x = {u ∈ R p : I 0 x ⊆ I 0 u and I C x ⊆ I 0 u }.
For the SVM, model identification boils down to finding the active set of the box constrained quadratic optimization problem after a finite number of iterations. We now turn to our identification result. To ensure model identification, we need the following (mild) assumption:

Assumption 5 (Locally C 2 ). For all j ∈ S x , g j is locally C 2 around x j , and f is locally C 2 around x .

It is satisfied for the Lasso and the dual SVM problem mentioned above, but also for sparse logistic regression or elastic net. The following theorem shows that the CD (Algorithm 1) has the model identification property with local constant step size 0 < γ j ≤ 1/L j :

Theorem 1 (Model identification of CD). Consider a solution x ∈ arg min x∈R p Φ(x) and S = S x . Suppose 1. Assumptions 1 to 5 hold.

2. The sequence (x (k) ) k≥0 generated by Algorithm 1 converges to x .

Then, Algorithm 1 identifies the model after a finite number of iterations, which means that there exists K > 0 such that for all k ≥ K, x (k)

S c = x S c .
This result implies that for k large enough, x (k) shares the support of x (potentially smaller).

Sketch of proof [Theorem 1]

• First we show that Assumptions 1 to 5 implies that g is partly smooth [START_REF] Lewis | Active sets, nonsmoothness, and sensitivity[END_REF] at x relative to the affine space x + T x .

• Then we show that for the CD Algorithm 1: dist ∂Φ(x (k) ), 0 → 0, when k → ∞, enabling us to apply Hare and Lewis (2004)[Thm. 5.3].

A full proof of Theorem 1 can be found in Appendix A. The first point is shown in appendix Lemma 1. We show the second point below: Proof As written in Algorithm 1, one update of coordinate descent reads:

1 γ j x (j-1,k) j -∇ j f x (j-1,k) - 1 γ j x (j,k) j ∈ ∂g j x (j,k) j 1 γ j x (k) j -∇ j f x (j-1,k) - 1 γ j x (k+1) j ∈ ∂g j x (k+1) j .
Since g is separable with non-empty subdifferential, the coordinate wise subdifferential of g is equal to the subdifferential of g, we then have

1 γ x (k) -∇ j f x (j-1,k) j∈[p] - 1 γ x (k+1) ∈ ∂g(x (k+1) ) , (9)
which leads to

1 γ x (k) -∇ j f x (j-1,k) j∈[p] - 1 γ x (k+1) + ∇f (x (k+1) ) ∈ ∂Φ(x (k+1) ) . ( 10
)
To prove support identification using Hare and Lewis (2004, Theorem 5.3), we need to bound the distance between ∂Φ(x (k+1) ) and 0, using Equation (10):

dist ∂Φ(x (k+1) ), 0

2 ≤ p j=1 x (k) j γ j -∇ j f (x (j-1,k) ) - x (k+1) j γ j + ∇ j f (x (k+1) ) 2 ≤ ||x (k) -x (k+1) || 2 γ -1 + p j=1 ∇ j f x (j-1,k) -∇ j f x (k+1) 2 ≤ ||x (k) -x (k+1) || 2 γ -1 + L 2 p j=1 x (j-1,k) -x (k+1) 2 ≤ ||x (k) -x (k+1) || 2 γ -1 + L 2 p j=1 p j ≥j x (k) j -x (k+1) j 2 →0 when k→∞ .
We thus have:

• dist ∂Φ(x (k+1) ), 0 → 0 • Φ(x (k) ) → Φ(x )
because Φ is prox-regular (since it is convex, see [START_REF] Poliquin | Prox-regular functions in variational analysis[END_REF] and subdifferentially continuous.

Then the conditions to apply Hare and Lewis (2004, Th. 5.3) are met and hence we have model identification after a finite number of iterations.

Comments on Theorem 1. It unifies several results found in the literature: [START_REF] Massias | Dual extrapolation for sparse generalized linear models[END_REF] showed model identification for the Lasso, solved with coordinate descent, but requiring uniqueness assumption. [START_REF] Nutini | Let's make block coordinate descent go fast: Faster greedy rules, message-passing, active-set complexity, and superlinear convergence[END_REF] showed some identification results under strong convexity assumption on f . Theorem 1 do not rely on any uniqueness, strong convexity, or local strong convexity hypothesis. Even if the solution of the optimization problem defined in Equation ( 1) is not unique, CD achieves model identification.

Local convergence rates

In this section, we prove local linear convergence of the CD Algorithm 1. After model identification, there exists a regime where the convergence towards a solution of Equation ( 1) is linear. Local linear convergence was already proved in various settings such as for ISTA and FISTA algorithms (i.e., with an 1 penalty, [START_REF] Tao | Local linear convergence of ISTA and FISTA on the LASSO problem[END_REF]) and then for the general Forward-Backward algorithm [START_REF] Liang | Local linear convergence of forward-backward under partial smoothness[END_REF].

Local linear convergence requires an additional assumption: restricted injectivity. It is classical for this type of analysis as it can be found in [START_REF] Liang | Activity Identification and Local Linear Convergence of Forward-Backward-type Methods[END_REF] and [START_REF] Poon | Trajectory of alternating direction method of multipliers and adaptive acceleration[END_REF]. Assumption 6. (Restricted injectivity) For a solution x ∈ arg min x∈R p Φ(x), the restricted Hessian to its generalized support S = S x is definite positive, i.e.,

∇ 2 S,S f (x ) 0 . ( 11
)
For the Lasso, Assumption 6 is a classical necessary condition to ensure uniqueness of the minimizer [START_REF] Fuchs | On sparse representations in arbitrary redundant bases[END_REF].

In order to study local linear convergence, we consider the fixed point iteration of a complete epoch (an epoch is a complete pass over all the coordinates). A full epoch of CD can be written:

x (k+1) = ψ(x (k) ) P p • . . . • P 1 (x (k) ) , (12)
where P j are coordinatewise sequential applications of the proximity operator P j : R p → R p :

x →             x 1 . . . x j-1 prox γj gj x j -γ j ∇ j f (x)
x j+1 . . .

x p            
.

Thanks to model identification (Theorem 1), we are able to prove that once the generalized support is correctly identified, there exists a regime where CD algorithm converges linearly towards the solution:

Theorem 2 (Local linear convergence). Consider a solution x ∈ arg min x∈R p Φ(x) and S = S x . Suppose 1. Assumptions 1 to 6 hold.

2. The sequence (x (k) ) k≥0 generated by Algorithm 1 converges to x .

3. The model has been identified i.e., there exists K ≥ 0 such as for all k ≥ K x (k)

S c = x S c .
Then (x (k) ) k≥K converges linearly towards x . More precisely, for any ν ∈ [ρ(J ψ S,S (x )), 1[, there exists K > 0 and a constant C such that for all k ≥ K, The complete proof of Theorem 2 can be found in Appendix B.

x (k) S -x S ≤ Cν (k-K) x (K) S -x S .

Sketch of proof [Theorem 2]

• A key element of the proof is to consider a full epoch of CD: it can be written as a fixed point iteration: x (k+1) = ψ(x (k) ) (see Equation ( 12)).

• We then show that the proximal operators, prox γj gj , evaluated at x jγ j ∇ j f (x ) are differentiable (for j ∈ [p]). Once stated, the differentiability of the proximal operator allows us to write the Taylor expansion of ψ:

x (k+1) -x = ψ(x k ) -ψ(x ) = J ψ(x ), x (k) -x + o( x (k) -x ).
• Capitalizing on model identification (Theorem 1) we start from x (k) S c = x S c and show the bound ρ(J ψ S,S (x )) < 1 on the spectral radius of the restricted Jacobian of ψ at x : J ψ S,S (x ).

• Finally, all the conditions are met to apply Polyak (1987, Th. 1, Sec. 2.1.2). The latter reference provides sufficient conditions for local linear convergence of sequences based non linear fixed point iterations.

Experiments

We now illustrate Theorems 1 and 2 on multiple datasets and estimators: the Lasso, the logistic regression and the SVM. In this section, we consider a design matrix A ∈ R n×p and a target y ∈ R n for regression (Lasso) and y ∈ {-1, 1} n for classification (logistic regression and support-vector machine). We used classical datasets from libsvm [START_REF] Chang | Libsvm: a library for support vector machines[END_REF] summarized in Table 1.

In Figures 1 to 3 the distance of the iterates to the optimum, x (k)x as a function of the number of iterations k is plotted as a solid blue line. The vertical red dashed line represents the iteration k where the model has been identified by CD (Algorithm 1) illustrating Theorem 1. The yellow dashed line represents the theoretical linear rate from Theorem 2. Theorem 2 gives the slope of the dashed yellow line, the (arbitrary) origin point of the theoretical rate line is chosen such that blue and yellow lines coincide at identification time, i.e., all lines intersect at this point. More precisely, if k denotes the iteration where model identification happens, the equation of the dashed yellow line is:

h(k) = x (k ) -x × ρ(J ψ S,S (x )) (k-k ) . ( 13
)
Once a solution x has been computed, one can calculate J ψ S,S (x ) and its spectral radius for each estimator.

For the experiments we used three different estimators that we detail here. Lasso. [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] The most famous estimator based on a nonsmooth optimization problem may be the Lasso. For a design matrix A ∈ R n×p and a target y ∈ R n it writes: arg min

x∈R p 1 2n ||Ax -y|| 2 + λ||x|| 1 . ( 14
)
The CD update for the Lasso is given by

x j ← ST γj λ x j -γ j A :,j (y -Ax) , (15)
where ST λ (x) = sign(x) • max(|x|λ, 0). The solution of Equation ( 14) is obtained using Algorithm 1 with constant stepsizes 1/γ j = ||A:,j || 2 n . Sparse logistic regression. The sparse logistic regression is an estimator for classification tasks. It is the solution of the following optimization problem, for a design matrix A ∈ R n×p and a target variable y ∈ {-1, 1} n , with σ(z)

1 1+e -z : arg min x∈R p - 1 n n i=1 log σ(y i x A i,: ) + λ x 1 . (16)
The CD update for the sparse logistic regression is

x j ← ST γj λ x j -γ j A :,j (y (σ(y Ax) -1)) . (17)
The constant stepsizes for the CD algorithm to solve Equation ( 16) are given by 1/γ j = ||A:,j || 2 4n . Support-vector machine. [START_REF] Boser | A training algorithm for optimal margin classifiers[END_REF] The support-vector machine (SVM) primal optimization problem is, for a design matrix A ∈ R n×p and Table 2: C values for SVM. dataset leukemia gisette rcv1 20news C value 10 1.5 10 -2 1.5 10 -2 5 10 -1 a target variable y ∈ {-1, 1} n : arg min

x∈R p 1 2 x 2 + C n i=1 max 1 -y i x A i,: , 0 . (18)
The SVM can be solved using the following dual optimization problem: arg min

w∈R n 1 2 w (y A)(y A) w - n i=1 w i subject to 0 ≤ w i ≤ C . (19)
The CD update for the SVM reads:

w i ← P [0,C] w i -γ i ((y A) i,: (y Aw) -1 ) , (20)
where P [0,C] (x) = min(max(0, x), C). The stepsizes of the CD algorithm to solve Equation ( 19) are given by 1/γ i = ||(y A) i,: || 2 . The values of the regularization parameter C for each dataset from Figure 3 are given in Table 2.

Comments on Figures 1 to 3. Finite time model identification and local linear convergence are illustrated on the Lasso, the sparse logistic regression and the SVM in Figures 1 to 3. As predicted by Theorem 1, the relative model is identified after a finite number of iterations. For the Lasso (Figure 1) and the sparse logistic regression (Figure 2), we observe that as the regularization parameter gets smaller, the number of iterations needed by the CD algorithm to identify the model increases. To our knowledge, this is a classical empirical observation, that is not backed up by theoretical results. After identification, the convergence towards a solution is linear as predicted by Theorem 2. The theoretical local speed of convergence provided by Theorem 2 seems like a sharp estimation of the true speed of convergence as illustrated by the three figures.

Note that on Figures 1 to 3 high values of λ (or small values of C) were required for the restricted injectivity Assumption 6 to hold. Indeed, despite its lack of theoretical foundation, it is empirically observed that, in general, the larger the value of λ, the smaller the cardinal of the generalized support: |S|. It makes the restricted injectivity Assumption 6: ∇ 2 S,S f (x ) 0 easier to be satisfied. For instance, for λ = λ max /20, the restricted injectivity Assumption 6 was not verified for a lot of datasets for the Lasso and the sparse logistic regression (Figures 1 and2). In the same vein, values of C for the SVM had to be chosen small enough, in order to make |S| not too large (Figure 3).

Note that finite time model identification is crucial to ensure local linear convergence, see for instance 20news dataset on Figure 3. However there exists very few quantitative theoretical results for the convergence speed of the model identification. [START_REF] Nutini | active-set complexity" of proximal gradient: How long does it take to find the sparsity pattern?[END_REF]; [START_REF] Sun | Are we there yet? manifold identification of gradient-related proximal methods[END_REF] tried to obtain some rates on the identification, quantifying "how much the problem is qualified", i.e., how much Assumption 4 is satisfied. But these theoretical results do not seem to explain fully the experimental results of the CD: in particular the identification speed of the model compared to other algorithms.

Limits. We would like to point out the limit of our analysis illustrated for the case of λ = λ max /15 for the sparse logistic regression and the rcv1 dataset in Figure 2. In this case, the solution may no longer be unique. The support gets larger and Assumption 6 is no longer met. In this case, the largest eigenvalue of J ψ S,S (x ) is exactly one, which leads to the constant rate observed in Figure 2. Despite the largest eigenvalue being exactly 1, a regime of locally linear convergence toward a (potentially non unique) minimizer is still observed. Linear convergence of non-strongly convex functions starts to be more and more understood [START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF]. Figure 2 with λ = λ max /15 for rcv1 suggests extensions of [START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF] could be possible in the nonsmooth case.

Conclusion and future work. In conclusion, we show finite time model identification for coordinate descent Algorithm 1 (Theorem 1). Thanks to this identification property we were able to show local linear rates of convergence (Theorem 2). These two theoretical results were illustrated on popular estimators (Lasso, sparse logistic regression and SVM dual) and popular machine learning datasets (Section 4).

A first natural extension of this paper would be to investigate block coordinate minimization: Theorem 1 could be extended for blocks under general partial smoothness assumption [START_REF] Hare | Identifying active constraints via partial smoothness and prox-regularity[END_REF]). However, it seems that Theorem 2 would require a more careful analysis. A second extension could be to show linear convergence without the restricted injectivity (Assumption 6), paving the way for a generalization of [START_REF] Necoara | Linear convergence of first order methods for non-strongly convex optimization[END_REF] as suggested by Figure 2.

A Proofs of model identification (Theorem 1)

Model identification often relies on the assumption that the nonsmooth function g is regular enough, or more precisely partly smooth. Loosely speaking, a partial smooth function behaves smoothly as it lies on the related model and sharply if we move normal to that model. Formally, we recall the definition of partly smooth functions restricted to the case of proper, lower semicontinuous and convex functions.

Definition 3 (Partial smoothness). Let g : R p → R be a proper closed convex function. g is said to be partly smooth at x relative to a set M ⊆ R n if there exists a neighbourhood U of x such that

• (Smoothness) M ∩ U is a C 2 -manifold and g restricted to M ∩ U is C 2 ,
• (Sharpness) The tangent space of M at x is the model tangent space T x where T x = Lin(∂g(x)) ⊥ ,

• (Continuity) The set valued mapping ∂g is continuous at x relative to M.

The class of partly smooth functions was first defined in [START_REF] Lewis | Active sets, nonsmoothness, and sensitivity[END_REF]. It encompasses a large number of known nonsmooth machine learning optimization penalties, such as the 1 -norm or box constraints to only name a few, see Vaiter et al. (2018, Section 2.1) for details. Interestingly, this framework enables powerful theoretical tools on model identification such as Hare and Lewis (2004)[Thm. 5.3]. For separable functions, next lemma gives an explicit link between the generalized support Definition 1 [START_REF] Sun | Are we there yet? manifold identification of gradient-related proximal methods[END_REF] and the framework of partial smooth functions [START_REF] Hare | Identifying active constraints via partial smoothness and prox-regularity[END_REF].

Lemma 1. Let x ∈ dom (g). If for every j ∈ S x , g j is locally C 2 around x j (Assumption 5), then g is partly smooth at x relative to x + T x .

Proof. We need to prove the three properties of the partial smoothness (Definition 3).

Smoothness. Let us write M x = x + T x the affine space directed by the model subspace and pointed by x . In particular, it is a C 2 -manifold.

For every j ∈ S x , g j is locally C 2 around x j , hence there exists a neighborhood U j of x j such that the restriction of f to U is twice continuously differentiable. For j ∈ S c

x , let's write U j = R. Take U = j∈[p] U j . This a neighborhood of x (it is open, and contains x ). Consider the restriction g |M x of g to M x . It is C 2 at each point of U since each coordinates (for j ∈ S x ) are C 2 around U j .

Sharpness. Since g is completly separable, we have that ∂g(x ) = ∂g 1 (x 1 )× . . . × ∂g p (x p ). Note that ∂g j (x j ) is a set valued mapping which is equal to the singleton {∇ j g(x j )} if g j is differentiable at x j or it is equal to an interval. The model tangent space T x of g at x is given by

T x = span(∂g(x )) ⊥ where span(∂g(x )) = aff(∂g(x )) -e x , (21)
with

e x = arg min e∈aff(∂g(x )) ||e|| , (22)
called the model vector.

In the particular case of separable functions, we have that aff (∂g(x )) = aff ∂g 1 (x 1 ) × . . . × ∂g p (x p ) = aff (∂g 1 (x 1 )) × . . . × aff ∂g p (x p ) .

In this case,

aff ∂g j (x j ) = {∇ j g(x j )} if j ∈ S x R otherwise and e x j = ∇ j g(x j ) if j ∈ S x 0 otherwise . (23)
Thus we have that span (∂g(x )) = aff (∂g(x ))e x = {x ∈ R p : ∀j ∈ S x , x j = 0} .

Then

T x = span (∂g(x )) ⊥ = {x ∈ R p : ∀j ∈ S c x , x j = 0} . (24)
Continuity. We are going to prove that ∂g is inner semicontinuous at x relative to M x , i.e., that for any sequence (x (k) ) of elements of M x converging to x and any η ∈ ∂g(x ), there exists a sequence of subgradients η (k) ∈ ∂g(x (k) ) converging to η.

Let x (k) be a sequence of elements of M x converging to x , or equivalently, let t (k) be a sequence of elements of T x converging to 0, and let η ∈ ∂g(x ).

For j ∈ S x , we choose η (k) j g j (x j +t (k) j ), using the smoothness property we have η ∈ ∂g(x (k) ). We have that η (k) ∈ ∂g(x k ) and η (k) converges towards η since g j is C 1 around x j for j ∈ S x , hence, g j (x j + t (k) j ) converges to g j (x j ) = ηj . Thus, it proves that g is partly smooth at x relative to x + T x .

The end of the proof of Theorem 1 is contained in Section 2.

B Proofs of local linear convergence (Theorem 2)

To simplify the notations in this section, S S x . Let us also write the element of S as follows: S = {j 1 , . . . , j |S| }. The first point of this proof is to write the CD algorithm as a fixed point iteration. A full epoch of CD can be written as x (k+1) = ψ(x (k) ) = P p • . . . • P 1 (x (k) ) .

(25) and

B j = M 1/2
:,j γ j ∂ x prox γj gj z j M 1/2 :,j .

Since only one coordinate change at each step, the chain rule leads to J ψ(x S ) = J P js (x S )J P js-1 (x S ) . . . J P j1 (x S ) = M -1/2 (Id -B js ) . . . (Id -B j1 )

A M 1/2
The next series of lemma will be usefull to prove that the spectral radius ρ J ψ(x S ) < 1.

Lemma 4. The matrix M defined in (38) is symmetric definite positive.

Proof. Using the non-expansivity of the prox, and the property ∂ x prox γj gj (z j ) > 0 for j ∈ S, diag(u) is a symmetric semidefinite matrix, so M is a sum of a symmetric definite positive matrix and a symmetric semidefinite matrix, hence M is symmetric definite positive.

Lemma 5. For all j ∈ S, the matrix B j defined in (40) has spectral norm bounded by 1, i.e., B j 2 ≤ 1.

Proof. B j is a rank one matrix which is the product of γ j ∂ x prox γj gj (z j )M 1/2 :,j and M 1/2 :,j

, its non-zeros eigenvalue is thus given by B j 2 = M 1/2 :,j γ j ∂ x prox γj gj (z j )M 1/2 :,j = γ j ∂ x prox γj gj (z j )M j,j = γ j ∂ x prox γj gj (z j )

      ∇ 2 j,j f (x ) 0≤ + 1 γ j ∂ x prox γj gj (z j ) - 1 γ j 0≤       . ( 41
)
By positivity of the two terms,

B j 2 = γ j ∂ x prox γj gj (z j ) ∇ 2 j,j f (x ) ≤Lj ≤ 1 γ j + 1 -∂ x prox γj gj (z j ) ≤ ∂ x prox γj gj (z j ) + 1 -∂ x prox γj gj (z j ) ≤ 1 . ( 42
)
Lemma 6. For all j ∈ S, B j / B j is an orthogonal projector onto Span(M 1/2 :,j ).

Figure 1 :

 1 Figure 1: Lasso, linear convergence. Distance to optimum, x (k)x , as a function of the number of iterations k, on 4 different datasets: leukemia, gisette, rcv1, and real-sim.

Figure 2 :

 2 Figure 2: Sparse logistic regression, linear convergence. Distance to optimum, x (k)x , as a function of the number of iterations k, on 4 different datasets: leukemia, gisette, rcv1, and real-sim.

Figure 3 :

 3 Figure 3: Support vector machine, linear convergence. Distance to optimum, x (k)x , as a function of the number of iterations k, on 4 different datasets: leukemia, gisette, rcv1 and 20news.

Table 1 :

 1 Characteristics of the datasets.

		Datasets #samples n #features p	density
		leukemia	38	7129	1
		gisette	6000	4955	1
		rcv1 real-sim 20news	20,242 72,309 5184	19,959 20,958 155,148	3.6 × 10 -3 2.4 × 10 -3 1.9 × 10 -3
		practical rate		theoretical rate	model identification
	0	50

Note that some local rates are shown in[START_REF] Tseng | Block-coordinate gradient descent method for linearly constrained nonsmooth separable optimization[END_REF] but under some strong hypothesis.

Acknowledgements

This work was partly supported by ANR GraVa ANR-18-CE40-0005 and ERC Starting Grant SLAB ERC-StG-676943.

We also define π x S c : R |S| → R p for all x S ∈ R |S| and all j ∈ S by

and for all j s ∈ S, Px S c js : R |S| → R |S| is the function defined for all x S ∈ R |S| and all j ∈ S by Px S c js (x S ) j =

x j if j = j s prox γj gj x jsγ js ∇ js f (π x S c (x S )) if j = j s .

(

Once the model is identified (Theorem 1), we have that there exists K ≥ 0 such that for all k ≥ K, we have that

When no confusion is possible, we denote by Pj the function Px S c js , hence still dependant on x S c . The following lemma shows that Pj is differentiable at the optimum.

Lemma 2. For all j ∈ S, Pj is differentiable at x S .

Proof. From Assumption 5, we know there exists a neighboorhood of x j denoted U such that, for j ∈ S, the restriction of g j to U is C 2 on U. In particular, it means that x j is a differentiable point of g j and given a pair (u, v)

we have 1 γj (vu) ∈ ∂g j (u) becomes

Let H(u) = (Id + g j )(u), since g j is twice differentiable at u, we have that

Thus, H : U → R is continuous and then H :

Using the implicit functions theorem, we have that there exists an open interval

Using (29) we thus have with the choice u = x j , v = x j -γ j ∇ j f (x ) that the map h coincides with prox γj gj on V and is differentiable at

For the sake of completness, we show that in fact prox γj gj is also differentiable on the complement of the generalized support at x j -∇ j f (x ).

Lemma 3. For all j ∈ S c , prox γj gj is constant around x j -∇ j f (x ). Moreover, the map x → prox γj gj (x j -∇ j f (x)) is differentiable at x with gradient 0.

Proof. Let ∂g j (x j ) = [a; b] and let z j = x j -∇ j f (x ), then combining the fixed point equation and Assumption 4 leads to:

Thus,

For all v ∈]γ j a+x j ; γ j b+x j [, we have 1 γj (v-x j ) ∈]a; b[= ri ∂g j (x j ) , i.e., prox γj gj (v) = x j . As f is C 2 in x , we have that x → prox γj gj (x j -∇ j f (x)) is differentiable at x with gradient being 0.

From Lemma 2, we have that Pj is differentiable at x S for all j ∈ S. Since x is an optimal point, the following fixed points equation holds:

The map ψ is then differentiable at x S since it is obtained as the composition of differentiable functions and that each function Pj is evaluated at a differentiable point (only one coordinate change at each step).

To compute, the Jacobian of Pj at x S , let us first notice that J P j (x S ) = e 1 . . . e j-1 v j e j+1 . . . e s ,

where v j = ∂ x prox γj gj z j e jγ j ∇ 2 j,: f (x ) and z j = x jγ j ∇ j f (x ). This matrix can be rewritten as J Pj (x S ) = Id |S| -e j e j + ∂ x prox γj gj z j e j e jγ j e j e j ∇ 2 f (x )

where

and u ∈ R |S| is defined for all j ∈ S by

Proof. It is clear that B j /||B j || is symmetric. We now prove that it is idempotent, i.e., (

Hence, B j /||B j || is an orthogonal projector.

Lemma 7. For all j ∈ S and for all

Proof. (44)

Combining Equations ( 43) and ( 44) leads to:

Lemma 8. The spectral norm of A is bounded by 1, i.e., (Id -B js ) . . .

we thus have for all j ∈ S, (Id -B j )x = x . One can thus successively apply Lemma 7 which leads to:

Moreover M 1/2 has full rank (see Lemma 4), thus x = 0 and (Id -B js ) . . . (Id -B j1 ) 2 < 1 .

From Lemma 8, ||A|| 2 < 1. Moreover A and J ψ(x S ) are similar matrices, then ρ(J ψ(x S )) = ρ(A) ≤ ||A|| 2 < 1.

To summarize, x S is the solution of a fixed point equation ψ(x S , x S c ) = x S . From Lemma 2, ψ(., x S c ) is differentiable at x S and the Jacobian at x S satifies the condition ρ(J ψ(x S )) < 1. Then all conditions are met to apply [START_REF] Polyak | Introduction to optimization. optimization software[END_REF][Theorem 1, Section 2.1.2] which proves local linear convergence.