Year: 2021

Abstract:
Many real-world applications of multi-agent reinforcement learning (RL), such as multi-robot navigation and decentralized control of cyber-physical systems, involve the cooperation of agents as a team with aligned objectives. We study multi-agent RL in the most basic cooperative setting — Markov teams — a class of Markov games where the cooperating agents share a common reward. We propose an algorithm in which each agent independently runs stage-based V-learning (a Q-learning style algorithm) to efficiently explore the unknown environment, while using a stochastic gradient descent (SGD) subroutine for policy updates. We show that the agents can learn an ε-approximate Nash equilibrium policy in at most ∝ Õ(1/ε) episodes. Our results advocate the use of a novel stage-based V-learning approach to create a stage-wise stationary environment. We also show that under certain smoothness assumptions of the team, our algorithm can achieve a nearly team-optimal Nash equilibrium. Simulation results corroborate our theoretical findings. One key feature of our algorithm is being decentralized, in the sense that each agent has access to only the state and its local actions, and is even oblivious to the presence of the other agents. Neither communication among teammates nor coordination by a central controller is required during learning. Hence, our algorithm can readily generalize to an arbitrary number of agents, without suffering from the exponential dependence on the number of agents.

heading: 1 INTRODUCTION
text: A variety of real-world sequential decision-making problems have been successfully addressed with reinforcement learning (RL), such as playing the game of Go (Silver et al., 2016), Poker (Brown & Sandholm, 2018), real-time strategy games (Vinyals et al., 2019), autonomous driving (ShalevShwartz et al., 2016), and robotics (Kober et al., 2013). Many of these applications can be cast as a multi-agent reinforcement learning (MARL) problem, where multiple agents are involved in an interactive environment. These successful applications have inspired a surging interest in understanding the theoretical aspects of MARL, which, however, is mostly devoted to the competitive setting, e.g., two-player zero-sum Markov games. This has left many problems in other important MARL settings, especially cooperative MARL, relatively open.
Cooperative MARL generally involves a team of agents collaborating on a common task. We roughly categorize cooperative MARL settings into three classes based on the information available to the agents. The first class is (1) centralized/joint learning (Boutilier, 1996; Claus & Boutilier, 1998), which assumes that there exists a single coordinator who can access the local information of all the agents, and makes decisions/learns policies jointly for all of them. This centralized training (though possibly decentralized execution) approach has become a common practice in empirical MARL (Oliehoek et al., 2008; Foerster et al., 2016; Lowe et al., 2017; Rashid et al., 2018; Son et al., 2019). Centralized learning essentially reduces the multi-agent problem to a single-agent one, but such a method does not scale well with the number of agents. In fact, one can show that the computation & sample complexity grows exponentially with the number of agents in the worst case. The second class is (2) learning with communications (Kar et al., 2013; Zhang et al., 2018; Dubey & Pentland, 2021), which typically assumes that the agents can communicate by sharing some local information with other agents (sometimes through a central agent). This approach gets around the centralized computation bottleneck, but instead suffers from additional communication overhead. It can also be
unrealistic in some real-world scenarios where communication may be expensive and/or unreliable, such as in unmanned aerial vehicle (UAV) field coverage (Pham et al., 2018).
Given the aforementioned limitations, we are more interested, in this paper, in a more practical setting: (3) decentralized learning1. We enforce a seemingly daunting restriction: Each agent only has access to its own local information (e.g., local actions and rewards), and can neither communicate with its teammates nor be coordinated by any central controller during learning. In fact, we require that each agent is completely unaware of the underlying structure of the game (whether it is cooperative or not), or even the presence of other agents. Among the three classes mentioned above, decentralized learning requires the weakest assumption on information availability, which makes it suitable for many practical multi-agent learning scenarios (Fudenberg et al., 1998; Wang & de Silva, 2008; Kuyer et al., 2008). In addition, decentralized learning is generally more scalable, and does not suffer from the exponential sample & computation complexity in the number of agents. Seemingly tempting, one might wonder how strong, if any, theoretical results can be established in this challenging decentralized setting. In the paper, we investigate the theoretical aspects of decentralized cooperative multi-agent RL in the non-asymptotic regime. A natural and fundamental theoretical question is:
Does any independent Q-learning style algorithm lead to Nash equilibria in cooperative MARL?
The answer to this question has remained elusive ever since it was raised in Lauer & Riedmiller (2000). To answer this question, we focus on arguably the most basic and fundamental cooperative MARL setting — Markov teams (Lauer & Riedmiller, 2000; Boutilier, 1996; Wang & Sandholm, 2002) – where the cooperating agents share a common team-reward function, and the state transition in the environment is affected by the joint actions of all agents. In fact, even in the basic Markov team setting, decentralized learning can be highly challenging. First, since both rewards and transitions are affected by the other agents, the environment becomes non-stationary from each agent’s own perspective when the other agents also update their policies. Hence, an agent needs to explore the unknown environment efficiently while keeping in mind that the information it gathers now will soon become outdated. This rules out many successful single-agent RL solutions, which critically rely on the assumption that the agent is learning in a stationary environment. Second, compared with RL in two-player zero-sum games, an additional challenge in teams is the existence of possibly multiple Nash equilibria (NE) with different values. The agents in teams hence have to additionally coordinate on which NE to converge to. Our contributions in the paper are summarized as follows.
Contributions. 1) As a warm-up, we start with decentralized matrix teams (Section 3), which constitute a special class of Markov teams with no state transitions. We show that a simple stochastic gradient descent (SGD) algorithm can provably find an ε-approximate NE policy using ∝ O(1/ε8) samples, which can be improved to ∝ Õ(1/ε4) by using variance reduction techniques. 2) For the general problem of decentralized Markov teams (Section 4), we propose an algorithm named StageBased V-learning with Stochastic Gradient Descent (V-learning SGD), where V-learning (a name first coined in Bai et al. (2020)) is a simple variant of Q-learning. Our algorithm can be executed in a decentralized way as we defined above. 3) We show in Section 5 that if all agents in the team run V-learning SGD, they can find an ε-approximate NE in at most Õ(1/ε4) · poly(S,H,Amax) episodes, where S is the number of states, Amax is the size of the largest action space, and H is the length of an episode. Our results rely on a novel investigation of stage-based V-learning to create a stage-wise stationary environment for the agents. 4) We further show that our algorithm can approach the best equilibrium, i.e., the team-optimality, in a specific class of problems named “smooth teams”. 5) In Section 6, we provide numerical results that corroborate our theoretical findings. To the best of our knowledge, our work appears to be the first to provide non-asymptotic guarantees for cooperative MARL with exploration, while having the additional advantage of being decentralized.
Related Work. A common mathematical framework of multi-agent RL is stochastic games (Shapley, 1953), which is also referred to as Markov games. Early attempts to learn NE in Markov games include Littman (1994; 2001); Hu & Wellman (2003); Hansen et al. (2013), but they either assume the transition kernel and rewards are known, or only yield asymptotic guarantees. Recently, various sample efficient methods have been proposed (Wei et al., 2017; Bai & Jin, 2020; Sidford et al., 2020;
1This setting has been studied under various names in the literature, including individual learning (Leslie & Collins, 2005), decentralized learning (Arslan & Yüksel, 2016), agnostic learning (Tian et al., 2021; Wei et al., 2021), and independent learning (Claus & Boutilier, 1998; Daskalakis et al., 2020). It also belongs to a more general category of teams/games with decentralized information structure (Ho, 1980; Nayyar et al., 2013a;b).
Xie et al., 2020; Bai et al., 2020; Liu et al., 2021; Zhao et al., 2021; Guo et al., 2021), mostly for learning in two-player zero-sum Markov games. Several works have investigated zero-sum games in a decentralized setting as we consider here (Daskalakis et al., 2020; Tian et al., 2021; Wei et al., 2021; Sayin et al., 2021), but these results do not carry over in any way to the decentralized cooperative/team setting. We refer the reader to Appendix A for a more detailed discussion on these related works.
Another line of research has considered RL in teams. Without enforcing a decentralized environment, Boutilier (1996) has proposed to coordinate the agents by letting them take actions in a lexicographic order. Wang & Sandholm (2002) have studied optimal adaptive learning that converges to the optimal NE in Markov teams. These two methods critically rely on communications among the agents (beforehand) and/or observing the teammates’ actions. In contrast, the distributed Q-learning algorithm in Lauer & Riedmiller (2000) is decentralized and coordination-free, which, however, only works for deterministic tasks, and has no non-asymptotic convergence guarantees. In fact, developing provable decentralized RL for cooperation in general stochastic environments is considered an important open problem in Lauer & Riedmiller (2000), which we address in this paper.
Arslan & Yüksel (2016) and Yongacoglu et al. (2019) have shown that decentralized Q-learning style algorithms can converge to NE in weakly acyclic games, which cover Markov teams as an important special case. Their decentralized settings are most similar to ours in that each agent is oblivious to the presence of the others. However, both of them require a coordinated exploration phase, and only yield asymptotic guarantees. In fact, it is not clear if their methods will converge in polynomial time. We significantly improve their results in this respect, by explicitly showing a polynomial sample complexity upper bound. Two contemporaneous works (Zhang et al., 2021b; Leonardos et al., 2021) have studied gradient play in Markov potential games,2 but they do not consider the aspect of exploring the unknown environment, one of the key issues in RL, and/or have to assume perfect knowledge of the environment. A concurrent work (Chang et al., 2021) has studied cooperative multi-player multi-armed bandits with information asymmetry, which parallels our matrix team setting. Nevertheless, (Chang et al., 2021) requires stronger conditions than our decentralized setting as their algorithm relies on playing a predetermined sequence of actions.

heading: 2 PRELIMINARIES
text: An N -player episodic Markov team is defined by a tuple (N , H,S, {Ai}Ni=1, r, P ), where (1) N = {1, 2, . . . , N} is the set of agents; (2) H ∈ N+ is the number of time steps in each episode; (3) S is the finite state space; (4) Ai is the finite action space for agent i ∈ N ; (5) r : [H]× S ×A → [0, 1] is the team reward function for all agents, where A = ×Ni=1Ai is the joint action (or action profile) space; and (6) P : [H]× S ×A → ∆(S) is the transition kernel. We remark that both the reward function and the state transition function depend on the joint actions of all the agents.
The agents interact in an unknown environment for K episodes. At each time step h ∈ [H], the agents observe the state sh ∈ S , and take actions aih ∈ Ai, i ∈ N simultaneously. All agents receive the same reward rh(sh, ah), where ah = (a1h, . . . , a N h ), and the environment transitions to the next state sh+1 ∼ Ph(·|sh, ah). We assume for simplicity that the initial state s1 of each episode is fixed. Note that the state transition here is general and not restricted to be deterministic, in contrast to Lauer & Riedmiller (2000). This makes decentralized learning considerably more challenging, as the agents cannot implicitly coordinate by enumerating/rehearsing all possible states. We focus on the decentralized setting, where each agent only observes the states, rewards, and its own actions, but not the actions of the other agents. We require that each agent is completely oblivious of the existence of the others, and is also not allowed to communicate. This decentralized information structure requires each agent to learn to make decisions based on only its local information.
For each agent i, a (Markov) policy πi : [H] × S → ∆(Ai) is a mapping from the time index and state space to a distribution over its own action space. Each agent seeks to find a policy that maximizes its own reward. A joint policy (or policy profile) π = (π1, . . . , πN ) induces a probability measure over the sequence of states and joint actions. For convenience, we use the superscript −i to denote the set of agents excluding agent i, i.e., N\{i}. For example, we can rewrite π = (πi, π−i)
2All our results in this paper immediately generalize to Markov potential games, which include Markov teams as a special case. For simplicity of notations, we illustrate our algorithm and results only in the Markov team setting. A short discussion on Markov potential games is included in Appendix B.
using this convention. For a policy profile π, and for any h ∈ [H], s ∈ S, and a ∈ A, we define the value function and the state-action value function (or Q-function) as follows:
V πh (s) def = Eπ [ H∑ h′=h rh′(sh′ , ah′) | sh = s ] , Qπh(s, a) def = Eπ [ H∑ h′=h rh′(sh′ , ah′) | sh = s, ah = a ] . (1)
For agent i, a policy πi? is a best response to π−i if V π i?,π−i h (s) = supπi V πi,π−i
h (s), for any step h ∈ [H] and state s ∈ S. A policy profile π = (πi, π−i) is a Nash equilibrium if πi is a best response to π−i for all i ∈ N . One can show that NE always exists in finite Markov teams (Boutilier, 1996). There can be multiple NE in a Markov team with possibly different values. One desirable equilibrium is the team-optimal one that achieves the highest value. Specifically, a policy profile π? is team-optimal if V π ?
h (s) = supπ V π h (s), for any step h ∈ [H] and state s ∈ S . It follows immediately
that a team-optimal policy profile is always a NE, but the converse is not true in general.
For general (non-Markov) policies, we can still define their values at step h = 1 in a sense similar to Equation (1). A best response πi? with respect to general policies π−i maximizes player i’s value at step 1, i.e., V π i?,π−i 1 (s1) = supπi V πi,π−i
1 (s1). The best response to the general policies of the opponents is not necessarily Markov. We also have an approximate notion of Nash equilibrium: Definition 1. (ε-approximate Nash equilibrium). For ε > 0, a general policy profile π = (πi, π−i) is an ε-approximate Nash equilibrium if V π i,π−i 1 (s1) ≥ supπi′ V πi ′ ,π−i
1 (s1)− ε, for all i ∈ N , where the supremum is taken over all general policies that are not necessarily Markov.
For notational convenience, in most parts of the paper we illustrate our algorithms and results for two-player teams, i.e., N = 2. It is straightforward to extend our results to the general N -player team as we defined above. With two players, we use A and B to denote the action spaces of players 1 and 2, respectively. Let S = |S|, A = |A|, B = |B|. Denote by Amax = max{A,B}, which is known to both players. We also rewrite the policies (π1, π2) as (µ, ν) for notational convenience.

heading: 3 WARM-UP: THE MATRIX TEAM
text: As a warm-up, we start with a simple matrix team problem, where two agents repeatedly interact in a single-stage team (Claus & Boutilier, 1998) with stochastic rewards. This is essentially a Markov team with a single state and H = 1. The intuitions and results in this section will become useful later in our analysis of the general Markov team problem.
A two-player finite matrix team can be described by a tuple (A,B, ρ), where A and B are the finite action spaces for agent 1 and agent 2, respectively, and the stochastic reward function ρ = (Pa,b : a ∈ A, b ∈ B) is a collection of reward distributions with support [0, 1]. The agents interact for T rounds. At each round t ∈ [T ], agent 1 (resp. agent 2) chooses an action at ∈ A (resp. bt ∈ B). Both agents receive the same reward Rt drawn i.i.d. from the reward distribution Pat,bt . We let r : A× B → [0, 1] be the expected value of the reward function, i.e., r(a, b) = ∫ 1 0 xdPa,b(x). We use µ ∈ ∆(A) and ν ∈ ∆(B) to denote the mixed strategies (or simply strategies) for agent 1 and agent 2, respectively. Slightly abusing notations, the value function for a strategy pair π = (µ, ν) is denoted by V (π) = V (µ, ν) def= Eµ×ν [r(a, b)], where a ∼ µ and b ∼ ν. In an unknown environment, a natural tactic for a self-interested agent is to adopt a no-regret learning algorithm to guide its action selection (Cesa-Bianchi & Lugosi, 2006). In our problem, we propose to let each agent run a simple no-regret learning algorithm, namely (projected) stochastic gradient descent (SGD). Our Stochastic Gradient Descent with Implicit Exploration (SGD-IX) algorithm for agent 1 is presented in Algorithm 1. The algorithm run by agent 2 is symmetric and hence not included. Our use of SGD is motivated by the simple observation that every first-order stationary point of the value function is also a NE for the matrix team (Lemma 4), and that SGD does lead to stationary points in non-convex optimization (Ghadimi & Lan, 2013). Our solution also relies on the fact that centralized SGD can be decoupled into a decentralized one (Lemma 5), as the joint strategy is a Cartesian product of individual strategies. For generality, we write our gradient descent step in Algorithm 1 as an online mirror descent (OMD) update. For more detailed descriptions and intuitions of Algorithm 1, we refer the reader to Appendix D. The following result shows that the output strategy pair of Algorithm 1 is close to a NE in expectation.
Algorithm 1: SGD-IX: Stochastic Gradient Descent with Implicit Exploration (for player 1) 1 Initialize: µ1,a ← 1/A, ∀a ∈ A, αt ← T−3/4,∀t ∈ [T ], and γ ← T−1/4. 2 for t← 1 to T do 3 Draw action at ∼ µt = (µt,1, . . . , µt,A); 4 Take action at, and observe the team reward Rt; 5 ∇̃µ,aLt ← 1−Rtµt,a+γ I{at = a},∀a ∈ A, and ∇̃µLt ← (∇̃µ,1Lt, . . . , ∇̃µ,ALt);
6 µt+1 ← argminµ∈∆(A) { 〈∇̃µLt, µ〉+ 1αtDω(µ, µt) } , where ω(µ) = 12 ‖µ‖ 2 2 , and
Dω (µ, µ ′) = ω(µ)− ω (µ′)− 〈∇ω (µ′) , µ− µ′〉;
7 Output µτ where τ is uniformly sampled from {1, . . . , T};
Theorem 1. Suppose that both agents run Algorithm 1 for T rounds. Let the two agents uniformly draw a time index τ from {1, . . . , T}, and denote by xτ = (µτ , ντ ) the output of the algorithm. For any ε > 0, if T = Ω(1/ε8) · poly(Amax), then (µτ , ντ ) is an ε-approximate Nash equilibrium in expectation, where poly(Amax) denotes a polynomial of Amax.
Note that to extract an approximate NE strategy, the agents need to sample the same time index τ . This requires some common source of randomness like a shared random seed. Such common randomness is also termed a correlation device, and is standard in decentralized learning (Bernstein et al., 2009; Arabneydi & Mahajan, 2015; Zhang et al., 2019). Note that the correlation device is never used during the learning process to coordinate the exploration, but is simply used to synchronize the selection of the strategies after they have already been generated by the learning algorithm. A common random seed is generally considered as a mild assumption and does not break the decentralized paradigm. In addition, we would like to remark that Theorem 1 is not an immediate consequence of, and, in fact, is fundamentally different from, the folklore result that the empirical frequency of plays of independent no-regret learning converges to the set of coarse correlated equilibria (CCE) in general-sum normal-form games (Hart & Mas-Colell, 2000). In teams, the set of NE can be a strict subset of CCE in general. Theorem 1 also crucially takes advantage of the first-order stationarity of NE in teams, a property that is clearly absent in the folklore result of general-sum games.
In the following theorem, we further show that the sample complexity of learning in matrix teams can be improved to Õ(1/ε4) by utilizing a variance reduction technique. The corresponding algorithm and its analysis are deferred to the appendices due to space limitations. Theorem 2. Let the agents run variance-reduced SGD (Algorithm 4 in Appendix E) for T rounds. For any ε > 0, if T = Ω̃(1/ε4) · poly(Amax), then (µτ , ντ ) is an ε-approximate NE in expectation.

heading: 4 A DECENTRALIZED COOPERATIVE MARL ALGORITHM
text: In this section, we introduce our algorithm Stage-Based V-learning with Stochastic Gradient Descent (V-learning SGD) for decentralized multi-agent RL in Markov teams. A well-known key challenge in MARL is that the environment becomes non-stationary from each agent’s own perspective when the other agents also update their policies (Busoniu et al., 2008; Zhang et al., 2021a). To address such a challenge, we propose to use a stage-based Q-learning style algorithm to create a stage-wise stationary environment for the agents, while invoking the SGD-IX subroutine (Algorithm 1) in each such stationary stage for policy updates. As an interesting side remark, stage-based Q-learning has also achieved near-optimal regret bounds in single-agent RL (Zhang et al., 2020b).
The algorithm run by agent 1 is presented in Algorithm 2. The algorithm for agent 2 is symmetric, by simply replacing the action space A with B. Each agent runs a stage-based V-learning algorithm independently, and maintains upper confidence bounds on the value functions to actively explore the unknown environment. For each step-state pair (h, s), we divide the visitations to this pair into multiple stages, where the lengths of the stages increase exponentially at a rate of (1 + 1/H) (Zhang et al., 2020b). Specifically, we let e1 = H , and ei+1 = b(1 + 1/H)eic, i ≥ 1 denote the lengths of the stages, and let the partial sums L def= { ∑j i=1 ei | j = 1, 2, 3, . . . } denote the set of ending times of the stages. For each (h, s) pair, we update our optimistic estimates V h(sh) of the value function at the end of each stage (that is, when the total number of visitations to (s, h) lies in the set
Algorithm 2: V-learning SGD: Stage-Based V-learning with Stochastic Gradient Descent 1 Initialize: V h(s)← H − h+ 1, Nh(s)← 0, Ňh(s)← 0, řh(s)← 0, v̌h(s)← 0, Ťh(s)← H ,
and µ1h(a | s)← 1/A, ∀h ∈ [H], s ∈ S, a ∈ A. 2 for episode k ← 1 to K do 3 Receive s1; 4 for step h← 1 to H do 5 Nh(sh)← Nh(sh) + 1, ň def = Ňh(sh)← Ňh(sh) + 1; 6 Take action ah ∼ µňh(· | sh), and observe reward rh and next state sh+1; 7 řh(sh)← řh(sh) + rh, v̌h(sh)← v̌h(sh) + V h+1(sh+1); 8 if Nh(sh) ∈ L then 9 //Entering a new stage
10 V h(sh)← řh(sh)ň + v̌h(sh) ň + bň, where bň ← 25H √ A3max/ň
1/4; 11 Ňh(sh)← 0, řh(sh)← 0, v̌h(sh)← 0, Ťh(sh)← ⌊ (1 + 1H )Ťh(sh) ⌋ ; 12 µ1h(a | sh)← 1/A, ∀a ∈ A; 13 αň ← ( Ťh(sh) )−3/4 , γ ← ( Ťh(sh) )−1/4 ; 14 ∇̃µ,aLň ← H−h+1−(rh+V h+1(sh+1))µňh(a|sh)+γ I{ah = a},∀a ∈ A; 15 ∇̃µLň ← (∇̃µ,1Lň, . . . , ∇̃µ,ALň); 16 µň+1h (· | sh)← argminµ∈∆(A) { 〈∇̃µLň, µ〉+ 1αňDω(µ, µ ň h(· | sh)) } , where
ω(µ) = 12 ‖µ‖ 2 2 , and Dω (µ, µ ′) = ω(µ)− ω (µ′)− 〈∇ω (µ′) , µ− µ′〉;
L), using samples only from this single stage. This way, our stage-based V-learning ensures that only the most recent O(1/H) fraction of the collected samples are used to calculate V h(sh), while the first 1−O(1/H) fraction is forgotten. Such a stage-based update framework in some sense mimics the celebrated optimistic Q-learning algorithm with a learning rate of αt = H+1H+t (Jin et al., 2018), which also roughly uses the last O(1/H) fraction of samples for value updates. However, we would like to specifically emphasize that the V-learning variant of Jin et al. (2018) does not directly apply to our problem due to some technical challenges that will be covered in Section 5. Unlike Jin et al. (2018), our stage-based V-learning ensures that the value function estimates are updated at a low frequency, a property also termed low-switching cost in single-agent RL (Bai et al., 2019), which is crucial to address the non-stationarity issue in our MARL problem.
Another key component of Algorithm 2 is the SGD subroutine. For each step h ∈ [H] and each state sh ∈ S, we treat the agents’ interactions at sh as a separate matrix team problem, similar to what we have defined in Section 3. The reward function for such a matrix team is the Q-function at the corresponding state, which is directly related to the optimistic value function estimates V h(sh) that we have constructed in the V-learning part. Since our stage-based V-learning only updates V h(sh) at the end of a stage, it creates a stage-wise stationary environment in the sense that the value of V h(sh) is fixed within each stage, which can be captured by our stationary problem formulation in Section 3. We can hence directly invoke Algorithm 1 and its analyses for the matrix team problem at each step-state pair (h, sh). Finally, we remark that both the V-learning and the SGD components of our algorithm are decentralized, which can be implemented using only the states observed and the local actions executed, without any communication or central coordination among the agents.

heading: 5 THEORETICAL ANALYSES
text: In this section, we present our main results on the sample complexity of V-learning SGD. We first introduce a few notations to facilitate the analysis. For a step h ∈ [H] of an episode k ∈ [K], we denote by skh the state that the agents observe at this time step. For any state s ∈ S, we let µkh(· | s) ∈ ∆(A) and νkh(· | s) ∈ ∆(B) be the strategies prescribed by Algorithm 2 to agents 1 and 2, respectively, at this step. Notice that such notations are well-defined for every s ∈ S even if s might not be the state skh that is actually visited. We further let µ k h = {µkh(· | s) : s ∈ S} and νkh = {νkh(· | s) : s ∈ S}. Let akh ∈ A and bkh ∈ B be the actual actions taken by the two agents.
Algorithm 3: Construction of the Output Policy µ̄kh 1 Input: The strategy trajectory {µkh} H,K h=1,k=1 specified by Algorithm 2. 2 Initialize: k′ ← k. 3 for step h′ ← h to H do 4 Receive sh′ ; 5 Take action ah′ ∼ µk ′
h′(· | sh′); 6 Uniformly sample i from {1, 2, . . . , Ňk′h′ (sh′)}; 7 Set k′ ← ľk′h′,i, where recall that ľk ′
h′,i is the index of the episode such that state sh′ was visited the i-th time (among the total Ňk ′
h′ (sh′) times) in the last stage;
For any s ∈ S, let Nkh (s), Ňkh (s), and V k
h(s) denote, respectively, the values of Nh(s), Ňh(s), and V h(s) at the beginning of the k-th episode. Note that it is proper to use the same notation to denote these values from the two agents’ perspectives, because both agents maintain the same estimates of these three terms as they can be calculated from the common observations (of the state-visitation).
Further, for a state skh, let ň k h denote the number of times that state s k h has been visited (at the h-th step) in the stage right before the current stage, and let ľkh,i denote the index of the episode that this state was visited the i-th time among the ňkh times. For notational convenience, we use ň to denote ňkh, and ľi to denote ľ k h,i, whenever h and k are clear from the context. With the new notations, the update rule in Line 10 of Algorithm 2 can be equivalently expressed as
V h(sh)← 1
ň ň∑ i=1 ( rh(sh, a ľi h , b ľi h ) + V ľi h+1(s ľi h+1) ) + bň. (2)
Similar to the situation of learning in zero-sum Markov games (Bai et al., 2020), the strategy trajectories {(µkh, νkh)} H,K h=1,k=1 themselves are not guaranteed to converge to an approximate NE. Instead, based on {(µkh, νkh)} H,K h=1,k=1, we extract an auxiliary set of policies {(µ̄kh, ν̄kh)} H,K h=1,k=1 that leads to an approximate NE. Our construction of the auxiliary policies, largely inspired by the “certified policies” (Bai et al., 2020) for learning in two-player zero-sum games, is formally presented in Algorithm 3. At the step h and state sh, µ̄kh first executes the strategy µ k h(· | sh). It then uniformly samples an episode index k′ from the total Ňkh (sh) episodes that the state sh was visited during the last stage. Finally, µ̄kh repeats a similar process at step h+ 1 and episode k
′, and so on. We define ν̄kh analogously. Compared with the certified policies in zero-sum games, our team problem has the additional challenge that the NE value is not necessarily unique, and we need to first specify which NE to compare with before we are able to define the notion of regret. Therefore, the analytical method of bounding the duality gap in zero-sum games (Bai et al., 2020) does not apply to our problem. We also remark that our construction of the auxiliary policies slightly simplifies the “certified policies” in the sense that we only need to uniformly sample an episode index at each step, while the certified policies rely on a rather complicated weighted sampling. Such simplification is a natural benefit of our stage-based V-learning that assigns uniform weights to all time steps of the same stage in history.
Clearly, the policies (µ̄kh, ν̄ k h) are non-Markov, and many notations defined for Markov policies in Section 2 no longer hold. Nonetheless, we can still define the value function and the Q-function starting from a step h ∈ [H] for a pair of non-Markov policies (µ̄kh, ν̄kh) in a similar way as in (1), because (µ̄kh, ν̄ k h) does not depend on the history before the h-th step. For notational convenience, we introduce the operators PhV (s, a, b) = Es′∼Ph(·|s,a,b)V (s′) for any value function V , and Dµ×νQ(s) = E(a,b)∼(µ,ν)Q(s, a, b). With these notations, the Bellman equations can be rewritten more succinctly as Qµ,νh (s, a, b) = ( rh + PhV µ,νh+1 ) (s, a, b), and V µ,νh (s) = (Dµh×νhQ µ,ν h ) (s) for any (s, a, b, h) ∈ S ×A× B × [H]. For each h ∈ [H], s ∈ S, define V ?,ν̄ k H+1
H+1 (s) = 0, and
V ?,ν̄kh h (s) def = 1
ň ň∑ i=1 max µ D µ×ν ľi h ( rh + PhV ?,ν̄ ľi h+1 h+1 ) (s), (3)
where the stages are partitioned in the same way as in the execution of Algorithm 2. In particular, if the visitation to a certain state s is in its first stage, we let V ?,ν̄ k h
h (s) = H − h+ 1 instead. We emphasize
that the maximum in (3) is taken with respect to each individual episode ľi, and hence V ?,ν̄kh h (s) can be considered as the value of the “dynamic best response” with respect to the non-Markov policy ν̄kh . Such a definition is stronger than the usual concept of “best fixed response in hindsight”, and is crucial to ensure that we obtain a NE in the end, instead of only a coarse correlated equilibrium (see also Remark 1 for detailed discussions). Define V µ̄ k h,?
h (s) analogously.
In the following, we start with an intermediate result, which states that the optimistic value function V k
h(s) is an upper bound of both V ?,ν̄kh h (s) and V µ̄kh,? h (s) in expectation.
Lemma 1. E[V kh(s)] ≥ E[V ?,ν̄kh h (s)] and E[V k h(s)] ≥ E[V µ̄kh,? h (s)],∀(s, h, k) ∈ S × [H]× [K].
The proof of Lemma 1 relies on our results from Section 3. Specifically, we can show that within each (stationary) stage, the agents face a matrix team problem that lasts for ň rounds, whose reward function is associated with the optimistic value function V k
h(s) at the given state s. Since in Algorithm 2, both agents essentially run an individual copy of Algorithm 1 for each state, our analysis exactly reduces to the investigation of Algorithm 1 in Section 3. We also remark that optimistic Q-learning with the celebrated learning rate of αt = H+1H+t (Jin et al., 2018) could fail at this step, because it induces a matrix team problem with weighted rewards where the weights vary over time and cannot be pre-computed. Such time-varying rewards are not compatible with our SGD framework, and as far as we know it is a challenging problem in stochastic non-convex optimization on its own. Interested readers are referred to Remarks 2 and 3 in Appendix G for detailed discussions of this technical challenge. We bypassed such a challenge by utilizing stage-based V-learning instead, which in essence assigns invariant and pre-computable weights to each time step in history. Our results hence advocate the use of stage-based V-learning in MARL, over the seminal work Jin et al. (2018), because the former creates a stage-wise stationary environment, and by nature provides a solution to the core issue of non-stationarity (Busoniu et al., 2008; Zhang et al., 2021a) in MARL.
Suppose that for every (h, k) ∈ [H] × [K], µ̄kh and ν̄kh use a common random seed to sample the episode index from {1, . . . , Ňkh (s)}. The value function for (µ̄kh, ν̄kh) can be written recursively as
V µ̄kh,ν̄ k h h (s) = 1
ň ň∑ i=1 D µ ľi h ×ν ľi h ( rh + PhV µ̄ ľi h+1 ,ν̄ ľi h+1 h+1 ) (s),
and V µ̄kH+1,ν̄ k H+1
H+1 (s) = 0. The following result shows that, in expectation (over the randomness of the reward sequence and the agents’ selected actions), the agents have no incentive to deviate from the policy pairs {(µ̄k1 , ν̄k1 )}Kk=1, up to an error term of the order O(K− 1 8 ).
Theorem 3. The auxiliary policies {(µ̄kh, ν̄kh)} H,K h=1,k=1 satisfy that
1
K K∑ k=1 E [ V ?,ν̄k1 1 (s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] ≤ O(H 17 8 S 1 8A 3 2 max/K 1 8 ), and
1
K K∑ k=1 E [ V µ̄k1 ,? 1 (s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] ≤ O(H 17 8 S 1 8A 3 2 max/K 1 8 ).
An immediate corollary is that, if we let the agents uniformly sample an episode index κ from {1, . . . ,K}, the policy pair (µ̄κ1 , ν̄κ1 ) we constructed is an approximate NE. Theorem 4. Suppose that the two agents run Algorithm 2 for K episodes with K = Ω(1/ε8) · poly(H,S,Amax), and uniformly sample an episode index κ from {1, . . . ,K}. Then, with probability at least 35 , the policy pair (µ̄ κ 1 , ν̄ κ 1 ) is an ε-approximate Nash equilibrium in expectation. A standard boosting technique (Mitzenmacher & Upfal, 2017) can increase the success probability to an arbitrary value 1− p for p ∈ (0, 25 ), by repeating the process for O(log 1 p ) times.
Remark 1. Theorem 4 states that we obtain an approximate NE with high probability for any realization of κ, rather than only in expectation over the randomness of κ. The former condition guarantees that each outcome (µ̄κ1 , ν̄ κ 1 ) of the sampling process is a Nash equilibrium with high probability; with the latter condition, we can only conclude that the uniform distribution over {(µ̄k1 , ν̄k1 )}Kk=1 is a coarse correlated equilibrium, a weaker solution concept, where the common random seed plays the role of the “trusted coordinator” in the language of correlated equilibria.
Our sample complexity bound depends on max{A,B} instead of A×B. This is the benefit from decentralized learning, and would not have been achieved by centralized approaches. Such an improvement would become more significant as the number of agents N increases. Our decentralized approach only depends on the largest single action space maxi∈N |Ai|, while the centralized methods would have an exponential dependence ΠNi=1|Ai|. Also note that our sample complexity bound in Theorem 4 holds in expectation. To obtain a standard high-probability result that holds with probability 1− p, one could either immediately apply Markov’s inequality and tolerate an additional O(1/p) factor of computation, or replace our SGD subroutine with one that has high-probability guarantees instead (Li & Orabona, 2020). We leave such improvements to our future work as they diverge from the main focus of this paper. Similar to Section 3, the sample complexity in Theorem 4 can be further improved to Õ(1/ε4) by incorporating variance reduction, which is sketched in Theorem 7 of Appendix G. We remark that we have not attempted to optimize the dependence on the other parameters H,S and Amax. There are many straightforward ways to obtain tighter bounds in these parameters, such as making the hyper-parameters αň and γ dependent on Amax. For completeness, a sample complexity lower bound in the order of Ω(1/ε2) · poly(H,S,Amax) is presented in Appendix G, which is achieved by a reduction to a single-agent RL problem. We leave the tightening of both the upper and lower bounds to our future work.
Finally, we show that our algorithm can nearly find the team-optimal NE in an important subclass of Markov teams named smooth teams. Our definition of a (λ, ρ)-smooth team, adapted from the definition of smooth games (Roughgarden, 2009; Radanovic et al., 2019), is formally introduced in Definition 5 of Appendix G. Define V ? to be the team-optimal value function, i.e., V ?h (s) = maxπ V π h (s) for any h ∈ [H], s ∈ S. The following theorem states that the output policies of Algorithm 2 converge to a λ/(1 + ρ) factor of team-optimality at a rate of O(K−1/8). Theorem 5. LetK = Ω(1/ε8)·poly(H,S,Amax). In a (λ, ρ)-smooth team, the value of the auxiliary policies {(µ̄kh, ν̄kh)} H,K h=1,k=1 satisfies
1
K K∑ k=1 E [ V µ̄k1 ,ν̄ k 1 1 (s1) ] ≥ λ 1 + ρ V ?1 (s1)− 1 1 + ρ O(K− 1 8 ).

heading: 6 SIMULATIONS
text: We empirically evaluate SGD-IX on a classic matrix team task, and V-learning SGD on a Markov team. Figure 1 illustrates the performances of the algorithms in terms of both the rewards and the L2 equilibrium gaps (which measure the L2 distance to a NE, formally defined in (39) of Appendix H). Our simulation results turn out to be more encouraging that what our theory suggests: Both the actual policy trajectories (“Current”) and the auxiliary policies (“Average”) converge to NE in the tasks we tested. Further, the rewards obtained approach those of a team-optimal oracle (“Centralized”), which suggests that our algorithm achieves the team-optimal NE very frequently in our simulations, even though our theory does not guarantee so in general. Detailed descriptions of the simulations are deferred to Appendix H due to space limitations.

heading: ETHICS STATEMENT
text: As a theory-oriented work, we do not believe that our research will cause any ethical issue, or put anyone at any disadvantage. In particular, we do not believe that our work will cause any discrimination / bias / fairness concerns, privacy and security issues, or legal compliance, and so on.

heading: REPRODUCIBILITY STATEMENT
text: For every theoretical claim in the paper, we included a clear statement and explanation of all necessary assumptions right before or after the theoretical claim. The complete proofs of all theoretical results in the paper can be found in Appendices E, F, and G. The source code used in the simulations are uploaded as supplementary materials along with the submission of the paper.

heading: Supplementary Materials for “Decentralized Cooperative

heading: Multi-Agent Reinforcement Learning with Exploration”

heading: A DETAILED DISCUSSIONS ON RELATED WORK
text: A common mathematical framework of multi-agent RL is stochastic games (Shapley, 1953), which is also referred to as Markov games. Early attempts to learn Nash equilibria in Markov games include Littman (1994; 2001); Hu & Wellman (2003); Hansen et al. (2013), but they either assume the transition kernel and rewards are known, or only yield asymptotic guarantees. More recently, various sample efficient methods have been proposed (Wei et al., 2017; Bai & Jin, 2020; Sidford et al., 2020; Xie et al., 2020; Bai et al., 2020; Liu et al., 2021; Zhao et al., 2021), mostly for learning in two-player zero-sum Markov games. Most notably, several works have investigated two-player zero-sum games in a decentralized environment: Daskalakis et al. (2020) have shown non-asymptotic convergence guarantees for independent policy gradient methods when the learning rates of the two agents follow a two-timescale rule. Tian et al. (2021) have studied online learning when the actions of the opponents are not observable, and have achieved the first sub-linear regret Õ(K 3 4 ) in the decentralized setting for K episodes. More recently, Wei et al. (2021) have proposed an Optimistic Gradient Descent Ascent algorithm with a slowly-learning critic, and have shown a strong finite-time last-iterate convergence result in the decentralized/agnostic environment. Overall, these works have mainly focused on two-player zero-sum games. These results do not carry over in any way to the decentralized cooperative/team setting that we consider in this paper.
Another line of research has considered RL in teams or cooperative games. Without enforcing a decentralized environment, Boutilier (1996) has proposed to coordinate the agents by letting them take actions in a lexicographic order. In a similar setting, Wang & Sandholm (2002) have studied optimal adaptive learning that converges to the optimal NE in Markov teams. Verbeeck et al. (2002) have presented an independent learning algorithm that achieves a Pareto optimal Nash equilibrium in common interest games with limited communication. These methods critically relied on communications among the agents (beforehand) or observing the teammates’ actions. In contrast, the distributed Q-learning algorithm in Lauer & Riedmiller (2000) is decentralized and coordinationfree, which, however, only works for deterministic tasks, and has no non-asymptotic guarantees. In fact, developing provable decentralized RL in teams with stochastic environment is considered as an open problem in Lauer & Riedmiller (2000), which we address in this paper.
Another relevant work (Arslan & Yüksel, 2016) has shown that decentralized Q-learning can converge to Nash equilibria in weakly acyclic games (WAGs), which cover Markov teams as an important special case. Their decentralized setting is most similar to ours in that each agent is completely oblivious to the presence of the others. Later, such a method has been improved in Yongacoglu et al. (2019) to achieve team-optimality. However, both of them require a coordinated exploration phase, and only yield asymptotic guarantees. In fact, to the best of our knowledge, it is not clear if their methods will converge in polynomial time. We significantly improve their results in this respect, by explicitly showing a polynomial sample complexity upper bound. Similarly, decentralized learning has also been studied in single-stage weakly acyclic games (Marden et al., 2009b) or potential games (a subclass of WAGs) (Marden et al., 2009a; Cohen et al., 2017). Two contemporaneous works (Zhang et al., 2021b; Leonardos et al., 2021) have studied gradient play in Markov potential games, which also cover Markov teams. Nonetheless, they do not consider the aspect of exploring the unknown environment, one of the key issues in RL, and/or have to assume perfect knowledge of the environment.
In general-sum normal-form games, a folklore result is that when the agents independently run noregret learning algorithms, their empirical frequency of plays converges to the set of coarse correlated equilibria (CCE) of the game (Hart & Mas-Colell, 2000). However, a CCE may suggest that the agents play obviously non-rational strategies. For example, Viossat & Zapechelnyuk (2013) have constructed an example where a CCE assigns positive probabilities only to strictly dominated strategies. On the other hand, given the PPAD completeness of finding a Nash equilibrium, convergence to NE seems hopeless in general. An impossibility result (Hart & Mas-Colell, 2003) has shown that uncoupled noregret learning does not converge to Nash equilibrium in general, due to the informational constraint
that the adjustment in an agent’s strategy does not depend on the reward functions of the others. Hence, convergence to Nash equilibria is guaranteed mostly in games with special reward structures, such as two-player zero-sum games (Freund & Schapire, 1999) and potential games (Kleinberg et al., 2009; Cohen et al., 2017).
Efficient exploration has also been widely studied in the literature of single-agent RL, see, e.g., Brafman & Tennenholtz (2002); Jaksch et al. (2010); Azar et al. (2017); Jin et al. (2018); Agarwal et al. (2020). For the tabular episodic setting, various methods (Azar et al., 2017; Zhang et al., 2020b; Menard et al., 2021) have achieved the sample complexity of Õ(H3SA/ε2), which matches the information-theoretical lower bound. When reduced to the bandit case, our problem is also related to the cooperative multi-armed bandit (MAB) problem (Lai et al., 2008; Avner & Mannor, 2014; Kalathil et al., 2014; Bubeck & Budzinski, 2020), originated from the literature of cognitive radio networks. The difference is that, in cooperative MAB, each agent is essentially interacting with an individual copy of the bandit, with an extra caution of action collisions; in our formulation, the reward function is defined on the Cartesian product of the action spaces, which allows the agents to be coupled in more general forms. A concurrent work (Chang et al., 2021) has studied cooperative multi-player multi-armed bandits with information asymmetry, which parallels our matrix team setting. Nevertheless, (Chang et al., 2021) requires stronger conditions than our decentralized setting as their algorithm relies on playing a predetermined sequence of actions.

heading: B DISCUSSIONS ON MARKOV POTENTIAL GAMES
text: An N -player episodic Markov potential game can be described by a tuple G = (N , H,S, {Ai}Ni=1, {ri}Ni=1, P ), where (1) N = {1, 2, . . . , N} is the set of agents; (2) H ∈ N+ is the number of time steps in each episode; (3) S is the finite state space; (4) Ai is the finite action space for agent i ∈ N ; (5) ri : [H] × S × A → [0, 1] is the reward function for agent i, where A = ×Ni=1Ai is the joint action space; and (6) P : [H] × S × A → ∆(S) is the transition kernel. We remark that the agents could have different reward functions, but the reward that an agent receives also depends on the joint actions of all the agents.
We define the space of policies for agent i by Πi = {πi : [H]×S → ∆(Ai)}. The joint policy space is denoted by Π = ×Ni=1Πi. For a joint policy π ∈ Π, and for any h ∈ [H], s ∈ S, and a ∈ A, we define the value function for agent i as:
V π,ih (s) def = Esh′+1∼Ph(·|sh′ ,ah′ )
[ H∑
h′=h
rih′(sh′ , ah′) | sh = s
] . (4)
We call G a Markov potential game if for every h ∈ [H] and s ∈ S, there exists a function Φsh : Π→ R, such that
Φsh(π i, π−i)− Φsh(πi ′ , π−i) = V (πi,π−i),i h (s)− V
(πi ′ ,π−i),i
h (s),
for all i ∈ N , πi, πi′ ∈ Πi and π−i ∈ Π−i. Such a definition extends the normal-form potential games (Monderer & Shapley, 1996) to state-dependent games. In words, Φsh can be considered as a potential function such that any unilateral change in an agent’s policy causes equal amount of change to the potential function and its own value function. It is known that every Markov potential game possesses at least one deterministic Nash equilibrium.
We can easily see that Markov teams are a special case of Markov potential games, where the potential function is simply the team value function. Since all the agents are correlated by a single potential function in Markov potential games, the process of learning Nash equilibria reduces to searching for the local optima of the potential function, which is exactly the core idea behind our design of Algorithm 1 for Markov teams. One can see that Lemmas 4 and 5 immediately generalize to Markov potential games, and hence the same sample complexity bounds from Theorems 1, 2 and 4 also hold in such a broader class of games.

heading: C AUXILIARY LEMMAS

heading: C.1 PROPERTIES OF THE BREGMAN DIVERGENCE
text: Let ω : X → R be a continuously differentiable mirror map that is 1-strongly convex on X , where X is a non-empty closed convex set. In this paper, we specifically take ω as the L2-induced regularizer, i.e., ω(x) = 12 ‖x‖ 2 2, although the following result holds more generally. Let Dω(x, y) = ω(x)− ω(y)− 〈x− y,∇ω(y)〉 be the ω-induced Bregman divergence. Lemma 2. (Beck, 2017, Section 9.2.1). The Bregman divergence satisfies the following properties:
(a) Three-point identity:
Dω(x, y) +Dω(y, z)−Dω(x, z) = 〈∇ω(z)−∇ω(y), x− y〉,∀x, y, z ∈ X .
(b) Suppose ϕ(x) is convex and z+ = argmin x∈X
{ ϕ(x) + 1αDω(x, z) } for some α > 0, we have
ϕ(x) + 1
α Dω(x, z) ≥ ϕ
( z+ ) + 1 α Dω ( z+, z ) + 1 α Dω ( x, z+ ) ,∀x ∈ X .

heading: C.2 PROPERTIES OF THE GRADIENT MAPPING
text: Lemma 3. (Agarwal et al., 2021, Proposition B.1). Let L : X → R be a β-smooth function. Define the gradient mapping as
Gη(x) = 1
η (ΠX (x+ η∇L(x))− x) .
The update rule for projected gradient ascent is x+ = x+ ηGη(x). If ‖Gη(x)‖2 ≤ ε, then
max x+δ∈X ,‖δ‖22≤1
δᵀ∇L(x+) ≤ ε(ηβ + 1).

heading: D AN EXTENDED DISCUSSION OF SECTION 3
text: We first introduce a few more notations for ease of presentation. Recall that r : A×B → [0, 1] is the expected value of the reward function, i.e., r(a, b) = ∫ 1 0 xdPa,b(x). We further let Q ∈ RA ×RB be the matrix of expected rewards. That is, given a mixed strategy pair π = (µ, ν), the corresponding expected reward for the players can be denoted by V (π) = V (µ, ν) = Eµ×ν [r(a, b)] = µᵀQν, where a ∼ µ and b ∼ ν. The function V is also called the value function. In certain circumstances, it is more convenient to work with costs instead of rewards. We let l : A × B → [0, 1] be the expected cost function such that l(a, b) = 1 − r(a, b),∀(a, b) ∈ A × B. Let M = 1A×B −Q be the matrix of expected costs, where 1A×B is an all-one matrix of the size A×B. Similar to the value function V , we define L(π) = L(µ, ν) def= Eµ×ν [l(a, b)] = µᵀMν as the corresponding loss function.
Our definition of an approximate Nash equilibrium (Definition 1) can be easily reduced to the matrix team setting as follows.
Definition 2. (ε-approximate NE for matrix teams). For ε > 0, a strategy pair (µ?, ν?) is an ε-approximate Nash equilibrium for a matrix team if V (µ?, ν?) ≥ V (µ, ν?)− ε, ∀µ ∈ ∆(A), and V (µ?, ν?) ≥ V (µ?, ν)− ε, ∀ν ∈ ∆(B).
In the following, we introduce the definition of an approximate stationary point.
Definition 3. (ε-approximate stationary point). A strategy pair π = (µ, ν) ∈ ∆(A) × ∆(B) is a (first-order) ε-approximate stationary point for the function V : ∆(A) ×∆(B) → R if for any δ1 ∈ RA, δ2 ∈ RB such that ‖δ1‖22 + ‖δ2‖ 2 2 ≤ 1, µ+ δ1 ∈ ∆(A), and ν + δ2 ∈ ∆(B), it holds that
δᵀ1∇µV (π) + δ ᵀ 2∇νV (π) ≤ ε.
Intuitively, π is a stationary point if the value function V cannot increase by more than ε along any direction that lies in the intersection of the strategy simplex and the neighborhood of π. Comparing the definitions of an approximate stationary point with an approximate NE easily leads to the following observation. Lemma 4. Let π? = (µ?, ν?) be an ε-approximate stationary point of the value function V for some ε > 0. Then, the strategy pair (µ?, ν?) is a √ 2ε-approximate Nash equilibrium of the corresponding matrix team.
All missing proofs of this section are deferred to Appendix E for clarity. With Lemma 4, learning a NE of the team reduces to finding a stationary point of the corresponding value function, or equivalently, the loss function. In what follows, we focus on the loss function L for reasons that will be clear soon. Numerous SGD methods have been proposed to find a stationary point in stochastic non-convex optimization (Ghadimi & Lan, 2013). For generality, we write our gradient step as an Online Mirror Descent (OMD) update with L2 regularization to search for the stationary point. The standard OMD algorithm performs the following update rule at each iteration:
xt+1 ← argmin x∈X
{ 〈∇̂L(xt), x〉+ 1
αt Dω(x, xt)
} , (5)
where X = ∆(A)×∆(B) is the feasible region, ∇̂L(xt) is an estimator of the subgradient of L at xt, αt > 0 is the stepsize, ω : X → R is a 1-strongly convex and continuously differentiable regularizer, and Dω (x, xt) = ω(x) − ω (xt) − 〈∇ω (xt) , x− xt〉 is the Bregman divergence induced by the regularizer ω.
The following nuance in the general OMD update requires some attention before we directly apply existing results. Since the optimization objective L is a function of π ∈ ∆(A) × ∆(B), existing methods typically proceed in a centralized way and manipulate π as a whole. In contrast, in our decentralized setting, the strategies µ and ν need to be updated separately and independently. However, such a nuance does not exist with SGD. Specifically, if we instantiate ω as an L2-induced regularizer, i.e., ω(x) = 12 ‖x‖ 2 2, then our OMD procedure is equivalent to projected SGD (Lattimore & Szepesvári, 2020, Example 28.2). In this case, the Bregman divergence simplifies to Dω(x, xt) = 1 2 ‖x− xt‖ 2 2. The update rule in (5) can hence be simplified to
xt+1 ← ΠX ( x− αt∇̂L(xt) ) ,
where ΠX (x) denotes the Euclidean projection of x onto X . We then point out the following simple fact, which states that conducting projected SGD on the optimization variable π = (µ, ν) as a whole is equivalent to running the same algorithm on each component of π independently.
Lemma 5. Let π = (µ, ν) be a strategy pair, and let ∇̂πL(π) = (∇̂µL(π), ∇̂νL(π)) be an arbitrary (not necessarily unbiased) gradient estimator ofL at π. Let π+ = Π∆(A)×∆(B)(π−α∇̂πL(π)) be the result of a centralized gradient update step with stepsize α > 0, and let µ+ = Π∆(A)(µ−α∇̂µL(π)) and ν+ = Π∆(B)(ν − α∇̂νL(π)) be the results of applying the gradient update to each agent’s strategy independently. Then, π+ = (µ+, ν+).
Therefore, it suffices to only analyze the centralized OMD algorithm with L2 regularization, as it is equivalent to decentralized SGD. Since our algorithm is essentially performing gradient descent instead of ascent, we choose the loss function L as our optimization objective.
We now turn to the gradient estimator ∇̂L(x). Existing stochastic non-convex optimization methods (Ghadimi & Lan, 2013) typically assume the existence of a stochastic oracle: Given any input x ∈ X , the oracle outputs a random vector ∇̂L(x, ξ)3 that is an unbiased estimator of the gradient and has bounded variance. More specifically, it is assumed that for any x ∈ X , Eξ[∇̂L(x, ξ)] ∈ ∂L(x), where ∂L(x) is the subdifferential set of L at x, and Eξ[||∇̂L(x, ξ)||22] ≤ ρ2 for some constant ρ > 0. The latter condition is also called the ρ-stochastic continuity of L (Lu, 2019). Such an oracle does not come for free in our problem, as the agents are faced with bandit feedback and can only observe
3ξ is a random variable that is used to explicitly show the randomness of the estimator. When there is no possibility of ambiguity, we suppress the ξ term and simply write the estimator as ∇̂L(x).
the reward associated with the selected action pair. To deal with this challenge, at the t-th round of the problem, we let agent 1 construct a biased gradient estimator ∇̃µLt in the form of
∇̃µ,iLt = 1−Rt µt,i + γ I{at = i},∀i ∈ A,
where µt,i is the probability that agent 1 takes action i ∈ A in round t, and γ > 0 is a constant specified in Algorithm 1. The biased gradient estimator ∇̃νLt for agent 2 is constructed similarly. We can see that the estimator we use tends to underestimate the actual cost. Such a technique is called Implicit eXploration (Neu, 2015) in the literature, and is widely used to prove high-probability regret bounds in multi-armed bandits, an objective that we do not pursue in this paper.
Our analysis of Algorithm 1 relies on showing that the biased estimator ∇̃Lt we constructed has bounded variance, and is upper bounded by the actual gradient in expectation, which satisfies the oracle requirements in a weaker sense. We can then analyze the non-asymptotic convergence behavior of our SGD algorithm to a stationary point. For generality, the proof we present in this paper follows the analysis of a general proximal point method (Davis & Drusvyatskiy, 2018; Zhang & He, 2018), by showing that the gradient of the Moreau envelope converges to 0, although simpler proofs for the convergence of SGD in non-convex optimization exist. Finally, we invoke Lemma 4 to conclude that the stationary point we find is close to a NE.

heading: E MISSING RESULTS IN SECTION 3 AND APPENDIX D

heading: E.1 PROOF OF LEMMA 4
text: Proof. Let π? = (µ?, ν?) be an ε-approximate stationary point of the value function V . For any strategy µ ∈ ∆(A),
V (µ, ν?)− V (µ?, ν?) = (µ− µ?)ᵀQν? = √
2 · 1√ 2
(µ− µ?)ᵀ∇µV (π?) ≤ √ 2ε,
where the last step follows from the ε-approximate stationarity of π? and ∥∥∥ 1√
2 (µ− µ?) ∥∥∥2 2 ≤ 1. A
similar argument can show that V (µ?, ν)− V (µ?, ν?) ≤ √
2ε,∀ν ∈ ∆(B). Therefore, (µ?, ν?) is a√ 2ε-approximate NE of the matrix team.

heading: E.2 PROOF OF LEMMA 5
text: Proof. Let π = (µ, ν). First, notice that π − α∇̂πL(π) = (µ − α∇̂Lµ(π), ν − α∇̂Lν(π)). Then, by the definition of the Euclidean projection,
π+ =Π∆(A)×∆(B) ( π − α∇̂πL(π) ) = argmin x∈∆(A)×∆(B) ∥∥∥π − α∇̂πL(π)− x∥∥∥2 2
= argmin x1∈∆(A),x2∈∆(B) (∥∥∥µ− α∇̂µL(π)− x1∥∥∥2 2 + ∥∥∥ν − α∇̂νL(π)− x2∥∥∥2 2 ) = argmin x1∈∆(A) ∥∥∥µ− α∇̂µL(π)− x1∥∥∥2 2 + argmin x2∈∆(B) ∥∥∥ν − α∇̂νL(π)− x2∥∥∥2 2
=Π∆(A) ( µ− α∇̂µL(π) ) + Π∆(B) ( ν − α∇̂νL(π) ) =µ+ + ν+.
This completes the proof of the lemma.

heading: E.3 PROOF OF LEMMA 7
text: Proof. Our proof relies on the smoothness of the function L and a gradient mapping property (Proposition B.1 in Agarwal et al. (2021), also reproduced as Lemma 3 in our Appendix C for
completeness). First, recall the definition that x̂ = argminz∈X {L(z) + λDω(z, x)}, and notice that the function L(z) + λDω(z, x) is (λ− 1)-strongly convex in z for any fixed x ∈ X . We know from Lemma 2 that L(x)− [L(x̂) + λDω(x̂, x)] ≥ (λ− 1)Dω(x, x̂). Therefore, we can derive that
L(x)− L(x̂)−Dω(x̂, x) =L(x)− (L(x̂) + λDω(x̂, x)) + (λ− 1)Dω(x̂, x) ≥(λ− 1)(Dω(x, x̂) +Dω(x̂, x)) =(λ− 1) ‖x− x̂‖22 .
It is easy to verify that L(x) is an (A+B)-smooth function, in the sense that ‖∇L(x)−∇L(y)‖2 ≤ (A+B) ‖x− y‖2. We can hence invoke Lemma 3 and conclude that x is an approximate stationary point, such that
max x+δ∈X ,‖δ‖22≤1
δᵀ∇L(x) ≤ 2(2Amax + λ) √ ε√
λ− 1 .

heading: E.4 PROOF OF THEOREM 1
text: Proof. From Lemma 5, we know that it suffices to analyze Algorithm 1 from a centralized perspective, because running projected gradient descent on the joint strategy x = (µ, ν) as a whole is equivalent to updating each agent’s strategy µ and ν individually. To show that xτ = (µτ , ντ ) is an approximate Nash equilibrium, we also know from Lemma 4 that the problem reduces to showing its approximate stationarity with respect to the loss function L.
For generality, the proof we present in this paper follows the analysis of a general proximal point method (Davis & Drusvyatskiy, 2018; Zhang & He, 2018), although many simpler proofs are given in the literature to show the convergence of SGD in non-convex optimization. To facilitate the analysis, we first introduce the following definition of the Bregman-Moreau envelope (Moreau, 1965; Drusvyatskiy & Paquette, 2019; Davis & Drusvyatskiy, 2018; Zhang & He, 2018).
Definition 4. (Bregman-Moreau envelope). For λ > 0, given a function L : X → R and a vector z ∈ X , the Bregman-Moreau envelope is defined as
Lλ(z) def = min x∈X
{ L(x) + 1
λ Dω(x, z)
} .
The corresponding Bregman proximal operator is
proxλ,L(z) def = argmin
x∈X
{ L(x) + 1
λ Dω(x, z)
} .
We start with the observation that L is a 1-weakly convex function, in the sense that L(x) + ω(x) is convex for ω(x) = 12 ‖x‖ 2 2. It then follows that for any 0 < λ < 1, the proximal operator is well-defined and unique.
Lemma 6. (Zhang & He, 2018, Lemma 2.2). Suppose L is 1-weakly convex on X , and let 0 < λ < 1. Then, for any vector z ∈ X , the function L(x) + 1λDω(x, z) is (λ
−1 − 1)-strongly convex in x, and the Bregman proximal operator proxλ,L(z) is uniquely defined.
For simplicity of notations, in the following, we use xt = (µt, νt) ∈ X to denote the strategy pair used by the two agents at round t ∈ [T ]. For any λ > 1 and x ∈ X , we let
x̂ def = prox1/λ,L(x),
and we know from Lemma 6 that x̂ is uniquely defined. By the definition of the Bregman envelope, L1/λ(xt+1) = L(x̂t+1) + λDω(x̂t+1, xt+1). The optimality of x̂t+1 implies that
L1/λ(xt+1) ≤ L(x̂t) + λDω(x̂t, xt+1). (6)
Let ∇̃Lt = (∇̃µLt, ∇̃νLt) ∈ RA+B be the biased gradient estimator constructed by concatenating the two players’ individual estimates in Line 5 of Algorithm 1. Further, let ∇L(xt) = (∇µL(xt),∇νL(xt)) = (Mνt,Mᵀµ) be the actual gradient of L at xt, a value that we do not have access to and is only used for analytical purposes. Our proof relies on a few properties of the Bregman divergence, which are standard results from the literature and are reproduced in Appendix C for completeness. Setting z = xt, z+ = xt+1, x = x̂t, α = αt, and ϕ(x) = 〈∇̃Lt, x〉 in Lemma 2(b) (Appendix C) leads to
αt 〈 ∇̃Lt, x̂t − xt+1 〉 ≥ Dω (x̂t, xt+1) +Dω (xt+1, xt)−Dω (x̂t, xt) . (7)
Rearranging and combining (6) and (7), L1/λ(xt+1) ≤ L(x̂t) + λDω(x̂t, xt+1)
≤L(x̂t) + λ ( αt 〈 ∇̃Lt, x̂t − xt+1 〉 −Dω (xt+1, xt) +Dω (x̂t, xt) ) =L1/λ(xt) + λ ( αt 〈 ∇̃Lt, x̂t − xt+1 〉 −Dω (xt+1, xt)
) =L1/λ(xt) + λαt 〈 ∇̃Lt, x̂t − xt 〉 + λ ( αt 〈 ∇̃Lt, xt − xt+1 〉 −Dω (xt+1, xt) ) , (8)
where the first equality is by the definition of the Bregman-Moreau envelope. The second term in (8) can be further decomposed into
〈∇̃Lt, x̂t − xt〉 = 〈∇L(xt), x̂t − xt〉︸ ︷︷ ︸ 1 + 〈∇̃Lt −∇L(xt), x̂t〉︸ ︷︷ ︸ 2 −〈∇̃Lt −∇L(xt), xt〉︸ ︷︷ ︸ 3 . (9)
In the following, we bound (the expected value of) each term in (9) separately. First, since L is 1-weakly convex on X , the first term in (9) is bounded by
1 = 〈∇L(xt), x̂t − xt〉 ≤ L(x̂t)− L(xt) + 1
2 ‖x̂t − xt‖22 = L(x̂t)− L(xt) +Dω(x̂t, xt). (10)
Second, by our construction of the biased estimator, for any i ∈ A, E[∇̃µ,iLt] = E [
1−Rt µt,i + γ
I{at = i} ]
≤E [
1−Rt µt,i
I{at = i} ] = ∑ j∈B P(bt = j) · E [ 1−Rt µt,i I{at = i} ]
= ∑ j∈B P(bt = j) · l(i, j) = ∇µ,iL(xt),
where here E[·] denotes the conditional expectation given the history of both agents up to the current round, and∇µ,iL(xt) denotes the i-th entry of the actual gradient ∇µL(xt). A similar result holds for agent 2, that is, E[∇̃ν,jLt] ≤ ∇ν,jL(xt),∀j ∈ B. Hence, we know that E[∇̃Lt −∇L(xt)] 0, where denotes element-wise comparison. Since x̂t lies in the the domain X = ∆(A)×∆(B), it is a non-negative vector, i.e., x̂t 0. We can therefore upper bound the expected value of the second term in (9) by
E[ 2 ] = E[〈∇̃Lt −∇L(xt), x̂t〉] ≤ 0. (11) Finally, recalling that xt = (µt, νt), we decompose the third term in (9) into two parts,
− 3 = −〈∇̃Lt −∇L(xt), xt〉 = 〈∇µL(xt)− ∇̃µLt, µt〉+ 〈∇νL(xt)− ∇̃νLt, νt〉. In the following, we show an upper bound of the first part, and the proof for the second part follows similarly. Using the definition of the inner product,
〈∇µL(xt)− ∇̃µLt, µt〉 = ∑ i∈A ( ∇µ,iL(xt)− ∇̃µ,iLt ) · µt,i (12)
where µt,i is the i-th component in the vector µt. Again, by taking the expectation of the biased gradient estimator,
E[∇̃µ,iLt] =E [
1−Rt µt,i + γ
I{at = i} ] = ∑ j∈B P(bt = j) · E [ 1−Rt µt,i + γ I{at = i} ]
= ∑ j∈B νt,j · l(i, j)µt,i µt,i + γ .
Plugging it back to (12), and recalling the definitions that ∇µL(xt) = Mνt and ∇µ,iL(xt) =∑ j∈B l(i, j)νt,j , we have that
E [ 〈∇µL(xt)− ∇̃µLt, µt〉 ] = ∑ i∈A ∑ j∈B l(i, j)νt,j − ∑ j∈B νt,j · l(i, j)µt,i µt,i + γ  · µt,i = ∑ i∈A ∑ j∈B l(i, j)νt,j ( 1− µt,i µt,i + γ ) · µt,i
= ∑ i∈A ∑ j∈B l(i, j)νt,j γ 1 + γ/µt,i
≤ ∑ i∈A ∑ j∈B l(i, j)νt,jγ
≤Aγ, where the first inequality is due to 1+γ/µt,i ≥ 1, and the second inequality uses the facts that l(i, j) ≤ 1,∀i ∈ A, j ∈ B and ∑ j∈B νt,j = 1. A similar argument shows that E [ 〈∇νL(xt)− ∇̃νLt, νt〉 ] ≤ Bγ. Combining the two parts, we conclude that the third term in (9) can be bounded by
−E[ 3 ] = −E[〈∇̃Lt −∇L(xt), xt〉] ≤ (A+B)γ ≤ 2Amaxγ. (13) From (9), we know that inequalities (10), (11), and (13) together can upper bound the expected value of the second term in (8). Finally, to bound the third term in (8), we apply Young’s inequality and obtain almost surely that
αt 〈 ∇̃Lt, xt − xt+1 〉 −Dω (xt+1, xt) = αt 〈 ∇̃Lt, xt − xt+1 〉 − 1
2 ‖xt+1 − xt‖22
≤1 2 α2t ∥∥∥∇̃Lt∥∥∥2 2 = 1 2 α2t ∑ i∈A (1−Rt)2I{at = i} (µt,i + γ)2 + ∑ j∈B (1−Rt)2I{bt = j} (νt,j + γ)2  ≤Amaxα 2 t
γ2 , (14)
where the last step holds because (µt,i + γ)2 ≥ γ2, and (1 − Rt)2I{at = i} ≤ 1 almost surely. Substituting (10), (11), (13), and (14) back to (8), we finally conclude that
E[L1/λ(xt+1)] ≤ E[L1/λ(xt) + λαt (L(x̂t)− L(xt) +Dω(x̂t, xt) + 2Amaxγ)] + Amaxλα
2 t
γ2 .
Rearranging, taking the sum from t = 1 to T , and telescoping leads to T∑ t=1 E[αt(L(xt)− L(x̂t)−Dω(x̂t, xt))]
≤ 1 λ E[L1/λ(x1)− L1/λ(xT+1)] + 2Amaxγ T∑ t=1 αt + Amax γ2 T∑ t=1 α2t
≤ 1 λ + 2Amaxγ T∑ t=1 αt + Amax γ2 T∑ t=1 α2t , (15)
where in the last step we used the fact that 0 ≤ L1/λ(x) ≤ 1,∀x ∈ X . In the following lemma, we show that bounding the LHS of (15) is sufficient to obtain an approximate stationary point.
Lemma 7. Let x̂ = prox1/λ,L(x) with λ > 1. For any ε > 0, if L(x)− L(x̂)−Dω(x̂, x) ≤ ε, then x is a 2(2Amax+λ) √ ε√
λ−1 -approximate stationary point of the value function V .
Dividing both sides of the inequality by ∑T t=1 αt, we obtain that∑T
t=1 E[αt(L(xt)− L(x̂t)−Dω(x̂t, xt))]∑T t=1 αt ≤ 1 + Amaxλγ2
∑T t=1 α 2 t
λ ∑T t=1 αt
+ 2Amaxγ. (16)
Recall that in Line 7 of Algorithm 1, the two agents uniformly sample a time index τ from {1, . . . , T} using a common random seed. Let xτ = (µτ , ντ ) be the output of the algorithm. Then, we can see that the LHS of (16) is equivalent to E[L(xτ )− L(x̂τ )−Dω(x̂τ , xτ )], because we assign equal values to the step sizes αt. Plugging back the parameter values αt = T−3/4 and γ = T−1/4, and letting λ = 2 yields,
E[L(xτ )− L(x̂τ )−Dω(x̂τ , xτ )] ≤ 4AmaxT−1/4. Together with Lemma 7, we can see that xτ is a 16 √ A3max/T
1/4-approximate stationary point of V in expectation. Finally, invoking Lemma 4, we can conclude that xτ is an 16 √ 2A3max/T
1/4- approximate Nash equilibrium in expectation. In other words, to obtain an ε-approximate NE, it suffices to use O(A12max/ε 8) samples.

heading: E.5 PROOF OF THEOREM 2
text: We first describe more explicitly the new variance-reduced algorithm that the agents are using, which has been skipped in the main text due to space limitations. Again, we only illustrate the algorithm for agent 1, as its counterpart for agent 2 is symmetric.
Similar to vanilla SGD (with no variance reduction) in Algorithm 1, agent 1 maintains a probability distribution µt over its action space for each time step t ∈ [T ]. Such a distribution is uniformly initialized, i.e., µ1,i = 1/A, ∀i ∈ A. At each time step t ∈ [T ], the agent first draws an action at according to a perturbation of the distribution µt = (µt,1, . . . , µt,A) to be specified later. It then executes the selected action at and observes the team reward Rt. Subsequently, the agent updates its probability distribution by following a variance-reduced gradient descent step.
The gradient update rule that we use and its convergence analysis are presented independently in Appendix F, as it can be considered as a solution to a standard constrained optimization problem on its own. Specifically, in Appendix F, we present a (projected) stochastic gradient descent algorithm that finds a stationary point in constrained smooth non-convex optimization at a rate of O(1/T 1/3). Our method is a straightforward extension of the STOchastic Recursive Momentum (STORM) algorithm proposed in Cutkosky & Orabona (2019), which uses a momentum-based approach to reduce the variances in SGD.
It is important to note that we rely on a non-adaptive variant of the STORM algorithm. In the original STORM (Algorithm 1 of Cutkosky & Orabona (2019)), the adaptive step size ηt is calculated in a data-dependent manner (using knowledge of ‖∇f(xt, ξt)‖). This is ordinarily a desired property, as it alleviates the need of manual parameter tuning, and also finds wide applications in heuristic methods like Adam (Kingma & Ba, 2014). However, such adaptivity breaks our decentralized paradigm, because the agents update the step sizes using the local samples, and hence can end up with different step size values. This breaks the equivalence (between centralized and decentralized SGD) that we established in Lemma 5, which critically relies on equal step size values among all agents. Consequently, we need to refer to a non-adaptive variant of STORM (Cutkosky & Orabona, 2019) that replaces the data-dependent values with a universal upper bound σ.
Finally, to prove Theorem 2, we only need to instantiate the various notations (Ft,∇ft, and ξt) of Appendix F in the context of the matrix team problem, and verify that all the assumptions in Appendix F are indeed satisfied. Then, we can complete the proof by invoking the convergence guarantee of Theorem 6 from Appendix F. It is easy to see that Ft represents the loss function L in Section 3, and the random variable ξt is associated with the action drawn at time step t. The notation ∇f deserves more commentary. To ensure that the bounded variance assumption from Appendix F is satisfied, we cannot let the agent draw an action according to the distribution µt, but instead according to a θ-greedy parameterization of µt for some θ > 0. Specifically, for a probability distribution µt at time t, the agent samples an action at according to µ̃t = (µ̃t,1, . . . , µ̃t,A), where: µ̃t,i = (1− θ)µt,i + θ/A, ∀i ∈ A. To put it in a different way, with probability 1 − θ, the agent randomly selects an action according to the distribution µt; with probability θ, the agent uniformly samples at from A. The gradient estimator is instantiated with respect to the θ-greedy parameterization as ∇f(µt, ξt) = (∇1f(µt, ξt), . . . ,∇Af(µt, ξt)), where
∇if(µt, ξt) = 1−Rt µ̃t,i I{at = i},∀i ∈ A.
Algorithm 4: STORM (Cutkosky & Orabona, 2019) with Projections 1 d1 ← ∇f(x1, ξ1); 2 for t← 1 to T do 3 ηt ← k(w+σ2t)1/3 ; 4 xt+1 ← ΠX (xt − ηtdt); 5 at+1 ← cη2t ; 6 dt+1 ← ∇f(xt+1, ξt+1) + (1− at+1)(dt −∇f(xt, ξt+1)); 7 Output xτ where τ is uniformly sampled from {1, . . . , T};
We can see that ∇f(µt, ξt) is an unbiased gradient estimator. Further, since µ̃t,i ≥ θ/A, ∀i ∈ A, we can also verify that∇f(µt, ξt) has a bounded variance. Simple calculations show that Assumption 1 in Appendix F holds under the following instantiation of parameters:
σ2 = 2Amax θ , and L = 4A2max θ2 .
We can then invoke Theorem 6 to conclude that
E [∥∥∥∥ 1ητ (x+τ+1 − xτ ) ∥∥∥∥2 ] ≤ Mw 1/3 Tk + Mσ2/3 T 2/3k ≤ C lnT T 2/3 ( Amax θ )11/3 ,
where xτ = (µτ , ντ ) and C is some absolute constant. To ensure E [∥∥∥ 1ητ (x+τ+1 − xτ )∥∥∥] ≤ ε, it
suffices that T = Ω̃ ((
A θ ) 11 2 · 1ε3 ) . If we further set θ = ε2/11, then we only need T = Ω̃ ( A11/2 ε4 ) .
From the gradient mapping property (Lemma 3), we know that ∥∥∥ 1ητ (x+τ+1 − xτ )∥∥∥ ≤ ε implies a 54εapproximate stationary point, where we used the condition that ητ ≤ 14L in Appendix F. Finally, we invoke Lemma 4 to conclude with an ε-approximate Nash equilibrium after T = Ω̃(1/ε4)·poly(Amax) rounds. This completes the proof of Theorem 2.

heading: F SGD WITH VARIANCE REDUCTION
text: In this section, we present a (projected) stochastic gradient descent method with variance reduction (Johnson & Zhang, 2013; Allen-Zhu & Hazan, 2016; Reddi et al., 2016). This method proceeds by adding a simple projection step to a non-adaptive variant of the STOchastic Recursive Momentum (STORM) algorithm proposed in Cutkosky & Orabona (2019).
We consider a generic stochastic non-convex optimization problem as follows: We are given an objective function F : RA → R, and our goal is to find a point x ∈ X ⊆ RA such that∇F (x) is close to 0, where X is the feasible region. We do not have accurate information about the function F , and can only access it through a stochastic sampling oracle f(·, ξ), where the random variable ξ represents the “randomness” of the oracle. Throughout this section, we make the following assumptions that are standard in the literature (Arjevani et al., 2019). All norms ‖·‖ in this section are L2 norms unless otherwise specified. Assumption 1. 1. We have access to a stream of random variables ξ1, . . . , ξT , such that the gra-
dient estimators are unbiased and have bounded variances: ∇Eξt [f(x, ξt)] = ∇F (x), and E[‖∇f(x, ξt)−∇F (x)‖2] ≤ σ2 for some σ > 0 for all t ∈ [T ] and x ∈ X .
2. The objective F has bounded initial sub-optimality and is L-smooth: F (x0)− infx∈X F (x) <∞, and ‖∇F (x)−∇F (y)‖ ≤ L · ‖x− y‖ ,∀x, y ∈ RA for some L > 0. The stochastic oracle also satisfies a mean-squared smoothness property for the same constant L: E[‖∇f(x, ξ)−∇f(y, ξ)‖2] ≤ L2 · ‖x− y‖2 ,∀x, y ∈ RA with probability 1.
Our method STORM with projections is formally presented in Algorithm 4. STORM uses a variant of momentum to reduce the variance of the gradients in SGD. It achieves an optimal convergence rate of O(1/T 1/3), which improves over the standard convergence rate O(1/T 1/4) of SGD with no
variance reduction. Before we present the convergence guarantee of Algorithm 4, we first introduce a few notations for ease of presentations.
For any t ∈ [T ], we break the update rule into two steps:
x̃t+1 def = xt − ηtdt, and xt+1 = ΠX (x̃t+1).
In addition, for each t ∈ [T ], we define x+t+1 def = ΠX (xt − ηt∇F (xt)) to be the next iterate updated using the full gradient ∇F (xt), a value we do not have access to. Define εt def = dt − ∇F (xt) to be the error in dt. The high-level procedure of our proof is to seek to upper bound the value
E [∑T
t=1 ∥∥∥ 1ηt (x+t+1 − xt)∥∥∥2], and then to invoke the gradient mapping property in Lemma 3 to conclude with a stationary point. This is in contrast with the unconstrained case, where Cutkosky & Orabona (2019) directly derive an upper bound of E [∑T t=1 ‖∇F (xt)‖ 2 ] . In the following, we start with a few technical lemmas.
Lemma 8. Suppose ηt ≤ 14L for all t ∈ [T ]. Then,
E[F (xt+1)− F (xt)] ≤ E [ − 3
16ηt
∥∥x+t+1 − xt∥∥2 + 7ηt8 ‖εt‖2 ] .
Proof. From the first-order optimality condition, we know that
〈x− xt+1, xt+1 − (xt − ηtdt)〉 ≥ 0,
for any x ∈ X . Taking x = xt leads to
〈xt − xt+1, xt+1 − xt〉+ 〈xt − xt+1, ηtdt〉 ≥ 0,
which in turn implies that
〈xt+1 − xt, dt〉 ≤ − 1
ηt ‖xt − xt+1‖2 . (17)
It follows that
〈∇F (xt), xt+1 − xt〉 = 〈dt − εt, xt+1 − xt〉
≤ − 1 ηt ‖xt − xt+1‖2 − 〈εt, xt+1 − xt〉
≤ − 1 ηt ‖xt − xt+1‖2 + ηt 2 ‖εt‖2 + 1 2ηt ‖xt+1 − xt‖2
=− 1 2ηt ‖xt − xt+1‖2 + ηt 2 ‖εt‖2 ,
where the first inequality uses (17), and the second inequality is due to Hölder’s inequality and Young’s inequality. From the smoothness of F ,
E[F (xt+1)] ≤E [ F (xt) + 〈∇F (xt), xt+1 − xt〉+ L
2 ‖xt+1 − xt‖2 ] ≤E [ F (xt)− 1
2ηt ‖xt − xt+1‖2 + ηt 2 ‖εt‖2 + L 2 ‖xt+1 − xt‖2 ] ≤E [ F (xt)− 3
8ηt ‖xt − xt+1‖2 + ηt 2 ‖εt‖2
] , (18)
where the last step uses ηt ≤ 14L . From the fact that ‖x+ y‖ 2 ≤ 2 ‖x‖2 + 2 ‖y‖2, we know∥∥x+t+1 − xt∥∥2 = ∥∥x+t+1 − xt+1 + xt+1 − xt∥∥2 ≤ 2∥∥x+t+1 − xt+1∥∥2 + 2 ‖xt+1 − xt‖2 .
Rearranging the terms leads to −‖xt − xt+1‖2 ≤ ∥∥x+t+1 − xt+1∥∥2 − 12 ∥∥x+t+1 − xt∥∥2 ≤‖ΠX (xt − ηt∇F (xt))−ΠX (xt − ηtdt)‖2 − 1
2 ∥∥x+t+1 − xt∥∥2 ≤‖(xt − ηt∇F (xt))− (xt − ηtdt)‖2 − 1
2 ∥∥x+t+1 − xt∥∥2 =η2t ‖εt‖
2 − 1 2 ∥∥x+t+1 − xt∥∥2 . (19) The second inequality uses the definition of x+t+1. The third step holds because the projection operator is non-expansive, i.e. ‖ΠX (x)−ΠX (y)‖ ≤ ‖x− y‖. Substituting (19) back to (18) leads to
E[F (xt+1)] ≤ E [ F (xt)− 3
16ηt
∥∥x+t+1 − xt∥∥2 + 7ηt8 ‖εt‖2 ] .
Rearranging the terms completes the proof.
Lemma 9. (Lemma 3 in Cutkosky & Orabona (2019)). For any t ∈ [T ], it holds that E [ (1− at)2
ηt−1 (∇f(xt, ξt)−∇F (xt)) · εt−1
] = 0,
E [ (1− at)2
ηt−1 (∇f(xt, ξt)−∇f(xt−1, ξt)−∇F (xt) +∇F (xt−1)) · εt−1
] = 0.
Lemma 10. (Adapted from Lemma 5 in Cutkosky & Orabona (2019)). With the notations in Algorithm 4, we have
E
[ ‖εt‖2
ηt−1
] ≤ E [ 2c2η3t−1σ 2 + (4L2η2t−1 + 1)(1− at)2
ηt−1 ‖εt−1‖2 +
4 (1− at)2 L2
ηt−1
∥∥x+t − xt−1∥∥2 ] .
Proof. First, observe that E [ ‖∇f(xt, ξt)−∇f(xt−1, ξt)−∇F (xt) +∇F (xt−1)‖2 ] ≤E [ ‖∇f(xt, ξt)−∇f(xt−1, ξt)‖2 + ‖∇F (xt)−∇F (xt−1)‖2
] − 2E [〈∇f(xt, ξt)−∇f(xt−1, ξt),∇F (xt)−∇F (xt−1)〉]
≤E [ ‖∇f(xt, ξt)−∇f(xt−1, ξt)‖2 + ‖∇F (xt)−∇F (xt−1)‖2 ] − 2E [E [〈∇f(xt, ξt)−∇f(xt−1, ξt),∇F (xt)−∇F (xt−1)〉 | ξ1, . . . , ξt−1]]
=E [ ‖∇f(xt, ξt)−∇f(xt−1, ξt)‖2 − ‖∇F (xt)−∇F (xt−1)‖2 ] ≤E [ ‖∇f(xt, ξt)−∇f(xt−1, ξt)‖2 ] . (20)
By the definition of εt, we have εt = dt−∇F (xt) = ∇f(xt, ξt)+(1−at)(dt−1−∇f(xt−1, ξt))− ∇F (xt). Therefore,
E
[ ‖εt‖2
ηt−1
] =E [ 1
ηt−1 ‖∇f(xt, ξt) + (1− at)(dt−1 −∇f(xt−1, ξt))−∇F (xt)‖2 ] =E [ 1
ηt−1 ‖at(∇f(xt, ξt)−∇F (xt)) + (1− at)(dt−1 −∇F (xt−1)) + (1− at)(∇f(xt, ξt)−∇f(xt−1, ξt)−∇F (xt) +∇F (xt−1))‖2 ]
≤E [ 2c2η3t−1 ‖∇f (xt, ξt)−∇F (xt)‖ 2 + 1
ηt−1 (1− at)2 ‖εt−1‖2
+ 2
ηt−1 (1− at)2 ‖∇f (xt, ξt)−∇f (xt−1, ξt)−∇F (xt) +∇F (xt−1)‖2
] ,
where in the last step we used Lemma 9 and the simple fact that ‖x+ y‖2 ≤ 2 ‖x‖2 +2 ‖y‖2. Further applying (20) and the assumption that E [ ‖∇f(xt, ξt)−∇F (xt)‖2 ] ≤ σ2 leads to
E
[ ‖εt‖2
ηt−1
] =E [ 2c2η3t−1σ 2 + (1− at)2
ηt−1 ‖εt−1‖2 +
2 (1− at)2
ηt−1 ‖∇f(xt, ξt)−∇f(xt−1, ξt)‖2
]
≤E [ 2c2η3t−1σ 2 + (1− at)2
ηt−1 ‖εt−1‖2 +
2 (1− at)2 L2
ηt−1 ‖xt − xt−1‖2
]
≤E [ 2c2η3t−1σ 2 + (1− at)2
ηt−1 ‖εt−1‖2 +
4 (1− at)2 L2
ηt−1
(∥∥xt − x+t ∥∥2 + ∥∥x+t − xt−1∥∥2) ]
≤E [ 2c2η3t−1σ 2 + (4L2η2t−1 + 1)(1− at)2
ηt−1 ‖εt−1‖2 +
4 (1− at)2 L2
ηt−1
∥∥x+t − xt−1∥∥2 ] .
The first inequality is due to the L-smoothness of the function f . The second inequality again uses the fact that ‖x+ y‖2 ≤ 2 ‖x‖2 + 2 ‖y‖2. The last step holds because of the non-expansiveness of the projection operator, that is,∥∥xt − x+t ∥∥ = ‖ΠX (xt−1 − ηt−1dt−1)−ΠX (xt−1 − ηt−1∇F (xt−1))‖ ≤ η2t−1 ‖εt−1‖2 . This completes the proof of the lemma.
Now, we are ready to present the convergence guarantee of Algorithm 4.
Theorem 6. (Adapted from Theorem 2 in Cutkosky & Orabona (2019)). Suppose the conditions in Assumption 1 are satisfied. For any b > 0, let k = bσ 2 3
L , c = 32L 2 + σ2/
( 7Lk3 ) =
L2 ( 32 + 1/ ( 7b3 )) , w = max ( (4Lk)3, 2σ2, ( ck 4L )3) = σ2 max ( (4b)3, 2, ( 32b+ 17b2 )3 /64 ) , and M = 16(F (x1)− F ?) + w 1/3σ2 2L2k + k3c2
L2 ln(T + 2). Then, the following convergence guarantee holds for Algorithm 4
E
[ 1
T T∑ t=1 ∥∥∥∥ 1ηt (x+t+1 − xt) ∥∥∥∥2 ] ≤ Mw 1/3 Tk + Mσ2/3 T 2/3k .
Proof. The proof is similar to that of Theorem 2 in Cutkosky & Orabona (2019), and we reproduce it here for convenience of the reader. First, define the Lyapunov function Φt = F (xt)+ 132L2ηt−1 ‖εt‖
2. From Lemma 10, we can derive that
E
[ ‖εt+1‖2
ηt − ‖εt‖
2
ηt−1
]
≤E [ 2c2η3t σ 2 + (4L2η2t + 1)(1− at+1)2
ηt ‖εt‖2 +
4 (1− at+1)2 L2
ηt
∥∥x+t+1 − xt∥∥2 − ‖εt‖2ηt−1 ]
≤E 2c2η3t σ2︸ ︷︷ ︸ At + ( (4L2η2t + 1)(1− at+1)2 ηt − 1 ηt−1 ) ‖εt‖2︸ ︷︷ ︸
Bt
+ 4 (1− at+1)2 L2
ηt ∥∥x+t+1 − xt∥∥2︸ ︷︷ ︸ Ct  . The first two terms At and Bt are exactly the same as in the proof of Theorem 2 in Cutkosky & Orabona (2019), and we refer to their results as follows:
T∑ t=1 At ≤ 2k3c2 ln(T + 2), and T∑ t=1 Bt ≤ −28L2 T∑ t=1 ηt ‖εt‖2 .
From w ≥ (4Lk)3, we know that ηt ≤ 14L . Further, since at+1 = cη 2 t , we have that at+1 ≤
ck 4Lw1/3 ≤ 1 for all t, and hence Ct ≤ 4L 2 ηt ∥∥x+t+1 − xt∥∥2. Putting it all together, we obtain 1
32L2 T∑ t=1
( ‖εt+1‖2
ηt − ‖εt‖
2
ηt−1
) ≤ k 3c2
16L2 ln(T + 2) + T∑ t=1 ( 1 8ηt ∥∥x+t+1 − xt∥∥2 − 7ηt8 ‖εt‖2 ) .
(21) From Lemma 8, we know that
E [Φt+1 − Φt] ≤ E [ − 3
16ηt ∥∥x+t+1 − xt∥∥2 + 7ηt8 ‖εt‖2 + 132L2ηt ‖εt+1‖2 − 132L2ηt−1 ‖εt‖2 ] .
Summing over t from 1 to T and then applying (21), we obtain E [ΦT+1 − Φ1] ≤ T∑ t=1 E [ − 3 16ηt ∥∥x+t+1 − xt∥∥2 + 7ηt8 ‖εt‖2 + 132L2ηt ‖εt+1‖2 − 132L2ηt−1 ‖εt‖2 ]
≤E
[ k3c2
16L2 ln(T + 2)− T∑ t=1 1 16ηt ∥∥x+t+1 − xt∥∥2 ] .
Rearranging the terms leads to
E [ T∑ t=1 1 ηt ∥∥x+t+1 − xt∥∥2 ] ≤E [ 16(Φ1 − ΦT+1) + k3c2 L2 ln(T + 2) ] ≤16(F (x1)− F ?) + 1
2L2η0 E[‖ε1‖2] +
k3c2
L2 ln(T + 2)
≤16(F (x1)− F ?) + w1/3σ2 2L2k + k3c2 L2 ln(T + 2),
where the last step holds due to the definition that η0 = kw1/3 . Since ηt is decreasing in t,
E [ T∑ t=1 1 ηt ∥∥x+t+1 − xt∥∥2 ] = E [ T∑ t=1 ηt ∥∥∥∥ 1ηt (x+t+1 − xt) ∥∥∥∥2 ] ≥ ηTE [ T∑ t=1 ∥∥∥∥ 1ηt (x+t+1 − xt) ∥∥∥∥2 ] .
Dividing both sides by TηT and recalling the definitionM = 16(F (x1)−F ?)+ w 1/3σ2 2L2k + k3c2
L2 ln(T+ 2), we obtain
E
[ 1
T T∑ t=1 ∥∥∥∥ 1ηt (x+t+1 − xt) ∥∥∥∥2 ] ≤ M TηT = M(w + σ2T )3 Tk ≤ Mw 1/3 Tk + Mσ2/3 T 2/3k ,
where in the last step we used the fact that (a+ b)1/3 ≤ a1/3 + b1/3.

heading: G MISSING RESULTS IN SECTION 5

heading: G.1 PROOF OF LEMMA 1
text: In the following, we give a proof for the first inequality. The second inequality can be proved via a similar argument. The desired result clearly holds for any state s that is in its first stage, due to our definition of V ?,ν̄ k h
h (s) for this special case. In the following, we only need to focus on the case where V h(s) has been updated at least once at the given state s before the k-th episode.
Our proof relies on induction on k ∈ [K]. First, the claim holds for k = 1 due to the aforementioned logic. For each step h ∈ [H] and s ∈ S, we consider the following two cases.
Case 1: V h(s) has just been updated in (the end of) episode k − 1. In this case,
V k h(s) = 1
ň ň∑ i=1 ( rh(s, a ľi h , b ľi h ) + V ľi h+1(s ľi h+1) ) + bň. (22)
By the definition of V ?,ν̄ k h
h (s) and the induction hypothesis,
V ?,ν̄kh h (s) = 1
ň ň∑ i=1 max µ D µ×ν ľih ( rh + PhV ?,ν̄ ľi h+1 h+1 ) (s)
≤ 1 ň ň∑ i=1 max µ D µ×ν ľih ( rh + PhV ľi h+1 ) (s), (23)
Further, we define the regret
Rň(s) = 1
ň ň∑ i=1 max µ D µ×ν ľih ( rh + PhV ľi h+1 ) (s)− 1 ň ň∑ i=1 D µ ľi h ×ν ľi h ( rh + PhV ľi h+1 ) (s). (24)
Since we only update the optimistic value function V h at the end of a stage, its value remain unchanged within each stage. Observe that in (24), the iterator ľi takes values of episode indices that belong to the same stage. Therefore, we know that V ľi h+1 is a constant value for all i ∈ {1, . . . , ň}. We can hence drop the index i and simply rewrite it as V ľ
h+1 when there is no ambiguity. With such formulation, one may quickly realize that bounding the regret Rň(s) reduces to analyzing a matrix team problem at each state s, as we have defined in Section 3. Specifically, the matrix team lasts for ň rounds. For a strategy pair (µ, ν) of the players, the expected reward function for this matrix team
is Dµ×ν(rh + PhV ľ
h+1)(s). Remark 2. It is worth mentioning that an additional step is required if one uses V-learning with the celebrated learning rate αt = H+1H+t (Jin et al., 2018) instead to update V h, which induces an update rule as follows:
V h (sh)← (1− αt)V h (sh) + αt ( rh (sh, ah, bh) + V h+1 (sh+1) + βt ) , (25)
where t is the number of times that sh has been visited, and βt is some bonus term. In this way, V h is updated every time the state sh is visited, and hence its value varies at each step. Drawing connections back to our approach, this will lead to a matrix team problem where the expected reward function is non-stationary over time. Our algorithm and analyses in Section 3 do not directly apply to such a non-stationary problem. However, in our experience, we were able to show that it is possible to carefully tune the order of the learning rate α in a way such that the level of non-stationarity in V h can be “endured” (guarantees sublinear regret in (24)) by a specifically-designed non-stationary variant of Algorithm 1. Therefore, in this step, the learning rate αt = H+1H+t requires more delicate treatments but does not completely fail.
Since in Algorithm 2, both players essentially run an individual copy of Algorithm 1 for each state, our analysis exactly reduces to the investigation of Algorithm 1 in the matrix team problem. One may also realize that taking the average over ň in (24) is equivalent to uniformly sample a time index τ from {1, . . . , ň} and then take the expectation over the randomness of such a sampling process. This is exactly how we constructed the approximate Nash equilibrium in Theorem 1. By applying the results (and proofs) of Theorem 1, we deduce that (µľτh , ν ľτ h ) is an 16H √ 2A3max/ň
1/4-approximate NE in expectation after running Algorithm 1 for ň rounds. The extra H factor comes from the additional nuance that in the problem formulation of Theorem 1, the reward is assumed to be bounded in [0, 1], while in the current matrix team problem, the rewards lie in [0, H]. The result above implies that for any µ†τ ∈ ∆(A) that is a best response of ν ľτ h , it holds that
E [ D µ†τ×ν ľτh (rh + PhV ľ h+1)(s)− Dµľτh ×ν ľτh (rh + PhV ľ h+1)(s) ] ≤ 16H √ 2A3max/ň 1/4,
where the expectation is taken over the randomness of τ , the rewards, and the algorithm itself. The above inequality is equivalent to
E[Rň] ≤ 16H √ 2A3max/ň 1/4, (26)
where the expectation is only with respect to the randomness of the rewards and the algorithm.
Remark 3. Again, we would like to discuss the alternative of using V-learning with the popular learning rate αt = H+1H+t (Jin et al., 2018). With such a learning rate, the update rule (25) of V h can be equivalently expressed as
V k
h(s) = α 0 tH + t∑ i=1 αit [ rh ( s, ak i h , b ki h ) + V ki h+1 ( sk i h+1 ) + βi ] ,
where t is the number of times that sh has been visited, ki is the index of the episode such that sh is visited the i-th time, and βi is some bonus term. The weights αit are given by
α0t = t∏ j=1 (1− αj) , and αit = αi t∏ j=i+1 (1− αj) ,∀1 ≤ i ≤ t.
If we redefine V µ̄ k h,? h (s) and V ?,ν̄kh h (s) to comply with this update rule, the regret in (24) would be rewritten as
Rt(s) = t∑ i=1 αit max µ D µ×νkih ( rh + PhV ki h+1 ) (s)− t∑ i=1 αitDµkih ×νkih ( rh + PhV ki h+1 ) (s).
For now, let us forget about the non-stationarity of V ki
h+1, and pretend that it can be considered as a constant term throughout the entire horizon by applying the method sketched in Remark 2. Even in this simplified scenario, for each fixed t, we still need to address a weighted matrix team problem where the reward at time i is assigned the weight αit. Further, as t varies, the weight α i t assigned to the same step i also changes over time. Unfortunately, such dynamically weighted rewards cannot be easily handled by Algorithm 1, as well as other stochastic non-convex optimization algorithms that we are aware of. During the execution of Algorithm 1, we are not able to modify the learning step size of a previous step online to cope with the changing weights. These weights also cannot be pre-computed, because it relies on knowing the total number of times that a certain state sh is visited during the entire horizon, which is impossible before seeing the output of the algorithm. Therefore, such an SGD framework is incompatible with the dynamic weights αit that change over t, and this technical difficulty prevents us from directly applying the popular learning rate αt = H+1H+t . Instead, we proposed to utilize stage-based V-learning, because it essentially assigns uniform weights to each time step in history. These weights only depend on the number of visitations to the state sh in a certain stage, which can be pre-computed from the recursive formula ei+1 = ⌊ (1 + 1H )ei ⌋ .
Further, let Fi be the σ-algebra generated by all the random variables before episode ľi. Then, we can see that {rh(s, aľih , b ľi h ) + V ľi h+1(s ľi h+1)}ňi=1 is a martingale with respect to {Fi}ňi=1. For any probability p ∈ (0, 1], let ι = log(2SKH/p). From the Azuma-Hoeffding inequality, it holds with probability at least 1− p/(2SHK) that
1
ň ň∑ i=1 D µ ľi h ×ν ľi h ( rh + PhV ľi h+1 ) (s)− 1 ň ň∑ i=1 ( rh(s, a ľi h , b ľi h ) + V ľi h+1(s ľi h+1) ) ≤ 2H √ ι/ň, (27)
where ι suppresses logarithmic terms. In fact, since we only need Lemma 1 to hold in expectation instead of with high probability, it suffices to replace (27) with an even simpler statement:
E
[ 1
ň ň∑ i=1 D µ ľi h ×ν ľi h ( rh + PhV ľi h+1 ) (s)− 1 ň ň∑ i=1 ( rh(s, a ľi h , b ľi h ) + V ľi h+1(s ľi h+1) )] = 0, (28)
which holds simply due to the martingale property. Finally, combining the results in (23), (24), (26), (28), we obtain that
E [ V ?,ν̄kh h (s) ] ≤E
[ 1
ň ň∑ i=1 max µ D µ×ν ľih ( rh + PhV ľi h+1 ) (s) ] ,
≤E
[ 1
ň ň∑ i=1 D µ ľi h ×ν ľi h ( rh + PhV ľi h+1 ) (s) ] + 16H √ 2A3max/ň 1/4
=E
[ 1
ň ň∑ i=1 ( rh(s, a ľi h , b ľi h ) + V ľi h+1(s ľi h+1) )] + 16H √ 2A3max/ň 1/4
≤E
[ 1
ň ň∑ i=1 ( rh(s, a ľi h , b ľi h ) + V ľi h+1(s ľi h+1) )] + bň
=E[V kh(s)],
where the second to last step is by the definition of bň = 25H √ A3max/ň
1/4. In the last step we used the formulation of V k
h(s) in (22).
Case 2: V h(s) was not updated in (the end of) episode k − 1. Since we have excluded the case that V h has never been updated, we are guaranteed that there exists an episode j such that V h(s) has been updated in the end of episode j − 1 most recently. In this case, E[V kh(s)] = E[V k−1 h (s)] = · · · = E[V jh(s)] ≥ E[V ?,ν̄jh h (s)], where the last step is by the induction hypothesis. Finally, observe that by our definition, the value of V ?,ν̄ i h
h (s) is a constant for all episode indices i that belong to the same stage. Since we know that episode j and episode k lie in the same stage, we can conclude that E[V ?,ν̄ k h
h (s)] = E[V ?,ν̄jh h (s)] ≤ E[V k h(s)].
Combining the two cases completes our proof.

heading: G.2 PROOF OF THEOREM 3
text: In the following, we present the proof for the first inequality. The second bound can be proved via a similar argument.
Before proceeding further, we first recall the definitions of several notations from Section 5 and define a few new ones. For a state skh, let n k h be the total number of episodes that this state has been visited (at the h-th step) prior to the current stage, and let lkh,i denote the index of the episode that this state was visited the i-th time among the total nkh times. Similarly, recall that ň k h denotes the number of visits to the state skh in the stage right before the current stage, and ľ k h,i denotes the i-th episode among the ňkh episodes. For simplicity, we use li and ľi to denote l k h,i and ľ k h,i, and ň to denote ň k h, whenever h and k are clear from the context.
From Lemma 1, we first obtain that K∑ k=1 E [ V ?,ν̄k1 1 (s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] ≤ K∑ k=1 E [ V k 1(s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] .
We hence only need to upper bound the RHS. For ease of exposition, we define the following notation:
δkh def = V k h(s k h)− V µ̄kh,ν̄ k h h (s k h),
The main idea of the proof is to upper bound ∑K k=1 δ k h by the next step ∑K k=1 δ k h+1, and then obtain a recursive formula. From the update rule of V k
h(s k h) in (2), we know that
V k
h(s k h) ≤ I[nkh = 0]H +
1
ň ň∑ i=1 ( rh(s, a ľi h , b ľi h ) + V ľi h+1(s ľi h+1) ) + bň,
where the I[nkh = 0] term counts for the event that the optimistic value function has never been updated for the given state.
Let Fi be the σ-algebra generated by all the random variables before the ľi-th episode. It can be seen
that {rh(s, aľih , b ľi h ) + V
µ̄ ľi h+1,ν̄ ľi h+1
h+1 (s ľi h+1)}ňi=1 is a martingale with respect to the filtration {Fi}ňi=1.
For any probability p ∈ (0, 1], let ι = log(2SKH/p). From the Azuma-Hoeffding inequality, it holds with probability at least 1− p/(2SHK) that
V µ̄kh,ν̄ k h h (s) = 1
ň ň∑ i=1 D µ ľi h ×ν ľi h ( rh + PhV µ̄ ľi h+1,ν̄ ľi h+1 h+1 ) (s)
≥ 1 ň ň∑ i=1 ( rh(s, a ľi h , b ľi h ) + V µ̄ ľi h+1,ν̄ ľi h+1 h+1 (s ľi h+1) ) − 2 √ H2ι/ň, (29)
where ι suppresses logarithmic terms. Substituting the two inequalities above into the definition of δkh, we have
δkh ≤I[nkh = 0]H + 1
ň ň∑ i=1 ( rh(s, a ľi h , b ľi h ) + V ľi h+1(s ľi h+1) ) + bň
− 1 ň ň∑ i=1 D µ ľi h ×ν ľi h ( rh + PhV µ̄ ľi h+1,ν̄ ľi h+1 h+1 ) (s)
≤I[nkh = 0]H + 1
ň ň∑ i=1 ( V ľi h+1(s ľi h+1)− V µ̄ ľi h+1,ν̄ ľi h+1 h+1 (s ľi h+1) ) + bň + 2 √ H2ι/ň
=I[nkh = 0]H + 1
ň ň∑ i=1 δľih+1 + bň + 2 √ H2ι/ň, (30)
where in the last step we used the definition of δľih+1. In fact, since we only need the results in Theorem 3 to hold in expectation rather than with high probability, we can simply replace (29) with
E [ V µ̄kh,ν̄ k h h (s) ] = 1
ň ň∑ i=1 E [ rh(s, a ľi h , b ľi h ) + V µ̄ ľi h+1,ν̄ ľi h+1 h+1 (s ľi h+1) ] ,
by invoking the martingale property. Consequently, we can also derive an expectation form of (30) as follows:
E [ δkh ] ≤ E [ I[nkh = 0]H + 1
ň ň∑ i=1 δľih+1 + bň
] . (31)
To find an upper bound of ∑K k=1 δ k h, we proceed to upper bound each term in the RHS of (31)
separately. First, notice that ∑K k=1 I [ nkh = 0 ] ≤ SH , because each fixed state-step pair (s, h)
contributes at most 1 to ∑K k=1 I [ nkh = 0 ] . Next, we turn to analyze the second term in the RHS of (31). Observe that
K∑ k=1 1 ňkh ňkh∑ i=1 δ ľkh,i h+1 = K∑ k=1 K∑ j=1 1 ňkh δjh+1 ňkh∑ i=1 1 [ ľkh,i = j ]
= K∑ j=1 δjh+1 K∑ k=1 1 ňkh ňkh∑ i=1 1 [ ľkh,i = j ] . (32)
For a fixed episode j, notice that ∑ňkh i=1 1[ľ k h,i = j] ≤ 1, and that ∑ňkh i=1 1[ľ k h,i = j] = 1 happens if and only if skh = s j h and (j, h) lies in the previous stage of (k, h) with respect to the state-step pair (skh, h). Define Kj def = {k ∈ [K] : ∑ňkh i=1 1[ľ k h,i = j] = 1}. We then know that all episode indices k ∈ Kj belong to the same stage, and hence these episodes have the same value of ňkh. That is, there
exists an integer Nj > 0, such that ňkh = Nj ,∀k ∈ Kj . Further, since the stages are partitioned in a way such that each stage is at most (1 + 1H ) times longer than the previous stage, we know that |Kj | ≤ (1 + 1H )Nj . Therefore, for every j, it holds that
K∑ k=1 1 ňkh ňkh∑ i=1 1 [ ľkh,i = j ] ≤ 1 + 1 H . (33)
Combining (32) and (33) leads to the following upper bound of the second term in (31):
K∑ k=1 1 ňkh ňkh∑ i=1 δ ľkh,i h+1 ≤ (1 + 1 H ) K∑ k=1 δkh+1. (34)
So far, we have obtained the following upper bound:
E [ K∑ k=1 δkh ] ≤ SH2 + (1 + 1 H )E [ K∑ k=1 δkh+1 ] + K∑ k=1 bňkh .
Iterating the above inequality over h = H,H − 1, . . . , 1 leads to
E [ K∑ k=1 δk1 ] ≤ O ( SH3 + H∑ h=1 K∑ k=1 (1 + 1 H )h−1bňkh ) , (35)
where we used the fact that (1 + 1H ) H ≤ e. In the following, we analyze the bonus term bňkh more carefully. Recall our definitions that e1 = H, ei+1 = ⌊ (1 + 1H )ei ⌋ , i ≥ 1, and bň =
25H √ A3max/ň 1/4. For any h ∈ [H],
K∑ k=1 (1 + 1 H )h−1bňkh ≤ K∑ k=1 (1 + 1 H )h−125H √√√√ A3max( Ňkh (s k h) )1/4
=25H √ A3max ∑ s∈S ∑ j≥1 (1 + 1 H )h−1e − 18 j K∑ k=1 I [ skh = s, Ň k h (s k h) = ej ] =25H √ A3max
∑ s∈S ∑ j≥1 (1 + 1 H )h−1w(s, j)e − 18 j ,
where we define w(s, j) def= ∑K k=1 I [ skh = s, Ň k h (s k h) = ej ] for any s ∈ S. If we further let w(s) def=∑
j≥1 w(s, j), we can see that ∑ s∈S w(s) = K. For each fixed state s, we now seek an upper bound of its corresponding j value, denoted as J in what follows. Since each stage is (1 + 1H ) times longer than its previous stage, we know that w(s, j) = ∑K k=1 I [ skh = s, Ň k h (s k h) = ej ] = ⌊ (1 + 1H )ej ⌋ for
any 1 ≤ j ≤ J . Since ∑J j=1 w(s, j) = w(s), we obtain that eJ ≤ (1 + 1 H )
J−1 ≤ 10 1+ 1H w(s) H by taking the sum of a geometric sequence. Therefore, by plugging in w(s, j) = ⌊ (1 + 1H )ej ⌋ ,
∑ j≥1 (1 + 1 H )h−1w(s, j)e − 18 j ≤ O  J∑ j=1 e 7 8 j  ≤ O (w(s) 78H 18) , where in the second step we again used the formula of the sum of a geometric sequence. Finally, using the fact that ∑ s∈S w(s) = K and applying the Cauchy-Schwartz inequality, we have
H∑ h=1 K∑ k=1 (1 + 1 H )h−1bňkh =O H2√A3max∑ s∈S ∑ j≥1 (1 + 1 H )h−1w(s, j)e − 18 j  ≤O ( K 7 8H 17 8 S 1 8A 3 2 max ) . (36)
Summarizing the results above leads to
E [ K∑ k=1 δk1 ] ≤ O ( SH3 +K 7 8H 17 8 S 1 8A 3 2 max ) .
When K is large enough such that K ≥ SH3/7, the second term becomes dominant, and we obtain the desired result:
K∑ k=1 E [ V ?,ν̄k1 1 (s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] ≤ K∑ k=1 E[δk1 ] ≤ O ( K 7 8H 17 8 S 1 8A 3 2 max ) .
This completes the proof of the theorem.

heading: G.3 PROOF OF THEOREM 4
text: The result essentially follows from the “regret to PAC guarantee” conversion given in Jin et al. (2018). Specifically, we know from Theorem 3 that
K∑ k=1 E [ V ?,ν̄k1 1 (s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] ≤ C1K 7 8H 17 8 S 1 8A 3 2 max (37)
for some C1 that is independent of K. It is important to notice that E [ V ?,ν̄k1 1 (s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] ≥ 0, (38)
for each individual k ∈ [K], because V ?,ν̄ k 1
1 (s1) is defined as the dynamic best response instead of the best fixed response in hindsight. Since the two agents uniformly draw the index κ from {1, 2, . . . ,K}, we can conclude from (37) and (38) that for any realization of κ,
E [ V ?,ν̄κ1 1 (s1)− V µ̄κ1 ,ν̄ κ 1 1 (s1) ] ≤ 5C1H 17 8 S 1 8A 3 2 max/K 1 8 ,
with probability at least 45 . To ensure 5C1H 17 8 S 1 8A 3 2 max/K 1 8 ≤ ε, it suffices to require that K ≥ 58C81H 17SA12max/ε 8. Similarly, if K ≥ 58C82H17SA12max/ε8 for some constant C2, then
E [ V µ̄κ1 ,? 1 (s1)− V µ̄κ1 ,ν̄ κ 1 1 (s1) ] ≤ 5C2H 17 8 S 1 8A 3 2 max/K 1 8 .
Therefore, if the agents run Algorithm 2 for at least K = Ω(H17SA12max/ε 8) episodes, the randomly drawn policy pair (µ̄κ1 , ν̄ κ 1 ) will be an ε-approximate NE in expectation with probability at least 3 5 by taking a union bound. To improve the success probability to an arbitrary value 1−p for p ∈ (0, 25 ), one can adopt the standard “boosting” technique in the literature of randomized algorithms (Mitzenmacher & Upfal, 2017) , by repeating the process for O(log 1p ) iterations, and then choosing the iteration with the smallest resulting regret. Finally, we remark that we have made no effort to optimize the sample complexity dependence on H,S, and Amax, which can certainly be improved.
G.4 V-LEARNING SGD WITH VARIANCE REDUCTION
Theorem 7. (Informal). For any ε > 0, replace the SGD-IX subroutine in Algorithm 2 with a variance-reduced SGD method (Algorithm 4), and set the bonus term to be bň = CHA6/11max (ln ň) 1/2
ň1/3ε1/3
for some absolute constant C instead. Suppose that the two agents run such an algorithm for K episodes with K = Ω̃(1/ε4) · poly(H,S,Amax), and uniformly sample an episode index κ from the set {1, . . . ,K}. Then, with probability at least 35 , the policy pair (µ̄ κ 1 , ν̄ κ 1 ) is an ε-approximate Nash equilibrium in expectation.
Proof sketch. The proof follows a similar procedure as in the proofs of Lemma 1 and Theorems 3 and 4. Due to their similarities, in the following, we only sketch some of the differences in the proofs. By setting the new bonus term to be bň = CHA6/11max (ln ň) 1/2
ň1/3ε1/3 , one can verify that the proof of Lemma 1
still holds. That is, the optimistic value function V k
h(s) (with a modified bonus term) is still a valid
upper bound of both V ?,ν̄ k h h (s) and V µ̄kh,? h (s) in expectation. Everything up to (35) in the proof of Theorem 3 also holds. Starting from (35), we further deduce that
H∑ h=1 K∑ k=1 (1 + 1 H )h−1bňkh ≤ Õ ( K 2 3 /ε 1 3 ) · poly(H,S,Amax),
by following a similar argument as when we derived (36). Plugging the above inequality back to the proof of Theorem 3, we obtain that
1
K K∑ k=1 E [ V ?,ν̄k1 1 (s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] ≤ Õ(K− 13 ε− 13 ) · poly(H,S,Amax), and
1
K K∑ k=1 E [ V µ̄k1 ,? 1 (s1)− V µ̄k1 ,ν̄ k 1 1 (s1) ] ≤ Õ(K− 13 ε− 13 ) · poly(H,S,Amax).
Finally, repeating the proof idea of Theorem 4 leads to the desired result.

heading: G.5 A SAMPLE COMPLEXITY LOWER BOUND
text: To obtain a sample complexity lower bound of the problem, we consider a simple instance where either A or B is a singleton, i.e., A = 1 or B = 1. Learning an approximate NE in such a Markov team reduces to finding a near-optimal policy in a single-agent RL problem. Applying the regret lower bound of single-agent RL yields the following result for RL in Markov teams.
Corollary 1. (Corollary of Jaksch et al. (2010)). For any algorithm, there exists a two-player Markov team problem that takes the algorithm at least Ω(H3SAmax/ε2) episodes to learn an ε-approximate Nash equilibrium.
We remark that such a lower bound might be very loose. Reducing to a single-agent RL problem evades the strategic learning behavior of the agents and the non-stationarity that such behavior causes to the environment, which in our opinion are the central difficulties of learning in Markov teams. To derive a tighter lower bound in our decentralized setting, one should also utilize the additional constraint that each agent only has access to its local information, a factor that Corollary 1 apparently does not take into account. It is hence unsurprising that when comparing Theorem 4 with Corollary 1, we can see an obvious gap in the parameter dependence.

heading: G.6 DECENTRALIZED MARL IN SMOOTH TEAMS
text: In the following, we address an important subclass of Markov teams, named smooth teams, using the V-learning SGD algorithm. More importantly, we show that our algorithm can achieve near team-optimality, i.e., find the best Nash equilibrium, in such problems.
Smooth games were first introduced in Roughgarden (2009) to study the Price of Anarchy (POA) in normal-form games. A large class of games are covered as examples of smooth games, including congestion games and many forms of auctions (Roughgarden, 2009; Syrgkanis & Tardos, 2013). The notion of smoothness was later extended to learning in normal-form games (Syrgkanis et al., 2015; Foster et al., 2016) and cooperative Markov games (Radanovic et al., 2019; Mao et al., 2020). This concept essentially ensures that the game has a bounded POA, and hence decentralized no-regret learning dynamics can possibly converge to near-optimality.
Define V ? to be the team-optimal value function, i.e., V ?h (s) = maxπ∈Π V π h (s) for any h ∈ [H], s ∈ S , where Π is the set of joint policies of the agents. We consider the following definition of a smooth Markov team:
Definition 5. (Adapted from Radanovic et al. (2019)). For λ ≥ 0 and 0 < ρ < 1, an N -player Markov team is (λ, ρ)-smooth if there exists a policy profile π? such that for any policy profile π = (πi, π−i):
V π i?,π−i
h (s) ≥ λ · V ? h (s)− ρ · V πh (s),∀i ∈ N , s ∈ S, h ∈ [H].
The (λ, ρ)-smoothness ensures that agent i continues doing well by playing its optimal policy even when the other agents are using slightly sub-optimal policies. It immediately follows that the resultant policies of Algorithm 2 converge to a λ/(1 +ρ) factor of the team-optimal NE at a rate of O(K−1/8).
Theorem 5. Let K = Ω(1/ε8) · poly(H,S,Amax). In a (λ, ρ)-smooth team, the value of running the auxiliary policies {(µ̄kh, ν̄kh)} H,K h=1,k=1 satisfies
1
K K∑ k=1 E [ V µ̄k1 ,ν̄ k 1 1 (s1) ] ≥ λ 1 + ρ V ?1 (s1)− 1 1 + ρ O(K− 1 8 ).
Proof. From Theorem 3, it holds that
1
K K∑ k=1 E [ V µ̄k1 ,ν̄ k 1 1 (s1) ] ≥ 1 K K∑ k=1 E [ V ?,ν̄k1 1 (s1) ] −O(K− 18 )
≥ 1 K K∑ k=1 ( λ · V ?1 (s1)− ρ · E [ V µ̄k1 ,ν̄ k 1 1 (s1) ]) −O(K− 18 ),
where the second step is simply by the definition of smoothness. Rearranging the terms leads to the desired result.
Equivalently, such a result also holds when the agents uniformly sample an episode index κ from {1, . . . ,K} and then run (µ̄κh, ν̄κh). In Theorem 5, we illustrated our results in the case of two agents, although one can easily verify that a similar result holds more generally for N -player smooth teams. Our approach significantly generalizes Radanovic et al. (2019) and Mao et al. (2020) in that we design natural update rules for both agents, who play symmetric roles in the self-play setting; the other two works only assign the algorithm to one agent, and have to assume that the policy of the other agent changes slowly.

heading: H SIMULATIONS
text: In this section, we demonstrate the empirical performances of our algorithms. We start by evaluating SGD-IX (Algorithm 1) on a classic matrix team problem. Then, we evaluate V-learning SGD (Algorithm 2) on a Markov team problem, and compare its performances with various benchmarks.

heading: H.1 MATRIX TEAMS
text: First, to evaluate Algorithm 1, we use a classic matrix team example from the literature (Claus & Boutilier, 1998; Lauer & Riedmiller, 2000). Its reward table is reproduced in Table 1, where agent 1 is the row player, and agent 2 is the column player, both being maximizers. The action spaces of the agents are A = {a0, a1, a2} and B = {b0, b1, b2}. There are three deterministic Nash equilibria in this team, among which two of them, (a0, b0) and (a2, b2), are team-optimal. It would be preferred that the agents not only learn a NE, but also settle on the same NE out of the two team-optimal ones.
We run our SGD-IX algorithm on this task for T = 5000 rounds, and we set the step size αt = 0.002 and the implicit exploration parameter γ = 0.001. We evaluate our algorithm in terms of both the rewards it obtained and its L2 equilibrium gap. Specifically, we define the L2 equilibrium gap as the L2 distance to a equilibrium point. For a pair of strategies (µ, ν) ∈ ∆(A)×∆(B), its L2 equilibrium gap is defined as:
Gap(µ, ν) def= ∥∥µ− µ†(ν)∥∥2 2 + ∥∥ν − ν†(µ)∥∥2 2 , (39)
where ν†(µ) (resp. µ†(ν)) is the best response with respect to µ (resp. ν), and ‖·‖2 is the L2 norm. The simulation results are presented in Figure 2. All results are averaged over 20 runs. Notice that we evaluate two sets of strategy trajectories: The “Current” strategy (µt, νt) is the strategy pair used by Algorithm 1 at round t, while the “Average” strategy (µ̄t, ν̄t) is the uniformly sampled strategy pair we obtain after running Algorithm 1 for t rounds. Drawing an analogy to the Markov team setting, “Current” denotes the actual policies {(µkh, νkh)} H,K h=1,k=1 used at the step (k, h) of Algorithm 2, while “Average” represents the auxiliary policies {(µ̄kh, ν̄kh)} H,K h=1,k=1 that we have constructed. In matrix teams, the auxiliary policy at round t simply reduces to uniformly drawing a random time index τ from {1, . . . , t} and running the strategy pair (µτ , ντ ).
Note that in Theorems 1 and 4 we only have theoretical guarantees for the auxiliary policies (“Average”) but not the actual policy trajectories (“Current”). Encouragingly, our simulation results show that the actual policy trajectories also enjoy convergence behavior to NE empirically, and in most cases the “Current” policies converge even faster. Specifically, from Figure 2(a), we can see that the equilibrium gap of both “Current” and “Average” converge to zero, indicating that they indeed find an equilibrium as the number of iterations increase. The convergence of “Average” slightly lags behind “Current” because in Algorithm 1, “Average” simply takes the mean value over the actual trajectories, which requires some time to reflect the convergence behavior. A more promising result is that from Figure 2(b), we can see that the rewards collected by “Current” and “Average” converge to values close to 9. This suggests that Algorithm 1 not only finds a NE in this task, but actually converges to a team-optimal equilibrium most of the time. It does not exactly reach the team-optimal value of 10 because it still converges to non-team-optimal NE at a rather low frequency.

heading: H.2 MARKOV TEAMS
text: We further evaluate Algorithm 2 on a Markov team task. Inspired by Yongacoglu et al. (2019), we construct a Markov team problem with two states S = {s0, s1}, where s0 is the “good state” and s1 is the “bad state”. Each agent has two candidate actions A = {a0, a1} or B = {b0, b1}. The reward function at each state is presented in Table 2. Specifically, at state s1, both agents get a reward of 0 no matter what actions they select, while at state s0, they will obtain a strictly positive reward if they either take the joint action (a0, b1) or the one (a1, b0). The state transition function is defined as follows:
Ph(s0 | s0 or s1, a0, b1) = 1− ε, Ph(s1 | s0 or s1, not (a0, b1)) = 1− ε, ∀h ∈ [H],
and all the other transitions happen with probability ε. Intuitively, no matter which state the agents are in, they will transition to the good state s0 with a high probability 1− ε at the next step as long as they select the action pair (a0, b1). All the other joint actions will lead to the bad state s1 with a high probability 1− ε. The task hence rewards the agents who learn to consistently play the action pair (a0, b1).
We run our V-learning SGD algorithm on this example for K = 5000 episodes, each episode containing H = 10 steps. We set the transition probability ε = 0.1, the step size αt = 0.0004, and the implicit exploration parameter γ = 0.002. We again evaluate our algorithm in terms of the reward and the equilibrium gap. The simulation results are presented in Figure 3, where all results are averaged over 20 runs. We would like to remind the reader that “Current” denotes the performances of the actual policy trajectory {(µkh, νkh)} H,K h=1,k=1, while “Average” represents the auxiliary policies {(µ̄kh, ν̄kh)} H,K h=1,k=1.
h=1,k=1, while “Average”
represents the auxiliary policies {(µ̄kh, ν̄kh)} H,K h=1,k=1. Shaded areas denote the standard deviations of the equilibrium gap or reward.
We compare our algorithm with two meaningful benchmarks. The first benchmark is a “Centralized” oracle. This oracle acts as a centralized coordinator that can control the actions of both agents. Such an oracle essentially converts the Markov team task into a single-agent RL problem. The (randomized) action space of the centralized agent is ∆(A × B), which is larger than the ∆(A) × ∆(B) space that we allow in our decentralized approach. Therefore, “Centralized” clearly upper bounds the performance that any decentralized learning algorithm can possibly achieve in this task. In our simulations, we implement “Centralized” by using a Hoeffding-based variant of a state-of-the-art single-agent RL algorithm UCB-ADVANTAGE (Zhang et al., 2020a). This algorithm has achieved the tightest sample complexity bound for single-agent RL in theory, and has also demonstrated remarkable empirical performances in practice (Mao et al., 2020). Such an algorithm could provide a strong performance upper bound in our task. The second benchmark we consider is the naïve “Independent” Q-learning. Specifically, we let each agent run a single-agent Q-learning algorithm independently, without being aware of the existence of the other agent or the structure of the Markov team. Each agent maintains an local optimistic Q-function, and takes greedy actions with respect to such optimistic estimates, without taking into account the other agents’ actions. Since the agents update their policies simultaneously, the stationarity assumption of the environment in single-agent RL quickly collapses, and the theoretical guarantees for single-agent Q-learning no longer hold. This is also reminiscent of the “independent learner” approach proposed in an early work (Claus & Boutilier, 1998) for learning in Markov teams. We believe such a benchmark could provide meaningful intuitions about the consequences of not taking care of the team structure in decentralized methods. In our simulations, we implement such a benchmark by letting each agent running a variant of the single-agent UCB-ADVANTAGE (Zhang et al., 2020a) algorithm independently, where the (randomized) action spaces of the agents are ∆(A) and ∆(B). We can observe from Figure 3(a) that the L2 equilibrium gaps of both “Current” and “Average” converge to zero, indicating that both of them find a Nash equilibrium. This is encouraging because our theoretical results in Theorem 4 only guarantee the performance of the auxiliary policies, while our simulation results suggest that the actual policy trajectories also converge to a NE. The actual policy sequence also converges faster than the auxiliary policies, for reasons similar to the matrix team simulations. From Figure 3(b), we can again see that our algorithm even finds to the team-optimal NE most of the time, as the rewards it obtained nearly match the “Centralized” oracle. On the other hand,
the “Independent” benchmark converges, albeit faster, to a clearly suboptimal value. This reiterates that the naïve idea of independent learning does not work well for multi-agent RL in general, and a careful treatment of the game/team structure (like our SGD subroutine) is necessary. Finally, the implemented V-learning SGD algorithm takes much fewer samples to find an approximate NE than our theoretical results suggested. This indicates that the theoretical bounds might be overly conservative, and our algorithm could converge much faster in practice.

Review summaries:
The authors establish a variety of novel theoretical results for Markov team problems. First of all, they present a novel non-asymptotic convergence bound for single-step Markov teams, using a no-regret independent learning variant of stochastic gradient descent (SGD-IX). Secondly, they extend SGD-IX with a stage-based design to arrive at a novel, decentralized MARL algorithm (Stage-based V-learning). For Stage-based V-learning, the authors additionally proceed to establish a novel non-asymptotic convergence bounds that do not depend on the size of the centralised joint policy, but instead just the largest individual action space. The authors also show how their algorithm sidesteps technical difficulties inherent to establishing convergence bounds for Optimistic Q-learning [Jin et al, 2018]. Using a simple matrix team task, as well as a Markov team task as an example, the authors demonstrate empirically that their algorithms can end up selecting favourable Nash equilibria.

Authors tackles the cooperative multi-agent reinforcement learning problem with stochastic dynamics. They introduced V-learning Stochastic Gradient Descent (SGD) with theoretical proofs of convergence to the optimal policy with desired e probability. The resulting sample complexity is O(1/e^4) * Poly (A,S,H), where H is the episode horizon. Authors also empirically tested the algorithm in a toy 2 agent matrix game showing convergence to the optimal policy.

This paper studies sample complexity bound for finding approximate Nash equilibrium in two-player Markov cooperative games using decentralized algorithms. The authors designed a stage-based V-learning algorithm, which achieves eps approximate Nash equilibrium in $1/eps^8$ episodes. Under the smoothness conditions, such equilibrium is a team-optimal Nash.