Roland Martin

Jérémie Giraud

Vitaliy Ogarko

Sébastien Chevrot

Stephen Beller

Pascal Gegout

Mark Jessell

Pascal Gégout

Three-dimensional gravity anomaly inversion in the Pyrenees using compressional seismic velocity model as structural similarity constraints

Keywords:

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

INTRODUCTION

During the last decade, seismic full waveform inversion has become a powerful and an increasingly popular technique for Earth imaging (Tarantola & Valette 1982a,b;[START_REF] Virieux | An overview of full-waveform inversion in exploration geophysics[END_REF][START_REF] Monteiller | A hybrid method to compute short-period synthetic seismograms of teleseismic body waves in a 3-D regional model[END_REF][START_REF] Monteiller | Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM-DSM hybrid method[END_REF]. However, with this method, density is retrieved with a lower resolution compared to seismic velocities or Lamé parameters, and is often not sufficiently reliable for interpretation. Recovery of density from seismic data alone is difficult and may require the utilisation of prior information [START_REF] Alemie | High-resolution three-term avo inversion by means of a trivariate cauchy probability distribution[END_REF][START_REF] Lines | Density and AVO[END_REF]. Gravity data inversion is the common way of retrieving densities but has a poor depth resolution and thus needs additional constraints. In this context, complementarity with seismic data has been recognized [START_REF] Lelievre | Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration[END_REF] as seismic tomography provides seismic velocity models that can be exploited jointly with gravity data for improved subsurface density imaging. Another motivation to perform gravity inversion constrained by seismic modeling is that densities inferred from compressional seismic velocity models through Nafe-Drake or linear scaling relationships may be affected by significant uncertainties that can exceed 10% [START_REF] Herceg | Sensitivity analysis of crustal correction for calculation of lithospheric mantle density from gravity data[END_REF].

One of our main motivations is thus to reduce significantly these uncertainties by recovering densities through gravity field data inversion constrained by seismic velocities, thereby preserving physical and structural similarities between seismic velocities and densities.

Another long-term goal is to refine existing density models at depth to get insights into the composition and thermal state of the lithosphere and to confront inverted models with existing geological models.

In this perspective, we use teleseismic FWI (TFWI) tomography as a priori information to constrain gravity inversion and to refine existing density models. Scaling laws provide useful information on density from crustal compressional velocity (V p ) models [START_REF] Brocher | Empirical relations between elastic wavespeeds and density in the Earth's crust[END_REF]) that can be used as prior information to constrain solutions of gravity inversion at lithospheric scale. In this context, we exploit the results of TFWI on data collected in the Western Pyrenees during the 2009-2014 PYROPE and IB-ERARRAY seismic experiments [START_REF] Chevrot | High-resolution imaging of the Pyrenees and Massif Central from the data of the PYROPE and IBERARRAY portable array deployments[END_REF][START_REF] Chevrot | The Pyrenean architecture as revealed by teleseismic P-to-S converted waves recorded along two dense transects[END_REF][START_REF] Chevrot | The non-cylindrical crustal architecture of the Pyrenees[END_REF] to constrain the inversion of gravity data provided by the BGI/Bureau Gravimétrique International (International Gravimetry Bureau). A V p seismic model has been obtained by full waveform of teleseismic P-waves recorded in the westernmost PYROPE transect [START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P-waves[END_REF]) in a 3D slice crossing the gravity anomaly of Mauleon's basin. This 3D slice was involving the seismic profile surrounded by temporary seismic stations to constrain laterally the TFWI. In a previous study in this region, where the strongest gravity anomalies in France are present [START_REF] Martin | A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data -Application to the western pyrenees[END_REF], we reconstructed density distributions with an unconstrained gravity inversion procedure but densities still remained poorly constrained at depth. A conclusion of that study was that further inversions were needed with additional constraints and prior information to provide robust images of deeper structures : without sufficient prior information, gravity inversion is heavily affected by non-uniqueness and the inverse problem suffers from rank deficiency. Here, we introduce supplementary structural similarity constraints in the gravity data cost function minimization.

We formulate these similarity constraints via cross gradients (e.g., [START_REF] Gallardo | Characterization of heterogeneous near-surface materials by joint 2d inversion of dc resistivity and seismic data[END_REF]) between density and V p models to obtain 3D density models that are consistent simultaneously with TFWI results and gravity data. In a similar work but more focused on offshore basin modelling, [START_REF] Welford | Crustal structure of Baffin Bay from constrained three-dimensional gravity inversion and deformable plate tectonic models[END_REF] performed 3D gravity modeling to obtain 3D density models constrained by structural models derived from seismic tomography. More recently, [START_REF] Darijani | Clustering and constrained inversion of seismic refraction and gravity data for overburden stripping: Application to uranium exploration in the athabasca basin, canada[END_REF] obtained 2D density models through a 2D gravity inversion constrained by a 2D seismic model at the near surface scale via a clustering/petrophysical constraint approach as in [START_REF] Lelievre | Joint inversion of seismic traveltimes and gravity data on unstructured grids with application to mineral exploration[END_REF]. In our study, we rely on a 3D model coming from TFWI and therefore we do not need such petrophysical/clustering constraints.

We introduce implicitly those constraints via a Nafe-Drake scaling relationship between density and V p .

Similar ideas have been proposed in other related works by [START_REF] Saunders | Constrained resistivity inversion using seismic data[END_REF] in the context of resistivity inversion constrained by seismic data, and by [START_REF] Brown | Seismically regularized controlled-source electromagnetic inversion[END_REF] and [START_REF] Yan | Two-dimensional magnetotelluric inversion using reflection seismic data as constraints and application in the COSC project[END_REF] for electromagnetic and magnetotelluric data inversion constrained by fixed seismic structures. Using linear relations between seismic velocites and densities or cross-gradient constraints between seismic velocities, densities and magnetic susceptibilities, [START_REF] Colombo | Geophysical modeling via simultaneous joint inversion of seismic, gravity, and electromagnetic data: Application to prestack depth imaging[END_REF], [START_REF] Stefano | Multiple-domain, simultaneous joint inversion of geophysical data with application to subsalt imaging[END_REF][START_REF] Moorkamp | A framework for 3-D joint inversion of MT, gravity and seismic refraction data[END_REF] performed joint inversion of seismic, gravity and electromagnetic data. In these earlier studies, seismic tomography relied essentially on ray-theory methods. [START_REF] Blom | Synthetic inversions for density using seismic and gravity data[END_REF] performed a joint full waveform/gravity data inversion for a global scale application, and found that additional seismic velocity-density correlations or structural constraints are needed to properly reconstruct density anomalies. [START_REF] Lin | The gramian method of joint inversion of the gravity gradiometry and seismic data[END_REF] used the joint seismic and gravity data inversion to image a synthetic salt-dome model. The method is constrained by a Gramian stabilizer applied to the gradient of seismic velocities and densities to increase the correlation between the different physical parameters of the inverse model or their spatial gradients. In Giraud et al. (2019b), gravity data inversions were constrained by geological structural uncertainty using entropy arguments.

This article is organized as follows. After introducing the geophysical and geological context of this study, we describe in a second part our inversion methodology and the optimization process. We then explain how to appropriately weight the different constraint terms involved in the cost function to be minimized and we show the models obtained after inversion. Then we present solutions of two kinds of gravity data inversions constrained by the V p model: a 2.5D gravity data inversion of a 1D Bouguer anomaly profile aligned with the seismic transect and a 3D inversion of the available gravity data in the 3D slice area. For the 2.5D gravity inversion, the 2.5D prior model is built by taking a 2D vertical section extracted from the 3D TFWI model and extended along either the direction perpendicular to the seismic profile or along the strike of the Pyrenees. As expected, we show that both the width and the strike of the anomalies influence the gravity response. However, the finite extent of the anomalies has a more dominant contribution than the orientation of the strike. All our 2.5D tests show that the lateral extension of the Mauleon's anomaly must be finite, with a width of about 30 km. We performed a 3D inversion to laterally localize the density anomalies, the 2.5D not being able to do it. The 2.5D inversions helped us to estimate the order of magnitude of the hyperparameters of regularization and to refine more quickly the range of the hyperparameters for the 3D inversions. Finally, we compare the different inverted models that we obtained in the 2.5D and 3D cases and propose a preliminary geological interpretation of our findings.

GEOLOGICAL AND GEOPHYSICAL CONTEXT

Geology of the area

The Pyrenees result from the collision between Iberia and Eurasia from late Cretaceous to early Miocene [START_REF] Lescoutre | Formation et réactivation du système de rift pyrénéo-cantabrique : héritage, segmentation et évolution thermique[END_REF][START_REF] Lescoutre | Role of rift-inheritance and segmentation for orogenic evolution: example from the Pyrenean-Cantabrian system[END_REF][START_REF] Gómez-Romeu | Role of rift structural inheritance in orogeny highlighted by the western pyrenees case-study[END_REF]. Positive complete Bouguer gravity anomalies characterize the northern flank of the Pyrenees, which are particularly strong in the Mauleon basin and close to the town of Saint-Gaudens (section A-A' in Figure 1-a). Understanding the nature and origin of these anomalies (from the upper mantle or from shallow crust structures) has been a long-lasting controversy, and was one of the main motivations of the PYROPE and IBERRARRAY seismic deployments from 2009 to 2014 [START_REF] Chevrot | High-resolution imaging of the Pyrenees and Massif Central from the data of the PYROPE and IBERARRAY portable array deployments[END_REF][START_REF] Chevrot | The Pyrenean architecture as revealed by teleseismic P-to-S converted waves recorded along two dense transects[END_REF], and of the later OROGEN deployments from 2015 to 2017 [START_REF] Chevrot | The non-cylindrical crustal architecture of the Pyrenees[END_REF]. Another motivation for studying the western part of the Pyrenees is that this region is seismically more active than the eastern part [START_REF] Souriau | A new synthesis of pyrenean seismicity and its tectonic implications[END_REF]. In particular it presents a pronounced concentration of seismicity beneath the northern Pyrenean foothills that is not observed in the eastern Pyrenees. This feature may be related to the strong structural asymmetry of the Pyrenees, with a continental subduction in the central and western Pyrenees that is absent in the eastern Pyrenees [START_REF] Chevrot | The non-cylindrical crustal architecture of the Pyrenees[END_REF]. Recently, FWI has been applied to teleseismic P waves recorded along the westernmost PYROPE profile, complemented with data recorded by nearby permanent stations. The resulting 3D V p model (Figure 2), which resolved the structures over a band of 60 km sourrounding the profile, revealed the presence of a shallow mantle body at 10 km depth below the Mauleon basin [START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P-waves[END_REF][START_REF] Martin | A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data -Application to the western pyrenees[END_REF][START_REF] Chevrot | The non-cylindrical crustal architecture of the Pyrenees[END_REF]). This mantle body, probably exhumed during the Cretaceous episode of rifting, explains the large positive Bouguer anomalies shown in Figure 1-b.

As mentioned in [START_REF] Martin | A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data -Application to the western pyrenees[END_REF], FWI can provide compressional and shear velocities as well as densities but in practice density is not really well reconstructed by FWI alone [START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P-waves[END_REF][START_REF] Beller | Lithospheric architecture of the South-Western Alps revealed by multiparameter teleseismic full-waveform inversion[END_REF]. Therefore, adding gravity data into the inversion can improve the density model reconstruction. Here, we aim at inverting a complete Bouguer gravity anomaly data set extracted from the BGI database under the constraint of a V p model obtained by inverting teleseismic P waveforms recorded by the westernmost PYROPE profile [START_REF] Chevrot | The Pyrenean architecture as revealed by teleseismic P-to-S converted waves recorded along two dense transects[END_REF].

Geophysical data and prior information for inversion

To perform gravity inversion constrained by compressional velocity models, we consider the dense (65062 points) complete Bouguer gravity anomaly dataset from the BGI database covering the French and Spanish continental regions (see the map of measured Bouguer gravity data in Figure 1a). These data have been used to obtain a 2 mn arc gridded global model of the complete Bouguer gravity anomalies [START_REF] Balmino | Spherical harmonic modelling to ultra high degree of Bouguer and isostatic anomalies[END_REF]Geodesist's-Handbook 2012) shown in Figure A1. For our study we prefered to use directly the measured data instead of the gridded data in order to make computations at the measurement points for maximum accuracy in the calculations. From this measured dataset, we extract a subset of 13000 points that are projected on a regular 2700 (135 × 20) points grid (340 km long × 50 km wide and rotated 45 o eastward) that includes the westernmost (A-A') 2.525 km resolution transect imaged by TFWI (See Figures 1b and1c). This seismic profile runs from coordinates (-2.615 • , 42.071 • ) to (0.32167 • , 44.4323 • ) and crosses the Ebro and Mauleon basins. Those 2700 gravity data points, located up to 25 km from the transect are used for 3D inversion and are also interpolated along this profile (136 points with a 2.525 km spacing) for 2.5D inversion. Along the profile, gravity anomalies show strong variations across the Pyrenees with anomalies, varying between -90 mgal in the south to +30 mgal in the north.

We now define the 3D a priori density model of the western Pyrenees that will be used in our Bouguer anomaly inversions. The density model is first constructed by extracting a V p volume around the seismic profile mentioned before. We scaled the a priori density model ρ using this V p model and using the Nafe-Drake law [START_REF] Nafe | Variation with depth in shallow and deep water marine sediments of porosity, density and the velocities of compressional and shear waves[END_REF][START_REF] Ludwig | Seismic Refraction, the Sea[END_REF][START_REF] Brocher | Empirical relations between elastic wavespeeds and density in the Earth's crust[END_REF] as in [START_REF] Martin | A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data -Application to the western pyrenees[END_REF]: 5). Furthermore, the computational mesh is distorted close to the surface to follow the topography.

ρ = 1.6612V p -0.4721V 2 p + 0.0671V 3 p -0.0043V 4 p + 0.000106V 5 p , for V p < 7.
The 2.5D models are built by extracting 2D V p and density sections from the 3D TFWI models along the transect and by replicating them invariantly to themselves over a 30 km (Figure 6-A) or a 44 km width . They are rotated with different angles θ with respect to the A-A' profile from 0 degree, which corresponds to the direction perpendicular to A-A' (line P-P') to 45 degrees which corresponds approximately to the strike of the Pyrenees (S-S' in Figure 4). In Figure 6, we show the influence of the angle of rotation and width of different 2.5D a priori models on the gravity response. For a 44 km width, the gravity anomalies are overestimated by a large amount for all angles of rotation (Figure 6-B), while for a 30 km width the gravity anomalies are in good agreement with the observed data. For a 30 km width and angles of 0 to 25 degrees, the amplitudes are very similar and are overestimated by 2-4 mgals while for 45 degrees they are underestimated by 2 mgals approximately.

This suggests that the strong anomaly produced by the shallow mantle body beneath the maximum gravity anomaly is oriented along a direction close to the strike of the Pyrenees (at an angle between 30 and 40 degrees from the P-P' transect). This width and those angles are close to those estimated ( 24 km width) at the the first order by [START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P-waves[END_REF]. However, note that for narrow density anomalies their orientation has a secondary effect compared to their lateral extent.

STRUCTURALLY CONSTRAINED INVERSION METHOD

Starting from a Nafe-Drake scaled density model as defined in the previous section, we explain now to perform the gravity data inversion under structural similarity constraint of densities by minimizing the cross-gradient between density and a fixed seismic velocity model obtained from TFWI.

Structural constraints

In our constrained inversion, structural correlations are taken into account between density anomalies ∆ρ and compressional velocites V p . We assume that at any point in the medium there exists, to first order, a spatially varying parameter β such that ∇∆ρ = β∇V p . Here, we consider the V p model obtained through TFWI as a suitable structural model describing the deep crustal architecture.

The structural constraints we apply here consist of minimizing the cross-gradient between density and velocity anomalies similarly to the joint resistivity-velocity inversion procedure used by [START_REF] Gallardo | Characterization of heterogeneous near-surface materials by joint 2d inversion of dc resistivity and seismic data[END_REF] and by many others since then [START_REF] Meju | Structural Coupling Approaches in Integrated Geophysical Imaging[END_REF]. Imposing a zero cross-gradient between two model distributions is equivalent to imposing local colinearity between the spatial gradients of the two models. The advantage of using the cross-gradient in our study is that we do not need to impose any predefined linear correlation function or any coefficient β in the inversion.

Inverse problem formulation and resolution

Let us now describe our methodology for gravity data inversion under seismic constraints. We minimize a cost function χ that involves the misfit between measured and computed complete Bouguer gravity anomalies, and regularize the problem by adding a damping of the norm of the density anomaly model ∆ρ and a structural constraint using the cross-gradient between density and V p models. The total cost function is thus

χ(m) =|| D obs -D cal || 2 l 2 +λ 2 || W (∆ρ -∆ρ p ) || 2 L 2 +α 2 || ∇∆ρ ∧ ∇V p || 2 L 2 , ( 3
)
where ∆ρ is the density anomaly, ∆ρ p the prior density anomaly, D obs and D cal the measured and theoretical Bouguer gravity data, and λ and α the regularization weights applied to the density model and the cross-gradients respectively. W is a depth weighting operator applied on the model and computed by integrating the sensitivity gravity kernel S over the whole mesh volume and for each observation point location [START_REF] Portniaguine | 3-D magnetic inversion with data compression and image focusing[END_REF]. For the density anomaly in cell k of the computational mesh, this matrix operator W is defined as

W kk = ( i=1,N obs S 2 i,k ) 1/2 , ( 4
)
where S ik is the element of the gravity kernel S linking the density anomaly ∆ρ k in cell k to the gravity anomaly

D i cal at location x i such that D i cal = S ik ∆ρ k .
The cross-gradient is discretized in space using a finite-difference scheme. In our Tomofast-x parallel inversion platform, we minimize Equation (3) using the LSQR algorithm of [START_REF] Paige | LSQR: An algorithm for sparse linear equations and sparse least squares[END_REF], which is implemented using a parallelized sparse matrix solver. At each inversion cycle (i.e., outer loop of LSQR inversion), the model is updated using an approximate solution using a fixed number of iterations in the inner loop of the LSQR algorithm [START_REF] Martin | Gravity inversion using wavelet-based compression on parallel hybrid CPU/GPU systems: Application to southwest Ghana[END_REF][START_REF] Martin | Parallel three-dimensional electrical capacitance data imaging using a nonlinear inversion algorithm and Lp norm-based model regularization[END_REF]Giraud et al. 2019a,b).

NUMERICAL RESULTS OF THE STRUCTURALLY CONSTRAINED INVERSION

Weighting of cost function terms in the 2.5D and full 3D inversion cases

To minimize the aggregate cost function (3) it is crucial to achieve optimum balancing between the different contributing terms. We determine the optimum weights λ and α following a L-curve approach (Hansen and O'Leary, 1993), which we extend to the two-parameter case by using a strategy similar in spirit to the Pareto optimization scheme [START_REF] Bijani | Physical-property-, lithologyand surface-geometry-based joint inversion using Pareto Multi-Objective Global Optimization[END_REF]).

Our analysis is based on the minimization of the data misfit cost function (term 1 in equation

(3), Figure 7-a), the model norm (term 2 in equation (3), Figure 7-d), cross-gradient norm (term 3 in equation ( 3), Figure 7-c) and the aggregate global cost function χ in equation ( 3), black solid curve in Figure 7-e). In Figure 7-b, we show the behavior of the cross-gradient and the model norms for a low data misfit of 10 -4 . We estimate the optimum (λ, α) values by identifying the inversion case fulfilling a series of conditions in terms of data misfit (Figure 7-a), model norm , and cross-gradient (Figure 7-c) values. The method has been applied to the four 2.5D a priori models rotated by 0, 20, 25 and 45 degrees respect to the P-P' direction. Similar inverted solutions have been obtained in the four cases and similar (λ, α) regularization parameters have been obtained. We describe now how those hyperparameters have been obtained for the 2.5D case (for 0 degrees for instance) and also for the 3D inversion and we present and compare the 2.5D and 3D results.

The process to determine λ and α is divided in three main stages. We first analyse a population of inverted models obtained from inversions using a range of values sampled on a discretized (λ, α) grid. To accelerate the search, we take a 2.5D a priori model averaged along the (A-A') transect that is extended invariantly to itself over 12 cells perpendicularly to the right and left of the profile. In our study, the (λ, α) grid consists of 11 x 11 values and is used to generate a collection of 121 inverted models. In a second stage, from this set of models, we then define a refined subrange of values thought to contain the optimum (λ, α) pair. In this refined subrange, here taken as [10 -8 , 5 10 -5 ] × [10 -5 , 10 -1 ], we refine the weights sampling around the optimal values. To this end, we perform 21 x 19 inversions on a refined set of (λ, α) pairs.

The estimate of the optimum set of weights (see Figure 7-a As shown in Figure A2, the data misfit function is well minimized in all cases, which suggests that the solution of the inverse problem is not unique. To obtain physical solutions at depth, as discussed in appendix B and observed in figure A3, it is thus crucial to perform the inversions considering all the terms of the total cost misfit function. We now focus our attention on the inverted solutions obtained under constraints. In Figure 9 we show the differences between inverted solutions obtained for the prior density model extracted from the TFWI model and extrapolated invariantly to itself in directions making an angle θ with the transect P-P' perpendicular to the A-A' profile. For θ=0, the model is invariant in the direction P-P' and for θ = 45 degrees it is invariant along the S-S' direction parallel to the strike of the Pyrenees passing through the highest Mauleon's gravity anomaly. As we can see, the differences are rather small, with maximum absolute differences and standard deviations of approximately 4 and 1.5 kg m 3 obtained for the 45 degrees case. The solutions are very close with a small but noticeable difference close to the surface. Again, these small differences result from the limited lateral extent of the anomalies as it was already the case for the responses of the four a priori models used in those 2.5D inversions.

For the 2.5D density modeling case, we first minimize the misfit cost function without the crossgradient term and we obtain a density distribution that is not well constrained at depth (Figure 8) as in [START_REF] Martin | A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data -Application to the western pyrenees[END_REF]. As already mentioned, instead of normalizing the sensitivity matrix by a depth weighting function, we use the integrated sensitivity technique of [START_REF] Portniaguine | 3-D magnetic inversion with data compression and image focusing[END_REF]. The differences between inverted and prior density models are shown in Figure 8. The solution is obviously not geologically relevant because the model perturbations are mainly concentrated close to the surface and decrease continuously with depth. To alleviate this, we take the same prior model and we now minimize the cost function with the cross-gradient constraint. This time, we obtain well constrained densities at depth with density perturbations that follow the spatial variations of the V p model (see the evolution of the cross-gradient in Figure 10). After 200 successive inversion cycles, the cross gradients are minimized in the first 60 km depth. In the first inversion cycles, the cross gradients are strong, but they progressively decrease and tend to very low values not only close to the surface but also at depth.

Note that the differences between a priori and inverted densities are as pronounced at depth as they are closer to surface . We note that this feature is not possible to achieve without geometrical/structural constraints. In Figure 7-e we can see the evolution of the different terms of the whole cost function that reaches a value of 10 -7 after 20 successive LSQR inversions (outer loop) with 300 iterations for each LSQR inversion. At the end of the inversion the cost function reaches a value of approximately 4 × 10 -8 and the data misfit has decreased by more than 80% (see Figure 7-e). The 2.5D inversion helps to tune the penalty parameters. However, it does not take into account correctly nor the localization of the lateral heterogeneities neither their orientation in the vicinity of the A-A' profile. It tends to concentrate the density perturbations along the central profile due to a different geometrical divergence of the gravity fields and also to decrease the density anomaly values below the location of strong gravity anomalies. It is thus always preferable to perform a full 3D inversion that allows to localize the density anomalies below strong gravity anomalies are present and also to better define their orientation.

Full 3D inversion case using topography

In order to have better images in the direction perpendicular to the profile, we perform gravity inversions for α and λ values (λ = 2 × 10 -7 and α = 10 -3 ) close to the optimal values obtained previously in the 2.5D case. As starting and a priori model we take the density model scaled with the 3D V p model obtained by FWI according to the Nafe-Drake law. The inner computational domain is the same as previously. However, a padding of 20 cells has been added to each outer vertical edge of the model to avoid edge effects. The 3D model has thus a total size of 175 × 60 × 30 grid points.

We perform 100 LSQR inversion cycles (outer loop) to achieve enough convergence of the whole inversion algorithm until reaching a total cost function value of 3.6 × 10 -7 . Gravity data are very well reproduced with point-to-point absolute errors lower than 0.4 mGal (4 × 10 -6 m/s 2 ), a total data misfit contribution of 5.46 × 10 -8 to the total misfit function, and a global relative error of 2.235 × 10 -3 .

Figure 12 shows the final model along the A-A' profile in a vertical section and also in a horizontal plane located at a 16.4 km depth (Figure 12-a), as well as in a section P-P' perpendicular to the profile (in the direction of the axis of the Pyrenees, Figure 12-d). Inverted densities are showing variations of +/-130 kg/m 3 slightly higher than to those obtained for the 2.5D inversion (+/-100 kg/m 3 ). Density anomalies close to the Moho are showing clear subduction of the Iberian plate beneath the European plate (negative density anomaly below 30 km depth) with part of the European Mantle rising up close to the surface (strong positive anomaly above the 30 km reference Moho depth). Indeed, in Figure 12, the new model emphasizes the subduction of the Iberian plate (represented by a strong negative density anomaly in blue) and the presence of exhumed mantle material (evidenced by a positive anomaly, high red values) reaching a 10 km depth. In the section perpendicular to the Pyrenean axis , the 3D inversion constrains the width of the exhumed mantle body to about 20 km in agreement with the 2.5D inversions [START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P-waves[END_REF]). However, the 3D full modeling provides a better reconstruction of the geometry of the anomalies in particular of the shallow mantle body, with better localization, orientation, and more accurate width. Furthermore, the exhumed anomaly can be located laterally, which is not possible by construction in the 2.5D configuration. This can be observed in Figure 12-a on the section located at a depth of 16.4 km, where the exhumed mantle anomaly is not exactly on the A-A' axis but is rotated with an angle of about 45 degrees with respect to P-P' direction.

The remaining differences between inverted and a priori models are mainly observed in the Iberian-European collisional domain. The initial cross-gradients are lower than those of the 2.5D case by an order of magnitude everywhere. Below 20 km depth, they are decreasing from about 4 × 10 -3 down to 1.5 × 10 -3 after inversion. At 5-10 km depth, cross-gradients are also reduced but remain important on top of the exhumed mantle body. In general, gravity data alone allow us to recover densities close to the surface but not at depth. But, thanks to the cross-gradient, we can also constrain deep density anomalies.

In Figure 13, a dispersion diagram between density and V p is shown for constrained and unconstrained solutions, and compared with the Nafe-Drake reference scaling. As in [START_REF] Martin | A high-order 3-D spectral-element method for the forward modelling and inversion of gravimetric data -Application to the western pyrenees[END_REF], the densities obtained without V p similarity constraints are overestimated and density values lower than 2550 kg/m 3 are almost absent close to the surface level and at depth. Moreover, adding the crossgradient constraint increases the amplitude of the variations of the densities compared to the densities obtained without cross-gradient constraints: it decreases further the densities where anomalies are negative close to the surface and also at depth until reaching the crust-mantle transition. This reduction has an amplitude ranging between -130 and +130 kg/m 3 . We now obtain lower densities in the crust and close to the surface in agreement with the results of [START_REF] Grandjean | Étude des structures crustales dans une portion de chaîne et de leurs relations avec les bassins sédimentaires, application aux Pyrénées occidentales[END_REF] in which densities are lower than 2500 kg/m 3 in the first 7 km in the basins.

GEOLOGICAL AND GEOMECHANICAL INTERPRETATION

Visually, the recovered density model bears important similarities with the tomographic model of [START_REF] Wang | The deep roots of the western Pyrenees revealed by full waveform inversion of teleseismic P-waves[END_REF], as can be seen in Figure 14. The model is also coherent with the recent collision models of Lescoutre ( 2019 We computed the Airy isostatic gravity anomalies by taking a reference Moho model obtained by assuming a 30 km compensation depth and topography interpolated form a ETOPO1 (1mn resolution).

Reference isostatic densities of 2670 kg/m 3 above this Moho and 3270 kg/m 3 beneath are used. In Figure 15, the observed isostatic anomalies are very similar to the isostatic anomalies computed in the final density model (i.e after inversion). The isostatic anomalies are reduced in the southern part of the profile compared to the Bouguer anomalies due to surface topography which predicts a deeper Moho (reaching around 36 km depth) that partly compensates the topography in that region, which is not the case in the northern part. The differences between observed and calculated anomalies are very small and reach maximum errors lower than 1 mGal (as can also been seen in Figure 12-d). This comparison demonstrates that it is equivalent to invert either Bouguer or isostatic anomalies.

CONCLUSIONS

We have performed 2.5D and full 3D gravity data inversions in the region showing the strongest gravity anomalies in the Pyrenees using structural constraints derived from a tomographic V p model to favor structural similarities between seismic velocity and density models. These constraints improved the density models at depth, by transferring structural information from TFWI to the density model, and providing structural details of the mantle wedge and the upper crust structures. Full 3D inversions with topography allow us to take into account lateral variations of densities and to localize heterogeneities at depth more accurately than with 2.5D inversions. 3D inversions also tend to better define their orientation (along the strike of the Pyrenees), whereas 2.5D inversions tend to decrease anomalies and concentrate them closer to the surface and along the central plane. The structural constraints adjust the densities not only close to the surface but also at depth close to the Moho with up to 10% density variations. However, 2.5D inversions allow us to tune quickly a range of penalty parameters of the cost function that can be used to estimate penalty hyperparameters for the full 3D inversions. These 3D inversions can be run in less than 20 minutes on 256 Intel Skylake 2.3GHz cores of the Olympe supercomputer. Inverting Bouguer or isostatic anomalies can be equivalently done and gives similar density models. They only differ by a shift generated by the removal of the reference Moho depth perturbation from the Bouguer anomalies.

In a future and second step, our methodology will be applied in 3D to the rest of the Pyrenees.

Future work will include joint cross-gradient inversion of gravity and seismic data at the scale of the whole Pyrenees. We will also consider adding constraints such as petrophysical and geological information as in [START_REF] Wehr | A three-dimensional model of the pyrenees and their foreland basins from geological and gravimetric data[END_REF]. In the present study, we performed inversions on medium sized supercomputing resources (50 to 200 processors clusters) but for the inversion of the whole Pyrenees area much more processors will be needed (at least 10 to 20 times more). In the mid-term, having access to the densities will provide key information on the mineralogic composition and thermal state of the lithosphere. To this end a full joint inversion seems to be a promising way to improve our capacity to image internal Earth structures by reducing model uncertainties. We can observe that the anomalies are overestimated when compared to the 30 km model width case with or without rotation of the model.

APPENDIX A: INVERSION ALGORITHM

A least-squares algorithm is used here to solve the matrix system as follows. Let us denote m the model vector containing the density anomalies ∆ρ. This model vector is obtained by solving the linear system of equation

     SW -1 λI d αW -1 ∇V p ∧ ∇      [W m] =      ∆g obs λW m p 0      , (A.1)
where m p = ∆ρ p is the a priori model and the gravity kernel S is depth-weighted by the operator W given in equation ( 4). The solution of (A.1) is obtained by an iterative LSQR algorithm [START_REF] Paige | LSQR: An algorithm for sparse linear equations and sparse least squares[END_REF]. At each i-th LSQR inversion cycle we solve the following linearized version of the system of equations (A.1)

For 1 ≤ i ≤ N iter max      SW -1 λI d αW -1 ∇V p ∧ ∇      [∆m i ] =      ∆g obs -∆g(m i-1 ) -λW (m i-1 -m p ) -α∇V p ∧ ∇m i-1      ∆m i = W -1 ∆m i m i = m i-1 + ∆m i End for , (A.2)
where N iter max is the maximum number of outer loop inversion cycles. Spatial gradients are computed using forward finite difference integration in the computation domain and using backward finite difference integration at the outer boundaries. Let us first consider the impact of an a priori model, a damping term on the model or a cross-gradient in the misfit function. As an example we consider the 2.5D model case (i.e for a 2D model defined as a Nafe-Drake scaled density model and extrapolated in the direction orthogonal to the A-A' profile).

As shown in Figures A2, the data misfit function is well minimized in all cases, which suggests that the solution of the inverse problem is not unique. Figure A3 shows the effective density models and the difference with the a priori model (scaled from the V p model using a Nafe-Drake law) or a density model of reference. In Figures 8 andA3, the inverted solutions are represented for inversions with and without starting models in the following cases :

• Case 1 : without damping (λ = 0) and without cross-gradient (α = 0)

• Case 2 : with damping (λ = 0) and without cross-gradient (α = 0)

• Case 3 : without damping (λ = 0) and with cross-gradient (α = 0)

where the non zero parameters λ (case 2) and α (case 3) are chosen for clarity as the optimal values (λ = 5 × 10 -7 and α = 2.3 × 10 -4 ) defined using the method described in the previous section.

Figure A2 shows the fit to gravity data obtained after inversion with or without starting model, with or without damping model regularization, and with or without cross-gradient constraints. Figure A3 shows the different density models obtained with the different types of inversions.

In Figures 8a-b and A3a-b, the case 1 is treated without or with a starting model. In the reconstructed model, the strongest density anomalies concentrate close to the surface where the highest gravity anomalies are present. This example illustrates that gravity inversion has a poor depth resolution. However, a better reconstruction of density anomalies close to the surface is achieved when a starting model is introduced (Figures A3c).

In case 2, densities are computed at depth with or without starting model and are exhibiting strong similarities (see Figures A3d-e). The solution is very close to the one obtained in case 1 when a starting model is introduced and no damping term is involved (see also Figure 8c). This result suggests that the damping term acts as a depth-dependent preconditionner of the matrix system (A.1) but the starting model does not prevent the inversion from concentrating the anomalies close to the surface.

In case 3 (Figure A3-f), the cross-gradient term adjusts the densities at depth to follow the variations of the V p model. But, as it is the case here, if no physical Nafe-Drake scaled starting model or if no damping term based on this prior model is introduced, these density variations are not following successfully the V p variations at depth, more specifically close to the crust-mantle transition (between 30 and 40 km depth). The inverted model obtained is smoothed too much everywhere in the computational domain. In Figure A4, an extreme case without data misfit minimization is considered while keeping the cross-gradient term only in the total misfit cost function. We can observe that the data distribution is reaching values closer to the real gravity data with a pronounced smoothing around the highest anomaly values.

In summary, these different tests suggest that more physical inverted physical models can be obtained during the inversion process by introducing a Nafe-Drake prior model as a starting model or into the damping model term, or both, and by applying simultaneously the cross-gradient term in order to distribute physically the density variations at depth, as can be seen in 8d-e. The simultaneous min-imization of the gravity data misfit, the model damping regularization term and the cross-gradient is necessary. Now, in the next section, we show how the different terms of the cost function are behaving when optimal λ and α are chosen.

  4 km/s (1) and ρ = ρ m = 3270kg/m 3 , for V p > 7.4 km/s (2) We define the density anomalies ∆ρ = ρ -ρ c in the crust (from the upper surface down to a reference Moho depth of 30 km) as the deviation from a reference crustal density ρ c = 2670 kg/m 3 , and ∆ρ = ρ -ρ m kg/m 3 beneath the Moho for a reference upper mantle density ρ m (at depths larger than 30 km). The 3D prior density model is shown Figure 3 with different cut sections at different depths (16.4 km, 50 km). Predicted Bouguer gravity anomalies (Figure 1-b) are relatively close to the measured Bouguer anomaly in the central and northern part of the profile. In our numerical tests, we have also observed that computed gravity responses are not influenced by density structures located farther than 25 km from the transect. As in Martin et al. (2017), a 340 km long × 50 km wide × 75 km deep 3D inner computational domain is considered and discretized over 136 × 20 × 30 prismatic bodies (2.525 km × 2.525 km × 2.533 km volume per cell). An external tapering is added to the inner computational domain (see Figure 5) in order to avoid edge effects. It involves 20 extra cells on each vertical outer interface of the computational domain. The 1 mn arc (≈ 1.78 km) ETOPO1 topography is added and interpolated at a resolution of 2.525 km resolution in longitude and latitude on top of the computational mesh. The gravity measurements are projected on this topography grid (Figures 4 and

  ) that balance the contribution of the different terms is obtained by identifying the case with lowest overall cost among the subset of inverted models honoring all the conditions (see Figure7-e). In these 2.5D inversions, we choose a damping model parameter λ = 3.5 × 10 -7 and a cross gradient weight α = 2.3 × 10 -4 . To ensure appropriate balancing between the different terms of the cost function, the optimum weights of the different terms have been inferred from the inflection points in the curves of data misfit curve versus model perturbations and cross-gradient norm (Figures7-b,c,d). We apply the L-curve analysis by performing many inversion runs with different λ and α values. An example is taken along the dashed lines shown in Figure7-e where the optimum α value is determined and fixed and an optimum λ value is determined accordingly. In a third stage, once the optimal values of the damping model parameter λ and of the cross gradient weight α have been found in 2.5D, we perform the inversions over the full 3D heterogeneous model with values of λ and α finely tuned around the values obtained in the previous 2.5D case.4.2 Cross-gradient constraints and their influence on density model reconstruction : the 2.5D case exampleTo show the impact of an a priori model and the different penalty terms of the misfit function on the inverted solutions, we performed different tests with or without damping or cross-gradient penalty in the misfit function. For the sake of clarity, the results of constrained and unconstrained inversions are described in appendix B for the 2.5D model example (i.e for an a priori 3D model defined as a Nafe-Drake scaled density model extrapolated in the direction P-P' orthogonal to the A-A' profile).

  ) (chapter III, Figure p.121, and Figure III-4, p. 131) and Lescoutre & Manatschal (2020); Gómez-Romeu et al. (2019). Our inversions show that the European mantle wedge and the Gavarnie nappe have strong imprint in the density model, with strong negative and positive density contrasts respectively (see Figure 12-c-left). The exhumed mantle body and the Gavarnie thrust show geometries in good agreement to the northern part of their transect passing through the Eastern part of Mauleon basin (profile A in Lescoutre (2019), Chapter III, Figure p.121). The southernmost part of our model is more similar to the geological profile passing through the west of Mauleon basin (profile B in Chapter III, p.121 Lescoutre (2019), reproduced in appendix A.5). Again, this simply results from the strong obliquity of the seismic transect with respect to the strike of the Pyrenees. Our 3D imaged section corresponds also to the profile passing through the Ebro basin and eastern Mauleon (section in of Figure III-4 of Lescoutre (2019), p. 131) and reproduced also in Figure A.5. In short, the AA' model profile (Figure 1) can be seen as an intermediate stage between these two south-north transects : one passing through the Eastern Mauleon basin and another one passing through Saint-Jean de Luz basin and Ebro basin at the west of Mauleon basin, with a central part where the geometry of the mantle body is in good agreement with Wang et al. (2016).

Figure 1 .

 1 Figure 1. (a) Measured Bouguer gravity anomaly data in Southwest France and North Spain. The western transect seismic profle A-A' crossing the Mauleon and Ebro basin in the Pyrenees close to Labour town (Chevrot et al. 2018) is shown in blue. (b-middle) Observed (black solid) and a priori (blue solid) Bouguer gravity anomalies close to the westernmost transect A-A'. (c) V p and Nafe-Drake scaled prior density model r ho along the A-A' profile.

Figure 2 .

 2 Figure 2. 3D V p model along the A-A' profile at 16.4 km (left column) and 50 km depth (right column). Sections perpendicular to the A-A' profile are shown too at x = 200 km (top), 226 km (middle, where the exhumed mantle is going up close to the surface) and 290 km.

Figure 3 .

 3 Figure 3. 3D Nafe-Drake scaled prior density model ρ along the A-A' profile at 16.4 km (left column) and km depth (right column). Sections perpendicular to the A-A' profile are shown too at x = 200 km (top), km (middle, where the exhumed mantle is going up close to the surface) and 290 km.

Figure 4 .

 4 Figure 4. Top view of the measured/not gridded data is shown (left) as well as the data interpolated (right) in the region of interest (inside the black box) containing the (A-A') seismic profile (red line). The interpolated data are gridded at 2.471 km resolution. The line segment P-P' perpendicular to the seismic profile and the Pyrenean strike segment S-S' passing through the Mauleon's gravity anomaly are shown. Both segments make an angle θ around 40 degrees.

Figure 5 .Figure 6 .

 56 Figure 5. Bouguer gravity anomaly interpolated on the topography in the computational slice of interest.ETOPO1 1 mn arc resolution topography interpolated at 2.471 km resolution is considered in the computational area of interest that includes the (A-A') seismic profile. The computational mesh used for data inversion takes into account the topography. The first layer of the computational mesh is distorted, follows the topography and is also shown on the figure. Lateral padding of the computational model is added, taking into account the topography outside the set of gravity data (in color).

Figure 7 .

 7 Figure 7. Data misfits are normalized by the sum of squared data. (a-Top) 3D representation of the three terms of the cost function. The optimal regularization parameters are indicated by the black dot. (b-Top) Section of this surface passing through this point. (c and d : Middle left and right) L curves to define the optimal model regularization and cross-gradient weights respectively corresponding to the optimal model damping and crossgradient terms shown in Figure (a). Refined subrange of weighting parameters (λ,α) varying in the set intervals [10 -8 ,5 10 -5 ] × [10 -5 ,10 -1 ] has been used. (e-Bottom) Evolution of data misfit, cross-gradient and model damping norms through inversion iterations.

Figure 8 .Figure 9 .

 89 Figure 8. Difference between inverted and prior density models in the 2.5D ('quasi 3D') case. (Top) Without cross-gradient constraint, and without (a) or with (b) starting model. (Middle) With cross-gradient constraint, and without (c) or with (d) starting model. Strong differences are appearing at depth, with variations exceeding + or -100 kg/m 3 , particularly when cross-gradient constraint is applied. Note also that differences are predominantly vertical. (Bottom-e) Observed and inverted gravity data with and without cross-gradient constraint.

Figure 10 .

 10 Figure 10. (Left) From top to bottom : cross-gradient evolution for 5, 10, 50, 100 and 200 inversion loops. Cross-gradient is decreasing as expected. (Right) Inverted effective densities with cross-gradient constraints.

Figure 11 .Figure 12 .

 1112 Figure 11. 3D inverted density model ρ along the A-A' profile at 16.4 km (left column) and 50 km depth (right column). Sections perpendicular to the A-A' profile are shown too at x = 200 km (top), 226 km (middle, where the exhumed mantle is going up close to the surface) and 290 km.

Figure 13 .Figure 14 .Figure 15 .

 131415 Figure 13. Dispersion diagram of ρ/V p properties around the initial ρ/V p Nafe-Drake scaling relation with and without seismic constraints. The correlation between inverted density and V p models obtained after gravity data inversion constrained by seismics is shown in red (blue for the non-constrained inversion case). When compared to the non-constrained densities, the seismic constrained densities are lowered close to the surface and also at depth, in the northern part of the section and in the diping part of the Iberian crust.

  APPENDIX B: INFLUENCE OF INITIAL MODELS AND PENALTY CONSTRAINTS ON DENSITY MODEL RECONSTRUCTION : DATA AND MODELS WITH OR WITHOUT STARTING MODELS, DAMPING REGULARIZATION AND CROSS-GRADIENT CONSTRAINTS.

Figure A2 .Figure A3 .Figure A4 .

 A2A3A4 Figure A2. Observed and inverted gravity data with (1) or without (0) starting model, with or without model damping regularization, with or without cross-gradient constraint.

ACKNOWLEDGMENTS

We thank the CALMIP supercomputing mesocentre of Toulouse, France, for their support through project P1138 on the Olympe supercomputer platform, and the TGCC/GENCI supercomputing center through project gen6351 on the CURIE supercomputer. We are grateful to the CNES granted projects 2017 and 2018-TOSCA/GET-GRAVI-GOCE for their financial support. This study was also supported by the Orogen research project, a tripartite academic-industry collaboration between CNRS, BRGM and Total. Jeremie Giraud was also supported by an Australian Government Research Training Program (RTP) Scholarship when part of this work was conducted.


Quentin Bertrand

Mathurin Massias

Anderson acceleration of coordinate descent

come   L'archive ouverte pluridisciplinaire

Introduction

Gradient descent is the workhorse of modern convex optimization [START_REF] Nesterov | Introductory lectures on convex optimization[END_REF][START_REF] Beck | First-Order Methods in Optimization[END_REF]. For composite problems, proximal gradient descent retains the nice properties enjoyed by the latter. In both techniques, inertial acceleration achieves optimal convergence rates [START_REF] Nesterov | A method for solving a convex programming problem with rate of convergence O(1/k 2 ). Soviet Math[END_REF][START_REF] Beck | A fast iterative shrinkagethresholding algorithm for linear inverse problems[END_REF].

Coordinate descent is a variant of gradient descent, which updates the iterates one coordinate at a time [START_REF] Tseng | A coordinate gradient descent method for nonsmooth separable minimization[END_REF][START_REF] Friedman | Regularization paths for generalized linear models via coordinate descent[END_REF]. Proximal coordinate descent has been applied to numerous Machine Learning problems (Shalev-Shwartz and [START_REF] Shalev-Shwartz | Accelerated proximal stochastic dual coordinate ascent for regularized loss minimization[END_REF][START_REF] Wright | Coordinate descent algorithms[END_REF][START_REF] Shi | A primer on coordinate descent algorithms[END_REF], in particular the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] or sparse logistic regression [START_REF] Ng | Feature selection, l1 vs. l2 regularization, and rotational invariance[END_REF]. It is used in preeminent packages such as scikitlearn [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF], glmnet (Friedman et al., ) -f (x * )

Figure 1: Suboptimality along iterations on the 1000 first features of the rcv1 dataset for a quadratic problem. GD: gradient descent, CD: coordinate descent.

2009), libsvm [START_REF] Fan | Liblinear: A library for large linear classification[END_REF] or lightning [START_REF] Blondel | Lightning: largescale linear classification, regression and ranking in python[END_REF]. On the theoretical side, inertial accelerated versions of coordinate descent [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF][START_REF] Lin | An accelerated proximal coordinate gradient method[END_REF][START_REF] Fercoq | Accelerated, parallel, and proximal coordinate descent[END_REF] achieves optimal rates. To obtain optimal rates, Anderson acceleration [START_REF] Anderson | Iterative procedures for nonlinear integral equations[END_REF]) is an alternative to inertia: it provides acceleration by exploiting the iterates' structure. This procedure has been known for a long time, under various names and variants [START_REF] Wynn | Acceleration techniques for iterated vector and matrix problems[END_REF][START_REF] Eddy | Extrapolating to the limit of a vector sequence[END_REF][START_REF] Smith | Extrapolation methods for vector sequences[END_REF], see [START_REF] Sidi | Vector extrapolation methods with applications[END_REF]; [START_REF] Brezinski | Shanks sequence transformations and anderson acceleration[END_REF] for reviews. Anderson acceleration enjoys accelerated rates on quadratic functions [START_REF] Golub | Chebyshev semi-iterative methods, successive overrelaxation iterative methods, and second order richardson iterative methods[END_REF], but theoretical guarantees in the nonquadratic case are weaker [START_REF] Scieur | Regularized nonlinear acceleration[END_REF]. Interestingly, numerical performances still show significant improvements on nonquadratic objectives. Anderson acceleration has been adapted to various algorithms such as Douglas-Rachford [START_REF] Fu | Anderson accelerated Douglas-Rachford splitting[END_REF], ADMM [START_REF] Poon | Trajectory of alternating direction method of multipliers and adaptive acceleration[END_REF] or proximal gradient descent [START_REF] Zhang | Globally convergent type-I Anderson acceleration for non-smooth fixed-point iterations[END_REF][START_REF] Mai | Anderson acceleration of proximal gradient methods[END_REF][START_REF] Poon | Geometry of first-order methods and adaptive acceleration[END_REF]. Among main benefits, the practical version of Algorithm 1 Offline Anderson extrapolation init: x (0) ∈ R p 1 for k = 1, . . . do 2 x (k) = T x (k-1) + b // regular linear iteration 3 U = [x (1)x (0) , . . . , x (k)x (k-1) ]

4 c = (U U ) -1 1 k /1 k (U U ) -1 1 k ∈ R k 5 x (k) e-off = k i=1 c i x (i)
// does not affect x (k) 6 return x

(k) e-off
Anderson acceleration is memory efficient, easy to implement, line search free, has a low cost per iteration and does not require knowledge of the strong convexity constant. Finally, it introduces a single additional parameter, which often does not require tuning (see Section 3.1).

In this work:

• We propose an Anderson acceleration scheme for coordinate descent, which, as visible on Figure 1, outperforms inertial and extrapolated gradient descent, as well as inertial coordinate descent, even reaching the first order optimal performance of conjugate gradient on this exemple.

• The acceleration is obtained eventhough the iteration matrix is not symmetric, a notable problem in the analysis of Anderson extrapolation.

• We empirically highlight that the proposed acceleration technique can generalize in the nonquadratic case (Algorithm 3) and significantly improve proximal coordinate descent algorithms (Section 3), which are state-of-the-art first order methods on the considered problems.

Notation The j-th line of the matrix A is A j: and its j-th column is A :j . The canonical basis vectors of R p are e j . The vector of size K with all one entries is 1 K . The spectral radius of the matrix A, ρ(A), is the largest eigenvalue modulus of A. The set of p by p symmetric semidefinite positive matrices is S p + . The condition number κ(A) of a matrix A is its largest singular value divided by its smallest. A positive definite matrix A induces the norm x A = √ x Ax. The proximity operator of the function g is prox g (x) = arg min y g(y) + 1 2 xy 2 .

Anderson extrapolation

Background

: x (0) ∈ R p 1 for k = 1, . . . do 2 x (k) = T x (k-1) + b // regular iteration 3 if k = 0 mod K then 4 U = [x (k-K+1) -x (k-K) , . . . , x (k) -x (k-1) ] 5 c = (U U ) -1 1 K /1 K (U U ) -1 1 K ∈ R K 6 x (k) e-on = K i=1 c i x (k-K+i) 7 x (k) = x (k) e-on
// base sequence changes 8 return x (k) iterations, that is:

x (k+1) = T x (k) + b , (1)
where the iteration matrix T ∈ R p×p has spectral radius ρ(T ) < 1. There exist two variants: offline and online, which we recall briefly.

Offline extrapolation (Algorithm 1), at iteration k, looks for a fixed point as an affine combination of the k first iterates: 1) , and solves for the coefficients c (k) ∈ R k as follows:

x (k) e-off = k 1 c (k) i x (i-
c (k) = arg min k 1 ci=1 k 1 c i x (i-1) -T k 1 c i x (i-1) -b 2 = arg min k 1 ci=1 k 1 c i x (i) -x (i-1) 2 = (U U ) -1 1 k /1 k (U U ) -1 1 k , (2)
where U = [x (1)x (0) , . . . , x (k)x (k-1) ] ∈ R p×k (and hence the objective rewrites U c 2 ). In practice, since -1) . The motivation for introducing the coefficients c (k) is discussed in more depth after Prop. 6 in [START_REF] Massias | Dual extrapolation for sparse generalized linear models[END_REF], and details about the closed-form solution can be found in Scieur et al. (2016, Lem. 2.4). In offline acceleration, more and more base iterates are used to produce the extrapolated point, but the extrapolation sequence does not affect the base sequence. This may not scale well since it requires solving larger and larger linear systems.

x (k) is available when c (k) is computed, one uses x (k) e = k 1 c (k) i x (i) instead of k 1 c (k) i x (i
A more practical variant is the online version (Algorithm 2), considered in this paper. The number of points to be extrapolated is fixed to K; x (1) , . . . , x (K) are computed normally with the fixed point iterations, but x (K) e is computed by extrapolating the iterates from x (1) to x (K) , and x (K) is taken equal to x (K) e . K normal iterates are then computed from x (K+1) to x (2K) then extrapolation is performed on these last K iterates, etc.

As we recall below, results on Anderson acceleration mainly concern fixed-point iterations with symmetric iteration matrices T , and results concerning non-symmetric iteration matrices are weaker [START_REF] Bollapragada | Nonlinear acceleration of momentum and primal-dual algorithms[END_REF]. Poon and Liang (2020, Thm 6.4) do not assume that T is symmetric, but only diagonalizable, which is still a strong requirement.

Proposition 1 (Symmetric T , Scieur 2019). Let the iteration matrix T be symmetric semi-definite positive, with spectral radius ρ = ρ(T ) < 1. Let x * be the limit of the sequence (x (k) ).

Let ζ = (1 - √ 1 -ρ)/(1 + √ 1 -ρ).
Then the iterates of offline Anderson acceleration satisfy, with B = (Id -T ) 2 :

x (k) e-off -x * B ≤ 2ζ k-1 1+ζ 2(k-1) x (0) -x * B , (3)
and thus those of online extrapolation satisfy:

x (k) e-on -x * B ≤ 2ζ K-1 1+ζ 2(K-1) k/K x (0) -x * B .
(4) [START_REF] Scieur | Regularized nonlinear acceleration[END_REF] showed that the offline version in Proposition 1 matches the accelerated rate of the conjugate gradient [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF]. As it states, gradient descent can be accelerated by Anderson extrapolation on quadratics.

Application to least squares The canonical application of Anderson extrapolation is gradient descent on least squares. Consider a quadratic problem, with b ∈ R p , H ∈ S p ++ such that 0 ≺ H L and L > 0:

x * = arg min

x∈R p 1 2 x Hx + b, x . (5)
A typical instance is overdetermined least squares with full-column rank design matrix A ∈ R n×p , and observations y ∈ R n , such that H = A A and b = -A y.

On Problem (5) gradient descent with step size 1/L reads:

x (k+1) = Id p -1 L H T GD ∈S p + x (k) + (-b/L b GD ) . (6)
Because they have this linear structure, iterates of gradient descent can benefit from Anderson acceleration, observing that the fixed point of x → T GD x + b GD solves (5), with T GD ∈ S p + . Anderson acceleration of gradient descent has therefore been well-studied beyond the scope of Machine Learning [START_REF] Pulay | Convergence acceleration of iterative sequences. the case of scf iteration[END_REF][START_REF] Eyert | A comparative study on methods for convergence acceleration of iterative vector sequences[END_REF]. However, on many Machine Learning problems, coordinate descent achieves far superior performance, and it is interesting to determine whether or not it can also benefit from Anderson extrapolation.

Linear iterations of coordinate descent

To apply Anderson acceleration to coordinate descent, we need to show that its iterates satisfy linear iterations as in ( 6). An epoch of cyclic coordinate descent for Problem (5) consists in updating the vector x one coordinate at a time, sequentially, i.e. for j = 1, . . . , p:

x j ← x j - 1 H jj (H j: x + b j ) , (7)
which can be rewritten, for j = 1, . . . , p:

x ← Id p - e j e j H jj H x - b j H jj e j . (8)
Thus, for primal iterates, as observed by Bertrand et al. (2020, Sec. A.3), one full pass (updating coordinates from 1 to p) leads to a linear iteration:

x (k+1) = T CD x (k) + b CD , (9)
with

T CD = Id p - epe p Hpp H . . . Id p - e1e 1
H11 H . Note that in the case of coordinate descent we write x (k) for the iterates after one pass of coordinate descent on all features, and not after each update (7). The iterates of coordinate therefore also have a fixed-point structure, but contrary to gradient descent, their iteration matrix T CD is not symmetric, which we address in Section 2.3.

Anderson extrapolation for nonsymmetric iteration matrices

Even on quadratics, Anderson acceleration with nonsymmetric iteration matrices is less developed, and the only results concerning its theoretical acceleration are recent and weaker than in the symmetric case.

Proposition 2 (Bollapragada et al. 2018, Thm 2.2). When T is not symmetric, and ρ(T ) < 1,

x (k) e-off -T x (k) e-off -b ≤ Id -ρ(T -Id) 2 P * (T )(x (1) -x (0) ) ,
where the unavailable polynomial P * minimizes P (T )(x (1)x (0) ) amongst all polynomials P of degree exactly k -1 whose coefficients sum to 1.

The quality of the bound (in particular, its eventual convergence to 0) crucially depends on P (T ) . Using the Crouzeix conjecture [START_REF] Crouzeix | Bounds for analytical functions of matrices[END_REF][START_REF] Bollapragada | Nonlinear acceleration of momentum and primal-dual algorithms[END_REF] managed to bound P (T ) , with P a polynomial:

P (T ) ≤ c max z∈W (T ) |P (z)| , (10)
with c ≥ 2 [START_REF] Crouzeix | Numerical range and functional calculus in hilbert space[END_REF][START_REF] Crouzeix | The numerical range is a (1+2)-spectral set[END_REF], and W (T ) the numerical range:

W (T ) {x * T x : x 2 = 1, x ∈ C p } . ( 11
)
Since there is no general formula for this bound, Bollapragada et al. ( 2018) used numerical bounds on

W (T 1 ) W (T 128 ) W (T 256 ) W (T 512 ) -2 -1 0 1 2 -2 -1 0 1 2 gina agnostic -2 -1 0 1 2 -1 0 1 rcv1 -2 -1 0 1 2 -1 0 1 real-sim -2-1 0 1 2 -2 -1 0 1 2 news20
Figure 2: Numerical range of T q as q varies; T is the iteration matrix of Ridge regression problems with conditioning κ = 10 3 , on 4 datasets. The black cross marks the (1, 0) point, which should lie outside the range for the theoretical bound to be useful.

W (T q ) to ensure convergence. Figure 2 displays the numerical range W (T q ) in the complex plane for q ∈ {1, 128, 256, 512}. In order to be able to apply the theoretical result from [START_REF] Bollapragada | Nonlinear acceleration of momentum and primal-dual algorithms[END_REF], one must chose q such that the point (1, 0) is not contained in W (T q ), and extrapolate x (0) , x (q) , x (2q) , . . . One can see on Figure 2 that large values of q are needed, unusable in practice: q = 512 is greater than the number of iterations needed to converge on some problems. Moreover, Anderson acceleration seems to provide speed up on coordinate descent even with q = 1 as we perform, which highlights the need for refined bounds for Anderson acceleration on nonsymmetric matrices.

We propose two means to fix this lack of theoretical results: to modify the algorithm in order to have a more amenable iteration matrix (Section 2.4), or to perform a simple cost function decrease check (Section 2.5).

Pseudo-symmetrization of T

A first idea to make coordinate descent theoretically amenable to extrapolation is to perform updates of coefficients from indices 1 to p, followed by a reversed pass from p to 1. This leads to an iteration matrix which is not symmetric either but friendlier: it writes with

T CD-sym H -1/2 SH 1/2 , ( 12
)
S = Id p -H 1/2 e1e 1 H11 H 1 2 × • • • × Id p -H 1 2 epe p Hpp H 1 2 × Id p -H 1 2 epe p Hpp H 1 2 × • • • × Id p -H 1 2 e1e 1 H11 H 1 2 . ( 13
)
S is symmetric, thus, S and T (which has the same eigenvalues as S), are diagonalisable with real eigenvalues. We call these iterations pseudo-symmetric, and show that this structure allows to preserve the guarantees of Anderson extrapolation.

Proposition 3 (Pseudosym. T = H -1/2 SH 1/2 ).

Let T be the iteration matrix of pseudo-symmetric coordinate descent: T = H -1/2 SH 1/2 , with S the symmetric semi-definite positive matrix of (12). Let x * be the limit of the sequence

(x (k) ). Let ζ = (1 - √ 1 -ρ)/(1 + √ 1 -ρ).
Then ρ = ρ(T ) = ρ(S) < 1 and the iterates of offline Anderson acceleration satisfy:

x (k) e-off -x * B ≤ κ(H) 2ζ k-1 1+ζ 2(k-1) x (0) -x * B , (14)
and thus those of online extrapolation satisfy:

x (k) e-on -x * B ≤ κ(H) 2ζ K-1 1+ζ 2(K-1) k/K x (0) -x * B . (15)
Proof of Proposition 3 can be found in Appendix B. Proposition 3 shows accelerated convergence rates for the offline Anderson acceleration, but a κ(H) appears in the rate of the online Anderson acceleration, meaning that K must be large enough that ζ K mitigates this effect. This factor however seems like a theoretical artefact of the proof, since we observed significant speed up of the online Anderson acceleration, even with bad conditioning of H (see Figure 3).

Algorithm 3 Online Anderson PCD (proposed) init: x (0) ∈ R p for k = 1, . . . do x = x (k-1) for j = 1, . . . p do xj = x j x j = prox λ L j gj (x j -A :j ∇f (Ax)/L j ) Ax += (x j -xj )A :j x (k) = x // regular iter. O(np) if k = 0 mod K then // extrapol., O(K 3 + pK 2 ) U = [x (k-K+1) -x (k-K) , . . . , x (k) -x (k-1) ] c = (U U ) -1 1 K /1 K (U U ) -1 1 K ∈ R K x e = K i=1 c i x (k-K+i) if f (Ax e ) + λg(x e ) ≤ f (x (k) ) + λg(x (k) ) then 13 x (k) = x e // guaranteed convergence return x (k)
Figure 3 illustrates the convergence speed of cyclic and pseudo-symmetric coordinate descent on the rcv1 dataset. Anderson acceleration provides speed up for both versions. Interestingly, on this quadratic problem, the non extrapolated pseudo-symmetric iterations perform poorly, worse than cyclic coordinate descent. However, the performances are reversed for their extrapolated counterparts: the pseudo-symmetrized version is better than the cyclic one (which has a nonsymmetric iteration matrix). Finally, Anderson extrapolation on the pseudo-symmetrized version even reaches the conjugate gradient performance.

Generalization to nonquadratic and proposed algorithm

After devising and illustrating an Anderson extrapolated coordinate descent procedure for a simple quadratic objective, our goal is to apply Anderson acceleration on problems where coordinate descent achieve state-of-the-art results, i.e., of the form:

min x∈R p f (Ax) + λg(x) := f (Ax) + λ p j=1 g j (x j ) , ( 16
)
where f : R n → R is convex, γ-smooth and g j 's are proper, closed and convex functions As examples, we allow g = 0, g = x 1 , g = 1 2 x 2 2 , g = x 1 + ρ 2λ x 2 . One pass of proximal coordinate descent from 1 to p can be seen as a nonlinear fixed point iteration:

x (k+1) = ψ(x (k) ) . ( 17
)
Proposition 4. If f is convex and smooth and C 2 , g j are convex smooth and C 2 , then ψ is differentiable and Therefore, iterations of proximal coordinate descent for this problem lead to noisy linear iterations. Proof of Proposition 4 can be found in Appendix B. Figure 4 shows the performance of Anderson extrapolation on a 2 -regularised logistic regression problem:

x (k+1) = Dψ(x * )(x (k) -x * ) + x * + o( x (k) -x * ) .
arg min x∈R p n i=1 log(1 + e -yiAi:x ) + λ 2 x 2 2 . (18)
One can see that despite the better theoretical properties of the pseudo-symmetrized coordinate descent, Anderson acceleration on coordinate descent seems to work better on the cyclic coordinate descent. We thus choose to apply Anderson extrapolation on the cyclic coordinate descent (Algorithm 3), while adding a step checking the decrease of the objective function in order to ensure convergence. Finally, we can also use Algorithm 3 in the non smooth case where g = • 1 , since coordinate descent achieves support identification when the solution is unique, after which the objective becomes differentiable. There is therefore a linear structure after a sufficient number of iterations (Massias et al., 2019, Prop. 10).

Experiments

An implementation relying on numpy, numba and cython [START_REF] Harris | Array programming with numpy[END_REF][START_REF] Lam | Numba: A llvm-based python jit compiler[END_REF][START_REF] Behnel | Cython: The best of both worlds[END_REF], with scripts to reproduce the figures, is available at https://mathurinm.github.io/ andersoncd We first show how we set the hyperparameters of Anderson extrapolation (Section 3.1). Then we show that Anderson extrapolation applied to proximal coordinate descent outperforms other first order algorithms on standard Machine Learning problems (Section 3.2).

Parameter setting

Anderson extrapolation relies on 2 hyperparameters: the number of extrapolated points K, and the amount of regularization eventually used when solving the linear system to obtain the coefficients c ∈ R K . Based on the conclusions of this section, we fix these parameters for all the subsequent experiments in Section 3.2: no regularization and K = 5.

Influence of the regularization. Scieur et al.

(2016) provided accelerated complexity rates for regularized Anderson extrapolation: a term λ reg c 2 is added to the objective of Equation ( 2). The closedform formula for the coefficients is then

(U U + λ reg Id K ) -1 1 K /1 K (U U + λ reg Id K ) -1 1 K .
However, similarly to Mai and Johansson (2019) and [START_REF] Poon | Geometry of first-order methods and adaptive acceleration[END_REF] we observed that regularizing the linear system does not seem necessary, and can even hurt the convergence speed. Figure 5 shows the influence of the regularization parameter on the convergence on the rcv1 dataset for a sparse logistic regression problem, with K = 5 and λ = λ max /30. The more the optimization problem is regularized, the more the convergence speed is deteriorated. We observed the same phenomenon when solving least squares problems (Figure 10, Appendix A.1). Thus we choose not to regularize when solving the linear system for the extrapolation coefficients. We simply check if the extrapolated point yields a lower objective function than the current regular iterate (see Algorithm 3).

Influence of K. Figure 6 shows the impact of K on the convergence speed. Although the performance depends on K, it seems that the dependency is loose, as for K ∈ {5, 10, 20} the acceleration is roughly the same. Therefore, we do not treat K as a parameter and fix it to K = 5. Every K iterations Anderson accelerated algorithms require to solve a K × K linear system. For K = 5 it is marginal compared to a gradient call: i.e., 5 3 + 5 2 p np in our settings.

Numerical comparison on Machine Learning problems

We compare multiple algorithms to solve popular Machine Learning problems: the Lasso, the elastic net, and sparse logistic regression (experiments on group lasso are in Appendix A.5). The compared algorithms are the following: proximal gradient descent (PGD, Combettes and Wajs 2005), Nesterov-like inertial PGD (FISTA, Beck and Teboulle 2009), Anderson accelerated PGD [START_REF] Mai | Anderson acceleration of proximal gradient methods[END_REF][START_REF] Poon | Geometry of first-order methods and adaptive acceleration[END_REF], proximal coordinate descent (PCD, Tseng and Yun 2009), inertial PCD [START_REF] Lin | An accelerated proximal coordinate gradient method[END_REF][START_REF] Fercoq | Accelerated, parallel, and proximal coordinate descent[END_REF], Anderson accelerated PCD (ours). We use datasets from libsvm [START_REF] Fan | Liblinear: A library for large linear classification[END_REF] and openml [START_REF] Feurer | Openml-python: an extensible python api for openml[END_REF] (Table 1), varying as much as possible to demonstrate the versatility of our approach. We also vary the convergence metric: we use suboptimality in the main paper, while graphs measuring the duality gaps are in Appendix A.

Lasso. Figure 7 shows the suboptimality f (x (k) )f (x * ) of the algorithms on the Lasso problem:

x * = arg min

x∈R p 1 2 y -Ax 2 + λ x 1 , (19)
as a function of the number of iterations for multiple datasets and values of λ. We parametrize λ as a fraction of λ max = A y ∞ , smallest regularization strength for which x * = 0. superiority of proximal coordinate descent over proximal gradient descent for Lasso problems on real-world datasets, and the benefits of extrapolation for coordinate descent. It shows that Anderson extrapolation can lead to a significative gain of performance. In particular Figure 7 shows that without restart, inertial coordinate descent [START_REF] Lin | An accelerated proximal coordinate gradient method[END_REF][START_REF] Fercoq | Accelerated, parallel, and proximal coordinate descent[END_REF] can slow down the convergence, despite its accelerated rate. Note that the smaller the value of λ, the harder the optimization: when λ decreases, more iterations are needed to reach a fixed suboptimality. The smaller λ is (i.e., the harder the problem), the more efficient Anderson extrapolation is.

As in [START_REF] Mai | Anderson acceleration of proximal gradient methods[END_REF], Anderson PGD performs well when p < n: Anderson PGD outperforms Anderson PCD on the gina agnostic dataset, where p < n. On all the other datasets, especially when p > n and when the values of λ are large, An-derson PCD outperforms Anderson PGD.

Other convergence metrics can be considered, since f (x * ) is unknown to the practitioner: for the Lasso, it is also common to use the duality gap as a stopping criterion [START_REF] Massias | Celer: a fast solver for the lasso with dual extrapolation[END_REF]. Thus, for completeness, we provide Figure 11 in appendix, which shows the duality gap as a function of the number of iterations. With this metric of convergence, Anderson PCD also significantly outperforms its competitors.

Elastic net. Anderson extrapolation is easy to extend to other estimators than the Lasso. Figure 8 (and Figure 12 in appendix) show the superiority of the Anderson extrapolation approach over proximal gradient descent and its accelerated version for the elastic net problem [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF]:

arg min x∈R p 1 2n y -Ax 2 + λ x 1 + ρ 2 x 2 2 . (20)
In particular, we observe that the more difficult the problem, the more useful the Anderson extrapolation: it is visible on Figures 8 and 12 that going from ρ = λ/10 to ρ = λ/100 lead to an increase in the number of iterations to achieve similar suboptimality for the classical coordinate descent, whereas the impact is more limited on the coordinate descent with Anderson extrapolation.

Finally, for a nonquadratic data-fit, here sparse logistic regression, we still demonstrate the applicability of extrapolated coordinate descent.

for multiple datasets and values of λ. We parametrize λ as a fraction of λ max = A y ∞ /2. As for the Lasso and the elastic net, the smaller the value of λ, the harder the problem and Anderson CD outperforms its competitors.

Conclusion

In this work, we have proposed to accelerate coordinate descent using Anderson extrapolation. We have exploited the fixed point itera-tions followed by coordinate descent iterates on multiple Machine Learning problems to improve their convergence speed. We have circumvented the nonsymmetricity of the iteration matrices by proposing a pseudo-symmetric version for which accelerated convergence rates have been derived. In practice, we have performed an extensive validation to demonstrate large benefits on multiple datasets and problems of interests. For future works, the excellent performance of Anderson extrapolation for cyclic coordinate descent calls for a more refined analysis of the known bounds, through a better analysis of the spectrum and numerical range of the iteration matrices.

A Additional experiments

In this section, we include the counterparts of Figures 7 to 9, but display the duality gap instead of the suboptimality. Indeed, since x * in not available in practice, the suboptimality cannot be used as a stopping criterion.

To create a dual feasible point, we use the classical technique of residual rescaling [START_REF]Sparse coding for machine learning, image processing and computer vision[END_REF].

A

A.5 Group Lasso

In this section we consider the group Lasso, with a design matrix A ∈ R n×p , a target y ∈ R n , and a partition G of [p] (elements of the partition being the disjoints groups):

arg min x∈R p 1 2 y -Ax 2 + λ g∈G x g , (22)
where for g ∈ G, x g ∈ R |g| is the subvector of x composed of coordinates in g. the group Lasso can be solved via Proximal Gradient Descent and by Block Coordinate Descent (BCD), the latter being amenable to Anderson Acceleration. As Figure 14 shows, the superiority of Anderson accelerated block coordinate descent is on par with the one observed on the problems studied above. Proof. Since x (i) = T x (i-1) + (x * -T x * ), c i (x (i)x (i-1) ) = c i (T x (i-1) + x * -T x *x (i-1) ) = (T -Id)c i (x (i-1)x * ) .

PGD

Hence, since

k 1 c i = 1, k i=1 c i (x (i) -x (i-1) ) = (T -Id) k i=1 c i x (i-1) -x * . ( 25
)
Lemma 6. For all c ∈ R k such that

Proof. In this proof, we denote by c * the solution of (2). We use the fact that for all c ∈ R k such that

c i T i-1 (x (0) -x * ) = k i=1 c i T i-1 (T -Id)(x (0) -x * ) = k i=1 c i T i-1 × (T -Id)(x (0) -x * ) ≤ H -1/2 k i=1 c i S i-1 H 1/2 × (T -Id)(x (0) -x * ) ≤ κ(H) k i=1 c i S i-1 × (T -Id)(x (0) -x * ) . (28)
Proof. We apply Lemma 6 by choosing c i equal to the Chebyshev weights c Cb i . Using the proof of Barré et al. (2020, Prop. B. 2), we have, with

ζ = 1- √ 1-ρ(T ) 1+ √ 1-ρ(T ) : k i=1 c Cb i S i-1 ≤ 2ζ k-1 1+ζ 2(k-1) . ( 29
)
Combined with Lemma 6 this concludes the proof:

(T -Id)(x ex * ) ≤ κ(H)

k i=1 c i S i-1 (T -Id)(x (0) -x * ) (30)
≤ κ(H) 2ζ k-1 1+ζ 2(k-1) (T -Id)(x (0)x * ) .

(31)

B.2 Proof of Proposition 4

Since g j are C 2 then prox gj are C 1 , see Gribonval and Nikolova (2020, Cor. 1.b). Moreover, f is C 2 and following [START_REF] Massias | Dual extrapolation for sparse generalized linear models[END_REF]; [START_REF] Klopfenstein | Model identification and local linear convergence of coordinate descent[END_REF] we have that:

ψ j : R p → R p x → prox gj            
x 1 . . .

x j-1 prox λgj /Lj x jγ j ∇ j f (x)

x j+1 . . .

x p             , (32)
is differentiable. Thus we have that the fixed point operator of coordinate descent:

ψ = ψ p •• • ••ψ 1 is differentiable.
Proposition 4 follows from the Taylor expansion of ψ in x * .

Figure 3 :

 3 Figure 3: OLS, rcv1. Suboptimality as a function of the number of gradient calls on the 5000 first columns of the dataset rcv1.

Figure 4 :

 4 Figure 4: 2 -regularised logistic regression, realsim. Suboptimality as a function of the number of gradient calls on the 2000 first features of the real-sim dataset, Tikhonov strength set so that κ = 10 5 .

Figure 5 :

 5 Figure 5: Influence of λ reg , sparse logitic regression, rcv1. Influence of the regularization amount when solving a sparse logistic regression using Anderson extrapolation with proximal coordinate descent (PCD) on the rcv1 dataset, K = 5, λ = λ max /30.

Figure 6 :

 6 Figure 6: Influence of K, quadratic, rcv1. Influence of the number of iterates K used to perform Anderson extrapolation with coordinate descent (CD) on a quadratic with the rcv1 dataset (1000 first columns).

FigureFigure 7 :

 7 Figure 7: Lasso, suboptimality. Suboptimality f (x (k) )f (x * ) as a function of the number of iterations for the Lasso on multiple datasets and values of λ.

Figure 8 :Figure 9 :

 89 Figure 8: Enet, suboptimality. Suboptimality f (x (k) )f (x * ) as a function of the number of iterations for the elastic net on Leukemia dataset, for multiple values of λ and ρ.

Figure 10 :Figure 11 :Figure 12 :Figure 13 :

 10111213 Figure 10: Influence of λ reg , quadratic, rcv1. Influence of the regularization amount when solving the linear system for Anderson extrapolation with coordinate descent (CD) on a quadratic with the rcv1 dataset (1000 first columns), K = 5.

Figure 14 :Lemma 5 .

 145 Figure14: Group Lasso, suboptimality. Suboptimality as a function of the number of iterations for the group Lasso on the Leukemia dataset, λ = λ max /100. Groups are artificially taken as consecutive blocks of 5 features.

  i (T -Id)(x (0)x * ) .

  i (x (i)x (i-1) ) = minc∈R k i ci=1 k i=1 c i (x (i)x (i-1) ) ≤ k i=1 c i (x (i)x (i-1) ) . (27)Then we use twice Lemma 5 for the left-hand and right-hand side of Equation (27). Using Lemma 5 with the c * i minimizing Equation (2) we have for all c i ∈ R such that k i=1 c i = 1 :(T -Id)(x ex * ) =

Table 1 :

 1 Datasets characteristics

	name	n	p	density
	gina agnostic	3468	970	1
	hiva agnostic	4229	1617	1
	leukemia	72	7129	1
	rcv1_train	20 242	19 960	3.7 10 -3
	real-sim	72 309	20 958	2.4 10 -3
	news20	19 996 632 983 6.1 10 -4