Paper deep dive
Cordial Learning: Distributed Training with Correlated Data
Sarah Shitrit, Ilai Bistritz
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:We consider a distributed learning task with agents that have correlated data. Specifically, the label of an agent depends on the input of other agents for the same sample, and these inputs are also correlated. Correlated data is the reality when agents share the same environment. Existing decentralized methods, such as federated learning, ignore the structure of the problem and perform poorly on correlated data. On the other hand, centralized approaches are infeasible due to privacy and communication constraints. We introduce cordial (correlated and distributed) learning to address this gap by sharing only low-dimensional outputs between the agents while training local models to extract informative signals from peers. This distributed learning induces a game in which the loss function of each agent depends on the models of others. Assuming a linear model, we prove that cordial learning converges with probability one to a globally optimal solution, despite the nonconvex global objective. Experiments on structured multi-digit MNIST tasks demonstrate that cordial learning remains highly effective even in highly nonlinear settings.
Tags
Links
- Source: https://arxiv.org/abs/2610.03330v1
- Canonical: https://arxiv.org/abs/2610.03330v1
Trouble viewing inline? Open PDF directly →
Full Text
115,249 characters extracted from source content.
Expand or collapse full text
Cordial Learning: Distributed Training with Correlated Data Sarah Shitrit sarahshitrit@mail.tau.ac.il Affiliation: School of Electrical and Computer Engineering Affiliation: Tel Aviv University Affiliation: Tel Aviv, Israel Ilai Bistritz ilaibistritz@tauex.tau.ac.il Affiliation: School of Industrial and Intelligent Systems Engineering and School of Electrical and Computer Engineering Affiliation: Tel Aviv University Affiliation: Tel Aviv, Israel Abstract We consider a distributed learning task with agents that have correlated data. Specifically, the label of an agent depends on the input of other agents for the same sample, and these inputs are also correlated. Correlated data is the reality when agents share the same environment. Existing decentralized methods, such as federated learning, ignore the structure of the problem and perform poorly on correlated data. On the other hand, centralized approaches are infeasible due to privacy and communication constraints. We introduce cordial (correlated and distributed) learning to address this gap by sharing only low-dimensional outputs between the agents while training local models to extract informative signals from peers. This distributed learning induces a game in which the loss function of each agent depends on the models of others. Assuming a linear model, we prove that cordial learning converges with probability one to a globally optimal solution, despite the nonconvex global objective. Experiments on structured multi-digit MNIST tasks demonstrate that cordial learning remains highly effective even in highly nonlinear settings. †shortheadings: Distributed Training with Correlated Data / Shitrit and Bistritz†firstpage: 1 keywords Federated Learning, Distributed Training, Game Theory, Collaborative Learning, Stochastic Optimization 1 Introduction When deployed, ML systems rarely operate in a vacuum. Their data comes from a shared environment with many ML systems. Thus, the data of many ML systems can be highly correlated. Current ML architectures are trained in solitude, ignoring any interaction they may encounter after deployment. In distributed learning, and federated learning (FL) in particular, all edge devices (agents) face the same learning task. The data collected by an edge device is valuable to all others since it extends their effective training dataset. FL then avoids sending explicit data to the server, maintaining user privacy, while still allowing agents to learn from each other’s data indirectly. When ML agents share an environment, one agent’s data is helpful not only because more data is better, but mostly because it shows the agent the part of the bigger picture that it cannot see alone. Technically, it means that the label of one agent depends on the input of others. Even if the agents have different ML tasks, so they cannot learn directly from others’ data, the data of others contains information that explains their label better. The current distributed learning paradigm cannot benefit from the correlated data of ML agents that share an environment after deployment. In this paper, we introduce a new paradigm we term cordial learning (CL), since it is correlated and distributed. Similar to FL, CL is concerned with edge devices that each have a private dataset. However, the l-th datapoint of all devices is correlated (input and label), for every index l. CL is the practical reality in the following example applications: Multi-Digit MNIST: Each classifier sees a handwritten digit as input. However, the digits of all agents are dependent (e.g., credit card numbers, error correction codes, or the same handwriting). For example, adjacent digits must add up to 10, or all digits tend to equal a common hidden digit. Remote Sensors: The agents are sensors that collect local environmental features (e.g., temperature, vibration, gas). The task is to classify whether a global phenomenon, such as a fire or earthquake, is occurring. Since all sensors in the same region are exposed to the same event, their data (inputs and labels) are correlated. However, each sensor only observes part of the full picture. Sharing model outputs, rather than raw data, allows the agents to benefit from the collective information while meeting the communication constraints. Medical Wearables in Shared Environments: Each agent monitors a different patient using wearable sensors that record local physiological signals (e.g., heart rate, skin temperature). The task is to detect the onset of a disease. When patients are nearby, such as on a hospital ward or during an outbreak, their health outcomes are correlated, as infection spreads. In this case, leveraging cross-patient correlations can improve prediction accuracy. However, direct data sharing is prohibited due to privacy constraints. A distributed learning approach that shares only model outputs can enable collaboration without compromising sensitive health data. Formation Control: A group of moving agents often has to maintain a target formation. Examples are platoons of autonomous vehicles or drones. Each vehicle collects local sensor data (e.g., LiDAR, radar, velocity) and must classify whether to accelerate, maintain speed, or decelerate. Since vehicles influence each other’s decisions, especially those of their neighbors, their data is correlated. However, each vehicle observes only partial information. Sharing embedded outputs rather than raw sensory data allows collaborative real-time decision-making while preserving autonomy. Our goal is to design a distributed training scheme that achieves close to optimal performance in this cordial learning scenario, while maintaining low communication overhead. We measure performance by the sum of loss functions of all agents. A centralized architecture would collect data from all agents and train one large model to fit their combined input to their combined labels. Then, the central server is needed in real-time, so all devices can send their input data to the server and receive the result back. Such a centralized architecture requires significant communication between a central server and the agents to collect all their data, both in training and in real-time. Then the central server must train a huge ML model with a large dataset to obtain all agent models. Perhaps even worse, collecting user data in the central server can violate user privacy. It also creates a single point of failure. Communication, latency, complexity, and privacy are known issues with any centralized architecture. In the context of cordial learning, the centralized architecture is also data-inefficient. The classification tasks are inherently distributed in the sense that each agent observes a different part of the whole problem. A centralized agent will waste a lot of data “understanding” the distributed nature of the problem (e.g., that there are N digits in the multi-digit example). In other words, since the problem is distributed, a data-efficient architecture will also have to be distributed. In CL, the learning tasks of the agents are intertwined. If two agents each observe a digit, knowing that both digits are even or odd together, then each agent only needs the other agent to detect the parity of its digit. However, being able to classify a digit is very much related to classifying its parity. On the other hand, only sending the decision is not enough information. If an agent believes that the digit is six, but digits nine and three together are more likely, then this information is lost. A similar tradeoff exists in remote sensors that detect fire - even if one sensor is sure that currently there is no fire in its region, some indicators can help sensors in more sensitive areas to lower their detection threshold. We propose using an embedded output as the basic communication unit between the agents. The embedding is in the dimension of the output, a short vector compared to the input vectors. This embedded output is similar to the softmax output of a neural network (N), but not the same. The softmax output of our agents is computed after incorporating the information from other agents, while the embedded output is processed before. In CL, agents maintain coordination layers dedicated to extracting useful information from the inputs they receive from other agents for the same sample. Thus, each agent simultaneously trains its N along with these additional layers. This introduces a unique coupling (or interaction) between the agents. The weights of the N of agent n affect not only the loss of agent n, but also the loss of all other agents, since it changes the information they receive from agent n. As our analysis reveals, this interaction requires careful design. The sum of loss functions is nonconvex even in the linear regression case. Thus, many suboptimal stationary points may exist. For the linear theoretical model, we show that if the learning of the “coordination layers” is slower than the learning of the embedding layers, then CL is guaranteed to converge to the optimal agent models that minimize the nonconvex sum of loss functions. This happens even though each agent locally optimizes only its loss function, ignoring the effect of its model on the loss of others. This makes CL a scalable, distributed training algorithm that can incorporate many devices, with little communication overhead. 2 Previous Work Our work belongs to the growing literature on on-device learning. When the datasets of different agents are independent, each agent has a disjoint dataset for the same classification task. This corresponds to the classic on-device learning setting addressed by federated learning (Kairouz et al., 2021; Mitra et al., 2021; Koloskova et al., 2022), where no inter-agent correlation exists. Reducing the communication of federated learning and distributed SGD has been recognized as crucial (Wang and Joshi, 2021). Unlike CL, FL is not designed to deal with correlated data. Additionally, FL requires that all devices have the same architecture so that the weights can be averaged. In CL, different devices can have different architectures. In Bistritz et al. (2020), we designed a distributed distillation scheme in which agents share softmax outputs on common unlabeled data to learn from each other’s private data without revealing it. Since softmax is a short vector, the resulting communication overhead is dramatically lower than with FL, and devices with different ML architectures can participate in the training. Distributed distillation does not deal with correlated data. In CL, agents share embedded outputs to learn from other agents about the same sample they are currently learning from. Game theory and learning have recently become intertwined. The literature on learning in games studies the interaction between players who learn to maximize their reward. Then, convergence to a Nash equilibrium (NE) is guaranteed for various distributed algorithms in many scenarios (Xu et al., 2025; Ratliff et al., 2016; Bistritz et al., 2022; Chasnov et al., 2020; Mazumdar et al., 2020). We study the opposite - how games can help in learning. Game theory has already been proven insightful for learning in generative adversarial networks (Goodfellow et al., 2020, GANs;). Unlike the zero-sum games studied in GANs, our players are cooperative and each is an individual ML system (Bistritz and Bambos, 2020). Distributed learning with linear models that observe local measurements with correlated noise was studied in Hong et al. (2022). The risk functions include a term that penalizes for deviating from neighboring models. In our scenario, players do not necessarily converge to the same model. On the contrary, players may face non-identically distributed data and need to take heterogeneous actions (e.g., the medical wearables example above). Performative prediction studies a game scenario where the actions of one learner affect both the risk scores and the distributions of the data of the others (Narang et al., 2023; Lin and Zrnic, 2024; Perdomo et al., 2020). The goal is to compute an equilibrium that predicts the behavior of selfish learners. In contrast, our agents are cooperative. Furthermore, given the distribution, the datapoints of different players in Narang et al. (2023) are independent, while ours are correlated. Finally, Narang et al. (2023) assumes that players can observe each other’s actions, which is unacceptable in our case since it would mean that players can observe the model parameters of everyone. CL resembles vertical federated learning (VFL) (Liu et al., 2024). In VFL, the different agents have different nonoverlapping features of the same sample. This means that the inputs of the agents are correlated. However, in VFL, there is one learning task and one label, which only the “active agent” possesses. In contrast, in CL, each agent has its own local learning task and its own label, and the labels are correlated between agents. In this sense, VFL is a special case of CL, where all labels are identical and only one agent orchestrates the learning. Multitask learning (MTL) (Maurer, 2006) is also concerned with intertwined learning tasks. However, MTL trains a centralized architecture and not distributed agents. Even more important is that MTL is trained for tasks that occur asynchronously in real-time. In contrast, CL trains distributed agents that learn from synchronous correlated data points in real-time. While successive classification shares the idea of interdependent models, it differs from our setting in several fundamental ways. In pipelines such as Faster R-CNN (Ren et al., 2015), the region proposal network and the final classifier share a common feature extraction backbone, introducing architectural coupling. However, the overall computation proceeds sequentially: region proposals are generated and then classified. In contrast, CL involves a network of distributed agents, each owning a private task and dataset. These agents do not share a backbone or data. Instead, they coordinate through shared latent representations generated during training. Information flows bidirectionally, and learning is formulated as a nonconvex game where each agent’s model parameters affect others. (a) (b) (c) Figure 1: Different agent configurations for a digit classification task. 3 Problem Formulation We consider N distributed agents, each on a different device. Agent n has a labeled training dataset n=l(n),yl(n)l=1DD_n= \ x_l (n ),y_l (n ) \_l=1^D, where l(n)∈ℝd x_l (n ) ^d is an input and yl(n)∈ℝmy_l (n ) ^m is a label (either for classification or regression). We assume that [l(n),yl(n)]n=1N [ x_l (n ),y_l (n ) ]_n=1^N is the l-th sample of a random vector [n,Yn]n=1N [ X_n,Y_n ]_n=1^N. The crux of our model is that (1,Y1),…,(N,YN) ( X_1,Y_1 ),…, ( X_N,Y_N ) are statistically dependent. The datasets of all agents are then D independent samples of the vector [n,Yn]n=1N [ X_n,Y_n ]_n=1^N. This can be thought of as one large feature vector l=(l(1),…,l(N)) x_l= ( x_l (1 ),…, x_l (N ) ) that is distributed across agents, which is reminiscent of the idea of vertical federated learning (Liu et al., 2024). However, in cordial learning, each agent has their own learning task with their own label vector. Ideally, to train the N agents, we would wish to minimize the following empirical risk L(1,…,N)=∑n=1N∑l=1Dℓ(n,[l(m),yl(m)]m=1N)L ( ω_1,…, ω_N )= _n=1^N _l=1^D ( ω_n, [ x_l (m ),y_l (m ) ]_m=1^N ) (1) for a loss function ℓ (e.g., MSE), where n ω_n is the model parameters of agent n. A centralized architecture to solve (1) treats the N agents as one big model =(1,…,N) ω= ( ω_1,…, ω_N ) where the input is l=(l(1),…,l(N)) x_l= ( x_l (1 ),…, x_l (N ) ) and the output is yl=(yl(1),…,yl(N))y_l= (y_l (1 ),…,y_l (N ) ). However, each agent n can only observe its local data (l(n),yl(n)) ( x_l (n ),y_l (n ) ) and not (l(m),yl(m)) ( x_l (m ),y_l (m ) ) for m≠nm≠ n. For this reason, the traditional approach is just to train each model separately, such that the n-th agent minimizes its loss function: L(n)(n)=∑l=1Dℓ(n,l(n),yl(n)).L (n ) ( ω_n )= _l=1^D ( ω_n, x_l (n ),y_l (n ) ). (2) 3.1 Linear Mean-Squared Error Model In general ML scenarios, the loss functions can be highly nonconvex and intractable. To capture the essence of the correlated data scenario and analyze the theoretical performance guarantees of cordial learning, we adopt a linear model. In this model, the random label vector of agent i is a linear mixture of the input of all agents: Yi=AiXi+∑j=1,j≠iNBjiXj+νiY_i=A_iX_i+ _j=1,j≠ i^NB^i_jX_j+ _i (3) where Ai∈ℝm×dA_i ^m× d and Bji∈ℝm×dB^i_j ^m× d are the unknown distribution parameters, and νi∈ℝm _i ^m is a zero-mean noise with variance ‖νi‖2=σ2E \\| _i\|^2 \=σ^2, that is independent of all Xi\X_i\. We assume without loss of generality that Xi=0E \X_i \=0 and XiXiT=IE \X_iX_i^T \=I for all i. As usual, each agent can preprocess the inputs to satisfy these assumptions. We also assume that ‖Xi‖4<∞E \\|X_i\|^4 \<∞ for all i. The naive approach to this problem is to ignore the correlation between the agents and train each of them separately, using the following loss function for agent i: Li(Mi)=‖Yi−MiXi‖2L_i(M_i)=E \ \|Y_i-M_iX_i \|^2 \ (4) where the norm is the Frobenius norm ‖A‖=Tr(AAT) \|A \|= (A^T) (throughout the paper). Such a training scheme is suboptimal. On the other extreme, we could train all agents in a centralized manner, which is encoded in the following loss function for agent i: Li(M1i,…,MNi)=‖Yi−∑j=1NMjiXj‖2L_i(M_1^i,…,M_N^i)=E \ \|Y_i- _j=1^NM_j^iX_j \|^2 \ (5) so the model parameters are given by the matrices MjiM_j^i for all j. This is equivalent to a large learning model that faces a prediction task with an mNmN-dimensional output. By construction, the centralized approach can achieve the minimal training loss. However, centralized training will require each agent to share all its data with all others. Then the central server will still be needed in real-time to collect the inputs from all devices, process them, and send back the results to the agents. The communication overhead of such sharing is infeasible in large networks, especially in real-time. Moreover, sharing explicit data violates user privacy. Last but not least, when the task structure is distributed (e.g., remote fire sensors or medical sensors), a distributed architecture is more data-efficient than a centralized one that wastes data to understand the distributed structure of the problem. Our goal is to bypass the issues of a centralized training scheme while improving performance over the naive isolated scheme. To that end, we design a distributed learning algorithm where each agent i processes the output of each other agent j≠ij≠ i using dedicated coordination layers Λji _j^i. The loss function for agent i is then given by: Li(,)=‖Yi−(MiXi+∑j=1,j≠iNΛjiMjXj)‖2=‖(Ai−Mi)Xi+∑j=1,j≠iN(Bji−ΛjiMj)Xj+νi‖2.L_i ( M, )=E \ \|Y_i- (M_iX_i+ _j=1,j≠ i^N _j^iM_jX_j ) \|^2 \=\\ E \ \| (A_i-M_i )X_i+ _j=1,j≠ i^N (B^i_j- _j^iM_j )X_j+ _i \|^2 \. (6) Therefore, each agent’s loss depends on the models of other agents. Our distributed training approach introduces a continuous game with the agents as players. The action of each player i is its local model MiM_i and Λji _j^i for each j≠ij≠ i. At each turn t, each agent will choose a new action MiM_i and Λjij\ _j^i\_j. The cost function of a player is just its loss function given in (6). The game induced by (6) presents a conflict between the players. Each player seeks to adjust MiM_i to minimize their own loss function. However, by doing so they change the loss values of others. In general, pushing MiM_i in the direction of Ai∗A_i^* (defined in Assumption 1 below) conflicts with pushing MiM_i such that ΛijMi=Bij _i^jM_i=B_i^j. Motivated by practical settings, we make the following assumption, which means that the agent model matrices are not arbitrary. Instead, their right processing (as modeled by MiM_i) can “filter” the input for all agents. For example, while one agent may only be interested in whether the digit the other sees is even, the first layers that are required to detect the digit from the background are helpful for all agents. Assumption 1. Define Ai∗=Ai+∑j≠iBjiCov(Xj,Xi)A_i^*=A_i+ _j≠ iB^i_j (X_j,X_i ) where Cov(Xj,Xi) (X_j,X_i ) is the covariance matrix of XjX_j and XiX_i. We assume there exists Λji∗ _j^i* such that Λji∗A∗j=Bij _j^i*A^*_j=B^i_j for all i≠ji≠ j. Ai∗A_i^* is the optimal linear estimator of YiY_i from XiX_i. Therefore, the assumption above means that the information that agent i can potentially extract from its input, given an optimal model, is enough to compute the effect of agent i on all other agents. Assumption 1 is much weaker than assuming that every improvement of MiM_i towards Ai∗A_i^* also improves the quality of the information received by other agents, i.e., pushes closer to ΛijMi=Bij _i^jM_i=B_i^j. Such an assumption would have made the game induced by (6) a potential game, with the sum of loss functions as the potential function (Héliou et al., 2017). Nevertheless, the cordial learning game is not a potential game. Often, updates of MiM_i by agent i to decrease its loss function in (6) increase the loss values for other agents. Although cordial learning can work without Assumption 1, its optimality is based on this assumption. Naturally, the more reasonable this assumption is in a given scenario, the more valuable the information the agent sends to its peers. Our goal is to design a distributed training algorithm that minimizes the global objective: minMii,Λjii≠j∑i=1N‖Yi−(MiXi+∑j=1,j≠iNΛjiMjXj)‖2. \M_i \_i, \ _j^i \_i≠ j _i=1^NE \ \|Y_i- (M_iX_i+ _j=1,j≠ i^N _j^iM_jX_j ) \|^2 \. (7) In the CL algorithm, agent i uses the local variables Zii=XiZ_i^i=X_i and Zji=Xj−Cov(Xj,Xi)XiZ_j^i=X_j-Cov(X_j,X_i)X_i for j≠ij≠ i. Thus, ZiiZ_i^i is uncorrelated with ZjiZ_j^i for every j≠ij≠ i. We then work with the residualized versions of (6) and (7), obtained by replacing XjX_j with ZjiZ_j^i in the loss function of agent i. Assumption 1 implies that Mi=Ai∗M_i=A_i^* and Λji=Λji∗ _j^i= _j^i* for all i and j≠ij≠ i constitute a globally optimal solution of (7) when using the ZjiZ_j^i variables, achieving the minimal total MSE Nσ2Nσ^2. However, the sum of the agents’ loss functions is nonconvex (see Appendix B). Each agent’s loss is only convex in its own parameters when the other agents’ embedding matrices are held fixed. Thus, convergence to a global minimum is not guaranteed even for centralized gradient descent over Mi,Λjii,j\M_i, _j^i\_i,j. This constitutes a major technical challenge in designing a distributed algorithm that optimizes the residualized version of (7). 4 Cordial Learning Algorithm and Performance We propose a novel algorithm, termed cordial learning (correlated and distributed learning), to deal with the correlated data scenario. In cordial learning (CL), agents only share their embedded outputs and never their explicit datapoints. Then, each agent’s model includes coordination layers to process the inputs from other agents. The pseudocode for CL is given in Algorithm 1. To keep it tractable, our theoretical model is linear (Maurer, 2006), but CL applies more broadly to any scenario in which training data is correlated between agents. We support this claim with experiments on deep neural networks (see Section 6). CL is distributed since each agent calculates the gradients over its local parameters only. Instead of optimizing the global objective (7) directly, each agent optimizes its local objective given in (6). The local loss function depends on Mi,Λ1i,…,ΛNiM_i, _1^i,…, _N^i for agent i, where MjM_j for j≠ij≠ i, the models of other agents, are treated as constants. Centralized gradient descent that optimizes (7) would require each agent i to know the gradient of the loss function of agent j with respect to MiM_i. CL bypasses that, as each agent runs a two-time-scale SGD to minimize its loss function (6) with respect to its parameters alone. The stochastic gradients of an agent employ information sent by others. In general, SGD does not converge in games. When it does converge, it is to a stationary point that, under mild conditions, is a Nash equilibrium (NE) (Mazumdar et al., 2020), which is typically socially inefficient (Roughgarden, 2015; Sessa et al., 2019). Surprisingly, our main theorem shows that in the game induced by the designed loss functions in (6), a two-time-scale SGD converges to a global minimum of the sum of loss functions. Algorithm 1 Cordial Learning Initialization: Let YiY_i be the m×1m× 1 label vector of agent i. Let XiX_i be the d×1d× 1 input vector of agent i. Mi,ΛjijM_i,\ _j^i\_j are initialized with random values for each i. Let η(t),ε(t),ρ(t)η(t), (t),ρ(t) be sequences that are tuned according to the conditions of Theorem 2. For t=0,…,Titer−1t=0,…,T_iter-1 training iterations: 1. Sample Batches: Sample two independent batches Bt0,Bt1B_t^0,B_t^1 of B data indices uniformly at random from a global index set shared across agents. Stack the local input vectors corresponding to these indices in XiBt0X_i^B_t^0 and XiBt1X_i^B_t^1, matrices of dimensions d×Bd× B. 2. Communication: Send the embedded outputs MiXiBt0M_iX_i^B_t^0 and MiXiBt1M_iX_i^B_t^1 to all other agents. 3. Construct Local Variables: For each sample ID l in Bt0B_t^0 and each j≠ij≠ i, let R~j,0i,l=1B−1∑k≠l,k∈Bt0XjkXikT R_j,0^i,l= 1B-1 _k≠ l,k∈ B_t^0X_j^kX_i^kT and R~j,1i=1B∑k∈Bt1XjkXikT R_j,1^i= 1B _k∈ B_t^1X_j^kX_i^kT. Compute MjZ~j,0l=MjXjl−MjR~j,0i,lXil,MjZ~j,1l=MjXjl−MjR~j,1iXil.M_j Z_j,0^l=M_jX_j^l-M_j R_j,0^i,lX_i^l, M_j Z_j,1^l=M_jX_j^l-M_j R_j,1^iX_i^l. (8) 4. Update MiM_i: Mi(t+1)=Mi(t)−ε(t)Mi(t)M_i(t+1)=M_i(t)- (t) g_M_i(t) (9) where Mi(t)=−2YiBt0(XiBt0)T+2Mi(t)XiBt0(XiBt0)T+2∑k=1,k≠iN(Πt[Λki(t)]Mk(t)Z~k,0Bt0(t))(XiBt0)T g_M_i(t)=-2Y_i^B_t^0(X_i^B_t^0)^T+2M_i(t)X_i^B_t^0(X_i^B_t^0)^T+\\ 2 _k=1,k≠ i^N ( _t [ _k^i(t) ]M_k(t) Z_k,0^B_t^0(t) )(X_i^B_t^0)^T (10) and Πt[Λ] _t[ ] is the projection into the set Λ∈ℝm×m∣‖Λ‖F≤ρ(t) \ ^m× m \| \|_F≤ρ(t) \. 5. Update Λji _j^i for all j≠ij≠ i: Λji(t+1)=Λji(t)−η(t)Λji(t) _j^i(t+1)= _j^i(t)-η(t) g_ _j^i(t) (11) where Λji(t)=2(Mi(t)XiBt0−YiBt0)(Mj(t)Z~j,1Bt0(t))T+2∑k≠iΛki(t)Mk(t)Z~k,0Bt0(t)(Mj(t)Z~j,1Bt0(t))T. g_ _j^i(t)=2 (M_i(t)X_i^B_t^0-Y_i^B_t^0 ) (M_j(t) Z_j,1^B_t^0(t) )^T+\\ 2 _k≠ i _k^i(t)M_k(t) Z_k,0^B_t^0(t) (M_j(t) Z_j,1^B_t^0(t) )^T. (12) End Figure 2: Cordial Learning information flow To deal with the correlated input data, each agent i subtracts the component that is correlated with XiX_i from all of its received data, resulting in local data variables Z~1i,…,Z~Ni Z_1^i,…, Z_N^i. The goal of these local variables is to make the expected gradient with respect to MiM_i independent of Λjij≠i\ ^i_j\_j≠ i and of Mjj≠i\M_j\_j≠ i (Lemma 3), and also reduce the variance of the stochastic gradients (Lemma 4 and Lemma 5). Let Cov(Xj,Xi) (X_j,X_i ) be the covariance matrix of XjX_j and XiX_i. Define Zii=XiZ_i^i=X_i and Zji=Xj−Cov(Xj,Xi)XiZ_j^i=X_j- (X_j,X_i )X_i, which is the version of XjX_j after subtracting the information it holds about XiX_i. We then have MjZji=MjXj−MjCov(Xj,Xi)Xi=MjXj−Cov(MjXj,Xi)XiM_jZ_j^i=M_jX_j-M_j (X_j,X_i )X_i=M_jX_j- (M_jX_j,X_i )X_i (13) which implies that agent i, which receives MjXjM_jX_j, can compute MjZjiM_jZ_j^i if it knows Cov(MjXj,Xi) (M_jX_j,X_i ). To run CL, MjZjiM_jZ_j^i is needed in computing the stochastic gradients Mi g_M_i and Λji g_ _j^i. Since Cov(Xj,Xi)Cov(X_j,X_i) is unknown, CL constructs two independent estimates of this covariance. At each iteration, two independent batches Bt0B_t^0 and Bt1B_t^1 are sampled. The batch Bt0B_t^0 is used to compute the stochastic gradient, while Bt1B_t^1 is used to construct an independent covariance estimate. For each sample l in Bt0B_t^0, define R~j,0i,l=1B−1∑k≠l,k∈Bt0XjkXikT,R~j,1i=1B∑k∈Bt1XjkXikT. R_j,0^i,l= 1B-1 _k≠ l,k∈ B_t^0X_j^kX_i^kT, R_j,1^i= 1B _k∈ B_t^1X_j^kX_i^kT. The corresponding local variables are Z~j,0i,l=Xjl−R~j,0i,lXil,Z~j,1i,l=Xjl−R~j,1iXil, Z_j,0^i,l=X_j^l- R_j,0^i,lX_i^l, Z_j,1^i,l=X_j^l- R_j,1^iX_i^l, where for j=ij=i, we set Z~i,0i,l=Xil Z^i,l_i,0=X^l_i. The covariance estimators are independent of the l-th sample, and independent of each other. CL uses Z~j,0i,l Z_j,0^i,l in the stochastic gradient with respect to MiM_i. In the stochastic gradient with respect to Λji _j^i, it uses Z~k,0i,l Z_k,0^i,l in the prediction term and the independent Z~j,1i,l Z_j,1^i,l in the multiplicative term. As Lemma 3 shows, the resulting stochastic gradients are unbiased. CL utilizes common randomness between the agents, which generate the correlated samples that correspond to the same sample ID. This can be implemented using traditional epochs where agents go over all datapoints in a pre-agreed order. The communication between agents in CL can be implemented with or without a central server. The agents communicate both in training and online after deployment, to extract information from the correlation between their online test (unlabeled) samples using the coordination layers they trained offline. At test time, agent i computes MjZji=MjXj−Cov(MjXj,Xi)XiM_jZ^i_j=M_jX_j-Cov(M_jX_j,X_i)X_i using an estimate of Cov(MjXj,Xi)Cov(M_jX_j,X_i) obtained from training samples with the final trained models. Our main theorem states that CL converges with probability 1 to a global minimum of the sum of loss functions of all agents. Theorem 2 makes a standard SGD assumption on the step sizes of Mi(t)M_i(t) and Λji(t) _j^i(t). On top of that, we assume that the update of the Λji(t) _j^i (t ) is slower than that of the Mi(t)M_i (t ), formalized by the requirement that ∑t=1∞η2(t)ε(t)<∞ _t=1^∞ η^2 (t ) (t )<∞. The projection of Λki(t) _k^i (t ) into the ρ(t)ρ (t ) radius ball is only done when calculating the stochastic gradient Mi(t) g_M_i(t). The coordination layers Λki(t) _k^i (t ) only affect the variance of Mi(t) g_M_i(t) and not its bias. Hence, this projection ensures that the variance of Mi(t) g_M_i(t) does not explode. This is guaranteed by the assumption that ∑t=1∞ε2(t)ρ2(t)<∞ _t=1^∞ ^2 (t )ρ^2 (t )<∞. Theorem 2. Consider N agents that run cordial learning (Algorithm 1), where the 2B2B samples of [Xn,Yn]n=1N[X_n,Y_n]_n=1^N in Bt0B_t^0 and Bt1B_t^1 are independent within each iteration and IID over t. Assume that Aii\A_i\_i and Bjii,j\B^i_j\_i,j satisfy Assumption 1. Let ε(t),η(t) (t),η(t) be step size sequences of the embedding layers and coordination layers, respectively, and ρ(t)ρ (t ) be the projection radius. Assume that ∑t=1∞ε(t)=∞ _t=1^∞ (t )=∞, ∑t=1∞η(t)=∞ _t=1^∞η (t )=∞, ρ(t)→∞ρ (t )→∞ and ∑t=1∞ε2(t)ρ2(t)<∞ _t=1^∞ ^2 (t )ρ^2 (t )<∞. Further assume that ∑t=1∞η2(t)ε(t)<∞ _t=1^∞ η^2 (t ) (t )<∞. Then, for all i, the loss function Li((t),(t))=‖Yi−Mi(t)Zii−∑k=1,k≠iNΛki(t)Mk(t)Zki‖2L_i ( M(t), (t) )=E \ \|Y_i-M_i(t)Z_i^i- _k=1,k≠ i^N _k^i(t)M_k(t)Z_k^i \|^2 \ (14) converges to its global minimum with probability 1. 5 Convergence Analysis In this section, we detail our convergence analysis that leads to the performance guarantees of Theorem 2. Throughout our analysis, we use the following filtration, which summarizes the past iterates of the algorithm: ℱt=σ(Mi(t)i=1N,Λji(t)i≠j,…,Mi(0)i=1N,Λji(0)i≠j).F_t=σ ( \M_i (t ) \_i=1^N, \ _j^i (t ) \_i≠ j,…, \M_i (0 ) \_i=1^N, \ _j^i (0 ) \_i≠ j ). (15) Lemma 3 below computes the stochastic gradients used in Algorithm 1. Naturally, the stochastic gradient of the loss with respect to the embedding parameters MiM_i depends on the models of other agents, and on the coordination parameters i ^i. However, the expected gradients of the loss function of agent i with respect to MiM_i depend only on MiM_i. This unique fact is crucial to our analysis. In particular, we can project Λji(t) _j^i (t ) into a ball with a controlled radius ρ(t)ρ (t ) without creating any bias in Mi g_M_i. Consequently, we can first show the convergence of Mi(t)M_i (t ), which is then used to prove convergence of Λji(t) _j^i (t ). Lemma 3. Let Ck,jC_k,j be the covariance matrix between ZkiZ^i_k and ZjiZ^i_j. Recall that in Algorithm 1, agent i uses, for each sample in the batch, the stochastic gradients Mi=−2YiZ~i,0iT+2MiZ~i,0iZ~i,0iT+2(∑j=1,j≠iNΠt[Λji(t)]MjZ~j,0i)Z~i,0iT g_M_i=-2Y_i Z_i,0^iT+2M_i Z_i,0^i Z_i,0^iT+2 ( _j=1,j≠ i^N _t [ _j^i (t ) ]M_j Z_j,0^i ) Z_i,0^iT (16) and for an agent j≠ij≠ i Λji=−2Yi(MjZ~j,1i)T+2MiZ~i,0i(MjZ~j,1i)T+2(∑k=1,k≠iNΛkiMkZ~k,0i)(MjZ~j,1i)T. g_ _j^i=-2Y_i(M_j Z_j,1^i)^T+2M_i Z_i,0^i(M_j Z_j,1^i)^T+2 ( _k=1,k≠ i^N _k^iM_k Z_k,0^i )(M_j Z_j,1^i)^T. (17) Then, the expectations of the stochastic gradients are given by Mi=2(Mi−Ai−∑j=1,j≠iNBjiCov(Xj,Xi))E\ g_M_i\=2 (M_i-A_i- _j=1,j≠ i^NB^i_j (X_j,X_i ) ) (18) and Λji=2∑k=1,k≠iN(ΛkiMk−Bki)Ck,jMjT.E\ g_ _j^i\=2 _k=1,k≠ i^N ( _k^iM_k-B^i_k )C_k,jM_j^T. (19) In particular, the fact that the expected gradient with respect to MiM_i is 2(Mi−Ai∗)2(M_i-A_i^*) suggests that in the absence of noise, MiM_i will converge to Ai∗A_i^* with no consideration to other agents. This would not be the case if each agent took a gradient step to minimize the global objective in (7). The gradient of (7) with respect to MiM_i has N−1N-1 other terms that push MiM_i to solve ΛijMi=Bij _i^jM_i=B^j_i. The next lemma bounds the variances of our stochastic gradients. We see that Λji _j^i converges to minimal error only since MjM_j converges to Aj∗A_j^* for other agents, which indicates that learning is coupled between agents. Nevertheless, the convergence of Λji _j^i also helps learning MiM_i, since the variance of the stochastic gradient of MiM_i tends to decrease as the error of Λji _j^i decreases (see (22) in the proof of Lemma 4). This shows that the learning of MiM_i and Λji _j^i is coupled. The decreasing variance of Mi g_M_i as Λji _j^i converges is the reason for working with the local Z~1i,…,Z~Ni Z_1^i,…, Z_N^i variables, which subtract the influence of XiX_i from the received variables. The bounds of Lemma 4 are worst case and therefore do not show the dependence on the convergence error of Λji _j^i. They also ignore the dependence on the variance of R~ji R_j^i that decreases in the batch size. The proof of Theorem 2 only needs the bounds as stated in Lemma 4. Performance in practice, of course, benefits a lot from the decreasing variances of the stochastic gradients. Lemma 4. Let Λji∗ _j^i* be such that Bij=Λji∗Aj∗B^i_j= _j^i*A_j^*. Recall the stochastic gradients Mi(t) g_M_i (t ) and Λji(t) g_ _j^i (t ) defined in (16) and (17), for the parameters at iteration t. Then for all i ‖Mi(t)‖2|ℱt=O(ρ2(t)(1+∑j=1N‖Mj(t)−Aj∗‖2))E\ \| g_M_i (t ) \|^2|F_t\=O (ρ^2 (t ) (1+ _j=1^N \|M_j (t )-A_j^* \|^2 ) ) (20) and for all i,ji,j, understanding that Λi=0 _i^i=0, ‖Λji(t)‖2|ℱt=O(∑l=1N(1+‖Λli(t)‖2)∑k=1N(1+‖Mk(t)−Ak∗‖2)2).E \ \| g_ _j^i (t ) \|^2|F_t \=O ( _l=1^N (1+ \| _l^i (t ) \|^2 ) _k=1^N (1+ \|M_k (t )-A_k^* \|^2 )^2 ). (21) Proof. Let R~ji=1B−1∑k≠lXjkXikT R_j^i= 1B-1 _k≠ lX_j^kX_i^kT be the estimated covariance matrix based on the batch samples. Let A~i=Ai+∑j≠iBjiR~ji A_i=A_i+ _j≠ iB^i_j R_j^i, which approximates Ai∗A_i^*. For sufficiently large constants K0,…,K10>0K_0,…,K_10>0 14‖Mi(t)‖2|ℱt=‖(Ai−Mi(t))XiXiT+∑j≠i(BjiXj−Πt[Λji(t)]Mj(t)Z~j,0i)XiT+νiXiT‖2|ℱt=(a)‖(Ai−Mi(t))XiXiT+∑j≠i(Bji(Z~j,0i+R~jiXi)−Πt[Λji(t)]Mj(t)Z~j,0i)XiT+νiXiT‖2|ℱt=‖(Ai+∑j≠iBjiR~ji−Mi(t))XiXiT+∑j≠i(Bji−Πt[Λji(t)]Mj(t))Z~j,0iXiT+νiXiT‖2|ℱt≤3‖νiXiT‖2+3‖Mi(t)−A~i‖2‖XiXiT‖2|ℱt+3N∑j≠i‖Bji−Πt[Λji(t)]Mj(t)‖2‖Z~j,0iXiT‖2|ℱt≤(b)3‖νi‖2‖Xi‖2+3‖Mi(t)−A~i‖2‖XiXiT‖2+3N∑j≠i‖Bji−Πt[Λji(t)]Mj(t)‖2‖Z~j,0iXiT‖2≤(c)K0σ2+K6‖Mi(t)−A~i‖2+K1∑j≠i(‖Bji−Πt[Λji(t)]Aj∗‖2+‖Πt[Λji(t)](Aj∗−Mj(t))‖2)≤K0σ2+2K6‖Mi(t)−Ai∗‖2+2K6‖A~i−Ai∗‖2+K1∑j≠i(‖Bji−Πt[Λji(t)]Aj∗‖2+‖Πt[Λji(t)]‖2‖Aj∗−Mj(t)‖2) 14E\ \| g_M_i (t ) \|^2|F_t\=\\ E \ \| (A_i-M_i (t ) )X_iX_i^T+ _j≠ i (B_j^iX_j- _t [ _j^i (t ) ]M_j (t ) Z_j,0^i )X_i^T+ _iX_i^T \|^2\ |\ F_t \ (a )=\\ E \ \| (A_i-M_i (t ) )X_iX_i^T+ _j≠ i (B_j^i ( Z_j,0^i+ R_j^iX_i )- _t [ _j^i (t ) ]M_j (t ) Z_j,0^i )X_i^T+ _iX_i^T \|^2\ |\ F_t \=\\ E \ \| (A_i+ _j≠ iB_j^i R_j^i-M_i (t ) )X_iX_i^T+ _j≠ i (B_j^i- _t [ _j^i (t ) ]M_j (t ) ) Z_j,0^iX_i^T+ _iX_i^T \|^2\ |\ F_t \≤\\ 3E \ \| _iX_i^T \|^2 \+3E \ \|M_i (t )- A_i \|^2 \|X_iX_i^T \|^2|F_t \\\ +3NE \ _j≠ i \|B_j^i- _t [ _j^i (t ) ]M_j (t ) \|^2 \| Z_j,0^iX_i^T \|^2|F_t \ (b )≤\\ 3E \ \| _i \|^2 \E \ \|X_i \|^2 \+3E \ \|M_i (t )- A_i \|^2 \E \ \|X_iX_i^T \|^2 \\\ +3N _j≠ i \|B_j^i- _t [ _j^i (t ) ]M_j (t ) \|^2E \ \| Z_j,0^iX_i^T \|^2 \ (c )≤\\ K_0σ^2+K_6E \ \|M_i (t )- A_i \|^2 \+K_1 _j≠ i ( \|B_j^i- _t [ _j^i (t ) ]A_j^* \|^2+ \| _t [ _j^i (t ) ] (A_j^*-M_j (t ) ) \|^2 )≤\\ K_0σ^2+2K_6 \|M_i (t )-A_i^* \|^2+2K_6E \ \| A_i-A_i^* \|^2 \\\ +K_1 _j≠ i ( \|B_j^i- _t [ _j^i (t ) ]A_j^* \|^2+ \| _t [ _j^i (t ) ] \|^2 \|A_j^*-M_j (t ) \|^2 ) (22) where (a) uses that Z~j,0i=Xj−R~jiXi Z_j,0^i=X_j- R_j^iX_i, (b) uses that νi _i is independent of XiX_i and that R~ji R^i_j is independent of the l-th sample, and (c) uses that ‖XiXiT‖2E \ \|X_iX_i^T \|^2 \ and ‖Z~j,0iXiT‖2E \ \| Z_j,0^iX_i^T \|^2 \ are bounded. Similarly, for the stochastic gradient with respect to Λji _j^i, 14‖Λji(t)‖2|ℱt=‖((Mi(t)−Ai)Xi+∑k≠i(Λki(t)Mk(t)Z~k,0i−BkiXk)−νi)Z~j,1iTMjT(t)‖2|ℱt≤3‖νiZ~j,1iTMjT(t)‖2+3‖(Mi(t)−A~i)XiZ~j,1iTMjT(t)‖2|ℱt+3N∑k≠i‖(Λki(t)Mk(t)−Bki)Z~k,0iZ~j,1iTMjT(t)‖2|ℱt≤3‖νiZ~j,1iTMjT(t)‖2+3‖Mi(t)−A~i‖2‖Mj(t)‖2‖XiZ~j,1iT‖2+3N∑k≠i‖Λki(t)Mk(t)−Bki‖2‖Mj(t)‖2‖Z~k,0iZ~j,1iT‖2≤(a)K0σ2+K2(1+‖Mi(t)−Ai∗‖2+‖A~i−Ai∗‖2)(‖Mj(t)−Aj∗‖2+‖Aj∗‖2)+NK3∑k≠i(‖Λki(t)Ak∗−Bki‖2+‖Λki(t)‖2‖Mk(t)−Ak∗‖2)(‖Mj(t)−Aj∗‖2+‖Aj∗‖2)≤K5(‖Mi(t)−Ai∗‖2+‖Mi(t)−Ai∗‖2‖Mj(t)−Aj∗‖2+‖Mj(t)−Aj∗‖2+σ2+1)+K6∑k≠i‖Λki(t)Ak∗−Bki‖2+K7‖Mj(t)−Aj∗‖2∑k≠i‖Λki(t)Ak∗−Bki‖2+K8‖Mj(t)−Aj∗‖2∑k≠i‖Λki(t)‖2‖Mk(t)−Ak∗‖2+K9∑k≠i‖Λki(t)‖2‖Mk(t)−Ak∗‖2≤K10∑l=1N(1+‖Λli(t)‖2)∑k=1N(1+‖Mk(t)−Ak∗‖2)2 14E \ \| g_ _j^i (t ) \|^2|F_t \=\\ E \ \| ( (M_i (t )-A_i )X_i+ _k≠ i ( _k^i (t )M_k (t ) Z_k,0^i-B_k^iX_k )- _i ) Z_j,1^iTM_j^T (t ) \|^2|F_t \≤\\ 3E \ \| _i Z_j,1^iTM_j^T (t ) \|^2 \+3E \ \| (M_i (t )- A_i )X_i Z_j,1^iTM_j^T (t ) \|^2|F_t \\\ +3N _k≠ iE \ \| ( _k^i (t )M_k (t )-B_k^i ) Z_k,0^i Z_j,1^iTM_j^T (t ) \|^2|F_t \≤\\ 3E \ \| _i Z_j,1^iTM_j^T (t ) \|^2 \+3E \ \|M_i (t )- A_i \|^2 \ \|M_j (t ) \|^2E \ \|X_i Z_j,1^iT \|^2 \\\ +3N _k≠ i \| _k^i (t )M_k (t )-B_k^i \|^2 \|M_j (t ) \|^2E \ \| Z_k,0^i Z_j,1^iT \|^2 \ (a )≤\\ K_0σ^2+K_2 (1+ \|M_i (t )-A_i^* \|^2+E \ \| A_i-A_i^* \|^2 \ ) ( \|M_j (t )-A_j^* \|^2+ \|A_j^* \|^2 )\\ +NK_3 _k≠ i ( \| _k^i (t )A_k^*-B_k^i \|^2+ \| _k^i (t ) \|^2 \|M_k (t )-A_k^* \|^2 ) ( \|M_j (t )-A_j^* \|^2+ \|A_j^* \|^2 )≤\\ K_5 ( \|M_i (t )-A_i^* \|^2+ \|M_i (t )-A_i^* \|^2 \|M_j (t )-A_j^* \|^2+ \|M_j (t )-A_j^* \|^2+σ^2+1 )\\ +K_6 _k≠ i \| _k^i (t )A_k^*-B_k^i \|^2+K_7 \|M_j (t )-A_j^* \|^2 _k≠ i \| _k^i (t )A_k^*-B_k^i \|^2\\ +K_8 \|M_j (t )-A_j^* \|^2 _k≠ i \| _k^i (t ) \|^2 \|M_k (t )-A_k^* \|^2+K_9 _k≠ i \| _k^i (t ) \|^2 \|M_k (t )-A_k^* \|^2≤\\ K_10 _l=1^N (1+ \| _l^i (t ) \|^2 ) _k=1^N (1+ \|M_k (t )-A_k^* \|^2 )^2 (23) where in (a) we used that ‖νiZ~j,1iT‖2,‖Z~k,0iZ~j,1iT‖2,‖XiZ~j,1iT‖2E \ \| _i Z_j,1^iT \|^2 \,E \ \| Z_k,0^i Z_j,1^iT \|^2 \,E \ \|X_i Z_j,1^iT \|^2 \ are bounded and 3‖νiZ~j,1iTMjT(t)‖2≤6‖νiZ~j,1iTAj∗T‖2+6‖νiZ~j,1iT‖2‖MjT(t)−Aj∗T‖2.3E \ \| _i Z_j,1^iTM_j^T (t ) \|^2 \≤ 6E \ \| _i Z_j,1^iTA_j^*T \|^2 \+6E \ \| _i Z_j,1^iT \|^2 \|M_j^T (t )-A_j^*T \|^2 \. (24) ∎ The following lemma is not necessary to prove Theorem 2. However, it sheds more light on how CL works. While Theorem 2 shows that the loss of each agent goes to its global minimum, Lemma 5 shows that each effective coordination error Λji(t)Aj∗−Bji _j^i(t)A_j^*-B_j^i converges to zero with probability one. Since Λji(t) ^i_j(t) is bounded and ρ(t)→∞ρ(t)→∞, we have Πt[Λji(t)]=Λji(t) _t [ ^i_j(t) ]= ^i_j(t) for all large t. Together with Mj(t)→Aj∗M_j(t)→ A^*_j (Theorem 2) and the bound (22) in the proof of Lemma 4, it implies that the variance of Mi(t) g_M_i(t) has a decreasing upper bound that reaches an error floor, determined by the label noise and covariance estimation. The assumptions of Lemma 5 essentially rule out redundant dimensions in the learning problem. If a nonzero linear combination of X1,…,XNX_1,…,X_N is constant, then some feature observed by one agent is perfectly determined by features observed by the others and can therefore be removed without losing information. Lemma 5. Assume that no nonzero linear combination of X1,…,XNX_1,…,X_N is a constant with probability 1. Then in CL, for each i,ji,j (i≠ji≠ j), ‖Λji(t)Aj∗−Bji‖→0 \| _j^i (t )A_j^*-B_j^i \|→ 0 with probability 1 as t→∞t→∞. Proof. Let fi(Λ1i,…,ΛNi)≜‖∑k=1,k≠iN(ΛkiAk∗−Bki)Zki‖2f_i ( _1^i,…, _N^i ) \ \| _k=1,k≠ i^N ( _k^iA_k^*-B_k^i )Z_k^i \|^2 \. Without loss of generality, let i=1i=1, and define Δ(t)=(Λ21(t)A2∗−B21⋯ΛN1(t)AN∗−BN1). (t)= pmatrix _2^1(t)A_2^*-B_2^1&·s& _N^1(t)A_N^*-B_N^1 pmatrix. Let C1C_1 be the block covariance matrix of (Z21,…,ZN1)(Z_2^1,…,Z_N^1), whose (k,j)(k,j) block is Ck,jC_k,j. Then f1(Λ11(t),…,ΛN1(t))=‖∑k=2N(Λk1(t)Ak∗−Bk1)Zk1‖2=Tr(Δ(t)C1ΔT(t)).f_1\! ( _1^1(t),…, _N^1(t) )=E \| _k=2^N ( _k^1(t)A_k^*-B_k^1 )Z_k^1 \|^2\\ =Tr ( (t)C_1 ^T(t) ). (25) Let Rj,kR_j,k be the covariance matrix of XjX_j and XkX_k. Define R1=(R1,1⏟D(R1,2⋯R1,N)⏟B(R2,1RN,1)⏟BT(R2,2⋯R2,N⋱RN,2⋯RN,N)⏟U).R_1= ( array[]c|c R_1,1_D& pmatrixR_1,2&·s&R_1,N pmatrix_B\\ pmatrixR_2,1\\ \\ R_N,1 pmatrix_B^T& pmatrixR_2,2&·s&R_2,N\\ & & \\ R_N,2&·s&R_N,N pmatrix_U array ). (26) The Schur complement of D in R1R_1 is R1/D=U−BTD−1BR_1/D=U-B^TD^-1B, which equals C1C_1 since D=ID=I. Let X be the vector that stacks X1,…,XNX_1,…,X_N. If R1R_1 is not full-rank, then there exists u≠u≠ 0 such that R1u=0R_1u=0. However, the variance of uTXu^TX is uTR1u=0u^TR_1u=0, so uTXu^TX is a constant with probability 1, contradicting the assumption. Therefore, R1R_1 is positive definite, implying that its Schur complement C1C_1 is also positive definite. Therefore, f1(Λ11(t),…,ΛN1(t)) f_1\! ( _1^1(t),…, _N^1(t) ) =Tr(Δ(t)C1ΔT(t)) =Tr ( (t)C_1 ^T(t) ) ≥λmin(C1)‖Δ(t)‖F2=λmin(C1)∑j=2N‖Λj1(t)Aj∗−Bj1‖F2. ≥ _ (C_1)\| (t)\|_F^2= _ (C_1) _j=2^N \| _j^1(t)A_j^*-B_j^1 \|_F^2. The proof of Theorem 2 shows that f1(Λ11(t),…,ΛN1(t))→0f_1\! ( _1^1(t),…, _N^1(t) )→ 0 with probability one. Since λmin(C1)>0 _ (C_1)>0, it follows that ‖Λj1(t)Aj∗−Bj1‖F⟶0 \| _j^1(t)A_j^*-B_j^1 \|_F 0 for every j≠1j≠ 1. The same argument applies to every agent i. ∎ 6 Experiments In this section, we test CL numerically. The first set of simulations uses our linear theoretical model, assumed in Theorem 2. The second set goes beyond the linear model and runs CL with deep neural networks that classify the digits of a composite MNIST image. 6.1 Linear Model Simulations Input Data Generation. We simulated N∈8,16,32N∈\8,16,32\ agents. Agent i has an input Xi∈ℝdX_i ^d and a label Yi∈ℝmY_i ^m, with d=2d=2 and m=1m=1. As assumed in Section 3.1, the inputs have zero mean and are whitened, so that XiXiT=IE\X_iX_i^T\=I. The inputs of different agents are correlated coordinate by coordinate. For each coordinate k=1,…,dk=1,…,d, the vector (X1,k,…,XN,k)(X_1,k,…,X_N,k) has a correlation matrix R∈ℝN×NR ^N× N, and different coordinates are independent, so that Cov(Xj,Xi)=RjiICov(X_j,X_i)=R_jiI. In every run we draw a new matrix R with unit diagonal and off-diagonal entries (Uji+Uij)/2(U_ji+U_ij)/2, where the entries of U are drawn uniformly from [−5/N,5/N][-5/N,5/N]. Such a matrix is not always positive semidefinite. In that case, we replace it by a nearby correlation matrix: we raise its eigenvalues to at least 10−1210^-12 and rescale it to a unit diagonal. This was needed in about 60% of the runs with N=8N=8, in less than 10% of the runs with N=16N=16, and in none with N=32N=32. We use two input distributions. In the first, the inputs are Gaussian with correlation matrix R. In the second, the inputs have Laplace marginals: we draw Gaussian inputs as before and map every coordinate g to F−1(Φ(g))F^-1( (g)), where Φ is the standard normal distribution function and F is that of the Laplace distribution with unit variance. The coordinates of each agent stay independent with unit variance, but their correlations across agents change. For Laplacian inputs, Cov(Xj,Xi)=R~jiICov(X_j,X_i)= R_jiI, where R~ji R_ji is the correlation of the transformed coordinates, which we compute from RjiR_ji by numerical integration. In the rest of this section, Cov(Xj,Xi)Cov(X_j,X_i) refers to this matrix for Laplacian inputs. Output Generation. The entries of Ai∈ℝm×dA_i ^m× d and Λi∗j∈ℝm×m ^i*_j ^m× m are drawn independently from (0,0.52)N(0,0.5^2) for Gaussian inputs and from a Laplace distribution with scale 0.40.4 for Laplacian inputs. To ensure that the fixed-point iteration below is a contraction, the coordination matrices of each agent i are scaled by si=min1,κγi,γi=∑j≠i∥Λji∗∥2∥Cov(Xj,Xi)∥2,s_i= \1, κ _i \, _i= _j≠ i \| ^i*_j \|_2 \|Cov(X_j,X_i) \|_2, with κ=0.8κ=0.8, where ∥⋅∥2\|·\|_2 is the spectral norm. We then compute Ai∗A^*_i as the fixed point of the Picard iteration Ai∗←Ai+∑j≠iΛji∗Aj∗Cov(Xj,Xi),A^*_i← A_i+ _j≠ i ^i*_jA^*_jCov(X_j,X_i), set Bij=Λi∗jA∗jB^i_j= ^i*_jA^*_j, and generate the labels as Yi=AiXi+∑j≠iBjiXjY_i=A_iX_i+ _j≠ iB^i_jX_j. By construction, Assumption 1 holds and Ai∗A^*_i is the optimal linear estimator of YiY_i from XiX_i. The labels are noiseless (σ2=0σ^2=0), so the global minimum of the loss (14) of every agent is 00. Training Procedure. The agents run Algorithm 1. At each iteration, they sample two independent batches Bt0B_t^0 and Bt1B_t^1 of B=10B=10 fresh samples, so the samples are independent across iterations, as Theorem 2 assumes. The stochastic gradients (10) and (12) are averaged over the batch. The initial values Mi(0)M_i(0) and Λji(0) ^i_j(0) are random, with entries drawn from the same distribution as the entries of AiA_i and of Λi∗j ^i*_j before scaling. The step sizes and the projection radius are ε(t)=Lm(t+t0)α,η(t)=Lλ(t+t0)β,ρ(t)=(t+1)0.05, (t)= L_m(t+t_0)^α, η(t)= L_λ(t+t_0)^β, ρ(t)=(t+1)^0.05, with α=0.6α=0.6, β=1β=1 and t0=100t_0=100, and with LmL_m and LλL_λ such that ε(0)=0.1 (0)=0.1 and η(0)=0.05η(0)=0.05. These sequences satisfy all the conditions of Theorem 2: ∑tε(t)=∑tη(t)=∞ _t (t)= _tη(t)=∞, ρ(t)→∞ρ(t)→∞, ∑tε2(t)ρ2(t)<∞ _t ^2(t)ρ^2(t)<∞ and ∑tη2(t)/ε(t)<∞ _tη^2(t)/ (t)<∞. Isolated agents minimize (4) by SGD with the same step size ε(t) (t), using the batch Bt0B_t^0. Each agent trains for 20,000 iterations. Evaluation. For cordial agents, we report the loss (14) of each agent, which is the objective of Theorem 2, averaged over the N agents. For isolated agents, we report the loss (4). The smallest possible isolated loss is ‖Yi−Ai∗Xi‖2E\\|Y_i-A^*_iX_i\|^2\, attained by Mi=Ai∗M_i=A^*_i. We compute both losses in closed form from the known input covariance, so the curves carry no evaluation noise. The residualization in (14) uses the true covariance Cov(Xj,Xi)Cov(X_j,X_i), while the agents themselves only use the estimates of Algorithm 1. Each configuration was repeated for 100 independent runs, each with new parameters, a new initialization and new samples. Results. Figure 3 shows the loss of cordial and isolated agents for all six configurations. Isolated agents converge to the best isolated predictor within about 100 iterations. Their median loss is at most 1.6%1.6\% above the best isolated loss from then on. Cordial agents start with a higher loss, since their coordination matrices are initialized at random, but it drops quickly. After 20 iterations, it is lower than that of the isolated agents in 94–100% of the runs of each configuration, and from iteration 140 on, it is lower in all 600 runs. The loss of cordial agents then keeps decreasing towards the global minimum. After 20,000 iterations, the median ratio between the loss of cordial and isolated agents ranges from 0.36%0.36\% (N=8N=8, Gaussian inputs) to 1.4%1.4\% (N=32N=32, Laplacian inputs), and it is below 2.4%2.4\% in every run. On the log-log scale, the loss of cordial agents decreases almost linearly, with a slope between −0.58-0.58 and −0.49-0.49 from iteration 2,000 to 20,000 in all configurations. The loss therefore decreases like a power of t. The loss at the last iteration is larger with more agents and with Laplacian inputs, while the rate of decrease is about the same in all configurations. Figure 3: Linear model: loss of cordial agents (14) and of isolated agents (4) for N=8,16,32N=8,16,32 agents with Gaussian (top) and Laplacian (bottom) inputs. The lines show the median over 100 runs of the loss averaged over the agents, and the shaded bands show the 10th to 90th percentiles across runs. The dashed line is the median loss of the best isolated predictor Mi=Ai∗M_i=A^*_i. The global minimum of the loss of cordial agents is 00. Both axes are logarithmic. 6.2 MNIST-Based Multi-Agent Simulations Input Data Generation. Each data sample is a horizontal concatenation of N MNIST digits, which forms a 28×28N28× 28N image. Agent i sees the i-th 28×2828× 28 digit and has to classify it. The digits of a sample are correlated, and the rule that generates them defines the dataset. Section 6.2.1 uses a deterministic rule (Pairwise-10), and Section 6.2.2 a probabilistic one (a common cause). We vary the amount of data by the number P of MNIST images per digit. The training composites are built from a pool of P MNIST training images of each digit, and the validation composites from a second pool of P training images of each digit, disjoint from the first. The test composites are built from a pool of P MNIST test images of each digit. Therefore, no image appears in two of the three sets. A composite takes a random image from the pool of each of its digits. Pixel values are scaled to [−1,1][-1,1]. We call one such random construction of the training, validation and test sets a draw of the images. Neural Network Architectures. Every agent runs the LeNet-5-style network of Figure 4. Its embedding layers are conv1 (6 filters of size 5×55× 5, padded to keep the 28×2828× 28 size), pool1 (2×22× 2 max pooling), conv2 (16 filters of size 5×55× 5), pool2, conv3 (120 filters of size 5×55× 5), fc1 (84 units) and fc2 (10 units). All layers use ReLU activations, and dropout follows every convolution, every pooling layer and fc1. The output hi∈ℝ10h_i ^10 of fc2, after its ReLU, is the embedded output of agent i, the analogue of MiXiM_iX_i. We compare three architectures: • Cordial Learning: Each agent sends hih_i to all others, both in training and at test time. During training, the communicated embedded outputs are computed with dropout, as in the forward pass. In its layer fc3, agent i computes zi=Wihi+bi+∑j≠iΛjihj+ciz_i=W_ih_i+b_i+ _j≠ i ^i_jh_j+c_i, where Wi,Λji∈ℝ10×10W_i, ^i_j ^10× 10 are weight matrices and bi,ci∈ℝ10b_i,c_i ^10 are biases, and its prediction is the softmax of fc4(ReLU(zi))(ReLU(z_i)), where fc4 has 10 units. The matrices Λji ^i_j and the bias cic_i form the coordination layer, which adds the message term ∑j≠iΛjihj+ci _j≠ i ^i_jh_j+c_i to ziz_i. It is initialized to zero. The received embedded outputs hjh_j are treated as constants, so the parameters of agent i are updated only by the gradient of its own loss, as in Algorithm 1. • Isolated Agents: The same network without the message term. The agents do not communicate. • Centralized Architecture: A single standard CNN that receives the entire 28×28N28× 28N image and predicts all the digits jointly. It has three convolutional layers with 5×55× 5 filters (padded to keep the size), each followed by ReLU and 2×22× 2 max pooling, a fully connected hidden layer, and an output layer with 10N10N units, with a softmax over each group of 10. Dropout follows every pooling layer and the hidden layer. The numbers of filters and hidden units are tuned from (6, 16, 120) filters with 128 hidden units, (16, 32, 64) with 128 or 256, and (32, 64, 128) with 256. It is trained on the sum of the losses of all agents. Training and Evaluation. All models are trained with Adam on batches of 65 composites, reshuffled every epoch. The loss is the sum over the agents of the cross-entropy of each agent’s digit. In epoch e, the learning rate of all layers except the coordination layers is εe=ε0/(1+e/t0)0.6 _e= _0/(1+e/t_0)^0.6. The coordination layers are not trained in the first K epochs, and after that their learning rate is ηe=η0/(1+(e−K)/t0) _e= _0/(1+(e-K)/t_0). Therefore, ηe/εe→0 _e/ _e→ 0 and ∑eηe2/εe<∞ _e _e^2/ _e<∞, as in Theorem 2. With Adam, however, these are nominal learning rates and not the step sizes of Theorem 2. The delay K lets the embeddings become informative before the agents start to rely on them. Isolated agents and the centralized network use the schedule εe _e for all their layers. The hyperparameters of every method were tuned with 4 agents on the validation set of the first draw of the images (see below), separately for every dataset and data size, with the same number of runs for all methods. The initial learning rate ε0 _0 is drawn log-uniformly from [5⋅10−4,5⋅10−2][5· 10^-4,5· 10^-2] ([10−4,10−2][10^-4,10^-2] for the centralized network), t0t_0 from 1,3,10,30,100,300\1,3,10,30,100,300\ (without the largest values for larger P), and the dropout rate from 0.1,0.25,0.5\0.1,0.25,0.5\. For cordial agents, η0 _0 is also drawn log-uniformly from [10−4,10−1][10^-4,10^-1], and K from a grid that scales with the number of epochs (e.g., 0,25,50,100,200,400\0,25,50,100,200,400\ for P=25P=25). For the centralized network, we also tune the numbers of filters and hidden units. For each method, we trained 40 random configurations, then 40 perturbations of the five best, and then retrained the three best with three more random seeds. For Pairwise-10 and for the common-cause dataset with p=0.6p=0.6, with P=25P=25 and P=500P=500, the candidates also include the configurations of an earlier random search, the same for all methods (12 and 10 per method, respectively, with epoch limits of 1,000 and 120). The configuration with the lowest validation loss, averaged over its four seeds, was selected. Tuning runs stop after 200, 80, and 15 epochs without an improvement of the validation loss for P=25P=25, 100, and 500, respectively, or at an epoch limit of 1,500, 600, and 180 epochs. With the selected hyperparameters, each method was trained with 10 random seeds for a fixed number of epochs: either the maximal number of epochs for tuning, or 1.5 times this number if one of the selected configurations reached its best validation loss in the last 10% of the tuning epochs. The number of agents N∈2,4,6N∈\2,4,6\ changes, but the hyperparameters tuned with 4 agents are kept. We keep the network with the lowest validation loss during training, which we evaluate every 5, 2, and 1 epochs for P=25,100P=25,100, and 500 (for cordial agents, from epoch K on), and evaluate it on the test set without dropout. We report the test loss per agent (the sum over agents of the mean test cross-entropy, divided by N) and the test accuracy averaged over the agents. For cordial agents, we also report the test accuracy when the message term is removed, which shows how much the agents rely on the received embedded outputs. The selected hyperparameters are listed in Table 1. With few images, the results depend on which images are drawn into the three pools. We therefore repeated the final runs of all methods on two further independent draws of the training, validation and test pools, with the same hyperparameters, numbers of epochs and random seeds, and report the results of all three draws together (30 runs per method). To compare two methods, we compare every run of one method with every run of the other on the same draw (100 pairs per draw) and report the number of pairs in which the first method has the lower test loss. The received embedded outputs are not residualized, so the coordination layers use the raw values hjh_j. Without residualization, the second batch Bt1B_t^1 is not needed. We also do not project the coordination layers. Cordial learning Isolated agents Centralized Dataset P ε0 _0 η0 _0 K t0t_0 Drop. ε0 _0 t0t_0 Drop. ε0 _0 t0t_0 Width Drop. Pairwise-10 25 0.0023 0.0052 50 300 0.5 0.0012 300 0.5 0.006 300 16-32-64; 128 0.5 100 0.0025 0.0071 100 100 0.5 0.0015 100 0.25 0.0041 100 16-32-64; 256 0.5 500 0.0011 0.018 10 10 0.25 0.001 30 0.25 0.0017 30 32-64-128; 256 0.5 Common, p=0.8p=0.8 25 0.0041 0.00017 100 100 0.5 0.0018 300 0.5 0.002 30 32-64-128; 256 0.5 100 0.0065 0.00019 0 30 0.25 0.0017 100 0.25 0.0025 100 32-64-128; 256 0.5 500 0.001 0.0001 0 30 0.25 0.0052 10 0.25 0.0015 30 32-64-128; 256 0.5 Common, p=0.6p=0.6 25 0.004 0.0013 400 30 0.5 0.0017 300 0.5 0.00075 10 32-64-128; 256 0.5 100 0.0044 0.00031 10 100 0.25 0.0014 30 0.25 0.0032 100 16-32-64; 128 0.5 500 0.0014 0.00028 0 30 0.25 0.0019 10 0.25 0.0012 30 32-64-128; 256 0.5 Table 1: Selected hyperparameters for Pairwise-10 and the common-cause datasets with p=0.8p=0.8 and p=0.6p=0.6 (tuned with 4 agents): initial learning rates ε0 _0 and η0 _0, decay constant t0t_0, start K of the coordination layers and dropout rate. The width of the centralized network gives the numbers of filters of its three convolutional layers and the number of hidden units. Figure 4: N Architecture 6.2.1 The Pairwise-10 Dataset In the Pairwise-10 dataset, the digits of a sample satisfy two rules. The digits at positions (1,2),(3,4),…(1,2),(3,4),… sum to 10, and so do the digits at positions (2,4),(4,6),…(2,4),(4,6),…. The valid tuples of digits are therefore (a,10−a)(a,10-a), (a,10−a,10−a,a)(a,10-a,10-a,a) and (a,10−a,10−a,a,a,10−a)(a,10-a,10-a,a,a,10-a) for N=2,4,6N=2,4,6, with a∈1,…,9a∈\1,…,9\. In particular, the digit of any single agent determines the digits of all other agents. Each of the three sets has P composites of every valid tuple, that is, 9P9P composites. We report P=25P=25, 100, and 500, from 225 to 4,500 composites per set. The final runs have 2,250, 600, and 180 epochs, respectively. Results. Figure 5 shows the test loss during training for P=25P=25, 100 and 500, and Table 2 lists the test loss and accuracy of the selected networks. CL is much better than isolated training in every setting. For P=25P=25, the median test loss per agent of isolated agents is 0.450.45, 0.300.30 and 0.300.30 for 2, 4 and 6 agents. Cordial learning reduces it to 0.230.23, 0.0140.014 and 0.00280.0028, and raises the median accuracy from 0.880.88, 0.910.91 and 0.910.91 to 0.940.94, 0.9960.996 and 1.0001.000. For P=500P=500, isolated agents reach an accuracy of 0.980.98 and a test loss of 0.060.06–0.070.07, and CL reduces the test loss to 0.0080.008, 0.00090.0009 and less than 10−410^-4. In every setting, every run of CL has a lower test loss than every run of isolated agents on the same draw, except for 2 agents with P=25P=25, where this holds for 297 of the 300 pairs of runs. The gain grows with the number of agents, which is expected, since every additional agent brings another view of the same information. The draw of the images matters for P=25P=25: the median test loss of isolated agents with 2 agents is 0.310.31, 0.470.47 and 0.490.49 on the three draws, and that of CL is 0.0930.093, 0.300.30 and 0.260.26. The advantage of CL holds on every draw. The centralized network shows what can be achieved when observing all digits. With P=500P=500, it reaches a test loss of 0.0020.002 for 2 agents and less than 0.00010.0001 for 4 and 6 agents. CL reaches almost the same accuracy (at least 0.9980.998), although the agents only exchange their 10-dimensional embedded outputs and never their inputs or gradients. The centralized network has the lower test loss in most pairs of runs. With less data, the centralized network often failed to train. With the selected configuration, 47 of its 180 runs with P=25P=25 and P=100P=100 did not learn at all: the network predicts the uniform distribution over the nine possible digits for every agent, with a test loss of ln9≈2.2 9≈ 2.2. This happened in 7, 12 and 15 of the 30 runs with P=25P=25 and 2, 4 and 6 agents, and in 0, 6 and 7 of the 30 runs with P=100P=100. This occurred even though none of the four tuning runs of the selected configuration failed. None of the 180 runs of CL with P≤100P≤ 100 failed, and their lowest accuracy is 0.900.90. For P=25P=25 and 2, 4 and 6 agents, the median test loss of CL is 0.230.23, 0.0140.014 and 0.00280.0028, and that of the centralized network is 0.300.30, 0.0270.027 and 1.11.1. Without the runs that did not learn, the median test loss of the centralized network is 0.230.23, 0.00020.0002 and less than 10−410^-4. The gain of CL comes from the received embedded outputs. When the message term is removed, the accuracy of cordial agents drops to between 0.650.65 and 0.950.95, below that of isolated agents in every setting. Cordial agents learn to combine their own image with the embedded outputs of the others, instead of classifying their own digit alone. This fits the structure of Pairwise-10, where a single received embedded output can carry as much information about the digit of agent i as its own image. Figure 5: Pairwise-10: test loss per agent during training of isolated agents, cordial agents and the centralized network. The lines show the median over 30 runs (10 random seeds on each of three draws of the images) and the shaded bands the interquartile range. The sharp drop in the loss of cordial agents is at epoch K, when the coordination layers start training. The vertical axis is logarithmic, and losses below 10−510^-5 (including losses that are exactly zero in single precision) are drawn at 10−510^-5. Test accuracy P N Method Runs messages removed Test loss per agent 25 2 Isolated 30 0.883 – 0.454 [0.397, 0.495] Cordial learning 30 0.944 0.760 0.232 [0.117, 0.302] Centralized 30 (7) 0.963 – 0.299 [0.146, 0.826] 4 Isolated 30 0.913 – 0.295 [0.274, 0.421] Cordial learning 30 0.996 0.684 0.014 [0.0087, 0.036] Centralized 30 (12) 0.996 – 0.027 [0.00013, 2.197] 6 Isolated 30 0.911 – 0.296 [0.191, 0.364] Cordial learning 30 1.000 0.651 0.0028 [0.00088, 0.0044] Centralized 30 (15) 0.547 – 1.128 [<10−4<10^-4, 2.197] 100 2 Isolated 30 0.960 – 0.182 [0.161, 0.205] Cordial learning 30 0.993 0.928 0.022 [0.016, 0.039] Centralized 30 0.998 – 0.015 [0.0013, 0.071] 4 Isolated 30 0.961 – 0.194 [0.159, 0.214] Cordial learning 30 0.999 0.903 0.0016 [0.00089, 0.0032] Centralized 30 (6) 0.999 – 0.0017 [<10−4<10^-4, 0.015] 6 Isolated 30 0.953 – 0.208 [0.189, 0.216] Cordial learning 30 1.000 0.882 0.00032 [0.00022, 0.00068] Centralized 30 (7) 1.000 – <10−4<10^-4 [<10−4<10^-4, 0.023] 500 2 Isolated 30 0.981 – 0.069 [0.059, 0.073] Cordial learning 30 0.998 0.954 0.0078 [0.0066, 0.011] Centralized 30 0.999 – 0.0024 [0.0017, 0.0036] 4 Isolated 30 0.982 – 0.063 [0.062, 0.067] Cordial learning 30 1.000 0.910 0.00085 [0.00029, 0.0018] Centralized 30 1.000 – <10−4<10^-4 [<10−4<10^-4, 0.0016] 6 Isolated 30 0.981 – 0.066 [0.064, 0.069] Cordial learning 30 1.000 0.916 <10−4<10^-4 [<10−4<10^-4, 0.00013] Centralized 30 1.000 – <10−4<10^-4 [<10−4<10^-4, <10−4<10^-4] Table 2: Pairwise-10: test accuracy (averaged over the agents) and test loss per agent of the network with the lowest validation loss, median over the runs, with the interquartile range of the loss in brackets. Runs is the number of runs: 30 for P=25P=25, 100 and 500 (three draws of the images, 10 random seeds each); the number of runs that did not learn (test loss ln9 9) is given in parentheses. For cordial agents, messages removed is the test accuracy when the message term is removed. 6.2.2 The Common-Cause Datasets In Pairwise-10, the digit of one agent determines the digits of all others. In the common-cause datasets, the dependence is probabilistic and much weaker. A hidden class z is drawn uniformly from 0,…,9\0,…,9\, and the digit of each agent is independently equal to z with probability p and uniform on the nine other digits otherwise. The digit of every agent is therefore uniform, and the digits of two agents are equal with probability p2+(1−p)2/9p^2+(1-p)^2/9. For p=0.1p=0.1 the digits are independent, and as p grows, the other agents carry more information about the digit of agent i. We report p=0.8p=0.8 and p=0.6p=0.6. As a reference for this information, the best prediction of the digit of an agent from the true digits of the other N−1N-1 agents alone has an accuracy of 0.640.64, 0.740.74 and 0.790.79 for N=2,4,6N=2,4,6 with p=0.8p=0.8, and 0.380.38, 0.450.45 and 0.540.54 with p=0.6p=0.6. Each of the three sets has 9P9P composites, the same number as in Pairwise-10. All hyperparameters were tuned separately for each value of p and P, as described above. The final runs have 1,500, 900, and 180 epochs for p=0.8p=0.8 and 2,250, 900, and 180 epochs for p=0.6p=0.6, for P=25P=25, 100, and 500. Results. Figures 6 and 7 show the test loss during training for p=0.8p=0.8 and p=0.6p=0.6, and Table 3 lists the test loss and accuracy of all methods. For p=0.8p=0.8 and P≥100P≥ 100, CL is better than isolated training in every setting and on every draw of the images. It has the lower test loss in 256 to 300 of the 300 pairs of runs for P=100P=100 and P=500P=500. Its median test loss is 18% to 35% lower, and its error rate is 15% to 28% lower than that of isolated agents. The gain is larger with 4 and 6 agents than with 2 agents. For P=500P=500 and 6 agents the test loss drops from 0.0680.068 to 0.0440.044 and the accuracy rises from 0.9820.982 to 0.9870.987. For P=25P=25, CL is better with 6 agents (274 of 300 pairs, a 9% lower loss); with 2 and 4 agents, the two methods are on par, where CL wins in 155 and 191 of 300 pairs, respectively. For p=0.6p=0.6, the less correlated scenario, the gain of CL depends on the amount of data. For P=500P=500, CL is better on every draw, in 246 to 282 of 300 pairs, with a 7% to 13% lower test loss. For P=100P=100, it is better with 4 agents (210 of 300 pairs) and on par with 2 and 6 agents. For P=25P=25, CL is worse than isolated training with 4 agents on every draw (78 of 300 pairs) and with 6 agents on two of the three draws (98 of 300 pairs), and on par with 2 agents. A likely reason is that the coordination layers overfit. With 225 training composites, the embedded output of a neighbor may identify its digit almost perfectly on the training images but not on new images, and the dependence between the digits is too weak to make up for these errors. Tuning selected the latest start of the coordination layers in its grid (K=400K=400 of 2,250 epochs), which suggests that on the validation set, the received embedded outputs help little. When the message term is removed, the accuracy of cordial agents drops in every setting with P≥100P≥ 100, and in 7 of these 12 settings it falls below that of isolated agents. For example, with p=0.8p=0.8, accuracy drops from 0.9850.985 to 0.9770.977 for P=500P=500 and 4 agents, against 0.9800.980 for isolated agents. Unlike in Pairwise-10, where the received embedded outputs can replace the agent’s own image, here they presumably correct the predictions on images whose digit is ambiguous. For P=500P=500, the centralized network is the best method and its test loss is 19% to 36% lower than that of CL. However, with P=25P=25, it is the worst method, with a test loss of 0.390.39 to 0.970.97, against 0.300.30 to 0.380.38 for isolated agents and 0.310.31 to 0.420.42 for CL. Its validation loss is lowest after a median of 96 to 242 epochs, while CL reaches its lowest validation loss near the end of training (after a median of 1,314 to 2,144 epochs), so the centralized network presumably overfits the 225 training composites. For P=100P=100, 18 of its 90 runs with p=0.8p=0.8 and 2 of its 90 runs with p=0.6p=0.6 did not learn at all and stayed at the uniform prediction (test loss ln10 10), as in Pairwise-10. CL has a lower median test loss than the centralized network in all settings with P=25P=25, for p=0.8p=0.8 with P=100P=100 and 4 or 6 agents, and for p=0.6p=0.6 with P=100P=100 and 6 agents. Unlike in Pairwise-10, this is not due to the runs that did not learn: without them, the median test loss of the centralized network for P=100P=100 is 0.160.16 and 0.300.30 for p=0.8p=0.8 with 4 and 6 agents, and 0.260.26 for p=0.6p=0.6 with 6 agents. Figure 6: Common-cause dataset with p=0.8p=0.8: test loss per agent during training of isolated agents, cordial agents and the centralized network. The lines show the median over 30 runs (10 random seeds on each of 3 draws) and the shaded bands the interquartile range. Figure 7: Common-cause dataset with p=0.6p=0.6: test loss per agent during training, as in Figure 6. Test loss per agent Test accuracy p P N Runs ISO CL Cen. ISO CL CL, removed Cen. CL << ISO 0.8 25 2 30 0.350 0.350 0.395 0.910 0.910 0.897 0.900 52% 4 30 0.298 0.310 0.827 0.914 0.911 0.891 0.771 64% 6 30 0.360 0.328 0.952 0.904 0.906 0.862 0.762 91% 100 2 30 0.203 0.166 0.122(1) 0.952 0.961 0.955 0.971 85% 4 30 0.196 0.147 0.175(8) 0.956 0.967 0.958 0.961 99% 6 30 0.172 0.119 0.319(9) 0.961 0.971 0.957 0.913 99.7% 500 2 30 0.080 0.061 0.039 0.980 0.983 0.978 0.989 97% 4 30 0.079 0.051 0.035 0.980 0.985 0.977 0.991 100% 6 30 0.068 0.044 0.036 0.982 0.987 0.978 0.990 100% 0.6 25 2 30 0.334 0.341 0.428 0.912 0.909 0.908 0.878 52% 4 30 0.381 0.422 0.826 0.892 0.876 0.878 0.759 26% 6 30 0.319 0.326 0.975 0.916 0.908 0.899 0.703 33% 100 2 30 0.208 0.213 0.134 0.953 0.958 0.957 0.968 43% 4 30 0.180 0.158 0.130(1) 0.959 0.965 0.963 0.971 70% 6 30 0.191 0.195 0.257(1) 0.952 0.954 0.950 0.936 53% 500 2 30 0.080 0.074 0.049 0.978 0.980 0.979 0.987 82% 4 30 0.072 0.067 0.051 0.979 0.982 0.979 0.987 93% 6 30 0.072 0.063 0.050 0.981 0.982 0.979 0.988 94% Table 3: Common-cause datasets: median test loss per agent and median test accuracy (averaged over the agents) of isolated agents (ISO), cordial learning (CL) and the centralized network (Cen.). Runs is the number of runs per method (three draws of the images for P=25P=25, 100 and 500). CL, removed is the accuracy of cordial agents when the message term is removed. CL << ISO is the fraction of pairs of runs on the same draw in which cordial learning has the lower test loss. A superscript gives the number of runs of the centralized network that did not learn (test loss was ln10 10). 7 Conclusions We studied a distributed learning scenario in which agents’ data are correlated. This generalizes the federated learning scenario where agents have independent data. To allow agents to learn from each other, we propose cordial learning, where agents send their embedded outputs to others. The agents train layers to learn how to extract valuable information from the embedded outputs of other agents. For a linear model where the label of an agent is a linear mixture of the inputs of all agents, under mild assumptions, we proved that cordial learning converges to models that minimize the sum of losses of all agents. Experiments on multi-digit MNIST with deep neural networks show that cordial learning performs well even beyond the theoretical model we analyzed. Our analysis focused on asymptotic convergence to the optimal models. The question of analyzing and optimizing the convergence rate remains open. The related question of proving optimal regret bounds in online learning settings is also open. Cordial learning is a first step towards studying distributed learning scenarios with correlated data. A challenging open question from this perspective is designing a cordial learning scheme where agents learn how to efficiently compress the information they send based on the loss functions of the receiving agents. Acknowledgments This research was supported by the Koret Foundation grant for Smart Cities and Digital Living 2030. Appendix A Proof of Lemma 3 We start by writing the stochastic loss: Li=(Yi−MiZ~i,0i−∑j=1,j≠iNΛjiMjZ~j,0i)T(Yi−MiZ~i,0i−∑j=1,j≠iNΛjiMjZ~j,0i)=YiTYi−YiTMiZ~i,0i−YiT∑j=1,j≠iNΛjiMjZ~j,0i−(MiZ~i,0i)TYi+(MiZ~i,0i)TMiZ~i,0i+(MiZ~i,0i)T∑j=1,j≠iNΛjiMjZ~j,0i−(∑j=1,j≠iNΛjiMjZ~j,0i)TYi+(∑j=1,j≠iNΛjiMjZ~j,0i)TMiZ~i,0i+(∑j=1,j≠iNΛjiMjZ~j,0i)T∑j=1,j≠iNΛjiMjZ~j,0i.L_i= (Y_i-M_i Z_i,0^i- _j=1,j≠ i^N _j^iM_j Z_j,0^i )^T (Y_i-M_i Z_i,0^i- _j=1,j≠ i^N _j^iM_j Z_j,0^i )=\\ Y_i^TY_i-Y_i^TM_i Z_i,0^i-Y_i^T _j=1,j≠ i^N _j^iM_j Z_j,0^i-(M_i Z_i,0^i)^TY_i+(M_i Z_i,0^i)^TM_i Z_i,0^i\\ +(M_i Z_i,0^i)^T _j=1,j≠ i^N _j^iM_j Z_j,0^i- ( _j=1,j≠ i^N _j^iM_j Z_j,0^i )^TY_i+ ( _j=1,j≠ i^N _j^iM_j Z_j,0^i )^TM_i Z_i,0^i\\ + ( _j=1,j≠ i^N _j^iM_j Z_j,0^i )^T _j=1,j≠ i^N _j^iM_j Z_j,0^i. (27) Thus the gradient of LiL_i with respect to MiM_i is ∂Li∂Mi ∂ L_i∂ M_i =∂Mi−2YiTMiZ~i,0i+(MiZ~i,0i)TMiZ~i,0i+2(∑j=1,j≠iNΛjiMjZ~j,0i)TMiZ~i,0i = ∂ M_i \-2Y_i^TM_i Z^i_i,0+(M_i Z^i_i,0)^TM_i Z^i_i,0+2 ( _j=1,j≠ i^N ^i_jM_j Z^i_j,0 )^TM_i Z^i_i,0 \ =−2YiZ~i,0iT+2MiZ~i,0iZ~i,0iT+2(∑j=1,j≠iNΛjiMjZ~j,0i)Z~i,0iT, =-2Y_i Z^iT_i,0+2M_i Z^i_i,0 Z^iT_i,0+2 ( _j=1,j≠ i^N ^i_jM_j Z^i_j,0 ) Z^iT_i,0, (28) where we used the derivative rules ∂(YTMX)∂M=YXT ∂(Y^TMX)∂ M=YX^T and ∂(XTATAX)∂A=2AXXT ∂(X^TA^TAX)∂ A=2AXX^T. The stochastic gradient Mi g_M_i in (16) is (28) with Λji ^i_j replaced by its projection Πt[Λji(t)] _t [ ^i_j(t) ]. The gradient of LiL_i with respect to Λji ^i_j is ∂Li∂Λji=−2(Yi−MiZ~i,0i−∑l=1,l≠iNΛliMlZ~l,0i)(MjZ~j,0i)T. ∂ L_i∂ ^i_j=-2 (Y_i-M_i Z^i_i,0- _l=1,l≠ i^N ^i_lM_l Z^i_l,0 ) (M_j Z^i_j,0 )^T. Algorithm 1 replaces Z~j,0i Z^i_j,0 in the last factor by Z~j,1i Z^i_j,1, which is computed from the independent batch Bt1B^1_t. This gives Λji=−2Yi(MjZ~j,1i)T+2MiZ~i,0i(MjZ~j,1i)T+2(∑l=1,l≠iNΛliMlZ~l,0i)(MjZ~j,1i)T. g_ ^i_j=-2Y_i(M_j Z^i_j,1)^T+2M_i Z^i_i,0(M_j Z^i_j,1)^T+2 ( _l=1,l≠ i^N ^i_lM_l Z^i_l,0 )(M_j Z^i_j,1)^T. (29) With this replacement, the two covariance estimates in each product are independent, which makes (29) unbiased (see (32)). The expected gradients are as follows: Mi=2(Mi−Ai)Z~i,0iZ~i,0iT+2∑j=1,j≠iN(ΛjiMjZ~j,0i−BjiXj)Z~i,0iT=2(Mi−Ai−∑j=1,j≠iNBjiCov(Xj,Xi))E \ g_M_i \=E \2 (M_i-A_i ) Z_i,0^i Z_i,0^iT+2 _j=1,j≠ i^N ( _j^iM_j Z_j,0^i-B_j^iX_j ) Z_i,0^iT \=\\ 2 (M_i-A_i- _j=1,j≠ i^NB_j^i (X_j,X_i ) ) (30) where we used Z~j,0iXiT=(Xj−R~j,0iXi)XiT=XjXiT−R~j,0iXiXiT=0,E \ Z_j,0^iX_i^T \=E \ (X_j- R_j,0^iX_i )X_i^T \=E\X_jX_i^T\-E\ R_j,0^i\E\X_iX_i^T\=0, and XjZ~i,0iT=Cov(Xj,Xi)E \X_j Z_i,0^iT \= (X_j,X_i ). Next, Λji \ g_ _j^i \ =(2(Mi−Ai)Xi+2∑k=1,k≠iN(ΛkiMkZ~k,0i−BkiXk))(Z~j,1i)TMjT =E \ (2(M_i-A_i)X_i+2 _k=1,k≠ i^N ( _k^iM_k Z_k,0^i-B_k^iX_k ) ) ( Z_j,1^i )^TM_j^T \ =2∑k=1,k≠iN(ΛkiMk−Bki)Ck,jMjT, =2 _k=1,k≠ i^N ( _k^iM_k-B_k^i )C_k,jM_j^T, (31) where we used Z~k,0i(Z~j,1i)T=(Xk−R~k,0iXi)(Xj−R~j,1iXi)T=XkXjT−XkXiT(R~j,1i)T−R~k,0iXiXjT+R~k,0iXiXiT(R~j,1i)T=Xk(Xj−R~j,1iXi)T=XkXjT−Cov(Xk,Xi)Cov(Xj,Xi)T=Ck,j.E \ Z_k,0^i ( Z_j,1^i )^T \=E \ (X_k- R_k,0^iX_i ) (X_j- R_j,1^iX_i )^T \=\\ E\X_kX_j^T\-E\X_kX_i^T\E\( R_j,1^i)^T\-E\ R_k,0^i\E\X_iX_j^T\+E\ R_k,0^i\E\X_iX_i^T\E\( R_j,1^i)^T\=\\ E \X_k (X_j- R_j,1^iX_i )^T \=E\X_kX_j^T\-Cov(X_k,X_i)Cov(X_j,X_i)^T=C_k,j. (32) Here we used that R~k,0i R_k,0^i is computed using the leave-one-out samples from Bt0B_t^0, while R~j,1i R_j,1^i is computed using the independent batch Bt1B_t^1. Hence, R~k,0i R_k,0^i, R~j,1i R_j,1^i, and the given sample are mutually independent. Moreover, R~k,0i=Cov(Xk,Xi),R~j,1i=Cov(Xj,Xi).E\ R_k,0^i\=Cov(X_k,X_i), \ R_j,1^i\=Cov(X_j,X_i). Appendix B Nonconvex Objective We show that the global objective in (7) is nonconvex, even in the two-agent scalar case. Let d=m=1d=m=1, and let X1X_1 and X2X_2 be zero-mean random variables satisfying [X12]=[X22]=1E[X_1^2]=E[X_2^2]=1 and [X1X2]=0E[X_1X_2]=0. Consider the noiseless labels Y1=X2Y_1=X_2 and Y2=X2Y_2=X_2. Then A1∗=0A_1^*=0 and A2∗=1A_2^*=1, so Assumption 1 is satisfied by taking Λ21∗=1 _2^1*=1 and Λ12∗=0 _1^2*=0. For two scalar agents, the global objective is F(M1,M2,Λ1,Λ2)= F(M_1,M_2, _1, _2)= [(Y1−M1X1−Λ1M2X2)2] \! [ (Y_1-M_1X_1- _1M_2X_2 )^2 ] +[(Y2−M2X2−Λ2M1X1)2] +E\! [ (Y_2-M_2X_2- _2M_1X_1 )^2 ] where Λ1=Λ21 _1= ^1_2 and Λ2=Λ12 _2= ^2_1. Next we restrict this objective to the affine subspace M1=0M_1=0 and Λ2=0 _2=0. On this subspace, h(M2,Λ1) h(M_2, _1) ≜F(0,M2,Λ1,0) F(0,M_2, _1,0) =[(X2−Λ1M2X2)2]+[(X2−M2X2)2] =E\! [ (X_2- _1M_2X_2 )^2 ]+E\! [ (X_2-M_2X_2 )^2 ] =(1−Λ1M2)2+(1−M2)2. =(1- _1M_2)^2+(1-M_2)^2. The Hessian of h is ∇2h(M2,Λ1)=(2+2Λ12−2+4Λ1M2−2+4Λ1M22M22).∇^2h(M_2, _1)= pmatrix2+2 _1^2&-2+4 _1M_2\\ -2+4 _1M_2&2M_2^2 pmatrix. At (M2,Λ1)=(0,0)(M_2, _1)=(0,0), ∇2h(0,0)=(2−2−20),∇^2h(0,0)= pmatrix2&-2\\ -2&0 pmatrix, whose determinant is −4-4. Hence, ∇2h(0,0)∇^2h(0,0) is indefinite, and h is not convex. Since the restriction of a convex function to any affine subspace must be convex, the full objective F(M1,M2,Λ1,Λ2)F(M_1,M_2, _1, _2) is not convex. Therefore, the global objective in (7) is nonconvex. Appendix C Proof of Theorem 2 We start by proving the convergence of Mi(t)M_i(t). We have ∑i=1N‖Ai∗−Mi(t+1)‖2∣ℱt=∑i=1N‖Ai∗−Mi(t)+ε(t)Mi(t)‖2∣ℱt=∑i=1N‖Ai∗−Mi(t)‖2+2ε(t)∑i=1N⟨Mi(t),Ai∗−Mi(t)⟩∣ℱt+ε2(t)∑i=1N‖Mi(t)‖2∣ℱt≤(a)(1+K1ε2(t)ρ2(t))∑i=1N‖Ai∗−Mi(t)‖2−4ε(t)∑i=1N‖Ai∗−Mi(t)‖2+K2ε2(t)ρ2(t)E \ _i=1^N \|A_i^*-M_i (t+1 ) \|^2 _t \=E \ _i=1^N \|A_i^*-M_i (t )+ (t ) g_M_i (t ) \|^2 _t \=\\ _i=1^N \|A_i^*-M_i (t ) \|^2+2 (t ) _i=1^NE \ g_M_i (t ),A_i^*-M_i (t ) _t \+ ^2 (t ) _i=1^NE \ \| g_M_i (t ) \|^2 _t \ (a )≤\\ (1+K_1 ^2 (t )ρ^2 (t ) ) _i=1^N \|A_i^*-M_i (t ) \|^2-4 (t ) _i=1^N \|A_i^*-M_i (t ) \|^2+K_2 ^2 (t )ρ^2 (t ) (33) where (a) uses Lemma 3 and Lemma 4. We also used that Mi(t)M_i(t) is measurable by ℱtF_t. Given the inequality above, we can apply the Robbins-Siegmund Theorem (Robbins and Siegmund, 1971) with (using their notation) zt=∑i=1N‖Ai∗−Mi(t)‖2z_t= _i=1^N \|A_i^*-M_i (t ) \|^2, βt=K1ε2(t)ρ2(t) _t=K_1 ^2 (t )ρ^2(t), ξt=K2ε2(t)ρ2(t) _t=K_2 ^2 (t )ρ^2(t) and ζt=4ε(t)∑i=1N‖Ai∗−Mi(t)‖2 _t=4 (t ) _i=1^N \|A_i^*-M_i (t ) \|^2. Hence, with probability 1, for all i: ∑t=1∞ε(t)∑i=1N‖Ai∗−Mi(t)‖2<∞ _t=1^∞ (t ) _i=1^N \|A_i^*-M_i (t ) \|^2<∞ (34) and ∑i=1N‖Ai∗−Mi(t)‖2 _i=1^N \|A_i^*-M_i (t ) \|^2 converges to some random variable, which must be zero with probability 1, since any positive probability to be positive would contradict (34), since ∑t=1∞ε(t)=∞ _t=1^∞ (t )=∞ by assumption. Now let fi(Λ1i,…,ΛNi)≜‖∑k=1,k≠iN(ΛkiAk∗−Bki)Zki‖2f_i ( _1^i,…, _N^i ) \ \| _k=1,k≠ i^N ( _k^iA_k^*-B_k^i )Z_k^i \|^2 \. We have for some constant K0>0K_0>0 ∑j≠i‖Λji(t+1)−Λji∗‖2∣ℱt=∑j≠i‖Λji(t)−η(t)Λji(t)−Λji∗‖2∣ℱt=∑j≠i‖Λji(t)−Λji∗‖2−2η(t)∑j≠i⟨Λji(t),Λji(t)−Λji∗⟩∣ℱt+η2(t)∑j≠i‖Λji(t)‖2∣ℱt≤(a)∑j≠i‖Λji(t)−Λji∗‖2−4η(t)∑j≠i⟨∑k≠i(Λki(t)Mk(t)−Bki)Ck,jMjT(t),Λji(t)−Λji∗⟩+K0η2(t)∑l=1N(1+‖Λli(t)‖2)∑k=1N(1+2‖Mk(t)−Ak∗‖2+‖Mk(t)−Ak∗‖4)E \ _j≠ i \| _j^i (t+1 )- _j^i* \|^2 _t \=E \ _j≠ i \| _j^i (t )-η (t ) g_ _j^i (t )- _j^i* \|^2 _t \=\\ _j≠ i \| _j^i (t )- _j^i* \|^2-2η (t )E \ _j≠ i g_ _j^i (t ), _j^i (t )- _j^i* _t \+η^2 (t )E \ _j≠ i \| g_ _j^i (t ) \|^2 _t \ (a )≤\\ _j≠ i \| _j^i (t )- _j^i* \|^2-4η (t ) _j≠ i _k≠ i ( _k^i (t )M_k (t )-B_k^i )C_k,jM_j^T (t ), _j^i (t )- _j^i* \\ +K_0η^2 (t ) _l=1^N (1+ \| _l^i (t ) \|^2 ) _k=1^N (1+2 \|M_k (t )-A_k^* \|^2+ \|M_k (t )-A_k^* \|^4 ) (35) where (a) follows from Lemma 3 and Lemma 4. To deal with the bias term in (35), we write: ∑j≠i⟨∑k≠i(Λki(t)Mk(t)−Bki)Ck,jMjT(t),Λji(t)−Λji∗⟩=∑j≠i⟨∑k≠i(Λki(t)Ak∗−Bki)Ck,jAj∗T,Λji(t)−Λji∗⟩+∑j≠i⟨∑k≠iBkiCk,j(Aj∗T−MjT(t)),Λji(t)−Λji∗⟩+∑j≠i⟨∑k≠iΛki(t)(Mk(t)Ck,jMjT(t)−Ak∗Ck,jAj∗T),Λji(t)−Λji∗⟩. _j≠ i _k≠ i ( _k^i (t )M_k (t )-B_k^i )C_k,jM_j^T (t ), _j^i (t )- _j^i* =\\ _j≠ i _k≠ i ( _k^i (t )A_k^*-B_k^i )C_k,jA_j^*T, _j^i (t )- _j^i* +\\ _j≠ i _k≠ iB_k^iC_k,j (A_j^*T-M_j^T (t ) ), _j^i (t )- _j^i* +\\ _j≠ i _k≠ i _k^i (t ) (M_k (t )C_k,jM_j^T (t )-A_k^*C_k,jA_j^*T ), _j^i (t )- _j^i* . (36) For the first term in (36), we use the convexity of fif_i to obtain ∑j≠i⟨∑k=1,k≠iN(Λki(t)Ak∗−Bki)Ck,jAj∗T,Λji(t)−Λji∗⟩=∑j≠i⟨12∇Λjifi(Λ1i(t),…,ΛNi(t)),Λji(t)−Λji∗⟩≥12fi(Λ1i(t),…,ΛNi(t))−12fi(Λ1i∗,…,ΛNi∗)=12fi(Λ1i(t),…,ΛNi(t)) _j≠ i _k=1,k≠ i^N ( _k^i (t )A_k^*-B_k^i )C_k,jA_j^*T, _j^i (t )- _j^i* =\\ _j≠ i 12 _ _j^if_i ( _1^i (t ),…, _N^i (t ) ), _j^i (t )- _j^i* ≥\\ 12f_i ( _1^i (t ),…, _N^i (t ) )- 12f_i ( _1^i*,…, _N^i* )= 12f_i ( _1^i (t ),…, _N^i (t ) ) (37) since fi(Λ1i∗,…,ΛNi∗)=0f_i ( _1^i*,…, _N^i* )=0 by definition of Λji∗ _j^i*. For the second term in (36) we use: −∑j≠i⟨∑k≠iBkiCk,j(Aj∗T−MjT(t)),Λji(t)−Λji∗⟩≤(a)∑j≠i∑k≠i‖BkiCk,j(Aj∗T−MjT(t))‖Λji(t)−Λji∗‖≤C5∑j≠i‖Mj(t)−Aj∗‖Λji(t)−Λji∗‖≤(b)C5∑j≠i‖Mj(t)−Aj∗‖Λji(t)−Λji∗‖2+C5∑j≠i‖Mj(t)−Aj∗‖- _j≠ i _k≠ iB_k^iC_k,j (A_j^*T-M_j^T (t ) ), _j^i (t )- _j^i* (a )≤\\ _j≠ i _k≠ i \|B_k^iC_k,j (A_j^*T-M_j^T (t ) ) \| \| _j^i (t )- _j^i* \|≤\\ C_5 _j≠ i \|M_j (t )-A_j^* \| \| _j^i (t )- _j^i* \| (b )≤\\ C_5 _j≠ i \|M_j (t )-A_j^* \| \| _j^i (t )- _j^i* \|^2+C_5 _j≠ i \|M_j (t )-A_j^* \| (38) where (a) (a ) follows from the Cauchy-Schwarz inequality and the triangle inequality, and (b) (b ) uses x≤x2+1x≤ x^2+1. For the third term in (36) we use: −∑j≠i⟨∑k≠iΛki(t)(Mk(t)Ck,jMjT(t)−Ak∗Ck,jAj∗T),Λji(t)−Λji∗⟩≤(a)∑j≠i‖Λji(t)−Λji∗‖∑k≠i‖Λki(t)(Mk(t)−Ak∗)Ck,jAj∗T‖+∑j≠i‖Λji(t)−Λji∗‖∑k≠i‖Λki(t)Mk(t)Ck,j(MjT(t)−Aj∗T)‖≤∑j≠i∑k≠i‖(Λki(t)−Λki∗)(Mk(t)−Ak∗)Ck,jAj∗T‖Λji(t)−Λji∗‖+∑j≠i∑k≠i‖Λki∗(Mk(t)−Ak∗)Ck,jAj∗T‖Λji(t)−Λji∗‖+∑j≠i∑k≠i‖(Λki(t)−Λki∗)Mk(t)Ck,j(MjT(t)−Aj∗T)‖Λji(t)−Λji∗‖+∑j≠i∑k≠i‖Λki∗Mk(t)Ck,j(MjT(t)−Aj∗T)‖Λji(t)−Λji∗‖≤(b)C1∑k=1N‖Mk(t)−Ak∗‖∑l≠i‖Λli(t)−Λli∗‖+C2∑k=1N‖Mk(t)−Ak∗‖∑l≠i‖Λli(t)−Λli∗‖2+C3∑k=1N‖Mk(t)−Ak∗‖2∑l≠i‖Λli(t)−Λli∗‖2+C4∑k=1N‖Mk(t)−Ak∗‖2∑l≠i‖Λli(t)−Λli∗‖≤C5∑k=1N‖Mk(t)−Ak∗‖∑l≠i‖Λli(t)−Λli∗‖2+C1∑k=1N‖Mk(t)−Ak∗‖+C6∑k=1N‖Mk(t)−Ak∗‖2∑l≠i‖Λli(t)−Λli∗‖2+C4∑k=1N‖Mk(t)−Ak∗‖2- _j≠ i _k≠ i _k^i (t ) (M_k (t )C_k,jM_j^T (t )-A_k^*C_k,jA_j^*T ), _j^i (t )- _j^i* (a )≤\\ _j≠ i \| _j^i (t )- _j^i* \| _k≠ i \| _k^i (t ) (M_k (t )-A_k^* )C_k,jA_j^*T \|\\ + _j≠ i \| _j^i (t )- _j^i* \| _k≠ i \| _k^i (t )M_k (t )C_k,j (M_j^T (t )-A_j^*T ) \|≤\\ _j≠ i _k≠ i \| ( _k^i (t )- _k^i* ) (M_k (t )-A_k^* )C_k,jA_j^*T \| \| _j^i (t )- _j^i* \|\\ + _j≠ i _k≠ i \| _k^i* (M_k (t )-A_k^* )C_k,jA_j^*T \| \| _j^i (t )- _j^i* \|\\ + _j≠ i _k≠ i \| ( _k^i (t )- _k^i* )M_k (t )C_k,j (M_j^T (t )-A_j^*T ) \| \| _j^i (t )- _j^i* \|\\ + _j≠ i _k≠ i \| _k^i*M_k (t )C_k,j (M_j^T (t )-A_j^*T ) \| \| _j^i (t )- _j^i* \| (b )≤\\ C_1 _k=1^N \|M_k (t )-A_k^* \| _l≠ i \| _l^i (t )- _l^i* \|+C_2 _k=1^N \|M_k (t )-A_k^* \| _l≠ i \| _l^i (t )- _l^i* \|^2+\\ C_3 _k=1^N \|M_k (t )-A_k^* \|^2 _l≠ i \| _l^i (t )- _l^i* \|^2+C_4 _k=1^N \|M_k (t )-A_k^* \|^2 _l≠ i \| _l^i (t )- _l^i* \|≤\\ C_5 _k=1^N \|M_k (t )-A_k^* \| _l≠ i \| _l^i (t )- _l^i* \|^2+C_1 _k=1^N \|M_k (t )-A_k^* \|+\\ C_6 _k=1^N \|M_k (t )-A_k^* \|^2 _l≠ i \| _l^i (t )- _l^i* \|^2+C_4 _k=1^N \|M_k (t )-A_k^* \|^2 (39) where (a) (a ) follows from the Cauchy-Schwarz inequality and the triangle inequality. In (b) we used that, for al=∥Λli(t)−Λli∗∥a_l=\| ^i_l(t)- ^i*_l\| (l≠il≠ i), ai=0a_i=0, bl=‖Ml(t)−Al∗‖b_l=\|M_l(t)-A^*_l\|, Sa=∑lalS_a= _la_l, Sb=∑lblS_b= _lb_l, Sa2=∑lal2S_a^2= _la_l^2 and Sb2=∑lbl2S_b^2= _lb_l^2, Cauchy-Schwarz and the L1-L2 norm inequality yield ∑k,lakalbk+∑k,lakblbk+∑k,lakbl+∑k,lakalbkbl=(Sa+Sb)∑kakbk+SaSb+(∑kakbk)2≤(Sa+Sb)Sa2Sb2+SaSb+Sa2Sb2≤N(Sa2+Sb2)Sa2Sb2+SaSb+Sa2Sb2≤NSa2Sb+NSaSb2+SaSb+Sa2Sb2. _k,la_ka_lb_k+ _k,la_kb_lb_k+ _k,la_kb_l+ _k,la_ka_lb_kb_l=\\ (S_a+S_b) _ka_kb_k+S_aS_b+ ( _ka_kb_k )^2≤\\ (S_a+S_b) S_a^2S_b^2+S_aS_b+S_a^2S_b^2≤ N ( S_a^2+ S_b^2 ) S_a^2S_b^2+S_aS_b+S_a^2S_b^2≤\\ NS_a^2S_b+ NS_aS_b^2+S_aS_b+S_a^2S_b^2. (40) Plugging (36), (37), (39) and (38) into (35), we obtain ∑j≠i‖Λji(t+1)−Λji∗‖2∣ℱt≤∑j≠i‖Λji(t)−Λji∗‖2−2η(t)fi(Λ1i(t),…,ΛNi(t))+C1∑k=1Nη(t)(‖Mk(t)−Ak∗‖+‖Mk(t)−Ak∗‖2+η(t)(‖Mk(t)−Ak∗‖4+1))∑j≠i‖Λji(t)−Λji∗‖2+C2∑k=1Nη(t)(‖Mk(t)−Ak∗‖+‖Mk(t)−Ak∗‖2+η(t)(‖Mk(t)−Ak∗‖4+1)).E \ _j≠ i \| _j^i (t+1 )- _j^i* \|^2 _t \≤\\ _j≠ i \| _j^i (t )- _j^i* \|^2-2η (t )f_i ( _1^i (t ),…, _N^i (t ) )\\ +C_1 _k=1^Nη (t ) ( \|M_k (t )-A_k^* \|+ \|M_k (t )-A_k^* \|^2+η (t ) ( \|M_k (t )-A_k^* \|^4+1 ) ) _j≠ i \| _j^i (t )- _j^i* \|^2\\ +C_2 _k=1^Nη (t ) ( \|M_k (t )-A_k^* \|+ \|M_k (t )-A_k^* \|^2+η (t ) ( \|M_k (t )-A_k^* \|^4+1 ) ). (41) Given (41), we can again apply the Robbins-Siegmund Theorem with zt=∑j≠i‖Λji(t)−Λji∗‖2z_t= _j≠ i \| _j^i (t )- _j^i* \|^2, βt=C1∑k=1Nη(t)(‖Mk(t)−Ak∗‖+‖Mk(t)−Ak∗‖2+η(t)+η(t)‖Mk(t)−Ak∗‖4) _t=C_1 _k=1^Nη (t ) ( \|M_k (t )-A_k^* \|+ \|M_k (t )-A_k^* \|^2+η (t )+η (t ) \|M_k (t )-A_k^* \|^4 ) (42) ξt=C2∑k=1Nη(t)(‖Mk(t)−Ak∗‖+‖Mk(t)−Ak∗‖2+η(t)+η(t)‖Mk(t)−Ak∗‖4) _t=C_2 _k=1^Nη (t ) ( \|M_k (t )-A_k^* \|+ \|M_k (t )-A_k^* \|^2+η (t )+η (t ) \|M_k (t )-A_k^* \|^4 ) (43) and ζt=2η(t)fi(Λ1i(t),…,ΛNi(t)). _t=2η (t )f_i ( _1^i (t ),…, _N^i (t ) ). We know from (34) that ∑t=1∞ε(t)‖Ai∗−Mi(t)‖2<∞ _t=1^∞ (t ) \|A_i^*-M_i (t ) \|^2<∞. Thus, using the assumption that ∑t=1∞η2(t)ε(t)<∞ _t=1^∞ η^2 (t ) (t )<∞: ∑t=1∞η(t)‖Ai∗−Mi(t)‖=∑t=1∞η(t)ε(t)ε(t)‖Ai∗−Mi(t)‖≤∑t=1∞η2(t)ε(t)+∑t=1∞ε(t)‖Ai∗−Mi(t)‖2<∞. _t=1^∞η (t ) \|A_i^*-M_i (t ) \|= _t=1^∞ η (t ) (t ) (t ) \|A_i^*-M_i (t ) \|≤\\ _t=1^∞ η^2 (t ) (t )+ _t=1^∞ (t ) \|A_i^*-M_i (t ) \|^2<∞. (44) Then, ∑t=1∞η(t)‖Ai∗−Mi(t)‖2 _t=1^∞η (t ) \|A_i^*-M_i (t ) \|^2 and ∑t=1∞η2(t)‖Ai∗−Mi(t)‖4 _t=1^∞η^2 (t ) \|A_i^*-M_i (t ) \|^4 are also finite since ‖Ai∗−Mi(t)‖ \|A_i^*-M_i (t ) \| converges to zero with probability 1, and ∑t=1∞η2(t)≤(suptε(t))∑t=1∞η2(t)ε(t)<∞ _t=1^∞η^2 (t )≤ ( _t (t ) ) _t=1^∞ η^2 (t ) (t )<∞. Hence, the Robbins-Siegmund Theorem implies that with probability 1, for all i: ∑t=1∞η(t)fi(Λ1i(t),…,ΛNi(t))<∞ _t=1^∞η (t )f_i ( _1^i (t ),…, _N^i (t ) )<∞ (45) and ∑j≠i‖Λji(t)−Λji∗‖2 _j≠ i \| _j^i (t )- _j^i* \|^2 converges to some random variable V as t→∞t→∞. In particular, i(t)t\ ^i(t)\_t is bounded. Define the closed set i≜i:fi(i)=0S_i \ ^i:f_i( ^i)=0 \ which is nonempty by Assumption 1. Let i0S_i^0 be a countable dense subset of iS_i. Applying the preceding Robbins-Siegmund argument to every fixed i∗∈i0 ^i* _i^0, and taking the intersection of the resulting probability-one events, shows that, with probability one, ‖i(t)−i∗‖ \| ^i(t)- ^i* \| converges for every i∗∈i0 ^i* _i^0. The same conclusion holds for every i∗∈i ^i* _i since for any i∗∈i ^i* _i and any sequence i∗,r∈i0 ^i*,r _i^0 converging to i∗ ^i*, |‖i(t)−i∗‖−‖i(t)−i∗,r‖|≤‖i∗−i∗,r‖, | \| ^i(t)- ^i* \|- \| ^i(t)- ^i*,r \| |≤ \| ^i*- ^i*,r \|, uniformly in t. Moreover, since ∑tη(t)=∞ _tη(t)=∞, (45) implies that lim inft→∞fi(i(t))=0. _t→∞f_i\! ( ^i(t) )=0. Consequently, there exists a subsequence trr\t_r\_r such that fi(i(tr))→0f_i\! ( ^i(t_r) )→ 0. By boundedness, this subsequence has a further subsequence, which we denote again by trr\t_r\_r, such that i(tr)→¯i ^i(t_r)→ ^i. Since fif_i is continuous, fi(¯i)=0f_i\! ( ^i )=0 and hence ¯i∈i ^i _i. As shown above, ‖i(t)−¯i‖\| ^i(t)- ^i\| has a limit. Along the subsequence trr\t_r\_r, it converges to zero. Therefore, its full-sequence limit is zero, and i(t)→¯i ^i(t)→ ^i with probability one. By continuity of fif_i, it follows that fi(i(t))→fi(¯i)=0f_i\! ( ^i(t) )→ f_i\! ( ^i )=0 with probability one. Together with Mk(t)→Ak∗M_k(t)→ A_k^* for every k, and the boundedness of Λki(t)t\ _k^i(t)\_t, this implies that the loss of each agent i, LiL_i, converges to its global minimum σ2σ^2 with probability one. References Bistritz and Bambos (2020) I. Bistritz and N. Bambos Cooperative multi-player bandit optimization. Advances in Neural Information Processing Systems 33. Cited by: §2. Bistritz et al. (2020) I. Bistritz, A. Mann, and N. Bambos Distributed distillation for on-device learning. Advances in Neural Information Processing Systems 33, p. 22593–22604. Cited by: §2. Bistritz et al. (2022) I. Bistritz, Z. Zhou, X. Chen, N. Bambos, and J. Blanchet No weighted-regret learning in adversarial bandits with delays. Journal of Machine Learning Research 23 (139), p. 1–43. Cited by: §2. Chasnov et al. (2020) B. Chasnov, L. Ratliff, E. Mazumdar, and S. Burden Convergence analysis of gradient-based learning in continuous games. In Uncertainty in Artificial Intelligence, p. 935–944. Cited by: §2. Goodfellow et al. (2020) I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio Generative adversarial networks. Communications of the ACM 63 (11), p. 139–144. Cited by: §2. Héliou et al. (2017) A. Héliou, J. Cohen, and P. Mertikopoulos Learning with bandit feedback in potential games. Advances in Neural Information Processing Systems 30. Cited by: §3.1. Hong et al. (2022) L. Hong, A. Garcia, and C. Eksin Distributed networked learning with correlated data. Automatica 137, p. 110134. Cited by: §2. Kairouz et al. (2021) P. Kairouz, H. B. McMahan, B. Avent, A. Bellet, M. Bennis, A. N. Bhagoji, K. Bonawitz, Z. Charles, G. Cormode, R. Cummings, et al. Advances and open problems in federated learning. Foundations and Trends® in Machine Learning 14 (1–2), p. 1–210. Cited by: §2. Koloskova et al. (2022) A. Koloskova, S. U. Stich, and M. Jaggi Sharper convergence guarantees for asynchronous SGD for distributed and federated learning. Advances in Neural Information Processing Systems 35, p. 17202–17215. Cited by: §2. Lin and Zrnic (2024) L. Lin and T. Zrnic Plug-in performative optimization. In Proceedings of the 41st International Conference on Machine Learning, R. Salakhutdinov, Z. Kolter, K. Heller, A. Weller, N. Oliver, J. Scarlett, and F. Berkenkamp (Eds.), Proceedings of Machine Learning Research, Vol. 235, p. 30546–30565. Cited by: §2. Liu et al. (2024) Y. Liu, Y. Kang, T. Zou, Y. Pu, Y. He, X. Ye, Y. Ouyang, Y. Zhang, and Q. Yang Vertical federated learning: concepts, advances, and challenges. IEEE Transactions on Knowledge and Data Engineering 36 (7), p. 3615–3634. Cited by: §2, §3. Maurer (2006) A. Maurer Bounds for linear multi-task learning. Journal of Machine Learning Research 7, p. 117–139. Cited by: §2, §4. Mazumdar et al. (2020) E. Mazumdar, L. J. Ratliff, and S. S. Sastry On gradient-based learning in continuous games. SIAM Journal on Mathematics of Data Science 2 (1), p. 103–131. Cited by: §2, §4. Mitra et al. (2021) A. Mitra, R. Jaafar, G. J. Pappas, and H. Hassani Linear convergence in federated learning: tackling client heterogeneity and sparse gradients. Advances in Neural Information Processing Systems 34, p. 14606–14619. Cited by: §2. Narang et al. (2023) A. Narang, E. Faulkner, D. Drusvyatskiy, M. Fazel, and L. J. Ratliff Multiplayer performative prediction: learning in decision-dependent games. Journal of Machine Learning Research 24 (202), p. 1–56. Cited by: §2. Perdomo et al. (2020) J. Perdomo, T. Zrnic, C. Mendler-Dünner, and M. Hardt Performative prediction. In International Conference on Machine Learning, p. 7599–7609. Cited by: §2. Ratliff et al. (2016) L. J. Ratliff, S. A. Burden, and S. S. Sastry On the characterization of local Nash equilibria in continuous games. IEEE Transactions on Automatic Control 61 (8), p. 2301–2307. Cited by: §2. Ren et al. (2015) S. Ren, K. He, R. Girshick, and J. Sun Faster R-CNN: towards real-time object detection with region proposal networks. Advances in Neural Information Processing Systems 28. Cited by: §2. Robbins and Siegmund (1971) H. Robbins and D. Siegmund A convergence theorem for non negative almost supermartingales and some applications. In Optimizing Methods in Statistics, p. 233–257. Cited by: Appendix C. Roughgarden (2015) T. Roughgarden Intrinsic robustness of the price of anarchy. Journal of the ACM (JACM) 62 (5), p. 1–42. Cited by: §4. Sessa et al. (2019) P. G. Sessa, M. Kamgarpour, and A. Krause Bounding inefficiency of equilibria in continuous actions games using submodularity and curvature. In The 22nd International Conference on Artificial Intelligence and Statistics, p. 2017–2027. Cited by: §4. Wang and Joshi (2021) J. Wang and G. Joshi Cooperative SGD: a unified framework for the design and analysis of local-update SGD algorithms. Journal of Machine Learning Research 22 (213), p. 1–50. Cited by: §2. Xu et al. (2025) Z. Xu, C. Yu, Y. Liang, Y. Wu, and Y. Wang Learning global Nash equilibrium in team competitive games with generalized fictitious cross-play. Journal of Machine Learning Research 26 (44), p. 1–30. Cited by: §2.