Paper deep dive
Efficient Gradient Methods for Distributed Saddle Problems
Ruichen Luo, Anton Rodomanov, Sebastian U. Stich
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 94%
Last extracted: 7/8/2026, 5:19:46 PM
Summary
This paper formalizes the distributed setting for Saddle Problems (SPs) and Variational Inequality Problems (VIPs), introducing a novel decoupled method (DM-SP) that achieves optimal communication cost within the zero-respecting framework. It improves upon the long-standing oracle complexity of the Extragradient method and establishes matching lower bounds for gradient-span algorithms, extending the results to multi-agent general-sum games.
Entities (8)
Relation Signals (6)
Ruichen Luo → affiliatedwith → Institute of Science and Technology Austria
confidence 99% · name=Ruichen Luo, affiliation=Institute of Science and Technology Austria
Saddle Problems (SPs) → generalizesto → Variational Inequality Problems (VIPs)
confidence 95% · Finally, we study the extension of distributed SP into Variational Inequality Problem (VIP), which generalizes two-player zero-sum games to multiplayer general-sum games.
Decoupled Method (DM-SP) → achievesoptimal → Communication Cost
confidence 94% · Our main result is a novel decoupled method that achieves optimal communication cost within the zero-respecting framework.
Decoupled Method (DM-SP) → improves → Extragradient Method (EG)
confidence 93% · Our method is based on a multi-stage reduction to the decoupled minimization of residual norms, which yields strict improvements over the best known communication cost for the class and the long-standing oracle cost of the Extragradient method.
Decoupled Method (DM-SP) → isoptimalwithin → Gradient-Span Algorithms
confidence 91% · Further, we show by a matching lower bound that our method is communication-optimal within the family of gradient-span algorithms.
Saddle Problems (SPs) → appliedin → Generative Adversarial Networks
confidence 90% · These problems have a wide array of modern applications, including the training of GANs (Goodfellow et al., 2014), robust optimization (Ben-Tal and Nemirovski, 2002), and equilibrium computation in game theory and multiagent systems
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:The distributed setting for Saddle Problems (SPs) has recently emerged as a framework for various modern applications in machine learning and multiagent systems. Despite its relevance, the theoretical foundations of this setting have not yet been thoroughly established. In this paper, we advance this research direction by formalizing the distributed setup for SPs and providing rigorous definitions of communication and computational costs. Our main result is a novel decoupled method that achieves optimal communication cost within the zero-respecting framework. Our method is based on a multi-stage reduction to the decoupled minimization of residual norms, which yields strict improvements over the best known communication cost for the class and the long-standing oracle cost of the Extragradient method. Further, we show by a matching lower bound that our method is communication-optimal within the family of gradient-span algorithms. Finally, we study the extension of distributed SP into Variational Inequality Problem (VIP), which generalizes two-player zero-sum games to multiplayer general-sum games. We show that our decoupled method achieves a new state-of-the-art communication complexity for this broader class.
Tags
Links
- Source: https://arxiv.org/abs/2605.18604v1
- Canonical: https://arxiv.org/abs/2605.18604v1
PDF not stored locally. Use the link above to view on the source site.
Full Text
163,567 characters extracted from source content.
Expand or collapse full text
title=Efficient Gradient Methods for Distributed Saddle Problems, date=May 15, 2026, name=Ruichen Luo, affiliation=Institute of Science and Technology Austria, email=rluo@ist.ac.at name=Anton Rodomanov, affiliation=CISPA Helmholtz Center, email=anton.rodomanov@cispa.de name=Sebastian U. Stich, affiliation=CISPA Helmholtz Center, email=stich@cispa.de The distributed setting for Saddle Problems (SPs) has recently emerged as a framework for various modern applications in machine learning and multiagent systems. Despite its relevance, the theoretical foundations of this setting have not yet been thoroughly established. In this paper, we advance this research direction by formalizing the distributed setup for SPs and providing rigorous definitions of communication and computational costs. Our main result is a novel decoupled method that achieves optimal communication cost within the zero-respecting framework. Our method is based on a multi-stage reduction to the decoupled minimization of residual norms, which yields strict improvements over the best known communication cost for the class and the long-standing oracle cost of the Extragradient method. Further, we show by a matching lower bound that our method is communication-optimal within the family of gradient-span algorithms. Finally, we study the extension of distributed SP into Variational Inequality Problem (VIP), which generalizes two-player zero-sum games to multiplayer general-sum games. We show that our decoupled method achieves a new state-of-the-art communication complexity for this broader class. convex optimization, distributed optimization, complexity bounds, saddle problems, variational inequalities, extra-gradient method, fast gradient method. authors thank Ali Zindari and Krishnendu Chatterjee for their helpful discussions and suggestions on this paper. RL acknowledges the support of ERC CoG 863818 (ForM-SMArt) and Austrian Science Fund (FWF) 10.55776/COE12. Gemini Pro 3.1 was used for part of the writing and analysis. 1 Introduction Motivation. Saddle problems (SPs) and their generalizations, variational inequality problems (VIPs), are of fundamental importance to optimization and game theory. These problems have a wide array of modern applications, including the training of GANs (Goodfellow et al., 2014), robust optimization (Ben-Tal and Nemirovski, 2002), and equilibrium computation in game theory and multiagent systems (von Neumann and Morgenstern, 1947; Rosen, 1965; Hu et al., 1998). The growing scale of modern problems—driven by applications in machine learning, complex game dynamics, and multiagent protocols—renders reliance on a single central processor increasingly impractical. Beyond scalability, a more crucial factor is that many applications are inherently distributed: agents are often geographically dispersed, driven by their own individual interests, and bound by privacy constraints that prohibit the sharing of raw data or utilities. Consequently, distributed computation has become an essential regime for these problems. This perspective underlies a growing body of work in large-scale learning, game-theoretic models, and multiagent systems (McMahan et al., 2017; Zhang et al., 2024; Conitzer and Sandholm, 2004; Nisan and Segal, 2006; Hart and Mansour, 2010; Yoon et al., 2025). In this work, we consider a natural setup where the decision variables and oracles of the SPs or VIPs are partitioned among distributed agents. For instance, in a classic saddle problem minmaxF(,) _x _yF(x,y), we consider one agent controls the minimizing variable x while another controls the maximizing variable y. This partition naturally models, for instance, the strategic autonomy of players in game theory, the interaction protocol in multiagent systems, and the physical separation of the generator and discriminator in GANs (Conitzer and Sandholm, 2004; Goodfellow et al., 2014). Since the decision variables are coupled within their utilities, these agents must coordinate to reach a mutual equilibrium. To do so, they form a communication network that allows them to exchange certain information, such as their current decision variables. Thus, this provides a natural distributed setup where the decision variables are separated among the different agents. While distributed optimization is well-established for finite-sum minimization and federated learning (Schmidt et al., 2017; McMahan et al., 2017), the literature on SPs has primarily focused on extending these data-distributed paradigms (Deng and Mahdavi, 2021; Beznosikov et al., 2025). In contrast, the study of distributed variables and oracles, which are essential to the multi-agent systems, remains a relatively new topic. Although a few recent works have touched upon this direction (Zhang et al., 2024; Zindari et al., 2025; Yoon et al., 2025; Yoon and Loizou, 2025), they predominantly focus on algorithms tailored to specific, favorable scenarios. Consequently, a fundamental gap persists: the lack of a theoretical framework for the general distributed settings of SPs and VIPs. Existing discussions regarding performance often remain at a vague conceptual level, lacking a rigorous formalization of the distributed environment itself. Specifically, there are no standardized definitions for communication and oracle costs in this context. Without such a foundation, it can be difficult to determine the performance limits (lower bounds) or to formally compare the efficiency of different protocols. To enable a rigorous analysis, it is essential to establish metrics that reflect the constraints of distributed multiagent systems, where network latency and bandwidth often dwarf local processing time. In this regime, the primary bottleneck is the communication cost (exchange rounds), while the computational cost (local gradient queries) is a secondary objective. Viewed through this lens, the Extragradient (EG) method (Tseng, 1995; Nemirovski, 2004) serves as the “gold standard” baseline, though the challenges associated with it differ by metric. Regarding computational complexity, consistently improving upon EG for general monotone problems has remained an elusive goal despite over two decades of research. Regarding communication complexity—a metric that has recently come into focus with the rise of distributed systems—EG similarly defines the current state-of-the-art. Surpassing this baseline in the general setting represents a new but critical open problem. This leads to the following research questions: • Formalization and Limits: How can we rigorously formalize the communication and oracle costs for distributed SPs? • Communication Efficiency: Can we design an algorithm that surpasses the state-of-the-art communication bounds for distributed SPs and VIPs? • Oracle Efficiency: Is it possible to consistently improve upon the long-standing oracle complexity of the EG method for general SPs? Contributions. We answer the aforementioned questions in the affirmative, which advances the current theory of distributed SPs and VIPs. • In Section 2, we formalize the distributed saddle-point problem, distributed methods, and their communication and oracle costs. We review EG and other algorithms, casting them as distributed methods to analyze their costs. • In Section 3, we start with a template DM-SP algorithm with a simple, one-loop communication protocol, which improves the state-of-the-art communication cost. • Continuing in Sections 4 and 5, after making the novel multi-stage reduction, we equip the template method with a concrete implementation, thereby consistently improving the (long-standing) oracle cost of EG for general SPs. • In Section 6, by drawing connection to classic convex minimization, we establish the lower bounds for both communication and oracle costs for distributed gradient-span algorithms. In particular, our result shows that our DM-SP algorithm is communication-optimal within the gradient-span algorithm family. • Finally, in Section 7, we extend the results to multi-agent settings. We propose DM-VIP and improve the state-of-the-art communication cost for the class of distributed VIPs. Table 1: Summary of algorithms and complexity results for distributed SPs. Method ℳ M Communication Cost Better Oracle?a Multi-Agentb EGc (θLxyDxDyϵ+LxDx2+LyDy2ϵ)O (θ L_xyD_xD_yε+ L_xD_x^2+L_yD_y^2ε ) No Yes DGDAc (log1ϵ)O ( 1ε ) (weakly coupled only) No Yes Cat-EGc ((LmaxD^xD^yϵ+LxD^x2+LyD^y2ϵ)log2(1ϵ))O ( ( L_ D_x D_yε+ L_x D_x^2+L_y D_y^2ε ) ^2 ( 1ε ) ) Maybe No Cat-Cat-DAGDAc (LxyD^xD^yϵlog3(1ϵ))O ( L_xy D_x D_yε ^3 ( 1ε ) ) Maybe No Lower Bound (Thm. 15) Ω(LxyDxDyϵ) ( L_xyD_xD_yε ) – – DM-SP (Thm. 1) (θLxyDxDyϵ)O (θ L_xyD_xD_yε ) Yes Yes a Indicates whether the method’s theoretical oracle cost outperforms the Extragradient (EG) baseline. b Indicates whether the method supports multi-agent extensions. c These methods are for the non-composite subclass SP∘ P_SP . Notations. Let [n]≜1,…,n[n] \1,…,n\, for any positive integer n. For any finite-dimensional real vector space ℰE, we denote its Euclidean norm by ∥⋅∥ℰ · _E and its dual norm by ∥⋅∥ℰ∗ · _E^*. Specifically, we equip the space ℰx=ℝnxE_x=R^n_x with the norm ∥x=⟨x,⟩1/2 _x= _xx,x ^1/2, where x:ℰx→ℰx∗P_x _x _x^* is a self-adjoint positive definite operator and the dual pairing ⟨ϕx,⟩ _x,x denotes ϕx() _x(x). We denote its corresponding dual norm by ∥⋅∥x∗ · _x^*. We assume analogous geometries for ℰy=ℝnyE_y=R^n_y, ℰi=ℝniE_i=R^n_i (i∈[K]i∈[K]), and ℰw=ℝnwE_w=R^n_w, associated with their respective operators yP_y, iP_i, and wP_w. For a function ψ:ℰ→ℝ∪+∞ψ ∪\+∞\, let domψdomψ denote its effective domain and ∂ψ()∂ψ(z) its subdifferential at ∈domψz ψ. Finally, for any set of vectors S, let spanSspanS denote its linear span. 2 Saddle problems with distributed oracles In the context of saddle problems, we consider two separate finite-dimensional real vector spaces, ℰx=(ℝnx,∥⋅∥x)E_x=(R^n_x, · _x) and ℰy=(ℝny,∥⋅∥y)E_y=(R^n_y, · _y). We are interested in solving composite Saddle Problems (SPs) of the following form: min∈domψxmax∈domψy[F(,)≜f(,)+ψx()−ψy()], _x _x\ _y _y\ [F(x,y) f(x,y)+ _x(x)- _y(y) ], (1) where ψx:ℰx→ℝ∪+∞ _x _x ∪\+∞\ and ψy:ℰy→ℝ∪+∞ _y _y ∪\+∞\ represent relatively simple local components (such as regularizers or indicator functions for constrained sets), and f(⋅,⋅)f(·,·) is a real-valued coupling function defined on an open set containing the domain Q≜domψx×domψyQ _x×dom _y. To simplify the notation, we denote the joint variable by ≜(,)∈Qz (x,y)∈ Q. 2.1 Distributed Methods A method is designed to solve a class of problems sharing a common structure. We begin by introducing the general notion of a problem class, and then formalize what it means for a method to be distributed. Problem class. A problem class, denoted by (ℱ,x,y,Δ,ϵ) P( F, O_x, O_y, ,ε), is a family of problems sharing • a function family ℱ F from which each problem instance is drawn; • two distributed oracles x O_x and y O_y, providing first-order information about an instance; • an accuracy measure Δ , quantifying the quality of a candidate solution; and • a target accuracy ϵ>0ε>0. A specific problem instance P∈P∈ P is determined by a function instance ℱ∈ℱF∈ F, the corresponding oracles (xf,yf)( O_x^f, O_y^f), the accuracy measure Δ , and the target accuracy ϵε. Solving P requires an algorithm to output a candidate solution ¯ z satisfying Δ(¯)≤ϵ ( z)≤ε, accessing only the distributed oracles. The specific problem class of interest in this paper is constructed in Section 2.2. The rest of this subsection formalizes what it means for an algorithm to be a distributed method for solving such a problem class, starting with an engineering description and then a mathematical formalization. Engineering description. We consider a distributed setup with two computational agents, Agent x and Agent y, each maintaining its decision variable (∈domψxx _x and ∈domψyy _y, respectively) in its own local memory. The two agents are fully distributed: neither has direct access to the other’s memory, and they communicate by exchanging messages over a shared channel. Agent x has local oracle access to its component ψx _x (for example, the ability to evaluate its proximal-point mapping) and to a problem-dependent first-order oracle x O_x whose concrete form is specified by the problem class; Agent y has the analogous access to ψy _y and y O_y. The agents proceed in discrete communication rounds. Within each round, each agent performs several local computational steps, each consisting of an oracle query at a chosen point, and the two agents exchange messages at the round boundary. After some number of rounds, the algorithm outputs a candidate solution ¯=(¯,¯) z=( x, y) that approximately solves Problem (1). Mathematical formalization. We adopt the framework of information-based complexity (Nemirovskij and Yudin, 1983). To keep the presentation general, we describe the algorithm in terms of abstract distributed oracles x O_x and y O_y; their concrete instantiation as partial-gradient oracles for SPs is given in Section 2.2. Suppose an algorithm ℳ M proceeds in T communication rounds, where T may be chosen adaptively. In each round t∈0,…,T−1t∈\0,…,T-1\, Agent x successively queries x O_x at τxt _x^t points xt,lz_x^t,l for l=0,…,τxt−1l=0,…, _x^t-1; symmetrically, Agent y queries y O_y at τyt _y^t points yt,lz_y^t,l. The number of local steps and the choice of query points are both decided by the agent based on its accumulated information. We track this via information sets, modeled as ordered sequences. Let Ixt,lI_x^t,l denote the oracle responses collected by Agent x prior to its (l+1)(l+1)-th query in round t. The information set is initialized empty, Ix0,0=∅I_x^0,0= , and each local query appends the corresponding oracle response, so that Ixt,τxtI_x^t, _x^t summarizes the round; Agent y is symmetric. Each agent’s local data, namely the initial point 0z^0 and its component ψx _x (or ψy _y), is treated as known a priori and is not part of the accumulating information set. At the round boundary, the agents exchange messages. From an engineering standpoint, each message is a deterministic function of the sender’s accumulated local information. Mathematically, no restriction is imposed on how a message is formed, and so without loss of generality we let each agent read the union of both agents’ information at the start of the next round: Ixt+1,0=Iyt+1,0=(Ixt,τxt,Iyt,τyt).I_x^t+1,0=I_y^t+1,0= (I_x^t, _x^t,I_y^t, _y^t ). (2) At the conclusion of round t, the algorithm produces a candidate solution ¯t+1=(¯t+1,¯t+1) z^t+1=( x^t+1, y^t+1) from the merged information set Ixt+1,0=Iyt+1,0I_x^t+1,0=I_y^t+1,0. Definition 1. An algorithm ℳ M is called a distributed method if, for every round t∈0,…,T−1t∈\0,…,T-1\, the following hold: 1. Each local query point xt,lz_x^t,l, l∈0,…,τxt−1l∈\0,…, _x^t-1\, is a deterministic function of Agent x’s current information set Ixt,lI_x^t,l; symmetrically, each yt,lz_y^t,l, l∈0,…,τyt−1l∈\0,…, _y^t-1\, is a deterministic function of Iyt,lI_y^t,l. 2. The candidate solution ¯t+1 z^t+1 is a deterministic function of the merged information set Ixt+1,0=Iyt+1,0I_x^t+1,0=I_y^t+1,0 defined in (2). While these deterministic mappings could in principle be randomized, we restrict attention to the deterministic case here for simplicity. Communication and oracle complexities. For a given problem instance P and target accuracy ϵ>0ε>0, the communication complexity of ℳ M on P, denoted TPℳT M_P, is the smallest integer k∈1,…,Tk∈\1,…,T\ such that the candidate solution ¯k z^k satisfies the target accuracy. The total numbers of local oracle queries made by Agent x and Agent y up to that point are Nx,Pℳ=∑t=0TPℳ−1τxtandNy,Pℳ=∑t=0TPℳ−1τyt.N_x,P M= _t=0^T M_P-1 _x^t N_y,P M= _t=0^T M_P-1 _y^t. The communication complexity and oracle complexity of ℳ M over a problem class P are defined by taking the supremum over all instances: Tℳ=supP∈TPℳandNℳ=supP∈(cxNx,Pℳ+cyNy,Pℳ),T M_ P= _P∈ PT M_P N_ P M= _P∈ P (c_xN_x,P M+c_yN_y,P M ), where cxc_x and cyc_y are fixed constants reflecting the computational costs of evaluating a single query to x O_x and y O_y, respectively. Because network communication typically forms the main bottleneck in distributed environments, we treat the communication complexity as the primary performance metric and the oracle complexity as a secondary measure of local computational effort. 2.2 Problem class SP P_SP We now specify the problem class of interest in this paper: composite saddle problems with distributed partial-gradient oracles. The four components of a problem class introduced in Section 2.1 are instantiated in turn below. Function family. We consider the function instances satisfying the following assumptions: (A1) For any fixed ∈domψyy _y, the function f(⋅,)f(·,y) is convex; and for any fixed ∈domψxx _x, the function f(,⋅)f(x,·) is concave. The functions ψx _x and ψy _y are proper, closed, and convex. (A2) Let Dx,Dy>0D_x,D_y>0 be distance parameters, and let 0=(0,0)∈Qz^0=(x^0,y^0)∈ Q be a given initial point. Relative to this initialization, Problem (1) has a saddle point (∗,∗)∈Q(x^*,y^*)∈ Q such that ∥0−∗∥x≤Dxand∥0−∗∥y≤Dy. ^0-x^* _x≤ D_x ^0-y^* _y≤ D_y. (A3) The function f(⋅,⋅)f(·,·) is continuously differentiable over Q. Moreover, its gradients are Lipschitz continuous. That is, with Lipschitz parameters Lx,Lxy,Ly>0L_x,L_xy,L_y>0, for all ,′∈domψxx,x _x and ,′∈domψyy,y _y, we have: ∥∇xf(′,′)−∇xf(,)∥x∗ _xf(x ,y )- _xf(x,y) _x^* ≤Lx∥′−∥x+Lxy∥′−∥y, ≤ L_x -x _x+L_xy -y _y, ∥∇yf(′,′)−∇yf(,)∥y∗ _yf(x ,y )- _yf(x,y) _y^* ≤Lxy∥′−∥x+Ly∥′−∥y. ≤ L_xy -x _x+L_y -y _y. Let ℱ F denote the function family consisting of all instances ℱ=(f,ψx,ψy,0)F=(f, _x, _y,z^0) that satisfy (A1), (A2) and (A3) for a fixed set of parameters (Lx,Lxy,Ly,Dx,Dy)(L_x,L_xy,L_y,D_x,D_y). The non-composite case ψx≡ψy≡0 _x≡ _y≡ 0 is included as a special instance. Distributed partial-gradient oracle. For SPs, the abstract oracles x O_x and y O_y are concretely realized as (deterministic) partial-gradient oracles: for a given function instance ℱ∈ℱF∈ F with coupling function f and any input point ∈Qz∈ Q, • Agent x queries x O_x, which returns xf()=∇xf() O_x^f(z)= _xf(z). • Agent y queries y O_y, which returns yf()=∇yf() O_y^f(z)= _yf(z). The oracles are strictly decoupled: each agent queries only its own oracle, with no access to the counterpart’s. As a concrete example, consider the objective f()=g()+fx()−fy()f(z)=g(z)+f_x(x)-f_y(y), where g()g(z) is a coupled global utility, while fx()f_x(x) and fy()f_y(y) are private utilities accessible only to Agents x and y, respectively. The partial-gradient oracles then take the form xf()=∇xg()+∇fx()yf()=∇yg()−∇fy(), O_x^f(z)= _xg(z)+∇ f_x(x) O_y^f(z)= _yg(z)-∇ f_y(y), for all ∈Qz∈ Q. Due to the distributed setting, Agent x is entirely blind to the private utility fyf_y and can only execute x O_x, and vice versa. Accuracy measure. To evaluate the quality of a candidate solution (¯,¯)∈Q( x, y)∈ Q, we rely on the restricted duality gap. Let ℬx≜∈ℰx∣∥0−∥x≤DxB_x \x _x ^0-x _x≤ D_x\ and ℬy≜∈ℰy∣∥0−∥y≤DyB_y \y _y ^0-y _y≤ D_y\ denote the balls of radii Dx,DyD_x,D_y around the initial points. Over the bounded domain ℬ≜ℬx×ℬyB _x×B_y, the duality gap is defined as Δ(¯,¯)≜max(,)∈ℬ∩Q[F(¯,)−F(,¯)]. ( x, y) _(x,y) ∩ Q [F( x,y)-F(x, y) ]. We say that a pair (¯,¯)∈Q( x, y)∈ Q is an ϵε-saddle point of Problem (1) if Δ(¯,¯)≤ϵ ( x, y)≤ε. Our goal is to design an algorithm that produces such an ϵε-saddle point for a given ϵ>0ε>0. We remark that for the classic problem of constrained optimization with bounded domains, one can enclose the constrained sets in the balls ℬxB_x and ℬyB_y with sufficiently large radius (e.g., the diameter of the constrained sets); the restricted saddle problem in form (1) then coincides with the original one. Problem class SP P_ SP. Assembling the function family ℱ F, the partial-gradient oracles (x,y)( O_x, O_y), the duality-gap accuracy measure Δ , and a target accuracy ϵ>0ε>0, we obtain the problem class of interest, denoted by SP(ℱ,x,y,Δ,ϵ) P_SP( F, O_x, O_y, ,ε), or for short SP P_SP. Solving an instance P∈SPP∈ P_SP requires an algorithm to output an ϵε-saddle point of ℱF utilizing the distributed oracles. To facilitate later discussion, we refer to the terms LxDx2L_xD_x^2 and LyDy2L_yD_y^2 as the diagonal conditioning, and the term LxyDxDyL_xyD_xD_y as the cross-coupled conditioning. 2.3 Existing algorithms from literature In this section, we review existing algorithms for solving SPs and analyze their communication and oracle costs within the distributed method framework. To keep the presentation concise, we summarize the methods and their limitations below, and defer their detailed algorithmic formulations, trajectories, and complexity propositions to Appendix A. Extragradient (EG). The classic EG method (Nemirovski, 2004; Juditsky et al., 2011) naturally fits our framework. Its distributed execution requires two communication rounds per iteration to evaluate coupled partial gradients at both the current and extrapolated points. It provides a robust and natural baseline for communication and oracle costs. Decoupled GDA (DGDA). The DGDA method (Zindari et al., 2025) attempts to reduce communication overhead by freezing the remote variable and taking multiple local gradient steps. While it achieves a fast logarithmic communication cost, its applications are highly restrictive: it only converges for weakly coupled strongly convex-strongly concave instances. For general problem class SP P_SP, the delayed remote variables cause the local updates to drift, leading the method to diverge. Catalyst acceleration. Using a Catalyst wrapper around EG (Cat-EG) (Lin et al., 2020; Yang et al., 2020; Lan and Li, 2026) accelerates the algorithm’s dependence on the diagonal conditioning. However, this comes with five significant caveats: (i) it requires a complicated, multi-loop communication protocol and careful parameter tuning; (i) it is highly sensitive to the inexactness of the diameter estimates D^x D_x and D^y D_y; (i) it introduces multiplicative logarithmic factors in the complexity; (iv) under certain conditioning, its theoretical complexity can be strictly worse than the unaccelerated EG baseline; and (v) it does not support extensions to multi-agent scenarios (cf. Section 7). Four-loop method. The Cat-Cat-DAGDA method (Wang and Li, 2020) applies double Catalyst wrappers around a decoupled accelerated GDA to further accelerate the cross-coupling term. Despite this theoretical improvement, it shares all five caveats of Cat-EG, introduces even more complicated nested loops into the communication protocol, and adds further logarithmic factors. Consequently, it serves primarily as a theoretical benchmark rather than a practical method in our setting. Other distributed stochastic gradient methods. Some recent papers (Zhang et al., 2024; Yoon et al., 2025; Yoon and Loizou, 2025) consider distributed SPs with stochastic gradient oracles. They propose different communication-efficient approaches; however, when applied to standard deterministic oracles considered in this paper, these methods fail to outperform EG. Consequently, as summarized in Table 1, the classic EG method remains a formidable baseline for SP P_SP, and improving its communication and oracle complexity remains a significant challenge. 3 Decoupled method for SPs When designing a communication-efficient method, the primary challenge is enabling distributed agents to compute local solutions independently despite the presence of cross-coupled functions. To address this, we propose a clean algorithmic template (or communication protocol) that reduces an SP into a sequence of coordinate-wise computational tasks. We highlight the key results and insights below, deferring the detailed derivation to Section 4. Assembled norm. Given parameters αx,αy>0 _x, _y>0 (to be specified later), we equip the joint space ℰ=ℰx×ℰyE=E_x×E_y with the assembled norm: ∥ℰ=⟨,⟩12=αx∥x2+αy∥y2for all ∈ℰ, _E= ,z 12= _x _x^2+ _y _y^2 all $ z $, (3) which corresponds to the block diagonal linear operator =αxx⊕αyyP= _xP_x _yP_y. Template DM-SP. Algorithm 1 outlines the Decoupled Method for Saddle Problems (DM-SP), which adapts the abstract framework of the Reduced-Operator Method (Nesterov, 2023) for distributed environments. The algorithm maintains a sequence of anchor points tv^t and proceeds iteratively. First, the agents decouple the joint problem by fixing the remote variable at the current anchor tv^t. This allows Agent x and Agent y to independently and concurrently solve their respective regularized local subproblems up to target accuracies δxt+1 _x^t+1 and δyt+1 _y^t+1 (Lines 4 and 5). Specifically, Agent x aims to approximately compute argmin∈domψx[f(,yt)+αxλt+12∥−xt∥x2+ψx()] *arg\,min_x _x[f(x,v^t_y)+ _x _t+12 -v^t_x _x^2+ _x(x)] by finding a point t+1x^t+1 whose regularized subgradient norm satisfies the exact mathematical bound specified in Line 4. Agent y symmetrically performs an approximate minimization for its corresponding objective −f(xt,)+αyλt+12∥−yt∥y2+ψy()-f(v^t_x,y)+ _y _t+12 -v^t_y _y^2+ _y(y). Following this local computation phase, the agents perform exactly two communication rounds to complete the iteration. In the first round (Line 6), the agents exchange their locally computed approximate solutions to assemble the joint intermediate point t+1=(t+1,t+1)z^t+1=(x^t+1,y^t+1). In the second round (Line 7), they use this assembled point to evaluate their local partial gradients, which they then exchange to form the full joint operator Vψ(t+1)V_ψ(z^t+1). Finally, using this assembled operator, the agents compute a closed-form step size at+1a_t+1, update the running ergodic average ¯t+1 z^t+1, and perform an extragradient-like step to generate the next anchor t+1v^t+1 (Lines 8 and 9). By structuring the method this way, DM-SP cleanly reduces the coupled SP into isolated coordinate-wise tasks with minimal communication overhead. Algorithm 1 DM-SP(f,(ψx,ψy),0,(λt)t≥1,(αx,αy))DM-SP(f,( _x, _y),z^0,( _t)_t≥ 1,( _x, _y)) 1:0=(x0,y0)=0v^0=(v^0_x,v^0_y)=z^0. 2:for t=0,1,…,T−1t=0,1,…,T-1 do 3: Let δxt+1=αxλt+12δ^t+1_x= _x _t+12 and δyt+1=αyλt+12δ^t+1_y= _y _t+12. 4: Agent x finds t+1x^t+1 and ψx′(t+1)∈∂ψx(t+1) _x (x^t+1)∈∂ _x(x^t+1) such that ∥∇xf(t+1,yt)+αxλt+1(t+1−xt)+ψx′(t+1)∥x∗≤δxt+1∥t+1−xt∥x. _xf(x^t+1,v_y^t)+ _x _t+1(x^t+1-v_x^t)+ _x (x^t+1) _x^*≤δ^t+1_x ^t+1-v^t_x _x. 5: Agent y finds t+1y^t+1 and ψy′(t+1)∈∂ψy(t+1) _y (y^t+1)∈∂ _y(y^t+1) such that ∥−∇yf(xt,t+1)+αyλt+1(t+1−yt)+ψy′(t+1)∥y∗≤δyt+1∥t+1−yt∥y. - _yf(v^t_x,y^t+1)+ _y _t+1(y^t+1-v_y^t)+ _y (y^t+1) _y^*≤δ^t+1_y ^t+1-v^t_y _y. 6: Exchange t+1x^t+1 and t+1y^t+1 to assemble t+1=(t+1,t+1)z^t+1=(x^t+1,y^t+1). 7: Calculate corresponding coordinates of Vψ(t+1)V_ψ(z^t+1), then exchange to assemble: Vψ(t+1)=(∇xf(t+1)+ψx′(t+1),−∇yf(t+1)+ψy′(t+1)).V_ψ(z^t+1)= ( _xf(z^t+1)+ _x (x^t+1),- _yf(z^t+1)+ _y (y^t+1) ). 8: Let at+1=2⟨Vψ(t+1),t−t+1⟩∥Vψ(t+1)∥ℰ∗2a_t+1= 2 V_ψ(z^t+1),v^t-z^t+1 V_ψ(z^t+1) _E^*^2 and generate solution ¯t+1=(∑i=1t+1ai)−1∑i=1t+1aii z^t+1= ( _i=1^t+1a_i )^-1 _i=1^t+1a_iz^i. 9: t+1=argmin∈Q[at+1⟨Vψ(t+1),⟩+12∥−t∥ℰ2]v^t+1= *arg\,min_v∈ Q\ [a_t+1 V_ψ(z^t+1),v + 12 -v^t _E^2 ]. 10:end for We refer to Algorithm 1 as a template method because we have not yet specified the implementations for the local computations in Lines 4 and 5. Provided that the inner local solvers in Lines 4 and 5 are standard gradient-based solvers, the template DM-SP procedure formally qualifies as a distributed method. Theorem 1. Consider the DM-SP template applied to SP P_SP, assuming its local trajectories satisfy Definition 1. With the parameter choices of αx=LxyDyDx _x= L_xyD_yD_x, αy=LxyDxDy _y= L_xyD_xD_y, and λt≡λ=2 _t≡λ=2, we have: TSPDM-SP≤2+4LxyDxDyϵ.T^DM-SP_ P_SP≤ 2+4 L_xyD_xD_yε. Remark 1 (Communication improvement). Theorem 1 shows that the communication cost of DM-SP depends only on the cross-coupled conditioning LxyDxDyL_xyD_xD_y, independent of the diagonal conditioning. In contrast, none of the existing methods reviewed in Section 2.3 has achieved this sharp communication guarantee. Specifically, the communication cost of the EG baseline is suboptimal due to its dependence on the diagonal conditioning. While advanced frameworks like Cat-Cat-DAGDA successfully isolate the communication cost from the diagonal conditioning, they suffer from highly complicated nested-loop designs and introduce poly-logarithmic overheads. Therefore, our DM-SP communication protocol represents a clear improvement over existing methods. Furthermore, as shown later in Section 6, our (LxyDxDyϵ)O ( L_xyD_xD_yε ) communication cost is minimax optimal within the family of distributed gradient-span algorithms. Remark 2 (Robustness to inexact distance estimates). Let us consider a practical scenario where the algorithm may not have the precise values of DxD_x and DyD_y in advance, but it has access to upper estimates D^x≥Dx D_x≥ D_x and D^y≥Dy D_y≥ D_y. Let θ≜DxD^yD^xDy+DyD^xD^yDx, θ D_x D_y D_xD_y+ D_y D_x D_yD_x, which quantifies the disproportionality between the true distance parameters and their estimates. Note that θ≥2θ≥ 2 with equality if and only if the estimates are proportional, i.e., D^x/Dx=D^y/Dy D_x/D_x= D_y/D_y. Now, consider the DM-SP template applied to SP P_SP. With the parameter choices of αx=LxyD^yD^x _x= L_xy D_y D_x, αy=LxyD^xD^y _y= L_xy D_x D_y, and λt≡λ=2 _t≡λ=2, we have: TSPDM-SP≤2+2θLxyDxDyϵ.T^DM-SP_ P_SP≤ 2+2θ L_xyD_xD_yε. In particular, θ provides a scale-invariant robustness compared to existing accelerated frameworks. As shown in Table 1, the communication complexities of Cat-EG and Cat-Cat-DAGDA scale directly with the product of the estimates, D^xD^y D_x D_y. Consequently, if both agents conservatively overestimate their domain sizes by a uniform factor c≫1c 1 (i.e., D^x=cDx D_x=cD_x and D^y=cDy D_y=cD_y), the communication cost of Catalyst-based methods inflates by a massive factor of c2c^2. For DM-SP, however, this uniform overestimation perfectly cancels out, yielding θ=2θ=2. 4 Novel multi-stage reduction and its building components In this section, we reveal the technical components of our DM-SP, which is built upon a novel multi-stage reduction. We first leverage the Reduced-Operator Method to reduce the problem to a Monteiro-Svaiter Subproblem (MSS). Then and crucially, we show that when this subproblem is weakly coupled, it can be solved by a Fully Decoupled Solver in one communication round. Consequently, the problem is further reduced to coordinate-wise Minimization of Residual Norms (MRNs). Finally, by exploiting the strong maximal monotonicity, the agents can apply existing accelerated methods to solve the MRNs to desired accuracy. 4.1 Preliminary: Variational inequality problems Let us first introduce the general notion of composite variational inequality problem (VIP) as the backbone of our problems. For any operator V(⋅):domψ→ℰ∗V(·) ψ ^* and any function ψ(⋅):ℰ→ℝ∪+∞ψ(·) ∪\+∞\, we say that ∗∈ℰz^* is a (strong) solution of the VIP of (V,ψ)(V,ψ) if ⟨V(∗),−∗⟩+ψ()≥ψ(∗), for all ∈domψ. V(z^*),z-z^* +ψ(z)≥ψ(z^*), for all z ψ. (4) Assumption for VIPs. Let us introduce the following assumption: (A1’) The function ψ is a (simple) proper closed convex function. The operator V is continuous and monotone over domψdomψ: that is, ⟨V(′)−V(),′−⟩≥0, for all ′,∈domψ. V(z )-V(z),z -z ≥ 0, for all z ,z ψ. There is another notion of a weak solution with the alternative formulation: ⟨V(),−∗⟩+ψ()≥ψ(∗), for all ∈domψ. V(z),z-z^* +ψ(z)≥ψ(z^*), for all z ψ. Under Assumption (A1’), the weak and strong solutions are equivalent. We refer the reader to standard texts such as Nemirovski (2004); Nesterov (2023) for a formal discussion of these two formulations. Moreover, under (A1’), a point ∗∈ℰz^* is the solution of (4) if and only if ∈V(∗)+∂ψ(∗)0∈ V(z^*)+∂ψ(z^*). Indeed, associated with the saddle problem given by (f,ψx,ψy)(f, _x, _y), let us consider the VIP given by (Vf,ψz)(V^f, _z), where Vf()≡(∇xf(),−∇yf()) and ψz()≡ψx()+ψy(), for all ∈Q.V^f(z)≡( _xf(z),- _yf(z)) and _z(z)≡ _x(x)+ _y(y), for all z∈ Q. (5) The associated VIP satisfies (A1’) whenever the saddle problem satisfies (A1). Hence, the solution of the VIP coincides with the saddle point (Nemirovski, 2004; Nesterov, 2023). 4.2 Reduced-operator method for VIPs Now, we introduce the Reduced-Operator Method (ROM) recently proposed in Nesterov (2023). In particular, we apply Nesterov’s general framework in a special way so as to reduce the VIP to a sequence of Monteiro-Svaiter Subproblems (MSSs) (Monteiro and Svaiter, 2013). The MSS asks to find a point for the regularized function such that the subgradient norm at this point is small compared to the distance from the initial point. Let us now define the MSS formally. Monteiro-Svaiter Subproblem. Given an operator V:domψ→ℰ∗V ψ ^*, a function ψ:ℰ→ℝ∪+∞ψ ∪\+∞\, a reference point ∈domψv ψ, and a real number λ>0λ>0, we say (+,ψ′(+))(z^+,ψ (z^+)) is a solution of the MSS if +∈domψz^+ ψ, ψ′(+)∈∂ψ(+)ψ (z^+)∈∂ψ(z^+), and ∥V(+)+ψ′(+)+λ(+−)∥ℰ∗≤λ∥+−∥ℰ. V(z^+)+ψ (z^+)+ (z^+-v) _E^*≤λ ^+-v _E. (6) We will discuss how to solve the MSSs later in Section 4.3. But for now, let us assume there exists a solver ℳMS(V,ψ,,λ)M^MS(V,ψ,v,λ) for the MSSs, which takes an MSS given by (V,ψ,,λ)(V,ψ,v,λ) and returns a solution of it. Built upon such a solver ℳMSM^MS, we now introduce ROM in Algorithm 2. At each iteration t: the solver ℳMSM^MS returns a solution (t+1,ψ′(t+1))(z^t+1,ψ (z^t+1)) for the MSS built at reference point tv^t; this solution is used as a midpoint to compute subgradient Vψ(t+1)V_ψ(z^t+1); then the ‘extragradient-type’ step is taken with the stepsize at+1a_t+1. Algorithm 2 ROM∥⋅∥ℰ(V,ψ,0,(λt)t≥1∣ℳMS)ROM_ · _E(V,ψ,z^0,( _t)_t≥ 1 ^MS) 1:A solver ℳMSM^MS for the MSSs. 2:0=0v^0=z^0. 3:for t=0,1,⋯t=0,1,·s do 4: (t+1,ψ′(t+1))=ℳMS(V,ψ,t,λt+1)(z^t+1,ψ (z^t+1))=M^MS(V,ψ,v^t, _t+1). 5: Vψ(t+1)=V(t+1)+ψ′(t+1)V_ψ(z^t+1)=V(z^t+1)+ψ (z^t+1). 6: at+1=2⟨Vψ(t+1),t−t+1⟩∥Vψ(t+1)∥ℰ∗2a_t+1= 2 V_ψ(z^t+1),v^t-z^t+1 V_ψ(z^t+1) _E^*^2. 7: t+1=argmin∈domψ[at+1⟨Vψ(t+1),⟩+12∥−t∥ℰ2]v^t+1= *arg\,min_v ψ [a_t+1 V_ψ(z^t+1),v + 12 -v^t _E^2 ]. 8:end for Next, let us show the convergence of ROM. Lemma 2. ROM (Algorithm 2) ensures for all ∈domψv ψ and for all T≥1T≥ 1, ∑t=0T−1at+1⟨Vψ(t+1),t+1−⟩≤12∥0−∥ℰ2−12∥T−∥ℰ2. _t=0^T-1a_t+1 V_ψ(z^t+1),z^t+1-v ≤ 12 ^0-v _E^2- 12 ^T-v _E^2. Moreover, we have at+1≥1λt+1a_t+1≥ 1 _t+1, for all t≥0t≥ 0. Proof. By the optimality of t+1v^t+1, we have for all ∈domψv ψ, at+1⟨Vψ(t+1),−t+1⟩+12∥t−∥ℰ2≥12∥t+1−∥ℰ2+12∥t−t+1∥ℰ2,a_t+1 V_ψ(z^t+1),v-v^t+1 + 12 ^t-v _E^2≥ 12 v^t+1-v _E^2+ 12 v^t-v^t+1 _E^2, and therefore, at+1⟨Vψ(t+1),−t+1⟩+12∥t−∥ℰ2 a_t+1 V_ψ(z^t+1),v-z^t+1 + 12 ^t-v _E^2 ≥at+1⟨Vψ(t+1),t+1−t+1⟩+12∥t+1−∥ℰ2+12∥t−t+1∥ℰ2 ≥ a_t+1 V_ψ(z^t+1),v^t+1-z^t+1 + 12 ^t+1-v _E^2+ 12 ^t-v^t+1 _E^2 =at+1⟨Vψ(t+1),t−t+1⟩+12∥t+1−∥ℰ2+at+1⟨Vψ(t+1),t+1−t⟩+12∥t−t+1∥ℰ2 =a_t+1 V_ψ(z^t+1),v^t-z^t+1 + 12 ^t+1-v _E^2+a_t+1 V_ψ(z^t+1),v^t+1-v^t + 12 ^t-v^t+1 _E^2 ≥at+1⟨Vψ(t+1),t−t+1⟩+12∥t+1−∥ℰ2−at+122∥Vψ(t+1)∥ℰ∗2 ≥ a_t+1 V_ψ(z^t+1),v^t-z^t+1 + 12 ^t+1-v _E^2- a_t+1^22 V_ψ(z^t+1) _E^*^2 ≥12∥t+1−∥ℰ2, ≥ 12 ^t+1-v _E^2, where the last inequality follows from the definition of at+1a_t+1 in Line 6 of Algorithm 2. Then, the desired bound follows from summing the above inequality over t from 0 to T−1T-1. Next, we show the lower bound for ata_t. For all t≥1t≥ 1, we have ⟨Vψ(t),t−1−t⟩−12λt∥Vψ(t)∥ℰ∗2 V_ψ(z^t),v^t-1-z^t - 12 _t V_ψ(z^t) _E^*^2 ≡λt2∥t−t−1∥ℰ2−12λt∥Vψ(t)+λt(t−t−1)∥ℰ∗2 ≡ _t2 ^t-v^t-1 _E^2- 12 _t V_ψ(z^t)+ _tP(z^t-v^t-1) _E^*^2 ≥0, ≥ 0, where the last inequality follows from Eq. 6. Therefore, we have at=2⟨Vψ(t),t−t⟩∥Vψ(t)∥ℰ∗2≥1λt.a_t= 2 V_ψ(z^t),v^t-z^t V_ψ(z^t) _E^*^2≥ 1 _t. ∎ 4.3 Fully decoupled solver for MSSs with weak couplings We now address the MSS introduced by the ROM in Section 4.2. Specifically, we are to deal with the MSS given by (Vf,ψz,,λ),(V^f, _z,v,λ), (7) where VfV^f and ψz _z are defined in Eq. 5, the reference point =(x,y)∈domψx×domψyv=(v_x,v_y) _x×dom _y, and the assembled norm ∥⋅∥ℰ · _E is associated with parameters αx _x and αy _y. We say that an MSS has a weak coupling if λ≥2L¯c≜2Lxyαxαy.λ≥ 2 L_ c 2L_xy _x _y. In this section, we will introduce a Fully Decoupled Solver (FDS), which reduces the weakly-coupled MSSs to coordinate-wise Minimization of Residual Norms (MRNs). Let us first define the problem of MRN. The MRN asks to find an approximate solution +w^+ of the VIP of (Vw,ψw)(V_w, _w) such that the residual norm is small compared to the distance from the initial point: Minimization of residual norm. Given an operator Vw:domψw→ℰw∗V_w _w _w^*, a function ψw:ℰw→ℝ∪+∞ _w _w ∪\+∞\, a reference point w∈domψwv_w _w, and an accuracy δ>0δ>0, we say (+,ψw′(+))(w^+, _w (w^+)) minimizes the residual norm to δ-relative distance accuracy, if +∈domψww^+ _w, ψw′(+)∈∂ψw(+) _w (w^+)∈∂ _w(w^+), and ∥Vw(+)+ψw′(+)∥w∗≤δ∥+−w∥w. V_w(w^+)+ _w (w^+) _w^*≤δ ^+-v_w _w. Let us, again, defer the discussion of solving MRNs to Section 4.4. But for now, let us assume there exist solvers ℳxMRN(Vx,ψ^x,x,δx)andℳyMRN(Vy,ψ^y,y,δy)M_x^MRN(V_x, ψ_x,v_x, _x) _y^MRN(V_y, ψ_y,v_y, _y) for the coordinate-wise MRNs in spaces ℰxE_x and ℰyE_y. In particular, these solvers take a coordinate-wise MRN problem and return a point and a subgradient satisfying the desired accuracy. A crucial step in our analysis relies on the Fully Decoupled Solver (FDS, Algorithm 3), which optimizes each decision variable independently. For an MSS with weak coupling, we establish in Lemma 3 that FDS returns a correct solution in a single round. Algorithm 3 FDS∥⋅∥ℰ((∇xf,−∇yf),(ψx,ψy),,λ∣(ℳxMRN,ℳyMRN))FDS_ · _E(( _xf,- _yf),( _x, _y),v,λ (M_x^MRN,M_y^MRN)) 1:Solvers ℳxMRNM_x^MRN and ℳyMRNM_y^MRN for the coordinate-wise MRNs. 2:ψ^x=ψx+αxλ2∥⋅−x∥x2 ψ_x= _x+ _xλ2 ·-v_x _x^2 and ψ^y=ψy+αyλ2∥⋅−y∥y2 ψ_y= _y+ _yλ2 ·-v_y _y^2. 3:Agent x and Agent y respectively compute (+,ψ^x′(+))=ℳxMRN(∇xf(⋅,y),ψ^x,x,δx) and (+,ψ^y′(+))=ℳyMRN(−∇yf(x,⋅),ψ^y,y,δy), gathered(x^+, ψ _x(x^+))=M_x^MRN ( _xf(·,v_y), ψ_x,v_x, _x ) and \\ (y^+, ψ _y(y^+))=M_y^MRN (- _yf(v_x,·), ψ_y,v_y, _y ), gathered where δx=αxλ2 _x= _xλ2 and δy=αyλ2 _y= _yλ2. 4:ψx′(+)=ψ^x′(+)−αxλx(+−x)ψ _x(x^+)= ψ _x(x^+)- _x _x(x^+-v_x) and ψy′(+)=ψ^y′(+)−αyλy(+−y)ψ _y(y^+)= ψ _y(y^+)- _y _y(y^+-v_y). 5:return (+,ψ′(+))(z^+,ψ (z^+)), where +=(+,+)z^+=(x^+,y^+) and ψ′(+)=(ψx′(+),ψy′(+))ψ (z^+)=(ψ _x(x^+),ψ _y(y^+)). Lemma 3. Consider the saddle problem given by (f,ψx,ψy)(f, _x, _y) which satisfies (A3). For λ≥2L¯cλ≥ 2 L_ c, FDS (Algorithm 3) returns a solution of the MSS given in Eq. 7. The correctness of the FDS for SPs can be implied as a direct consequence of the correctness of a more general version of FDS for VIPs, which will be introduced later in Appendix C. Therefore, let us defer this proof to Lemma 22 in Appendix C. 4.4 Minimization of residual norms We arrive at the last building component, the Minimization of Residual Norms (MRNs). Assumptions for MRNs. Let us introduce the following assumptions: (Â1) The function ψw _w is a (simple) proper closed convex function. The operator VwV_w is monotone over domψwdom _w. (Â2) The set-valued operator Vw+∂ψwV_w+∂ _w is μ-strongly maximally monotone over domψwdom _w. That is, ⟨Vw(′)+ψw′(′)−Vw()−ψw′(),′−⟩≥μ∥′−∥w2, V_w(w )+ _w (w )-V_w(w)- _w (w),w -w ≥μ -w _w^2, for all ′,∈domψww ,w _w, ψw′(′)∈∂ψw(′)ψ _w(w )∈∂ _w(w ), and ψw′()∈∂ψw()ψ _w(w)∈∂ _w(w), (Â3) The operator Vw()V_w(w) is L-Lipschitz continuous over ∈domψww _w. (Â4) The operator Vw=∇fwV_w=∇ f_w, where fwf_w is a continuously differentiable function defined on an open set containing domψwdom _w. The theoretical guarantee provided in the literature is usually based on the distance-to-solution accuracy (cf. Definition 2). We show in Lemma 4 that, under strong maximal monotonicity, the relative distance accuracy required in this paper can be implied from distance-to-solution accuracy. Definition 2 (Distance-to-solution accuracy). We say that (+,ψw′(+))(w^+, _w (w^+)) satisfies ξ-distance-to-solution accuracy if +∈domψww^+ _w, ψw′(+)∈∂ψw(+) _w (w^+)∈∂ _w(w^+), and ∥Vw(+)+ψw′(+)∥w∗≤ξ∥w−~∥w V_w(w^+)+ _w (w^+) _w^*≤ξ v_w- w _w for some ~ w in the solution set of the VIP of (Vw,ψw)(V_w, _w). Lemma 4. Consider the MRN problem given by (Vw,ψw,w,δ)(V_w, _w,v_w,δ) which satisfies (Â2). Let ξ≤μδμ+δξ≤ μδμ+δ. If (+,ψw′(+))(w^+, _w (w^+)) satisfies ξ-distance-to-solution accuracy, then (+,ψw′(+))(w^+, _w (w^+)) is a solution of the MRN problem. Proof. In view of the triangle inequality and the μ-strong maximal monotonicity of Vw+ψwV_w+ _w, we have ∥w−~∥w≤∥+−w∥w+∥+−~∥w _w- w _w≤ ^+-v_w _w+ ^+- w _w ≤∥+−w∥w+1μ∥Vw(+)+ψw′(+)∥w∗ ≤ ^+-v_w _w+ 1μ V_w(w^+)+ _w (w^+) _w^* ≤∥+−w∥w+ξμ∥w−~∥w. ≤ ^+-v_w _w+ ξμ _w- w _w. Then, we have ∥w−~∥w≤μ−ξ∥+−w∥w. _w- w _w≤ μ-ξ ^+-v_w _w. Therefore, we have ∥Vw(+)+ψw′(+)∥w∗≤ξ∥w−~∥w≤μξμ−ξ∥+−w∥w≤δ, V_w(w^+)+ _w (w^+) _w^*≤ξ _w- w _w≤ μξμ-ξ ^+-v_w _w≤δ, where the last inequality follows from the assignment ξ≤μδμ+δξ≤ μδμ+δ. ∎ We will leverage efficient existing solvers for the MRN problems. In particular, we are to deal with the specific MRNs given in Line 3 in Algorithm 3, where the corresponding coordinate-wise operators are gradients of smooth convex functions. Therefore, we can leverage the existing accelerated gradient methods from the literature. Let us apply, for instance, the Accumulative Regularization Method (ARM) from Lan et al. (2023), whose detailed pseudocode is presented in Algorithm 4 in Appendix B for completeness. Let us denote this algorithm as ARM(∇fw,ψw,w,ξ∣L),ARM\ (∇ f_w, _w,v_w,ξ L), (8) which takes an MRN instance of interest, has knowledge of the parameter L in (Â3), and returns a solution that satisfies ξ-distance-to-solution accuracy. Now, we state the oracle complexity of MRN with respect to the distance-to-solution accuracy. The original result of Lan et al. (2023) is given in projected gradient norm, which can be converted to the subgradient norm considered in this paper. We defer the detailed proof to Appendix B. Lemma 5 (Lan et al. 2023). Assume (Â1), (Â3), (Â4), and that the solution set of the VIP of (Vw,ψw)(V_w, _w) is non-empty. Let (+,ψw′(+))=ARM(∇fw,ψw,w,ξ∣L).(w^+, _w (w^+))=ARM\ (∇ f_w, _w,v_w,ξ L). Then, ARM takes no more than 34⋅3L2ξ34· 3L2ξ queries to ∇fw(⋅)∇ f_w(·) and ensures that (+,ψw′(+))(w^+, _w (w^+)) satisfies ξ-distance-to-solution accuracy. 5 Decoupled method for SPs: Concrete implementation We are now back to considering the original SPs in Eq. 1. Let us combine the technical components in Section 4 and present the final, implementable version of DM-SP. Implementation of DM-SP. We use the ARM solver in Eq. 8 for Minimization of Residual Norms: ℳxMRN(Vx,ψ^x,x,δx)≜ARM(Vx,ψ^x,x,2δx3∣Lx),ℳyMRN(Vy,ψ^y,y,δy)≜ARM(Vy,ψ^y,y,2δy3∣Ly). gatheredM^MRN_x(V_x, ψ_x,v_x, _x) (V_x, ψ_x,v_x, 2 _x3 L_x ),\\ M^MRN_y(V_y, ψ_y,v_y, _y) (V_y, ψ_y,v_y, 2 _y3 L_y ). gathered (9) Consider the assembled norm ∥⋅∥ℰ · _E with parameters αx _x and αy _y. Then, for any Monteiro-Svaiter Subproblem given by (Vf,ψz,,λ)(V^f, _z,v,λ), we leverage the solver FDS-ARM(Vf,ψz,,λ)=FDS∥⋅∥ℰ(Vf,ψz,,λ∣(ℳxMRN,ℳyMRN))FDS-ARM(V^f, _z,v,λ)=FDS_ · _E (V^f, _z,v,λ (M^MRN_x,M^MRN_y) ) Finally, we obtain the concrete algorithm DM-SP as follows: ROM∥⋅∥ℰ(Vf,ψz,0,(λt)t≥1∣FDS-ARM).ROM_ · _E (V^f, _z,z^0,( _t)_t≥ 1 -ARM ). (10) Combining Lemmas 2, 3, 4 and 5, we are ready to prove the main convergence lemma for SPs in Lemma 6. Lemma 6. Consider DM-SP with the implementation in Eq. 10, applied to SP P_SP. Under (A2), (A1) and (A3), for λt+1≡λ≥2Lxyαxαy _t+1≡λ≥ 2L_xy _x _y, the algorithm takes no more than 2T2T communication rounds, no more than T⋅(1+349Lx2αxλ)T· (1+34 9L_x2 _xλ ) queries to ∇xf _xf, and no more than T⋅(1+349Ly2αyλ)T· (1+34 9L_y2 _yλ ) queries to −∇yf- _yf, and obtains an ϵε-saddle point ¯T, z^T, where T=⌈αxλDx2+αyλDy22ϵ⌉.T= _xλ D_x^2+ _yλ D_y^22ε . Proof. By (A1), we have Δ(¯T)≤(∑t=0T−1at+1)−1max∈ℬ∩Q[∑t=0T−1at+1⟨Vψ(t+1),t+1−⟩]. ( z^T)≤ ( _t=0^T-1a_t+1 )^-1 _z ∩ Q\ [ _t=0^T-1a_t+1 V_ψ(z^t+1),z^t+1-z ]. Further, with λ≥2L¯cλ≥ 2 L_ c, by Lemmas 2 and 3, we have Δ(¯T)≤(∑t=0T−1at+1)−1max∈ℬ∩Q[∑t=0T−1at+1⟨Vψ(t+1),t+1−⟩] ( z^T)≤ ( _t=0^T-1a_t+1 )^-1 _z ∩ Q\ [ _t=0^T-1a_t+1 V_ψ(z^t+1),z^t+1-z ] ≤(∑t=0T−1at+1)−1max∈ℬ∩Q[αx2∥0−∥x2+αy2∥0−∥y2] ≤ ( _t=0^T-1a_t+1 )^-1 _z ∩ Q\ [ _x2 ^0-x _x^2+ _y2 ^0-y _y^2 ] ≤(∑t=0T−11λt+1)−1⋅12(αxDx2+αyDy2)≤ϵ, ≤ ( _t=0^T-1 1 _t+1 )^-1· 12( _xD_x^2+ _yD_y^2)≤ε, where the last inequality follows from the assignments of (λt)t≥1( _t)_t≥ 1 and T. Therefore, the number of communication rounds is bounded by 2T2T. Now we count the number of gradient queries. By Lemma 4, ARM always returns the solution with the required relative distance accuracy; and in view of Lemma 5, it takes no more than 349Lx2αxλ34 9L_x2 _xλ gradient queries to ∇xf _xf and no more than 349Ly2αyλ34 9L_y2 _yλ gradient queries to −∇yf- _yf. Therefore, the numbers of gradient queries to ∇xf _xf and −∇yf- _yf are bounded by T⋅(1+349Lx2αxλ)T· (1+34 9L_x2 _xλ ) and T⋅(1+349Ly2αyλ)T· (1+34 9L_y2 _yλ ), respectively. ∎ Finally, we conclude with the following guarantee on the oracle cost. Theorem 7. Consider DM-SP with the implementation in Eq. 10, applied to SP P_SP. With the same choices of αx _x, αy _y, and (λt)t≥0( _t)_t≥ 0 as in Theorem 1, we have: NSPDM-SP=(cx+cy)2LxyDxDyϵ+102(LxyDxDyϵ)12(cx(LxDx2ϵ)12+cy(LyDy2ϵ)12).N_ P_SP^DM-SP=(c_x+c_y) 2L_xyD_xD_yε+102 ( L_xyD_xD_yε ) 12 (c_x ( L_xD_x^2ε ) 12+c_y ( L_yD_y^2ε ) 12 ). As detailed in Section 2.3, alternative distributed algorithms either fail to outperform EG in general for SP P_SP or suffer from certain theoretical and practical caveats. Consequently, the EG method remains a crucial baseline for oracle complexity, which we now compare against. Remark 3 (Oracle improvement). Let us compare the oracle cost of DM-SP against the EG baseline, which requires NSP∘EG=(cx+cy)⋅(LxyDxDyϵ+LxDx2ϵ+LyDy2ϵ)N^EG_ P _SP=(c_x+c_y)· ( L_xyD_xD_yε+ L_xD_x^2ε+ L_yD_y^2ε ) oracle costs (Proposition 17). We conclude that the computational cost of DM-SP is consistently bounded by that of EG. Furthermore, it yields a substantial improvement when LxDx2+LyDy2≫LxyDxDy+LxyDxDy⋅(cxcx+cyLxDx2+cycx+cyLyDy2).L_xD_x^2+L_yD_y^2 L_xyD_xD_y+ L_xyD_xD_y· ( c_xc_x+c_y L_xD_x^2+ c_yc_x+c_y L_yD_y^2 ). For instance, assuming uniform oracle costs (cx=cyc_x=c_y), this improvement occurs when the diagonal conditioning dominates the cross-coupled conditioning, i.e., LxDx2+LyDy2≫LxyDxDyL_xD_x^2+L_yD_y^2 L_xyD_xD_y. To our knowledge, DM-SP is the first method to consistently improve upon the EG oracle cost for SP P_SP. 6 Lower complexity bounds for SPs In this section, we establish lower complexity bounds for distributed SPs. In particular, we focus on algorithms in the distributed gradient-span framework (Nesterov, 2004). This restriction may not be absolutely necessary, and we might be able to avoid it via more complicated reasoning. However, it naturally holds for the majority of practical algorithms as detailed in Appendix A. Hence, we consider this as a crucial first step towards establishing information-theoretic lower bounds for distributed SPs. To establish the lower bounds, we analyze a subclass of problems, denoted by SP∘⊂SP P_SP ⊂ P_SP, where the local regularizing components are identically zero (i.e., ψx≡0 _x≡ 0 and ψy≡0 _y≡ 0). Establishing lower bounds on this unconstrained, smooth subclass yields lower bounds for the broader class SP P_SP. Because the local functions are zero, the agents’ updates rely exclusively on the partial gradients of the coupling function f, and no subgradients are involved. For the lower bound analysis, we consider the case where the messages exchanged between rounds allow each agent to reconstruct the other’s historical gradients (which can be seen as the worst case for lower bound analysis). Thus, Agent x utilizes its own gradients up to the current local step, and Agent y’s gradients up to the end of the previous round. 6.1 Distributed gradient-span framework Building upon the information-based complexity established in Section 2.1, the gradient-span assumption imposes a structural restriction on the deterministic mappings. We assume that the iterates are constructed by taking linear combinations of the historically observed preconditioned partial gradients. Algorithm trajectories and histories. To describe the available gradients at any given step, we use a more specific notation for the local query trajectories. In each round t∈0,…,T−1t∈\0,…,T-1\, Agent x and Agent y generate local query trajectories of lengths τxt _x^t and τyt _y^t. The query points consist of the agent’s local variable and its delayed approximation of the remote variable: Z^xt=xt,l=(t,l,^t,l)l=0τxt−1andZ^yt=yt,l=(^t,l,t,l)l=0τyt−1. Z_x^t= \z_x^t,l=(x^t,l, y^t,l) \_l=0 _x^t-1 Z_y^t= \z_y^t,l=( x^t,l,y^t,l) \_l=0 _y^t-1. Let Zxt−1=⋃i=0t−1Z^xiZ_x^t-1= _i=0^t-1 Z_x^i and Zyt−1=⋃i=0t−1Z^yiZ_y^t-1= _i=0^t-1 Z_y^i denote the accumulated sets of query points from all prior rounds (with Zx−1=Zy−1=∅Z_x^-1=Z_y^-1= ). Within round t, the sets of queries made up to local step l are denoted by Zxt,l=Zxt−1∪xt,ii=0l−1Z_x^t,l=Z_x^t-1∪\z_x^t,i\_i=0^l-1 and Zyt,l=Zyt−1∪yt,ii=0l−1Z_y^t,l=Z_y^t-1∪\z_y^t,i\_i=0^l-1. Finally, the accumulated queries up to the end of round t are Zxt=Zxt,τxtZ_x^t=Z_x^t, _x^t and Zyt=Zyt,τytZ_y^t=Z_y^t, _y^t. Assumption 1 (Local variables). For any instance P∈SP∘P∈ P_SP and all rounds t∈0,…,T−1t∈\0,…,T-1\, the local variables queried by Agent x and Agent y satisfy, for all respective local steps l: t,l ^t,l ∈0+x−1span∇xf()∣∈Zxt,l, ^0+P_x^-1span \ _xf(z) ∈ Z_x^t,l \, t,l ^t,l ∈0+y−1span∇yf()∣∈Zyt,l. ^0+P_y^-1span \ _yf(z) ∈ Z_y^t,l \. Assumption 2 (Remote variables). For any instance P∈SP∘P∈ P_SP and all rounds t∈0,…,T−1t∈\0,…,T-1\, the delayed approximations of the remote variables queried by Agent x and Agent y satisfy, for all respective local steps l: ^t,l y^t,l ∈0+y−1span∇yf()∣∈Zyt−1, ^0+P_y^-1span \ _yf(z) ∈ Z_y^t-1 \, ^t,l x^t,l ∈0+x−1span∇xf()∣∈Zxt−1. ^0+P_x^-1span \ _xf(z) ∈ Z_x^t-1 \. Assumption 3 (Candidate solutions). For any instance P∈SP∘P∈ P_SP and all rounds t∈0,…,T−1t∈\0,…,T-1\, the candidate solutions ¯t+1=(¯t+1,¯t+1) z^t+1=( x^t+1, y^t+1) generated at the end of round t satisfy ¯t+1 x^t+1 ∈0+x−1span∇xf()∣∈Zxt, ^0+P_x^-1span \ _xf(z) ∈ Z_x^t \, ¯t+1 y^t+1 ∈0+y−1span∇yf()∣∈Zyt. ^0+P_y^-1span \ _yf(z) ∈ Z_y^t \. Definition 3. An algorithm is called gradient-span if Assumptions 1, 2 and 3 are all satisfied. These assumptions specify that updates are confined to the available gradient subspaces. The impact of delayed communication is modeled by restricting the cross-variable approximations ( y and x) to the remote gradients from Zyt−1Z_y^t-1 and Zxt−1Z_x^t-1, respectively. This gradient-span family includes the majority of practical algorithms, such as EG, Decoupled-GDA, Cat-EG, and Cat-Cat-DAGDA. A more detailed discussion is provided in Appendix A. Finally, we remark that our proposed algorithm, DM-SP, also operates within this framework. Specifically, each iteration t of DM-SP requires exactly two communication rounds. In the first round, the local solvers perform multiple steps to compute the approximate solutions t+1x^t+1 and t+1y^t+1, generating local trajectories denoted by Z^x2t Z_x^2t and Z^y2t Z_y^2t. In the second round, the agents evaluate partial gradients at the assembled point t+1z^t+1, yielding singleton trajectories Z^x2t+1=Z^y2t+1=t+1 Z_x^2t+1= Z_y^2t+1=\z^t+1\. 6.2 Lower complexity bounds Now we provide detailed proofs for the lower complexity bounds of distributed gradient-span algorithms applied to SP P_SP. We focus on a subclass of unconstrained SPs with ψx=ψy=0 _x= _y=0, denoted SP∘ P_SP . We assume the initial points are 0=x^0=0 and 0=y^0=0. Lower bounds established for this subclass hold for the general class SP P_SP. Let us consider the spaces ℰx=ℝnxE_x=R^n_x and ℰy=ℝnyE_y=R^n_y equipped with the standard Euclidean inner product ⟨⋅,⋅⟩ ·,· and the corresponding Euclidean norm ∥⋅∥ · . The cases with general preconditioned norms ∥⋅∥x · _x and ∥⋅∥y · _y can be proven similarly. Function subfamilies. To establish the overall lower bound, we decompose the general function family into three distinct subfamilies: Definition 4. Let ℱx=ℱ(Lx,0,0,Dx,0) F_x= F(L_x,0,0,D_x,0) denote the set of functions Fx(,)=12∥x−x∥2,F_x(x,y)= 12 _xx-b_x ^2, where the matrix x∈ℝnx×nxA_x ^n_x× n_x satisfies ∥x∥2≤Lx _x ^2≤ L_x, and the vector x∈ℝnxb_x ^n_x is such that the linear system x=xA_xx=b_x has a solution ∗∈ℰxx^* _x satisfying ∥∗∥≤Dx ^* ≤ D_x. Definition 5. Let ℱy=ℱ(0,0,Ly,0,Dy) F_y= F(0,0,L_y,0,D_y) denote the set of functions Fy(,)=−12∥y−y∥2,F_y(x,y)=- 12 _yy-b_y ^2, where the matrix y∈ℝny×nyA_y ^n_y× n_y satisfies ∥y∥2≤Ly _y ^2≤ L_y, and the vector y∈ℝnyb_y ^n_y is such that the linear system y=yA_yy=b_y has a solution ∗∈ℰyy^* _y satisfying ∥∗∥≤Dy ^* ≤ D_y. Definition 6. Let ℱxy=ℱ(0,Lxy,0,Dx,Dy) F_xy= F(0,L_xy,0,D_x,D_y) denote the set of functions Fxy(,)=⟨xy−xy,⟩,F_xy(x,y)= _xyx-b_xy,y , where the matrix xy∈ℝny×nxA_xy ^n_y× n_x satisfies ∥xy∥≤Lxy _xy ≤ L_xy, and the vector xy∈ℝnyb_xy ^n_y is such that the linear system xy=xyA_xyx=b_xy has a solution ∗∈ℰxx^* _x satisfying ∥∗∥≤Dx ^* ≤ D_x. Distributed oracles. We consider the pairs of distributed oracles (x,y)( O_x, O_y) that return the partial gradients of the functions in these subfamilies: • For Fx∈ℱxF_x∈ F_x: x(,)=∇xFx(,)=x⊤(x−x) O_x(x,y)= _xF_x(x,y)=A_x (A_xx-b_x) and y(,)=−∇yFx(,)= O_y(x,y)=- _yF_x(x,y)=0. • For Fy∈ℱyF_y∈ F_y: x(,)=∇xFy(,)= O_x(x,y)= _xF_y(x,y)=0 and y(,)=−∇yFy(,)=y⊤(y−y) O_y(x,y)=- _yF_y(x,y)=A_y (A_yy-b_y). • For Fxy∈ℱxyF_xy∈ F_xy: x(,)=∇xFxy(,)=xy⊤ O_x(x,y)= _xF_xy(x,y)=A_xy y and y(,)=−∇yFxy(,)=xy−xy O_y(x,y)=- _yF_xy(x,y)=b_xy-A_xyx. Accuracy measures. For any approximate solution ¯=(¯,¯)∈ℰx×ℰy z=( x, y) _x×E_y, the restricted primal-dual gaps over the bounded sets evaluate to: • For Fx∈ℱxF_x∈ F_x: ΔFx(¯)=Fx(¯,¯)−min∈ℬxFx(,¯)=12∥x¯−x∥2 _F_x( z)=F_x( x, y)- _x _xF_x(x, y)= 12 _x x-b_x ^2. • For Fy∈ℱyF_y∈ F_y: ΔFy(¯)=max∈ℬyFy(¯,)−Fy(¯,¯)=12∥y¯−y∥2 _F_y( z)= _y _yF_y( x,y)-F_y( x, y)= 12 _y y-b_y ^2. • For Fxy∈ℱxyF_xy∈ F_xy: ΔFxy(¯)=max∈ℬy⟨xy¯−xy,⟩−min∈ℬx⟨xy−xy,¯⟩≥Dy∥xy¯−xy∥ _F_xy( z)= _y _y _xy x-b_xy,y - _x _x _xyx-b_xy, y ≥ D_y _xy x-b_xy , where the inequality holds since min∈ℬx⟨xy−xy,¯⟩≤⟨xy∗−xy,¯⟩=0 _x _x _xyx-b_xy, y ≤ _xyx^*-b_xy, y =0. Problem subclasses. These constructions yield three distinct problem subclasses: • x=(ℱx,x,y,ϵ) P_x= ( F_x, O_x, O_y,ε ), • y=(ℱy,x,y,ϵ) P_y= ( F_y, O_x, O_y,ε ), and • xy=(ℱxy,x,y,ϵ) P_xy= ( F_xy, O_x, O_y,ε ). The worst-case complexity for general SPs is bounded below by the maximum complexity among these three subclasses. For the rest of this section, we present our lower bound proofs by establishing a connection from distributed SPs to convex minimization. In Section 6.2.1, we recall and properly rescale Nesterov’s construction for unconstrained convex optimization. We analyze the three problem subclasses separately in Sections 6.2.2 and 6.2.3, mapping the lower bound for each subclass to this core convex minimization problem. Finally, we combine the results to establish the overall lower bound for distributed SPs in Section 6.2.4. 6.2.1 The “worst function in the world” (with proper rescaling) While the original formulation (Theorem 2.1.7 by Nesterov (2004)) establishes a lower bound with respect to an arbitrary initial distance ∥0−∗∥ ^0-v^* , this is insufficient for the lower bound proof in this paper. We are working on a lower bound analysis with a more specific problem subclass, where the algorithms may be designed with prior knowledge of a bound D on the distance to the optimum. To this end, we present a refined version of Nesterov’s proposition by properly rescaling Nesterov’s original construction to ensure that the distance to the optimum is bounded by D, while maintaining the identical lower bound on the function value. For any given matrix ∈ℝm×nA ^m× n and vector ∈ℝmb ^m, we define the k-th Krylov subspace as follows: ℋk(,)≜span⊤,(⊤)⊤,…,(⊤)k−1⊤.H^k(A,b) \A b,(A A)A b,…,(A A)^k-1A b \. Proposition 8. Let L>0L>0, D>0D>0, and integer 1≤k≤minm−1,n−121≤ k≤ \m-1,n\-12. Then, there exist a matrix =(L,k)∈ℝm×nA=A(L,k) ^m× n with ∥≤L ≤ L and a vector =(L,D,k)∈ℝmb=b(L,D,k) ^m, such that min∈ℋk(,)12∥−∥2≥3L2D232(k+1)2, _v ^k(A,b)\ 12 -b ^2≥ 3L^2D^232(k+1)^2, and the linear system =Av=b has a solution ∗∈ℝnv^* ^n satisfying ∥∗∥≤D ^* ≤ D. Proof. Let p=2k+1p=2k+1. By the condition k≤minm−1,n−12k≤ \m-1,n\-12, we have p+1≤mp+1≤ m and p≤np≤ n. Let p∈ℝp×pM_p ^p× p be the symmetric tridiagonal matrix defined as: p=[2−10⋯0−12−1⋯00−12⋱⋮⋱−100⋯−12].M_p= bmatrix2&-1&0&·s&0\\ -1&2&-1&·s&0\\ 0&-1&2& & \\ & & & &-1\\ 0&0&·s&-1&2 bmatrix. Let p∈ℝ(p+1)×pB_p ^(p+1)× p be the matrix such that p⊤p=pB_p B_p=M_p, defined as: p=[100⋯0−110⋯00−11⋱⋮⋱000⋯−1100⋯0−1].B_p= bmatrix1&0&0&·s&0\\ -1&1&0&·s&0\\ 0&-1&1& & \\ & & & &0\\ 0&0&·s&-1&1\\ 0&0&·s&0&-1 bmatrix. We define the matrix ∈ℝm×nA ^m× n as the block matrix: =L2[p(p+1)×(n−p)(m−p−1)×p(m−p−1)×(n−p)].A= L2 bmatrixB_p&0_(p+1)×(n-p)\\ 0_(m-p-1)× p&0_(m-p-1)×(n-p) bmatrix. The matrix ⊤∈ℝn×nA A ^n× n is given by the block-diagonal matrix: ⊤=L24[p].A A= L^24 bmatrixM_p&0\\ 0&0 bmatrix. The spectral norm of A is bounded as ∥=λmax(⊤)=L2λmax(p)≤L = _ (A A)= L2 _ (M_p)≤ L. Let γ=D6(p+1)p(2p+1)γ=D 6(p+1)p(2p+1). We define the vector ∈ℝp+1u ^p+1 by its coordinates: u1=p+1,andui=−1p+1for 2≤i≤p+1.u_1= pp+1, u_i=- 1p+1 2≤ i≤ p+1. By the structure of p⊤B_p , we have p⊤=1(p)∈ℝpB_p u=e_1^(p) ^p. Its squared norm evaluates to: ∥2=(p+1)2+p(−1p+1)2=p2+p(p+1)2=p+1. ^2= ( pp+1 )^2+p (- 1p+1 )^2= p^2+p(p+1)^2= pp+1. We define ∈ℝmb ^m as the block vector: =γL2[m−p−1].b=γ L2 bmatrixu\\ 0_m-p-1 bmatrix. The linear system =Av=b has a solution ∗∈ℝnv^* ^n given by the block vector: ∗=γ[p−11(p)n−p].v^*=γ bmatrixM_p^-1e_1^(p)\\ 0_n-p bmatrix. By the structure of p−1M_p^-1, the coordinates of ∗v^* are vi∗=γp+1−ip+1v_i^*=γ p+1-ip+1 for 1≤i≤p1≤ i≤ p, and 0 otherwise. Its squared norm evaluates to: ∥∗∥2=γ2(p+1)2∑j=1pj2=γ2p(2p+1)6(p+1)=D2. ^* ^2= γ^2(p+1)^2 _j=1^pj^2=γ^2 p(2p+1)6(p+1)=D^2. The Krylov subspace is defined as: ℋk(,)=span⊤,(⊤)⊤,…,(⊤)k−1⊤.H^k(A,b)=span\A b,(A A)A b,…,(A A)^k-1A b\. Since p⊤=1(p)B_p u=e_1^(p), the initial vector evaluates to ⊤=γL241(n)A b=γ L^24e_1^(n). Successive multiplication of 1(n)e_1^(n) by ⊤A A expands the non-zero support by one standard basis vector at a time. Thus, the subspace spans the first k standard basis vectors: ℋk(,)=span1(n),2(n),…,k(n).H^k(A,b)=span\e_1^(n),e_2^(n),…,e_k^(n)\. For any ∈ℋk(,)v ^k(A,b), its non-zero support is confined to the first k coordinates. Let k∈ℝkv_k ^k denote these first k coordinates, such that the first p coordinates of v are p=(k⊤,p−k⊤)⊤v_p=(v_k ,0_p-k ) . Let k∈ℝk×kM_k ^k× k be the leading principal submatrix of pM_p. The squared residual norm for ∈ℋk(,)v ^k(A,b) evaluates to: 12∥−∥2 12 -b ^2 =L28∥pp−γ∥2 = L^28 _pv_p- ^2 =L28(⟨p,p⊤pp⟩−2γ⟨p,p⊤⟩+γ2⟨,⟩) = L^28 ( _p,B_p B_pv_p -2γ _p,B_p u +γ^2 ,u ) =L28(⟨k,kk⟩−2γ⟨k,1(k)⟩+γ2p+1). = L^28 ( _k,M_kv_k -2γ _k,e_1^(k) +γ^2 pp+1 ). Denoting v1=⟨k,1(k)⟩v_1= _k,e_1^(k) , the norm simplifies to: 12∥−∥2=L28(⟨k,kk⟩−2γv1+γ2p+1). 12 -b ^2= L^28 ( _k,M_kv_k -2γ v_1+γ^2 pp+1 ). Minimizing this residual norm over k∈ℝkv_k ^k yields the optimal solution k∗=γk−11(k)v_k^*= _k^-1e_1^(k), equivalently characterized by kk∗=γ1(k)M_kv_k^*= _1^(k). By the structure of k−1M_k^-1, the coordinates of k∗v_k^* are vk,i∗=γk+1−ik+1v_k,i^*=γ k+1-ik+1 for 1≤i≤k1≤ i≤ k. Substituting k=k∗v_k=v_k^* into the kv_k-dependent terms (and noting v1=⟨k∗,1(k)⟩=vk,1∗v_1= _k^*,e_1^(k) =v_k,1^*), the minimum value evaluates to: (⟨k,kk⟩−2γv1)|k=k∗ ( _k,M_kv_k -2γ v_1 ) |_v_k=v_k^* =⟨k∗,kk∗⟩−2γvk,1∗ = _k^*,M_kv_k^* -2γ v_k,1^* =⟨k∗,γ1(k)⟩−2γvk,1∗ = _k^*, _1^(k) -2γ v_k,1^* =γvk,1∗−2γvk,1∗=−γvk,1∗=−γ2k+1. =γ v_k,1^*-2γ v_k,1^*=-γ v_k,1^*=-γ^2 kk+1. Substituting these values gives the minimum residual norm over the Krylov subspace: min∈ℋk(,)12∥−∥2 _v ^k(A,b) 12 -b ^2 =L2γ28(p+1−k+1) = L^2γ^28 ( pp+1- kk+1 ) =L2γ28(2k+12k+2−k+1) = L^2γ^28 ( 2k+12k+2- kk+1 ) =L2γ216(k+1) = L^2γ^216(k+1) =L2D216(k+1)6(2k+2)(2k+1)(4k+3) = L^2D^216(k+1) 6(2k+2)(2k+1)(4k+3) =3L2D24(8k2+10k+3). = 3L^2D^24(8k^2+10k+3). Since 8k2+10k+3≤8(k+1)28k^2+10k+3≤ 8(k+1)^2, this residual norm is lower bounded by 3L2D232(k+1)2 3L^2D^232(k+1)^2. ∎ 6.2.2 Quadratic subclasses x P_x and y P_y In this section, we apply Proposition 8 to establish lower bounds for the quadratic subclasses x P_x and y P_y. It should be noted that to establish lower bounds for quadratic subclasses, only the assumptions on local variables (Assumption 1) and candidate solutions (Assumption 3) are used. For any instance Fx∈xF_x∈ P_x of the form Fx()F_x(x), the partial gradient with respect to y is identically zero, and the gradient ∇xFx() _xF_x(x) depends only on the local variable x. As a result, the communication between agents provides no additional information, and the restriction on Agent x reduces to the standard gradient-span condition for single-node unconstrained convex optimization. That is, for this subclass, Assumptions 1 and 3 become: t,l ^t,l ∈0+x−1span∇xFx(i,j)∣0≤i<t, 0≤j≤τxi−1 or i=t, 0≤j≤l−1, ^0+P_x^-1span \ _xF_x(x^i,j) 0≤ i<t,\ 0≤ j≤ _x^i-1 or i=t,\ 0≤ j≤ l-1 \, ¯t+1 x^t+1 ∈0+x−1span∇xFx(i,j)∣0≤i≤t, 0≤j≤τxi−1. ^0+P_x^-1span \ _xF_x(x^i,j) 0≤ i≤ t,\ 0≤ j≤ _x^i-1 \. Notice that the gradient span sequences are naturally confined to the standard Krylov subspaces. Specifically, for x∈ℝnx×nxA_x ^n_x× n_x and x∈ℝnxb_x ^n_x, the subspace satisfies the algebraic progression: (x⊤x)ℋm(x,x)⊆ℋm+1(x,x)andx⊤x∈ℋm+1(x,x).(A_x A_x)H^m(A_x,b_x) ^m+1(A_x,b_x) _x b_x ^m+1(A_x,b_x). Proposition 9. Let ℳ M be a distributed method satisfying Assumption 1 and Assumption 3, operating on an instance Fx∈xF_x∈ P_x of the form Fx()=12∥x−x∥2F_x(x)= 12 _xx-b_x ^2. Let ¯ x be the candidate solution generated by Agent x after evaluating NxN_x partial gradients. Then, ¯∈ℋNx(x,x) x ^N_x(A_x,b_x). Proof. For Fx∈xF_x∈ P_x, we have ψx=ψy=0 _x= _y=0. Let 1,…,Nxz_1,…,z_N_x be the sequence of query points evaluated by Agent x, where i=(i,i)z_i=(x_i,y_i). Since 0=x^0=0, Assumption 1 and Assumption 3 require that each local query point m+1x_m+1 and the candidate solution ¯ x reside in the span of historical gradients. Let Sm≜span∇xFx(i)∣1≤i≤m⊆ℰxS_m \ _xF_x(z_i) 1≤ i≤ m\ _x. We have m+1∈Smx_m+1∈ S_m and ¯∈SNx x∈ S_N_x. We show Sm⊆ℋm(x,x)S_m ^m(A_x,b_x) by induction. The base case m=0m=0 holds since S0=⊆ℋ0(x,x)S_0=\0\ ^0(A_x,b_x). Assume Sm⊆ℋm(x,x)S_m ^m(A_x,b_x) for some m≥0m≥ 0. The gradient at m+1x_m+1 evaluates to ∇xFx(m+1)=x⊤xm+1−x⊤x _xF_x(x_m+1)=A_x A_xx_m+1-A_x b_x. Since m+1∈Sm⊆ℋm(x,x)x_m+1∈ S_m ^m(A_x,b_x), applying the algebraic progression properties yields: ∇xFx(m+1)∈(x⊤x)ℋm(x,x)−x⊤x⊆ℋm+1(x,x). _xF_x(x_m+1)∈(A_x A_x)H^m(A_x,b_x)-A_x b_x ^m+1(A_x,b_x). Thus, Sm+1=Sm+span∇xFx(m+1)⊆ℋm+1(x,x)S_m+1=S_m+span\ _xF_x(x_m+1)\ ^m+1(A_x,b_x). By induction, SNx⊆ℋNx(x,x)S_N_x ^N_x(A_x,b_x), which implies ¯∈ℋNx(x,x) x ^N_x(A_x,b_x). ∎ Analogously, for any instance Fy∈yF_y∈ P_y evaluated by Agent y, of the form Fy(,)=−12∥y−y∥2F_y(x,y)=- 12 _yy-b_y ^2, the candidate solution generated after NyN_y queries satisfies ¯∈ℋNy(y,y) y ^N_y(A_y,b_y). Theorem 10. Let ℳ M be a distributed method satisfying Assumption 1 and Assumption 3 for problem class SP P_SP, where nx≥23LxDx232ϵ+2n_x≥ 2 3L_xD_x^232ε+2 and ny≥23LyDy232ϵ+2n_y≥ 2 3L_yD_y^232ε+2. Then, we have Nxℳ≥cx(3LxDx232ϵ−1),N M_ P_x≥ c_x ( 3L_xD_x^232ε-1 ), and analogously, Nyℳ≥cy(3LyDy232ϵ−1).N M_ P_y≥ c_y ( 3L_yD_y^232ε-1 ). Proof. We prove the bound for x P_x. Let K=⌊3LxDx232ϵ⌋K= 3L_xD_x^232ε . Since nx≥2K+2n_x≥ 2K+2, we apply Proposition 8 to obtain K=(Lx,K)A_K=A( L_x,K) and K=(Lx,Dx,K)b_K=b( L_x,D_x,K) constructed in ℝnx×nxR^n_x× n_x and ℝnxR^n_x. We set x=KA_x=A_K, x=Kb_x=b_K, and define FK(,)=12∥x−x∥2F_K(x,y)= 12 _xx-b_x ^2. By Proposition 8, ∥x∥2≤Lx _x ^2≤ L_x, and the linear system x=xA_xx=b_x has a solution ∗x^* satisfying ∥∗∥≤Dx ^* ≤ D_x. Hence, FK∈xF_K∈ P_x. When ℳ M is applied to PKP_K, suppose it generates an ϵε-saddle point ¯ z utilizing NxN_x queries to x O_x. If Nx≥KN_x≥ K, the bound Nx≥3LxDx232ϵ−1N_x≥ 3L_xD_x^232ε-1 holds. If Nx<KN_x<K, Proposition 9 implies ¯∈ℋNx(x,x) x ^N_x(A_x,b_x). Bounding the restricted duality gap via Proposition 8 yields: ϵ≥ΔFK(¯)=12∥x¯−x∥2≥min∈ℋNx(x,x)12∥x−x∥2≥3LxDx232(Nx+1)2.ε≥ _F_K( z)= 12 _x x-b_x ^2≥ _v ^N_x(A_x,b_x) 12 _xv-b_x ^2≥ 3L_xD_x^232(N_x+1)^2. Rearranging gives Nx≥3LxDx232ϵ−1N_x≥ 3L_xD_x^232ε-1. Multiplying by cxc_x yields Nxℳ≥cx(3LxDx232ϵ−1)N M_ P_x≥ c_x ( 3L_xD_x^232ε-1 ). The proof for y P_y is symmetric, setting y=(Ly,K)A_y=A( L_y,K) and y=(Ly,Dy,K)b_y=b( L_y,D_y,K) for K=⌊3LyDy232ϵ⌋K= 3L_yD_y^232ε . ∎ 6.2.3 Bilinear subclass xy P_xy In this section, we establish the lower bounds for the bilinear subclass xy P_xy. It should be noted that to establish lower bounds for the bilinear subclass, only the assumptions on remote variables (Assumption 2) and candidate solutions (Assumption 3) are used. Consider any function from xy P_xy of the form F(,)=⟨−,⟩F(x,y)= -b,y , where ∈ℝny×nxA ^n_y× n_x and ∈ℝnyb ^n_y. The partial gradients are ∇xF(,)=⊤ _xF(x,y)=A y and ∇yF(,)=− _yF(x,y)=Ax-b. For any ¯=(¯,¯)∈ℰx×ℰy z=( x, y) _x×E_y, the restricted primal-dual gap evaluates to: ΔF(¯)≥max∈ℬy⟨¯−,⟩=Dy∥¯−∥. _F( z)≥ _y _y\ x-b,y =D_y x-b . We define the coupled Krylov subspaces in ℰyE_y and ℰxE_x as: ℋyk(,)≜span,(⊤),…,(⊤)k−1,H^k_y(A,b) \b,(AA )b,…,(AA )^k-1b \, ℋxk(,)≜span⊤,(⊤)⊤,…,(⊤)k−1⊤≡ℋk(,).H^k_x(A,b) \A b,(A A)A b,…,(A A)^k-1A b \ ^k(A,b). By definition, these subspaces satisfy the alternating properties: ℋxm(,)+span=ℋym+1(,)and⊤ℋym(,)=ℋxm(,)⊆ℋxm+1(,).AH^m_x(A,b)+span\b\=H^m+1_y(A,b) H^m_y(A,b)=H^m_x(A,b) ^m+1_x(A,b). Proposition 11. Let ℳ M be a distributed method satisfying Assumption 2 and Assumption 3, operating on F(,)=⟨−,⟩F(x,y)= -b,y . Let ¯ x be the candidate solution generated after T communication rounds. Then, ¯∈ℋx⌈(T−1)/2⌉(,) x (T-1)/2 _x(A,b). Proof. Let Sxt≜span∇xF()∣∈Zxt⊆ℰxS_x^t \ _xF(z) ∈ Z_x^t\ _x and Syt≜span−∇yF()∣∈Zyt⊆ℰyS_y^t \- _yF(z) ∈ Z_y^t\ _y. We show by induction that Sxt⊆ℋx⌈t/2⌉(,)S_x^t t/2 _x(A,b) and Syt⊆ℋy⌊t/2⌋+1(,)S_y^t t/2 +1_y(A,b). The base case t=−1t=-1 holds since Sx−1=⊆ℋx0(,)S_x^-1=\0\ ^0_x(A,b) and Sy−1=⊆ℋy0(,)S_y^-1=\0\ ^0_y(A,b). Assume the claim holds for round t−1t-1. In round t, Agent x queries points using remote variables y. By Assumption 2, ^∈Syt−1 y∈ S_y^t-1. Applying the alternating properties, the evaluated gradient satisfies: ∇xF(,^)=⊤^∈⊤ℋy⌊(t−1)/2⌋+1(,)=ℋx⌊(t−1)/2⌋+1(,)=ℋx⌈t/2⌉(,). _xF(x, y)=A y H (t-1)/2 +1_y(A,b)=H (t-1)/2 +1_x(A,b)=H t/2 _x(A,b). Thus, Sxt⊆ℋx⌈t/2⌉(,)S_x^t t/2 _x(A,b). Similarly, Agent y queries points using remote variables x. By Assumption 2, ^∈Sxt−1 x∈ S_x^t-1. The evaluated gradient satisfies: −∇yF(^,)=−^∈span+ℋx⌈(t−1)/2⌉(,)=ℋy⌈(t−1)/2⌉+1(,)=ℋy⌊t/2⌋+1(,).- _yF( x,y)=b-A x \b\+AH (t-1)/2 _x(A,b)=H (t-1)/2 +1_y(A,b)=H t/2 +1_y(A,b). Thus, Syt⊆ℋy⌊t/2⌋+1(,)S_y^t t/2 +1_y(A,b). By induction, the claim holds for all t. By Assumption 3, the candidate solution satisfies ¯∈SxT−1⊆ℋx⌈(T−1)/2⌉(,) x∈ S_x^T-1 (T-1)/2 _x(A,b). ∎ Theorem 12. Let ℳ M be a distributed method satisfying Assumption 2 and Assumption 3 for problem class xy P_xy, where minnx,ny≥2LxyDxDy3ϵ+2 \n_x,n_y\≥ 2L_xyD_xD_y3ε+2. Then, we have Txyℳ≥2LxyDxDy3ϵ−2.T M_ P_xy≥ 2L_xyD_xD_y3ε-2. Proof. Let K=⌊LxyDxDy3ϵ⌋K= L_xyD_xD_y3ε . Since nx,ny≥2K+2n_x,n_y≥ 2K+2, we apply Proposition 8 to obtain K=(Lxy,K)∈ℝny×nxA_K=A(L_xy,K) ^n_y× n_x and K=(Lxy,Dx,K)∈ℝnyb_K=b(L_xy,D_x,K) ^n_y. We set xy=KA_xy=A_K and xy=Kb_xy=b_K, and define FK(,)=⟨xy−xy,⟩F_K(x,y)= _xyx-b_xy,y . By Proposition 8, ∥xy∥=∥K∥≤Lxy _xy = _K ≤ L_xy. The linear system xy=xyA_xyx=b_xy has a solution ∗x^* satisfying ∥∗∥≤Dx ^* ≤ D_x. Thus, FK∈xyF_K∈ P_xy. When ℳ M is applied to PKP_K, suppose it generates an ϵε-saddle point (¯,¯)( x, y) after T communication rounds. If T≥2KT≥ 2K, then T≥2LxyDxDy3ϵ−2T≥ 2L_xyD_xD_y3ε-2 holds. If T<2KT<2K, Proposition 11 implies ¯∈ℋxk(xy,xy) x ^k_x(A_xy,b_xy), where k=⌈(T−1)/2⌉k= (T-1)/2 . Since ℋxk(xy,xy)≡ℋk(K,K)H^k_x(A_xy,b_xy) ^k(A_K,b_K), the restricted duality gap on PKP_K evaluates to: ΔFK(¯,¯)≥Dy∥K¯−K∥. _F_K( x, y)≥ D_y _K x-b_K . Applying Proposition 8 yields: ϵ≥ΔFK(¯,¯)≥Dymin∈ℋk(K,K)∥K−K∥≥Dy3Lxy2Dx22(8k2+10k+3)≥LxyDxDy3(k+1).ε≥ _F_K( x, y)≥ D_y _v ^k(A_K,b_K) _Kv-b_K ≥ D_y 3L_xy^2D_x^22(8k^2+10k+3)≥ L_xyD_xD_y3(k+1). Rearranging gives k≥LxyDxDy3ϵ−1k≥ L_xyD_xD_y3ε-1. Because k=⌈(T−1)/2⌉k= (T-1)/2 , we have T≥2k≥2LxyDxDy3ϵ−2T≥ 2k≥ 2L_xyD_xD_y3ε-2. ∎ The oracle lower bounds can be proved using identical subspace confinement arguments. Rather than repeating the proof, we provide the main proposition and theorem here. Proposition 13. Let ℳ M be a distributed method satisfying Assumption 2 and Assumption 3, operating on F(,)=⟨−,⟩F(x,y)= -b,y . Let (¯,¯)( x, y) be the solution generated after NxN_x and NyN_y queries to x O_x and y O_y, respectively. Then, ¯∈ℋxk(,) x ^k_x(A,b), where k=min(Nx,Ny)k= (N_x,N_y). Theorem 14. Let ℳ M be a distributed method satisfying Assumption 2 and Assumption 3 for problem class xy P_xy, where minnx,ny≥2LxyDxDy3ϵ+2 \n_x,n_y\≥ 2L_xyD_xD_y3ε+2. Then, we have Nxyℳ≥(cx+cy)(LxyDxDy3ϵ−1).N M_ P_xy≥(c_x+c_y) ( L_xyD_xD_y3ε-1 ). 6.2.4 Lower complexity bounds Finally, we derive the lower bound for SP P_SP by assembling the lower bounds obtained from the three subclasses. Theorem 15. Let ℳ M be a distributed gradient-span algorithm for SP P_SP. Suppose the dimensions of ℰxE_x and ℰyE_y are sufficiently large such that nx≥2LxyDxDy3ϵ+3LxDx28ϵ+2n_x≥ 2L_xyD_xD_y3ε+ 3L_xD_x^28ε+2 and ny≥2LxyDxDy3ϵ+3LyDy28ϵ+2n_y≥ 2L_xyD_xD_y3ε+ 3L_yD_y^28ε+2. We have: TSPℳ≥2LxyDxDy3ϵ−2,NSPℳ≥cx+cy9LxyDxDyϵ+cx33LxDx232ϵ+cy33LyDy232ϵ−2cx+2cy3. gatheredT M_ P_SP≥ 2L_xyD_xD_y3ε-2,\\ N M_ P_SP≥ c_x+c_y9 L_xyD_xD_yε+ c_x3 3L_xD_x^232ε+ c_y3 3L_yD_y^232ε- 2c_x+2c_y3. gathered Proof of Theorem 15. The general class of convex-concave saddle point problems SP P_SP contains the three unregularized subclasses constructed in the previous sections: the x-quadratic subclass x P_x, the y-quadratic subclass y P_y, and the bilinear subclass xy P_xy. The worst-case complexity for an algorithm operating over the entire class SP P_SP is bounded from below by the maximum of the complexities required for these individual subclasses. Because a gradient-span algorithm satisfies all three assumptions (Assumptions 1, 2 and 3), we can assemble these results. By Theorem 12, the communication complexity over the class is bounded by the communication complexity of the bilinear subclass: TSPℳ≥Txyℳ≥2LxyDxDy3ϵ−2.T M_ P_SP≥ T M_ P_xy≥ 2L_xyD_xD_y3ε-2. For the computational complexity, we combine the independent lower bounds established in Theorem 10 and Theorem 12. The total computational complexity is bounded by the maximum of the three individual requirements: NSPℳ≥maxNxyℳ,Nxℳ,Nyℳ.N M_ P_SP≥ \N M_ P_xy,N M_ P_x,N M_ P_y \. Using the algebraic property maxa,b,c≥13(a+b+c) \a,b,c\≥ 13(a+b+c), we obtain the lower bound: NSPℳ N M_ P_SP ≥13[Nxyℳ+Nxℳ+Nyℳ] ≥ 13 [N M_ P_xy+N M_ P_x+N M_ P_y ] ≥13[(cx+cy)(LxyDxDy3ϵ−1)+cx(3LxDx232ϵ−1)+cy(3LyDy232ϵ−1)] ≥ 13 [(c_x+c_y) ( L_xyD_xD_y3ε-1 )+c_x ( 3L_xD_x^232ε-1 )+c_y ( 3L_yD_y^232ε-1 ) ] =cx+cy9LxyDxDyϵ+cx33LxDx232ϵ+cy33LyDy232ϵ−2cx+2cy3. = c_x+c_y9 L_xyD_xD_yε+ c_x3 3L_xD_x^232ε+ c_y3 3L_yD_y^232ε- 2c_x+2c_y3. This establishes the stated lower bounds for the general problem class and concludes the proof. ∎ Remark 4 (Communication optimality). Theorem 15 confirms that the communication lower bound for distributed SPs depends only on the cross-coupled conditioning. More importantly, the communication lower bound perfectly matches our upper bound proven in Theorem 1 up to a constant. Therefore, we conclude that our DM-SP is a communication-optimal algorithm within the gradient-span framework. We note, however, that a gap remains between the achieved oracle costs and the theoretical oracle lower bound established in Theorem 15. Yet this is a known open question even for non-distributed SPs. 7 Variational inequality problems with distributed oracles Motivation. Thus far, we have studied SPs, which naturally model two-player zero-sum games. To capture more complex multiagent interactions (such as equilibrium computation in multiplayer general-sum games, network routing, and multiagent reinforcement learning), we extend our algorithmic framework to the broader class of monotone Variational Inequality Problems (VIPs) with separable composite terms and distributed oracles. We briefly outline the problem class and our results here, deferring the detailed presentation to Appendix C. Problem class. We consider a distributed multiagent setting with a star communication network over a product space ℰ=ℰ1×⋯×ℰKE=E_1×·s×E_K, where a joint decision variable is partitioned among K autonomous agents as =(1,…,K)z=(z_1,…,z_K). Let us consider the problem class VIP P_VIP as follows: • Operators and local components: The problem is governed by a joint monotone operator V()=(V1(),⋯,VK())V(z)=(V_1(z),·s,V_K(z)) and a separable local composite function ψ()=∑i=1Kψi(i)ψ(z)= _i=1^K _i(z_i) defined on Q=domψ1×⋯×domψKQ=dom _1×·s×dom _K. We assume the solution set is bounded by local distance parameters Di>0D_i>0 for each agent. Furthermore, the operator satisfies block-wise Lipschitz continuity: for any fixed −jz_-j, the mapping Vi(j;−j)V_i(z_j;z_-j) is LijL_ij-Lipschitz continuous with respect to jz_j. • Distributed oracles: Each Agent i∈[K]i∈[K] controls its local variable i∈domψiz_i _i, has access to its private function ψi _i and a partial oracle i()=Vi() O_i(z)=V_i(z). • Accuracy measure: The goal is to find an ϵε-approximate solution ¯∈Q z∈ Q such that the restricted gap Δ(¯)≜sup∈ℬ∩Q[⟨V(),¯−⟩+ψ(¯)−ψ()] ( z) _z ∩ Q[ V(z), z-z +ψ( z)-ψ(z)] satisfies Δ(¯)≤ϵ ( z)≤ε, where the bounded domain ℬB is defined by the initial point 0z^0 and distance parameters DiD_i, i∈[K]i∈[K]. Conditionings. Let L¯ij≜maxLij,Lji L_ij \L_ij,L_ji\. Similar to SPs, we hereby define the cross-coupled conditioning, denoted by ∑i∈[K]Ai _i∈[K]A_i with Ai≜Di∑j≠iL¯ijDjA_i D_i _j≠ i L_ijD_j, that quantifies the cross-dependencies between the agents over the network. The diagonal conditioning, denoted by ∑i∈[K]Bi _i∈[K]B_i with Bi≜LiiDi2B_i L_iD_i^2, measures the self-dependency within a single agent’s domain. New state-of-the-art communication cost. In the multiagent setting, Extragradient (EG) remains the state-of-the-art method, requiring TVIPEG=(∑i∈[K]Ai+Biϵ)T^EG_ P_VIP=O ( _i∈[K] A_i+B_iε ) communication cost, which depends on both the cross-coupled and diagonal conditionings. By extending our decoupled template to distributed VIPs, we propose DM-VIP, and establish a much better communication cost. Theorem 16. Consider the DM-VIP algorithm for problem class VIP P_VIP. We have TVIPDM-VIP=(∑i∈[K]Ai/ϵ).T^DM-VIP_ P_VIP=O ( _i∈[K]A_i/ε ). Remark 5. The communication cost in Theorem 16 completely drops the dependence on the diagonal conditioning. Consequently, we have substantially improved the state-of-the-art EG cost when the diagonal conditioning dominates, i.e., when ∑iBi≫∑iAi _iB_i _iA_i. 8 Conclusion and limitations This paper studies communication and oracle costs in distributed SPs and VIPs. For the class of SPs, we settle the communication complexity in the distributed setup within gradient-span framework, and consistently improve the long-standing oracle cost of EG method. For the class of distributed VIPs, we improve the state-of-the-art communication cost. The following directions are not addressed in this paper and are left for future work: (a) closing the gap of oracle costs; (b) showing lower bound for non-zero-sum games; and (c) showing information-theoretic lower bounds for randomized methods. References A. Ben-Tal and A. Nemirovski (2002) Robust optimization–methodology and applications. Mathematical programming 92 (3), p. 453–480. Cited by: §1. A. Beznosikov, V. Samokhin, and A. Gasnikov (2025) Distributed saddle point problems: lower bounds, near-optimal and robust algorithms. Optimization Methods and Software, p. 1–18. Cited by: §1. R. I. Boţ and E. Chenchene (2024) Extra-gradient method with flexible anchoring: strong convergence and fast residual decay. arXiv preprint arXiv:2410.14369. Cited by: §C.4.2. V. Conitzer and T. Sandholm (2004) Communication complexity as a lower bound for learning in games. In Proceedings of the twenty-first international conference on Machine learning, p. 24. Cited by: §1, §1. Y. Deng and M. Mahdavi (2021) Local stochastic gradient descent ascent: convergence analysis and communication efficiency. In International Conference on Artificial Intelligence and Statistics, p. 1387–1395. Cited by: §1. I. J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio (2014) Generative adversarial nets. Advances in neural information processing systems 27. Cited by: §1, §1. S. Hart and Y. Mansour (2010) How long to equilibrium? the communication complexity of uncoupled equilibrium procedures. Games and Economic Behavior 69 (1), p. 107–126. Cited by: §1. J. Hu, M. P. Wellman, et al. (1998) Multiagent reinforcement learning: theoretical framework and an algorithm.. In ICML, Vol. 98, p. 242–250. Cited by: §1. A. Juditsky, A. Nemirovski, et al. (2011) First order methods for nonsmooth convex large-scale optimization, i: utilizing problems structure. Optimization for Machine Learning 30 (9), p. 149–183. Cited by: Appendix C, §2.3, Proposition 17, Proposition 20. G. Lan and Y. Li (2026) A novel catalyst scheme for stochastic minimax optimization. Mathematical Programming, p. 1–49. Cited by: §2.3, Proposition 18. G. Lan, Y. Ouyang, and Z. Zhang (2023) Optimal and parameter-free gradient minimization methods for convex and nonconvex optimization. arXiv preprint arXiv:2310.12139. Cited by: Appendix B, §4.4, §4.4, Lemma 5. T. Lin, C. Jin, and M. I. Jordan (2020) Near-optimal algorithms for minimax optimization. In Conference on learning theory, p. 2738–2779. Cited by: §2.3. B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas (2017) Communication-efficient learning of deep networks from decentralized data. In Artificial intelligence and statistics, p. 1273–1282. Cited by: §1, §1. R. D. Monteiro and B. F. Svaiter (2013) An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods. SIAM Journal on Optimization 23 (2), p. 1092–1125. Cited by: §4.2. A. Nemirovski (2004) Prox-method with rate of convergence O(1/t)O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems. SIAM Journal on Optimization 15 (1), p. 229–251. Cited by: Appendix C, §1, §2.3, §4.1, §4.1. A. S. Nemirovskij and D. B. Yudin (1983) Problem complexity and method efficiency in optimization. Wiley-Interscience. Cited by: §2.1. Y. E. Nesterov (2004) Introductory lectures on convex optimization - A basic course. Applied Optimization, Vol. 87, Springer. External Links: Link, Document, ISBN 978-1-4613-4691-3 Cited by: §6.2.1, §6. Y. Nesterov (2023) High-order reduced-gradient methods for composite variational inequalities. arXiv preprint arXiv:2311.15154. Cited by: §3, §4.1, §4.1, §4.2. N. Nisan and I. Segal (2006) The communication requirements of efficient allocations and supporting prices. Journal of Economic Theory 129 (1), p. 192–224. Cited by: §1. J. B. Rosen (1965) Existence and uniqueness of equilibrium points for concave n-person games. Econometrica: Journal of the Econometric Society, p. 520–534. Cited by: §1. M. Schmidt, N. Le Roux, and F. Bach (2017) Minimizing finite sums with the stochastic average gradient. Mathematical Programming 162 (1), p. 83–112. Cited by: §1. P. Tseng (1995) On linear convergence of iterative methods for the variational inequality problem. Journal of Computational and Applied Mathematics 60 (1-2), p. 237–252. Cited by: §1. J. von Neumann and O. Morgenstern (1947) Theory of games and economic behavior. Princeton university press. Cited by: §1. Y. Wang and J. Li (2020) Improved algorithms for convex-concave minimax optimization. Advances in Neural Information Processing Systems 33, p. 4800–4810. Cited by: Appendix A, §2.3, Proposition 19. J. Yang, S. Zhang, N. Kiyavash, and N. He (2020) A catalyst framework for minimax optimization. Advances in Neural Information Processing Systems 33, p. 5667–5678. Cited by: §2.3, Proposition 18. T. Yoon, S. Choudhury, and N. Loizou (2025) Multiplayer federated learning: reaching equilibrium with less communication. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, External Links: Link Cited by: §1, §1, §2.3. T. Yoon and N. Loizou (2025) PEARL-prox: proximal algorithm for resolving player drift in multiplayer federated learning. In OPT 2025: Optimization for Machine Learning, Cited by: §1, §2.3. S. Zhang, S. Choudhury, S. U. Stich, and N. Loizou (2024) Communication-efficient gradient descent-ascent methods for distributed variational inequalities: unified analysis and local updates. In The Twelfth International Conference on Learning Representations, Cited by: §1, §1, §2.3. A. Zindari, P. Yazdkhasti, A. Rodomanov, T. Chavdarova, and S. U. Stich (2025) Decoupled SGDA for games with intermittent strategy communication. In Forty-second International Conference on Machine Learning, External Links: Link Cited by: Appendix A, §1, §2.3. Appendix A Detailed review of existing algorithms In this section, we provide the detailed formulations, distributed trajectories, and formal complexity results for the algorithms summarized in Section 2.3. We will also show that all the algorithms described below indeed satisfy the gradient-span assumptions in Definition 3. We recall that SP∘ P_SP denotes the subclass of SP P_SP where the local components ψx _x and ψy _y are zero. In this setting, the problem reduces to finding a saddle point of a smooth convex-concave function f. We note that many algorithms discussed below do not handle general composite functions and only deal with problem instances from this non-composite subclass. As in Remark 2, let us consider a practical scenario where the algorithms may not have the precise values of DxD_x and DyD_y in advance, but they have access to upper estimates D^x≥Dx D_x≥ D_x and D^y≥Dy D_y≥ D_y. Let θ≜DxD^yD^xDy+DyD^xD^yDx, θ D_x D_y D_xD_y+ D_y D_x D_yD_x, which quantifies the disproportionality between the true distance parameters and their estimates. Extragradient (EG). Suppose we apply EG to a non-composite problem instance P=(f,0,0,0)∈SP∘P=(f,0,0,z^0)∈ P _SP. The method iteratively maintains and updates two sequences of variables k,k∈ℰx×ℰyv^k,z^k _x×E_y. Initializing with 0=0v^0=z^0, the updates for iteration k=0,1,…,K−1k=0,1,…,K-1 are given by: k+1 ^k+1 =k−ηk−1Vf(k), =v^k-η^kP^-1V^f(v^k), k+1 ^k+1 =k−ηk−1Vf(k+1), =v^k-η^kP^-1V^f(z^k+1), where ηk>0η^k>0 is the step size at iteration k. Although the EG method is conventionally formulated as a K-iteration loop, its distributed execution requires 2K2K communication rounds. In each iteration, the agents must synchronize twice to evaluate the coupled partial gradients at kv^k and at k+1z^k+1. Thus, the EG method corresponds to a distributed gradient-span algorithm with T=2KT=2K rounds and local step lengths τxt=τyt=1 _x^t= _y^t=1, where the trajectories are given by: Z^x2k=Z^y2k=kandZ^x2k+1=Z^y2k+1=k+1, Z_x^2k= Z_y^2k=\v^k\ Z_x^2k+1= Z_y^2k+1=\z^k+1\, for k=0,…,K−1k=0,…,K-1. To align with Definition 1, the method produces a candidate solution ¯t+1 z^t+1 after each round t∈0,…,T−1t∈\0,…,T-1\. Because each iteration requires two communication rounds, the algorithm effectively updates its output only every two rounds. Specifically, after completing iteration k (which corresponds to round t=2k+1t=2k+1), it outputs the ergodic average ¯2k+2=1∑i=0kηi(∑i=0kηii+1) z^2k+2= 1 _i=0^kη^i ( _i=0^kη^iz^i+1 ). During the intermediate rounds (at t=2kt=2k), it simply retains the previous solution by setting ¯2k+1=¯2k z^2k+1= z^2k (where ¯0=0 z^0=z^0). Since the trajectories and the candidate solutions are formed by linear combinations of the evaluated gradients, they satisfy the span conditions in Definition 3. The complexity of the EG method provides a natural baseline. Translating the classic results into our distributed complexity measures yields the following upper bounds. Proposition 17 (Juditsky et al. 2011, Eq. (6.21)). Consider the EG method applied to SP∘ P _SP. With the parameter choices of αx=LxD^x+LxyD^yD^x _x= L_x D_x+L_xy D_y D_x, αy=LyD^y+LxyD^xD^y _y= L_y D_y+L_xy D_x D_y, and ηk≡1η^k≡ 1, we have TSP∘EG≤θ⋅LxyDxDyϵ+LxDx2ϵ+LyDy2ϵ,NSP∘EG≤(cx+cy)⋅(θ⋅LxyDxDyϵ+LxDx2ϵ+LyDy2ϵ). gatheredT^EG_ P _SP≤θ· L_xyD_xD_yε+ L_xD_x^2ε+ L_yD_y^2ε,\\ N^EG_ P _SP≤(c_x+c_y)· (θ· L_xyD_xD_yε+ L_xD_x^2ε+ L_yD_y^2ε ). gathered Decoupled GDA. Let us apply DGDA to an instance P∈SP∘P∈ P_SP with a fixed local trajectory length τ. The algorithm maintains local iterates t,lx^t,l and t,ly^t,l for round t=0,…,T−1t=0,…,T-1 and local step l=0,…,τl=0,…,τ. At the beginning of round t, the agents synchronize by exchanging their latest local iterates. Specifically, Agent x receives ^t y^t and Agent y receives ^t x^t, defined as: ^t≜0if t=0t−1,τif t>0and^t≜0if t=0t−1,τif t>0. x^t casesx^0&if t=0\\ x^t-1,τ&if t>0 cases y^t casesy^0&if t=0\\ y^t-1,τ&if t>0 cases. The agents then initialize their local variables for the current round as t,0=^tx^t,0= x^t and t,0=^ty^t,0= y^t. With the remote variables firmly fixed, the agents execute τ local gradient steps. For l=0,…,τ−1l=0,…,τ-1, the local updates are given by: t,l+1 ^t,l+1 =t,l−ηxx−1∇xf(t,l,^t), =x^t,l- _xP_x^-1 _xf(x^t,l, y^t), t,l+1 ^t,l+1 =t,l+ηyy−1∇yf(^t,t,l), =y^t,l+ _yP_y^-1 _yf( x^t,y^t,l), where ηx,ηy>0 _x, _y>0 are the local step sizes. The DGDA algorithm yields the following trajectories: Z^xt=(t,l,^t)l=0τ−1andZ^yt=(^t,t,l)l=0τ−1. Z_x^t= \(x^t,l, y^t) \_l=0^τ-1 Z_y^t= \( x^t,y^t,l) \_l=0^τ-1. After each round t∈0,…,T−1t∈\0,…,T-1\, the method returns the updated local iterates as the candidate solution ¯t+1=(t,τ,t,τ) z^t+1=(x^t,τ,y^t,τ). Because the trajectories and the candidate solutions are formed entirely by linear combinations of the evaluated gradients, DGDA is indeed a distributed gradient-span algorithm by Definition 3. While DGDA lies perfectly in our framework, its theoretical guarantees are highly restrictive. The algorithm is only proven to converge for restricted strongly convex-strongly concave problem instances where the cross-coupling between the variables is sufficiently weak [Zindari et al., 2025]. In this narrowly defined regime, DGDA achieves a logarithmic communication complexity of (log1ϵ)O ( 1ε ), which is a clear improvement over the EG baseline. However, for general problem instances in SP∘ P_SP with stronger coupling, the delayed remote variables cause the local updates to drift, ultimately leading the method to diverge. Consequently, DGDA fails to provide any meaningful complexity guarantee for the general problem class SP∘ P_SP under consideration. Catalyst acceleration. Catalyst methods first add small regularization terms to the objective: f(,)+ϵ4D^x2∥−0∥x2−ϵ4D^y2∥−0∥y2f(x,y)+ ε4 D_x^2 -x^0 _x^2- ε4 D_y^2 -y^0 _y^2, reducing the problem to a strongly convex-strongly concave one. Then, Catalyst introduces an outer loop, indexed by k=0,1,…,K−1k=0,1,…,K-1, designed to balance the conditioning between the two variables. For instance, when the conditioning of x is worse (i.e. LxD^x2≥LyD^y2L_x D_x^2≥ L_y D_y^2), the method carefully maintains an extrapolation sequence (~k)k=0K−1( x^k)_k=0^K-1 and, in each outer iteration, adds a proximal term λx2∥−~k∥x2 _x2 - x^k _x^2 to the objective. Conversely, if the conditioning of y is worse, the outer loop would instead maintain an extrapolation sequence for y and add a corresponding regularization term for y. An inner base algorithm (EG in this case) is then deployed to solve this regularized subproblem to a specified accuracy. To simplify the notation, let Lmax=maxLx,Ly,LxyL_ = \L_x,L_y,L_xy\. Proposition 18 (Yang et al. [2020], Lan and Li [2026]). Consider the Catalyst framework equipped with the EG method as the inner solver, denoted by Cat-EG, applied to SP∘ P_SP . Then, we have: TSP∘Cat-EG=((LmaxD^xD^yϵ+LxD^x2ϵ+LyD^y2ϵ)log2(1ϵ)),NSP∘Cat-EG=(cx+cy)⋅((LmaxD^xD^yϵ+LxD^x2ϵ+LyD^y2ϵ)log2(1ϵ)). gatheredT^Cat-EG_ P_SP =O ( ( L_ D_x D_yε+ L_x D_x^2ε+ L_y D_y^2ε ) ^2 ( 1ε ) ),\\ N^Cat-EG_ P_SP =(c_x+c_y)·O ( ( L_ D_x D_yε+ L_x D_x^2ε+ L_y D_y^2ε ) ^2 ( 1ε ) ). gathered Now, let us explain the caveats we mentioned in Section 2.3 regarding the Cat-EG method in more detail. To explain the second caveat of Cat-EG, its sensitivity to inexact diameter estimates, we compare it with the EG baseline and our proposed method. Because Catalyst uses D^x D_x and D^y D_y to set the initial regularization, its complexity scales directly with these estimates rather than the true distances DxD_x and DyD_y. Specifically, the diagonal terms in EG or DM-SP depend strictly on the true distances, whereas in Cat-EG they scale with D^x D_x and D^y D_y. Similarly, for the cross-coupled term, EG or DM-SP depends on the true distances multiplied by the proportionality ratio θ. If the estimates are loose but proportional (e.g., D^x=cDx D_x=cD_x and D^y=cDy D_y=cD_y), θ remains 22, leaving the complexity unaffected by the overestimation factor c. In contrast, the Cat-EG coupled term scales with D^xD^y D_x D_y, meaning any overestimation of D^x≫Dx D_x D_x or D^y≫Dy D_y D_y directly inflates the bound. Consequently, EG and our proposed DM-SP method are much more robust to inexact distance estimates, provided the estimates are roughly proportional. Furthermore, we note that the fourth caveat of Cat-EG (performing worse than unaccelerated EG) can be easily verified. For instance, consider a problem instance where Lx=106L_x=10^6, Lxy=1L_xy=1, Ly=10−6L_y=10^-6, Dx=10−3D_x=10^-3, and Dy=103D_y=10^3. Under this conditioning, the theoretical upper bound of Cat-EG significantly exceeds that of standard EG. Four-loop method. The Cat-Cat-DAGDA method [Wang and Li, 2020] first adds small (ϵ)O(ε) regularizations to both x and y to ensure the objective is strongly convex-strongly concave. It then executes four nested loops, which justifies our naming convention: (i) The first loop is a Catalyst outer loop adding a proximal regularization term to the primal variable; (i) The second loop is another Catalyst outer loop adding a proximal regularization term to the dual variable; (i) The third loop manages communication by exchanging and freezing the remote variables, identical to DGDA; and (iv) The fourth loop performs local computations, but unlike the standard gradient steps in DGDA, it employs an accelerated gradient method (hence DAGDA) to solve the inner subproblems. While the fourth loop of the method is a local computation loop, the first three loops all require communication rounds. Let us state the complexity results of this method without going into the tedious details of the algorithm. Proposition 19 (Wang and Li [2020]). Consider the Cat-Cat-DAGDA method applied to SP∘ P_SP . Then, we have: TSP∘Cat-Cat-DAGDA=(LxyD^xD^yϵlog3(1ϵ)),NSP∘Cat-Cat-DAGDA=(cx+cy)⋅((LmaxLxyD^xD^yϵ+LxD^x2ϵ+LyD^y2ϵ)log4(1ϵ)). gatheredT^Cat-Cat-DAGDA_ P_SP =O ( L_xy D_x D_yε ^3 ( 1ε ) ),\\ N^Cat-Cat-DAGDA_ P_SP =(c_x+c_y)·O ( ( L_ L_xy D_x D_yε+ L_x D_x^2ε+ L_y D_y^2ε ) ^4 ( 1ε ) ). gathered Appendix B MRN solver: Accumulative Regularization Method Algorithm 4 ARM(∇fw,ψw,w,ξ∣L)ARM (∇ f_w, _w,v_w,ξ L ) 1:Set τ=2+max0,⌈log4(3L2ξ)⌉τ=2+ \0, _4 ( 3L2ξ ) \ and σ(0)=0σ^(0)=0. 2:Set σ(k)=4k−32ξ3σ^(k)=4^k-3 2ξ3 and Nk=⌈16Lσ(k)⌉N_k= 16 Lσ^(k) for k=1,…,τk=1,…,τ. 3:Initialize ¯(0)=(0)=w w^(0)=w^(0)=v_w. 4:for k=1,…,τk=1,…,τ do 5: γ(k)=1−σ(k−1)σ(k)γ^(k)=1- σ^(k-1)σ^(k) 6: ¯(k)=(1−γ(k))¯(k−1)+γ(k)(k−1) w^(k)=(1-γ^(k)) w^(k-1)+γ^(k)w^(k-1) 7: % Begin Inner Subroutine: Nesterov’s Accelerated Gradient Method 8: Initialize k(0)=(k−1)x_k^(0)=w^(k-1), k(0)=(k−1)y_k^(0)=w^(k-1), t0=1t_0=1, and Lk=L+σ(k)L_k=L+σ^(k). 9: for i=0,…,Nk−1i=0,…,N_k-1 do 10: ∇fw(k)(k(i))=∇fw(k(i))+σ(k)(k(i)−¯(k))∇ f_w^(k)(y_k^(i))=∇ f_w(y_k^(i))+σ^(k)(y_k^(i)- w^(k)) 11: k(i+1)=argmin∈domψw⟨∇fw(k)(k(i)),⟩+Lk2∥−k(i)∥w2+ψw()x_k^(i+1)= *arg\,min_w _w\ \ ∇ f_w^(k)(y_k^(i)),w + L_k2 -y_k^(i) _w^2+ _w(w) \ 12: ti+1=1+1+4ti22t_i+1= 1+ 1+4t_i^22 13: k(i+1)=k(i+1)+ti−1ti+1(k(i+1)−k(i))y_k^(i+1)=x_k^(i+1)+ t_i-1t_i+1(x_k^(i+1)-x_k^(i)) 14: end for 15: % End Inner Subroutine 16: (k)=k(Nk)w^(k)=x_k^(N_k) 17: ψw′((k))=−∇fw(k)(k(Nk−1))−Lk(k(Nk)−k(Nk−1)) _w (w^(k))=-∇ f_w^(k)(y_k^(N_k-1))-L_k(x_k^(N_k)-y_k^(N_k-1)) 18:end for 19:return ((τ),ψw′((τ)))(w^(τ), _w (w^(τ))) Proof of Lemma 5. It was originally shown in [Lan et al., 2023, Theorem 3.1] that ARM takes no more than 343L2ξ34 3L2ξ gradient queries and obtains ′∈domψww _w, such that there exists ~∈domψw w _w with ∇fw(~)∈−∂ψw(~)∇ f_w( w)∈-∂ _w( w), and ∥2L(+−′)∥w≤23ξ∥w−~∥w, where +=argmin∈domψw[⟨∇fw(′),⟩+ψw()+L∥−′∥w2]. 2L(w^+-w ) _w≤ 23ξ _w- w _w, where w^+= *arg\,min_w _w\ [ ∇ f_w(w ),w + _w(w)+L -w _w^2 ]. By the optimality of +w^+, there exists ψw′(+)∈∂ψw(+)ψ _w(w^+)∈∂ _w(w^+) such that ∇fw(′)+ψw′(+)+2L(+−′)=.∇ f_w(w )+ψ _w(w^+)+2LP(w^+-w )=0. Then, we have ∥∇fw(+)+ψw′(+)∥w∗ ∇ f_w(w^+)+ _w (w^+) _w^* ≤∥∇fw(′)+ψw′(+)∥w∗+∥∇fw(+)−∇fw(′)∥w∗ ≤ ∇ f_w(w )+ψ _w(w^+) _w^*+ ∇ f_w(w^+)-∇ f_w(w ) _w^* =∥2L(+−′)∥w∗+∥∇fw(+)−∇fw(′)∥w∗ = 2LP(w^+-w ) _w^*+ ∇ f_w(w^+)-∇ f_w(w ) _w^* ≤3L∥+−′∥w ≤ 3L ^+-w _w ≤ξ∥w−~∥w. ≤ξ _w- w _w. ∎ Appendix C Monotone composite variational inequality problems In this section, we study variational inequality problems (VIPs) [Nemirovski, 2004, Juditsky et al., 2011], a generalization of SPs that captures, for instance, multiplayer general-sum games. C.1 Problem formulation VIPs (with separable composite terms). Let us consider the VIP in Eq. 4, where ℰ=ℰ1×⋯×ℰKE=E_1×·s×E_K is the direct product of K finite-dimensional real vector spaces. For all i∈[K]i∈[K]: let the mapping Vi:domψ→ℰi∗V_i ψ _i^*, and let the function ψi:ℰi→ℝ∪+∞ _i _i ∪\+∞\. We consider the decomposition of V()=(V1(),⋯,VK())V(z)= (V_1(z),·s,V_K(z) ) and ψ()=ψ1(1)+⋯+ψK(K)ψ(z)= _1(z_1)+·s+ _K(z_K), for all =(1,⋯,K)∈ℰz=(z_1,·s,z_K) . Moreover, we denote domψ=domψ1×⋯×domψK≜Qdomψ=dom _1×·s×dom _K Q. Assumptions for VIPs. Let us make the following assumptions: (A2’) Let 0=(10,⋯,K0)∈Qz^0=(z^0_1,·s,z^0_K)∈ Q be a given point. There exists ∗=(1∗,⋯,K∗)∈Qz^*=(z_1^*,·s,z_K^*)∈ Q in the solution set of the VIP of (V,ψ)(V,ψ), such that for all i∈[K]i∈[K]: i∗∈ℬiz_i^* _i, where ℬi≜i∈ℰi∣∥i0−i∥i≤DiB_i \z_i _i _i^0-z_i _i≤ D_i\ and Di>0D_i>0 is a given distance. (A3’) The operator Vi(j;−j)V_i(z_j;z_-j) is LijL_ij-Lipschitz continuous in j∈domψjz_j _j for any fixed −j∈domψ1×⋯×domψj−1×domψj+1×⋯×domψKz_-j _1×·s×dom _j-1×dom _j+1×·s×dom _K. 111For all =(1,⋯,K)∈Qz=(z_1,·s,z_K)∈ Q, we use the following notations for simplicity: (j;−j)≜(z_j;z_-j) and −j≜(1,⋯,j−1,j+1,⋯,K)z_-j (z_1,·s,z_j-1,z_j+1,·s,z_K). Let the operator family ℱVIP F_VIP be comprised of all the operators with initialization points ((Vi)i∈[K],0)∈VIP ((V_i)_i∈[K],z^0 )∈ P_VIP, such that Assumptions (Â1), (Â3) and (A2’) are satisfied. Notations. To simplify the notations, let us denote ≜(1,⋯,K)∈Qz (z_1,·s,z_K)∈ Q in the context of VIPs. Let us denote L¯ij≜maxLij,Lji,Ai≜Di(∑j∈[K]∖iL¯ijDj), and Bi≜L¯iiDi2, for all i,j∈[K]. L_ij \L_ij,L_ji\,\ A_i D_i ( _j∈[K] \i\ L_ijD_j ), and B_i L_iD_i^2, for all i,j∈[K]. We refer to ∑iAi _iA_i as the cross-coupled conditioning and ∑iBi _iB_i as the diagonal conditioning, and we say that the diagonal conditioning dominates when ∑iBi≫∑iAi _iB_i _iA_i. C.2 Communication and computational costs Distributed oracles. We consider a distributed setting with K agents very similar to the one in Sections 2.2 and 2.1. For all i∈[K]i∈[K]: Agent i controls decision variable i∈domψiz_i _i, has direct access to the function ψi _i, and has access to the oracle i()=Vi() O_i(z)=V_i(z) for ∈Qz∈ Q. We consider per query to i O_i costs ci≥0c_i≥ 0, i∈[K]i∈[K]. Accuracy measure. We consider the following accuracy measure for VIPs: Δ(¯)≜sup∈ℬ∩Q⟨V(),¯−⟩+ψ(¯)−ψ(), for all ∈Q, ( z) _z ∩ Q V(z), z-z +ψ( z)-ψ(z), for all z∈ Q, where ℬ≜ℬ1×⋯×ℬKB _1×·s×B_K. We say that a point ¯∈Q z∈ Q is an ϵε-approximate solution of the VIP if Δ(¯)≤ϵ. ( z)≤ε. Our goal is to find such an ϵε-approximate solution for any ϵ>0ε>0. Problem class. We formally define the overall problem class, denoted by VIP(ℱVIP,(i)i∈[K],ϵ) P_VIP ( F_VIP,( O_i)_i∈[K],ε ), or for short VIP P_VIP. A specific problem instance P∈VIPP∈ P_VIP is constructed by drawing an operator instance (with initial point) V from ℱVIP F_VIP, equipping it with the distributed oracles (i)i∈[K]( O_i)_i∈[K], and specifying a target accuracy ϵ>0ε>0. Solving the instance P requires an algorithm to output an ϵε-approximate solution of V utilizing the distributed oracles. Distributed algorithms for VIPs, communication and oracle costs. To provide the formal definitions of distributed algorithms for the problem class VIP P_VIP, we generalize the information-based framework to K agents. Every query point, message, and output is generated as a deterministic mapping of the information available to the agent at that step. Suppose an algorithm ℳ M proceeds in T rounds. In each round t∈0,…,T−1t∈\0,…,T-1\, each Agent i∈[K]i∈[K] executes multiple local computational steps to generate local query points denoted by it,l=(i,1t,l,…,i,it,l,…,i,Kt,l)for l∈0,…,τit−1,z_i^t,l=(z_i,1^t,l,…,z_i,i^t,l,…,z_i,K^t,l) l∈\0,…, _i^t-1\, where i,it,lz_i,i^t,l is the local variable updated by Agent i, and i,jt,lz_i,j^t,l (for j≠ij≠ i) is the delayed approximation of Agent j’s variable utilized by Agent i. Let Iit,lI_i^t,l denote the accumulated information sequence available to Agent i prior to making its (l+1)(l+1)-th local oracle query in round t. The base case at initialization is Ii0,0=(ψi,0)I_i^0,0=( _i,z^0). During the local computational steps l∈0,…,τit−1l∈\0,…, _i^t-1\, the information sequence of Agent i updates sequentially by appending the newly acquired oracle response for ViV_i: Iit,l+1=(Iit,l,Vi(it,l)).I_i^t,l+1= (I_i^t,l,V_i(z_i^t,l) ). After the local steps in round t, the agents exchange messages. Let Mj→itM_j→ i^t denote the message sent from Agent j to Agent i. The information sequence available to Agent i at the beginning of round t+1t+1 appends the received messages to its prior local history: Iit+1,0=(Iit,τit,Mj→itj≠i).I_i^t+1,0= (I_i^t, _i^t,\M_j→ i^t\_j≠ i ). Definition 7. An algorithm ℳ M is called a distributed algorithm for problem class VIP P_VIP if, when applied to any instance P∈VIPP∈ P_VIP, its execution satisfies the following conditions for all t∈0,…,T−1t∈\0,…,T-1\ and agents i∈[K]i∈[K]: 1. Local Computation: The query points are determined entirely by the locally available information. For all local steps l∈0,…,τit−1l∈\0,…, _i^t-1\, there exists a deterministic mapping function it,lA_i^t,l such that: it,l=it,l(Iit,l).z_i^t,l=A_i^t,l (I_i^t,l ). 2. Communication: The messages exchanged are produced by deterministic mappings of the sender’s local information. For any j≠ij≠ i, there exists a mapping function ℳi→jtM_i→ j^t such that: Mi→jt=ℳi→jt(Iit,τit).M_i→ j^t=M_i→ j^t (I_i^t, _i^t ). 3. Candidate solution: The candidate solutions ¯t+1=(¯1t+1,…,¯Kt+1) z^t+1=( z_1^t+1,…, z_K^t+1) are constructed from the respective agents’ updated information sets. There exists a deterministic mapping function ¯it+1 A_i^t+1 such that: ¯it+1=¯it+1(Iit+1,0). z_i^t+1= A_i^t+1 (I_i^t+1,0 ). For a given instance P∈VIPP∈ P_VIP and a target accuracy ϵ>0ε>0, we define the communication cost required by a distributed algorithm ℳ M on P, denoted by TPℳT M_P, as the smallest integer k∈1,…,Tk∈\1,…,T\ such that the candidate solution ¯k z^k satisfies the target accuracy ϵε. The total number of local oracle queries evaluated by Agent i for instance P up to this point is given by the cumulative number of local steps taken, denoted by Ni,Pℳ=∑r=0TPℳ−1τirN_i,P M= _r=0^T M_P-1 _i^r. Let cic_i denote the computational cost of evaluating a single partial oracle ViV_i. The communication cost and oracle cost of algorithm ℳ M over the entire problem class VIP P_VIP are defined by taking the supremum over all instances: TVIPℳ=supP∈VIPTPℳandNVIPℳ=supP∈VIP(∑i∈[K]ciNi,Pℳ).T M_ P_VIP= _P∈ P_VIPT M_P N_ P_VIP M= _P∈ P_VIP\ ( _i∈[K]c_iN_i,P M ). We first state the classic results of the EG method in Proposition 20, which remains the state-of-the-art communication complexity bound. Proposition 20 (Juditsky et al. 2011, Eq. (6.21)). For any target accuracy ϵ>0ε>0, the communication cost of EG is bounded by (∑i∈[K]Ai+Biϵ),O ( _i∈[K] A_i+B_iε ), and the computational cost of EG is bounded by ((∑i∈[K]ci)(∑i∈[K]Aiϵ)+(∑i∈[K]ci)(∑i∈[K]Biϵ)).O ( ( _i∈[K]c_i ) ( _i∈[K] A_iε )+ ( _i∈[K]c_i ) ( _i∈[K] B_iε ) ). C.3 Decoupled method for variational inequality problems Now, we present our DM-VIP method, which extends the DM-SP into multiplayer general-sum games. Assembled norm. Given parameters αi>0 _i>0 for all i∈[K]i∈[K] (to be specified later), we equip the joint space ℰ=ℰ1×⋯×ℰKE=E_1×…×E_K with the assembled norm: ∥ℰ=⟨,⟩12=∑i=1Kαi∥i∥i2for all ∈ℰ, _E= ,z 12= _i=1^K _i _i _i^2 all $ z $, (11) which corresponds to the block diagonal linear operator =α11⊕⋯⊕αKKP= _1P_1 … _KP_K. Accordingly, we equip the dual space ℰ∗=ℰ1∗×⋯×ℰK∗E^*=E_1^*×…×E_K^* with the corresponding dual norm: ∥ℰ∗=⟨,−1⟩12=∑i=1Kαi−1∥i∥i∗2for all ∈ℰ∗. _E^*= ,P^-1g 12= _i=1^K _i^-1 _i _i^*^2 all $ g ^*$. (12) Template DM-VIP. To extend our decoupled framework to block composite variational inequality problems (VIPs), we first define the coupled conditioning constant, which characterizes the interaction between the K distinct blocks: L¯c≜maxj∈[K][(αjDj)−1∑i∈[K]∖jL¯ij(∑l∈[K]∖iL¯ilDl)αi]. L_ c _j∈[K] [( _jD_j)^-1 _i∈[K] \j\ L_ij ( _l∈[K] \i\ L_ilD_l ) _i ]. (13) Algorithm 5 outlines the Decoupled Method for block composite VIPs (DM-VIP), generalizing the DM-SP procedure. The algorithm maintains a sequence of anchor points tv^t and orchestrates iterative updates among K agents over a distributed network. At the start of each iteration, the agents decouple the joint problem by fixing their remote variables to the current anchor components −itv^t_-i. This allows each Agent i to independently and concurrently solve its regularized local subproblem. Specifically, each agent invokes an internal solver ℳiMRNM^MRN_i to minimize the local residual norm (MRN) up to a target accuracy δit+1=αiλt+12 _i^t+1= _i _t+12 (Line 5). Following this local computation phase, the K agents perform exactly two communication rounds to complete the iteration. In the first round (Line 6), the agents broadcast their locally computed approximate solutions it+1z_i^t+1 to assemble the joint intermediate point t+1z^t+1. In the second round (Line 7), they use this assembled point to evaluate their local partial operators Vi(t+1)V_i(z^t+1), which they subsequently exchange to form the full joint operator Vψ(t+1)V_ψ(z^t+1). Finally, using this assembled joint operator, the agents compute a closed-form reduced step size at+1a_t+1, update the running ergodic average ¯t+1 z^t+1, and perform a joint extragradient-like step to generate the next anchor t+1v^t+1 (Lines 8 and 9). Algorithm 5 DM-VIP(K,(Vi)i∈[K],(ψi)i∈[K],0,(λt)t≥1,(αi)i∈[K]∣(ℳiMRN)i∈[K])DM-VIP (K,(V_i)_i∈[K],( _i)_i∈[K],z^0,( _t)_t≥ 1,( _i)_i∈[K] (M^MRN_i)_i∈[K] ) 1:A local solver ℳiMRNM^MRN_i for the minimization of residual norms. 2:0=(10,…,K0)=0v^0=(v^0_1,…,v^0_K)=z^0. 3:for t=0,1,…,T−1t=0,1,…,T-1 do 4: Let δit+1=αiλt+12δ^t+1_i= _i _t+12 for all i∈[K]i∈[K]. 5: Concurrently for all i∈[K]i∈[K], Agent i computes (it+1,ψi′(it+1))=ℳiMRN(Vi(⋅;−it),ψi+αiλt+12∥⋅−it∥i2,it,δit+1). (z^t+1_i, _i (z^t+1_i) )=M^MRN_i (V_i(·;v^t_-i),\, _i+ _i _t+12 ·-v^t_i _i^2,\,v^t_i,\,δ^t+1_i ). 6: All agents exchange it+1z^t+1_i to assemble the joint point t+1=(1t+1,…,Kt+1)z^t+1=(z^t+1_1,…,z^t+1_K). 7: Agents locally compute Vi(t+1)V_i(z^t+1) and exchange to assemble the full operator: Vψ(t+1)=V(t+1)+(ψ1′(1t+1),…,ψK′(Kt+1)).V_ψ(z^t+1)=V(z^t+1)+ ( _1 (z^t+1_1),…, _K (z^t+1_K) ). 8: Let at+1=2⟨Vψ(t+1),t−t+1⟩∥Vψ(t+1)∥ℰ∗2a_t+1= 2 V_ψ(z^t+1),v^t-z^t+1 V_ψ(z^t+1) _E^*^2 and generate solution ¯t+1=(∑i=1t+1ai)−1∑i=1t+1aii z^t+1= ( _i=1^t+1a_i )^-1 _i=1^t+1a_iz^i. 9: t+1=argmin∈Q[at+1⟨Vψ(t+1),⟩+12∥−t∥ℰ2]v^t+1= *arg\,min_v∈ Q\ [a_t+1 V_ψ(z^t+1),v + 12 -v^t _E^2 ]. 10:end for Similar to the saddle point setting, we refer to Algorithm 5 as a template method because we have abstracted the exact implementation of the inner solvers. For the theoretical guarantees presented below, we merely assume the existence of an algorithmic oracle ℳiMRN(V^i,ψ^i,i,δi)M^MRN_i( V_i, ψ_i,v_i, _i) capable of taking an MRN instance in the space ℰiE_i and successfully returning a sufficiently accurate solution. We defer the detailed implementation of these local solvers to Eq. 15 at the end of this section. For now, let us proceed with the communication complexity of the template DM-VIP method. Theorem 21. Consider the DM-VIP template applied to VIP P_VIP. With the parameter choices of αi=∑j∈[K]∖iL¯ijDjDi _i= _j∈[K] \i\ L_ijD_jD_i for all i∈[K]i∈[K], and λt≡λ≥2L¯c _t≡λ≥ 2 L_ c, we have: TVIPDM-VIP≤2+∑i,j∈[K]i≠j2L¯ijDiDjϵ.T^DM-VIP_ P_VIP≤ 2+ _ subarrayci,j∈[K]\\ i≠ j subarray 2 L_ijD_iD_jε. For the distributed VIP considered in this paper, the classic EG method represents the best known communication cost. Now, let us compare our communication guarantee with that of EG. Remark 6 (Improved communication). The classic EG method takes 1ϵ∑i∈[K](Ai+Bi) 1ε _i∈[K](A_i+B_i) communication rounds (cf. Proposition 20). Our communication cost in Theorem 21 is (∑i∈[K]Ai/ϵ)O ( _i∈[K]A_i/ε ), which is consistently no worse compared to that of EG and is substantially faster when the “diagonal conditioning” dominates—i.e., ∑i∈[K]Bi≫∑i∈[K]Ai. _i∈[K]B_i _i∈[K]A_i. To our knowledge, DM-VIP improves the state-of-the-art communication cost for distributed VIPs. C.4 Detailed proofs C.4.1 Proof for FDS Let us provide the detailed pseudocode of FDS for VIPs in Algorithm 6. Then, we prove the correctness of the solution returned by FDS. Algorithm 6 FDS∥⋅∥ℰ((Vi)i∈[K],(ψi)i∈[K],,λ∣(ℳiMRN)i∈[K])FDS_ · _E ((V_i)_i∈[K],( _i)_i∈[K],v,λ (M^MRN_i)_i∈[K] ) 1:Solver ℳiMRNM^MRN_i for the minimization of residual norms, for all i∈[K]i∈[K]. 2:for i∈[K]i∈[K] do 3: δi=αiλ2 _i= _iλ2. 4: ψ^i=ψi+αiλ2∥⋅−i∥i2 ψ_i= _i+ _iλ2 ·-v_i _i^2. 5: (i+,ψi′(i+))=ℳiMRN(Vi(⋅;−i),ψ^i,i,δi)(z_i^+,ψ _i(z_i^+))=M^MRN_i (V_i(·;v_-i), ψ_i,v_i, _i ). 6: ψi′(i+)=ψ^i′(i+)−αiλi(i+−i) _i (z_i^+)= ψ_i (z_i^+)- _i _i(z_i^+-v_i). 7:end for 8:return (+,ψ′(+))(z^+,ψ (z^+)), where +=(1+,⋯,K+)z^+=(z^+_1,·s,z^+_K) and ψ′(+)=(ψ1′(1+),⋯,ψK′(K+))ψ (z^+)=(ψ _1(z_1^+),·s,ψ _K(z_K^+)). Lemma 22. Under (A3’), for λ≥2L¯cλ≥ 2 L_ c, FDS (Algorithm 6) returns the correct solution of the MS subproblem given by (V,ψ,,λ)(V,ψ,v,λ). Proof of Lemma 22. For all i∈[K]i∈[K], by (A3’) and then by the relative distance accuracy, we have ∥Vi(+)+ψi′(i+)+αiλi(i+−i)∥i∗ V_i(z^+)+ _i (z^+_i)+ _i _i(z^+_i-v_i) _i^* (14) ≤∥Vi(i+;−i)+ψi′(+)+αiλi(i+−i)∥i∗+∑j∈[K]∖iLij∥j+−j∥j ≤ V_i(z^+_i;v_-i)+ _i (z^+)+ _i _i(z^+_i-v_i) _i^*+ _j∈[K] \i\L_ij ^+_j-v_j _j ≤δi∥i+−i∥i+∑j∈[K]∖iLij∥j+−j∥j. ≤ _i ^+_i-v_i _i+ _j∈[K] \i\L_ij ^+_j-v_j _j. Finally, we assemble the norms: ∥V(+)+ψ′(+)+λ(+−)∥ℰ∗2 V(z^+)+ψ (z^+)+ (z^+-v) _E^*^2 =∑i∈[K]αi−1∥Vi(+)+ψi′(i+)+αiλi(i+−i)∥i∗2 = _i∈[K] _i^-1 V_i(z^+)+ _i (z^+_i)+ _i _i(z^+_i-v_i) _i^*^2 ≤(14)∑i∈[K]αi−1(δi∥i+−i∥i+∑j∈[K]∖iLij∥j+−j∥j)2 ( eq:FDS-bound-on-block-coordinate-i)≤ _i∈[K] _i^-1 ( _i ^+_i-v_i _i+ _j∈[K] \i\L_ij ^+_j-v_j _j )^2 ≤∑i∈[K]2αi−1[δi2∥i+−i∥i2+(∑j∈[K]∖iLij∥j+−j∥j)2] ≤ _i∈[K]2 _i^-1 [ _i^2 ^+_i-v_i _i^2+ ( _j∈[K] \i\L_ij ^+_j-v_j _j )^2 ] =λ22∑i∈[K](αi∥i+−i∥i2)+2∑i∈[K][αi−1(∑j∈[K]∖iLij∥j+−j∥j)2] = λ^22 _i∈[K] ( _i ^+_i-v_i _i^2 )+2 _i∈[K] [ _i^-1 ( _j∈[K] \i\L_ij ^+_j-v_j _j )^2 ] ≤λ22∥+−∥ℰ2+2∑i∈[K][αi−1(∑l∈[K]∖iLilDl)(∑j∈[K]∖iLijDj∥j+−j∥j2)] ≤ λ^22 ^+-v _E^2+2 _i∈[K] [ _i^-1 ( _l∈[K] \i\L_ilD_l ) ( _j∈[K] \i\ L_ijD_j ^+_j-v_j _j^2 ) ] =λ22∥+−∥ℰ2+2∑j∈[K][∥j+−j∥j2Dj(∑i∈[K]∖jLij(∑l∈[K]∖iLilDl)αi)] = λ^22 ^+-v _E^2+2 _j∈[K] [ z^+_j-v_j _j^2D_j ( _i∈[K] \j\ L_ij ( _l∈[K] \i\L_ilD_l ) _i ) ] ≤(13)λ22∥+−∥ℰ2+2∥+−∥ℰ2L¯c2 ( eq:coupled-conditioning)≤ λ^22 ^+-v _E^2+2 ^+-v _E^2 L_ c^2 ≤λ22∥+−∥ℰ2+λ22∥+−∥ℰ2=λ2∥+−∥ℰ2. ≤ λ^22 ^+-v _E^2+ λ^22 ^+-v _E^2=λ^2 ^+-v _E^2. ∎ C.4.2 Proof for MRN Our algorithm is built upon Lemma 23, the proof of which can be found in [Boţ and Chenchene, 2024, Corollary 2.4]. Lemma 23. Assume (Â1), (Â3), and that the solution set of the VIP of (Vw,ψw)(V_w, _w) is non-empty. Then, there exists an algorithm, denoted by (+,ψw′(+))=FEGM(Vw,ψw,w,ξ∣L),(w^+, _w (w^+))=FEGM\ (V_w, _w,v_w,ξ L), which takes no more than C0⋅LξC_0· Lξ operator queries and returns (+,ψw′(+))(w^+, _w (w^+)) that satisfies ξ-distance-to-solution accuracy, where C0>0C_0>0 is some fixed constant. C.4.3 Concrete implementation We are now back to considering the VIPs. Let us use FEGM in Lemma 23 for the minimization of residual norms: ℳiFEGM(V^i,ψ^i,i,δi)≜FEGM(V^i,ψ^i,i,2δi3∣Lii).M^FEGM_i( V_i, ψ_i,v_i, _i) ( V_i, ψ_i,v_i, 2 _i3 L_i). Then, for any Monteiro-Svaiter Subproblem given by (V,ψ,,λ)(V,ψ,v,λ), we leverage the solver FDS-FEGM(V,ψ,,λ)=FDS∥⋅∥ℰ(V,ψ,,λ∣(ℳiFEGM)i∈[K]).FDS-FEGM(V,ψ,v,λ)=FDS_ · _E (V,ψ,v,λ (M^FEGM_i)_i∈[K] ). Finally, we obtain the concrete algorithm DM-VIP as follows: ROM∥⋅∥ℰ((Vi)i∈[K],(ψi)i∈[K],0,(λt)t≥1∣FDS-FEGM). _ · _E ((V_i)_i∈[K],( _i)_i∈[K],z^0,( _t)_t≥ 1 -FEGM ). (15) Combining Lemmas 2, 22, 4 and 23, with the implementation in Eq. 15, we conclude that Lemma 24 holds for the constant C0C_0 from Lemma 23. We include the complete proof below. Lemma 24. Under (A2’) and (A3’), for λt+1≡λ≥2L¯c, _t+1≡λ≥ 2 L_ c, DM-VIP (Algorithm 5) with the implementation in Eq. 15 takes no more than 2T2T communication rounds and no more than T⋅(1+C0⋅3Liiαiλ)T· (1+C_0· 3L_i _iλ ) queries to ViV_i, for all i∈[K]i∈[K], and obtains an ϵε-approximate solution ¯T, z^T, where T=⌈∑iαiλDi22ϵ⌉T= _i _iλ D_i^22ε and C0>0C_0>0 is some fixed constant. Proof of Lemma 24. By (A1’), we have Δ(¯T)≤(∑t=0T−1at+1)−1max∈ℬ∩Q[∑t=0T−1at+1⟨Vψ(t+1),t+1−⟩]. ( z^T)≤ ( _t=0^T-1a_t+1 )^-1 _z ∩ Q\ [ _t=0^T-1a_t+1 V_ψ(z^t+1),z^t+1-z ]. Further, with λ≥2L¯cλ≥ 2 L_ c, by Lemmas 2 and 22, we have Δ(¯T)≤(∑t=0T−1at+1)−1max∈ℬ∩Q[∑t=0T−1at+1⟨Vψ(t+1),t+1−⟩] ( z^T)≤ ( _t=0^T-1a_t+1 )^-1 _z ∩ Q\ [ _t=0^T-1a_t+1 V_ψ(z^t+1),z^t+1-z ] ≤(∑t=0T−1at+1)−1[∑i∈[K](αi2maxi∈ℬi∩domψi∥i0−i∥i2)] ≤ ( _t=0^T-1a_t+1 )^-1 [ _i∈[K] ( _i2 _z_i _i _i\ ^0_i-z_i _i^2 ) ] ≤(∑t=0T−11λt+1)−1⋅12∑i∈[K]αiDi2≤ϵ, ≤ ( _t=0^T-1 1 _t+1 )^-1· 12 _i∈[K] _iD_i^2≤ε, where the last inequality follows from the assignments of (λt)t≥1( _t)_t≥ 1 and T. Therefore, the number of communication rounds is bounded by 2T2T. Now we count the number of gradient queries. By Lemma 4, FEGM always returns the solution with the required relative distance accuracy; and in view of Lemma 23, it takes no more than C0⋅3LiiαiλC_0· 3L_i _iλ gradient queries to ViV_i, for all i∈[K]i∈[K]. Therefore, the numbers of queries to ViV_i are bounded by T⋅(1+C0⋅3Liiαiλ)T· (1+C_0· 3L_i _iλ ), for all i∈[K]i∈[K]. ∎ Remark 7 (Oracle comparison). Under the same choice of parameters as in Theorem 21, the oracle cost of DM-VIP is bounded by 2ϵ(∑i∈[K]ci)(∑i∈[K]Ai)+3C0ϵ(∑i∈[K]BiciAi)(∑i∈[K]Ai). 2ε ( _i∈[K]c_i ) ( _i∈[K]A_i )+ 3C_0ε ( _i∈[K] B_ic_iA_i ) ( _i∈[K]A_i ). (16) Compared to the computational cost of EG, which is given by 1ϵ(∑i∈[K]ci)(∑i∈[K]Ai)+1ϵ(∑i∈[K]ci)(∑i∈[K]Bi), 1ε ( _i∈[K]c_i ) ( _i∈[K]A_i )+ 1ε ( _i∈[K]c_i ) ( _i∈[K]B_i ), our computational cost in Eq. 16 differs in the second term. Consequently, our DM-VIP may offer an advantage or disadvantage depending on the relative conditioning of AiA_i, BiB_i, and cic_i for i∈[K]i∈[K].