Conceptio › Archive › arXiv CS

arXiv CSopen access

Efficient Gradient Methods for Distributed Saddle Problems

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

Open Source ↗Direct PDF ↓

cloud distributed-computing parallel-computing

distributed computing, parallel computing, cloud

Efficient Gradient Methods for Distributed Saddle Problems∗ Ruichen Luo†

Anton Rodomanov‡

Sebastian U. Stich§

arXiv:2605.18604v1 [math.OC] 18 May 2026

May 15, 2026

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. Keywords: convex optimization, distributed optimization, complexity bounds, saddle problems, variational inequalities, extra-gradient method, fast gradient method.

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 ∗

The 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. † Institute of Science and Technology Austria. E-mail: [email protected]. ‡ CISPA Helmholtz Center. E-mail: [email protected]. § CISPA Helmholtz Center. E-mail: [email protected].

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 minx maxy F (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 stateof-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 gradientspan 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.1: Summary of algorithms and complexity results for distributed SPs. Method M

Communication Cost

EGc

Lxy Dx Dy + ϵ ϵ O log 1ϵ (weakly coupled only) q 2 +L D̂ 2 Lx D̂x Lmax D̂x D̂y y y log2 1ϵ O + ϵ ϵ L D̂ D̂ O xy ϵ x y log3 1ϵ O θ

DGDA

Cat-EG

Cat-Cat-DAGDA

Lower Bound (Thm. 15) Ω DM-SP (Thm. 1)

2 2 +Ly Dy Lx Dx

Lxy Dx Dy ϵ L D D θ xy ϵ x y

Better Oracle?a

Multi-Agentb

Yes

Maybe

–

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 P ◦ . SP

Notations. Let [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 Ex = Rnx with the norm ∥x∥x = ⟨Px x, x⟩1/2 , where Px : Ex → Ex∗ is a self-adjoint positive definite operator and the dual pairing ⟨ϕx , x⟩ denotes ϕx (x). We denote its corresponding dual norm by ∥·∥x∗ . We assume analogous geometries for Ey = Rny , Ei = Rni (i ∈ [K]), and Ew = Rnw , associated with their respective operators Py , Pi , and Pw . For a function ψ : E → R ∪ {+∞}, let dom ψ denote its effective domain and ∂ψ(z) its subdifferential at z ∈ dom ψ. Finally, for any set of vectors S, let span S denote its linear span.

Saddle problems with distributed oracles

In the context of saddle problems, we consider two separate finite-dimensional real vector spaces, Ex = (Rnx , ∥·∥x ) and Ey = (Rny , ∥·∥y ). We are interested in solving composite Saddle Problems (SPs) of the following form: min max F (x, y) ≜ f (x, y) + ψx (x) − ψy (y) , (2.1) x∈dom ψx y∈dom ψy

where ψx : Ex → R ∪ {+∞} and ψy : Ey → R ∪ {+∞} represent relatively simple local components (such as regularizers or indicator functions for constrained sets), and f (·, ·) is a real-valued coupling function defined on an open set containing the domain Q ≜ dom ψx ×dom ψy . To simplify the notation, we denote the joint variable by z ≜ (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 P(F , Ox , Oy , ∆, ϵ), is a family of problems sharing • a function family F from which each problem instance is drawn; • two distributed oracles Ox and Oy , providing first-order information about an instance; • an accuracy measure ∆, quantifying the quality of a candidate solution; and • a target accuracy ϵ > 0. A specific problem instance P ∈ P is determined by a function instance F ∈ F , the corresponding oracles (Oxf , Oyf ), 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 (x ∈ dom ψx and y ∈ dom ψ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 (for example, the ability to evaluate its proximal-point mapping) and to a problem-dependent first-order oracle Ox whose concrete form is specified by the problem class; Agent y has the analogous access to ψy and Oy . 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 (2.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 Ox and Oy ; 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 − 1}, Agent x successively queries Ox at t,l t t τxt points zt,l x for l = 0, . . . , τx − 1; symmetrically, Agent y queries Oy at τy points zy . 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,l denote the oracle responses collected by Agent x prior to its (l+1)-th query in round t. The information set is initialized empty, Ix0,0 = ∅, and each local query appends the corresponding t,τ t oracle response, so that Ix x summarizes the round; Agent y is symmetric. Each agent’s local data, namely the initial point z0 and its component ψx (or ψ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: t t,τ t Ixt+1,0 = Iyt+1,0 = Ixt,τx , Iy y . (2.2) At the conclusion of round t, the algorithm produces a candidate solution z̄t+1 = (x̄t+1 , ȳt+1 ) from the merged information set Ixt+1,0 = Iyt+1,0 . Definition 1. An algorithm M is called a distributed method if, for every round t ∈ {0, . . . , T − 1}, the following hold: t 1. Each local query point zt,l x , l ∈ {0, . . . , τx − 1}, is a deterministic function of Agent x’s t,l t current information set Ix ; symmetrically, each zt,l y , l ∈ {0, . . . , τy − 1}, is a determint,l istic function of Iy . 2. The candidate solution z̄t+1 is a deterministic function of the merged information set Ixt+1,0 = Iyt+1,0 defined in (2.2). While these deterministic mappings could in principle be randomized, we restrict attention to the deterministic case here for simplicity. 5

Communication and oracle complexities. For a given problem instance P and target accuracy ϵ > 0, the communication complexity of M on P , denoted TPM , is the smallest integer k ∈ {1, . . . , T } such that the candidate solution 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 TPM −1 TPM −1 X X t M M τyt . τx and Ny,P = Nx,P = t=0

t=0

The communication complexity and oracle complexity of M over a problem class P are defined by taking the supremum over all instances: M M M M = sup cx Nx,P + cy Ny,P , TP = sup TPM and NP P ∈P

P ∈P

where cx and cy are fixed constants reflecting the computational costs of evaluating a single query to Ox and Oy , 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 PSP

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 y ∈ dom ψy , the function f (·, y) is convex; and for any fixed x ∈ dom ψx , the function f (x, ·) is concave. The functions ψx and ψy are proper, closed, and convex. (A2) Let Dx , Dy > 0 be distance parameters, and let z0 = (x0 , y0 ) ∈ Q be a given initial point. Relative to this initialization, Problem (2.1) has a saddle point (x∗ , y∗ ) ∈ Q such that ∥x0 − x∗ ∥x ≤ Dx and ∥y0 − y∗ ∥y ≤ Dy . (A3) The function f (·, ·) is continuously differentiable over Q. Moreover, its gradients are Lipschitz continuous. That is, with Lipschitz parameters Lx , Lxy , Ly > 0, for all x, x′ ∈ dom ψx and y, y′ ∈ dom ψy , we have: ∥∇x f (x′ , y′ ) − ∇x f (x, y)∥x∗ ≤ Lx ∥x′ − x∥x + Lxy ∥y′ − y∥y , ∥∇y f (x′ , y′ ) − ∇y f (x, y)∥y∗ ≤ Lxy ∥x′ − x∥x + Ly ∥y′ − y∥y . Let F denote the function family consisting of all instances F = (f, ψx , ψy , z0 ) that satisfy (A1) to (A3) for a fixed set of parameters (Lx , Lxy , Ly , Dx , Dy ). The non-composite case ψx ≡ ψy ≡ 0 is included as a special instance.

Distributed partial-gradient oracle. For SPs, the abstract oracles Ox and Oy are concretely realized as (deterministic) partial-gradient oracles: for a given function instance F ∈ F with coupling function f and any input point z ∈ Q, • Agent x queries Ox , which returns Oxf (z) = ∇x f (z). • Agent y queries Oy , which returns Oyf (z) = ∇y f (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 (z) = g(z) + fx (x) − fy (y), where g(z) is a coupled global utility, while fx (x) and fy (y) are private utilities accessible only to Agents x and y, respectively. The partial-gradient oracles then take the form Oxf (z) = ∇x g(z) + ∇fx (x)

Oyf (z) = ∇y g(z) − ∇fy (y),

for all z ∈ Q. Due to the distributed setting, Agent x is entirely blind to the private utility fy and can only execute Ox , and vice versa. Accuracy measure. To evaluate the quality of a candidate solution (x̄, ȳ) ∈ Q, we rely on the restricted duality gap. Let Bx ≜ {x ∈ Ex | ∥x0 − x∥x ≤ Dx } and By ≜ {y ∈ Ey | ∥y0 − y∥y ≤ Dy } denote the balls of radii Dx , Dy around the initial points. Over the bounded domain B ≜ Bx × By , the duality gap is defined as ∆(x̄, ȳ) ≜ max F (x̄, y) − F (x, ȳ) . (x,y)∈B∩Q

We say that a pair (x̄, ȳ) ∈ Q is an ϵ-saddle point of Problem (2.1) if ∆(x̄, ȳ) ≤ ϵ. Our goal is to design an algorithm that produces such an ϵ-saddle point for a given ϵ > 0. We remark that for the classic problem of constrained optimization with bounded domains, one can enclose the constrained sets in the balls Bx and By with sufficiently large radius (e.g., the diameter of the constrained sets); the restricted saddle problem in form (2.1) then coincides with the original one. Problem class PSP . Assembling the function family F , the partial-gradient oracles (Ox , Oy ), the duality-gap accuracy measure ∆, and a target accuracy ϵ > 0, we obtain the problem class of interest, denoted by PSP (F , Ox , Oy , ∆, ϵ), or for short PSP . Solving an instance P ∈ PSP requires an algorithm to output an ϵ-saddle point of F utilizing the distributed oracles. To facilitate later discussion, we refer to the terms Lx Dx2 and Ly Dy2 as the diagonal conditioning, and the term Lxy Dx Dy 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. 7

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 PSP , 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; (ii) it is highly sensitive to the inexactness of the diameter estimates D̂x and D̂y ; (iii) 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.1, the classic EG method remains a formidable baseline for PSP , and improving its communication and oracle complexity remains a significant challenge.

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 crosscoupled 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 (to be specified later), we equip the joint space E = Ex × Ey with the assembled norm: q 1 ∥z∥E = ⟨Pz, z⟩ 2 = αx ∥x∥2x + αy ∥y∥2y for all z ∈ E, (3.1) which corresponds to the block diagonal linear operator P = αx Px ⊕ αy Py . Template DM-SP. Algorithm 3.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 vt and proceeds iteratively. First, the agents decouple the joint problem by fixing the remote variable at the current anchor vt . This allows Agent x and Agent y to independently and concurrently solve their respective regularized local subproblems up to target accuracies δxt+1 and δyt+1 (Lines 4 and 5). Specifically, Agent x aims to approximately compute arg minx∈dom ψx [f (x, vyt ) + αx λt+1 ∥x − vxt ∥2x + ψx (x)] by finding a point xt+1 whose regularized subgradient norm 2 satisfies the exact mathematical bound specified in Line 4. Agent y symmetrically perα λ forms an approximate minimization for its corresponding objective −f (vxt , y)+ y 2t+1 ∥y− 2 t vy ∥y + ψ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 zt+1 = (xt+1 , yt+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ψ (zt+1 ). Finally, using this assembled operator, the agents compute a closed-form step size at+1 , update the running ergodic average z̄t+1 , and perform an extragradient-like step to generate the next anchor vt+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. We refer to Algorithm 3.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 PSP , assuming its local trajectories L Dy Lxy Dx satisfy Definition 1. With the parameter choices of αx = xy Dx , αy = Dy , and λt ≡ λ = 2, we have: Lxy Dx Dy DM-SP TP ≤2+4 . SP ϵ Remark 1 (Communication improvement). Theorem 1 shows that the communication cost of DM-SP depends only on the cross-coupled conditioning Lxy Dx Dy , 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 9

Algorithm 3.1 DM-SP(f, (ψx , ψy ), z0 , (λt )t≥1 , (αx , αy )) 0 , v0 ) = z0 . 1: v0 = (vx y

2: for t = 0, 1, . . . , T − 1 do 3: 4:

α λ

Let δxt+1 = αx λ2t+1 and δyt+1 = y 2t+1 . Agent x finds xt+1 and ψx′ (xt+1 ) ∈ ∂ψx (xt+1 ) such that ∥∇x f (xt+1 , vyt ) + αx λt+1 (xt+1 − vxt ) + ψx′ (xt+1 )∥x∗ ≤ δxt+1 ∥xt+1 − vxt ∥x .

Agent y finds yt+1 and ψy′ (yt+1 ) ∈ ∂ψy (yt+1 ) such that ∥−∇y f (vxt , yt+1 ) + αy λt+1 (yt+1 − vyt ) + ψy′ (yt+1 )∥y∗ ≤ δyt+1 ∥yt+1 − vyt ∥y .

6: 7:

8: 9:

Exchange xt+1 and yt+1 to assemble zt+1 = (xt+1 , yt+1 ). Calculate corresponding coordinates of Vψ (zt+1 ), then exchange to assemble: Vψ (zt+1 ) = ∇x f (zt+1 ) + ψx′ (xt+1 ), −∇y f (zt+1 ) + ψy′ (yt+1 ) . 2⟨Vψ (zt+1 ),vt −zt+1 ⟩ and generate solution z̄t+1 = ∥Vψ (zt+1 )∥2E ∗ vt+1 = arg minv∈Q at+1 ⟨Vψ (zt+1 ), v⟩ + 12 ∥v − vt ∥2E .

Let at+1 =

Pt+1

i=1 ai

−1 Pt+1

i i=1 ai z .

10: end for

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 L D D later in Section 6, our O xy ϵ x 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 Dx and Dy in advance, but it has access to upper estimates D̂x ≥ Dx and D̂y ≥ Dy . Let θ ≜

Dx D̂y

Dy D̂x

, D̂x Dy D̂y Dx which quantifies the disproportionality between the true distance parameters and their estimates. Note that θ ≥ 2 with equality if and only if the estimates are proportional, i.e., D̂x /Dx = D̂y /Dy . Now, consider the DM-SP template applied to PSP . With the parameter choices of αx =

Lxy D̂y L D̂ , αy = xy x , and λt ≡ λ = 2, we have: D̂x D̂y DM-SP TP ≤ 2 + 2θ SP

Lxy Dx Dy . ϵ

In particular, θ provides a scale-invariant robustness compared to existing accelerated frameworks. As shown in Table 1.1, the communication complexities of Cat-EG and CatCat-DAGDA scale directly with the product of the estimates, D̂x D̂y . Consequently, if both agents conservatively overestimate their domain sizes by a uniform factor c ≫ 1 (i.e., D̂x = cDx and D̂y = cDy ), the communication cost of Catalyst-based methods inflates by a massive factor of c2 . For DM-SP, however, this uniform overestimation perfectly cancels out, yielding θ = 2. 10

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 coordinatewise 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 ψ → E ∗ and any function ψ(·) : E → R ∪ {+∞}, we say that z∗ ∈ E is a (strong) solution of the VIP of (V, ψ) if ⟨V (z∗ ), z − z∗ ⟩ + ψ(z) ≥ ψ(z∗ ), for all z ∈ dom ψ.

(4.1)

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 ψ: that is, ⟨V (z′ ) − V (z), z′ − z⟩ ≥ 0, for all z′ , z ∈ dom ψ. There is another notion of a weak solution with the alternative formulation: ⟨V (z), z − z∗ ⟩ + ψ(z) ≥ ψ(z∗ ), for all z ∈ dom ψ. 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∗ ∈ E is the solution of (4.1) if and only if 0 ∈ V (z∗ ) + ∂ψ(z∗ ). Indeed, associated with the saddle problem given by (f, ψx , ψy ), let us consider the VIP given by (V f , ψz ), where V f (z) ≡ (∇x f (z), −∇y f (z)) and ψz (z) ≡ ψx (x) + ψy (y), for all z ∈ Q.

(4.2)

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. 11

Monteiro-Svaiter Subproblem. Given an operator V : dom ψ → E ∗ , a function ψ : E → R ∪ {+∞}, a reference point v ∈ dom ψ, and a real number λ > 0, we say (z+ , ψ ′ (z+ )) is a solution of the MSS if z+ ∈ dom ψ, ψ ′ (z+ ) ∈ ∂ψ(z+ ), and ∥V (z+ ) + ψ ′ (z+ ) + λP(z+ − v)∥E ∗ ≤ λ∥z+ − v∥E .

(4.3)

We will discuss how to solve the MSSs later in Section 4.3. But for now, let us assume there exists a solver MMS (V, ψ, v, λ) for the MSSs, which takes an MSS given by (V, ψ, v, λ) and returns a solution of it. Built upon such a solver MMS , we now introduce ROM in Algorithm 4.1. At each iteration t: the solver MMS returns a solution (zt+1 , ψ ′ (zt+1 )) for the MSS built at reference point vt ; this solution is used as a midpoint to compute subgradient Vψ (zt+1 ); then the ‘extragradient-type’ step is taken with the stepsize at+1 . Algorithm 4.1 ROM∥·∥E (V, ψ, z0 , (λt )t≥1 | MMS ) Require: A solver MMS for the MSSs. 1: v0 = z0 . 2: for t = 0, 1, · · · do 3: (zt+1 , ψ ′ (zt+1 )) = MMS (V, ψ, vt , λt+1 ). 4: Vψ (zt+1 ) = V (zt+1 ) + ψ ′ (zt+1 ). 5:

at+1 =

2⟨Vψ (zt+1 ),vt −zt+1 ⟩ . ∥Vψ (zt+1 )∥2E ∗

vt+1 = arg minv∈dom ψ at+1 ⟨Vψ (zt+1 ), v⟩ + 12 ∥v − vt ∥2E . 7: end for

Next, let us show the convergence of ROM. Lemma 2. ROM (Algorithm 4.1) ensures for all v ∈ dom ψ and for all T ≥ 1, T −1 X t=0

1 1 at+1 ⟨Vψ (zt+1 ), zt+1 − v⟩ ≤ ∥v0 − v∥2E − ∥vT − v∥2E . 2 2

1 Moreover, we have at+1 ≥ λt+1 , for all t ≥ 0.

Proof. By the optimality of vt+1 , we have for all v ∈ dom ψ, 1 1 1 at+1 ⟨Vψ (zt+1 ), v − vt+1 ⟩ + ∥vt − v∥2E ≥ ∥vt+1 − v∥2E + ∥vt − vt+1 ∥2E , 2 2 2

and therefore, 1 at+1 ⟨Vψ (zt+1 ), v − zt+1 ⟩ + ∥vt − v∥2E 2 1 1 ≥ at+1 ⟨Vψ (zt+1 ), vt+1 − zt+1 ⟩ + ∥vt+1 − v∥2E + ∥vt − vt+1 ∥2E 2 2 1 1 t+1 2 t+1 t t+1 − v∥E + at+1 ⟨Vψ (zt+1 ), vt+1 − vt ⟩ + ∥vt − vt+1 ∥2E = at+1 ⟨Vψ (z ), v − z ⟩ + ∥v 2 2 2 a 1 ≥ at+1 ⟨Vψ (zt+1 ), vt − zt+1 ⟩ + ∥vt+1 − v∥2E − t+1 ∥Vψ (zt+1 )∥2E ∗ 2 2 1 t+1 ≥ ∥v − v∥2E , 2 where the last inequality follows from the definition of at+1 in Line 5 of Algorithm 4.1. Then, the desired bound follows from summing the above inequality over t from 0 to T − 1. Next, we show the lower bound for at . For all t ≥ 1, we have ⟨Vψ (zt ), vt−1 − zt ⟩ −

1 ∥Vψ (zt )∥2E ∗ 2λt

λt t 1 ∥z − vt−1 ∥2E − ∥Vψ (zt ) + λt P(zt − vt−1 )∥2E ∗ 2 2λt ≥ 0, ≡

where the last inequality follows from Eq. (4.3). Therefore, we have at =

4.3

2⟨Vψ (zt ), vt − zt ⟩ 1 ≥ . λt ∥Vψ (zt )∥2E ∗

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 (V f , ψz , v, λ), (4.4) where V f and ψz are defined in Eq. (4.2), the reference point v = (vx , vy ) ∈ dom ψx × dom ψy , and the assembled norm ∥·∥E is associated with parameters αx and αy . 2Lxy We say that an MSS has a weak coupling if λ ≥ 2L̄c ≜ √ . In this section, we αx αy 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 ) 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 → Ew∗ , a function ψw : Ew → R ∪ {+∞}, a reference point vw ∈ dom ψw , and an accuracy δ > 0, we ′ (w+ )) minimizes the residual norm to δ-relative distance accuracy, if say (w+ , ψw + ′ (w+ ) ∈ ∂ψ (w+ ), and w ∈ dom ψw , ψw w ′ ∥Vw (w+ ) + ψw (w+ )∥w∗ ≤ δ∥w+ − vw ∥w .

Let us, again, defer the discussion of solving MRNs to Section 4.4. But for now, let us assume there exist solvers MMRN (Vx , ψ̂x , vx , δx ) x

and

MMRN (Vy , ψ̂y , vy , δy ) y

for the coordinate-wise MRNs in spaces Ex and Ey . 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 4.2), 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 4.2 FDS∥·∥E ((∇x f, −∇y f ), (ψx , ψy ), v, λ | (MMRN , MMRN )) x y Require: Solvers MMRN and MMRN for the coordinate-wise MRNs. x y αy λ αx λ 2 1: ψ̂x = ψx + 2 ∥· − vx ∥x and ψ̂y = ψy + 2 ∥· − vy ∥2y . 2: Agent x and Agent y respectively compute (x+ , ψ̂x′ (x+ )) = MMRN ∇x f (·, vy ), ψ̂x , vx , δx and x (y+ , ψ̂y′ (y+ )) = MMRN −∇y f (vx , ·), ψ̂y , vy , δy , y α λ

where δx = α2x λ and δy = 2y . 3: ψx′ (x+ ) = ψ̂x′ (x+ ) − αx λPx (x+ − vx ) and ψy′ (y+ ) = ψ̂y′ (y+ ) − αy λPy (y+ − vy ). 4: return (z+ , ψ ′ (z+ )), where z+ = (x+ , y+ ) and ψ ′ (z+ ) = (ψx′ (x+ ), ψy′ (y+ )).

Lemma 3. Consider the saddle problem given by (f, ψx , ψy ) which satisfies (A3). For λ ≥ 2L̄c , FDS (Algorithm 4.2) returns a solution of the MSS given in Eq. (4.4). 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: (Â1) The function ψw is a (simple) proper closed convex function. The operator Vw is monotone over dom ψw . 14

(Â2) The set-valued operator Vw + ∂ψw is µ-strongly maximally monotone over dom ψw . That is, ′ ′ Vw (w′ ) + ψw (w′ ) − Vw (w) − ψw (w), w′ − w ≥ µ∥w′ − w∥2w , ′ (w′ ) ∈ ∂ψ (w′ ), and ψ ′ (w) ∈ ∂ψ (w), for all w′ , w ∈ dom ψw , ψw w w w (Â3) The operator Vw (w) is L-Lipschitz continuous over w ∈ dom ψw . (Â4) The operator Vw = ∇fw , where fw is a continuously differentiable function defined on an open set containing dom ψw . The theoretical guarantee provided in the literature is usually based on the distanceto-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. ′ (w+ )) satisfies ξDefinition 2 (Distance-to-solution accuracy). We say that (w+ , ψw ′ (w+ ) ∈ ∂ψ (w+ ), and ∥V (w+ ) + ψ ′ (w+ )∥ ∗ ≤ distance-to-solution accuracy if w+ ∈ dom ψw , ψw w w w w ξ∥vw − w̃∥w for some w̃ in the solution set of the VIP of (Vw , ψw ).

Lemma 4. Consider the MRN problem given by (Vw , ψw , vw , δ) which satisfies (Â2). Let µδ ′ (w+ )) satisfies ξ-distance-to-solution accuracy, then (w+ , ψ ′ (w+ )) ξ ≤ µ+δ . If (w+ , ψw w is a solution of the MRN problem. Proof. In view of the triangle inequality and the µ-strong maximal monotonicity of Vw + ψw , we have ∥vw − w̃∥w ≤ ∥w+ − vw ∥w + ∥w+ − w̃∥w 1 ′ ≤ ∥w+ − vw ∥w + ∥Vw (w+ ) + ψw (w+ )∥w∗ µ ξ ≤ ∥w+ − vw ∥w + ∥vw − w̃∥w . µ Then, we have ∥vw − w̃∥w ≤

µ ∥w+ − vw ∥w . µ−ξ

Therefore, we have ′ ∥Vw (w+ ) + ψw (w+ )∥w∗ ≤ ξ∥vw − w̃∥w ≤

µξ ∥w+ − vw ∥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 2 in Algorithm 4.2, 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 B.1 in Appendix B for completeness. Let us denote this algorithm as ARM (∇fw , ψw , vw , ξ | L), 15

(4.5)

which takes an MRN instance of interest, has knowledge of the parameter L in (Â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 (Â1), (Â3), (Â4), and that the solution set of the VIP of (Vw , ψw ) is non-empty. Let ′ (w+ , ψw (w+ )) = ARM (∇fw , ψw , vw , ξ | L). q + ′ + Then, ARM takes no more than 34· 3L 2ξ queries to ∇fw (·) and ensures that (w , ψw (w )) satisfies ξ-distance-to-solution accuracy.

Decoupled method for SPs: Concrete implementation

We are now back to considering the original SPs in Eq. (2.1). Let us combine the technical components in Section 4 and present the final, implementable version of DM-SP. Implementation of DM-SP. of Residual Norms:

We use the ARM solver in Eq. (4.5) for Minimization

2δx | Lx , 3 2δ y | Ly . MMRN (Vy , ψ̂y , vy , δy ) ≜ ARM Vy , ψ̂y , vy , y 3

MMRN (Vx , ψ̂x , vx , δx ) ≜ ARM Vx , ψ̂x , vx , x

(5.1)

Consider the assembled norm ∥·∥E with parameters αx and αy . Then, for any MonteiroSvaiter Subproblem given by (V f , ψz , v, λ), we leverage the solver FDS-ARM(V f , ψz , v, λ) = FDS∥·∥E V f , ψz , v, λ | (MMRN , MMRN ) x y Finally, we obtain the concrete algorithm DM-SP as follows: ROM∥·∥E V f , ψz , z0 , (λt )t≥1 | FDS-ARM .

(5.2)

Combining Lemmas 2 to 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. (5.2), applied to PSP . 2L Under (A1) to (A3), for λt+1 ≡ λ ≥ √αxxyαy , the algorithm takes no more than 2T communication rounds, no more than r 9Lx T · 1 + 34 2αx λ

queries to ∇x f , and no more than s T · 1 + 34

9Ly 2αy λ

queries to −∇y f , and obtains an ϵ-saddle point z̄T , where l αx λD2 + αy λD2 m x y T = . 2ϵ Proof. By (A1), we have T

∆(z̄ ) ≤

−1 TX

at+1

−1

TX −1

max

z∈B∩Q

t=0

t+1

at+1 ⟨Vψ (z

t+1

), z

− z⟩ .

t=0

Further, with λ ≥ 2L̄c , by Lemmas 2 and 3, we have ∆(z̄T ) ≤

−1 TX

at+1

t=0

≤

−1 TX t=0 −1 TX t=0

at+1

−1

TX −1

max

z∈B∩Q

max

z∈B∩Q

hα

at+1 ⟨Vψ (zt+1 ), zt+1 − z⟩

t=0

∥x0 − x∥2x +

i αy 0 ∥y − y∥2y 2

1 −1 1 · (αx Dx2 + αy Dy2 ) ≤ ϵ, λt+1 2

9Lx more than 34 2α gradient queries to ∇x f and no more than 34 2αyyλ gradient queries xλ to −∇y f . Therefore, the numbers q of gradient queries to ∇x f and −∇y f are bounded by q 9L 9Lx T · 1 + 34 2αx λ and T · 1 + 34 2αyyλ , respectively.

Finally, we conclude with the following guarantee on the oracle cost. Theorem 7. Consider DM-SP with the implementation in Eq. (5.2), applied to PSP . With the same choices of αx , αy , and (λt )t≥0 as in Theorem 1, we have: DM-SP NP = (cx + cy ) SP

L D D 1 L D2 1 Ly D2 1 2Lxy Dx Dy xy x y 2 y 2 x x 2 + 102 cx + cy . ϵ ϵ ϵ ϵ

As detailed in Section 2.3, alternative distributed algorithms either fail to outperform EG in general for PSP 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 EG NP ◦ = (cx + cy ) · SP

L D D Ly Dy2 Lx Dx2 xy x y + + ϵ ϵ ϵ

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 c p p cy q x Lx Dx2 + Ly Dy2 ≫ Lxy Dx Dy + Lxy Dx Dy · Lx Dx2 + Ly Dy2 . cx + cy cx + cy For instance, assuming uniform oracle costs (cx = cy ), this improvement occurs when the diagonal conditioning dominates the cross-coupled conditioning, i.e., Lx Dx2 + Ly Dy2 ≫ p Lxy Dx Dy . To our knowledge, DM-SP is the first method to consistently improve upon the EG oracle cost for PSP .

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 PSP PSP , where the local regularizing components are identically zero (i.e., ψx ≡ 0 and ψy ≡ 0). Establishing lower bounds on this unconstrained, smooth subclass yields lower bounds for the broader class PSP . 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 gradientspan 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 −1}, Agent x and Agent y generate local query trajectories of lengths τxt and 18

τyt . The query points consist of the agent’s local variable and its delayed approximation of the remote variable: t t t,l t,l τy −1 t,l t,l τx −1 and Ẑyt = zt,l Ẑxt = zt,l y = (x̂ , y ) l=0 . x = (x , ŷ ) l=0 S St−1 i i t−1 = Let Zxt−1 = t−1 i=0 Ẑx and Zy i=0 Ẑy denote the accumulated sets of query points −1 −1 from all prior rounds (with Zx = Zy = ∅). Within round t, the sets of queries made up t,l l−1 t−1 ∪ {zt,i }l−1 . Finally, to local step l are denoted by Zxt,l = Zxt−1 ∪ {zt,i x }i=0 and Zy = Zy y i=0 t,τ t

t,τ t

the accumulated queries up to the end of round t are Zxt = Zx x and Zyt = Zy y . ◦ and all rounds t ∈ {0, . . . , T − Assumption 1 (Local variables). For any instance P ∈ PSP 1}, the local variables queried by Agent x and Agent y satisfy, for all respective local steps l: t,l xt,l ∈ x0 + P−1 x span ∇x f (z) | z ∈ Zx , t,l yt,l ∈ y0 + P−1 y span ∇y f (z) | z ∈ Zy . ◦ and all rounds t ∈ Assumption 2 (Remote variables). For any instance P ∈ PSP {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−1 ŷt,l ∈ y0 + P−1 , y span ∇y f (z) | z ∈ Zy t−1 x̂t,l ∈ x0 + P−1 . x span ∇x f (z) | z ∈ Zx ◦ and all rounds t ∈ Assumption 3 (Candidate solutions). For any instance P ∈ PSP {0, . . . , T − 1}, the candidate solutions z̄t+1 = (x̄t+1 , ȳt+1 ) generated at the end of round t satisfy t x̄t+1 ∈ x0 + P−1 x span ∇x f (z) | z ∈ Zx , t ȳt+1 ∈ y0 + P−1 y span ∇y f (z) | z ∈ Zy .

Definition 3. An algorithm is called gradient-span if Assumptions 1 to 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−1 and Zxt−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 xt+1 and yt+1 , generating local trajectories denoted by Ẑx2t and Ẑy2t . In the second round, the agents evaluate partial gradients at the assembled point zt+1 , yielding singleton trajectories Ẑx2t+1 = Ẑy2t+1 = {zt+1 }.

6.2

Lower complexity bounds

Now we provide detailed proofs for the lower complexity bounds of distributed gradientspan algorithms applied to PSP . We focus on a subclass of unconstrained SPs with ◦ . We assume the initial points are x0 = 0 and y0 = 0. Lower ψx = ψy = 0, denoted PSP bounds established for this subclass hold for the general class PSP . Let us consider the spaces Ex = Rnx and Ey = Rny equipped with the standard Euclidean inner product ⟨·, ·⟩ and the corresponding Euclidean norm ∥·∥. The cases with general preconditioned norms ∥·∥x and ∥·∥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 Fx = F (Lx , 0, 0, Dx , 0) denote the set of functions 1 Fx (x, y) = ∥Ax x − bx ∥2 , 2 where the matrix Ax ∈ Rnx ×nx satisfies ∥Ax ∥2 ≤ Lx , and the vector bx ∈ Rnx is such that the linear system Ax x = bx has a solution x∗ ∈ Ex satisfying ∥x∗ ∥ ≤ Dx . Definition 5. Let Fy = F (0, 0, Ly , 0, Dy ) denote the set of functions 1 Fy (x, y) = − ∥Ay y − by ∥2 , 2 where the matrix Ay ∈ Rny ×ny satisfies ∥Ay ∥2 ≤ Ly , and the vector by ∈ Rny is such that the linear system Ay y = by has a solution y∗ ∈ Ey satisfying ∥y∗ ∥ ≤ Dy . Definition 6. Let Fxy = F (0, Lxy , 0, Dx , Dy ) denote the set of functions Fxy (x, y) = ⟨Axy x − bxy , y⟩, where the matrix Axy ∈ Rny ×nx satisfies ∥Axy ∥ ≤ Lxy , and the vector bxy ∈ Rny is such that the linear system Axy x = bxy has a solution x∗ ∈ Ex satisfying ∥x∗ ∥ ≤ Dx . Distributed oracles. We consider the pairs of distributed oracles (Ox , Oy ) that return the partial gradients of the functions in these subfamilies: • For Fx ∈ Fx : Ox (x, y) = ∇x Fx (x, y) = A⊤ x (Ax x−bx ) and Oy (x, y) = −∇y Fx (x, y) = 0. • For Fy ∈ Fy : Ox (x, y) = ∇x Fy (x, y) = 0 and Oy (x, y) = −∇y Fy (x, y) = A⊤ y (Ay y− by ). • For Fxy ∈ Fxy : Ox (x, y) = ∇x Fxy (x, y) = A⊤ xy y and Oy (x, y) = −∇y Fxy (x, y) = bxy − Axy x.

Accuracy measures. For any approximate solution z̄ = (x̄, ȳ) ∈ Ex ×Ey , the restricted primal-dual gaps over the bounded sets evaluate to: • For Fx ∈ Fx : ∆Fx (z̄) = Fx (x̄, ȳ) − minx∈Bx Fx (x, ȳ) = 12 ∥Ax x̄ − bx ∥2 . • For Fy ∈ Fy : ∆Fy (z̄) = maxy∈By Fy (x̄, y) − Fy (x̄, ȳ) = 12 ∥Ay ȳ − by ∥2 . • For Fxy ∈ Fxy : ∆Fxy (z̄) = maxy∈By ⟨Axy x̄ − bxy , y⟩ − minx∈Bx ⟨Axy x − bxy , ȳ⟩ ≥ Dy ∥Axy x̄ − bxy ∥, where the inequality holds since minx∈Bx ⟨Axy x − bxy , ȳ⟩ ≤ ⟨Axy x∗ − bxy , ȳ⟩ = 0. Problem subclasses.

These constructions yield three distinct problem subclasses: • Px = Fx , Ox , Oy , ϵ , • Py = Fy , Ox , Oy , ϵ , and • Pxy = Fxy , Ox , Oy , ϵ .

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 ∥v0 − 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 A ∈ Rm×n and vector b ∈ Rm , we define the k-th Krylov subspace as follows: Hk (A, b) ≜ span A⊤ b, (A⊤ A)A⊤ b, . . . , (A⊤ A)k−1 A⊤ b . Proposition 8. Let L > 0, D > 0, and integer 1 ≤ k ≤ min{m−1,n}−1 . Then, there exist 2 a matrix A = A(L, k) ∈ Rm×n with ∥A∥ ≤ L and a vector b = b(L, D, k) ∈ Rm , such that 1 3L2 D2 min ∥Av − b∥2 ≥ , 32(k + 1)2 v∈Hk (A,b) 2 and the linear system Av = b has a solution v∗ ∈ Rn satisfying ∥v∗ ∥ ≤ D. 21

Proof. Let p = 2k + 1. By the condition k ≤ min{m−1,n}−1 , we have p + 1 ≤ m and p ≤ n. 2 p×p Let Mp ∈ R be the symmetric tridiagonal matrix defined as:   2 −1 0 · · · 0 −1 2 −1 · · · 0     ..  ..  . .  Mp =  0 −1 2 .   . . . . .. .. . . −1  .. 0 0 · · · −1 2 Let Bp ∈ R(p+1)×p be the matrix such that B⊤ p Bp = Mp , defined as: 

1 0 0 −1 1 0    0 −1 1 Bp =   . .. ..  .. . .  0 0 ··· 0 0 ···

··· ··· .. .

 0 0  ..  .  .  .. . 0  −1 1  0 −1

We define the matrix A ∈ Rm×n as the block matrix: L Bp 0(p+1)×(n−p) A= . 2 0(m−p−1)×p 0(m−p−1)×(n−p) The matrix A⊤ A ∈ Rn×n is given by the block-diagonal matrix: L2 Mp 0 . A⊤ A = 0 0 4 p p The spectral norm of A is bounded as ∥A∥ = λmax (A⊤ A) = L2 λmax (Mp ) ≤ L. q 6(p+1) Let γ = D p(2p+1) . We define the vector u ∈ Rp+1 by its coordinates: u1 =

p , p+1

and

ui = −

1 p+1

for 2 ≤ i ≤ p + 1.

(p)

⊤ p By the structure of B⊤ p , we have Bp u = e1 ∈ R . Its squared norm evaluates to:

∥u∥2 =

p 2 1 2 p2 + p p +p − = = . p+1 p+1 (p + 1)2 p+1

We define b ∈ Rm as the block vector: L u . b=γ 2 0m−p−1 The linear system Av = b has a solution v∗ ∈ Rn given by the block vector: " # (p) M−1 e ∗ p 1 v =γ . 0n−p 22

p+1−i ∗ ∗ By the structure of M−1 p , the coordinates of v are vi = γ p+1 for 1 ≤ i ≤ p, and 0 otherwise. Its squared norm evaluates to: p

X γ2 p(2p + 1) ∥v ∥ = j2 = γ2 = D2 . 2 (p + 1) 6(p + 1) ∗ 2

j=1

The Krylov subspace is defined as: Hk (A, b) = span{A⊤ b, (A⊤ A)A⊤ b, . . . , (A⊤ A)k−1 A⊤ b}. (p)

(n)

L ⊤ Since B⊤ p u = e1 , the initial vector evaluates to A b = γ 4 e1 . Successive multipli(n)

cation of e1 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: (n)

(n)

Hk (A, b) = span{e1 , e2 , . . . , ek }. For any v ∈ Hk (A, b), its non-zero support is confined to the first k coordinates. Let vk ∈ Rk denote these first k coordinates, such that the first p coordinates of v are ⊤ k×k be the leading principal submatrix of M . The vp = (vk⊤ , 0⊤ p p−k ) . Let Mk ∈ R squared residual norm for v ∈ Hk (A, b) evaluates to: 1 L2 ∥Av − b∥2 = ∥Bp vp − γu∥2 2 8 L2 ⊤ 2 = ⟨vp , B⊤ B v ⟩ − 2γ⟨v , B u⟩ + γ ⟨u, u⟩ p p p p p 8 2 p L (k) ⟨vk , Mk vk ⟩ − 2γ⟨vk , e1 ⟩ + γ 2 . = 8 p+1 (k)

Denoting v1 = ⟨vk , e1 ⟩, the norm simplifies to: 1 L2 p ∥Av − b∥2 = ⟨vk , Mk vk ⟩ − 2γv1 + γ 2 . 2 8 p+1 (k)

Minimizing this residual norm over vk ∈ Rk yields the optimal solution vk∗ = γM−1 k e1 , (k) −1 ∗ equivalently characterized by Mk vk = γe1 . By the structure of Mk , the coordinates ∗ = γ k+1−i for 1 ≤ i ≤ k. Substituting v = v∗ into the v -dependent terms of vk∗ are vk,i k k k k+1 (k)

∗ ), the minimum value evaluates to: (and noting v1 = ⟨vk∗ , e1 ⟩ = vk,1

⟨vk , Mk vk ⟩ − 2γv1

vk =vk∗

∗ = ⟨vk∗ , Mk vk∗ ⟩ − 2γvk,1 (k)

∗ = ⟨vk∗ , γe1 ⟩ − 2γvk,1 ∗ ∗ ∗ = γvk,1 − 2γvk,1 = −γvk,1 = −γ 2

k . k+1

Substituting these values gives the minimum residual norm over the Krylov subspace: 1 L2 γ 2 p k ∥Av − b∥2 = − 8 p+1 k+1 v∈Hk (A,b) 2 k L2 γ 2 2k + 1 − = 8 2k + 2 k + 1 2 2 L γ = 16(k + 1) 6(2k + 2) L2 D 2 = 16(k + 1) (2k + 1)(4k + 3) 3L2 D2 = . 4(8k 2 + 10k + 3) min

3L D Since 8k 2 + 10k + 3 ≤ 8(k + 1)2 , this residual norm is lower bounded by 32(k+1) 2.

6.2.2

Quadratic subclasses Px and Py

In this section, we apply Proposition 8 to establish lower bounds for the quadratic subclasses Px and Py . 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 ∈ Px of the form Fx (x), the partial gradient with respect to y is identically zero, and the gradient ∇x Fx (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 singlenode unconstrained convex optimization. That is, for this subclass, Assumptions 1 and 3 become: i,j i xt,l ∈ x0 + P−1 x span ∇x Fx (x ) | 0 ≤ i < t, 0 ≤ j ≤ τx − 1 or i = t, 0 ≤ j ≤ l − 1 , i,j i x̄t+1 ∈ x0 + P−1 x span ∇x Fx (x ) | 0 ≤ i ≤ t, 0 ≤ j ≤ τx − 1 . Notice that the gradient span sequences are naturally confined to the standard Krylov subspaces. Specifically, for Ax ∈ Rnx ×nx and bx ∈ Rnx , the subspace satisfies the algebraic progression: m m+1 (A⊤ (Ax , bx ) x Ax )H (Ax , bx ) ⊆ H

and

m+1 A⊤ (Ax , bx ). x bx ∈ H

Proposition 9. Let M be a distributed method satisfying Assumption 1 and Assumption 3, operating on an instance Fx ∈ Px of the form Fx (x) = 12 ∥Ax x − bx ∥2 . Let x̄ be the candidate solution generated by Agent x after evaluating Nx partial gradients. Then, x̄ ∈ HNx (Ax , bx ). Proof. For Fx ∈ Px , we have ψx = ψy = 0. Let z1 , . . . , zNx be the sequence of query points evaluated by Agent x, where zi = (xi , yi ). Since x0 = 0, Assumption 1 and Assumption 3 require that each local query point xm+1 and the candidate solution x̄ reside in the span of historical gradients. Let Sm ≜ span{∇x Fx (zi ) | 1 ≤ i ≤ m} ⊆ Ex . We have xm+1 ∈ Sm and x̄ ∈ SNx . 24

We show Sm ⊆ Hm (Ax , bx ) by induction. The base case m = 0 holds since S0 = {0} ⊆ H0 (Ax , bx ). Assume Sm ⊆ Hm (Ax , bx ) for some m ≥ 0. The gradient at xm+1 evaluates to ⊤ m ∇x Fx (xm+1 ) = A⊤ x Ax xm+1 − Ax bx . Since xm+1 ∈ Sm ⊆ H (Ax , bx ), applying the algebraic progression properties yields: m ⊤ m+1 ∇x Fx (xm+1 ) ∈ (A⊤ (Ax , bx ). x Ax )H (Ax , bx ) − Ax bx ⊆ H

Thus, Sm+1 = Sm + span{∇x Fx (xm+1 )} ⊆ Hm+1 (Ax , bx ). By induction, SNx ⊆ HNx (Ax , bx ), which implies x̄ ∈ HNx (Ax , bx ). Analogously, for any instance Fy ∈ Py evaluated by Agent y, of the form Fy (x, y) = − 12 ∥Ay y−by ∥2 , the candidate solution generated after Ny queries satisfies ȳ ∈ HNy (Ay , by ). Theorem 10. Let M be a distributedqmethod satisfying Assumption 1 and Assumption 3 q for problem class PSP , where nx ≥ 2

3Lx Dx2 32ϵ + 2 and ny ≥ 2

M NP ≥ cx x

3Ly Dy2 32ϵ + 2. Then, we have

3Lx Dx2 −1 , 32ϵ

and analogously, s

3Ly Dy2 −1 . 32ϵ q 2 x Dx Proof. We prove the bound for Px . Let K = ⌊ 3L32ϵ ⌋. Since nx ≥ 2K + 2, we apply √ √ Proposition 8 to obtain AK = A( Lx , K) and bK = b( Lx , Dx , K) constructed in Rnx ×nx and Rnx . We set Ax = AK , bx = bK , and define FK (x, y) = 12 ∥Ax x − bx ∥2 . By Proposition 8, ∥Ax ∥2 ≤ Lx , and the linear system Ax x = bx has a solution x∗ satisfying ∥x∗ ∥ ≤ Dx . Hence, FK ∈ Px . When M is applied to PK , suppose itqgenerates an ϵ-saddle point z̄ utilizing Nx M NP ≥ cy y

x Dx queries to Ox . If Nx ≥ K, the bound Nx ≥ 3L32ϵ − 1 holds. If Nx < K, Proposition 9 N x implies x̄ ∈ H (Ax , bx ). Bounding the restricted duality gap via Proposition 8 yields:

1 1 3Lx Dx2 ∥Ax v − bx ∥2 ≥ ϵ ≥ ∆FK (z̄) = ∥Ax x̄ − bx ∥2 ≥ min . 2 32(Nx + 1)2 v∈HNx (Ax ,bx ) 2 q q 2 3Lx Dx2 x Dx M ≥c Rearranging gives Nx ≥ 3L32ϵ − 1. Multiplying by cx yields NP −1 . x 32ϵ x p p The proof for Py is symmetric, setting Ay = A( Ly , K) and by = b( Ly , Dy , K) q 3Ly Dy2 for K = ⌊ 32ϵ ⌋.

6.2.3

Bilinear subclass Pxy

In this section, we establish the lower bounds for the bilinear subclass Pxy . 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 Pxy of the form F (x, y) = ⟨Ax−b, y⟩, where A ∈ Rny ×nx and b ∈ Rny . The partial gradients are ∇x F (x, y) = A⊤ y and ∇y F (x, y) = Ax − b. For any z̄ = (x̄, ȳ) ∈ Ex × Ey , the restricted primal-dual gap evaluates to: ∆F (z̄) ≥ max ⟨Ax̄ − b, y⟩ = Dy ∥Ax̄ − b∥. y∈By

We define the coupled Krylov subspaces in Ey and Ex as: Hyk (A, b) ≜ span b, (AA⊤ )b, . . . , (AA⊤ )k−1 b , Hxk (A, b) ≜ span A⊤ b, (A⊤ A)A⊤ b, . . . , (A⊤ A)k−1 A⊤ b ≡ Hk (A, b). By definition, these subspaces satisfy the alternating properties: AHxm (A, b) + span{b} = Hym+1 (A, b)

and

A⊤ Hym (A, b) = Hxm (A, b) ⊆ Hxm+1 (A, b).

Proposition 11. Let M be a distributed method satisfying Assumption 2 and Assumption 3, operating on F (x, y) = ⟨Ax − b, y⟩. Let x̄ be the candidate solution generated ⌈(T −1)/2⌉ after T communication rounds. Then, x̄ ∈ Hx (A, b). Proof. Let Sxt ≜ span{∇x F (z) | z ∈ Zxt } ⊆ Ex and Syt ≜ span{−∇y F (z) | z ∈ Zyt } ⊆ Ey . ⌈t/2⌉

⌊t/2⌋+1

We show by induction that Sxt ⊆ Hx (A, b) and Syt ⊆ Hy (A, b). The base case t = −1 holds since Sx−1 = {0} ⊆ Hx0 (A, b) and Sy−1 = {0} ⊆ Hy0 (A, b). Assume the claim holds for round t − 1. In round t, Agent x queries points using remote variables ŷ. By Assumption 2, ŷ ∈ Syt−1 . Applying the alternating properties, the evaluated gradient satisfies: ∇x F (x, ŷ) = A⊤ ŷ ∈ A⊤ Hy⌊(t−1)/2⌋+1 (A, b) = Hx⌊(t−1)/2⌋+1 (A, b) = Hx⌈t/2⌉ (A, b). ⌈t/2⌉

Thus, Sxt ⊆ Hx (A, b). Similarly, Agent y queries points using remote variables x̂. By Assumption 2, x̂ ∈ Sxt−1 . The evaluated gradient satisfies: −∇y F (x̂, y) = b−Ax̂ ∈ span{b}+AHx⌈(t−1)/2⌉ (A, b) = Hy⌈(t−1)/2⌉+1 (A, b) = Hy⌊t/2⌋+1 (A, b). ⌊t/2⌋+1

Thus, Syt ⊆ Hy

(A, b). By induction, the claim holds for all t. By Assumption 3, ⌈(T −1)/2⌉

the candidate solution satisfies x̄ ∈ SxT −1 ⊆ Hx

(A, b).

Theorem 12. Let M be a distributed method satisfying Assumption 2 and Assumption 3 2L D D for problem class Pxy , where min{nx , ny } ≥ xy3ϵ x y + 2. Then, we have M TP ≥ xy

2Lxy Dx Dy − 2. 3ϵ

D D

Proof. Let K = ⌊ xy 3ϵx y ⌋. Since nx , ny ≥ 2K + 2, we apply Proposition 8 to obtain AK = A(Lxy , K) ∈ Rny ×nx and bK = b(Lxy , Dx , K) ∈ Rny . We set Axy = AK and bxy = bK , and define FK (x, y) = ⟨Axy x − bxy , y⟩. By Proposition 8, ∥Axy ∥ = ∥AK ∥ ≤ Lxy . The linear system Axy x = bxy has a solution x∗ satisfying ∥x∗ ∥ ≤ Dx . Thus, FK ∈ Pxy . When M is applied to PK , suppose it generates an ϵ-saddle point (x̄, ȳ) after T com2L D D munication rounds. If T ≥ 2K, then T ≥ xy3ϵ x y − 2 holds. If T < 2K, Proposition 11 implies x̄ ∈ Hxk (Axy , bxy ), where k = ⌈(T − 1)/2⌉. Since Hxk (Axy , bxy ) ≡ Hk (AK , bK ), the restricted duality gap on PK evaluates to: ∆FK (x̄, ȳ) ≥ Dy ∥AK x̄ − bK ∥. Applying Proposition 8 yields: s ϵ ≥ ∆FK (x̄, ȳ) ≥ Dy

min

∥AK v − bK ∥ ≥ Dy

v∈Hk (AK ,bK )

Rearranging gives k ≥ 2Lxy Dx Dy − 2. 3ϵ

Lxy Dx Dy − 1. 3ϵ

3L2xy Dx2 Lxy Dx Dy ≥ . 2(8k 2 + 10k + 3) 3(k + 1)

Because k = ⌈(T − 1)/2⌉, we have T ≥ 2k ≥

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 (x, y) = ⟨Ax − b, y⟩. Let (x̄, ȳ) be the solution generated after Nx and Ny queries to Ox and Oy , respectively. Then, x̄ ∈ Hxk (A, b), where k = min(Nx , Ny ). Theorem 14. Let M be a distributed method satisfying Assumption 2 and Assumption 3 2L D D for problem class Pxy , where min{nx , ny } ≥ xy3ϵ x y + 2. Then, we have Lxy Dx Dy M NPxy ≥ (cx + cy ) −1 . 3ϵ 6.2.4

Lower complexity bounds

Finally, we derive the lower bound for PSP by assembling the lower bounds obtained from the three subclasses. Theorem 15. Let M be a distributed gradient-span algorithm for PSP q. Suppose the

dimensions of Ex and Ey are sufficiently large such that nx ≥ q 3Ly Dy2 2Lxy Dx Dy ny ≥ + + 2. We have: 3ϵ 8ϵ

2Lxy Dx Dy + 3ϵ

3Lx Dx2 + 2 and 8ϵ

2Lxy Dx Dy − 2, 3ϵ s r 2 3Ly Dy2 2cx + 2cy c + c L D D c c 3L D x y xy x y y x x x M NP ≥ + + − . SP 9 ϵ 3 32ϵ 3 32ϵ 3 M TP ≥ SP

Proof of Theorem 15. The general class of convex-concave saddle point problems PSP contains the three unregularized subclasses constructed in the previous sections: the xquadratic subclass Px , the y-quadratic subclass Py , and the bilinear subclass Pxy . The worst-case complexity for an algorithm operating over the entire class PSP 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 to 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: 2Lxy Dx Dy − 2. 3ϵ 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: o n M M M M , N , N NP ≥ max N Py . Px Pxy SP M M TP ≥ TP ≥ xy SP

Using the algebraic property max{a, b, c} ≥ 31 (a + b + c), we obtain the lower bound: i 1h M M M M ≥ NP N + N + N Px Py SP 3 Pxy s " # r 2 3Ly Dy2 Lxy Dx Dy 1 3Lx Dx ≥ (cx + cy ) − 1 + cx − 1 + cy −1 3 3ϵ 32ϵ 32ϵ s r 2 cx + cy Lxy Dx Dy cx 3Lx Dx cy 3Ly Dy2 2cx + 2cy = + + − . 9 ϵ 3 32ϵ 3 32ϵ 3 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 communicationoptimal 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.

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. 28

Problem class. We consider a distributed multiagent setting with a star communication network over a product space E = E1 × · · · × EK , where a joint decision variable is partitioned among K autonomous agents as z = (z1 , . . . , zK ). Let us consider the problem class PVIP as follows: • Operators and local components: The problem is governed by a joint monotone operator V (z) = (V1 (z), · · · , VK (z)) and a separable local composite function ψ(z) = P K i=1 ψi (zi ) defined on Q = dom ψ1 × · · · × dom ψK . We assume the solution set is bounded by local distance parameters Di > 0 for each agent. Furthermore, the operator satisfies block-wise Lipschitz continuity: for any fixed z−j , the mapping Vi (zj ; z−j ) is Lij -Lipschitz continuous with respect to zj . • Distributed oracles: Each Agent i ∈ [K] controls its local variable zi ∈ dom ψi , has access to its private function ψi and a partial oracle Oi (z) = Vi (z). • Accuracy measure: The goal is to find an ϵ-approximate solution z̄ ∈ Q such that the restricted gap ∆(z̄) ≜ supz∈B∩Q [⟨V (z), z̄ − z⟩ + ψ(z̄) − ψ(z)] satisfies ∆(z̄) ≤ ϵ, where the bounded domain B is defined by the initial point z0 and distance parameters Di , i ∈ [K]. Conditionings. Let L̄ij ≜ max{L we hereby define the crossPij , Lji }. Similar to SPs,P coupled conditioning, denoted by A with A ≜ D i i i∈[K] i j̸=i L̄ij Dj , that quantifies the cross-dependencies between the agents over the network. The diagonal conditioning, P B with Bi ≜ Lii Di2 , measures the self-dependency within a single denoted by i i∈[K] agent’s domain. New state-of-the-art communication cost. In the multiagent setting, P ExtragraAi +Bi EG dient (EG) remains the state-of-the-art method, requiring TP = O i∈[K] ϵ VIP 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 PVIP . We have X DM-VIP TP = O A /ϵ . i VIP i∈[K]

Remark 5. The communication cost in Theorem 16 completely drops the dependence on the diagonal conditioning. Consequently, we have substantially improved state-ofP the P the-art EG cost when the diagonal conditioning dominates, i.e., when i Bi ≫ i Ai .

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 29

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 Aharon Ben-Tal and Arkadi Nemirovski. Robust optimization–methodology and applications. Mathematical programming, 92(3):453–480, 2002. 1 Aleksandr Beznosikov, Valentin Samokhin, and Alexander Gasnikov. Distributed saddle point problems: lower bounds, near-optimal and robust algorithms. Optimization Methods and Software, pages 1–18, 2025. 1 Radu Ioan Boţ and Enis Chenchene. Extra-gradient method with flexible anchoring: Strong convergence and fast residual decay. arXiv preprint arXiv:2410.14369, 2024. C.4.2 Vincent Conitzer and Tuomas Sandholm. Communication complexity as a lower bound for learning in games. In Proceedings of the twenty-first international conference on Machine learning, page 24, 2004. 1 Yuyang Deng and Mehrdad Mahdavi. Local stochastic gradient descent ascent: Convergence analysis and communication efficiency. In International Conference on Artificial Intelligence and Statistics, pages 1387–1395. PMLR, 2021. 1 Ian J Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. Advances in neural information processing systems, 27, 2014. 1 Sergiu Hart and Yishay Mansour. How long to equilibrium? the communication complexity of uncoupled equilibrium procedures. Games and Economic Behavior, 69(1): 107–126, 2010. 1 Junling Hu, Michael P Wellman, et al. Multiagent reinforcement learning: theoretical framework and an algorithm. In ICML, volume 98, pages 242–250, 1998. 1 Anatoli Juditsky, Arkadi Nemirovski, et al. First order methods for nonsmooth convex large-scale optimization, ii: utilizing problems structure. Optimization for Machine Learning, 30(9):149–183, 2011. 2.3, 17, C, 20 Guanghui Lan and Yan Li. A novel catalyst scheme for stochastic minimax optimization. Mathematical Programming, pages 1–49, 2026. 2.3, 18 Guanghui Lan, Yuyuan Ouyang, and Zhe Zhang. Optimal and parameter-free gradient minimization methods for convex and nonconvex optimization. arXiv preprint arXiv:2310.12139, 2023. 4.4, 4.4, 5, B Tianyi Lin, Chi Jin, and Michael I Jordan. Near-optimal algorithms for minimax optimization. In Conference on learning theory, pages 2738–2779. PMLR, 2020. 2.3

Brendan McMahan, Eider Moore, Daniel Ramage, Seth Hampson, and Blaise Aguera y Arcas. Communication-efficient learning of deep networks from decentralized data. In Artificial intelligence and statistics, pages 1273–1282. PMLR, 2017. 1 Renato DC Monteiro and Benar Fux Svaiter. An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods. SIAM Journal on Optimization, 23(2):1092–1125, 2013. 4.2 Arkadi Nemirovski. Prox-method with rate of convergence O(1/t) for variational inequalities with Lipschitz continuous monotone operators and smooth convex-concave saddle point problems. SIAM Journal on Optimization, 15(1):229–251, 2004. 1, 2.3, 4.1, 4.1, C Arkadij Semenovič Nemirovskij and David Borisovich Yudin. Problem complexity and method efficiency in optimization. Wiley-Interscience, 1983. 2.1 Yurii Nesterov. High-order reduced-gradient methods for composite variational inequalities. arXiv preprint arXiv:2311.15154, 2023. 3, 4.1, 4.1, 4.2 Yurii E. Nesterov. Introductory Lectures on Convex Optimization - A Basic Course, volume 87 of Applied Optimization. Springer, 2004. ISBN 978-1-4613-4691-3. doi: 10.1007/978-1-4419-8853-9. URL https://doi.org/10.1007/978-1-4419-8853-9. 6, 6.2.1 Noam Nisan and Ilya Segal. The communication requirements of efficient allocations and supporting prices. Journal of Economic Theory, 129(1):192–224, 2006. 1 J Ben Rosen. Existence and uniqueness of equilibrium points for concave n-person games. Econometrica: Journal of the Econometric Society, pages 520–534, 1965. 1 Mark Schmidt, Nicolas Le Roux, and Francis Bach. Minimizing finite sums with the stochastic average gradient. Mathematical Programming, 162(1):83–112, 2017. 1 Paul Tseng. On linear convergence of iterative methods for the variational inequality problem. Journal of Computational and Applied Mathematics, 60(1-2):237–252, 1995. 1 John von Neumann and Oskar Morgenstern. Theory of games and economic behavior. Princeton university press, 1947. 1 Yuanhao Wang and Jian Li. Improved algorithms for convex-concave minimax optimization. Advances in Neural Information Processing Systems, 33:4800–4810, 2020. 2.3, A, 19 Junchi Yang, Siqi Zhang, Negar Kiyavash, and Niao He. A catalyst framework for minimax optimization. Advances in Neural Information Processing Systems, 33:5667–5678, 2020. 2.3, 18 TaeHo Yoon and Nicolas Loizou. Pearl-prox: Proximal algorithm for resolving player drift in multiplayer federated learning. In OPT 2025: Optimization for Machine Learning, 2025. 1, 2.3 31

TaeHo Yoon, Sayantan Choudhury, and Nicolas Loizou. Multiplayer federated learning: Reaching equilibrium with less communication. In The Thirty-ninth Annual Conference on Neural Information Processing Systems, 2025. URL https://openreview.net/ forum?id=9JX8XrTVEz. 1, 2.3 Siqi Zhang, Sayantan Choudhury, Sebastian U Stich, and Nicolas Loizou. Communicationefficient gradient descent-ascent methods for distributed variational inequalities: Unified analysis and local updates. In The Twelfth International Conference on Learning Representations, 2024. 1, 2.3 Ali Zindari, Parham Yazdkhasti, Anton Rodomanov, Tatjana Chavdarova, and Sebastian U Stich. Decoupled SGDA for games with intermittent strategy communication. In Forty-second International Conference on Machine Learning, 2025. URL https://openreview.net/forum?id=ZYkFTSEZ6k. 1, 2.3, 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. ◦ denotes the subclass of P We recall that PSP SP where the local components ψx and ψ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 Dx and Dy in advance, but they have access to upper estimates D̂x ≥ Dx and D̂y ≥ Dy . Let θ≜

Dx D̂y D̂x Dy

Dy D̂x D̂y Dx

which quantifies the disproportionality between the true distance parameters and their estimates. Extragradient (EG). Suppose we apply EG to a non-composite problem instance ◦ . The method iteratively maintains and updates two sequences P = (f, 0, 0, z0 ) ∈ PSP k k of variables v , z ∈ Ex × Ey . Initializing with v0 = z0 , the updates for iteration k = 0, 1, . . . , K − 1 are given by: zk+1 = vk − η k P−1 V f (vk ), vk+1 = vk − η k P−1 V f (zk+1 ), where η 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 2K communication rounds. In each iteration, the agents must synchronize twice to evaluate the coupled partial gradients at vk and at zk+1 . Thus, the EG method corresponds to a distributed gradient-span algorithm with T = 2K rounds and local step lengths τxt = τyt = 1, where the trajectories are given by: Ẑx2k = Ẑy2k = {vk }

and Ẑx2k+1 = Ẑy2k+1 = {zk+1 },

for k = 0, . . . , K − 1. To align with Definition 1, the method produces a candidate solution z̄t+1 after each round t ∈ {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 to round t = 2k + 1), it outputs the ergodic Pk corresponds i zi+1 . During the intermediate rounds (at t = 2k), it average z̄2k+2 = Pk1 i η i=0 i=0 η

simply retains the previous solution by setting z̄2k+1 = z̄2k (where z̄0 = z0 ). Since the trajectories and the candidate solutions are formed by linear combinations of the evaluated gradients, they satisfy the span conditions in Definition 3. 33

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 L D̂ +L

D̂

L D̂ +L

xy y ◦ . With the parameter choices of α = x x PSP , αy = y y xy x D̂x D̂y we have Ly Dy2 Lxy Dx Dy Lx Dx2 EG TP + + , ◦ ≤ θ · SP ϵ ϵ ϵ L D D Ly Dy2 Lx Dx2 xy x y EG NP + + . ◦ ≤ (cx + cy ) · θ · SP ϵ ϵ ϵ

D̂x

, and η k ≡ 1,

◦ with a fixed local Decoupled GDA. Let us apply DGDA to an instance P ∈ PSP trajectory length τ . The algorithm maintains local iterates xt,l and yt,l for round t = 0, . . . , T −1 and local step l = 0, . . . , τ . At the beginning of round t, the agents synchronize by exchanging their latest local iterates. Specifically, Agent x receives ŷt and Agent y receives x̂t , defined as: ( ( 0 x if t = 0 y0 if t = 0 t x̂t ≜ and ŷ ≜ . t−1,τ t−1,τ x if t > 0 y if t > 0

The agents then initialize their local variables for the current round as xt,0 = x̂t and yt,0 = ŷt . With the remote variables firmly fixed, the agents execute τ local gradient steps. For l = 0, . . . , τ − 1, the local updates are given by: t,l t xt,l+1 = xt,l − ηx P−1 x ∇x f (x , ŷ ), t t,l yt,l+1 = yt,l + ηy P−1 y ∇y f (x̂ , y ),

where ηx , ηy > 0 are the local step sizes. The DGDA algorithm yields the following trajectories: τ −1 τ −1 Ẑxt = (xt,l , ŷt ) l=0 and Ẑyt = (x̂t , yt,l ) l=0 . After each round t ∈ {0, . . . , T − 1}, the method returns the updated local iterates as the candidate solution z̄t+1 = (xt,τ , yt,τ ). 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 convexstrongly 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 O log 1ϵ , which is a clear improvement over the ◦ with stronger coupling, the EG baseline. However, for general problem instances in PSP delayed remote variables cause the local updates to drift, ultimately leading the method to diverge. Consequently, DGDA fails to provide any meaningful complexity guarantee ◦ under consideration. for the general problem class PSP

Catalyst acceleration. Catalyst methods first add small regularization terms to the objective: f (x, y) + ϵ 2 ∥x − x0 ∥2x − ϵ 2 ∥y − y0 ∥2y , reducing the problem to a strongly 4D̂x

4D̂y

convex-strongly concave one. Then, Catalyst introduces an outer loop, indexed by k = 0, 1, . . . , K − 1, designed to balance the conditioning between the two variables. For instance, when the conditioning of x is worse (i.e. Lx D̂x2 ≥ Ly D̂y2 ), the method carefully maintains an extrapolation sequence (x̃k )K−1 k=0 and, in each outer iteration, adds a proximal λx k 2 term 2 ∥x − x̃ ∥x 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 = max{Lx , Ly , Lxy }. 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 PSP Then, we have: s s 2 Ly D̂y2 L D̂ D̂ L D̂ max x y x x Cat-EG 2 1 TP ◦ =O + + log , SP ϵ ϵ ϵ ϵ s s 2 Ly D̂y2 Lmax D̂x D̂y Lx D̂x Cat-EG 2 1 NP ◦ = (cx + cy ) · O + + log . SP ϵ ϵ ϵ ϵ 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 and D̂y to set the initial regularization, its complexity scales directly with these estimates rather than the true distances Dx and Dy . Specifically, the diagonal terms in EG or DM-SP depend strictly on the true distances, whereas in Cat-EG they scale with D̂x and 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 and D̂y = cDy ), θ remains 2, leaving the complexity unaffected by the overestimation factor c. In contrast, the Cat-EG coupled term scales with D̂x D̂y , meaning any overestimation of D̂x ≫ Dx or D̂y ≫ Dy 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 = 106 , Lxy = 1, Ly = 10−6 , Dx = 10−3 , and Dy = 103 . 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 convexstrongly 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; (ii) The second loop is another Catalyst outer loop adding a proximal 35

regularization term to the dual variable; (iii) 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 ◦ . Then, we have: to PSP 1 Lxy D̂x D̂y , log3 SP ϵ ϵ s s p 2 Ly D̂y2 Lmax Lxy D̂x D̂y Lx D̂x Cat-Cat-DAGDA 4 1 NP ◦ = (cx + cy ) · O log . + + SP ϵ ϵ ϵ ϵ Cat-Cat-DAGDA TP =O ◦

MRN solver: Accumulative Regularization Method

Algorithm B.1 ARM ∇fw , ψw , vw , ξ | L 1: Set τ = 2 + max 0, log4 3L and σ (0) = 0. 2ξ l q m 2ξ L 2: Set σ (k) = 4k−3 3 and Nk = 16 for k = 1, . . . , τ . σ (k) 3: Initialize w̄(0) = w(0) = vw . 4: for k = 1, . . . , τ do 5: 6:

7: 8: 9: 10: 11: 12: 13:

(k−1)

γ (k) = 1 − σσ(k) w̄(k) = (1 − γ (k) )w̄(k−1) + γ (k) w(k−1) % Begin Inner Subroutine: Nesterov’s Accelerated Gradient Method (0) (0) Initialize xk = w(k−1) , yk = w(k−1) , t0 = 1, and Lk = L + σ (k) . for i = 0, . . . , Nk − 1 do (k) (i) (i) (i) ∇fw (yk ) = ∇fw (yk ) + σ (k) (yk − w̄(k) ) (i+1) (k) (i) (i) xk = arg minw∈dom ψw ⟨∇fw (yk ), w⟩ + L2k ∥w − yk ∥2w + ψw (w) √ 1+ 1+4t2i ti+1 = 2 (i+1)

(i+1)

−1 yk = xk + ttii+1 (xk end for % End Inner Subroutine

(i)

− xk )

(N )

w(k) = xk k ′ (w(k) ) = −∇f (k) (y(Nk −1) ) − L (x(Nk ) − y(Nk −1) ) 15: ψw w k k k k 16: end for ′ (w(τ ) )) 17: return (w(τ ) , ψw 14:

Proof of Lemma 5. It was q originally shown in [Lan et al., 2023, Theorem 3.1] that ARM takes no more than 34

3L ′ 2ξ gradient queries and obtains w ∈ dom ψw , such that there

exists w̃ ∈ dom ψw with ∇fw (w̃) ∈ −∂ψw (w̃), and 2 ∥2L(w+ −w′ )∥w ≤ ξ∥vw −w̃∥w , where w+ = arg min ⟨∇fw (w′ ), w⟩+ψw (w)+L∥w−w′ ∥2w . 3 w∈dom ψw ′ (w+ ) ∈ ∂ψ (w+ ) such that By the optimality of w+ , there exists ψw w ′ ∇fw (w′ ) + ψw (w+ ) + 2LP(w+ − w′ ) = 0.

Then, we have ′ ∥∇fw (w+ ) + ψw (w+ )∥w∗ ′ ≤ ∥∇fw (w′ ) + ψw (w+ )∥w∗ + ∥∇fw (w+ ) − ∇fw (w′ )∥w∗

= ∥2LP(w+ − w′ )∥w∗ + ∥∇fw (w+ ) − ∇fw (w′ )∥w∗ ≤ 3L∥w+ − w′ ∥w ≤ ξ∥vw − w̃∥w .

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 generalsum games.

C.1

Problem formulation

VIPs (with separable composite terms). Let us consider the VIP in Eq. (4.1), where E = E1 × · · · × EK is the direct product of K finite-dimensional real vector spaces. For all i ∈ [K]: let the mapping Vi : dom ψ → Ei∗ , and let the function ψi : Ei → R∪{+∞}. We consider the decomposition of V (z) = V1 (z), · · · , VK (z) and ψ(z) = ψ1 (z1 ) + · · · + ψK (zK ), for all z = (z1 , · · · , zK ) ∈ E. Moreover, we denote dom ψ = dom ψ1 × · · · × dom ψK ≜ Q. Assumptions for VIPs.

Let us make the following assumptions:

(A2’) Let z0 = (z01 , · · · , z0K ) ∈ Q be a given point. There exists z∗ = (z∗1 , · · · , z∗K ) ∈ Q in the solution set of the VIP of (V, ψ), such that for all i ∈ [K]: z∗i ∈ Bi , where Bi ≜ {zi ∈ Ei | ∥z0i − zi ∥i ≤ Di } and Di > 0 is a given distance. (A3’) The operator Vi (zj ; z−j ) is Lij -Lipschitz continuous in zj ∈ dom ψj for any fixed z−j ∈ dom ψ1 × · · · × dom ψj−1 × dom ψj+1 × · · · × dom ψK . 1 Let the operator family FVIP be comprised of all the operators with initialization points 0 (Vi )i∈[K] , z ∈ PVIP , such that Assumptions (Â1) to (Â3) are satisfied. 1

For all z = (z1 , · · · , zK ) ∈ Q, we use the following notations for simplicity: (zj ; z−j ) ≜ z and z−j ≜ (z1 , · · · , zj−1 , zj+1 , · · · , zK ).

Notations. To simplify the notations, let us denote z ≜ (z1 , · · · , zK ) ∈ Q in the context of VIPs. Let us denote X L̄ij ≜ max{Lij , Lji }, Ai ≜ Di L̄ij Dj , and Bi ≜ L̄ii Di2 , for all i, j ∈ [K]. j∈[K]\{i}

P P We refer to i Ai as the cross-coupled conditioning and i Bi asPthe diagonal P conditioning, and we say that the diagonal conditioning dominates when i Bi ≫ i Ai .

C.2

Communication and computational costs

Distributed oracles. We consider a distributed setting with K agents very similar to the one in Sections 2.1 and 2.2. For all i ∈ [K]: Agent i controls decision variable zi ∈ dom ψi , has direct access to the function ψi , and has access to the oracle Oi (z) = Vi (z) for z ∈ Q. We consider per query to Oi costs ci ≥ 0, i ∈ [K]. Accuracy measure.

We consider the following accuracy measure for VIPs:

∆(z̄) ≜ sup ⟨V (z), z̄ − z⟩ + ψ(z̄) − ψ(z), for all z ∈ Q, z∈B∩Q

where B ≜ B1 × · · · × BK . We say that a point z̄ ∈ Q is an ϵ-approximate solution of the VIP if ∆(z̄) ≤ ϵ. Our goal is to find such an ϵ-approximate solution for any ϵ > 0. Problem class. We formally define the overall problem class, denoted by PVIP FVIP , (Oi )i∈[K] , ϵ , or for short PVIP . A specific problem instance P ∈ PVIP is constructed by drawing an operator instance (with initial point) V from FVIP , equipping it with the distributed oracles (Oi )i∈[K] , and specifying a target accuracy ϵ > 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 PVIP , 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 − 1}, each Agent i ∈ [K] executes multiple local computational steps to generate local query points denoted by t,l t,l t,l zt,l i = (zi,1 , . . . , zi,i , . . . , zi,K )

for l ∈ {0, . . . , τit − 1},

t,l where zt,l i,i is the local variable updated by Agent i, and zi,j (for j ̸= i) is the delayed approximation of Agent j’s variable utilized by Agent i. Let Iit,l denote the accumulated information sequence available to Agent i prior to making its (l + 1)-th local oracle query in round t. The base case at initialization is Ii0,0 = (ψi , z0 ). During the local computational steps l ∈ {0, . . . , τit − 1}, the information

sequence of Agent i updates sequentially by appending the newly acquired oracle response for Vi : Iit,l+1 = Iit,l , Vi (zt,l i ) . t After the local steps in round t, the agents exchange messages. Let Mj→i denote the message sent from Agent j to Agent i. The information sequence available to Agent i at the beginning of round t + 1 appends the received messages to its prior local history:

t,τ t t Iit+1,0 = Ii i , {Mj→i }j̸=i . Definition 7. An algorithm M is called a distributed algorithm for problem class PVIP if, when applied to any instance P ∈ PVIP , its execution satisfies the following conditions for all t ∈ {0, . . . , T − 1} and agents i ∈ [K]: 1. Local Computation: The query points are determined entirely by the locally available information. For all local steps l ∈ {0, . . . , τit − 1}, there exists a deterministic mapping function At,l i such that: t,l t,l zt,l i = Ai Ii .

2. Communication: The messages exchanged are produced by deterministic mappings of the sender’s local information. For any j ̸= i, there exists a mapping function Mti→j such that: t,τ t t Mi→j = Mti→j Ii i . t+1 3. Candidate solution: The candidate solutions z̄t+1 = (z̄t+1 1 , . . . , z̄K ) are constructed from the respective agents’ updated information sets. There exists a deterministic mapping function Āt+1 such that: i

Iit+1,0 . z̄t+1 = Āt+1 i i For a given instance P ∈ PVIP and a target accuracy ϵ > 0, we define the communication cost required by a distributed algorithm M on P , denoted by TPM , as the smallest integer k ∈ {1, . . . , T } such that the candidate solution z̄k satisfies the target accuracy ϵ. The total number of local oracle queries evaluated by Agent i for instance P up to M = this point is given by the cumulative number of local steps taken, denoted by Ni,P PTPM −1 r τi . r=0 Let ci denote the computational cost of evaluating a single partial oracle Vi . The communication cost and oracle cost of algorithm M over the entire problem class PVIP are defined by taking the supremum over all instances: X M M M M TP = sup T and N = sup c N i i,P . P PVIP VIP P ∈PVIP

P ∈PVIP

i∈[K]

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, the communication cost of EG is bounded by O

X Ai + B i ϵ

i∈[K]

and the computational cost of EG is bounded by O

X i∈[K]

C.3

X X Bi X Ai ci + . ϵ ϵ i∈[K]

i∈[K]

Decoupled method for variational inequality problems

Now, we present our DM-VIP method, which extends the DM-SP into multiplayer generalsum games. Assembled norm. Given parameters αi > 0 for all i ∈ [K] (to be specified later), we equip the joint space E = E1 × · · · × EK with the assembled norm: v uK uX 1 ∥z∥E = ⟨Pz, z⟩ 2 = t αi ∥zi ∥2i for all z ∈ E, (C.1) i=1

which corresponds to the block diagonal linear operator P = α1 P1 ⊕ · · · ⊕ αK PK . Ac∗ with the corresponding dual norm: cordingly, we equip the dual space E ∗ = E1∗ × · · · × EK v uK uX 1 −1 2 αi−1 ∥gi ∥2i∗ for all g ∈ E ∗ . (C.2) ∥g∥E ∗ = ⟨g, P g⟩ = t i=1

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: v P u X L̄ij u l∈[K]\{i} L̄il Dl −1 L̄c ≜ t max (αj Dj ) . (C.3) αi j∈[K] i∈[K]\{j}

Algorithm C.1 outlines the Decoupled Method for block composite VIPs (DM-VIP), generalizing the DM-SP procedure. The algorithm maintains a sequence of anchor points vt 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 t . This allows each Agent i to remote variables to the current anchor components v−i independently and concurrently solve its regularized local subproblem. Specifically, each agent invokes an internal solver MMRN to minimize the local residual norm (MRN) up i αi λt+1 t+1 to a target accuracy δi = 2 (Line 4). Following this local computation phase, the K agents perform exactly two communication rounds to complete the iteration. In the first round (Line 5), the agents broadcast 40

their locally computed approximate solutions zt+1 to assemble the joint intermediate i t+1 point z . In the second round (Line 6), they use this assembled point to evaluate their local partial operators Vi (zt+1 ), which they subsequently exchange to form the full joint operator Vψ (zt+1 ). Finally, using this assembled joint operator, the agents compute a closed-form reduced step size at+1 , update the running ergodic average z̄t+1 , and perform a joint extragradientlike step to generate the next anchor vt+1 (Lines 7 and 8). Algorithm C.1 DM-VIP K, (Vi )i∈[K] , (ψi )i∈[K] , z0 , (λt )t≥1 , (αi )i∈[K] | (MMRN )i∈[K] i

Require: A local solver MMRN for the minimization of residual norms. i 0 ) = z0 . 1: v0 = (v10 , . . . , vK 2: for t = 0, 1, . . . , T − 1 do 3: Let δit+1 = αi λ2t+1 for all i ∈ [K]. 4: Concurrently for all i ∈ [K], Agent i computes αi λt+1 ′ t+1 MRN t zt+1 ∥· − vit ∥2i , vit , δit+1 . Vi (·; v−i ), ψi + i , ψi (zi ) = Mi 2 5: 6:

7: 8:

t+1 All agents exchange zt+1 to assemble the joint point zt+1 = (zt+1 1 , . . . , zK ). i t+1 Agents locally compute Vi (z ) and exchange to assemble the full operator: t+1 ′ Vψ (zt+1 ) = V (zt+1 ) + ψ1′ (zt+1 1 ), . . . , ψK (zK ) . 2⟨Vψ (zt+1 ),vt −zt+1 ⟩ and generate solution z̄t+1 = ∥Vψ (zt+1 )∥2E ∗ vt+1 = arg minv∈Q at+1 ⟨Vψ (zt+1 ), v⟩ + 12 ∥v − vt ∥2E .

Let at+1 =

Pt+1

i=1 ai

−1 Pt+1

i i=1 ai z .

9: end for

Similar to the saddle point setting, we refer to Algorithm C.1 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 MMRN (V̂i , ψ̂i , vi , δi ) capable of taking an MRN instance in the space Ei and suci cessfully returning a sufficiently accurate solution. We defer the detailed implementation of these local solvers to Eq. (C.5) at the end of this section. For now, let us proceed with the communication complexity of the template DM-VIP method. Theorem 21. P Consider the DM-VIP template applied to PVIP . With the parameter choices of αi =

j∈[K]\{i} L̄ij Dj

for all i ∈ [K], and λt ≡ λ ≥ 2L̄c , we have:

DM-VIP TP ≤2+ VIP

X 2L̄ij Di Dj . ϵ

i,j∈[K] i̸=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 X (Ai + Bi ) ϵ i∈[K]

communication rounds (cf. Proposition 20). Our communication cost in Theorem 21 is P O A /ϵ , which is consistently no worse compared to that of EG and is substani i∈[K] tially faster when the “diagonal conditioning” dominates—i.e., X X Ai . Bi ≫ i∈[K]

i∈[K]

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 C.2. Then, we prove the correctness of the solution returned by FDS. Algorithm C.2 FDS∥·∥E (Vi )i∈[K] , (ψi )i∈[K] , v, λ | (MMRN )i∈[K] i

Require: Solver MMRN for the minimization of residual norms, for all i ∈ [K]. i 1: for i ∈ [K] do 2: δi = α2i λ . 3: ψ̂i = ψi + α2i λ ∥· − vi ∥2i . ′ + MRN V (·; v ), ψ̂ , v , δ . 4: (z+ i −i i i i i , ψi (zi )) = Mi + ′ + 5: ψi′ (z+ i ) = ψ̂i (zi ) − αi λPi (zi − vi ). 6: end for + + + ′ (z+ )). 7: return (z+ , ψ ′ (z+ )), where z+ = (z1 , · · · , zK ) and ψ ′ (z+ ) = (ψ1′ (z1 ), · · · , ψK K Lemma 22. Under (A3’), for λ ≥ 2L̄c , FDS (Algorithm C.2) returns the correct solution of the MS subproblem given by (V, ψ, v, λ). Proof of Lemma 22. For all i ∈ [K], by (A3’) and then by the relative distance accuracy, we have + ∥Vi (z+ ) + ψi′ (z+ i ) + αi λPi (zi − vi )∥i∗ + ′ + ≤ ∥Vi (z+ i ; v−i ) + ψi (z ) + αi λPi (zi − vi )∥i∗ +

X j∈[K]\{i}

≤ δi ∥z+ i − v i ∥i +

Lij ∥z+ j − v j ∥j .

j∈[K]\{i}

Lij ∥z+ j − v j ∥j (C.4)

Finally, we assemble the norms: ∥V (z+ ) + ψ ′ (z+ ) + λP(z+ − v)∥2E ∗ X + 2 αi−1 ∥Vi (z+ ) + ψi′ (z+ = i ) + αi λPi (zi − vi )∥i∗ i∈[K] (C.4)

≤

αi−1 δi ∥z+ i − vi ∥i +

i∈[K]

≤

= ≤

2 λ2 2

Lij ∥z+ j − v j ∥j

2αi−1

2 δi2 ∥z+ i − v i ∥i +

Lij ∥z+ j − v j ∥j

2 αi ∥z+ i − v i ∥i + 2

i∈[K]

αi−1

i∈[K]

X X + 2 αi−1 ∥z − v∥E + 2 i∈[K]

≤

C.4.2

Lil Dl

Lij + ∥zj − vj ∥2j Dj j∈[K]\{i} P Lij Lil Dl

l∈[K]\{i}

j∈[K]

Lij ∥z+ j − vj ∥j

j∈[K]\{i}

2 X X ∥z+ λ2 + j − v j ∥j 2 ∥z − v∥E + 2 = 2 Dj

≤

j∈[K]\{i}

(C.3) λ2

j∈[K]\{i}

i∈[K]

λ2

l∈[K]\{i}

αi

i∈[K]\{j}

∥z+ − v∥2E + 2∥z+ − v∥2E L̄2c

λ2 + λ2 ∥z − v∥2E + ∥z+ − v∥2E = λ2 ∥z+ − v∥2E . 2 2

Proof for MRN

Our algorithm is built upon Lemma 23, the proof of which can be found in [Boţ and Chenchene, 2024, Corollary 2.4]. Lemma 23. Assume (Â1), (Â3), and that the solution set of the VIP of (Vw , ψw ) is non′ (w+ )) = FEGM (V , ψ , v , ξ | empty. Then, there exists an algorithm, denoted by (w+ , ψw w w w L ′ (w+ )) that satL), which takes no more than C0 · ξ operator queries and returns (w+ , ψw isfies ξ-distance-to-solution accuracy, where C0 > 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: MFEGM (V̂i , ψ̂i , vi , δi ) ≜ FEGM(V̂i , ψ̂i , vi , i

2δi | Lii ). 3

Then, for any Monteiro-Svaiter Subproblem given by (V, ψ, v, λ), we leverage the solver FDS-FEGM(V, ψ, v, λ) = FDS∥·∥E V, ψ, v, λ | (MFEGM )i∈[K] . i Finally, we obtain the concrete algorithm DM-VIP as follows:

ROM∥·∥E (Vi )i∈[K] , (ψi )i∈[K] , z0 , (λt )t≥1 | FDS-FEGM .

(C.5)

Combining Lemmas 2, 4, 22 and 23, with the implementation in Eq. (C.5), we conclude that Lemma 24 holds for the constant C0 from Lemma 23. We include the complete proof below. Lemma 24. Under (A2’) and (A3’), for λt+1 ≡ λ ≥ 2L̄c , DM-VIP (Algorithm C.1) with the implementation in Eq. (C.5) takes no more than 2T communication rounds and no more than T · 1 + C0 ·

3Lii αi λ

queries to Vi , for all i ∈ [K], and obtains an ϵ-approximate solution z̄T , where P i αi λDi2 T = 2ϵ and C0 > 0 is some fixed constant. Proof of Lemma 24. By (A1’), we have ∆(z̄T ) ≤

−1 TX

at+1

−1

TX −1

max

z∈B∩Q

t=0

at+1 ⟨Vψ (zt+1 ), zt+1 − z⟩ .

t=0

Further, with λ ≥ 2L̄c , by Lemmas 2 and 22, we have T

∆(z̄ ) ≤

−1 TX

at+1

t=0

≤

−1 TX t=0

at+1

−1

max

TX −1

z∈B∩Q

t+1

at+1 ⟨Vψ (z

t+1

), z

− z⟩

t=0

−1 h X α i i max ∥z0i − zi ∥2i 2 zi ∈Bi ∩dom ψi i∈[K]

−1 TX 1 −1 1 X · αi Di2 ≤ ϵ, ≤ λt+1 2 t=0

i∈[K]

where the last inequality follows from the assignments of (λt )t≥1 and T . Therefore, the number of communication rounds is bounded by 2T . 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 ii takes no more than C0 · 3L αi λ gradient queries to Vi, for all i ∈ [K]. Therefore, the numbers ii of queries to Vi are bounded by T · 1 + C0 · 3L αi λ , for all 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 X X 3C0 X Bi ci X Ai + ci Ai . ϵ ϵ Ai i∈[K]

i∈[K]

(C.6)

i∈[K]

Compared to the computational cost of EG, which is given by 1 X X 1 X X Ai + Bi , ci ci ϵ ϵ i∈[K]

i∈[K]

our computational cost in Eq. (C.6) differs in the second term. Consequently, our DM-VIP may offer an advantage or disadvantage depending on the relative conditioning of Ai , Bi , and ci for i ∈ [K].

Record · ID 200445 · SHA-256 5e13044bf624743e
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.