Model-Based Reinforcement Learning for Offline Zero-Sum Markov Games

1. Introduction

Multiagent reinforcement learning (MARL), a subfield of reinforcement learning (RL) that involves multiple individuals interacting/competing with each other in a shared environment, has garnered widespread recent interest, partly sparked by its capability of achieving superhuman performance in game playing and autonomous driving (Shalev-Shwartz et al. 2016, Baker et al. 2019, Berner et al. 2019, Brown and Sandholm 2019, Jaderberg et al. 2019, Vinyals et al. 2019). The coexistence of multiple players—whose own welfare might come at the expense of other parties involved—makes MARL inherently more intricate than the single-agent counterpart.

A standard framework to describe the environment and dynamics in competitive MARL is Markov games (MGs), which are generally attributed to Shapley (1953) (originally referred to as stochastic games). Given the conflicting needs of the players, a standard goal in Markov games is to seek some sort of steady-state solutions with the Nash equilibrium (NE) being arguably the most prominent one. Whereas computational intractability has been observed when calculating NEs in general-sum MGs and/or MGs with more than two players (Daskalakis et al. 2009, Daskalakis 2013), an assortment of tractable algorithms have been put forward to solve two-player, zero-sum Markov games. On this front, a large strand of recent works revolves around developing sample- and computation-efficient paradigms (Bai et al. 2020, Xie et al. 2020, Zhang et al. 2020a, Liu et al. 2021, Tian et al. 2021). What is particularly noteworthy here is the recent progress in overcoming the so-called curse of multiple agents (Bai et al. 2020, Jin et al. 2021a, Li et al. 2022b); that is, although the total number of joint actions exhibits exponential scaling in the number of agents, learnability of Nash equilibria becomes plausible even when the sample size scales only linearly with the maximum cardinality of the individual action spaces. See also Jin et al. (2021a), Song et al. (2021), Mao and Başar (2022), and Daskalakis et al. (2023) for similar accomplishments in learning coarse correlated equilibria in multiplayer general-sum MGs.

The aforementioned works permit online data collection either via active exploration of the environment or through sampling access to a simulator. Nevertheless, the fact that real-time data acquisition might be unaffordable or unavailable—for example, it could be time-consuming, costly, and/or unsafe in healthcare and autonomous driving—constitutes a major hurdle for widespread adoption of these online algorithms. This practical consideration inspires a recent flurry of studies collectively referred to as offline RL or batch RL (Kumar et al. 2020, Levine et al. 2020) with the aim of learning based on a historical data set of logged interactions.

1.1. Data Coverage for Offline Markov Games

The feasibility and efficiency of offline RL are largely governed by the coverage of the offline data in hand. On one hand, if the available data set covers all state–action pairs adequately, then there is sufficient information to guarantee learnability; on the other hand, full data coverage imposes an overly stringent requirement that is rarely fulfilled in practice and is oftentimes wasteful in terms of data efficiency. Consequently, a recurring theme in offline RL gravitates around the quest for algorithms that work under minimal data coverage. Encouragingly, the recent advancement on this frontier (e.g., Rashidinejad et al. 2021, Xie et al. 2021) uncovers the sufficiency of single-policy data coverage in single-agent RL; namely, offline RL becomes information-theoretically feasible as soon as the historical data covers the part of the state–action space reachable by a single target policy.

Unfortunately, single-policy coverage is provably insufficient when it comes to Markov games with negative evidence observed in Cui and Du (2022b). Instead, a sort of unilateral coverage—that is, a condition that requires the data to cover not only the target policy pair but also any unilateral deviation from it—seems necessary to ensure learnability of Nash equilibria in two-player, zero-sum MGs. Employing the so-called unilateral concentrability coefficient $C^{\star }$ to quantify such unilateral coverage as well as the degree of distribution shift (which we define shortly in Assumption 1), Cui and Du (2022b) demonstrate how to find ε-Nash solutions in a finite-horizon, two-player, zero-sum MG once the number of sample rollouts exceeds

Here, S is the number of shared states; A and B represent, respectively, the number of actions of the max-player and the min-player; H stands for the horizon length; and the notation $\tilde{O}(·)$ denotes the order-wise scaling with all logarithmic dependency hidden.

Despite being an intriguing polynomial sample complexity bound, a shortfall of (1) lies in its unfavorable scaling with AB (i.e., the total number of joint actions), which is substantially larger than the total number of individual actions A + B. Whether it is possible to alleviate this curse of multiple agents in a two-player, zero-sum Markov game—and, if so, how to accomplish it—is the key question to be investigated in the current paper.

1.2. An Overview of Main Results

The objective of this paper is to design a sample-efficient, offline RL algorithm for learning Nash equilibria in two-player, zero-sum Markov games, ideally breaking the curse of multiple agents. Focusing on γ-discounted, infinite-horizon MGs, we propose a model-based paradigm—called the value iteration with lower confidence bounds for zero-sum Markov game (VI-LCB-Game)—that is capable of learning an ε-approximate Nash equilibrium with sample complexity

where $C_{\mathsf{\text{clipped}}}^{\star }$ is the so-called clipped unilateral concentrability coefficient (formalized in Assumption 2) and always satisfies $C_{\mathsf{\text{clipped}}}^{\star }\le C^{\star }$. Our result strengthens prior theory in Cui and Du (2022b) by a factor of $\min \{A,B\}$ (if we view the horizon length H in finite-horizon MGs and the effective horizon $\frac{1}{1-\gamma }$ in the infinite-horizon counterpart as equivalence). To demonstrate that this bound is essentially unimprovable, we develop a matching minimax lower bound (up to some logarithmic factor), thus settling this problem. Our algorithm is a pessimistic variant of value iteration with carefully designed Bernstein-style penalties, which requires neither sample splitting nor sophisticated schemes such as reference-advantage decomposition. The fact that our sample complexity result holds for the full ε-range (i.e., any $\epsilon \in \left(0,\frac{1}{1-\gamma }\right]$) unveils that sample efficiency is achieved without incurring any burn-in cost.

Finally, when finalizing the current paper, we became aware of an independent study (Cui and Du 2022a; posted to arXiv on June 1, 2022) that also manages to overcome the curse of multiple agents in a two-player, zero-sum Markov game, on which we elaborate toward the end of Section 3.

1.3. Notation

Before proceeding, let us introduce some notation that is used throughout. With slight abuse of notation, we use P to denote a probability transition kernel and the associated probability transition matrix exchangeably. We also use the notation μ exchangeably for a probability distribution and its associated probability vector (and we often do not specify whether μ is a row vector or column vector as long as it is clear from the context). For any two vectors $x={[x_i]}_{i=1}^n$ and $y={[y_i]}_{i=1}^n$, we use $x°y={[x_i{y}_i]}_{i=1}^n$ to denote their Hadamard product, and we also define $x^2={[x_i^2]}_{i=1}^n$ in an entry-wise fashion. For a finite set $\mathcal{S}=\{1,\dots ,S\}$, we let $\mathrm{\Delta }(\mathcal{S})≔\{x\in {\mathbb{R}}^S|1^{\top }x=1,x\ge 0\}$ represent the probability simplex over the set $\mathcal{S}$.

2. Problem Formulation

In this section, we introduce the background of zero-sum Markov games, followed by a description of the offline data set.

2.1. Preliminaries

2.1.1. Zero-Sum, Two-Player Markov Games.

Consider a discounted, infinite-horizon, zero-sum MG (Shapley 1953, Littman 1994) as represented by the tuple $\mathcal{M}\mathcal{G}=(\mathcal{S},\mathcal{A},\mathcal{B},P,r,\gamma )$. Here, $\mathcal{S}=\{1,\dots ,S\}$ is the shared state space; $\mathcal{A}=\{1,\dots ,A\}$ (respectively, $\mathcal{B}=\{1,\dots ,B\}$) is the action space of the max-player (min-player); $P:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to \mathrm{\Delta }(\mathcal{S})$ is the (a priori unknown) probability transition kernel, where $P(s\prime |s,a,b)$ denotes the probability of transitioning from state s to state $s\prime$ if the max-player executes action a and the min-player chooses action b; $r:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to [0,1]$ is the reward function such that $r(s,a,b)$ indicates the immediate reward observed by both players in state s when the max-player takes action a and the min-player takes action b; and $\gamma \in (0,1)$ is the discount factor with $\frac{1}{1-\gamma }$ commonly referred to as the effective horizon. Throughout this paper, we primarily focus on the scenario in which S, A, B, and $\frac{1}{1-\gamma }$ could all be large. Additionally, for notational simplicity, we define the vector $P_{s,a,b}\in {\mathbb{R}}^{1\times S}$ as $P_{s,a,b}≔P(·|s,a,b)$ for any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$.

2.1.2. Policy, Value Function, Q-Function, and Occupancy Distribution.

Let $\mu :\mathcal{S}\to \mathrm{\Delta }(\mathcal{A})$ and $\nu :\mathcal{S}\to \mathrm{\Delta }(\mathcal{B})$ be (possibly random) stationary policies of the max-player and the min-player, respectively. In particular, $\mu (·|s)\in \mathrm{\Delta }(\mathcal{A})$ ($\nu (·|s)\in \mathrm{\Delta }(\mathcal{B})$) specifies the action selection probability of the max-player (min-player) in state s. The value function $V^{\mu ,\nu }:\mathcal{S}\to \mathbb{R}$ for a given product policy $\mu \times \nu$ is defined as

where the expectation is taken with respect to the randomness of the trajectory ${\{(s_t,a_t,b_t)\}}_{t\ge 0}$ induced by the product policy $\mu \times \nu$ (i.e., for any $t\ge 0$, the players take $a_t\sim \mu (·|s_t)$ and $b_t\sim \nu (·|s_t)$ independently conditional on the past) and the probability transition kernel P (i.e., $s_{t+1}\sim P(·|s_t,a_t,b_t)$ for $t\ge 0$). Similarly, we can define the Q-function $Q^{\mu ,\nu }:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to \mathbb{R}$ for a given product policy $\mu \times \nu$ as follows:where the actions are drawn from $\mu \times \nu$ except for the initial time step (namely, for any $t\ge 1$, we execute $a_t\sim \mu (·|s_t)$ and $b_t\sim \nu (·|s_t)$ independently conditional on the past). Additionally, for any state distribution $\rho \in \mathrm{\Delta }(\mathcal{S})$, we introduce the following notation tailored to the weighted value function of policy pair $(\mu ,\nu )$:

Moreover, we define the discounted occupancy measures associated with an initial state distribution $\rho \in \mathrm{\Delta }(\mathcal{S})$ and the product policy $\mu \times \nu$ as follows:

where the sample trajectory ${\{(s_t,a_t,b_t)\}}_{t\ge 0}$ is initialized with $s_0\sim \rho$ and then induced by the product policy $\mu \times \nu$ and the transition kernel P as before. It is clearly seen from the preceding definition that

2.1.3. Nash Equilibrium.

In general, the two players have conflicting goals with the max-player aimed at maximizing the value function and the min-player minimizing the value function. As a result, a standard compromise in Markov games becomes finding an NE. To be precise, a policy pair $({\mu }^{\star },{\nu }^{\star })$ is said to be a Nash equilibrium if no player can benefit from unilaterally changing the player’s own policy given the opponent’s policy (Nash 1951); that is, for any policies $\mu :\mathcal{S}\to \mathrm{\Delta }(\mathcal{A})$ and $\nu :\mathcal{S}\to \mathrm{\Delta }(\mathcal{B})$, one has

As is well-known (Shapley 1953), there exists at least one Nash equilibrium $({\mu }^{\star },{\nu }^{\star })$ in the discounted, two-player, zero-sum Markov game, and every NE results in the same value function:

In addition, when the max-player’s policy μ is fixed, it is clearly seen that the MG reduces to a single-agent Markov decision process (MDP). In light of this, we define, for each $s\in \mathcal{S}$,

each of which corresponds to the optimal value function of one player with the opponent’s policy frozen. Moreover, for any policy pair $(\mu ,\nu )$, the following weak duality property always holds:

In this paper, our goal can be posed as calculating a policy pair $(\widehat{\mu },\widehat{\nu })$ such that

where $\rho \in \mathrm{\Delta }(\mathcal{S})$ is some prescribed initial state distribution and $\epsilon \in \left(0,\frac{1}{1-\gamma }\right]$ denotes the target accuracy level. The gap $V^{\star ,\widehat{\nu }}(\rho )-V^{\widehat{\mu },\star }(\rho )$ is often referred to as the duality gap of $(\widehat{\mu },\widehat{\nu })$ in the rest of the present paper.

2.2. Offline Data Set (Batch Data Set)

Suppose that we have access to a historical data set containing a batch of N sample transitions $\mathcal{D}={\{(s_i,a_i,b_i,s_i^{\prime })\}}_{1\le i\le N}$, which are generated independently from a distribution $d_{\mathsf{b}}\in \mathrm{\Delta }(\mathcal{S}\times \mathcal{A}\times \mathcal{B})$ and the probability transition kernel P, namely,

The goal is to learn an approximate Nash equilibrium on the basis of this historical data set.

In general, the data distribution $d_{\mathsf{b}}$ might deviate from the one generated by a Nash equilibrium $({\mu }^{\star },{\nu }^{\star })$. As a result, whether reliable learning is feasible depends heavily upon the quality of the historical data. To quantify the quality of the data distribution, Cui and Du (2022b) introduce the following unilateral concentrability condition.

Assumption 1

(Unilateral Concentrability). Suppose that the following quantity

$$C^{\star }≔\max \left\{\underset{\mu ,s,a,b}{\sup }\frac{d^{\mu ,{\nu }^{\star }}\left(s,a,b;\rho \right)}{d_{\mathsf{b}}\left(s,a,b\right)},\underset{\nu ,s,a,b}{\sup }\frac{d^{{\mu }^{\star },\nu }\left(s,a,b;\rho \right)}{d_{\mathsf{b}}\left(s,a,b\right)}\right\}$$(6)

is finite, where we define $0/0=0$ by convention. This quantity $C^{\star }$ is termed the unilateral concentrability coefficient.

In words, this quantity $C^{\star }$ employs certain density ratios to measure the distribution mismatch between the target distribution and the data distribution in hand. On the one hand, Assumption 1 is substantially weaker than the type of uniform coverage requirement (which imposes a uniform bound on the density ratio $\frac{d^{\mu ,\nu }\left(s,a,b;\rho \right)}{d_{\mathsf{b}}\left(s,a,b\right)}$ over all $(\mu ,\nu )$ simultaneously) as (6) freezes the policy of one side, exhausting over all policies of the other side. On the other hand, Assumption 1 remains more stringent than a single-policy coverage requirement (which only requires the data set to cover the part of the state–action space reachable by a given policy pair $({\mu }^{\star },{\nu }^{\star })$) as (6) requires the data to cover those state–action pairs reachable by any unilateral deviation from the target policy pair $({\mu }^{\star },{\nu }^{\star })$. As posited by Cui and Du (2022b), unilateral coverage (i.e., a finite $C^{\star }<\infty$) is necessary for learning Nash equilibria in Markov games, which stands in sharp contrast to the single-agent case in which single-policy concentrability suffices for finding the optimal policy (Rashidinejad et al. 2021, Xie et al. 2021, Li et al. 2024b).

In this paper, we introduce a modified assumption that might give rise to slightly improved sample complexity bounds.

Assumption 2

(Clipped Unilateral Concentrability). Suppose that the following quantity

$$C_{\mathsf{\text{clipped}}}^{\star }≔\max \left\{\underset{\mu ,s,a,b}{\sup }\frac{\min \left\{d^{\mu ,{\nu }^{\star }}\left(s,a,b;\rho \right),\frac{1}{S\left(A+B\right)}\right\}}{d_{\mathsf{b}}\left(s,a,b\right)},\underset{\nu ,s,a,b}{\sup }\frac{\min \left\{d^{{\mu }^{\star },\nu }\left(s,a,b;\rho \right),\frac{1}{S\left(A+B\right)}\right\}}{d_{\mathsf{b}}\left(s,a,b\right)}\right\}$$(7)

is finite, where we define $0/0=0$ by convention. This quantity $C_{\mathsf{\text{clipped}}}^{\star }$ is termed the clipped unilateral concentrability coefficient.

In a nutshell, when $d^{\mu ,{\nu }^{\star }}(s,a,b;\rho )$ or $d^{{\mu }^{\star },\nu }(s,a,b;\rho )$ is reasonably large (i.e., larger than $\frac{1}{S\left(A+B\right)}$), Assumption 2 no longer requires the data distribution $d_{\mathsf{b}}$ to scale proportionally with $d^{\mu ,{\nu }^{\star }}$ or $d^{{\mu }^{\star },\nu }$, thus resulting in a (slight) relaxation of Assumption 1. Comparing (7) with (6) immediately reveals that

holds all the time. Further, it is straightforward to verify that $C^{\star }\ge \max \{A,B\}$; in comparison, $C_{\mathsf{\text{clipped}}}^{\star }$ can be as small as $\frac{2AB}{S(A+B)}$ as shown in our lower bound construction in Online Appendix EC.2.

3. Algorithm and Main Theory

In this section, we propose a pessimistic, model-based, offline algorithm—called VI-LCB-Game—to solve the two-player, zero-sum Markov games. The proposed algorithm is then shown to achieve minimax-optimal sample complexity in finding an approximate Nash equilibrium of the Markov game given offline data.

3.1. Algorithm Design

3.1.1. The Empirical Markov Game.

With the offline data set ${\{(s_i,a_i,b_i,s_i^{\prime })\}}_{1\le i\le N}$ in hand, we can readily construct an empirical Markov game. To do so, we first compute the sample size

for each $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$. The empirical transition kernel $\widehat{P}:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to \mathrm{\Delta }(\mathcal{S})$ is then constructed as follows:for any $s\prime \in \mathcal{S}$ and any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$. Throughout this paper, we often let ${\widehat{P}}_{s,a,b}\in {\mathbb{R}}^{1\times S}$ abbreviate $\widehat{P}(·|s,a,b)$. In addition, the empirical reward function $\widehat{r}:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to \mathbb{R}$ is taken to befor any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$. Armed with these components, we arrive at an empirical, zero-sum Markov game, denoted by $\widehat{\mathcal{M}\mathcal{G}}=(\mathcal{S},\mathcal{A},\mathcal{B},\widehat{P},\widehat{r},\gamma )$.

3.1.2. Pessimistic Bellman Operators.

Recall that the classic Bellman operator $\mathcal{T}:{\mathbb{R}}^{\mathit{\text{SAB}}}\to {\mathbb{R}}^{\mathit{\text{SAB}}}$ is defined such that (Shapley 1953, Lagoudakis and Parr 2002), for any $Q:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to \mathbb{R}$,

where $V:\mathcal{S}\to \mathbb{R}$ is the value function associated with the input Q, that is,

Note, however, that we are in need of modified versions of the Bellman operator in order to accommodate the offline setting. In this paper, we introduce the pessimistic Bellman operator ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}$ (${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}$) for the max-player (min-player) as follows: for every $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$,

where V is again defined in (10). The additional term $\beta (s,a,b;V)$ is incorporated into the operators in order to implement pessimism; informally, we anticipate this penalty term to help ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}$ (${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}$) produce a conservative estimate of the Q-function from the max-player’s (min-player’s) viewpoint. Here and throughout, we choose this term based on Bernstein-style concentration bounds; specifically, we takefor some sufficiently large constant $C_{\mathsf{b}}>0$, where $1-\delta$ denotes the target success probability, and the empirical variance term is defined as

It is well-known that the classic Bellman operator $\mathcal{T}$ satisfies the γ-contraction property, which guarantees fast global convergence of classic value iteration. As it turns out, the pessimistic Bellman operators introduced also enjoy the γ-contraction property in the sense that

see Lemma 1 for precise statements.

3.1.3. Pessimistic Value Iteration with Bernstein-Style Penalty.

With the pessimistic Bellman operators in place, we are positioned to present the proposed paradigm. Our algorithm maintains the Q-function iterates $\{Q_{\mathsf{\text{pe}},t}^{-}\}$, the policy iterates $\{{\mu }_t^{-}\}$ and $\{{\nu }_t^{-}\}$, and the value function iterates $\{V_{\mathsf{\text{pe}},t}^{-}\}$ from the max-player’s perspective; at the same time, it also maintains an analogous group of iterates $\{Q_{\mathsf{\text{pe}},t}^{+}\}$, $\{{\mu }_t^{+}\}$ and $\{{\nu }_t^{+}\}$, and $\{V_{\mathsf{\text{pe}},t}^{+}\}$ from the min-player’s perspective. The updates of the two groups of iterates are carried out in a completely decoupled manner except when determining the final output.

In what follows, let us describe the update rules from the max-player’s perspective. For notational simplicity, we write $\mu (s)≔\mu (·|s)\in \mathrm{\Delta }(\mathcal{A})$ and $\nu (s)≔\nu (·|s)\in \mathrm{\Delta }(\mathcal{B})$ whenever it is clear from the context. In each round $t=1,2,\dots ,$ we carry out the following update rules:

1. Updating Q-function estimates: Run a pessimistic variant of value iteration to yield

   $$Q_{\mathsf{\text{pe}},t}^{-}={\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}\left(Q_{\mathsf{\text{pe}},t-1}^{-}\right).$$(15)

   The γ-contraction property (14) helps ensure sufficient progress made in each iteration of this update rule.

2. Updating policy estimates: We then adjust the policies based on the updated Q-function estimates (15). Specifically, for each $s\in \mathcal{S}$, we compute the Nash equilibrium $({\mu }_t^{-}(s),{\nu }_t^{-}(s))\in \mathrm{\Delta }(\mathcal{A})\times \mathrm{\Delta }(\mathcal{B})$ of the zero-sum matrix game with payoff matrix $Q_{\mathsf{\text{pe}},t}^{-}(s,·,·)$. It is worth noting that there is a host of methods for efficiently calculating the NE of a zero-sum matrix game; prominent examples include linear programming and no-regret learning (Raghavan 1994, Freund and Schapire 1999, Rakhlin and Sridharan 2013, Roughgarden 2016).

3. Policy evaluation: for each $s\in \mathcal{S}$, update the value function estimates based on the updated policies $({\mu }_t^{-}(s),{\nu }_t^{-}(s))$ as follows:

   $$V_{\mathsf{\text{pe}},t}^{-}(s)=\underset{a\sim {\mu }_t^{-}(s),b\sim {\nu }_t^{-}(s)}{\mathbb{E}}\left[Q_{\mathsf{\text{pe}},t}^{-}\left(s,a,b\right)\right].$$

The updates for $\{Q_{\mathsf{\text{pe}},t}^{+}\},\hspace{0.17em}\{{\mu }_t^{+}\}$, and $\{{\nu }_t^{+}\}$ from the min-player’s perspective are carried out in an analogous and completely independent manner; see Algorithm 1 for details.

3.1.4. Final Output.

By running these update rules for $T=\lceil \frac{\log (N/(1-\gamma ))}{\log (1/\gamma )}\rceil$ iterations, we arrive at the Q-function estimates

in addition to two sets of policy estimates

The final policy estimate of the algorithm is then chosen to be

The full algorithm is summarized in Algorithm 1.

Algorithm 1

(VI-LCB-Game)

Initialization: set $Q_{\mathsf{\text{pe}},0}^{-}(s,a,b)=0$ and $Q_{\mathsf{\text{pe}},0}^{+}(s,a,b)=\frac{1}{1-\gamma }$ for all $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$; set $T=\lceil \frac{\log (N/(1-\gamma ))}{\log (1/\gamma )}\rceil$.

Compute the empirical transition kernel $\widehat{P}$ as (8) and the empirical reward function $\widehat{r}$ as (9).

For: $t=1,\dots ,T$ do

• Update

  $$\begin{array}{ll}{Q}_{\mathsf{\text{pe}},t}^{-}\left(s,a,b\right)\hfill & ={\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}\left(Q_{\mathsf{\text{pe}},t-1}^{-}\right)=\max \left\{\widehat{r}\left(s,a,b\right)+\gamma {\widehat{P}}_{s,a,b}{V}_{\mathsf{\text{pe}},t-1}^{-}-\beta \left(s,a,b;V_{\mathsf{\text{pe}},t-1}^{-}\right),\hspace{0.17em}0\right\},\hfill \\ Q_{\mathsf{\text{pe}},t}^{+}\left(s,a,b\right)\hfill & ={\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}\left(Q_{\mathsf{\text{pe}},t-1}^{+}\right)=\min \left\{\widehat{r}\left(s,a,b\right)+\gamma {\widehat{P}}_{s,a,b}{V}_{\mathsf{\text{pe}},t-1}^{+}+\beta \left(s,a,b;V_{\mathsf{\text{pe}},t-1}^{+}\right),\hspace{0.17em}\frac{1}{1-\gamma }\right\},\hfill \end{array}$$
  where
  $$\beta \left(s,a,b;V\right)=\min \left\{\max \left\{\sqrt{\frac{C_{\mathsf{b}} \log \frac{N}{\delta }}{N\left(s,a,b\right)}{\mathsf{\text{Var}}}_{{\widehat{P}}_{s,a,b}}\left(V\right)},\frac{2C_{\mathsf{b}} \log \frac{N}{\delta }}{\left(1-\gamma \right)N\left(s,a,b\right)}\right\},\frac{1}{1-\gamma }\right\}+\frac{4}{N}$$
  for some sufficiently large constant $C_{\mathsf{b}}>0$ with ${\mathsf{\text{Var}}}_{{\widehat{P}}_{s,a,b}}(V)$ defined in (13).

• For each $s\in \mathcal{S}$, compute

  $$\begin{array}{ll}({\mu }_t^{-}(s),{\nu }_t^{-}(s))\hfill & =\mathsf{\text{MatrixNash}}(Q_{\mathsf{\text{pe}},t}^{-}\left(s,·,·\right)),\hfill \\ ({\mu }_t^{+}(s),{\nu }_t^{+}(s))\hfill & =\mathsf{\text{MatrixNash}}(Q_{\mathsf{\text{pe}},t}^{+}\left(s,·,·\right)),\hfill \end{array}$$
  where, for any matrix $M\in {\mathbb{R}}^{A\times B}$, the function $\mathsf{\text{MatrixNash}}(M)$ returns a solution $(\widehat{w},\widehat{z})$ to the minimax program ${\max }_{w\in \mathrm{\Delta }(\mathcal{A})} {\min }_{z\in \mathrm{\Delta }(\mathcal{B})}\hspace{0.17em}{w}^{\top }Mz$.

• For each $s\in \mathcal{S}$, update

  $$\begin{array}{ll}{V}_{\mathsf{\text{pe}},t}^{-}(s)\hfill & =\underset{a\sim {\mu }_t^{-}(s),b\sim {\nu }_t^{-}(s)}{\mathbb{E}}\left[Q_{\mathsf{\text{pe}},t}^{-}\left(s,a,b\right)\right],\hfill \\ V_{\mathsf{\text{pe}},t}^{+}(s)\hfill & =\underset{a\sim {\mu }_t^{+}(s),b\sim {\nu }_t^{+}(s)}{\mathbb{E}}\left[Q_{\mathsf{\text{pe}},t}^{+}\left(s,a,b\right)\right].\hfill \end{array}$$

 Output: the policy pair $(\widehat{\mu },\widehat{\nu })$, where $\widehat{\mu }={\{{\mu }_T^{-}(s)\}}_{s\in \mathcal{S}}$ and $\widehat{\nu }={\{{\nu }_T^{+}(s)\}}_{s\in \mathcal{S}}$.

3.2. Theoretical Guarantees

Our main result is to uncover the intriguing sample efficiency of the proposed model-based algorithm. This is formally stated as follows with the proof postponed to Section 6.

Theorem 1.

Consider any initial state distribution $\rho \in \mathrm{\Delta }(\mathcal{S})$ and suppose that Assumption 2 holds. Assume that $1/2\le \gamma <1$ and consider any $\delta \in (0,1)$ and $\epsilon \in \left(0,\frac{1}{1-\gamma }\right]$. Then, with probability exceeding $1-\delta$, the policy pair $(\widehat{\mu },\widehat{\nu })$ returned by Algorithm 1 satisfies

$$V^{\widehat{\mu },\star }\left(\rho \right)-\epsilon \le V^{\star }\left(\rho \right)\le V^{\star ,\widehat{\nu }}\left(\rho \right)+\epsilon ,$$

as long as the sample size exceeds

$$N\ge c_1\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{{\left(1-\gamma \right)}^3{\epsilon }^2}\log \frac{N}{\delta }$$

for some sufficiently large constant $c_1>0$.

Remark 1.

Our result and analysis are inspired by prior works that show that model-based RL achieves, in multiple settings, sample efficiency without the need of variance reduction (Agarwal et al. 2020; Li et al. 2024a, b). The proof of this sample complexity bound entails several key analysis ingredients: (i) a leave-one-out analysis argument that proves effective in decoupling complicated statistical dependency and (ii) a careful self-bounding trick (i.e., upper bounding a certain quantity by a contraction of itself in addition to some other error terms) to derive a sharp control of the target duality gap. See Section 6 for details. Although techniques such as leave-one-out analysis are used in some prior RL literature (Agarwal et al. 2020; Li et al. 2024a, b), as far as we know, our work applies this technique for the first time to multiagent reinforcement learning. It has been observed that extending the algorithmic or analysis ideas in single-agent RL to the multiagent counterpart often leads to suboptimal sample complexity bounds that scale linearly in the total number of joint actions AB (Zhang et al. 2020a, Cui and Du 2022b). In contrast, our analysis framework leads to an optimal sample complexity bound that scales linearly in the total number of individual actions A + B.

Remark 2.

A line of recent works focuses on instance-optimality of RL algorithms (Khamaru et al. 2021a, b; Mou et al. 2022). However, it remains challenging to establish instance-dependent bounds for multiagent RL even in two-player, zero-sum Markov games because of the difficulties arising from offline data and multiagent settings. Unlike RL with a generative model (simulator) that can generate independent samples for all state–action pairs, offline RL suffers from substantially more challenges, such as distribution shift and limited data coverage, making it more difficult to derive instance-dependent error bounds. In addition, the prior literature Khamaru et al. (2021a) that establishes instance optimality of variance-reduced Q-learning algorithms for the optimal value estimation problem requires one of the following two conditions: the optimal policy is unique or a meaningful sample complexity bound that depends on an optimality gap can be obtained. However, neither condition has a direct analog in zero-sum Markov games; this is because the Nash equilibrium in a zero-sum Markov game is not unique in general, and there is no well-defined analog of optimality gap for zero-sum Markov games. Detailed discussion on the challenges and difficulty of extending our analysis to develop instance-dependent error bounds can be found in Section 6.5.

The sample complexity needed for Algorithm 1 to compute a policy pair with ε-duality gap is at most

$$\tilde{O}\left(\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{{\left(1-\gamma \right)}^3{\epsilon }^2}\right),$$(18)

which accommodates any target accuracy within the range $\left(0,\frac{1}{1-\gamma }\right]$. In addition to linear dependency on $C_{\mathsf{\text{clipped}}}^{\star }$, the sample complexity bound (18) scales linearly (as opposed to quadratically) with the aggregate size A + B of the individual action spaces. It is noteworthy that our algorithm is a fairly straightforward implementation of the model-based approach (except that the pessimism principle is incorporated) and does not require either sample splitting or sophisticated schemes such as variance reduction (Zhang et al. 2020a, b; Li et al. 2021; Xie et al. 2021; Yan et al. 2023).

As it turns out, the preceding sample complexity theory for Algorithm 1 matches the minimax lower limit modulo some logarithmic term as asserted by the following theorem. This minimax lower bound—whose proof is postponed to Online Appendix EC.2—is inspired by prior lower bound theory for single-agent MDPs (e.g., Azar et al. 2013, Li et al. 2024b) and might shed light on how to establish lower bounds for other game-theoretic settings.

Theorem 2.

Consider any $S\ge 2,\hspace{0.17em}A\ge 2$, $B\ge 2,\hspace{0.17em}\gamma \in \left[\frac{2}{3},1\right)$, and $C_{\mathsf{\text{clipped}}}^{\star }\ge \frac{2AB}{S(A+B)}$, and define the set

$$\begin{array}{ll}\mathsf{\text{MG}}\left(C_{\mathsf{\text{clipped}}}^{\star }\right)\hfill & ≔\left\{\right\{\mathcal{M}\mathcal{G},\rho ,d_{\mathsf{b}}\}\hspace{0.17em}|\hspace{0.17em}|\mathcal{S}|=S,\hspace{0.17em}|\mathcal{A}|=A,\hspace{0.17em}|\mathcal{B}|=B,\hfill \\ \hfill & \hspace{1em}\hspace{1em}\rho \in \mathrm{\Delta }\left(\mathcal{S}\right),\hspace{0.17em}{d}_{\mathsf{b}}\in \mathrm{\Delta }\left(\mathcal{S}\times \mathcal{A}\times \mathcal{B}\right),\hspace{0.17em}\exists an NE \left({\mu }^{\star },{\nu }^{\star }\right) of \mathcal{M}\mathcal{G} \mathit{\text{such}} \mathit{\text{that}} \hfill \\ \hfill & \max \left\{\underset{\mu ,s,a,b}{\sup }\frac{\min \left\{d^{\mu ,{\nu }^{\star }}\left(s,a,b;\rho \right),\frac{1}{S\left(A+B\right)}\right\}}{d_{\mathsf{b}}\left(s,a,b\right)},\underset{\nu ,s,a,b}{\sup }\frac{\min \left\{d^{{\mu }^{\star },\nu }\left(s,a,b;\rho \right),\frac{1}{S\left(A+B\right)}\right\}}{d_{\mathsf{b}}\left(s,a,b\right)}\}=C_{\mathsf{\text{clipped}}}^{\star }\right\}.\hfill \end{array}$$

Then, there exist some universal constants $c_2,c_{\epsilon }>0$ such that, for any $\epsilon \in \left(0,\frac{1}{c_{\epsilon }(1-\gamma )\log (A+B)}\right]$, if the sample size obeys

$$N<\frac{c_{2}S\left(A+B\right)C_{\mathsf{\text{clipped}}}^{\star }}{{\left(1-\gamma \right)}^3{\epsilon }^2 \log \left(A+B\right)},$$

then one necessarily has

$$\underset{(\widehat{\mu },\widehat{\nu })}{\inf } \underset{\{\mathcal{M}\mathcal{G},\rho ,d_{\mathsf{b}}\}\in \mathsf{\text{MG}}(C_{\mathsf{\text{clipped}}}^{\star })}{\sup }\mathbb{E}\left[V^{\star ,\widehat{\nu }}\left(\rho \right)-V^{\widehat{\mu },\star }\left(\rho \right)\right]\ge \epsilon .$$

Here, the infimum is taken over all estimators $(\widehat{\mu },\widehat{\nu })$ for the Nash equilibrium based on the batch data set $\mathcal{D}={\{(s_i,a_i,b_i,s_i^{\prime })\}}_{i=1}^n$ generated according to (5).

Remark 3.

The target we are estimating is the NE of a zero-sum MG, which is more challenging than standard statistical estimation problems in the sense that (i) NE is not unique in general and (ii) the error metric is a duality gap. It is challenging to use standard proof frameworks such as Fano’s and Le Cam’s methods to derive a meaningful lower bound for this problem. To overcome this challenge, we construct a family of hard Markov game instances indexed by a binary parameter $\theta \in {\{0,1\}}^{\max \{A,B\}}$ and then put a prior distribution over this set and compute the posterior probability of failure to differentiate each entry of θ. These steps taken together carefully allow us to compute the desired minimax risk.

As a direct implication of Theorem 2, if the total number of samples in the offline data set obeys

$$N<\frac{c_{2}S\left(A+B\right)C_{\mathsf{\text{clipped}}}^{\star }}{{\left(1-\gamma \right)}^3{\epsilon }^2 \log \left(A+B\right)},$$

then one can construct a hard Markov game instance such that no algorithm whatsoever can reach a duality gap below ε. This, taken collectively (18), unveils, up to some logarithmic factor, the minimax statistical limit for finding NEs based on offline data.

Our theory makes remarkable improvement upon prior art, which can be seen through comparisons with the most relevant prior work (Cui and Du 2022b) (even though the focus therein is finite-horizon, zero-sum MGs). On a high level, Cui and Du (2022b) propose an algorithm that combines pessimistic value iteration with variance reduction (also called reference-advantage decomposition; Zhang et al. 2020b), which provably finds an ε-Nash policy pair using

$$\tilde{O}\left(\frac{C^{\star }SAB{H}^3}{{\epsilon }^2}\right)$$(19)

sample trajectories provided that $\epsilon \le 1/H$. Here, H stands for the horizon length of the finite-horizon Markov game, and $C^{\star }$ is the unilateral concentrability coefficient tailored to the finite-horizon setting. Despite the difference between discounted infinite- and finite-horizon settings, our algorithm design and theory achieve several improvements upon Cui and Du (2022b):

• Perhaps most importantly, our result scales linearly in the total number of individual actions A + B (as opposed to the number of joint actions AB as in Cui and Du 2022b), which manages to alleviate the curse of multiple agents in two-player, zero-sum Markov games.

• Our theory accommodates the full ε-range $\left(0,\frac{1}{1-\gamma }\right]$, which is much wider than the range $(0,1/H]$ covered by Cui and Du (2022b) (if we view the effective horizon $\frac{1}{1-\gamma }$ in the infinite-horizon case and the horizon length H in the finite-horizon counterpart as equivalence).

• The algorithm design herein is substantially simpler than Cui and Du (2022b): it neither requires sample splitting to decouple statistical dependency, nor relies on reference-advantage decomposition techniques to sharpen the horizon dependency.

When we were finalizing the present manuscript, we became aware of the independent work Cui and Du (2022a) proposing a different offline algorithm—based on incorporation of strategy-wise lower confidence bounds—that improved the prior art as well. When it comes to two-player, zero-sum Markov games with finite horizon and nonstationary transition kernels, Cui and Du (2022a, algorithm 1) provably yields an ε-Nash policy pair using

$$\tilde{O}\left(\frac{C^{\star }S\left(A+B\right)H^4}{{\epsilon }^2}\right)$$(20)

sample trajectories, each containing H samples. This bound (20) is at least a factor of H above the minimax limit. It is worth noting that Cui and Du (2022a) is able to accommodate offline, multiagent, general-sum MGs although the algorithm proposed therein becomes computationally intractable when going beyond two-player, zero-sum MGs.

4. Related Works

4.1. Offline RL and Pessimism Principle

The principle of pessimism in the face of uncertainty, namely, being conservative in value estimation of those state–action pairs that have been under-covered, has been adopted extensively in recent development of offline RL. A highly incomplete list includes Kumar et al. (2020), Kidambi et al. (2020), Yu et al. (2020; 2021a, b), Yin et al. (2021a, c), Rashidinejad et al. (2021), Jin et al. (2021b), Xie et al. (2021), Liu et al. (2020), Zhang et al. (2021c), Chang et al. (2021), Yin and Wang (2021), Uehara and Sun (2021), Munos (2003, 2007), Zanette et al. (2021), Yan et al. (2023), Li et al. (2022a, 2024b), Shi et al. (2022), Cui and Du (2022b), Zhong et al. (2022), Lu et al. (2022), Wang et al. (2022), and Xu and Liang (2022), which unveils the efficacy of the pessimism principle in both model-based and model-free approaches. Among this body of prior works, the ones that are most related to the current paper are Cui and Du (2022a, b), and Zhong et al. (2022), both of which focus on episodic, finite-horizon, zero-sum Markov games with two players. More specifically, Cui and Du (2022b) demonstrate that a unilateral concentrability condition is necessary for learning NEs in offline settings and propose a pessimistic value iteration with reference-advantage decomposition to enable sample efficiency. Zhong et al. (2022) propose a pessimistic minimax value iteration algorithm, which achieves appealing sample complexity in the presence of linear function representation and was recently improved by Xiong et al. (2022). In the concurrent work, Cui and Du (2022a) propose a different pessimistic algorithm that designs lower confidence bounds for policy pairs instead of state–action pairs; for two-player, zero-sum MGs, their algorithm is capable of achieving a sample complexity proportional to A + B. In the single-agent, offline RL setting, Rashidinejad et al. (2021), Yan et al. (2023), and Li et al. (2024b) study offline RL for infinite-horizon MDPs, and Jin et al. (2021b), Xie et al. (2021), Shi et al. (2022), and Li et al. (2024b) look at the finite-horizon episodic counterpart, all of which operate upon some single-policy concentrability assumptions. Among these works, Li et al. (2024b) and Yan et al. (2023) achieved minimax-optimal sample complexity $\tilde{O}\left(\frac{S{C}^{\star }}{{(1-\gamma )}^3{\epsilon }^2}\right)$ for discounted, infinite-horizon MDPs by means of model-based and model-free algorithms, respectively; similar results have been established for finite-horizon MDPs as well (Xie et al. 2021; Yin et al. 2021b, c; Shi et al. 2022; Li et al. 2024b).

4.2. Multiagent RL and Markov Games

The concept of Markov games—also under the name of stochastic games—dates back to Shapley (1953), which has become a central framework to model competitive multiagent decision making. A large strand of prior works studies how to efficiently solve Markov games when perfect model description is available (Littman 1994, 2001; Hu and Wellman 2003; Hansen et al. 2013; Perolat et al. 2015; Daskalakis et al. 2020, 2023; Cen et al. 2021; Wei et al. 2021; Zhao et al. 2021; Chen et al. 2022; Mao and Başar 2022). Recent years have witnessed much activity in studying the sample efficiency of learning Nash equilibria in zero-sum Markov games, covering multiple different types of sampling schemes; for instance, Wei et al. (2017), Xie et al. (2020), Bai et al. (2020), Bai and Jin (2020), Liu et al. (2021), Jin et al. (2021a), Song et al. (2021), Mao and Başar (2022), Daskalakis et al. (2023), Tian et al. (2021), and Chen et al. (2022) focus on the online explorative environments, whereas Zhang et al. (2020a) pays attention to the scenario that assumes sampling access to a generative model. Whereas the majority of these works exhibits a sample complexity that scales at least as $\tilde{O}(\mathit{\text{SAB}})$ in order to learn an approximate NE, the recent work Bai et al. (2020) proposes a V-learning algorithm attaining a sample complexity that scales linearly with $S(A+B)$, thus matching the minimax-optimal lower bound up to a factor of H². When a generative model is available, Li et al. (2022b) further develops an algorithm that learns ε-Nash using $\tilde{O}\left(\frac{H^{4}S(A+B)}{{\epsilon }^2}\right)$ samples, which attains the minimax lower bound for nonstationary, finite-horizon MGs. The setting of general-sum, multiplayer Markov games is much more challenging given that learning Nash equilibria is known to be PPAD-complete (Daskalakis et al. 2009, Daskalakis 2013). Shifting attention to more tractable solution concepts, Jin et al. (2021a), Daskalakis et al. (2023), Mao and Başar (2022), and Song et al. (2021) propose algorithms that provably learn (coarse) correlated equilibria with sample complexities that scale linearly with ${\max }_i A_i$ (where A_i is the number of actions of the ith player), thereby breaking the curse of multiagents. Additionally, there are also several works investigating the turn-based setting in which the two players take actions in turn; see Sidford et al. (2020), Cui and Yang (2021), Jia et al. (2019), and Jin et al. (2022). Moreover, another two works Zhang et al. (2021b) and Abe and Kaneko (2020) study offline sampling oracles under uniform coverage requirements (which are clearly more stringent than the unilateral concentrability assumption). The interested readers are also referred to Zhang et al. (2021a) and Yang and Wang (2020) for an overview of recent development.

4.3. Model-Based RL

The method proposed in the current paper falls under the category of model-based algorithms, which decouple model estimation and policy learning (planning). The model-based approach is extensively studied in the single-agent setting, including the online exploration setting (Azar et al. 2017, Zhang et al. 2023), the case with a generative model (Azar et al. 2013, Agarwal et al. 2020, Jin and Sidford 2021, Wang et al. 2021, Li et al. 2024a), the offline RL setting (Xie et al. 2021, Li et al. 2024b), and turn-based Markov games (Cui and Yang 2021). Encouragingly, the model-based approach is capable of attaining minimax-optimal sample complexities in a variety of settings (e.g., Azar et al. 2017, Agarwal et al. 2020, Zhang et al. 2023, Li et al. 2024b), sometimes even without incurring any burn-in cost (Cui and Yang 2021; Zhang et al. 2023; Li et al. 2024a, b). The method proposed in Cui and Du (2022b) also exhibits the flavor of a model-based algorithm although an additional variance reduction scheme is incorporated in order to optimize the horizon dependency.

5. Additional Notation

Let us collect a set of additional notations that are used in the analysis. First of all, for any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$, any vector $V\in {\mathbb{R}}^S$, and any probability transition kernel $P:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to \mathrm{\Delta }(\mathcal{S})$, we define

where $P_{s,a,b}$ abbreviates $P(·|s,a,b)$ as usual. When the max-player’s policy μ is fixed, the Markov game reduces to a (single-agent) MDP for the min-player. For any MDP, it is known that there exists at least one policy that simultaneously maximizes the value function (Q-function) for all states (state–action pairs) (Bertsekas 2017). In light of this, when the policy μ of the max-player is frozen, we denote by ${\nu }_{\mathsf{\text{br}}}(\mu )$ the optimal policy of the min-player, which is often referred to as the best response of the min-player when the max-player adopts policy μ. Similarly, we can define the best response of the max-player when the min-player adopts policy ν, which we denote by ${\mu }_{\mathsf{\text{br}}}(\nu )$. These allow one to definefor all $s\in \mathcal{S}$, andfor all $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$. Note that the definitions of $V^{\mu ,\star }$ and $V^{\star ,\nu }$ here are consistent with the ones in Section 2.

6. Proof of Theorem 1

Toward proving Theorem 1, we first state a slightly stronger result as follows.

Theorem 3.

Consider any initial state distribution $\rho \in \mathrm{\Delta }(\mathcal{S})$ and suppose that Assumption 2 holds. Assume that $1/2\le \gamma <1$. Then, with probability exceeding $1-\delta$, the policy pair $(\widehat{\mu },\widehat{\nu })$ returned by Algorithm 1 satisfies

$$V^{\star }\left(\rho \right)-V^{\widehat{\mu },\star }\left(\rho \right)\le c_0\sqrt{\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{{\left(1-\gamma \right)}^{3}N}\log \frac{N}{\delta }}+c_0\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{{\left(1-\gamma \right)}^{2}N}\log \frac{N}{\delta },$$(22a)

$$V^{\star ,\widehat{\nu }}\left(\rho \right)-V^{\star }\left(\rho \right)\le c_0\sqrt{\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{{\left(1-\gamma \right)}^{3}N}\log \frac{N}{\delta }}+c_0\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{{\left(1-\gamma \right)}^{2}N}\log \frac{N}{\delta }$$(22b)

for some sufficiently large constant $c_0>0$. As an immediate consequence, the duality gap of $(\widehat{\mu },\widehat{\nu })$ obeys, with probability at least $1-\delta$, that

$$V^{\star ,\widehat{\nu }}\left(\rho \right)-V^{\widehat{\mu },\star }\left(\rho \right)\le 2c_0\sqrt{\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{{\left(1-\gamma \right)}^{3}N}\log \frac{N}{\delta }}+2c_0\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{{\left(1-\gamma \right)}^{2}N}\log \frac{N}{\delta }.$$(23)

As can be straightforwardly verified, Theorem 1 is a direct consequence of Theorem 3 (by taking the right-hand side of (23) to be no larger than ε).

The remainder of this section is, thus, dedicated to establishing Theorem 3. Before proceeding, let us now take a moment to provide a brief road map of the proof.

1. We first show in Section 6.2 that the pessimistic Bellman operators ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}$ and ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}$ introduced in (11) are both monotone and γ-contractive and admit unique fixed points $Q_{\mathsf{\text{pe}},t}^{-\star }$ and $Q_{\mathsf{\text{pe}},t}^{-\star }$, respectively. These properties reveal that the pessimistic value iterations ${\{Q_{\mathsf{\text{pe}},t}^{-}\}}_{1\le t\le T}$ (${\{Q_{\mathsf{\text{pe}},t}^{+}\}}_{1\le t\le T}$) in Algorithm 1 converge to $Q_{\mathsf{\text{pe}},t}^{-\star }$ ($Q_{\mathsf{\text{pe}},t}^{-\star }$) at a geometric rate, and therefore, it suffices to analyze the fixed points $Q_{\mathsf{\text{pe}},t}^{-\star }$ and $Q_{\mathsf{\text{pe}},t}^{+\star }$.

2. Next, we show Bernstein-style concentration bounds for random quantities such as $({\widehat{P}}_{s,a,b}-P_{s,a,b})V_{\mathsf{\text{pe}},t}^{-\star }$ and $({\widehat{P}}_{s,a,b}-P_{s,a,b})V_{\mathsf{\text{pe}},t}^{+\star }$ in Section 6.3, in which $V_{\mathsf{\text{pe}},t}^{-\star }$ and $V_{\mathsf{\text{pe}},t}^{+\star }$ are the value functions associated with $Q_{\mathsf{\text{pe}},t}^{-\star }$ and $Q_{\mathsf{\text{pe}},t}^{+\star }$. Because of the complicated statistical dependency between ${\widehat{P}}_{s,a,b}$ and $V_{\mathsf{\text{pe}},t}^{-\star }$, we use a leave-one-out argument to establish this concentration result in Lemma 2.

3. Finally, based on the aforementioned results, we derive error bounds for $V^{\star }(\rho )-V^{\widehat{\mu },\star }(\rho )$ and $V^{\star ,\widehat{\nu }}(\rho )-V^{\star }(\rho )$ in Section 6.4. Our analysis makes use of a self-bounding trick, which allows one to derive sharp estimation error bounds that turn out to be minimax-optimal.

6.1. Preliminary Facts

Before continuing, we collect several preliminary facts that are useful throughout.

1. For any $Q_1,Q_2:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to \mathbb{R}$, we have

   $${\Vert V_1-V_2\Vert }_{\infty }\le {\Vert Q_1-Q_2\Vert }_{\infty },$$(24)
   where V₁ (V₂) denotes the value function associated with Q₁ (Q₂); see (10) for the precise definition.

2. For any $V_1,V_2:\mathcal{S}\to \left[0,\frac{1}{1-\gamma }\right]$, any probability transition kernel $P:\mathcal{S}\times \mathcal{A}\times \mathcal{B}\to \mathrm{\Delta }(\mathcal{S})$, and any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$, we have

   $$|{\mathsf{\text{Var}}}_{P_{s,a,b}}\left(V_1\right)-{\mathsf{\text{Var}}}_{P_{s,a,b}}\left(V_2\right)|\le \frac{4}{1-\gamma }{\Vert V_1-V_2\Vert }_{\infty },$$(25)
   where ${\mathsf{\text{Var}}}_{P,s,a,b}(V)$ is defined in (21).

3. As a consequence, we also know that, for any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$ and any $V_1,V_2:\mathcal{S}\to \left[0,\frac{1}{1-\gamma }\right]$, the corresponding penalty terms (cf. (12)) obey

   $$|\beta \left(s,a,b;V_1\right)-\beta \left(s,a,b;V_2\right)|\le 2{\Vert V_1-V_2\Vert }_{\infty }.$$(26)

The proof of the preceding results can be found in Online Appendix EC.1.1.

6.2. Step 1: Key Properties of Pessimistic Bellman Operators

Recall the definition of the pessimistic Bellman operators ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}$ and ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}$ introduced in (11). The following lemma gathers a couple of key properties of these two operators.

Lemma 1.

The following properties hold true:

• (Monotonicity) For any $Q_1\ge Q_2$, we have ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}(Q_1)\ge {\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}(Q_2)$ and ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}(Q_1)\ge {\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}(Q_2)$.

• (Contraction) Both operators are γ-contractive in the ${\mathcal{\ell }}_{\infty }$ sense, that is,

  $$\Vert {\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}\left(Q_1\right)-{\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}\left(Q_2\right){\Vert }_{\infty }\le \gamma \Vert Q_1-Q_2{\Vert }_{\infty },\hspace{1em}\hspace{1em}\Vert {\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}\left(Q_1\right)-{\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}\left(Q_2\right){\Vert }_{\infty }\le \gamma {\Vert Q_1-Q_2\Vert }_{\infty }$$
  for any Q₁ and Q₂.

• (Uniqueness of fixed points) ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}$ (${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{+}$) has a unique fixed point $Q_{\mathsf{\text{pe}}}^{-\star }$ ($Q_{\mathsf{\text{pe}}}^{+\star }$), which also satisfies $0\le Q_{\mathsf{\text{pe}}}^{\star -}(s,a,b)\le \frac{1}{1-\gamma }$ $(0\le Q_{\mathsf{\text{pe}}}^{+\star }(s,a,b)\le \frac{1}{1-\gamma })$ for any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$.

Proof.

See Online Appendix EC.1.2. □

Next, we make note of several immediate consequences of Lemma 1. Here and throughout, $V_{\mathsf{\text{pe}}}^{-\star }$ and $V_{\mathsf{\text{pe}}}^{+\star }$ are defined to be the value functions (see (10)) associated with $Q_{\mathsf{\text{pe}}}^{-\star }$ and $Q_{\mathsf{\text{pe}}}^{+\star }$, respectively.

• First, the preceding lemma implies that

  $$Q_{\mathsf{\text{pe}},t}^{-}\le Q_{\mathsf{\text{pe}}}^{-\star }\hspace{1em}(\forall t\ge 0)\hspace{1em} \text{hence}\hspace{1em}\hspace{1em}{Q}_{\mathsf{\text{pe}}}^{-}\le Q_{\mathsf{\text{pe}}}^{-\star }.$$(27)

  To see this, we first note that $Q_{\mathsf{\text{pe}},0}^{-}=0\le Q_{\mathsf{\text{pe}}}^{-\star }$. Next, suppose that $Q_{\mathsf{\text{pe}},t}^{-}\le Q_{\mathsf{\text{pe}}}^{-\star }$ for some iteration $t\ge 0$; then, the monotonicity of ${\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}$ (cf. Lemma 1) tells us that

  $$Q_{\mathsf{\text{pe}},t+1}^{-}={\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}\left(Q_{\mathsf{\text{pe}},t}^{-}\right)\le {\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}(Q_{\mathsf{\text{pe}}}^{-\star })=Q_{\mathsf{\text{pe}}}^{-\star },$$
  from which (27) follows.

• In addition, the γ-contraction property in Lemma 1 leads to

  $${\Vert V_{\mathsf{\text{pe}}}^{-}-V_{\mathsf{\text{pe}}}^{-\star }\Vert }_{\infty }\le {\Vert Q_{\mathsf{\text{pe}}}^{-}-Q_{\mathsf{\text{pe}}}^{-\star }\Vert }_{\infty }\le \frac{1}{N},$$(28)
  and, to justify this, observe that
  $$\begin{array}{ll}{\Vert Q_{\mathsf{\text{pe}},t}^{-}-Q_{\mathsf{\text{pe}}}^{-\star }\Vert }_{\infty }\hfill & ={\Vert {\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}(Q_{\mathsf{\text{pe}},t-1}^{-})-{\widehat{\mathcal{T}}}_{\mathsf{\text{pe}}}^{-}(Q_{\mathsf{\text{pe}}}^{-\star })\Vert }_{\infty }\le \gamma {\Vert Q_{\mathsf{\text{pe}},t-1}^{-}-Q_{\mathsf{\text{pe}}}^{-\star }\Vert }_{\infty }\hfill \\ \hfill & \le \cdots \le {\gamma }^t{\Vert Q_{\mathsf{\text{pe}},0}^{-}-Q_{\mathsf{\text{pe}}}^{-\star }\Vert }_{\infty }\le \frac{{\gamma }^t}{1-\gamma },\hfill \end{array}$$
  which, together with $T=\lceil \frac{\log (N/(1-\gamma ))}{\log (1/\gamma )}\rceil$ and (24), gives
  $${\Vert V_{\mathsf{\text{pe}}}^{-}-V_{\mathsf{\text{pe}}}^{-\star }\Vert }_{\infty }\le \Vert Q_{\mathsf{\text{pe}}}^{-}-Q_{\mathsf{\text{pe}}}^{-\star }{\Vert }_{\infty }=\Vert Q_{\mathsf{\text{pe}},T}^{-}-Q_{\mathsf{\text{pe}}}^{-\star }{\Vert }_{\infty }\le \frac{{\gamma }^T}{1-\gamma }\le \frac{1}{N}.$$

• A similar argument also yields

  $$Q_{\mathsf{\text{pe}}}^{+}\ge Q_{\mathsf{\text{pe}}}^{+\star },\hspace{1em}\Vert Q_{\mathsf{\text{pe}}}^{+}-Q_{\mathsf{\text{pe}}}^{+\star }{\Vert }_{\infty }\le 1/N,\hspace{1em}\Vert V_{\mathsf{\text{pe}}}^{+}-V_{\mathsf{\text{pe}}}^{+\star }{\Vert }_{\infty }\le 1/N.$$(29)

6.3. Step 2: Decoupling Statistical Dependency and Establishing Pessimism

To proceed, we rely on the following theorem to quantify the difference between $\widehat{P}$ and P when projected onto a value function direction.

Lemma 2.

For any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$ satisfying $N(s,a,b)\ge 1$ with probability exceeding $1-\delta$,

$$|({\widehat{P}}_{s,a,b}-P_{s,a,b})\tilde{V}|\le \tilde{c}\sqrt{\frac{1}{N\left(s,a,b\right)}{\mathsf{\text{Var}}}_{{\widehat{P}}_{s,a,b}}(\tilde{V})\log \frac{N}{\delta }}+\tilde{c}\frac{\log \frac{N}{\delta }}{\left(1-\gamma \right)N\left(s,a,b\right)}$$(30)

for some sufficiently large constant $\tilde{c}>0$, and

$${\mathsf{\text{Var}}}_{{\widehat{P}}_{s,a,b}}(\tilde{V})\le 2{\mathsf{\text{Var}}}_{P_{s,a,b}}(\tilde{V})+O\left(\frac{1}{{\left(1-\gamma \right)}^{2}N\left(s,a,b\right)}\log \frac{N}{\delta }\right)$$(31)

hold simultaneously for all $\tilde{V}\in {\mathbb{R}}^S$ satisfying $0\le \tilde{V}\le \frac{1}{1-\gamma }1$ and $\min \{{\Vert \tilde{V}-V_{\mathsf{\text{pe}}}^{-\star }\Vert }_{\infty },{\Vert \tilde{V}-V_{\mathsf{\text{pe}}}^{+\star }\Vert }_{\infty }\}\le 1/N$.

The proof is deferred to Online Appendix EC.1.3.

In words, the first result (30) delivers some Bernstein-type concentration bound, whereas the second result (31) guarantees that the empirical variance estimate (i.e., the plug-in estimate) is close to the true variance. It is worth noting that Lemma 2 does not require $\tilde{V}$ to be statistically independent from ${\widehat{P}}_{s,a,b}$, which is particularly crucial when coping with the complicated statistical dependency of our iterative algorithm. The proof of Lemma 2 is established upon a leave-one-out analysis argument (see, e.g., Chen et al. 2019a, b; Agarwal et al. 2020; Ma et al. 2020; Chen et al. 2021; Li et al. 2024a, b) that helps decouple statistical dependency; see details in Online Appendix EC.1.3. Armed with Lemma 2, we can readily see that

holds for any $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$ satisfying $N(s,a,b)\ge 1$ and any $\tilde{V}$ satisfying the conditions in Lemma 2. In turn, this important fact allows one to justify that $Q_{\mathsf{\text{pe}}}^{-}$ ($Q_{\mathsf{\text{pe}}}^{+}$) is indeed an upper (lower) bound on $Q^{\widehat{\mu },\star }$ ($Q^{\star ,\widehat{\nu }}$) as formalized subsequently.

Lemma 3.

With probability exceeding $1-\delta$, it holds that

$$Q_{\mathsf{\text{pe}}}^{-}\le Q^{\widehat{\mu },\star },\hspace{1em}\hspace{1em}{Q}_{\mathsf{\text{pe}}}^{+}\ge Q^{\star ,\widehat{\nu }},\hspace{1em}{V}_{\mathsf{\text{pe}}}^{-}\le V^{\widehat{\mu },\star }\hspace{1em}\mathit{\text{and}}\hspace{1em}\hspace{1em}{V}_{\mathsf{\text{pe}}}^{+}\ge V^{\star ,\widehat{\nu }}.$$

The proof is provided in Online Appendix EC.1.4.

This lemma makes clear a key rationale for the principle of pessimism: we want the Q-function estimates to be always conservative uniformly over all entries.

6.4. Step 3: Bounding $V^{\star }(\rho )-V^{\widehat{\mu },\star }(\rho )$ and $V^{\star ,\widehat{\nu }}(\rho )-V^{\star }(\rho )$

Before proceeding to bound $V^{\star }-V^{\widehat{\mu },\star }$, we first develop a lower bound on $V_{\mathsf{\text{pe}}}^{-}$ given that $V^{\widehat{\mu },\star }$ is lower bounded by $V_{\mathsf{\text{pe}}}^{-}$ (according to Lemma 3). Toward this end, we invoke the definition of $V_{\mathsf{\text{pe}}}^{-}$ to reach

where we set the policy of the max-player to be ${\mu }^{\star }$ on the right-hand side of the preceding equation. Clearly, there exists a deterministic policy ${\nu }_0:\mathcal{S}\to \mathrm{\Delta }(\mathcal{B})$ such thatfor any $s\in \mathcal{S}$; for instance, one can simply set, for any $s\in \mathcal{S}$,with ${\mathrm{𝟙}}_{b_s}$ denoting a probability vector that is nonzero only in b_s. This deterministic policy ν₀ helps us lower bound $V_{\mathsf{\text{pe}}}^{-}$ as accomplished in the following lemma. Here and as follows, we define two vectors $r^{{\mu }^{\star },{\nu }_0},{\beta }^{{\mu }^{\star },{\nu }_0}\in {\mathbb{R}}^S$ and a probability transition kernel $P^{{\mu }^{\star },{\nu }_0}:\mathcal{S}\to \mathrm{\Delta }(\mathcal{S})$ restricted to ${\mu }^{\star }$ and ν₀ such that, for any $s,s\prime \in \mathcal{S}$,

Lemma 4.

With probability exceeding $1-\delta$, we have

$$V_{\mathsf{\text{pe}}}^{-}\ge r^{{\mu }^{\star },{\nu }_0}+\gamma P^{{\mu }^{\star },{\nu }_0}{V}_{\mathsf{\text{pe}}}^{-}-2{\beta }^{{\mu }^{\star },{\nu }_0}.$$(37)

The proof is deferred to Online Appendix EC.1.5.

In addition, we can invoke Lemma 3 and the fact that $V^{\star }=V^{{\mu }^{\star },{\nu }^{\star }}=V^{{\mu }^{\star },\star }$ to reach

which motivates us to look at $V^{{\mu }^{\star },{\nu }_0}-V_{\mathsf{\text{pe}}}^{-}$. Toward this, we note that the Bellman equation tells us that

Taking (37) and (39) collectively yields

thus resulting in a self-bounding type of relation. Applying (40) recursively, we arrive at thatholds for all positive integers n. Letting $n\to \infty$ and recalling that the vector $d^{{\mu }^{\star },{\nu }_0}≔{[d^{{\mu }^{\star },{\nu }_0}(s;\rho )]}_{s\in \mathcal{S}}$ obeys (see (2))we arrive atwhere the last line makes use of the fact that ${\Vert {\rho }^{\top }{(P^{{\mu }^{\star },{\nu }_0})}^n\Vert }_1=1$ for any $n\ge 1$, and hence, ${\gamma }^n{\rho }^{\top }{(P^{{\mu }^{\star },{\nu }_0})}^n\to 0$ as $n\to \infty$ when $\gamma <1$.

In order to further control (42), we resort to the following lemma for bounding ${(d^{{\mu }^{\star },{\nu }_0})}^{\top }{\beta }^{{\mu }^{\star },{\nu }_0}$, whose proof can be found in Online Appendix EC.1.6.

Lemma 5.

There exists some large enough universal constant $c_6>0$ such that

$${(d^{{\mu }^{\star },{\nu }_0})}^{\top }{\beta }^{{\mu }^{\star },{\nu }_0}\le c_6\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{\left(1-\gamma \right)N}\log \frac{N}{\delta }+c_6\sqrt{\frac{C_{\mathsf{\text{clipped}}}^{\star }S\left(A+B\right)}{N\left(1-\gamma \right)}\log \frac{N}{\delta }}.$$

This is with probability exceeding $1-\delta$.

To finish, taking (38), (42) and Lemma 5 together gives

This completes the proof for Claim (22a). The proof for the other claim (22b) follows from an almost identical argument and is, hence, omitted.

6.5. Discussion: Instance-Dependent Statistical Bounds?

Thus far, we have presented the proof of Theorem 1 that concerns the minimax optimality of the model-based algorithm. Note that a recent line of work has attempted to move beyond minimax-optimal statistical guarantees and pursue more refined instance-optimal (or locally minimax) performance guarantees (Khamaru et al. 2021a, b, Mou et al. 2022). Here, we take a moment to discuss the challenges that need to be overcome in order to extend our analysis in an instance-optimal fashion.

• A crucial step that allows us to obtain error bounds that scale linearly with A + B instead of the ones that scale linearly with AB in Cui and Du (2022b) is the introduction of an auxiliary policy ν₀ in (34). This allows us to upper bound the error with

  $${\rho }^{\top }(V^{{\mu }^{\star },{\nu }_0}-V_{\mathsf{\text{pe}}}^{-})\le \frac{2}{1-\gamma }{(d^{{\mu }^{\star },{\nu }_0})}^{\top }{\beta }^{{\mu }^{\star },{\nu }_0};$$
  see (42). Whereas this facilitates our analysis, the terms ${\beta }^{{\mu }^{\star },{\nu }_0}$ (defined in (36)) and $d^{{\mu }^{\star },{\nu }_0}$ (defined in (41)) both depend on the auxiliary policy ν₀. Because of the complicated dependency between ν₀ (which is determined by a random function $V_{\mathsf{\text{pe}}}^{-}$) and the model parameters, it remains quite challenging to connect these error terms with instance-dependent quantities (i.e., model parameters) without losing optimality.

• Note that we might be able to resolve this issue in a coarse way, for example, by taking the supremum over all possible policy ν:

  $${\rho }^{\top }(V^{{\mu }^{\star },{\nu }_0}-V_{\mathsf{\text{pe}}}^{-})\le \frac{2}{1-\gamma }\underset{\nu \in \mathrm{\Delta }(\mathcal{B})}{\sup }{(d^{{\mu }^{\star },\nu })}^{\top }{\beta }^{{\mu }^{\star },\nu }.$$(43)

  Let us assume for the moment that this could work (despite the potential suboptimality of this error bound) and see what this lead to. By checking the proof of Lemma 5, we can see that, in order to upper bound (43) in an instance-optimal manner, it is important to relate ${\mathsf{\text{Var}}}_{P_{s,a,b}}(V_{\mathsf{\text{pe}}}^{-})$ to model parameters for all $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$. Ideally, we can replace $V_{\mathsf{\text{pe}}}^{-}$ with the value function associated with the Nash equilibrium $V^{\star }$. However, based on our current analysis framework, we can only show that $V^{\star }(\rho )$ and $V_{\mathsf{\text{pe}}}^{-}(\rho )$ are close, which is insufficient to guarantee the closeness of ${\mathsf{\text{Var}}}_{P_{s,a,b}}(V_{\mathsf{\text{pe}}}^{-})$ and ${\mathsf{\text{Var}}}_{P_{s,a,b}}(V^{\star })$ for all $(s,a,b)\in \mathcal{S}\times \mathcal{A}\times \mathcal{B}$.

• As we briefly mention in Remark 2, prior literature Khamaru et al. (2021a), which establishes instance optimality for the optimal value estimation problem in single-agent RL under a generative model, requires one of the following two conditions: the optimal policy is unique or a sample complexity bound that depends on an optimality gap

  $$\mathrm{\Delta }≔\underset{\pi \in \mathbf{\Pi }\backslash {\mathbf{\Pi }}^{\star }}{\text{min}}{\Vert Q^{\star }-r-\gamma P^{\pi }{Q}^{\star }\Vert }_{\infty },$$(44)
  where $Q^{\star }$ is the optimal Q-function, $\mathbf{\Pi }$ is the set of deterministic policies, and ${\mathbf{\Pi }}^{\star }$ is the set of optimal (deterministic) policies. However, neither condition has a direct analog in two-player, zero-sum Markov games: for the first one, this is because the Nash equilibrium of a zero-sum Markov game is not unique in general; for the second one, this is because the Nash equilibrium policy pair is usually random, and it is not clear how to define a nonzero optimality gap such as (44).

In view of these challenges, our current analysis framework remains incapable of deriving instance-optimal performance guarantees. Accomplishing instance-optimal results for zero-sum Markov games might require substantially more refined analysis techniques, and we leave this important direction to future investigation.

7. Discussion

In the present paper, we propose a model-based offline algorithm, which leverages the principle of pessimism in solving two-player, zero-sum Markov games on the basis of past data. In order to find an ε-approximate Nash equilibrium of the Markov game, our algorithm requires no more than $\tilde{O}\left(\frac{S(A+B)C^{\star }}{{(1-\gamma )}^3{\epsilon }^2}\right)$ samples, and this sample complexity bound is provably minimax optimal for the entire range of target accuracy level $\epsilon \in \left(0,\frac{1}{1-\gamma }\right]$. Our theory improves upon prior sample complexity bounds in Cui and Yang (2021) in terms of the dependency on the size of the action space. Another appealing feature is the simplicity of our algorithm, which does not require complicated variance reduction schemes and is, hence, easier to implement and interpret. Moving forward, there are a couple of interesting directions that are worthy of future investigation. For instance, one natural extension is to explore whether the current algorithmic idea and analysis extend to multiagent, general-sum Markov games with the goal of learning other solution concepts of equilibria such as coarse correlated equilibria (given that finding Nash equilibria in general-sum games is PPAD-complete). Another topic of interest is to design model-free algorithms for offline NE learning in zero-sum or general-sum Markov games. Furthermore, the current paper focuses attention on tabular Markov games, and it would be of great interest to design sample-efficient, offline, multiagent algorithms in the presence of function approximation.