\begin{document}

\title{DR-DSGD: A Distributionally Robust Decentralized\ Learning Algorithm over Graphs}

\begin{abstract}
In this paper, we propose to solve a regularized distributionally robust learning problem in the decentralized setting, taking into account the data distribution shift. By adding a Kullback-Liebler regularization function to the robust min-max optimization problem, the learning problem can be reduced to a modified robust minimization problem and solved efficiently. Leveraging the newly formulated optimization problem, we propose a robust version of Decentralized Stochastic Gradient Descent (DSGD), coined Distributionally Robust Decentralized Stochastic Gradient Descent (DR-DSGD). Under some mild assumptions and provided that the regularization parameter is larger than one, we theoretically prove that DR-DSGD achieves a convergence rate of  $\mathcal{O}\left(1/\sqrt{KT} + K/T\right)$, where $K$ is the number of devices and $T$ is the number of iterations. Simulation results show that our proposed algorithm can improve the worst distribution test accuracy by up to $10\%$. Moreover, DR-DSGD is more communication-efficient than DSGD since it requires fewer communication rounds (up to $20$ times less) to achieve the same worst distribution test accuracy target. Furthermore, the conducted experiments reveal that DR-DSGD results in a fairer performance across devices in terms of test accuracy.
\end{abstract}

\section{Introduction}\label{introduction}
Federated learning (FL) is a learning framework that allows the training of a model across multiple devices under the orchestration of a parameter server (PS). Unlike the traditional way of training ML models, where the individual data of the devices are shared with the PS, FL hides the raw data, so it can significantly reduce the communication cost. When combined with a privacy-preservation mechanism, it further ensures privacy. Training using FedAvg presents several challenges that need to be tackled. It often fails to address the data heterogeneity issue. Though many enhanced variants of federated averaging (FedAvg) \citep{li2020federated,liang2019variance,karimireddy2020scaffold,karimireddy2020mime,wang2020tackling,Mitra2021,horvath2022fedshuffle} were proposed to solve this issue in the mean risk minimization setting, local data distributions might differ greatly from the average distribution. Therefore, a considerable drop in the global model performance on local data is seen, suggesting that the mean risk minimization may not be the right objective. Another major issue in FL is fairness. In many cases, the resultant learning models are biased or unfair in the sense that they discriminate against certain device groups \citep{hardt2016equality}. Finally, FL relies on the existence of a PS to collect and distribute model parameters, which is not always feasible or even accessible to devices that are located far away. 

Even though several FL algorithms  \citep{kairouz2019, yang2019, li2020} have been proposed for the distributed learning problem, FedAvg-style methods \citep{fedavg,li2020federated,wang2020tackling} remain the state-of-the-art algorithm. Specifically, FedAvg entails performing one or multiple local iterations at each device before communicating with the PS, which in turn performs periodic averaging. However, because FedAvg is based on the empirical risk minimization (ERM) to solve the distributed learning problem, i.e. FedAvg minimizes the empirical distribution of the local losses, its performance deteriorates when the local data are distributed non-identically across devices. While the ERM formulation assumes that all local data come from the same distribution, local data distributions might significantly diverge from the average distribution. As a result, even though the global model has a good average test accuracy, its performance locally drops when the local data are heterogeneous. In fact, increasing the diversity of local data distributions has been shown to reduce the generalization ability of the global model derived by solving the distributed learning problem using FedAvg \citep{Li2020On, sahu2018convergence, zhao2018federated}. 

While there are many definitions of robustness, the focus of our work is on distributionally robustness, that is, being robust to the distribution shift of the local data distributions. We consider a distributionally robust learning perspective by seeking the best solution for the worst-case distribution. Another key focus of this work is to investigate the fairness of the performance across the different devices participating in the learning. Fairness aims to reduce the difference in performance on the local datasets to ensure that the model performance is uniform across the devices participating in the learning process. In the FL context, achieving fair performance among devices is a critical challenge. In fact, existing techniques in FL, such as FedAvg \citep{fedavg} lead to non-uniform performance across the network, especially for large networks, since they favour or hurt the model performance on certain devices. While the average performance is high, these techniques do not ensure a uniform performance across devices. 

PS-based learning (star topology) incurs a significant bottleneck in terms of communication latency, scalability, bandwidth, and fault tolerance. Decentralized topologies circumvent these limitations and hence have significantly greater scalability to larger datasets and systems. In fact, while the communication cost increases with the number of devices in the PS-based topology, it is generally constant (in a ring or torus topology), or a slowly increasing function in the number of devices since decentralizing learning only requires on-device computation and local communication with neighboring devices without the need of a PS. Several works investigated the decentralizing learning problem \citep{yuan2016convergence,zeng2018nonconvex,wang2019matcha,wei2012distributed,shi2014linear,benissaid2021, duchi2011dual}; however, while interesting none of these works has considered solving the decentralized learning problem in a distributionally robust manner.
\paragraph{Summary of Contributions.}
The main contributions of this paper are summarized as follows
\begin{itemize}
    \item We propose a \textbf{distributionally robust learning algorithm}, dubbed Distributionally Robust Decentralized Stochastic Gradient Descent (DR-DSGD), that solves the learning problem in a decentralized manner while being robust to data distribution shift. To the best of our knowledge, our framework is the first to solve the regularized distributionally robust optimization problem in a decentralized topology\footnote{After acceptance of the paper, the work \citep{zecchin2022communication} was brought to our attention, where the authors propose to solve a decentralized distributionally robust optimization problem using a single loop gradient descent/ascent with compressed consensus scheme.}.
    \item Under some mild assumptions and provided that the regularization parameter is larger than one, we prove that DR-DSGD achieves a fast convergence rate of $\mathcal{O}\left(\frac{1}{\sqrt{KT}} + \frac{K}{T}\right)$, where $K$ is the number of devices and $T$ is the number of iterations, as shown in Corollary \ref{corollary}. Note that, unlike existing FL frameworks that rely on the \textbf{unbiasedness} of the stochastic gradients, our analysis is more challenging and different from the traditional analyses for decentralized SGD since it involves \textbf{biased} stochastic gradients stemming from the compositional nature of the reformulated loss function.
    \item We demonstrate the robustness of our approach compared to vanilla decentralized SGD via numerical simulations. It is shown that DR-DSGD leads to an improvement of up to $10\%$ in the worst distribution test accuracy while achieving a reduction of up to $20$ times less in terms of communication rounds. 
    \item Furthermore, we show by simulations that DR-DSGD leads to a fairer performance across devices in terms of test accuracy. In fact, our proposed algorithm reduces the variance of test accuracies across all devices by up to $60\%$ while maintaining the same average accuracy.
\end{itemize}
\paragraph{Paper Organization.}
The remainder of this paper is organized as follows. In Section \ref{probForm}, we briefly describe the problem formulation and show the difference between the ERM and distributionally robust optimization (DRO) formulation. Then, we present our proposed framework, \emph{DR-DSGD}, for solving the decentralized learning problem in a distributionally robust manner in Section \ref{propsol}. In Section \ref{convergence}, we prove the convergence of DR-DSGD theoretically under some mild conditions. Section \ref{eval} validates the performance of DR-DSGD by simulations and shows the robustness of our proposed approach compared to DSGD. Finally, the paper concludes with some final remarks in Section \ref{SecConc}. The details of the proofs of our results are deferred to the appendices.

\section{Related Works}
\paragraph{Robust Federated Learning.} 
Recent robust FL algorithms \citep{mohri2019agnostic,reisizadeh2020robust,deng2021distributionally, fedboost} have been proposed for the learning problem in the PS-based topology. Instead of minimizing the loss with respect to the average distribution among the data distributions from local clients, the authors in \citep{mohri2019agnostic} proposed agnostic federated learning (AFL), which optimizes the global model for a target distribution formed by any mixture of the devices' distributions. Specifically, AFL casts the FL problem into a min-max optimization problem and finds the worst loss over all possible convex combinations of the devices' distributions. \citet{reisizadeh2020robust} proposed FedRobust, a variant of local stochastic gradient descent ascent (SGDA), aiming to learn a model for the worst-case affine shift by assuming that a device's data distribution is an affine transformation of a global one. However, FedRobust requires each client to have enough data to estimate the local worst-case shift; otherwise, the global model performance on the worst distribution deteriorates. The authors in \citep{deng2021distributionally} proposed a distributionally robust federated averaging (DRFA) algorithm with reduced communication. Instead of using the ERM formulation, the authors adopt a DRO objective by formulating a distributed learning problem to minimize a distributionally robust empirical loss, while periodically averaging the local models as done in FedAvg \citep{fedavg}. Using the Bregman Divergence as the loss function, the authors in \citep{fedboost} proposed FedBoost, a communication-efficient FL algorithm based on learning the optimal mixture weights on an ensemble of pre-trained models by communicating only a subset of the models to any device. The work of \citep{pillutla2019robust} tackles the problem of robustness to corrupted updates by applying robust aggregation based on the geometric median. Robustness to data and model poisoning attacks was also investigated by \citep{li2021ditto}. DRO has been applied in multi-regional power systems to improve reliability by considering the uncertainty of wind power distributions and constructing a multi-objective function that maintains a trade-off between the operational cost and risk \citep{Lipeng2020, Hu2021}. However, none of these works provides any convergence guarantees. Furthermore, our analysis is different from these works since it involves biased stochastic gradients.
\paragraph{Fairness in Federated Learning.} Recently, there has been a growing interest in developing FL algorithms that guarantee fairness across devices \citep{mohri2019agnostic, qffl,term}. Inspired by works in fair resource allocation for wireless networks, the authors in \citep{qffl} proposed $q$-FFL, an FL algorithm that addresses fairness issues by minimizing an average reweighted loss parameterized by $q$. The proposed algorithm assigns larger weights to devices with higher losses to achieve a uniform performance across devices. Tilted empirical risk minimization (TERM), proposed in \citep{term}, has a similar goal as  $q$-FFL, i.e. to achieve fairer accuracy distributions among devices while ensuring similar average performance. 
\paragraph{Decentralized Learning.}
Decentralized optimization finds applications in various areas including wireless sensor networks \citep{mihaylov2009decentralized,avci2018wireless,soret2021learning}, networked multi-agent systems \citep{inalhan2002decentralized,ren2007information,johansson2008distributed}, as well as the area of smart grid implementations \citep{kekatos2012distributed}. Several popular algorithms based on gradient descent \citep{yuan2016convergence,zeng2018nonconvex,wang2019matcha}, alternating direction method of multipliers (ADMM) \citep{wei2012distributed,shi2014linear,benissaid2021}, or dual averaging \citep{duchi2011dual} have been proposed to tackle the decentralized learning problem. 
\section{Notations \& Problem Formulation}\label{probForm} 
\subsection{Notations}
Throughout the whole paper, we use bold font for vectors and matrices. The notation $\nabla f$ stands for the gradient of the function $f$, and $\mathbb{E}[\cdot]$ denotes the expectation operator. The symbols $\|\cdot\|$, and $\|\cdot\|_F$  denote the $\ell_2$-norm of a vector, and the Frobenius norm of a matrix, respectively. For a positive integer number $n$, we write $[n] \triangleq \{1, 2,\dots, n\}$. The set of vectors of size $K$ with all entries being positive is denoted by $\mathbb{R}_{+}^K$. The notations $\bm{0}$ and $\bm{1}$ denote a vector with all entries equal to zero, or one, respectively (its size is to be understood from the context). Furthermore, we define the matrices: $\bm{I}$ the identity matrix and $\bm{J} = \frac{1}{K} \bm{1} \bm{1}^T$. For a square matrix $\bm{A}$, $\text{Tr}(\bm{A})$ is the trace of $\bm{A}$, i.e. the sum of elements on the main diagonal. Finally, for the limiting behavior of functions, $f=\mathcal{O}(g)$ means that $f$ is bounded above up to a constant factor by $g$, asymptotically. 
\subsection{Problem Formulation}
We consider a connected network consisting of a set $\mathcal{V}$ of $K$ devices. Each device $i \in [K]$ has its data distribution $\mathcal{D}_i$ supported on domain $\Xi_i:= (\mathcal{X}_i,\mathcal{Y}_i)$. The connectivity among devices is represented as an undirected connected communication graph $\mathcal{G}$ having the set $\mathcal{E} \subseteq \mathcal{V} \times \mathcal{V}$ of edges, as illustrated in Fig. \ref{decentralized}. The set of neighbors of device $i$ is defined as $\mathcal{N}_i = \{ j | (i, j) \in \mathcal{E} \}$ whose cardinality is $|\mathcal{N}_i| = d_{i}$. Note that $(i, j) \in \mathcal{E}$ if and only if devices $i$ and $j$ are connected by a communication link; in other words, these devices can exchange information directly. All devices collaborate to solve the optimization problem given by
\begin{align}\label{eq0}
\min_{\bm{\Theta} \in  \mathbb{R}^{d}} \sum_{i=1}^K \frac{n_i}{n} f_i(\bm{\Theta}) \quad \text{where} \quad f_i(\bm{\Theta}) = \mathbb{E}_{\xi_i \sim \mathcal{D}_i}[\ell(\bm{\Theta};\xi_i)],
\end{align}
where $\xi_i:=(x_i,y_i)$ denotes the set of features $x_i$ and labels $y_i$ of device $i$. The function $\ell(\bm{\Theta}, \xi_i)$ is the cost of predicting $y_i$ from $x_i$, where $\bm{\Theta}$ denotes the model parameters, e.g., the weights/biases of a neural network. Here, $n_i$ denotes the number of training examples drawn from $\mathcal{D}_i$ and $n = \sum_{i=1}^K n_i$ is the total number of examples. 

Without loss of generality, we assume in the remainder that all devices have the same number of samples, and therefore $n_i/n = 1/K$, $\forall i \in [K]$. In this case, problem (\ref{eq0}) writes as
\begin{align}\label{eq00}
&\min_{\bm{\Theta} \in  \mathbb{R}^{d}} \frac{1}{K} \sum_{i=1}^K f_i(\bm{\Theta}).
\end{align}

\begin{figure}
\caption{Illustration of a graph topology ($K=8$) and the interactions between the devices ($D_1-D_8$).}
\label{decentralized}
\end{figure}

\begin{algorithm}[b]
\renewcommand{\algorithmiccomment}[1]{#1}
	\caption{\textsc{Vanilla Decentralized SGD (DSGD)}}
	\label{dsgd}
	\begin{algorithmic}[1]

		\FOR[{{\it in parallel for all devices $i \in [K]$}}]{$t$ \textbf{in} $0,\dots, T-1$}
		\STATE Sample a mini-batch of size $\{\xi_j^{t}\}_{j=1}^B$, compute gradient $\bm{g}_i(\bm{\theta}_i^{t}) := \frac{1}{B} \sum_{j=1}^B \nabla \ell(\bm{\theta}_i^{t}, \xi_j^{t})$
		\STATE $\bm{\theta}_i^{t + \frac{1}{2}} := \bm{\theta}_i^{t} - \eta \bm{g}_{i}(\bm{\theta}_i^{t})$
		\STATE Send $\bm{\theta}_i^{t + \frac{1}{2}}$ to neighbors
		\STATE $\bm{\theta}_i^{t + 1} := \sum_{j = 1}^{K} W_{ij} \bm{\theta}_j^{t + \frac{1}{2}}$
		\ENDFOR
	\end{algorithmic}
\end{algorithm}

One way to solve (\ref{eq0}) in a decentralized way is to use vanilla decentralized SGD (DSGD) \citep{yuan2016convergence}. As shown in Algorithm \ref{dsgd}, each device in DSGD performs two steps: (i) a local stochastic gradient update (Line 3) using the learning rate $\eta$, and (ii) a consensus operation in which it averages its model with its neighbors' models (Lines 4-5) using the weights of the connectivity (mixing) matrix of the network $\bm{W} = [W_{ij}] \in \mathbb{R}^{K\times K}$. The mixing matrix $\bm{W}$ is often assumed to be a symmetric ($\bm{W} =\bm{W}^T$) and doubly stochastic ($\bm{W} \bm{1} = \bm{1}, \bm{1}^T \bm{W}  =\bm{1}^T$) matrix, such that $W_{ij} \in [0,1]$, and if $(i,j) \notin \mathcal{E}$, then $W_{ij} = 0$. Assuming $\bm{W}$ to be symmetric and doubly stochastic is crucial
to ensure that the devices achieve consensus in terms of converging to the same stationary point.

While formulating problem (\ref{eq00}), we assume that the target distribution is given by
\begin{align}
\overline{\mathcal{D}} = \frac{1}{K} \sum_{i=1}^K \mathcal{D}_i.  
\end{align}
However, the heterogeneity of local data owned by the devices involved in the learning presents a significant challenge in the FL setting. In fact, models resulting from solving (\ref{eq00}) lack robustness to distribution shifts and are vulnerable to adversarial attacks \citep{bhagoji2019analyzing}.
This is mainly due to the fact that the target distribution may be significantly different from $\overline{\mathcal{D}}$ in practice. 
\section{Proposed Solution}\label{propsol}
Our aim is to learn a global model $\bm{\Theta}$ from heterogeneous data coming from these possibly non-identical data distributions of the devices in a decentralized manner. To account for heterogeneous data distribution across devices, the authors in \citep{mohri2019agnostic} proposed agnostic FL, where the target distribution is given by
\begin{align}
\mathcal{D}_{\bm{\lambda}} = \sum_{i=1}^K \lambda_i \mathcal{D}_i,
\end{align}
where the weighting vector $\bm{\lambda}$ belongs to the $K$-dimensional simplex, $\Delta  = \{ \bm{\lambda} = (\lambda_1, \dots, \lambda_K)^T \in \mathbb{R}^K_+ : \sum_{i=1}^K \lambda_i =1 \}$. Note that this target distribution is more general than $\overline{\mathcal{D}}$ and it reduces to $\overline{\mathcal{D}}$ when $\lambda_i = 1/K, \forall i \in [K]$.

Unlike $\overline{\mathcal{D}}$ which gives equal weight to all distributions $\{\mathcal{D}_i\}_{i=1}^K$ during training, $\mathcal{D}_{\bm{\lambda}}$ is rather a mixture of devices' distributions, where the
unknown mixture weight $\bm{\lambda}$ is learned during training and not known a priori. In this case, the  distributionally robust empirical loss problem is given by the following min-max optimization problem
\begin{align} \label{eq:3}
 \min_{\bm{\Theta} \in \mathbb{R}^d} \max_{\bm{\lambda} \in \Delta} \sum_{i=1}^K \lambda_i f_i(\bm{\Theta}).
\end{align}
Although several distributed algorithms \citep{mohri2019agnostic,reisizadeh2020robust,deng2021distributionally, fedboost} have been proposed for (\ref{eq:3}), solving it in a decentralized fashion (in the absence of a PS) is a challenging task. Interestingly, when introducing a regularization term in (\ref{eq:3}) and by appropriately choosing the regularization function, the min-max optimization problem can be reduced to a robust minimization problem that can be solved in a decentralized manner, as shown later on. Specifically, the regularized version of problem (\ref{eq:3}) can be written as follows
\begin{align} \label{eq:4}
 \min_{\bm{\Theta} \in \mathbb{R}^d} \max_{\bm{\lambda} \in \Delta} \sum_{i=1}^K \lambda_i f_i(\bm{\Theta}) - \mu \phi(\bm{\lambda},1/K),
\end{align}
where $\mu >0$ is a regularization parameter, and $\phi(\bm{\lambda},1/K)$
is a divergence measure between $\{\lambda_i\}_{i=1}^K$ and the uniform probability that assigns the same weight $1/K$ to every device's distribution. The function $\phi$ can be seen as a penalty that ensures that the weight $\lambda_i$ is not far away from $1/K$. Note that $\mu \rightarrow 0$ gives the distributionally robust problem (\ref{eq:3}) while when $\mu \rightarrow \infty$, we recover the standard empirical risk minimization problem (\ref{eq00}). 

\begin{algorithm}[t!]
\renewcommand{\algorithmiccomment}[1]{#1}
	\caption{\textsc{Distributionally Robust Decentralized SGD (DR-DSGD)}}
	\label{drdsgd}
	\begin{algorithmic}[1]

		\FOR[{{\it in parallel for all devices $i \in [K]$}}]{$t$ \textbf{in} $0,\dots, T-1$}
		\STATE Sample a mini-batch of size $\{\xi_j^{t}\}_{j=1}^B$, compute the gradient $\bm{g}_i(\bm{\theta}_i^{t}) := \frac{1}{B} \sum_{j=1}^B \nabla \ell(\bm{\theta}_i^{t}, \xi_j^{t})$ and \\ $h(\bm{\theta}_i^{t};\mu) :=  \exp\big(\frac{1}{B} \sum_{j=1}^B \ell(\bm{\theta}_i^t, \xi_j^t)/\mu\big)$
		\STATE $\bm{\theta}_i^{t + \frac{1}{2}} := \bm{\theta}_i^{t} - \eta \times \frac{h_i(\bm{\theta}_i^{t};\mu)}{\mu} \times \bm{g}_{i}(\bm{\theta}_i^{t})$
		\STATE Send $\bm{\theta}_i^{t + \frac{1}{2}}$ to neighbors
		\STATE $\bm{\theta}_i^{t + 1} := \sum_{j = 1}^{K} W_{ij} \bm{\theta}_j^{t + \frac{1}{2}}$
		\ENDFOR
	\end{algorithmic}
\end{algorithm}

Although different choices of the $\phi$-divergence can be considered in (\ref{eq:4}), the robust optimization community has been particularly interested in the Kullback–Leibler (KL) divergence owing to its simplified formulation \citep{esfahani2018data}. In fact, when we consider the KL divergence, i.e. $\phi(\bm{\lambda},1/K)= \sum_{i=1}^K \lambda_i\log(\lambda_i K)$, then by exactly maximizing over $\bm{\lambda} \in \Delta$, the min-max problem, given in (\ref{eq:4}), is shown to be equivalent to \citep[Lemma 1]{huang2021compositional}
\begin{align} \label{eq:5}
\min_{\bm{\Theta} \in \mathbb{R}^d} \mu\log\left( \frac{1}{K}\sum_{i=1}^K \exp\left(f_i(\bm{\Theta})/\mu\right)\right).
\end{align}
Since $\log(\cdot)$ is a monotonically increasing function, then instead of solving (\ref{eq:5}), we simply solve the following problem
\begin{align} \label{eq:6}
 \min_{\bm{\Theta} \in \mathbb{R}^d} F(\bm{\Theta}) \triangleq \frac{1}{K} \sum_{i=1}^K F_i(\bm{\Theta}),
\end{align}
where $F_i(\bm{\Theta}) = \exp\big(f_i(\bm{\Theta})/\mu\big), \forall i \in [K]$. 
\begin{remark}
    Note that problem (\ref{eq:4}) is equivalent to problem (\ref{eq:3}) as $\mu \rightarrow 0$. The proximity of the solutions of (\ref{eq:4}) to (\ref{eq:3}) in the \textbf{convex} case is discussed in Appendix \ref{Appendix:prox}. In the remainder of this work, we focus on solving  the regularized distributional robust optimization problem \citep{mohri2019agnostic, deng2021distributionally} given in (\ref{eq:4}).
\end{remark}
Although any decentralized learning algorithm can be used to solve (\ref{eq:6}), the focus of this paper is to propose a distributionally robust implementation of DSGD. Our framework, coined Distributionally Robust Decentralized SGD (DR-DSGD), follows similar steps as DSGD with the main difference in the local update step
\begin{align}\label{update}
\bm{\theta}_i^{t + 1} = \sum_{i = 1}^{K} W_{ij} \left( \bm{\theta}_i^{t} - \frac{\eta}{\mu}  h(\bm{\theta}_i^{t};\mu) \bm{g}_{i}(\bm{\theta}_i^{t})\right),
\end{align}
where $h(\bm{\theta}_i^{t};\mu) :=  \exp\big(\frac{1}{B} \sum_{j=1}^B \ell(\bm{\theta}_i^t, \xi_j^t)/\mu\big)$ and $B$ is the mini-batch size $\{\xi_j^t\}_{j=1}^B$. Introducing the term $h(\bm{\theta}_i^t;\mu)/\mu$ in line 3 of Algorithm \ref{drdsgd} makes the algorithm more robust to the heterogeneous setting and ensures fairness across devices, as will be shown in the numerical simulations section.
\section{Convergence Analysis}\label{convergence}
This section provides a theoretical analysis of the convergence rate of the DR-DSGD algorithm. Before stating the main results of the paper, we make the following assumptions. \\
\textbf{Assumption 1. (Smoothness)} There exist constants $L_F$ and $L_1$, such that $\forall \bm{\theta_1}, \bm{\theta_2} \in \mathbb{R}^d$, we have
\begin{align}
&\|\nabla F(\bm{\theta_1}) - \nabla F(\bm{\theta_2}) \| \leq L_F \|\bm{\theta_1}-\bm{\theta_2}\|, \\
&\mathbb{E}[\|\bm{g}_i(\bm{\theta_1}) - \bm{g}_i(\bm{\theta_2}) \|] \leq L_1 \|\bm{\theta_1}-\bm{\theta_2}\|. \label{ff1} 
\end{align}

\textbf{Assumption 2. (Gradient Boundedness)}
The gradient of $f_i(\cdot)$ is bounded, i.e. there exits $G_1$ such that $\forall \bm{\theta} \in \mathbb{R}^d$, we have
\begin{align}
    & \|\nabla f_i(\bm{\theta})\| \leq G_1. \label{as000}
\end{align}
\textbf{Assumption 3. (Variance Boundedness)}
The variances of stochastic gradient $g_i(\cdot)$ and the function $l_i(\cdot,\xi)$ are bounded, i.e. there exit positive scalars $\sigma_1$ and $\sigma_2$ such that $\forall \bm{\theta} \in \mathbb{R}^d$, we have
\begin{align}
    & \mathbb{E}\left[|\ell(\bm{\theta}, \xi_i)-f_i(\bm{\theta})|^2\right] \leq \sigma_1^2, \\
    & \mathbb{E}\left[\|\bm{g}_i(\bm{\theta})-\nabla f_i(\bm{\theta})\|^2\right] \leq \sigma_2^2.
\end{align}
\textbf{Assumption 4. (Function Boundedness)} The function $F(\cdot)$ is lower bounded, i.e. $F_{inf} = \underset{\bm{\theta} \in \mathbb{R}^d}{\inf} F(\bm{\theta})$ such that $F_{inf} > - \infty$ and the functions $\ell(\cdot, \xi_i)$ are bounded, i.e. there exists $M>0$ such that $\forall \bm{\theta} \in \mathbb{R}^d$, we have
\begin{align}\label{as1}
    \mathbb{E}\left[|\ell(\theta, \xi_i)|\right] \leq M.
\end{align}
\textbf{Assumption 5. (Spectral Norm)} The spectral norm, defined as
$\rho = \|\mathbb{E}\left[\bm{W}^T\bm{W}\right]-\bm{J}\|$, is assumed to be less than 1.

Assumptions 1-5 are key assumptions that are often used in the context of distributed and compositional optimization \citep{ wang2016accelerating, wang2017stochastic, Li2020On, huang2021compositional}. Note that (\ref{as1}) can be fulfilled by imposing loss clipping \citep{xu2006robust,wu2007robust,yang2010relaxed} to most commonly-used loss functions. Furthermore, although the categorical cross-entropy function is not bounded upwards, it will only take on large values if the predictions are wrong. In fact, under some mild assumptions, the categorical cross-entropy will be around the entropy of a $M$-class uniform distribution, which is $\log(M)$ (refer to Appendix \ref{Appendix:bound} for a more in-depth discussion). Considering a relatively moderate number of classes, e.g., $M=100$, we get $\log(M) < 5$. Thus, the cross-entropy will generally be relatively small in practice, and the assumption will hold in practice.
\begin{remark}
It is worth mentioning that since the exponential function is convex, and non-decreasing function, then $\{\exp(f_i(\cdot))\}_{i=1}^K$ is convex when $\{f_i(\cdot)\}_{i=1}^K$ is convex. In this case, the convergence of the proposed method follows directly from existing results in the literature \citep{yuan2016convergence}.
\end{remark}
In the remainder of this section, we focus on the non-convex setting and we start by introducing the following matrices 
\begin{align}
    \bm{\theta}^{t} &= \left[\bm{\theta}_1^{t},\dots,\bm{\theta}_K^{t}\right],\\
    \nabla \bm{F}^{t} &= \left[\nabla \bm{F}_1(\bm{\theta}_1^{t}),\dots,\nabla \bm{F}_K(\bm{\theta}_K^{t})\right], \\
    \bm{U}^{t} &= \frac{1}{\mu}\left[ h(\bm{\theta}_1^{t};\mu) \bm{g}_1(\bm{\theta}_1^t),\dots,  h(\bm{\theta}_K^{t};\mu) \bm{g}_K(\bm{\theta}_K^t)\right].
\end{align}
\begin{remark}
In the non-convex case, the convergence analysis is more challenging than in the case of DSGD. In fact, due to the compositional nature of the local loss functions, the stochastic gradients are biased, i.e.
\begin{align}
\mathbb{E}\left[\exp\big(\ell(\bm{\theta}_i^t, \xi_i^t)/\mu\big) \bm{g}_{i}(\bm{\theta}_i^{t})\right] \neq \exp\big(f_i(\bm{\theta}_i^t)/\mu\big) \nabla f_{i}(\bm{\theta}_i^{t}).  
\end{align}
\end{remark}
To proceed with the analysis, we start by writing the matrix form of the update rule (\ref{update}) as
\begin{align}\label{newupdate}
    \bm{\theta}^{t+1} = \left(\bm{\theta}^{t} - \eta \bm{U}^{t}\right)\bm{W}.
\end{align}
Multiplying both sides of the update rule (\ref{newupdate}) by $\bm{1}/K$, we get
\begin{align}
    \bar{\bm{\theta}}^{t+1} = \bar{\bm{\theta}}^{t} - \frac{\eta}{ K}\bm{U}^{t}\bm{1},
\end{align}
where $\bar{\bm{\theta}}^{t}$ is the averaged iterate across devices defined as $\bar{\bm{\theta}}^{t} = \frac{1}{K} \sum_{i=1}^K \bm{\theta}_i^t$. Now, we are in a position to introduce our first Lemma.
\begin{lemma}\label{lemma1}
For any matrix $\bm{A} \in \mathbb{R}^{d\times K}$, we have
\begin{align}
    \mathbb{E}\left[\|\bm{A} \left( \bm{W}^n- \bm{J}\right)\|_F^2\right] \leq \rho^n \|\bm{A}\|_F^2.
\end{align}
\end{lemma}
\begin{proof}
The details of the proof can be found in Appendix \ref{proofLemma1}. 
\end{proof}
The following lemma provides some of the inequalities that we need for the proof of the main result.
\begin{lemma}\label{lemma3}
From Assumption 4, we have that
\begin{itemize}
    \item[(a)] There exists a constant $L_2(\mu)$ such that $\forall y_1, y_2 \in \mathcal{Y}$, we have
    \begin{align}
        \mathbb{E}[|\exp\left(y_1\right)\!-\! \exp\left(y_2\right) |] \leq L_2(\mu) |y_1-y_2|,
    \end{align}
    where $\mathcal{Y} \triangleq \{y = \ell(\bm{\theta}, \xi_i)/\mu \text{ such that } \bm{\theta} \in \mathbb{R}^d\}$ is the range of functions $\{\ell(\bm{\theta}, \xi_i)/\mu\}$.
    \item[(b)] There exists a constant $G_2(\mu)$ such that $\forall \theta \in \mathbb{R}^d$, we have
    \begin{align}
        \mathbb{E}[|\exp(\ell(\bm{\theta}, \xi_i)/\mu)|] \leq G_2(\mu). \label{ff0}
    \end{align}
    \item[(c)] There exists a constant $\sigma_3(\mu)$ such that $\forall \theta \in \mathbb{R}^d$, we have
    \begin{align}
        \mathbb{E}\left[|\exp(\ell(\bm{\theta}, \xi_i)/\mu)-\exp(f_i(\bm{\theta})/\mu)|^2\right] \leq (\sigma_3(\mu))^2. \label{ff3}
    \end{align}
\end{itemize}
\end{lemma}
\begin{proof}
The details of the proof can be found in Appendix \ref{proofLemma3}. 
\end{proof}
Using Lemmas \ref{lemma1} and \ref{lemma3}, we can now state Lemma \ref{lemma2}, which gives an upper bound on the discrepancies among the local models. 
\begin{lemma}\label{lemma2}
Let $\eta$ satisfy $\eta L_{\mu} < \frac{\mu (1-\sqrt{\rho})}{4 G_{\mu} \sqrt{\rho}}$. Provided that all local models are initiated at the same point, the discrepancies among the local models $\mathbb{E}\left[\|\bm{\theta}^{t}(\bm{I}\!-\!\bm{J})\|_F^2\right]$ can be upper bounded by 
\begin{align}
\frac{1}{K T} \sum_{t=1}^T \mathbb{E}\left[\|\bm{\theta}^t(\bm{I}-\bm{J})\|_F^2\right]\leq \frac{2 \eta^2 \rho G_{\mu}^2 [8\sigma_{\mu}^2(L_{\mu}^2+1)+ G_{\mu}^2]}{\mu^2 (1-\gamma_{\mu}) (1-\sqrt{\rho})^2},
\end{align}
where $\sigma_{\mu}=\max\{\sigma_1, \sigma_2,\sigma_3(\mu)\}$, $G_{\mu}=\max\{G_1,G_2(\mu)\}$, $L_{\mu} = \max\{L_F, L_1,L_2(\mu)\}$, and $\gamma_{\mu} = \frac{16 \eta^2 \rho L_{\mu}^2 G_{\mu}^2}{\mu^2 (1-\sqrt{\rho})^2}$.
\end{lemma}
\begin{proof}
The proof is deferred to Appendix \ref{proofLemma2}. In a nutshell, the proof uses the update rule of DR-DSGD, given in (\ref{newupdate}), the special property of the mixing matrix, and Lemma \ref{lemma1}.
\end{proof}
Next, we present the main theorem that states the 
convergence of our proposed algorithm.
\begin{theorem}\label{thm1} 
Let $\eta$ satisfy $\eta L_{\mu} < \min\{\frac{\mu (1-\sqrt{\rho})}{8 G_{\mu} \sqrt{\rho}},1\}$ and provided that all local models are initiated at the same point, then the averaged gradient norm is upper bounded as follows
\begin{align}\label{con_rate}
      \frac{1}{T} \sum_{t=1}^T \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right]\leq \frac{2 (F(\bar{\bm{\theta}}^{1})-F_{inf})}{\eta T}+ \frac{2 G_{\mu}^2 \sigma_{\mu}^2(L_{\mu}^2 + \mu^2)}{\mu^4 B}+\frac{16 \rho \eta^2 L_{\mu}^2  G_{\mu}^2 (G_{\mu}^2+\mu^2) [8\sigma_{\mu}^2(L_{\mu}^2+1)+ G_{\mu}^2]}{3 \mu^8 (1-\sqrt{\rho})^2}.
\end{align}
\end{theorem}
\begin{proof}
The proof of Theorem \ref{thm1} is detailed in Appendix \ref{proofthm1}. 
\end{proof}
The first term of the right-hand side of (\ref{con_rate}) is neither affected by the graph parameters nor by the regularization parameter $\mu$ and it exhibits a linear speedup in $T$. This term exists also in the convergence analysis of DSGD \citep{wang2019matcha,lian2017can}. The graph topology affects the third term via the value of the spectral norm $\rho$. In fact, a smaller value of $\rho$ makes the third term smaller. Finally, the regularization parameter $\mu$ has an impact on both the second and third terms. A  weakness of the derived bound is that as $\mu$ gets smaller, both the second and third terms grow larger. In the limit case, when $\mu \rightarrow 0$, both terms grow to infinity. In the next corollary, we limit ourselves to the case when $\mu \geq 1$ and show that under this setting DR-DSGD is guaranteed to converge.
\begin{remark}
Our analysis is still valid in the case where the mixing matrix changes at every iteration. Similar theoretical guarantees hold provided that the matrices $\{\bm{W}^t\}_{t=1}^T$ are independent and identically distributed and their spectral norm $\rho^t < 1$.
\end{remark}
When $\mu \geq 1$, it can be seen from the proof of Lemma \ref{lemma3}, that we can find constants, $L_2$, $G_2$, and $\sigma_3$, that do not depend on $\mu$. In this case, (\ref{con_rate}) can be written as
\begin{align}
    \frac{1}{T} \sum_{t=1}^T \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right]\leq \frac{2 (F(\bar{\bm{\theta}}^{1})-F_{inf})}{\eta T}+ \frac{2 G^2 \sigma^2(L^2 + 1)}{B}+\frac{16 \rho \eta^2 L^2  G^2 (G^2+1) [8\sigma^2(L^2+1)+ G_2]}{3 (1-\sqrt{\rho})^2}.
\end{align}
  Furthermore, if the learning rate $\eta$ and the  mini-batch size $B$ are chosen properly, we obtain the following corollary.
\begin{corollary}\label{corollary}
If $\mu \geq 1$ and if we choose $\eta =  \frac{1}{2L + \sqrt{T/K}}$ and $B = \sqrt{KT}$, then we have
\begin{align}
    \frac{1}{T} \sum_{t=1}^T \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right] = \mathcal{O}\left(\frac{1}{\sqrt{KT}} + \frac{K}{T}\right).
\end{align}
\end{corollary}
We show empirically in the next section that considering a regularization parameter such that $\mu \geq 1$ still leads to significant gains in terms of the worst distribution test accuracy and improves fairness across devices compared to DSGD.
\section{Experiments} \label{eval}
In this section, we validate our theoretical results and show the communication-efficiency, robustness and fairness of our proposed approach, DR-DSGD, compared to its non-robust counterpart DSGD. 
\subsection{Simulation Settings} 
For our experiments, we consider the image classification task using two main datasets: \textsc{Fashion MNIST} \citep{fashionmnist} and \textsc{CIFAR10} \citep{krizhevsky2009learning}. We implement DR-DSGD and DSGD algorithms using PyTorch. For \textsc{Fashion MNIST}, we use an MLP model with ReLU activations having two hidden layers with 128 and 64 neurons, respectively. For the \textsc{CIFAR10} dataset, we use a CNN model composed of three convolutional layers followed by two fully connected layers, each having 500 neurons. For each dataset, we distribute the data across the $K$ devices in a pathological non-IID way, as in \citep{fedavg}, to mimic an actual decentralized learning setup. More specifically, we first order the samples according to the labels and divide data into shards of equal sizes. Finally, we assign each device the same number of chunks. This will ensure a pathological non-IID partitioning of the data, as most devices will only have access to certain classes and not all of them. Unless explicitly stated, we choose the learning rate $\eta = \sqrt{K/T}$ and the mini-batch size $B=\sqrt{KT}$.

For the graph generation, we generate randomly a network consisting of $K$ devices with a connectivity ratio $p$ using the networkx package \citep{networkx}. The parameter $p$ measures the sparsity of the graph. While smaller values of $p$ lead to a sparser graph, the generated graph becomes denser as $p$ approaches $1$. We use the Metropolis weights to construct the mixing matrix $\bm{W}$ as follow
\begin{align*}
    W_{ij}=\left\{
         \begin{array}{ll}
           1/\left(1+\max\{d_i,d_j\}\right),     &\text{ if } (j,i)\in \mathcal{E},\\
           0,                              &\text{ if } (j,i)\notin \mathcal{E} \text{ and } j\neq i,\\
           1-\sum_{l\in\mathcal{N}_i} W_{il},&\text{ if } j=i,
         \end{array}
       \right.
\end{align*}
Unless otherwise stated, the graphs used in our experiments are of the Erdős-Rényi type.
\begin{figure*}
\caption{Performance comparison between DR-DSGD and DSGD in terms of: (a) average test accuracy, (b) worst test accuracy, and (c) STDEV of test accuracy for \textsc{Fashion MNIST}
\label{fig1:a}
\label{fig1:b}
\label{fig1:c}
\label{fig1}
\end{figure*}

\begin{figure*}
\caption{Performance comparison between DR-DSGD and DSGD in terms of: (a) average test accuracy, (b) worst test accuracy, and (c) STDEV of test accuracy for \textsc{CIFAR10}
\label{fig2:a}
\label{fig2:b}
\label{fig2:c}
\label{fig2}
\end{figure*}

\subsection{Robustness \& Communication-Efficiency} 
In this section, we consider $K=10$ devices and $\mu=6$. For \textsc{Fashion MNIST}, we consider a value of $p=0.3$ while we take $p=0.5$ for \textsc{CIFAR10}. For each experiment, we report both the average test accuracy, the worst distribution test accuracy, and their corresponding one standard error shaded area based on five runs. The worst distribution test accuracy is defined as the worst of all test accuracies. The performance comparison between DR-DSGD and DSGD for \textsc{Fashion MNIST} and \textsc{CIFAR10} dataset is reported in Figs. \ref{fig1} and \ref{fig2}, respectively. Both experiments show that while DR-DSGD achieves almost the same average test accuracy as DSGD, it significantly outperforms DSGD in terms of the worst distribution test accuracy. For the gap between both algorithms, the improvement, in terms of the worst distribution test accuracy, is of the order of $7\%$ for \textsc{Fashion MNIST} and $10\%$ for \textsc{CIFAR10} datasets, respectively. This is mainly due to the fact that while the ERM objective sacrifices the worst-case performance for the average one, DRO aims to lower the variance while maintaining good average performance across the different devices. Not only our proposed algorithm achieves better performance than DSGD, but it is also more communication-efficient. In fact, for the same metric requirement, DR-DSGD requires fewer communication rounds than DSGD. Since our approach exponentially increases the weight of high training loss devices, it converges much faster than DSGD. For instance, in the experiment using the \textsc{Fashion MNIST} dataset, DR-DSGD requires $10\times$ fewer iterations than DSGD to achieve $70\%$ worst distribution test accuracy. Finally, in each experiment, we plot the standard deviation (STDEV) of the different devices' test accuracies for both algorithms. We can see from both Figs. \ref{fig1:c} and \ref{fig2:c} that DR-DSGD has a smaller STDEV compared to DSGD, which reflects that DR-DSGD promotes more fairness among the devices. 

\subsection{Fairness} 
From this section on, we consider $K=25$. To investigate the fairness of the performance across the devices, we run the experiments on \textsc{Fashion MNIST} and \textsc{CIFAR10} datasets reporting the final test accuracy on each device in the case when $\mu=9$. In Figs. \ref{fig3:a} and \ref{fig3:b}, we plot the worst test accuracy distribution across devices. We note that DR-DSGD results in a more concentrated distribution in both experiments, hence a fairer test accuracy distribution with lower variance. For instance, DR-DSGD reduces the variance of accuracies across all devices by $60\%$ on average for the \textsc{Fashion MNIST} experiment while keeping almost the same average accuracy.

\begin{figure}
\caption{Performance comparison between DR-DSGD and DSGD in terms of worst test accuracy distribution for: (a) \textsc{Fashion MNIST}
\label{fig3:a}
\label{fig3:b}
\label{fig3}
\end{figure}
\subsection{Tradeoff Between Fairness \& Average Test Accuracy} In this section, we show how $\mu$ controls the trade-off between fairness and average test accuracy. To this end, we report, in Table \ref{table:1}, the average, and worst $10\%$ test accuracy, as well as the STDEV based on five runs for $T=300$ and for different values of $\mu$ for both datasets. As expected, higher values of $\mu$ give more weight to the regularization term; hence driving the values of $\lambda_i$ closer to the average weight $1/K$. Therefore, as $\mu$ increases, the average test accuracy increases but the worst ($10\%$) test accuracy decreases. Conversely, the worst test accuracy increases as the value of $\mu$ decreases at the cost of a drop in the average test accuracy. Furthermore, the STDEV decreases for smaller values of $\mu$ ensuring a fairer performance across devices. 

\subsection{Impact of Graph Topology} 
In this section, we start by inspecting the effect of the graph sparsity on the performance of DR-DSGD and DSGD in terms of the worst distribution test accuracy by considering different connectivity ratios $p \in \{0.3, 0.45, 0.6\}$ and for $\mu=6$. The results are reported in Fig \ref{fig4} for both datasets for $\mu=6$. We can see that as the graph becomes denser, i.e. as $p$ increases, the performance of both algorithms improves in terms of the worst test distribution. Nonetheless, it is clear that DR-DSGD outperforms DSGD for three different values of $p$ for both datasets. Next, we explore the performance of both algorithms on several other types of graph-types, specifically geometric (Fig. \ref{fig5:a}), ring (Fig. \ref{fig5:b}) and grid graphs (Fig. \ref{fig5:c}). The first row represents the graph topology, while the second represents the worst distribution test accuracy as a function of the number of communication rounds for \textsc{Fashion MNIST}. We can see that DR-DSGD outperforms DSGD for the three graphs considered in Fig. \ref{fig5} by achieving a higher worst distribution test accuracy. Furthermore, we note that both algorithms converge faster as the graph becomes denser. For instance, both algorithms require fewer communication rounds when the graph topology is geometric (Fig. \ref{fig5:a}) compared to the ring graph (Fig. \ref{fig5:b}).

\begin{table}[t]
\caption{Statistics of the test accuracy distribution for different values of $\mu$.}\label{table:1}
\centering
\begin{tabular}{ l | c | ccc}
    \toprule
    \multirow{2}{*}{\textbf{Dataset}} & \multirow{2}{*}{$\mu$} & \textbf{Average} & \textbf{\small{Worst 10\%}} & \textbf{STDEV} \\
    & & (\%) & (\%) &  \\
    \hline
    \multirow{3}{*}{\textsc{FMNIST}} & $\mu=2$ & 71.5 {\scriptsize $\pm$ 1.3}  & \textbf{49.1} {\scriptsize $\pm$ 2.4} &  \textbf{11.4} {\scriptsize $\pm$ 0.3} \\
    & $\mu=3$ & 72.3 {\scriptsize $\pm$ 1.1}  & 48.8 {\scriptsize $\pm$ 3}  &  11.8 {\scriptsize $\pm$ 1.5} \\
    & $\mu=5$ & \textbf{73.4} {\scriptsize $\pm$ 2.4}  & 44.5 {\scriptsize $\pm$ 4.4}  &  13.4  {\scriptsize $\pm$ 2.1} \\
    \hline
    \multirow{3}{*}{\textsc{CIFAR10}} & $\mu=2$ & 57.2 {\scriptsize $\pm$ 2.4}  & \textbf{50.9} {\scriptsize $\pm$ 1.5} &  \textbf{10.3} {\scriptsize $\pm$ 0.6} \\
    & $\mu=3$ & 59.83 {\scriptsize $\pm$ 1.6}  & 48.9 {\scriptsize $\pm$ 1.8}  &  10.9 {\scriptsize $\pm$ 0.4} \\
    & $\mu=5$ & \textbf{61} {\scriptsize $\pm$ 1.4}  & 44.9 {\scriptsize $\pm$ 1.9}  &  11.2  {\scriptsize $\pm$ 1.3} \\
    \bottomrule
\end{tabular}
\end{table}

\begin{figure}
\caption{Performance comparison between DR-DSGD and DSGD in terms of worst test accuracy distribution for different values of $p$ for: (a) \textsc{Fashion MNIST}
\label{fig4:a}
\label{fig4:b}
\label{fig4}
\end{figure}

\begin{figure*}
\caption{Performance comparison between DR-DSGD and DSGD in terms of worst test accuracy for: (a) geometric graph, (b) ring graph, and (c) grid graph for \textsc{Fashion MNIST}
\label{fig5:a}
\label{fig5:b}
\label{fig5:c}
\label{fig5:d}
\label{fig5:e}
\label{fig5:f}
\label{fig5}
\end{figure*}
\subsection{Limitations}\label{lim}
In this section, we highlight some of the limitations of our approach, while outlining several potential directions for future work.
\begin{itemize}
    \item \textbf{Solving the unregularized distributionally robust problem:} Instead of solving (\ref{eq:3}), we propose DR-DSGD to solve its regularized form given in (\ref{eq:3}). This is mainly motivated by the fact that solving the min-max problem in (\ref{eq:3}) in a decentralized fashion is a challenging task. A weakness of the convergence analysis is that the derived convergence rate is obtained only for values of $\mu$ larger than one. For $\mu \rightarrow 0$, the bound, derived in Theorem \ref{thm1}, grows to infinity.
    \item \textbf{Relaxing some of the assumptions:} As pointed out in \citep{khaled2020tighter}, the bounded gradient assumption has been criticised in the FL literature. Adopting a more reasonable assumption is left for future work.
\end{itemize}
\section{Conclusion}\label{SecConc}
This paper proposes a distributionally robust decentralized algorithm, DR-DSGD, that builds upon the decentralized stochastic gradient descent (DSGD) algorithm. The proposed framework is the first to solve the distributionally robust learning problem over graphs. Simulation results indicate that our proposed algorithm is more robust across heterogeneous data distributions while being more communication-efficient than its non-robust counterpart, DSGD. Furthermore, the proposed approach ensures fairer performance across all devices compared to DSGD. 
\section*{Broader Impact Statement}
This paper proposes a distributionally-robust decentralized learning algorithm. Our approach minimizes the maximum loss among the worst-case distributions across devices' data. As a result, even if the data distribution across devices is significantly heterogeneous, DR-DSGD guarantees a notion of fairness across devices. More specifically, the proposed algorithm ensures that the trained model has almost similar performance for all devices, rather than just performing well on a few of them while doing poorly on other devices.
\section*{Acknowledgments} 
This work was supported in part by the Academy of Finland
6G Flagship under grant No. 318927, in part by project SMARTER, in part
by projects EU-ICT IntellIoT under grant No. 957218, EUCHISTERA LearningEdge, CONNECT, Infotech-NOOR, and NEGEIN. We would like also to thank the
anonymous reviewers for their insightful comments that led to improving the manuscript.

\appendix
\section{Appendix}
\subsection{Worst-case Bound on Categorical Cross-entropy}\label{Appendix:bound}
Let's first examine the behavior of a randomly initialized network. With random weights, the many units/layers will usually compound to result in the network outputting approximately uniform predictions. In a classification problem with $M$ classes, we will get probabilities of around $1/M$ for each category. In fact, the cross-entropy for a single data point is defined as
\begin{align}
    CE = -\sum_{m=1}^M y_m \log(\hat{y}_m),
\end{align}
where $y$ are the true probabilities (labels), $\hat{y}$ are the predictions. With hard labels (i.e. one-hot encoded), only a single $y_m$ is 1, and all others are 0. Thus, $CE$ reduces to $-\log(\hat{y}_m)$, where $m$ is now the correct class. With a randomly initialized network, we have $\hat{y}_m \sim 1/M$, therefore, we get $-\log(1/M) = \log(M)$. Since the training objective is usually to reduce cross-entropy, we can think of $\log(M)$ as a worst-case bound.
\subsection{Proximity of the Solutions of (\ref{eq:4}) to (\ref{eq:3})}\label{Appendix:prox}
In this appendix, we discuss briefly the Proximity of the solutions of (\ref{eq:4}) to (\ref{eq:3}) in the \textbf{convex} case. According to \citep{qian2019robust}, problem (\ref{eq:4}) can be re-written using the duality theory as
\begin{align}
    \min_{\Theta \in \mathbb{R}^d} \max_{\lambda \in \Delta} \sum_{i=1}^K \lambda_i f_i(\Theta) \\
    s.t. ~ \phi(\lambda,1/K) \leq \tau 
\end{align}
where to every $\mu$, we associate a $\tau$ value. From \citep{namkoong2016stochastic, duchi2016statistics}, for \textbf{convex} loss functions, we have that
\begin{align}
    \max_{\lambda \in \mathcal{P}_{\rho,K}} \sum_{i=1}^K \lambda_i \ell_i(\Theta) = \frac{1}{K} \sum_{i=1}^K  \ell_i(\Theta) + \sqrt{\tau Var_{P_0}[\ell(\Theta, \xi)]} + o_{P_0}(K^{-\frac{1}{2}}),
\end{align}
where $\ell(\Theta) = [\ell_1(\Theta), \dots, \ell_K(\Theta)]^T \in \mathbb{R}^K$ is the vector of losses, $\mathcal{P}_0$ is the empirical probability distribution,  and the ambiguity set $\mathcal{P}_{\rho,K}$ is defined as $\mathcal{P}_{\rho,K} = \{\lambda \in \mathbb{R}^K \sum_{i=1}^K \lambda_i = 1, \lambda \geq 0, \phi(\lambda,1/K) \leq \tau \}$. 
Note that problem (\ref{eq:3}) is equivalent to (\ref{eq:4}) when $\tau = 0$. Thus, we can link the objective of the problem (\ref{eq:4}) to (\ref{eq:3}) as
\begin{align}
    \min_{\Theta \in \mathbb{R}^d} \max_{\lambda \in \mathcal{P}_{\rho,K}} \sum_{i=1}^K \lambda_i f_i(\Theta) = \min_{\Theta \in \mathbb{R}^d} \max_{\lambda \in \mathcal{P}_{0,K}} \sum_{i=1}^K \lambda_i f_i(\Theta) +  \sqrt{\tau Var_{P_0}[\ell(\Theta, \xi)]}.
\end{align}
\subsection{Basic identities and inequalities}\label{ineq}
We start by summarizing the main identities and inequalities used in the proof. Let $\{\bm{a}_s\}_{s=1}^S$ be a sequence of vectors in $\mathbb{R}^d$, $b_1$ and $b_2$ two scalars, $\bm{C}_1$ and $\bm{C}_2$ two matrices, and $\epsilon>0$, then we have
\begin{align}\label{id1}
    \left\|\sum_{s=1}^S \bm{a}_s\right\|^2 \leq S \sum_{s=1}^S \|\bm{a}_s\|^2.
\end{align}
\begin{align}\label{id2}
    2 \langle \bm{a}_1, \bm{a}_2 \rangle = \|\bm{a}_1\|^2 + \|\bm{a}_2\|^2-\|\bm{a}_1-\bm{a}_2\|^2.
\end{align}
\begin{align}\label{id3}
    2 b_1 b_2 \leq \frac{b_1^2}{\epsilon} + \epsilon b_2^2, \forall \epsilon > 0.
\end{align}
\begin{align}\label{id4}
    \text{(Cauchy-Schwarz)} \quad |\text{Tr}\{\bm{C}_1 \bm{C}_2\}| \leq \|\bm{C}_1\|_F \|\bm{C}_2\|_F.
\end{align}
\subsection{Proof of Lemma \ref{lemma1}}\label{proofLemma1}
Let $\bm{a}_i^T$ denote the i$^{th}$ row vector of matrix $\bm{A}$ and $\bm{e}_i$ the i$^{th}$ vector of the canonical basis of $\mathbb{R}^K$, then we can write
\begin{align}\label{lemeq1}
    \nonumber &\mathbb{E}\left[\|\bm{A} \left( \bm{W}^n- \bm{J}\right)\|_F^2\right]\\
    \nonumber &= \sum_{i=1}^K \left\|\bm{a}_i^T\left(\bm{W}^n \bm{e}_i - \frac{\bm{1}}{K}\right)\right\|^2\\
    &\leq \sum_{i=1}^K \|\bm{a}_i^T\|^2 \left\|\bm{W}^n \bm{e}_i - \frac{\bm{1}}{K}\right\|^2.
\end{align}
From \citep[Lemma 5]{lian2017can}, we have
\begin{align}\label{lemeq2}
\left\|\bm{W}^n \bm{e}_i - \frac{\bm{1}}{K}\right\|^2 \leq \rho^n.    
\end{align}
Replacing \eqref{lemeq2} in \eqref{lemeq1}, we get
\begin{align}
\mathbb{E}\left[\|\bm{A} \left( \bm{W}^n- \bm{J}\right)\|_F^2\right] \leq \rho^n \|\bm{A}\|_F^2.    
\end{align}
Hence, the proof is completed.
\subsection{Proof of Lemma \ref{lemma3}}\label{proofLemma3}
From Assumption 4, we have that the function $\theta \mapsto \ell(\theta, \xi)$ is bounded, then $\theta \mapsto \exp(\ell(\theta, \xi))$ is $C^1$ on a compact, and as a consequence, there exists $L_e > 0$ such that
\begin{align}\label{le22}
    \mathbb{E}[|\exp\left(z_1\right)\!-\! \exp\left(z_2\right) |] \leq L_e |z_1-z_2|,
\end{align}
where $z_i = \ell(\theta_i, \xi)$. Next, let $z_1 > z_2$ and define $t_i = \exp(z_i)$, then we can write
\begin{align}
    \mathbb{E}[t_1^{\frac{1}{\mu}}- t_2^{\frac{1}{\mu}}] = \frac{1}{\mu} \mathbb{E}\left[\int_{t_2}^{t_1} x^{\frac{1}{\mu}-1} dx\right] \leq \frac{\exp\left((\frac{1}{\mu}-1)M\right)}{\mu}  \mathbb{E}[t_1-t_2] = \frac{\exp\left((\frac{1}{\mu}-1)M\right)}{\mu} \mathbb{E}\left[\exp\left(z_1\right)-\exp\left(z_2\right)\right].
\end{align}
Using (\ref{le22}), we get
\begin{align}
    \mathbb{E}[t_1^{\frac{1}{\mu}}- t_2^{\frac{1}{\mu}}] \leq \exp\left((\frac{1}{\mu}-1)M\right) L_e \mathbb{E}\left[\frac{z_1}{\mu}-\frac{z_2}{\mu}\right]
\end{align}
A similar result can be obtained if we considered $z_1 < z_2$, hence, we can write
\begin{align}
    \mathbb{E}\left[|t_1^{\frac{1}{\mu}}- t_2^{\frac{1}{\mu}}|\right] \leq \exp\left((\frac{1}{\mu}-1)M\right) L_e \mathbb{E}\left[\left|\frac{z_1}{\mu}-\frac{z_2}{\mu}\right|\right]. 
\end{align}
Using the definition of $\mathcal{Y}$, we can write
\begin{align}
    \mathbb{E}\left[|\exp(y_1) - \exp(y_2)|\right] \leq L_2(\mu) |y_1-y_2|,
\end{align}
where $L_2(\mu) := \exp\left((\frac{1}{\mu}-1)M\right) L_e$ and hence the proof of statement $(a)$.\\
For statement $(b)$, using Hoeffding's lemma, we can write
\begin{align}
    \mathbb{E}\left[\exp(\ell(\theta, \xi_i)/\mu) \right]\leq \exp\left(\frac{\mathbb{E}\left[\ell(\theta, \xi_i)\right]}{\mu} + \frac{M^2}{2 \mu^2}\right) \leq \exp\left(\frac{M}{\mu} + \frac{M^2}{2 \mu^2}\right).
\end{align}
By choosing $G_2(\mu) := \exp\left(\frac{M}{\mu} + \frac{M^2}{2 \mu^2}\right)$, we prove statement $(b)$. \\
Finally, for statement $(c)$, we can write
\begin{align}
    \mathbb{E}\left[|\exp(\ell(\bm{\theta}, \xi_i)/\mu)-\exp(f_i(\bm{\theta})/\mu)|^2\right] \leq 2 \left( \mathbb{E}\left[\exp(2 \ell(\bm{\theta}, \xi_i)/\mu)\right] +  \exp(2 f_i(\bm{\theta})/\mu)\right).
\end{align}
Using Assumption 4 and Hoeffding's lemma, we get
\begin{align}
    \mathbb{E}\left[|\exp(\ell(\bm{\theta}, \xi_i)/\mu)-\exp(f_i(\bm{\theta})/\mu)|^2\right] \leq (\sigma_{3}(\mu))^2,
\end{align}
where $(\sigma_{3}(\mu))^2 := 2 \left( \exp\left(\frac{2M}{\mu} + \frac{M^2}{\mu^2}\right) + \exp\left(\frac{2M}{\mu}\right) \right).$
\subsection{Proof of Lemma \ref{lemma2}}\label{proofLemma2}
Using the update rule (\ref{newupdate}) and the identity $\bm{W} \bm{J} = \bm{J} \bm{W} = \bm{J}$, we can write
\begin{align}\label{lem2eq1}
    \bm{\theta}^t (\bm{I}-\bm{J}) = \left(\bm{\theta}^{t-1}-\eta \bm{U}^{t-1}\right) \bm{W} (\bm{I}-\bm{J}) = \bm{\theta}^{t-1} (\bm{I}-\bm{J}) \bm{W} - \eta \bm{U}^{t-1} \bm{W} (\bm{I}-\bm{J}).
\end{align}
Writing (\ref{lem2eq1}) recursively, we get
\begin{align}
    \bm{\theta}^t (\bm{I}-\bm{J}) = \bm{\theta}^{0} (\bm{I}-\bm{J}) \bm{W}^{t} - \eta \sum_{\tau=0}^{t-1} \bm{U}^{\tau} \left(\bm{W}^{t-\tau}-\bm{J}\right) = - \eta \sum_{\tau=0}^{t-1} \bm{U}^{\tau} \left(\bm{W}^{t-\tau}-\bm{J}\right),
\end{align}
where we used the fact that all local models are initiated at the same point, i.e. $\bm{\theta}^{0} (\bm{I}-\bm{J}) \bm{W}^{t} = \bm{0}$.\\
Thus, we can write
\begin{align}\label{lem2eq2}
    \nonumber & \frac{1}{KT} \sum_{t=1}^T\mathbb{E}\left[\|\bm{\theta}^{t}(\bm{I}-\bm{J})\|_F^2 \right]\\
    \nonumber & = \frac{\eta^2}{KT} \sum_{t=1}^T \mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \bm{U}^{\tau} \left(\bm{W}^{t-\tau}-\bm{J}\right)\right\|_F^2\right] \\
    \nonumber & = \frac{\eta^2}{KT} \sum_{t=1}^T \mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \left(\bm{U}^{\tau}-\nabla \bm{F}^\tau + \nabla \bm{F}^\tau \right) \left(\bm{W}^{t-\tau}-\bm{J}\right)\right\|_F^2\right] \\
    & \leq \frac{2 \eta^2}{KT} \sum_{t=1}^T \mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \left(\bm{U}^{\tau}-\nabla \bm{F}^\tau \right) \left(\bm{W}^{t-\tau}-\bm{J}\right)\right\|_F^2\right] + \frac{2 \eta^2}{KT} \sum_{t=1}^T \mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \nabla \bm{F}^\tau \left(\bm{W}^{t-\tau}-\bm{J}\right)\right\|_F^2\right],
\end{align}
where we used (\ref{id1}) (for $S = 2$) in the last inequality. Let $\bm{B}_{\tau,t} = \bm{W}^{t-\tau}-\bm{J}$. We start by examining the first term of (\ref{lem2eq2}) by writing
\begin{align}\label{lem2eq3}
\nonumber & \mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \left(\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right) \bm{B}_{\tau,t}\right\|_F^2\right]\\
\nonumber & = \sum_{\tau=0}^{t-1} \mathbb{E}\left[\left\|\left(\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right) \bm{B}_{\tau,t}\right\|_F^2\right] + \sum_{\tau=0}^{t-1} \sum_{\substack{\tau'=0,\\ \tau' \neq \tau}}^{t-1} \mathbb{E}\left[\text{Tr}\{\bm{B}_{\tau,t}^T (\bm{U}^{\tau}-\nabla \bm{F}^{\tau})^T (\bm{U}^{\tau'}-\nabla \bm{F}^{\tau'}) \bm{B}_{\tau',t}\}\right]\\
\nonumber & \leq \sum_{\tau=0}^{t-1} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2 \left\|\bm{B}_{\tau,t}\right\|_F^2\right]+ \sum_{\tau=0}^{t-1} \sum_{\substack{\tau'=0,\\ \tau' \neq \tau}}^{t-1} \mathbb{E}\left[\left\|\left(\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right) \bm{B}_{\tau,t} \right\|_F \|(\bm{U}^{\tau'}-\nabla \bm{F}^{\tau'}) \bm{B}_{\tau',t}\|_F\right]\\
& \leq \sum_{\tau=0}^{t-1} \rho^{t-\tau} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2 \right]+\sum_{\tau=0}^{t-1}\sum_{\substack{\tau'=0,\\ \tau' \neq \tau}}^{t-1} \left\{ \frac{\rho^{t-\tau}}{2\epsilon} \mathbb{E}\left[\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau}\|_F^2\right]+\frac{\epsilon \rho^{t-\tau'}}{2} \mathbb{E}\left[\left\|\bm{U}^{\tau'}-\nabla \bm{F}^{\tau'} \right\|_F^2\right]\right\},
\end{align}
where we have used Lemma \ref{lemma1} and inequalities (\ref{id3}) and (\ref{id4}). Setting $\epsilon = \rho^{\frac{\tau'-\tau}{2}}$, we can further write (\ref{lem2eq3}) as
\begin{align}\label{lem2main0}
\nonumber & \mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \left(\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right) \bm{B}_{\tau,t}\right\|_F^2\right]\\   
\nonumber & \leq \sum_{\tau=0}^{t-1} \rho^{t-\tau} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2 \right]+\sum_{\tau=0}^{t-1} \sum_{\substack{\tau'=0,\\ \tau' \neq \tau}}^{t-1} \frac{\rho^{t-\frac{\tau+\tau'}{2}}}{2} \left( \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^\tau \right\|_F^2 + \left\| \bm{U}^{\tau'}-\nabla \bm{F}^{\tau'}\right\|_F^2\right]\right)\\
\nonumber & \leq \sum_{\tau=0}^{t-1} \rho^{t-\tau} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2 \right]+\sum_{\tau=0}^{t-1} \sum_{\substack{\tau'=0,\\ \tau' \neq \tau}}^{t-1} \rho^{t-\frac{\tau+\tau'}{2}} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2\right]\\
\nonumber & \leq \sum_{\tau=0}^{t-1} \rho^{t-\tau} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2 \right]+\sum_{\tau=0}^{t-1}  \rho^{\frac{t-\tau}{2}} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2\right] \sum_{\substack{\tau'=0,\\ \tau' \neq \tau}}^{t-1} \rho^{\frac{t-\tau'}{2}}\\
\nonumber &= \sum_{\tau=0}^{t-1} \rho^{t-\tau} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2 \right]+\sum_{\tau=0}^{t-1}  \rho^{\frac{t-\tau}{2}} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2\right] \left(\sum_{\tau'=0}^{t-1} \rho^{\frac{t-\tau'}{2}}-\rho^{\frac{t-\tau}{2}}\right)\\
\nonumber& \leq \sum_{\tau=0}^{t-1}  \rho^{\frac{t-\tau}{2}} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2\right] \sum_{\tau'=0}^{t-1} \rho^{\frac{t-\tau'}{2}}\\
& \leq \frac{\sqrt{\rho}}{1-\sqrt{\rho}} \sum_{\tau=0}^{t-1} \rho^{\frac{t-\tau}{2}} \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2\right],
\end{align}
where in the last inequality, we used $\sum_{\tau=0}^{t-1} \rho^{\frac{t-\tau}{2}} = \sqrt{\rho}^t +  \sqrt{\rho}^{t-1} + \dots + \sqrt{\rho} \leq \frac{\sqrt{\rho}}{1-\sqrt{\rho}}$. Now, let's focus on finding an upper bound for the term $\mathbb{E}\left[\left\|\bm{U}^{\tau}\!-\!\nabla \bm{F}^{\tau} \right\|_F^2\right]$. To this end, we start by writing
\begin{align}\label{lemm2eq4}
\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2 = \frac{1}{\mu^2}\sum_{i=1}^K \left\|h(\bm{\theta}_i^{\tau};\mu) \bm{g}_i(\bm{\theta}_i^{\tau}) - \exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right) \nabla f_i(\bm{\theta}_i^{\tau}) \right\|^2.
\end{align}
Next, we can write the following
\begin{align}\label{lemm2eq5}
\nonumber &h(\bm{\theta}_i^{\tau};\mu) \bm{g}_i(\bm{\theta}_i^{\tau})-\exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right) \nabla f_i(\bm{\theta}_i^{\tau})\\
\nonumber &=h(\bm{\theta}_i^{\tau};\mu) \bm{g}_i(\bm{\theta}_i^{\tau})-h(\bm{\theta}_i^{\tau};\mu) \bm{g}_i(\bar{\bm{\theta}}^{\tau})+h(\bm{\theta}_i^{\tau};\mu) \bm{g}_i(\bar{\bm{\theta}}^{\tau})-h(\bm{\theta}_i^{\tau};\mu) \nabla f_i(\bar{\bm{\theta}}^{\tau})+h(\bm{\theta}_i^{\tau};\mu) \nabla f_i(\bar{\bm{\theta}}^{\tau})\\
&-\exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right) \nabla f_i(\bar{\bm{\theta}}^{\tau})+\exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right) \nabla f_i(\bar{\bm{\theta}}^{\tau})-\exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right) \nabla f_i(\bm{\theta}_i^{\tau}).
\end{align}
Using the decomposition (\ref{lemm2eq5}) in (\ref{lemm2eq4}) and taking the expected value while using the inequality (\ref{id1}) (for $S=4$), we get 
\begin{align}\label{lem2main1}
\nonumber & \mathbb{E}\left[\left\|\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right\|_F^2\right]\\    
\nonumber &\leq\frac{4}{\mu^2}\sum_{i=1}^K\mathbb{E}\left[\left\|h(\bm{\theta}_i^{\tau};\mu) \left(\bm{g}_i(\bm{\theta}_i^{\tau})- \bm{g}_i(\bar{\bm{\theta}}^{\tau})\right)\right\|^2\right]+\frac{4}{\mu^2}\sum_{i=1}^K\mathbb{E}\left[\left\|h(\bm{\theta}_i^{\tau};\mu) \left(\bm{g}_i(\bar{\bm{\theta}}^{\tau})-\nabla f_i(\bar{\bm{\theta}}^i)\right)\right\|^2\right]\\
\nonumber &+\frac{4}{\mu^2}\sum_{i=1}^K\mathbb{E}\left[\left\|\nabla f_i(\bar{\bm{\theta}}^{\tau})\left(h(\bm{\theta}_i^{\tau};\mu)-\exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right) \right)\right\|^2\right]\\
\nonumber &+\frac{4}{\mu^2}\sum_{i=1}^K\mathbb{E}\left[\left\|\exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right)\left(\nabla f_i(\bar{\bm{\theta}}^{\tau})-\nabla f_i(\bm{\theta}_i^{\tau})\right)\right\|^2\right]\\
\nonumber &\leq \frac{4 (G_2(\mu))^2}{\mu^2} \sum_{i=1}^K \mathbb{E}\left[\|\bm{g}_i(\bm{\theta}_i^{\tau})- \bm{g}_i(\bar{\bm{\theta}}^{\tau})\|^2\right] +\frac{4 (G_2(\mu))^2}{\mu^2}\sum_{i=1}^K \mathbb{E}\left[\|\bm{g}_i(\bar{\bm{\theta}}^{\tau})- \nabla f_i(\bar{\bm{\theta}}^{\tau})\|^2\right] \\
\nonumber &+\frac{4G_1^2}{\mu^2}\sum_{i=1}^K \mathbb{E}\left[\left|h(\bm{\theta}_i^{\tau};\mu)-\exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right)\right|^2\right]+\frac{4 (G_2(\mu))^2}{\mu^2}\sum_{i=1}^K \mathbb{E}\left[\|\nabla f_i(\bar{\bm{\theta}}^{\tau})- \nabla f_i(\bm{\theta}_i^{\tau})\|^2\right]\\
\nonumber &\leq \frac{8 G_{\mu}^2 L_{\mu}^2}{\mu^2} \sum_{i=1}^K \mathbb{E}\left[\|\bar{\bm{\theta}}^{\tau}-\bm{\theta}_i^{\tau}\|^2\right] + \frac{8 G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2+1) K}{\mu^2 B} \\
&\leq \frac{8 G_{\mu}^2 L_{\mu}^2}{\mu^2} \mathbb{E}\left[\|\bm{\theta}^{\tau}(\bm{I}-\bm{J})\|_F^2\right] + \frac{8 G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2+1) K}{\mu^2},
\end{align}
where we used Assumptions 1-3 and defined the constants $\sigma_{\mu} =\max\{\sigma_1, \sigma_2,\sigma_3(\mu)\}$, $G_{\mu}=\max\{G_1,G_2(\mu)\}$ and $L_{\mu} = \max\{L_F, L_1, L_2(\mu)\}$. Going back to (\ref{lem2main0}), and using (\ref{lem2main1}), we can write
\begin{align}\label{lem2main3}
\nonumber &\frac{1}{KT} \sum_{t=1}^T \mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \left(\bm{U}^{\tau}-\nabla \bm{F}^{\tau} \right) \bm{B}_{\tau,t}\right\|_F^2\right]\\   
\nonumber & \leq \frac{8 \sqrt{\rho} G_{\mu}^2 L_{\mu}^2}{\mu^2 (1-\sqrt{\rho})} \frac{1}{KT} \sum_{t=1}^T \sum_{\tau=0}^{t-1} \rho^{\frac{t-\tau}{2}} \mathbb{E}\left[\|\bm{\theta}^{\tau}(\bm{I}-\bm{J})\|_F^2\right]+ \frac{8 \sqrt{\rho} G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2+1)}{\mu^2 (1-\sqrt{\rho}) T} \sum_{t=1}^T \sum_{\tau=0}^{t-1} \rho^{\frac{t-\tau}{2}}\\
\nonumber & \leq \frac{8 \sqrt{\rho} G_{\mu}^2 L_{\mu}^2}{\mu^2 (1-\sqrt{\rho})} \frac{1}{KT} \sum_{t=1}^T \mathbb{E}\left[\|\bm{\theta}^t(\bm{I}-\bm{J})\|_F^2\right]  \sum_{\tau=0}^{T-t} \rho^{\frac{\tau}{2}} + \frac{8 \sqrt{\rho} G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2+1)}{\mu^2 (1-\sqrt{\rho})  T} \sum_{t=1}^T \sum_{\tau=0}^{T-t} \rho^{\frac{\tau}{2}}\\
& \leq \frac{8 \rho G_{\mu}^2 L_{\mu}^2}{\mu^2 (1-\sqrt{\rho})^2} \frac{1}{KT} \sum_{t=1}^T \mathbb{E}\left[\|\bm{\theta}^t(\bm{I}-\bm{J})\|_F^2\right] + \frac{8 \rho G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2+1)}{\mu^2 (1-\sqrt{\rho})^2}.
\end{align}
Now, we focus on the second term of (\ref{lem2eq2}). Following similar steps as when bounding the first term of (\ref{lem2eq2}), we get
\begin{align}\label{001}
\mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \nabla \bm{F}^{\tau} \bm{B}_{\tau,t}\right\|_F^2\right] \leq \frac{\sqrt{\rho}}{1-\sqrt{\rho}} \sum_{\tau=0}^{t-1} \rho^{\frac{t-\tau}{2}} \mathbb{E}\left[\left\|\nabla \bm{F}^{\tau}\right\|_F^2\right].
\end{align}
Next, we look for an upper bound for the term $\mathbb{E}\left[\left\|\nabla \bm{F}^{\tau}\right\|_F^2\right]$. To this end, we start by writing
\begin{align}\label{002}
\nonumber &\mathbb{E}\left[\left\|\nabla \bm{F}^{\tau}\right\|_F^2\right]\\
\nonumber &= \frac{1}{\mu^2} \sum_{i=1}^K \mathbb{E}\left[\left\|\exp\left(\frac{f_i(\bm{\theta}_i^{\tau})}{\mu}\right)\nabla f_i(\bm{\theta}_i^{\tau})\right\|^2\right]\\
\nonumber &\leq \frac{(G_2(\mu))^2}{\mu^2} \mathbb{E}\left[\left\|\nabla f_i(\bm{\theta}_i^{\tau})\right\|^2\right]\\
\nonumber &\leq \frac{G1^2 (G_2(\mu))^2 K}{\mu^2} \\
&\leq \frac{G_{\mu}^4 K}{\mu^2}.
\end{align}
Therefore, we get
\begin{align}\label{lem2main2}
\frac{1}{KT} \sum_{t=1}^T \mathbb{E}\left[\left\|\sum_{\tau=0}^{t-1} \nabla \bm{F}^{\tau} \bm{B}_{\tau,t}\right\|_F^2\right] \leq \frac{G_{\mu}^4 \rho}{\mu^2 (1-\sqrt{\rho})^2}.
\end{align}
Next, using (\ref{lem2eq2}), (\ref{lem2main3}), and (\ref{lem2main2}), we can write
\begin{align}
\frac{1}{K T} \sum_{t=1}^T \mathbb{E}\left[\|\bm{\theta}^t(\bm{I}-\bm{J})\|_F^2\right] \leq \frac{16 \eta^2 \rho L_{\mu}^2 G_{\mu}^2}{\mu^2 (1-\sqrt{\rho})^2} \frac{1}{K T} \sum_{t=1}^T \mathbb{E}\left[\|\bm{\theta}^t(\bm{I}-\bm{J})\|_F^2\right]+\frac{2 \eta^2 \rho G_{\mu}^2 [8\sigma_{\mu}^2(L_{\mu}^2+1)+ G_{\mu}^2]}{\mu^2(1-\sqrt{\rho})^2}.
\end{align}
Let $\gamma_{\mu} = \frac{16 \eta^2 \rho L_{\mu}^2 G_{\mu}^2}{\mu^2 (1-\sqrt{\rho})^2}$, then we obtain
\begin{align}
\frac{1}{K T} \sum_{t=1}^T \mathbb{E}\left[\|\bm{\theta}^t(\bm{I}-\bm{J})\|_F^2\right]\leq \frac{2 \eta^2 \rho G_{\mu}^2 [8\sigma_{\mu}^2(L_{\mu}^2+1)+ G_{\mu}^2]}{\mu^2 (1-\gamma_{\mu}) (1-\sqrt{\rho})^2}, 
\end{align}
which concludes the proof of Lemma \ref{lemma2}.
\subsection{Proof of Theorem \ref{thm1}}\label{proofthm1}
Since the objective function $F(\cdot)$ has a Lipschitz gradient, we can write
\begin{align}
    F(\bar{\bm{\theta}}^{t+1}) - F(\bar{\bm{\theta}}^{t})
    \leq \langle {\nabla F(\bar{\bm{\theta}}^{t}}), {\bar{\bm{\theta}}^{t+1}}-\bar{\bm{\theta}}^{t} \rangle + \frac{L_{\mu}}{2} \|\bar{\bm{\theta}}^{t+1}-\bar{\bm{\theta}}^{t}\|^2. 
\end{align}
Plugging the update rule $\bar{\bm{\theta}}^{t+1} = \bar{\bm{\theta}}^{t} - \eta \bm{U}^{t}\bm{1}/K$, we have
\begin{align}\label{eqqq0}
    F(\bar{\bm{\theta}}^{t+1}) - F(\bar{\bm{\theta}}^{t})
    &\leq -\eta \langle \nabla F(\bar{\bm{\theta}}^{t}), \frac{\bm{U}^{t}\bm{1}}{K} \rangle + \frac{\eta^2L_{\mu}}{2}\left\|\frac{\bm{U}^{t}\bm{1}}{K}\right\|^2.
\end{align}
Using the identity (\ref{id2}), we can write the first term of the left hand-side of (\ref{eqqq0}) as
\begin{align}\label{eqq1}
\langle \nabla F(\bar{\bm{\theta}}^{t}), \frac{\bm{U}^{t}\bm{1}}{K} \rangle = \frac{1}{2} \left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2 + \left\|\frac{\bm{U}^{t}\bm{1}}{K}\right\|^2 - \left\|F(\bar{\bm{\theta}}^{t})-\frac{\bm{U}^{t}\bm{1}}{K} \right\|^2\right].
\end{align}
Next, we look for an upper bound for the third term of (\ref{eqq1}). To this end, we start by writing
\begin{align}
    \nonumber &\nabla F(\bar{\bm{\theta}}^{t}) - \frac{\bm{U}^{t}\bm{1}}{K}\\
    \nonumber & = \frac{1}{\mu K} \sum_{i=1}^K \left\{ \exp\left(\frac{f_i(\bar{\bm{\theta}}^{t})}{\mu}\right)\nabla f_i(\bar{\bm{\theta}}^{t})-h(\bm{\theta}_i^{t};\mu) g_i(\bm{\theta}_i^{t})\right\}\\
    \nonumber & = \frac{1}{\mu K} \sum_{i=1}^K \left\{ \exp\left(\frac{f_i(\bar{\bm{\theta}}^{t})}{\mu}\right)\nabla f_i(\bar{\bm{\theta}}^{t})- h(\bar{\bm{\theta}}^{t};\mu) \nabla f_i(\bar{\bm{\theta}}^{t})\right\} + \frac{1}{\mu K} \sum_{i=1}^K \left\{ h(\bar{\bm{\theta}}^{t};\mu) \nabla f_i(\bar{\bm{\theta}}^{t})-h(\bm{\theta}_i^{t};\mu) \nabla f_i(\bar{\bm{\theta}}^{t})\right\}\\
    &+ \frac{1}{\mu K} \sum_{i=1}^K \left\{ h(\bm{\theta}_i^{t};\mu) \nabla f_i(\bar{\bm{\theta}}^{t})-h(\bm{\theta}_i^{t};\mu) g_i(\bar{\bm{\theta}}^{t})\right\}+ \frac{1}{\mu K} \sum_{i=1}^K \left\{ h(\bm{\theta}_i^{t};\mu) g_i(\bar{\bm{\theta}}^{t})-h(\bm{\theta}_i^{t};\mu) g_i(\bm{\theta}_i^{t})\right\}.
\end{align}
Using the inequality (\ref{id1}), we get
\begin{align}
 \nonumber &\left\|\nabla F(\bar{\bm{\theta}}^{t}) - \frac{\bm{U}^{t}\bm{1}}{K}\right\|^2\\
    \nonumber &\leq\frac{4}{\mu^2 K} \sum_{i=1}^K \left\|\left(\exp\left(\frac{f_i(\bar{\bm{\theta}}^{t})}{\mu}\right) -h(\bar{\bm{\theta}}^{t};\mu)\right) \nabla f_i(\bar{\bm{\theta}}^{t})\right\|^2+\frac{4}{\mu^2 K} \sum_{i=1}^K \left\|\left(h(\bar{\bm{\theta}}^{t};\mu) -h(\bm{\theta}_i^{t};\mu)\right) \nabla f_i(\bar{\bm{\theta}}^{t})\right\|^2\\
    &+\frac{4}{\mu^2 K} \sum_{i=1}^K \left\|h(\bm{\theta}_i^{t};\mu) \left(\nabla f_i(\bar{\bm{\theta}}^{t})- g_i(\bar{\bm{\theta}}^{t})\right)\right\|^2
    +\frac{4}{\mu^2 K} \sum_{i=1}^K \left\|h(\bm{\theta}_i^{t};\mu) \left( g_i(\bar{\bm{\theta}}^{t})- g_i(\bm{\theta}_i^{t})\right)\right\|^2.
\end{align}
Using assumption 2 and taking the expected value from both sides, we can write
\begin{align}
     \nonumber &\mathbb{E}\left[\left\|\nabla F(\bar{\bm{\theta}}^{t}) - \frac{\bm{U}^{t}\bm{1}}{K}\right\|^2\right]\\
     \nonumber &\leq \frac{4 G_1^2}{\mu^2 K}  \sum_{i=1}^K \underbrace{\mathbb{E}\left[\left|\exp\left(\frac{f_i(\bar{\bm{\theta}}^{t})}{\mu}\right)-h(\bar{\bm{\theta}}^{t};\mu) \right|^2\right]}_\text{$T_1$}+\frac{4 G_1^2}{\mu^2 K} \sum_{i=1}^K \underbrace{\mathbb{E}\left[\left|h(\bar{\bm{\theta}}^{t};\mu) -h(\bm{\theta}_i^{t};\mu)\right|^2\right]}_\text{$T_2$}\\
     &+\frac{4 (G_2(\mu))^2}{\mu^2 K} \sum_{i=1}^K \underbrace{\mathbb{E}\left[\left\|\nabla f_i(\bar{\bm{\theta}}^{t})-g_i(\bar{\bm{\theta}}^{t})\right\|^2\right]}_\text{$T_3$}+\frac{4 (G_2(\mu))^2}{\mu^2 K} \sum_{i=1}^K \underbrace{\mathbb{E}\left[\left\|g_i(\bar{\bm{\theta}}^{t})- g_i(\bm{\theta}_i^{t})\right\|^2\right]}_\text{$T_4$}.
\end{align}
We start by bounding the term $T_1$
\begin{align}
    \nonumber T_1 &= \mathbb{E}\left[\left|\exp\left(\frac{f_i(\bar{\bm{\theta}}^{t})}{\mu}\right)-h(\bar{\bm{\theta}}^{t};\mu) \right|^2\right]\\
    \nonumber & \stackrel{(\ref{as1})}{\leq} L_{\mu}^2 \mathbb{E}\left[\left|\frac{1}{B} \sum_{i=1}^B \frac{\ell(\bar{\bm{\theta}}^{t}, \xi_i^t)}{\mu}- \frac{f_i(\bar{\bm{\theta}}^{t})}{\mu}\right|^2\right]\\
    \nonumber &= \frac{L_{\mu}^2}{B^2 \mu^2}  \mathbb{E}\left[\left|\sum_{i=1}^B \left(\ell(\bar{\bm{\theta}}^{t}, \xi_i^t)- f_i(\bar{\bm{\theta}}^{t}) \right)\right|^2\right]\\
    \nonumber & = \frac{L_{\mu}^2}{B^2 \mu^2} \mathbb{E}\left[\left|\sum_{i=1}^B X_i^t 
    \right|^2\right] \\
    & = \frac{L_{\mu}^2}{B^2 \mu^2} \left( \sum_{i=1}^B \mathbb{E}\left[\left| X_i^t 
    \right|^2\right] + \sum_{j \neq i} \mathbb{E}\left[\langle X_i^t, X_j^t 
    \rangle\right] \right),
\end{align}
where $X_i^t = \ell(\bar{\bm{\theta}}^{t}, \xi_i^t)- f_i(\bar{\bm{\theta}}^{t})$. Since $\xi_i^t$ and $\xi_j^t$ are independent for $j \neq i$, then $\mathbb{E}\left[\langle X_i^t, X_j^t 
    \rangle\right] = 0$. Then, using (\ref{ff3}), we can write that
    \begin{align}
        T_1  = \frac{L_{\mu}^2}{B^2 \mu^2} \sum_{i=1}^B \mathbb{E}\left[\left| X_i^t 
    \right|^2\right] \leq \frac{L_{\mu}^2}{B^2 \mu^2} (B \sigma_{\mu}^2) = \frac{L_{\mu}^2 \sigma_{\mu}^2}{B \mu^2}.
    \end{align}
Similarly, we can bound the term $T_3$ as $T_3 \leq \frac{\sigma_{\mu}^2}{B}$. Next, we bound the term $T_2$ as
\begin{align}
    \nonumber T_2 &= \mathbb{E}\left[\left|h(\bar{\bm{\theta}}^{t};\mu) -h(\bm{\theta}_i^{t};\mu)\right|^2\right]\\
    \nonumber & = \mathbb{E}\left[\left|\exp\left(\frac{1}{B} \sum_{j=1}^B \frac{\ell(\bar{\bm{\theta}}^{t}, \xi_j^t)}{\mu}\right)-\exp\left(\frac{1}{B} \sum_{j=1}^B \frac{\ell(\bm{\theta}_i^{t}, \xi_j^t)}{\mu}\right)\right|^2\right]\\
    \nonumber & \stackrel{(\ref{as1})}{\leq} L^2 \mathbb{E}\left[\left|\frac{1}{B} \sum_{j=1}^B \frac{\ell(\bar{\bm{\theta}}^{t}, \xi_j^t)}{\mu}-\frac{1}{B} \sum_{j=1}^B \frac{\ell(\bm{\theta}_i^{t}, \xi_j^t)}{\mu}\right|^2\right] \\
    \nonumber & = \frac{L_{\mu}^2}{B^2 \mu^2} \mathbb{E}\left[\left|\sum_{j=1}^B ( \ell(\bar{\bm{\theta}}^{t}, \xi_j^t)- \ell(\bm{\theta}_i^{t}, \xi_j^t))\right|^2\right]\\
     \nonumber & \leq \frac{L_{\mu}^2}{B \mu^2} \sum_{j=1}^B  \mathbb{E}\left[\left|\ell(\bar{\bm{\theta}}^{t}, \xi_j^t)- \ell(\bm{\theta}_i^{t}, \xi_j^t)\right|^2\right] \\
    & \stackrel{\mathrm{(a)}}{\leq} \frac{L_{\mu}^2 G_{\mu}^2}{\mu^2}  \mathbb{E}\left[\left\|\bar{\bm{\theta}}^{t}-\bm{\theta}_i^{t}\right\|^2\right],
\end{align}
where we have used the fact that (\ref{as000}) gives that $\ell(\bar{\bm{\theta}}^{t}, \xi_j^t)$ is $G$-Lipschitz continuous in $\mathrm{(a)}$ and (\ref{ff0}) in the last inequality. Finally, using (\ref{ff1}), we can bound $T_4$ as 
$T_4 \leq L_{\mu}^2 \mathbb{E}\left[\|\bar{\bm{\theta}}^{t}-\bm{\theta}_i^{t}\|^2\right]$. Using the bounds on the terms $T_1$-$T_4$, we can write
\begin{align}\label{eqq2}
     \mathbb{E}\left[\left\|\nabla F(\bar{\bm{\theta}}^{t}) - \frac{\bm{U}^{t}\bm{1}}{K}\right\|^2\right]\leq \frac{4 G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2 + \mu^2)}{\mu^4 B} +  \frac{4 G_{\mu}^2 (G_{\mu}^2+\mu^2) L_{\mu}^2}{\mu^4 K} \sum_{i=1}^K \mathbb{E}\left[\|\bar{\bm{\theta}}^{t}-\bm{\theta}_i^{t}\|^2\right].
\end{align}
Replacing (\ref{eqq2}) in (\ref{eqq1}), we obtain
\begin{align}
\mathbb{E}\left[\langle \nabla F(\bar{\bm{\theta}}^{t}), \frac{\bm{U}^{t}\bm{1}}{K} \rangle\right]\geq \frac{1}{2} \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right] + \frac{1}{2} \mathbb{E}\left[\left\|\frac{\bm{U}^{t}\bm{1}}{K}\right\|^2\right]-\frac{2 G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2 + \mu^2)}{\mu^4 B}-\frac{2 G_{\mu}^2 (G_{\mu}^2+\mu^2) L_{\mu}^2}{\mu^4 K}  \mathbb{E}\left[\|\bm{\theta}^{t}(\bm{I}-\bm{J})\|_F^2\right].
\end{align}
Going back to (\ref{eqqq0}), we can write
\begin{align}
    \nonumber &\mathbb{E}\left[ F(\bar{\bm{\theta}}^{t+1}) - F(\bar{\bm{\theta}}^{t})\right]\\
    &\leq - \frac{\eta}{2} \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right] + \frac{\eta}{2} \left(\eta L_{\mu}-1\right) \mathbb{E}\left[\left\|\frac{\bm{U}^{t}\bm{1}}{K}\right\|^2\right] +\frac{2 G_{\mu}^2 \sigma_{\mu}^2 \eta (L_{\mu}^2 + \mu^2)}{\mu^4 B}+\frac{4 G_{\mu}^2 L_{\mu}^2 \eta (G_{\mu}^2+\mu^2)}{\mu^4 K}  \mathbb{E}\left[\|\bm{\theta}^{t}(\bm{I}-\bm{J})\|_F^2\right].
\end{align}
Setting the learning rate $\eta$ such that $\eta L_{\mu} \leq 1$, and taking the average over $t \in [1,T]$, we get
\begin{align}\label{main0}
    \frac{\mathbb{E}\left[ F(\bar{\bm{\theta}}^{T}) - F(\bar{\bm{\theta}}^{1})\right]}{T}\leq - \frac{\eta}{2 T} \sum_{t=1}^T \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right] + \frac{2 \eta G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2 + \mu^2) }{\mu^4 B} + \frac{2 \eta G_{\mu}^2 L_{\mu}^2 (G_{\mu}^2+\mu^2)}{\mu^4 K T} \sum_{t=1}^T \mathbb{E}\left[\|\bm{\theta}^{t}(\bm{I}\!-\!\bm{J})\|_F^2\right].
\end{align}
Re-arranging the terms and using assumption 5, we can write
\begin{align}\label{main1}
    \frac{1}{T} \sum_{t=1}^T \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right] \leq \frac{2 (F(\bar{\bm{\theta}}^{1})-F_{inf})}{\eta T} + \frac{2 G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2 + \mu^2) }{\mu^4 B} + \frac{2 G_{\mu}^2 L_{\mu}^2 (G_{\mu}^2+\mu^2)}{\mu^4 K T} \sum_{t=1}^T \mathbb{E}\left[\|\bm{\theta}^{t}(\bm{I}-\bm{J})\|_F^2\right].
\end{align}
Using Lemma \ref{lemma2} in (\ref{main1}), we get
\begin{align}
    \frac{1}{T} \sum_{t=1}^T \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right] \leq \frac{2 (F(\bar{\bm{\theta}}^{1})-F_{inf})}{\eta T} +  \frac{2 G_{\mu}^2 \sigma_{\mu}^2 (L_{\mu}^2 + \mu^2) }{\mu^4 B} + \frac{4 \eta^2 L_{\mu}^2 \rho G_{\mu}^2 (G_{\mu}^2+\mu^2) [8\sigma_{\mu}^2(L_{\mu}^2+1)+ G_{\mu}^2]}{\mu^8 (1-\gamma) (1-\sqrt{\rho})^2}.
\end{align}
Furthermore, choosing $\eta$ such that $\eta L_{\mu} < \frac{\mu (1-\sqrt{\rho})}{8 G_{\mu} \sqrt{\rho}}$ ensures that $\gamma_{\mu} < \frac{1}{4}$, then we can write
\begin{align}
   \frac{1}{T} \sum_{t=1}^T \mathbb{E}\left[\|\nabla F(\bar{\bm{\theta}}^{t})\|^2\right] \leq \frac{2 (F(\bar{\bm{\theta}}^{1})-F_{inf})}{\eta T}+ \frac{2 G_{\mu}^2 \sigma_{\mu}^2(L_{\mu}^2 + \mu^2)}{\mu^4 B}+\frac{16 \eta^2 L_{\mu}^2 \rho G_{\mu}^2 (G_{\mu}^2+\mu^2) [8\sigma_{\mu}^2(L_{\mu}^2+1)+ G^2]}{3 \mu^8 (1-\sqrt{\rho})^2},
\end{align}
which finalizes the proof.

\end{document}