# Two-Timescale Q-Learning with Function Approximation
in Zero-Sum Stochastic Games

Zaiwei Chen$1,*$, Kaiqing Zhang$2$, Eric Mazumdar$1,\dagger$, Asuman
Ozdaglar$3$, Adam Wierman$1,\ddagger$
$1$CMS, Caltech<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
$2$ECE, University of Maryland, College Park<EMAIL_ADDRESS>
$3$EECS, MIT<EMAIL_ADDRESS>
###### Abstract

We consider two-player zero-sum stochastic games and propose a two-timescale
$Q$-learning algorithm with function approximation that is payoff-based,
convergent, rational, and symmetric between the two players. In two-timescale
$Q$-learning, the fast-timescale iterates are updated in spirit to the
stochastic gradient descent and the slow-timescale iterates (which we use to
compute the policies) are updated by taking a convex combination between its
previous iterate and the latest fast-timescale iterate. Introducing the slow
timescale as well as its update equation marks as our main algorithmic
novelty. In the special case of linear function approximation, we establish,
to the best of our knowledge, the first last-iterate finite-sample bound for
payoff-based independent learning dynamics of these types. The result implies
a polynomial sample complexity to find a Nash equilibrium in such stochastic
games.

To establish the results, we model our proposed algorithm as a two-timescale
stochastic approximation and derive the finite-sample bound through a
Lyapunov-based approach. The key novelty lies in constructing a valid Lyapunov
function to capture the evolution of the slow-timescale iterates.
Specifically, through a change of variable, we show that the update equation
of the slow-timescale iterates resembles the classical smoothed best-response
dynamics, where the regularized Nash gap serves as a valid Lyapunov function.
This insight enables us to construct a valid Lyapunov function via a
generalized variant of the Moreau envelope of the regularized Nash gap. The
construction of our Lyapunov function might be of broad independent interest
in studying the behavior of stochastic approximation algorithms.

## 1 Introduction

Recent years have seen remarkable successes of reinforcement learning (RL)
across many application domains, such as the game of Go [98], autonomous
driving [96], and robotics [61]. Among existing learning dynamics,
$Q$-learning [109] and its variants, such as the celebrated deep $Q$-network
(DQN) [81], are among the most popular. Although $Q$-learning was originally
proposed for single-agent RL, since many real-world applications by nature
involve the decision-making and learning of multiple agents in a shared
environment, it is natural to explore the potential of $Q$-learning in multi-
agent systems [21, 121]. However, in contrast to the single-agent setting
(where $Q$-learning has a solid theoretical foundation [106, 42, 31, 70]), due
to the interactions among agents in multi-agent systems, it is still an active
research topic to theoretically investigate whether the agents will jointly
achieve any performance guarantees if each of them independently implements
$Q$-learning (or its variants) [68, 93].

One of the key challenges in the deployment of $Q$-learning to real-world
applications is the curse of dimensionality [14, 87]. Specifically, the state
space in practical applications can be massively or even infinitely large,
making it computationally inefficient, if not intractable, to estimate and
store the $Q$-functions in learning. A common way to deal with this curse has
been to employ _function approximation_ for the $Q$-function. The high-level
idea is to approximate the solution of the RL problem from a chosen function
class, in which each function can be specified by a parameter (or weight) with
a much smaller dimension than the size of the state-action space of the
underlying RL problem. A representative example of $Q$-learning with function
approximation is the DQN, in which neural networks are used as a means for
approximating the $Q$-functions. Unfortunately, when function approximation is
employed, the “deadly triad”111The deadly triad refers to the combination of
off-policy sampling, function approximation, and bootstrapping [103]. [103]
would occur in $Q$-learning, resulting in the algorithm to diverge even under
linear function approximation, as demonstrated by several counterexamples in
the literature [3, 28]. To overcome the deadly triad, variants of $Q$-learning
have been proposed, such as Greedy-GQ [74], $Q$-learning with a target network
[28, 123], fitted $Q$-iteration [83], and Zap $Q$-learning [27], etc. However,
the results are mostly limited to the single-agent RL setting. Therefore, it
is natural to ask:

Is it possible to make $Q$-learning with function approximation provably work
for multi-agent systems?

In this work, we provide an answer to the above question in the benchmark
setting of two-player zero-sum stochastic games [97], where both players aim
to optimize their strategies in a competitive environment with evolving
states. Specifically, we propose a two-timescale $Q$-learning dynamics that
can incorporate function approximation. In particular, our learning dynamics
enjoy several desired properties in multi-agent reinforcement learning (MARL),
i.e., the learning dynamics are _independent_ (require no coordination between
the players [86]), _rational_ (each player will converge to the best-response
to the opponent if the opponent plays an asymptotically stationary policy
[19]), and _payoff-based_ (each player can only observe their realized payoff
as well as the system state during learning [44, 76]). On the theoretical
side, we establish the last-iterate finite-sample convergence guarantees in
the case when the function approximator has a linear structure. We detail our
contributions in the next subsection.

### 1.1 Main Contributions

##### The Learning Dynamics.

From the perspective of player $i\in\\{1,2\\}$, our two-timescale $Q$-learning
has two loops and maintains $4$ sets of parameters (or weights) $w^{i}_{t,k}$,
$\theta^{i}_{t,k}$, $\bar{w}_{t}^{i}$, and $\bar{\theta}_{t}^{i}$ that we use
to approximate the local $q$-functions, where $k$ is the inner-loop iteration
index and $t$ is the outer-loop iteration index. Similarly to the DQN in the
single-agent setting, the two sets of parameters $\bar{w}_{t}^{i}$ and
$\bar{\theta}_{t}^{i}$ are called the _target networks_ , which are fixed in
the inner loop and are synchronized to the last iterates $w_{t,K}^{i}$ and
$\theta_{t,K}^{i}$ from the inner loop, respectively, where $K$ is the number
of iterations of the inner loop. As will be illustrated in Section 3.2, such
an update rule for $\bar{w}_{t}^{i}$ and $\bar{\theta}_{t}^{i}$ can be
interpreted as an approximate version of the _minimax value iteration_ , a
well-known equilibrium-computation algorithm for zero-sum stochastic games
[97].

The parameters $w^{i}_{t,k}$ and $\theta^{i}_{t,k}$ are updated in the inner
loop in a two-timescale fashion, and is designed to find the values of a
sequence of induced zero-sum matrix games, which are needed for implementing
the minimax value iteration in the outer loop. Specifically, $w_{t,k}^{i}$ is
the fast-timescale parameter and is updated according to the projected
stochastic gradient descent (SGD) for minimizing a variant of the Bellman
error, which reduces to the temporal-difference (TD)-learning when using
linear function approximation [102, 107]. The slow-timescale iterate
$\theta_{t,k}^{i}$ is used to compute the policies for playing the game, which
we propose to update by performing the convex combination
$\theta_{t,k+1}^{i}=(1-\beta_{k})\theta_{t,k}^{i}+\beta_{k}w_{t,k}^{i}$ with
stepsize $\beta_{k}>0$. Importantly, we show in Section 5 that the seemingly
simple update equation for $\theta_{t,k}^{i}$ can be interpreted as a
generalization of the smoothed best-response dynamics (classical learning
dynamics for solving zero-sum matrix games) implemented in the _parameter
space_ of the function approximation class.

Notably, introducing the two-timescale structure and the update equation of
the slow-timescale iterate $\theta_{t,k}^{i}$ is essential for connecting our
algorithm to the classical smoothed best-response dynamics, which also
motivates our Lyapunov-based analysis framework, and thus presents our main
algorithmic novelty.

##### The Last-Iterate Finite-Sample Analysis.

When using linear function approximation, under certain assumptions, we
establish the last-iterate finite-sample guarantee of our proposed learning
dynamics measured by the Nash gap, which implies a polynomial sample
complexity to find a Nash equilibrium of the underlying stochastic game. To
the best of our knowledge, this is the first result that shows last-iterate
convergence bounds for payoff-based independent learning dynamics under
function approximation for stochastic games. See Section 1.2 for a detailed
comparison with the existing literature.

##### A New Lyapunov Function for Analysis.

To establish the results, we model our proposed algorithm as a two-timescale
stochastic approximation, and the overall framework for our analysis is a
Lyapunov-based approach, which has been widely used to study stochastic
approximation algorithms in optimization and machine learning [5, 18, 65, 101,
31, 32]. The main challenge, as in all of the Lyapunov-based approaches, is
the construction of a valid Lyapunov function to capture the stability of the
iterates. As illustrated at the beginning of this subsection, one of our main
algorithmic novelties is to propose to update the slow-timescale iterate
$\theta_{t,k}^{i}$ (which we use to compute the policies) according to the
convex combination of the previous $\theta_{t,k}^{i}$ and the updated
$w_{t,k}^{i}$, as
$\theta_{t,k+1}^{i}=(1-\beta_{k})\theta_{t,k}^{i}+\beta_{k}w_{t,k}^{i}$.
Although such an update equation is simple to implement, it cannot be
interpreted as any standard RL algorithms such as policy iteration or policy
gradient [104], which makes studying $\\{\theta_{t,k}^{i}\\}$ particularly
challenging as there is no off-the-shelf Lyapunov function.

To overcome this challenge, firstly, we show that the update equation for the
local $q$-function associated with the parameter $\theta_{t,k}^{i}$ (after a
change of variable) mimics the update equation for the policies in the
smoothed best-response dynamics. This observation motivates us to use the
entropy-regularized Nash gap (ERNG) (which was known to be a valid Lyapunov
function for studying the smoothed best-response dynamics for matrix games
[52]) as a starting point for our Lyapunov function construction. However,
several issues prevent us from directly using the ERNG to study
$\theta_{t,k}^{i}$. For example, since the parameter of our function
approximator in general lives in a real-vector space (beyond the simplexes),
the entropy function (which is the regularizer in ERNG) is not even directly
well-defined in the parameter space. In addition, the change of variable that
enabled us to connect the update equation for $\theta_{t,k}^{i}$ to the
smoothed best-response dynamics requires the invertibility of the reward
matrix at each state, which may not hold in general. To resolve these issues,
we construct our Lyapunov function as the infimal convolution between the ERNG
and the square of a properly defined seminorm, which can also be interpreted
as a generalization of the Moreau envelope. The construction process is
detailed in Section 5, and the roadmap behind it might be of independent
interest.

### 1.2 Related Work

##### Single-Agent $Q$-Learning (with Function Approximation).

$Q$-learning, first proposed in [109], has become one of the most classical
and empirically popular algorithms in RL. As a result, the theoretical
understanding of $Q$-learning has attracted a great deal of attention. While
early literature has focused on the _asymptotic_ convergence guarantees [106,
54, 17], the more recent ones have shared the common interest in establishing
_finite-sample_ guarantees [42, 10, 11, 69, 71, 31, 55]. To overcome the curse
of dimensionality, $Q$-learning is usually incorporated with function
approximation in practical applications, a typical example of which is the DQN
[80]. However, since $Q$-learning is inherently a bootstrapped off-policy
learning algorithm (as our goal is to estimate the state-action value function
of an optimal policy while using samples collected from a potentially sub-
optimal one), under function approximation, the so-called phenomenon of
_deadly triad_ [103] will occur, and the algorithm can be unstable and easily
diverge [3, 28]. Overcoming the deadly triad in $Q$-learning (empirically and
also theoretically) has been an active research topic in the RL community, see
e.g., [77, 33, 49, 66, 126, 114, 22, 24, 123, 1, 119, 74, 105, 117, 58, 111,
125, 83, 40, 28, 78]. Among existing approaches, the introduction of a target
network in $Q$-learning [80, 123, 28] has proven to be an effective idea to
improve the stability of $Q$-learning, which is also used in the present work.

##### Finite-Sample Analysis.

We here focus on the most related works on finite-sample analysis through the
lens of stochastic approximation [90, 12, 16, 13, 63], which is an iterative
framework for solving root-finding problems, and has been widely applied in
continuous optimization (to solve for the zeros of the gradient operator) [65]
and RL (to solve the Bellman equation) [14]. Specifically, for linear
stochastic approximation, finite-sample analysis has been conducted in [15,
101, 64, 82, 62], and for nonlinear stochastic approximation, finite-sample
analysis has been conducted in [42, 71, 30, 89, 31, 116, 126], motivated by
various RL algorithms, such as $Q$-learning [109] and SARSA [92]. In addition,
stochastic gradient descent/ascent, as a special case of stochastic
approximation, has been extensively studied in the literature; see [18, 8, 65]
and the references therein for the finite-sample analysis of stochastic
gradient algorithms.

##### Best-Response-Based Dynamics in Stochastic Games.

Learning in games has been one of the core topics in game theory, and the
_best-response-based_ dynamics has served as one of the most classical ones
[48]. A prominent example of such dynamics is the (smoothed) fictitious play
[20, 91, 47], wherein each player responds with the (smoothed) best action it
can choose given its belief of the opponents’ strategies. Such best-response-
type dynamics provide a way to justify _equilibria_ as the long-run outcome of
bounded rational players optimizing myopically in games [48]. Despite their
long history, the literature on understanding the behavior of best-response-
type dynamics in _stochastic games_ [97, 43] – a fundamental framework of MARL
[72] and games in dynamic environments – has not become rich until recently
[93, 94, 95, 6, 7, 86, 29, 23]. The results, however, either concern the
tabular setting without function approximation, or have established only
_asymptotic_ convergence guarantees. Note that, in [68], the authors proposed
an individual $Q$-learning algorithm for zero-sum matrix games and established
the asymptotic convergence, where they connected the ordinary differential
equation (ODE) associated with their $Q$-function to the smoothed best-
response dynamics via the chain rule of calculus. This is similar in spirit to
our change of variable in the analysis.

The existing work most related to ours is [32], where the authors proposed a
learning algorithm that combines value iteration with smoothed best-response
dynamics, and conducted last-iterate finite-sample analysis. We now highlight
the main differences between this work and [32]. Algorithmically, since [32]
considers the tabular setting while we consider the function approximation
setting, the learning algorithms are largely different as we cannot directly
operate in the policy space, but have to operate in the parameter space of the
function approximation class. Technically, while both our work and [32] have
used Lyapunov-based arguments to establish finite-sample guarantees, the
Lyapunov functions, especially the one we use to study the slow-timescale
iterates $\\{\theta_{t,k}^{i}\\}$ for computing the policies, are largely
different from those used in [32]. In particular, since [32] directly uses the
smoothed best-response dynamics for their policy update equation, their
Lyapunov function for studying the policies is the ERNG, which is not the case
in our setting. The construction of a valid Lyapunov function for
$\theta_{t,k}^{i}$ presents one of our main technical novelties, as
illustrated in Section 1.1.

##### Independent Learning in Games.

Independent learning is a desired property in learning in games, where each
player makes decisions using local information, without communication or
coordination with other players. It has been extensively studied in normal-
form/matrix game settings in the forms of (smoothed) fictitious play [20, 91,
47], and in general, no-regret learning [25]. In stochastic games, provable
(independent) learning has received increasing attention recently [2, 37, 120,
93, 86, 110, 94, 95, 6, 7, 56, 100, 75, 122, 39, 38, 41, 99, 23, 32], as well
as the hardness results of independent learning [73, 46]. However, these works
all focused on the _tabular_ setting without resorting to function
approximation. Several recent works that have incorporated function
approximation in learning in stochastic games [113, 53, 124, 57, 26, 85, 45,
35] are mostly focused on the online exploration setting with regret
guarantees, a more challenging setting, without last-iterate convergence
guarantees. More importantly, most of these algorithms are either asymmetric
or require coordination among the agents, making them not independent learning
algorithms. Very recently, [108, 36] have both studied independent function
approximation for finite-horizon episodic stochastic games, in the online
exploration setting. However, they used the regret bound as the performance
metric, which does not imply last-iterate finite-sample guarantees and is thus
not comparable to our results.

##### Outline.

The rest of this paper is organized as follows. In Section 2, we present the
problem formulation, which includes the model of a zero-sum stochastic game
and the notation we will use throughout the paper. In Section 3, we introduce
our learning dynamics, i.e., the two-timescale $Q$-learning with function
approximation, its connection to the celebrated DQN algorithm, and how it can
be viewed as an approximate version of the minimax value iteration for zero-
sum stochastic games. In Section 4, we present our main theoretical result –
the finite-sample guarantee of our proposed learning dynamics, followed by a
detailed roadmap of the construction of a valid Lyapunov function to study the
slow-timescale iterates in Section 5. Finally, we conclude the paper in
Section 6. The complete proof of all theoretical results can be found in the
appendix.

## 2 Problem Formulation

An infinite-horizon two-player zero-sum stochastic game is defined by a tuple
$(\mathcal{S},\mathcal{A}^{1},\mathcal{A}^{2},p,R_{1},R_{2},\gamma)$, where
$\mathcal{S}$ is a finite state space, $\mathcal{A}^{i}$, $i\in\\{1,2\\}$, is
the finite action space of player $i$,
$p:\mathcal{S}\times\mathcal{A}^{1}\times\mathcal{A}^{2}\mapsto\Delta(\mathcal{S})$
is the state transition kernel (where $\Delta(\mathcal{S})$ stands for the
probability simplex supported on $\mathcal{S}$),
$R_{i}:\mathcal{S}\times\mathcal{A}^{i}\times\mathcal{A}^{-i}\mapsto\mathbb{R}$,
$i\in\\{1,2\\}$, is the reward function of player $i$, and $\gamma\in[0,1)$ is
the discount factor. Here, we use $-i$ to denote the index of the opponent of
player $i$. In zero-sum stochastic games, the reward functions satisfy
$\sum_{i=1,2}R_{i}(s,a^{i},a^{-i})=0$ for all $(s,a^{1},a^{2})$. Throughout,
we assume that $\max_{s,a^{1},a^{2}}|R_{1}(s,a^{1},a^{2})|\leq 1$, which is
without loss of generality because we are working with a finite stochastic
game222Our results can be generalized to a stochastic game with an infinite
state space but a finite action space. We focus on the finite setting for ease
of presentation..

At time step $k\geq 0$, denote the current state of the environment by
$\mathcal{S}_{k}$. For $i\in\\{1,2\\}$, player $i$ plays the game with an
action $A_{k}^{i}\sim\pi^{i}(\cdot\mid S_{k})$, where $\pi^{i}$ is called a
policy and can be viewed as a mapping from the state space $\mathcal{S}$ to
the set of probability distributions supported on the action space
$\mathcal{A}^{i}$. The environment then transitions to a new state $S_{k+1}$
based on the transition kernel $p(\cdot\mid S_{k},A_{k}^{1},A_{k}^{2})$, and
returns player $i$ a one-stage reward $R_{i}(S_{k},A_{k}^{i},A_{k}^{-i})$. The
process is then repeated.

Given a joint policy $\pi=(\pi^{1},\pi^{2})$ and $i\in\\{1,2\\}$, the local
state-action function
$q_{\pi}^{i}:\mathcal{S}\times\mathcal{A}^{i}\mapsto\mathbb{R}$, the global
value function $v_{\pi}^{i}:\mathcal{S}\mapsto\mathbb{R}$, and the expected
value function $U^{i}(\pi^{i},\pi^{-i})\in\mathbb{R}$ are defined as

$\displaystyle
q_{\pi}^{i}(s,a^{i})=\mathbb{E}_{\pi}\left[\sum_{k=0}^{\infty}\gamma^{k}R_{i}(S_{k},A_{k}^{i},A_{k}^{-i})\;\middle|\;S_{0}=s,A_{0}^{i}=a^{i}\right],\quad\forall\;(s,a^{i})\in\mathcal{S}\times\mathcal{A}^{i},$

$v_{\pi}^{i}(s)=\mathbb{E}_{A^{i}\sim\pi^{i}(\cdot|s)}[q_{\pi}^{i}(s,A^{i})]$
for all $s\in\mathcal{S}$, and
$U^{i}(\pi^{i},\pi^{-i})=\mathbb{E}_{S\sim\rho_{o}(\cdot)}[v_{\pi}^{i}(S)]$,
respectively, where we use $\mathbb{E}_{\pi}[\,\cdot\,]$ to indicate that the
actions are chosen according to the joint policy $\pi=(\pi^{1},\pi^{2})$, and
$\rho_{o}(\cdot)$ is an arbitrary initial distribution over the states.

Unlike in the single-agent setting, there may not exist a universal optimal
policy for each player. Therefore, we use the Nash gap
$\text{NG}(\cdot,\cdot)$ to measure the performance of a joint policy
$\pi=(\pi^{1},\pi^{2})$, which is defined as

$\displaystyle\text{NG}\left(\pi^{1},\pi^{2}\right)=\sum_{i=1,2}\big{\\{}\max_{\hat{\pi}^{i}}U^{i}(\hat{\pi}^{i},\pi^{-i})-U^{i}(\pi^{i},\pi^{-i})\big{\\}}.$

Note that $\text{NG}\left(\pi^{1},\pi^{2}\right)=0$ if and only if
$\max_{\hat{\pi}^{i}}U^{i}(\hat{\pi}^{i},\pi^{-i})=U^{i}(\pi^{i},\pi^{-i})$
for $i\in\\{1,2\\}$, which implies that the joint policy
$\pi=(\pi^{1},\pi^{2})$ is a Nash equilibrium [84].

##### Notation.

Given $i\in\\{1,2\\}$, for any
$x\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}||\mathcal{A}^{-i}|}$, we use
$x(s)$ to denote an $|\mathcal{A}^{i}|\times|\mathcal{A}^{-i}|$ matrix with
its $(a^{i},a^{-i})$’s entry being $x(s,a^{i},a^{-i})$. Similarly, for any
$y\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}|}$, we use $y(s)$ to denote an
$|\mathcal{A}^{i}|$ – dimensional vector with its $a^{i}$-th entry being
$y(s,a^{i})$. For example, for any $s\in\mathcal{S}$, $\pi^{i}(s)$ is a vector
living in the probability simplex $\Delta(\mathcal{A}^{i})$, and
$q_{\pi}^{i}(s)$ is a real-valued vector in $\mathbb{R}^{|\mathcal{A}^{i}|}$.

## 3 The Learning Dynamics

In this section, we present our two-timescale $Q$-learning with function
approximation. Before that, we formally introduce the function approximation
architecture.

##### Function Approximation.

Given $i\in\\{1,2\\}$, let
$q_{w^{i}}^{i}\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}|}$ be the
approximation of the local $q$-function of player $i$, where
$w^{i}\in\mathbb{R}^{d_{i}}$ is the parameter (or weight) and $d_{i}$ is a
positive integer that is usually smaller than
$|\mathcal{S}||\mathcal{A}^{i}|$. We assume that $q_{w^{i}}^{i}$ is
differentiable with respect to $w^{i}$. As an example of function
approximation, consider a neural network with the dimension of the input
vector being $d_{i}$ and the dimension of the output vector being
$|\mathcal{S}||\mathcal{A}^{i}|$. With the parametrization of the local
$q$-function introduced above, we consider using the following softmax policy
class:

$\displaystyle\Pi^{i}=\left\\{\pi_{\theta^{i}}^{i}(a^{i}\mid
s)=\frac{\exp(q_{\theta^{i}}^{i}(s,a^{i})/\tau)}{\sum_{\tilde{a}^{i}\in\mathcal{A}^{i}}\exp(q_{\theta^{i}}^{i}(s,\tilde{a}^{i})/\tau)},\,\forall\,(s,a^{i})\,\middle|\,\theta^{i}\in\mathbb{R}^{d_{i}}\right\\},$

where $\tau>0$ is called the temperature. For simplicity of notation, given
$\tau>0$, let
$\sigma_{\tau}:\mathbb{R}^{|\mathcal{A}^{i}|}\mapsto\mathbb{R}^{|\mathcal{A}^{i}|}$
be the softmax operator with temperature $\tau$ defined as
$[\sigma_{\tau}(x^{i})](a^{i})=\exp(x^{i}(a^{i})/\tau)/\sum_{\tilde{a}^{i}\in\mathcal{A}^{i}}\exp(x^{i}(\tilde{a}^{i})/\tau)$
for any $a^{i}\in\mathcal{A}^{i}$ and
$x^{i}\in\mathbb{R}^{|\mathcal{A}^{i}|}$.

##### Linear Function Approximation.

A special case of function approximation is when $q_{w^{i}}^{i}$ is linear in
$w^{i}$, which is called linear function approximation [104, 107].
Specifically, for $i\in\\{1,2\\}$, player $i$ is associated with a feature
matrix
$\Phi^{i}\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}|\times\mathbb{R}^{d_{i}}}$,
and the local $q$-function is parametrized as $q_{w^{i}}^{i}=\Phi^{i}w^{i}$.
Note that the tabular setting corresponds to $\Phi^{i}=I^{i}$, where $I^{i}$
is the identity matrix with dimension $|\mathcal{S}||\mathcal{A}^{i}|$. We
assume without loss of generality that $\Phi^{i}$ has linearly independent
columns and is normalized so that $\max_{s,a^{i}}\|\phi^{i}(s,a^{i})\|_{2}\leq
1$ [14], where $\phi^{i}(s,a^{i})\in\mathbb{R}^{d_{i}}$ denotes the
$(s,a^{i})$-th row of $\Phi^{i}$. Moreover, for any $s\in\mathcal{S}$, let
$\Phi^{i}_{s}\in\mathbb{R}^{|\mathcal{A}^{i}|\times d_{i}}$ be the “sub-
feature matrix” of $\Phi^{i}$ at state $s$, that is, the $a^{i}$-th row of
$\Phi_{s}^{i}$ is given by $\phi^{i}(s,a^{i})$ for all
$a^{i}\in\mathcal{A}^{i}$. In this case, the softmax policy class of player
$i\in\\{1,2\\}$ can be compactly written as
$\Pi^{i}=\\{\pi_{\theta^{i}}^{i}(s)=\sigma_{\tau}(\Phi_{s}^{i}\theta^{i}),\,\,\forall\,s\in\mathcal{S}\mid\theta^{i}\in\mathbb{R}^{d_{i}}\\}$,
where we recall our notation that
$\pi_{\theta^{i}}^{i}(s)\in\Delta(\mathcal{A}^{i})$.

### 3.1 Two-Timescale $Q$-Learning with Function Approximation

For $i\in\\{1,2\\}$, our two-timescale $Q$-learning with function
approximation is presented in Algorithm 1. Note that Algorithm 1 has two loops
and maintains $4$ sets of iterates $\\{w_{t,k}^{i}\\}$,
$\\{\theta_{t,k}^{i}\\}$, $\\{\bar{w}_{t}^{i}\\}$, and
$\\{\bar{\theta}_{t}\\}$. In Algorithm 1 Line $6$, we use $\lceil\cdot\rceil$
to represent the truncation operator with truncation radius $r=1/(1-\gamma)$.
That is, for any $x\in\mathbb{R}$, we have $\lceil x\rceil=x$ when
$x\in[-r,-r]$, $\lceil x\rceil=r$ when $x>r$, and $\lceil x\rceil=-r$ when
$x<-r$. For a vector-valued $x$, the operator $\lceil\cdot\rceil$ operates
component-wisely. In Algorithm 1 Line 7, we use
$\text{Proj}_{M}^{i}:\mathbb{R}^{d_{i}}\mapsto\mathbb{R}^{d_{i}}$ to represent
the projection operator onto the $\ell_{2}$-norm ball with radius $M$ centered
at the origin, that is,
$\text{Proj}_{M}^{i}(x)=\arg\min_{y\in\mathbb{R}^{d_{i}},\|y\|_{2}\leq
M}\|x-y\|_{2}$, where the projection radius $M\geq 0$ is a tunable parameter.
Both the truncation and the projection are introduced for the stability of our
learning dynamics, which are common in the existing literature even in the
single-agent setting [28, 83, 81, 126].

Algorithm 1 Two-Timescale $Q$-Learning with Function Approximation

1: Input: Integers $K$ and $T$, initializations
$w_{0,0}^{i}=\theta_{0,0}^{i}=\bar{w}_{0}^{i}=\bar{\theta}_{0}^{i}=0$, and
$S_{0}\in\mathcal{S}$

2: for $t=0,1,\cdots,T-1$ do

3: for $k=0,1,\cdots,K-1$ do

4:
$\theta^{i}_{t,k+1}=\theta_{t,k}^{i}+\beta_{k}(w_{t,k}^{i}-\theta_{t,k}^{i})$
and compute
$\pi_{t,k+1}^{i}(S_{k})=\sigma_{\tau}(q_{\theta_{t,k+1}^{i}}^{i}(S_{k}))$

5: Play $A_{k}^{i}\sim\pi_{t,k+1}^{i}(\cdot\mid S_{k})$ (against $A_{k}^{-i}$)
and observe $S_{k+1}\sim p(\cdot\mid S_{k},A_{k}^{i},A_{k}^{-i})$

6:
$\Delta_{t,k}^{i}=R_{i}(S_{k},A_{k}^{i},A_{k}^{-i})+\gamma\sigma_{\tau}(q_{\bar{\theta}_{t}^{i}}^{i}(S_{k+1}))^{\top}\lceil
q_{\bar{w}_{t}^{i}}^{i}(S_{k+1})\rceil-q_{w_{t,k}^{i}}^{i}(S_{k},A_{k}^{i})$

7:
$w_{t,k+1}^{i}=\text{Proj}_{M}^{i}\big{(}w_{t,k}^{i}+\alpha_{k}\nabla\big{[}q_{w_{t,k}^{i}}^{i}(S_{k},A_{k}^{i})\big{]}\Delta_{t,k}^{i}\big{)}$

8: end for

9: $\bar{w}_{t+1}^{i}=w_{t,K}^{i}$ and
$\bar{\theta}_{t+1}^{i}=\theta_{t,K}^{i}$

10: $w_{t+1,0}^{i}=w_{t,K}^{i}$, $\theta_{t+1,0}^{i}=\theta_{t,K}^{i}$, and
$S_{0}=S_{K}$

11: end for

Before diving into the details of Algorithm 1, we see that the learning
dynamics are payoff-based and independent because agent $i$ only needs to
observe the current state of the environment and the realized payoff to
implement the learning dynamics (cf. Algorithm 1 Line $5$). In addition, the
learning dynamics are symmetric between the two players and are driven by a
single trajectory of Markovian samples. Therefore, Algorithm 1 can be
implemented in an online manner without the need to access a generative model
or constantly reset the system.

##### Algorithm Details.

For a fixed outer-loop index $t$, consider the update equations in the inner
loop of Algorithm 1. Note that $\bar{w}_{t}^{i}$ and $\bar{\theta}_{t}$ are
fixed in the inner loop. The iterate $\theta_{t+1,k}^{i}$ is obtained by
simply taking a convex combination (with parameter $\beta_{k}$) between its
previous iterate $\theta_{t,k}^{i}$ and $w_{t,k}^{i}$ (cf. Algorithm 1 Line
$4$). The parameter $\beta_{k}$ is also called the stepsize, or learning rate.
The iterate $w_{t,k}^{i}$ is updated according to a projected SGD (cf.
Algorithm 1 Line $7$) for minimizing a variant of the Bellman error, which
will be elaborated in detail in Section 3.2. In the special case of linear
function approximation, since $\nabla
q_{w_{t,k}^{i}}^{i}(S_{k},A_{k}^{i})=\phi^{i}(S_{k},A_{k}^{i})$, Algorithm 1
Line $7$ is in spirit to the TD-learning in RL [102], where $\Delta_{t,k}^{i}$
is known as the temperal difference and $\alpha_{k}$ is the stepsize. We
require $\beta_{k}\ll\alpha_{k}$ so that $\theta_{t,k}^{i}$ evolves at a much
slower rate compared to that of $w_{t,k}^{i}$, which is why we call Algorithm
1 two-timescale $Q$-learning.

In the outer loop of Algorithm 1, we update $\bar{w}_{t}^{i}$ and
$\bar{\theta}_{t}^{i}$ by synchronizing them to the values of $w_{t,K}^{i}$
and $\theta_{t,K}^{i}$ from the last inner loop. Following the terminology
used in the existing literature [81], we call $\bar{w}_{t}^{i}$ and
$\bar{\theta}_{t}^{i}$ the target network parameters. Note that in single-
agent $Q$-learning, there is only one set of target network parameters, which
corresponds to $\\{\bar{w}_{t}^{i}\\}$ in Algorithm 1. Here, we have another
set of target network parameters $\\{\bar{\theta}_{t}\\}$, which we use to
compute the policies and are designed for the setting of zero-sum stochastic
games. In Algorithm 1 Line $10$, we set the initializations of the iterates
for the next inner loop to be their corresponding last iterates of the
previous inner loop, ensuring that Algorithm 1 follows a single trajectory of
Markovian samples.

##### Connection to the DQN.

Suppose that we use the deep neural net as a means for function approximation.
Then the main difference between Algorithm 1 and the DQN [81] is the two-
timescale structure. In fact, if we do not choose $\beta_{k}\ll\alpha_{k}$ but
instead do the opposite by setting $\beta_{k}\equiv 1$. Since every time we
simply assign the fast-timescale iterate $w_{t,k}^{i}$ to the slow-timescale
iterate $\theta_{t,k}^{i}$ (cf. Algorithm 1 Line $4$), there are essentially
only two sets of iterates $\\{w_{t,k}^{i}\\}$ and $\bar{w}_{t}^{i}$. In this
case, Algorithm 1 reduces to the DQN in [81], with modular changes such as
replacing the softmax by the hardmax and the use of experience replay for
efficient sampling.

######  Remark.

In light of the discussion in the previous paragraph, the main algorithmic
novelty of Algorithm 1 is essentially that we introduced another set of
parameters $\theta_{t,k}^{i}$ (along with its target network
$\bar{\theta}_{t}^{i}$) to compute the policies instead of directly taking
softmax based on the local $q$-function associated with $w_{t,k}^{i}$.
Although the update equation for $\theta_{t,k}^{i}$ is a simple convex
combination, we illustrate in Section 5 that such an update equation resembles
the smoothed best-response dynamics in zero-sum games, which enables us to
show the last-iterate convergence of Algorithm 1 in the special case of linear
function approximation.

### 3.2 The Viewpoint of Algorithm 1 as the Minimax Value Iteration

In single-agent $Q$-learning with target networks, it was shown in the
literature that the algorithm mimics the value iteration [83, 28]. In this
section, we also elaborate on how to view Algorithm 1 as an approximate
version of the minimax value iteration, which is a well-known algorithm for
solving zero-sum stochastic games. This also serves as a starting point for
our analysis as will be presented in Section 5. We begin by introducing the
minimax value iteration to provide context.

##### The Minimax Value Iteration.

Given $i\in\\{1,2\\}$, let
$\mathcal{B}^{i}:\mathbb{R}^{|\mathcal{S}|}\mapsto\mathbb{R}^{|\mathcal{S}|}$
be the minimax Bellman operator defined as
$\mathcal{B}^{i}(v^{i})(s)=\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}$
for all $s\in\mathcal{S}$, where the operator
$\mathcal{T}^{i}:\mathbb{R}^{|\mathcal{S}|}\mapsto\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}||\mathcal{A}^{-i}|}$
is defined as

$\displaystyle\mathcal{T}^{i}(v^{i})(s,a^{i},a^{-i})=R_{i}(s,a^{i},a^{-i})+\gamma\mathbb{E}\left[v^{i}(S_{1})\mid
S_{0}=s,A_{0}^{i}=a^{i},A_{0}^{-i}=a^{-i}\right]$ (1)

for all $v^{i}\in\mathbb{R}^{|\mathcal{S}|}$ and $(s,a^{i},a^{-i})$. Here, we
recall our notation that $\mathcal{T}^{i}(v^{i})(s)$ is an
$|\mathcal{A}^{i}|\times|\mathcal{A}^{-i}|$ matrix with its
$(a^{i},a^{-i})$-th entry being $\mathcal{T}^{i}(v^{i})(s,a^{i},a^{-i})$. In
the minimax value iteration, Player $i$ initializes
$v_{0}^{i}\in\mathbb{R}^{|\mathcal{S}|}$ arbitrarily and iteratively updates
$v_{t}^{i}$ according to

$\displaystyle
v_{t+1}^{i}(s)=\mathcal{B}^{i}(v_{t}^{i})(s),\quad\forall\,s\in\mathcal{S}.$
(2)

It is well-known that the operator $\mathcal{B}^{i}(\cdot)$ is a contraction
mapping [72], hence having a unique fixed-point
$v_{*}^{i}\in\mathbb{R}^{|\mathcal{S}|}$ [4]. As a result, the minimax value
iteration converges geometrically fast to $v_{*}^{i}$. In addition, given
$i\in\\{1,2\\}$, a joint policy $\pi=(\pi^{1},\pi^{2})$ satisfying
$\pi^{i}(s)\mathcal{T}^{i}(v_{*}^{i})(s)\pi^{-i}(s)=\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{*}^{i})(s)\mu^{-i}$
is a Nash equilibrium of the zero-sum stochastic game. Although the minimax
value iteration is provably efficient, it is not a payoff-based independent
learning algorithm, as each agent needs to know the underlying model of the
stochastic game (i.e., the transition kernel and the reward function) to
implement the algorithm. Next, we illustrate why Algorithm 1 can be viewed as
an approximate version of the minimax value iteration, while only using the
realized payoffs for implementation.

We begin by reformulating Algorithm 1. Observe that the target network
parameters $\bar{\theta}_{t}^{i}$ and $\bar{w}_{t}^{i}$ are only used in
Algorithm 1 Line $7$ and are synchronized to $\theta_{t,K}^{i}$ and
$w_{t,K}^{i}$ in the outer loop. Therefore, letting
$v_{t+1}^{i}(s)=\sigma_{\tau}(q_{\theta_{t,K}^{i}}^{i}(s))^{\top}\lceil
q_{\theta_{t,K}^{i}}^{i}(s)\rceil$ for all $s\in\mathcal{S}$ and $t$, we can
rewrite Algorithm 1 in an equivalent form in Algorithm 2. Note that Algorithm
2 is not meant for implementation but only used for the analysis.

Algorithm 2 Equivalent Formulation of Algorithm 1

1: Input: Integers $K$ and $T$, initializations
$w_{0,0}^{i}=\theta_{0,0}^{i}=0$, $v_{0}^{i}=0$, and $S_{0}\in\mathcal{S}$

2: for $t=0,1,\cdots,T-1$ do

3: for $k=0,1,\cdots,K-1$ do

4:
$\theta^{i}_{t,k+1}=\theta_{t,k}^{i}+\beta_{k}(w_{t,k}^{i}-\theta_{t,k}^{i})$
and $\pi_{t,k+1}^{i}(s)=\sigma_{\tau}(q_{\theta_{t,k+1}^{i}}^{i}(s))$ for all
$s\in\mathcal{S}$

5: Play $A_{k}^{i}\sim\pi_{t,k+1}^{i}(\cdot\mid S_{k})$ (against $A_{k}^{-i}$)
and observe $S_{k+1}\sim p(\cdot\mid S_{k},A_{k}^{i},A_{k}^{-i})$

6:
$w_{t,k+1}^{i}=\text{Proj}_{M}^{i}\big{(}w_{t,k}^{i}+\alpha_{k}\nabla\big{[}q_{w_{t,k}^{i}}^{i}(S_{k},A_{k}^{i})\big{]}(R_{i}(S_{k},A_{k}^{i},A_{k}^{-i})+\gamma
v_{t}^{i}(S_{k+1})-q_{w_{t,k}^{i}}^{i}(S_{k},A_{k}^{i}))\big{)}$

7: end for

8: $v_{t+1}^{i}(s)=\pi_{t,K}^{i}(s)^{\top}\lceil q_{w_{t,K}^{i}}^{i}(s)\rceil$
for all $s\in\mathcal{S}$

9: $w_{t+1,0}^{i}=w_{t,K}^{i}$, $\theta_{t+1,0}^{i}=\theta_{t,K}^{i}$, and
$S_{0}=S_{K}$

10: end for

To interpret Algorithm 2 as the minimax value iteration, in view of Eq. (2)
and Algorithm 2 Line $8$, we would like to have

$\displaystyle q_{w_{t,K}^{i}}^{i}(s)\approx\;$
$\displaystyle\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s),\quad\forall\,s\in\mathcal{S},$
(3)
$\displaystyle\pi^{i}_{t,K}(s)^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)\approx\;$
$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\mu^{-i},\quad\forall\,s\in\mathcal{S}.$
(4)

Note that Eq. (4) means that the joint policy (at state $s$)
$(\pi_{t,K}^{i}(s),\pi_{t,K}^{-i}(s))$ is an approximate Nash equilibrium of
the matrix game with the payoff matrix being $\mathcal{T}^{i}(v_{t}^{i})(s)$.

##### The Fast-Timescale Iterates.

We first elaborate on how Eq. (3) can be achieved with Algorithm 2 Line $6$,
which updates the fast-timescale parameter $w_{t,k}^{i}$. Given
$i\in\\{1,2\\}$ and $t\geq 0$, for any
$v_{t}^{i}\in\mathbb{R}^{|\mathcal{S}|}$ and joint policy
$(\pi_{t}^{1},\pi_{t}^{2})$ that are not time-varying in $k$, consider
minimizing the objective function

$\displaystyle J(w^{i})=\;$
$\displaystyle\sum_{s\in\mathcal{S}}\mu_{\pi_{t}}(s)\sum_{a^{i}\in\mathcal{A}^{i}}\pi_{t}^{i}(a^{i}|s)\left(\sum_{a^{-i}\in\mathcal{A}^{-i}}\mathcal{T}^{i}(v_{t}^{i})(s,a^{i},a^{-i})\pi_{t}^{-i}(a^{-i}|s)-q_{w^{i}}^{i}(s,a^{i})\right)^{2}$

using the projected SGD, where $\mu_{\pi_{t}}\in\mathbb{R}^{|\mathcal{S}|}$
represents the stationary distribution of the Markov chain $\\{S_{k}\\}$
induced by the joint policy $(\pi_{t}^{1},\pi_{t}^{2})$, assuming the uniform
ergodicity of $\\{S_{k}\\}$. Specifically, with a sample
$(S_{k},A_{k}^{1},A_{k}^{2},S_{k+1})$ such that
$S_{k}\sim\mu_{\pi_{t}}(\cdot)$, $A_{k}^{i}\sim\pi_{t}^{i}(\cdot|s)$ for
$i\in\\{1,2\\}$, and $S_{k+1}\sim p(\cdot\mid S_{k},A_{k}^{1},A_{k}^{2})$, in
view of the explicit expression of the operator $\mathcal{T}^{i}(\cdot)$ (cf.
Eq. (1)), the projected SGD for minimizing $J(\cdot)$ is given as

$\displaystyle
w_{t,k+1}^{i}=\text{Proj}_{M}^{i}\big{(}w_{t,k}^{i}+\alpha_{k}\nabla\big{[}q_{w_{t,k}^{i}}^{i}(S_{k},A_{k}^{i})\big{]}(R_{i}(S_{k},A_{k}^{i},A_{k}^{-i})+\gamma
v_{t}^{i}(S_{k+1})-q_{w_{t,k}^{i}}^{i}(S_{k},A_{k}^{i}))\big{)},$ (5)

where $M$ is the projection radius. When using linear function approximation,
the objective function $J(\cdot)$ is quadratic (hence strongly convex) in
$w^{i}$. Therefore, with appropriately chosen stepsizes and a large enough
$K$, the output $w_{t,K}^{i}$ of the projected SGD in Eq. (5) is guaranteed to
converge to the global optimal solution of $J(\cdot)$. In addition, when the
function approximation class has enough representation power, we would expect
the $q$-function associated with $w_{t,K}^{i}$ to satisfy
$q_{w_{t,K}^{i}}^{i}(s)\approx\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t}^{-i}(s)$
for all $s\in\mathcal{S}$.

Although we motivated the projected SGD update (5) when the joint policy
$(\pi_{t}^{1},\pi_{t}^{2})$ is fixed in updating $w_{t,k}^{i}$, the reasoning
here can be easily extended to the inner loop of Algorithm 2, where
$\pi_{t,k}^{i}$ is time-varying in $k$. Specifically, in the inner loop of
Algorithm 2, observe that the parameters $\\{\theta_{t,k}^{i}\\}$ of the
policies $\\{\pi_{t,k}^{i}\\}$ are updated in a slower timescale compared to
the parameters $\\{w_{t,k}^{i}\\}$ of the local $q$-functions
$\\{q_{w_{t,k}^{i}}^{i}\\}$. Therefore, from the perspective of $w_{t,k}^{i}$,
the policy $\pi_{t,k}^{i}$ is as if it were stationary in $k$. As a result, we
would expect
$q_{w_{t,K}^{i}}^{i}(s)\approx\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)$,
as desired in Eq. (3). To make the claim rigorous, we need to carefully
analyze the two-timescale stochastic approximation in the inner loop of
Algorithm 2, which is done in Appendix A.4.

##### The Slow-Timescale Iterates.

Now that we have elaborated that the fast-timescale iterate in the inner loop
of Algorithm 2 provides us with a $q_{w_{t,K}^{i}}^{i}$ such that
$q_{w_{t,K}^{i}}^{i}(s)\approx\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)$,
it remains to demonstrate that the slow-timescale iterates
$\\{\theta_{t,k}^{i}\\}$ provide us with a joint policy that satisfies Eq.
(4). The problem here is essentially to solve a matrix game (for each state
$s$) with the payoff matrix $\mathcal{T}^{i}(v_{t}^{i})(s)$. To achieve that,
one of the most classical and promising learning dynamics is the smoothed
best-response dynamics [47, 51, 32], which maintain a policy $\pi_{t,k}^{i}$,
and update it iteratively towards the smoothed best response to the opponent’s
latest policy according to

$\displaystyle\pi_{t,k+1}^{i}(s)=\pi_{t,k}^{i}(s)+\beta_{k}(\sigma_{\tau}(\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,k}^{-i}(s))-\pi_{t,k}^{i}(s)),\quad\forall\;s\in\mathcal{S},$
(6)

where $\beta_{k}$ is the stepsize and $\sigma_{\tau}(\cdot)$ is the softmax
operator. Suppose that both players follow the smoothed best-response
dynamics, they are guaranteed to find the Nash distribution [68, 52], which
approximates a Nash equilibrium as the temperature $\tau$ goes to zero. There
are two challenges in implementing Eq. (6). First, it requires the quantity
$\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,k}^{-i}(s)$ (which involves the
opponent’s policy and the underlying model), and second, it operates directly
in the policy space instead of in the parameter space, which violates the
purpose of using function approximation. The first challenge has been overcome
by introducing the fast-timescale iterate $\\{w_{t,k}^{i}\\}$ to estimate the
quantity $\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,k}^{-i}(s)$ via the projected
SGD. To overcome the second challenge, we need to extend the idea of the
smoothed best-response dynamics to the parameter space of the function
approximation class.

In view of Eq. (6), the key insight here is that, in each iteration, the
policy is updated by taking a step towards the smoothed best-response to the
opponent’s latest policy. Inspired by this idea, in Algorithm 2 Line $4$, we
propose to update the policy parameter $\theta_{t,k+1}^{i}$ by taking a convex
combination between the previous policy parameter $\theta_{t,k}^{i}$ and the
fast-timescale iterate $w_{t,k}^{i}$. Recall that the $q$-function associated
with $w_{t,k}^{i}$ is constructed as an estimator of the marginalized payoff
to the opponent $\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,k}^{-i}(s)$. Therefore,
this is in the same spirit as the smoothed best-response dynamics in the
policy space as, in each iteration, the parameter of the policy is updated by
taking a step towards the parameter of the smoothed best response, which is
estimated through the local $q$-function associated with iterate
$w_{t,k}^{i}$.

######  Remark.

Although Algorithm 2 Line $4$ is simple, intuitive, and mimics the smoothed
best-response dynamics in the parameter space, rigorously analyzing it is
particularly challenging. This is because, when using softmax policy, even
under linear function approximation, a linear update in the parameter space
does not imply a linear update in the policy space as in the smoothed best-
response dynamics in Eq. (6). Therefore, in a Lyapunov-based approach, there
is no existing Lyapunov function to study the update equation for the slow-
timescale iterates, the construction of which presents our main technical
novelty. This will be illustrated in detail in Section 5.

## 4 Last-Iterate Finite-Sample Analysis

In the rest of this paper, we consider using linear function approximation.
Denote
$\Pi=\\{(\pi^{1},\pi^{2})\mid\pi^{i}\in\Pi^{i},\,\forall\,i\in\\{1,2\\}\\}$ as
the set of joint policies that are representable using our softmax policy
class. Before presenting our main theoretical results, we first state our
assumptions.

######  Assumption 4.1.

For any joint policy $(\pi^{1},\pi^{2})\in\Pi$, the Markov chain $\\{S_{k}\\}$
induced by $(\pi^{1},\pi^{2})$ is irreducible and aperiodic. In addition,

  1. (1)

there exist $C>0$ and $\rho\in(0,1)$ such that
$\sup_{\pi\in\Pi}\max_{s\in\mathcal{S}}\left\|P_{\pi}(s,\cdot)-\mu_{\pi}(\cdot)\right\|_{\text{TV}}\leq
C\rho^{k}$ for all $k\geq 0$, where
$P_{\pi}\in\mathbb{R}^{|\mathcal{S}|\times|\mathcal{S}|}$ and
$\mu_{\pi}\in\Delta(\mathcal{S})$ are the transition matrix and the stationary
distribution of the Markov chain $\\{S_{k}\\}$ induced by $\pi$, respectively;

  2. (2)

it holds that
$\lambda:=\min_{i\in\\{1,2\\}}\inf_{\pi\in\Pi}\lambda_{\min}((\Phi^{i})^{\top}D_{\pi}^{i}\Phi^{i})>0$,
where
$D_{\pi}^{i}\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}|\times|\mathcal{S}||\mathcal{A}^{i}|}$
is a diagonal matrix with diagonal components
$\\{\mu_{\pi}(s)\pi^{i}(a^{i}|s)\\}_{(s,a^{i})\in\mathcal{S}\times\mathcal{A}^{i}}$
and $\lambda_{\min}(\cdot)$ returns the smallest eigenvalue of a symmetric
matrix.

Assumption 4.1 is standard in the existing literature studying RL algorithms
under time-varying policies [126, 60, 34, 29, 115, 112, 88], which ensures
that all policies encountered from the algorithm trajectory have enough
exploration. Since we use softmax policies and both the fast and the slow-
timescale iterates are uniformly bounded (due to the projection operator
$\text{Proj}_{M}^{i}(\cdot)$), there exists $\ell>0$ such that
$\pi_{t,k}^{i}(a^{i}\mid s)\geq\ell$ for all $(s,a^{i})$, $i\in\\{1,2\\}$, and
$t,k$. Therefore, with some additional work, it can be shown that Assumption
4.1 is satisfied if the following (weaker) assumption is imposed.

######  Assumption 4.2.

There exists a joint policy $(\pi_{b}^{1},\pi_{b}^{2})$ such that the induced
Markov chain $\\{S_{k}\\}$ is irreducible and aperiodic.

The implication from Assumption 4.2 to Assumption 4.1 was formally established
in [32, Lemma 4.1]. However, in that case, the convergence rate of the
learning dynamics involves problem-dependent constants that are implicit
functions of the underlying transition probabilities of the stochastic game.
In this work, for ease of exposition, we directly impose Assumption 4.1. Given
a precision level $\delta>0$, we define

$\displaystyle z_{\delta}=\min\left\\{k\geq
0\,:\,\sup_{\pi\in\Pi}\max_{s\in\mathcal{S}}\left\|P_{\pi}^{k}(s,\cdot)-\mu_{\pi}(\cdot)\right\|_{\text{TV}}\leq\delta\right\\},$
(7)

which can be viewed as the uniform mixing time of the Markov chain
$\\{S_{k}\\}$ induced by any joint policy from the policy class $\Pi$. Under
Assumption 4.1, it is easy to see that
$z_{\delta}=\mathcal{O}(\log(1/\delta))$.

To state the next assumption, we need to define a variant of the inherent
Bellman error. For any joint policy $(\pi^{1},\pi^{2})\in\Pi$ and
$i\in\\{1,2\\}$, let
$\bar{\mathcal{H}}_{\pi}^{i}:\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}|}\mapsto\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}|}$
be an operator defined as

$\displaystyle[\bar{\mathcal{H}}_{\pi}^{i}(q^{i})](s,a^{i})=\sum_{a^{-i}\in\mathcal{A}^{-i}}\left(R_{i}(s,a^{i},a^{-i})+\gamma\mathbb{E}\left[\pi^{i}(S_{1})^{\top}q^{i}(S_{1})\,\middle|\,S_{0}=s,A_{0}^{i}=a^{i}\right]\right)\pi^{-i}(a^{-i}\mid
s)$

for all $(s,a^{i})\in\mathcal{S}\times\mathcal{A}^{i}$ and
$q^{i}\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}|}$, where we recall our
notation that $\pi^{i}(s)\in\Delta(\mathcal{A}^{i})$ and
$q^{i}(s)\in\mathbb{R}^{|\mathcal{A}^{i}|}$ for all $s\in\mathcal{S}$. Define

$\displaystyle\mathcal{E}=\sum_{i=1,2}\sup_{\pi\in\Pi,\tilde{w}^{i}\in\mathbb{R}^{d_{i}}}\inf_{w^{i}\in\mathbb{R}^{d_{i}}}\left\|\Phi^{i}w^{i}-\bar{\mathcal{H}}_{\pi}^{i}(\lceil\Phi^{i}\tilde{w}^{i}\rceil)\right\|_{\infty}.$
(8)

######  Assumption 4.3.

It holds that $\mathcal{E}=0$.

In the special case of tabular setting, since
$\Phi^{i}=I_{|\mathcal{S}||\mathcal{A}^{i}|}$ is the identity matrix,
Assumption 4.3 is automatically satisfied. More generally, Assumption 4.3
means that the function approximation class is rich enough, and can be viewed
as an extension of the Bellman completeness assumption from the single-agent
setting [83, 118, 1] to the setting of solving zero-sum stochastic games with
independent learning.

Next, we state the condition for choosing the stepsizes.

######  Condition 4.1.

We choose $\tau\leq 1/(1-\gamma)$, $\alpha_{0}\leq 1/\lambda$, and
$c_{\alpha,\beta}:=\beta_{k}/\alpha_{k}\leq\min(\frac{1}{64L_{p}A_{\max}},\lambda)$,
where $A_{\max}=\max\\{|\mathcal{A}^{1}|,|\mathcal{A}^{2}|\\}$ and
$L_{p}:=\frac{\log(\rho/C)}{\log(\rho)}+\frac{1}{1-\rho}$.

Now, we present the last-iterate finite-sample analysis of Algorithm 1 when
using constant stepsizes $\alpha_{k}\equiv\alpha$ and
$\beta_{k}\equiv\beta=c_{\alpha,\beta}\alpha$, where
$c_{\alpha,\beta}\in(0,1)$ is the stepsize ratio. The results for using
diminishing stepsizes are straightforward extensions.

######  Theorem 4.1.

Suppose that we use linear function approximation, Assumptions 4.1 and 4.3 are
satisfied, the stepsizes verify Condition 4.1, and the projection radius
$M\geq\lambda^{-1/2}(1-\gamma)^{-1}$, Then, when both players follow the
learning dynamics presented in Algorithm 1, we have for any $T\geq 0$ and
$K\geq z_{\beta}$ that

$\displaystyle\mathbb{E}[\text{NG}(\pi_{T,K}^{1},\pi_{T,K}^{2})]\lesssim\;$
$\displaystyle\underbrace{\frac{A_{\max}^{2}T\gamma^{T-1}}{\tau^{2}(1-\gamma)^{4}}}_{\mathcal{E}_{1}}+\underbrace{\frac{\tau\log(A_{\max})}{(1-\gamma)^{2}}}_{\mathcal{E}_{2}}+\underbrace{\frac{A_{\max}^{3}K\left(1-\beta/2\right)^{\frac{K-z_{\beta}-1}{2}}}{\tau^{3}\lambda(1-\gamma)^{7}}}_{\mathcal{E}_{3,1}}$
$\displaystyle+\underbrace{\frac{A_{\max}^{2}z_{\beta}\alpha^{\frac{1}{2}}}{\tau^{\frac{5}{2}}\lambda^{\frac{5}{4}}(1-\gamma)^{\frac{13}{2}}}+\frac{A_{\max}^{3}\beta}{\tau^{3}\lambda^{2}(1-\gamma)^{7}\alpha}+\frac{A_{\max}^{3}z_{\beta}^{2}\alpha}{\tau^{4}\lambda^{5/2}(1-\gamma)^{7}}+\frac{A_{\max}^{5}\beta^{2}}{\tau^{5}\lambda^{4}(1-\gamma)^{8}\alpha^{2}}}_{\mathcal{E}_{3,2}},$
(9)

where $z_{\beta}$ is the mixing time with precision $\beta$ (see Eq. (7)), and
we use $a\lesssim b$ to mean that there exists an absolute constant $c>0$ such
that $a\leq cb$. As a result, given $\epsilon>0$, to achieve
$\mathbb{E}[\text{NG}(\pi_{T,K}^{1},\pi_{T,K}^{2})]\leq\epsilon$, the sample
complexity is
$\tilde{\mathcal{O}}(\text{poly}(1/\epsilon.1/\lambda,1/(1-\gamma),A_{\max})$.

To the best of our knowledge, Theorem 4.1 presents the first last-iterate
finite-sample analysis of payoff-based independent learning that successfully
incorperates function approximation. We next provide an interpretation for
each term on the RHS of Eq. (4.1).

The term $\mathcal{E}_{1}$ represents the error in the outer loop. As
illustrated in Section 3.2, the outer loop of Algorithm 1 is an approximation
of the minimax value iteration. Therefore, since the minimax value iteration
converges geometrically fast, the first term also converges geometrically in
$T$. The term $\mathcal{E}_{2}$ represents the error due to using the softmax
policy class with a fixed temperature $\tau$. To see why it arises, consider a
special case where the stochastic game has a unique Nash equilibrium, which
consists of a pair of deterministic policies. In this case, due to the
deterministic nature of the Nash equilibrium and the stochastic nature of the
softmax policies, we cannot make the Nash gap vanish with a fixed temperature
$\tau$. As a result, the error term $\mathcal{E}_{2}$ is proportional to
$\tau$, and can be made arbitrarily small by using a small enough $\tau$. The
terms $\mathcal{E}_{3,1}$ and $\mathcal{E}_{3,2}$ on the RHS of Eq. (4.1)
represent the convergence error in the inner loop of Algorithm 1, which is a
two-timescale stochastic approximation algorithm. The term $\mathcal{E}_{3,1}$
goes to zero geometrically as $K$ goes to infinity and the terms in
$\mathcal{E}_{3,2}$ are constants that are functions of the stepsize for the
fast-timescale and the ratio between the two stepsizes.

## 5 Construction of the Lyapunov Function for the Slow-Timescale Iterates

Throughout this section, we will work with Algorithm 2 and linear function
approximation. The overall framework for our proof is a Lyapunov-based
approach. As illustrated in Section 3.2, the outer loop of Algorithm 2 is
designed to approximate the minimax value iteration, and the inner loop is
designed to solve an induced matrix game for each state $s\in\mathcal{S}$
using a two-timescale stochastic approximation algorithm, where the fast-
timescale iterates are updated using TD-learning to estimate the marginalized
payoff to the opponent, and the slow-timescale iterates are updated in spirit
to the smoothed best-response dynamics, albeit in the parameter space of the
function approximation class. Inspired by [32], we construct a Lyapunov
function for each set of the iterates (i.e., $\\{v_{t}^{i}\\}$,
$\\{w_{t,k}^{i}\\}$, and $\\{\theta_{t,k}^{i}\\}$), derive a Lyapunov drift
inequality for each, and finally combine them together to establish the
finite-sample bounds. The main challenge here (see the last remark of Section
3.2 for more details) is to construct a valid Lyapunov function to study the
slow-timescale iterates $\\{\theta_{t,k}^{i}\\}$, which will be the focus of
this section.

Since we are concerning the inner loop of Algorithm 2, we will omit the
subscript $t$. To begin with, due to the fact that $\Phi^{i}$ has linearly
independent columns, Algorithm 2 Line $4$ can be equivalently written as

$\displaystyle q^{i}_{\theta_{k+1}^{i}}(s)=\;$ $\displaystyle
q^{i}_{\theta_{k}^{i}}(s)+\beta_{k}\left(q^{i}_{w_{k}^{i}}(s)-q^{i}_{\theta_{k}^{i}}(s)\right)$
$\displaystyle=\;$ $\displaystyle
q^{i}_{\theta_{k}^{i}}(s)+\beta_{k}\left(\mathcal{T}^{i}(v^{i})(s)\sigma_{\tau}(q^{-i}_{\theta_{k}^{-i}}(s))-q^{i}_{\theta_{k}^{i}}(s)\right)+\beta_{k}\left(q^{i}_{w_{k}^{i}}(s)-\mathcal{T}^{i}(v^{i})(s)\sigma_{\tau}(q^{i}_{\theta_{k}^{-i}}(s)\right)$
(10)

for all $s\in\mathcal{S}$. Recall that we designed the update equation for the
fast-timescale iterates $\\{w_{k}^{i}\\}$ precisely to estimate
$\mathcal{T}^{i}(v^{i})(s)\sigma_{\tau}(q_{\theta_{k}^{-i}}^{-i}(s))=\mathcal{T}^{i}(v^{i})(s)\pi_{k}^{-i}(s)$.
Therefore, the last term in the RHS of Eq. (10) corresponds to the estimation
error of the fast-timescale iterate $w_{k}^{i}$, and can be controlled by
studying Algorithm 2 Line $6$. As a result, to understand the evolution of the
slow-timescale iterates $\\{\theta_{k}^{i}\\}$, it essentially reduces to
studying the following update equation:

$x_{0}^{i}\in\mathbb{R}^{|\mathcal{A}^{i}|},\quad
x_{k+1}^{i}=x_{k}^{i}+\beta_{k}(X_{i}\sigma_{\tau}(x_{k}^{-i})-x_{k}^{i}),\quad
i\in\\{1,2\\},$ (11)

where $X_{i}\in\mathbb{R}^{|\mathcal{A}^{i}|\times|\mathcal{A}^{-i}|}$ for
$i\in\\{1,2\\}$. To make the connection between Eqs. (10) and (11), fix
$s\in\mathcal{S}$, identify $x_{k}^{i}=q_{\theta_{k}^{i}}^{i}(s)$ and
$X_{i}=\mathcal{T}^{i}(v^{i})(s)$ for $i\in\\{1,2\\}$. For the purpose of
studying the slow-timescale iterates, we assume that
$X_{1}+X_{2}^{\top}=0$333In Algorithm 2, we do not necessarily have
$\mathcal{T}^{1}(v_{t}^{1})(s)+(\mathcal{T}^{2}(v_{t}^{2})(s))^{\top}=0$
because $v_{t}^{1}$ and $v_{t}^{2}$ are independently maintained by players
$1$ and $2$, and do not sum up to zero in general. Such a non-zero sum
structure of the induced matrix game in the inner loop has to be taken into
account (see Appendix A), but is not the focus of this section.. Our goal is
to show that $(x_{k}^{1},x_{k}^{2})$ of Eq. (11) converges to
$(x_{*}^{1},x_{*}^{2})$ that satisfies
$x_{*}^{i}=X_{i}\sigma_{\tau}(x_{*}^{-i})$ for $i\in\\{1,2\\}$, in which case
the induced policy pair $(\sigma_{\tau}(x_{*}^{1}),\sigma_{\tau}(x_{*}^{2}))$
is an approximate Nash equilibrium (or the Nash distribution [67]) of the
matrix game with payoff matrices $X_{1}$ and $X_{2}$. In the rest of this
section, we present the roadmap of the construction of our Lyapunov function
to study $(x_{k}^{1},x_{k}^{2})$ generated by Eq. (11).

### 5.1 The Change of Variable

As illustrated in Section 3.2, we view Algorithm 2 Line $4$ as an extension of
the smoothed best-response dynamics implemented in the parameter space of the
function approximation class. To formally make the connection, we state the
smoothed best-response dynamics for solving the zero-sum matrix game with
payoff matrices $X_{1}$ and $X_{2}$ in the following:

$\displaystyle\hat{\pi}_{k+1}^{i}=\hat{\pi}_{k}^{i}+\beta_{k}(\sigma_{\tau}(X_{i}\hat{\pi}_{k}^{-i})-\hat{\pi}_{k}^{i}),\quad
i\in\\{1,2\\}.$ (12)

Next, we present a change of variable to connect Eq. (12) with Eq. (11), which
relies on the following assumption.

######  Assumption 5.1.

The payoff matrix $X_{1}$ is invertible.

######  Remark.

Since $X_{1}+X_{2}^{\top}=0$, $X_{1}$ being invertible also implies $X_{2}$
being invertible. We want to emphasize that Assumption 5.1 is made only for
the illustration of the construction of our Lyapunov function. In fact,
neither Theorem 4.1 nor the final form of our Lyapunov function for studying
$\\{\theta_{k}^{i}\\}$ requires this assumption.

Under Assumption 5.1, by multiplying $X_{i}^{-1}$ on both sides of Eq. (11)
and denoting $\zeta_{k}^{-i}=X_{i}^{-1}x_{k}^{i}$ for $i\in\\{1,2\\}$, we can
write Eq. (11) equivalently as

$\displaystyle\zeta_{k+1}^{i}=\zeta_{k}^{i}+\beta_{k}(\sigma_{\tau}(X_{i}\zeta_{k}^{-i})-\zeta_{k}^{i}),\quad
i\in\\{1,2\\}.$ (13)

Observe that Eq. (13) seems to be identical to the smoothed best-response
dynamics presented in Eq. (12). However, we want to highlight one major
difference. Note that the iterate $\zeta_{k}^{i}$ is obtained from a change of
variable. Therefore, $\zeta_{k}^{i}$ is not necessarily a policy (or a
probability distribution) as in the smoothed best-response dynamics. In view
of Eq. (13), one may ask that, if we initialized $\zeta_{0}^{i}$ in the
probability simplex $\Delta(\mathcal{A}^{i})$, we would be able to show that
$\zeta_{k}^{i}\in\Delta(\mathcal{A}^{i})$ for all $k\geq 0$ using an induction
argument. However, this is also not possible because of the following two
reasons. (i) Recall from Eqs. (10) and (11) that $X_{i}$ corresponds to
$\mathcal{T}^{i}(v^{i})(s)$ and $x_{0}^{i}$ corresponds to
$q_{\theta_{0}^{i}}^{i}(s)$. Therefore, in order to ensure that
$\zeta_{0}^{-i}=X_{i}^{-1}x_{0}^{i}=\mathcal{T}^{i}(v^{i})(s)q_{\theta_{0}^{i}}^{i}(s)\in\Delta(\mathcal{A}^{i})$,
we would require access to the underlying model parameters to compute
$\mathcal{T}^{i}(v^{i})(s)$. This is not possible in payoff-based independent
learning. (ii) Recall that Eq. (11) is a simplified version of Eq. (10) by
ignoring the estimation error of the fast-timescale iterates
$\\{w_{k}^{i}\\}$, which is stochastic by nature. Therefore, even if we
managed to initialize $\theta_{0}^{i}$ such that
$\zeta_{0}^{i}\in\Delta(\mathcal{A}^{i})$, the stochastic error in $w_{k}^{i}$
may drive $\zeta_{k}^{i}$ out of the probability simplex.

In summary, by performing a change of variable, the resulting update equation
(cf. Eq. (13)) is identical in form to the smoothed best-response dynamics.
However, there are two issues:

  1. (1)

Issue 1. the iterates are not necessarily probability distributions;

  2. (2)

Issue 2. the change of variable requires the invetibility of $X_{1}$, which
may not hold in general.

In the next two subsections, we will start with the Lyapunov function for
studying the smoothed best-response dynamics, i.e., the regularized Nash gap,
and work our way back to construct the Lyapunov function for studying
$(x_{k}^{1},x_{k}^{2})$ generated by Eq. (11), while resolving the two issues
stated above.

### 5.2 The Moreau Envelope for a Smooth Domain Extension

In the existing literature, to study the smooth best-response dynamics (12),
the regularized Nash gap
$\text{RNG}_{\tau}:\Delta(\mathcal{A}^{1})\times\Delta(\mathcal{A}^{2})\mapsto\mathbb{R}$
defined as

$\displaystyle\text{RNG}_{\tau}(\hat{\pi}^{1},\hat{\pi}^{2})=\sum_{i=1,2}\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\\{(\mu^{i}-\hat{\pi}^{i})^{\top}X_{i}\hat{\pi}^{-i}+\tau\nu(\mu^{i})-\tau\nu(\hat{\pi}^{i})\\},\quad\forall\,(\hat{\pi}^{1},\hat{\pi}^{2}),$
(14)

has been shown to be a valid Lyapunov function [52], where
$\nu(\mu^{i})=-\sum_{a^{i}}\mu^{i}(a^{i})\log(\mu^{i}(a^{i}))$ is the entropy
function.

Naturally, we want to use $\text{RNG}_{\tau}(\cdot,\cdot)$ to study
$(\zeta_{k}^{1},\zeta_{k}^{2})$ generated by Eq. (13). However, since
$\zeta_{k}^{i}$ is not necessarily a probability distribution, the entropy
function of $\zeta_{k}^{i}$ is not even well-defined, which prevents us from
directly using $\text{RNG}_{\tau}(\cdot,\cdot)$. Therefore, we need to extend
the domain of the regularized Nash gap from the joint probability simplex
$\Delta(\mathcal{A}^{1})\times\Delta(\mathcal{A}^{2})$ to the joint real
vector space
$\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}$.

To provide a starting point, we first present a naive way to extend the domain
of the regularized Nash gap. For $i\in\\{1,2\\}$, let
$\bar{V}^{i}:\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}\mapsto\mathbb{R}$
be defined as

$\displaystyle\bar{V}^{i}(\zeta^{1},\zeta^{2})=\begin{dcases}\max_{u^{i}\in\Delta(\mathcal{A}^{i})}\\{(u^{i}-\zeta^{i})^{\top}X_{i}\zeta^{-i}+\tau\nu(u^{i})-\tau\nu(\zeta^{i})\\},&\zeta^{i}\in\Delta(\mathcal{A}^{i}),\\\
+\,\infty,&\text{Otherwise},\end{dcases}$ (15)

and define
$\bar{V}:\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}\mapsto\mathbb{R}$
as $\bar{V}=\bar{V}^{1}+\bar{V}^{2}$. Note that the function value of
$\bar{V}$ was directly assigned to infinity if the variable
$(\zeta^{1},\zeta^{2})$ is not in the joint probability simplex
$\Delta(\mathcal{A}^{1})\times\Delta(\mathcal{A}^{2})$. However, the function
$\bar{V}(\cdot,\cdot)$ fails as a Lyapunov function due to its non-smoothness.
What we truly need is to smoothly extend the domain of the regularized Nash
gap. To achieve that, we use the Moreau envelope as a means for the
construction, which is defined in the following.

######  Definition 5.1 (The Moreau Envelope).

Let $h:\mathbb{R}^{d}\mapsto\mathbb{R}$ be a closed and convex function and
$c>0$. The Moreau envelope of $h(\cdot)$ is defined to be the function
$M_{h}^{c}(x)=\inf_{u\in\mathbb{R}^{d}}\\{h(u)+\frac{1}{2c}\|x-u\|_{2}^{2}\\}$.

######  Remark.

The Moreau envelope has many nice properties, one of which is the smoothness
[8]. It was previously used in [50, 9] to study convex but non-smooth
optimization problems, and used in [31, 59] to construct smooth Lyapunov
functions.

Next, we use the Moreau envelope to smoothly extend the domain of the
regularized Nash gap. We start with the function $\bar{V}^{i}(\cdot,\cdot)$
defined in Eq. (15). For any $i\in\\{1,2\\}$ and $\mu>0$, let
$\hat{V}^{i}:\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}\mapsto\mathbb{R}$
be defined as

$\displaystyle\hat{V}^{i}(\zeta^{1},\zeta^{2})=\;$
$\displaystyle\inf_{\bar{u}^{i}\in\mathbb{R}^{|\mathcal{A}^{i}|}}\left\\{\bar{V}^{i}(\bar{u}^{i},\zeta^{-i})+\frac{1}{2\mu}\|\zeta^{i}-\bar{u}^{i}\|_{2}^{2}\right\\}$
$\displaystyle=\;$
$\displaystyle\min_{\bar{u}^{i}\in\Delta(\mathcal{A}^{i})}\left\\{\max_{u^{i}\in\Delta(\mathcal{A}^{i})}\\{(u^{i}-\bar{u}^{i})^{\top}X_{i}\zeta^{-i}+\tau\nu(u^{i})-\tau\nu(\bar{u}^{i})\\}+\frac{1}{2\mu}\|\zeta^{i}-\bar{u}^{i}\|_{2}^{2}\right\\},$
(16)

where Eq. (16) follows from the definition of $\bar{V}^{i}(\cdot,\cdot)$ and
the Weierstrass extreme value theorem. Let
$\hat{V}:\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}\mapsto\mathbb{R}$
be defined as

$\displaystyle\hat{V}(\zeta^{1},\zeta^{2})=\;$
$\displaystyle\sum_{i=1,2}\hat{V}^{i}(\zeta^{1},\zeta^{2})$ $\displaystyle=\;$
$\displaystyle\sum_{i=1,2}\min_{\bar{u}^{i}\in\Delta(\mathcal{A}^{i})}\left\\{\max_{u^{i}\in\Delta(\mathcal{A}^{i})}\\{(u^{i}-\bar{u}^{i})^{\top}X_{i}\zeta^{-i}+\tau\nu(u^{i})-\tau\nu(\bar{u}^{i})\\}+\frac{1}{2\mu}\|\zeta^{i}-\bar{u}^{i}\|_{2}^{2}\right\\},$
(17)

which successfully extends the domain of the regularized Nash gap from the
joint probability simplex
$\Delta(\mathcal{A}^{1})\times\Delta(\mathcal{A}^{2})$ to the joint real
vector space
$\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}$. In fact,
one can show (following the same line of analysis as in Appendix A.4.1) that
$\hat{V}(\cdot,\cdot)$ is a valid Lyapunov function for
$(\zeta_{k}^{1},\zeta_{k}^{2})$ from Eq. (13).

### 5.3 Overcoming the Invertibility Issue by Generalizing the Moreau Envelope

Now that we have constructed a Lyapunov function
$\hat{V}(\zeta_{1},\zeta^{2})$ for $(\zeta_{k}^{1},\zeta_{k}^{2})$ generated
by Eq. (13), it remains to perform the inverse change of variable to finish
constructing the Lyapunov function to study $(x_{k}^{1},x_{k}^{2})$ generated
by Eq. (11). Since $\zeta_{k}^{-i}=X_{i}^{-1}x_{k}^{i}$ for $i\in\\{1,2\\}$,
let
$\tilde{V}:\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}\mapsto\mathbb{R}$
be defined as

$\displaystyle\tilde{V}(x^{1},x^{2})=\;$
$\displaystyle\hat{V}(X_{2}^{-1}x_{2},X_{1}^{-1}x_{1})$ $\displaystyle=\;$
$\displaystyle\sum_{i=1,2}\min_{\bar{u}^{i}\in\Delta(\mathcal{A}^{i})}\max_{u^{i}\in\Delta(\mathcal{A}^{i})}\left\\{(u^{i}-\bar{u}^{i})^{\top}x^{i}+\tau\nu(u^{i})-\tau\nu(\bar{u}^{i})+\frac{1}{2\mu}\|X_{-i}^{-1}x^{-i}-\bar{u}^{i}\|_{2}^{2}\right\\},$
(18)

which would serve as a valid Lyapunov function to study
$(x_{k}^{1},x_{k}^{2})$ from Eq. (11) if Assumption 5.1 were satisfied, that
is, $X_{1}$ (or equivalently $X_{2}$) were indeed invertible. However, this
may not be true in general. For example, when players $1$ and $2$ have
different numbers of actions, the matrix $X_{1}$ (or $X_{2}$) is not even a
square matrix, hence can not be invertible. Therefore, the last challenge we
need to overcome is to remove the invertibility assumption.

Note that $\\{X_{i}^{-1}\\}_{i\in\\{1,2\\}}$ only appear inside the norm
square on the RHS of Eq. (18), which arises due to using the Moreau envelope
as a means for smooth approximation. Tracing back to Definition 5.1, the
Moreau envelope is nothing but an infimal convolution between the function
$h(\cdot)$ and the $\ell_{2}$-norm square function. We next generalize the
definition of the Moreau envelope by replacing $\|\cdot\|_{2}$ with a properly
defined seminorm, the advantage of which is that, after performing the inverse
change of variable as in Eq. (18), the inverse of $X_{i}$ vanishes as if we
never needed the invertibility of $X_{i}$.

Given a matrix $W\in\mathbb{R}^{d_{2}\times d_{1}}$, let $\|\cdot\|_{W}$ be a
function in $\mathbb{R}^{d_{1}}$ defined as
$\|x\|_{W}=\sqrt{x^{\top}W^{\top}Wx}$. Note that, unless $W$ has independent
columns, $\|\cdot\|_{W}$ is a seminorm instead of a norm.

######  Definition 5.2 (The Generalized Moreau Envelope).

Let $h:\mathbb{R}^{d_{1}}\mapsto\mathbb{R}$ be a closed and convex function,
$W\in\mathbb{R}^{d_{2}\times d_{1}}$ be a matrix, and $c>0$. The generalized
Moreau envelope of $h(\cdot)$ is defined to be the function
$M_{h}^{W,c}(x)=\inf_{u\in\mathbb{R}^{d}}\\{h(u)+\frac{1}{2c}\|x-u\|_{W}^{2}\\}$.

Now, starting from Eq. (16), instead of using the standard Moreau envelope
(cf. Definition 5.1) as a means for smoothly extending the domain of the
regularized Nash gap, we use the generalized Moreau envelope with $W=X_{-i}$
for defining $\hat{V}^{i}(\cdot,\cdot)$. In this case, after performing the
inverse change of variable, we arrive at the final Lyapunov function for
studying $(x_{k}^{1},x_{k}^{2})$ generated by Eq. (11):

$\displaystyle V(x^{1},x^{2})=\;$
$\displaystyle\sum_{i=1,2}\min_{\bar{u}^{i}\in\Delta(\mathcal{A}^{i})}\max_{u^{i}\in\Delta(\mathcal{A}^{i})}\left\\{(u^{i}-\bar{u}^{i})^{\top}x^{i}+\tau\nu(u^{i})-\tau\nu(\bar{u}^{i})+\frac{1}{2\mu}\|X_{-i}^{-1}x^{-i}-\bar{u}^{i}\|_{X_{-i}}^{2}\right\\}$
$\displaystyle=\;$
$\displaystyle\sum_{i=1,2}\min_{\bar{u}^{i}\in\Delta(\mathcal{A}^{i})}\max_{u^{i}\in\Delta(\mathcal{A}^{i})}\left\\{(u^{i}-\bar{u}^{i})^{\top}x^{i}+\tau\nu(u^{i})-\tau\nu(\bar{u}^{i})+\frac{1}{2\mu}\|x^{-i}-X_{-i}\bar{u}^{i}\|_{2}^{2}\right\\}.$
(19)

Importantly, while we used the invertibility of $X_{1}$ to present the roadmap
in constructing our Lyapunov function, since $X_{1}^{-1}$ (or $X_{2}^{-1}$)
does not appear in the definition of $V(\cdot,\cdot)$; see Eq. (19), it is as
if we never needed such invertibility.

It is easy to see that $V(x^{1},x^{2})$ is non-negative by definition. In
addition, we have $V(x^{1},x^{2})=0$ if and only if
$x^{i}=X_{i}\sigma_{\tau}(x^{-i})$ for $i\in\\{1,2\\}$, as desired. To see
this, observe that when $x^{i}=X_{i}\sigma_{\tau}(x^{-i})$ for
$i\in\\{1,2\\}$, $\sigma_{\tau}(x^{i})$ is simultaneously the optimal solution
of both
${\max}_{u^{i}\in\Delta(\mathcal{A}^{i})}\left\\{(u^{i})^{\top}x^{i}+\tau\nu(u^{i})\right\\}$
and
${\min}_{\bar{u}^{i}\in\Delta(\mathcal{A}^{i})}\\{-(\bar{u}^{i})^{\top}x^{i}-\tau\nu(\bar{u}^{i})+\frac{1}{2\mu}\|x^{-i}-X_{-i}\bar{u}^{i}\|_{2}^{2}\\}$
for $i\in\\{1,2\\}$. The following proposition shows that
$(x_{k}^{1},x_{k}^{2})$ updated according to Eq. (11) produces a negative
drift with respect to $V(\cdot,\cdot)$.

######  Proposition 5.1 (Proof in Appendix A).

Consider $(x_{k}^{1},x_{k}^{2})$ generated by Eq. (11). By choosing
$\mu=\tau/64$, we have for all $k\geq 0$ that

$\displaystyle V(x_{k+1}^{1},x_{k+1}^{2})\leq\;$
$\displaystyle\left(1-\frac{\beta_{k}}{2}\right)V(x_{k}^{1},x_{k}^{2})+\frac{138\max_{i\in\\{1,2\\}}(\|x_{0}^{i}\|_{2}+2\|X_{i}\|_{1,2})}{\tau}\beta_{k}^{2},$
(20)

where $\|X_{i}\|_{1,2}=\max_{\|x\|_{1}=1}\|X_{i}x\|_{2}$ is an induced matrix
norm.

The first term on the RHS of Eq. (20) represents the negative drift of
$(x_{k}^{1},x_{k}^{2})$ with respect to the Lyapunov function
$V(\cdot,\cdot)$, and the second term is a higher order error term, which is
due to the fact that Eq. (11) is a discrete update equation. To this end, we
have successfully constructed a Lyapunov function to study
$(x_{k}^{1},x_{k}^{2})$ generated by Eq. (11), and therefore the slow-
timescale iterates $\\{\theta_{t,k}^{i}\\}$ of Algorithm 2.

## 6 Conclusion

In this work, we consider two-player zero-sum stochastic games and develop a
two-timescale $Q$-learning that is payoff-based and uses function
approximation. In the special case of the linear function approximation, our
proposed learning dynamics provably enjoy last-iterate finite-sample
guarantees. To establish the result, we use a Lyapunov approach. The key
novelty lies in the construction of a valid Lyapunov function for the slow-
timescale iterates, which involves a change of variable and the use of a
generalized variant of the Moreau envelope of the regularized Nash gap.

We recognize $3$ future directions of this work. Firstly, we would like to
remove the projection operator $\text{Proj}_{M}^{i}(\cdot)$ in Algorithm 1 as
the projection radius depends on unknown parameters of the stochastic game,
which is not ideal from a practical viewpoint. Secondly, while this is the
first result that establishes last-iterate finite-sample bounds for payoff-
based independent learning dynamics under function approximation, we believe
the convergence rate is not tight, and improving the rate is another future
direction. Finally, we would like to see if the algorithmic ideas and
technical tools developed in this work, especially the Lyapunov function
construction process illustrated in Section 5, can be applied to other MARL
settings.

## References

  * Agarwal et al., [2021] Agarwal, N., Chaudhuri, S., Jain, P., Nagaraj, D. M., and Netrapalli, P. (2021).  Online target $Q$-learning with reverse experience replay: Efficiently finding the optimal policy for linear MDPs.  In International Conference on Learning Representations.
  * Bai et al., [2020] Bai, Y., Jin, C., and Yu, T. (2020).  Near-optimal reinforcement learning with self-play.  Advances in Neural Information Processing Systems, 33.
  * Baird, [1995] Baird, L. (1995).  Residual algorithms: Reinforcement learning with function approximation.  In Machine Learning Proceedings 1995, pages 30–37. Elsevier.
  * Banach, [1922] Banach, S. (1922).  Sur les opérations dans les ensembles abstraits et leur application aux équations intégrales.  Fund. math, 3(1):133–181.
  * Bansal and Gupta, [2019] Bansal, N. and Gupta, A. (2019).  Potential-function proofs for gradient methods.  Theory of Computing, 15(1):1–32.
  * [6] Baudin, L. and Laraki, R. (2022a).  Fictitious play and best-response dynamics in identical interest and zero-sum stochastic games.  In International Conference on Machine Learning, pages 1664–1690. PMLR.
  * [7] Baudin, L. and Laraki, R. (2022b).  Smooth Fictitious Play in Stochastic Games with Perturbed Payoffs and Unknown Transitions.  In Advances in Neural Information Processing Systems.
  * Beck, [2017] Beck, A. (2017).  First-order methods in optimization, volume 25.  SIAM.
  * Beck and Teboulle, [2012] Beck, A. and Teboulle, M. (2012).  Smoothing and first order methods: A unified framework.  SIAM Journal on Optimization, 22(2):557–580.
  * Beck and Srikant, [2012] Beck, C. L. and Srikant, R. (2012).  Error bounds for constant step-size $Q$-learning.  Systems & control letters, 61(12):1203–1208.
  * Beck and Srikant, [2013] Beck, C. L. and Srikant, R. (2013).  Improved upper bounds on the expected error in constant step-size $Q$-learning.  In 2013 American Control Conference, pages 1926–1931. IEEE.
  * Benaïm, [2006] Benaïm, M. (2006).  Dynamics of stochastic approximation algorithms.  In Seminaire de Probabilites XXXIII, pages 1–68. Springer.
  * Benveniste et al., [2012] Benveniste, A., Métivier, M., and Priouret, P. (2012).  Adaptive algorithms and stochastic approximations, volume 22.  Springer Science & Business Media.
  * Bertsekas and Tsitsiklis, [1996] Bertsekas, D. P. and Tsitsiklis, J. N. (1996).  Neuro-dynamic programming.  Athena Scientific.
  * Bhandari et al., [2018] Bhandari, J., Russo, D., and Singal, R. (2018).  A Finite Time Analysis of Temporal Difference Learning With Linear Function Approximation.  In Conference On Learning Theory, pages 1691–1692.
  * Borkar, [2009] Borkar, V. S. (2009).  Stochastic approximation: a dynamical systems viewpoint, volume 48.  Springer.
  * Borkar and Meyn, [2000] Borkar, V. S. and Meyn, S. P. (2000).  The ODE method for convergence of stochastic approximation and reinforcement learning.  SIAM Journal on Control and Optimization, 38(2):447–469.
  * Bottou et al., [2018] Bottou, L., Curtis, F. E., and Nocedal, J. (2018).  Optimization methods for large-scale machine learning.  Siam Review, 60(2):223–311.
  * Bowling and Veloso, [2001] Bowling, M. and Veloso, M. (2001).  Rational and convergent learning in stochastic games.  In Proceedings of the 17th International Joint Conference on Artificial Intelligence - Volume 2, IJCAI’01, page 1021–1026, San Francisco, CA, USA. Morgan Kaufmann Publishers Inc.
  * Brown, [1951] Brown, G. W. (1951).  Iterative solution of games by fictitious play.  Activity Analysis of Production and Allocation, 13(1):374–376.
  * Busoniu et al., [2008] Busoniu, L., Babuska, R., De Schutter, B., et al. (2008).  A comprehensive survey of multiagent reinforcement learning.  IEEE Transactions on Systems, Man, and Cybernetics, Part C, 38(2):156–172.
  * [22] Cai, Q., Yang, Z., Lee, J. D., and Wang, Z. (2023a).  Neural temporal difference and $Q$-learning provably converge to global optima.  Mathematics of Operations Research.
  * [23] Cai, Y., Luo, H., Wei, C.-Y., and Zheng, W. (2023b).  Uncoupled and convergent learning in two-player zero-sum Markov games.  arXiv preprint arXiv:2303.02738.
  * Carvalho et al., [2020] Carvalho, D., Melo, F. S., and Santos, P. (2020).  A new convergent variant of $Q$-learning with linear function approximation.  Advances in Neural Information Processing Systems, 33.
  * Cesa-Bianchi and Lugosi, [2006] Cesa-Bianchi, N. and Lugosi, G. (2006).  Prediction, Learning, and Games.  Cambridge University Press.
  * [26] Chen, F., Mei, S., and Bai, Y. (2022a).  Unified algorithms for RL with decision-estimation coefficients: No-regret, PAC, and reward-free learning.  arXiv preprint arXiv:2209.11745.
  * [27] Chen, S., Devraj, A. M., Lu, F., Busic, A., and Meyn, S. (2020a).  Zap $Q$-learning with nonlinear function approximation.  Advances in Neural Information Processing Systems, 33:16879–16890.
  * [28] Chen, Z., Clarke, J. P., and Maguluri, S. T. (2022b).  Target Network and Truncation Overcome The Deadly triad in $Q$-Learning.  Preprint arXiv:2203.02628.
  * [29] Chen, Z., Ma, S., and Zhou, Y. (2021a).  Sample efficient stochastic policy extragradient algorithm for zero-sum markov game.  In International Conference on Learning Representations.
  * [30] Chen, Z., Maguluri, S. T., Shakkottai, S., and Shanmugam, K. (2020b).  Finite-Sample Analysis of Contractive Stochastic Approximation Using Smooth Convex Envelopes.  Advances in Neural Information Processing Systems, 33.
  * [31] Chen, Z., Maguluri, S. T., Shakkottai, S., and Shanmugam, K. (2021b).  A Lyapunov theory for finite-sample guarantees of asynchronous $Q$-learning and TD-learning variants.  Preprint arXiv:2102.01567.
  * Chen et al., [2023] Chen, Z., Zhang, K., Mazumdar, E., Ozdaglar, A., and Wierman, A. (2023).  A Finite-Sample Analysis of Payoff-Based Independent Learning in Zero-Sum Stochastic Games.  Preprint arXiv:2303.03100.
  * Chen et al., [2019] Chen, Z., Zhang, S., Doan, T. T., Clarke, J.-P., and Maguluri, S. T. (2019).  Finite-Sample Analysis of Nonlinear Stochastic Approximation with Applications in Reinforcement Learning.  Preprint arXiv:1905.11425.
  * [34] Chen, Z., Zhou, Y., Chen, R.-R., and Zou, S. (2022c).  Sample and communication-efficient decentralized actor-critic algorithms with finite-time analysis.  In International Conference on Machine Learning, pages 3794–3834. PMLR.
  * Cisneros-Velarde and Koyejo, [2023] Cisneros-Velarde, P. and Koyejo, S. (2023).  Finite-sample guarantees for nash Q-learning with linear function approximation.  arXiv preprint arXiv:2303.00177.
  * Cui et al., [2023] Cui, Q., Zhang, K., and Du, S. S. (2023).  Breaking the curse of multiagents in a large state space: RL in Markov games with independent linear function approximation.  arXiv preprint arXiv:2302.03673.
  * Daskalakis et al., [2020] Daskalakis, C., Foster, D. J., and Golowich, N. (2020).  Independent policy gradient methods for competitive reinforcement learning.  Advances in neural information processing systems, 33:5527–5540.
  * Daskalakis et al., [2022] Daskalakis, C., Golowich, N., and Zhang, K. (2022).  The complexity of Markov equilibrium in stochastic games.  Preprint arXiv:2204.03991.
  * Ding et al., [2022] Ding, D., Wei, C.-Y., Zhang, K., and Jovanovic, M. (2022).  Independent policy gradient for large-scale markov potential games: Sharper rates, function approximation, and game-agnostic convergence.  In International Conference on Machine Learning, pages 5166–5220. PMLR.
  * Du et al., [2019] Du, S. S., Luo, Y., Wang, R., and Zhang, H. (2019).  Provably efficient $Q$-learning with function approximation via distribution shift error checking oracle.  In Proceedings of the 33rd International Conference on Neural Information Processing Systems, pages 8060–8070.
  * Erez et al., [2022] Erez, L., Lancewicki, T., Sherman, U., Koren, T., and Mansour, Y. (2022).  Regret minimization and convergence to equilibria in general-sum Markov games.  Preprint arXiv:2207.14211.
  * Even-Dar and Mansour, [2003] Even-Dar, E. and Mansour, Y. (2003).  Learning rates for $Q$-learning.  Journal of Machine Learning Research, 5(Dec):1–25.
  * Fink, [1964] Fink, A. M. (1964).  Equilibrium in a stochastic $n$-person game.  Journal of science of the hiroshima university, series ai (mathematics), 28(1):89–93.
  * Foster and Young, [2006] Foster, D. and Young, H. P. (2006).  Regret testing: Learning to play Nash equilibrium without knowing you have an opponent.  Theoretical Economics, 1(3):341–367.
  * [45] Foster, D. J., Foster, D. P., Golowich, N., and Rakhlin, A. (2023a).  On the complexity of multi-agent decision making: From learning in games to partial monitoring.  arXiv preprint arXiv:2305.00684.
  * [46] Foster, D. J., Golowich, N., and Kakade, S. M. (2023b).  Hardness of independent learning and sparse equilibrium computation in Markov games.  arXiv preprint arXiv:2303.12287.
  * Fudenberg and Kreps, [1993] Fudenberg, D. and Kreps, D. (1993).  Learning mixed equilibria.  Games and Economic Behavior, 5:320–367.
  * Fudenberg and Levine, [1998] Fudenberg, D. and Levine, D. K. (1998).  The Theory of Learning in Games, volume 2.  MIT press.
  * Gao et al., [2021] Gao, Z., Ma, Q., Başar, T., and Birge, J. R. (2021).  Finite-Sample Analysis of Decentralized $Q$-Learning for Stochastic Games.  Preprint arXiv:2112.07859.
  * Guzmán and Nemirovski, [2015] Guzmán, C. and Nemirovski, A. (2015).  On lower complexity bounds for large-scale smooth convex optimization.  Journal of Complexity, 31(1):1–14.
  * Harris, [1998] Harris, C. (1998).  On the rate of convergence of continuous-time fictitious play.  Games and Economic Behavior, 22(2):238–259.
  * Hofbauer and Hopkins, [2005] Hofbauer, J. and Hopkins, E. (2005).  Learning in perturbed asymmetric games.  Games and Economic Behavior, 52(1):133–152.
  * Huang et al., [2022] Huang, B., Lee, J. D., Wang, Z., and Yang, Z. (2022).  Towards general function approximation in zero-sum Markov games.  In International Conference on Learning Representations.
  * Jaakkola et al., [1994] Jaakkola, T., Jordan, M. I., and Singh, S. P. (1994).  Convergence of stochastic iterative dynamic programming algorithms.  In Advances in neural information processing systems, pages 703–710.
  * Jin et al., [2018] Jin, C., Allen-Zhu, Z., Bubeck, S., and Jordan, M. I. (2018).  Is $Q$-learning provably efficient?  In Proceedings of the 32nd International Conference on Neural Information Processing Systems, pages 4868–4878.
  * Jin et al., [2021] Jin, C., Liu, Q., Wang, Y., and Yu, T. (2021).  V-learning – A simple, efficient, decentralized algorithm for multiagent RL.  Preprint arXiv:2110.14555.
  * Jin et al., [2022] Jin, C., Liu, Q., and Yu, T. (2022).  The power of exploiter: Provable multi-agent RL in large state spaces.  In International Conference on Machine Learning, pages 10251–10279. PMLR.
  * Jin et al., [2020] Jin, C., Yang, Z., Wang, Z., and Jordan, M. I. (2020).  Provably efficient reinforcement learning with linear function approximation.  In Conference on Learning Theory, pages 2137–2143. PMLR.
  * Kalogiannis et al., [2022] Kalogiannis, F., Anagnostides, I., Panageas, I., Vlatakis-Gkaragkounis, E.-V., Chatziafratis, V., and Stavroulakis, S. A. (2022).  Efficiently computing Nash equilibria in adversarial team Markov games.  In The Eleventh International Conference on Learning Representations.
  * Khodadadian et al., [2022] Khodadadian, S., Doan, T. T., Romberg, J., and Maguluri, S. T. (2022).  Finite sample analysis of two-time-scale natural actor-critic algorithm.  IEEE Transactions on Automatic Control.
  * Kober et al., [2013] Kober, J., Bagnell, J. A., and Peters, J. (2013).  Reinforcement learning in robotics: A survey.  The International Journal of Robotics Research, 32(11):1238–1274.
  * Konda and Tsitsiklis, [2004] Konda, V. R. and Tsitsiklis, J. N. (2004).  Convergence rate of linear two-time-scale stochastic approximation.  The Annals of Applied Probability, 14(2):796 – 819.
  * Kushner and Clark, [2012] Kushner, H. J. and Clark, D. S. (2012).  Stochastic approximation methods for constrained and unconstrained systems, volume 26.  Springer Science & Business Media.
  * Lakshminarayanan and Szepesvari, [2018] Lakshminarayanan, C. and Szepesvari, C. (2018).  Linear Stochastic Approximation: How Far Does Constant Step-Size and Iterate Averaging Go?  In International Conference on Artificial Intelligence and Statistics, pages 1347–1355.
  * Lan, [2020] Lan, G. (2020).  First-order and Stochastic Optimization Methods for Machine Learning.  Springer.
  * Lee and He, [2020] Lee, D. and He, N. (2020).  A unified switching system perspective and convergence analysis of $Q$-learning algorithms.  Advances in Neural Information Processing Systems, 33:15556–15567.
  * Leslie and Collins, [2003] Leslie, D. S. and Collins, E. (2003).  Convergent multiple-timescales reinforcement learning algorithms in normal form games.  The Annals of Applied Probability, 13(4):1231–1251.
  * Leslie and Collins, [2005] Leslie, D. S. and Collins, E. J. (2005).  Individual $Q$-learning in normal form games.  SIAM Journal on Control and Optimization, 44(2):495–514.
  * [69] Li, G., Cai, C., Chen, Y., Gu, Y., Wei, Y., and Chi, Y. (2021a).  Tightening the dependence on horizon in the sample complexity of $Q$-learning.  In International Conference on Machine Learning, pages 6296–6306. PMLR.
  * [70] Li, G., Shi, L., Chen, Y., Gu, Y., and Chi, Y. (2021b).  Breaking the sample complexity barrier to regret-optimal model-free reinforcement learning.  Advances in Neural Information Processing Systems, 34.
  * Li et al., [2020] Li, G., Wei, Y., Chi, Y., Gu, Y., and Chen, Y. (2020).  Sample Complexity of Asynchronous $Q$-Learning: Sharper Analysis and Variance Reduction.  In Advances in Neural Information Processing Systems, volume 33, pages 7031–7043. Curran Associates, Inc.
  * Littman, [1994] Littman, M. L. (1994).  Markov games as a framework for multi-agent reinforcement learning.  In Proceedings of the Eleventh International Conference on International Conference on Machine Learning, pages 157–163.
  * Liu et al., [2022] Liu, Q., Wang, Y., and Jin, C. (2022).  Learning Markov games with adversarial opponents: Efficient algorithms and fundamental limits.  In International Conference on Machine Learning, pages 14036–14053. PMLR.
  * Maei et al., [2010] Maei, H. R., Szepesvári, C., Bhatnagar, S., and Sutton, R. S. (2010).  Toward off-policy learning control with function approximation.  In Proceedings of the 27th International Conference on Machine Learning (ICML-10), pages 719–726.
  * Mao et al., [2022] Mao, W., Yang, L., Zhang, K., and Başar, T. (2022).  On improving model-free algorithms for decentralized multi-agent reinforcement learning.  In International Conference on Machine Learning, pages 15007–15049. PMLR.
  * Marden and Shamma, [2012] Marden, J. R. and Shamma, J. S. (2012).  Revisiting log-linear learning: Asynchrony, completeness and payoff-based implementation.  Games and Economic Behavior, 75(2):788–808.
  * Melo et al., [2008] Melo, F. S., Meyn, S. P., and Ribeiro, M. I. (2008).  An analysis of reinforcement learning with function approximation.  In Proceedings of the 25th international conference on Machine learning, pages 664–671.
  * Meyn, [2023] Meyn, S. (2023).  Stability of $Q$-learning through design and optimism.  arXiv preprint arXiv:2307.02632.
  * Mitrophanov, [2005] Mitrophanov, A. Y. (2005).  Sensitivity and convergence of uniformly ergodic markov chains.  Journal of Applied Probability, 42(4):1003–1014.
  * Mnih et al., [2013] Mnih, V., Kavukcuoglu, K., Silver, D., Graves, A., Antonoglou, I., Wierstra, D., and Riedmiller, M. (2013).  Playing atari with deep reinforcement learning.  Preprint arXiv:1312.5602.
  * Mnih et al., [2015] Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A. A., Veness, J., Bellemare, M. G., Graves, A., Riedmiller, M., Fidjeland, A. K., Ostrovski, G., et al. (2015).  Human-level control through deep reinforcement learning.  nature, 518(7540):529–533.
  * Mou et al., [2020] Mou, W., Li, C. J., Wainwright, M. J., Bartlett, P. L., and Jordan, M. I. (2020).  On linear stochastic approximation: Fine-grained Polyak-Ruppert and non-asymptotic concentration.  In Conference on Learning Theory, pages 2947–2997. PMLR.
  * Munos and Szepesvári, [2008] Munos, R. and Szepesvári, C. (2008).  Finite-Time Bounds for Fitted Value Iteration.  Journal of Machine Learning Research, 9(5).
  * Nash Jr, [1950] Nash Jr, J. F. (1950).  Equilibrium points in $n$-person games.  Proceedings of the national academy of sciences, 36(1):48–49.
  * Ni et al., [2022] Ni, C., Song, Y., Zhang, X., Jin, C., and Wang, M. (2022).  Representation learning for general-sum low-rank Markov games.  arXiv preprint arXiv:2210.16976.
  * Ozdaglar et al., [2021] Ozdaglar, A., Sayin, M. O., and Zhang, K. (2021).  Independent learning in stochastic games.  arXiv preprint arXiv:2111.11743.
  * Powell, [2007] Powell, W. B. (2007).  Approximate Dynamic Programming: Solving the curses of dimensionality, volume 703.  John Wiley & Sons.
  * Qiu et al., [2021] Qiu, S., Yang, Z., Ye, J., and Wang, Z. (2021).  On finite-time convergence of actor-critic algorithm.  IEEE Journal on Selected Areas in Information Theory, 2(2):652–664.
  * Qu and Wierman, [2020] Qu, G. and Wierman, A. (2020).  Finite-Time Analysis of Asynchronous Stochastic Approximation and $Q$-Learning.  In Conference on Learning Theory, pages 3185–3205. PMLR.
  * Robbins and Monro, [1951] Robbins, H. and Monro, S. (1951).  A stochastic approximation method.  The Annals of Mathematical Statistics, pages 400–407.
  * Robinson, [1951] Robinson, J. (1951).  An iterative method of solving a game.  Annals of Mathematics, pages 296–301.
  * Rummery and Niranjan, [1994] Rummery, G. A. and Niranjan, M. (1994).  On-line $Q$-learning using connectionist systems, volume 37.  University of Cambridge, Department of Engineering Cambridge, UK.
  * Sayin et al., [2021] Sayin, M., Zhang, K., Leslie, D., Basar, T., and Ozdaglar, A. (2021).  Decentralized $Q$-learning in zero-sum Markov games.  Advances in Neural Information Processing Systems, 34:18320–18334.
  * [94] Sayin, M. O., Parise, F., and Ozdaglar, A. (2022a).  Fictitious play in zero-sum stochastic games.  SIAM Journal on Control and Optimization, 60(4):2095–2114.
  * [95] Sayin, M. O., Zhang, K., and Ozdaglar, A. (2022b).  Fictitious Play in Markov Games with Single Controller.  In Proceedings of the 23rd ACM Conference on Economics and Computation, pages 919–936.
  * Shalev-Shwartz et al., [2016] Shalev-Shwartz, S., Shammah, S., and Shashua, A. (2016).  Safe, multi-agent, reinforcement learning for autonomous driving.  Preprint arXiv:1610.03295.
  * Shapley, [1953] Shapley, L. S. (1953).  Stochastic games.  Proceedings of the National Academy of Sciences, 39(10):1095–1100.
  * Silver et al., [2016] Silver, D., Huang, A., Maddison, C. J., Guez, A., Sifre, L., Van Den Driessche, G., Schrittwieser, J., Antonoglou, I., Panneershelvam, V., Lanctot, M., et al. (2016).  Mastering the game of go with deep neural networks and tree search.  nature, 529(7587):484.
  * Song et al., [2023] Song, Z., Lee, J. D., and Yang, Z. (2023).  Can we find Nash equilibria at a linear rate in Markov games?  arXiv preprint arXiv:2303.03095.
  * Song et al., [2022] Song, Z., Mei, S., and Bai, Y. (2022).  When can we learn general-sum markov games with a large number of players sample-efficiently?  In International Conference on Learning Representations.
  * Srikant and Ying, [2019] Srikant, R. and Ying, L. (2019).  Finite-Time Error Bounds For Linear Stochastic Approximation and TD Learning.  In Conference on Learning Theory, pages 2803–2830.
  * Sutton, [1988] Sutton, R. S. (1988).  Learning to predict by the methods of temporal differences.  Machine learning, 3(1):9–44.
  * Sutton, [2015] Sutton, R. S. (2015).  Introduction to reinforcement learning with function approximation.  In Tutorial at the Conference on Neural Information Processing Systems, page 33.
  * Sutton and Barto, [2018] Sutton, R. S. and Barto, A. G. (2018).  Reinforcement learning: An introduction.  MIT press.
  * Szepesvári and Munos, [2005] Szepesvári, C. and Munos, R. (2005).  Finite time bounds for sampling based fitted value iteration.  In Proceedings of the 22nd international conference on Machine learning, pages 880–887.
  * Tsitsiklis, [1994] Tsitsiklis, J. N. (1994).  Asynchronous stochastic approximation and $Q$-learning.  Machine learning, 16(3):185–202.
  * Tsitsiklis and Van Roy, [1997] Tsitsiklis, J. N. and Van Roy, B. (1997).  An analysis of temporal-difference learning with function approximation.  IEEE Transactions on Automatic Control, 42(5):674–690.
  * Wang et al., [2023] Wang, Y., Liu, Q., Bai, Y., and Jin, C. (2023).  Breaking the curse of multiagency: Provably efficient decentralized multi-agent RL with function approximation.  arXiv preprint arXiv:2302.06606.
  * Watkins and Dayan, [1992] Watkins, C. J. and Dayan, P. (1992).  $Q$-learning.  Machine learning, 8(3-4):279–292.
  * Wei et al., [2021] Wei, C.-Y., Lee, C.-W., Zhang, M., and Luo, H. (2021).  Last-iterate convergence of decentralized optimistic gradient descent/ascent in infinite-horizon competitive Markov games.  In Conference on Learning Theory, pages 4259–4299. PMLR.
  * Weisz et al., [2021] Weisz, G., Amortila, P., and Szepesvári, C. (2021).  Exponential lower bounds for planning in MDPs with linearly-realizable optimal action-value functions.  In Algorithmic Learning Theory, pages 1237–1264. PMLR.
  * Wu et al., [2020] Wu, Y. F., Zhang, W., Xu, P., and Gu, Q. (2020).  A finite-time analysis of two time-scale actor-critic methods.  Advances in Neural Information Processing Systems, 33:17617–17628.
  * Xie et al., [2020] Xie, Q., Chen, Y., Wang, Z., and Yang, Z. (2020).  Learning zero-sum simultaneous-move Markov games using function approximation and correlated equilibrium.  In Conference on Learning Theory, pages 3674–3682. PMLR.
  * Xu and Gu, [2020] Xu, P. and Gu, Q. (2020).  A finite-time analysis of $Q$-learning with neural network function approximation.  In International Conference on Machine Learning, pages 10555–10565. PMLR.
  * Xu and Liang, [2021] Xu, T. and Liang, Y. (2021).  Sample complexity bounds for two timescale value-based reinforcement learning algorithms.  In International Conference on Artificial Intelligence and Statistics, pages 811–819. PMLR.
  * Yan et al., [2022] Yan, Y., Li, G., Chen, Y., and Fan, J. (2022).  The efficacy of pessimism in asynchronous Q-learning.  Preprint arXiv:2203.07368.
  * Yang and Wang, [2020] Yang, L. and Wang, M. (2020).  Reinforcement learning in feature space: Matrix bandit, kernels, and regret bound.  In International Conference on Machine Learning, pages 10746–10756. PMLR.
  * Zanette et al., [2020] Zanette, A., Lazaric, A., Kochenderfer, M., and Brunskill, E. (2020).  Learning near optimal policies with low inherent Bellman error.  In International Conference on Machine Learning, pages 10978–10989. PMLR.
  * Zanette and Wainwright, [2022] Zanette, A. and Wainwright, M. (2022).  Stabilizing $Q$-learning with linear architectures for provable efficient learning.  In International Conference on Machine Learning, pages 25920–25954. PMLR.
  * Zhang et al., [2020] Zhang, K., Kakade, S., Başar, T., and Yang, L. (2020).  Model-based multi-agent RL in zero-sum Markov games with near-optimal sample complexity.  Advances in Neural Information Processing Systems, 33:1166–1178.
  * [121] Zhang, K., Yang, Z., and Başar, T. (2021a).  Multi-agent reinforcement learning: A selective overview of theories and algorithms.  Handbook of Reinforcement Learning and Control, pages 321–384.
  * Zhang et al., [2022] Zhang, R., Liu, Q., Wang, H., Xiong, C., Li, N., and Bai, Y. (2022).  Policy optimization for Markov games: Unified framework and faster convergence.  In Advances in Neural Information Processing Systems.
  * [123] Zhang, S., Yao, H., and Whiteson, S. (2021b).  Breaking the Deadly Triad with a Target Network.  In Proceedings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pages 12621–12631. PMLR.
  * Zhao et al., [2021] Zhao, Y., Tian, Y., Lee, J. D., and Du, S. S. (2021).  Provably Efficient Policy Optimization for Two-Player Zero-Sum Markov Games.  Preprint arXiv:2102.08903.
  * Zhou et al., [2021] Zhou, D., Gu, Q., and Szepesvari, C. (2021).  Nearly minimax optimal reinforcement learning for linear mixture Markov decision processes.  In Conference on Learning Theory, pages 4532–4576. PMLR.
  * Zou et al., [2019] Zou, S., Xu, T., and Liang, Y. (2019).  Finite-sample analysis for SARSA with linear function approximation.  In Advances in Neural Information Processing Systems, pages 8668–8678.

Appendices

##  Appendix A Proof of Theorem 4.1

We begin by rewriting Algorithm 3 in the following when using linear function
approximation, where we denote $\tilde{q}_{t,k}^{i}=\Phi^{i}\theta_{t,k}^{i}$
and $q_{t,k}^{i}=\Phi^{i}w_{t,k}^{i}$ for all $t,k$. Other notation that will
be used in the proof is summarized in Appendix A.1. Our proof is divided into
$4$ steps. In Appendix A.2, we bound the Nash gap by the Lyapunov functions.
In Appendix A.3, we analyze the outer loop, which is an approximation of the
minimax value iteration. In Appendix A.4, we analyze the inner loop, where the
slow-timescale iterates $\\{\theta_{t,k}^{i}\\}$ are studied in Appendix A.4.1
using the Lyapunov function constructed in Section 5, and the fast-timescale
iterates $\\{w_{t,k}^{i}\\}$ are studied in Appendix A.4.3. Finally, in
Appendix A.5, we establish the finite-sample bounds. The proof of all
supporting lemmas is provided in Appendix A.6.

Algorithm 3 Algorithm 1 Under Linear Function Approximation

1: Input: Integers $K$ and $T$, initializations
$\tilde{q}_{0,0}^{i}=q_{0,0}^{i}=0$, $v_{0}^{i}=0$, and
$\pi_{0,0}^{i}(s)=\text{Unif}(\mathcal{A}^{i})$ for all $s\in\mathcal{S}$

2: for $t=0,1,\cdots,T-1$ do

3: for $k=0,1,\cdots,K-1$ do

4:
$\tilde{q}^{i}_{t,k+1}=\tilde{q}_{t,k}^{i}+\beta_{k}(q_{t,k}^{i}-\tilde{q}_{t,k}^{i})$
and compute
$\pi_{t,k+1}^{i}(S_{k})=\sigma_{\tau}(\tilde{q}_{t,k+1}^{i}(S_{k}))$

5: Play $A_{k}^{i}\sim\pi_{k+1}^{i}(\cdot\mid S_{k})$ and observe $S_{k+1}\sim
p(\cdot\mid S_{k},A_{k}^{i},A_{k}^{-i})$

6:
$w_{t,k+1}^{i}=\text{Proj}_{M}^{i}\left[w_{t,k}^{i}+\alpha_{k}\phi^{i}(S_{k},A_{k}^{i})\left(R_{i}(S_{k},A_{k}^{i},A_{k}^{-i})+\gamma
v_{t}^{i}(S_{k+1})-\phi^{i}(S_{k},A_{k}^{i})^{\top}w_{t,k}^{i}\right)\right]$

7: end for

8: $v_{t+1}^{i}(s)=\pi_{t,K}^{i}(s)^{\top}\lceil q_{t,K}^{i}(s)\rceil$ for all
$s\in\mathcal{S}$

9: $S_{0}=S_{K}$, $q_{t+1,0}^{i}=q_{t,K}^{i}$, and
$\tilde{q}_{t+1,0}^{i}=\tilde{q}_{t,K}^{i}$

10: end for

### A.1 Notation

We begin by summarizing the notation we are going to use here.

  1. (1)

Given non-negative integers $k_{1}\leq k_{2}$, we denote
$\alpha_{k_{1},k_{2}}=\sum_{k=k_{1}}^{k_{2}}\alpha_{k}$ and
$\beta_{k_{1},k_{2}}=\sum_{k=k_{1}}^{k_{2}}\beta_{k}$.

  2. (2)

Given a pair of matrices
$\\{X^{i}\in\mathbb{R}^{|\mathcal{A}^{i}|\times|\mathcal{A}^{-i}|}\\}_{i\in\\{1,2\\}}$,
we define a function
$V_{X}:\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}\mapsto\mathbb{{R}}$
as

$\displaystyle
V_{X}(x^{1},x^{2})=\sum_{i=1,2}\max_{u^{i}\in\Delta(\mathcal{A}^{i})}\min_{\bar{u}^{i}\in\Delta(\mathcal{A}^{i})}\left\\{(u^{i}-\bar{u}^{i})^{\top}x^{i}+\tau\nu(u^{i})-\tau\nu(\bar{u}^{i})+\frac{1}{2\mu}\|x^{-i}-X_{-i}\bar{u}^{i}\|_{2}^{2}\right\\}$

for all $x^{1}\in\mathbb{R}^{|\mathcal{A}^{1}|}$ and
$x^{2}\in\mathbb{R}^{|\mathcal{A}^{2}|}$, where $\nu(\cdot)$ is the entropy
function defined as
$\nu(x^{i})=-\sum_{a^{i}\in\mathcal{A}^{i}}x^{i}(a^{i})\log(x^{i}(a^{i}))$ for
$i\in\\{1,2\\}$. The function $V_{X}(\cdot,\cdot)$ serves as our Lyapunov
function to analyze the evolution of the policy parameter
$\\{\theta_{t,k}^{i}\\}$. A sequence of properties related to
$V_{X}(\cdot,\cdot)$ is established in Appendix A.4.2.

  3. (3)

Given a pair of $v$-functions $(v^{1},v^{2})$ and a state $s\in\mathcal{S}$,
we define
$V_{v,s}:\mathbb{R}^{|\mathcal{A}^{1}|}\times\mathbb{R}^{|\mathcal{A}^{2}|}\mapsto\mathbb{{R}}$
as $V_{v,s}(x^{1},x^{2})=V_{X}(x^{1},x^{2})$ with
$X_{i}=\mathcal{T}^{i}(v^{i})(s)$ for $i\in\\{1,2\\}$.

  4. (4)

For any joint policy $(\pi^{1},\pi^{2})$ and state $s\in\mathcal{S}$, given
$i\in\\{1,2\\}$, let
$v^{i}_{*,\pi^{-i}}(s)=\max_{\hat{\pi}^{i}}v^{i}_{\hat{\pi}^{i},\pi^{-i}}(s)$,
$v^{i}_{\pi^{i},*}=\min_{\hat{\pi}^{-i}}v^{i}_{\pi^{i},\hat{\pi}^{-i}}(s)$,
$v^{-i}_{\pi^{-i},*}(s)=\min_{\hat{\pi}^{i}}v^{-i}_{\pi^{-i},\hat{\pi}^{i}}(s)$,
and
$v^{-i}_{*,\pi^{i}}(s)=\max_{\hat{\pi}^{-i}}v^{-i}_{\hat{\pi}^{-i},\hat{\pi}^{i}}(s)$.
Note that we have $v^{i}_{*,\pi^{-i}}+v^{-i}_{\pi^{-i},*}=0$ and
$v^{i}_{\pi^{i},*}+v^{-i}_{*,\pi^{i}}=0$ for $i\in\\{1,2\\}$.

  5. (5)

For any $t,k$ and $i\in\\{1,2\\}$, we define
$q_{t,k}^{i},\tilde{q}_{t,k}^{i}$, and
$\bar{q}_{t,k}^{i}\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{i}|}$ as

$\displaystyle
q_{t,k}^{i}=\Phi^{i}w_{t,k}^{i},\quad\tilde{q}_{t,k}^{i}=\Phi^{i}\theta_{t,k}^{i},\quad\bar{q}_{t,k}^{i}(s)=\mathcal{T}^{i}(v^{i})(s)\pi_{t,k}^{-i}(s),\quad\forall\,s\in\mathcal{S}.$

Moreover, we let $\bar{w}_{t,k}^{i}\in\mathbb{R}^{d_{i}}$ be such that
$\bar{q}_{t,k}^{i}=\Phi^{i}\bar{w}_{t,k}^{i}$. Note that $\bar{w}_{t,k}^{i}$
exists and is unique under Assumption 4.3.

  6. (6)

For any $t,k\geq 0$, define
$\mathcal{L}_{v}(t)=\sum_{i=1,2}\|v_{t}^{i}-v_{*}^{i}\|_{\infty}$,
$\mathcal{L}_{\text{sum}}(t)=\|v_{t}^{1}+v_{t}^{2}\|_{\infty}$,
$\mathcal{L}_{\theta}(t,k)=\max_{s\in\mathcal{S}}V_{v_{t},s}(\tilde{q}_{t,k}^{1}(s),\tilde{q}_{t,k}^{2}(s))$,
and
$\mathcal{L}_{w}(t,k)=\sum_{i=1,2}\|w_{t,k}^{i}-\bar{w}_{t,k}^{i}\|_{2}^{2}$.

### A.2 Bounding the Nash Gap

######  Lemma A.1.

It holds for any joint policy $\pi=(\pi^{1},\pi^{2})$ that

$\displaystyle\text{NG}\left(\pi^{1},\pi^{2}\right)\leq
2\sum_{i=1,2}\left\|v^{i}_{\pi^{i},*}-v^{i}_{*}\right\|_{\infty}.$

###### Proof of Lemma A.1.

By definition of the expected value function, given $i\in\\{1,2\\}$, we have

$\displaystyle\max_{\hat{\pi}^{i}}U^{i}(\hat{\pi}^{i},\pi^{-i})=\;$
$\displaystyle\max_{\hat{\pi}^{i}}\mathbb{E}_{S\sim\rho_{o}(\cdot)}\left[v_{\hat{\pi}^{i},\pi^{-i}}^{i}(S)\right]$
$\displaystyle\leq\;$
$\displaystyle\mathbb{E}_{S\sim\rho_{o}(\cdot)}\left[\max_{\hat{\pi}^{i}}v_{\hat{\pi}^{i},\pi^{-i}}^{i}(S)\right]$
(Jensen’s inequality) $\displaystyle=\;$
$\displaystyle\mathbb{E}_{S\sim\rho_{o}(\cdot)}\left[v_{*,\pi^{-i}}^{i}(S)\right],$

where the last line follows from $v_{*,\pi^{-i}}^{i}(s)\geq
v_{\hat{\pi}^{i},\pi^{-i}}^{i}(s)$ for any policy $\hat{\pi}^{i}$ and
$s\in\mathcal{S}$. Therefore, we have

$\displaystyle\text{NG}\left(\pi^{1},\pi^{2}\right)=\;$
$\displaystyle\sum_{i=1,2}\left[\max_{\hat{\pi}^{i}}U^{i}(\hat{\pi}^{i},\pi^{-i})-U^{i}(\pi^{i},\pi^{-i})\right]$
$\displaystyle\leq\;$ $\displaystyle\sum_{i=1,2}\mathbb{E}_{S\sim
p_{o}}\left[v^{i}_{*,\pi^{-i}}(S)-v^{i}_{\pi^{i},\pi^{-i}}(S)\right]$
$\displaystyle\leq\;$
$\displaystyle\sum_{i=1,2}\left\|v^{i}_{*,\pi^{-i}}-v^{i}_{\pi^{i},\pi^{-i}}\right\|_{\infty}.$
(21)

To proceed, observe that for any $i\in\\{1,2\\}$ and $s\in\mathcal{S}$, we
have

$\displaystyle\left|v^{i}_{*,\pi^{-i}}(s)-v^{i}_{\pi^{i},\pi^{-i}}(s)\right|=\;$
$\displaystyle v^{i}_{*,\pi^{-i}}(s)-v^{i}_{\pi^{i},\pi^{-i}}(s)$
($v^{i}_{*,\pi^{-i}}(s)\geq v^{i}_{\pi^{i},\pi^{-i}}(s)$ for all $s$)
$\displaystyle\leq\;$ $\displaystyle
v^{i}_{*,\pi^{-i}}(s)-v^{i}_{\pi^{i},*}(s)$ ($v^{i}_{\pi^{i},\pi^{-i}}(s)\geq
v^{i}_{\pi^{i},*}(s)$ for all $s$) $\displaystyle=\;$
$\displaystyle-v^{-i}_{\pi^{-i},*}(s)-v^{i}_{\pi^{i},*}(s)$
($v^{i}_{*,\pi^{-i}}+v^{-i}_{\pi^{-i},*}=0$) $\displaystyle=\;$ $\displaystyle
v_{*}^{-i}(s)-v^{-i}_{\pi^{-i},*}(s)+v_{*}^{i}(s)-v^{i}_{\pi^{i},*}(s)$
($v_{*}^{i}+v_{*}^{-i}=0$) $\displaystyle\leq\;$
$\displaystyle\sum_{i=1,2}\left\|v^{i}_{\pi^{i},*}-v^{i}_{*}\right\|_{\infty}.$

Since the RHS of the previous inequality does not depend on $s$, we have for
$i\in\\{1,2\\}$ that

$\displaystyle\left\|v^{i}_{*,\pi^{-i}}-v^{i}_{\pi^{i},\pi^{-i}}\right\|_{\infty}\leq\sum_{i=1,2}\left\|v^{i}_{\pi^{i},*}-v^{i}_{*}\right\|_{\infty}.$

The result follows from using the previous inequality in Eq. (A.2). ∎

######  Lemma A.2.

For any $\tilde{q}^{1}\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{1}|}$ and
$\tilde{q}^{2}\in\mathbb{R}^{|\mathcal{S}||\mathcal{A}^{2}|}$, let
$\pi^{i}(s)=\sigma_{\tau}(\tilde{q}^{i}(s))$ for all $s\in\mathcal{S}$ and
$i\in\\{1,2\\}$. Then, we have for any
$v^{1},v^{2}\in\mathbb{R}^{|\mathcal{S}|}$ that

$\displaystyle\sum_{i=1,2}\|v^{i}_{\pi^{i},*}-v^{i}_{*}\|_{\infty}\leq\;$
$\displaystyle\frac{2}{1-\gamma}\left(\sum_{i=1,2}\|v^{i}-v^{i}_{*}\|_{\infty}+2\|v^{1}+v^{2}\|_{\infty}\right.$
$\displaystyle+\left.\frac{9A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\max_{s\in\mathcal{S}}V_{v,s}(\tilde{q}^{1}(s),\tilde{q}^{2}(s))+6\tau\log(A_{\max})+\mu\right).$

###### Proof of Lemma A.2.

Since $v^{-i}_{*}-v^{-i}_{\pi^{-i},*}=v^{i}_{*,\pi^{-i}}-v^{i}_{*}$ for
$i\in\\{1,2\\}$, it is enough to bound
$\|v^{i}_{*,\pi^{-i}}-v^{i}_{*}\|_{\infty}$. Given $i\in\\{1,2\\}$, for any
$s\in\mathcal{S}$, we have

$\displaystyle 0\leq\;$ $\displaystyle v^{i}_{*,\pi^{-i}}(s)-v^{i}_{*}(s)$
$\displaystyle=\;$
$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i}_{*,\pi^{-i}})(s)\pi^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{*}^{i})(s)\mu^{-i}$
(Bellman optimality equations) $\displaystyle\leq\;$
$\displaystyle\left|\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i}_{*,\pi^{-i}})(s)\pi^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i}_{*})(s)\pi^{-i}(s)\right|$
$\displaystyle+\left|\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i}_{*})(s)\pi^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)\right|$
$\displaystyle+\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}$
$\displaystyle+\left|\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{*}^{i})(s)\mu^{-i}\right|$
(Triangle inequality) $\displaystyle\leq\;$
$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\left|(\mu^{i})^{\top}(\mathcal{T}^{i}(v^{i}_{*,\pi^{-i}})(s)-\mathcal{T}^{i}(v^{i}_{*})(s))\pi^{-i}(s)\right|$
$\displaystyle+\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\left|(\mu^{i})^{\top}(\mathcal{T}^{i}(v^{i}_{*})(s)-\mathcal{T}^{i}(v^{i})(s))\pi^{-i}(s)\right|$
$\displaystyle+\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}$
$\displaystyle+\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\max_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}\left|(\mu^{i})^{\top}(\mathcal{T}^{i}(v^{i})(s)-\mathcal{T}^{i}(v_{*}^{i})(s))\mu^{-i}\right|.$

Using Lemma A.19 (2), we have from the previous inequality that

$\displaystyle 0\leq\;$ $\displaystyle v^{i}_{*,\pi^{-i}}(s)-v^{i}_{*}(s)$
$\displaystyle\leq\;$
$\displaystyle\gamma\|v^{i}_{*,\pi^{-i}}-v^{i}_{*}\|_{\infty}+2\gamma\|v^{i}-v^{i}_{*}\|_{\infty}$
$\displaystyle+\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}.$
(22)

It remains to bound the last term on the RHS of Eq. (22). Observe that

$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}$
$\displaystyle\leq\;$
$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)-\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\pi^{i}(s))^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}$
$\displaystyle=\;$
$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)+\max_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{-i})^{\top}\mathcal{T}^{-i}(v^{-i})(s)\pi^{i}(s)$
$\displaystyle-\max_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{-i})^{\top}\mathcal{T}^{-i}(v^{-i})(s)\pi^{i}(s)-\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\pi^{i}(s))^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}$
$\displaystyle=\;$
$\displaystyle\sum_{i=1,2}\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)+\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\pi^{i}(s))^{\top}\mathcal{T}^{i}(-v^{-i})(s)\mu^{-i}$
$\displaystyle-\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\pi^{i}(s))^{\top}\mathcal{T}^{i}(v^{i})(s)\mu^{-i}$
$\displaystyle\leq\;$
$\displaystyle\sum_{i=1,2}\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)+\max_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}\left|(\pi^{i}(s))^{\top}(\mathcal{T}^{i}(-v^{-i})(s)-\mathcal{T}^{i}(v^{i})(s))\mu^{-i}\right|$
$\displaystyle\leq\;$
$\displaystyle\sum_{i=1,2}\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v^{i})(s)\pi^{-i}(s)+\gamma\|v^{1}+v^{2}\|_{\infty}$
$\displaystyle\leq\;$
$\displaystyle\frac{9A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}V_{v,s}(\tilde{q}^{1}(s),\tilde{q}^{2}(s))+2\|v^{1}+v^{2}\|_{\infty}+6\tau\log(A_{\max})+\mu$
(23)

where the last line follows from Lemma A.19 (3). Using the previous inequality
in Eq. (22), we have

$\displaystyle 0\leq v^{i}_{*,\pi^{-i}}(s)-v^{i}_{*}(s)\leq\;$
$\displaystyle\gamma\|v^{i}_{*,\pi^{-i}}-v^{i}_{*}\|_{\infty}+2\|v^{i}-v^{i}_{*}\|_{\infty}+2\|v^{1}+v^{2}\|_{\infty}$
$\displaystyle+\frac{9A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\max_{s\in\mathcal{S}}V_{v,s}(\tilde{q}^{1}(s),\tilde{q}^{2}(s))+6\tau\log(A_{\max})+\mu.$

Since the RHS of the previous inequality does not depend on $s$, we have

$\displaystyle\|v^{i}_{*,\pi^{-i}}-v^{i}_{*}\|_{\infty}\leq\;$
$\displaystyle\gamma\|v^{i}_{*,\pi^{-i}}-v^{i}_{*}\|_{\infty}+2\|v^{i}-v^{i}_{*}\|_{\infty}+2\|v^{1}+v^{2}\|_{\infty}$
$\displaystyle+\frac{9A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\max_{s\in\mathcal{S}}V_{v,s}(\tilde{q}^{1}(s),\tilde{q}^{2}(s))+6\tau\log(A_{\max})+\mu,$

which by rearranging terms implies

$\displaystyle\|v^{i}_{*,\pi^{-i}}-v^{i}_{*}\|_{\infty}\leq\;$
$\displaystyle\frac{1}{1-\gamma}\left(2\|v^{i}-v^{i}_{*}\|_{\infty}+2\|v^{1}+v^{2}\|_{\infty}\right.$
$\displaystyle+\left.\frac{9A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\max_{s\in\mathcal{S}}V_{v,s}(\tilde{q}^{1}(s),\tilde{q}^{2}(s))+6\tau\log(A_{\max})+\mu\right).$

Summing up the previous inequality for $i\in\\{1,2\\}$, we obtain

$\displaystyle\sum_{i=1,2}\|v^{i}_{*,\pi^{-i}}-v^{i}_{*}\|_{\infty}\leq\;$
$\displaystyle\frac{2}{1-\gamma}\left(\sum_{i=1,2}\|v^{i}-v^{i}_{*}\|_{\infty}+2\|v^{1}+v^{2}\|_{\infty}\right.$
$\displaystyle+\left.\frac{9A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\max_{s\in\mathcal{S}}V_{v,s}(\tilde{q}^{1}(s),\tilde{q}^{2}(s))+6\tau\log(A_{\max})+\mu\right).$

∎

Combing Lemmas A.1 and A.2 and using our definition of the Lyapunov functions
$\mathcal{L}_{v}(\cdot)$, $\mathcal{L}_{\text{sum}}(\cdot)$, and
$\mathcal{L}_{\theta}(\cdot,\cdot)$ to simplify the notation, we have the
following result.

######  Lemma A.3.

It holds that

$\displaystyle\text{NG}\left(\pi_{T,K}^{1},\pi_{T,K}^{2}\right)\leq\frac{4}{1-\gamma}\left(\mathcal{L}_{v}(T)+2\mathcal{L}_{\text{sum}}(T)+\frac{9A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\mathcal{L}_{\theta}(T,K)+6\tau\log(A_{\max})+\mu\right).$

### A.3 Analysis of the Outer Loop: the $v$-Function Update

######  Lemma A.4.

It holds for all $t\geq 0$ that

$\displaystyle\mathcal{L}_{v}(t+1)\leq\gamma\mathcal{L}_{v}(t)+\frac{17A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\mathcal{L}_{\theta}(t,K)+2\mathcal{L}_{\text{sum}}(t)+2\mathcal{L}_{w}(t,K)^{1/2}+12\tau\log(A_{\max})+\mu.$

###### Proof of Lemma A.4.

For any $t\geq 0$ and $s\in\mathcal{S}$, we have

$\displaystyle\left|v_{t+1}^{i}(s)-v_{*}^{i}(s)\right|=\;$
$\displaystyle\left|(\pi_{t,K}^{i}(s))^{\top}\lceil q_{t,K}^{i}(s)\rceil-
v_{*}^{i}(s)\right|$ (Algorithm 3 Line $8$) $\displaystyle\leq\;$
$\displaystyle\left|\mathcal{B}^{i}(v_{t}^{i})(s)-\mathcal{B}^{i}(v_{*}^{i})(s)\right|+\left|(\pi_{t,K}^{i}(s))^{\top}\lceil
q_{t,K}^{i}(s)\rceil-\mathcal{B}^{i}(v_{t}^{i})(s)\right|$ ($v_{*}^{i}$ is the
unique fixed point of $\mathcal{B}^{i}(\cdot)$) $\displaystyle\leq\;$
$\displaystyle\gamma\|v_{t}^{i}-v_{*}^{i}\|_{\infty}+\left|(\pi_{t,K}^{i}(s))^{\top}\lceil
q_{t,K}^{i}(s)\rceil-\mathcal{B}^{i}(v_{t}^{i})(s)\right|,$ (24)

where the last line follows from the contraction property of the minimax
Bellman operator $\mathcal{B}^{i}(\cdot)$. Next, we bound the second term on
the RHS of the previous inequality. Observe that

$\displaystyle\left|(\pi_{t,K}^{i}(s))^{\top}\lceil
q_{t,K}^{i}(s)\rceil-\mathcal{B}^{i}(v_{t}^{i})(s)\right|$
$\displaystyle\leq\;$ $\displaystyle\left|(\pi_{t,K}^{i}(s))^{\top}\lceil
q_{t,K}^{i}(s)\rceil-(\pi_{t,K}^{i}(s))^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)\right|$
$\displaystyle+\left|(\pi_{t,K}^{i}(s))^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)-(\pi_{t,K}^{i}(s))^{\top}\tilde{q}_{t,K}^{i}(s)\right|$
$\displaystyle+\left|(\pi_{t,K}^{i}(s))^{\top}\tilde{q}_{t,K}^{i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\tilde{q}_{t,K}^{i}(s)\right|$
$\displaystyle+\left|\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\tilde{q}_{t,K}^{i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)\right|$
$\displaystyle+\left|\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\mu^{-i}\right|$
$\displaystyle\leq\;$
$\displaystyle\left\|q_{t,K}^{i}(s)-\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)\right\|_{\infty}+2\left\|\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)-\tilde{q}_{t,K}^{i}(s)\right\|_{\infty}$
$\displaystyle+\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i}-\pi_{t,K}^{i}(s))^{\top}\tilde{q}_{t,K}^{i}(s)$
$\displaystyle+\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\mu^{-i}.$
(25)

We next bound each term on the RHS of the previous inequality. For the first
term, we have

$\displaystyle\left\|q_{t,K}^{i}(s)-\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)\right\|_{\infty}^{2}=\;$
$\displaystyle\left\|q_{t,K}^{i}(s)-\bar{q}_{t,K}^{i}(s)\right\|_{\infty}^{2}$
$\displaystyle\leq\;$
$\displaystyle\left\|q_{t,K}^{i}-\bar{q}_{t,K}^{i}\right\|_{\infty}^{2}$
$\displaystyle=\;$
$\displaystyle\left\|\Phi^{i}(w_{t,K}^{i}-\bar{w}_{t,K}^{i})\right\|_{\infty}^{2}$
$\displaystyle\leq\;$
$\displaystyle\|\Phi^{i}\|_{2,\infty}\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}^{2}$
$\displaystyle\leq\;$
$\displaystyle\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}^{2},$ (26)

where the last line follows from
$\|\Phi^{i}\|_{2,\infty}=\max_{s,a^{i}}\|\phi^{i}(s,a^{i})\|_{2}\leq 1$. For
the second term on the RHS of Eq. (25), we have

$\displaystyle
2\left\|\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)-\tilde{q}_{t,K}^{i}(s)\right\|_{\infty}\leq\;$
$\displaystyle
2\left\|\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)-\tilde{q}_{t,K}^{i}(s)\right\|_{2}$
$\displaystyle\leq\;$
$\displaystyle\frac{4\sqrt{2}A_{\max}}{\sqrt{\tau}(1-\gamma)}V_{v_{t},s}^{1/2}(\tilde{q}_{t,K}^{1}(s),\tilde{q}_{t,K}^{2}(s))$
(Lemma A.12 and
$\|\mathcal{T}^{i}(v_{t}^{i})(s)\|_{2}\leq\frac{A_{\max}}{1-\gamma}$)
$\displaystyle\leq\;$
$\displaystyle\tau+\frac{8A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}V_{v_{t},s}(\tilde{q}_{t,K}^{1}(s),\tilde{q}_{t,K}^{2}(s)),$
(27)

where the last inequality follows from $a^{2}+b^{2}\geq 2ab$. For the third
term on the RHS of Eq. (25), since

$\displaystyle\pi_{t,K}^{i}(s)=\sigma_{\tau}(\tilde{q}_{t,K}^{i}(s))={\arg\max}_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\\{(\mu^{i})^{\top}\tilde{q}_{t,K}^{i}(s)+\tau\nu(\mu^{i})\\},$
(Algorithm 3 Line $4$)

we have

$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i}-\pi_{t,K}^{i}(s))^{\top}\tilde{q}_{t,K}^{i}(s)$
$\displaystyle\leq\;$
$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\left\\{(\mu^{i}-\pi_{t,K}^{i}(s))^{\top}\tilde{q}_{t,K}^{i}(s)+\tau\nu(\mu^{i})-\tau\nu(\pi_{t,K}^{i}(s))\right\\}+2\tau\log(A_{\max})$
$\displaystyle=\;$ $\displaystyle 2\tau\log(A_{\max}).$ (28)

For the last term on the RHS of Eq. (25), we have by Eq. (23) that

$\displaystyle\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)-\max_{\mu^{i}\in\Delta(\mathcal{A}^{i})}\min_{\mu^{-i}\in\Delta(\mathcal{A}^{-i})}(\mu^{i})^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\mu^{-i}$
$\displaystyle\leq\;$
$\displaystyle\frac{9A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}V_{v_{t},s}(\tilde{q}_{t,K}^{1}(s),\tilde{q}_{t,K}^{2}(s))+2\|v_{t}^{1}+v_{t}^{2}\|_{\infty}+6\tau\log(A_{\max})+\mu.$
(29)

Using the bounds we obtained in Eqs. (26), (27), (28), and (29) together in
Eq. (25), we have

$\displaystyle\left|(\pi_{t,K}^{i}(s))^{\top}\lceil\tilde{q}_{t,K}^{i}(s)\rceil-\mathcal{B}^{i}(v_{t}^{i})(s)\right|\leq\;$
$\displaystyle\frac{17A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}V_{v_{t},s}(\tilde{q}_{t,K}^{1}(s),\tilde{q}_{t,K}^{2}(s))+2\|v_{t}^{1}+v_{t}^{2}\|_{\infty}$
$\displaystyle+\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}+12\tau\log(A_{\max})+\mu.$

Using the previous inequality in Eq. (24), we have

$\displaystyle\|v_{t+1}^{i}-v_{*}^{i}\|_{\infty}\leq\;$
$\displaystyle\gamma\|v_{t}^{i}-v_{*}^{i}\|_{\infty}+\frac{17A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\mathcal{L}_{\theta}(t,K)+2\mathcal{L}_{\text{sum}}(t)+\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}$
$\displaystyle+12\tau\log(A_{\max})+\mu.$

Summing up the resulting inequality for $i\in\\{1,2\\}$, we obtain

$\displaystyle\mathcal{L}_{v}(t+1)\leq\;$
$\displaystyle\gamma\mathcal{L}_{v}(t)+\frac{17A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\mathcal{L}_{\theta}(t,K)+2\mathcal{L}_{\text{sum}}(t)+\sum_{i=1,2}\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}+12\tau\log(A_{\max})+\mu$
$\displaystyle\leq\;$
$\displaystyle\gamma\mathcal{L}_{v}(t)+\frac{17A_{\max}^{2}}{\tau^{2}(1-\gamma)^{2}}\mathcal{L}_{\theta}(t,K)+2\mathcal{L}_{\text{sum}}(t)+2\mathcal{L}_{w}(t,K)^{1/2}+12\tau\log(A_{\max})+\mu,$

where the last line follows from

$\displaystyle\sum_{i=1,2}\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}=\sum_{i=1,2}(\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}^{2})^{1/2}\leq
2\bigg{(}\sum_{i=1,2}\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}^{2}\bigg{)}^{1/2}.$
(30)

∎

######  Lemma A.5.

It holds for all $t\geq 0$ that

$\displaystyle\mathcal{L}_{\text{sum}}(t+1)\leq\gamma\mathcal{L}_{\text{sum}}(t)+2\mathcal{L}_{w}(t,K)^{1/2}.$

###### Proof of Lemma A.5.

For any $t\geq 0$ and $s\in\mathcal{S}$, we have

$\displaystyle\left|\sum_{i=1,2}v_{t+1}^{i}(s)\right|=\;$
$\displaystyle\left|\sum_{i=1,2}(\pi_{t,K}^{i}(s))^{\top}\lceil
q_{t,K}^{i}(s)\rceil\right|$ (Algorithm 3 Line $8$) $\displaystyle=\;$
$\displaystyle\left|\sum_{i=1,2}(\pi_{t,K}^{i}(s))^{\top}\lceil
q_{t,K}^{i}(s)\rceil-(\pi_{t,K}^{i}(s))^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)\right|$
$\displaystyle+\left|\sum_{i=1,2}(\pi_{t,K}^{i}(s))^{\top}\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)\right|$
$\displaystyle\leq\;$
$\displaystyle\gamma\|v_{t}^{1}+v_{t}^{2}\|_{\infty}+\sum_{i=1,2}\left\|q_{t,K}^{i}(s)-\mathcal{T}^{i}(v_{t}^{i})(s)\pi_{t,K}^{-i}(s)\right\|_{\infty}$
(Lemma A.19) $\displaystyle\leq\;$
$\displaystyle\gamma\|v_{t}^{1}+v_{t}^{2}\|_{\infty}+\sum_{i=1,2}\|w_{t,K}^{i}-\bar{w}_{t,K}^{i}\|_{2}$
(Eq. (26)) $\displaystyle\leq\;$
$\displaystyle\gamma\|v_{t}^{1}+v_{t}^{2}\|_{\infty}+2\mathcal{L}_{w}(t,K)^{1/2}.$

where the last line follows from Eq. (30). Since the RHS of the previous
inequality does not depend on $s$, we have the desired result. ∎

### A.4 Analysis of the Inner Loop

In this section, we will focus on analyzing the update equations of the two-
timescale iterates $\\{w_{t,k}^{i}\\}$ and $\\{\theta_{t,k}^{i}\\}$ updated in
the inner loop of Algorithm 3. We begin by writing only the inner loop of
Algorithm 3 in Algorithm 4, where we omit the subscript $t$. The results
derived for Algorithm 4 can be easily combined with the analysis of the outer
loop in previous sections using the tower property of conditional expectations
and the Markov property.

Algorithm 4 The Inner Loop of Algorithm 3

1: Input: Integer $K$, initializations $w_{0}^{i},\theta_{0}^{i}$ satisfying
$\|w_{0}^{i}\|_{2}\leq M$, $\|\theta_{0}^{i}\|_{2}\leq M$, and a value
function $v^{i}$ (satisfying $\|v^{i}\|_{\infty}\leq 1/(1-\gamma)$) from the
outer loop.

2: for $k=0,1,\cdots,K-1$ do

3:
$\tilde{q}^{i}_{k+1}=\tilde{q}_{k}^{i}+\beta_{k}(q_{k}^{i}-\tilde{q}_{k}^{i})$
and compute $\pi_{k+1}^{i}(S_{k})=\sigma_{\tau}(\tilde{q}_{k+1}^{i}(S_{k}))$

4: Play $A_{k}^{i}\sim\pi_{k+1}^{i}(\cdot\mid S_{k})$ (against $A_{k}^{-i}$),
and observe $S_{k+1}\sim p(\cdot\mid S_{k},A_{k}^{i},A_{k}^{-i})$

5:
$w_{k+1}^{i}=\text{Proj}_{M}^{i}[w_{k}^{i}+\alpha_{k}\phi^{i}(S_{k},A_{k}^{i})(R_{i}(S_{k},A_{k}^{i},A_{k}^{-i})+\gamma
v^{i}(S_{k+1})-\phi^{i}(S_{k},A_{k}^{i})^{\top}w_{k}^{i})]$

6: end for

#### A.4.1 Analysis of the Slow-Timescale Iterates

Note that Algorithm 4 Line $3$ is equivalent to

$\displaystyle\tilde{q}^{i}_{k+1}(s)=\tilde{q}_{k}^{i}(s)+\beta_{k}(\mathcal{T}^{i}(v^{i})(s)\sigma_{\tau}(\tilde{q}_{k}^{-i}(s))-\tilde{q}_{k}^{i}(s))+\beta_{k}(q_{k}^{i}(s)-\mathcal{T}^{i}(v^{i})(s)\sigma_{\tau}(\tilde{q}_{k}^{-i}(s))),\;\;\forall\,s\in\mathcal{S}.$
(31)

The last term on the RHS of the previous inequality represents the estimation
error of the fast-timescale iterates $\\{w_{k}^{i}\\}$, which will be analyzed
in Appendix A.4.3. Here we focus on the evolution of the $q$-function
$\tilde{q}_{k}^{i}$ associated with the slow-timescale iterate
$\theta_{k}^{i}$, which we use to compute the policies. In view of Eq. (31),
it is enough to study the following joint update equation:

$\displaystyle x_{k+1}=\;$ $\displaystyle
x_{k}+\beta_{k}(B\sigma_{\tau}(y_{k})-x_{k}+\mathcal{E}_{k,x}),$ (32)
$\displaystyle y_{k+1}=\;$ $\displaystyle
y_{k}+\beta_{k}(A\sigma_{\tau}(x_{k})-y_{k}+\mathcal{E}_{k,y}),$ (33)

where $x_{k},\mathcal{E}_{k,x}\in\mathbb{R}^{n}$,
$y_{k},\mathcal{E}_{k,y}\in\mathbb{R}^{m}$, $A\in\mathbb{R}^{m\times n}$, and
$B\in\mathbb{R}^{n\times m}$. Note that we do not assume $A+B^{\top}=0$. We
next use a Lyapunov approach to study the behavior of $(x_{k},y_{k})$. Let

$\displaystyle V_{1}(x,y)=\;$
$\displaystyle\max_{u^{\prime}\in\Delta^{m}}\min_{u\in\Delta^{m}}\left\\{(u^{\prime}-u)^{\top}y+\tau\nu(u^{\prime})-\tau\nu(u)+\frac{1}{2\mu}\|x-Bu\|_{2}^{2}\right\\},$
$\displaystyle V_{2}(x,y)=\;$
$\displaystyle\max_{v^{\prime}\in\Delta^{n}}\min_{v\in\Delta^{n}}\left\\{(v^{\prime}-v)^{\top}x+\tau\nu(v^{\prime})-\tau\nu(v)+\frac{1}{2\mu}\|y-Av\|_{2}^{2}\right\\},$

where $\mu>0$ is a tunable parameter. Let
$V:\mathbb{R}^{n}\times\mathbb{R}^{m}\mapsto\mathbb{R}$ be the Lyapunov
function defined as $V(x,y)=V_{1}(x,y)+V_{2}(x,y)$. For simplicity of
notation, we denote

$\displaystyle p(\mu,x,y)=\;$
$\displaystyle{\arg\min}_{u\in\Delta^{m}}\left\\{-u^{\top}y-\tau\nu(u)+\frac{1}{2\mu}\|x-Bu\|_{2}^{2}\right\\},$
$\displaystyle q(\mu,x,y)=\;$
$\displaystyle{\arg\min}_{v\in\Delta^{n}}\left\\{-v^{\top}x-\tau\nu(v)+\frac{1}{2\mu}\|y-Av\|_{2}^{2}\right\\},$

both of which are well-defined due to strong convexity. The following
proposition shows that the joint update equations (32) and (33) produce a
negative drift with respect to the Lyapunov function $V(\cdot,\cdot)$.

######  Proposition A.1.

Suppose that there exists $L_{b}\geq 0$ such that
$\|B\sigma_{\tau}(y_{k})-x_{k}+\mathcal{E}_{k,x}\|_{2}\leq L_{b}$ and
$\|A\sigma_{\tau}(x_{k})-y_{k}+\mathcal{E}_{k,y}\|_{2}\leq L_{b}$ for all
$k\geq 0$. Then, by choosing $\mu=\tau/64$, we have for all $k\geq 0$ that

$\displaystyle V(x_{k+1},y_{k+1})\leq\;$
$\displaystyle\left(1-\frac{\beta_{k}}{2}\right)V(x_{k},y_{k})+\frac{520\beta_{k}}{\tau}(\|\mathcal{E}_{k,x}\|_{2}^{2}+\|\mathcal{E}_{k,y}\|_{2}^{2})$
$\displaystyle+4\beta_{k}\max_{i,j}|A(i,j)+B(j,i)|+\frac{138L_{b}\beta_{k}^{2}}{\tau}.$

######  Remark.

By letting $A=X_{1}$, $B=X_{2}$, $\mathcal{E}_{k,x}=0$, and
$\mathcal{E}_{k,y}=0$, Proposition A.1 implies Proposition 5.1.

###### Proof of Proposition A.1.

A sequence of properties regarding the Lyapunov function $V(\cdot)$ is
provided in Appendix A.4.2. For any $k\geq 0$, we have by the smoothness
property of the Lyapunov function (cf. Lemma A.8) that

$\displaystyle V(x_{k+1},y_{k+1})-V(x_{k},y_{k})\leq\;$
$\displaystyle\underbrace{\langle\nabla_{x}V(x_{k},y_{k}),x_{k+1}-x_{k}\rangle}_{E_{1}}+\underbrace{\frac{L_{V}}{2}\|x_{k+1}-x_{k}\|_{2}^{2}}_{E_{2}}$
$\displaystyle+\underbrace{\langle\nabla_{y}V(x_{k},y_{k}),y_{k+1}-y_{k}\rangle}_{E_{3}}+\underbrace{\frac{L_{V}}{2}\|y_{k+1}-y_{k}\|_{2}^{2}}_{E_{4}}.$
(34)

We next bound all the terms on the RHS of the previous inequality. For the
term $E_{2}$, we have

$\displaystyle
E_{2}=\frac{L_{V}}{2}\|x_{k+1}-x_{k}\|_{2}^{2}=\frac{L_{V}\beta_{k}^{2}}{2}\|B\sigma_{\tau}(y_{k})-x_{k}+\mathcal{E}_{k,x}\|_{2}^{2}\leq\frac{L_{V}L_{b}^{2}\beta_{k}^{2}}{2}.$

Similarly, we also have $E_{4}\leq\frac{L_{V}L_{b}^{2}\beta_{k}^{2}}{2}$.
Next, we consider the term $E_{1}$. Using Eq. (32) and Lemma A.7, we have

$\displaystyle E_{1}=\;$
$\displaystyle\langle\nabla_{x}V(x_{k},y_{k}),x_{k+1}-x_{k}\rangle$
$\displaystyle=\;$ $\displaystyle\beta_{k}\underbrace{\frac{1}{\mu}\langle
x_{k}-Bp(\mu,x_{k},y_{k}),B\sigma_{\tau}(y_{k})-x_{k}\rangle}_{E_{1,1}}+\beta_{k}\underbrace{\frac{1}{\mu}\langle
x_{k}-Bp(\mu,x_{k},y_{k}),\mathcal{E}_{k,x}\rangle}_{E_{1,2}}$
$\displaystyle+\beta_{k}\underbrace{\langle\sigma_{\tau}(x_{k})-q(\mu,x_{k},y_{k}),B\sigma_{\tau}(y_{k})-x_{k}\rangle}_{E_{1,3}}+\beta_{k}\underbrace{\langle\sigma_{\tau}(x_{k})-q(\mu,x_{k},y_{k}),\mathcal{E}_{k,x}\rangle}_{E_{1,4}}$
(35)

We next bound the terms $\\{E_{1,i}\\}_{1\leq i\leq 4}$. For the term
$E_{1,1}$, we have

$\displaystyle E_{1,1}=\;$ $\displaystyle\frac{1}{\mu}\langle
x_{k}-Bp(\mu,x_{k},y_{k}),B\sigma_{\tau}(y_{k})-Bp(\mu,x_{k},y_{k})+Bp(\mu,x_{k},y_{k})-x_{k}\rangle$
$\displaystyle=\;$ $\displaystyle\frac{1}{\mu}\langle
B^{\top}(x_{k}-Bp(\mu,x_{k},y_{k})),\sigma_{\tau}(y_{k})-p(\mu,x_{k},y_{k})\rangle-\frac{1}{\mu}\|x_{k}-Bp(\mu,x_{k},y_{k})\|_{2}^{2}$
$\displaystyle=\;$ $\displaystyle\langle-
y_{k}-\tau\nabla\nu(p(\mu,x_{k},y_{k})),\sigma_{\tau}(y_{k})-p(\mu,x_{k},y_{k})\rangle-\frac{1}{\mu}\|x_{k}-Bp(\mu,x_{k},y_{k})\|_{2}^{2}$
(Lemma A.9) $\displaystyle=\;$ $\displaystyle-