Kanemichi Muraoka

Miharu Eguchi

Osamu Ishitani
email: ishitani@chem.titech.ac.jp

François Cheviré

Kazuhiko Maeda
email: maedak@chem.titech.ac.jp

Selective CO2 reduction into

Keywords: Artificial photosynthesis, Heterogeneous photocatalysis, Mixed anion compounds, Solar fuels

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Catalytic conversion of anthropogenic CO 2 into energy rich chemicals with the aid of visible light energy could be a potential means of addressing the shortage of carbon resources and global warming issues [1][2][3][4][5][6][7][8]. Hybrid materials that consist of visible-light-absorbing semiconductors and (photo)catalytically active metal complexes have been reported to function as good photocatalysts for CO 2 reduction under visible light [2,7,8]. For example, a binuclear Ru(II) complex having a photosensitizer and catalyst units in one molecule (RuRu', see Fig. 1) works as an efficient homogeneous photocatalyst for CO 2 reduction, but requires a relatively strong electron donor for operation because of its low photooxidation ability. Combining the RuRu' with a visible-light-absorbing semiconductor that has strong oxidation ability enables one to accomplish CO 2 reduction even with the use of a weak electron donor [7]. In some photoelectrochemical systems constructed from a similar binuclear metal complex and a semiconductor, sacrificial reagent-free CO 2 reduction has been achieved [9][10][11]. As these hybrid systems are based on the combination of photocatalysis of semiconductors and molecules, it is important to develop both of the two in harmony. Focusing on the semiconductor component, various semiconductor materials including TaON [12], Ta 3 N 5 [13], GaN:ZnO [14], and Li 2 LaTa 2 O 6 N [15], C/N polymers [16] and Pb 2 Ti 2 O 5.4 F 1.2 [17] have been reported to become active semiconductor components for the Z-scheme CO 2 reduction with the aid of RuRu'. It has been also shown that for efficient CO 2 reduction, the potential of the conduction band minimum (CBM) of a semiconductor applied to the Z-scheme system has to be more negative than -1.3 V (vs Ag/AgNO 3 ), which is equivalent to the excited-state oxidation potential of the photosensitizer unit in RuRu', because suppression of back electron transfer from RuRu' to the semiconductor is essential [7]. Although several non-oxide materials have been reported to work as the semiconductor component in the hybrid photocatalyst with RuRu' for visible-light-driven CO 2 reduction, it is still important to search for a new semiconductor material for the application.

In heterogeneous photocatalysis for water splitting, rare earth-based semiconductor materials are known to show unique photocatalytic activities, depending on the rare-earth element. Ln-containing oxide semiconductors have been reported to show photocatalytic activities for water splitting under UV irradiation [18][19][20][21][22][23][24][25]. For visible light photocatalysis, Domen et al. have reported that oxysulfides Ln 2 Ti 2 S 2 O 5 showed activities for H 2 and O evolution half reactions in the presence of an electron donor and an acceptor, respectively [26]. Some Ln-containing transition metal oxynitrides have been shown to act as photocatalysts for the water reduction and/or oxidation half reactions under visible light [27][28][29]. Visible-light-driven overall water splitting has been accomplished using LaMg x Ta 1- [30]. Thus, there are a lot of examples of Ln-based water splitting photocatalysts. Recently, rare-earth species have been reported as effective modifiers for CO reduction over Ag/Ga 2 O 3 photocatalyst under UV irradiation [31]. However, there are very few that describe visible-light-driven Z-scheme CO 2 reduction using a Ln-containing semiconductor photocatalyst.

x O 1+3x N 2-3x (x ≥ 1/3)
We have previously reported Y-Ta oxynitride with a defect fluorite structure served as an effective semiconductor component in Z-scheme CO 2 reduction to produced formate, with confirmation by isotope tracer experiment using 13 CO 2 [32]. In this work, we examined Ln-Ta oxynitrides (LnTaO x N y ; Ln = Nd, Sm, Gd, Tb, Dy and Ho) as the new candidates for the Z-scheme CO 2 reduction system. The photocatalytic activities were discussed on the basis of the results of physicochemical analyses.

Experimental

Materials and reagents

Nd(NO 3 ) Gases of NH 3 (>99.9995%) and CO 2 (>99.9995%) were purchased from Sumitomo Seika Chemicals Co. and Taiyo Nippon Sanso Co., respectively.

MeCN was distilled over P 2 O 5 twice, and then distilled over CaH 2 prior to use. TEOA was distilled under reduced pressure (<1 Torr) and then stored under N 2 prior to use.

Synthesis of Ln-Ta oxynitrides LnTaO x N y

Ln-Ta oxides were synthesized by the polymerized complex (PC) method reported by Kakihana [33]. 0.014 mol of TaCl 5 was dissolved in 100 mL of methanol. Subsequently, 0.21 mol of CA was added to the methanol solution with continuous stirring at room temperature to produce a metal-CA complex. A 0.014 mol quantity of Ln(NO 3 ) 3 hydrate was then added to the solution and the mixture was magnetically stirred to make a transparent solution, after which 0.84 mol of EG was added. The resulting mixture was gradually heated to approximately 473 K to remove the methanol and to accelerate the esterification reaction between CA and EG. Upon continued heating, the solution began to gel, eventually transitioning to a glassy, transparent, brownish resin. This resin was thermally decomposed in a heating mantle at 673 K to produce a grey solid, and carbon-based compounds were removed by calcination on an Al 2 O 3 plate at 873 K for 2 h to yield a powder. This amorphous Ln-Ta oxide powder was nitrided by heating at 1123 K under a flow of gaseous NH 3 in a tubular furnace [34]. After nitridation, the resulting powder was washed six times with distilled water and dried at 343 K overnight.

Modification with Ag nanoparticles

The as-prepared LnTaO x N y powder was modified with Ag nanoparticles by an impregnation method. 10 mL of an aqueous AgNO 3 solution (containing 1.5 wt% Ag vs. LnTaO x N y ) was slowly added dropwise into 10 mL of H 2 O, in which 100 mg of LnTaO x N y was suspended, with continuous stirring. H 2 O was subsequently removed from the mixture by distillation under reduced pressure in a rotary evaporator. Finally, the as-obtained powder was placed in an alumina boat and heated under a flow of H 2 (20 mL min -1 ) at 473 K for 1 h in a tube furnace.

Our previous studies have shown that the loaded Ag nanoparticles on a visible-lightabsorbing semiconductor, which are in the form of nanoparticles of 3-30 nm in size, work as electron traps to facilitate interfacial electron transfer from the semiconductor to the photosensitizing unit of RuRu' during Z-scheme CO 2 reduction [16,32,35]. The detailed results and discussion can be found elsewhere.

Adsorption of RuRu' onto Ag-loaded LnTaO x N y

The as-prepared Ag/LnTaO x N y powder (30 mg) was dispersed in 15 mL of MeCN containing RuRuʹ (3.0 µmol g -1 vs. Ag/LnTaO x N y ). RuRu' was synthesized according to our previous report [12]. The suspension was stirred at room temperature in the dark overnight to allow for adsorption/desorption equilibration. The resulting powders were collected by filtration and washed with 5 mL of MeCN or methanol. UV-vis absorption measurements revealed that RuRu' quantitatively adsorbed on the surface of Ag/LnTaO x N y . After washing, the powders were dried at room temperature for 3 h under reduced pressure. In our previous studies, this procedure allows one to couple RuRu' on the surface of a semiconductor while keeping the structure of RuRu' [12,16,17]. It has been confirmed that both Ag and RuRu' were necessary to obtaining appreciable formate production [13,16,32,35]. The as-obtained hybrid material was employed as the photocatalyst for the CO 2 reduction reaction.

Characterization of materials

The prepared materials were characterized using X-ray diffraction (XRD; MiniFlex600, Rigaku, Cu K α radiation), UV-vis diffuse reflectance spectroscopy (DRS; V-565, Jasco), and transmission electron microscopy (TEM; H-7650, Hitachi High-Technologies). The Brunauer-Emmett-Teller (BET) surface area was measured using a nitrogen adsorption instrument (BELSOEP-mini II, MicrotracBEL) at liquid-nitrogen temperature.

Mott-Schottky plot measurements

LnTaO x N y electrodes were prepared by electrophoretic deposition on a fluorine-doped tin oxide (FTO) substrate. A pair of FTO glasses were immersed in 50 mL of acetone containing 100 mg of LnTaO x N y powder and 10 mg of I 2 . The FTO electrodes were connected to a direct current power supply (GW Instek PSW 80-13.5), to apply a bias of 50 V for 1 min. The LnTaO x N y -coated FTO electrode was finally heated at 673 K for 1 h under a N 2 flow in tube furnace.

Mott-Schottky plots were measured in MeCN/TEOA (4/1 v/v) mixed solution containing a 0.1 M Et 4 NBF 4 as an electrolyte at room temperature with a potentiostat (BAS ALS/CHI-760e) in a three-electrode system. LnTaO x N y /FTO electrodes were used as working electrodes. A Pt wire and a Ag/AgNO 3 electrode were used as the counter and reference electrodes, respectively. Mott-Schottky plots were recorded at a frequency of Hz.

Photocatalytic CO 2 reduction reaction

8 mg of photocatalyst powder was dispersed in 4 mL of a MeCN/TEOA mixed solution (4:1 v/v) in an 8 mL test tube. The suspension was subsequently purged with CO 2 gas for min prior to irradiation. The light source employed in this work was a 400 W high pressure Hg lamp (SEN) with an aqueous NaNO 2 solution filter to allow for visible light irradiation (λ > 400 nm). After 15 h of visible light irradiation, the liquid phase of the reaction mixture was analyzed using a capillary electrophoresis system (Otsuka Electronics Co. CAPI-3300) and the gaseous phase of the reaction cell was analyzed with a gas chromatograph (GL science, GC323) employing a thermal conductivity detector (TCD), an activated carbon column and argon as the carrier gas. The selectivity toward formate production during the CO 2 reduction reaction was calculated on the basis of the ratio of formate generated to the total amount of reduction products (i.e., formate, CO, and H 2 ). Blank tests were performed using the same photocatalyst suspension but without irradiation. Note that the amount of formate generation reported in this work was calculated by subtracting the blank portion.

The turnover number (TON) for formate production was calculated as:

TON = [Formate]/[RuRu'],
where [Formate] and [RuRu'] indicate the amounts of formate produced and RuRu' introduced in the hybrid photocatalyst, respectively.

Results and discussion

Synthesis and characterization of LnTaO x N y

It is known that the Ln-Ta-O-N system contains several phases such as perovskites, pyrochlores and defect fluorites, depending on the Ln element [34,36]. It has been suggested that LnTaO x N y with smaller Ln 3+ cations prefer to possess defect fluorite structure [34]. The obtained phase is also dependent on the synthesis route. As a crystalline phase, pyrochlore structure has been observed for Ln = Ce → Gd, but the use of soft chemistry route for the precursor synthesis enables to stabilize the defect fluorite type for a wide range of Ln [37]. Fig. 2 shows the XRD patterns of the as-synthesized LnTaO x N y materials. According to the previous literatures, we deduced that the as-synthesized LnTaO x N y had defect fluorite structure as the main phase, except for the Nd specimen that contained several different phases (see Fig. S1). Considering the larger ionic size of Nd 3+ , the Nd specimen may be pyrochlore (or a mixture with fluorite). As the atomic number increased, the positions of the diffraction peaks successively shifted to higher 2θ angles. This shift is reasonable considering that the ionic radii of Ln 3+ become smaller with increasing the atomic number due to the lanthanide contraction. A similar peak shift with respect to Ln 3+ ions has been observed in Ln 2 Ti 2 S 2 O 5 [26]. On the other hand, the full width at half maximum (FWHM) of the diffraction peaks were irregularly changed with respect to Ln 3+ in LnTaO x N y . As listed in Table 1, the FWHM of the 110 reflection is relatively small for Ln = Nd, Tb and Ho compared to others, suggesting higher crystallinity of LnTaO x N y (Ln = Nd, Tb and Ho). Fig. 3 shows TEM images of the same LnTaO x N y materials. All of them consisted of aggregated particles with primary particle sizes that range from several tens of nm to 50 nm, forming larger secondary particles. The specific surface areas of the synthesized LnTaO x N y were determined to be 14-15 m 2 g -1 (Table 1), which is very close to the value of an Y-Ta oxynitride reported previously (16 m 2 g -1 ) [32]. Therefore, no significant impact of the lanthanide element on the morphology and the specific surface area was observed.

Optical absorption properties and band-edge potentials

UV-vis DRS of LnTaO x N y are shown in Figs. 4 andS2, along with their oxide precursors for comparison. The absorption band edges of the precursor oxides appeared at around 300 nm, consistent with the report by Machida et al. on crystalized LnTaO 4 (Ln = Nd and Sm) [20]. Some sharp peaks appearing at wavelengths longer than 300 nm are attributed to f-f transitions derived from Ln 3+ ions [19][20][21].

Nitridation of these oxides with dry NH 3 at 1123 K to yield oxynitrides resulted in a significant red-shift of absorption edge caused by the reduction of the band gap, which arises from the N-2p state at the upper part of the valence band in LnTaO x N y , as observed in many oxynitride materials [38]. As the atomic number of Ln elements increased, the absorption band edge of LnTaO x N y tended to shift to longer wavelengths (Fig. S2), except for the Nd derivative that contained some impurity phases (Fig. S1). An identical trend has been observed in Ln 2 Ti 2 S 2 O 5 , in which the O-2p/S-3p/Ln-4f and S-3p/Ln-4f hybridized orbitals in the valence and conduction bands are claimed to cause the red-shift of absorption edge [26].

It should be also noted that LnTaO x N y (Ln = Nd, Ho) had characteristic absorption peaks derived from the f-f transitions, indicating that the Ln-4f states existed in the band gap of LnTaO x N.

Flat-band potential (E FB ) is an important parameter of a semiconductor that influences metal-complex/semiconductor hybrid photocatalysts [14]. To determine the E FB values of the as-prepared LnTaO x N y , Mott-Schottky plot measurements were conducted. As shown in Fig. 5, all of the Mott-Schottky plots have positive slopes, indicating n-type semiconducting character of LnTaO x N y . This should be reasonable, considering the fact that the LnTaO x N y compounds having defect fluorite structure contain a considerable number of anionic defects; the formation of anionic vacancies is accompanied by the generation of electrons, which is the origin of n-type semiconductivity. The E FB values were obtained by extrapolation of the linear portion to the x-axis intercept, and were listed in Table 1. The result showed that the values were a little fluctuant with respect to the Ln element, but ranged from -1.5 to -1.3 V vs. Ag/AgNO 3 . The conduction band minimum (CBM) of an n-type semiconductor lies negatively ca. 0.1-0.3 V from the E FB , depending on the conductivity [39]. Therefore, the CBM values of the LnTaO x N y could be estimated to be ca. -1.6±0.2 V, which is similar to or slightly more negative than that reported in Y-Ta oxynitride with the same defect fluorite structure (-1.40 V) [32] and is also negative enough to allow for the forward electron transfer from the semiconductor to the excited state of RuRu', while suppressing the back electron transfer from the excited RuRu' [7].

Potential / V vs. Ag/AgNO 3 Tb

Potential / V vs. Ag/AgNO 3 Dy 3.0

Potential / V vs. Ag/AgNO 3 Ho Fig. 5. Mott-Schottky plots for the as-synthesized LnTaO x N y .

Photocatalytic activity

The as-synthesized LnTaO x N y materials were tested as the semiconductor components of the hybrid CO 2 photoreduction system with RuRu'. CO 2 reduction reaction was conducted in a mixed solution of MeCN and TEOA under visible-light irradiation (λ >400 nm), in which TEOA acted as an electron donor to scavenge holes in the valence band of LnTaO x N y . As shown in Fig. 6, photocatalytic activities for the CO 2 reduction reaction were strongly dependent on the Ln element. In the cases of Ln = Nd and Sm, no formate was detected. On the other hand, other materials (Ln = Gd, Tb, Dy, and Ho) catalytically produced formate, with TONs greater than 1. Among the LnTaO x N y materials tested, the Tb derivative showed the highest activity, which was almost the same level as Y-Ta oxynitride reported before [32]. These active LnTaO x N y materials gave very high selectivity to formate (>99%), with negligible production of other products (e.g., CO or H 2 ).

Factors affecting photocatalytic activity

Thus, it was shown that the photocatalytic activity for the Z-scheme CO 2 reduction using LnTaO x N y was strongly influenced by the Ln element in LnTaO x N y . To accomplish the reaction, electrons and holes photogenerated in LnTaO x N y have to be consumed by reduction and oxidation processes, respectively, while minimizing recombination. It has been reported that the electron-hole recombination in a semiconductor photocatalyst can be minimized by improving the crystallinity of the semiconductor [40,41]. In the present work, the importance of crystallinity could be seen; SmTaO x N y and GdTaO x N y having lower crystallinity (i.e., greater FWHM, see Table 1) showed relatively low (or negligible) photocatalytic activities (Fig. 6). The inactivity of the Nd derivative is most likely due to the low quality of the material, as it contained some unreacted Nd impurities (Fig. S1); therefore, electron-hole recombination would be significant.

In addition to crystallinity of semiconductor, it appears that there are some other factors affecting the photocatalytic activity. As listed in Table 1, LnTaO x N y (Ln = Tb, Dy and Ho) had similar crystallinity to each other, but showed different activities (Fig. 6). While TbTaO x N y and HoTaO x N y had almost the same crystallinity and E FB , the former showed three times higher activity than the latter. UV-vis DRS (Fig. 4) showed an obvious contribution of Ln-4f states in the band gap of the lower activity specimen (Ln = Ho), in contrast to TbTaO x N y that showed the highest photocatalytic activity. The existence of Ln-4f states in the band gap of a semiconductor may lead to lower photocatalytic activity because such states can work as effective traps of photogenerated charge carriers, as reported earlier [20].

In the case of NdTaO x N y , as compared to TbTaO x N y , it is likely that the E FB (in other words, CBM potential) is another important factor that determines the photocatalytic performance of LnTaO x N y , in addition to the negative impact of the existence of Ln-4f states in the band gap. In fact, the E FB of NdTaO x N y was 0.15 V more positive than that of TbTaO x N y , indicating that the reactivity of the photoexcited electron is lower in the former than in the latter. Although one may think that the difference of 0.1-0.2 V is small, it has been reported that such a small difference (0.1-0.2 V) in the CBM potential of semiconductors causes a significant difference in photocatalytic activity for the Z-scheme CO 2 reduction using a visible-light-absorbing semiconductor and RuRu' [14].

On the basis of the discussion above, high crystallinity, more negative conduction band potential and the absence of Ln-4f states in the band gap are important factors for LnTaO x N y semiconductors showing high photocatalytic activity. Nevertheless, the absolute activity of the present hybrid materials (even for the best-performing one with TbTaO x N y ) was still very low, as compared to similar hybrid materials constructed with other non-oxide semiconductors such as TaON or Ta 3 N 5 . The synthesized TbTaO x N y seems to work under longer wavelength visible light (e.g., >500 nm) owing to its extended absorption band (~600 nm), but the TbTaO x N y -based hybrid did not work under >500 nm irradiation. This means that improvement of the preparation method for LnTaO x N y is necessary, and will be investigated as part of our future works.

Conclusions

Hybrid photocatalysts were developed using visible-light-absorbing LnTaO x N y (Ln = Nd, Sm, Gd, Tb, Dy and Ho) semiconductors and a Ru(II) binuclear complex for CO 2 reduction under visible light (λ > 400 nm). Selective reduction of CO 2 into formate was achieved with hybrids consisting of LnTaO x N y (Ln = Gd, Tb, Dy and Ho) in MeCN containing TEOA as an electron donor. The highest activity was obtained with the Tb derivative. The use of LnTaO x N y with high crystallinity and more negative conduction band potential, while avoiding participation of Ln-4f states in the band gap structure, was found to be essential to attaining high photocatalytic activity. We believe that the knowledge obtained in this work will be useful for designing a new semiconductor material not only for the Zscheme CO 2 reduction with a supramolecular complex but also for a system combined with a molecular catalyst, because semiconductor photocatalysis is common in both systems.

Supporting Information. Additional characterization data. This material is available free of charge via the Internet. Ln-Ta oxynitrides with defect fluorite structure were found to serve as effective semiconductor components in visible-light-driven Z-scheme CO 2 reduction with a binuclear Ru(II) complex, giving high selectivity to formate (>99%).

Graphical abstract

Fig. 1 .

 1 Fig. 1. Hybrid photocatalyst constructed from a visible-light-absorbing semiconductor and a binuclear Ru(II) complex (RuRu').

Fig. 2 .

 2 Fig. 2. XRD patterns for the as-synthesized LnTaO x N y .

Fig. 3 .Fig. 4 .

 34 Fig. 3. TEM images for the as-synthesized LnTaO x N y .

Fig. 6 .

 6 Fig. 6. Photocatalytic activities for CO 2 reduction over RuRu'/Ag/LnTaO x N y under visible light (λ >400 nm). Reaction conditions: catalyst, 8.0 mg; reactant solution, mixed solution of MeCN and TEOA (4:1 v/v, 4.0 mL); light source, 400 W high-pressure Hg lamp with a K 2 CrO 4 aqueous solution filter. Reaction time: 15 h. The dotted red line indicates the level of TON of 1 with respect to RuRu'. No reaction took place in the absence of either LnTaO x N y or RuRu'.

  3

Table 1 .

 1 FWHM, specific surface areas of the as-synthesized LnTaO x N y .

	Ln FWHM of the 110 reflection Specific surface area (m 2 g -1 ) Flat-band potential (V vs. Ag/AgNO 3 )
	Nd 0.23	14	-1.44
	Sm 0.35	15	-1.38
	Gd 0.37	15	-1.44
	Tb 0.27	14	-1.47
	Dy 0.28	14	-1.32
	Ho 0.28	14	-1.48

ACKNOWLEDGMENTS

This work was supported by a Grant-in-Aid for Scientific Research on Innovative Area "Mixed Anion (Projects JP16K21724 and JP16H06441)" from the Japan Society for the Promotion of Science (JSPS) and a CREST program (Project JPMJCR13L1) from the Japan Science and Technology Agency (JST). K. Muraoka wishes to acknowledge support in the form of a JSPS Fellowship for Young Scientists (Project JP17J06914).


Tadié Jm
email: adel.maamar@chu-rennes.fr

Reply to: return of spontaneous circulation after an out-of-hospital cardiac arrest: an acute brain injury like others?

Keywords: Cardiac arrest, Ultrasonography Transcranial Doppler, prognosis, critically ill patients, cerebral monitoring

Word count: 370. Despite its non-invasive nature, Transcranial Doppler Sonography (TCD) remains underused at the return of spontaneous circulation (ROSC) phase after out-of-hospital cardiac arrest (OHCA). This is probably due to the low interrater reliability of some TCD measurements and the conflicting results in previous studies [START_REF] Cardim | A comparison of non-invasive versus invasive measures of intracranial pressure in hypoxic ischaemic brain injury after cardiac arrest[END_REF][START_REF] Doepp | Duplex sonography of cerebral blood flow after cardiac arrest-A prospective observational study[END_REF]. In our study, we found that patients with poor neurological outcomes had lower diastolic flow velocity (DFV) and higher pulsatility index (PI) at early stage after ROSC than patients with good neurological outcomes [START_REF] Rafi | Doppler sonography of cerebral blood flow for early prognostication after out-of-hospital cardiac arrest: DOTAC study[END_REF]. Therefore, we assumed that TCD measurements may allow to estimate cerebral blood flow (CBF), cerebral vascular resistance and therefore to optimize cerebral perfusion [START_REF] Edouard | Non-invasive assessment of cerebral perfusion pressure in brain injured patients with moderate intracranial hypertension[END_REF].

Nevertheless, TCD measurements provide only a surrogate of the cerebral blood flow (CBF) and does not allow to affirm if the cerebral autoregulation is preserved or not [START_REF] Sundgreen | Autoregulation of cerebral blood flow in patients resuscitated from cardiac arrest[END_REF].

As Engrand et al. have pointed out, the relationship between TCD and CBF, and the interpretation of an elevated PI are complex and not completely understood.

J o u r n a l P r e -p r o o f

Indeed, several processes would occur following a recovered cardiac arrest: elevation of the CBF, metabolic disorders, imbalance of the balance supply and need for tissue oxygen, reperfusion lesions, all contributing to cell death [START_REF] Sinha | Monitoring the Brain After Cardiac Arrest: a New Era[END_REF].

As mentioned in our study, the CBF alteration profile in post cardiac arrest has been studied in animals, and shows four phases [START_REF] Iordanova | Alterations in Cerebral Blood Flow after Resuscitation from Cardiac Arrest[END_REF]: absence of flow during cardiac arrest, hyperemia in the first 30 minutes following ROSC, hypoperfusion in the following 6 hours and then last fixed phase with either a return to normal, or persistence of prolonged hypoperfusion, or hyperemia. We hypothesize that TCD measurement should be focused on the early phase of ROSC and may allow the development of individualized targeted therapies to improve cerebral perfusion at this phase. Different targeted therapies have been tested in animals which would improve the neurological prognosis by reducing early hyperemia of the subcortical areas [START_REF] Shaik | 20-Hydroxyeicosatetraenoic Acid Inhibition by HET0016 Offers Neuroprotection, Decreases Edema, and Increases Cortical Cerebral Blood Flow in a Pediatric Asphyxial Cardiac Arrest Model in Rats[END_REF][START_REF] Leonov | Hypertension with hemodilution prevents multifocal cerebral hypoperfusion after cardiac arrest in dogs[END_REF].

We believe that TCD should complement other neuromonitoring tools. The tools for monitoring cerebral hemodynamics are widely used for traumatic brain injury [START_REF] Bouzat | Transcranial Doppler to screen on admission patients with mild to moderate traumatic brain injury[END_REF], but are poorly developed in the management of the patient in post-ROSC.

Finally, results of our study strongly encourage us to carry out further investigations to understand the different mechanisms involved at early stage of ROSC in CBF modifications, as recommended for others acute brain injury situations.

Conflicts of interest statement:

On behalf of all authors, the corresponding author states that there are no conflicts of interest.

J o u r n a l P r e -p r o o f


Francis Bach
email: francis.bach@inria.fr

On the Effectiveness of Richardson Extrapolation in Machine Learning

Richardson extrapolation is a classical technique from numerical analysis that can improve the approximation error of an estimation method by combining linearly several estimates obtained from different values of one of its hyperparameters, without the need to know in details the inner structure of the original estimation method. The main goal of this paper is to study when Richardson extrapolation can be used within machine learning, beyond the existing applications to step-size adaptations in stochastic gradient descent. We identify two situations where Richardson interpolation can be useful: (1) when the hyperparameter is the number of iterations of an existing iterative optimization algorithm, with applications to averaged gradient descent and Frank-Wolfe algorithms (where we obtain asymptotically rates of O(1/k 2 ) on polytopes, where k is the number of iterations), and (2) when it is a regularization parameter, with applications to Nesterov smoothing techniques for minimizing non-smooth functions (where we obtain asymptotically rates close to O(1/k 2 ) for non-smooth functions), and ridge regression. In all these cases, we show that extrapolation techniques come with no significant loss in performance, but with sometimes strong gains, and we provide theoretical justifications based on asymptotic developments for such gains, as well as empirical illustrations on classical problems from machine learning.

Introduction

Many machine learning methods can be cast as looking for approximations of some ideal quantity which cannot be readily computed from the data at hand: this ideal quantity can be the predictor learned from infinite data, or an iterative algorithm run for infinitely many iterations. Taking theirs roots in optimization and more generally numerical analysis, many accelerations techniques have been developed to tighten these approximations with as few changes as possible to the original method.

While some acceleration techniques add some simple modifications to a known algorithm, such as Nesterov acceleration for the gradient descent method [START_REF] Yurii | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF], extrapolation techniques are techniques that do not need to know the fine inner structure of the method to be accelerated. These methods are only based on the observations of solutions of the original method.

Extrapolation techniques work on the vector-valued output x t ∈ R d of the original method that depends on some controllable real-valued quantity t, which can be the number of iterations or some regularization parameter, and more generally any parameter that controls both the running time and the approximation error of the algorithm. When t tends to t ∞ (which is typically 0 or +∞),

x t 2x t -x 2t x 2t x * x * + t∆ x * + 2t∆
Figure 1: Illustration of Richardson extrapolation for t ∞ = 0 and x t = x * + t∆ + O(t 2 ). Iterates (in black) with their first-order expansions (in red). The deviations (represented by circles) are of order O(t 2 ). Adapted from [START_REF] Dieuleveut | Bridging the gap between constant step size stochastic gradient descent and Markov chains[END_REF].

we will assume that x t has an asymptotic expansion of the form

x t = x * + g t + O(h t ),
where x * is the desired output, g t ∈ R d is the asymptotic equivalent of x t -x * , and h t = o( g t ).

The key question in extrapolation is the following: from the knowledge of x t for potentially several t's, how can we better approximate x * , without the full knowledge of g t ?

For exponentially converging algorithms, there exist several "non-linear" schemes that combine linearly several values of x t with weights that depend non-linearly on the iterates, such as Aitken's ∆ 2 process [START_REF] Craig | On Bernoulli's numerical solution of algebraic equations[END_REF] or Anderson acceleration [START_REF] Donald G Anderson | Iterative procedures for nonlinear integral equations[END_REF], which has recently been shown to provide significant acceleration to linearly convergent gradient-based algorithms [START_REF] Scieur | Regularized nonlinear acceleration[END_REF]. In this paper, we consider dependence in powers of t, where Richardson extrapolation excels (see, e.g., [START_REF] Fry | The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam[END_REF][START_REF] Joyce | Survey of extrapolation processes in numerical analysis[END_REF][START_REF] Gautschi | Numerical Analysis[END_REF].

We thus assume that g t = t α • ∆, and h t = t β is a power of t such that h t = o( g t ), where α ∈ R is known but ∆ ∈ R d is unknown, that is

x t = x * + t α ∆ + O(t β ).
In all our cases, α = -1 when t ∞ = +∞ and α = 1 when t = 0. Richardson extrapolation is simply combining two iterates with different values of t so that the zero-th order term x * is preserved, while the first-order term cancels, for example:

2x t -x 2 1/α t = 2(x * + t α ∆ + O(t β )) -(x * + 2t α ∆ + O(t β )) = x * + O(t β ).
See an illustration in Figure 1 for α = 1, β = 2, and t ∞ = 0. Note that: (a) the choice of 2 1/α = 1 as a multiplicative factor is arbitrary and chosen for its simplicity when |α| = 1, and (b) Richardson extrapolation can be used with m + 1 iterates to remove the first m terms in a asymptotic expansion, where the powers of the expansion are known and not the associated vector-valued constants (see examples in Section 3).

The main goal of this paper is to study when Richardson extrapolation can be used within machine learning, beyond the existing explicit applications to step-size adaptations in stochastic gradient descent [START_REF] Durmus | Stochastic gradient Richardson-Romberg Markov chain Monte Carlo[END_REF][START_REF] Dieuleveut | Bridging the gap between constant step size stochastic gradient descent and Markov chains[END_REF], and implicit wide-spread use in "tail-averaging" (see more details in Section 2.1).

We identify two situations where Richardson interpolation can be useful:

x 0

x 1

x 2

x 0

x 1

x 2

∆

Figure 2: Left: Oscillating convergence, where Richardson extrapolation does not lead to any gain. Right: non-oscillating convergence, with a main direction ∆ (in red dotted), where Richardson extrapolation can be beneficial if the oscillations orthogonal to the direction ∆ are negligible compared to convergence along the direction ∆.

• t = k is the number of iterations of an existing iterative optimization algorithm converging to x * , where then α = -1 and t ∞ = +∞, and Richardson extrapolation considers, for k even, x

k = 2x k -x k/2 . We consider in Section 2, averaged gradient descent and Frank-Wolfe algorithms (where we obtain asymptotically rates of O(1/k 2 ) on polytopes, where k is the number of iterations).

• t = λ is a regularization parameter, where then α = 1 and t ∞ = 0, and Richardson extrapolation considers x

(1) λ = 2x λ -x 2λ . We consider in Section 3, Nesterov smoothing techniques for minimizing non-smooth functions (where we obtain asymptotically rates close to O(1/k 2 ) for non-smooth functions), and ridge regression (where we obtain estimators with lower bias).

As we will show, extrapolation techniques come with no significant loss in performance, but with sometimes strong gains, and the goal of this paper is to provide theoretical justifications for such gains, as well as empirical illustrations on classical problems from machine learning. Note that we aim for the simplest asymptotic results (most can be made non-asymptotic with extra assumptions).

Extrapolation on the Number of Iterations

In this section, we consider extrapolation based on the number of iterations k, that is, for the simplest case x

(1)

k = 2x k -x k/2
, for optimization algorithms aimed at minimizing on R d If x k is converging to a minimizer x * , then so is x k/2 , and thus also x

(1) k = 2x k -x k/2 ; moreover, we have x (1) k -x * 2 2 x k -x * 2 + x k/2 -x * 2
, so even if there are no cancellations, performance is never significantly deteriorated (the risk is essentially to lose half of the iterations).

The potential gains depend on the way x k converges to x * . The existence of a convergence rate of the form f

(x k ) -f (x * ) = O(1/k) or O(1/k 2 )
is not enough, as Richardson extrapolation requires a specific direction of asymptotic convergence. As illustrated in Figure 2, some algorithms are oscillating around their solutions, while some converge with a specific direction. Only the latter ones can be accelerated with Richardson extrapolation, while the former ones are good candidates for Anderson acceleration [START_REF] Donald G Anderson | Iterative procedures for nonlinear integral equations[END_REF][START_REF] Scieur | Regularized nonlinear acceleration[END_REF].

We now consider three algorithms: (1) averaged gradient descent, where extrapolation is at its best, as it transforms an O(1/t 2 ) convergence rate into an exponential one, (2) accelerated gradient descent, where extrapolation does not bring anything, and (3) Frank-Wolfe algorithms, where the situation is mixed (sometimes it helps, sometimes it does not).

Averaged gradient descent

We consider the usual gradient descent algorithm

x k = x k-1 -γf (x k-1 ),
where γ 0 is a step-size, with Polyak-Ruppert averaging [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF][START_REF] Ruppert | Efficient estimations from a slowly convergent Robbins-Monro process[END_REF]:

xk = 1 k k-1 i=0 x i .
Averaging is key to robustness to potential noise of the gradients [START_REF] Polyak | Acceleration of stochastic approximation by averaging[END_REF][START_REF] Nemirovski | Robust stochastic approximation approach to stochastic programming[END_REF]. However it comes with the unintended consequence of losing the exponential forgetting of initial conditions for strongly convex problems [START_REF] Bach | Non-asymptotic analysis of stochastic approximation algorithms for machine learning[END_REF].

A common way to restore exponential convergence (up to the noise level in the stochastic case) is to consider "tail-averaging", that is, to replace xk by the average of only the latest k/2 iterates [START_REF] Jain | Parallelizing stochastic gradient descent for least squares regression: Mini-batching, averaging, and model misspecification[END_REF]. As shown below for k even, this corresponds exactly to Richardson extrapolation:

2 k k-1 i=k/2 x i = 2 k k-1 i=0 x i - 2 k k/2-1 i=0 x i = 2x k -xk/2 .
While [START_REF] Jain | Parallelizing stochastic gradient descent for least squares regression: Mini-batching, averaging, and model misspecification[END_REF] focuses on a non-asymptotic analysis for stochastic problems for leastsquares regression, we now provide an asymptotic analysis for general convex objective functions and non-stochastic problems (see proof in Appendix B).

Proposition 1 Assume f convex, three-times differentiable with Hessian eigenvalues bounded by L, bounded third-order derivatives, and a unique minimizer

x * ∈ R d such that f (x * ) is positive definite. If γ 1/L, then xk = x * + 1 k ∆ + O(exp(-kλ)),
where ∆ = ∞ i=0 (x i -x * ) and λ is proportional to γλ min (f (x * )).

Note that: (a) before Richardson extrapolation, the asymptotic convergence rate will be of the order O(1/k 2 ), which is better than the usual O(1/k) upper-bound for the rate of gradient descent, but with a stronger assumption that in fact leads to exponential convergence before averaging, (b) while ∆ has a simple expression, it cannot be computed in practice, (c) that Richardson extrapolation leads to an exponentially convergent algorithm from an algorithm converging asymptotically in O(1/k 2 ) for functions values, and (d) that in the presence of noise in the gradients, the exponential convergence would only be up to the noise level. See Figure 3 (left plot) for an illustration with noisy gradients.

Accelerated gradient descent

In the section above, we considered averaged gradient descent, which is asymptotically converging as O(1/k 2 ) and on which Richardson extrapolation could be used with strong gains. Is it possible also for the accelerated gradient descent [START_REF] Yurii | A method for solving the convex programming problem with convergence rate O(1/k 2 )[END_REF], which has a (non-asymptotic) convergence rate of O(1/k 2 ) for convex functions?

It turns out that the behavior of the iterates of accelerated gradient descent is exactly of the form depicted in the left plot of Figure 2: that is, the iterates x k oscillate around the optimum, as can be seen from the spectral analysis for quadratic problems, in continuous time [START_REF] Su | A differential equation for modeling Nesterov's accelerated gradient method: theory and insights[END_REF] or discrete time [START_REF] Flammarion | From averaging to acceleration, there is only a step-size[END_REF]. Richardson extrapolation is of no help, but is not degrading performance too much. See Figure 3 (right plot) for an illustration.

f (x) = 1 n n i=1 log(1 + exp(-b i x a i )), with (a i , b i ) ∈ R d × {-1, 1}.
The covariance matrix of (Gaussian) inputs has eigenvalues 1/j, j = 1, . . . , d; the lowest eigenvalue is 1/400 and therefore we can see the effect of strong convexity starting between k = 100 and 1000 iterations; moreover, for the regular averaged recursion, the line in the log-log plot has slope -2. Right: accelerated gradient descent on a quadratic optimization problem in dimension d = 1000 and a Hessian whose eigenvalues are 1/j 2 , j = 1, . . . , d; with such eigenvalues, the local linear convergence is not observed and we have a line of slope -2.

Frank-Wolfe algorithms

We now consider Frank-Wolfe algorithms (also known as conditional gradient algorithms) for minimizing a smooth convex function f on a compact convex set K. These algorithms are dedicated to situations where one can easily minimize linear functions on K (see, e.g., Jaggi, 2013, and references therein). The algorithm has the following form:

xk ∈ arg min x∈K f (x k-1 ) + f (x k-1 ) (x -x k-1 ) x k = (1 -ρ k )x k-1 + ρ k xk .
That is, the first order Taylor expansion of f at x k-1 is minimized, ending up typically in an extreme point xk or K and a convex combination of x k-1 and xk is considered. While some form of line search can be used to find ρ k , we consider so-called "open loop" schemes where ρ k = 1/k or ρ k = 2/(k + 1) [START_REF] Dunn | Conditional gradient algorithms with open loop step size rules[END_REF][START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF].

In terms of function values, these two variants are known to converge at respective rates O(log(k)/k) and O(1/k). Moreover, as illustrated in Figure 4, they are known to zig-zag towards the optimal point. Avoiding this phenomenon can be done in several ways, for examples through optimizing over all convex combinations of the xi 's for i k [START_REF] Balder | Simplicial decomposition in nonlinear programming algorithms[END_REF], or through so-called "away steps" [START_REF] Guélat | Some comments on Wolfe's away step[END_REF][START_REF] Lacoste-Julien | On the global linear convergence of Frank-Wolfe optimization variants[END_REF]. In this section, we consider Richardson extrapolation and assume for simplicity that K is a polytope (which is a typical use case for Frank-Wolfe algorithms). Note here that we are considering asymptotic convergence rates, and even without extrapolation (but with a local strong-convexity assumption), we can beat the O(1/k) rates for the step-size ρ k = 2/(k + 1).1 Asymptotic expansion. In order to provide the proposition below (see Appendix C for a proof) that characterizes the zig-zagging phenomenon, we assume regularity properties similar to Section 2.1

x * K y 2 y 1 x 0 x 1 x 2 ∆ Figure 4
: Frank-Wolfe algorithm zigzagging. Starting from x 0 , the algorithm always moves towards one of the extreme points of K, with an average direction of ∆.

0 2 4 log 10 (k) -6 -4 -2 0 log 10 [ f(x)-f * ] step 1/k x k 2x k -x k/2 0 2 4 log 10 (k) -6 -4 -2 0 log 10 [ f(x)-f * ]
step 2/(k+1)

x k 2x k -x k/2
Figure 5: Frank-Wolfe for the "constrained logistic Lasso", that is, min

x∈R d f (x) such that x 1 c, with f (x) = 1 n n i=1 log(1+exp(-b i x a i ))
. We consider n = 400 observations in dimension d = 400, sampled from a standard normal distribution, and with a constraint on the 1 -norm. Left: step size 1/k with slopes -1 (blue) and -2 (red). Right: step size 2/(k + 1), with slope approximately -2 for the two curves.

0 2 4 log 10 (k) -8 -6 -4 -2 0 2 log 10 [ f(x)-f * ] step 1/k x k 2x k -x k/2 0 2 4 log 10 (k) -8 -6 -4 -2 0 2 log 10 [ f(x)-f * ]
step 2/(k+1)

x k 2x k -x k/2
Figure 6: Frank-Wolfe for the dual of robust regression. We consider the dual of absolute loss regression with n = 400 observations in dimension d = 200, sampled from a standard normal distribution, with a squared Euclidean norm penalty. The primal problem is inf

y∈R d 1 n b -Ax 1 + λ 2 y 2 2 while the dual problem is sup x ∞ 1 -f (x), with f (x) = -1 n x b + 1 2n 2 λ
x AA x. Left: step size 1/k with slopes -1 (blue) and -2 (red). Right: step size 2/(k + 1), with slope -2 for the two curves.

and that the unique minimizer is "in the middle of a face" of K, which is often referred to as constraint qualification in optimization [START_REF] Nocedal | Numerical Optimization[END_REF].

Proposition 2 Assume f convex, three-times differentiable with bounded third-order derivatives in the polytope K, and with a unique minimizer x * ∈ R d such that f (x * ) is positive definite. Moreover, we assume x * is strictly in a (m -1)-dimensional face of K which is the convex hull K * of extreme points y 1 , . . . , y m ∈ R d , and for which min y∈K f (x * ) y is attained only by elements of K * . Then:

• For ρ k = 1/k, x k = x * + 1 k ∆ 1 + O(1/k 2 ). This implies f (x k ) -f (x * ) = 1 k ∆ 1 f (x * ) + O(1/k 2 ) and f (2x k -x k/2 ) -f (x * ) = O(1/k 2 ). • For ρ k = 2/(k + 1), x k = x * + 1 k(k+1) ∆ 2 + O(1/k 2 ). This implies f (x k ) -f (x * ) = O(1/k 2 ) and f (2x k -x k/2 ) -f (x * ) = O(1/k 2 ).
The two vectors ∆ 1 and ∆ 2 , are orthogonal (for the dot-product defined by f (x * )) to the span of all y i -x * , i = 1, . . . , m

We now discuss the consequences of the proposition above.

Step-size ρ k = 1/k. Although it leads to a worse performance than the step-size ρ k = 2/(k + 1) both in theory (extra logarithmic term) and in practice, we consider it as (1) this is the "historical" step-size [START_REF] Dunn | Conditional gradient algorithms with open loop step size rules[END_REF] with an interesting behavior in our set-up, and (2) the corresponding dual algorithm is the subgradient method with plain averaging (see, e.g., [START_REF] Bach | Duality between subgradient and conditional gradient methods[END_REF], which is sometimes preferred in online learning [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF].

As shown in Prop. 2, Richardson extrapolation allows to go from an O(1/k) to an O(1/k 2 ) convergence rate. In the left plots of Figure 5 and Figure 6, we can observe the benefits of Richardson extrapolation on two optimization problems with the step-size ρ k = 1/k. Note that: (a) asymptotically, there is provably no extra logarithmic factor like we have for the non-asymptotic convergence rate, and (b) that the Richardson extrapolated iterate may not be within K, but is O(1/k) away from it (in our simulations, we simply make the iterate feasible by rescaling).

Step-size ρ k = 2/(k + 1). As shown in Prop. 2, we already get an performance of O(1/k 2 ) without extrapolation (which is a new asymptotic result for Frank-Wolfe algorithm on polytopes), and Richardson extrapolation does not lead to any acceleration. In the right plots of Figure 5 and Figure 6, we indeed see no benefits (but no strong degradation).

Note that here (and for both step-sizes), higher order Richardson would not lead to further cancellation as within the span of the supporting face, we have an oscillating behavior similar to the left plot of Figure 2. Moreover, although we do not have a proof, the closed loop algorithm exhibits the same behavior as the step ρ k = 1/k, both with and without extrapolation, which is consistent with the analysis of [START_REF] Michael | A tight upper bound on the rate of convergence of Frank-Wolfe algorithm[END_REF]. It would also be interesting to consider the benefits for Richardson extrapolation for strongly-convex sets [START_REF] Garber | Faster rates for the Frank-Wolfe method over strongly-convex sets[END_REF].

Extrapolation on the Regularization Parameter

In this section, we explore the application of Richardson extrapolation to regularization methods. In a nutshell, regularization allows to make an estimation problem more stable (less subject to variations for statistical problems) or the algorithm faster (for optimization problems). However, regularization adds a bias that needs to be removed. In this section, we apply Richardson extrapolation to the regularization parameter to reduce this bias. We consider two applications where we can provably show some benefits: (a) smoothing for non-smooth optimization in Section 3.1, and (b) ridge regression in Section 3.2.

Other applications include the classical use within integration methods, where the technique is often called Richardson-Romberg extrapolation [START_REF] Gautschi | Numerical Analysis[END_REF], and for bias removal in constantstep-size stochastic gradient descent for sampling [START_REF] Durmus | Stochastic gradient Richardson-Romberg Markov chain Monte Carlo[END_REF] and optimization [START_REF] Dieuleveut | Bridging the gap between constant step size stochastic gradient descent and Markov chains[END_REF].

Smoothing non-smooth problems

We consider the minimization of a convex function of the form f (x) = h(x)+g(x), where h is smooth and g is non-smooth. These optimization problems are ubiquitous in machine learning and signal processing, where the lack of smoothness can come from (a) non-smooth losses such as max-margin losses used in support vector machines and more generally structured output classification [START_REF] Taskar | Learning structured prediction models: A large margin approach[END_REF][START_REF] Tsochantaridis | Large margin methods for structured and interdependent output variables[END_REF], and (b) sparsity-inducing regularizers (see, e.g., Bach et al., 2012, and references therein). While many algorithms can be used to deal with this non-smoothness, we consider a classical smoothing technique below.

Nesterov smoothing. In this section, we consider the smoothing approach of [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF] where the non-smooth term is "smoothed" into g λ , where λ is a regularization parameter, and accelerated gradient descent is used to minimize h + g λ .

A typical way of smoothing the function g is to add λ times a strongly convex regularizer to the Fenchel conjugate of g (see an example below); as shown by [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF], this leads to a function g λ which has a smoothness constant (defined as the maximum of the largest eigenvalues of all Hessians) proportional to 1/λ, with a uniform error of λ between g and g λ . Given that accelerated gradient descent leads to an iterate with excess function values proportional to 1/(λk 2 ) after k iterations, with the choice of λ ∝ 1/k, this leads to an excess in function values proportional to 1/k, which improves on the subgradient method which converges in O(1/ √ k).

Richardson extrapolation. If we denote by x λ the minimizer of h + g λ and x * the global minimizer of f = h + g, if we can show that

x λ = x * + λ∆ + O(λ 2 ), then x (1) λ = 2x λ -x 2λ = O(λ 2 ) and we can expand f (x (1) λ ) = f (x * ) + O(λ 2 ), which is better than the O(λ) approximation without extrapolation.
Then, with λ ∝ k -2/3 , to balance the two terms 1/(λk 2 ) and λ 2 , we get an overall convergence rate for the non-smooth problem of k -4/3 . We now make this formal for the special (but still quite generic) case of polyhedral functions g, and also consider m-step Richardson extrapolation, which leads to a convergence rate arbitrarily close to O(1/k 2 ).

Polyhedral functions. We consider a polyhedral function of the form

g(x) = max i∈{1,...,m} a i x -b i = max i∈{1,...,m} (Ax -b) i ,
where A ∈ R m×d and b ∈ R m . This form includes traditional regularizers such as the 1 -norm, the ∞ -norm, grouped 1 -∞ -norms [START_REF] Negahban | Joint support recovery under high-dimensional scaling: Benefits and perils of 1,∞ -regularization[END_REF], or more general sparsityinducing norms [START_REF] Bach | Structured sparsity-inducing norms through submodular functions[END_REF].

We consider the smoothing of g as:

g λ (x) = max η∈∆m η (Ax -b) -λϕ(η),
for some strongly convex function ϕ, typically, the negative entropy

m i=1 η i log η i -η i , or 1 2 η 2 2 .
For our asymptotic expansion, we also need a form a constraint qualification (see proof in Appendix D).

Proposition 3 Assume h convex, three-times differentiable with bounded third-order derivatives g convex, and a unique minimizer x * ∈ R d of h + g such that h (x * ) is positive definite. Assume there exists η * ∈ R m in the simplex ∆ m (defined as the vectors in R m with non-negative components summing to one), such that, for the support I ⊂ {1, . . . , m} of η * (that is, the set of non-zeros), h (x * ) + A η * = 0, and max i∈{1,...,m} (Ax * -b) i is only attained for all i ∈ I.

Assume moreover the submatrix A I obtained by taking the the rows of A indexed by I has full rank. We denote by x λ a minimizer of h(x) + g λ (x), and η λ the corresponding dual variable. Then:

x λ = x * + λ∆ + O(λ 2 ), with ∆ = h (x * ) -1 A I A I h (x * ) -1 A I -1 (η * ) I
for the quadratic penalty, and a similar expression for the entropic penalty.

The proposition above implies that (see detailed proof in Appendix D for details) the smoothing technique asymptotically adds a bias of order λ

f (x λ ) = f (x * ) + O(λ),
where we recover the usual upper-bound in O(λ), confirming the result from [START_REF] Nesterov | Smooth minimization of non-smooth functions[END_REF]. The key other consequence is that

f (x (1) λ ) = f (2x λ -x 2λ ) = f (x * ) + O(λ 2 ),
which shows the benefits of Richardson extrapolation.

Multiple-step Richardson extrapolation. Given that one-step Richardson extrapolation allows to go from a bias of O(λ) to O(λ 2 ), a natural extension is to consider m-step Richardson extrapolation [START_REF] Gautschi | Numerical Analysis[END_REF], that is, a combination of m + 1 iterates:

x (m) λ = m+1 i=1 α (m) i x iλ ,
where the weights α (m) i are such that the first m powers in the Taylor expansion of x λ cancel out. This can be done by solving the linear system based on the following equations:

m+1 i=1 α (m) i = 1 (1) ∀j ∈ {1, . . . , m}, m+1 i=1 α (m) i i j = 0. ( 2
)
Using the same technique as Pagès (2007, Lemma 3.1), this is a Vandermonde system with a closed form solution (see proof in Appendix E.3):

α (m) i = (-1) i-1 m + 1 i .
We show in Appendix D the following proposition:

Proposition 4 On top of assumptions from Prop. 3, assume h is (m + 2)-times differentiable with bounded derivatives. Then

f (x λ ) = f (x * ) + O(λ m+1 ).
Thus, within the smoothing technique, if we consider λ ∝ 1/k 2/(m+2) , to balance the terms 1/(λk 2 ) and λ m+1 , we get an error for the non-smooth problem of 1/k 2(m+1)/(m+2) , which can get arbitrarily close to 1/k 2 when m gets large. The downsides are that (a) the constants in front of the asymptotic equivalent may blow up (a classical problem in high-order expansions), and (b) m-step extrapolation requires running the algorithm m times (this can be down in parallel). In our experiments below, 3-step extrapolation already brings in most of the benefits.

Experiments. We consider the penalized Lasso problem:

min x∈R d 1 2n n i=1 (b i -x a i ) 2 + λ x 1 , where (a i , b i ) ∈ R d × R for i = 1, .
. . , n, for d = 100 and n = 100, and with input data distributed as a standard normal vector. We use either a dual entropic penalty or a dual quadratic penalty for smoothing each component |x j | of the 1 -norm of x. Plots for the quadratic penalty are presented here in Figure 7, while plots for the entropic penalty are presented in Figure 9 in Appendix D with the same conclusions.

In the left plot of Figure 7, we illustrate the dependence of f (x λ ) -f (x * ) on λ for Richardson extrapolation with various orders, while in the right plot of Figure 7, we study the effect of extrapolation to solve the non-smooth problem. For a series of regularization parameters equal to 2 i for i between -18 and 1 (sampled every 1/5), we run accelerated gradient descent on h + g λ and we plot the value of f (x) -f (x * ) for the various estimates, where for each number of iterations, we minimize over the regularization parameter. This is an oracle version of varying λ as a function of the number of iterations (a detailed evaluation where λ depends on k could also be carried out). In Figure 7, we plot the excess function values as a function of the number of iterations, taking into account that m-step Richardson extrapolation requires m-times more iterations. We see that we get a strong improvement approaching 1/k 2 .

-15 -10 From non-linear programming to linear programming. When we use the entropic penalty, the smoothing framework is generally applicable in most non-linear programming problems (see, e.g., [START_REF] Cominetti | Stable exponential-penalty algorithm with superlinear convergence[END_REF]. It is interesting to note that typically when applying the entropic penalty, the deviation to the global optimizer is going to zero exponentially in -1/λ for some of the components (see a proof for our particular case in Appendix D), but not for the corresponding dual problem (which is our primal problem).

Another classical instance of entropic regularization in machine learning leads to the Sinkhorn algorithm for computing optimal transport plans [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. For that problem, the entropic penalty is put directly on the original problem, and the deviation in estimating the optimal transport plan can be shown to be asymptotically exponential in -1/λ [START_REF] Cominetti | Asymptotic analysis of the exponential penalty trajectory in linear programming[END_REF], and thus there Richardson extrapolation is not helpful (unless one wants to estimate the Kantorovich dual potentials).

Improving bias in ridge regression

We consider the ridge regression problem, that is, w λ is the unique minimizer of min

w∈R d 1 2n y -Φw 2 2 + λ 2 w 2 2 ,
where Φ ∈ R n×d is a feature vector and y ∈ R n a vector of responses [START_REF] Friedman | The Elements of Statistical Learning[END_REF]. The solution may be obtained in closed form by solving the normal equations, as

w λ = (Φ Φ + nλI) -1 Φ y.
The regularization term λ 2 w 2 2 is added to avoid overfitting and control the variability of w λ due to the randomness in the training data (the higher the λ, the more control); however, it does create a bias that goes down as λ goes to zero. Richardson extrapolation can be used to reduce this bias. We thus consider w

(1) λ = 2w λ -w 2λ and more generally

w (m) λ = m i=0 α (m) i w iλ ,
with the same weights as defined in Eq. (1) and Eq. ( 2). In order to compute w (m) λ , either m ridge regression problems can be solved, or a closed-form spectral formula can be used based on a single singular value decomposition of the kernel matrix (see Appendix E.3 for details).

Theoretical analysis. Following [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF], we assume for simplicity that Φ is deterministic and that y = z + ε, where ε has zero mean and covariance matrix σ 2 I. We consider the in-sample error of ŷλ = Φw λ = K(K + nλI) -1 y = Ĥλ y, where K = ΦΦ is the usual kernel matrix, and Ĥλ is the smoothing matrix, which is equal to I for very small λ and equal to zero for very large λ. We consider the so-called "in-sample generalization error", that is, we want to minimize

1 n E ŷλ -z 2 2 = bias( Ĥλ ) + variance( Ĥλ ), where bias( Ĥλ ) = 1 n ( Ĥλ -I)z 2 2 variance( Ĥλ ) = σ 2 n tr Ĥ2 λ .
The bias term is increasing in λ, while the variance term in decreasing in λ, and there is thus a trade-off between these two terms. To find the optimal λ, assumptions need to be made on the problem, regarding the eigenvalues of K, and the components of z in the eigenbasis of K. That is, following the notations of Bach (2013, Section 4.3), we assume that the eigenvalues of K are Θ(nµ i ) (that is bounded from above and below by constants times nµ i ) and the coordinates of z in the eigenbasis of K are Θ( √ nν i ). The precise trade-off depends on the rates at which µ i and ν i decay to zero. A classical situation is µ i ∼ i -2β and ν i ∼ i -2δ , where β > 1/2 and δ > 1/2 (to ensure finite energy). As detailed in Appendix E:

• the variance term is equivalent to σ 2 n λ -1/2β and does not depend on z or δ; • the bias term depends on both δ and β: for signals which are not too smooth (i.e., not too fast decay of ν i , and thus small δ), that is if δ < 2β + 1/2, then the bias term is equivalent to λ (2δ-1)/2β and we can thus find the optimal λ as (σ 2 /n) β/δ leading to predictive performance of (σ 2 /n) 1-1/2δ , which happens to be optimal [START_REF] Caponnetto | Optimal rates for the regularized least-squares algorithm[END_REF]. However, when δ > 2β + 1/2, a phenomenon called "saturation" occurs, and the bias term is equivalent to λ 2 (independent of δ), and the optimized predictive performance is (σ 2 /n) 1-1/(4β+1) , which is not optimal anymore.

As shown in Appendix E, by reducing the bias, with m-step Richardson interpolation, we can show that the variance term is bounded by a constant times the usual one, while the bias term is equivalent to λ (2δ-1)/2β for a wider range of δ, that is, δ < 2(m+ 1)β + 1/2, which recovers for m = 0 the non-extrapolated estimate. This leads to optimal statistical performance for a wider range of problems.

Experiments. As an illustration, we consider a ridge regression problem with data uniformly sampled on the unit sphere in dimension d = 40 with n = 200 observations, and y generated as a linear function of the input plus some noise. We consider the rotation invariant kernel equal to the expectation k(x, x ) = E τ (1 + x x )σ(τ x)σ(τ x ), for τ uniform on the sphere. This is equal to, for σ(α) = 1 α>0 (see [START_REF] Cho | Kernel methods for deep learning[END_REF][START_REF] Bach | Breaking the curse of dimensionality with convex neural networks[END_REF]: ))

k(x, y) ∝ (1 + x x ) π -arccos(x x ) .
Figure 8: Left: regularization path for the classical iterate w λ , one step of Richardson w

(1)

λ and 3 steps w

(3) λ Right: optimal error as a function of the order of the Richardson step.

When the number of observations n tends to infinity, the eigenvalues of 1 n K are known to converge to the eigenvalues of a certain infinite-dimensional operator [START_REF] Koltchinskii | Random matrix approximation of spectra of integral operators[END_REF]. As shown by [START_REF] Bach | Breaking the curse of dimensionality with convex neural networks[END_REF], the corresponding eigenvalues of the kernel matrix decay as i -1-1/d . We consider z generated as a linear function so that ν i has a finite number of non zero components in the eigenbasis of K.

In the left plot of Figure 8, we consider 1-step and 3-step Richardson extrapolation and plot the generalization error (averaged over 10 replications) as a function of the regularization parameter: we can see that as expected, (a) with extrapolation the curves move to the right (we can use a larger λ for a similar performance, which is advantageous as iterative algorithms are typically faster), and (b) the minimal error is smaller (which is true here because we learn a smooth function). In the right plot of Figure 8, we study the effect of increasing the order m of extrapolation, showing that the larger the better with some saturation. With m infinite, there will be overfitting as the corresponding spectral filter is non-stable, but this happens very slowly (see Appendix E.3 for details).

Conclusion

In this paper, we presented various applications of Richardson extrapolation to machine learning optimization and estimation problems, each time with an asymptotic analysis showing the potential benefits. For example, when using the number of iterations of an iterative algorithm to perform extrapolation, we can accelerate Frank-Wolfe algorithms to have an asymptotic rates of O(1/k 2 ) on polytopes for the step-size ρ k = 1/k and locally strongly-convex functions (this is achieved without extrapolation for the step-size ρ k = 2/(k + 1)). When extrapolating based on the regularization parameter, we can accelerate Nesterov smoothing technique to have asymptotic rates close to O(1/k 2 ).

We also highlighted situations where Richardson extrapolation does not bring any benefits (but does not degrade performance much), namely when applied to accelerated gradient descent or the Sinkhorn algorithm for optimal transport.

The analysis in this paper can be extended in a number of ways: (1) while the paper has focused on asymptotic analysis for simplicity, non-asymptotic analysis could be carried out to study more finely when acceleration starts, (2) we have focused on deterministic optimization algorithms, and extensions to stochastic algorithms could be derived, along the lines of the work of [START_REF] Dieuleveut | Bridging the gap between constant step size stochastic gradient descent and Markov chains[END_REF], (3) we have primarily focused on convex optimization algorithms but non-convex extensions, like done by [START_REF] Scieur | Nonlinear acceleration of deep neural networks[END_REF] for Anderson acceleration, could also lead to acceleration

A Preliminary considerations

We first start with lemmas that we will need in subsequent proofs; the second one shows strong convexity on a level set once we assume that the Hessian at optimum is positive definite.

Lemma 1 Assume f convex, three-times differentiable with bounded third-order derivatives, and a point x * such that f (x * ) is positive definite. Then there exists c > 0 such that for any x ∈ R d .

1 2 (x -x * ) f (x * )(x -x * ) c ⇒ f (x) 1 2 f (x * ). Proof Since f (x * ) is positive definite, λ min (f (x * )) > 0, and 1 2 (x -x * ) f (x * )(x -x * ) c implies that x -x * 2 2 2c λmin(f (x * ))
. Thus, since third order derivatives of f are bounded, if

1 2 (x -x * ) f (x * )(x -x * ) c, we have f (x) -f (x * )
op Ac for some constant A, and thus

f (x) f (x * ) -AcI f (x * ) - Ac λ min (f (x * )) f (x * ) 1 2 f (x * ) if Ac 1 2 λ min (f (x * )).
Lemma 2 Assume f convex, three-times differentiable with bounded third-order derivatives, and that x * ∈ R d is a minimizer of f on R d , such that f (x * ) is positive definite. Then x * is the unique minimizer of f and there exists c > 0 such that for any

x ∈ R d , f (x) -f (x * ) c ⇒ f (x) 1 2 f (x * ).
Proof Using the above Lemma 1, then there exists c > 0 such that

1 2 (x -x * ) f (x * )(x -x * ) c ⇒ f (x) 1 2 f (x * ).
We now show that f (x)-f (x * ) c/2 implies that we stay in the region where 1 2 (x-x * ) f (x * )(xx * ) c. Indeed, using the Taylor formula with integral remainder, for 1 2 (x -x * ) f (x * )(x -x * ) c (where the lower bound on Hessians above holds), we have

f (x) -f (x * ) = 0 + 1 0 (x -x * ) f (x * + t(x -x * ))(x -x * )(1 -t)dt 1 2 1 0 (1 -t)dt (x -x * ) f (x * )(x -x * ) = 1 4 (x -x * ) f (x * )(x -x * ) and thus f (x) -f (x * ) c/2 implies 1 2 (x -x * ) f (x * )(x -x * ) c
, and we get the desired result.

B Proof of Prop. 1 (averaged gradient descent)

In this particular case of unconstrained gradient descent, f

(x k ) -f (x * ) 1 γk x 0 -x * 2
, as soon as γ 1/L [START_REF] Nesterov | Introductory Lectures on Convex Optimization: A Basic Course[END_REF]. This implies from Lemma 2 in Appendix A, that for k larger than some k 0 , all iterates are such that f (x k ) 1 2 f (x * ), and thus, after that k, we are in the strongly convex case, where x k -x * c ρ k , where c , ρ depend on the lowest eigenvalue µ of f (x * ) as c = x k0 -x * ρ -k0 and ρ = (1 -γµ/2).

Thus, with xk =

1 k k-1 i=0 x i , k(x k -x * ) tends to ∞ i=0 (x i -x * ) when k → +∞ (since the series is convergent), and k(x k -x * ) - ∞ i=0 (x i -x * ) = - ∞ i=k (x i -x * ), leading to k(x k -x * ) - ∞ i=0 (x i -x * ) = O(ρ k ), and thus, with ∆ = ∞ i=0 (x i -x * ) (which is hard to compute a priori), xk = x * + 1 k ∆ + O(ρ k ).
C Proof of Prop. 2 (Frank-Wolfe)

Preliminary remarks. We consider the step-sizes ρ k = 1 k and ρ k = 2 k+1 , for which the respective convergence rates for f (x k ) -f (x * ) are of the form c k and c log k k , for constants c depending on the smoothness of f and the diameter of the compact set (see, e.g., [START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF]. When running Frank-Wolfe (with any of the classical versions with open-loop step-sizes), we thus have f

(x k ) -f (x * ) = O((log k) β /k) with β ∈ {0, 1}.
Because of affine invariance of the Frank-Wolfe algorithm (and because f (x * ) is invertible), we can assume without loss of generality that f (x * ) = I. Moreover, the constraint qualification implies that f (x * ) = 0, thus, using Taylor expansion with integral remainder like in Lemma 2, if x -x * 2 is small enough then

f (x) -f (x * ) f (x * ) (x -x * ) + 1 4 x -x * 2 2 .
Since x * is the minimizer of f on K and x ∈ K, then f (x * ) (x -x * ) 0 and thus, we have that

x k -x * = O((log k) β/2 / √ k)
and (from Lemma 1), that f is locally strongly-convex.

Analysis of Frank-Wolfe. The assumption which is made implies that there exists σ > 0 such that the ball of center x * and radius σ intersected with the affine hull of y 1 , . . . , y m is included in the convex hull of y 1 , . . . , y m , as well α ∈ (0, 1) such that if cos(f (x * ), z) 1 -α, then min y∈K z y is attained only by elements of the convex hull of y 1 , . . . , y m . Thus, for k large enough, that is greater than some k 0 (which can be quantified from σ, α and other quantities), all elements of arg min y∈K f (x k-1 ) y are in the convex hull of y 1 , . . . , y m , that is, only the correct extreme points are selected. Denoting Π the orthogonal projection on the span of y 1 -x * , . . . , y m -x * , we have:

x k = (1 -ρ k )x k-1 + ρ k xk where xk ∈ arg min y∈K f (x k-1 ) y,
and thus, subtracting x * :

x k -x * = (1 -ρ k )(x k-1 -x * ) + ρ k (x k -x * ),
leading to, using the projections Π and I -Π:

Π(x k -x * ) = (1 -ρ k )Π(x k-1 -x * ) + ρ k (x k -x * ), and, because xk -x * is in the span of y 1 -x * , . . . , y m -x * , (I -Π)(x k -x * ) = (1 -ρ k )(I -Π)(x k-1 -x * ).
We now consider these two terms separately.

Convergence of (I -Π)(x k -x * ). For ρ k = 2 k+1 , we have, in closed form for k k 0 :

(I -Π)(x k -x * ) = k -1 k + 1 (I -Π)(x k-1 -x * ) = k 0 (k 0 + 1) k(k + 1) (I -Π)(x k0 -x * ).
For ρ k = 1 k , we have:

(I -Π)(x k -x * ) = k -1 k (I -Π)(x k-1 -x * ) = k 0 k (I -Π)(x k0 -x * ).
Convergence of Π(x k -x * ). We now look at the convergence of Π(x k -x * ). We have

Π(x k -x * ) 2 = (1 -ρ k )Π(x k-1 -x * ) + ρ k (x k -x * ) 2 = (1 -ρ k ) 2 Π(x k-1 -x * ) 2 + ρ 2 k xk -x * 2 + 2(1 -ρ k )ρ k (x k -x * ) Π(x k-1 -x * ).
Because of the ball assumption (that is the existence of a ball around x * that is contained in the supporting face of K), we have

f (x k-1 ) (x k-1 -xk ) = max y∈K f (x k-1 ) (x k-1 -y) f (x k-1 ) (x k-1 -x * ) + Πf (x k-1 ) σ, which leads to f (x k-1 ) (x * -xk ) Πf (x k-1 ) σ.
Moreover, using a Taylor expansion with f (x * ) = I, we have f

(x k-1 ) = f (x * ) + f (x * )(x k-1 - x * ) + O((log k) β /k) = f (x * ) + (x k-1 -x * ) + O((log k) β /k). Moreover, by optimality of x * , f (x * ) (x k -x * ) = 0. Therefore we have, for ρ k = O(1/k): Πx k -x * 2 (1 -ρ k ) 2 Πx k-1 -x * 2 + ρ 2 k diam(K) 2 -2(1 -ρ k )ρ k Πf (x k-1 ) σ + O((log k) β /k 2 ) (1 -ρ k ) 2 Πx k-1 -x * 2 -2(1 -ρ k )ρ k Πx k-1 -x * σ + O((log k) β /k 2 ).
For ρ k = 1/k, we get

k 2 Πx k -x * 2 (k -1) 2 Πx k-1 -x * 2 -2 (k -1) Πx k-1 -x * 2 σ + O(log k).
For ρ k = 2/(k + 1), we get

k 2 Πx k -x * 2 (k -1) 2 Πx k-1 -x * 2 -2 (k -1) Πx k-1 -x * 2 σ + O(1).
We then can use the following simple lemma on sequences 2 : if u k 0 such that u 0 = 0 and

u k u k-1 -A √ u k-1 + Bv k-1 for (v k ) non-decreasing and positive, then u k B 2 A 2 v 2 k + Bv k . This leads to a bound in O(1) for k 2 Πx k -x * 2 for ρ k = 2/(k + 1), and in O((log k) 2 ) for ρ k = 1/k.
Note that this is similar to the proof of the convergence of Frank-Wolfe algorithms to an interior point of the feasible set, for which we have a rate of convergence of f [START_REF] Chen | Super-samples from kernel herding[END_REF].

(x k ) -f (x * ) = O(1/k 2 )
2 The proof is by induction. This is true for k = 0. If this is true for k -1, then either (a) u k-1

B 2 A 2 v k-1 , then u k u k-1 B 2 A 2 v 2 k + Bv k because v k-1 v k , or (b) u k-1 B 2 A 2 v k-1 and then u k u k-1 + Bv k-1 B 2 A 2 v 2 k-1 + Bv k-1 B 2 A 2 v 2 k + Bv k .
Putting things together. We then have for ρ k = 1 k ,

x k = x * + k 0 k (I -Π)x k0 + O((log k) 2 /k 2 ), which leads to f (x k ) -f (x * ) = k0 k f (x * ) (I -Π)x k0 + O((log k) 2 /k 2
), which can be put back into the original bound, to obtain the bound without the logarithmic factor (since we have now replaced O(log(k)/k) by O((log k) 2 /k 2 ) = o(1/k) in the start of the proof). The dependence in k can here lead to an acceleration.

For

ρ k = 2 k+1 , x k = x * + k 0 (k 0 + 1) k(k + 1) (I -Π)x k0 + O(1/k 2 ), which leads to f (x k ) -f (x * ) = k0(k0+1) k(k+1) f (x * ) (I -Π)x k0 + O(1/k 2 ).
The dependence in k does not lead to an acceleration.

Note that the terms in O(1/k 2 ) are not amenable to Richardson extrapolation because they are oscillating.

D Proofs of Prop. 3 and Prop. 4 (Nesterov smoothing)

We denote by x λ the minimizer of h(x) + g λ (x), and η λ the corresponding dual variable. The dual variable η λ is unique because of the strong-convexity of ϕ, while the primal variable is unique due to the same reasoning as in Appendix A (when λ tends to zero, x λ has to be close to x * , and in the neighborhood of x * , h is strongly-convex).

The primal problem is min

x∈R d h(x) + λϕ * Ax -b λ ,
while the dual problem is max η∈∆m -λϕ(η) -h * (-A η) -η b.

The primal and dual solutions x λ and η λ are related through duality for ϕ, that is,

η λ = ∂ϕ * Ax λ -b λ ,
and through duality for h, that is,

x λ = ∂h * (-A η λ ).
Since we consider functions ϕ which are uniformly bounded, then we know that f (x λ )-f (x * ) = O(λ) and thus because the function is locally strongly convex, we have:

x λ -x * = O(λ 1/2 ).

D.1 Quadratic penalty

With a quadratic penalty, and small enough λ, the solution η λ will have the same sparsity pattern, using standard techniques from active set methods [START_REF] Nocedal | Numerical Optimization[END_REF]-that is, show that the dual solution constrained to the same active set leads to a globally optimal primal/dual solution. Moreover, because, once restricted to I, the dual function is locally strongly convex, and because ϕ is bounded, the deviation in dual function values is less than O(λ) and thus η λ -η * = O(λ 1/2 ) (which is a bound we are going to improve below).

The optimality conditions for the dual problem become (stationarity with respect to η I ):

0 = -λ(η λ ) I + A I ∂h * (-A I (η λ ) I ) -b I .
Since h is twice differentiable at x * , h (x * ) is invertible, and x * ∈ ∂h * (-A η * ), by the implicit function theorem, h * is twice differentiable at -A η * , and its Hessian is h (x * ) -1 . We can thus further expand the optimality condition above as ). Since we know already that η λ -η * = O(λ 1/2 ), this leads to η λ -η * = O(λ), which in turn leads to

0 = -λ(η λ ) I + A I ∂h * (-A I (η * ) I ) -A I ∂ 2 h * (-A I (η * ) I )A I ((η λ ) I -(η * ) I ) -b I + O( η λ -η * 2 ). Since ∂ 2 h * (-A I (η * ) I ) = h (x * ) -
(η λ ) I = (η * ) I -λ A I h (x * ) -1 A I -1 (η * ) I + O(λ 2 ),
which is the desired expansion for the dual variable. We then get:

x λ = ∂h * (-A η λ ) = x * + λh (x * ) -1 A I A I h (x * ) -1 A I -1 (η * ) I + O(λ 2 ) = x * + λ∆ + O(λ 2 ).
Thus, 2x λ -x 2λ = O(λ 2 ), and

h(x λ ) + g(x λ ) = h(x * ) + g(x * ) + λh (x * ) ∆ + g(x * + λ∆) -g(x * ) + O(λ 2 ) = h(x * ) + g(x * ) + g(x * + λ∆) -g(x * ) + λh (x * ) ∆ + O(λ 2 ).
Since by optimality h (x * ) ∈ -∂g(x * ), the term g(x * + λ∆) -g(x * ) + λh (x * ) ∆ resembles a Taylor expansion of g at x * , but in general, we cannot have a term in O(λ 2 ) because of the non-smoothness of g. For example, for the 1 -norm, we get g(x * + λ∆) -g(x * ) + λf (x * ) ∆ = λ ∆ I c 1 + λf (x * ) ∆, is not zero, and only O(λ).

For Richardson extrapolation, we get:

h(2x λ -x 2λ ) + g(2x λ -x 2λ ) = h(x * ) + g(x * ) + O(λ 2 ),
and thus an improvement from O(λ) to O(λ 2 ).

D.2 Entropic penalty

For ϕ(η) = m i=1 η i log η i -η i , we have ϕ (η) i = log η i , and we cannot use anymore thet fact that η λ has the same sparsity pattern as η * , since all components of η λ are non zero. However, since the entropy is bounded over the simplex, we still have η λ -η * = O(λ 1/2 ), and from the same reasoning as for the quadratic penalty, x λ -x * = O(λ 1/2 ).

Thus, by writing primal-dual optimality conditions, we get:

-λ log η λ + Ax λ -b = 0 h (x λ ) + A η λ = 0.
This implies that for i / ∈ I,

log η i ∼ (Ax * -b) i -sup j (Ax * -b j ),
which is strictly negative by assumption. Thus (η I c ) λ = O(ρ -1/λ ) for a certain ρ ∈ (0, 1). We now show that (η I ) λ = (η * ) I + λ∆ + O(λ 2 ). From the optimality conditions, we get:

c ) λ = O(ρ -1/λ ), this leads to -λ(log(η) * ) I -A I h (x * ) -1 A I ((η λ ) I -(η * ) I ) = O( (η λ ) I -(η * ) I 2 )+O(λ (η λ ) I -(η * ) I )+O(ρ -1/λ ),
which in turn leads to

(η λ ) I = (η * ) I -λ A I f (x * ) -1 A I -1 (log(η * )) I + O(λ 2 ),
which is the desired expansion for the dual variable. We then get:

x λ = ∂h * (-A η λ ) = x * + λh (x * ) -1 A I A I h (x * ) -1 A I -1 (log(η * )) I + O(λ 2 ) = x * + λ∆ + O(λ 2 ).
Note that following [START_REF] Cominetti | Stable exponential-penalty algorithm with superlinear convergence[END_REF], we could get a single proof for all ϕ(η) = n i=1 ψ(η i ), with replacing log η * by ψ (η * ), with a condition to ensure that the zero variables lead to vanishing terms.

D.3 Higher-order expansions (Prop. 4)

In order to use m-step Richardson extrapolation, we need to have a bound of the form

x λ = x * + m i=1 ∆ i λ i + O(λ m+1 ).
We only consider for simplicity the quadratic penalty, where we simply need an expansion of (η λ ) I (minor modifications would lead to a proof for the entropic penalty since there η I c is exponentially small).

The expansion can be obtained from the implicit function theorem applied to the equation 0 = -λ(η λ ) I + A I ∂h * (-A I (η λ ) I ) -b I , which is of the form H((η λ ) I , λ) = 0, where H has highorder derivatives as long as h * is sufficiently differentiable, and the partial derivative with respect to (η λ ) I is an invertible matrix. The high-order differentiability of h * around -A η * comes from the implicit function theorem, applied to the definition of the gradient of the Fenchel conjugate.

E Ridge regression E.1 Standard extrapolation

We have ŷλ = K(K + nλI) -1 y = Ĥλ y, and thus 2ŷ λ -ŷ2λ = [(2 Ĥλ -Ĥ2λ )]y, and we can compute explicitly

2 Ĥλ -Ĥ2λ -I = 2K(K + nλI) -1 -K(K + 2nλI) -1 -I = 2nλ (K + 2nλI) -1 -(K + nλI) -1 = -2(nλ) 2 (K + nλI) -1 (K + 2nλI) -1 . Thus, bias(2 Ĥλ -Ĥ2λ ) 4n 3 λ 4 z (K + nλI) -4 z variance(2 Ĥλ -Ĥ2λ ) σ 2 n 4 tr[K(K + nλI) -1 ] 2 + tr[K(K + 2nλI) -1 ] 2 5σ 2 n tr[K(K + nλI) -1 ] 2 .
We have bias(2 Ĥλ -Ĥ2λ ) 4bias( Ĥλ ) and variance(2 Ĥλ -Ĥ2λ ) 5variance( Ĥλ ), so the two terms never incur more than a constant factor. However, the bias can be much improved. Following Bach (2013, Section 4.3), if the eigenvalues of K are Θ(nµ i ) and the coordinates of z in the eigenbasis of K are Θ( √ nν i ), we can compute equivalents (up to constant terms) of the bias and variance terms, for different types of decays. See Table 1 with m = 0. Since the optimal predictive performance is ( σ 2 n ) 1-1/(2δ) , the only potential gains to go from m = 0 (no extrapolation) to m > 0 (extrapolation) occur when (ν i ) has a fast decay (that is, first, third and fifth lines).

For the first line in Table 1, we will show that the bias term for Richardson extrapolation is of the order λ 4 if 2δ > 8β + 1, and equal to λ (2δ-1)/2β when 2δ < 8β + 1. We will thus increase the regime of validity of the bias term that leads to optimal performance, from 2δ < 4β+1 to 2δ < 8β+1. More generally, as we show below in Appendix E.2, the bias for m-step extrapolation is proportional to n 2m+1 λ 2m+2 z (K + nλI) -2m-2 z, and we bound it directly following closely the computations of Bach (2013, Appendix C.2):

n 2m+1 λ 2m+2 z (K + nλI) -2m-2 z = n 2m+2 λ 2m+2 n i=1 ν i (nµ i + nλ) 2m+2 = λ 2m+2 n i=1 ν i (µ i + λ) 2m+2 = λ 2m+2 n i=1 i -2δ (i -2β + λ) 4 2λ 2m+2 n 1 t -2δ (t -2β + λ) 2m+2 dt = 2λ 2m+2 n 1 t 4(m+1)β-2δ (1 + λt 2β ) 2m+2 dt. If 2δ -4(m + 1)β > 1, then we have an upper bound of 2λ 2m+2 n 1 t 4(m+1)β-2δ dt = O(λ 2m+2 ). If 2δ -4(m + 1)β < 1, then we can further bound 2λ 2m+2 n 1 t 4(m+1)β-2δ (1 + λt 2β ) 2m+2 dt = 2λ 2m+2 λn 2β λ [(u/λ) 1/2β ] 4(m+1)β-2δ+1 (1 + u) 2m+2 1 2β du with the change of variable u = λt 2β , = 2λ 2m+2-(2m+2)+δ/β-1/(2β) λn 2β λ u (4(m+1)β-2δ+1))/((4(m+1)β) (1 + u) 2m+2 1 2β du = O(λ (2δ-1)/(2β) ),
because the integral in convergent. The bias term for the third and fifth lines of Table 1 needs modifications in exactly the same way the corresponding proof from Bach (2013, Appendix C.2).

Table 1: Variance, bias, optimal regularization parameter, corresponding prediction performance, for several decays of eigenvalues and signal coefficients, and m-th order Richardson extrapolation (we always assume δ > 1/2, β > 1/2, ρ > 0, κ > 0, to make the series summable). All entries are functions of i, n or λ and are only asymptotically bounded below and above, i.e., correspond to the asymptotic notation Θ(•). Adapted from [START_REF] Bach | Sharp analysis of low-rank kernel matrix approximations[END_REF], changes due to potential Richardson extrapolation are highlighted in red.

In order to find the optimal regularization parameter, we minimize with respect to λ, that leads to λ 2(m+1)+1/2β ∝ (σ 2 /n), leading to an optimal prediction performance proportional to (σ 2 /n) τ , with τ = 2(m+1)

2(m+1)+ 1 2β = 1 - 1 2β 2(m+1)+ 1 2β = 1 - 1 4(m+1)β+1 .

E.2 Multiple extrapolation steps

Using m steps of Richardson interpolation, we prove that we can get this regime as long as 2δ < 4(m + 1)β +1, and thus with no limit if m is large enough.

We thus consider Ĥ(m) = (K + nλI) -1/2 I + nλ(i -1)(K + λI) -1 -1 (K + nλI) -1/2 = (K + nλI) -1/2 m-1 j=0 (-1) j nλ(K + nλ(i -1)I) -1 j (K + nλI) -1/2 +(-1) m (K + nλI) -1/2 nλ(i -1)(K + nλI) -1 m I + nλ(K + λ(i -1)I) -1 -1 (K + nλI) -1/2 = m-1 j=0 (-1) j (nλ) j (K + nλI) -j-1 (i -1) j + (-1) m (nλ) m (K + nλI) -m (i -1) m (K + nλiI) -1 , where we have used the sum of a geometric series, with A = nλ(i -1)(K + λI) Similarly to the case m = 1, ( Ĥ(m) λ ) 2 is upper-bounded by a sum of matrices which are all less than K 2 (K + nλI) -2 (for the order between symmetric matrices), leading to a bound on the variance term as:

variance( Ĥ(m) λ ) = σ 2 tr[ Ĥ(m) λ ] 2 m σ 2 n tr[K(K + nλI) -1 ] 2 .
Here m and m are constants that could be explicitly computed. As shown in the section, these constants have to diverge when m tends to +∞.

E.3 Explicit expression

Expression for α (m) i

. We first give a proof for the explicit expression for α (m) i

. One approach is to solve Vandermonde matrices like done by [START_REF] Pagès | Multi-step Richardson-Romberg extrapolation: remarks on variance control and complexity[END_REF] in a similar context, but given the conjecture, we can simply check that it satisfies Eq. (1) and Eq. (2).

For Eq. (1), we have:

m+1 i=1 (-1) i-1 m + 1 i = 1 - m+1 i=0
(-1) i m + 1 i = 0, using the binomial formula. For Eq. ( 2), we have:

m+1 i=1 (-1) i-1 m + 1 i i = m+1 i=1
(-1) i-1 (m + 1) m i -1 = (m + 1) × 0 = 0, and, more generally for any j ∈ {1, . . . , m}, m+1 i=j (-1) i-1 m + 1 i i(i -1) • • • (i -j + 1) = m+1 i=j (-1) i-1 (m + 1)m • • • (m -j + 2) m -j + 1 i -j = 0, also using the binomial formula. This shows that m+1 i=j (-1) i-1 m+1 i i j = 0 for all j ∈ {1, . . . , m}, which finishes the proof of the formula for α where s : R n×n → R n×n is a spectral function defined on symmetric matrices from a function (note the classical overloaded notation) s : R → R, by keeping eigenvectors unchanged and applying s to eigenvalues (Golub and Van Loan, 1989, Chapter 11). We have, using a representation by an integral and the binomial formula:

s(µ) = 1 - m+1 i=1 i(-1) i-1 m + 1 i 1 µ + i = 1 - m+1 i=1
(m + 1)(-1) i-1 m i -1

1 0 t µ+i-1 dt = 1 -(m + 1) 1 0 t µ m+1 i=1 (-1) i-1 m i -1 t i-1 dt = 1 -(m + 1) 1 0 t µ (1 -t) m dt = 1 -(m + 1) Γ(1 + µ)Γ(1 + m) Γ(2 + m + µ) ,
using the expression of the Beta function in term of the Gamma function Γ [START_REF] Abramowitz | Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables[END_REF]. We can simply the expression as follows:

s(µ) = 1 - (m + 1)! (µ + 1)(µ + 2) • • • (µ + m + 1)
.

This provides a new closed-form expression for Richardson extrapolation, as well as it provides an equivalent when m tends to +∞, as s(µ) ∼ 1 -Γ(1+µ) m µ . Therefore, when m → +∞, and µ > 0, then s(µ) tends to one (which implies that the constants m and m cannot remain bounded). Therefore, the variance term converges to σ 2 , but very slowly: the method does not blow up, but does not learn either.

Figure 3 :

 3 Figure 3: Left: averaged gradient descent on a logistic regression problem in dimension d = 400, and with n = 4000 observations min x∈R d f (x) = 1

Figure 7 :

 7 Figure7: Richardson extrapolation for Nesterov smoothing on a penalized Lasso problem, with regularization by the quadratic penalty. Left: dependence of f (x λ ) -f (x * ) on λ, for Richardson extrapolation of order m, we indeed recover a slope of m + 1 in the log-log plot. Right: Optimization error vs. number of iterations; where we go from a slope of -1 (blue curves) to improved slopes of -4/3 (red curve) and -8/5 (black curve). See text for details.

  1 and b I = A I ∂h * (-A I (η * ) I ) (because of optimality of x * and η * ), this leads to 0 = -λ(η λ ) I + A I h (x * ) -1 A I ((η λ ) I -(η * ) I ) + O( η λ -η * 2

0=Figure 9 :

 9 Figure9: Richardson extrapolation for Nesterov smoothing on a penalized Lasso problem, with regularization by the entropic penalty. Left: dependence of f (x λ ) -f (x * ) on λ, for Richardson extrapolation of order m, we indeed recover a slope of m + 1 in the log-log plot. Right: Optimization error vs. number of iterations; where we go from a slope of -1 (blue curves) to improved slopes of -4/3 (red curve) and -8/5 (black curve). See text for details.

  j A j + (-1) m (I + A) -1 A m . j (nλ) j (K + nλI) -j-1 (i -1) j m (nλ) m (K + nλI) -m (i -1) m (K + nλiI) -1 .The first term is exactly zero by definition of the Richardson weights α 1) m (K + nλI) -m (K + nλiI) -1 .Thus (nλ) -m-1 (K + nλI) -m/2-1/2 Ĥ(m) λ -I (K + nλI) -m/2-1/2has an operator norm bounded by a constant that depends on m (and not on other quantities like λ or n). This allows to bound the bias as biasz 2 m n 2m+1 λ 2m+2 z (K + nλI) -2m-2 z.

.

  Expression for the smoothing matrix. We can now provide an explicit expression for the extrapolated smoothing matrix. We have: nλiI) -1 = s(K/(nλ)),

Note that: (a) the lower bound with dependence O(1/k) from Canon and Cullum (1968) only applies to Frank-Wolfe algorithms with line-search, and (b) that our bounds are local and the constants have to depend on dimension so as to not contradict the lower bound from[START_REF] Jaggi | Revisiting Frank-Wolfe: Projection-free sparse convex optimization[END_REF].

Acknowledgements

This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute). We also acknowledge support the European Research Council (grant SEQUOIA 724063).

(µ i ) (ν i ) var. bias optimal λ pred. perf. condition