Does Worst-Performing Agent Lead the Pack? Analyzing Agent Dynamics in Unified Distributed SGD

Authors: jie  hu
yi-ting  ma
do  young
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-02 03:46:33

Abstract: Distributed learning is essential to train machine learning algorithms across *heterogeneous* agents while maintaining data privacy. We conduct an asymptotic analysis of Unified Distributed SGD (UD-SGD), exploring a variety of communication patterns, including decentralized SGD and local SGD within Federated Learning (FL), as well as the increasing communication interval in the FL setting. In this study, we assess how different sampling strategies, such as *i.i.d.* sampling, shuffling, and Markovian sampling, affect the convergence speed of UD-SGD by considering the impact of agent dynamics on the limiting covariance matrix as described in the Central Limit Theorem (CLT). Our findings not only support existing theories on linear speedup and asymptotic network independence, but also theoretically and empirically show how efficient sampling strategies employed by individual agents contribute to overall convergence in UD-SGD. Simulations reveal that a few agents using highly efficient sampling can achieve or surpass the performance of the majority employing moderately improved strategies, providing new insights beyond traditional analyses focusing on the worst-performing agent.

Main paper: 
 introduction: 
 Distributed learning deals with the training of models across multiple agents over a communication network in a distributed manner, while addressing the challenges of privacy, scalability, and highdimensional data [11, 55] . Each agent i ∈ [N ] holds a private dataset X i and an agent-specified loss function F i : R d × X i → R that depends on the model parameter θ ∈ R d and a data point X ∈ X i . The goal is then to find a local minima θ * of the objective function f (θ) ≜1 N N i=1 f i (θ), where agent i's loss function f i (θ) ≜ E X∼Di [F i (θ, X)] and D i represents the target distribution of data for agent i. 1 Each agent i can locally compute the gradient ∇F i (θ, X) ∈ R d w.r.t. θ for every sampled data point X ∈ X i . Due to the distributed nature, {D i } i∈[N ] and {X i } i∈ [N ] are not necessarily identically distributed over [N ] so that the minima of each local function f i (θ) can be far away from L. This is particularly relevant in decentralized training data, e.g., Federated Learning (FL) with heterogeneous data across data centers or devices [81, 31] .In this paper, we focus on Unified Distributed SGD (UD-SGD), where each agent i ∈ [N ] updates its model parameter θ i n+1 in a two-step process:Local update:EQUATIONAggregation:EQUATIONwhere γ n denotes the step size, X i n is the data sampled by agent i at time n (i.e., agent dynamics), and W n = [w n (i, j)] i,j∈ [N ] represents the doubly-stochastic communication matrix satisfying w n (i, j) ≥ 0 and 1 T W n = 1 T , W n 1 = 1. In the special case of N = 1, (1) simplifies to the vanilla SGD where W n = 1 for all n. UD-SGD covers a wide range of distributed algorithms, e.g., decentralized SGD (DSGD) [71, 80, 61, 68] , distributed SGD with changing topology (DSGD-CT) [24, 43] , local SGD (LSGD) in FL [55, 76] , and its variant aimed at reducing communication costs (LSGD-RC) [51] . 1 : GD-SGD algorithm with a communication network of N = 5 agents, each holding potentially distinct datasets; e.g., agent j (in blue) samples X j i.i.d. and agent i (in red) samples X i via Markovian trajectory.Versatile Communication Patterns {W n }: For visualization, we depict the scenarios of UD-SGD (1) in Figure 1 . In DSGD, each agent (node) in the graph communicates with its neighbors after each SGD computation via W n , representing the underlying network topology. As a special case, central server-based aggregation, forming a fully connected network, translates W n into a rank-1 matrix W n = 11 T /N . To minimize communication expenses, FL variants allow each agent to perform multiple SGD steps before aggregation [55, 67, 76] , resulting in a communication interval of length K and a consistent pattern W n = W for n = mK, ∀m ∈ N, and W n = I N otherwise. In particular, i) W = 11 T /N corresponds to LSGD with full agent participation (LSGD-FP) [76, 42, 51] ; ii) W is a random matrix generated by partial agent participation (LSGD-PP) [55, 18, 74] ; iii) W is generated by Metropolis-Hasting algorithm in decentralized setting, e.g., hybrid LSGD (HLSGD) [37, 32] and decentralized FL (DFL) [46, 77, 16] . We defer further discussion of W to Appendix F.1.Markovian vs i.i.d. Sampling: Agents typically employ i.i.d. or Markovian sampling, as illustrated in the bottom brown box of Figure 1 . In cases where agents have full access to their data, DSGD with i.i.d sampling has been extensively studied [60, 43, 61, 47] . In FL, many application-oriented LSGD variants have been investigated [51, 18, 77, 32, 37, 53] . However, these works solely focus on i.i.d. sampling, restricting their applicability to Markovian sampling scenarios.Markovian sampling, which has received increased attention in limited settings (see Table 1 ), is vital where agents lack independent data access. For instance, in statistical applications, agents with an unknown a priori distribution often use Markovian sampling over i.i.d. sampling [40, 63] . In HLSGD across device-to-device (D2D) networks [32, 37] , random walks reduce communication costs compared to the frequent aggregations required by Gossip algorithms [38, 28, 4] . For single-agent scenarios, vanilla SGD with Markovian noise, as applied in a D2D network, has shown improved communication efficiency and privacy [69, 28, 35] . In contrast, for agents with full data access, Markov Chain Monte Carlo (MCMC) methods can be more efficient than i.i.d. sampling, especially in high-dimensional spaces with constraints [27, 40] , where acceptance-rejection methods [12] lead to computational inefficiency (e.g., wasted samples) due to multiple rejections before obtaining a sample that satisfies constraints [26, 69] . In addition, shuffling methods can be considered as high-order Markov chains [38] , which achieves faster convergence than i.i.d. sampling [1, 79, 78] .Limitations of Non-Asymptotic Analysis on Agent's Sampling Strategy: Recent studies on the non-asymptotic behavior of DSGD and LSGD variants under Markovian sampling, as summarized in Table 1 , have made significant strides. However, these works often fall short in accurately revealing the statistical influence of each agent dynamics {X i n } on the performance of UD-SGD. For instance, [71, 68] proposed the error bound O( 1/ log 2 (1/ρ) n 1-a ), where a ∈ (0.5, 1] and ρ denotes the identical mixing rate for all agents, overlooking agent heterogeneity in sampling strategy. A similar assumption to ρ is also evident in [42] . More recent contributions from [80, 72] have attempted to relax these constraints by considering a finite-time bound of O(τ2 mix /(n + 1)), where τ mix is the mixing time of the slowest agent. This approach, however, inherently focuses on the worst-performing agent, neglecting how other agents with faster mixing rates might positively influence the system. 2 Such an analysis fails to capture the collective impact of other agents on the overall system performance, Table 1 : Comparison of recent works in distributed learning: We classify the communication patterns into seven categories, i.e., DSGD, DSGD-CT, LSGD-FP, LSGD-PP, LSGD-RC, HLSGD and DFL. We mark 'UD-SGD' when all aforementioned patterns are included and the detailed discussion on {W n } is referred to Appendix F.1. Abbreviations: 'Asym.' = 'Asymptotic', 'D.A.B' = 'Differentiating Agent Behavior', 'L.S.' = 'Linear Speedup', 'A.N.I.' = 'Asymptotic Network Independence'. agent is challenging due to privacy concerns or sporadic unreachability. Since agents in distributed learning have the freedom to choose their sampling strategies, it's vital to understand how each agent's improved sampling approach contributes to the overall convergence speed of the UD-SGD algorithm. This understanding is key to enhancing system performance, particularly in large-scale machine learning scenarios where agent heterogeneity is a defining feature.Rationale for Asymptotic Analysis: Recent trends in convergence analysis have leaned towards non-asymptotic methods, yet it's crucial to recognize the complementary role of asymptotic analysis for a better understanding of convergence behaviors, as highlighted in [9, 56, 25, 39] . For vanilla SGD, [59, 17] emphasized that central limit theorem (CLT) is far less asymptotic than it may appear under both i.i.d. and Markovian sampling. Notably, the limiting covariance matrix, a key statistical feature in vanilla SGD's CLT, also prominently features in high-probability bound [59] , explicit finite-time bound [17] and 1-Wasserstein distance in the non-asymptotic CLT [66] . [38] further underscored this by numerically showing that the limiting covariance matrix provides a more precise depiction of convergence than the mixing rates often used in finite-time upper bounds [26, 69] . Moreover, they argued that finite-time analysis may not suitably apply to certain efficient high-order Markov chains, due to the lack of comparative mixing-rate metrics.Our Contributions: We present an asymptotic analysis of the UD-SGD algorithm (1) under heterogeneous agent dynamics {X i n } and a large family of communication patterns {W n }. Specifically, • Under appropriate assumptions, all agents performing (1) asymptotically reach the consensus and find θ* : ∀i ∈ [N ], θ n ≜ 1 N N i=1 θ in denotes the average model parameter among all agents, we haveEQUATIONMoreover, we derive the CLT of UD-SGD in the form ofEQUATIONOur framework addresses technical challenges in quantifying consensus error under various communication patterns and slowly increasing communication interval. This shows a substantial extension compared to previous studies [58, 43, 51] , particularly in regulating the growth of communication stationary distribution. This naturally incurs a maximum operator, and thus convergence is strongly influenced by the slowest mixing rate, i.e., the worst-performing agent.intervals (Assumption 2.3-ii) and in proving the scaled consensus error's boundedness (Lemma B.1). Furthermore, we reformulate UD-SGD as a stochastic approximation-like iteration and tackle the Markovian noise term using the Poisson equation, a technique previously confined only to vanilla SGD with Markovian sampling [17, 38, 52] . The key here is to devise the noise decomposition that separates the consensus error among all agents from the error caused by the bias from the Markov chain, which aligns with the target distribution only asymptotically at infinity, not at finite times.• In analyzing (3), we derive the exact form of V as 1N 2 N i=1 V i .Here, V i is the limiting covariance matrix of agent i, which depends mainly on its sampling strategy {X i n }. This allows us to show that improving individual agents' sampling strategy can reduce the covariance in CLT, which in turn implies a smaller mean-square error (MSE) for large time n. This is a significant advancement over previous finite-sample bounds that only account for the worst-performing agent and do not fully capture the effect of individual agent dynamics on overall system performance. Our CLT result (3) also treats recent findings in [38] as a very special case with N = 1, where the relationship therein between the sampling efficiency of the Markov chain and the limiting covariance matrix in the CLT of vanilla SGD, can carry over to our UD-SGD.• We demonstrate that our analysis supports recent findings from studies such as [42] , which exhibited linear speedup scaling with the number of agents under LSGD-FP with Markovian sampling; and [62, 61] , which examined the notion of 'asymptotic network independence' for DSGD with i.i.d. sampling, where the convergence of the algorithm (1) at large time n depends solely on the left eigenvector of W n ( 1 N 1 considered in this paper) rather than the specific communication network topology encoded in W n , but now under Markovian sampling. We extend these findings in view of CLT to a broader range of communication patterns {W n } and general sampling strategies {X i n }. • We conduct numerical experiments using logistic regression and neural network training with several choices of agents' sampling strategies, including a recently proposed one via nonlinear Markov chain [25] . Our results uncover a key phenomenon: a handful of compliant agents adopting highly efficient sampling strategies can match or exceed the performance of the majority using moderately improved strategies. This finding is crucial for practical optimization in large-scale learning systems, moving beyond the current literature that only considers the worst-performing agent in more restrictive settings. 
preliminaries: 
 Basic Notations: We use ∥v∥ to indicate the Euclidean norm of a vector v ∈ R d and ∥M∥ to indicate the spectral norm of a matrix M ∈ R d×d . The identity matrix of dimension d is denoted by I d , and the all-one (resp. all-zero) vector of dimension N is denoted by 1 (resp. 0). Let J ≜ 11 T /N . The diagonal matrix with the entries of v on the main diagonal is written as diag(v). We also use '⪰' for Loewner ordering such that A ⪰ B is equivalent to x T (A -B)x ≥ 0 for any x ∈ R d .Asymptotic Covariance Matrix: Asymptotic variance is a widely used metric for evaluating the second-order properties of Markov chains associated with a scalar-valued test function in the MCMC literature, e.g., Chapter 6.3 [12] , and asymptotic covariance matrix is its multivariate version for a vector-valued function. Specifically, we consider a finite, irreducible, aperiodic and positive recurrent (ergodic) Markov chain {X n } n≥0 with transition matrix P and stationary distribution π, and the estimator μn (g) ≜ 1 n n-1 s=0 g(X s ) for any vector-valued function g : [N ] → R d . According to the ergodic theorem [12, 13] , we have lim n→∞ μn (g) = E π (g) a.s.. As defined in [13, 38] , the asymptotic covariance matrix Σ X (g) for a vector-valued function g(•) is given byEQUATIONwhere∆ n ≜ n-1 s=0 (g(X s ) -E π (g)).By following the algebraic manipulations in [12, Theorem 6.3.7] for asymptotic variance (univariate version), we can rewrite (4) in a matrix form such thatEQUATIONwhereG ≜ [g(1), • • • , g(N )] T ∈ R N ×d and Z ≜ [I N -P + 1π T ] -1. This matrix form explicitly shows the dependence on the transition matrix P and its stationary distribution π, and will be utilized in our Theorem 3.3.The UD-SGD in (1) can be expressed in a compact iterative form, i.e., we haveEQUATIONat each time n, where each agent i samples according to its own Markovian trajectory {X i n } n≥0 with stationary distribution π i such that E X∼πi [F i (θ, X)] = f i (θ). Let K l denote the communication interval between the (l -1)-th and l-th aggregation among N agents, and n l ≜ l m=1 K m be the time instance for the l-th aggregation. We also define τ n ≜ min l {l : n l ≥ n} as the index of the upcoming aggregation at time n such that K τn indicates the communication interval for the τ n -th aggregation, or more precisely, the length of the communication interval that includes the time index n. The communication pattern follows that W n = I n if n ̸ = n l and W n = W otherwise for l ≥ 1, where the examples of W will be discussed in Appendix F.1. Note that i) when K l = 1, (6) reduces to DSGD; ii) when K l = K > 1, (6) becomes the local SGD in FL. iii) When K l increases with l, we recover some choices of K l studied in [51] beyond LSGD-RC with i.i.d. sampling. This increasing communication interval aims to further reduce the frequency of aggregation among agents for lower communication costs, but now under a Markovian sampling setting and a wider range of communication patterns. We below state the assumptions needed for the main theoretical results. Assumption 2.1 (Regularity of the gradient). For each i ∈ [N ] and X ∈ X i , the function F i (θ, X) is L-smooth in terms of θ, i.e., for anyEQUATIONIn addition, we assume that the objective function f is twice continuously differentiable and µ-strongly convex only around the local minima θ * ∈ L, i.e.,EQUATIONAssumption 2.1 imposes the regularity conditions on the gradient ∇F i (•, X) and Hessian matrix of the objective function f (•), as is commonly assumed in [10, 45, 29, 38] . Note that (7) requires per-sample Lipschitzness of ∇F i and is stronger than the Lipschitzness of its expected version ∇f i , which is commonly assumed under i.i.d sampling setting [73, 50, 30] . However, we remark that this is in line with previous work on DSGD and LSGD-FP under Markovian sampling as well [71, 42, 80] , because ∇F i (θ, X) is no longer the unbiased stochastic version of ∇f i (θ) and the effect of {X i n } has to be taken into account in the analysis. The local strong convexity at the minimizer is commonly assumed to analyze the convergence of the algorithm under both asymptotic and non-asymptotic analysis [10, 29, 38, 45, 52, 80] . Assumption 2.2 (Ergodicity of Markovian sampling). {X i n } n≥0 is an ergodic Markov chain with stationary distribution π i such that E X∼πi [F i (θ, X)] = f i (θ), and is independent from {X j n } n≥0 , j ̸ = i. The ergodicity of the underlying Markov chains, as stated in Assumption 2.2, is commonly assumed in the literature [26, 69, 80, 42, 38] . This assumption ensures the asymptotic unbiasedness of the loss function F i (θ, •), which takes i.i.d. sampling as a special case. Assumption 2.3 (Decreasing step size and slowly increasing communication interval). i) For bounded communication interval K τn ≤ K, ∀n, we assume the polynomial step size γ n = 1/n a and a ∈ (0.5, 1]; Or ii) If K τn → ∞ as n → ∞, we assume γ n = 1/n and define η n = γ n K L+1 τn , where the sequence{K l } l≥0 satisfies n η 2 n < ∞, K τn = o(γ -1/2(L+1) n), and lim l→∞ η n l +1 /η n l+1 +1 = 1.In Assumption 2.3, the polynomial step size γ n is standard in the literature and it has the property 17, 38] . Inspired by [51] , we introduce η n to control the step size within each l-th communication interval with length K l to restrict the growth of K l . Specifically,n γ n = ∞, n γ 2 n < ∞ [n η 2 n < ∞ and K τn = o(γ -1/2(L+1) n) ensure that η n → 0 and K τn does not increase too fast in n. lim l→∞ η n l +1 /η n l+1 +1 = 1 sets the restriction on the increment from n l to n l+1 . Several practical forms of K l suggested by [51] , including K l ∼ log(l) and K l ∼ log log(l), also satisfy Assumption 2.3-ii). We defer to Appendix A the mathematical verification of these two types of K l , together with the practical implications of increasing communication interval K l . Remark 1. In Assumption 2.3, we incorporate an increasing communication interval along with a step size γ n = 1/n. This complements the choice of step size γ n in [51, Assumption 3.3] , where γ n = 1/n a for a ∈ (0.5, 1). It is important to note, however, that the increasing communication interval specified in [51, Assumption 3.2] is applicable only in i.i.d sampling. Under the Markovian sampling framework, the expression ∇F i (θ, X)-∇f i (θ) loses its unbiased and Martingale difference properties. Consequently, the Martingale CLT application as utilized by [51] does not directly extend to Markovian sampling. To address this, we adapted techniques from [29, 58] to accommodate the increasing communication interval within the Markovian sampling setting and various communication patterns. This adaptation necessitates γ n = 1/n, a specification not covered in [51] . Exploring more general forms of K l that could relax this assumption is outside the scope of our current study. Assumption 2.4 (Stability on model parameter). We assume sup n ∥θ i n ∥ < ∞ almost surely ∀i ∈ [N ]. Assumption 2.4 claims that the sequence of {θ i n } always remains in a path-dependent compact set. It is to ensure the stability of the algorithm that serves the purpose of analyzing the convergence, which is often assumed under the asymptotic analysis of vanilla SGD with Markovian noise [23, 29, 52] . As mentioned in [58, 70] , checking Assumption 2.4 is challenging and requires case-by-case analysis, even under i.i.d. sampling. Only recently the stability of SGD under Markovian sampling has been studied in [9] , but the result for UD-SGD remains unknown in the literature. Thus, we analyze each agent's sampling strategy in the asymptotic regime under this stability condition. Assumption 2.5 (Contraction property of communication matrix). i). {W n } n≥0 is independent of the sampling strategy {X i n } n≥0 for all i ∈ [N ] and is assumed to be doubly-stochastic for all n; ii). At each aggregation step n l , W n l is independently generated from some distribution P n l such that∥E W∼Pn l [W T W]-J∥ ≤ C 1 < 1 for some constant C 1 .The doubly-stochasticity of W n in Assumption 2.5-i) is widely assumed in the literature [54, 24, 43, 80] . Assumption 2.5-ii) is a contraction property to ensure that agents employing UD-SGD will asymptotically achieve the consensus, which is also common in [7, 24, 80] . Examples of W that satisfy Assumption 2.5-ii), e.g., Metropolis-Hasting matrix, partial agent participation in FL, are deferred to Appendix F.1 due to space constraint. 
asymptotic analysis of ud-sgd: 
 Almost Sure Convergence: Let θ n ≜ 1 N N i=1 θ in represent the consensus among all the agents at time n, we establish the asymptotic consensus of the local parameters θ i n , as stated in Lemma 3.1. Lemma 3.1. Under Assumptions 2.1, 2.3, 2.4 and 2.5, the consensus error θ i n -θ n diminishes to zero at the rate specified below: Almost surely, for every agent i ∈ [N ],EQUATIONLemma 3.1 indicates that all agents asymptotically reach consensus at a rate ofO(γ n ) (or O(η n )).This finding extends the scope of [58, Proposition 1], incorporating considerations for Markovian sampling, FL settings, and increasing communication interval K l . The proof, detailed in Appendix B, primarily tackles the challenge of establishing the boundedness of the sequences {γ-1 n (θ i n -θ n )} (or {η -1 n (θ i n -θ n )}) almost surely for all i ∈ [N ]. This is specifically analyzed in Lemma B.1. Next, with additional Assumption 2.2, we are able to obtain the almost sure convergence of θ n to θ * ∈ L. Theorem 3.2. Under Assumptions 2.1 -2.5, the consensus θ n converges to L almost surely, i.e.,EQUATION) Theorem 3.2 is achieved by decomposing the Markovian noise term ∇F i (θ i n , X i n ) -∇f i (θ i n ), using the Poisson equation technique as discussed in [6, 29, 17] , into a Martingale difference noise term, along with additional noise terms. We then reformulate ( 6) into an iteration akin to stochastic approximation, as depicted in (56) . The subsequent step involves verifying the conditions on these noise terms under our stated assumptions. Crucially, this theorem also establishes that UD-SGD ensures an almost sure convergence of each agent to a local minimum θ * ∈ L, even in scenarios where the communication interval K l gradually increases, in accordance with Assumption 2.3-ii). The detailed proof of this theorem is provided in Appendix C.U i ≜ Σ X i (∇F i (θ * , •))represent the asymptotic covariance matrix (defined in ( 5)) associated with each agent i ∈ [N ], given their sampling strategy {X i n } and function∇F i (θ * , •). Define U ≜ 1 N 2 N i=1 U i .We assume the polynomial step-size γ n ∼ γ ⋆ /n a , a ∈ (0.5, 1] and γ ⋆ > 0. In the case of a = 1, we further assume γ ⋆ > 1/2µ, where µ is defined in (8) . For notational simplicity, and without loss of generality, our remaining CLT result is stated while conditioning on the event that {θ n → θ * } for some θ * ∈ L. Theorem 3.3. Let Assumptions 2.1 -2.5 hold. Then,EQUATIONwhere the limiting covariance matrix V is in the form ofEQUATIONHere, we haveM = -H if a ∈ (0.5, 1), or M = I d /2γ ⋆ -H if a = 1, where H is defined in (8) .EQUATIONThe proof, presented in Appendix D, addresses the technical challenges in deriving the CLT for UD-SGD, specifically the second-order conditions in decomposing the Markovian noise term, which is not present in the i.i.d. sampling case [58, 43, 51] . We decompose∇F i (θ n , X i n )-∇f i (θ n ) into three parts in (48) using Poisson equation: e i n+1 , ν i n+1 , ξ i n+1 .The consensus error θ i n -θ n embedded in noise terms e i n+1 and ξ i n+1 is a new factor, whose characteristics have been quantified in our Lemma 3.1 but are not present in the single-agent scenario analyzed as an application of stochastic approximation in [22, 29] . The specifics of this analysis are expanded upon in Appendices D.1 to D.3. We require γ ⋆ > 1/2µ for a = 1 to ensure that the largest eigenvalue of M is negative, as this is a necessary condition for the existence of V in (12) (otherwise integration diverges). In the case where there is only one agent (N = 1), V and V ′ reduce to the matrices specified in the CLT result of vanilla SGD [29, 38, 52] . In addition, for a special case of constant communication interval in Assumption 2.3-i) and i.i.d. sampling as shown in Table 1 , we recover the CLT of LSGD-RC in [51] . See Appendix E for detailed discussions.n } for large time n, i.e., E[∥θ n -θ * ∥ 2 ]= d i=1 e T i E[(θ n -θ * )(θ n -θ * ) T ]e i ≈ γ n d i=1 e T i Ve i = γ n Tr(V), where e i is the d-dimensional vector of all zeros except 1 at the i-th entry. This indicates that a smaller limiting covariance matrix V, according to the Loewner order, results in a smaller trace of V and consequently in a reduced MSE for large n. Consideration for smaller V will be presented in the next section, where agents have the opportunity to improve their sampling strategies. Remark 2. Studies by [62, 61] have shown that in DSGD with a fixed doubly-stochastic matrix W, the influence of communication topology diminishes after a transient period. Our Theorem 3.3 extends these findings to Markovian sampling and a broader spectrum of communication patterns as in Table 1 . This extension is based on the fact that the consensus error, impacted primarily by the communication pattern, decreases faster than the CLT scale O( √ γ n ) and is thus not the dominant factor in the asymptotic regime, as suggested by Lemma 3.1. Remark 3. Recent studies have highlighted linear speedup with increasing number of agents N in the dominant term of their finite-sample error bounds under DSGD-CT with i.i.d. sampling [43] and LSGD-FC with Markovian sampling [42] . However, our Theorem 3.3 demonstrates this phenomenon under more diverse communication patterns and Markovian sampling in Table 1 via the leading term V in our CLT. Specifically, it scales with 1/N , i.e. V = V/N , where V = 1 N N i=1 V i denotes the average limiting covariance matrices across all N agents andV i = ∞ 0 e Mt U i e Mt dt, suggesting that the MSE E[∥θ n -θ * ∥ 2 ] will be improved by 1/N . A similar argument also applies to V ′ in (13), i.e., V ′ = V′ /N , where V′ = 1 N N i=1 V ′ i and V ′ i = H -1 U i H -1 .Impact of Agent's Sampling Strategy: In the literature, the mixing time-based technique has been widely used in the non-asymptotic analysis in SGD, DSGD and various LSGD variants in FL [26, 69, 68, 80, 42] , i.e., for each agent i ∈ [N ] and some constant C,EQUATIONwhere ρ i is the mixing rate of the underlying Markov chain. However, typical non-asymptotic analyses often rely on ρ ≜ max i ρ i among N agents, i.e., the worst-performing agent in their finite-time bounds [80, 72] , or assume an identical mixing rate across all N agents [42, 68] .In contrast, Remark 3 highlights that each agent holds its own limiting covariance matrices V i and V ′ i , which are predominantly governed by the matrix U i , capturing the agent's sampling strategy {X i n } and contributing equally to the overall performance of UD-SGD. For each agent i, denote by U X i and U Y i the asymptotic covariance matrices associated with two candidate sampling strategies {X i n } and {Y i n }, respectively. Let V X and V Y be the limiting covariance matrices of the distributed system in (12) , where agent i employs {X i n } and {Y i n }, respectively, while keeping other agents' sampling strategies unchanged. Then, we have the following result. Corollary 3.4. For agent i, if there exist two sampling strategies {X i n } n≥0 and {Y i n } n≥0 such thatU X i ⪰ U Y i , we have V X ⪰ V Y . Corollary 3.4 directly follows from the definition of Loewner ordering, and Loewner ordering being closed under addition (i.e., A ⪰ B implies A+C ⪰ B+C). It demonstrates that even a single agent improves its sampling strategy from {X i n } to {Y i n }, it leads to an overall reduction in V (in terms of Loewner ordering), thereby decreasing the MSE and benefiting the entire group of N agents. The subsequent question arises: How do we identify an improved sampling strategy {Y i n } over the baseline {X i n }? This question has been partially addressed by [57, 48, 38] , which qualitatively investigates the 'efficiency ordering' of two sampling strategies. In particular, [38, Theorem 3.6 (i)] shows that sampling strategy {Y n } is more efficient than {X n } if and only if Σ X (g) ⪰ Σ Y (g) for any vector-valued function g(•) ∈ R d . Consequently, in the UD-SGD framework, employing a more efficient sampling strategy {Y i n } over the baseline {X i n } by agent i leads toΣ X i (∇F i (θ * , •)) ⪰ Σ Y i (∇F i (θ * , •)), thus satisfying U X i ⪰ U Y i .This finding, as per Corollary 3.4, implies an overall improvement in UD-SGD.For illustration purposes, we list a few examples where two competing sampling strategies follow efficiency ordering: i) When an agent has complete access to the entire dataset (e.g., deep learning), shuffling techniques like single shuffling and random reshuffling are more efficient than i.i.d. sampling [38, 78] ; ii) When an agent works with a graph-like data structure and employs a random walk, e.g., agent i in Figure 1 , using non-backtracking random walk (NBRW) is more efficient than simple random walk (SRW) [48] . iii) A recently proposed self-repellent random walk (SRRW) is shown to achieve near-zero sampling variance, indicating even higher sampling efficiency than NBRW and SRW [25] . 3 This random-walk-based sampling finds a particular application in large-scale FL within D2D networks (e.g., mobile networks, wireless sensor networks), where each agent acts as an edge server or access point, gathering information from the local D2D network [37, 32] . Employing a random walk over local D2D network for each agent constitutes the sampling strategy. Theorem 3.3 and Corollary 3.4 not only qualitatively compare these sampling strategies but also allow for a quantitative assessment of the overall system enhancement. Since every agent contributes equally to the limiting covariance matrix V of the distributed system as in Remark 3, a key application scenario is to encourage a subset of compliant agents to adopt highly efficient strategies like SRRW, potentially yielding better performance than universally upgrading to slightly improved strategies like NBRW. This approach, more feasible and impactful in large-scale machine learning scenarios where some agents cannot freely modify their sampling strategies, is a unique aspect of our framework not addressed in previous works focusing on the worst-performing agent [80, 42, 68, 72] . 
experiments: 
 In this section, we empirically evaluate the effect of agents' sampling strategies under various communication patterns in UD-SGD. We consider the L 2 -regularized binary classification problemEQUATION)where the feature vector x i,j and its corresponding label y i,j are held by agent i, with a penalty parameter κ set to 1. We use the ijcnn1 dataset [14] with 22 features in each data point and 50k data points in total, which is evenly distributed to two groups with 50 agents each (N = 100 agents NBRW or SRRW with reweighting to sample its local dataset with uniform weight. In this simulation, we assume that agents can only communicate through a communication network using the DSGD algorithm. This scenario with heterogeneous agents, as depicted in Figure 1 , is of great interest in large-scale machine learning [37, 32] . In addition, we employ a decreasing step size γ n = 1/n in our UD-SGD framework (1) because it is typically used for the strongly convex objective function and is tested to have the fastest convergence in this simulation setup. Due to space constraints, we defer detailed simulation setup, including the introduction of SRW, NBRW, and SRRW, to Appendix G.1.The simulation results are obtained through 120 independent trials. In Figure 2 (a), we assume that the first group of agents perform either i.i.d. sampling or shuffling method, while the other group of agents all change their sampling strategies from baseline SRW to NBRW and SRRW, as shown in the legend. This plot shows that improved sampling strategy leads to overall convergence speedup since NBRW and SRRW are more efficient than SRW [38, 25] . Furthermore, it illustrates that SRRW is significantly more efficient than NBRW in this simulation setup, i.e., SRRW ≫ NBRW > SRW in terms of sampling efficiency. While keeping the second group of agents unchanged, we can see that shuffling method outperforms i.i.d. sampling with smaller asymptotic MSE. However, shuffling method may not perform perfectly for small time n due to slow mixing behavior in the initial period, which is also observed in the single-agent scenario in [65, 1, 38] . The error bar therein also indicates that the random-walk sampling strategy has a significant impact on the overall system performance and SRRW has smaller variance than NBRW and SRW.In Figure 2 (b), we let the first group of agents perform i.i.d. sampling while only changing a portion of agents in the second group to upgrade from SRW to SRRW, e.g., 30 SRW 20 SRRW in the legend means that there are 30 agents using SRW while the rest 20 agents in the second group upgrade to SRRW. We observe that more agents willing to upgrade from SRW to SRRW lead to smaller asymptotic MSE, as predicted by Theorem 3.3 and Remark 3. This improvement in MSE reduction doesn't scale linearly with more agents adopting SRRW because each agent holds its own dataset that are not necessarily identical, resulting in different individual limiting covariance matricesV i ̸ = V j .While maintaining i.i.d. sampling for the first group of agents, we compare the performance when the second group of agents in Figure 2 (c) employ NBRW or SRRW. Remarkably, the case with only 10 agents out of 50 agents in the second group adopting far more efficient sampling strategy (40 SRW, 10 SRRW) through incentives or compliance already produces a smaller MSE than all 50 agents using slightly better strategy (50 NBRW). The performance gap becomes even more pronounced when 20 agents upgrade from SRW to SRRW (30 SRW, 20 SRRW). We show that the performance of a distributed system can be improved significantly when a small proportion of agents adopt highly efficient sampling strategies. Furthermore, in Appendix G.2, we simulate an image classification task with CIFAR-10 dataset [44] by training a 5-layer CNN and ResNet-18 model collaboratively through a 10-agent network. The result is illustrated in Figure 3 , where SRRW outperforms NBRW and SRW as expected. In summary, we find that upgrading even a small portion of agents to efficient sampling strategies (e.g., shuffling method, NBRW, SRRW under different dataset structures) improves system performance in UD-SGD. These results are consistent in binary and image classification tasks, underscoring that every agent matters in distributed learning. 
conclusion: 
 In this work, we develop an UD-SGD framework that establishes the CLT of various distributed algorithms with Markovian sampling. We overcome technical challenges such as quantifying consensus error under very general communication patterns and decomposing Markovian noise through the Poisson equation, which extends the analysis beyond the single-agent scenario. We demonstrate that even if only a few agents optimize their sampling strategies, the entire distributed system will benefit with a smaller limiting covariance in the CLT, suggesting a reduced MSE. This finding challenges the current established upper bounds where the worst-performing agent leads the pack. Future studies could pivot towards developing fine-grained finite-time bounds to individually characterize each agent's behavior, and theoretically analyze the effect of SRRW in UD-SGD.[79] Chulhee Yun, Suvrit Sra, and Ali Jadbabaie. Open problem: Can single-shuffle sgd be better than reshuffling sgd and gd? In Proceedings of Thirty Fourth Conference on Learning Theory, volume 134, pages 4653-4658. PMLR, Aug 2021.[80] Sihan Zeng, Thinh T Doan, and Justin Romberg. Finite-time convergence rates of decentralized stochastic approximation with applications in multi-agent and multi-task learning. IEEE Transactions on Automatic Control, 2022.[81] Yue Zhao, Meng Li, Liangzhen Lai, Naveen Suda, Damon Civin, and Vikas Chandra. Federated learning with non-iid data. arXiv preprint arXiv:1806.00582, 2018. 
 
 Appendix: A.1 Suitable choices of K l When wet let K l ∼ log(l) (resp. K l ∼ log log(l)), as suggested by [51] , it trivially satisfiesK τn = o(γ -1/2(L+1) n ) = o(n 1/2(L+1) ) since by definition K τn < K n ∼ log(n) (resp. log log(n)), and log(n) = o(n ϵ ) (resp. log log(n) = o(n ϵ )) for any ϵ > 0. Besides, n η 2 n = n γ 2 n K 2(L+1) τn ≲ n n -2 n 2(L+1)ϵ = n n 2(L+1)ϵ-2 .To ensure n η 2 n < ∞, it is sufficient to have 2(L + 1)ϵ -2 < -1, or equivalently, ϵ < 1/2(L + 1). Since ϵ can be arbitrarily small to satisfy the condition, n η 2 n < ∞ is satisfied. When K l ∼ log(l), we can rewrite the last condition asEQUATIONwhere we have n l ∼ log(l!) such that K l+1 /n l = log(l + 1)/ log(l!) → 0 and log(l + 1)/ log(l) → 1 as l → ∞, which leads to lim n→∞ η n l +1 /η n l+1 +1 = 1. Similarly, for K l ∼ log log(l), we have n l ∼ log( l s=1 log(s)) such that K l+1 /n l ∼ log log(l + 1)/ log log( l s=1 log(s)) → 0 and log log(l + 1)/ log log(l) → 1 as l → ∞, which also leads to lim n→∞ η n l +1 /η n l+1 +1 = 1.In this assumption, we allow the number of local iterations to go to infinity asymptotically. In distributed learning environments such as mobile, IoT, and wireless sensor networks, where nodes are often constrained by battery life, increasing communication interval in Assumption 2.3-ii) plays a crucial role in balancing energy costs with communication effectiveness. It allows agents to communicate more frequently early on, leading to a faster initial convergence to the neighborhood of θ * . Then, we slow down the communication frequency between agents to conserve energy, leveraging the diminishing returns on accuracy improvements from additional communications.Consider the scenario where devices across multiple clusters collaborate on a distributed optimization task, utilizing local datasets. Devices within each cluster form a communication network that allows a virtual agent to perform a heterogeneous Markov chain trajectory via random walk, or an i.i.d. sequence in a complete graph with self-loops, depending on the application context. Each cluster features an edge server that supports the exchange of model estimates with neighboring clusters. By performing K local updates before uploading these to the cluster's edge server, the model benefits from reduced communication overhead. As the frequency of updates between devices and edge servers decreases -optimized by gradually increasing K -we effectively lower communication costs, particularly as the model estimation θ n is close to θ * . B Proof of Lemma 3.1Let J ⊥ ≜ I N -J ∈ R N ×N and J ⊥ ≜ J ⊥ ⊗ I d ∈ R N d×N d , where ⊗ is the Kronecker product. Let Θ n = [(θ 1 n ) T , • • • , (θ N n ) T ] T ∈ R N d .Then, motivated by [58] , we define a sequenceϕ n ≜ η -1 n+1 J ⊥ Θ n ∈ R N d in the increasing communication interval case (resp. ϕ n ≜ γ -1 n+1 J ⊥ Θ n in the bounded communication interval case), where η n+1 is defined in Assumption 2.3-ii). J ⊥ Θ n = Θ n -1N (11 T ⊗ I d )Θ n represents the consensus error of the model. We first give the following lemma that shows the pathwise boundedness of ϕ n . Lemma B.1. Let Assumptions 2.1, 2.3, 2.4 and 2.5 hold. For any compact setΩ ⊂ R N d , the sequence ϕ n satisfies sup n E[∥ϕ n ∥ 2 1 ∩ j≤n-1 {Θj ∈Ω} ] < ∞.n ≥ 0, E[∥J ⊥ Θ n ∥ 2 1 ∩ j≤n-1 {Θj ∈Ω} ] = η 2 n+1 E[∥ϕ n ∥ 2 1 ∩ j≤n-1 {Θj ∈Ω} ] ≤ Cη 2n+1 for some constant C that depends on C 1 and Ω. Along with Assumption 2.4 such that ∥J ⊥ Θ n ∥ is always bounded per each trajectory, it means∥J ⊥ Θ n ∥1 ∩ j≤n-1 {Θj ∈Ω} = O(η n ) a.s.Let {Ω m } m≥0 be a sequence of increasing compact subset of R N d such that m Ω m = R N d . Then, we know that for any m ≥ 0,EQUATION(17) indicates either one of the following two cases:• there exists some trajectory-dependent index m ′ such that each trajectory {Θ n } n≥0 is always within the compact set Ω m ′ , i.e., 1 ∩ j≤n {Θj ∈Ω m ′ } = 1 (satisfied by the construction of increasing compact sets {Ω m } m≥0 and Assumption 2.4), and we have∥J ⊥ Θ n ∥ = O(η n ) such that lim n→∞ J ⊥ Θ n = 0;• Θ n will escape the compact set Ω m eventually for any m ≥ 0 in finite time such that 1 ∩ j≤n-1 {Θj ∈Ωm} = 0 when n is large enough.We can see the second case contradicts Assumption 2.4 because we assume every trajectory {Θ n } n≥0 is within some compact set. Therefore, (17) for any m ≥ 0 is equivalent to showing ∥J ⊥ Θ n ∥ = O(η n ) and lim n→∞ J ⊥ Θ n = 0. Under Assumption 2.3-i) we can obtain similar result ∥J ⊥ Θ n ∥ = O(γ n ) by following the same steps as above, which completes the proof of Lemma 3.1.Proof of Lemma B.1. We begin by rewriting (6) in the matrix form,EQUATIONwhere X n ≜ (X 1 n , X 2 n , • • • , X N n ) and ∇F(Θ n , X n ) ≜ [∇F 1 (θ 1 n , X 1 n ) T , • • • , ∇F N (θ N n , X N n ) T ] T ∈ R N d . Recall θ n ≜ 1 N N i=1 θ i n ∈ R d and we have [θ T n , • • • , θ T n ] T = 1 N (11 T ⊗ I d )Θ n ∈ R N d .EQUATIONwhere the equality comes from 1 N (11 18) and ( 19), we haveT ⊗ I d )W n = 1 N (11 T W n ⊗ I d ) = 1 N (11 T ⊗ I d ). With (EQUATIONwhere the second equality comes from W n -1 N (11EQUATION)Define the filtration{F n } n≥0 as F n ≜ σ{Θ 0 , X 0 , W 0 , Θ 1 , X 1 , W 1 , • • • , X n-1 , W n-1 , Θ n , X n }.Recursively computing (21) w.r.t the time interval [n l , n l+1 ] givesEQUATIONwhere is the backward multiplier, the second equality comes fromEQUATIONwhere the first inequality comes fromJ T ⊥ J ⊥ = J ⊥ and η k+2 ≥ η n l+1 +1 for k ∈ [n l , n l+1 -1]. Then, we analyze the norm of the gradient ∥∇F(Θ k , X k )∥ in the second term on the RHS of ( 23) conditioned on F n l . By Assumption 2.4, we assume Θ n l is within some compact set Ω at timen l such that sup i∈[N ],X i ∈Xi ∇F i (θ i n l , X i ) ≤ C Ω for some constant C Ω . For n = n l + 1 and any X ∈ X 1 × X 2 × • • • × X N , we have ∥∇F(Θ n l +1 , X)∥ ≤ ∥∇F(Θ n l +1 , X) -∇F(Θ n l , X)∥ + ∥∇F(Θ n l , X)∥. Considering ∥∇F(Θ n l , X)∥, we have sup X ∥∇F(Θ n l , X)∥ 2 ≤ N i=1 sup X i ∈Xi ∥∇F i (θ i n l , X i )∥ 2 ≤ N C 2 Ω such that ∥∇F(Θ n l , X)∥ ≤ √ N C Ω .In addition, we haveEQUATION)For n = n l + 2 and any X, we have∥∇F(Θ n l +2 , X)∥ ≤ ∥∇F(Θ n l +2 , X) -∇F(Θ n l +1 , X)∥ + ∥∇F(Θ n l +1 , X)∥.Similar to the steps in ( 24), we haveEQUATION)Then, ∥∇F(Θ n l +2 , X)∥ ≤ (1 + γ n l +2 L) sup X ∥∇F(Θ n l +1 , X)∥ and, together with (25), we haveEQUATIONBy induction,∥∇F(Θ n l +m , X)∥ ≤ m s=1 (1 + γ n l +s L) √ N C Ω for m ∈ [1, K l+1 -1].The next step is to analyze the growth rate of EQUATIONwhere the last inequality comes fromK l+1 -1 m=0 m L < K l+1 (K l+1 -1) L < K L+1 l+1 .We can see the sum of the norm of the gradients are bounded by √ N e L C Ω , which only depends on the compact set Ω at time n = n l .Let δ 1 ∈ (C 1 , 1). Since from Assumption 2.3-ii), lim l→∞ η n l +1 /η n l+1 +1 = 1, there exists some large enough l 0 such that ( ηn l +1 ηn l+1 +1 ) 2 C 1 < δ 1 < δ 2 := (δ 1 + 1)/2 < 1 for any l > l 0 . Note that δ 1 depends only on C 1 and is independent of F n . Then, let CΩ := √ N e L C Ω , we can rewrite (23) asEQUATIONwhereM Ω satisfies M Ω > 8 C2 Ω /(1 -δ 1 ) + δ 1 C2Ω , which is derived from rearranging (29) asM Ω ≥ (δ 1 -δ 2 )∥ϕ n l ∥ 2 + 2δ 1 CΩ ∥ϕ n l ∥ + δ 1 C2Ω and upper bounding the RHS. Upon noting that 1 ∩ j≤n l {Θj ∈Ω} ≤ 1 ∩ j≤n l-1 {Θj ∈Ω} , we obtainEQUATIONThe induction leads toE[∥ϕ n l+1 ∥ 2 1 ∩ j≤n l {Θj ∈Ω} ] ≤ δ n l+1 -n l 0 2 E[∥ϕ n l 0 ∥ 2 1 ∩ j≤n l 0 -1 {Θj ∈Ω} ] + M/(1 -δ 2 ) < ∞ for any l ≥ l 0 .Besides, for m ∈ (n l , n l+1 ), by following the above steps (23) applied to (21) , we haveEQUATIONBy (28) we already show that ∥n l+1 -1 k=n l ∇F(Θ k ,X k ) K L+1 l+1 ∥ < ∞ conditioned on F n l . Therefore, E[∥ϕ m ∥ 2 1 ∩ j≤n l {Θj ∈Ω} ] < ∞ for m ∈ (n l , n l+1). This completes the boundedness analysis ofE[∥ϕ n ∥ 2 1 ∩ j≤n-1 {Θj ∈Ω} ].Case 2 (Bounded communication interval K τn ≤ K): In this case, we do not need the auxiliary step size η n and can directly work on γ n = 1/n a for a ∈ (0.5, 1]. Similar to (20) , we haveEQUATIONand let b n ≜ γ n /γ n+1 , dividing both sides of above equation by γ n+2 givesEQUATIONThen, by following the similar steps in ( 22) and ( 23), we obtainEQUATIONAlso similar to ( 25) -( 28), we can bound the sum of the norm of the gradients asEQUATIONNow that K l is bounded above by K,k s=n l (1 + γ s+1 L) ≤ e L k s=n lγs+1 < e L K-1 s=0 γs+1 := C K . Then, we further bound (35) asEQUATIONThe subsequent proof is basically a replication of ( 29) - (31) and is therefore omitted.We focus on analyzing the convergence property of θ, which is obtained by left multiplying (18) withEQUATIONwhere the second equality comes from W n being doubly stochastic and 1 N (1T ⊗ I d )W n = 1 N (1 T W n ⊗ I d ) = 1 N (1 T ⊗ I d ).For self-contained purpose, we first give the almost sure convergence result for the stochastic approximation that will be used in our proof.Theorem C.1 (Theorem 2 [23] ). Consider the stochastic approximation in the form ofEQUATIONAssume that C1. w.p.1, the closure of {θ n } n≥0 is a compact subset of R d ; C2. {γ n } is a decreasing sequence of positive number such that n γ n = ∞; C3. w.p.1, lim p→∞ p n=1 γ n e n exists and is finite. Moreover, lim n→∞ r n = 0.C4. vector-valued function h is continuous on R d and there exists a continuously differentiable functionV : R d → R such that ⟨∇V (θ), h(θ)⟩ ≤ 0 for all θ ∈ R d . Besides, the interior of V (L) is empty where L ≜ {θ ∈ R d : ⟨∇V (θ), h(θ)⟩ = 0}.Then, w.p.1, lim sup n d(θ n , L) = 0.We can rewrite (37) asEQUATIONand work on the converging behavior of the third and fourth term. By definition of function ∇f (•), we haveEQUATIONBy the Lipschitz continuity of function ∇F i (•, X) in ( 7), we haveEQUATIONwhere the second inequality comes from the Cauchy-Schwartz inequality. In Appendix B, we have shown lim n J ⊥ Θ n = 0 almost surely such that lim n→∞ r n = 0 almost surely.Next, we further decompose the fourth term in (39) . For an ergodic transition matrix P and a function v associated with the same state space X , define the operator P k v(x) ≜ y∈X P k (x, y)v(y) for the k-step transition probability P k (x, y). Denote by P 1 , • • • , P N the underlying transition matrices of all N agents with corresponding stationary distribution π 1 , • • • , π N . Then, for every functionEQUATIONThe solution of the Poisson equation ( 42) has been studied in the literature, e.g., [17, 38] . For self-contained purpose, we derive the closed-form m θ i (x) from scratch. First of all, we can obtain function m θ i (x) in the recursive form as follows,EQUATIONIt is not hard to check that (43) satisfies (42) . Note that by induction we getEQUATIONUsing the Lipschitzness property of P i m θ (X) in Lemma C.2, we haveEQUATIONIn Appendix B we have shownlim n→∞ J ⊥ Θ n = 0 almost surely. Moreover, ∥ 1 N (1 T ⊗ I d )∇F(Θ n , X n )∥ is bounded per sample path. Therefore, lim n→∞ 1 N N i=1 ∥ξ i n+1 ∥ = 0 such that lim n→∞ 1 N N i=1 ξ i n+1 = 0 almost surely.To sum up, we decompose (39) intoEQUATIONNow that lim p→∞ p n=11 N N i=1 γ n e i n and lim p→∞ p n=0 1 N N i=1 γ n+1 (ν i n -ν i n+1) converge and are finite, lim n→∞ r n = 0, lim n→∞ 1 N N i=1 ξ i n = 0 , all the conditions of C3 in Theorem C.1 are satisfied. Additionally, Assumption 2.4 corresponds to C1, Assumption 2.3 meets C2, and C4 is automatically satisfied when we choose the lyapunov functionV (θ) = f (θ). Therefore, lim sup n inf θ * ∈L ∥θ n -θ * ∥ = 0.To obtain Theorem 3.3, we need to utilize the existing CLT result for general SA in Theorem D.1 and check all the necessary conditions therein. Theorem D.1 (Theorem 2.1 [29] ). Consider the stochastic approximation iteration (38), assume C1. Let θ * be the root of function h, i.e., h(θ * ) = 0, and assume lim n→∞ θ n = θ * . Moreover, assume the mean field h is twice continuously differentiable in a neighborhood of θ * , and the Jacobian H ≜ ∇h(θ * ) is Hurwitz, i.e., the largest real part of its eigenvalues B < 0;C2. The step size n γ n = ∞, n γ 2 n < ∞, and either (i ). log(γ n-1 /γ n ) = o(γ n ), or (ii). log(γ n-1 /γ n ) ∼ γ n /γ ⋆ for some γ ⋆ > 1/2|B|; C3. sup n ∥θ i n ∥ < ∞EQUATIONC5. Let r n = r(1)EQUATIONn , r n is F n -adapted, andEQUATIONThen,EQUATIONwhereEQUATIONNote that the matrix U in the condition C4(b) of Theorem D.1 was assumed to be positive definite in the original Theorem 2.1 [29] . It was only to ensure that the solution V to the Lyapunov equation ( 60) is positive definite, which was only used for the stability of the related autonomous linear ODE (e.g., Theorem 3.16 [15] or Theorem 2.2.3 [36] ). However, in this paper, we do not need strict positive definite matrix V. Therefore, we extend U to be positive semi-definite such that V is also positive semi-definite (see Lemma D.2 for the closed form of matrix V). Such kind of extension does not change any of the proof steps in [29] .Our Assumption 2.1 corresponds to C1 by letting function h(θ) = -∇f (θ) therein. We can also let γ ⋆ in Theorem 3.3 large enough to satisfy C2. The typical form of step size, also indicated in [29] , is polynomial step size γ n ∼ γ ⋆ /n a for a ∈ (0.5, 1]. Note that a ∈ (0.5, 1) satisfies C2 (i) and a = 1 satisfies C2 (ii). Assumption 2.4 corresponds to C3. 4To check condition C4, we need to analyze the Martingale difference sequence {e i n }. Recalle i n+1 = m θ i n (X i n ) -P i m θ i n (X i n-1) such that there exists a constant C,EQUATIONSince ∥m θ i n (Y )∥ < ∞ almost surely by Assumption 2.4 and X i is a finite state space, at all time n, we haveEQUATION)and there exists another constant C ′ such that by definition of P i m θ i n (X i n-1 ), we haveEQUATIONTherefore, E[∥e i n+1 ∥ 2+τ |F n ] < ∞ a.s. for all n and C4.(a) is satisfied.We now turn to C4.(b). Note that for any i ̸ = j, we have E[e i n+1 (e j n+1 ) T |F n ] = 0 due to the independence between agent i and j, and E[e i n+1 |F n ] = 0. Then, we haveT |F n ] = E[e i n+1 |F n ] • E[(e j n+1 )EQUATIONThe analysis of E[e i n+1 (e i n+1 ) T |F n ] is inspired by Section 4 [29] and Section 4.3.3 [22] , where they constructed another Poisson equation to further decompose the noise terms therein. 5 Here, expanding By the definition (66) and Assumption 2.4, for a sample path, sup n ∥G i (θ i n , x)∥ < ∞ for any x ∈ X i , as well as sup n ∥g i (θ i n )∥ < ∞, which leads to sup n ∥φ i θ i n (x)∥ < ∞ for any x ∈ X i because of (68) .Then, we have sup n ∥D EQUATIONwhere∥φ i θ i 0 (X i 0 )∥ + ∥φ i θ i n (X i n )∥ ≤ D 1for a given sample path, D 2 is the Lipschitz constant of φ i θ i (x), and ∥∇F i (x i , X i )∥ ≤ D Ω for any x i ∈ Ω and X i ∈ X i . Then,γ n n k=1 D (2,b) i,k ≤ D 2 D Ω γ n n k=1 γ k + γ n D 1 → 0 as n → ∞ (75) because γ n n k=1 γ k = O(n 1-2a) by assumption 2.3. Therefore, the second condition of C4 is satisfied.We now analyze condition C5. The decreasing rate of each term in (56) has been proved in Appendix C. Specifically, by assumption 2.4, there exists a compact subset for a given sample path, and• we have shown that r(A) n = O(η n ) a.s., which implies r (A) n = o( √ γ n ) a.s. • For 1 N N i=1 ξ i n , in the case of increasing communication interval, 1 N N i=1 ξ i n = O(γ n + η n ), by Assumption 2.3-ii), we know (γ n + η n )/ √ γ n = √ γ n + √ γ n K L+1 τn = o(1) such that ∥ 1 N N i=1 ξ i n ∥ = o( √ γ n ) almost surely.On the other hand, in the case of bounded communication interval, 1N N i=1 ξ i n = O(γ n ) such that ∥ 1 N N i=1 ξ i n ∥ = o( √ γ n ) a.s. • Since sup n ∥ν i n ∥ < ∞ almost surely, we have sup p ∥ 1 N N i=1 p k=0 (ν i k -ν i k+1 )∥ = sup p ∥ 1 N N i=1 (ν i 0 -ν i p+1 )∥ < ∞ almost surely. Then, √ γ p ∥ 1 N N i=1 p k=0 (ν i k - ν i k+1 )∥ = O( √ γ p ) leads to √ γ p ∥ 1 N N i=1 p k=0 (ν i k -ν i k+1 )∥ = o(1) a.s. Let r (1) n ≜ r (A) n + 1 N N i=1 ξ i n and r (2) n ≜ 1 N N i=1 (ν i k -ν i k+1 ).From above, we can see that C5 in Theorem D.1 is satisfied and we show that all the conditions in Theorem D.1 have been satisfied.Lastly, we need to analyze the closed-form expression of U as in C4 (b) of Theorem D.1. Recall that U = 1 N 2 N i=1 U i and U i = g i (θ * ) in (69) . We now give the exact form of function g i (θ * ) as follows:EQUATIONwhere the second equality comes from the recursive form of m θ i (x) in (45) , and that the process {X n } n≥0 is in its stationary regime, i.e., X 0 ∼ π i from the beginning. The last equality comes from rewriting Cov(∇F i (θ * , X i ), ∇F i (θ * , X j )) in a matrix form. Note that g i (θ * ) is exactly the asymptotic covariance matrix of the underlying Markov chain {X i n } n≥0 associated with the test function ∇F i (θ * , •). By utilizing the following lemma, we can obtain the explicit form of V as defined in (60) . Lemma D.2 (Lemma D.2.2 [41] ). If all the eigenvalues of matrix M have negative real part, then for every positive semi-definite matrix U there exists a unique positive semi-definite matrix V satisfying U + MV + VM T = 0. The explicit solution V is given asEQUATIONWe now consider the CLT result of Polyak-Ruppert averaging θn = 1 n n-1 k=0 θ k . The steps follow similar way by verifying that the conditions in the related CLT of Polyak-Ruppert averaging for the stochastic approximation are satisfied. The additional assumption is given below. agent holds a disjoint local dataset (non-overlapping data points for every agent), while we distribute the ijcnn1 dataset [14] with more varied distribution among 100 agents by leveraging Dirichlet distribution with the default alpha value of 0.5 in Figure 2 (d) -2(f).Moreover, we assume that all agents are distributed over a communication network. In order to create this network among 100 agents and the graph-like dataset structure held by each agent, we utilize connected sub-graphs from the real-world graph Facebook in SNAP [49] . All 100 agents collaborate together to generate a deterministic communication matrix W = [W ij ] with Metropolis Hasting algorithm of the following form: For i, j ∈ [N ], we haveW ij = min 1 di , 1 dj if agent j is the neighbor of agent i, 0 otherwise, W ii = 1 - j∈[N ] W ij ,where d i represents the degree of agent i in the graph. The communication interval K is set to 1, as is the usual choice in DSGD [71, 80, 61, 68] .For the first group of agents, we assume they have full access to their datasets, thus performing i.i.d. sampling or single shuffling. In particular,• i.i.d. sampling employed by agent i means that the data point X i n is independently and uniformly sampled from its dataset X i at each time n.• Single shuffling, by its name, only shuffles the dataset once and adheres to that specific order throughout the training process.On the other hand, within the second group of agents, we assume that they hold graph-like datasets. Now, we introduce simple random walk (SRW), non-backtracking random walk (NBRW), and self-repellent random walk (SRRW) in order:• SRW refers to the walker that chooses one of the neighboring nodes uniformly at random.• NBRW, as studied in [2, 48, 5] , is a variation of SRW, which selects one of the neighbors uniformly at random, with the exception of the one visited in the last step. • SRRW, recently proposed by [25] , is designed with a nonlinear transition kernel K[x] ∈ [0, 1] N ×N of the following form:EQUATIONwhere matrix P = [P ij ] is the transition kernel of the baseline Markov chain and µ = [µ i ] is its corresponding stationary distribution. Additionally, α denotes the force of self repellence, and larger α leads to stronger force of self repellence, thus higher sampling efficiency [25, Corollary 4.3]. Moreover, vector x ∈ R N is in the interior of probability simplex, representing the empirical distribution, in other words, the visit frequency of each node in the graph. The update rule of this empirical distribution is in the following form:EQUATIONwhere β n ≜ (n + 1) -b is the step size of SRRW iterates. b = 1 was original proposed in [25] and is recently extended to b ∈ (0.5, 1) in [39] . In this simulation, we use SRW as the baseline Markov chain of SRRW, and in turn µ is proportional to the degree distribution. We also assume x 0 = 1/N , i.e., each node has been visited once, and choose the step size β n = (n + 1) -0.8 , force of self repellence α = 20 according to the suggestion in [39, Section 4] .Since SRW, NBRW, and SRRW all admit the stationary distribution that is proportional to degree distribution, in order to obtain unfirom target in (15) , we need to reweight the gradient computed by each agent i in order to maintain an asymptotic unbiased gradient. Thus, agent i should modify its SGD update from (1a) to the following: EQUATIONIn this part, we perform the image classification task through a 5-layer neural network, where the CIFAR10 dataset [44] with 50k image data is evenly distributed to 10 agents. Each agent possesses 5k images, which are further divided into 200 batches, each batch with 25 images.The Convolutional Neural Network (CNN) model used in this simulation encompasses:• Two convolutional layers (i.e., nn.Conv2d (3, 6, 5) and nn.Conv2d(6, 16, 5)), each followed by ReLU activation functions to introduce non-linearity and max pooling (i.e., nn.MaxPool2d(2, 2)) to reduce spatial dimensions.• Three fully connected (linear) layers, concluding with a softmax output to handle the multi-class classification problem.Similar to the simulation setup in Section 4, among the 10 participating agents, five have unrestricted access to their respective data allocations, enabling them to utilize the shuffling method to iterate through their batches. The other five agents are designed to simulate limited data access scenarios. Their data access is structured using five distinct graph topologies extracted from the SNAP dataset collection [49] , each graph simulating a unique communication pattern among the batches (nodes) of data. Within these topologies, agents adopt one of three random walk strategies -SRW, NBRW, and SRRW, all with importance reweighting -to sample the batches for training.Local model training is conducted for five epochs at each agent before aggregating the updates at a central server -a process repeated for a total of 200 communication rounds. Each epoch consists of a full traversal of the local dataset of agents in the first group, or 200 batches sampled for training in the second group of agents. To mimic realistic conditions, we also introduce partial agent participation where only 40% of agents are selected randomly to transmit their updates in each round, reflecting the intermittent communication in real-world FL deployments. Lastly, the selection of the step size β n for SRRW iterates (90) is a critical aspect of our experiments. In this simulation, we experiment with various values of b ∈ {0.501, 0.6, 0.7, 0.8, 0.9} to determine the most advantageous setting for maximizing the benefits of the SRRW strategy. Based on our findings, the best choice for the SRRW step size is b = 0.501, in other words, β n = (n + 1) -0.501 .The simulation result is quantified by averaging the training loss across ten repeated trials for each configuration. As depicted in Figure 3 (a), the training results are consistent with our previous findings in Figure 2 (a) in the context of the FL framework and the training of neural networks: the use of a more effective sampling approach, even for a portion of the agents, results in significant enhancements in the overall training of the model, and this improvement is further enhanced through the highly efficient sampling strategy SRRW.In Figure 3 (b) and 3(c), we perform image classification experiments in the FL setting with partial client participation. Only 4 random agents will participate in the training process at each aggregation phase. In Figure 3 (b), we fix the sampling strategy (shuffling, SRRW with α = 10) and test the linear speedup effect for the 5-layer CNN model by duplicating 10 agents to N agents with N ∈ {10, 20, 30, 40}, keeping the same participation ratio 0.4. As can be seen from the plot, the training loss is inversely proportional to the number of agents, i.e., at 200 rounds, the training loss is 0.52 for 10 agents, 0.23 for 20 agents, 0.18 for 30 agents, and 0.12 for 40 agents. In Figure 3 (c), we extend the current simulation from 5-layer CNN model to ResNet-18 model [34] in order to numerically test the performance of different sampling strategies in a more complex neural network training task. By fixing the shuffling in the first group of agents, we observe that improving Markovian sampling from SRW to NBRW, then to SRRW, gives accelerated training process with smaller training loss.Our study provides crucial insights into the identification of nuanced agents' sampling behaviors in UD-SGD, where improving each agent's sampling strategy speeds up overall system performance without additional computational burden except the additional storage for the visit counts used for sampling their datasets. Our UD-SGD is scalable in terms of larger datasets as the sampling strategy (i.e., random walk) utilized by each agent leverages only local information for its dataset. However, our work has two limitations that should be acknowledged.1. Assumption 2.4 posits that the parameter trajectory {θ n } is almost surely bounded, which is a strong assumption. This is crucial for guaranteeing the well-defined nature of all related quantities. Some mechanisms such as projections onto a compact subset [45, Chapter 5.1], or truncation-based method with expanding compact subsets can do the trick to ensure that the iteration is always bounded [3] . As mentioned in our discussion after Assumption 2.4, only recently the stability of SGD under Markovian sampling has been guaranteed to hold for some class of objective function f [9] , while the discussion on stability issue under multi-agent scenario with Markovian sampling remains an open problem and we do not pursue to remove this stability assumption in this paper. 2. Asymptotic analysis: The main results of our work, i.e., almost sure convergence and central limit theorem in distributed optimization, are based on asymptotic analysis and might not accurately represent the finite-sample performance of each contributing agent in the system. The state-of-the-art finite-sample analysis in the literature only focuses on the worst-performing agent that cannot capture the nuanced differences between agents' sampling strategies, with the explanation detailed in Footnote 2. Thus, a finite-sample error bound that distinguish every agent's dynamics is still unknown and regarded as a future direction.Throughout the paper we don't impose convexity assumption on f (θ). For convex f (θ), L is the global minima. For non-convex f (θ), L represents the collection of local minima, which is of great interest in neural network training for sufficiently good performance[20,19]. With an additional condition such as the Polyak-Lojasiewicz inequality, non-convex f (θ) is ensured to have a unique minima[1,75,78].Although improving the finite-time upper bound to distinguish each agent may not be the focus of the aforementioned works, their analyses require every Markov chain to be close to some neighborhood of itsNote that SRRW is a nonlinear Markov chain that depends on the relative visit counts of each node in the graph. While its application in single-agent optimization has been studied in[39], expanding the theoretical examination of SRRW to multi-agent scenarios is beyond the scope of this paper. However, we can still numerically evaluate the performance of UD-SGD with multiple agents on general communication matrices using SRRW as a highly efficient sampling strategy in SectionTheorem D.1 is slightly modified in terms of condition C3, which is mentioned as a special case in Section 2.2[29]. For the sake of mathematical simplicity, we stick to condition C3 in the proof.However, we note that[29,22] considered the Lipschitz continuity of function F i θ i (x) defined in(68) as an assumption instead of a conclusion, where we give a detailed proof for this. We also obtain matrix Ui in an explicit form, which coincides with the definition of asymptotic covariance matrix and was not simplified in[29]. The discussion on the improvement of Ui is outlined in Section 3.2, which was not the focus of[29,22] and was not covered therein.In Appendix F.2, we will discuss the mathematical equivalence between our client sampling strategy and the commonly used one in the FL literature[50,30]. 


Decision: All reviewers agree that this is a good paper. It introduces theoretical advancements in a generalized distributed learning framework, offering a refined analysis that captures the behavior of individual agents. The new asymptotic analysis effectively explains linear speedup with the number of agents and demonstrates asymptotic network independence. The findings also show that upgrading even a small subset of agents can positively impact the entire system, a result validated by the experiments. Consequently, the paper is recommended for acceptance.

Reviews: review: paper_summary: paper_summary: This paper conducts an asymptotic analysis of Unified Distributed SGD (UD-SGD), which has a generalized communication patterns (modelled with a doubly stochastic communication matrix). The paper investigates several different sampling strategies, such as i.i.d. sampling, shuffling, and Markovian sampling
strength_weakness: strength_weakness: Empirical results are quite sparse. Not that this is a big problem for a theory paper, but maybe some additional studies would make it more interesting. I have some suggestions:
- Are all of the results obtained using the same communication matrix (referring to line 982 in appendix G1)? Since the theory is generic for all double stochastic W maybe the authors could present the average result over 5 random W?
- How was local data allocated to each client? The text said each client holds 500 data points from a bigger dataset -- are they iid? If so, it might be interesting to investigate non iid local data (using Dirichlet distribution, which is a common setting in FL)
- How does client dropout (also a common setting in FL) affect the theory and empirical results?
questions: questions: I have no further question.
limitations: limitations: No negative societal impact.
paper_summary: paper_summary: This paper studies the asymptotic convergence behavior of federated learning under the UD-SGD framework with Markovian data. The authors establish a new central limit theorem that considers the strategy of every agent, which goes beyond the existing bounds that only focus on the worst-performing agent. Their theory emphasizes the importance of every individual and also explains linear speedup and asymptotic network independence.
strength_weakness: strength_weakness: 1. As the author mentioned, Assumption 2.4 seems strong and finite-sample analysis is preferred,
2. As far as I could see in the experiments, the author only implemented DSGD but the UD-SGD framework has more algorithms.
3. The linear speedup scaling with the number of agents and asymptotic network independence are only mentioned in theory but do not have empirical justification.
4. The experiments for neural networks are somewhat toy.
questions: questions: 1. What is the main technical challenge when utilizing Poisson equation to prove Theorem 3.3? Specifically, what is new compared with the analysis in reference [23] and [30]?
2. Do the authors believe the current asymptotic analysis technique could be extended to (federated) reinforcement learning? I saw some finite-sample analyses for federated RL in your references like [72].
limitations: limitations: Yes.
paper_summary: paper_summary: This paper provides an asymptotic analysis of Unified Distributed SGD (UD-SGD) under heterogeneous agent dynamics and a large family of communication topology. It shows that under certain assumptions: i) regularity of the gradient, **ii) Ergodicity of Markovian sampling**, iii) decreasing step size and intervals not increasing so fast, iv) stability on model parameter v) contraction property of douly-stochastic communication matrix, UD-SGD algorithm guarantees that every agent parameter $\theta_n$ converges to the same minimum value $\theta^*$ with some variance matrix $V$ and the mean value of all agents $\overline{\theta}=\sum_{i=1}^n\theta_i$ also converges to $\theta_i$ with some variance matrix $V'$.
strength_weakness: strength_weakness: 1. In line 26, the authors claim that $\mathcal{L}$ represents the collection of local minima of objective function $f$. In line 246, the authors claim that $\theta^*\in \mathcal{L}$. Since when $f$ is non-convex, $\mathcal{L}$ do not have a single element and $\theta^*$ could have a lot of selections. How could we ensure that the convergence of consensus $\theta_n$ is unique? And how to choose $\theta^*$? Will that be the global minima?
2. The asymptotic analysis Theorem 3.3 is mainly based on Assumption 2.2, which assumes that the dynamic agent $\{X_{n}^i\}$ is an ergodic Markov chain. Could the author give some sampling strategies that satisfy this property and are now widely used in DSGD? What are the advantages of Markov sampling compared with the i.i.d sampling method?
3. In figure 2 of the experimental section, the authors only compare the performance of different Markov sampling methods. Could the author also compare with the iid sampling and sampling shuffle strategies?
questions: questions: Same as the weakness.
limitations: limitations: Same as the weakness.

Comments: title: Thank you for the detailed response
comment: I don't find many weaknesses to begin with, and I'm happy to maintain my positive score.
Thanks to the authors for their detailed response. They answer all my concerns and I would like to keep my score.
Thank you for the detailed response. I do not have further questions now.