Paper deep dive
Decentralized Diffusion Policy Learning for Enhanced Exploration in Cooperative Multi-agent Reinforcement Learning
Yuyang Zhang, Haldun Balim, Na Li
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 7/8/2026, 11:29:50 AM
Summary
This paper addresses exploration limitations in cooperative multi-agent reinforcement learning (MARL) caused by projecting energy-based policies onto unimodal Gaussian distributions. It proposes Decentralized Diffusion Policy Learning (DDPL), which uses denoising diffusion probabilistic models to capture multi-modal action distributions for enhanced exploration. DDPL is trained online via Importance Sampling Score Matching (ISSM) and demonstrates improved performance across continuous-action MARL benchmarks.
Entities (7)
Relation Signals (6)
DDPL → uses → DDPM
confidence 95% · we propose decentralized diffusion policy learning (DDPL), which parameterizes each agent's policy with a denoising diffusion probabilistic model
DDPM → captures → Multi-modal Action Distributions
confidence 90% · an expressive generative model that captures multi-modal action distributions for enhanced exploration.
ISSM → enables → Online Training
confidence 90% · DDPL enables efficient online training of diffusion policies via importance sampling score matching (ISSM)
DDPL → evaluatedon → MPE, MuJoCo, IsaacLab, SMAX
confidence 90% · We evaluate DDPL on representative continuous-action MARL benchmarks, including multi-agent particle environment, multi-agent MuJoCo, IsaacLab, and JAX-reimplemented StarCraft multi-agent challenge
Gaussian Policy → hinders → exploration
confidence 90% · the limited expressiveness of Gaussian policies severely hinders exploration in DecSPG, and this limitation worsens as the number of agents grows.
DecSPG → uses → Gaussian Policy
confidence 90% · In practice, however, such energy-based policies are intractable to maintain and are commonly projected onto the Gaussian policy class.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Cooperative multi-agent reinforcement learning (MARL) involves complex agent interactions and requires effective exploration strategies. A prominent class of MARL algorithms, decentralized softmax policy gradient (DecSPG), addresses this through energy-based policy updates. In practice, however, such energy-based policies are intractable to maintain and are commonly projected onto the Gaussian policy class. In this work, we show that the limited expressiveness of Gaussian policies severely hinders exploration in DecSPG, and this limitation worsens as the number of agents grows. To address this issue, we propose decentralized diffusion policy learning (DDPL), which parameterizes each agent's policy with a denoising diffusion probabilistic model, an expressive generative model that captures multi-modal action distributions for enhanced exploration. DDPL enables efficient online training of diffusion policies via importance sampling score matching (ISSM), a novel training method with theoretical guarantee. We evaluate DDPL on representative continuous-action MARL benchmarks, including multi-agent particle environment, multi-agent MuJoCo, IsaacLab, and JAX-reimplemented StarCraft multi-agent challenge, and observe consistently improved performance.
Tags
Links
- Source: https://arxiv.org/abs/2605.07101v1
- Canonical: https://arxiv.org/abs/2605.07101v1
Trouble viewing inline? Open PDF directly →
Full Text
123,719 characters extracted from source content.
Expand or collapse full text
Decentralized Diffusion Policy Learning for Enhanced Exploration in Cooperative Multi-agent Reinforcement Learning Yuyang Zhang1,2 Haldun Balim1 Na Li1 1SEAS, Harvard University 2Kempner Institute, Harvard University yuyangzhang@g, hbalim@fas, nali@seas.harvard.edu Abstract Cooperative multi-agent reinforcement learning (MARL) involves complex agent interactions and requires effective exploration strategies. A prominent class of MARL algorithms, decentralized softmax policy gradient (DecSPG), addresses this through energy-based policy updates. In practice, however, such energy-based policies are intractable to maintain and are commonly projected onto the Gaussian policy class. In this work, we show that the limited expressiveness of Gaussian policies severely hinders exploration in DecSPG, and this limitation worsens as the number of agents grows. To address this issue, we propose decentralized diffusion policy learning (DDPL), which parameterizes each agent’s policy with a denoising diffusion probabilistic model, an expressive generative model that captures multi-modal action distributions for enhanced exploration. DDPL enables efficient online training of diffusion policies via importance sampling score matching (ISSM), a novel training method with theoretical guarantee. We evaluate DDPL on representative continuous-action MARL benchmarks, including multi-agent particle environment, multi-agent MuJoCo, IsaacLab, and JAX-reimplemented StarCraft multi-agent challenge, and observe consistently improved performance. 1 Introduction Multi-agent reinforcement learning (MARL) has found broad applicability across a range of domains (Gronauer and Diepold, 2022), including traffic networks (Chu et al., 2019), smart grids (Chen et al., 2021), smart buildings (Yu et al., 2020), and robotics (Orr and Dutta, 2023; Yang and Gu, 2004). In MARL, agents learn policies by interacting with the environment and adapting to the evolving behavior of other agents, with the goal of maximizing their rewards. Compared to single agent reinforcement learning, MARL poses additional challenges for exploration in the high-dimensional joint action space. In this paper, we focus on cooperative MARL, where agents share a common global reward and seek to maximize it collectively. In this setting, effective exploration is essential: agents must try diverse coordination patterns to discover high-reward joint actions, while premature commitment to a single pattern can trap the team in a suboptimal equilibrium. A prominent class of algorithms for this setting is decentralized softmax policy gradient (DecSPG) (Zhang et al., 2022) and its variants (Sun et al., 2023; Fox et al., 2022; Yang et al., 2024; Sahitaj et al., 2025; Liu et al., 2024; Zhong et al., 2024) These methods iteratively perform decentralized energy-based policy updates, which naturally generate multimodal action distributions for diverse exploration and provably converge to Nash equilibria when policy updates are exact. In practice, however, such energy-based action distributions are notoriously difficult to represent exactly (Song and Kingma, 2021). Most practical algorithms therefore project these energy-based policies onto simple policy classes, most commonly Gaussian policies (Liu et al., 2024; Zhong et al., 2024), favoring computational efficiency at the expense of expressiveness. Different from the original energy-based target policies, Gaussian policies lack expressiveness and often fail to capture multimodal action distributions. We find that this restriction can significantly hinder exploration and cause algorithms to converge prematurely to suboptimal equilibrium. The failure mode is particularly severe in multi-agent settings, where the number of viable coordination patterns grows with the number of agents, and the cost of missing them compounds accordingly. This leads to a central challenge when applying the DecSPG-based algorithms in practice: designing policy representations that are both expressive enough for effective exploration and complex coordination and tractable enough for online training and sampling. Fortunately, diffusion models offer a natural remedy (Ho et al., 2020; Austin et al., 2021). They are expressive generative models capable of representing complex, multimodal distributions while remaining efficient to train and sample from. This presents a great opportunity to parameterize reinforcement learning policies using diffusion models, especially in cooperative MARL where exploration with expressive policies is critical. The main obstacle to deploying diffusion policies in DecSPG-based algorithms is the lack of efficient online training algorithms. Standard training algorithms, e.g., score matching and its variants (Song et al., 2021; Lipman et al., 2022; Rombach et al., 2022), are primarily designed for offline settings and typically require samples from the target distribution, which are rarely available in online MARL. Most existing work either proposes alternative training methods that are computationally heavy (Wang et al., 2024; Zhang et al., 2025; Ding et al., 2025) or requires additional information, e.g., the gradient of the Q function, that is nontrivial to estimate (Psenka et al., 2023; Jain et al., 2024; Yang et al., 2023). Drawing on ideas from loss reweighting (Ma et al., 2025a; Dong et al., 2025), we derive an online variant of the score-matching loss that trains diffusion policies without requiring samples from the target policy distribution. Contributions. This paper formally shows that the limited expressiveness of standard policy parameterizations is a key obstacle to effective exploration in cooperative MARL, and develops a diffusion-based approach to overcome this limitation. Specifically, • First (Section˜3), we identify an explicit failure mode of unimodal policy classes and motivate the need for expressive, multimodal policy representations. We construct a representative environment in which Gaussian approximations to the energy-based policy updates of DecSPG severely restrict exploration and lead to convergence to suboptimal policies. Under mild conditions, we prove that the success probability of the resulting algorithm is at most 0.45n0.45^n, decaying exponentially to zero as the number of agents increases. • Then (Section˜4), we propose decentralized diffusion policy learning (DDPL), a DecSPG-based algorithm that parameterizes decentralized agent policies with denoising diffusion probabilistic models (DDPMs). To enable online training, we derive an importance-sampling-based score matching objective that moves the diffusion policy toward the target energy-based distribution without requiring samples from that distribution. We also provide theoretical sample-complexity guarantees for the proposed training method. • Finally (Section˜5), we evaluate DDPL on representative continuous-action MARL benchmarks, including the multi-agent particle environment (MPE), multi-agent MuJoCo (MaMuJoCo), IsaacLab bi-shadow-hand, and JAX-reimplemented StarCraft Multi-Agent Challenge (SMAX) (Rutherford et al., 2024), demonstrating improved performance and sample efficiency. Notations. For integer n, we use [n][n] to denote 1,2,⋯,n\1,2,·s,n\. Given agent set [n][n] and integer i, we use −i-i to denote agents [n]\i[n] \i\. We use 0,1n\0,1\^n to denote all binary vectors in ℝnR^n, i.e., 0,1n=x=(x1,⋯,xn)∈ℝn:xi∈0,1\0,1\^n=\x=(x_1,·s,x_n) ^n:x_i∈\0,1\\. We use n 1_n to denote the n-dimensional vector (1,⋯,1)(1,·s,1). We use ≲ to hide dependencies on absolute constants and log factors. 2 Preliminaries 2.1 Cooperative Multi-Agent Reinforcement Learning In cooperative MARL, agents collectively maximize a global reward. A standard model for this setting is the identical-interest Markov game with global state and decentralized policies (Matignon et al., 2012). The formal definition is as follows. Identical-Interest Markov Games. An identical-interest Markov game is defined by the tuple ℳ=(,ρ0,ii=1n,r,P,γ)M=(S,ρ^0,\A_i\_i=1^n,r,P,γ). Here S is the global state space and ρ0∈Δ()ρ^0∈ (S) is the initial state distribution. Each agent i∈[n]i∈[n] has an individual action space iA_i. We write the joint action as ≔(1,⋯,n)∈1×⋯×n≔ a ( a_1,·s, a_n) _1×·s×A_n . The reward function r:×↦[rmin,rmax]r:S×A [r_ ,r_ ] assigns a common reward to all agents, and the transition kernel P:×↦Δ()P:S×A (S) specifies the next-state distribution. The scalar γ∈(0,1)γ∈(0,1) is the discount factor. A policy π:↦Δ()π:S (A) is called decentralized if it factorizes across agents as π(|)=∏i=1nπi(i|)π( a| s)= _i=1^n _i( a_i| s). We denote such a policy by π=π1×⋯×πnπ= _1×·s× _n. For a decentralized policy π, we define the action-value function, value function, and objective as Qπ(,)≔π[∑t=0∞γt(r(t,t))|0=,0=],Vπ()≔∼π(⋅|)[Qπ(,)]Q^π( s, a) _π [ _t=0^∞γ^t (r( s^t, a^t) )| s^0= s, a^0= a ], V^π( s) _ a π (·| s ) [Q^π( s, a) ], and J(π)≔∼ρ0[Vπ()]J (π ) _ s ρ^0 [V^π( s) ], respectively. For each agent i, we further define the averaged action-value function, which evaluates the effect of agent i’s local action i a_i averaged over the actions of all other agents −i a_-i under the decentralized policy π: Q¯iπ(,i)≔−i∼π−i(⋅|)[Qπ(,(i,−i))]. Q_i^π( s, a_i) _ a_-i _-i(·| s) [Q^π( s,( a_i, a_-i)) ]. Nash Equilibrium (NE). A policy π⋆=π1⋆×⋯×πn⋆π =π _1×·s×π _n is a Nash equilibrium if J(πi⋆,π−i⋆)≥J(πi′,π−i⋆)J( _i , _-i )≥ J( _i , _-i ) for all i∈[n],πi′i∈[n], _i . Furthermore, we define the Nash Equilibrium Gap (NE-gap) of a policy π=π1×⋯×πnπ= _1×·s× _n as NE-gap(π)≔maxi∈[n]NE-gapi(π)NE-gap(π) _i∈[n]NE-gap_i(π), where NE-gapi(π)≔maxπi′J(πi′,π−i)−J(πi,π−i)NE-gap_i(π) _ _i J( _i , _-i)-J( _i, _-i). Policy π is called an ϵε-NE if NE-gap(π)≤ϵNE-gap(π)≤ε. 2.2 Denoising Diffusion Probabilistic Model The denoising diffusion probabilistic model (DDPM) is an expressive generative model for learning complex distributions. The forward process corrupts the data from the data distribution 0∼p0 a_0 p_0 towards the noise distribution (0,I)N(0,I) by adding Gaussian noise, i.e., τ=1−βττ−1+βτϵ a_τ= 1- _τ a_τ-1+ _τε. Here βτ=1H\ _τ\_τ=1^H is the variance schedule and ϵ∼(0,I)ε (0,I). With ατ=1−βτ _τ=1- _τ, α¯τ=∏h=1ταh α_τ= _h=1^τ _h, we can express τ a_τ as τ=α¯τ0+1−α¯τϵ a_τ= α_τ a_0+ 1- α_τε. We denote the marginal distribution of τ a_τ as pτp_τ and the conditional distribution of τ|0 a_τ| a_0 as pτ|0p_τ|0. The reverse process denoises samples from H∼(0,I) a_H (0,I) towards data distribution p0p_0: τ−1=1ατ(τ+βτsτθ(τ))+στϵ a_τ-1= 1 _τ ( a_τ+ _τs_τ^θ( a_τ) )+ _τε. Here στ _τ is a hyperparameter and we let στ2=βτ/ατ _τ^2= _τ/ _τ in this paper. sτθ(τ)s_τ^θ( a_τ) is the score network with parameters θ, which can be trained via denoising score matching (Ho et al., 2020) by minimizing ℒmarginal(θ)=∑τ=1Hτ∼pτβτ2στ2ατ‖sτθ(τ)−∇τlogpτ(τ)‖2L^marginal(θ)= _τ=1^HE_ a_τ p_τ _τ^2 _τ^2 _τ \|s_τ^θ ( a_τ )- _ a_τ p_τ( a_τ) \|^2. It is shown in Ma et al. (2025b) that minimizing ℒmarginalL^marginal is equivalent to minimizing the following conditional loss that is more tractable: ℒ(θ)=0∼p0∑τ=1Hβτ2στ2ατϵτ∼(0,I)‖sτθ(α¯τ0+1−α¯τϵτ)−ϵτ1−α¯τ‖2. splitL(θ)=&E_ a_0 p_0 _τ=1^H _τ^2σ^2_τ _τE_ _τ (0,I) \|s_τ^θ ( α_τ a_0+ 1- α_τ _τ )- _τ 1- α_τ \|^2. split (1) 3 Comparing Multi-Modal Policies with Gaussian Policies in MARL We first present decentralized softmax policy gradient (DecSPG) (Zhang et al., 2022), which maintains decentralized energy-based policies and converges to NEs under exact policy updates. In practice, the exact updates are intractable when action space is large or even continuous, and the policies are commonly projected onto the Gaussian policy class. In this section, we investigate how this projection affects performance through concrete examples. We show that Gaussian policies cause severe under-exploration, and this failure mode amplifies exponentially with the number of agents. 3.1 Decentralized Softmax Policy Gradient Algorithm˜1 iteratively estimates each agent’s averaged Q function Q^iπk−1 Q_i^π^k-1 and updates each agent’s policy towards the energy-based target πi⋆,k∝πik−1exp(η1−γQ^iπk−1)π ,k_i _i^k-1 ( η1-γ Q_i^π^k-1). Each agent thereby captures the influence of other agents’ policies π−ik−1 _-i^k-1 through the averaged Q function and gradually reweights toward actions with higher rewards under π−ik−1 _-i^k-1. This energy-based target is generally intractable, and the standard approach is to project it onto a parametric class Πi _i. Algorithm 1 Decentralized Softmax Policy Gradient (Dec-SPG) (Zhang et al., 2022) 1: Init: Initial policy π0=π10×⋯×πn0π^0= _1^0×·s× _n^0, learning rate η, batch size B; 2: for epoch k∈[K]k∈[K] do 3: Take B steps with policy πk−1π^k-1 and store to buffer =,,r,′D=\ s, a,r, s \; 4: Estimate average Q function Q¯iπk−1 Q_i^π^k-1 for each agent i with data in D, denoted by Q^iπk−1 Q_i^π^k-1; 5: Update policy πik(⋅|) _i^k(·| s) for each agent i: πik(⋅|)←argminπi∈ΠiKL(πi∥πi⋆,k(⋅|)),πi⋆,k(⋅|)∝πik−1(⋅|)exp(η1−γQ^iπk−1(,⋅)). split _i^k(·| s)← _ _i∈ _iKL ( _i\| _i ,k(·| s) ), _i ,k(·| s) _i^k-1(·| s) ( η1-γ Q_i^π^k-1( s,·) ). split When the Q estimation and policy updates are exact, Algorithm˜1 provably converges to a Nash equilibrium given finite state and action spaces (Theorem˜C.1). The average squared NE-gap decays at rate (n/K)O(n/K), where K and n are the number of epochs and agents, respectively. This guarantee, however, requires the policy update to be exact, i.e., πik=πi⋆,k _i^k= _i ,k, which is rarely satisfied in practice due to the policy projection onto Πi _i. 3.2 Multi-Modal Policies Improve Exploration We now investigate how the policy class Πi _i affects the performance of Algorithm˜1, comparing all Gaussian policies against all policies. Although Gaussians are popular for their efficiency, we show that their unimodal nature causes insufficient exploration, breaking the guarantee of Algorithm˜1. We first build intuition with a single-agent MDP, then extend to an identical-interest game where the failure amplifies exponentially with the number of agents. 3.2.1 A Single-Agent Example We construct MDP ℳsaM_sa with three states =0,1,2S=\ s_0, s_1, s_2\, fixed initial state 0 s_0, action space 1=ℝA_1=R, and discount factor γ=0.1γ=0.1. The reward function r(,)r( s, a) is defined as: r(,)=(=0)rstart()+(=1)c1+(=2)0 splitr( s, a)=I ( s= s_0 )r_start( a)+I ( s= s_1 )c_1+I ( s= s_2 )0 split with rstart()=−ν2+c0exp(−α(+β)2)+c0exp(−α(−β)2)r_start( a)=-ν a^2+c_0 (-α( a+β)^2)+c_0 (-α( a-β)^2) and parameters ν=10,α=250,β=0.5,c0=2.7ν=10,α=250,β=0.5,c_0=2.7, and c1=20000c_1=20000. The transition probability is defined as: ℙ(1|0,)=pstart(),ℙ(2|0,)=1−ℙ(1|0,),ℙ(2|1,)=ℙ(2|2,)=1. ( s_1| s_0, a)=p_start( a), ( s_2| s_0, a)=1-P( s_1| s_0, a), ( s_2| s_1, a )=P ( s_2| s_2, a )=1. where pstart()=cpexp(−α(+β)2)⋅(|+β|≤1/4)p_start( a)=c_p (-α( a+β)^2)·I (| a+β|≤ 1/4 ) with cp=0.005c_p=0.005. Intuitively, this MDP rewards actions near ±β±β and grants a small probability of transitioning to ‘goal’ state 1 s_1 for actions near −β-β. Otherwise the environment transitions to the absorbing state 2 s_2 with zero reward. Learning an optimal policy requires extensive exploration in 0 s_0. Figure 1: Policies during training with Gaussian policies and all policies in the single-agent example. We then run Algorithm˜1 with batch size B=200B=200, η=0.01(1−γ)η=0.01(1-γ), and Πi _i being either Gaussian policies or all policies. We plot representative learned policies (Figure˜1) and the average reward of the policies (Figure˜2, left) during training. Looking at the trained policies (Figure˜1), we observe: • Initially (Epoch 20), both policy classes lack action samples near ±β±β, and Q^(0,)≈−ν2 Q( s_0, a)≈-ν a^2. This leads to a zero-mean Gaussian policy proportional to exp(−η/(1−γ)⋅kν2) (-η/(1-γ)· kν a^2 ) in epoch k. • Later (Epoch 100), more samples near ±β±β accumulate and Q^(0,) Q( s_0, a) captures reward terms c0exp(−α(±β)2)c_0 (-α( a±β)^2 ). Therefore, with all policies, the energy-based reweighting term exp(ηQ^/(1−γ)) (η Q/(1-γ)) gradually increases action probability near =±β a=±β, resulting in a multi-modal exploration distribution. In contrast, with Gaussian policies, the updated policy is stuck near 0. The reward signal near =±β a=±β in exp(ηQ^/(1−γ)) (η Q/(1-γ)) is too weak to take effect in the Gaussian policy projection. Therefore, the policy continues to shrink its variance due to the −ν2-ν a^2 term in Q Q, leading to insufficient exploration near ±β±β. • Finally (Epoch 130 & 150), with more samples near −β-β, Algorithm˜1 with all policies reaches state 1 s_1 and converges to optimal action =−β a=-β in state 0 s_0. However, Algorithm˜1 with Gaussian policies continues to output policy centered at 0. Figure 2: Expected return during training for Algorithm˜1 with all policies and Gaussian policies. Left: single-agent. Right: multi-agent with n=2n=2. The following theorem confirms this failure mode theoretically. It shows that starting from a zero-mean Gaussian with modest variance, similar to the first plot in Figure˜1, the Gaussian-projected policy update keeps the zero mean and monotonically shrinks the variance. The policy converges to a delta distribution around 0, leading to a low probability of reaching 1 s_1 throughout the entire algorithm. The assumption that Q^1π Q_1^π is accurate for reward rsa¬r _sa isolates the failure mode of Gaussian policies by only considering the estimation error at state 1 s_1. This is a natural condition before 1 s_1 is reached. In practice, other estimation errors compound the issue. Theorem 3.1. Consider Algorithm˜1 with batch size B=32B=32, learning rate η=0.01(1−γ)η=0.01(1-γ), and Π1 _1 being Gaussian policy class. Consider any epoch k before 1 s_1 is first reached and any policy π1k=(μk,(σk)2) _1^k=N(μ^k,(σ^k)^2) with μk=0μ^k=0. Suppose Q^1π1k Q_1^π^k_1 is accurate for ℳsaM_sa with reward rsa¬(,)=(=0)rstart()r _sa( s, a)=I ( s= s_0 )r_start( a). Performing policy update (line 5, Algorithm˜1) with gradient descent yields a Gaussian policy π1k+1=(μk+1,(σk+1)2)π^k+1_1=N(μ^k+1,(σ^k+1)^2) with μk+1=0,1(σk+1)2≥1(σk)2+0.128. splitμ^k+1=0, 1 (σ^k+1 )^2≥ 1 (σ^k )^2+0.128. split Moreover, starting from π10=(0,0.1) _1^0=N(0,0.1), the probability of reaching state 1 s_1 is at most 0.450.45.∎ 3.2.2 A Multi-Agent Example We now extend the above example to an identical-interest Markov game ℳmaM_ma with agents [n][n], state space =−n∪0,1nS=\ -1_n\∪\0,1\^n, fixed initial state start=−n s_start= -1_n, and individual action space i=ℝA_i=R. The reward function is defined as: r(,)=(=−n)∑i=1nrstart(i)+(≠±n)c1∑i=1n[]i+(=n)(c2+c3exp(−⊺)), splitr( s, a)=I ( s=- 1_n ) _i=1^nr_start( a_i)+I ( s≠± 1_n )\,c_1 _i=1^n[ s]_i+I ( s= 1_n ) (c_2+c_3 (- a a ) ), split with parameters c1=20000,c2=30000c_1=20000,c_2=30000, and c3=10000c_3=10000. The transition probability is defined as: ℙ([]i=1|=−n,)=pstart(i),ℙ([]i=0|=−n,)=1−pstart(i),ℙ(|,)=1,∀≠−n,∀. split&P ([ s]_i=1| s= -1_n, a )=p_start( a_i), ([ s]_i=0| s= -1_n, a )=1-p_start( a_i),\\ &P ( 0| s, a )=1, ∀ s≠ -1_n,∀ a. split ℳmaM_ma embeds n copies of the single-agent MDP at the initial state −n -1_n. Each agent i independently transitions coordinate []i[ s]_i to 11 with probability pstart(i)p_start( a_i). The ‘goal’ state =n s= 1_n is reached only when all agents succeed. Therefore, each agent learns independently before =n s= 1_n is reached, which we will show amplifies the single-agent failure mode exponentially in n. Proposition 3.2. Consider Algorithm˜1 with batch size B=32B=32, learning rate η=0.01(1−γ)η=0.01(1-γ), and Πi _i being the Gaussian policy class. For every i∈[n]i∈[n] and every epoch k where []i=1[ s]_i=1 has not been reached, suppose Q^iπk Q_i^π^k is accurate for ℳmaM_ma with reward rma¬(,)=(=−n)∑j=1nrstart(j)+(≠−n)c1∑j≠i[]jr_ma ( s, a)=I ( s=- 1_n ) _j=1^nr_start( a_j)+I ( s≠- 1_n )\,c_1 _j≠ i[ s]_j. Then starting from πi0=(0,0.1)i∈[n]\ _i^0=N(0,0.1)\_i∈[n], the probability of reaching state =n s= 1_n using Algorithm˜1 and Gaussian policies is at most 0.45n0.45^n.∎ Here, assuming Q^iπk Q_i^π^k being accurate for reward rma¬r_ma isolates the failure mode by only considering the estimation error in reward terms c1[]ic_1[ s]_i and (c2+c3exp(−⊺))(c_2+c_3 (- a a )), which is a natural condition before []i[ s]_i is reached. In practice, other estimation errors compound the issue. Under this assumption, the success probability of reaching =n s= 1_n decays exponentially in the number of agents, indicating that Gaussian policies severely hurt the performance of the algorithm in the multi-agent setting. We confirm this empirically with n=2n=2 in Figure˜2 (right): Algorithm˜1 with Gaussian policies stays near zero return throughout training, while Algorithm˜1 with all policies discovers the high-reward equilibrium and converges to near-optimal performance. 4 Decentralized Diffusion Policy Learning To perform the energy-based policy updates in Algorithm˜1 efficiently and accurately, we parameterize each agent’s policy using a DDPM, an expressive generative model capable of representing complex multi-modal distributions. We then develop a provably efficient online training algorithm for diffusion models based on score matching, and present decentralized diffusion policy learning (DDPL), an online MARL algorithm with decentralized diffusion policies. Without loss of generality, we consider bounded action space i⊆∈ℝd:‖≤cbA_i \ a ^d: \| a \|≤ c_b\ for each agent i throughout this section. Many popular MARL benchmarks, including MPE, MaMuJoCo, IsaacLab, SMAX, and real-world applications involve such bounded action spaces. 4.1 Importance Sampling Score Matching (ISSM) Algorithm Design. Consider policy πiθi _i _i parameterized by a DDPM with score network si,τθis_i,τ _i and parameter θi∈Θ _i∈ . Given state s, sampling from i∼πiθi(⋅|) a_i _i _i(·| s) is a two-step procedure. First, we sample i,0 a_i,0 via the H-step reverse process i,τ−1=1ατ(i,τ+βτsi,τθi(i,τ|))+στϵ a_i,τ-1= 1 _τ( a_i,τ+ _τs_i,τ _i( a_i,τ| s))+ _τε and the noise schedule in LABEL:eq:scheduler. Second, we project i,0 a_i,0 onto iA_i: i=Proji(i,0)=argmin∈i‖−i,0‖ a_i=Proj_A_i( a_i,0)= _ a _i \| a- a_i,0 \|. In epoch k, we aim to train si,τθi(⋅|)s_i,τ _i(·| s) so that distribution πik(⋅|)=πiθi(⋅|) _i^k(·| s)= _i _i(·| s) satisfies πik(⋅|)∝πik−1(⋅|)exp(η1−γQ^iπk−1(,⋅)) _i^k(·| s) _i^k-1(·| s) ( η1-γ Q_i^π^k-1( s,·)) for every state ∈ s . For simplicity, we consider the score network training for a specific agent i and state ∈ s . A core challenge in applying standard score matching is that minimizing the denoising score matching loss ℒi,(θi)L_i, s( _i) (Equation˜1) requires samples from the target distribution i∼πi⋆,k(⋅|) a_i _i ,k(·| s). Such samples are unavailable before the policy update. To tackle this issue, we adopt the idea of importance sampling by drawing samples from an easy-to-sample distribution p~(⋅|) p(·| s), which will be specified later. With samples i(j)∼p~(⋅|)j=1N\ a_i^(j) p(·| s)\_j=1^N, we compute the following importance sampling score matching loss: ℒi,is(θi)≔1N∑j=1Nπi⋆,k(i(j)|)p~(i(j)|)⋅∑τ=1Hβτ2στ2ατϵτ∼(0,I)∥si,τθi(α¯τi(j)+1−α¯τϵτ|)−ϵτ1−α¯τ∥2. split&L_i, s^is( _i) 1N _j=1^N _i ,k( a_i^(j)| s) p( a_i^(j)| s)· _τ=1^H _τ^2 _τ^2 _τE_ _τ (0,I) \|s_i,τ _i ( . α_τ a_i^(j)+ 1- α_τ _τ | s )- _τ 1- α_τ \|^2. split (2) In expectation, ℒi,is(θi)L_i, s^is( _i) is equivalent to the original score matching loss ℒi,(θi)L_i, s( _i) (Section˜B.1). As a result, one only needs to minimize the ℒi,is(θi)L_i, s^is( _i) within Θ and get minimizer θ^i≔argminθi∈Θℒi,is(θi) θ_i _ _i∈ L_i, s^is( _i). Actions i∼πiθ^i(⋅|) a_i _i θ_i(·| s) can then be sampled using the reverse process with score network si,τθ^i(⋅|)s_i,τ θ_i(·| s), followed by a projection onto action space iA_i. Theoretical Guarantee. We now show that the importance sampling score matching is provably efficient. Due to space limit, we defer Section˜B.2 to the appendix, which poses standard and mild assumptions on the expressiveness, complexity, and Lipschitzness of score function class ΘS_ . Theorem 4.1. Consider score function class ΘS_ satisfying Section˜B.2 with constants c1c_1 and c2c_2. Consider target distribution πi⋆,k(⋅|) _i ,k(·| s) and sampling distribution p~(⋅|) p(·| s), both supported on iA_i. Let cw≔maxw()c_w _ aw( a), where w()≔πi⋆,k(|)/p~(|)w( a) π ,k_i( a| s)/ p( a| s). If sample size N≳log(dNlogH/δ)N (d N H/δ), then πiθ^i(⋅|) _i θ_i(·| s) satisfies the following with probability at least 1−δ1-δ DTV2(πi⋆,k(⋅|)∥πiθ^i(⋅|))≲d2log6H2+dexp(34D4(πi⋆,k(⋅|)∥p~(⋅|)))N. splitD_TV^2 ( _i ,k(·| s)\| _i θ_i(·| s) ) & d^2 ^6HH^2+ d ( 34D_4 ( _i ,k(·| s)\| p(·| s) ) ) N. split Here d is action dimension, H is diffusion denoising steps, and D4(p∥q)=13lnx∼q[(p(x)/q(x))4]D_4(p\|q)= 13 _x q [(p(x)/q(x))^4 ] is the fourth-order Renyi divergence. ≲ and ≳ absorbs constants c1,c2,cb,cwc_1,c_2,c_b,c_w, and log(H,N,1/δ) (H,N,1/δ). ∎ The above theorem bounds error of the learned policy with two terms. The first term is the DDPM discretization error which decays with the number of denoising steps H. The second term is the statistical error which is proportional to 1/N1/ N where N is the sample size from sampling distribution p~(⋅|) p(·| s). We can therefore guarantee the accuracy of the learned target distribution πiθ^i(⋅|) _i θ_i(·| s) with large denoising steps H and number of samples N. The one-step policy update error in Theorem˜4.1 also decays with D4(πi⋆,k(⋅|)∥p~(⋅|))D_4( _i ,k(·| s)\| p(·| s)). The closer p~(⋅|) p(·| s) is to πi⋆,k(⋅|) _i ,k(·| s), the fewer samples we need for an accurate update. This motivates choosing πik−1(⋅|) _i^k-1(·| s) as the sampling distribution. Theoretical guarantee of Algorithm˜1 (Theorem˜C.1) suggests a decreasing distance between the current policy πik−1(⋅|) _i^k-1(·| s) and the target policy πi⋆,k(⋅|) _i ,k(·| s), since the sequence of policies πi⋆,k(⋅|)\ _i ,k(·| s)\ converges. Therefore, with p~(⋅|)=πik−1(⋅|) p(·| s)= _i^k-1(·| s), we can expect accurate policy update with less and less samples during training. With this choice, the importance sampling weight in LABEL:eq:flow_matching_is simplifies to πi⋆,k(i(j)|)/p~(i(j)|)=exp(η1−γQ^iπk−1(,i(j))) _i ,k( a_i^(j)| s)/ p( a_i^(j)| s)= ( η1-γ Q_i^π^k-1( s, a_i^(j))). We further quantify how the one-step policy update error propagates through epochs and influences NE-gap of the learned policies. We develop Theorem˜C.2 building on previous literature Theorem˜C.1. Specifically, if diffusion updates are ϵis _is-accurate in TV distance, the squared average NE-gap, NE-gap2(πk)/KNE-gap^2(π^k)/K, scales as (n/K+n2ϵis)O(n/K+n^2 _is) where n is the number of agents and K is the total number of episodes. Compared to the original rate in Theorem˜C.1, Theorem˜C.2 only incurs an additive error linear in ϵis _is. Since NE-gap measures the deviation of a policy from a Nash equilibrium, this bound implies that the learned policies are approximate Nash equilibria as K grows and ϵis _is shrinks. Due to space constraints, we defer the formal statement and its proof to Appendix C. 4.2 The Full Decentralized Diffusion Policy Learning Algorithm The previous subsection establishes per-state updates for Algorithm˜1 with decentralized diffusion policies. To obtain a more practical algorithm, we extend Algorithm˜1 to large state spaces with neural network function approximation and replace exact ISSM loss minimizations with gradient steps, yielding Algorithm˜2. In epoch k, each agent i maintains a Q-network Q^iψik(,i) Q_i _i^k( s, a_i) and a diffusion policy πiθik(i|) _i _i^k( a_i| s) with score network si,τθik(i|)s_i,τ _i^k( a_i| s), both taking state s as input to generalize across the state space. The agent first estimates the Q function (line 5), then updates the policy via one gradient step on the importance sampling score matching loss with πiθik−1 _i _i^k-1 as sampling distribution (line 6). Both updates are decentralized because they only require local actions i a_i. Algorithm 2 Decentralized Diffusion Policy Learning (DDPL) 1: Init: Policy parameters θ0=θi0i∈[n]θ^0=\ _i^0\_i∈[n], Q parameters ψ0=ψi0i∈[n]ψ^0=\ _i^0\_i∈[n], learning rate η, rollout batch size BrB_r, update batch size BuB_u, smoothing factor ξ. 2: for epoch k∈[K]k∈[K] do 3: Rollout the current policy πθk−1π^θ^k-1 for BrB_r steps and store to buffer D; 4: For each agent i, sample batch ℬi←(,i,r,′)B_i←\( s, a_i,r, s )\ of size BuB_u from D and next actions ′,i(j)∼πiθik−1j∈[N]\ a_ s ,i^(j) _i _i^k-1\_j∈[N] for all ′∈ℬi s _i; 5: Q update. For each agent i, minimize ℒiQ(ψik−1)←∑(,i,r,′)∈ℬi‖Q^iψik−1(,i)−r−γN∑j=1NQ^iψ¯ik−1(′,′,i(j))‖2,L_i^Q( _i^k-1)← _( s, a_i,r, s ) _i \| Q_i _i^k-1( s, a_i)-r- γN _j=1^N Q_i ψ_i^k-1( s , a_ s ,i^(j)) \|^2, and update ψik←Adam(ψik−1,∇ℒiQ) _i^k ( _i^k-1, _i^Q), ψ¯ik←ξψik+(1−ξ)ψ¯ik−1 ψ_i^k←ξ _i^k+(1-ξ) ψ_i^k-1. 6: Policy update. For each agent i, minimize the importance sampling score matching loss ℒ′,iis(θik−1)←1N∑j=1Nexp(η1−γQ^iψik(′,′,i(j)))l′,iθik−1(′,i(j)),where l,iθ()≔∑τ=1Hβτ2στ2ατϵτ∼(0,I)∥si,τθ(α¯τ+1−α¯τϵτ|)−ϵτ1−α¯τ∥2. split&L_ s ,i^is( _i^k-1)← 1N _j=1^N ( η1-γ Q_i _i^k ( s , a_ s ,i^(j) ) )l _i^k-1_ s ,i( a_ s ,i^(j)),\\ &where l_ s,i^θ( a) _τ=1^H _τ^2 _τ^2 _τE_ _τ (0,I) \|s_i,τ^θ ( . α_τ a+ 1- α_τ _τ | s )- _τ 1- α_τ \|^2. split and update θik←Adam(θik−1,∑′∇ℒ′,iis) _i^k ( _i^k-1, _ s _ s ,i^is). 5 Simulations Experiment setup. We implement DDPL (Algorithm˜2) in JAX and evaluate on 8 continuous action tasks including 4 MaMuJoCo, 1 MPE, 1 IsaacLab, and 2 SMAX tasks. Hyperparameters are listed in Table˜3, among which diffusion policy update learning rate η is chosen after light tuning and all other hyperparameters follow those in the codebase (Wang et al., 2024) without tuning. We train DDPL for 5×1065× 10^6, 1×1071× 10^7, and 5×1075× 10^7 environment steps on MPE, MaMuJoCo/IsaacLab, and SMAX tasks, respectively. Training was performed on RTX 5090 GPUs, and the training time ranges from 1 to 8 hours per seed. For each environment, we computed the average episode returns over 50 trajectories. Its mean and std. are then computed across 3 random seeds. We report the mean and std. of average episode returns during training in Figure˜3, where curves show means and shaded areas show std., and report the highest mean during training in Table˜1. Our training results are benchmarked against two on-policy baselines, HAPPO Liu et al. (2024) and MAPPO Yu et al. (2022), and one off-policy baseline, HASAC Liu et al. (2024). Main results. From Table˜1, Algorithm˜2 attains the best mean return on 7 of 8 tasks, with a performance improvement of 25.78%25.78\% on IsaacLab-ShadowHandOver, 132.56%132.56\% on 3s_vs_5z, and 155.84%155.84\% on 6h_vs_8z over the best baseline. Moreover, the shapes of the learning curves suggest that DDPL’s gains come from sufficient exploration rather than early commitment to suboptimal solutions. For example, on 3s_vs_5z, Algorithm˜2 obtains lower returns than the baselines in the initial exploration stage. After exploration, the performance then improves rapidly and converges near the maximum return of 22. We attribute the early lag to broader exploration in the action space, which is the cost paid for discovering higher-reward equilibria that unimodal baselines never reach. Similar dynamics appear on the IsaacLab and MaMuJoCo tasks. Figure 3: Average episode return during training for Algorithm˜2 and baselines. Curves and shaded areas show the mean and std. over three random seeds, respectively. Environment DDPL HASAC HAPPO MAPPO Ant-v2-2x4 6650.32±94.18 6650.32± 94.18 6301.93±422.646301.93± 422.64 5881.92±87.505881.92± 87.50 5215.60±922.765215.60± 922.76 Ant-v2-4x2 6460.28±242.05 6460.28± 242.05 6408.54±161.96 6408.54± 161.96 6176.43±22.576176.43± 22.57 6011.51±195.106011.51± 195.10 Ant-v2-8x1 6633.27±121.71 6633.27± 121.71 6204.69±153.576204.69± 153.57 5949.69±45.615949.69± 45.61 4977.60±706.164977.60± 706.16 HumanoidStandup-v2-17x1 172373.97±24918.25 172373.97± 24918.25 160698.78±2993.16160698.78± 2993.16 167491.90±10143.52167491.90± 10143.52 116959.76±11131.91116959.76± 11131.91 Simple-Reference-v2 −10.35±0.52-10.35± 0.52 −9.59±0.35 -9.59± 0.35 −15.89±1.31-15.89± 1.31 −14.20±1.92-14.20± 1.92 ShadowHandOver 669.25±44.82 669.25± 44.82 532.07±26.87532.07± 26.87 18.07±1.6518.07± 1.65 16.68±0.6016.68± 0.60 3s_vs_5z 2.00±0.00 2.00± 0.00 0.86±0.060.86± 0.06 0.59±0.030.59± 0.03 0.42±0.010.42± 0.01 6h_vs_8z 1.97±0.02 1.97± 0.02 0.77±0.030.77± 0.03 0.71±0.030.71± 0.03 0.65±0.020.65± 0.02 Table 1: Highest mean of the episode return with standard deviation over 3 random seeds. Ablation Study: diffusion vs. Gaussian policies. To provide further evidence on the mechanism of DDPL’s performance improvement, we probe the diffusion policies from DDPL during training and compare them with Gaussian policies from DDPL with the Gaussian policy class. Figure˜4 (left) shows the average return, and Figure˜4 (right) shows the action distribution at a representative state for two action dimensions from two agents. The Gaussian variant converges early to a suboptimal equilibrium with episode return ≈0.65≈ 0.65, and subsequent updates rarely changes its unimodal distribution. Contrarily, DDPL with diffusion policies maintains multi-modal action distributions for exploration in early training, discovers higher-reward modes the Gaussian variant never visits, and converges to a qualitatively different equilibrium with return ≈2≈ 2. This result matches our theoretical analysis in Section˜3, verifying that the multimodality of diffusion policies are essential to more effective exploration and improved performance. Sensitivity analysis. We test DDPL’s sensitivity to two diffusion hyperparameters on 3s_vs_5z: the number of denoising steps H and the noise schedule. Figure˜5 shows that final return is stable across H∈10,20,30H∈\10,20,30\ (with cosine schedule) and across linear, cosine, and VP schedules (with H=20H=20). This indicates that DDPL’s performance does not rely on careful tuning of these hyperparameters. Training and inference time. We measure the wall-clock time of Algorithm˜2, averaged over 10000 training and inference steps on IsaacLab Bi-ShadowHandOver with an AMD EPYC 9275F CPU (24 cores, 48 threads) and an NVIDIA RTX 5090 GPU (32GB memory). We compare with DACER (Wang et al., 2024) and QVPO (Ding et al., 2024), two single-agent diffusion-based algorithms applied to the environment in a centralized manner, as well as DDPL with the diffusion policy class replaced by Gaussian policy class. All algorithms are implemented in JAX with inference and training jitted, and all diffusion policies use 20 denoising steps. Results are shown in Figure˜5. The inference time of DDPL, which remains below 1 ms, is moderate. The training overhead is 39% smaller than DACER and 81% faster than QVPO, which we attribute to ISSM’s efficient online training. Figure 4: LEFT: the average episode return of DDPL with diffusion and Gaussian policies; RIGHT: the action distribution at a representative state for two action dimensions from two agents. Figure 5: Sensitivity to denoising steps H (left) and noise schedule (right) in SMAX-3s_vs_5z. Algorithm Train (ms) Inf. (ms) DDPL 0.820 3.684 DACER 0.814 6.045 QVPO 0.777 19.720 DDPL-Gaussian 0.374 2.599 Table 2: Training and inference time in IsaacLab Bi-ShadowHand. 6 Conclusions & Limitations This paper studies exploration in cooperative MARL through the lens of policy expressiveness. We show that a common implementation choice of projecting energy-based policies into Gaussian policy class can substantially hinder exploration and lead to degraded performance. The failure mode becomes more severe as the number of agents grows. To address this, we propose decentralized diffusion policy learning (DDPL), which parameterizes decentralized policies with diffusion models and trains them online via a novel importance sampling score matching (ISSM) objective with provable sample complexity and NE-gap guarantees. Empirically, DDPL improves performance and sample efficiency on representative continuous-action MARL benchmarks (MPE, MaMuJoCo, IsaacLab, and SMAX), highlighting the benefit of multi-modal diffusion policies for enhanced exploration and discovering higher-reward equilibria. Limitations. DDPL parameterizes fully decentralized policies. Although they require minimal communication, they do not exploit network MARL settings where low-cost communication with neighboring agents is encouraged. Extending the framework to such partially decentralized policies is a natural next step. Also, the paper focuses on the identical-interest setting. Extending to general-sum or competitive settings is left to future work. References J. Austin, D. D. Johnson, J. Ho, D. Tarlow, and R. van den Berg (2021) Structured denoising diffusion models in discrete state-spaces. In Advances in Neural Information Processing Systems, M. Ranzato, A. Beygelzimer, Y. Dauphin, P.S. Liang, and J. W. Vaughan (Eds.), Vol. 34, p. 17981–17993. Cited by: §1. P. L. Bartlett, D. J. Foster, and M. J. Telgarsky (2017) Spectrally-normalized margin bounds for neural networks. Advances in neural information processing systems 30. Cited by: §B.2. D. Chen, K. Chen, Z. Li, T. Chu, R. Yao, F. Qiu, and K. Lin (2021) Powernet: multi-agent deep reinforcement learning for scalable powergrid control. IEEE Transactions on Power Systems 37 (2), p. 1007–1017. Cited by: §1. S. Chen, S. Chewi, J. Li, Y. Li, A. Salim, and A. R. Zhang (2022) Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. arXiv preprint arXiv:2209.11215. Cited by: §B.2. T. Chu, J. Wang, L. Codecà, and Z. Li (2019) Multi-agent deep reinforcement learning for large-scale traffic signal control. IEEE transactions on intelligent transportation systems 21 (3), p. 1086–1095. Cited by: §1. S. Ding, K. Hu, Z. Zhang, K. Ren, W. Zhang, J. Yu, J. Wang, and Y. Shi (2024) Diffusion-based reinforcement learning via q-weighted variational policy optimization. Advances in Neural Information Processing Systems 37, p. 53945–53968. Cited by: §5. S. Ding, K. Hu, S. Zhong, H. Luo, W. Zhang, J. Wang, J. Wang, and Y. Shi (2025) GenPO: generative diffusion models meet on-policy reinforcement learning. arXiv preprint arXiv:2505.18763. Cited by: §1. X. Dong, J. Cheng, and X. S. Zhang (2025) Maximum entropy reinforcement learning with diffusion policy. arXiv preprint arXiv:2502.11612. Cited by: §1. B. L. Edelman, S. Goel, S. Kakade, and C. Zhang (2022) Inductive biases and variable creation in self-attention mechanisms. In International Conference on Machine Learning, p. 5793–5831. Cited by: §B.2. R. Fox, S. M. Mcaleer, W. Overman, and I. Panageas (2022) Independent natural policy gradient always converges in markov potential games. In International Conference on Artificial Intelligence and Statistics, p. 4414–4425. Cited by: §1. S. Gronauer and K. Diepold (2022) Multi-agent deep reinforcement learning: a survey. Artificial Intelligence Review 55 (2), p. 895–943. Cited by: §1. J. Ho, A. Jain, and P. Abbeel (2020) Denoising diffusion probabilistic models. In Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (Eds.), Vol. 33, p. 6840–6851. Cited by: §B.2, §1, §2.2. V. Jain, T. Akhound-Sadegh, and S. Ravanbakhsh (2024) Sampling from energy-based policies using diffusion. arXiv preprint arXiv:2410.01312. Cited by: §1. H. Kim, G. Papamakarios, and A. Mnih (2021) The lipschitz constant of self-attention. In International Conference on Machine Learning, p. 5562–5571. Cited by: §B.2. [15] G. Li and Y. Yan O (d/t) convergence theory for diffusion probabilistic models under minimal assumptions. In The Thirteenth International Conference on Learning Representations, Cited by: §B.2, §B.2. Y. Lipman, R. T. Chen, H. Ben-Hamu, M. Nickel, and M. Le (2022) Flow matching for generative modeling. arXiv preprint arXiv:2210.02747. Cited by: §1. J. Liu, Y. Zhong, S. Hu, H. Fu, Q. Fu, X. Chang, and Y. Yang (2024) Maximum entropy heterogeneous-agent reinforcement learning. In ICLR, Cited by: §1, §1, §5. H. Ma, T. Chen, K. Wang, N. Li, and B. Dai (2025a) Efficient online reinforcement learning for diffusion policy. arXiv preprint arXiv:2502.00361. Cited by: §1. H. Ma, T. Chen, K. Wang, N. Li, and B. Dai (2025b) Efficient online reinforcement learning for diffusion policy. In Proceedings of the 42nd International Conference on Machine Learning, A. Singh, M. Fazel, D. Hsu, S. Lacoste-Julien, F. Berkenkamp, T. Maharaj, K. Wagstaff, and J. Zhu (Eds.), Proceedings of Machine Learning Research, Vol. 267, p. 41837–41853. Cited by: §2.2. L. Matignon, G. J. Laurent, and N. Le Fort-Piat (2012) Independent reinforcement learners in cooperative markov games: a survey regarding coordination problems. The Knowledge Engineering Review 27 (1), p. 1–31. Cited by: §2.1. J. Orr and A. Dutta (2023) Multi-agent deep reinforcement learning for multi-robot applications: a survey. Sensors 23 (7), p. 3625. Cited by: §1. M. Psenka, A. Escontrela, P. Abbeel, and Y. Ma (2023) Learning a diffusion model policy from rewards via q-score matching. arXiv preprint arXiv:2312.11752. Cited by: §1. R. Rombach, A. Blattmann, D. Lorenz, P. Esser, and B. Ommer (2022) High-resolution image synthesis with latent diffusion models. In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, p. 10684–10695. Cited by: §1. A. Rutherford, B. Ellis, M. Gallici, J. Cook, A. Lupu, G. Ingvarsson, T. Willi, R. Hammond, A. Khan, C. S. de Witt, et al. (2024) Jaxmarl: multi-agent rl environments and algorithms in jax. Advances in Neural Information Processing Systems 37, p. 50925–50951. Cited by: 3rd item. R. Sahitaj, P. Sasnauskas, Y. Yalın, D. Mandal, and G. Radanović (2025) Independent learning in performative markov potential games. arXiv preprint arXiv:2504.20593. Cited by: §1. Y. Song, C. Durkan, I. Murray, and S. Ermon (2021) Maximum likelihood training of score-based diffusion models. Advances in neural information processing systems 34, p. 1415–1428. Cited by: §1. Y. Song and D. P. Kingma (2021) How to train your energy-based models. arXiv preprint arXiv:2101.03288. Cited by: §1. Y. Sun, T. Liu, R. Zhou, P. Kumar, and S. Shahrampour (2023) Provably fast convergence of independent natural policy gradient for markov potential games. Advances in neural information processing systems 36, p. 43951–43971. Cited by: §1. S. Vempala and A. Wibisono (2019) Rapid convergence of the unadjusted langevin algorithm: isoperimetry suffices. Advances in neural information processing systems 32. Cited by: §B.2.2, §B.2.2. A. Virmaux and K. Scaman (2018) Lipschitz regularity of deep neural networks: analysis and efficient estimation. Advances in neural information processing systems 31. Cited by: §B.2. Y. Wang, L. Wang, Y. Jiang, W. Zou, T. Liu, X. Song, W. Wang, L. Xiao, J. Wu, J. Duan, et al. (2024) Diffusion actor-critic with entropy regulator. Advances in Neural Information Processing Systems 37, p. 54183–54204. Cited by: §1, §5, §5. E. Yang and D. Gu (2004) Multiagent reinforcement learning for multi-robot systems: a survey. Technical report tech. rep. Cited by: §1. L. Yang, Z. Huang, F. Lei, Y. Zhong, Y. Yang, C. Fang, S. Wen, B. Zhou, and Z. Lin (2023) Policy representation via diffusion probability model for reinforcement learning. arXiv preprint arXiv:2305.13122. Cited by: §1. T. Yang, S. Cen, Y. Wei, Y. Chen, and Y. Chi (2024) Federated natural policy gradient and actor critic methods for multi-task reinforcement learning. Advances in Neural Information Processing Systems 37, p. 121304–121375. Cited by: §1. C. Yu, A. Velu, E. Vinitsky, J. Gao, Y. Wang, A. Bayen, and Y. Wu (2022) The surprising effectiveness of ppo in cooperative multi-agent games. Advances in neural information processing systems 35, p. 24611–24624. Cited by: §5. L. Yu, Y. Sun, Z. Xu, C. Shen, D. Yue, T. Jiang, and X. Guan (2020) Multi-agent deep reinforcement learning for hvac control in commercial buildings. IEEE Transactions on Smart Grid 12 (1), p. 407–419. Cited by: §1. R. Zhang, J. Mei, B. Dai, D. Schuurmans, and N. Li (2022) On the global convergence rates of decentralized softmax gradient play in markov potential games. Advances in Neural Information Processing Systems 35, p. 1923–1935. Cited by: Appendix C, Appendix C, Theorem C.1, §1, §3, Algorithm 1. Y. Zhang, S. Yu, T. Zhang, M. Guang, H. Hui, K. Long, Y. Wang, C. Yu, and W. Ding (2025) SAC flow: sample-efficient reinforcement learning of flow-based policies via velocity-reparameterized sequential modeling. arXiv preprint arXiv:2509.25756. Cited by: §1. Y. Zhong, J. G. Kuba, X. Feng, S. Hu, J. Ji, and Y. Yang (2024) Heterogeneous-agent reinforcement learning. Journal of Machine Learning Research 25 (32), p. 1–67. Cited by: §1, §1. Z. Zhu, F. Locatello, and V. Cevher (2023) Sample complexity bounds for score-matching: causal discovery and generative modeling. Advances in Neural Information Processing Systems 36, p. 3325–3337. Cited by: §B.2. Appendix A Theoretical Results for Section˜3.2 Proof of Theorem˜3.1. For simplicity, we omit all subscripts i in this proof. Step 1. We first calculate the KL divergence between any Gaussian policy π(⋅|0)=(μ,σ2)π(·| s_0)=N(μ,σ^2) and π~k+1(⋅|0)=πk(⋅|0)exp(η1−γQ^πk(0,⋅))/Z π^k+1(·| s_0)=π^k(·| s_0) ( η1-γ Q^π^k( s_0,·) )/Z, where Z is the normalizing factor. ℒ(μ,σ2)≔KL(π(⋅|0)∥π~k+1(⋅|0))=∼π(⋅|0)(logπ(|0)π~k+1(|0))=∼π(⋅|0)[logπ(|0)πk(|0)]−η1−γ∼π(⋅|0)[Q^πk(0,)]+logZ=log(σkσ)+σ2+(μ−μk)22(σk)2−12−η1−γ∼π(⋅|0)[Q^πk(0,)]+logZ. splitL(μ,σ^2) &KL (π(·| s_0)\| π^k+1(·| s_0) )=E_ a π(·| s_0) ( π( a| s_0) π^k+1( a| s_0) )\\ =&E_ a π(·| s_0) [ π( a| s_0)π^k( a| s_0) ]- η1-γE_ a π(·| s_0) [ Q^π^k( s_0, a) ]+ Z\\ =& ( σ^kσ )+ σ^2+ (μ-μ^k )^22(σ^k)^2- 12- η1-γE_ a π(·| s_0) [ Q^π^k( s_0, a) ]+ Z. split (3) Since the environment does not reach 1 s_1 in epochs [k][k], and Q^1π1k Q_1^π^k_1 is accurate for ℳsaM_sa with reward rsa¬(,)=(=0)rstart()r _sa( s, a)=I ( s= s_0 )r_start( a), we have Q^πk(0,)=rstart() Q^π^k( s_0, a)=r_start( a). We can then derive the following closed-form solution for the second last term: ∼π(⋅|0)[Q^πk(0,)]=∼π(⋅|0)[−ν2+c0exp(−α(+β)2)+c0exp(−α(−β)2)]=−ν(μ2+σ2)+c01+2ασ2exp(−α(μ−β)21+2ασ2)+c01+2ασ2exp(−α(μ+β)21+2ασ2). split&E_ a π(·| s_0) [ Q^π^k( s_0, a) ]\\ =&E_ a π(·| s_0) [-ν a^2+c_0 (-α( a+β)^2 )+c_0 (-α( a-β)^2 ) ]\\ =&-ν (μ^2+σ^2 )+ c_0 1+2ασ^2 (- α(μ-β)^21+2ασ^2 )+ c_0 1+2ασ^2 (- α(μ+β)^21+2ασ^2 ). split Step 2. The updated mean and covariance pair (μk+1,(σk+1)2)(μ^k+1,(σ^k+1)^2) is the first order stationary point of loss function ℒ(μ,σ2)L(μ,σ^2). To analyze μk+1μ^k+1, we first derive the gradient with respect to μ and evaluate it at μ=μk=0μ=μ^k=0. ∂μℒ(μ,σ2)=μ−μk(σk)2−η1−γ∂μ∼π(⋅|0)[Q^πk(0,)]. split& ∂μL(μ,σ^2)= μ-μ^k(σ^k)^2- η1-γ ∂μE_ a π(·| s_0) [ Q^π^k( s_0, a) ]. split (4) For the last term, we have: ∂μ∼π(⋅|0)[Q^πk(0,)]|μ=μk=−2νμ+2c01+2ασ2α(β−μ)1+2ασ2exp(−α(μ−β)21+2ασ2)+2c01+2ασ2α(−β−μ)1+2ασ2exp(−α(μ+β)21+2ασ2)|μ=μk=0. split& . ∂μE_ a π(·| s_0) [ Q^π^k( s_0, a) ] |_μ=μ^k\\ =&-2νμ+ 2c_0 1+2ασ^2 α (β-μ )1+2ασ^2 (- α(μ-β)^21+2ασ^2 ) .+ 2c_0 1+2ασ^2 α (-β-μ )1+2ασ^2 (- α(μ+β)^21+2ασ^2 ) |_μ=μ^k\\ =&0. split Substituting back, we have: ∂μℒ(μ,σ2)|μ=μk=0. split& . ∂μL(μ,σ^2) |_μ=μ^k=0. split (5) Therefore, starting from μk=0μ^k=0, we have μk+1=0μ^k+1=0. Step 3. Next, we derive ∂ℒ(μ,σ2)/∂σ2 (μ,σ^2)/∂σ^2: ∂σ2ℒ(μ,σ2)=−12σ2+12(σk)2−η1−γ∂σ2∼π(⋅|0)[Q^πk(0,)]. split ∂σ^2L(μ,σ^2)=&- 12σ^2+ 12(σ^k)^2- η1-γ ∂σ^2E_ a π(·| s_0) [ Q^π^k( s_0, a) ]. split (6) For the last term, we have: ∂σ2∼π(⋅|0)[Q^πk(0,)]=−ν+c0α(2α(μ−β)2−1−2ασ2)(1+2ασ2)5/2exp(−α(μ−β)21+2ασ2)+c0α(2α(μ+β)2−1−2ασ2)(1+2ασ2)5/2exp(−α(μ+β)21+2ασ2). split ∂σ^2E_ a π(·| s_0) [ Q^π^k( s_0, a) ]=&-ν+c_0 α (2α(μ-β)^2-1-2ασ^2 ) (1+2ασ^2 )^5/2 (- α(μ-β)^21+2ασ^2 )\\ &+c_0 α (2α(μ+β)^2-1-2ασ^2 ) (1+2ασ^2 )^5/2 (- α(μ+β)^21+2ασ^2 ). split (7) Since μ stays at 0, the above simplifies to ∂σ2∼π(⋅|0)[Q^πk(0,)]=−ν+2c0α(2αβ2−1−2ασ2)(1+2ασ2)5/2exp(−αβ21+2ασ2) split& ∂σ^2E_ a π(·| s_0) [ Q^π^k( s_0, a) ]=-ν\\ &+2c_0 α (2αβ^2-1-2ασ^2 ) (1+2ασ^2 )^5/2 (- αβ^21+2ασ^2 )\\ split (8) Substituting back into Equation˜6 and setting the gradient to 0, we have: 12σ2=12(σk)2−η1−γ∂σ2∼π(⋅|0)[Q^πk(0,)]=12(σk)2+η1−γ(ν−2c0α(2αβ2−1−2ασ2)(1+2ασ2)5/2exp(−αβ21+2ασ2)). split 12σ^2=& 12(σ^k)^2- η1-γ ∂σ^2E_ a π(·| s_0) [ Q^π^k( s_0, a) ]\\ =& 12(σ^k)^2+ η1-γ (ν-2c_0 α (2αβ^2-1-2ασ^2 ) (1+2ασ^2 )^5/2 (- αβ^21+2ασ^2 ) ). split (9) With η=0.01(1−γ),γ=0.1,ν=10,c0=2.7,α=250η=0.01(1-γ),γ=0.1,ν=10,c_0=2.7,α=250, and β=0.5β=0.5, we have 1σ2=1(σk)2+0.02(10−1350124−500σ2(1+500σ2)5/2⋅exp(−62.51+500σ2))≥1(σk)2+0.02(10−3.6)≥1(σk)2+0.128. split 1σ^2=& 1(σ^k)^2+0.02 (10-1350 124-500σ^2 (1+500σ^2 )^5/2· (- 62.51+500σ^2 ) )\\ ≥& 1(σ^k)^2+0.02 (10-3.6 )≥ 1(σ^k)^2+0.128. split (10) Since the above inequality holds for all stationary points σ, we know that 1(σk)2≥1(σk−1)2+0.128≥⋯≥1(σ0)2+0.128k. split 1 (σ^k )^2≥ 1 (σ^k-1 )^2+0.128≥·s≥ 1 (σ^0 )^2+0.128k. split (11) Step 4. Finally, we lower bound the probability that the algorithm never reaches 1 s_1. For any epoch k before 1 s_1 is reached for the first time, consider policy πk=(0,(σk)2)π^k=N(0,(σ^k)^2) with (σk)2(σ^k)^2 satisfying Equation˜11. We upper bound the probability of reaching state 1 s_1 with a single action sample from πkπ^k, denoted by pkp_k. pk=∫−β−14−β+1412π(σk)2exp(−22(σk)2)⋅cpexp(−α(+β)2)d=(i)C′2π(σ′)2∫−β−14−β+14exp(−(−μ′)22(σ′)2)d≤C′Φ(−β+14−μ′σ′) splitp_k=& _-β- 14^-β+ 14 1 2π(σ^k)^2 (- a^22 (σ^k )^2 )· c_p (-α( a+β)^2 )d a\\ (i)=& C 2π(σ )^2 _-β- 14^-β+ 14 (- ( a-μ )^22 (σ )^2 )d a\\ ≤&C ( -β+ 14-μ σ )\\ split (12) In (i)(i) we have defined μ′=−2αβ(σk)21+2α(σk)2,(σ′)2=(σk)21+2α(σk)2,C′=cp1+2α(σk)2exp(−αβ21+2α(σk)2). splitμ =&- 2αβ(σ^k)^21+2α(σ^k)^2, (σ )^2= (σ^k )^21+2α(σ^k)^2, C = c_p 1+2α(σ^k)^2 (- αβ^21+2α(σ^k)^2 ). split (13) For k≥150000k≥ 150000, we have that (σk)2≤1/(10+0.128k)≤1/19210(σ^k)^2≤ 1/(10+0.128k)≤ 1/19210 and that −μ′≤0.05-μ ≤ 0.05. This leads to pk≤C′Φ(−15σ′)≤C′Φ(−151+2α(σk)2(σk)2)≤C′Φ(−150.128k+510)≤5C′exp(−(0.00256k+10.2))2π(0.128k+510). splitp_k≤&C (- 15σ )≤ C (- 15 1+2α (σ^k )^2 (σ^k )^2 )≤ C (- 15 0.128k+510 )\\ ≤& 5C (- (0.00256k+10.2 ) ) 2π (0.128k+510 ). split (14) In the last line, we have used Φ(−x)≤exp(−x2/2)/(2πx) (-x)≤ (-x^2/2)/( 2πx). This is because Φ(−x)=12π∫x∞exp(−t22)dt≤12π∫x∞txexp(−t22)dt=12πxexp(−t2/2)|∞x=12πxexp(−x2/2) (-x)= 1 2π _x^∞ (- t^22 )dt≤ 1 2π _x^∞ tx (- t^22 )dt= 1 2πx (-t^2/2 )|^x_∞= 1 2πx (-x^2/2 ) for x>0x>0. We now use the above inequality to lower bound the probability of not reaching state 1 s_1 in all epochs, denoted by p¯∞ p_∞. p¯∞≥∏k=0∞(1−pk)32=∏k=0∞exp(ln(1−pk))32≥(i)exp(−32∑k=0∞pk1−pk)=exp(−32∑k=0149999pk1−pk)exp(−32∑k=150000∞pk1−pk)≥exp(−32∑k=0149999pk1−pk)exp(−320.9999∑k=150000∞pk) split p_∞≥& _k=0^∞ (1-p_k )^32= _k=0^∞ ( (1-p_k ) )^32 (i)≥ (-32 _k=0^∞ p_k1-p_k )\\ =& (-32 _k=0^149999 p_k1-p_k ) (-32 _k=150000^∞ p_k1-p_k )\\ ≥& (-32 _k=0^149999 p_k1-p_k ) (- 320.9999 _k=150000^∞p_k )\\ split (15) In (i)(i), we have used ln(1−x)≥−x/(1−x) (1-x)≥-x/(1-x) for x∈(0,1)x∈(0,1). We bound ∑k=150000∞pk _k=150000^∞p_k as follows. ∑k=150000∞pk≤∫149999∞5C′exp(−(0.00256x+10.2))2π(0.128x+510)dx≤5C′exp(−10.2)1020π∫149999∞exp(−0.00256x)dx≤10−100. split _k=150000^∞p_k≤& _149999^∞ 5C (- (0.00256x+10.2 ) ) 2π (0.128x+510 )dx\\ ≤& 5C (-10.2 ) 1020π _149999^∞ (-0.00256x )dx≤ 10^-100. split (16) Therefore, p¯∞≥exp(−32∑k=0149999pk1−pk)exp(−320.9999∑k=150000∞pk)≥0.55. split p_∞≥& (-32 _k=0^149999 p_k1-p_k ) (- 320.9999 _k=150000^∞p_k )≥ 0.55. split (17) ∎ Proof of Proposition 3.2. We let EiE_i be the event that []i=1[ s]_i=1 at least once when executing Algorithm˜1, and let kik_i denote the first epoch when []i=1[ s]_i=1. Then ℙ(visit =n at least once)≤ℙ(⋂i∈[n]Ei). split&P (visit s= 1_n at least once ) ( _i∈[n]E_i ). split (18) Notice that event EicE_i^c, i.e., the agents never visit []i=1[ s]_i=1, only depends on the individual policy πik(⋅|start) _i^k(·| s_start) for epoch k≤kik≤ k_i. Moreover, the individual policy πik _i^k with k≤kik≤ k_i can be expressed recursively as πik(⋅|start)=argminπiKL(πi∥π~ik),π~ik(⋅)∝πik−1(⋅|0)exp(ηQ^iπk−1(0,⋅)(1−γ)). split& _i^k(·| s_start)= _ _iKL ( _i\| π_i^k ), π_i^k(·) _i^k-1(·| s_0) ( η Q_i^π^k-1( s_0,·)(1-γ) ). split (19) Since Q^iπk−1 Q_i^π^k-1 is accurate for ℳmaM_ma reward function rma¬(,)=(=−n)∑j=1nrstart(j)+(≠−n)c1∑j≠i[]jr_ma ( s, a)=I ( s=- 1_n ) _j=1^nr_start( a_j)+I ( s≠- 1_n )\,c_1 _j≠ i[ s]_j, policy πi _i only contributes through term rstart(j)r_start( a_j). Therefore, Q^iπk−1(0,⋅)=rstart(i)+g(π−i) Q_i^π^k-1( s_0,·)=r_start( a_i)+g( _-i) for some function g. Substituting back gives: π~ik(⋅)∝πik−1(⋅|0)exp(ηQ^iπk−1(0,⋅)(1−γ))∝πik−1(⋅|0)exp(ηrstart(i)(1−γ)). split& π_i^k(·) _i^k-1(·| s_0) ( η Q_i^π^k-1( s_0,·)(1-γ) ) _i^k-1(·| s_0) ( η r_start( a_i)(1-γ) ). split (20) Therefore, the policy update of agent i is independent of other agents for k≤kik≤ k_i. Thus, events Eici∈[n]\E_i^c\_i∈[n] remain independent of each other, and ℙ(Eic)≥0.55P(E_i^c)≥ 0.55 by Theorem˜3.1. As a result, events Eii∈[n]\E_i\_i∈[n] are also independent and ℙ(Ei)≤0.45P(E_i)≤ 0.45 for all i∈[n]i∈[n]. This gives ℙ(visit =n at least once)≤ℙ(⋂i∈[n]Ei)=∏i∈[n]ℙ(Ei)≤0.45n. split&P (visit s= 1_n at least once ) ( _i∈[n]E_i )= _i∈[n]P (E_i )≤ 0.45^n. split (21) ∎ Appendix B Theoretical Results for Section˜4 B.1 Importance Sampling Score Matching Loss To show the equivalence, we first rewrite the importance sampling loss as follows ℒi,is(θi)≔1N∑j=1Nπi⋆,k(i(j)|)p~(i(j)|)liθi(i(j)|),whereliθi(i(j)|)≔∑τ=1Hβτ2στ2ατϵτ∼(0,I)∥si,τθi(α¯τi(j)+1−α¯τϵτ|)−ϵτ1−α¯τ∥2. split&L_i, s^is( _i) 1N _j=1^N _i ,k( a_i^(j)| s) p( a_i^(j)| s)l _i_i( a_i^(j)| s),\\ &where l _i_i( a_i^(j)| s) _τ=1^H _τ^2 _τ^2 _τE_ _τ (0,I) \|s_i,τ _i ( . α_τ a_i^(j)+ 1- α_τ _τ | s )- _τ 1- α_τ \|^2. split (22) In expectation, ℒi,is(θi)L_i, s^is( _i) is equivalent to the original denoising score matching loss ℒi,(θi)L_i, s( _i) because i(j)∼p~(⋅|)[ℒi,is(θi)]=1N∑j=1N∫p~(i(j)|)πi⋆,k(i(j)|)p~(i(j)|)⋅liθi(i(j)|)di(j)=∫πi⋆,k(i|)liθi(i|)di=i∼πi⋆,k(⋅|)[liθi(i|)]=ℒi,(θi). splitE_ a_i^(j) p(·| s) [L_i, s^is( _i) ]=& 1N _j=1^N p( a_i^(j)| s) _i ,k( a_i^(j)| s) p( a_i^(j)| s)· l _i_i( a_i^(j)| s)d a_i^(j)\\ =& _i ,k( a_i| s)l _i_i( a_i| s)d a_i=E_ a_i _i ,k(·| s) [l _i_i( a_i| s) ]\\ =&L_i, s( _i). split (23) B.2 Proof of Theorem˜4.1 The variance schedule considered in this theorem is commonly studied in diffusion model literature [Li and Yan, ], and is defined for sufficiently large constants c′,c′c ,c : β1=1Hc′,βτ+1=c′logHminβ1(1+c′logH)τ,1,ατ=1−βτ,α¯τ=∏h=1ταh,στ2=βτατ. split& _1= 1H^c , _τ+1= c HH \ _1 (1+ c HH )^τ,1 \,\\ & _τ=1- _τ, α_τ= _h=1^τ _h, _τ^2= _τ _τ. split (24) To proceed, we first introduce the following assumptions on the score function class Θ≔(τ,)↦sτθ():τ∈[H],∈i,θ∈ΘS_ \(τ, a) s^θ_τ( a):τ∈[H], a _i,θ∈ \. Assumption B.1. (i). Expressiveness. There exists θ⋆∈Θθ ∈ such that sτθ⋆(τ)=∇τpτ(τ)s^θ _τ( a_τ)= _ a_τp_τ( a_τ). (i). Complexity. The κ-covering number of ΘS_ 111The κ-covering number of function class ΘS_ is the cardinality of the smallest parameter subset Θ~⊂Θ ⊂ such that for every θ∈Θθ∈ , there exists θ~∈Θ~ θ∈ such that ‖sτθ()−sτθ~()‖≤κ\|s^θ_τ( a)-s θ_τ( a)\|≤κ for all τ∈[H]τ∈[H] and a in the support. , denoted by (Θ,κ)C(S_ ,κ), satisfies log(Θ,κ)≤c1log(c1/κ) (S_ ,κ)≤ c_1 (c_1/κ) for constant c1≥1c_1≥ 1 and all κ>0κ>0. (i). Linear growth. ‖sτθ()‖≤c2(1+‖) \|s_τ^θ( a) \|≤ c_2(1+ \| a \|) for constant c2≥1c_2≥ 1, all θ∈Θθ∈ , τ∈[H]τ∈[H], and ∈ℝd a ^d. Section˜B.2 is mild and standard in function approximation literature. Part (i) on expressiveness is satisfied by common neural network architectures, which are universal function approximators. Parts (i) and (i) are satisfied by various popular neural network architectures including feedforward networks, convolutional networks, and transformers, where (i) is established in related literature ([Bartlett et al., 2017, Edelman et al., 2022]) and (i) follows by their Lipschitzness ([Virmaux and Scaman, 2018, Kim et al., 2021]). Moreover, the true score function sτθ⋆(τ)s^θ _τ( a_τ) is often assumed to be Lipschitz in the rich diffusion literature [Chen et al., 2022, Zhu et al., 2023], implying the linear growth of the true score function. This further justifies the validity of part (i). We also note that Part (i) is mainly for simplicity of the proof. We can easily extend it to ϵε-realizability, which assumes the existence of θ⋆∈Θθ ∈ such that sτθ⋆(τ)s^θ _τ( a_τ) is only ϵε-accurate. This will lead to an additional error term of size ϵε in the final squared TV distance bound. Proof. Recall that the target distribution is πi⋆,k(⋅|) _i ,k(·| s), the true score function is si,τθ⋆(⋅|)s_i,τ^θ (·| s), and the sampling distribution is p~(⋅|) p(·| s). The learned action distribution with or without projection onto the bounded action space iA_i are denoted by πiθ^i(⋅|)π θ_i_i(·| s) or πiθ^i,noclip(⋅|) _i θ_i,noclip(·| s), respectively. For simplicity, we omit subscript i and the dependence on state s and epoch k, reducing the above notations to π⋆(⋅)π (·), sτθ⋆(⋅)s_τ^θ (·), p~(⋅) p(·), πθ^(⋅)π θ(·), and πθ^,noclip(⋅)π θ,noclip(·). Moreover, we define importance weight w()w( a), per-sample loss lθ()l^θ( a), and per-sample weighted loss fθ()f^θ( a) as follows. All are supported on i⊆:‖≤cbA_i \ a: \| a \|≤ c_b\. fθ()≔w()lθ(),w()≔π⋆()p~(),lθ()≔∑τ=1Hβτ2στ2ατϵτ∼(0,I)‖sτθ(α¯τ+1−α¯τϵτ)−ϵτ1−α¯τ‖2=∑τ=1Hβτϵτ∼(0,I)‖sτθ(α¯τ+1−α¯τϵτ)−ϵτ1−α¯τ‖2 splitf^θ( a) &w( a)l^θ( a), w( a) π ( a) p( a),\\ l^θ( a) & _τ=1^H _τ^2 _τ^2 _τE_ _τ (0,I) \|s_τ^θ ( α_τ a+ 1- α_τ _τ )- _τ 1- α_τ \|^2\\ =& _τ=1^H _τE_ _τ (0,I) \|s_τ^θ ( α_τ a+ 1- α_τ _τ )- _τ 1- α_τ \|^2\\ split (25) The loss can then be defined as follows for (j)∼iidp~j∈[N]\ a^(j) iid p\_j∈[N]. ℒ(θ)≔∼p~[fθ()]=∼π⋆[lθ()],ℒ^(θ)≔1N∑j=1Nfθ((j)). split&L(θ) _ a p[f^θ( a)]=E_ a π [l^θ( a)], L(θ) 1N _j=1^Nf^θ( a^(j)). split (26) We will analyze the solution to the optimization problem θ^≔argminθ∈Θℒ^(θ) θ _θ∈ L(θ) and the resulting projected policy πθ^(⋅)π θ(·). The proof proceeds in two steps. In Step 1, we establish a concentration bound on ‖ℒ^(θ)−ℒ(θ)‖\| L(θ)-L(θ)\| for all θ∈Θθ∈ , and utilize the above concentration to upper bound the population loss difference ℒ(θ^)−ℒ(θ⋆)L( θ)-L(θ ). In Step 2, we connect this loss difference to the TV distance DTV(π⋆(⋅),πθ^(⋅))D_TV(π (·),π θ(·)). Step 1: Concentration of ℒ^(θ)−ℒ(θ) L(θ)-L(θ). We start from a 1/N1/ N-cover of function class ℱΘ=fθ(⋅):θ∈ΘF_ =\f^θ(·):θ∈ \, denoted by ℱΘ~F_ . By Lemma˜B.2, we know that |ℱΘ~|≤(3(c′+c′)c1c2cw(1+cb+2d)NlogH)c1|F_ |≤(3(c +c )c_1c_2c_w(1+c_b+2 d) N H)^c_1. We first bound the mean and variance of fθ~()f θ( a) for any θ~∈Θ~ θ∈ and ∈¯ a∈ A. For its mean, we have the following bound by Lemma˜B.3, ‖fθ~()−fθ~()‖≤max∈¯‖fθ~()‖≤cwmax∈¯‖lθ~()‖≤3(c′+c′)c22(1+cb+2d)2logH⏟κ1cw. split\|f θ( a)-Ef θ( a)\|≤& _ a∈ A\|f θ( a)\|≤ c_w _ a∈ A\|l θ( a)\|\\ ≤& 3(c +c )c_2^2(1+c_b+2 d)^2 H_ _1c_w. split (27) Here the first inequality is because fθ~()>0f θ( a)>0. Moreover, for its variance, we have ∼p~(‖fθ~()−fθ~()‖2)≤∼p~(‖fθ~()‖2)≤(i)∼p~(w()4)∼p~(lθ~()4)≤(ii)κ12∼p~(w()4). splitE_ a p ( \|f θ( a)-Ef θ( a) \|^2 )≤&E_ a p ( \|f θ( a) \|^2 ) (i)≤ E_ a p(w( a)^4) E_ a p(l θ( a)^4)\\ (i)≤& _1^2 E_ a p(w( a)^4). split (28) Here (i)(i) is by Cauchy-Schwarz and (ii)(i) is by Lemma˜B.3. Therefore, by the Bernstein’s inequality, we know that with probability at least 1−δ′1-δ for any δ′∈(0,1/|ℱΘ~|)δ ∈(0,1/|F_ |), ‖1N∑j=1Nfθ~((j))−fθ~()‖≤2κ12log(2/δ′)∼p~(w()4)4N⏟ℐ1+2κ1cwlog(2/δ′)31N⏟ℐ2. split \| 1N _j=1^Nf θ( a^(j))-Ef θ( a) \|≤& 2 _1^2 (2/δ ) [4]E_ a p(w( a)^4) N_I_1+ 2 _1c_w (2/δ )3 1N_I_2. split (29) Since N≥2cw2log(2/δ′)/9N≥ 2c_w^2 (2/δ )/9 and (w()4)≥(w())4≥1E(w( a)^4) (w( a))^4≥ 1, we know that ℐ1≥ℐ2I_1 _2. Therefore, ‖1N∑j=1Nfθ~((j))−fθ~()‖≤2ℐ1=22κ1log2δ′∼p~(w()4)4N. split \| 1N _j=1^Nf θ( a^(j))-Ef θ( a) \|≤ 2I_1=&2 2 _1 2δ [4]E_ a p(w( a)^4) N. split (30) With a union bound, Equation˜30 holds for all θ~∈Θ~ θ∈ with probability at least 1−|ℱΘ~|δ′=1−δ1-|F_ |δ =1-δ. Here we have let δ=|ℱΘ~|δ′∈(0,1)δ=|F_ |δ ∈(0,1). Finally, we extend the above bound on θ~∈Θ~ θ∈ to any θ∈Θθ∈ . By the definition of 1/N1/ N cover ℱΘ~F_ , for any θ∈Θθ∈ , we can find θ~∈Θ~ θ∈ such that ‖fθ()−fθ~()‖≤1/N,∀∈A¯ \|f^θ( a)-f θ( a) \|≤ 1/ N,~∀ a∈ A. Therefore, the following holds for any θ∈Θθ∈ with probability at least 1−δ1-δ, ‖ℒ^(θ)−ℒ(θ)‖= \| L(θ)-L(θ) \|= ‖1N∑j=1Nfθ((j))−fθ()‖≤‖1N∑j=1Nfθ~((j))−fθ~()‖+2N \| 1N _j=1^Nf^θ( a^(j))-Ef^θ( a) \|≤ \| 1N _j=1^Nf θ( a^(j))-Ef θ( a) \|+ 2 N ≤ ≤ (22κ1log2|ℱΘ~|δ+2)∼p~(w()4)4N. (2 2 _1 2|F_ |δ+2 ) [4]E_ a p(w( a)^4) N. By definition of θ θ, we know that ℒ^(θ^)−ℒ^(θ⋆)≤0 L( θ)- L(θ )≤ 0. Therefore, with probability at least 1−δ1-δ, ℒ(θ^)−ℒ(θ⋆)≤ℒ(θ^)−ℒ^(θ^)+ℒ^(θ^)−ℒ^(θ⋆)+ℒ^(θ⋆)−ℒ(θ⋆)≤ℒ(θ^)−ℒ^(θ^)+ℒ^(θ⋆)−ℒ(θ⋆)≤2maxθ∈Θ‖ℒ(θ)−ℒ^(θ)‖≤(42κ1log2|ℱΘ~|δ+4)∼p~(w()4)4N. splitL( θ)-L(θ )≤&L( θ)- L( θ)+ L( θ)- L(θ )+ L(θ )-L(θ )\\ ≤&L( θ)- L( θ)+ L(θ )-L(θ )≤ 2 _θ∈ \|L(θ)- L(θ) \|\\ ≤& (4 2 _1 2|F_ |δ+4 ) [4]E_ a p(w( a)^4) N. split (31) Recall that ℒmarginal(θ)=∑τ=1Hβτ2στ2ατ(τ)‖sτθ((τ))−sτθ⋆((τ))‖2=∑τ=1H(1−ατ)(τ)‖sτθ((τ))−sτθ⋆((τ))‖2. splitL^marginal(θ)=& _τ=1^H _τ^2 _τ^2 _τE_ a(τ) \|s_τ^θ ( a(τ) )-s_τ^θ ( a(τ) ) \|^2\\ =& _τ=1^H(1- _τ)E_ a(τ) \|s_τ^θ ( a(τ) )-s_τ^θ ( a(τ) ) \|^2. split (32) By Proposition 3.1 in Ho et al. [2020], we know that ℒmarginal(θ)=ℒ(θ)+C3L^marginal(θ)=L(θ)+C_3 for constant C3C_3 independent of θ. Therefore, with probability at least 1−δ1-δ, ℒmarginal(θ^)≤ℒmarginal(θ^)−ℒmarginal(θ⋆)=ℒ(θ^)−ℒ(θ⋆)≤(42κ1log2|ℱΘ~|δ+4)∼p~(w()4)4N≤(42κ1c1log2κ1c1cwNδ+4)∼p~(w()4)4N. splitL^marginal( θ)≤&L^marginal( θ)-L^marginal(θ )=L( θ)-L(θ )\\ ≤& (4 2 _1 2|F_ |δ+4 ) [4]E_ a p(w( a)^4) N\\ ≤& (4 2 _1 c_1 2 _1c_1c_w Nδ+4 ) [4]E_ a p(w( a)^4) N. split (33) Here the first inequality is because ℒmarginal(θ⋆)=0L^marginal(θ )=0 by Section˜B.2 (i), and the last inequality is by |ℱΘ~|≤(3(c′+c′)c1c2cw(1+cb+2d)NlogH)c1≤(2κ1c1cwN/δ)c1|F_ |≤(3(c +c )c_1c_2c_w(1+c_b+2 d) N H)^c_1≤(2 _1c_1c_w N/δ)^c_1. Step 2: Connecting to DTV(π⋆(⋅),πθ^(⋅)) D_TV(π (·),π θ(·)). First, recall that πθ^,noclipπ θ,noclip and π⋆π are policies induced by the reverse diffusion denoising process at τ=0τ=0 with score functions sθ^s θ and sθ⋆s^θ , respectively. We further introduce policies π1θ^,noclipπ θ,noclip_1 and π1⋆π _1 to represent the policies induced by the diffusion denoising process at diffusion step τ=1τ=1. Combining Equation (20-23 (i)) and the last inequality in [Li and Yan, , Section 4], we have the following for some constant κ2≥1 _2≥ 1 DTV2(π1⋆,π1θ^,noclip)≤κ2d2log6H2+κ2∑τ=2H1−ατ2⋅(τ)[‖sτθ^((τ))−sτθ⋆((τ))‖2]. splitD_TV^2 (π _1,π θ,noclip_1 )≤& _2 d^2 ^6HH^2+ _2 _τ=2^H 1- _τ2·E_ a(τ) [ \|s_τ θ ( a(τ) )-s_τ^θ ( a(τ) ) \|^2 ]. split (34) Utilizing the one-step denoising TV bound Lemma˜B.5 with σ2=β1/α1,β=β1,α=α1σ^2= _1/ _1,β= _1,α= _1, DTV2(π0⋆,π0θ^,noclip)≤β12(1)[‖s1θ^((1))−s1⋆((1))‖2]+2DTV2(π1⋆,π1θ^,noclip)+β1c22dα1(α1/β1−c2)≤(i)3κ2d2log6H2+κ2∑τ=1H(1−ατ)⋅(τ)[‖sτθ^((τ))−sτθ⋆((τ))‖2]=3κ2d2log6H2+κ2ℒmarginal(θ^). splitD_TV^2 (π _0,π θ,noclip_0 )≤& _12E_ a(1) [ \|s_1 θ ( a(1) )-s_1 ( a(1) ) \|^2 ]+2D_TV^2 (π _1,π θ,noclip_1 )\\ & 20.00003pt+ _1c_2^2d _1 ( _1/ _1-c_2 )\\ (i)≤&3 _2 d^2 ^6HH^2+ _2 _τ=1^H (1- _τ )·E_ a(τ) [ \|s_τ θ ( a(τ) )-s_τ^θ ( a(τ) ) \|^2 ]\\ =&3 _2 d^2 ^6HH^2+ _2L^marginal( θ). split (35) Here (i)(i) is because c22β1d/(α1(α1/β1−c2))≤κ2d2log6H/H2c_2^2 _1d/( _1( _1/ _1-c_2))≤ _2d^2 ^6H/H^2 for sufficiently large c′c . Finally, we connect back to the projected policy π0θ _0 θ. By definition, applying the projection Proji()=argmin′∈i‖′−‖Proj_A_i( a)= _ a _i \| a - a \| to samples from π0θ^,noclip(⋅)π θ,noclip_0(·) yields π0θ^(⋅)π θ_0(·). Moreover, applying Proji(⋅)Proj_A_i(·) to π⋆(⋅)π (·) doesn’t change the sample distribution since π⋆(⋅)π (·) is only defined on iA_i. Therefore, by the data processing inequality, we have DTV(π0⋆,π0θ^)≤DTV(π0⋆,π0θ^,noclip)D_TV(π _0,π θ_0)≤ D_TV(π _0,π θ,noclip_0). Combining with Equation˜35 gives the following for N≥2c1cw2log(2κ1c1cwN/δ)/9≥2cw2log(2|ℱΘ~|/δ)/9N≥ 2c_1c_w^2 (2 _1c_1c_w N/δ)/9≥ 2c_w^2 (2|F_ |/δ)/9: DTV2(π0⋆,π0θ^)≤3κ2d2log6H2+κ2ℒmarginal(θ^)≤3κ2d2log6H2+κ2(42κ1c1log2κ1c1cwNδ+4)∼p~(w()4)4N≲d2log6H2+d∼p~(w()4)4N=d2log6H2+dexp(34D4(π⋆∥p~))N. split&D_TV^2 (π _0,π θ_0 )≤ 3 _2 d^2 ^6HH^2+ _2L^marginal( θ)\\ ≤&3 _2 d^2 ^6HH^2+ _2 (4 2 _1 c_1 2 _1c_1c_w Nδ+4 ) [4]E_ a p(w( a)^4) N\\ & d^2 ^6HH^2+ d [4]E_ a p(w( a)^4) N\\ =& d^2 ^6HH^2+ d ( 34D_4 (π \| p ) ) N. split (36) Here in the second last line, the dependencies on constants and polylog factors κ1=3(c′+c′)c22(1+cb+2d)2logH _1=3(c +c )c_2^2(1+c_b+2 d)^2 H, κ2 _2, c1log(2κ1c1cwN/δ) c_1 (2 _1c_1c_w N/δ) are absorbed using ≲ . In the last line, we used the definition of fourth-order Renyi divergence D4(p||q)=13lnx∼q[(p(x)/q(x))4]D_4(p||q)= 13 _x q [(p(x)/q(x))^4 ]. ∎ B.2.1 Supporting Details Lemma B.2 (Covering and boundedness of ℱΘF_ ). Consider the setting of Theorem˜4.1 and notations in Equation˜25. The κ-covering number of ℱΘ=fθ(⋅)=w(⋅)lθ(⋅):θ∈ΘF_ =\f^θ(·)=w(·)l^θ(·):θ∈ \, denoted by (ℱΘ,κ)C(F_ ,κ), satisfies log(ℱΘ,κ)≤c1log3(c′+c′)c1c2cw(1+cb+2d)logHκ. (F_ ,κ)≤ c_1 3(c +c )c_1c_2c_w(1+c_b+2 d) Hκ. (37) Proof. Consider any κ>0κ>0 and any κ-cover of ΘS_ , denoted by Θ~S_ . Then by definition, for any θ∈Θθ∈ , there exists θ~∈Θ~ θ∈ such that maxτ,‖sτθ()−sτθ~()‖≤κ _τ, a\|s_τ^θ ( a )-s_τ θ ( a )\|≤κ. We now bound ‖fθ()−fθ~()‖\|f^θ( a)-f θ( a)\| for any a in the support of f, i.e., :‖≤cb\ a: \| a \|≤ c_b\, as follows. ‖fθ()−fθ~()‖=‖w()lθ()−w()lθ~()‖≤∥cw∑τ=1Hβτϵτ∼(0,I)[∥sτθ(α¯τ+1−α¯τϵτ)−ϵτ1−α¯τ∥2−∥sτθ~(α¯τ+1−α¯τϵτ)−ϵτ1−α¯τ∥2]∥=cw∑τ=1Hβτϵτ∼(0,I)‖sτθ(α¯τ+1−α¯τϵτ)−sτθ~(α¯τ+1−α¯τϵτ)‖⋅‖sτθ(α¯τ+1−α¯τϵτ)+sτθ~(α¯τ+1−α¯τϵτ)−2ϵτ1−α¯τ‖≤cwκ⋅∑τ=1Hβτϵτ∼(0,I)[∥sτθ(α¯τ+1−α¯τϵτ)∥+∥sτθ~(α¯τ+1−α¯τϵτ)∥+2‖ϵτ‖1−α¯τ]≤(i)2cwκ⋅∑τ=1Hβτ(c2(1+cb+d)+d1−α¯τ)≤2c2cw(1+cb+2d)κ⋅∑τ=1Hβτ1−α¯τ≤(ii)3(c′+c′)c2cw(1+cb+2d)logH⋅κ. split& \|f^θ( a)-f θ( a) \|= \|w( a)l^θ( a)-w( a)l θ( a) \|\\ ≤& \|c_w _τ=1^H _τE_ _τ (0,I) [ \|s_τ^θ ( α_τ a+ 1- α_τ _τ )- _τ 1- α_τ \|^2 . .\\ & . .- \|s_τ θ ( α_τ a+ 1- α_τ _τ )- _τ 1- α_τ \|^2 ] \|\\ =&c_w _τ=1^H _τE_ _τ (0,I) \|s_τ^θ ( α_τ a+ 1- α_τ _τ )-s_τ θ ( α_τ a+ 1- α_τ _τ ) \|\\ &· \|s_τ^θ ( α_τ a+ 1- α_τ _τ )+s_τ θ ( α_τ a+ 1- α_τ _τ )- 2 _τ 1- α_τ \|\\ ≤&c_wκ· _τ=1^H _τE_ _τ (0,I) [ \|s_τ^θ ( α_τ a+ 1- α_τ _τ ) \| .\\ & .+ \|s_τ θ ( α_τ a+ 1- α_τ _τ ) \|+ 2 \| _τ \| 1- α_τ ]\\ (i)≤&2c_wκ· _τ=1^H _τ (c_2 (1+c_b+ d )+ d1- α_τ )\\ ≤&2c_2c_w (1+c_b+2 d )κ· _τ=1^H _τ1- α_τ\\ (i)≤&3(c +c )c_2c_w (1+c_b+2 d ) H·κ. split (38) Here (i)(i) is by Lemma˜B.3 and (ii)(i) is by Lemma˜B.4. Therefore, ℱΘ~F_ is a κ′κ -cover of ℱΘF_ with κ′=3(c′+c′)c2cw(1+cb+2d)logH⋅κ =3(c +c )c_2c_w (1+c_b+2 d ) H·κ. We can then conclude that the κ′κ -log covering number of ℱΘF_ is upper bounded by: log(ℱΘ,κ′)≤log|ℱΘ~|=c1logc1κ=c1log3(c′+c′)c1c2cw(1+cb+2d)logHκ′. split (F_ ,κ )≤ |F_ |=c_1 c_1κ=c_1 3(c +c )c_1c_2c_w(1+c_b+2 d) Hκ . split (39) ∎ Lemma B.3. Suppose H≥3c′logH≥ 3c H. Consider score function sτθ(⋅)s_τ^θ(·) satisfying Section˜B.2 (i) with constant c2c_2 and loss lθ(⋅)l^θ(·) defined in Equation˜25. They satisfy the following for all θ, τ∈[H]τ∈[H], and a with ‖≤cb \| a \|≤ c_b: ϵ∼(0,I)‖sτθ(α¯τ+1−α¯τϵ)‖2≤c22(1+cb+d)2, _ε (0,I) \|s_τ^θ ( α_τ a+ 1- α_τε ) \|^2≤ c_2^2 (1+c_b+ d )^2, (40a) ‖lθ()‖≤3(c′+c′)(c22(1+cb+d)2+d)logH. \|l^θ( a) \|≤ 3(c +c ) (c_2^2 (1+c_b+ d )^2+d ) H. (40b) Proof. We first prove Equation˜40a as follows: ϵ∼(0,I)‖sτθ(α¯τ+1−α¯τϵ)‖2≤(i)c22ϵ∼(0,I)[(1+α¯τ‖+1−α¯τ‖ϵ‖)2]≤(ii)c22(1+cbα¯τ+d1−α¯τ)2≤c22(1+cb+d)2. splitE_ε (0,I) \|s_τ^θ ( α_τ a+ 1- α_τε ) \|^2 (i)≤&c_2^2E_ε (0,I) [ (1+ α_τ \| a \|+ 1- α_τ \|ε \| )^2 ]\\ (i)≤&c_2^2 (1+c_b α_τ+ d 1- α_τ )^2\\ ≤&c_2^2 (1+c_b+ d )^2. split (41) Here (i)(i) is by Section˜B.2 (i); (ii)(i) is by expanding the square and using ‖ϵ‖≤‖ϵ‖2=dE \|ε \|≤ E\|ε\|^2= d. To prove Equation˜40b, recall the definition of lθ()l^θ( a): lθ()≔∑τ=1Hβτϵτ∼(0,I)‖sτθ(α¯τ+1−α¯τϵτ)−ϵτ1−α¯τ‖2. splitl^θ( a) _τ=1^H _τE_ _τ (0,I) \|s_τ^θ ( α_τ a+ 1- α_τ _τ )- _τ 1- α_τ \|^2. split (42) Therefore, ‖lθ()‖≤2∑τ=1Hβτ⋅ϵτ∼(0,I)[‖sτθ(α¯τ+1−α¯τϵτ)‖2+‖ϵτ‖21−α¯τ]≤(i)2∑τ=1Hβτ⋅(c22(1+cb+d)2+d1−α¯τ)≤2(c22(1+cb+d)2+d)∑τ=1Hβτ1−α¯τ≤(ii)3(c′+c′)(c22(1+cb+d)2+d)logH. split \|l^θ( a) \|≤&2 _τ=1^H _τ·E_ _τ (0,I) [ \|s_τ^θ ( α_τ a+ 1- α_τ _τ ) \|^2+ \| _τ \|^21- α_τ ]\\ (i)≤&2 _τ=1^H _τ· (c_2^2 (1+c_b+ d )^2+ d1- α_τ )\\ ≤&2 (c_2^2 (1+c_b+ d )^2+d ) _τ=1^H _τ1- α_τ\\ (i)≤&3(c +c ) (c_2^2 (1+c_b+ d )^2+d ) H. split (43) Here (i)(i) is by Equation˜41 and (ii)(i) is by Lemma˜B.4. ∎ Lemma B.4. Suppose H≥3c′logH≥ 3c H. Then the schedule specified in LABEL:eq:scheduler satisfies ∑τ=1Hβτατ(1−α¯τ)≤32(c′+c′)logH split _τ=1^H _τ _τ(1- α_τ)≤ 32(c +c ) H split (44) Proof. We first notice that βτατ(1−α¯τ)=1−ατατ(1−α¯τ)=−1−α¯τ−11−α¯τ+1ατ=1−1−α¯τ−11−α¯τ+1−ατατ≤log(1−α¯τ1−α¯τ−1)+βτατ,τ≥2. split _τ _τ(1- α_τ)=& 1- _τ _τ(1- α_τ)=- 1- α_τ-11- α_τ+ 1 _τ=1- 1- α_τ-11- α_τ+ 1- _τ _τ\\ ≤& ( 1- α_τ1- α_τ-1 )+ _τ _τ, τ≥ 2. split (45) For τ=1τ=1, we have βτατ(1−α¯τ)=1−1−α¯τ−11−α¯τ+1−ατατ=1+βτατ _τ _τ(1- α_τ)=1- 1- α_τ-11- α_τ+ 1- _τ _τ=1+ _τ _τ. Summing up τ∈[H]τ∈[H] gives ∑τ=1Hβτατ(1−α¯τ)≤1+∑τ=2Hlog(1−α¯τ1−α¯τ−1)+∑τ=1Hβτατ≤1+log(1−α¯H1−α¯1)+∑τ=1Hβτminτατ≤(i)1+c′logH+c′logH1−c′logH/H≤(ii)32(c′+c′)logH. split _τ=1^H _τ _τ(1- α_τ)≤&1+ _τ=2^H ( 1- α_τ1- α_τ-1 )+ _τ=1^H _τ _τ≤ 1+ ( 1- α_H1- α_1 )+ _τ=1^H _τ _τ _τ\\ (i)≤&1+c H+ c H1-c H/H (i)≤ 32(c +c ) H. split (46) Here in (i)(i) we have used (1−α¯H)/(1−α¯1)≤1/β1=Hc′(1- α_H)/(1- α_1)≤ 1/ _1=H^c , ∑τ=1Hβτ≤∑τ=1Hc′logH/H≤c′logH _τ=1^H _τ≤ _τ=1^Hc H/H≤ c H, and ατ≥1−c′logH/H _τ≥ 1-c H/H, and (ii)(i) follows by H≥3c′logH≥ 3c H. ∎ B.2.2 One-step Denoising Bound Lemma B.5 (One-step Denoising KL/TV Bound). Consider random variable X0∈ℝdX_0 ^d with score function s0⋆(x0)=∇x0logℙX0(x0)s_0 (x_0)= _x_0 _X_0(x_0) and one-step forward process X1=αX0+βϵx, splitX_1= αX_0+ β _x, split (47) where ϵx∼(0,I) _x (0,I), α,βα,β are positive scalars satisfying α+β=1α+β=1. Let s1⋆(x1)=∇x1logℙX1(x1)s_1 (x_1)= _x_1 _X_1(x_1) be the score function of X1X_1. Consider random variable Y1∈ℝdY_1 ^d and one-step backward process Y0=(Y1+βs1(Y1))/α+σϵy, splitY_0= (Y_1+β s_1(Y_1) )/ α+σ _y, split (48) where ϵy∼(0,I) _y (0,I), σ>0σ>0, and s1(y1)s_1(y_1) is the approximate score function of X1X_1. Suppose s0⋆(⋅)s_0 (·) is c2c_2-Lipschitz and suppose α/β>c2α/β>c_2. We have: DTV2(ℙX0∥ℙY0)≤A+2DTV2(ℙX1∥ℙY1),DKL(ℙX0∥ℙY0)≤A+DKL(ℙX1∥ℙY1), splitD_TV^2(P_X_0\|P_Y_0)≤&A+2D_TV^2(P_X_1\|P_Y_1), D_KL(P_X_0\|P_Y_0)≤ A+D_KL(P_X_1\|P_Y_1), split (49) where A=β22σ2αx1∼X1‖s1⋆(x1)−s1(x1)‖2+σ2dα/β−c2(c22+(1σ2−αβ)2)A= β^22σ^2αE_x_1 X_1 \|s_1 (x_1)-s_1(x_1) \|^2+ σ^2dα/β-c_2(c_2^2+( 1σ^2- αβ)^2). Proof. For notational simplicity, we will use ℙU|v(u)P_U|v(u) to abbreviate conditional probability ℙU|V(u|v)P_U|V(u|v) in this proof. We first prove the TV bound. Step 1. Connecting DTV(ℙX∥ℙY) D_TV(P_X_0\|P_Y_0) to the score difference. Define auxiliary random variable X~0=(X1+βs1(X1))/α+σϵ X_0= (X_1+β s_1(X_1) )/ α+σε for ϵ∼(0,I)ε (0,I). By the triangle inequality, we have that DTV2(ℙX0∥ℙY0)≤2DTV2(ℙX0∥ℙX~0)+2DTV2(ℙX~0∥ℙY0). splitD_TV^2(P_X_0\|P_Y_0)≤&2D_TV^2(P_X_0\|P_ X_0)+2D_TV^2(P_ X_0\|P_Y_0). split (50) Since X~0 X_0 and Y0Y_0 are generated by the same backward process (Equation˜48) from distributions X1X_1 and Y1Y_1, respectively, by the data processing inequality, we know that DTV(ℙX~0∥ℙY0)≤DTV(ℙX1∥ℙY1)D_TV(P_ X_0\|P_Y_0)≤ D_TV(P_X_1\|P_Y_1). Substituting back into Equation˜50 gives DTV2(ℙX0∥ℙY0)≤2DTV2(ℙX0∥ℙX~0)+2DTV2(ℙX1∥ℙY1). splitD_TV^2(P_X_0\|P_Y_0)≤&2D_TV^2(P_X_0\|P_ X_0)+2D_TV^2(P_X_1\|P_Y_1). split (51) For the first term, we have that DTV(ℙX0∥ℙX~0)=12∫|∫(ℙX0|x1(x0)−ℙX~0|x1(x0))ℙX1(x1)dx1|dx0≤12∫|ℙX0|x1(x0)−ℙX~0|x1(x0)|dx0ℙX1(x1)dx1=X1[DTV(ℙX0|x1∥ℙX~0|x1)]. splitD_TV(P_X_0\|P_ X_0)=& 12 | (P_X_0|x_1(x_0)-P_ X_0|x_1(x_0) )P_X_1(x_1)dx_1 |dx_0\\ ≤& 12 |P_X_0|x_1(x_0)-P_ X_0|x_1(x_0) |dx_0P_X_1(x_1)dx_1\\ =&E_X_1 [D_TV (P_X_0|x_1\|P_ X_0|x_1 ) ]. split (52) Substituting back gives DTV2(ℙX0∥ℙY0)≤2X12[DTV(ℙX0|x1∥ℙX~0|x1)]+2DTV2(ℙX1∥ℙY1)≤2X1[DTV2(ℙX0|x1∥ℙX~0|x1)]+2DTV2(ℙX1∥ℙY1)≤X1[DKL(ℙX0|x1∥ℙX~0|x1)]+2DTV2(ℙX1∥ℙY1). splitD_TV^2(P_X_0\|P_Y_0)≤&2E_X_1^2 [D_TV (P_X_0|x_1\|P_ X_0|x_1 ) ]+2D_TV^2(P_X_1\|P_Y_1)\\ ≤&2E_X_1 [D_TV^2 (P_X_0|x_1\|P_ X_0|x_1 ) ]+2D_TV^2(P_X_1\|P_Y_1)\\ ≤&E_X_1 [D_KL (P_X_0|x_1\|P_ X_0|x_1 ) ]+2D_TV^2(P_X_1\|P_Y_1). split (53) Here the last line is by Pinsker’s inequality. Since X~0|x1∼((x1+βs1(x1))/α,σ2I) X_0|x_1 ((x_1+β s_1(x_1))/ α,σ^2I), its negative log density −logℙX~0|x1- _ X_0|x_1 is 1/σ21/σ^2-strongly convex. Therefore, by [Vempala and Wibisono, 2019, Equation (4)], DKL(ℙX0|x1∥ℙX~0|x1)≤σ22x0∼ℙX0|x1‖∇x0logℙX0|x1(x0)−∇x0logℙX~0|x1(x0)‖2 splitD_KL(P_X_0|x_1\|P_ X_0|x_1)≤& σ^22E_x_0 _X_0|x_1 \| _x_0 _X_0|x_1(x_0)- _x_0 _ X_0|x_1(x_0) \|^2 split (54) By definition of X0|x1X_0|x_1 and X~0|x1 X_0|x_1, we can write the log probabilities explicitly: ∇x0logℙX0|x1(x0)=∇x0log(ℙX1|x0(x1)ℙX0(x0)ℙX1(x1))=∇x0(logℙX1|x0(x1)+logℙX0(x0))=−αβ(x0−x1α)+s0⋆(x0),∇x0logℙX~0|x1(x0)=−1σ2(x0−x1+βs1(x1)α). split _x_0 _X_0|x_1(x_0)=& _x_0 ( P_X_1|x_0(x_1)P_X_0(x_0)P_X_1(x_1) )= _x_0 ( _X_1|x_0(x_1)+ _X_0(x_0) )\\ =&- αβ (x_0- x_1 α )+s_0 (x_0),\\ _x_0 _ X_0|x_1(x_0)=&- 1σ^2 (x_0- x_1+β s_1(x_1) α ). split (55) Substituting back gives DKL(ℙX0|x1∥ℙX~0|x1)≤σ22x0∼ℙX0|x1‖s0⋆(x0)−βσ2αs1(x1)+(1σ2−αβ)(x0−x1α)⏟ℐ‖2. split&D_KL(P_X_0|x_1\|P_ X_0|x_1)≤ σ^22E_x_0 _X_0|x_1 \| s_0 (x_0)- βσ^2 αs_1(x_1)+ ( 1σ^2- αβ ) (x_0- x_1 α )_I \|^2. split (56) Finally, combining the above inequality with Equation˜53 gives DTV2(ℙX0∥ℙY0)≤X1[DKL(ℙX0|x1∥ℙX~0|x1)]+2DTV2(ℙX1∥ℙY1)≤σ22x1∼ℙX1x0∼ℙX0|x1‖ℐ‖2+2DTV2(ℙX1∥ℙY1). splitD_TV^2(P_X_0\|P_Y_0)≤&E_X_1 [D_KL (P_X_0|x_1\|P_ X_0|x_1 ) ]+2D_TV^2(P_X_1\|P_Y_1)\\ ≤& σ^22E_x_1 _X_1E_x_0 _X_0|x_1 \|I \|^2+2D_TV^2(P_X_1\|P_Y_1). split (57) Step 2. Upper bounding score difference term ℐ I. Before further derivation, we show ℙX0|x1(x0)P_X_0|x_1(x_0) is strongly log-concave. Note that −∇x02logℙX0|x1(x0)=−∇x0(s0⋆(x0)−αβ(x0−x1α))=−∇x0s0⋆(x0)+αβI⪰(αβ−c2)I. split- _x_0^2 _X_0|x_1(x_0)=&- _x_0 (s_0 (x_0)- αβ (x_0- x_1 α ) )=- _x_0s_0 (x_0)+ αβI\\ & ( αβ-c_2 )I. split (58) Here the last line is because s0⋆(⋅)s_0 (·) is c2c_2-Lipschitz. Since α/β>c2α/β>c_2, ℙX0|x1P_X_0|x_1 is (α/β−c2)(α/β-c_2)-strongly-log-concave. Utilizing this fact, we will now upper bound x0∼ℙX0|x1‖ℐ‖2≤‖x0∼ℙX0|x1[ℐ]‖2+tr(Varx0∼ℙX0|x1[ℐ]). splitE_x_0 _X_0|x_1 \|I \|^2≤\|E_x_0 _X_0|x_1[I]\|^2+tr (Var_x_0 _X_0|x_1 [I ] ). split (59) Firstly, we derive x0∼ℙX0|x1[ℐ]E_x_0 _X_0|x_1[I]. Since ℙX0|x1P_X_0|x_1 is strongly log-concave, it decays exponentially to 0 as ‖x0‖→∞ \|x_0 \|→∞. Therefore, by the divergence theorem, for all i∈[d]i∈[d] ∫div(ℙX0|x1(x0)ei)dx0=∮∞ℙX0|x1(x0)ei⋅n^dS=0. split (P_X_0|x_1(x_0)e_i )dx_0= _∞P_X_0|x_1(x_0)e_i· ndS=0. split (60) Therefore, 0=∫div(ℙX0|x1(x0)ei)dx0=∫∇x0ℙX0|x1(x0)⋅eidx0+∫ℙX0|x1(x0)div(ei)dx0=∫∇ℙX0|x1(x0)⋅eidx0. split0=& (P_X_0|x_1(x_0)e_i )dx_0= _x_0P_X_0|x_1(x_0)· e_idx_0+ _X_0|x_1(x_0)div(e_i)dx_0\\ =& _X_0|x_1(x_0)· e_idx_0. split (61) Since the above holds for all i∈[d]i∈[d], we know that 0=∫∇x0ℙX0|x1(x0)dx0=X0|x1[∇x0logℙX0|x1(x0)]=X0|x1[−αβ(x0−x1α)+s0⋆(x0)]. split0=& _x_0P_X_0|x_1(x_0)dx_0=E_X_0|x_1 [ _x_0 _X_0|x_1(x_0) ]\\ =&E_X_0|x_1 [- αβ (x_0- x_1 α )+s_0 (x_0) ]. split (62) Here the last line is by Equation˜55. Therefore, combining with the definition of ℐI gives x0∼ℙX0|x1[ℐ]=x0∼ℙX0|x1[1σ2(x0−x1α)−βσ2αs1(x1)]=1σ2(x1+βs1⋆(x1)α−x1α)−βσ2αs1(x1)=βσ2α(s1⋆(x1)−s1(x1)). splitE_x_0 _X_0|x_1 [I ]=&E_x_0 _X_0|x_1 [ 1σ^2 (x_0- x_1 α )- βσ^2 αs_1(x_1) ]\\ =& 1σ^2 ( x_1+β s_1 (x_1) α- x_1 α )- βσ^2 αs_1(x_1)\\ =& βσ^2 α (s_1 (x_1)-s_1(x_1) ). split (63) Here the second line is by the Tweedie’s identity. Secondly, we bound tr(Varx0∼ℙX0|x1[ℐ])tr(Var_x_0 _X_0|x_1 [I ]). By the definition of ℐI, tr(Varx0∼ℙX0|x1[ℐ])=∑i=1dVarx0∼ℙX0|x1[s0⋆(x0)+(1σ2−αβ)x0]i. splittr(Var_x_0 _X_0|x_1 [I ])=& _i=1^dVar_x_0 _X_0|x_1 [s_0 (x_0)+ ( 1σ^2- αβ )x_0 ]_i. split (64) Since ℙX0|x1P_X_0|x_1 is (α/β−c2)(α/β-c_2)-strongly-log-concave, it satisfies the Poincare inequality [Vempala and Wibisono, 2019, Equation (24)] with constant 1α/β−c2 1α/β-c_2, leading to: tr(Varx0∼ℙX0|x1[ℐ])=∑i=1dVarx0∼ℙX0|x1[s0⋆(x0)+(1σ2−αβ)x0]i≤1α/β−c2∑i=1dx0∼ℙX0|x1[∥∇[s0⋆(x0)+(1σ2−αβ)x0]i∥2]≤2dα/β−c2(c22+(1σ2−αβ)2). splittr(Var_x_0 _X_0|x_1 [I ])=& _i=1^dVar_x_0 _X_0|x_1 [s_0 (x_0)+ ( 1σ^2- αβ )x_0 ]_i\\ ≤& 1α/β-c_2 _i=1^dE_x_0 _X_0|x_1 [ \|∇ [s_0 (x_0)+ ( 1σ^2- αβ )x_0 ]_i \|^2 ]\\ ≤& 2dα/β-c_2 (c_2^2+ ( 1σ^2- αβ )^2 ). split (65) Combining Equations˜59, 63 and 64 gives x0∼ℙX0|x1‖ℐ‖2≤β2σ4α‖s1⋆(x1)−s1(x1)‖2+2βdα−c2β(c22+(1σ2−αβ)2). splitE_x_0 _X_0|x_1 \|I \|^2≤ β^2σ^4α \|s_1 (x_1)-s_1(x_1) \|^2+ 2β dα-c_2β (c_2^2+ ( 1σ^2- αβ )^2 ). split (66) Step 3. Conclusion. Substituting the above inequality back into Equation˜57 gives DTV2(ℙX0∥ℙY0)≤β22σ2αx1∼X1‖s1⋆(x1)−s1(x1)‖2+σ2βdα−c2β(c22+(1σ2−αβ)2)+2DTV2(ℙX1∥ℙY1). splitD_TV^2(P_X_0\|P_Y_0)≤& β^22σ^2αE_x_1 X_1 \|s_1 (x_1)-s_1(x_1) \|^2\\ &+ σ^2β dα-c_2β (c_2^2+ ( 1σ^2- αβ )^2 )+2D_TV^2(P_X_1\|P_Y_1). split (67) Similarly, one can derive a bound on KL divergence: DKL(ℙX0∥ℙY0)≤X1[DKL(ℙX0|x1∥ℙX~0|x1)]+DKL(ℙX1∥ℙY1)≤σ22x1∼ℙX1x0∼ℙX0|x1‖ℐ‖2+DKL(ℙX1∥ℙY1)≤β22σ2αx1∼X1‖s1⋆(x1)−s1(x1)‖2+σ2βdα−c2β(c22+(1σ2−αβ)2)+DKL(ℙX1∥ℙY1). split&D_KL(P_X_0\|P_Y_0) _X_1 [D_KL (P_X_0|x_1\|P_ X_0|x_1 ) ]+D_KL(P_X_1\|P_Y_1)\\ ≤& σ^22E_x_1 _X_1E_x_0 _X_0|x_1 \|I \|^2+D_KL(P_X_1\|P_Y_1)\\ ≤& β^22σ^2αE_x_1 X_1 \|s_1 (x_1)-s_1(x_1) \|^2+ σ^2β dα-c_2β (c_2^2+ ( 1σ^2- αβ )^2 )+D_KL(P_X_1\|P_Y_1). split (68) ∎ Appendix C Proof of Theorem˜C.2 Theorem C.1 (Theorem 5 in Zhang et al. [2022]). Suppose state space S and action spaces ii∈[n]\A_i\_i∈[n] are finite. Suppose all stationary points of the return function J(π)J(π) are isolated. Suppose the estimated averaged Q functions Q^iπk Q_i^π^k are accurate, i.e., Q^iπk=Q¯iπk Q_i^π^k= Q_i^π^k. Suppose the updated policies πik=πi⋆,k _i^k= _i ,k. Then for learning rate η=(1−γ)2/(2n(rmax−rmin))η=(1-γ)^2/(2n(r_ -r_ )), ∑k=0K−1NE-gap(πk)2K≤6M(rmax−rmin)2c(1−γ)3nK. split _k=0^K-1NE-gap (π^k )^2K≤ 6M(r_ -r_ )^2c(1-γ)^3 nK. split Here related constants are defined as follows: rmax≔max,r(,),rmin=min,r(,),M≔supπmax1dπ(),c=infkcπk,dπ()≔(1−γ)0∼ρ0[∑t=0∞γtℙ(t=|0)],cπ≔min∑i⋆∈argmaxiQ¯iπ(,i)πi(i⋆|). splitr_ & _ s, ar( s, a), r_ = _ s, ar( s, a), M _π _ s 1d^π( s), c= _kc_π^k,\\ d^π( s) &(1-γ)E_ s_0 ρ^0 [ _t=0^∞γ^tP ( s_t= s| s_0 ) ], c_π _ s _ a_i ∈ _ a_i Q_i^π ( s, a_i ) _i ( a_i | s ). split Intuitively, the quadratic mean of the NE-gaps converges to 0 at a rate of (n/K)O( n/K), where n is the number of agents and K is the number of epochs. The output policy is therefore an approximate NE. Theorem C.2. Consider the setting of Theorem˜C.1 and Algorithm˜1. Suppose in every epoch k, the learned policy πikπ^k_i of all agents i∈[n]i∈[n] satisfies DTV(πik(⋅|),πi⋆,k(⋅|))≤ϵisD_TV( _i^k(·| s), _i ,k(·| s))≤ _is for all i∈[n],k∈[K],∈i∈[n],k∈[K], s , then ∑k=0K−1NE-gap(πk)2K≤6M(rmax−rmin)2c(1−γ)3nK+12M(rmax−rmin)2c(1−γ)4n2ϵis, split _k=0^K-1NE-gap (π^k )^2K≤ 6M(r_ -r_ )^2c(1-γ)^3 nK+ 12M(r_ -r_ )^2c(1-γ)^4n^2 _is, split where M and c are defined in Theorem˜C.1. ∎ Proof of Theorem˜C.2. We let πik(⋅|)π^k_i(·| s) denote the approximate policy learned by importance sampling score matching with N samples from πik−1(⋅|)π^k-1_i(·| s), and let π~ik(⋅|) π^k_i(·| s) denote the exact softmax policy update from πk−1π^k-1, i.e., π~ik(⋅|)∝πik−1(⋅|)exp(η1−γQ¯iπk−1(,⋅)) π_i^k(·| s) _i^k-1(·| s) ( η1-γ Q_i^π^k-1( s,·) ). By assumption, DTV(π~ik(⋅|),πik(⋅|))≤ϵisD_TV ( π_i^k(·| s), _i^k(·| s) )≤ _is for all i∈[n],k∈[K],∈i∈[n],k∈[K], s . Step 1: Perturbation bound on J. We first bound ‖J(π)−J(π~)‖ \|J(π)-J( π) \| for any two joint policies π=(π1,⋯,πn)π=( _1,·s, _n) and π~=(π~1,⋯,π~n) π=( π_1,·s, π_n) satisfying DTV(πi(⋅|),π~i(⋅|))≤ϵD_TV( _i(·| s), π_i(·| s))≤ε for all i,i, s. To do this, we introduce auxiliary policies νi≔(π1,…,πi,π~i+1,…,π~n),∀i∈[n]ν^i ( _1,…, _i,\, π_i+1,…, π_n ),∀ i∈[n] and ν0≔π~ν^0 π. By telescoping, we know that J(π)−J(π~)=∑i=1n(J(νi)−J(νi−1)).J(π)-J( π)= _i=1^n (J(ν^i)-J(ν^i-1) ). (69) Since intermediate policies νiν^i and νi−1ν^i-1 only differ in the i-th component, we can apply the performance difference lemma (Lemma 8 of Zhang et al. [2022]) and get: J(νi)−J(νi−1)=11−γ∑,idνi()πi(i|)A¯iνi−1(,i).J(ν^i)-J(ν^i-1)= 11-γ _ s, a_id^ν^i( s) _i( a_i| s)\, A^ν^i-1_i( s, a_i). (70) Here A¯iπ(,i):=∑−iπ−i(−i|)Aiπ(,i,−i) A^π_i( s, a_i):= _ a_-i _-i( a_-i| s)\,A^π_i( s, a_i, a_-i) and Aiπ(,i,−i)=Qiπ(,i,−i)−Viπ()A^π_i( s, a_i, a_-i)=Q^π_i( s, a_i, a_-i)-V^π_i( s). By the definition of A¯iπ(,i) A^π_i( s, a_i), we know that ∑iπi(i|)A¯iπ(,i)=0 _ a_i _i( a_i| s)\, A^π_i( s, a_i)=0. Since νi−1ν^i-1 has π~i π_i as the i-th component, ∑iπ~i(i|)A¯iνi−1(,i)=0 _ a_i π_i( a_i| s)\, A^ν^i-1_i( s, a_i)=0. Substituting this equation into Equation˜70 leads to J(νi)−J(νi−1)=11−γ∑,idνi()(πi(i|)−π~i(i|))A¯iνi−1(,i)J(ν^i)-J(ν^i-1)= 11-γ _ s, a_id^ν^i( s) ( _i( a_i| s)- π_i( a_i| s) ) A^ν^i-1_i( s, a_i) (71) Taking the absolute values gives: ‖J(νi)−J(νi−1)‖≤ \|J(ν^i)-J(ν^i-1) \|≤ 11−γ∑dνi()∥πi(⋅|)−π~i(⋅|)∥1⋅maxi|A¯iνi−1(,i)| 11-γ _ sd^ν^i( s)\, \| _i(·| s)- π_i(·| s) \|_1· _ a_i | A^ν^i-1_i( s, a_i) | (72) ≤ ≤ 11−γ⋅∑dνi()2DTV(πi(⋅|),π~i(⋅|))⋅rmax−rmin1−γ 11-γ· _ sd^ν^i( s)2D_TV ( _i(·| s), π_i(·| s) )· r_ -r_ 1-γ (73) ≤ ≤ 2(rmax−rmin)(1−γ)2ϵ. 2(r_ -r_ )(1-γ)^2ε. (74) Substituting into (69) gives ‖J(π)−J(π~)‖≤2n(rmax−rmin)(1−γ)2ϵ. \|J(π)-J( π) \|≤ 2n(r_ -r_ )(1-γ)^2ε. (75) Step 2: Reward improvement from importance sampling score matching. By Lemma 20 and Lemma 21 in Zhang et al. [2022], since π~k π^k is the exact softmax update from πk−1π^k-1 with accurate Q functions, we have J(π~k)−J(πk−1)≥Lemma 201η∑i=1n∑dπ~k()log(∑iπik−1(i|)exp(η1−γA¯iπk−1(,i)))≥Lemma 21cη3MNE-gap(πk−1)2, splitJ( π^k)-J(π^k-1) Lemma 20≥& 1η _i=1^n _ sd π^k( s) ( _ a_i _i^k-1( a_i| s) ( η1-γ A^π^k-1_i( s, a_i) ) )\\ Lemma 21≥& cη3MNE-gap (π^k-1 )^2, split (76) where constants c and M are defined in Theorem˜C.1. Combining with Step 1 gives: J(πk)−J(πk−1)=J(πk)−J(π~k)+J(π~k)−J(πk−1)≥cη3MNE-gap(πk−1)2−2n(rmax−rmin)(1−γ)2ϵis. splitJ(π^k)-J(π^k-1)=&J(π^k)-J( π^k)+J( π^k)-J(π^k-1)\\ ≥& cη3MNE-gap (π^k-1 )^2- 2n(r_ -r_ )(1-γ)^2 _is. split (77) Step 3: Connecting reward improvement to Nash gap. Summing from k=1k=1 to K gives: J(πK)−J(π0)≥cη3M∑k=0K−1NE-gap(πk)2−2nK(rmax−rmin)(1−γ)2ϵis. splitJ(π^K)-J(π^0)≥ cη3M _k=0^K-1NE-gap (π^k )^2- 2nK(r_ -r_ )(1-γ)^2 _is. split (78) Since J(πK)−J(π0)≤(rmax−rmin)/(1−γ)J(π^K)-J(π^0)≤(r_ -r_ )/(1-γ), we have: 1K∑k=0K−1NE-gap(πk)2≤3M(rmax−rmin)cη(1−γ)1K+6Mn(rmax−rmin)cη(1−γ)2ϵis. split 1K _k=0^K-1NE-gap (π^k )^2≤& 3M(r_ -r_ )cη(1-γ) 1K+ 6Mn(r_ -r_ )cη(1-γ)^2 _is. split (79) Substituting η=(1−γ)2/(2n(rmax−rmin))η=(1-γ)^2/ (2n(r_ -r_ ) ) gives: 1K∑k=0K−1NE-gap(πk)2≤6Mn(rmax−rmin)2c(1−γ)31K+12Mn2(rmax−rmin)2c(1−γ)4ϵis. split 1K _k=0^K-1NE-gap (π^k )^2≤& 6Mn(r_ -r_ )^2c(1-γ)^3 1K+ 12Mn^2(r_ -r_ )^2c(1-γ)^4 _is. split (80) This completes the proof. ∎ Appendix D Experiment Details Name Value Name Value Rollout batch size BrB_r 20 Diffusion noise schedule linear Update batch size B 256 Diffusion noise schedule start 0.001 Smoothing parameter ξ 0.005 Diffusion noise schedule end 0.999 Score/Q network (MLP) size 3x256 Diffusion steps 20 Score/Q network activation mish Adam learning rate linear(1e-4,5e-5) Diffusion policy update lr η 1 Table 3: Hyperparameters for Section˜5