Solving Zero-Sum Markov Games with Continous State via Spectral Dynamic Embedding

Authors: chenhao  zhou
zebang  shen
chao  zhang
hanbin  zhao
hui  qian
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-14 16:45:22

Abstract: In this paper, we propose a provably efficient natural policy gradient algorithm called Spectral Dynamic Embedding Policy Optimization (\SDEPO) for two-player zero-sum stochastic Markov games with continuous state space and finite action space.
  In the policy evaluation procedure of our algorithm, a novel kernel embedding method is employed to construct a finite-dimensional linear approximations to the state-action value function.
  We explicitly analyze the approximation error in policy evaluation, and show that \SDEPO\  achieves an $\tilde{O}(\frac{1}{(1-\gamma)^3\epsilon})$ last-iterate convergence to the $\epsilon-$optimal Nash equilibrium, which is independent of the cardinality of the state space.
  The complexity result matches the best-known results for global convergence of policy gradient algorithms for single agent setting.
  Moreover, we also propose a practical variant of \SDEPO\ to deal with continuous action space and empirical results demonstrate the practical superiority of the proposed method.

Main paper: 
 introduction: 
 Two-player zero-sum stochastic Markov games (2p0s-MGs) has been the focus of research across a range of research communities. In this problem, two players select their actions based on the current state simultaneously and independently. Player one aims to maximize the return based on the reward provided by the environment, while player two aims to minimize it. For 2p0s-MGs with finite state space, tabular methods [Alacaoglu et al., 2022 , Bai and Jin, 2020 , Daskalakis et al., 2020 , Wei et al., 2021 , Zhao et al., 2022] represent the state-action value function with tables, which results in a sample complexity depending on the cardinality of the state spaces. To deal with 2p0s-MGs with complex state space, researchers recently employ function approximations of the state-action value function to deal with large-scale discrete/continuous state space, including linear function approximations [Xie et al., 2020 , Chen et al., 2022] , kernel function approximations [Junchi Li et al., 2022 , Qiu et al., 2021] and general function classes such as neural networks [Jin et al., 2022 , Huang et al., 2021] . Basically, these methods use samples to construct the function approximations and update the policies with value iteration. Theoretical analyses show that these methods possess sample complexities independent of the state space's cardinality. Note that these methods fail to explicitly utilize the dynamics information of the underlying environments and only achieve a sub-optimal sample complexity Õ(ϵ -2 ) to find an ϵ-optimal Nash equilibrium for 2p0s-MGs with known system dynamics.In this paper, we introduce a spectral dynamic embedding method, which explicitly uses the dynamics information to approximate the state-action value functions, and propose an efficient naturalpolicy-gradient-type algorithm, called Spectral Dynamic Embedding Policy optimization (SDEPO), Table 1 : Comparison of policy optimization methods for finding an ϵ-optimal NE of two-player zero-sum episodic Markov games in terms of the duality gap. Here, H refers to the horizon length and d is the dimension of their features. For simplicity, we ignore the problem-dependent constant.for 2p0s-MGs. In particular, the spectral dynamic embedding method directly constructs truncated linear approximations to the transition dynamics of a Markov game in a kernel space, and implements dynamic programming to calculate the state-action value function approximation. The superiority of spectral dynamic embedding has been justified in single agent setting [Ren et al., 2022 [Ren et al., , 2023]] . We leverage two kernel feature generation methods for the truncated approximation, namely random feature generation and Nyström feature generation, and analyze the approximation error of these two methods during policy evaluation. Our contributions lie in the following folds.1. We present a truncated kernel-based linearization method for the state-action value function approximation in two-player zero-sum Markov games with continuous state space. With the random/Nyström feature generation, this method automatically generates truncated kernel representation from system dynamics, bypassing the difficulty of kernel feature decision existed in other kernel approximation methods. By integrating the acquired kernel features into the temporal difference learning process, we estimate the state-action value functions through the least square policy evaluation. Leveraging the obtained value functions, policy improvement can be achieved through the natural policy gradient descent/ascent approach. 2. We establish rigorous analysis of the approximation error in the truncated value-function approximation and the statistical error induced in policy improvement procedure with finite samples. Our theoretical analysis demonstrates that SDEPO achieves a near-optimal Õ( 1 (1-γ) 3 ϵ ) last-iterate convergence to the ϵ-optimal Nash equilibrium for 2p0s-MGs with continuous state space and finite action space, where γ represents the discounted factor. This complexity result matches the best-known results for the policy gradient algorithm to achieve global convergence in the single agent setting. Moreover, we also propose a practical variant of SDEPO to deal with continuous action space and empirical results demonstrate the practical superior performance of the proposed method.1.1 Related works RL methods for 2p0s-MGs. There is a large body of literature on MARL for two-player MGs. Alacaoglu et al. [2022] , Daskalakis et al. [2020] , Zhao et al. [2022] focus on the tabular setting, i.e., the state-action space can be represented by a table with moderate size To address the challenge of continuous state spaces, many researchers have designed algorithms based on function approximation. Xie et al. [2020] , Chen et al. [2022] investigated methods based on linear function approximation, using a set of linear features to represent the state transition function and reward function. Jin et al. [2022] , Huang et al. [2021] employed general function approximation for MGs with low multi-agent Bellman eluder dimension and MGs with a finite minimax Eluder dimension. Although these algorithms yield strong theoretical guarantees, they are computationally inefficient. Junchi Li et al. [2022] , Qiu et al. [2021] studied learning MGs with kernel approximation. Their approachs assume that there are a set of (possibly infinite-dimensional) kernel features that span the transition function or value function space. However, finding a good set of kernel features is a very challenging task, making their assumption difficult to satisfy. Additionally, infinite-dimensional kernel features are infeasible to compute, so their method requires finite-dimensional approximation, and the errors caused by finite-dimensional approximation currently lack analysis. Tabular methods for 2p0s-MGs where system dynamics are known. As shown in Table 1 , a parallel line of research aims to solve 2p0s-MGs with system dynamics known. For the infinite-horizon discounted setting, Wei et al. [2021] proposed an optimistic gradient descent ascent (OGDA) method which achieves a last-iterate convergence at an O |S| 3(1-γ) 9 ϵ 2 iteration complexity. Cen et al. [2021] established linear last-iterate convergence of entropy-regularized OMWU. For the finite-horizon episodic setting, Zhang et al. [2022] showed that the modified optimistic Follow-The-Regularized-Leader method finds an ϵ-optimal NE in O H 4 ϵ iterations and Cen et al. [2023] proposed a singleloop policy optimization algorithm that implies the last-iterate convergence with an iteration complexity of O H 3 ϵ . However, their methods are all confined to tabular setting, whereas our method can handle problems in continuous state spaces. RL with Function Approximation. Function approximation in single-agent RL has been extensively studied in recent years to achieve a better sample complexity that depends on the complexity of function approximators rather than the size of the state/action space. One line of work studies RL with linear function approximation [Yang and Wang, 2019, Jin et al., 2020] . Typically, these methods assume the optimal value function can be well approximated by linear functions, and achieve polynomial sample efficiency guarantees related to feature dimension under certain regularity conditions. Another line of works studied the MDPs with general nonlinear function approximations [Jiang et al., 2017 , Jin et al., 2021] . Jiang et al. [2017] , Jin et al. [2021] present algorithms with PAC guarantees for problems with low Bellman rank and low BE dimension, respectively. We note that MGs are inherently more complex than MDPs due to their min-max nature and it is generally difficult to directly extend these results to the dual-player dynamic setting of MGs. 
background and preliminaries: 
 In this section, we present the necessary definitions that will be adopted throughout the paper. In Section 2.1, we formally describe the setup for two-player zero-sum stochastic Markov games with simultaneous moves. In Section 2.2, We briefly introduce the background knowledge about positive definite kernels and their decompositions.A two-player zero-sum stochastic Markov games with simultaneous moves is defined by the tuple (S, A 1 , A 2 , r, P, γ), where S is the state space, A i is the set of actions that player i ∈ {1, 2} can take, r is the reward function, P is the transition function and γ ∈ [0, 1) is the discounted factor. At each step t, given the state s ∈ S, P1 and P2 take actions a ∈ A 1 and b ∈ A 2 , respectively, and then both receive the reward r(s, a, b). The system then shifts to a new state s ′ ∼ P(•|s, a, b) according to the transition kernel. Throughout this paper, we assume for simplicity that S = R d , A 1 = A 2 = A and that the rewards r(s, a, b) are deterministic functions of the tuple (s, a, b) taking value in [-1, 1]; generalization to the setting with A 1 ̸ = A 2 and stochastic rewards is straightforward.Denote by ∆ ≡ ∆(A) the probability simplex over the action space A. A stochastic policy of P1 is a sequence of functions π := (π t : S → ∆) t . At each step t and state s ∈ S, P1 takes an action sampled from the distribution π t (s) over A. Similarly, a stochastic policy of P2 is given by the sequence π := (π t : S → ∆) t .For a fixed pair of policies (π, π) for both players, the value and Q (a.k.a. action-value) functions for the above game can be defined as following:V π,π (s) : = E ∞ t=0 γ t r(s t , a t , b t )|s 0 = s , Q π,π (s, a, b) : = E ∞ t=0 γ t r(s t , a t , b t )|s 0 = s, a 0 = a, b 0 = b ,where the expectation is overa t ∼ π t (•|s t ), b t ∼ π t (•|s t ) and s t+1 ∼ P(•|s t , a t , b t ).In the zero-sum setting, for a given initial state s 0 , P1 seeks to maximize V π,π (s 0 ) whereas P2 aims to minimize it. Accordingly, we introduce the value and Q functions when P1 plays the best response to a fixed policy π of P2:V * ,π (s) = max π V π,π (s) and Q * ,π (s, a, b) = max π Q π,π (s, a, b).Similarly, when P2 plays the best response to P1's policy π, we defineV π, * (s) = min π V π,π (s) and Q π, * (s, a, b) = min π Q π,π (s, a, b).A Nash Equilibrium (NE) of the game is a pair of stochastic policies (π * , π * ) that are the best response to each other, which we write as V π * ,π * (s) = V * ,π * (s) = V π * , * (s). NE always exists for discounted two-player zero-sum Markov Games [Filar and Vrieze, 2012] . Correspondingly, let a, b) denote the values of the NE. In practice, we always seek to find an ϵ-optimal Nash equilibrium, which is a pair of stochastic policies (π, π) that satisfiesV * (s) := V π * ,π * (s) and Q * (s, a, b) := Q π * ,π * (s,E s∼µ0 max π ′ ,π ′ V π ′ ,π (s) -V π,π ′ (s) ≤ ϵ. for a initial state distribution µ 0 .We are interested in finding a one-sided ϵ-optimal Nash equilibrium, similar to Alacaoglu et al. [2022] , Zhao et al. [2022] , Daskalakis et al. [2020] , Zhang et al. [2019] . In particular, for the initial state distribution µ 0 , we seek π out such thatE s∼µ0 max π V π,π out (s) -V * (s) ≤ ϵ.We assume that the transition function satisfies the following assumption. Assumption 1. For each (s t , a t , b t ) ∈ S × A × A, we assume thats t+1 = f (s t , a t , b t ) + ϵ t , where ϵ t ∼ N (0, σ 2 I d ).The function f : S × A × A → S describes the general dynamics and {ϵ t } ∞ t=0 are independent Gaussian noises. In other words,P(s t+1 |s t , a t , b t ) ∝ exp - ∥f (s t , a t , b t ) -s t+1 ∥ 2 2 2σ 2 (1)To efficiently represent the continuous state space of 2p0s-MGs which can not be handled by traditional tabular methods, it is common to embed the continuous state space into a kernel space. A widely used kernel is the positive definite (PD) kernel.Definition 1 ((Positive-Definite) Kernel [Mohri, 2018] ). A symmetric function k : X × X → R is said to be a positive definite kernel if for any {x 1 , . . . , x m } ⊂ X , the matrixK = [k(x i , x j )] ij ∈ R m×m is symmetric positive-definite.PD kernels admit many decompositions, such as Bochner decomposition [Devinatz, 1953] , Mercer decomposition [Mercer, 1909] , Canonical decomposition [Stochel, 1992] and Kolmogorov decomposition [Ghaemi et al., 2021] . Among these decompositions, Bochner decomposition [Devinatz, 1953] and Mercer decomposition [Mercer, 1909] have recently draw significant attention since they lead to efficient, low-dimensional approximations of the kernel, reducing computational complexity [Liu et al., 2021] . Definition 2 (Bochner decomposition [Rudin, 2017] and Mercer decomposition [Mercer, 1909] ).Let X ⊂ R d be a compact domain, µ a strictly positive Borel measure on X , and k(x, x ′ ) = G(x -x ′ ) a bounded continuous shift-invariant positive definite kernel. Then, k(x, x ′ ) admits Bochner decomposition, i.e. there exists a non-negative measure ω, such thatEQUATIONand Mercer decomposition, i.e. there exists a countable orthonormal basisEQUATIONwhere the convergence is absolute and uniform for all (x, x ′ ) ∈ X × X . Without loss of generality, we assumeσ 1 ≥ σ 2 ≥ • • • > 0.It can be verified that the Gaussian kernel, k(x, x ′ ) = exp(- ∥x-x ′ ∥ 2 2σ 2), meets the conditions in the above definition, and it admits both a Bochner decomposition and a Mercer decomposition [Ren et al., 2023] . Note that according to the Bochner/Mercer decomposition, the kernel can be represented with infinite basis. Practically, one can construct finite-dimensional approximations to the positive-definite kernel based on the Bochner/Mercer decompositions with random/Nyström features [Rahimi and Recht, 2008, Williams and Seeger, 2000] , respectively. Recently, Ren et al. [2023] utilize a finite-dimensional approximation to represent the environment in stochastic nonlinear control problems and showed its superior performance, motivating our work. 
policy optimization with spectral dynamic embedding: 
 In this section, we begin by introducing spectral dynamic embedding to represent the environment of the Markov game. This approach allows us to express the Q-function for any policy pair with infinite dimensional features. Subsequently, we develop finite-dimensional approximated features for computational tractability. With these features, each player conducts least square policy evaluation to estimate the Q-function of current policy pair based on the generated features and then improve the policy by natural policy gradient based on the estimated Q-function. This leads to Spectral Dynamic Embedding Policy Optimization (SDEPO) in Algorithm 3.By interpreting the state transition function (1) as a Gaussian kernel, we can decompose the transition function of a Markov game using Bochner decomposition and Mercer decomposition, as detailed below. Lemma 1 (Spectral Dynamic Embedding). Consider any α ∈ [0, 1). Denotek α (x, x ′ ) = exp -(1-α 2 )∥x-x ′ ∥ 2 2σ 2for any 0 ≤ α < 1. We can decompose k α (x, x ′ ) using Bochner decomposition and Mercer decomposition.Letψ ω (s, a, b) = g α (f (s, a, b)) α d cos ω ⊤ f (s, a, b) √ 1 -α 2 , sin ω ⊤ f (s, a, b) √ 1 -α 2 , χ ω (s ′ ) = p α (s ′ )[cos( 1 -α 2 ω ⊤ s ′ ), sin( 1 -α 2 ω ⊤ s ′ )] ⊤ , where g α (f (s, a, b)) := exp α 2 ∥f (s,a,b)∥ 2 2(1-α 2 )σ 2, ω ∼ N (0, σ -2 I d ), and p α (s ′ ) =α d (2πσ 2 ) d/2 exp -∥αs ′ ∥ 2 2σ 2is a Gaussian distribution for s ′ with standard deviation σ α . Using Bochner decomposition, we haveP (s ′ |s, a, b) =E ω∼N (0,σ -2 I d ) ψ ω (s, a, b) ⊤ χ ω (s ′ ) := ⟨ψ ω (s, a, b), χ ω (s ′ )⟩ N (0,σ -2 I d ) .(4)Let µ be a strictly positive Borel measure on S. By Mercer's theorem, k α admits a decomposition for any (x, x ′ ) ∈ X [Fasshauer, 2011] :EQUATION)The tunable parameter α offers advantages for theoretical analysis and can also be leveraged to enhance empirical performance.Let [m] where ω i ∼ N (0, σ -2 I d ) and construct the featureϕ ω (•) = [ψ ω (•), r(•)] and ϕ M (•) = [ψ M (•), r(•)].EQUATION) set ϕ(s, a, b) = [ψ rf (s, a, b), r(s, a, b)].Algorithm 2 Nyström Features Generation Data: Transition Model s ′ = f (s, a, b) + ε where ε ∼ N (0, σ 2 I d ),(nNys) i,j = k α (x i , x j ).Perform eigendecomposition on the Gram matrix K (nnys) U = ΛU , with λ 1 ≥ • • • ≥ λ nnys denoting the corresponding eigenvalues. Construct the featureEQUATION) set ϕ(s, a, b) = [ψ nys (s, a, b), r(s, a, b)].Lemma 2. For any policy pair, there exist weights {θ π,π ω } (where ω ∼ N (0, σ -2 I d )) and θ π,π M ∈ ℓ 2 such that the corresponding value function satisfiesQ π,π (s, a, b) = ϕ ω (s, a, b), θ π,π ω N (0,σ -2 I d ) = ϕ M (s, a, b), θ π,π M ℓ2We provide the proof of this section in Appendix A.Basically, a desirable feature embedding method should provide a good approximation to the underlying dynamics of MGs with a modest feature dimension. However, the dimension of both Bochner and Mercer spectral dynamic embedding is infinite, which is computationally intractable, motivating us to investigate the finite-dimensional embeddings. We construct finite-dimensional truncations of Bochner and Mercer embedding using random feature [Ren et al., 2022] and Nyström feature [Williams and Seeger, 2000] , respectively, to provide efficient finite linear approximations to represent the transition kernel. As show in Algorithm 1 and Algorithm 2, random feature is the Monte-Carlo approximation for the Bochner embedding and Nyström feature approximates the subspace spanned by the top eigenfunctions of the Mercer embedding via eigendecomposition of an empirical Gram matrix. A detailed derivation of the Nyström method is provided in Appendix B. We analyze the approximation error due to using the finite-dimensional basis in Appendix D. The finite-dimensional embedding is a crucial component of our algorithm, Spectral Dynamic Embedding Policy Optimization (SDEPO), which is given as Algorithm 3.After generating the finite-dimensional embedding, SDEPO is divided in two stages. It finds an approximate solution π k in Stage 1, which is then utilized in Stage 2 to derive an approximate solution π k . In each stage, there are two main components: policy evaluation and policy improvement. Least square policy evaluation is conducted for estimating the state-action value function of current policy pair upon the generated finite-dimensional truncation featuresQ πt,π t (s, a, b) = ϕ(s, a, b) ⊤ w πt,π t .In the policy improvement procedure, based on the approximated value function, the natural policy gradient method is used to adjust the policy of each player iteratively in an alternating fashion. Initializeθ 0 = θ 0 = 0 and π 0 (•|s) = π 0 (•|s) = U nif (A) for all s ∈ S. for k = 0, 1, • • • , K do Stage 1 Initialize θ k,0 = θ k , θ k,0 = θ k and π k,0 (•|s) = π k (•|s), π k,0 (•|s) = π k (•|s) for all s ∈ S. for t = 0, 1, • • • , T -1 do Sample i.i.d. {s i , a i , b i , s ′ i , a ′ i , b ′ i } i∈[n] with policy pair π k-1 , π k-1 , where s ′ i = f (s i , a i , b i ) + ε .Initialize ŵk,t,0 = 0.EQUATION) Output π k = 1 T T t=1 π k,t . Stage 2 Initialize θ ′ k,0 = 0 and π ′ k,0 (•|s) = U nif (A). for t = 0, 1, • • • , T -1 do Sample i.i.d. {s i , a i , b i , s ′ i , a ′ i , b ′ i } i∈[n] with policy pair π ′ k,t , π k ,where s ′ i = f (s i , a i , b i ) + ε.Initialize ζk,t,0 = 0 for all s ∈ S.EQUATIONOutputπ k = π ′ k, t, where t ∈ [T ] is selected uniformly at random. π t+1 is updated as π t+1 (•|s) = arg max π(•|s)∈∆(A) π(•|s), E b∼π t (•|s) [ϕ(s, •, b)] ⊤ w πt,π t + 1 η KL (π(•|s)||π t (•|s))yielding the following closed-form solution,π t+1 (a|s) ∝ π t (a|s) exp E b∼π t (•|s) [ϕ (s, a, b)] ⊤ ηw πt,π t = exp E b∼π t (•|s) [ϕ (s, a, b)] ⊤ θ t+1where θ t+1 = t i=0 ηw πi,π i . π is updated by similar update rule. Note that in SDEPO, each player is required to have knowledge of their opponent's strategy, which is crucial for effective policy optimization and is a common requirement in numerous existing algorithms for Markov games. [Xie et al., 2020 , Chen et al., 2022 , Junchi Li et al., 2022 , Alacaoglu et al., 2022]  
theoretical results: 
 The major difficulty in analyzing the convergence of SDEPO is that the policies of both players evolve at each iteration, leading to non-stationary in the environment for each player. Additionally, there exists a certain level of estimation error in the transition and Q-function for each player. As a consequence, the vanilla proof strategy of convergence of policy optimization used in majority of the literature is no longer applicable. In this section, we provide rigorous investigation of the impact of the approximation error for policy evaluation and derive the convergence of policy optimization. We first specify the commonly used assumptions [Yu and Bertsekas, 2008 , Jin et al., 2020 , Agarwal et al., 2021 , Abbasi-Yadkori et al., 2019 , Ren et al., 2023] , under which we derive our theoretical results below. Assumption 2 (Regularity Condition for Dynamics). For the dynamic function f , there exists aconstant c f , such that ∥f (s, a, b)∥ ≤ c f for all s ∈ S, a ∈ A, b ∈ A.Assumption 3 (Regularity Condition for Stationary Distribution). The stationary distribution ν π,π for all policy pair (π, π) has full support, and satisfies the following conditions with Υ 1 , Υ 2 > 0:λ min E νπ,π ϕ(s, a, b)ϕ(s, a, b) ⊤ ≥ Υ 1 , λ min E νπ,π ϕ(s, a, b) ϕ(s, a, b) -γE νπ,π ϕ(s ′ , a ′ , b ′ ) ⊤ ≥ Υ 2 . Assumption 4 (Regularity Condition for Feature). The features {ϕ ω (s, a, b)} (s,a,b)∈S×A×A and {ϕ ω (s, a, b)} (s,a,b)∈S×A×A are linearly independent.First, we have the following bound on the error for policy evaluation. For convenience, we focus solely on the error for least square policy evaluation in Stage 1. A similar analysis can be conducted for policy evaluation error in Stage 2. We decompose the error for policy evaluation into two parts, one is the approximation error due to the finite number of features, and one is the statistical error due to the finite number of samples. Combine the approximation error and the statistical error, we have the following bound on the error for policy evaluation. Detailed proof can be found in Appendix D.Theorem 1. Let L = Θ(log n). Denote gα := sup s,a,b gα(f (s,a,b)) α dand Qπ,π Φ,L = ϕ ⊤ ŵL . With probability at least 1 -δ, we have that for the random features,EQUATION)and for the Nyström features,EQUATION)As shown in Theorem D.1, Nyström method improves the approximation error from O(m -1 ) to O(n -1 nys ) with a mild assumption to make in the kernel literature (cf. Theorem 2 in [Belkin, 2018] ). Then, we provide the error analysis for policy optimization. We remark that in the following proof we assume the action space |A| is finite for simplicity. The following assumption assumes that the selection probability of each action is positive under each state, which is a commonly used assumption in the analysis of policy gradient type methods [Alacaoglu et al., 2022 , Lan, 2023] . 4Assumption 5. There exists a constant c such that, for any policy iterate pair π k , π k , for any state action tuple s, a, b, it holds that π k (a|s) ≥ c > 0, π k (b|s) ≥ c > 0. We now present the result for policy optimization. We use Perolat et al. [2015] 's error propagation framework, which needs the results in each stage. The detailed proof of it can be found in Appendix E. First, we analyze the error in the Stage 1 of Algorithm 3.Lemma 3. Denote Q k-1 = Q π k-1 ,π k-1 , π s Q s π s = E a∼π(•|s),b∼π(•|s) [Q(s, a, b)]. Let Assump- tion 5 hold. In Stage 1 of Algorithm 3, E max π πQ k-1 π k -min π π k Q k-1 π ≤ 1 ηT log 1 c + 2∥ Qπ k-1 ,π k-1 -Q π k-1 ,π k-1 ∥ νπ k-1 ,π k-1 + η 2(1 -γ) 2 .This result shows how the error in the approximation of the Q-function and the learning rate η impact the optimization performance in Stage 1. Next, we turn to Stage 2 and provide an upper bound on the one-sided error. Lemma 4. Let Assumption 5 hold and µ 0 be a state distribution. In Stage 2 of Algorithm 3,1 T T t=1 E s∼µ0 V π * k ,π k (s) -V π ′ k,t ,π k (s) ≤ 1 1 -γ log 1 c ηT + η (1 -γ) 2 + η∥ Qπ ′ k,t ,π k -Q π ′ k,t ,π k ∥ ν π ′ k,t ,π k .Building on the results from Stage 1 and Stage 2, along with the policy evaluation error, and using the error propagation framework from Perolat et al. [2015] , we can now derive the overall convergence result for the policy optimization process. Specifically, we obtain the following proposition, which quantifies the iteration complexity required to reach a near-optimal solution:Proposition 1. Let L = Θ(log n) and T = log 1/c log n ( 1 (1-γ) 2 + gα (1-γ) 2 √ m + g6 α m 3 (1-γ)Υ 2 1 Υ2 √ n ) for random features and T = log 1/c log n ( 1 (1-γ) 2 + gα (1-γ) 2 nnys + g6 α m 3 (1-γ)Υ 2 1 Υ2 √ n ) for Nyström features. Iteration complexity to get E s∼µ0 max π V π,π k (s) -V * (s) ≤ ϵ is O( 1 (1-γ) 3 ϵ ).Proposition 1 shows the iteration complexity to one-sided NE and it can be directly extended to establish a two-sided NE by applying the algorithm with the roles switched [Zhao et al., 2022] . 
the practical variant of sdepo: 
 Many practical 2p0s-MGs not only have continuous state spaces but also continuous action spaces, such as such as real-time strategy games [Vinyals et al., 2019 , Berner et al., 2019] and robust policy optimization [Pinto et al., 2017] . Note that the natural policy gradient update ( 31) and ( 12 We restrict the players' policies π ψ and π ψ to Gaussians with the reparametrization trick, i.e., a t = f ψ (ϵ t ; s t ) and b t = f ψ (ϵ ′ t ; s t ) where ϵ t and ϵ ′ t are input noise vectors, sampled from a Gaussian. The policy parameters can be learned by minimizingJ(ψ) = E st∼D,ϵt,ϵ ′ t ∼N log π ψ (f ψ (ϵ t ; s t )|s t ) -Q θ (s t , π ψ (f ψ (ϵ t ; s t )|s t ), π ψ (f ψ (ϵ ′ t ; s t )|s t )), andJ(ψ) = E st∼D,ϵt,ϵ ′ t ∼N log π ψ (f ψ (ϵ ′ t ; s t )|s t ) -Q θ (s t , π ψ (f ψ (ϵ t ; s t )|s t ), π ψ (f ψ (ϵ ′ t ; s t )|s t )), where D is the replay buffer, and π ψ and π ψ are defined implicitly in terms of f ψ and f ψ , respectively. Note that tabular methods [Alacaoglu et al., 2022 , Bai and Jin, 2020 , Daskalakis et al., 2020 , Wei et al., 2021 , Zhao et al., 2022] can discretize the state/action spaces to handle applications with continuous state/action spaces. However, directly applying such methods often incurs the curse of dimension and requires an excessive amount of computation resources even for a small size problem. On the other hand, although there exist theoretical-guaranteed methods for 2p0s-MGs with continuous state and action space, they all involve a computational intractable subroutine, i.e., Qiu et al. [2021] , Junchi Li et al. [2022] need to solve a difficult 'find_ne'/'find_cce' subroutine, and Jin et al. [2022] , Huang et al. [2021] have to tackle a comprehensive constrained optimization problem. 
numerical verification: 
 In this section, we present two experiments to evaluate our methods. The first experiment focuses on a simple zero-sum Markov game featuring a continuous state space and a finite action space, aiming to validate the convergence of SDEPO. The second experiment adapts a multi-agent scenario inspired by the simple push [Lowe et al., 2017] , where both the state and action spaces are continuous, to assess the effectiveness of SDEPO-NN. In the first experiment, we designed a simple zero-sum Markov game with a continuous state and finite action space (S = R, |A| = 5). The state space is partitioned into 42 distinct intervals: one interval for (-∞, -10), 40 intervals evenly spaced by 0.5 units in the range [-10, 10), and one interval for (10, ∞). In the i-th interval, the transition dynamics are defined by P (s, a, b) = f (s, a, b) + ϵ, where ϵ ∼ N (0, 1), and f (s, a, b) = ϵ i,a,b , with ϵ i,a,b ∼ Unif(-10.5, 10.5). The reward function is r(s, a, b) = ϵ ′ i,a,b , where ϵ ′ i,a,b ∼ Unif(-1, 1). The initial state distribution is assumed to be uniform over [-10.5, 10.5 ]. We ran SDEPO for 120 iterations, and measured the convergence of π by metrics in Proposition 1. As shown in Figure 1 , SDEPO with random features and Nyström features both converge after 60 iterations. We discretized the state space of this environment and compared it with OFTRL [Zhang et al., 2022] , a tabular method where the environment is known. We adopted the parameter settings recommended in [Zhang et al., 2022] and adjusted the environment to a 100-horizon setting. As shown in Figure 1 , our method demonstrated superior convergence in this environment. This likely stems from the fact that OFTRL operates on the discretized state space, whereas our method computes on the original state space. Next, we conduct experiments on an adapted version of simple push [Lowe et al., 2017] , wherein both the state and action spaces are continuous. This problem consists of two agents and one landmark. Each agent receives a reward for proximity to the landmark while ensuring the other agent remains distant. Thus, agents must learn to stay close to the landmark and simultaneously push the other agent away. At each time step, a noise ϵ ∼ N (0, σ 2 I d ) is added to the state. We implemented SDEPO-NN with random features (SDEPO-NN rf ) and Nyström features (SDEPO-NN nys ), comparing them against methods where Q functions do not utilize spectral dynamical embedding (NPG-NN). 2 : Comparison of winning rates between NPG-NN, SDEPO-NN rf , and SDEPO-NN nys in Simple Push with σ = 0.1/0.01. The results before and after / correspond to σ = 0.1 and 0.01, respectively. 
conclusion: 
 In this paper, we propose a provably efficient natural policy gradient algorithm for two-player zerosum stochastic Markov games with continuous state. We analyze the approximation error and convergence of the algorithm. To deal with continuous action spaces, a practical variant is provided and demonstrates superior performance. A possible direction is to extend our methods to independent learning setting. 
experiments compute resources: 
 Question: For each experiment, does the paper provide sufficient information on the computer resources (type of compute workers, memory, time of execution) needed to reproduce the experiments? Answer: [Yes] Justification: Section 6 Guidelines:• The answer NA means that the paper does not include experiments.• The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage. • The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute. • The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper). 9. Code Of Ethics Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics https://neurips.cc/public/EthicsGuidelines? Answer: [Yes] Justification: Section 6 Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 10. Broader Impacts Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed? Answer: [Yes] Justification: Section 7 Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations. • The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA] Justification: the paper poses no such risks. Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort. 12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected? Answer: [Yes] Justification: We cite the original paper that produced the environment. Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators. 13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? 
 
 Appendix: In this section, we derivate spectral dynamic embedding and then prove the Q-function for arbitrary policy pair can be linearly represented by the feature functions {ϕω (s, a, b)} with ω ∼ N 0, σ 2 I d and ϕ M (s, a, b) ∈ ℓ 2 . Lemma A.1 (i.e. Lemma 1). Consider any α ∈ [0, 1). Denote k α (x, x ′ ) = exp -(1-α 2 )∥x-x ′ ∥ 2 2σ 2for any 0 ≤ α < 1. We can decompose k α (x, x ′ ) using Bochner decomposition and Mercer decomposition. Letψ ω (s, a, b) = g α (f (s, a, b)) α d cos ω ⊤ f (s, a, b) √ 1 -α 2 , sin ω ⊤ f (s, a, b) √ 1 -α 2 , χ ω (s ′ ) = p α (s ′ )[cos( 1 -α 2 ω ⊤ s ′ ), sin( 1 -α 2 ω ⊤ s ′ )] ⊤ ,whereg α (f (s, a, b)) := exp α 2 ∥f (s,a,b)∥ 2 2(1-α 2 )σ 2 , ω ∼ N (0, σ -2 I d ), and p α (s ′ ) = α d (2πσ 2 ) d/2 exp -∥αs ′ ∥ 2 2σ 2is a Gaussian distribution for s ′ with standard deviation σ α . Using Bochner decomposition, we have 15) By Mercer's theorem, k α admits a decomposition for any (x, x ′ ) ∈ X :P (s ′ |s, a, b) =E ω∼N (0,σ -2 I d ) ψ ω (s, a, b) ⊤ χ ω (s ′ ) := ⟨ψ ω (s, a, b), χ ω (s ′ )⟩ N (0,σ -2 I d ) , (EQUATION)We denoteEQUATION)Proof. For Bochner decomposition, we first notice that ∀α ∈ (0, 1), we haveEQUATION)The factorization of the transition P (s ′ |s, a, b) in ( 15) can thus be derived from the property of the Gaussian kernel by applying Bochner decomposition in Lemma 2 to the second term in (18).For the Mercer decomposition, the proof is analogous, where here we apply Mercer's theorem to decompose the middle term in terms of the kernel k αk α s ′ , f (s, a, b) 1 -α 2 = exp   - (1 -α 2 ) s ′ -f (s,a,b) 1-α 2 2 2σ 2    = exp - (1 -α 2 )s ′ -f (s, a, b) 2 2σ 2 (1 -α 2 ) . Let ϕ ω (•) = [ψ ω (•), r(•)] and ϕ M (•) = [ψ M (•), r(•)].Here we prove the Q-function for arbitrary policy pair can be linearly represented by the feature functions {ϕ ω (s, a, b)} with ω ∼ N 0, σ 2 I d and ϕ M (s, a, b) ∈ ℓ 2 .Lemma A.2 (i.e. Lemma 2). For any policy, there exist weights {θ π,π ω } (where ω ∼ N (0, σ -2 I d )) and θ π,π M ∈ ℓ 2 such that the corresponding value function satisfiesQ π,π (s, a, b) = ϕ ω (s, a, b), θ π,π ω N (0,σ -2 I d ) = ϕ M (s, a, b), θ π,π M ℓ2 Proof. Denote µ ω (•) = [χ ω (•), 0] ⊤ , µ M (•) = [1, χ M (•)] ⊤ , θ r = [0, 0, 1] ⊤ , θ r,M = [1, 0, 0, . . . ] ⊤ ∈ ℓ 2 .Our claim can be verified easily by applying the decompositions to Bellman recursion:Q π,π (s, a, b) = r(s, a, b) + γE P V π,π (s ′ ) (19) = ϕ ω (s, a, b), θ r + γ S µ ω (s ′ )V π,π (s ′ )ds ′ θ π,π ω N (0,σ -2 I d ) = ϕ M (s, a, b), θ r,M +γ S µ M (s ′ )V π,π (s ′ )ds ′ θ π,π M ℓ2 ,where the second equation comes from the Bochner decomposition and the third equation comes from the Mercer decomposition.Consider a bounded and continuous positive definite kernel k(x, y) defined on a compact space S, along with any probability measure µ on S, e.g. the µ nys considered in Algorithm 3. Mercer's theorem guarantees the existence of eigenvalues {σ j } ∞ j=1 and orthonormal eigenvectors {e j } ∞ j=1 ⊂ L 2 (µ) such that for any j ∈ N and any x ∈ S, the following holds: S k(x, y)e j (y)dµ(y) = σ j e j (x).(20)The Nyström method provides an approximation of the Mercer eigendecomposition of k(x, y) = ∞ j=1 σ j e j (x)e j (y) in the formkn m (x, y) = m i=1 σi êi (x)ê i (y)for some positive integer m. Here, the pairs where the (σ i , êi ) are determined through a numerical approximation of the eigenfunction problem outlined (20). We will now describe the procedure in detail.Recall that n nys refers to the number of samples used to construct the Nyström features. For simplicity, in the sections of the Appendix that address Nyström approximation results, we denote n := n nys , unless otherwise noted. It is important not to confuse this n with the n used to represent the number of samples for statistical learning of the Q-function. Suppose we draw n independent samples X n = {x s } n s=1 from S according to the distribution µ. For any eigenfunction e j , the eigenfunction problem from (20) can be numerically approximated by:EQUATION)EQUATIONLet K (n) ∈ R n×n denote the Gram matrix where (K n ) rs = k(x r , x s ), and let the eigendecomposition of K (n) be K (n) U = ΛU with U orthogonal and Λ diagonal. To satisfy (21a), the condition must hold for all x ∈ X n . Specifically, for any r ∈ [n], we should expectEQUATIONfor any eigenfunction e j (•) and its corresponding eigenvalue σ j of k(•, •). For any eigenvector u i of the Gram matrix K (n) , we have the relation1 n n s=1 k(x r , x s )(u i ) s = λ i n (u i ) r .It is therefore natural to extend u i ∈ R n to an eigenfunction ẽi by setting ẽi (x s ) = (u i ) s for any x s ∈ X n , and for other x, defineẽi (x) := 1 λ i n s=1 k(x, x s )(u i ) s ,with associated eigenvalue λi n . To satisfy the orthonormality condition from (21b), since1 n n s=1 ẽi (x s ) 2 = 1 n n s=1 (u i ) 2 s = 1 n , we scale êi (x) := √ nẽ i (x) = √ n λ i n s=1 k(x, x s )(u i ) s .Using these scaled eigenfunctions, we define the Nyström approximation kn of the original kernel k asEQUATION)where k(x, X n ) is a row vector with components k(x, x s ), and k(X n , y) is a column vector with components k(x s , y). Based on (23), for any m ≤ n, we then define the rank-m Nyström approximation asEQUATION)which is the rank-m Nyström kernel approximation used in the paper. Since it can be written as kn ,m (x, y) = φ nys (x) ⊤ φ nys (y), whereEQUATIONwe can view φ nys (•) ∈ R m as the rank-m Nyström features corresponding to the Nyström approximation.Lemma C.3 (Performance difference lemma. See [Alacaoglu et al., 2022] ). For any policies π 1 , π 2 , π and any state s 0 ,V π1,π (s 0 ) -V π2,π (s 0 ) = 1 1 -γ E s ′ ∼d π 1 ,π s 0 ⟨E b∼π(•|s) Q π2,π (s, •, b), π 1 (•|s) -π 2 (•|s)⟩ where d π1,π2 s0 (s) = (1 -γ)∞ t=0 γ t P π1,π2 (s t = s|s 0 ), where P π1,π2 (s t = s|s 0 ) denote the probability that s t = s after starting at s 0 and following the policies π 1 and π 2 .We provide a brief outline of our overall proof strategy for policy evaluation. For convenience, we focus solely on the error for least square policy evaluation in Stage 1. A similar analysis can be conducted for policy evaluation error in Stage 2.We decompose the error into two parts, one is the approximation error due to the limitation of our basis (i.e., finite m in Line 1 of Algorithm 1 and Line 4 of Algorithm 2), and one is the statistical error due to the finite number of samples we use (i.e., finite n in Line 27 of Algorithm 3). For notational simplicity, we omit π, π and use ν to denote the stationary distribution corresponding to π, π in this section. For the ease of presentation, we omit the polynomial dependency on c f and focus on the dependency of other terms of interest.We start by deriving a bound on the approximation error when representing the Q-function for a given policy pair using an imperfect, finite-dimensional feature set. In such cases, the best possible approximation, denoted as Q, is the solution to a projected Bellman equation (cf. Yu and Bertsekas [2008] ), which is defined as follows: given any (possibly finite-dimensional) feature map Φ := {ϕ(s, a, b)} (s,a,b)∈S×A×A , the approximation QΦ π,π is defined by: QΦ π,π = Πν, Φ(r + P π,π QΦ π,π ),where ν represents the stationary distribution under π, π, and the operator P π,π is given by(P π,π f )(s, a, b) = E (s ′ ,a ′ ,b ′ )∼P (s,a,b)×π×π f (s ′ , a ′ , b ′ ).(27) and Π ν,Φ is the projection operator as defined asEQUATION)Our focus on Qπ Φ is motivated by the fact that the least-squares policy evaluation step in Algorithm 3 (see equation ( 11)) recovers Qπ Φ if the number of samples, n, goes to infinity. The effect of a finite sample size n on the statistical error will be discussed later.Next, we present our bound on the approximation error for the random and Nyström features. We begin with the random feature approach, which . To do so, we first need the following technical result, which relies on the following technical lemma, adapted from Lemma 1 in Rahimi and Recht [2008] . Lemma D.4 (cf. Lemma 1 from Rahimi and Recht [2008] ). Let p be a distribution on a space Ω, and consider a mapping ϕ(x; ω) ∈ R ℓ . Supposef * (x) = Ω p(ω)β(ω) ⊤ ϕ(x; ω)dω,for some vector β(ω) ∈ R ℓ where sup x,ω β(ω) ⊤ ϕ(x; ω) ≤ C for some C > 0. Consider {ω i } k i=1 drawn iid from p, and denote the sample average of f * as f(x) = 1 K K k=1 β(ω k ) ⊤ ϕ(x; ω k ).Then, for any δ > 0, with probability at least 1 -δ over the random draws of{ω i } k i=1 , X f (x) -f * (x) 2 dµ(x) ≤ C √ K 1 + 2 log 1 δ .Based on Lemma D.4, we derive the approximation error with random features. Proposition D.1 (Q-Approximation error with random features). We define the feature map Φ rf for random features as followsEQUATION)where ψ rf (s, a, b) is defined in (7) in Algorithm 3. Then, for any δ > 0, with probability at least 1 -δ, we have thatEQUATION)where ∥ • ∥ ν is the L 2 norm defined as ∥f ∥ ν = f 2 dν, and Qπ,π Φ rf is defined in (26).Proof. By leveraging the contraction property and the results from Yu and Bertsekas [2008] , we establish the following bound:EQUATION)where Π ν,Φ rf is defined as in (28). By definition Π ν,Φ rf is contractive under ∥ • ∥ ν . Note thatQ π,π (s, a, b) = r(s, a, b) + γE ω∼N (0,σ -2 I d ) ϕ ω (s, a, b) ⊤ S µ ω (s ′ )V π,π (s ′ )ds ′ .Using Hölder's inequality, along with the fact that ∥V π ∥ ∞ = O((1 -γ) -1 ), and denoting β π,π ω:= S µ ω (s ′ )V π,π (s ′ )ds ′ , we have for every ω that β π,π ω ∞ = O((1 -γ) -1 ). Additionally, recalling that ϕ ω (s, a, b) = [ψ ω (s, a, b), r(s, a, b)], we have sup ((s,a,b),ω) |ϕ ω (s, a, b) ⊤ β π,π ω | ≤ sup ((s,a,b),ω) |ψ ω (s, a, b)| 1 β π,π ω ∞ since the coordinate in β π,π ω corresponding to the reward part of ϕ ω (s, a, b) is 0. Noting sup ((s,a,b),ω) |ψ ω (s, a, b)| 1 ≤ 2 sup (s,a,b) g α (f (s, a, b)) α d , we conclude that sup ((s,a,b),ω) |ϕ ω (s, a, b) ⊤ β π,π ω | ≤ sup ((s,a,b),ω) |ψ ω (s, a, b)| 1 β π,π ω ∞ = O sup (s,a,b) g α (f (s, a, b)) α d (1 -γ)Applying Lemma D.4, we obtain the following bound:EQUATION)Substituting ( 32) into (31) completes the proof.As shown above, the approximation error for random features decreases at a rate of O m -1/2 with high probability as the number of random features m increases.Next, we turn to the approximation error for Nyström features. To establish this result, we first present the following key finding. Lemma D.5 (cf. Lemma 8 from Ren et al. [2023] ). Consider the following Mercer decomposition (on S) of k α (•, •):EQUATION)where {e i } ∞ i=1 forms a countable orthonormal basis for L 2 (µ nys ) with corresponding values{σ i } ∞ i=1 . Let X nnys = {x i } nnys i=1 be an i.i.d n nys -point sample from µ nys . In addition, let λ 1 ≥ λ 2 ≥ • • • ≥ λ nnys denotethe eigenvalues of the (unnormalized) Gram matrix K (nnys) in its eigendecomposition K (nnys) U = ΛU where U ⊤ U = U ⊤ = I. Suppose that σ j , λ j /n nys ≲ exp(-βj 1/h ) for some β > 0 and h > 0. Suppose that n nys ≥ 3 and that ⌊(2log n nys ) h /β h ⌋ ≤ m ≤ n nys .Consider the rank-m kernel approximation knnys α,m constructed using Nyström features, defined as follows:EQUATION)where φ nys (•) ∈ R m and is defined asEQUATION)where u i denotes the i-th column of U , and k α (X nnys , •) denotes a n nys -dimensional vector where(k α (X nnys , •)) ℓ = k α (x ℓ , •); for details on how φ nys (•) is derived, see Appendix B.Then, for any δ > 0, with probability at least 1 -δ,EQUATION)Leveraging Lemma D.5, we can now derive the following result for the approximation error associated with Nyström features.Proposition D.2 (Q-Approximation error with Nyström features). Suppose all the assumptions in Lemma D.5 hold. Suppose also that we pick the sampling distribution µ nys such thatEQUATION)We define the feature map Φ nys for the Nyström features as follows:EQUATION)where ψ nys (s, a, b) ∈ R m is defined in (38) in Algorithm 3. Then, for any δ > 0, with probability at least 1 -δ,Q π,π -Qπ,π nys ν ≤ Õ γg α (1 -γ) 2 n nys ≤ Õ γg α (1 -γ) 2 m , where ∥ • ∥ ν is the L 2 norm defined as ∥f ∥ ν = f 2 dν, Qπ,πΦnys is defined in (26), and gα := sup s,a,b gα(f (s,a,b)) α d .Proof. Following the analysis at the beginning of the proof of Proposition D.1, and using the contraction property along with results from Yu and Bertsekas [2008] , we have the following bound:EQUATION)where Π ν,Φnys is defined as in (28).Next, we need to bound the term Q π,π -Π ν,Φnys Q π,π ν . Starting from (34), and recalling the definition of φ nys (•) from ( 35), we express our Nyström approximation of the kernelEQUATION)Since we haveψ nys (s, a, b) = g α (f (s, a, b))φ nys f (s,a,b) 1-α 2, where ψ nys (•) is defined as in (8), this suggests the following Nyström-based approximation for the Q-function:Qπ,π nys (s, a, b) := r(s, a, b) + γψ nys (s, a, b) ⊤ S φ nys (s ′ )V π,π (s ′ )p α (s ′ )ds ′ α d .It is important to note that this approximation is a valid solution to the objective:arg min f ∈span(Φnys) E ν Q π,π (s, a, b) -f (s, a, b) 2 .Thus, we haveEQUATION)Using the one-to-one correspondence between Reproducing Kernel Hilbert Spaces (RKHS) and PD kernels, we can express the kernel function k α as an inner product in the corresponding RKHS, H kα . Specifically, we havek α s ′ , f (s, a, b) 1 -α 2 = ⟨k α f (s, a, b) 1 -α 2 , • , k α (s ′ , •)⟩ H kα knnys m s ′ , f (s, a, b) 1 -α 2 = k α ( f (s, a, b) 1 -α 2 , X n ) m i=1 1 λ i u i u ⊤ i k α (X n , s ′ ) = k α ( f (s, a, b) 1 -α 2 , •), k α (X n , •) H kα m i=1 1 λ i u i u ⊤ i k α (X n , s ′ ) = ⟨k α f (s, a, b) 1 -α 2 , • , knnys m (s ′ , )⟩ H kαBy choosing µ nys = p α (x), continuing from (41), we haveT 1 = S knnys m s ′ , f (s, a, b) 1 -α 2 -k α s ′ , f (s, a, b) 1 -α 2 p α (s ′ )ds ′ ν = S ⟨k α f (s, a, b) 1 -α 2 , • , knnys m (s ′ , )⟩ H kα -⟨k α f (s, a, b) 1 -α 2 , • , k α (s ′ , •)⟩ H kα p α (s ′ )ds ′ ν = S k α f (s, a, b) 1-α 2 , • , knnys m (s ′ , •)-k α (s ′ , •) H kα dp α (s ′ ) ν ≤ S k α f (s, a, b) 1 -α 2 , f (s, a, b) 1-α 2 (k α - knnys m )(s ′ , s ′ )dp α (s ′ ) ν ≤ S (k α - knnys m )(s ′ , s ′ )dp α (s ′ ) νTo move from the second-to-last line to the last line, we used the fact that k α (•, •) ≤ 1. Next, applying Lemma D.5 and using the decay assumption on the eigenvalues from both the Mercer expansion and the empirical Gram matrix, we obtain a probabilistic bound. Specifically, for any δ > 0, with probability at least 1 -δ, we haveT 1 ≤ S (k α - knnys m )(s ′ , s ′ )µ(s ′ )ds ′ = Õ 1 n nys .Thus, by combining ( 41) with the bound on T 1 , we conclude the desired result, demonstrating that the error introduced by the Nyström approximation decays at the rate Õ(1/n nys ).Remark D.1. As demonstrated by the proposition above, a key advantage of the Nyström method is its ability to reduce the approximation error to O( n -1 nys ), where n nys represents the number of samples used to construct the Gram matrix for generating Nyström features. This rate of O( n -1 nys ) consistently outperforms the O(1/m) where m is the number of features, since SDEPO is designed to select m ≤ n nys . The only requirement for this improvement is that the eigenvalues of the empirical Gram matrix and Mercer expansion meet certain decay assumptions, which is a standard assumption to make in the kernel literature (cf. Theorem 2 in Belkin [2018] ).We now present the bound on the statistical error arising from using a finite number of samples, n. The result holds for both Nyström and random features.Proposition D.3. For each policy pair {π, π} encountered in the algorithm, let Qπ,π Φ,L denote the policy given byQπ,π Φ,L (s, a, b) = ϕ(s, a, b) ⊤ ŵL, where ŵL is defined as in (11), and Φ can either be the reward concatenated with the Nyström or random features, i.e. either Φ nys or Φ rf . Then, for sufficiently large n, there exists an universal constant C > 0 independent of m, n, L and (1 -γ) -1 , such that with probability at least 1 -δ, we haveEQUATION)where we recall Qπ,π Φ is defined as in (26).Proof. To simplify the presentation, let us define Φ as the concatenation of ϕ(s, a, b) across all (s, a, b) ∈ S × A × A, and define the operator P π,π as(P π,π f )(s, a, b) = E (s ′ ,a ′ )∼P (s,a,b)×π×π f (s ′ , a ′ , b ′ ).Additionally, define w as the solution to the equation:w = E ν ϕ(s, a, b)ϕ(s, a, b) ⊤ -1 E ν ϕ(s, a, b) r(s, a, b) + γE (s ′ ,a ′ ,b ′ )∼P (s,a,b)×π×π ϕ(s ′ , a ′ , b ′ ) ⊤ w , and let Q(s, a, b) = ϕ(s, a, b) ⊤ w.It is straightforward to observe that w is the fixed point of the population (i.e., n → ∞) projected least square update (11). Furthermore, we can express w asw = E ν ϕ(s, a, b) ϕ(s, a, b) -γE (s ′ ,a ′ ,b ′ )∼P (s,a,b)×π×π ϕ(s ′ , a ′ , b ′ ) ⊤ -1 E ν [ϕ(s, a, b)r(s, a, b)] .For the sake of brevity, we will omit the subscript when ν refers to the Lebesgue measure. We define Πν and P π,π as the empirical counterparts of Π ν and P π,π , respectively. Under the update rule (11), we have the following relations:Φ ŵt+1 = Πν (r + γ P π,π Φ ŵt ), Φ w =Π ν (r + γP π,π Φ w), which lead to Φ( w -ŵt+1 ) = (Π ν -Πν )r + γ(Π ν P π,π )Φ( w -ŵt ) + γ(Π ν P π,π -Πν P π,π )Φ ŵt .Applying the triangle inequality, we can bound this as:∥Φ( w -ŵt+1 )∥ ν ≤ γ ∥Φ( w -ŵt )∥ ν + Π ν -Πν r ν + γ Π ν P π,π -Πν P π,π Φ ŵt ν ,where we use the contractivity under ∥ • ∥ ν . Telescoping over t, we have∥Φ( w -ŵL )∥ ν ≤γ T ∥Φ( w -ŵ0 )∥ ν + 1 1 -γ Π ν -Πν r ν + γ 1 -γ max l∈[L] Π ν P π,π -Πν P π,π Φ ŵt ν .The following proof follows the same reasoning as Appendix E in Ren et al. [2023] .Proposition D.3 shows that the statistical error of the linear components can be broken down into two parts. The first part, γ L Qπ,π ν , arises from initializing with ŵ0 = 0 and decreases as the number of least-squares policy evaluation iterations L increases. The second part,C g6 α m 3 polylog(m,L/δ) (1-γ)Υ 2 1 Υ2√ n , represents the statistical error from using a finite sample size and diminishes as n increases. By choosing L = Θ(log n), we can balance both parts, resulting in an estimation error that decreases at a rate of O(n -1/2 ) with high probability.Combine the approximation error in Proposition D.1 and Proposition D.2 with the statistical error from Proposition D.3, we obtain the following bound on the total error for least-squares policy evaluation: Theorem D.1 (i.e. Theorem D.1). Let L = Θ(log n). With probability at least 1 -δ, we have that for the random features Φ rf ,EQUATION)and for the Nyström featuresQ π,π - Qπ,π Φnys,L ν = Õ      gα (1 -γ) 2 n nys approx. error + g6 α m 3 (1 -γ)Υ 2 1 Υ 2 √ n stat. error      .Proof. This result follows directly from the triangle inequality, i.e.,Q π,π - Qπ,π Φ,L ν ≤ Q π,π - Qπ,π Φ ν + Qπ,π Φ - Qπ,π Φ,L ν .The first term, representing the approximation error, is bounded by Proposition D.1 and Proposition D.2 for random and Nyström features, respectively. The second term, capturing the statistical error, is bounded by Proposition D.3.Theorem D.1 provides the estimation error bound for the Q-function using least-squares policy evaluation, applicable to both random and Nyström features. The result reveals a key tradeoff between the approximation error and statistical error. For random features, increasing the number of features m enhances the ability to approximate the original infinite-dimensional function space, as reflected by the Õ 1 √ m approximation error term on the RHS of (43). However, this improvement comes at the cost of needing more learning samples n during policy evaluation to effectively train the weights, as indicated by the Õ m 3 √ n statistical error term in (43). For Nyström features, this tradeoff may be better especially since the approximation term scales as Õ(1/n nys ) ≤ Õ(1/m) while thestatistical error term remains the same as with random features. As a result, the Nyström method offers a better approximation error rate with additional mild assumption.In this section, we denoteQ k-1 = Q π k-1 ,π k-1 , and use the notation π s Q s π s = E a∼π(•|s),b∼π(•|s) [Q(s, a, b)].To establish the convergence of the policy optimization procedure, we follow the error propagation framework from Perolat et al. [2015] , which consists of two key stages:Stage 1: Identify an approximate solution π k such thatE max π π s Q s k-1 π s k -min π max π π s Q s k-1 π s = ϵ k 1 (s),where the expectation is taken over the randomness of the algorithm used to generate π k .In our analysis, we bound the stronger quantity, which is called duality gapE max π π s Q s k-1 π s k -min π π s k Q s k-1 π s ≥ ϵ k 1 (s),where the bound follows from the definition of Nash equilibrium, asmin π π k Q k-1 π ≤ min π max π πQ k-1 π.Stage 2: Identify an approximate solution π k such thatV π k ,π * k (s) -EV π k ,π k (s) = ϵ k 2 (s), where the expectation is over the randomness of the algorithm used to generate π k .Following the analysis in Perolat et al. [2015] , we conclude that there exists a constant C, such thatEQUATION)First, we analyze the error in the Stage 1 of Algorithm 3. Lemma E.6 (i.e. Lemma 3). Let Assumption 5 hold. In Stage 1 of Algorithm 3, , a, b) ]. First by definition of π k and the standard formulation for the duality gap, we have for all sE max π πQ k-1 π k -min π π k Q k-1 π ≤ 1 ηT log 1 c + 2∥ Qπ k-1 ,π k-1 -Q π k-1 ,π k-1 ∥ νπ k-1 ,π k-1 + η 2(1 -γ) 2 . Proof. Denote Q = Q k-1 . Recall the notation π s Q s π s = E a∼π(•|s),b∼π(•|s) [Q(sπ s Q s π k -π s k Q s π s = 1 T T t=1 ⟨E b∼π k,t (•|s) Q(s, a, •), π(•|s)⟩- 1 T T t=1 ⟨E a∼π k,t (•|s) Q(s, •, b), π(•|s)⟩ = 1 T T t=1 ⟨E b∼π k,t (•|s) Q(s, a, •), π(•|s)-π k,t (•|s)⟩ -⟨E a∼π k,t (•|s) Q(s, •, b), π(•|s)-π k,t (•|s)⟩.From the update rule of π k,t+1 , it holds for all k, s, π that⟨∇KL(π k,t+1 , π k,t ) + ηE a∼π k,t [ Qπ k-1 ,π k-1 (s, a, •)], π(•|s) -π k,t+1 (•|s)⟩ ≥ 0.Applying the three-point identity yields:KL(π(•|s), π k,t+1 ( • |s)) ≤ KL(π(•|s), π k,t (•|s)) -KL(π k,t+1 (•|s), π k,t (•|s)) +η⟨E a∼π k-1 [ Qπ k-1 ,π k-1 (s, a, •) -Q π k-1 ,π k-1 (s, a, •)], π(•|s) -π k,t+1 (•|s)⟩ +η⟨E a∼π k-1 [Q k-1 (s, a, •)], π(•|s) -π k,t+1 (•|s)⟩.We can bound the inner products using Cauchy-Schwarz, Young's and Pinsker's inequalities,η⟨E a∼π k-1 [ Qπ k-1 ,π k-1 (s, a, •) -Q π k-1 ,π k-1 (s, a, •)], π(•|s) -π k,t+1 (•|s)⟩ ≤ 2η∥ Qπ k-1 ,π k-1 -Q π k-1 ,π k-1 ∥ νπ k-1 ,π k-1 η⟨E a∼π k-1 [Q π k-1 ,π k-1 (s, a, •)], π(•|s) -π k,t+1 (•|s)⟩ =η⟨E a∼π k-1 [Q π k-1 ,π k-1 (s, a, •)], π(•|s) -π k,t (•|s)⟩ +η⟨E a∼π k-1 [Q π k-1 ,π k-1 (s, a, •)], π k,t (•|s) -π k,t+1 (•|s)⟩ ≤η⟨E a∼π k-1 [Q π k-1 ,π k-1 (s, a, •)], π(•|s) -π k,t (•|s)⟩ + η 2 2(1 -γ) 2 + KL(π k,t+1 (•|s), π k,t (•|s))Combining these estimates gives us:⟨E a∼π k-1 [Q π k-1 ,π k-1 (s, a, •)], π k,t (•|s) -π(•|s)⟩ + 1 η KL(π(•|s), π k,t+1 (•|s)) ≤ 1 η KL(π(•|s), π k,t (•|s)) + 2∥ Qπ k-1 ,π k-1 -Q π k-1 ,π k-1 ∥ νπ k-1 ,π k-1 + η 2 2(1 -γ) 2Summing these inequalities results in:1 T T t=1 ⟨E a∼π k,t [Q π k,t ,π k,t (s, a, •)], π k,t (•|s) -π(•|s)⟩ ≤ 1 ηT KL(π, π k,0 ) + 2∥ Qπ k-1 ,π k-1 -Q π k-1 ,π k-1 ∥ νπ k-1 ,π k-1 + η 2 2(1 -γ) 2 ≤ 1 ηT log 1 c + 2∥ Qπ k-1 ,π k-1 -Q π k-1 ,π k-1 ∥ νπ k-1 ,π k-1 + η 2(1 -γ) 2This result shows how the error in the approximation of the Q-function and the learning rate η impact the optimization performance in Stage 1. Next, we turn to Stage 2 and provide an upper bound on the one-sided error. Lemma E.7 (i.e. Lemma 4). Let Assumption 5 hold and µ 0 be a state distribution. In Stage 2 of Algorithm 3,1 T T t=1 E s∼µ0 V π * k ,π k (s) -V π ′ k,t ,π k (s) ≤ 1 1 -γ log 1 c ηT + η (1 -γ) 2 + η∥ Qπ ′ k,t ,π k -Q π ′ k,t ,π k ∥ ν π ′ k,t ,π k .Proof. By the update rule for π ′ k,t+1 , for all s and π, we have thatKL(π(•|s), π ′ k,t+1 (•|s)) ≤ KL(π(•|s), π ′ k,t (•|s)) -η⟨E b∼π k [ Qπ ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t+1 (•|s)⟩ -KL(π ′ k,t+1 (•|s), π ′ k,t (•|s)) = KL(π(•|s), π ′ k,t (•|s)) -η⟨E b∼π k [ Qπ ′ k,t ,π k (s, •, b) -Q π ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t+1 (•|s)⟩ -η⟨E b∼π k [Q π ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t+1 (•|s)⟩ -KL(π ′ k,t+1 (•|s), π ′ k,t (•|s))We can bound the inner products using Cauchy-Schwarz, Young's and Pinsker's inequalitiesη⟨E b∼π k [ Qπ ′ k,t ,π k (s, •, b) -Q π ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t+1 (•|s)⟩ ≤ O(η∥ Qπ ′ k,t ,π k -Q π ′ k,t ,π k ∥ ν π ′ k,t ,π k) andη⟨E b∼π k [Q π ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t+1 (•|s)⟩ = -η⟨E b∼π k [Q π ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t (•|s)⟩ -η⟨E b∼π k [Q π ′ k,t ,π k (s, •, b)], π ′ k,t (•|s) -π ′ k,t+1 (•|s)⟩ ≤ -η⟨E b∼π k [Q π ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t (•|s)⟩ + η 2 ∥Q π ′ k,t ,π k (s, a, b)∥ 2 ∞ 2 + KL(π ′ k,t+1 (•|s), π ′ k,t (•|s))Consequently, we have that⟨E b∼π k [Q π ′ k,t ,π k (s, •, b)], π(•|s)-π ′ k,t (•|s)⟩ + 1 η KL(π(•|s), π ′ k,t+1 (•|s))| ≤ 1 η KL(π, π ′ k,t ) + O(η(∥Q π ′ k,t ,π k (s, a, b)∥ 2 ∞ + ∥ Qπ ′ k,t ,π k -Q π ′ k,t ,π k ∥ ν π ′ k,t ,π k))Summing the inequality, we get1 T T t=0 ⟨E b∼π k [Q π ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t (•|s)⟩ (45) ≤ 1 ηT KL(π, π ′ k,t ) + η (1 -γ) 2 + η∥ Qπ ′ k,t ,π k -Q π ′ k,t ,π k ∥ ν π ′ k,t ,π k (46)By the performance difference lemma, we obtainV π * k ,π k (s 0 ) -V π ′ k,t ,π k (s 0 ) = 1 1 -γ E s∼d π * k ,π k s 0 ⟨E b∼π k (•|s) Q π ′ k,t ,π k (s, •, b), π * k (•|s) -π ′ k,t (•|s)⟩.Proposition E.4. Let Assumption 5 hold and µ 0 be a state distribution,E s∼µ0 max π V π,π K (s) -V * (s) ≤O CK (1 -γ) 2 ( log 1 c ηT + η (1 -γ) 2 + η∥ Qπ k-1 ,π k-1 -Q π k-1 ,π k-1 ∥ νπ k-1 ,π k-1 + η∥ Qπ ′ k,t ,π k -Q π ′ k,t ,π k ∥ ν π ′ k,t ,π k ) + γ K 1 -γ .where C is a problem-dependent constant.Proof. Inserting the results of Lemma E.6,Lemma E.7 to (44) gives the result.Let L = Θ(log n) and T = log 1/c log n ( 1 (1-γ) 2 + gα (1-γ) 2 √ m + g6 α m 3 (1-γ)Υ 2 1 Υ2 √ n ) for random fea- tures and T = log 1/c log n ( 1 (1-γ) 2 + gα (1-γ) 2 nnys + g6 α m 3 (1-γ)Υ 2 1 Υ2√ n ) for Nyström features. Combining the error of policy evaluation in Theorem D.1, we have that the iteration complexity to getE s∼µ0 max π V π,π k (s) -V * (s) ≤ ϵ is O( 1 (1-γ) 3 ϵ ).In this section, we are going to use β-greedy policy to avoid Assumption 5. Let us define the modified policies with greedy explorationπ = (1 -β)π + βU nif (A), π = (1 -β)π + βU nif (A)Now we are going to sample with the π, π and the algorithm will read as Algorithm 4.Lemma E.8. In Stage 1 of Algorithm 4,E max π πQ k-1 π k -min π π k Q k-1 π ≤ η 2(1 -γ) 2 + 32γβ (1 -γ) 2 + 1 ηT log 1 |A| + 2η∥ Qπ k-1 ,π k-1 -Q∥ ν πk-1 ,π k-1 . Proof. Denote Q = Q k-1 , Q = Q πk-1 ,π k-1 .First by definition, it holds for all sπ s Q s π s k -π s k Q s π s = 1 T T t=1 ⟨E b∼π k,t (•|s) Q(s, •, b), π(•|s) -π k,t (•|s)⟩ -⟨E a∼π k,t (•|s) Q(s, a, •), π(•|s) -π k,t (•|s)⟩ = 1 T T t=1 ⟨E b∼π k,t (•|s) Q(s, •, b), π(•|s) -π k,t (•|s)⟩ -⟨E a∼π k,t (•|s) Q(s, a, •), π(•|s) -π k,t (•|s)⟩ + 1 T T t=1 ⟨E b∼π k,t (•|s) [Q(s, •, b) -Q(s, •, b)], π(•|s) -π k,t (•|s)⟩ -⟨E a∼π k,t (•|s) [Q(s, a, •) -Q(s, a, •)], π(•|s) -π k,t (•|s)⟩For the error terms note⟨E b∼π k,t (•|s) [Q(s, •, b) -Q(s, •, b)], π(•|s) -π t (•|s)⟩ ≤ 2∥E b∼π k,t (•|s) Q(s, •, b) -Q(s, •, b)∥ ∞ ≤ 2γ max a,b |E s ′ ∼P (•|s,a,b) [V π k-1 ,π k-1 (s ′ ) -V πk-1 ,π k-1 (s ′ )]| ≤ 2γ∥V π k-1 ,π k-1 -V πk-1 ,π k-1 ∥ ∞ ≤ 16γβ (1 -γ) 2 ,where the last step is due to the Lipschitzness of the value function due to performance difference lemma and that the policies πk-1 , πk-1 and π k-1 , π k-1 differ at most by β.By three point identity, we haveKL(π(•|s), π k,t+1 (•|s)) ≤ KL(π(•|s), π k,t (•|s)) + η⟨E a∼ πk-1 [ Qπ k-1 ,π k-1 (s, a, •) -Q(s, a, •)], π(•|s) -π k,t+1 (•|s)⟩ + η⟨E a∼ πk,t [ Q(s, a, •)], π(•|s) -π k,t+1 (•|s)⟩ -KL(π k,t+1 (•|s), π k,t (•|s)).Again we use Cauchy-Schwarz and Youngss inequalitiesEQUATIONEQUATIONOutputπ k = 1 T T t=1 π k,t . Stage 2: Initialize θ ′ k,0 = 0 and π ′ k,0 (•|s) = U nif (A). for t = 0, 1, • • • , T -1 do Sample i.i.d. {(s i , a i , b i , s ′ i ), a ′ i , b ′ i } i∈[n] with policy pair πt k , πk , where s ′ i = f (s i , a i , b i ) + ε.Initialize ζk,t,0 = 0.EQUATION) Output π k = π ′ k,t, where t ∈ [T ] is selected uniformly at random.η⟨E a∼ πk,t [ Q(s, a, •)], π(•|s) -π k,t+1 (•|s)⟩ =η⟨E a∼ πk,t [ Q(s, a, •)], π(•|s) -π k,t (•|s)⟩ +η⟨E a∼ πk,t [ Q(s, a, •)], π k,t (•|s) -π k,t+1 (•|s)⟩ ≤η⟨E a∼ πk,t [ Q(s, a, •)], π(•|s) -π k,t (•|s)⟩ + η 2 2(1 -γ) 2 + KL(π k+1 (•|s), π k (•|s))Using these estimations, we have that⟨E a∼ πk,t [ Qπ k,t ,π k,t (s, a, •)], π k,t (•|s) -π(•|s)⟩ + 1 η KL(π(•|s), π k,t+1 (•|s)) ≤ 1 η KL(π(•|s), π k,t (•|s)) + η 2(1 -γ) 2 + 2η∥ Qπ k-1 ,π k-1 -Q∥ ν πk-1 ,π k-1Summing the inequality, we get1 T T t=1 ⟨E a∼ πk,t [ Qπ k,t ,π k,t (s, a, •)], π k,t (•|s) -π(•|s)⟩ ≤ 1 ηT KL(π, π k,0 ) + 2η∥ Qπ k-1 ,π k-1 -Q∥ ν πk-1 ,π k-1 ≤ 1 ηT log 1 |A| + η 2(1 -γ) 2 + 2η∥ Qπ k-1 ,π k-1 -Q∥ ν πk-1 ,π k-1 Thus E max π πQ k-1 π k -min π π k Q k-1 π ≤ η 2(1 -γ) 2 + 32γβ (1 -γ) 2 + 1 ηT log 1 |A| + 2η∥ Qπ k-1 ,π k-1 -Q∥ ν πk-1 ,π k-1 Lemma E.9. In Stage 2 of Algorithm 4, E 1 T T t=0 V π * k ,π k (s 0 ) -V π ′ k,t ,π k (s 0 ) ≤O    log |A| (1 -γ)ηT + η (1 -γ) 3 + γβ η(1 -γ) 3 + ∥ Qπ ′ k,t ,π k -Q∥ ν π ′ k,t ,π k 1 -γ    . Proof. Denote Q = Q π′ k,t ,π k .By the update rule of the algorithm,KL(π(•|s), π ′ k,t+1 (•|s)) ≤ KL(π(•|s), π ′ k,t (•|s)) -η⟨E b∼π k [ Qπ ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t+1 (•|s)⟩ -KL(π ′ k,t+1 (•|s), π ′ k,t (•|s)) = KL(π(•|s), π ′ k,t+1 (•|s)) -KL(π ′ k,t+1 (•|s), π ′ k,t (•|s)) -η⟨E b∼π k [ Qπ ′ k,t ,π k (s, •, b) -Q(s, •, b)], π(•|s) -π ′ k,t+1 (•|s)⟩ -η⟨E b∼π k [ Q(s, •, b)], π ′ k,t (•|s) -π ′ k,t+1 (•|s)⟩ -η⟨E b∼π k [Q π′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t (•|s)⟩ -η⟨E b∼π k [ Q(s, •, b) -Q π′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t (•|s)⟩We bound the inner products similarly η⟨E b∼π k [ Qπ ′ k,t ,π k (s, •, b) -Q(s, •, b)], π(•|s) -π ′ k,t+1 (•|s)⟩ ≤ O(η∥ Qπ ′ k,t ,π k -Q∥ ν π′ k,t ,π k ) η⟨E b∼π k [ Q(s, •, b)], π ′ k,t (•|s) -π ′ k,t+1 (•|s)⟩ ≤ η 2 2(1 -γ) 2 + KL(π ′ k,t+1 (•|s), π ′ k,t (•|s)) and η⟨E b∼π k [ Q(s, •, b) -Q π ′ k,t ,π k (s, •, b)], π(•|s) -π ′ k,t (•|s)⟩ ≤ 16γβ (1 -γ) 2{(s i , a i , b i , s ′ i , a ′ i , b ′ i )} to replay buffer.for each gradient step do Update θ,θ by minimizing eqs. ( 55) and ( 57), respectively.Update ψ, ψ by minimizing eqs. ( 61) and ( 62), respectively.Consequently, we have thatEQUATIONBy the performance difference lemmaV π * k ,π k (s 0 ) -V π ′ k,t ,π k (s 0 ) = 1 1 -γ E s∼d π * k ,π k s 0 ⟨E b∼π k (•|s) Q π ′ k,t ,π k (s, •, b), π * k (•|s) -π ′ k,t (•|s)⟩.Summing the inequality, we getEQUATIONProposition E.5. For the output of player2,E s∼µ0 max π V π,π K (s) -V * (s) ≤O CK (1 -γ) 2 ( log 1 c ηT + η (1 -γ) 2 + γβ η(1 -γ) 2 + η∥ Qπ ′ k,t ,π k -Q π ′ k,t ,π k ∥ ν π ′ k,t ,π k +η∥ Qπ k-1 ,π k-1 -Q∥ ν πk-1 ,π k-1 ) + γ K 1 -γ .where C is a problem-dependent constant.Proof. Combining Lemma E.8,Lemma E.9 and (44) gives the result.To deal with continuous action space, we propose a practical variant of SDEPO, named SDEPO-NN, to deal with continuous (or large-scale discrete ) action space. SDEPO-NN utilizes neural networks in policy π and the state-action value function approximation Q, detailed in Algorithm F.• It is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified. • For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates). • If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.Given a probability measure µ defined on X ⊆ R d , we denote L2(µ) as the set of functions f :X → R such that ∫ X (f (x)) 2 dµ(x) < ∞.The ℓ2 norm is the usual norm for square-summable infinite-dimensional vectors v indexed by the positive integers, such that ∥v∥2 ℓ 2 = ∑ ∞ i=1 v 2 i .This assumption can be removed if we use a β-greedy mixed version of π k , π k . As the β-greedy mixed policies are shifted with the original ones, we provide the convergence analyses of SDEPO with β-greedy exploration in Appendix E.1.The ℓ2 norm is the usual norm for square-summable infinite-dimensional vectors v indexed by the positive integers, such that ∥v∥2 ℓ 2 = ∑ ∞ i=1 v 2 i . 


Decision: The paper introduces the Spectral Dynamic Embedding Policy Optimization (SDEPO) algorithm to address two-player zero-sum Markov games with continuous state spaces, providing both theoretical guarantees and practical performance improvements. The approach leverages spectral embedding and offers a convergence rate that matches the best-known single-agent RL results, with a clear motivation and solid theoretical foundation. Overall, the paper is considered technically solid with the recommendation leaning towards acceptance.

Reviews: review: paper_summary: paper_summary: The authors introduce an innovative approach to solving 2p0s-MGs with continuous state spaces, providing both theoretical guarantees and practical improvements over existing methods. The SDEPO algorithm and its variants offer efficient and scalable solutions for complex Markov games, potentially applicable to various domains in reinforcement learning.
strength_weakness: strength_weakness: 1. The spectral embedding methods can be computationally intensive in practice due to the complexity of handling spectral dynamic embeddings.
2. Why were these specific feature generation methods chosen? Is the proposed method sensitive to feature generation methods?
3. The experiments are somewhat limited, expanding the empirical section to include more complex and diverse scenarios would significantly strengthen the paper.
questions: questions: Please see the weakness part.
limitations: limitations: Yes.
paper_summary: paper_summary: This paper proposes a new algorithm named Spectral Dynamic Embedding Policy Optimization (SDEPO) to solve the zero-sum Markov games with continous state and finite actions. The convergence analysis indicates that the proposed method achieves the best-known sample complexity as the case of finite-state space; this paper is the first theoretical result in handling the continuous state space with known dynamic and infinite horizon.
strength_weakness: strength_weakness: 1. The Assumption 1 is not reasonable. It says, whatever the state $s$ and the action $(a,b)$ are, the agent can move to any state $s'$ with a positive probability. Please correct me if I am wrong.

2. I am confused about what is new in the Spectral Dynamic Embedding method. It seems that both Bochner and Mercer theorem are well-known. This paper simply applies them to represent the transition probability and the reward using some kernels. Then everything is the same as traditional method in RL. 

3. A mild comment on Assumption 3: Since the optimal policy might be deterministic, it means that $\pi(a|s)$ is likely to be zero for some $a$. During the training, the policy $\pi_k$ will tend to the optimal policy; the mass at non-optimal action will also approach to $0$. It means if  $\underline{c}$ is larger than $\epsilon$, then $\pi_k$ will never converge to the optimal action in the sense of $L_\infty$ norm. From my understanding, the author needs to set $\underline{c}$ to be $\epsilon$ and it won't affect the complexity.
questions: questions: Q1:  Can the author justify the use of Assumption 1? It seems to be unrealistic in RL. 

Q2: This paper seems to be a simple combination of linear MDP + policy gradient. What is the novel part of this work? I feel hard to consider representing $\mathbb{P}$ and $r$ using kernel methods as a new thing.
limitations: limitations: NA
paper_summary: paper_summary: This paper studies the two-player zero-sum stochastic Markov games (2p0s-MGs) with large scale or continuous state spaces. These problems have a large cardinality and function approximation methods are needed. The paper consider a spectral dynamic embedding method and proposed SDEPO. This methods utilized the transition dynamics in the construction of the state-space value function. SDEPO is able to converge with order $1/\epsilon$, which matches the optimal rate in single agent RL.
Theorems are provided for the last iterate convergence of the SDEPO algorithm. The effectiveness of the algorithm has been verified in games against baseline methods.

Generally the paper is well structured, but section 3-4 should be better explained while section 5 and 6 focused on the "practical algorithms" is stretched quite far from the analytical results in the previous sections.
strength_weakness: strength_weakness: The only evaluation for the proposed algorithm is a rate of success in playing games with baseline algorithms. To the reviewer, this seems to be very limited. For a submission with strong theoretical focus, current result fails to validate the convergence properties of the proposed algorithm.

The sections for main results are not very well-written and is a bit difficult to read, more explanation would be appreciated. Although this could be due to the page limit.

Assumption 5 is in fact quite strong, a brief discussion on the impact and reasoning should be provided.

Section 5 and 6 seems a bit rushed and is intended to bring out the neural networks, the prior sections discussed the setting with tabular actions, where in these sections the action space is seen as continuous and more algorithms have been added, with no analytical results. I suggest the authors focusing on the existing setting with better presentation, explanation and more experiments.

Another problem this paper did not address is what are the current existing algorithms involving dynamics and function approximation in the single agent setting. The single agent RL with function approximation literature should be somewhat addressed in general.
questions: questions: What is the effect of truncation w.r.t. the later regularity condition assumptions? Does a more limiting truncation negatively impact these assumptions on the problem?

What is the reason for the consideration of one-sided $\epsilon$−optimal Nash equilibrium? The author stated that many existing works also consider this, but an explanation would be appreciated.
limitations: limitations: The authors did not fully address the limitations of this paper. The authors mentioned that the algorithm is non-independent and relies on some assumptions. There are some other limitations of this work, which has been raised in the previous sections.

Comments: comment: Dear Reviewer [sWhd],

Thank you again for your constructive feedback. The end of the discussion period is close. We would be grateful if we could hear your feedback regarding our answers to the reviews. Here, we clarify the following part of our work again.

[Novelty of this paper] The primary goal of our paper is to propose a provably efficient method for two-player zero-sum stochastic Markov games with a continuous state space. Achieving Nash equilibria of Markov games is inherently more complex than identifying optimal policies of MDPs due to the min-max nature of Markov games and the non-stationary environment that each player faces. Our method first generates the $m$ feature embeddings through Spectral Dynamic Embedding method, and then evaluates the Q functions based on the embeddings and applies natural gradient policy improvements for each player based on the evaluated Q functions. Unlike linear Markov games, where the underlying representation is assumed to be known, we need to choose a proper feature embedding method and a suitable policy optimization strategy. We prove that the two parts can interact well to achieve the convergence of policies to a Nash equilibrium with a desirable overall iteration complexity. Specifically, we theoretically analyze the policy evaluation error of each player’s Q-function based on the embeddings and bound the error of the one-side Nash equilibrium explicitly with the policy evaluation error. Note that even the kernel method and natural gradient method seems common in the RL literature, the merit of our method lies in identifying effective embeddings to enhance policy optimization and analyze the impact of evaluation error on the convergence to NE. 

[The rationality of Assumption 1] Assumption 1 assumes that next state is determined by certain deterministic function with an additional Gaussian noise. Such assumption is actually realistic (satisfied in many RL tasks) and is commonly assumed in both theory and practice in the RL literature. Your main concern of assumption 1 is that it requires the transition probability to remain consistently positive, attributed to the Gaussian noise. Nevertheless, due to the $3-\sigma$ rule, most of next state’s probability mass would concentrate within a limited area. Moreover, our results can be generalized to compactly supported transition randomness which admits certain kernelized representation (such as truncated Gaussian). We will add more discussion about this in our revision.


We sincerely hope to engage in further discussions with you. Thank you for your time and consideration!
We sincerely thank you for recognizing our efforts in rebuttal. We greatly appreciate your decision to increase the score. We will incorporate the discussion on the theoretical analyses and newly updated numerical results in our revision.
We thank you for your review and appreciate your time reviewing our paper. The end of the discussion period is close. We would be grateful if we could hear your feedback regarding our answers to the reviews. We are happy to address any remaining points during the remaining discussion period. Thanks a lot in advance!
Thank you for the detailed explanations provided by the authors and the newly updated numerical results. The authors have addressed my questions and I have decided to increase my score.