We sincerely thank the reviewer for their feedback and valuable suggestions. We address each point carefully and explicitly below:

Regarding Numerical Illustrations (Weaknesses and Q2)

We would like to note that this work is primarily theoretical, and the full empirical evaluation of robust TD learning is left as a future work for application papers. That said, because our guarantees hinge on the one‑step semi‑norm contraction properties (characterized in Lemma 3.1 and Theorem 3.2), we now include numerical evaluations that directly verify strict contraction across the settings studied. These results empirically support the key structural claim used by our analysis. Due to space constraint, we refer to our additions in the response to Reviewer rBB6 (Regarding Numerical Validations).

Regarding Discounted Robust RL Methods (Q1)

We make the following comparisons between our work and the prior robust discounted RL works. In discounted RL, contraction comes for free from the discount factor $\beta<1$ in a (weighted) sup-norm [1], and standard stochastic fixed point iterations yields rates accordingly. In contrast, average-reward has no inherent contraction and our analysis builds a new semi-norm under which the robust average-reward operator is a contraction via ergodicity. The algorithmic analysis is therefore fundamentally different from discounted works, and we will add a short remark clarifying this correspondence and the different constants.

Regarding Robust Level $\delta$ (Q3)

Thanks for bringing up this important point. We note that the bound for contamination (Eq. 138) and Wasserstein (Eq. 141) do not have a dependance over $\delta$, and thus the limit as $\delta\to 0$ is the same. Thus, the result do not degrade with increasing $\delta$. However, we note that contraction factor depends on $\delta$ since the maximum spectral radius in the uncertainty set will increase with $\delta$. However, there is an additional factor of $(1+\frac{1}{\delta})$ for the TV uncertainty set (see Eq. 128, 139), which is a loose bound and is left as a future work to make this tighter. We will mention the discussion in the final version.

Regarding Lower Bounds (Q4)

We note that since the non-robust serves as a special case for the robust policy evaluation, the lower bound of non-robust can serve as a hard example for the robust setting. Nevertheless, obtaining tight lower bounds that also track $\delta$ for robust evaluation remains an interesting open problem. We will mention the above in the revised version.

Regarding Typos and the Missing Definition of $\Delta$ (Q5-Q7 and Paper Formatting Concerns)

Thanks for identifying the typos, and we will fix them in the revised draft. In particular, we will define the simplex $\Delta(\mathcal S)$ as the probability simplex over the state space and use this notation consistently.

References

[1] Weighted Sup-Norm Contractions in Dynamic Programming: A Review and Some New Applications, 2012.

We sincerely thank the reviewer for their feedback and valuable suggestions. We address each point carefully and explicitly below:

Regarding Optimality and Lower Bounds (Weakness 1)

We agree our rates are order‑optimal only in $\epsilon$ ($\tilde O(\epsilon^{-2})$). Regarding other factors (e.g., $|\mathcal S|, |\mathcal A|$ and structural constants), we will clarify that our bounds have a linear $|\mathcal S||\mathcal A|$ factor from per–state–action oracle calls, and also depend on the contraction gap $1-\gamma$, the span bound on $V$, and an additional $\log(1/\epsilon)$ from truncation in the TV and Wasserstein cases.

Regarding the lower bounds, we note that classical regret lower bound for finding an optimal policy in non-robust average reward reinforcement learning in [1] is $\Omega(\sqrt{|S| |A| T D})$ where $D$ denotes the diameter of the MDP. However, this does not directly translate into a policy‑evaluation sample‑complexity lower bound of $\Omega(\sqrt{S A D} \epsilon^{-2})$. This is because the conversion from regret in control to estimation error in evaluation is non‑trivial and depends on additional structure (e.g., mixing time/diameter and whether a generative model is available). To the best of our knowledge, there have been no works specifically characterizing the exact lower bound for policy evaluation in (robust) average reward reinforcement learning with the focus on $|S|$ and $|A|$. We will therefore not claim tightness in $|\mathcal S| |\mathcal A|$ and treat sharpening these dependencies as open. Regarding the lower-bound justifying the claim that the dependence on $\epsilon$, we note that we provide the standard mean estimation as a hard example and demonstrate the sample complexity of $\Omega(\epsilon^{-2})$ in Appendix F.2. We will mention the above in our revised draft.

Regarding Section 5 (Weakness 2)

We note that Section 5 is a bit concise due to page constraints, and we will expand it in the final version due to availability of an additional page. We have provided the detailed analysis for this Section in Appendix F, some of the key details can be moved to Section 5 with additional space availability in the final version. We would add the discussion of the additional logarithmic terms in the sample complexity results for TV and Wasserstein distance uncertainty sets. These are due to the bias of the estimators for the support functions $\hat{\sigma}(V)$ in Algorithm 1. In contrast, the estimator of the support function in the contamination uncertainty set is unbiased (as shown in Eq. (19)), and thus the corresponding sample complexity for policy evaluation not have the logarithmic terms.

We also agree that a brief discussion on how these policy evaluation methods could be fit in within a procedure to find the optimal policy is important, specifically, our algorithm can be used to estimate the robust policy sub-gradients, which can be then utilized to perform policy optimizations accordingly, we will add the above in the revised draft.

Regarding Rectangularity (Q1)

Thanks for raising up this point. We note that our algorithms do focus on contamination, TV, and Wasserstein ambiguity around a nominal model $P$. We will clarify the presentation by first fixing $P$, and then defining the $(s,a)$‑rectangular set as $\mathcal U=\prod_{s,a}\mathcal U_{s,a}$ where $\mathcal U_{s,a}$ is, respectively, a contamination/TV/Wasserstein neighborhood of $P(\cdot\mid s,a)$. In our algorithms, the nominal generative model first supplies samples from $P$, and then we compute the support function $\min_{q\in\mathcal U_{s,a}}\langle q, V\rangle$ by decoupling across $(s,a)$, this procedure is achievable due to rectangularity assumption. Thus, our $(s,a)$‑rectangular assumption makes Algorithm 1 implementable (closed‑form for contamination; truncated‑MLMC for TV/W) and is standard in robust DP to preserve Bellman‑style decomposability. We will add the above to the revised draft.

Regarding Motivations for Different Uncertainty Sets (Q2)

We provide the motivations for the different types of uncertainty sets studied in this work as follows:

• Contamination ($\delta$). Models outliers or rare faults [2]: with prob. $\delta$, the next‑state law can be arbitrary (e.g., sensor glitch or unmodeled jump to a bad region). Leads to a simple worst‑case expectation combining $P$ with the worst state value.
• Total variation (TV radius). Models categorical misspecification or discretization error [3]: bounded $\ell_1$ mass shift per $(s,a)$; captures sparse but possibly large reallocation of probability mass.
• Wasserstein (metric radius). Models smooth model drift when states have a geometry [4] (e.g., position/velocity): mass can move to nearby states with cost proportional to distance; realistic for learned simulators and dynamics linearization.

We will add a compact paragraph with these examples in the revised draft

Regarding Typos (Q3)

Thanks for identifying the typos, we agree with the review on the definition of $E$. The typos will be fixed in the revised draft.

Regarding Model-Based Methods (Q4)

We agree that our focus is model‑free setting. We note that the model-based counterpart is not straightforward, and there is a new paper [5] in ICML 2025 (which came online after our work was submitted), that studies robust model-based average reward setting with TV, chi-square, and contamination uncertainty sets. We will add this to Related Work in the final version, while note that the approach of model-based that is used in this paper which is reduction from discounted to average does not work well in model-free settings due to higher dependence on $1/(1-\gamma)$ in the model-free setups.

Regarding $g$ and Theorem 2.1 (Q5)

We agree $g^\pi$ is a scalar independent of the starting state under a standard unichain/ergodicity assumption. Our Theorem 2.1 indeed assumes ergodicity later in the section. Thanks for pointing this out. we will move that assumption next to the theorem statement to avoid ambiguity in the revised draft.

Regarding Line 18 in Algorithm 2 (Q6)

Thanks for pointing out this detail. In our algorithm, we use uniform averaging to control the span/semi‑norm uniformly over states. The analysis also goes through with stationary‑distribution weighting ($d_{P}$); we will add a remark noting this variant and that it can reduce variance in practice without changing the rate.

References

[1] Near-optimal regret bounds for reinforcement learning, 2008.

[2] A General Decision Theory for Huber’s $\epsilon$-Contamination Model, 2017.

[3] Partial Policy Iteration for $L_1$-Robust Markov Decision Processes, 2021.

[4] First-Order Methods for Wasserstein Distributionally Robust MDPs, 2021.

[5] A Reduction Framework for Distributionally Robust Reinforcement Learning under Average Reward, 2025.

We sincerely thank the reviewer for their feedback and valuable suggestions. We address each point carefully and explicitly below:

Regarding Numerical Validations (Weakness 1)

We like to note that this work is primarily theoretical, and the full empirical evaluation of robust TD learning is left as a future work for application papers. That said, because our guarantees hinge on the one‑step semi‑norm contraction properties (characterized in Lemma 3.1 and Theorem 3.2), we now include numerical evaluations that directly verify strict contraction across the settings studied. These results empirically support the key structural claim used by our analysis.

Evaluation of Lemma 3.1

Lemma C.2 is the technical backbone for Lemma 3.1, as Lemma C.2 constructs the fixed‑kernel seminorm and provides the one-step contraction for a given $P$. So we perform numerical evaluations on Lemma C.2 to demonstrate the one-step contraction property for Lemma 3.1.

For an ergodic Markov matrix $P$ with stationary $d$, we follow the steps in Appendix C.1 and construct $||\cdot||_P$ in Eq. (36) with $\alpha=\min(0.99,(1+\rho(Q))/2)$ and $\epsilon=0.25(1-\alpha)$ where $Q=P-\mathbf ed^\top$. Then the one‑step contraction factor would be $\beta=\alpha+\epsilon<1$.

To generate ergodic matrices with dimension $n$, let $I_n$ be the identity matrix and $S_n$ be the cyclic shift matrix $(S_ne_i=e_{i+1 \bmod n})$. We provide the following 4 examples:

• $P_1 = 0.5I_5 + 0.5S_5$
• $P_2 = 0.6I_6 + 0.4S_6$
• $P_3 = 0.55I_7 + 0.45S_7$
• $P_4 = 0.6I_8 + 0.3S_8 + 0.1S_8^2$

We generate 1000 random unit vectors $x$ and calculate each of the ratio $\frac{||P_i x||_P}{||x||_P}$. We present the empirical results as follows:

Evaluation of Theorem 3.2

Lemma C.4 is the key step for Theorem 3.2 as Lemma C.4 proves a uniform one‑step contraction across all $P$ in the uncertainty. So we perform numerical evaluations on Lemma C.4 to demonstrate the one-step contraction property for Theorem 3.2 under contamination, TV, and Wasserstein-1 distance uncertainty set. We select $P_1$, $P_2$, $P_3$, and $P_4$ as 4 examples of the nominal model.

To perform numerical approximation of $\lVert\cdot\rVert_{\mathcal P}$ defined in Eq. (50), we approximate $\lVert\cdot\rVert_{\mathcal P}$ by (i) discretizing the uncertainty set and (ii) using a finite‑product to approximate the extremal‑norm. First, we sample a family $\{P^{(i)}\}_{i=1}^m\subset\mathcal P$

of size $m$ and form $Q_i:=P^{(i)}-\mathbf e d_{P^{(i)}}^\top$ with $\hat r=\max_i \rho(Q_i)$. We set $\alpha=\min(0.99,(1+\hat r)/2)$ and $\epsilon\in(0,1-\alpha)$. To approximate the extremal norm, we build a library of scaled products of the $Q_i$’s up to maximum length $K$ for each $k=0,1,\dots,K$ we draw a set of products $M_{k,j}=Q_{i_k}\cdots Q_{i_1}$; the count of such draws at each $k$ is the products per length (denoted as samples_per_k in the table). This defines the surrogate $||z||_{\rm ext}^{(K)}$

=$\max_{0\le k\le K,j}\ \alpha^{-k}||M_{k,j}z||_2$. We then set

$||x||_{\mathcal P}^{(K)}$

$=\max_{i=1,\dots,m} ||Q_i x||_{\rm ext}^{(K)}$

$+\epsilon\min_{c}||x-c \mathbf e||_{\mathrm{ext}}^{(K)},$

We generate 50 random unit vectors $x$ for each of the sampled uncertainty matrix inside

$\{P^{(i)}\}_{i=1}^m\subset\mathcal P$

, and calculate each of the ratio of

$||P^{(i)} x||_{\mathcal{P}}$

and $||x||_{\mathcal{P}}$

We present the empirical results as follows:

Contamination

TV

Wasserstein-1

Interpretations

Note that for the above tables, max span ratio denotes the empirical largest one‑step span contraction coefficient over the sampled families. This value equals to $1$ for all the settings, meaning no strict one‑step contraction in the span.

The quantities ratio_min / ratio_median / ratio_p90 / ratio_max summarize the empirical one‑step ratios computed under the constructed semi-norms over all sampled $P$ and random unit directions $x$. They report, respectively, the minimum, median, 90th percentile, and maximum observed value across those tests. In every table we have ratio_max < 1, so we observe empirical one‑step contraction under our semi-norms even when span does not contract. Moreover, ratio_max $\le$ $\gamma$ (robust case) or $\le$ $\beta$ (non-robust case), which is consistent with the corresponding theoretical contraction factor guaranteed by our constructions.

The above numerical results and the code for the experiments would be released in the revised manuscript.

Regarding Assumption 2.3 (Weakness 2)

Thanks for raising this question. We provide Remark 2.4 after this Assumption 2.3 showing how it relaxes to ergodicity of the nominal model only when the uncertainty radius is small. The assumption of ergodicity for nominal model aligns with common practice in the non‑robust average‑reward RL literature (such as [1]-[3]). Additionally, In our numerical validations, we also assume only the nominal model is ergodic and still observe strict contraction for reasonably large radii, in our case, $\delta=0.15$.

Regarding Algorithm 1's Sample Complexity (Q1)

We would like to clarify that in Algorithm. 1, the expected sample cost of a single call to the robust‑expectation oracle is $O(N_{\max})$ by Theorem 4.1, and in Appendix F4 we set $N_{\max}={O}(\log(1/\epsilon))$ to control bias/variance. Hence we note that Algorithm 1 does not require exponentially many i.i.d. samples; the overall policy‑evaluation complexity remains $\tilde O(\epsilon^{-2})$. Due to space constraint, we refer to our response to Reviewer uoff (Regarding MLMC) for a more detailed explanation of why we do not need exponentially many i.i.d. samples. We will mention the above in the revised draft to prevent this misunderstanding.

Regarding Typos and Definition of $\mathbf e$ (Q2-Q4)

Thanks for identifying the typos, we agree with the review on the definition of $\sigma(V)$ and the duplicated reference, and they will be fixed in the revised draft. We note that $\bf e$ is initially defined in line 144, we will also re-emphasize this definition in Lemma 3.1 in the revised draft.

Regarding Limitations

We agree that our tabular analysis and the use of a generative model are simplifying assumptions. Our goal here is to establish the first finite‑sample guarantees for policy evaluation in robust average‑reward RL. To the best of our knowledge, the most closely related works cited in Appendix A also study tabular problems and adopt a generative‑model abstraction. This abstraction is particularly natural in the robust setting, where the motivation is assuming access to the simulated environment. We have noted these limitations in the conclusion. We view these as promising directions but outside the scope of this first finite‑sample study.

References

[1] Finite Sample Analysis of Average-Reward TD Learning and Q-Learning, 2021.

[2] Learning Infinite-horizon Average-reward MDPs with Linear Function Approximation, 2021.

[3] A sharper global convergence analysis for average reward reinforcement learning via an actor-critic approach, 2025.

Dear Reviewer rBB6,

Thank you again for your thoughtful and constructive comments. We wanted to follow up to see if you had any additional questions or concerns following our rebuttal. We’d be glad to provide further clarification during the discussion phase if needed.

Regards,

Authors

We sincerely thank the reviewer for their feedback and valuable suggestions. We address each point carefully and explicitly below:

Regarding MLMC (Q1, Q5 and Weaknesses)

Thanks for bringing up comparisons between our work and the two Blanchet papers. We note that [1] obtains $\tilde O(\epsilon^{-2})$ using unbiased and untruncated MLMC under a smooth‑functional assumption (differentiable with Hölder‑continuous derivative, as mentioned in their Theorem 1), which enables second‑order multilevel cancellation. [2] studies discounted robust RL with KL divergence uncertainty sets, where these conditions hold, so their unbiased‑MLMC analysis in [1] carries through. In contrast, we focus on TV and Wasserstein distance uncertainty sets, and the support functions in our setting are Lipschitz with no uniform differentiability guarantees. Thus, the analysis of the unbiased MLMC in [1] and [2] is no longer applicable. We therefore adopt a different framework where we apply the truncated MLMC method and still yielding the $\tilde O(\epsilon^{-2})$ sample complexity.

We note that we characterize the sample complexity, bias and variance properties of our truncated MLMC approach (Algorithm. 1) in our Theorem 4.1-4.4. In particular, we note that due to truncation, the expected number of samples to perform Algorithm. 1 is $O(N_{\rm max})$ instead of $O(2^{N_{\rm max}})$ with $N_{\rm max}$ being the maximum truncation level. Thus, we agree with the reviewer that if one uses bias to be $2^{-N/2}$, the level is indeed $O(\log(\epsilon^{-2}))$, which is equivalent to $O(\log(\epsilon^{-1}))$ (what we wrote in the draft). However, this will use $O(\log(\epsilon^{-2}))$ samples instead of $O(\epsilon^{-2})$ within one iteration with constant order probability, which still leads to our $\tilde{O}(\epsilon^{-2})$ sample complexity result.

Regarding Model-Based Methods (Weaknesses)

We agree that our focus is model‑free since the problem of policy evaluation is essentially a model-free approach by its nature, and we will add emphasis on this in the revised draft. We note that the problem is not straightforward even for model-based settings, and indeed there is a new paper [3] in ICML 2025 (which came online after our work was submitted), that studies robust model-based average reward setting with TV, chi-square, and contamination uncertainty sets. We will add this to Related Work in the final version, while note that the approach of model-based that is used in this paper which is reduction from discounted to average does not work well in model-free settings due to higher dependence on $1/(1-\gamma)$ in the model-free setups.

Regarding the Semi‑Norm Contraction and $\gamma$ (Q2-Q4 and Weaknesses)

Thanks a lot for the great suggestion on providing examples to demonstrate the contraction properties in Lemma 3.1 and Theorem 3.2. Due to space constraint, we refer to our additions in the response to Reviewer rBB6 (Regarding Numerical Validations), where we provide detailed numerical validations and demonstrate the semi-norm contraction properties.

Regarding the interpretation of $\gamma$. Note that in Lemma C.3, by setting $\epsilon \rightarrow 0$ and $\alpha \rightarrow r^\star$, our construction allows $\gamma$ to be made arbitrarily close to $r^\star = \sup_{P\in\mathcal P}\rho(Q_P)$ defined in Eq. (41). Since $r^\star$ represents the second‑largest eigenvalue among all the $P\in\mathcal P$, this semi-norm is constructive and $\gamma$ tracks the spectrum that matters for stability. We note that in the reversible (or lazy‑reversible) case, $r^\star$ directly controls mixing: for each $P$, where we have so $t_{\mathrm{mix}}=\tilde O\big(1/(1-r^\star)\big)$ (shown in Theorem 12.4 of [4]). For non‑reversible chains, the standard replacement is the pseudo‑spectral gap via multiplicative reversibilization; we will add a short remark connecting $r^\star$ to that quantity for completeness.

We also would like to note that the kernel of a semi-norm is the set of all vectors in the vector space that map to zero under the semi-norm, we will add this explanation in the revised draft to improve clarity.

Regarding Future Work (Q6 and Weaknesses)

We agree with the reviewer that the policy optimization and control problem is an important direction. However, policy optimization is based on policy evaluation, which is the focus of this paper. Thus, one cannot obtain policy optimization results without having a complete analysis of policy evaluation. Further, we note that the current policy optimization papers in robust average-reward setup [5] assume knowledge of exact sub-gradient, while our method provides actual sampling techniques with finite-time analysis for policy evaluation. Thus, our work can potentially be used in their analysis to replace the knowledge of exact sub-gradient. However, a detailed analysis is left as a future work, and will be mentioned in the paper.

Regarding the Mapping $P\mapsto\rho(Q_P)$ (Q7)

We note that the continuity of the mapping holds under the assumptions in the paper. On finite ergodic chains, the stationary distribution $d_P$ depends continuously on $P$; hence $Q_P=P-\mathbf e d_P^\top$ is continuous in $P$. The spectral radius is a continuous function of the matrix entries (Appendix D in [6]), so the composition $P\mapsto\rho(Q_P)$ is continuous on the ergodic region. We will add the detailed proof on this in the final version.

References

[1] Unbiased multilevel Monte Carlo: Stochastic optimization, steady-state simulation, quantiles, and other applications, 2019.

[2] A Finite Sample Complexity Bound for Distributionally Robust Q-learning, 2023.

[3] A Reduction Framework for Distributionally Robust Reinforcement Learning under Average Reward, 2025.

[4] Markov Chains and Mixing Times, 2017.

[5] Policy Optimization for Robust Average Reward MDPs, 2024.

[6] Matrix Analysis, 2012.

Thank you for your feedback and the corrections of the details of the MLMC approaches. We are glad to hear that our response addressed your major concerns.

Additionally, thanks for bringing up [01] and [02]. We note that according to NeurIPS policy, both [01] and [02] are considered concurrent to this paper (and were posted around the NeurIPS main paper deadline). Nevertheless, we agree that [02] is a model-based method while [01] utilizes a model-free value iteration approach. Specifically, [01] studies robust average reward RL under contamination and $l_p$ norm uncertainty sets, leveraging a value-iteration style Halpern iteration (RHI). Since [01] utilizes properties of non-expansiveness (along with some properties in the uncertainty sets) instead of contraction, [01] does not have the contraction factor in the sample complexity results. We also note that RHI is analyzed under the assumption that every stochastic increment it adds to the running Bellman estimate is unbiased. This is guaranteed by the recursive sampling procedure. For Contamination and $l_p$ norm uncertainty, the robust-Bellman update can be expressed as one next-state expectation plus a deterministic norm term, so a single draw gives an increment whose conditional mean is exactly the true difference, guaranteeing unbiasedness. However, for a Wasserstein ball, the correction is a global optimal-transport cost that no finite set of local samples can reproduce on average, so any stochastic estimate is inevitably biased. Thus, one limitation of [01] is that RHI fails under Wasserstein uncertainty sets, while our work does not have this limitation.

We will add a discussion of the above concurrent works in the revised draft.