domain: sde stat

text:
Hervé Douville
email: herve.douville@meteo.fr

Katharine M Willett

A drier than expected future, supported by near-surface relative humidity observations

Teaser (one-sentence summary): Climate models underestimate both recent and future nearsurface air drying over land, even when

Revised Manuscript: Confidential Template revised February 2021 Almost ten years ago, a commentary in Science (1) warned that changes in precipitation do not tell the whole story of hydrological changes caused by human-induced global warming: "Many regions will get more rain, but it appears that few will get enough to keep pace with the growing evaporative demand." This suggestion was grounded in the predicted decrease in annual mean near-surface relative humidity (RH) over most land regions. It was heeded by the first working group (WG1) of the Intergovernmental Panel on Climate Change (IPCC) which released last year its sixth Assessment Report (AR6). Chapter 8 pointed out that water cycle changes result not only from changes in precipitation patterns, but also from an increase in the atmospheric evaporative demand due to greater warming over land than over the ocean and, consequently, a reduction in continental near-surface RH [START_REF] Douville | IPCC AR6 WG1[END_REF][START_REF] Allan | Advances in understanding large-scale responses of the water cycle to climate change[END_REF][START_REF] Byrne | Understanding Decreases in Land Relative Humidity with Global Warming: Conceptual Model and GCM Simulations[END_REF]). Yet, both observed and simulated changes in evapotranspiration, and in precipitation minus evaporation (P-E), are still very uncertain, even when globally averaged over land [START_REF] Douville | IPCC AR6 WG1[END_REF][START_REF] Robertson | Reconciling Land-Ocean Moisture Transport Variability in Reanalyses with P-ET in Observationally Driven Land Surface Models[END_REF][START_REF] Elbaum | Uncertainty in projected changes in precipitation minus evaporation: Dominant role of dynamic circulation changes and weak role for thermodynamic changes[END_REF]. Structural and parametric model uncertainties still dominate the overall uncertainty in climate projections [START_REF] Douville | IPCC AR6 WG1[END_REF][START_REF] Lehner | Partitioning climate projection uncertainty with multiple large ensembles and CMIP5/6[END_REF] and the gap between understanding and predicting water cycle changes has not been bridged from AR5 to AR6.

It is thus urgent to combine our improved theoretical understanding [START_REF] Allan | Advances in understanding large-scale responses of the water cycle to climate change[END_REF] with better constrained projections to provide more reliable climate information to decision makers. As the climate warms, both theoretical (8) and modelling [START_REF] O'gorman | How closely do changes in surface and column water vapor follow Clausius-Clapeyron scaling in climate change simulations?[END_REF] studies predict that RH will remain approximately constant at the global scale, leading to an increase in specific humidity, a positive water vapor feedback, and an overall intensification of the global water cycle [START_REF] Douville | IPCC AR6 WG1[END_REF]. An early formal attribution study identified a human-induced significant global-scale increase in surface specific humidity, with relative humidity remaining approximately constant [START_REF] Willett | Attribution of observed surface humidity changes to human influence[END_REF]. Yet, the constant near-surface RH hypothesis seems to be no more valid over land, at least since the early 21 st century [START_REF] Simmons | Low-frequency variations in surface atmospheric humidity, temperature, and precipitation: Inferences from reanalyses and monthly gridded observational data sets[END_REF]. Global projections from the 5 th phase of the coupled model project (CMIP5) showed uncertain changes with contrasted geographical patterns [START_REF] Byrne | Understanding Decreases in Land Relative Humidity with Global Warming: Conceptual Model and GCM Simulations[END_REF]. In the highest emission scenario, they revealed a strong land-ocean contrast in the response of near-surface RH to global warming, with a possible but model-dependent decrease over land. Yet, there were at least two studies suggesting that most CMIP5 models underestimate the projected drying over land given their limited ability to capture the observed negative trends in RH (i.e., a rising water vapor deficit) since the late 20 th century [START_REF] Douville | Midlatitude summer drying: An underestimated threat in CMIP5 models?[END_REF][START_REF] Dunn | Comparison of land surface humidity between observations and CMIP5 models[END_REF].

In the present study, we use a recent update of quality-controlled observations of surface temperature (HadCRUT5; 14) and RH (HadISDH; 15) to constrain global projections from the state-of-the-art CMIP6 climate models. Compared to CMIP5, CMIP6 models do not show any improvement in their ability to converge on the response of near-surface RH to global warming (Fig. S1), even in the highest emission scenario where such uncertainties cannot be attributed to internal climate variability [START_REF] Lehner | Partitioning climate projection uncertainty with multiple large ensembles and CMIP5/6[END_REF]. Over most land areas, the range of model responses varies from limited wetting to substantial drying. Only a few regions, particularly the Mediterranean basin, South Africa and Amazonia, are subject to a robust drying, whose magnitude is however still highly model-dependent. Such uncertainties remain a major obstacle for adaptation strategies, but can be potentially narrowed through the use of observational constraints. For this purpose, we use the KCC (Kriging for Climate Change) Bayesian statistical method inspired from data assimilation [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF]. This package has been recently enriched to combine two observational constraints [START_REF] Qasmi | Reducing uncertainty in local temperature projections[END_REF][START_REF] Douville | Global warming at near-constant tropospheric relative humidity is supported by observations[END_REF]. The method relies on a probabilistic assessment in which a CMIP ensemble of GCMs is taken as the "prior" distribution of projected anomalies which is then narrowed with observations. This is a strong but common assumption which neglects the possibility that the models are all wrong, especially if they share common flaws in their structure or historical 2

Revised Manuscript: Confidential Template revised February 2021 forcings, While this may lead to overconfident constrained projections, we will show that consistent results are obtained using either CMIP6 or the previous-generation CMIP5 models. Moreover, KCC also accounts for observational uncertainties and may be thus less sensitive to this assumption than previously proposed Bayesian methods. More details about the KCC method can be found in the Materials and Methods section of the Supplementary Materials.

Results

We first illustrate the method using the global-mean land surface response of RH (hereafter LSRH) projected by a subset of thirty-two CMIP6 models with only one realization of the historical experiment and the SSP5-8.5 high emission scenario (Fig. 1). It should be here emphasized that KCC aims at constraining the forced model response, after discarding the contribution of internal climate variability. In Fig. 1a, we first check that the HadCRUT5 observed estimates (200 members) of the historical (1850-2021) warming allow us to constrain the projected increase in the global-mean near-surface air temperature (GSAT). A 33% narrowing of the 5-95% confidence interval is obtained at the end of the 21 st century as well as a downward shift of the best estimate, in good agreement with a previous assessment based on a larger subset of CMIP6 models [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF].

Our results however show that the same HadCRUT5 observations poorly constrain the projections of the LSRH response. Even the sign of the constrained LSRH response remains uncertain across the 21 st century (Fig. 1b), and the confidence interval of the posterior is only narrowed by 17% at the end of the century compared to the prior. This narrowing is mostly the result of an upward shift of the lower 5% percentile. This limited constraint is consistent with the fact that the projected global warming does not tell the whole story about regional climate change [START_REF] Shepherd | Storyline approach to the construction of regional climate change information[END_REF]. Even aggregated over global land, RH is not strongly constrained by HadCRUT5 observations despite the theoretical expectation that trends in continental humidity are directly linked to the ocean warming [START_REF] Byrne | Trends in continental temperature and humidity directly linked to ocean warming[END_REF]. Similar results are found (not shown) when replacing HadCRUT5 by global sea surface temperature observations (HadISST1). This result confirms that there may be a gap between our physical understanding of climate change and our ability to deliver accurate projections [START_REF] Douville | Water remains a blind spot in climate change policies[END_REF].

The RH record from the HadISDH dataset, although being available only since 1973, provides another potential constraint on the projected LSRH anomalies (Fig. 1c). As for GSAT, we here take advantage of the estimated observational uncertainty estimate in each 5° by 5° grid cell of HadISDH to generate 200 estimates of LSRH (cf. Supplementary Materials). Compared to the prior, the HadISDH observations allow us to narrow by 22% the posterior distribution of LRSH projections at the end of the 21 st century. Beyond the reduced inter-model spread and unlike the GSAT observational constraint, the observational constraint here rules out a possible future increase in LSRH, with an ensemble mean response that is close to -3%. Combining the HadCRUT5 and HadISDH constraints has a limited effect on the ensemble mean response, but reduces the inter-model spread by 37% in 2100 (Fig. 1d). This substantial narrowing of the model uncertainties highlights the power of the KCC method and the complementarity of the proposed observational constraints, although the obtained reduction of the inter-model spread is not fully additive. Template revised February 2021

Similar conclusions can be drawn from the intermediate emission scenario (SSP2-4.5, Fig. S2) and further reinforce our confidence in the results. The constrained ensemble-mean drying found in the high-emission scenario is greater (-4% in 2100) when averaging RH in the northern midlatitudes (Fig. S3), and even stronger (-5% in 2100) when focusing on the June to September season (Fig. S4). In this region and season, where and when a large fraction of the global food production is concentrated, the land surface drying is inevitable, all the more so as greenhouse gas emissions are unabated, and potentially associated with considerable ecological and agricultural impacts [START_REF]WG2 Summary for Policymakers[END_REF]. Compared to LSRH, the boreal summer mid-latitude drying is even less constrained by the observed global warming. Yet, the combined effect of the HadCRUT5 and HadISDH observational constraints is still substantial as evidenced by a 30% reduction in the inter-model spread. In all cases, none of the available CMIP6 models is able to capture the magnitude of the early 21 st century drying. As discussed in the first KCC study [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF], this remark highlights one possible limitation common to all Bayesian statistical methods: the fact that the prior distribution may not provide a comprehensive sampling of uncertainty.

Yet, and unlike most alternative methods, the KCC method allows recent and future climate change to be constrained in a consistent manner. If we now focus on a subset of nine CMIP6 models that have provided at least three realizations of the historical simulations with individual radiative forcings (cf. Supplementary Materials), the double observational constraint provides an unequivocal demonstration that human emissions of greenhouse gases have been responsible for the recent continental near-surface drying (Fig. 2). This drying has been however countered by both anthropogenic and natural (volcanic) aerosols over much of the 20 th century, which may explain why it has not been observed before the 1990s. This finding is consistent with the latest IPCC report [START_REF] Douville | IPCC AR6 WG1[END_REF] and with previous evidence of a strong anthropogenic aerosol influence on the global water cycle across the 20th century [START_REF] Wu | Anthropogenic impact on Earth's hydrological cycle[END_REF][START_REF] Bonfils | Human influence on joint changes in temperature, rainfall and continental aridity[END_REF].

As most CMIP6 models used in Fig. 1, the selected subset used in Fig. 2 tends to overestimate the observed global warming, but to underestimate the observed continental drying (Fig. 2a).

Although the prior average drying varies depending on which model ensemble is chosen, their results are very close once the observational constraints are applied. Yet, the structural LSRH vs GSAT dependence across the nine CMIP6 models do not allow the KCC method to produce a scaled drying that is fully consistent with the observations. Looking at the best estimates, the posterior drying remains close to -0.5%/K while observations rather suggest a threefold -1.5%/K scaling (Fig. 2a). Similar results are found when all available CMIP6 models shown in Fig. 1 are considered (Fig. S5), thereby indicating that the underestimation of the scaled drying is not an artefact of the model subset used in Fig. 2. This mismatch highlights the need for large ensembles of independent models within CMIP and the importance of accounting for observational errors in KCC.

Compared to the prior ensemble mean response (dashed lines), the posterior response (solid lines) is not much different for non-GHG forcings, but suggests an overestimated sensitivity to the GHG forcing over the period covered by the HadISDH dataset (Fig. 2b). At first glance, this result is not consistent with the underestimated response found in Fig. 1 but would be consistent if Fig. 1 had been produced with the subset of only nine models available for detectionattribution (not shown). In the northern mid-latitudes (Fig. S6), the observed recent drying is underestimated, in agreement with the KCC results obtained across the 21 st century (Fig. S3d). Not surprisingly, the influence of the anthropogenic aerosols is more pronounced than at the global scale. Yet, the response to both natural and anthropogenic aerosols is not much influenced Template revised February 2021 by observations, that mostly constrain the response to the GHG forcing. This finding is consistent with a substantial narrowing of the CMIP6 projections using KCC given the increasing influence of GHGs across the 21 st century in both the intermediate and high-emission scenarios.

Discussion

Our results show that the continental drying projected by the CMIP6 models is highly modeldependent and that the multi-model ensemble mean is not necessarily the most likely response.

Constraining the projections with both HadCRUT5 and HadISDH observations allows us to narrow model uncertainties but also to predict a stronger than simulated decrease in RH, especially in the northern mid-latitudes. In contrast, the ensemble-mean continental drying conditional to the observed historical warming only is weaker than simulated and the associated inter-model spread is not much reduced compared to the prior distribution. This finding highlights a potential drawback of the current mitigation policy framings, based on global warming levels alone [START_REF] Hausfather | Climate simulations: recognize the 'hot model' problem[END_REF]. Other policy-relevant metrics could be considered, not only because of the uncertain hydrological consequences at a given global warming level [START_REF] Douville | IPCC AR6 WG1[END_REF], but also given the plausible vulnerability of the mitigation options compatible with the Paris Agreement to the forthcoming water cycle changes [START_REF] Douville | Water remains a blind spot in climate change policies[END_REF].

As an illustration, Fig. 3 provides a storyline of plausible climate changes associated with the top ten CMIP6 models projecting a strong continental drying across the 21 st century (Fig. S7). The SSP5-8.5 high-emission scenario is used again to illustrate this storyline approach. Compared to the full set of CMIP6 models (Fig. S8), the selected models show as expected a stronger recent drying, but which is also more realistic with respect to the best guess estimate derived from the HadISDH multi-member dataset (Fig. S5). Although they are thus more credible than the "hot models" [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF][START_REF] Hausfather | Climate simulations: recognize the 'hot model' problem[END_REF], they also project a stronger surface warming (Fig. 3b), especially over the continental areas that exhibit a stronger drying (Fig. 3a). Significant differences in the projected precipitation anomalies are found over both land and ocean surfaces, including a lesser increase in the boreal mid-latitudes as well as contrasted patterns in the tropics (Fig. 3c).

More surprisingly, total evapotranspiration anomalies show an overall weaker increase over land but a stronger increase over the ocean (Fig. 3d). This land-sea contrast is intriguing since it is sometimes assumed that surface evaporation will increase proportionally to the control evaporation under warmer conditions ( 6). Yet, possible changes in both soil moisture and plant stomatal conductance can invalidate this hypothesis over land. The weaker increase in land surface evapotranspiration despite a stronger near-surface drying (Fig. 3a) indeed suggests a key role for land surface processes. This hypothesis is consistent with a detailed analysis of the RH response to an abrupt quadrupling of atmospheric CO 2 in the CNRM-CM6-1 global climate model, where a direct CO 2 effect on plant transpiration was identified as a key factor for understanding a stronger continental drying than in the previous-generation CNRM-CM5 model [START_REF] Douville | Drivers of the enhanced decline of land surface relative humidity in CNRM-CM6-1[END_REF]. Yet, a key question that could not be investigated in this study based on atmosphere-only simulations was the extent to which the land surface processes can influence the projected global warming.

Revised Manuscript: Confidential Template revised February 2021

In the understanding of changes in climate sensitivity between successive generations of global climate models, land surface processes have been rarely mentioned as one of the possible causes [START_REF] Gettelman | High climate sensitivity in the Community Earth System Model Version 2 (CESM2)[END_REF]. Yet, our current understanding of climate variability at shorter timescales highlight that land surface anomalies can have remote and circumglobal effects [START_REF] Teng | Circumglobal Response to Prescribed Soil Moisture over North America[END_REF]. Idealized deforestation experiments also reveal non-local effects, including a global ocean cooling in the CMIP6 models where deforestation first leads to a widespread and significant cooling above the deforested areas [START_REF] Boysen | Global climate response to idealized deforestation in CMIP6 models[END_REF]. Such a possible land surface influence on global mean temperature is also consistent with Fig. 3b that shows an enhanced near-surface warming over the ocean in the strong-drying storyline. Such an interpretation remains however highly speculative given the strong overlap between the ten hottest and ten driest models within CMIP6 (Fig. S7). We just note here that both storylines are not mutually exclusive and that land surface feedbacks can contribute to amplify global warming.

Looking again at all available CMIP6 models (Fig. 4), there is a clear emergent relationship between near-surface drying and near-surface warming globally averaged over land (Fig. 4a). Such a relationship does not imply a causality, but is consistent with a possible amplification of global warming. Note that this emergent relationship would be much stronger if the four configurations of the ECMWF model (i.e., the only models that show an unrealistic increase in LSRH according to KCC) were discarded from the CMIP6 ensemble. Even without these four models, the linear regression fit would show a non-zero offset that is consistent with a fast adjustment (i.e., drying) induced by increased CO 2 levels [START_REF] Douville | Drivers of the enhanced decline of land surface relative humidity in CNRM-CM6-1[END_REF]. The possible outliers highlight the well-known weakness of emergent constraints that are not necessarily robust across successive generations of models [START_REF] Sanderson | The potential for structural errors in emergent constraints[END_REF]. In contrast, the KCC method is based on direct observational constraints and was shown to be relatively robust [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF][START_REF] Douville | Global warming at near-constant tropospheric relative humidity is supported by observations[END_REF]. This is further supported by the present study which show similar conclusions for CMIP6 (Fig. 2) and CMIP5 (Fig. S9) models.

An emergent constraint based on model biases in LSRH (Fig. S10a) would be also consistent with the KCC results regarding the possible underestimation of the projected drying by a significant fraction of the available CMIP6 models. Yet, the link between model climatologies and model responses is much weaker across CMIP5 models (Fig. S10b). This result agrees with the growing evidence that emergent constraints based on model biases may lead to overconfident projections [START_REF] Sanderson | The potential for structural errors in emergent constraints[END_REF]. It does not contradict our previous conclusions and our "strong-drying" storyline since KCC is based on projected anomalies relative to a reference climatology and thus does not explicitly consider model errors in mean climate. Yet, it may suggest that the calibration of coupled models should not only be based on radiative equilibrium at the top of atmosphere and resulting SSTs, but also on the water budget at the land surface and resulting soil moisture.

Coming back to Fig. 4, there is overall no clear link between LSRH and precipitation over land (Fig. 4b). In contrast, there is a strong emergent relationship between LSRH and land surface evapotranspiration (Fig. 4c). Since local evaporation does not strongly control local RH (because of moisture advection), Fig. 4c may be interpreted as the signature of a soil moisture and/or stomatal influence on both evapotranspiration and RH. Moreover, a stronger near-surface drying is generally associated with a strong ratio between annual runoff and annual precipitation (Fig. 4d), thereby further suggesting that changes in the land surface water budget may influence LSRH. Again, this emergent relationship would be stronger if the multiple configurations of the ECMWF model were not considered. These four models had been already noticed in Fig. 4a where they contribute to weaken the relationship between projected changes in RH vs Tair over land. Our results therefore support the claim that potential outliers may lead to spurious emergent constraints and overconfident projections across a limited set of models [START_REF] Sanderson | The potential for structural errors in emergent constraints[END_REF]. Yet, our more robust KCC method and the relationships shown in Fig. 4 also advocate for a stronger attention paid to the representation of land surface processes in global Earth System models [START_REF] Van Den Hurk | LS3MIP (v1.0) contribution to CMIP6: the Land Surface, Snow and Soil moisture Model Intercomparison Project -aims, setup and expected outcome[END_REF], and highlights the potential need for a stronger structural diversity therein.

MATERIAL AND METHODS

HadISDH and HadCRUT5 observations

HadISDH is a global gridded monthly mean surface humidity dataset. Quality controlled and homogenised / bias adjusted monthly mean anomalies (relative to a 1991-2020 base period) are provided alongside uncertainty estimates (observation and gridbox sampling). Actual values, climatological mean and standard deviation are also provided. The dataset begins in January 1973 and is updated annually. Here we use the 1973-2021 timeseries and focus over land only. HadISDH uses simultaneous subdaily temperature and dew point temperature data from >4500 quality controlled HadISD stations that have sufficiently long records. All humidity variables, including RH, are calculated at hourly resolution and monthly means are created. Monthly means are then homogenised using a pairwise algorithm to detect and adjust for features within the data that do not appear to be of climate origin. Stations with very large (> 15 % for RH) adjustments applied are removed. Measurement, climatological and homogeneity adjustment uncertainty is estimated for each month. Climatological averages over the 1991 to 2020 period are calculated and monthly mean climate anomalies obtained. These anomalies (in addition to climatological mean and standard deviation, actual values and uncertainty components) are then averaged over 5° by 5° gridboxes. Given the uneven distribution of stations over time and space, sampling uncertainty is estimated for each gridbox month. Although the coverage of land areas is not complete even at a 5° resolution (e.g., some limited gaps in the tropics), we consider that the whole timeseries are representative of the global land domain. This is indeed consistent with our previous finding [START_REF] Douville | Drivers of the enhanced decline of land surface relative humidity in CNRM-CM6-1[END_REF] that applying a HadISDH mask to the ERA5 reanalysis does not alter the annual mean timeseries compared to the full-coverage ERA5 global land values (cf. their Fig. 1). Regarding observational uncertainties, HadISDH does not directly provide multiple members as for the HadCRUT5 data sets. Yet, we have implemented a simple method to generate 200 estimates of the 1973-2021 annual mean timeseries. We have assumed no strong seasonality in the observational uncertainties so that the annual mean values of the monthly uncertainty estimates have been used to generate 200 members of annual mean values assuming a normal distribution. So doing, we keep the steep reduction in observational uncertainties across the 1973-2021 period, with much larger uncertainties in the early record.

HadCRUT5 is one of the main datasets used to monitor global and regional surface temperature variability and trends and is here only briefly described (more details can be found at https://www.metoffice.gov.uk/hadobs/hadcrut5/). It is a global surface temperature product that combines gridded land surface weather station air temperatures from CRUTEM5 with gridded ship and buoy SSTs from HadSST4. The observations have been quality controlled and bias adjusted to remove non-climatic artefacts. HadCRUT5 is a departure from its predecessors in that two different versions are offered. Both versions include monthly anomalies over 1850present on a 5° by 5° lat-lon grid. The first version has no interpolation of missing grid box Revised Manuscript: Confidential Template revised February 2021 values. As such there are wide coverage gaps in the polar regions and in the interiors of some continents like Africa and South America. The second version is a spatially complete, infilled or "analyzed" version with (almost) no coverage gaps and is the most suitable version to use for direct comparisons with climate model outputs (as in the present study). Both versions include a 200-member ensemble that samples the distribution of the systematic observational uncertainty as well as the analysis uncertainties.

CMIP6 and CMIP5 models

We make use of a large set of global climate models from both CMIP6 and CMIP5. We took all models providing at least one historical simulation and the corresponding SSP5-8.5 highemission scenario for CMIP6 (RCP8.5 for CMIP5) for both tas and hurs monthly mean variables, corresponding to near-surface air temperature and relative humidity respectively. As a result, we considered thirty-two CMIP6 models (slightly less for the intermediate SSP2-4.5 scenario, cf. missing models in Fig. S7

): ACCESS-CM2, ACCESS-ESM1-5, CESM2-WACCM, CMCC-CM2-SR5, CMCC-ESM2, CNRM-CM6-1-HR, CNRM-CM6-1, CNRM-ESM2-1, CanESM5, EC-Earth3-CC, EC-Earth3-Veg-LR, EC-Earth3-Veg, EC-Earth3, FGOALS-f3-L, FGOALS-g3, FIO-ESM-2-0, GFDL-CM4, GFDL-ESM4, GISS-E2-1-G, HadGEM3-GC31-LL, HadGEM3-GC31-MM, INM-CM4-8, INM-CM5-0, IPSL-CM6A-LR, KACE-1-0-G, MIROC- ES2L, MIROC6, MPI-ESM1-2-HR, MPI-ESM1-2-LR, MRI-ESM2-0, NorESM2-MM, UKESM1-0-LL.
Alternatively, we also used a subset of nine CMIP6 models for the detectionand-attribution analysis in Fig. 2 (ACCESS-ESM-1-5, CanESM5, CNRM-CM6-1, FGOALS-g3, GFDL-ESM4, HadGEM3-GC31-LL, IPSL-CM6A-LR, MIROC6, MRI-ESM2-0). These nine models were selected since they provide at least three realizations for the following four experiments: hist-ALL (historical simulations with all natural and anthropogenic forcings), hist-GHG (simulations driven by evolving GHG concentrations only), hist-AER (simulations driven by evolving anthropogenic aerosol loadings only), and hist-NAT (simulations driven by evolving solar and volcanic forcings only). Historical simulations and the corresponding RCP8.5 highemission scenario from only twenty CMIP5 models (bcc-csm1-1-m, BNU-ESM, CanESM2, CCSM4, CESM1-CAM5, CNRM-CM5, CSIRO-Mk3-6-0, GISS-E2-H, GISS-E2-H-CC, GISS-E2-R, GISS-E2-R-CC, HadGEM2-ES, inmcm4, MIROC5, MIROC-ESM, MIROC-ESM-CHEM, MRI-CGCM3, MRI-ESM1, NorESM1-M, NorESM1-ME) have been also used. All model outputs have been retrieved from the IPSL ESPRI platform (https://mesocentre.ipsl.fr/plate-forme-physique/). Annual mean anomalies have been estimated relative to the 1995-2014 baseline reference period, as in the latest IPCC report. A bilinear regridding to a common horizontal grid (T127 grid of the medium-resolution CNRM-CM6-1 model) has been used for the global mapping of the multi-model ensemble statistics (e.g., Fig. S1 and Fig 3), but the GSAT and global land average RH timeseries have been calculated from raw model outputs.

KCC method

The observational constraint method, called Kriging for Climate Change (KCC), has been previously applied to global and regional warming [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF][START_REF] Qasmi | Reducing uncertainty in local temperature projections[END_REF], as well as to global total precipitable water [START_REF] Douville | Global warming at near-constant tropospheric relative humidity is supported by observations[END_REF]. It can deal with at least two different variables [START_REF] Qasmi | Reducing uncertainty in local temperature projections[END_REF][START_REF] Douville | Global warming at near-constant tropospheric relative humidity is supported by observations[END_REF], and there is no need that related observations cover the same period of time. KCC consists of three steps. First, the forced response of each climate model is estimated over the whole 1850-2100 period (after concatenation of historical simulations with corresponding 21 st century projections). In order to Template revised February 2021 also get attribution statements, the responses to ALL (all forcings), NAT (natural forcings only) and GHG forcings are estimated separately. Second, the sample of the forced responses from available climate models is used as a prior of the real-world forced response, assuming that "models are statistically indistinguishable from the truth". Third, observations are used to derive a posterior distribution of the past and future forced response given observations. This Bayesian method can be summarised using the following equation:

y = Hx + ε,
where y is the timeseries of observations (a vector), x is the timeseries of the forced response (a vector), H is an observational operator (matrix), ε is the random noise associated with internal variability and measurement errors (a vector), and ε ∼ N(0, Σ y ), where N stands for the multivariate Gaussian distribution. In the case of two variables (as in the present study), x and y are simply built as the concatenation of the related timeseries. Climate models are used to construct a prior on x: π(x) = N(μ x , Σ x ). Then the posterior distribution given observations y can be derived as p(x|y) = N(μ p , Σ p ). Remarkably, μ p and Σ p are available in closed-form expressions.

In the following, we mostly assess the forced response of annual and global mean land surface relative humidity (LSRH), as well as the response to specific subsets of radiative forcings (attribution). These forced responses are then constrained by GMST (a good surrogate for GSAT) and/or by observations of LSRH. The same technique can be easily applied to selected regions (e.g., the northern mid-latitudes) and seasons (e.g., June to September). Therefore, we consider the following CMIP vector:

x = (T all , RH all , RH ghg , RH nat ), where each element is an entire 1850-2100 timeseries of the forced response, T and RH stand for GSAT and LSRH, respectively. "all", "ghg" or "nat" are the subsets of external forcings considered. Similarly, we define an observed vector as: y = (T obs , RH obs ), i.e., only observed timeseries are used in y. The length of these time-series varies: 1850-2021 for GMST, 1973-2021 only for LSRH. As a result, x is a very long vector, and all attribution or projection diagnoses presented below can be derived from the posterior distribution p(x|y). Following [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF], μ x and Σ x are estimated as the sample mean and covariance of the forced responses. Σ y requires statistical modelling of internal variability and measurement errors, and we use a mix of auto-regressive processes of order 1 (AR1) to model internal climate variability. The intrinsic variance of both GMST and LSRH is derived from observations after subtracting the multi-model mean estimate of the forced response of GSAT and LSRH respectively. We also assume a dependence between GMST and LSRH internal variability, by accounting for the correlation between the two residuals in Σ y . Such a dependence is however limited (correlation of 0. 1 but for a subset of twenty CMIP5 models in the RCP8.5 highemission scenario. The narrowing of model uncertainties is less pronounced than for the CMIP6 models in the SSP5-8.5 scenario, which may be due to the much smaller subset of independent models (only few modelling centers which provide multiple configurations of the same model). Yet, the effect of the KCC method on the ensemble mean response is very consistent between the CMIP5 and CMIP6 ensembles.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Constrained versus unconstrained annual forced changes in global-mean nearsurface air temperature (GSAT) and global-mean land near-surface relative humidity (LSRH). Mean (solid lines) and 5-95% range (shading) of the prior and posterior distributions of the annual forced GSAT and LSRH response to both natural and anthropogenic radiative forcings in historical simulations and SSP5-8.5 projections from 32 CMIP6 models: a) HadCRUT-only observational constraint on GSAT; b) HadCRUT-only observational constraint on LSRH; c) HadISDH-only observational constraint on LSRH; d) HadCRUT and HadISDH observational constraints on LSRH. After the constraint, the 5-95% interval at the end of the 21 st century is reduced by 33%, 17%, 22% and 37% in panel a, b, c, and d respectively. The mean observed anomalies are show as black (gray) filled circles when they are (not) used for constraining the model response.

Fig. 3 .Fig. 4 . 5 Fig. S2 .

 345S2 Fig. 3. Strong drying storylines. Difference in the ensemble mean anomalies (2081-2100 minus 1995-2014 climatology) between a subset of ten CMIP6 models showing the strongest globalmean continental drying in the SSP8-5.85 scenario and the full set of available CMIP6 models for the same scenario: a) near-surface relative humidity (%); b) near-surface temperature (K); c) total precipitation (mm/day); d) total surface evaporation (mm/day). Stippling denotes the regions where the difference is significant at the 5% level.

Fig. S3 .

 S3 Fig. S3. Same as in Fig. 1 but for LSRH averaged over the northern mid-latitude (30°N-60°N) land areas only. After the constraint, the 5-95% interval at the end of the 21 st century is reduced by 33%, 12%, 20% and 31% in panel a, b, c, and d respectively.

Fig. S4 .

 S4 Fig. S4. Same as in Fig. 1 but for LSRH averaged over the northern mid-latitude (30°N-60°N) land areas and from June to September only. After the constraint, the 5-95% interval at the end of the 21 st century is reduced by 33%, 11%, 18% and 30% in panel a, b, c, and d respectively.

Fig. S5 .

 S5 Fig.S5. Same as in Fig.2abut using all available CMIP6 models used in Fig.1. Note that observational uncertainties are slightly different between the two figures since the median forced model response is used to assess the observed internal variability in the KCC method.

Fig. S6 .

 S6 Fig. S6. Same as in Fig. 2 but for the northern mid-latitude land areas only.

Fig. S7 .

 S7 Fig. S7. Sorted GSAT (a) and LSRH (b) anomalies across thirty-two CMIP6 models. Anomalies are shown for two scenarios (SSP2-4.5 and SSP5-8.5) over the 2081-2100 period relative to the 1995-2014 baseline and are sorted based on the SSP5-8.5 high-emission scenario (in red). Models highlighted with a star symbol are those used in Fig. 2. The top 10 driest models in panel b are those used in Fig.3.

Fig. S8 .

 S8 Fig. S8. Ensemble mean anomalies (2081-2100 minus 1995-2014 climatology) from thirty-two CMIP6 models in the SSP8-5.85 high-emission scenario: a) near-surface relative humidity (%); b) near-surface temperature (K); c) total precipitation (mm/day); d) total surface evaporation (mm/day). This supplementary figure provides the background for assessing the strong drying storyline shown in Fig. 3.

Fig. S9 .

 S9 Fig.S9. Same as in Fig.1but for a subset of twenty CMIP5 models in the RCP8.5 highemission scenario. The narrowing of model uncertainties is less pronounced than for the CMIP6 models in the SSP5-8.5 scenario, which may be due to the much smaller subset of independent models (only few modelling centers which provide multiple configurations of the same model). Yet, the effect of the KCC method on the ensemble mean response is very consistent between the CMIP5 and CMIP6 ensembles.

Fig. S10 .

 S10 Fig. S10. Emergent constraint on LSRH projections. Scatterplots of LSRH anomalies as a function of the LSRH present-day (1995-2014) climatologies across all available a) CMIP6 models in the SSP-5.8.5 high-emission scenario, b) CMIP5 models in the corresponding RCP8.5 scenario. All anomalies are estimated as differences between the 2081-2100 and 1995-2014 climatologies. Each symbol represents one model. Blue/red bars denote +/-1 standard deviation of interannual variability. A linear regression fit is drawn as a black solid line (R 2 is the fraction of explained variance). The dashed vertical black line highlights the observed climatology according to the HadISDH dataset.

Acknowledgements:

The authors are grateful to all global modelling centers which have contributed to CMIP6, as well as to the people in charge of the ESGF archive (https://esgfnode.llnl.gov/projects/cmip6/) and of the ESPRI platform (https://mesocentre.ipsl.fr/plate-formephysique/). Thanks are also due to A. Ribes and S. Qasmi for their guidance regarding the implementation of the KCC method, and to A. Ribes and C. Delire for their useful comments about the preliminary manuscript. Funding: The study has been partly supported by the European Union's Horizon 2020 research and innovation program under the EUCP (Grant Agreement 776613) and the CONSTRAIN (Grant Agreement 820829) projects. Kate Willett was supported by the Met Office Hadley Centre Climate Programme funded by BEIS and Defra.

Data and materials availability: The HadCRUT5 observations are available at https://www.metoffice.gov.uk/hadobs/hadcrut5 / . The HadISDH observations are available at https://www.metoffice.gov.uk/hadobs/hadisdh / . All CMIP6 datasets are available at https://esgfnode.llnl.gov/search/cmip6 / . The full KCC package is available on gitlab at https://gitlab.com/saidqasmi/KCC/ and is also available under a GNU General Public License, version 3 (GPLv3) at https://doi.org/10.5281/zenodo.5233947. The more specific R scripts and data files needed to plot Fig. 1 and 2 with KCC are freely accessible from zenodo at: https://doi.org/10.5281/zenodo.7994731. Other codes for data curation, as well as drawing Fig. 3 and Fig. 4 are basic applications of CliMAF (https://climaf.readthedocs.io/en/master/ ) but can be also easily reproduced with your own favorite graphic tools. All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at https://science.org/...

Author contributions:

Data curation: HD, KW Conceptualization: HD Methodology: HD Calculation and vizualization: HD Writing: HD, KW

Competing interests: Authors declare that they have no competing interests.

Supplementary Materials for


Mohamad Mestoukirdi

Matteo Zecchin

David Gesbert

Qianrui Li

Nicolas Gresset

User-Centric Federated Learning

Keywords: Personalized federated learning, distributed optimization, user-centric aggregation, statistical learning theory

Data heterogeneity across participating devices poses one of the main challenges in federated learning as it has been shown to greatly hamper its convergence time and generalization capabilities. In this work, we address this limitation by enabling personalization using multiple user-centric aggregation rules at the parameter server. Our approach potentially produces a personalized model for each user at the cost of some extra downlink communication overhead. To strike a trade-off between personalization and communication efficiency, we propose a broadcast protocol that limits the number of personalized streams while retaining the essential advantages of our learning scheme. Through simulation results, our approach is shown to enjoy higher personalization capabilities, faster convergence, and better communication efficiency compared to other competing baseline solutions.

I. INTRODUCTION

Federated learning [START_REF] Mcmahan | Communication-efficient learning of deep networks from decentralized data[END_REF] has seen great success, being able to solve distributed learning problems in a communicationefficient and privacy-preserving manner. Specifically, federated learning provides to clients (e.g. smartphones, IoT devices, and organizations) the possibility of collaboratively train a model under the orchestration of a parameter server (PS) by iteratively aggregating locally optimized models and without off-loading local data [START_REF] Kairouz | Advances and open problems in federated learning[END_REF]. The original aggregation policy was implemented by Federated Averaging (FedAvg) [START_REF] Mcmahan | Communication-efficient learning of deep networks from decentralized data[END_REF], has been devised under the assumption that clients' local datasets are statistically identical, an assumption that is hardly met in practice. In fact, clients typically store datasets that are statistically heterogeneous and different in size [START_REF] Sattler | Clustered federated learning: Model-agnostic distributed multitask optimization under privacy constraints[END_REF], and are mainly interested in learning models that generalize well over their local data distribution through collaboration. Generally speaking, FedAvg exhibits slow convergence and poor generalization capabilities in such non-IID setting [START_REF] Li | Federated optimization in heterogeneous networks[END_REF]. To address these limitations, a large body of literature deals with personalization as a technique to reduce the detrimental effect of non-IID data. A straightforward solution consists in producing adapted models at a device scale by local fine-tuning procedures. Borrowing ideas from Model Agnostic Meta-Learning (MAML) [START_REF] Finn | Model-agnostic metalearning for fast adaptation of deep networks[END_REF], federated learning can be exploited in order to find a launch model that can be later personalized at each device using few gradient iterations † Equal contribution. M. Mestoukirdi, M. Zecchin and D. Gesbert are with the Communication Systems Department, EURECOM, Sophia-Antipolis, France. Emails: {mestouki, zecchin, gesbert}@eurecom.fr. Q.Li and N. Gresset are with Mitsubishi Electric R&D Centre Europe. Emails: {Q.Li, N.Gresset}@fr.merce.mee.com. The work of M. Zecchin is funded by the Marie Sklodowska Curie action WINDMILL (grant No. 813999). The work of M. Mestoukirdi is funded by Mitsubishi Electric R&D Centre Europe. [START_REF] Fallah | Personalized federated learning with theoretical guarantees: A model-agnostic metalearning approach[END_REF], [START_REF] Jiang | Improving federated learning personalization via model agnostic meta learning[END_REF]. Alternatively, local adaptation can be obtained by tuning only the last layer of a globally trained model [START_REF] Manoj Ghuhan Arivazhagan | Federated learning with personalization layers[END_REF] or by interpolating between a global model and locally trained ones [START_REF] Deng | Adaptive personalized federated learning[END_REF], [START_REF] Hanzely | Federated learning of a mixture of global and local models[END_REF]. However, these methods can fail at producing models with an acceptable generalization performance even for synthetic datasets [START_REF] Caldas | Leaf: A benchmark for federated settings[END_REF]. Adaptation can also be obtained leveraging user data similarity to personalize the training procedure. For instance, a Mixture of Experts formulation has been considered to learn a personalized mixing of the outputs of a commonly trained set of models [START_REF] Reisser | Federated mixture of experts[END_REF]. Similarly, [START_REF] Marfoq | Federated multi-task learning under a mixture of distributions[END_REF] proposed a distributed Expectation-Maximization (EM) algorithm concurrently converges to a set of shared hypotheses and a personalized linear combination of them at each device. Furthermore, [START_REF] Zhang | Personalized federated learning with first order model optimization[END_REF] proposed a personalized aggregation rule at the user side based on the validation accuracy of the locally trained models at the different devices. In order to be applicable, these techniques need to strike a good balance between communication overhead and the amount of personalization in the system. In fact, if on one hand, the expressiveness of the mixture is proportional to the number of mixed components; on the other, the communication load is linear in this quantity. Clustered Federated Learning (CFL) measures the similarity among the model updates during the optimization process in order to lump together users in homogeneous groups. For example, [START_REF] Briggs | Federated learning with hierarchical clustering of local updates to improve training on noniid data[END_REF], [START_REF] Sattler | Clustered federated learning: Model-agnostic distributed multitask optimization under privacy constraints[END_REF] proposed a hierarchical strategy in which the original set of users is gradually divided into smaller groups and, for each group, the federated learning algorithm is branched in a new decoupled optimization problem.

In this work, we propose a different approach to achieve personalization by allowing multiple user-centric aggregation strategies at the PS. The mixing strategies account for the existence of heterogeneous clients in the system and exploit estimates of the statistical similarity among clients that are obtained at the beginning of the federated learning procedure. Furthermore, the number of distinct aggregation rules -also termed personalized streams -can be fixed in order to strike a good trade-off between communication and learning efficiency.

In particular, the contributions of this work are:

1) We propose a user-centric aggregation rule to personalize users' local models. This rule exploits a novel similarity score that quantifies the discrepancy between individual user data distributions. Different from previous algorithms based on user clustering [START_REF] Sattler | Clustered federated learning: Model-agnostic distributed multitask optimization under privacy constraints[END_REF], [START_REF] Briggs | Federated learning with hierarchical clustering of local updates to improve training on noniid data[END_REF], our approach enables collaboration across all the nodes during training and, as a result, outperforms the above techniques, especially when clear clusters of users do not exist.

Conversely to [START_REF] Marfoq | Federated multi-task learning under a mixture of distributions[END_REF], personalization is performed at the PS and therefore without the additional cost of transmitting multiple models to each client. 2) Leveraging results from domain adaptation theory, we provide an upper bound on the risk w.r.t. the local data distribution of the personalized models obtained by our aggregation strategy. The result is used to obtain insights on how to determine the degree of collaboration among the devices. 3) We propose a heuristic strategy to compute the mixing coefficients for the personalized aggregation without accounting for the communication overhead. Then, in order to limit the communication burden introduced by the personalized aggregation, we propose to limit the number of personalized streams using the centroids obtained by clustering the mixing coefficient vectors. 4) We provide simulation results for different scenarios and demonstrate that our approach exhibits faster convergence, higher personalization capabilities, and communication efficiency compared to other popular baseline algorithms.

II. LEARNING WITH HETEROGENEOUS DATA SOURCES

In this section, we provide theoretical guarantees for learners that combine data from heterogeneous data distributions. The set-up mirrors the one of personalized federated learning and the results are instrumental to derive our user-centric aggregation rule. In the following, we limit our analysis to the discrepancy distance (4) but it can be readily extended to other divergences [START_REF] Mestoukirdi | User-centric federated learning: Trading off wireless resources for personalization[END_REF].

In the federated learning setting, the weighted combination of the empirical loss terms of the collaborating devices represents the customary training objective. Namely, in a distributed system with m nodes, each endowed with a dataset D i of n i IID samples from a local distribution P i , the goal is to find a predictor f : X → Ŷ from a hypothesis class F that minimizes

L(f, w) = m i=1 w i n i (x,y)∈Di (f (x), y) (1)
where : Ŷ × Y → R + is a loss function and w = (w 1 , . . . , w m ) is a weighting scheme. In case of identically distributed local datasets, the typical weighting vector is w =

1 i ni (n 1 , . . . , n m )
, the relative fraction of data points stored at each device. This particular choice minimizes the variance of the aggregated empirical risk, which is also an unbiased estimate of the local risk at each node in this scenario. However, in the case of heterogeneous local distributions, the minimizer of w-weighted risk may transfer poorly to certain devices whose target distribution differs from the mixture P w = m i=1 w i P i . Furthermore, it may not exists a single weighting strategy that yields a universal predictor with satisfactory performance for all participating devices. To address the above limitation of a universal model, personalized federated learning allows adapting the learned solution at each device. In order to better understand the potential benefits and drawbacks coming from the collaboration with statistically similar but not identical devices, let us consider the point of view of a generic node i that has the freedom of choosing the degree of collaboration with the other devices in the distributed system. Namely, identifying the degree of collaboration between node i and the rest of users by the weighting vector w i = (w i,1 , . . . , w i,m ) (where w i,j defines how much node i relies on data from user j), we define the personalized objective for user i

L(f, w i ) = m j=1 w i,j n j (x,y)∈Dj (f (x), y) (2)
and the resulting personalized model

f wi = arg min f ∈F L(f, w i ). (3)
We now seek an answer to: "What's the proper choice of w i in order to obtain a personalized model f wi that performs well on the target distribution P i ?". This question is deeply tied to the problem of domain adaptation, in which the goal is to successfully aggregate multiple data sources in order to produce a model that transfers positively to a different and possibly unknown target domain. In our context, the dataset D i is made of data points drawn from the target distribution P i and the other devices' datasets provide samples from the sources {P j } j =i . Leveraging results from domain adaptation theory [START_REF] Ben-David | A theory of learning from different domains[END_REF], we provide learning guarantees on the performance of the personalized model f wi to gauge the effect of collaboration that we later use to devise the weights for the user-centric aggregation rules.

In order to avoid negative transfer, it is crucial to upper bound the performance of the predictor w.r.t. to the target task. The discrepancy distance introduced in [START_REF] Mansour | Domain adaptation: Learning bounds and algorithms[END_REF] provides a measure of similarity between learning tasks that can be used to this end. For a functional class F : X → Ŷ and two distributions P, Q on X , the discrepancy distance is defined as

d F (P, Q) = sup f,f ∈F |E x∼P [ (f, f )] -E x∼Q [ (f, f )]| (4)
where we streamlined notation denoting f (x) by f . For bounded and symmetric loss functions that satisfy the triangular inequality, the previous quantity allows obtaining the following inequality

E (x,y)∼P [ (f, y)] ≤ E (x,y)∼Q [ (f, y)] + d F (P, Q) + λ where λ = inf f ∈F E (x,y)∼P [ (f, y)] + E (x,y)∼Q [ (f, y)] .
We can exploit the inequality to obtain the following risk guarantee for f wi w.r.t the true minimzer f * of the risk for the distribution P i .

Theorem 1. For a symmetric and B-bounded range loss function that satisfies the triangular inequality, w.p. 1 -δ Fig. 1: Personalized Federated Learning with user-centric aggregates at round t.

the predictor f wi satisfies

E (x,y)∼Pi [ ( f wi , y)] -E (x,y)∼Pi [ (f * , y)] ≤ B m j=1 w 2 i,j n j 2d i n i log e i n i d + log 2 δ + 2 m j=1 w i,j d F (P i , P j ) + 2λ
where d is the VC-dimension of the function space resulting from the composition of F and and

λ = arg min f ∈F E (x,y)∼Pi [ (f, y)] + E (x,y)∼P w i [ (f, y)] .
Sketch of Proof. The proof works by bounding the population risk of (3) w.r.t. the local measure P i and, subsequently, the estimation error of the weighted empirical risk minimizer. Full details are provided in [START_REF] Mestoukirdi | User-centric federated learning: Trading off wireless resources for personalization[END_REF].

The theorem highlights that a fruitful collaboration should strike a balance between the bias terms due to dissimilarity between the local distribution and the risk estimation gains provided by the data points of other nodes. Analytically minimizing the upper bounds seems an appealing solution; however, the divergence terms are difficult to compute, especially under the privacy constraints that federated learning imposes. For this reason, in the following, we consider a heuristic method based on the similarity of the readily available users' model updates to estimate the collaboration coefficients.

III. USER-CENTRIC AGGREGATION

For a suitable hypothesis class parametrized by θ ∈ R d , federated learning approaches use an iterative procedure to minimize the aggregate loss (1) with w = 1 i ni (n 1 , . . . , n m ). At each round t, the PS broadcasts the parameter vector θ t-1 and then combines the locally optimized models by the clients {θ t-1 i } m i=1 according to the following aggregation rule

θ t ← m i=1 n i m j=1 n j θ t-1 i .
As mentioned in Sec. II, this aggregation rule has two shortcomings: it does not take into account the data heterogeneity across users, and it is bounded to produce a single solution.

For this reason, we propose a user-centric model aggregation scheme that takes into account the data heterogeneity across the different nodes participating in training and aims at neutralizing the bias induced by a universal model. Our proposal generalizes the naïve aggregation of FedAvg, by assigning a unique set of mixing coefficients w i to each user i, and consequently, a user-specific model aggregation at the PS side. Namely, at the PS side, the following set of user-centric aggregation steps are performed

θ t i ← m j=1 w i,j θ t-1/2 j for i = 1, . . . , m (5)
where now, θ

t-1/2 j
is the locally optimized model at node j starting from θ t-1 j , and θ t i is the user-centric aggregated model for user i at communication round t.

As we elaborate next, the mixing coefficients are heuristically defined based on a distribution similarity metric and the dataset size ratios. These coefficients are calculated before the start of federated training. The similarity score we propose is designed to favor collaboration among similar users and takes into account the relative dataset sizes, as more intelligence can be harvested from clients with larger data availability. Using these user-centric aggregation rules, each node ends up with its own personalized model that yields better generalization for the local data distribution. It is worth noting that the usercentric aggregation rule does not produce a minimizer of the user-centric aggregate loss given by (2). At each round, the PS aggregates model updates computed starting from a different set of parameters. Nonetheless, we find it to be a good approximation of the true update since personalized models for similar data sources tend to propagate in a close neighborhood. The aggregation in [START_REF] Zhang | Personalized federated learning with first order model optimization[END_REF] capitalizes on the same intuition.

A. Computing the collaboration coefficients

Computing the discrepancy distance (4) can be challenging in high-dimension, especially under the communication and privacy constraints imposed by federated learning. For this reason, we propose to compute the mixing coefficient based on the relative dataset sizes and the distribution similarity metric given by

∆ i,j ( θ) = 1 n i (x,y)∈Di ∇ (f θ , y) - 1 n j (x,y)∈Dj ∇ (f θ , y) 2 ≈ E z∼Pi ∇ (f θ , y) -E z∼Pj ∇ (f θ , y) 2
where the quality of the approximation depends on the number of samples n i and n j . The mixing coefficients for user i are then set to the following normalized exponential function

w i,j = nj ni e -1 2σ i σ j ∆i,j ( θ) m j =1 n j ni e -1 2σ i σ j ∆ i,j ( θ) for j = 1, . . . , m. (6)
The mixture coefficients are calculated at the PS during a special round prior to federated training. During this round, the PS broadcasts a common model denoted θ to the users, which compute the full gradient on their local datasets. At the same time, each node i locally estimates the value σ 2 i partitioning the local data in K batches {D k i } K k=1 of size n k and computing

σ 2 i = 1 K K k=1 1 n k (x,y)∈D k i ∇ (f θ , y) - 1 n i (x,y)∈Di ∇ (f θ , y) 2 (7) where σ 2
i is an estimate of the gradient variance computed over local datasets D k i sampled from the same target distribution. Once all the necessary quantities are computed, they are uploaded to the PS, which proceeds to calculate the mixture coefficients and initiates the federated training using the custom aggregation scheme given by [START_REF] Finn | Model-agnostic metalearning for fast adaptation of deep networks[END_REF]. Note that the proposed heuristic embodies the intuition provided by Th. 1. In fact, in the case of homogeneous users, it falls back to the standard FedAvg aggregation rule, while in the case of node i has an infinite amount of data it degenerates to the local learning rule which is optimal in that case.

B. Reducing the communication load

A full-fledged personalization by the means of the usercentric aggregation rule ( 5) would introduce a m-fold increase in communication load during the downlink phase as the original broadcast transmission is replaced by unicast ones. Although from a learning perspective the user-centric learning scheme is beneficial, it is also possible to consider overall system performance from a learning-communication trade-off point of view. The intuition is that, for small discrepancies between the user data distributions, the same model transfer positively to statistically similar devices. In order to strike a suitable tradeoff between learning accuracy and communication overhead we hereby propose to adaptively limit the number of personalized downlink streams. In particular, for a number of personalized models m t , we run a k-means clustering scheme with k = m t over the set of collaboration vectors {w i } m i=1 and we select the centroids { ŵi } mt i=1 to implement the m t personalized streams. We then proceed to replace the unicast transmission with group broadcast ones, in which all users belonging to the same cluster c receive the same personalized model ŵc . Choosing the right value for the number of personalized streams is critical in order to save communication bandwidth but at the same time obtain satisfactory personalization capabilities. It can be experimentally shown that clustering quality indicators such as the Silhouette score over the user-centric weights can be used to guide the search for the suitable number of streams m t [START_REF] Mestoukirdi | User-centric federated learning: Trading off wireless resources for personalization[END_REF].

IV. EXPERIMENTS

We now provide a series of experiments to showcase the personalization capabilities and communication efficiency of the proposed algorithm.

A. Set-up

In our simulation we consider a handwritten character/digit recognition task using the EMNIST dataset [START_REF] Cohen | Emnist: Extending mnist to handwritten letters[END_REF] and an image classification task using the CIFAR-10 dataset [START_REF] Krizhevsky | Learning multiple layers of features from tiny images[END_REF]. Data heterogeneity is induced by splitting and transforming the dataset in a different fashion across the group of devices. In particular, we analyze three different scenarios:

• Character/digit recognition with user-dependent label shift in which 10k EMNIST data points are split across 20 users according to their labels. The label distribution follows a Dirichlet distribution with parameter 0.4, as in [START_REF] Marfoq | Federated multi-task learning under a mixture of distributions[END_REF], [START_REF] Wang | Tackling the objective inconsistency problem in heterogeneous federated optimization[END_REF]. • Character/digit recognition with user-dependent label shift and covariate shift in which 100k samples from the EMNIST dataset are partitioned across 100 users each with a different label distribution, as in the previous scenario. Additionally, users are clustered in 4 group and at each group images are rotated of {0 • , 90 • , 180 • , 270 • } respectively. • Image classification with user-dependent concept shift in which the CIFAR-10 dataset is distributed across 20 users which are grouped in 4 clusters, for each group we apply a different random label permutation. For each scenario, we aim at solving the task at hand by leveraging the distributed and heterogeneous datasets. We compare our algorithm against four different baselines: FedAvg, local learning, CFL [START_REF] Sattler | Clustered federated learning: Model-agnostic distributed multitask optimization under privacy constraints[END_REF] and FedFomo [START_REF] Zhang | Personalized federated learning with first order model optimization[END_REF]. In all scenarios and for all algorithms, we train a LeNet-5 convolutional neural network [START_REF] Lecun | Gradientbased learning applied to document recognition[END_REF] using a stochastic gradient descent optimizer with a fixed learning rate η = 0.1 and momentum β = 0.9.

B. Personalization performance

We now report the average accuracy over 5 trials attained by the different approaches. We also study the personalization performance of our algorithm when we restrain the overall number of personalized streams, namely the number of personalized models that are concurrently learned. In Fig. 2a we report the average validation accuracy in the EMNIST label shift scenario. We first notice that in the case of label shift, harvesting intelligence from the datasets of other users amounts to a large performance gain compared to the localized learning strategy. This indicates that data heterogeneity is moderate and collaboration is fruitful. Nonetheless, personalization can still provide gains compared to FedAvg. Our solution yields a validation accuracy which is increasing in the number of personalized streams. Allowing maximum personalization, namely a different model at each user, we obtain a 3% gain in the average accuracy compared to FedAvg. CFL is not able to transfer intelligence among different groups of users and attains performance similar to the FedAvg. This behavior showcases the importance of soft clustering compared to the hard one for the task at hand. We find that FedFOMO, despite excelling in the case of strong statistical heterogeneity, fails to harvest intelligence in the label shift scenario. In Fig. 2b we report the personalization performance for the second scenario. In this case, we also consider the oracle baseline, which corresponds to running 4 different FedAvg instances, one for each cluster of users, as if the 4 groups of users were known beforehand. Different from the previous scenario, the additional shift in the covariate space renders personalization necessary in order to attain satisfactory performance. In fact, the oracle training largely outperforms FedAvg. Furthermore, as expected, our algorithm matches the oracle final performance when the number of personalized streams is 4 or more. Also, CLF and FedFOMO are able to correctly identify the 4 clusters. However, the former exhibits slower convergence due to the hierarchical clustering over time while the latter plateaus to a lower average accuracy level. We turn now to the more challenging CIFAR-10 image classification task. In Fig. 2c we report the average accuracy of the proposed solution for a varying number of personalized streams, the baselines, and the oracle solution. As expected, the label permutation renders collaboration extremely detrimental as the different learning tasks are conflicting. As a result, local learning provides better accuracy than FedAvg. On the other hand, personalization can still leverage data among clusters and provide gains also in this case. Our algorithm matches the oracle performance for a suitable number of personalized streams. This scenario is particularly suitable for hard clustering, which isolates conflicting data distributions. As a result, CFL matches the proposed solution. FedFOMO promptly detects clusters and therefore quickly converges, but it attains lower average accuracy compared to the proposed solution.

The performance reported so far is averaged over users and therefore fails to capture the existence of outliers performing worse than average. In order to assess the fairness of the training procedure, in Table I we report the worst user performance in the federated system. The proposed approach produces models with the highest worst case in all three scenarios.

C. Communication Efficiency

Personalization comes at the cost of increased communication load in the downlink transmission from the PS to the federated user. In order to compare the algorithm convergence time, we parametrize the distributed system using two parameters. We define by ρ = T ul T dl the ratio between model transmission time in uplink (UL) and downlink (DL). Typical values of ρ in wireless communication systems are in the [START_REF] Kairouz | Advances and open problems in federated learning[END_REF][START_REF] Li | Federated optimization in heterogeneous networks[END_REF] range because of the larger transmitting power of the base station compared to the edge devices. Furthermore, to account for unreliable computing devices, we model the random computing time T i at each user i by a shifted exponential r.v. with a cumulative distribution function P [T i > t] = 1 -1(t ≥ T min ) 1 -e -µ(t-Tmin)

where T min representing the minimum possible computing time and 1/µ being the average additional delay due to random computation impairments. Therefore, for a population of m devices, we then have

T comp = E [max{T 1 , . . . , T m }] = T min + H m µ
where H m is the m-th harmonic number. To study the communication efficiency we consider the simulation scenario with the EMNIST dataset with label and covariate shift. In Fig. 3 we report the time evolution of the validation accuracy in 3 different systems. A wireless systems with slow UL ρ = 4 and unreliable nodes T min = T dl = 1 µ , a wireless system with fast uplink ρ = 2 and reliable nodes T min = T dl , 1 µ = 0 and a wired system ρ = 1 (symmetric UL and DL) with reliable nodes T min = T dl , 1 µ = 0. The increased DL cost is negligible for wireless systems with strongly asymmetric UL/DL rates and in these cases, the proposed approach largely outperforms the baselines. In the case of more balanced UL and DL transmission times ρ = [START_REF] Mcmahan | Communication-efficient learning of deep networks from decentralized data[END_REF][START_REF] Kairouz | Advances and open problems in federated learning[END_REF] and reliable nodes, it

Fig. 2 :

 2 Fig. 2: Evolution of the average validation accuracy in the three simulation scenarios.

TABLE I :

 I Worst user performance averaged over 5 experiments.

		Local FedAvg Oracle CFL [3] FedFOMO [14]	Proposed
	EMNIST label shift	58.8	68.9	-	70.3	70.0	73.2 for k = 20
	EMNIST covariate and label shift	56.0	67.5	77.4	76.1	73.6	76.4 for k = 4
	CIFAR concept shift	35.7	19.6	49.1	48.6	45.5	49.1 for k = 4

becomes instead necessary to properly choose the number of personalized streams in order to render the solution practical. Nonetheless, the proposed approach remains the best even in this case for k = 4. Note that FedFOMO incurs a large communication cost as personalized aggregation is performed at the client-side.

V. CONCLUSION

In this work, we presented a novel federated learning algorithm that exploits multiple user-centric aggregation rules to produce personalized models. The aggregation rules are based on user-specific mixture coefficients that can be computed during one communication round prior to federated training. Additionally, in order to limit the communication burden of personalization, we propose a simple strategy to effectively limit the number of personalized streams. We experimentally study the performance of the proposed solution across different tasks. Overall, our solution yields personalized models with higher testing accuracy while at the same time being more communication-efficient compared to the competing baselines.


Eunjeong Jeong
email: eunjeong.jeong@eurecom.fr

Matteo Zecchin
email: matteo.zecchin@eurecom.fr

Marios Kountouris
email: marios.kountouris@eurecom.fr

Asynchronous Decentralized Learning over Unreliable Wireless Networks

Keywords: asynchronous decentralized learning, over-theair computation, device-to-device communication

Decentralized learning enables edge users to collaboratively train models by exchanging information via deviceto-device communication, yet prior works have been limited to wireless networks with fixed topologies and reliable workers. In this work, we propose an asynchronous decentralized stochastic gradient descent (DSGD) algorithm, which is robust to the inherent computation and communication failures occurring at the wireless network edge. We theoretically analyze its performance and establish a non-asymptotic convergence guarantee. Experimental results corroborate our analysis, demonstrating the benefits of asynchronicity and outdated gradient information reuse in decentralized learning over unreliable wireless networks.

I. INTRODUCTION

Distributed learning algorithms empower devices in wireless networks to collaboratively optimize the model parameters by alternating between local optimization and communication phases. Leveraging the aggregated computational power available at the wireless network edge in a communication efficient [START_REF] Mcmahan | Communication-efficient learning of deep networks from decentralized data[END_REF] and privacy preserving manner [START_REF] Yan | Distributed autonomous online learning: Regrets and intrinsic privacy-preserving properties[END_REF], distributed learning is considered to be a key technology enabler for future intelligent networks. A promising paradigm, which enables collaborative learning among edge devices communicating in a peer-to-peer (server-less) manner, is decentralized learning [START_REF] Tsitsiklis | Distributed asynchronous deterministic and stochastic gradient optimization algorithms[END_REF]. Differently from federated learning, decentralized algorithms do not require a star topology with a central parameter server, thus being more flexible with respect to the underlying connectivity [START_REF] Koloskova | A unified theory of decentralized sgd with changing topology and local updates[END_REF]. This feature renders decentralized learning particularly appealing for future wireless networks with device-to-device communication. Several decentralized learning schemes over wireless networks have been proposed and analyzed [START_REF] Xing | Decentralized federated learning via sgd over wireless d2d networks[END_REF]- [START_REF] Shi | Over-the-air decentralized federated learning[END_REF], highlighting the key role of over-the-air computation (AirComp) [START_REF] Amiri | Machine learning at the wireless edge: Distributed stochastic gradient descent over-the-air[END_REF] for low-latency training at the edge. Prior works have mainly considered wireless networks of reliable workers communicating in a fixed topology throughout the entire training procedure. Nevertheless, these assumptions are hardly met in practical systems, in which communication links can be intermittent or blocked, and devices may become temporarily unavailable due to computation impairments or energy saving reasons. Asynchronous distributed training has been shown to mitigate the effect of stragglers (slow workers) [START_REF] Dutta | Slow and stale gradients can win the race[END_REF]- [START_REF] Adikari | Decentralized optimization with non-identical sampling in presence of stragglers[END_REF]. However, harnessing the potential benefits of asynchronism in decentralized learning over unreliable wireless networks remains elusive.

In this paper, we propose an asynchronous implementation of decentralized stochastic gradient descent (DSGD) as a means to address the inherent communication and computation impairments of heterogeneous wireless networks. In particular, we study decentralized learning over a wireless network with a random time-varying communication topology, comprising unreliable devices that can become stragglers at any point of the learning process. To account for communication impairments, we propose a consensus strategy based on time-varying mixing matrices determined by the instantaneous network state. At the same time, we design the learning rates at the edge devices in such a way so as to preserve the stationary point of the original network objective in spite of the devices' heterogeneous computational capabilities. Finally, we provide a non-asymptotic convergence guarantee for the proposed algorithm, demonstrating that decentralized learning is possible even when outdated information from slow devices is used to locally train the models. Experimental results confirm our analysis and show that reusing stale gradient information can speed up convergence of asynchronous DSGD.

II. SYSTEM MODEL

We consider a network consisting of m wireless edge devices, in which each node i is endowed with a local loss function f i : R d → R and local parameter estimate θ i ∈ R d . The network objective consists in minimizing the aggregate network loss subject to a consensus constraint minimize θ1,...,θm

f (θ 1 , . . . , θ m ) := 1 m m i=1 f i (θ i ) (1)
s.t. θ 1 = θ 2 = • • • = θ m .
This corresponds to the distributed empirical risk minimization problem whenever f i is a loss term over a local dataset. In the following, f (θ) denotes the network objective

f (θ 1 , . . . , θ m ) θ1=•••=θm=θ and θ = 1/m m i=1 θ i .
To solve (1), we consider a DSGD algorithm according to which devices alternate between a local optimization based on gradient information (computation phase) and a communication phase.

A. Computation model

To locally optimize the model estimate θ i , we assume that each device can query a stochastic oracle satisfying the following properties. Assumption 1. At each node i, the gradient oracle g i (θ) satisfies the following properties for all θ ∈ R d

• E[g i (θ)] = ∇ θ f i (θ) (unbiasedness) • E g i (θ) -∇ θ f i (θ) 2 ≤ σ 2 (bounded variance) • E g i (θ) ≤ G 2 (bounded magnitude).
We admit the existence of straggling nodes and that a random subset of devices can become inactive or postpone local optimization procedures, e.g., due to computation impairments or energy constraints. As a result, devices may join the communication phase and disseminate a model that has been updated using gradient information computed using previous model estimates, or a model that has not been updated at all from the previous iteration(s). Formally, at every optimization round t, the local update rule is

θ (t+ 1 2 ) i = θ (t) i , if device i is straggler at round t θ (t) i -η t i g i (θ (t-τi)
), otherwise where η t i is a local learning rate and the delay τ i ≥ 0 accounts for the staleness of the gradient information at device i.

B. Communication model

The channel between any pair of device i and j follows a Rayleigh fading model. At every communication iteration t, devices can exchange information according to a connectivity graph G (t) = (V, E (t) ), where V = {1, 2, . . . , m} indices the network nodes and (i, j) ∈ E (t) if devices i and j can communicate during round t. We consider symmetric communication links; therefore the communication graph is undirected. While the connectivity graph is assumed to remain fixed within the optimization iteration, it may vary across optimization iterations due to deep fading, blockage, and/or synchronization failures.

III. ASYNCHRONOUS DECENTRALIZED SGD

The proposed asynchronous DSGD procedure, which takes into account both computation and communication failures, is detailed in Algorithm 1.

At the beginning of each training iteration t, non straggling devices update the local estimate θ (t) i according to (2) using a potentially outdated gradient information. Subsequently, based on the current connectivity graph G (t) = (V, E (t) ), devices agree on a symmetric and doubly stochastic mixing matrix W (t) using a Metropolis-Hastings weighting scheme [START_REF] Xiao | Distributed average consensus with timevarying Metropolis weights[END_REF]. The weights are very simple to compute and are amenable for distributed implementation. In particular, each device requires only knowledge of the degrees of its neighbors to determine the weights on its adjacent edges.

After that, it follows a communication phase in which devices exchange the updated estimates and employ a gossip scheme based on W (t) . To leverage AirComp capabilities, Algorithm 1 Asynchronous Decentralized SGD

Input: θ (0) i = 0 ∈ R d Output: θ(T ) 1: for t in [0, T ] do 2:
for each non straggling devices do Determine matrix

W (t) based on G (t) 6: for s in [1, S t ] do 7:
if s ≡ 0 (mod 2) then 8:

# Broadcast phase

9:
for each device i scheduled in slot s do 10:

Device i transmits [START_REF] Ozfatura | Decentralized SGD with over-theair computation[END_REF] 11:

Each device j ∈ N (t) i receives [START_REF] Xing | Federated learning over wireless device-to-device networks: Algorithms and convergence analysis[END_REF] 12:

Each device j ∈ N (t) i estimates [START_REF] Shi | Over-the-air decentralized federated learning[END_REF] 13:

end for for each device do 24:

model consensus as in [START_REF] Amiri | Machine learning at the wireless edge: Distributed stochastic gradient descent over-the-air[END_REF] 25:

end for 26: end for devices employ analog transmission together with the scheduling scheme proposed in [START_REF] Xing | Federated learning over wireless device-to-device networks: Algorithms and convergence analysis[END_REF]. Accordingly, the communication phase is divided into multiple pairs of communication slots. Each pair consists of an AirComp slot and a broadcast slot as illustrated in Fig. 1. During the AirComp slot s, the star center i receives the superposition of the signals transmitted by its neighboring devices

N (t) (i) = {j ∈ V : (i, j) ∈ E (t) }.
In particular, each scheduled node j ∈ N (t) (i) transmits to the star center i

x (s,t) j = γ (s,t) i h (s,t) i,j w (t) i,j θ (t+ 1 2 ) j (3)
where h (s,t) i,j ∈ C d is the channel coefficient between user i and j during slot s, γ (s,t) i ∈ R is a power alignment coefficient, and w (t) i,j is the (i, j) entry of the mixing matrix W . The star center i receives the aggregated signal

y (s,t) i = j∈N (i) h (s,t) i,j x (s,t) j + z (s,t) i (4)
where z (s,t) i ∼ N (0, σ w 1 d ) is a noise vector, and estimates the aggregated model as

ŷ(s,t) i = y (s,t) i γ (s,t) i = j∈N (i) w (t) i,j θ (t+ 1 2 ) j + z (s,t) i γ (s,t) i . (5)
On the other hand, during a broadcast slot s, scheduled node i transmits using a power scaling factor α (s,t) i the signal

x (s,t) i = α (s,t) i θ (t+ 1 2 ) i (6)
and all neighboring devices j ∈ N (t) (i) receive

y (s,t) j = h (s,t) j,i x (s,t) i + z (s,t) j (7)
and estimate the updated model as

ŷ(s,t) j = w (t) j,i y (s,t) j α (s,t) i h (s,t) j,i = w (t) j,i θ (t+ 1 2 ) i + z (s,t) j √ α i h j,i . (8
θ (t+1) i = (1 -ζ)θ (t+ 1 2 ) i + ζ    m j=1 w (t) i,j θ (t+ 1 2 ) j + ñi (t)    (9)
where ñ(t) i ∼ N (0, σ(t) w,i 1 d ) is a noise vector term that accounts for the aggregation of noise components during AirComp and broadcast transmissions at device i during communication phase t.

IV. CONVERGENCE ANALYSIS

In this section, we study the effect of communication and computation failures on the asynchronous DGSD procedure and prove its convergence.

A. Effect of Communication Failures

Communication impairments amount for a random connectivity graph with an edge set that differs at each different optimization iteration. From an algorithmic perspective, random communication impairments result in DSGD with stochastic mixing matrices. A particular class of stochastic mixing matrices are those that satisfy the expected consensus property.

Definition 1 (Expected Consensus Rate [START_REF] Koloskova | A unified theory of decentralized sgd with changing topology and local updates[END_REF]). A random matrix W ∈ R m×m is said to satisfy the expected consensus with rate p if for any

X ∈ R d×m E W W X -X 2 F ≤ (1 -p) X -X 2 F
where X = X 11 T m and the expectation is w.r.t. the random matrix W . Lemma 1. If the event that the connectivity graph G (t) is connected at round t has a probability q > 0 and the Metropolis-Hastings weighting is used to generated the mixing W (t) , the expected consensus rate is satisfied with rate p = qδ > 0, with δ being the expected consensus rate in case of a connected topology.

Proof. See Appendix A.

If the expected consensus is satisfied, it is then possible to establish a convergent behavior for the estimates generated by the proposed algorithm.

Lemma 2 (Consensus inequality). Under Assumption 1, after T iterations, decentralized SGD with a constant learning rate η and consensus step size ζ satisfies

m i=1 θ (T ) -θ(T ) 2 ≤η 2 12mG 2 (pζ) 2 + ζ 2 p m i=1 σ 2 w,i where σ 2 w,i = max T t=0 E ñ(t) i 2 .
Proof. See Appendix B.

Overall, communication failures amount to a reduced expected consensus rate compared to the scenario with perfect communication. At the same time, dropping users that are delayed and are unable to synchronize and perform Air-Comp, renders the communication protocol more flexible. For instance, in Fig. 2, we consider a network of nine nodes organized according to different topologies and show the evolution of the average spectral gap of the mixing matrix with Metropolis-Hastings weights, whenever devices not satisfying a certain delay constraint are dropped. As expected, stricter delay requirements result in sparser effective communication graphs and mixing matrices with smaller spectral gaps.

B. Effect of Computation Failures

Random computation impairments make the group of devices that effectively update the model parameter vary over time. To account for this in the analysis, we introduce a virtual learning rate that is zero in case of failed computation. Namely, the learning rate at device i during computation round t becomes

η(t) i = 0, if i is straggler at round t η (t) i , otherwise where η (t) i
is a specified learning rate value in case of successful computation. Furthermore, to ensure that the procedure converges to stationary points of the network objective even when edge devices have different computing capabilities, the expected learning rates have to be equalized. In particular, if

E[η (t)
i ] = η, ∀i, we have that stationary points are maintained in expectation, namely

m i=1 E[η (t) i ]∇f i (θ) = 0 =⇒ ∇f (θ) = 0.
Finally, the existence of straggling devices introduces asynchronicity in the decentralized optimization procedure. In particular, a device i that fails at completing the gradient computation at a given optimization iteration is allowed to apply the result in a later one, without discarding the computation results. While we do not specify the delay distribution, we rather introduce the following assumption regarding the staleness of gradients.

Assumption 2. For all iteration t, there exists a constant γ ≤ 1 such that

E ∇f ( θ(t) ) - m i=1 ∇f i (θ (t-τi) i ) m 2 ≤ γE ∇f ( θ(t) ) 2 + L 2 m i=1 E θ (t) i -θ(t) 2 m .
The above assumption is similar to the one in [START_REF] Dutta | Slow and stale gradients can win the race[END_REF] with an additional consensus error term. Note that the value of γ is proportional to the staleness of the gradients and in case of perfect synchronization (γ = 0) the bound amounts to a standard consensus error term.

C. Convergence Guarantee

In this subsection, we demonstrate the convergence of the decentralized optimization procedure to a stationary point of the problem (1).

Theorem 1. Consider a network of unreliable communicating devices in which the expected consensus rate is satisfied with constant p and each device can be a straggler with probability ρ i < 1. If Assumptions 1 and 2 are satisfied, asynchronous DSGD with constant learning rate

η i = min j (1 - ρ j )/( √ 4LT (1 -ρ i ))
and consensus rate ζ = 1/T 3/8 satisfies the following stationary condition

1 T T t=1 ∇f ( θ(t) ) 2 ≤ 8 √ L(f ( θ(T ) ) -f * ) γ ρ min √ T + 3G 2 L T 1/4 p 2 γ + L 4T σ 2 mγ min j (1 -ρ j ) + m i=1 σ 2 w,i mγ 2L 2 γ pT 3/8 + 4L √ L mT 1/4 ρ min where γ = 1 -γ, ρ min = min j (1 -ρ j ) and f * = min θ∈R d f (θ). Proof. See Appendix C.
The above theorem establishes a vanishing bound on the stationarity of the returned solution, which involves quantities related to both communication and computation impairments. In particular, the constant of the slowest vanishing terms T -1/4 contains the term p related to random connectivity, as well as γ and ρ min due to stragglers.

V. NUMERICAL RESULTS

The effectiveness of the proposed asynchronous DSGD scheme is assessed using a network of m = 15 devices that collaboratively optimize the parameters of a convolutional neural network (CNN) for image classification with Fashion-MNIST. Gradients are calculated using batches of 16 data samples and the performance is evaluated using a test set of 500 images. We model the channel gain between each device pair as Rayleigh fading and we assume a shifted exponential computation time at each device, i.e., T comp = T min +Exp(µ) with T min = 0.25s and µ = 1. In Fig. 3, nodes communicate only when the channel is in favorable conditions, i.e., when the channel gain exceeds a certain minimum threshold h min . This allows to save energy; however, while higher threshold values result into lower average energy consumption, they also produce mixing matrices with smaller consensus rate, thus increasing the convergence time.

To study the effect of computation impairments, our proposed asynchronous learning algorithm is compared with: (i) synchronous DSGD, which waits for all devices to finish their computations; and (ii) synchronous DSGD with a delay barrier T max , which discards computation from users that violate the maximum computing time. Compared to the latter, our asynchronous procedure allows for slow devices to reuse stale gradient computations during later iterations. In Fig. 4, we plot the evolution of the test accuracy of the aforementioned algorithms under two different values of T max . For a moderate delay constraint T max = E[T comp ], asynchronous DSGD and synchronous DSGD with delay barrier perform similarly as the fraction of slow users is modest. Nonetheless, imposing a delay constraint and discarding slow devices greatly reduces the training time compared to the synchronous DSGD case. On the other hand, for a stringent delay requirement, T max = 4 5 E[T comp ], reusing stale gradients turns out to be beneficial and the proposed asynchronous DSGD attains higher accuracy faster compared to the synchronous DSGD with a delay barrier.

VI. CONCLUSION

In this work, we have proposed and analyzed an asynchronous implementation of DSGD, which enables decentralized optimization over realistic wireless networks with unreliable communication and heterogeneous devices in terms of computation capabilities. We have studied the effect of both communication and computation failures on the training performance and proved non-asymptotic convergence guarantees for the proposed algorithm. The main takeaway is that reusing outdated gradient information from slow devices is beneficial in asynchronous decentralized learning.

APPENDIX

A. Proof of Lemma 1

Define the event E (t) := {G (t) is connected} and its complementary event Ē(t) . Whenever the Metropolis-Hasting weights are obtained from a connected graph, the resulting mixing matrix W (t) has a consensus rate greater than zero. Therefore, there exists δ > 0 such that

E W (t) |E (t) W (t) X -X 2 F ≤ (1 -δ) W (t) X -X 2 F It follows that, for any X ∈ R d×m E W (t) W (t) X -X 2 F =qE W (t) |E (t) W (t) X -X 2 F + (1 -q)E W (t) | Ē(t) X -X 2 F ≤q(1 -δ) W (t) X -X 2 F + (1 -q) X -X 2 F
where we have lower bounded the consensus rate by zero in case of disconnected topologies. Grouping terms and having assumed q > 0, we obtain that the expected consensus is satisfied with rate (1 -qδ) > 0.

B. Proof of Lemma 2

Similarly to [START_REF] Xing | Federated learning over wireless device-to-device networks: Algorithms and convergence analysis[END_REF], [START_REF] Koloskova | Decentralized deep learning with arbitrary communication compression[END_REF] we establish the following recursive inequality m i=1 E θ (t) -θ(t) and then solving the recursion we obtain the final expression.

C. Proof of Theorem 1

We denote stale gradients by g i ( θ(t) i ) = g i (θ (t-τi) i

). According to the update rule, at each iteration t + 1, we have

E[f ( θt+1 )] = E f θt - 1 m m i=1

Fig. 1 .

 1 Fig. 1. An example of the timeline for one training iteration composed of alternate Broadcast and AirComp slots.

  ) At the end of the communication phase, each node i obtains the new estimate θ (t+1) i combining all received signals and using a consensus with step size ζ ∈ (0, 1]

Fig. 2 .

 2 Fig.2. Average spectral gap under different delay constraints for mesh, ring, and two-dimensional torus topologies with 9 nodes. Each link is associated to a completion time ∼ Exp(1) and is dropped if it exceeds the delay tolerance value.

h min = 2 Fig. 3 .

 23 Fig. 3. Test accuracy versus time under different channel gain thresholds. Smaller thresholds result in larger average consensus rates and therefore in faster convergence.

Fig. 4 .

 4 Fig. 4. Test accuracy for the asynchronous, synchronous with delay barrier, and synchronous schemes under two different values of Tmax.

  For an L-smooth objective function, we have

								Defining γ = (1 -γ), telescoping and taking expectations we
								obtain
	E[f ( θ(t+1) )] ≤ f ( θ(t) ) -+ L 2m 2 E m i=1 η(t) m 1 i g i ( i=1 m θ(t) i )) ∇f ( θ(t) ), E[η :=T1 2 :=T2 + L 2m 2 ζ 2 i=1 (t) i g i ( m E	θ(t) i ))] ñ(t)	1 T	T t=1	∇f ( θ(t) )	2	≤2 + f ( θ0 ) -f ( θT ) ηT γ 1 T T t=1 L 2 m + E θ σ 2 4ηγ mT (t) i -θ(t) 2 mγ i=1 + 1 T T t=1 m 2 Lζ 2 i=1 ηm 2 γ E ñ(t) i .
								Defining σ 2 w,i = max T t=0 E	ñ(t) i	2	and bounding the consen-
								sus term by Lemma 2, we obtain
								1 T	T t=1	∇f ( θ(t) )	2	≤2	f ( θ0 ) -f ( θT ) ηT γ
								+	L 2 mγ	η 2 12mG 2 (pζ) 2 + ζ	2 p	m i=1	σ 2 w,i
								+	σ 2 4ηγ mT	+	Lζ 2 ηm 2 γ	m i=1	σ 2 w,i .
							2
				∇f i (	θ(t) i )	.
	Summing T 1 and T 2 we obtain
	T 1 + T 2 ≤ -+	η 2 ηL 2 (1 -γ) ∇f ( θ(t) ) i=1 2m m θ (t) i -θ(t) 2 2	+	σ 2 8mT
	-	η 4m 2	m i=1	∇f i (	θ(t)

η(t) i g i ( θ(t) i ) + ζ ñ(t) i

where the expectation is w.r.t. the stochastic gradients, the communication noise Ξ (t) , and the computation and commu-nication failures at iteration t + 1. i )

2

.

The final result is obtained setting η = 1 √ 4LT and ζ = 1 T 3/8 .

* These authors contributed equally to this work. The work of E. Jeong is supported from a Huawei France-funded Chair on Future Wireless Networks. The work of M. Zecchin is funded by the Marie Sklodowska Curie action WINDMILL (grant No. 813999).