Hierarchical Multilevel Monte Carlo for Order-Optimal Neural Actor-Critic in Average-Reward CMDPs Ankur Naskar1 , Vaneet Aggarwal2 , 1
Computer Science and Automation, Indian Institute of Science, India 2 School of Industrial Engineering, Purdue University, USA [email protected], [email protected]
arXiv:2607.28390v1 [cs.LG] 30 Jul 2026
Abstract Constrained Markov Decision Processes (CMDPs) provide a natural framework for reinforcement learning in safety-critical applications, where agents maximize long-term reward while satisfying long-term constraints. Although primal-dual actorcritic methods with linear critics are well understood, extending order-optimal convergence guarantees to neural critics in average-reward CMDPs has remained open. The main challenge is a fundamental bias-cost tradeoff in neural critic estimation: under Neural Tangent Kernel (NTK) analysis, reducing critic bias requires substantially increasing critic optimization cost, preventing order-optimal convergence in the primal-dual framework. We resolve this bottleneck by introducing a hierarchical Multilevel Monte Carlo (MLMC) neural critic that performs debiasing simultaneously across trajectory sampling and critic optimization. The resulting estimator attains the bias of a long critic optimization run with only logarithmic expected sample cost. Building on this estimator, we develop a primal-dual Natural Actor-Critic algorithm that achieves both an optimality gap and a constraint violation of order Õ(T −1/2 ). This establishes the first order-optimal convergence guarantees for infinite-horizon average-reward CMDPs with general policy parameterization and neural critics, while eliminating the need to know the underlying mixing time. Our results are novel even in the unconstrained setting.
1
Introduction
Reinforcement learning (RL) has emerged as a powerful framework for sequential decision-making, with recent successes largely driven by deep neural networks. Many applications require agents not only to maximize long-term performance but also to satisfy resource or operational constraints, including transportation (Al-Abbasi, Ghosh, and Aggarwal 2019; Haliem et al. 2021), communication networks (Panju et al. 2021), robotics (Chen, Lan, and Aggarwal 2023; Gonzalez et al. 2023), and healthcare (Tamboli et al. 2024). Such problems are naturally modeled as constrained Markov Decision Processes (CMDPs), where the objective is to maximize the long-run average reward while satisfying long-term safety or resource constraints. Compared with discounted formulations, the average-reward criterion directly captures steady-state performance and is therefore particularly suitable for continuing decision-making tasks. Primal-dual policy-gradient and actor-critic methods have become the dominant framework for solving CMDPs due
to their scalability and compatibility with function approximation. In these methods, the actor updates a parameterized policy using gradient information estimated by a critic, while a dual variable enforces long-term constraints. For linear critics, recent advances have established order-optimal convergence guarantees of Õ(T −1/2 ) under general policy parameterization (Xu et al. 2025; Mondal and Aggarwal 2024), matching the minimax sample complexity lower bound for average-reward reinforcement learning (Jin and Sidford 2021). Despite these advances, extending order-optimal guarantees to neural critics remains an important open problem. Neural networks provide substantially greater representational power than linear critics and constitute the standard function approximator in modern reinforcement learning. Recent analyses based on the Neural Tangent Kernel (NTK) have enabled finite-sample guarantees for neural actor-critic methods by showing that, within a neighborhood of the network initialization, neural critics behave approximately as linear models. These techniques have led to order-optimal convergence for unconstrained discounted MDPs (Ganesh et al. 2025). However, for average-reward CMDPs, the best existing result achieves only a Õ(T −1/4 ) convergence rate (Satheesh et al. 2026), leaving a substantial gap from the optimal Õ(T −1/2 ) rate attainable with linear critics (Xu et al. 2025). The difficulty is not merely the use of neural networks, but a fundamental bias-cost tradeoff introduced by neural critic estimation. Under the NTK framework, the projection required to maintain the linear approximation causes the squared critic bias to decay at the same rate as its meansquared error (MSE). Consequently, reducing the critic bias necessitates running the critic optimization for substantially more iterations. In contrast, the primal-dual actor-critic algorithm operates within a nested-loop architecture, where each outer primal-dual update can only afford a limited amount of critic computation while maintaining the overall sample budget. As a result, the critic MSE remains of order Θ̃(1) (this implies Õ(1) squared critic bias under NTK) whereas orderoptimal convergence requires the squared critic bias to decay at a rate Õ(T −1 ). This creates a fundamental bias-cost bottleneck that prevents existing neural actor-critic algorithms from attaining order-optimal convergence.
In this work, we resolve this bottleneck by introducing a hierarchical Multilevel Monte Carlo (MLMC) neural critic. Existing MLMC techniques randomize trajectory lengths to remove sampling bias under Markovian observations. Such estimators inherit the bias of a long trajectory while incurring only a logarithmic sampling cost. However, they do not address the optimization bias arising from finite neural critic training. Our key insight is to introduce a second MLMC layer that randomizes the optimization horizon itself, thereby debiasing critic optimization in addition to trajectory sampling. Our main contributions are summarized as follows. • Resolving the Neural Critic Bias-Cost Tradeoff: We identify the fundamental bottleneck preventing orderoptimal convergence of primal-dual actor-critic methods with neural critics. To overcome this challenge, we introduce a novel hierarchical MLMC neural critic that performs multilayer debiasing across both trajectory sampling and critic optimization. Our estimator inherits the bias of a vanilla critic run for T iterations while requiring only O(log T ) expected sample cost. • First Order-Optimal Guarantees with Neural Critics: Building upon the proposed hierarchical MLMC neural critic, we establish global convergence and constraintviolation rates of Õ(T −1/2 ). To the best of our knowledge, this is the first primal-dual Natural Actor-Critic algorithm with neural critics and general policy parameterization to achieve order-optimal convergence for infinite-horizon CMDPs. Our results are novel even for the unconstrained average-reward setting. • Mixing-Time-Free Algorithm: Unlike existing methods that achieve order-optimal convergence with neural critics, our algorithm requires no knowledge of the mixing time of the underlying CMDP, making it considerably more practical for large-scale reinforcement learning problems.
1.1
Related Works
Average-Reward Constrained MDPs: CMDPs, introduced by (Altman 2021), provide the standard framework for safe reinforcement learning. Recent work has focused on primal-dual policy gradient methods for high-dimensional CMDPs (Ding et al. 2020; Paternain et al. 2019). Initial finite-sample guarantees were established for tabular and linear-function-approximation settings (Xu, Liang, and Lan 2021; Chen, Jain, and Luo 2022; Bai, Bedi, and Aggarwal 2023), while (Agarwal, Bai, and Aggarwal 2022) and (Ghosh, Zhou, and Shroff 2023) achieved optimal convergence with zero constraint violation in model-based and linear modelfree settings, respectively. For general policy parameterizations, (Bai, Mondal, and Aggarwal 2024a) established a suboptimal rate, whereas (Xu et al. 2025) obtained the optimal Õ(T −1/2 ) convergence rate. Actor-Critic Methods with Neural Critics: Early finitesample analyses of actor-critic methods assumed linear critics (Khodadadian, Chen, and Maguluri 2021; Cayci, He, and Srikant 2024). More recent works leverage Neural Tangent
Kernel (NTK) theory (Jacot, Gabriel, and Hongler 2018) to analyze neural critics in discounted MDPs (Fu, Yang, and Wang 2021; Tian, Olshevsky, and Paschalidis 2023; Cayci, He, and Srikant 2024; Gaur et al. 2024; Ganesh et al. 2025). In particular, (Ganesh et al. 2025) established a Õ(T −1/2 ) convergence rate for unconstrained discounted MDPs. For average-reward constrained MDPs, (Satheesh et al. 2026) established a Õ(T −1/4 ) convergence rate, which remains the state-of-the-art. MLMC-Based Estimation: Existing analyses under Markovian sampling often rely on sample thinning (Gaur et al. 2024; Ganesh, Mondal, and Aggarwal 2025a), which requires knowledge of the mixing time. Recent works have instead leveraged Multilevel Monte Carlo (MLMC) methods (Blanchet and Glynn 2015; Beznosikov et al. 2023; Patel et al. 2024; Ganesh, Mondal, and Aggarwal 2025b) to construct bias-controlled estimators via randomized trajectory lengths, thereby eliminating the need for data dropping and mixing-time estimates. In the CMDP setting, (Xu et al. 2025) used MLMC to mitigate the bias-cost tradeoff in primal-dual natural actor-critic algorithms with linear critics.
2
Problem Setup
We consider an infinite-horizon average-reward constrained Markov Decision Process (CMDP) M := (S, A, P, r, c, ρ), where S and A are the state and action spaces, respectively, P : S × A → ∆(S) denotes the transition kernel,1 r : S × A → [0, 1] and c : S × A → [−1, 1] are reward and cost functions, respectively, and ρ ∈ ∆(S) denotes the initial state distribution. A stationary policy π : S → ∆(A) generates a trajectory (st , at )∞ t=0 , where s0 ∼ ρ, at ∼ π(·|st ), and st+1 ∼ P(·|st , at ). Throughout this work, we make the following assumption on the underlying CMDP. Assumption 1 (Ergodicity). For every policy π, the Markov chain (st ) induced by π is irreducible and aperiodic. Under Assumption 1, for every policy π, the Markov chain (st ) has a unique stationary distribution dπ ∈ (0, 1)S , independent of ρ, such P that (dπ )⊤ P π = (dπ )⊤ , where for states s, s′ , P π (s, s′ ) := a π(a|s) P(s′ |s, a). The following is an important quantity associated with the CMDP M. π Definition 1 (Mixing Time). The mixing time τmix of M under policy π is defined as n o π τmix := min t > 0 : max (P π )t (s, ·) − dπ TV ≤ 1/4 , s
where ∥ · ∥TV denotes the total variation norm. Further, the π uniform mixing time of M is defined as τmix := supπ τmix . The mixing time quantifies the convergence rate of (st ) generated by any policy π to its stationary distribution dπ . Assumption 1 implies τmix < ∞ (Satheesh et al. 2026). Given a policy π, we define its average reward Jrπ and average constraint cost Jcπ as "T −1 # X 1 π Ju := lim E u(st , at ) s0 ∼ ρ , (1) T →∞ T t=0 1
∆(U ) denotes the probability simplex over a finite set U.
Reference Ganesh, Mondal, and Aggarwal (2025b) Xu et al. (2025) Xu et al. (2025) Ganesh et al. (2025) Satheesh et al. (2026) Our Work
Neural Critic ✗ ✗ ✗ ✓ ✓ ✓
Convergence Rate e −1/2 ) O(T e −1/2+ε ) O(T e −1/2 ) O(T e −1/2 ) O(T e −1/4 ) O(T e −1/2 ) O(T
Constraint Violation N/A e O(T −1/2+ε ) e −1/2 ) O(T N/A e −1/4 ) O(T e −1/2 ) O(T
Unknown Mixing Time ✓ ✓ ✗ ✗ ✓ ✓
Table 1: Comparison with related works under general policy parameterizations. The table indicates the type of critic employed (✗ for linear vs ✓ for neural), the convergence and constraint-violation rates achieved (with “N/A” denoting unconstrained MDPs), and whether knowledge of the mixing time is required. for each u ∈ {r, c}, where the expectation is taken over trajectories (st , at )∞ t=0 generated by policy π. Under the ergodicity assumption, we have Juπ = E(s,a)∼ν π [u(s, a)], where ν π denotes the stationary occupancy measure defined as ν π (s, a) := dπ (s)π(a|s), for each state-action pair (s, a). A common practice, when dealing with a large or continuous state space, is to parameterize the policy as πθ , where θ ∈ Θ ⊂ Rd and d ≪ |S||A|. Hereafter, we adopt the shorthand Juθ = Juπθ and ν θ = ν πθ . Further, we let Eθ [·] denote the expectation w.r.t. (s, a, s′ , a′ ) satisfying (s, a) ∼ ν θ , s′ ∼ P(·|s, a), and a′ ∼ πθ (·|s′ ). Given these quantities, our objective is to solve the following constrained optimization: max Jrθ θ∈Θ
s.t. Jcθ ≥ 0.
(2)
To ensure that (2) admits a feasible solution, we make the following assumption (Wei, Liu, and Ying 2022; Bai, Mondal, and Aggarwal 2024a). Assumption 2 (Slater Condition). There exists δ ∈ (0, 1) and θ̄ ∈ Θ such that Jcθ̄ ≥ δ.
2.1
The primal-dual approach reformulates the constrained optimization (2) as the saddle-point problem θ∈Θ λ≥0
(3)
where L is called the Lagrangian, and λ is referred to as the dual variable. A common gradient-based scheme for (3) is Natural Policy Gradient (NPG), defined by the updates θk+1 = θk + αF (θk )† ∇θ L(θk , λk ), λk+1 = Π[0,2/δ] λk − βJcθk . Here, F (θ) is the Fisher information matrix defined as i h ⊤ F (θ) := Eθ ∇θ ln πθ (a|s) (∇θ ln πθ (a|s)) ,
(4)
(5)
† denotes the Moore-Penrose pseudoinverse, α and β are the primal and dual stepsizes, respectively, δ is the constant from Assumption 2, and ΠA denotes projection onto a set A ⊆ R. The NPG vector wk∗ := F (θk )† ∇L(θk , λk ) is obtained by minimizing the function f (θk , λk , w) :=
1 ⊤ w F (θk )w − w⊤ ∇θ L(θk , λk ). 2
∇w f (θk , λk , w) = F (θk )w − ∇θ Jrθk + λk ∇θ Jcθk , (7) which in turn requires computing ∇θ Juθ , u ∈ {r, c}. By the policy gradient theorem (Sutton et al. 1999), for each u ∈ {r, c}, we have ∇θ Juθ = Eθ Aθu (s, a) ∇θ ln πθ (a|s) , (8) where Aθu ∈ R|S||A| denotes the advantage function defined as Aθu (s, a) := Qθu (s, a) − Ea′ ∼πθ (·|s) Qθu (s, a′ ) , (9) and Qθu ∈ R|S||A| is a solution to the Bellman equation Q(s, a) = u(s, a) − Juθ + Eθ [Q(s′ , a′ )] . (10) Such a function is unique (up to additive constants) and is called the state-action value function for the policy πθ . Consequently, implementing (4) hinges on estimating (Juθk , Qθuk ) for u ∈ {r, c}. We now describe the critic estimator used for this purpose.
2.2
Primal-Dual Natural Policy Gradient
max min L(θ, λ) := Jrθ + λJcθ ,
This is achieved by applying gradient descent with
(6)
Neural Critic Estimation
Given a policy πθ , the average reward/cost Juθ can be estimated with relative ease. However, estimating the function Qθu is often computationally challenging. In this work, we employ a neural critic for this task. Let ϕ : S × A → Rd be a feature map such that ∥ϕ(s, a)∥ ≤ 1, for every state(1) action pair (s, a). Given matrices Wu ∈ Rm×d and (ℓ) m×m {Wu }L , we let the vector ℓ=2 ⊂ R h i⊤ (L) (1) ζu := Vec⊤ Wu , . . . , Vec⊤ Wu serve as our critic parameter.2 These matrices consist of trainable weights in an L-layer feedforward neural network defined as 1 x(ℓ) σ Wu(ℓ) x(ℓ−1) (s, a) , ℓ ∈ [L], (11) u (s, a) := √ u m (0)
where xu (s, a) := ϕ(s, a) and σ is an element-wise activation function. For each u ∈ {r, c}, the output of this neural network lets us approximate Qθu (s, a) as 1 Q(s, a, ζu ) := √ b⊤ x(L) (s, a), (12) m u u 2
Vec(·) denotes vectorization by stacking columns.
where bu ∈ Rm is a fixed vector. Throughout the paper, we make the following assumption. Assumption 3 (Smooth Activation). The activation function σ is L1 -Lipschitz and L2 -smooth. Remark 1. Assumption 3 is satisfied by commonly used smooth approximations of ReLU such as Sigmoid, ELU, and GeLU (Hendrycks and Gimpel 2016; Clevert, Unterthiner, and Hochreiter 2015). For subsequent analysis, we rely on the Neural Tangent Kernel (NTK) regime, which requires restricting the critic parameter ζu to a ball B of radius R around the initialization ζ0 . This lets us approximate the neural network output Q(·, ·, ζ) by the linear function class n o FB,m := Qlin (·, ·, ζ) : ζ ∈ B , (13) where Qlin (·, ·, ζ0 ), is defined as Qlin (·, ·, ζ) := Q(·, ·, ζ0 ) + (ζ − ζ0 )⊤ ∇ζ Q(·, ·, ζ0 ). (14)
3
Algorithm
We now present our hierarchical MLMC-based primaldual Natural Actor-Critic algorithm (HiMLMC-PD-NAC) for solving (2), with its pseudocode provided in Algorithm 1. As the transition kernel and occupancy measures are unknown, we rely on sample-based estimates. Our algorithm updates the primal-dual pair (θk , λk ) via the update rule θk+1 = θk + αwk λk+1 = Π[0,2/δ] (λk − βηc,k ) ,
⊤ ξu,k,MLMC = ηu,k,MLMC , (ζu,k,MLMC )⊤ i h k k k := ξu,1 + 1{2Hk ≤Hmax } 2Hk ξu,h − ξ u,hk /2 , k
(17)
k where for H ∈ {1, hk , hk /2}, ξu,H denotes the output from the vanilla-critic subroutine after H iterations.
Inner MLMC layer: The vanilla-critic subroutine seeks to solve the optimization h i 1 2 min Ru (θk , η) := Eθk (η − u(s, a)) η 2 h 2 i 1 min Eu (θk , ζ) := Eθk Q(s, a, ζ) − Qθuk (s, a) . ζ 2 We optimize these objectives together using Stochastic Gradient Descent (SGD). At each iteration 0 ≤ h < hk , we es- k k timate the joint gradient ∇η Ru (θk , ηu,h ), ∇ζ E(θk , ζu,h ) and update the vector k ⊤ k k ξu,h = ηu,h , (ζu,h )⊤ . To reduce gradient bias, we employ a second MLMC construction. Let Uk,h ∼ Geom(1/2) and define ℓk,h := 1 + 1{2Uk,h ≤Tmax } 2Uk,h − 1 . (18) ℓ
(15)
where α, β, and δ are as in (4), whereas wk and ηc,k are approximations to the NPG vector wk∗ and the averageconstraint-cost Jcθk associated with policy πθk , respectively. The update (15) is run for K > 0 iterations, yielding the output {(θk , λk )}k<K .
3.1
and the MLMC critic parameter
Critic Estimator
For each outer-loop iteration 0 ≤ k < K, we approximate the critic (Juθk , Qθuk ) by an estimator (ηu,k,MLMC , Qu,k,MLMC ). The critic bias is controlled using a hierarchical two-layer MLMC scheme with outer- and inner-layer truncation parameters Hmax and Tmax , respectively. The outer MLMC layer operates on several outputs from a vanilla-critic subroutine, with each vanilla output corresponding to a different number of iterations. Specifically, after hk iterations, k the vanilla subroutine produces the output ξu,h . These k vanilla outputs are then combined by the outer layer to form (ηu,k,MLMC , Qu,k,MLMC ). The inner MLMC layer is used for gradient estimation within the vanilla subroutine. Outer MLMC layer: Sample Hk ∼ Geom(1/2) and let hk := 1 + 1{2Hk ≤Hmax } (2Hk − 1). Then, we define the critic estimation function k Qu,k,MLMC (·, ·) := Q ·, ·, ζu,1 + 1{2Hk ≤Hmax } h i k k × 2Hk Q ·, ·, ζu,h − Q ·, ·, ζu,h (16) k k /2
k,h k,h Given a trajectory T k,h := (sk,h t , at )t=0 generated by the policy πθk , we define the naive gradient-estimator as PT −1 k,h k,h k,h k,h v̂u,T := T1 t=0 v̂uk,h (sk,h t , at , st+1 , at+1 ), with
ˆ η Rk,h (s, a) c ∇ u v̂uk,h (s, a, s′ , a′ ) := ˆ γ k,h , ∇ζ Eu (s, a, s′ , a′ )
(19)
where cγ > 0 is a scaling parameter, and k ˆ η Rk,h ∇ u (s, a) := ηu,h − u(s, a), ˆ ζ Euk,h (s, a, s′ , a′ ) := ∇ζ Q(s, a, ζ0 ) ∇
k k k × ηu,h − u(s, a) + Q(s, a, ζu,h ) − Q(s′ , a′ , ζu,h ) . Then, the corresponding MLMC-based gradient estimator is defined as k,h k,h + 1{2Uk,h ≤Tmax } v̂u,MLMC := v̂u,1 × 2Uk,h v̂ k,hUk,h − v̂ k,hUk,h −1 . u,2
u,2
(20)
Finally, we run the following update rule for H iterations: k,h k k ξu,h+1 = ΠB ξu,h − γξ,H v̂u,MLMC , (21) where γξ,H is a stepsize dependent on H, while ΠB projects only the ζ portion onto the NTK ball B. After completing H iterations, with each H ∈ {1, hk /2, hk }, we have the quantities required in (17).
3.2
NPG Estimator
For each outer iteration 0 ≤ k < K, this subroutine estimates the NPG vector wk∗ associated with the current policy πθk . As discussed in Section 2.1, the NPG vector is obtained by minimizing f (θk , λk , w) defined in (6). To employ an MLMC-based SGD scheme, we sample Uk,h ∼ Geom(1/2), k,h ℓk,h and generate a trajectory T k,h := (sk,h t , at )t=0 using policy πθk , where ℓk,h is as defined in (18). We define the naive T -sample estimators for F (θk ) and ∇θ Juθk as T −1
F̂Tk,h :=
1 X k,h k,h k,h F̂ (st , at ), T t=0 T −1
and
k,h ĝu,T :=
1 X k,h k,h k,h k,h k,h k,h ĝ (st , at , āt , st+1 , at+1 ), T t=0 u
respectively, where āk,h ∼ πθk (·|sk,h t t ) and ⊤
F̂ k,h (s, a) := ∇θ ln πθk (a|s) (∇θ ln πθk (a|s)) , ĝuk,h (s, a, ā, s′ , a′ ) := ∇θ ln πθk (a|s) u(s, a) − ηu,k (22) + Qu,k,MLMC (s′ , a′ ) − Qu,k,MLMC (s, ā) . Then, the corresponding MLMC estimators are defined as k,h X̂MLMC := X̂1k,h + 1{2Uk,h ≤Tmax } k,h × 2Uk,h X̂ k,h − X̂ Uk,h −1 , Uk,h 2
2
(23)
for each X̂ k,h ∈ {F̂ k,h , ĝuk,h }. Finally, with these estimates, we run the following SGD update k wh+1 = whk h i k,h k,h k,h − γw F̂MLMC whk − ĝr,MLMC + λk ĝc,MLMC .
(24)
k After Hw iterations, we set wk := wH as the NPG estimate w for policy πθk .
Remark 2 (Expected Sample Cost). Since Hk and Uk,h are independent, we have E[hk ℓk,h ] = E[hk ] · E[ℓk,h ] = O(⌊log2 Hmax ⌋ · ⌊log2 Tmax ⌋). Thus, for each outer iteration, our expected critic-estimation cost is O(ln Hmax · ln Tmax ). Similarly, our NPG-estimation cost is O(Hw · ln Tmax ).
4
Assumptions and Main Results
Before we present our finite-time convergence bounds in Theorem 1, we state the following assumptions. Assumption 4 (Score Function). There exist G1 , G2 > 0, such that for all θ, θ′ ∈ Θ, and state-action pair (s, a), ∥∇θ ln πθ (a|s)∥ ≤ G1
and
∥∇θ ln πθ (a|s) − ∇θ ln πθ′ (a|s)∥ ≤ G2 ∥θ − θ′ ∥. Assumption 5 (Fisher Non-Degeneracy). There exists µF > 0, such that for all θ ∈ Θ, F (θ) ⪰ µF I.
Assumption 6 (Policy Parametrization Error). There exists εbias > 0 such that for all θ ∈ Θ and λ ∈ [0, 2/δ], we have min Lν π∗ (θ, λ, w) ≤ εbias , ∗
w∈Rd π∗
where π and ν denote an optimal stationary policy3 and its occupancy measure, respectively, while the function Lν (θ, λ, w) is called the transferred compatible function and is defined as h 2 i Eν w⊤ ∇θ ln πθ (a|s) − Aθr (s, a) − λAθc (s, a) . Remark 3. Assumptions 4–6 are standard in the PG literature (Liu et al. 2020; Papini et al. 2018; Xu, Gao, and Gu 2020; Fatkhullin et al. 2023; Ding et al. 2025). Assumption 4 requires the score function to be bounded and Lipschitz, and is satisfied by Gaussian policies with bounded variance and softmax policies. Assumption 5 ensures that f (θ, λ, w) is strongly convex and admits a unique minimizer. Assumption 6 characterizes the richness of the policy class. Assumption 7 (Critic Approximation Error). There exists a finite εapp > 0, such that for every u ∈ {r, c}, and θ ∈ Θ, we have h 2 i Eθ Qθu (s, a) − ΠFB,m Qθu (s, a) ≤ εapp , where ΠFB,m denotes projection onto the linear function class defined in (13). Remark 4. Assumption 7 ensures that the linear function class FB,m approximates the state-action value function with bounded error. Under NTK, it is equivalent to (Gaur et al. 2024, Assumption 2.5) and (Cayci, He, and Srikant 2024, Assumption 2). We make one final assumption commonly used for average-reward MDPs (Ganesh, Mondal, and Aggarwal 2025b; Xu et al. 2025; Satheesh et al. 2026). Assumption 8 (Smooth Objectives). The average-reward Jrθ and the average-constraint-cost Jcθ are L-smooth in the policy parameter θ. Assumption 8 is satisfied by several commonly used policy parameterizations, including the softmax policy. Under Assumptions 1–8, we have the following result. Theorem 1 (Global Convergence and Constraint ViolaT tion). Choose the vanilla-critic stepsizes γξ,H = 2µln , with ϕH suitable µϕ > 0, and the NPG stepsize γw = µ2FlnHTw . Set the network width m = Hmax , and the NTK radius R = Θ(ln T ). Set the primal-dual stepsizes α = β = √1K , and horizon lengths Tmax = Hmax = K = Θ(T ), and Hw = Θ(ln T ). Then, we have K−1 √ 1 X θk 1 √ π∗ Jr − E Jr = Õ √ + εbias + εapp , K T k=0 K−1 √ 1 X θk 1 √ − E Jc = Õ √ + εbias + εapp , K T k=0 where π ∗ is as defined in Assumption 6. 3 ∗ π is a solution to the constrained optimization: maxπ Jrπ subject to Jcπ ≥ 0.
Remark 5 (Optimal Global Convergence). Theorem 1 shows that, under the prescribed parameter√ choices, the optimality gap decays to zero at a rate of Õ(1/ T ), modulo the policy parametrization error εbias and the neural function approximation error εapp . Since the total expected sample cost is K(Hw + ln Hmax ) ln Tmax = Θ̃(T ), our convergence rate is optimal. Remark 6 (Constraint Violation Rate). The expected constraint√ violation rate also decays to zero at the rate of Õ(1/ T ) (up to the parametrization and neural critic approximation errors), matching the existing state-of-theart (Xu et al. 2025). Remark 7 (Mixing-Time Knowledge). The parameter choices in Theorem 1 require no knowledge of the mixing time of the underlying CMDP.
5
Proof Sketch
We now outline the proof of Theorem 1. Our proof consists of the following key steps. 1. Use the dual update in (15) to express the optimality gap and constraint violation as a combination of the Lagrangian error and critic bounds. 2. Bound the Lagrangian error using the primal update in (15) in terms of the NPG bias and MSE. 3. Exploit the linear structure of the NPG update (24) together with MLMC bounds for the Fisher matrix and policy-gradient estimates to control the NPG bias and MSE in terms of the critic parameter. 4. Bound the critic bias and MSE: (a) Bound the vanilla critic: decompose the vanilla critic update (21) into a perturbed linear recursion, then control the linear components using the inner MLMC layer (20), and the nonlinear perturbation via the NTKregime. (b) Bound the MLMC critic: combine the vanilla-critic bounds with the outer MLMC layer (17). Before we elaborate on these steps, we introduce the following notation: let Ek [·] := E[·|Fk ], where Fk denotes the σ-algebra generated by all samples till the (k − 1)-th outer loop iteration. Step 1. Optimality Gap and Constraint Violation: We borrow from the proof strategy of (Ganesh, Mondal, and Aggarwal 2025b), later adapted for primal-dual methods with a linear critic by (Xu et al. 2025). Let L(π ∗ , λ) be defined as ∗ ∗ Jrπ + λJcπ , where π ∗ is an optimal stationary policy. Using the dual update in (15) and the definition of L, we can bound the optimality gap and constraint violation as given by the following two lemmas.
Lemma 1 (Optimality Gap). Under Assumption 2, we have K−1 1 X π∗ E Jr − Jrθk K k=0
K−1
=O
1 Xh 2 ∗ βEk ξc,k − ξc,k,MLMC K k=0 i ∗ + E Ek ξc,k,MLMC − ξc,k ! K−1 1 X ∗ +β+ E L(π , λk ) − L(θk , λk ) . K k=0
Lemma 2 (Constraint Violation Rate). Under Assumption 2, we have K−1 K−1 1 X θk 1 Xh 2 ∗ − βE ξc,k,MLMC − ξc,k E Jc = O K K k=0 k=0 i ∗ + E Ek ξc,k,MLMC − ξc,k K−1 1 1 X ∗ +β+ E L(π , λk ) − L(θk , λk ) . + βK K k=0
The detailed proofs of Lemmas 1 and 2 are deferred to Appendices A and B, respectively. Step 2. Lagrangian Error: In light of Lemmas 1 and 2, our goal for the remainder of this section is to bound the Lagrangian error. Under Assumption 8, we have the following result. Lemma 3. Let (θk , λk ) be updated by (15). Under Assumptions 4–6, the following holds for 0 ≤ k < K : K−1 √ 1 X 1 E L(π ∗ , λk ) − L(θk , λk ) = O εbias + α + K αK k=0 K−1 K−1 α X 1 X E∥Ek [wk ] − wk∗ ∥ + E∥wk − wk∗ ∥2 . + K K k=0
k=0
The proof of Lemma 3 is given in Appendix C. Step 3. NPG Estimator Analysis: The NPG update (24) is a stochastic linear recursion with the mean drift (∇θ Jrθk + λk ∇θ Jcθk ) − F (θk )w. Therefore, standard arguments on stochastic linear recursions allow us to bound the NPG bias and MSE in terms of the biases and MSEs of the Fisher and policy gradient estimates. For the Fisher matrix, the bias arises from replacing the occupancy measure ν θk with finite-length trajectories generated by πθk . Standard MLMC arguments yield a squared bias of order O(τmix /Tmax ) and MSE of order O(ln Tmax ). In contrast, the advantage function (9) contributes additional bias to the PG k,h estimator ĝu,MLMC through the neural critic approximation. The NTK regime lets us replace the neural critic Q(·, ·, ζ) by its first-order Taylor approximation Qlin (·, ·, ζ) ∈ FB,m . Let h i⊤ ∗ ∗ ∗ ∗ ∗ := Juθk and ζu,k is a ξu,k := ηu,k , (ζu,k )⊤ , where ηu,k
minimum-norm solution to the following projected Bellman equation h Eθk ∇ζ Q(s, a, ζ0 ) Juθk − u(s, a) i + Qlin (s, a, ζ) − Qlin (s′ , a′ , ζ) = 0. (25) Then, we have the following result. Theorem 2 (NPG bounds). Let Assumptions 1– 7 hold and choose γw = µ2FlnHTw . Then, the following bounds hold: 1 τmix 1 2 Ek [wk ] − wk∗ = Õ + + + εapp 2 T Tmax m X 2 2 ∗ + δQ , + Ek [ξu,k,MLMC ] − ξu,k u∈{r,c}
1 1 + τmix + + εapp T2 m X 2 2 ∗ , + σQ + Ek ξu,k,MLMC − ξu,k
Ek wk − wk∗
2
= Õ
u∈{r,c}
where P 2 2 := u Ek [Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC )] , δQ P 2 2 := u Ek Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC ) . σQ Inside the NTK ball, the neural critic Q(·, ·, ζ) decomposes into the linear critic Qlin (·, ·, ζ) and a nonlinear remainder. This nonlinear remainder yields the Õ(1/m) term in Theorem 2. Whereas, εapp reflects the approximation capacity of the linear critic. Consequently, the analysis is reduced to controlling the critic bias and MSE. Deferring the proof of Theorem 1 to Appendix D, we now establish the critic bounds in the following step. Step 4. Critic Estimator Analysis: Recall the hierarchical MLMC scheme introduced in Section 3.1. Since the outer layer combines outputs of the vanilla-critic subroutine corresponding to different iteration budgets, as in (17), we begin k by analyzing the vanilla-critic ξu,H generated by (21) with an arbitrary iteration budget H. The NTK projection ΠB introduces an additional error in k . This projection error is controlled the bias analysis of ξu,H by the probability with which iterates escape the NTK ball, which, in turn, is bounded in terms of the vanilla-critic MSE (see (Ganesh et al. 2025, Appendix A.2)). To bound the vanilla-critic MSE, we decompose (21) as a perturbed linear recursion. Under the NTK regime, the linear critic Qlin (·, ·, ζ) approximates the neural critic up to an error Õ(m−1/2 ) (see Lemma 6). By leveraging this linear approximation, we can write (21) as: k,h k k k ξu,h+1 = ΠB ξu,h + γξ,H B̂u,MLMC − Âk,h MLMC ξu,h k,h + γξ,H Ṽu,MLMC , k,h where the linear component (B̂u,MLMC − Âk,h MLMC ξ) is obtained by replacing Q(·, ·, ζ) with Qlin (·, ·, ζ) in (19) and
using the inner-layer MLMC construction defined in (20). Collecting the remainder yields the nonlinear component k,h k,h Ṽu,MLMC . Under the NTK regime, we have ∥Ṽu,MLMC ∥ = −1/2 Õ(m ), while standard MLMC bounds apply on the linear component. This results in the following bound for the vanilla critic (detailed proof in Appendix E). k Lemma 4 (Vanilla-Critic Bounds). Let ξu,H be generated by (21) after H iterations and let Assumptions 1–7 hold. T Choose γξ,H = 2µln , for sufficiently small µϕ , and let cγ ϕH be sufficiently large. Then, for all u ∈ {r, c}, we have τmix 1 τmix τmix 2 k ∗ Ek ξu,H = Õ + − ξu,k + + . T2 Tmax H m Remark 8 (Bias-Cost Tradeoff). Owing to the NTK projection ΠB , our squared critic bias has the same Õ(H −1 ) rate as the MSE. In contrast, a linear critic yields squared bias −1 of order Õ(Tmax ) (Xu et al. 2025, Lemma 4.8). However, our primal-dual scheme requires a large outer-loop size K to control constraint violations, making a large H prohibitively k for gradient estimation expensive. Hence, directly using ξu,H introduces a fundamental bias-cost tradeoff. We overcome the above bias-cost tradeoff by exploiting the outer MLMC layer defined by (16) and (17). Recall that the iteration budget for the vanilla subroutine is set to hk := 1 + 1{2Hk ≤Hmax } (2Hk − 1) with Hk ∼ Geom(1/2). A telescoping-sum argument then shows that k E ξu,MLMC = E ξu,H and E hk = ⌊log2 Hmax ⌋ . max Therefore, while the MLMC estimator inherits the bias of a vanilla critic updated for Hmax iterations, the expected number of vanilla iterations is only O(ln Hmax ). Using a similar argument for the MSE, we prove the following bounds in Appendix E. Theorem 3 (MLMC-Critic Bounds). Consider the setting of Lemma 4. Then, the critic function Qu,k,MLMC satisfies 2 2 = Õ (Hmax /m) . Moreover, the = Õ (1/m) and σQ δQ critic parameter ξu,k,MLMC satisfies: ∗ Ek ξu,k,MLMC − ξu,k
2
= Õ(Hmax /m) , τmix τmix τmix 2 ∗ Ek [ξu,k,MLMC ] − ξu,k = Õ . + + Tmax Hmax m Setting m = Hmax results in a critic MSE of order Õ(1) −1/2 −1/2 and a critic bias of order Õ Tmax + Hmax . Finally, in Appendix F, we combine Lemmas 1–3 with Theorems 2 and 3, and apply the prescribed parameter choices to complete the proof of Theorem 1.
6
Conclusion
In this work, we study infinite-horizon average-reward CMDPs under general policy parameterization and neural critic estimation. We propose a primal-dual natural actorcritic algorithm equipped with a novel hierarchical MLMCbased critic and establish optimal convergence guarantees for neural critics in the CMDP literature.
References Agarwal, M.; Bai, Q.; and Aggarwal, V. 2022. Concave Utility Reinforcement Learning with Zero-Constraint Violations. Transactions on Machine Learning Research. Al-Abbasi, A. O.; Ghosh, A.; and Aggarwal, V. 2019. Deeppool: Distributed model-free algorithm for ride-sharing using deep reinforcement learning. IEEE Transactions on Intelligent Transportation Systems, 20(12): 4714–4727. Altman, E. 2021. Constrained Markov decision processes. Routledge. Bai, Q.; Bedi, A. S.; and Aggarwal, V. 2023. Achieving zero constraint violation for constrained reinforcement learning via conservative natural policy gradient primal-dual algorithm. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, 6737–6744. Bai, Q.; Mondal, W. U.; and Aggarwal, V. 2024a. Learning general parameterized policies for infinite horizon average reward constrained MDPs via primal-dual policy gradient algorithm. Advances in Neural Information Processing Systems, 37: 108566–108599. Bai, Q.; Mondal, W. U.; and Aggarwal, V. 2024b. Regret analysis of policy gradient algorithm for infinite horizon average reward markov decision processes. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 38, 10980–10988. Beznosikov, A.; Samsonov, S.; Sheshukova, M.; Gasnikov, A.; Naumov, A.; and Moulines, E. 2023. First order methods with markovian noise: from acceleration to variational inequalities. Advances in Neural Information Processing Systems, 36: 44820–44835. Blanchet, J. H.; and Glynn, P. W. 2015. Unbiased Monte Carlo for optimization and functions of expectations via multi-level randomization. In 2015 Winter Simulation Conference (WSC), 3656–3667. IEEE. Cayci, S.; He, N.; and Srikant, R. 2024. Finite-Time Analysis of Entropy-Regularized Neural Natural Actor-Critic Algorithm. Transactions on Machine Learning Research, 2024. Chen, J.; Lan, T.; and Aggarwal, V. 2023. Option-Aware Adversarial Inverse Reinforcement Learning for Robotic Control. In Proceedings IEEE International Conference on Robotics and Automation. Chen, L.; Jain, R.; and Luo, H. 2022. Learning infinitehorizon average-reward Markov decision process with constraints. In International Conference on Machine Learning, 3246–3270. Clevert, D.-A.; Unterthiner, T.; and Hochreiter, S. 2015. Fast and accurate deep network learning by exponential linear units (elus). arXiv preprint arXiv:1511.07289, 4(5): 11. Ding, D.; Zhang, K.; Basar, T.; and Jovanovic, M. 2020. Natural policy gradient primal-dual method for constrained markov decision processes. Advances in Neural Information Processing Systems, 33: 8378–8390. Ding, D.; Zhang, K.; Duan, J.; Basar, T.; and Jovanovic, M. R. 2025. Convergence and sample complexity of natural policy gradient primal-dual methods for constrained mdps. Journal of Machine Learning Research, 26(256): 1–76.
Dorfman, R.; and Levy, K. Y. 2022. Adapting to mixing time in stochastic optimization with markovian data. In International Conference on Machine Learning, 5429–5446. PMLR. Fatkhullin, I.; Barakat, A.; Kireeva, A.; and He, N. 2023. Stochastic policy gradient methods: Improved sample complexity for fisher-non-degenerate policies. In International Conference on Machine Learning, 9827–9869. Fu, Z.; Yang, Z.; and Wang, Z. 2021. Single-Timescale Actor-Critic Provably Finds Globally Optimal Policy. In International Conference on Learning Representations. Ganesh, S.; Chen, J.; Mondal, W. U.; and Aggarwal, V. 2025. Order-optimal global convergence for actor-critic with general policy and neural critic parametrization. In The 41st Conference on Uncertainty in Artificial Intelligence. Ganesh, S.; Mondal, W. U.; and Aggarwal, V. 2025a. OrderOptimal Regret with Novel Policy Gradient Approaches in Infinite Horizon Average Reward MDPs. In The 28th International Conference on Artificial Intelligence and Statistics. Ganesh, S.; Mondal, W. U.; and Aggarwal, V. 2025b. A Sharper Global Convergence Analysis for Average Reward Reinforcement Learning via an Actor-Critic Approach. In Proceedings of the 42nd International Conference on Machine Learning, 18206–18227. Gaur, M.; Bedi, A.; Wang, D.; and Aggarwal, V. 2024. Closing the Gap: Achieving Global Convergence (Last Iterate) of Actor-Critic under Markovian Sampling with Neural Network Parametrization. In Forty-first International Conference on Machine Learning. Ghosh, A.; Zhou, X.; and Shroff, N. 2023. Achieving sublinear regret in infinite horizon average reward constrained mdp with linear function approximation. In The Eleventh International Conference on Learning Representations. Gonzalez, G.; Balakuntala, M.; Agarwal, M.; Low, T.; Knoth, B.; Kirkpatrick, A. W.; McKee, J.; Hager, G.; Aggarwal, V.; Xue, Y.; et al. 2023. Asap: A semi-autonomous precise system for telesurgery during communication delays. IEEE Transactions on Medical Robotics and Bionics, 5(1): 66–78. Haliem, M.; Mani, G.; Aggarwal, V.; and Bhargava, B. 2021. A distributed model-free ride-sharing approach for joint matching, pricing, and dispatching using deep reinforcement learning. IEEE Transactions on Intelligent Transportation Systems, 22(12): 7931–7942. Hendrycks, D.; and Gimpel, K. 2016. Gaussian error linear units (gelus). arXiv preprint arXiv:1606.08415. Jacot, A.; Gabriel, F.; and Hongler, C. 2018. Neural tangent kernel: Convergence and generalization in neural networks. Advances in neural information processing systems, 31. Jin, Y.; and Sidford, A. 2021. Towards tight bounds on the sample complexity of average-reward MDPs. In International Conference on Machine Learning, 5055–5064. Khodadadian, S.; Chen, Z.; and Maguluri, S. T. 2021. Finitesample analysis of off-policy natural actor-critic algorithm. In International Conference on Machine Learning, 5420– 5431.
Liu, Y.; Zhang, K.; Basar, T.; and Yin, W. 2020. An improved analysis of (variance-reduced) policy gradient and natural policy gradient methods. Advances in Neural Information Processing Systems, 33: 7624–7636. Mondal, W. U.; and Aggarwal, V. 2024. Sample-efficient constrained reinforcement learning with general parameterization. Advances in Neural Information Processing Systems, 37: 68380–68405. Panju, M.; Raghu, R.; Sharma, V.; Aggarwal, V.; and Ramachandran, R. 2021. Queueing theoretic models for uncoded and coded multicast wireless networks with caches. IEEE Transactions on Wireless Communications, 21(2): 1257–1271. Papini, M.; Binaghi, D.; Canonaco, G.; Pirotta, M.; and Restelli, M. 2018. Stochastic variance-reduced policy gradient. In International conference on machine learning, 4026– 4035. Patel, B.; Suttle, W. A.; Koppel, A.; Aggarwal, V.; Sadler, B. M.; Manocha, D.; and Bedi, A. 2024. Towards Global Optimality for Practical Average Reward Reinforcement Learning without Mixing Time Oracles. In Proceedings of the 41st International Conference on Machine Learning, 39889– 39907. Paternain, S.; Chamon, L.; Calvo-Fullana, M.; and Ribeiro, A. 2019. Constrained reinforcement learning has zero duality gap. Advances in Neural Information Processing Systems, 32. Satheesh, A.; Barman, P. K.; Mondal, W. U.; and Aggarwal, V. 2026. Global Convergence of Average Reward Constrained MDPs with Neural Critic and General Policy Parameterization. In Proceedings of the 42nd Conference on Uncertainty in Artificial Intelligence (UAI). Sutton, R. S.; McAllester, D.; Singh, S.; and Mansour, Y. 1999. Policy gradient methods for reinforcement learning with function approximation. Advances in neural information processing systems, 12. Tamboli, D.; Chen, J.; Jotheeswaran, K. P.; Yu, D.; and Aggarwal, V. 2024. Reinforced sequential decision-making for sepsis treatment: The posnegdm framework with mortality classifier and transformer. IEEE Journal of Biomedical and Health Informatics, 28(5): 3114–3122. Tian, H.; Olshevsky, A.; and Paschalidis, Y. 2023. Convergence of actor-critic with multi-layer neural networks. Advances in neural information processing systems, 36: 9279– 9321. Wei, H.; Liu, X.; and Ying, L. 2022. A provably-efficient model-free algorithm for infinite-horizon average-reward constrained markov decision processes. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 36, 3868–3876. Xu, P.; Gao, F.; and Gu, Q. 2020. Sample Efficient Policy Gradient Methods with Recursive Variance Reduction. In International Conference on Learning Representations. Xu, T.; Liang, Y.; and Lan, G. 2021. Crpo: A new approach for safe reinforcement learning with convergence guarantee. In International Conference on Machine Learning, 11480– 11491.
Xu, Y.; Ganesh, S.; Mondal, W. U.; Bai, Q.; and Aggarwal, V. 2025. Global Convergence for Average Reward Constrained MDPs with Primal-Dual Actor Critic Algorithm. In The Thirty-ninth Annual Conference on Neural Information Processing Systems.
Algorithm 1: HiMLMC-PD-NAC Require: θ0 , λ0 = 0, α, β, γw , γξ (·), K, Hw , Hmax , Tmax 1: for k = 0, . . . , K − 1 do 2: for u ∈ {r, c} do 3: Compute critic (Qu,k , ηu,k ) ← Algorithm 2 (θk , u, Hmax , Tmax ) 4: end for 5: Compute NPG estimate wk ← Algorithm 3 (θk , λk , Qr,k , Qc,k , Hw , Tmax ) 6: Update primal variable θk+1 ← θk + αwk 7: Update dual variable λk+1 ← Π[0,2/δ] (λk − βηc,k ) 8: end for 9: return {θk , λk }K−1 k=0
Algorithm 2: Hierarchical MLMC Critic Require: θk , u, Hmax , Tmax 1: Sample Hk ∼ Geom(1/2) 2: hk = 1 + (2Hk − 1)1{2Hk ≤Hmax } k,(0)
k,(1)
k,(2)
3: Initialize ξu,0 , ξu,0 , ξu,0 4: for h = 0, . . . , hk − 1 do 5: Sample Uk,h ∼ Geom(1/2) and compute ℓk,h using (18) ℓ
k,h k,h Generate trajectory T k,h := (sk,h t , at )t=0 using current policy πθk , k,h Construct MLMC critic gradient vbu,MLMC using (20) if h < 1 then k,(0) k,(1) k,(2) Update ξu,h , ξu,h , and ξu,h , using (21) with stepsizes γξ (1), γξ (hk /2), γξ (hk ) end if if h < hk /2 then k,(1) k,(2) 12: Update ξu,h and ξu,h using (21) with stepsizes γξ (hk /2) and γξ (hk ) 13: end if k,(2) 14: Update ξu,h (21) with stepsize γξ (hk ) 15: end for 16: Compute critic function Qu,k,MLMC (·, ·) and parameter ξu,k,MLMC using (16) and (17) 17: return (Qu,k,MLMC (·, ·), ηu,k,MLMC )
6: 7: 8: 9: 10: 11:
Algorithm 3: MLMC NPG Estimator Require: θk , λk , Qr,k , Qc,k , Hw , Tmax 1: Initialize w0k 2: for h = 0, . . . , Hw − 1 do 3: Sample Uk,h ∼ Geom(1/2) and compute ℓk,h using (18) k,h ℓk,h 4: Generate trajectory T k,h := (sk,h t , at )t=0 using current policy πθk , kh kh 5: Construct MLMC Fisher estimate FbMLMC and MLMC PG estimate gbu,MLMC using (23) k 6: Update NPG estimate wh using (24) 7: end for k 8: return wH w
A
Using the definition of L(θ, λ) and L(π ∗ , λ), we have
Bounding Optimality Gap: Proof of Lemma 1
K−1 K−1 ∗ 1 X 1 X π∗ E Jr − Jrθk = − E λk Jcπ − Jcθk K K
Using the fact that λ0 = 0, we have for all K > 0, λ2K =
K−1 X
k=0
λ2k+1 − λ2k
k=0
+
k=0 K−1 (a) X
=
Π[0,2/δ] [λk − βηc,k,MLMC ]
2
− Π[0,2/δ] (λk )
2
k=0
X (b) K−1 ≤
1 K
K−1 X
E L(π ∗ , λk ) − L(θk , λk ) .
k=0
We combine the above relation with (26) to complete the proof of Lemma 1.
2
(λk − βηc,k,MLMC ) − λ2k
k=0
≤ β2
K−1 X
2 ηc,k,MLMC − 2β
k=0 (c)
≤ 2β 2
K−1 X
K−1 X
ηc,k,MLMC − Jcθk
2
+ 2β 2
K−1 X
Jcθk
2
We have from the dual update in (15), that
k=0
λk Jcθk − ηc,k,MLMC
+ 2β
∗ λk Jcπ − Jcθk ,
k=0
where (a) uses from the dual update in (15), (b) follows since any projection is non-expansive, and (c) follows since ∗ Jcπ ≥ 0, and the fact that c(s, a) ≤ 1, for all s, a implies ηc,k ≤ 1. Since λK ≥ 0, dividing both sides by 2βK and taking conditional expectation gives us K−1
−
1 X π∗ λk Jc − Jcθk K k=0
K−1
≤β+
1 X λk Jcθk − Ek ηc,k,MLMC K k=0
K−1
+
β X 2 Ek Jcθk − ηc,k,MLMC 2K
≤ β+
K−1 2 X θk Jc − Ek ηc,k,MLMC δK
Summing both sides over 0 ≤ k < K, multiplying with β/K, and taking expectation gives us
k=0
+ (b)
β 2K
2 ≤ β+ δK
K−1 X
Ek Jcθk − ηc,k,MLMC
2
k=0 K−1 X
∗ ξc,k − Ek ξc,k,MLMC
K−1 1 X 2 E Jcθk λk − K δ
k=0
k=0
K−1
K−1
β X 2 ∗ + Ek ξc,k − ξc,k,MLMC , 2K
2
≤ 2β 2 Ek Jcθk − ηc,k,MLMC 2 + 2β λk − Jcθk − Ek ηc,k,MLMC δ 2 θk − 2βJc λk − + 2β 2 (Jcθk )2 δ (b) 2 ∗ = O β 2 Ek ξc,k − ξc,k,MLMC 2 θk ∗ 2 λk − + β ξc,k − Ek ξc,k,MLMC − βJc +β , δ where (a) follows from non-expansiveness of the projection ∗ Π[0,2/δ] , and (b) follows by using the definitions of ξc,k and θk ξc,k,MLMC , and using λk ≤ 2/δ and Jc ≤ 1.
k=0
(a)
2 2 2 2 Ek λk+1 − − λk − δ δ 2 2 (a) 2 2 ≤ Ek λk − β ηc,k,MLMC − − λk − δ δ 2 2 = β 2 Ek ηc,k,MLMC − 2βEk ηc,k,MLMC λk − δ
k=0 K−1 X
Bounding Constraint Violation Rate: Proof of Lemma 2
k=0
K−1 X k=0
+ 2β
B
λk ηc,k,MLMC
(26)
k=0
where (a) follows since |Jcθk | ≤ 1, 0 ≤ λk ≤ 2/δ, and (b) ∗ follows from the definitions of ξc,k and ξc,k,MLMC .
1 β X 2 ∗ + E ξc,k − ξc,k,MLMC βK K k=0 ! K−1 X 1 ∗ + E ξc,k − Ek ξc,k,MLMC . (27) K
=O β+
k=0
C
Now, note that i 2 K−1 X ∗ Jcθk Jrπ − Jrθk − δ k=0 k=0 K−1 K−1 i h X ∗ 2 (a) X π∗ θk π θk ≤ Jr − J r + Jc − Jc δ k=0 k=0 K−1 K−1 X 2 (b) X ∗ Jcθk . = [L(π , λk ) − L(θk , λk )] + λk − δ
K−1 Xh
k=0
k=0
∗
where (a) uses the fact that Jcπ ≥ 0, while (b) uses the definition of the Lagrangian L. Dividing both sides by K, taking expectation, and substituting the last term using (27) then gives us ∗
Jrπ −
K−1 K−1 1 X θk 2 X θk E Jr + E −Jc K δK k=0
=O β+
k=0
1 β + βK K
1 + K
K−1 X
K−1 X
∗ − ξc,k,MLMC E ξc,k
2
k=0
k=0
k=0
We define a policy π̄K that picks a policy from {πθk }k<K , uniformly at random. Then, (28) implies 2 ∗ Jrπ − E Jrπ̄K + E − Jcπ̄K δ K−1 1 β X 2 ∗ =O β+ E ξc,k − ξc,k,MLMC + βK K k=0
1 K
K−1 X
∗ − Ek ξc,k,MLMC E ξc,k
k=0
! K−1 1 X ∗ + L(π , λk ) − L(θk , λk ) . K k=0
By (Xu et al. 2025, Lemma G.6), we conclude K−1 1 X E − Jcπ̄K = E − Jcθk K k=0
K−1
=O
δ δβ X 2 ∗ + δβ + E ξc,k − ξc,k,MLMC βK K k=0
+
δ K
∗
∗
Let expectation w.r.t stationary distributions dπ or ν π be denoted by E∗ [·]. Then, we have for for 0 ≤ k < K, E∗ KL (π ∗ (·|s)||πθk (·|s)) − KL π ∗ (·|s)||πθk+1 (·|s) πθ (a|s) = E∗ ln k+1 πθk (a|s) (a)
≥ E∗ ⟨∇θ ln πθk (a|s), θk+1 − θk ⟩ −
G2 ∥θk+1 − θk ∥2 2
G2 α2 ∥wk ∥2 2 G2 α2 ≥ αE∗ ⟨∇θ ln πθk (a|s), wk ⟩ − ∥wk ∥2 2 G2 α2 (c) = α [L(π ∗ , λk ) − L(θk , λk )] − ∥wk ∥2 2 + αE∗ ⟨∇θ ln πθk (a|s), wk ⟩ h ∗ i h ∗ i − α Jrπ − Jrθk − αλk Jcπ − Jcθk (b)
= αE∗ ⟨∇θ ln πθk (a|s), wk ⟩ −
G2 α 2 ∥wk ∥2 2 + αE∗ ⟨∇θ ln πθk (a|s), wk − wk∗ ⟩ + αE∗ ⟨∇θ ln πθk (a|s), wk∗ ⟩ − Aθrk (s, a) − λk Aθck (s, a) (d)
∗ − Ek ξc,k,MLMC E ξc,k
! K−1 1 X ∗ + L(π , λk ) − L(θk , λk ) . (28) K
+
Lagrange Error Analysis: Proof of Lemma 3
K−1 X
∗ E ξc,k − Ek ξc,k,MLMC
k=0
! K−1 δ X ∗ + L(π , λk ) − L(θk , λk ) . K k=0
This completes the proof of Lemma 2.
= α [L(π ∗ , λk ) − L(θk , λk )] −
√ G2 α2 ∥wk ∥2 − α εbias , 2 + αE∗ ⟨∇θ ln πθk (a|s), wk − wk∗ ⟩ , (e)
≥ α [L(π ∗ , λk ) − L(θk , λk )] −
where (a) follows from Assumption 4, (b) follows from (15), (c) follows from (3), (d) follows from Lemma 5, and (e) follows √ from combining Assumption 6 with the fact that x 7→ x is concave and Assumption 6. Now, rearranging terms and taking the conditional expectation w.r.t Fk gives √ 1 Ek [L(π ∗ , λk ) − L(θk , λk )] − εbias − α × E∗ KL (π ∗ (·|s)||πθk (·|s)) − KL π ∗ (·|s)||πθk+1 (·|s) (a)
G2 α ∥wk ∥2 2 (b) ≤ G1 ∥Ek [wk ] − wk∗ ∥ + G2 α E∥wk − wk∗ ∥2 + E∥wk∗ ∥2 ≤ − ⟨Eπ∗ [∇θ ln πθk (a|s)] , Ek [wk ] − wk∗ ⟩ +
(c)
2 ≤ G1 ∥Ek [wk ] − wk∗ ∥ + G2 αE∥wk − wk∗ ∥2 + O(τmix α),
where (a) follows since θk is Fk -measurable and wk is inde∗ pendent of ν π , (b) uses the inequality (a + b2 ) ≤ 2a2 + 2b2 , while (c) follows as we use Assumptions 2, 4, and 5 to get wk∗ = F (θk )† ∇θ L(θk , λk ) 1 2 ≤ ∇θ Jrθk + ∇θ Jcθk , µF δ and use Assumption 1 to conclude ∇θ Juθk ≤ G1 Eπθk Aθuk (s, a) = O(τmix ).
Finally, the result follows by summing both sides over 0 ≤ k < K, dividing by K, using the fact that the KL-divergence is non-negative, and taking the expectation.
D
NPG Estimator Analysis: Proof of Theorem 2
Throughout the subsequent sections, we let Ek,0 [·] denote the conditional expectation Ek [·]. Further, for h > 0, we let Ek,h [·] denote the conditional expectation w.r.t Fk,h , i.e., all samples obtained upto the k-th outer-loop iteration and the subsequent (h − 1) inner-loop iterations. Recall from (24), that the NPG update is given by the following linear stochastic recursion: k wh+1 = whk h i k,h k,h k,h + γw ĝr,MLMC + λk ĝc,MLMC − F̂MLMC whk , (29) k,h k,h k,h where ĝr,MLMC +λk ĝc,MLMC and F̂MLMC are MLMC-based estimators for the Lagrangian gradient ∇θ L(θk , λk ) and the Fisher matrix F (θk ), respectively. Note that, under Assumption 5, wk∗ is the unique vector satisfying F (θk )wk∗ = ∇θ L(θk , λk ). We claim that the NPG update satisfies the conditions needed in Lemma 7, which bounds the convergence rate achieved by a general linear stochastic recursion. Note that, by Assumptions 4 and 5, for all θ ∈ Θ and λ ∈ [0, 2/δ],
∥∇θ L(θ, λ)∥ ≤
2 4G1 τmix . 1+ δ
Next, we bound the bias and MSE of MLMC estik,h k,h mators ĝu,MLMC and F̂MLMC . Equation (22) shows the PT −1 kh naive estimator for F (θk ) as T1 t=0 F̂ k (skh t , at ) = ⊤ ∇θ ln πθk (a|s) (∇θ ln πθk (a|s)) . From Assumption 4, we have i h Eθk F̂ kh (s, a) − F (θk ) = 0 and
ĝuk,h (s, a, ā, s′ , a′ ) ≤ G1 u(s, a) + ηu,k,MLMC + 2 Qu,k,MLMC = O ηu,k,MLMC + Qlin (·, ·, ζu,k,MLMC ) + Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC ) = O ηu,k,MLMC + ζu,k,MLMC + Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC ) . Thus, the single-sample MSE for the naive PG estimator is 2
k,h k,h k,h k,h θk Ek ĝuk,h (sk,h t , at , āt , st+1 , at+1 ) − ∇θ Ju 2 ∗ 2 = Õ σQ + ξu,k,MLMC − ξu,k ,
where we define 2 σQ := Ek Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC )
µF ≤ ∥F (θ)∥ ≤ G21 and
Here, Qu,k,MLMC denotes the critic estimator for πθk defined in (17). Recall from (14) that Qlin be the first-order Taylor approximation to Q(·, ·, ζ). Then, we have
kh F̂ kh (skh t , at ) − F (θk )
2
≤ 4G41 ,
kh Thus, by Lemma 8, we we have the bias and MSE of F̂MLMC as h i 2 kh δF2 := Ekh F̂MLMC − F (θk ) = O (τmix /Tmax ) , (30) 2 kh σF2 := Ekh F̂MLMC − F (θk ) = O (τmix ln Tmax ) . k,h Likewise, we bound the MSE and bias in ĝu,MLMC , recall θk from (22) that the naive estimator for ∇θ Ju is given by PT −1 k,h k,h k,h k,h k,h k,h 1 t=0 ĝu (st , at , āt , st+1 , at+1 ), where T ĝuk,h (s, a, ā, s′ , a′ ) = u(s, a) − ηu,k + Qu,k,MLMC (s′ , a′ ) − Qu,k,MLMC (s, ā) ∇θ ln πθ (a|s).
2
.
Now, we bound the single-sample bias for the naive PG estimator. Note that, Eθk ĝuk,h (s, a, ā, s′ , a′ ) = Eθk ∇θ ln πθk (s, a) h × u(s, a) − Juθk + Qθuk (s′ , a′ ) − Qθuk (s, ā) ∗ + Juθk − ηu,k,MLMC + Qθk (s, ā) − Qlin (s, ā, ζu,k ) ∗ − Qθk (s′ , a′ ) − Qlin (s′ , a′ , ζu,k ) ′ ′ ∗ + Qlin (s , a , ζu,k,MLMC ) − Qlin (s′ , a′ , ζu,k ) ∗ − Qlin (s, ā, ζu,k,MLMC ) − Qlin (s, ā, ζu,k ) + Qu,k,MLMC (s, ā) − Qlin (s, ā, ζu,k,MLMC ) i − Qu,k,MLMC (s , a ) − Qlin (s , a , ζu,k,MLMC ) . (31) ′
′
′
′
Applying the NTK bounds from Lemma 6 on (31) and using (8), we have the single-sample PG bias as Eθk
ĝuk,h (s, a, ā, s′ , a′ )
2
− ∇ θ J θk
1 2 ∗ = Õ εapp + + ξu,k,MLMC − ξu,k m + Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC )
2
and 2 k,h ′ ′ Ek Eθk ĝu (s, a, ā, s , a ) − ∇θ J θk 1 2 ∗ = Õ εapp + + Ek ξu,k,MLMC − ξu,k m 2 . + Ek Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC )
Consequently, applying Lemma 8, we get the following MSE k,h and bias for ĝu,MLMC : σg2 :=
X
k,h − ∇θ Juθk Ek,h ĝu,MLMC
2
u
X 1 2 2 ∗ + σQ ξu,k,MLMC − ξu,k , + = Õ εapp + m u (32) and δg2 :=
X
2 k,h − ∇θ Juθk Ek,h ĝu,MLMC
u
X 1 2 ∗ 2 ξu,k,MLMC − ξu,k + σQ + , = Õ εapp + m u X 2 k,h δ̄g2 := Ek ĝu,MLMC − ∇θ Juθk
u
X 1 2 ∗ 2 Ek ξu,k,MLMC − ξu,k = Õ εapp + + δQ + , m u (33) where the sum is over u ∈ {r, c} and 2 2 . δQ := Ek Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC ) Using these bounds, we are ready to apply Lemma 7 to the NPG update (29). Choosing the NPG stepsize γw := µ2FlnHTw , we have for T sufficiently large: 1 1 ∗ 2 2 2 2 2 Ek wk − wk = Õ + σ + σ g + δF + δg , T2 Hw F 1 2 + δF2 + δ̄g2 Ek [wk ] − wk∗ = Õ T2 1 2 2 2 2 2 + δP σ + σ g + δF + δg . Hw F Finally, substituting the bias and MSE terms in the above relations with (30), (32), and (33) completes the proof of Theorem 2.
Q(·, ·, ζ) with the linear critic Qlin (·, ·, ζ) in (19), we decompose the naive single-sample critic-gradient estimator as h i k v̂ k,h (z, z ′ ) := Âk,h (z, z ′ ) ξu,h − B̂uk,h (z, z ′ ) − Ṽuk,h (z, z ′ ), where Âk,h (z, z ′ ) := B̂uk,h (z, z ′ ) :=
Critic Estimator Analysis: Proofs of Lemma 4 and Theorem 3
Proof of Lemma 4. We begin by writing the vanilla-critic update as a perturbed linear recursion. Let z, z ′ denote the state action pairs (s, a), (s′ , a′ ), respectively. Further, let ψ(z) denote ∇ζ Q(z, ζ0 ). Then, replacing the neural critic
cγ ψz
0 ⊤ , ψ(z) [ψ(z) − ψ(z ′ )]
cγ u(z) , (u(z) − Q(z, ζ0 ) + Q(z ′ , ζ0 )) ψ(z)
and Ṽuk,h (z, z ′ ) collects the remainder terms. This decomposition then translates to the following decomposition of the MLMC critic-gradient estimator in (20): h i k,h k,h k,h k v̂MLMC := Âk,h (35) MLMC ξu,h − B̂u,MLMC − Ṽu,MLMC , k,h k,h where Âk,h MLMC , B̂u,MLMC , and Ṽu,MLMC are obtained by replacing v̂uk,h in (20) with Âk,h , B̂uk,h , and Ṽuk,h , respectively. With the decomposition (35), the vanilla-critic update (21) can be written as the following perturbed linear recursion:
k,h k k k ξu,h+1 = ΠB ξu,h + γξ,H B̂u,MLMC − Âk,h MLMC ξu,h k,h + γξ,H Ṽu,MLMC . (36) ∗ Recall that ξu,k is the unique solution to A(θk )ξ = Bu (θk ), where
A(θk ) := Eθk [Âk,h (z, z ′ )], and
Bu (θk ) := Eθk [B̂uk,h (z, z ′ )].
(37)
k,h k,h Let Âk,h MLMC , B̂u,MLMC , and Ṽu,MLMC satisfy the following MSE bounds 2 2 Ek,h ∥Âk,h MLMC − A(θk )∥ ≤ σA , k,h 2 Ek,h ∥B̂u,MLMC − Bu (θk )∥2 ≤ σB ,
and
k,h ∥2 ≤ σV2 , ∥Ṽu,MLMC
(38)
and the following bias bounds 2 2 ∥Ek,h [Âk,h MLMC ] − A(θk )∥ ≤ δA ,
and
E
(34)
k,h 2 ∥Ek,h [B̂u,MLMC ] − Bu (θk )∥2 ≤ δB .
(39)
If we define the auxiliary linear recursion χku,h+1 k,h k = ΠB χku,h + γξ,H B̂u,MLMC − Âk,h , MLMC χu,h
(40)
k where χku,0 = ξu,0 , then subtracting (36) from (40) gives us k ξu,h+1 − χku,h+1 (a)
≤ =
(a)
=
2
2 k k,h I − γξ,H Âk,h ξu,h − χku,h − γξ,H Ṽu,MLMC MLMC k I − γξ,H A(θk ) ξu,h − χku,h 2 k,h + γξ,H A(θk ) − Âk,h MLMC − γξ,H Ṽu,MLMC
k I − γξ,H A(θk ) ξu,h − χku,h
2
2
2
k,h 2 2 A(θk ) − Âk,h + γξ,H Ṽu,MLMC + γξ,H MLMC E D k,h 2 − 2γξ,H A(θk ) − Âk,h MLMC , Ṽu,MLMC D + 2γξ,H A(θk ) − Âk,h MLMC , k E × I − γξ,H A(θk ) ξu,h − χku,h , (41)
where (a) follows from the non-expansiveness of ΠB . Taking conditional expectation on both sides, we have 2
k Ek,h ξu,h+1 − χku,h+1 k 2 (a) = I − γξ,H A(θk ) ξu,h − χku,h 2
k,h 2 2 + γξ,H Ek,h A(θk ) − Âk,h + γξ,H Ek,h Ṽu,MLMC MLMC D E k,h 2 − 2γξ,H Ek,h A(θk ) − Âk,h MLMC , Ṽu,MLMC D + 2γξ,H A(θk ) − Ek,h [Âk,h MLMC ], k E × I − γξ,H A(θk ) ξu,h − χku,h (b)
≤ (1 − µϕ γξ,H )
2
2
2
2 2 σA + σV2 + 2γξ,H
Hence, to conclude Lemma 4, it remains to bound the MSEs k,h k,h and biases of Âk,h MLMC , B̂u,MLMC , and Ṽu,MLMC . By Lemma 6 and (37), the corresponding single-sample biases satisfy Eθk [Ṽuk,h (z, z ′ )] = Õ(R m−1/2 ), Eθk [Âk,h (z, z ′ )] − A(θk ) = 0, Eθk [B̂uk,h (z, z ′ )] − Bu (θk ) = 0,
(46)
while under NTK regime, the single-sample MSEs satisfy
2
= O(c2γ ), (47)
2 σA = O(τmix ln Tmax ),
2 −1 δA = O(τmix Tmax ),
2 σB = O(τmix ln Tmax ),
2 −1 δB = O(τmix Tmax ).
(48)
Finally, we substitute the above bounds into (44) and (45), and use the inequality
k ∗ Ek ξu,H − ξu,k
2
k ≤ exp (−µϕ γξ,H ) ξu,h − χku,h 2 2 + O γξ,H σA + σV2 + γξ,H R δA ,
(45)
Applying Lemma 8, we conclude σV2 = Õ(τmix m−1 ) and
k + 2γξ,H (1 − µϕ γξ,H ) ξu,h − χku,h δA (d)
2
∗ Ek χku,H − ξu,k 1 2 2 2 2 = Õ σA + σB + δA + δB + σV2 . H
B̂uk,h (z, z ′ ) − Eθk [B̂uk,h (z, z ′ )] = Õ(c2γ ).
D + 2γξ,H A(θk ) − Ek,h [Âk,h MLMC ], k E × I − γξ,H A(θk ) ξu,h − χku,h k ≤ (1 − µϕ γξ,H ) ξu,h − χku,h
∗ Taking the NTK radius R sufficiently large so that ξu,k lies in the NTK ball, we apply the same recipe used to obtain (44) to the auxiliary update (40). This gives us
Âk,h (z, z ′ ) − Eθk [Âk,h (z, z ′ )]
2 A(θk ) − Âk,h MLMC 2 k,h 2 + 2γξ,H Ek,h Ṽu,MLMC
2 + 2γξ,H Ek,h
(c)
T Taking R = Θ(ln T ) and γξ,H = 2µln , and using the fact ϕH k k that ξu,0 = χu,0 , 1 2 k 2 2 Ek ξu,H − χku,H = Õ σA + σV2 + δA , (44) H
Ṽuk,h (z, z ′ ) − Eθk [Ṽuk,h (z, z ′ )] = Õ(R m−1/2 ),
2
k ξu,h − χku,h
k (39). Finally, (d) follows since ξu,H , χku,H lie in the NTK ball and by using the inequality 1 − x < e−x . Then, applying (42) recursively and following the arguments in the proof of (Xu et al. 2025, Theorem B.1), we get 2 2 k k Ek ξu,H − χku,H = O Ek ξu,0 − χku,0 e−µϕ Hγξ,H 2 2 + γξ,H σA + σV2 + R δA . (43)
2
k ≤ 2Ek ξu,H − χku,H
(42)
k where (a) follows since ξu,h , χku,h are Fk,h -measurable, (b)
follows since assuming ker(A(θk )) = {0}, (this holds if the vectors {∇ζ Q(s, a, ζ0 )} do not span the all ones vector in R|S||A| ) we can show that for cγ sufficiently large, ξ ⊤ A(θk )ξ ≥ µϕ ∥ξ∥2 (see (Satheesh et al. 2026, Lemma 17)), and using Cauchy-Schwarz inequality, and (c) follows by choosing γξ,H < µϕ /2 and using the bounds (38) and
2
∗ + 2Ek χ∗u,k − ξu,k
2
.
(49)
This completes the proof of Lemma 4. Next, we conclude this section with the proof of Theorem 3. Proof of Theorem 3. Recall from (16) and (17), that k Qu,k,MLMC (·, ·) := Q ·, ·, ζu,1 h i k k + 1{2Hk ≤Hmax } 2Hk Q ·, ·, ζu,2 − Q ·, ·, ζu,2 Hk Hk −1
where (a) follows from (52), (b) follows since ∥ · ∥ is convex, and (c) follows from Lemma 4. Now, we bound the critic MSE. Using a telescoping sum as earlier, we have the critic parameter MSE as
and ξu,k,MLMC k k k , := ξu,1 + 1{2Hk ≤Hmax } 2Hk ξu,2 Hk − ξ Hk −1 u,2
∗ Ek ξu,k,MLMC − ξu,k
where Hk ∼ Geom(1/2). Then, we have k Ek Qu,k,MLMC (·, ·) − Ek Q ·, ·, ζu,1 (a)
2 ⌊log2 Hmax ⌋
(a)
k ∗ ≤ 2Ek ξu,1 − ξu,k
⌊log2 Hmax ⌋
X
=
2 k k × Ek ξu,2 h − ξu,2h−1
h=1 k ∗ ≤ 2Ek ξu,1 − ξu,k
⌊log2 Hmax ⌋
=
h
2 ·2
⌊log2 Hmax ⌋
(b)
h h i i k k − Ek Q ·, ·, ζu,2 × Ek Q ·, ·, ζu,2 h h−1 X
−h
×
h=1 2 k ∗ 2Ek ξu,2h − ξu,k + 2Ek
×
h i k k = Ek Q ·, ·, ζu,2 − Ek Q ·, ·, ζu,1 , (50) ⌊log2 Hmax ⌋
(c)
h=1
where (a) uses (16) and (b) follows since Hk ∼ Geom(1/2). Hence, h i k Ek Qu,k,MLMC (·, ·) = Ek Q ·, ·, ζu,2 . (51) ⌊log2 Hmax ⌋ A similar telescoping argument gives us h i k Ek ξu,k,MLMC = Ek ξu,2 ⌊log2 Hmax ⌋ .
X
k ∗ + = 2Ek ξu,1 − ξu,k
k k Ek Q ·, ·, ζu,2 − Ek Q ·, ·, ζu,2 h h−1
=
22h+1 · 2−h k ∗ ξu,2 h−1 − ξu,k
2h+2
h=1 2 2 k ∗ k ∗ Ek ξu,2h − ξu,k + Ek ξu,2 h−1 − ξu,k
⌊log2 Hmax ⌋ h X 2 Hmax = Õ , = Õ 1 + m m
where (a) uses (17), (b) follows since Hk ∼ Geom(1/2), and (c) follows from Lemma 4. Similar arguments show the following: 2
(a)
(52)
(53)
h=1
Ek Qu,k,MLMC − Qlin (·, ·, ζu,k,MLMC )
2 Thus, we have the critic-function bias δQ h i 2 = Ek Qu,k,MLMC (·, ·) − Qlin (·, ·, ζu,k,MLMC )
2
k ≤ 2Ek Q(·, ·, ζu,1 ) − Qlin (·, ·, ζu,k,MLMC ) ⌊log2 Hmax ⌋
X
+
h=1
2 · 4h 2 k k Ek Q(·, ·, ζu,2 h ) − Q(·, ·, ζu,2h−1 ) 2h
k = 2Ek Q(·, ·, ζu,1 ) − Qlin (·, ·, ζu,k,MLMC )
2 k k = Ek Q ·, ·, ζu,2 − Qlin ·, ·, ζu,2 ⌊log2 Hmax ⌋ ⌊log2 Hmax ⌋
(a)
2
⌊log2 Hmax ⌋
X
+
(b)
k k 2h+1 Ek Q(·, ·, ζu,2 h ) − Q(·, ·, ζu,2h−1 )
2
h=1
= Õ (1/m) , (b)
where (a) follows from (51) and (52), together with the fact that Qlin (·, ·, ζ) is linear in ζ. Whereas, (b) follows from the NTK bounds in Lemma 6 since ζu,2⌊log2 Hmax ⌋ is contained within the NTK ball B. Further, the critic-parameter bias satisfies 2 ∗ Ek ξu,k,MLMC − ξu,k 2 k (a) ∗ = Ek ξu,2 ⌊log2 Hmax ⌋ − ξu,k (b)
2
⌊log2 Hmax ⌋
h h i i k k × Ek Q ·, ·, ζu,2 − Ek Q ·, ·, ζu,2 h h−1 ⌊log2 Hmax ⌋
X
+
h=1
X
2 · 4h · P(Hk = h)
h=1
2h · P(Hk = h)
(b)
X
+
2
k ∗ ≤ Ek ξu,2 ⌊log2 Hmax ⌋ − ξu,k τmix τmix τmix 1 (c) + + ⌊log H ⌋ + √ = Õ T2 Tmax m 2 2 max τmix τmix τmix = Õ + + , Tmax Hmax m
k ≤ 2Ek Qlin (·, ·, ζu,1 ) − Qlin (·, ·, ζu,k,MLMC )
2
⌊log2 Hmax ⌋
X
+
k k 2h+1 Ek Qlin (·, ·, ζu,2 h ) − Qlin (·, ·, ζu,2h−1 )
2
h=1
+ Õ (c)
1 m
⌊log2 Hmax ⌋ X 1+ 2h+1 h=1
k = Õ Ek ζu,1 − ζu,k,MLMC
2
+
Hmax m
⌊log2 Hmax ⌋
X
+
k k 2h+1 Ek ζu,2 h − ζu,2h−1
2
h=1
Hmax + = Õ m
(d)
⌊log2 Hmax ⌋
X h=1
h 2 Õ , m
(54)
where (a) is obtained from (16) and using Hk ∼ Geom(1/2), k k k to obtain (b), for each ζ ∈ ζu,1 , ζu,2 h , ζu,2h−1 , we replace Q(·, ·, ζ) with Qlin (·, ·, ζ) and bound ∥Q(·, ·, ζ)−Qlin (·, ·, ζ)∥ using the NTK bounds in Lemma 6, (c) follows since Qlin (·, ·, ζ) is Lipschitz in ζ, and (d) uses (53) and Lemma 4. Therefore, we conclude that Hmax 2 σQ = Õ . m
We finally combine the bounds derived in Theorems 2 and 3, and Lemmas 3, 2, and 1 to prove Theorem 1. We choose m = Hmax and Hmax = Tmax = T, as prescribed in Theorem 1. Then, Theorem 3 gives the following critic bounds:
and
and
2 δQ = Õ (1/T ) ,
(55)
2
∗ Ek ξu,k,MLMC − ξu,k = Õ (1) , 2 ∗ Ek ξu,k,MLMC − ξu,k = Õ (τmix /T ) .
(56)
Plugging the bounds from (55) and (56) into Theorem 2 and taking Hw = Θ(ln T ), we get 2 Ek wk − wk∗ = Õ τmix + εapp , (57) τmix 2 . Ek [wk ] − wk∗ = Õ εapp + T Next, we substitute the above NPG √ bias and MSE bounds into Lemma 3 and choose α = 1/ K and K = Θ(T ) to get the following Lagrange error: K−1 1 X E L(π ∗ , λk ) − L(θk , λk ) K k=0 √ τmix √ 1 √ 2 =O εbias + εapp + τmix α + + √ αK T 2 √ τmix √ = Õ εbias + εapp + √ . (58) T
Finally, we plug the√ bounds from (58) and (56) into Lemma 1 and choose β = 1/ T . This gives us the following averaged error: K−1 2 √ τmix 1 X θk √ π∗ εbias + εapp + √ . Jr − E Jr = Õ K T k=0 Likewise, plugging (58) and (56) into Lemma 2 gives the constraint violation rate: K−1 2 √ 1 X θk τmix √ − E Jc = Õ εbias + εapp + √ . K T k=0 This concludes the proof of Theorem 1.
The following lemma is the standard performance difference lemma (Bai, Mondal, and Aggarwal 2024b, Lemma 4) for average-reward MDPs. Lemma 5 (Performance Difference Lemma). For any two policies πθ and πθ′ , and u ∈ {r, c}, we have ′
′
Juπ − Juπ ≤ E(s,a)∼ν π [Aθu (s, a)], ′
Proof of Theorem 1
2 σQ = Õ (1)
Technical lemmas
where Aθu is the advantage function of policy πθ′ as defined in (9).
This completes the proof of Theorem 3.
F
G
The following lemma, proved in (Satheesh et al. 2026, Lemma 15), collects several useful bounds for neural critics under the NTK regime. Lemma 6 (NTK bounds). Choose R = O(ln T ) and let k ξu,h ∈ B, for all 0 ≤ h < hk . Then, the following holds with probability at least (1 − ε − 2Le−Cm ) : p (a) ∥∇ζ Q(·, ·, ζ)∥ ≤ C1 ; |Q(·, ·, ζ)| = O ln(hk /ε) ; √ ln(hk /ε) √ (b) ∥Q(·, ·, ζ) − Qlin (·, ·, ζ)∥ = O ; m √ ln(hk /ε) √ ; (c) ∥∇ζ Q(·, ·, ζ) − ∇ζ Q(·, ·, ζ0 )∥ = O m √ ln(hk /ε) √ (d) ∥Ṽuk,h (·, ·, ·, ·)∥ = O . m The following lemma, proved in (Ganesh, Mondal, and Aggarwal 2025b, Theorem 2), characterizes the bias and MSE of a stochastic linear recursion. Lemma 7 (Stochastic Linear Recursion). Consider the linear stochastic recursion xh+1 = xh + γ q̂h − P̂h xh , d where (q̂h , P̂h )H−1 h=0 are noisy estimates for (q, P ) ∈ R × Rd×d . The target is an x∗ such that P x∗ = q. Suppose that (q̂h , P̂h ) satisfies:
Eh ∥P̂h − P ∥2 ≤ σP2 ,
Eh ∥q̂h − q∥2 ≤ σq2 ,
∥Eh [P̂h ] − P ∥2 ≤ δP2 ,
∥Eh [q̂h ] − q∥2 ≤ δq2 .
Further, let ∥E[q̂h ] − q∥2 ≤ δ̄q2 ,
λP ≤ ∥P ∥ ≤ ΛP ,
∥q̂h ∥ ≤ Λq ,
and
γ≤
λP . (24σP2 + 8Λ2P )
Then, after H iterations, we have E∥xH − x∗ ∥2 ≤ ∥x0 − x∗ ∥2 e−γλP H + O (γR0 + R1 ) ∥E[xH ] − x∗ ∥2 ≤ ∥x0 − x∗ ∥2 e−β̄λP H + δP2 + O δP2 (γR0 + R1 ) + R̄1 , where R0 = O σP2 + σq2 , R1 = O δP2 + δq2 , and R̄1 = O δP2 + δ̄q2 .
The following lemma is adapted for biased estimates from (Dorfman and Levy 2022, Lemma 3.1), and summarizes the computational cost, bias, and MSE guarantees of MLMC estimation. Lemma 8 (MLMC bounds). Let (Zt ) be a timehomogeneous ergodic Markov chain with a unique invariant distribution dZ , and a mixing-time τmix . Assume that F (x, Z) is an estimate of F (x). Let for all t ≥ 0, ∥EdZ [F (x, Z)] − F (x)∥2 ≤ δ 2 and
E∥F (x, Zt ) − EdZ [F (x, Z)]∥2 ≤ σ 2 .
Further, let Q ∼ Geom(1/2). Then, the MLMC estimator defined as gMLMC := g1 + 12Q ≤Tmax · 2Q (g2Q − g2Q−1 ) , Pj−1 where gj := j −1 t=0 F (x, Zt ), satisfies the following bounds. (a) E[gMLMC ] = E[g2⌊log2 Tmax ⌋ ]; (b) E∥F (x) − gMLMC ∥2 = O σ 2 τmix log2 Tmax + δ 2 ; −1 (c) ∥F (x) − E[gMLMC ]∥2 = O σ 2 τmix Tmax + δ2 .