Thank you for your time and detailed review. Please see our responses below.

A discussion on the effects of scaling state and action space might be nice to include

The bound in Theorem 5.1 only depends on the size of the action space, and only scales logarithmically with this number. Importantly, the bound does not depend directly on the size of the state space at all, only on the dimension of the feature map $d$. We will be sure to add this discussion in the revision.

Insight into how the current work may be extended to the general non-linear function approximation setting would make the paper comprehensive

The main obstacle in generalizing our algorithm in its exact form to more general function approximation is the confidence ellipsoid, $\lVert \xi \rVert_{\hat{\Lambda}} \leq \beta$ . This is based on martingale concentration techniques from the linear MDP literature (Lemma B.1). These techniques allow us to very precisely quantify the level of uncertainty in the dataset. It is not yet clear how to do this for more general function approximation like neural networks.

However a promising avenue may be to adapt techniques from papers such as [Xie et al., 2021] and [Cheng et al., 2022]. Specifically, one can consider removing the pessimism parameter $\xi$ and considering a Lagrangian relaxation of the critic's problem. Specifically, in episode $k$ we may ask the critic to find a function $q_k\in \mathcal{Q}$ minimizing $J + \frac{\lambda}{N} \sum_{i = 1}^N (q(s_i, a_i) + J - r(s_i, a_i) - E_{a \sim \pi_k}[q(s_i', a)])^2,$ where $\mathcal{Q}$ is some suitable class for function approximation, $\lambda$ is a Lagrange multiplier which can be tuned to control the degree of pessimism, and $\frac{1}{N} \sum_{i = 1}^N (q(s_i, a_i) + J - r(s_i, a_i) - E_{a \sim \pi_k}[q(s_i', a)])^2$ is the average square TD error. The main challenge we see here, which is not present in the discounted setting, is the presence of the additional parameter $J$ meant to estimate the average reward of policy $\pi_k$. Ideally, $J$ should depend on $q$ in such a way such that $J(q_k)$ is a pessimistic estimate of the average reward of policy $\pi_k$. We leave this interesting problem for future work.

While weaker than uniform ergodicity, the unichain assumption may not always hold. As mentioned in the section on future directions, it would be interesting to see how this affects algorithm design and analysis.

There are two main difficulties we come across without the unichain assumption. The first is that, in this case, $J^{\pi}$ not be constant and could depend the initial state. The formulation of the optimization problem (6) may need to change to reflect this. The second is that, without Assumption 3.1, the true Q-functions may be unbounded. The stability of our algorithm relies on the members of the function class $\mathcal{Q}(B_w)$ being bounded. If the true Q functions are unbounded, but the elements of $\mathcal{Q}(B_w)$ are, it may be unreasonable to expect the completeness and approximate realizability assumptions to hold. We will add the discussion in the revision.

Thank you very much for your time and detailed review. Please see our responses to your questions below.

Clarity

There are a few typos in the proof. Thank you for pointing these out!

The $w_k^*$ appearing in the proof of Lemma B.3 is a typo. It should be $\bar{w}_k$ throughout the proof. We have corrected this.

The second inequality in line 610 should read $\frac{1}{K} \sum_{k = 1}^{K} J^{\pi} - \hat{J_k} + \kappa_{B_w, B_{\theta}} \leq 2B_w \sqrt{\frac{\log A}{K}} + 2\beta \lVert\phi^{\mu^{\pi}}\rVert_{\hat{\Lambda}^{-1}} + \varepsilon_{B_w, B_{\theta}} + \kappa_{B_w, B_{\theta}}.$ We have corrected this as well.

We will also update section 4.2 to be clear that the policy update follows mirror ascent with KL divergence regularization.

Originality

However, it seems that the idea of proof procedure (Lemma 6.2 and 6.3) is similar to the FOPO algorithms' proof in [Wei et al, 2021]. The only difference is that [Wei et al, 2021] discusses online problems with Optimism, and this paper discusses offline problems with pessimism.

Although our analysis is partly inspired by [Wei et al., 2021], we would like to point out some notable differences below.

1. Unlike [Wei et al, 2021], we do not require the perfect linear MDP assumption. Specifically, we only assume the approximate realizability and restricted closedness properties. Additional work is required for the analysis to handle this level of generality. For instance, in the proof of Lemma 6.2, we cannot assume that $w_k^{\star}$ parameterizes the true Q-function $q^{\pi_k}$ exactly. So the function $q_k^{\star}(s, a) = \phi(s, a)^{\top}w_k^{\star}$ need not be a solution to the Bellman equation. This results in the additional term $\hat{\Lambda}^{-1} \sum_{i = 1}^{N} \phi(s_i, a_i) \Delta_{k}(s_i,a_i)$ in equation (17), where $\Delta_k(s, a) = q_k^{\star}(s, a) + J_k^{\star} - T^{\pi_k}q_k^{\star}(s, a)$ is the irreducible Bellman error due to realizability holding only approximately. When arguing that $w_k^{\star}, J_k^{\star}$ is feasible for the optimization problem (6) with high probability, this term needs to be accounted for. This is why the additional factor of $\kappa_{\mathcal{Q}(B_w), \Pi} \sqrt{N}$ appears in the definition of the confidence parameter $\beta$. A similar error term also appears in Lemma 6.3 due to the restricted closedness assumption.

2. Because our algorithm is policy based, our analysis is built on the extended performance difference lemma (Lemma 6.1). This differs from the proof structure in [Wei et al., 2021] since their algorithm takes greedy actions. More specifically, the way we decompose the sub-optimality in step 3 of the proof sketch is quite different. We decompose the sub-optimality into an optimization term and a bias term. Our Lemma 6.3 focuses on controlling the bias in the estimate of the Q-function for each policy $\pi_k$ encountered during the algorithm's execution. Such a statement does not appear in [Wei et al., 2021] because thier algorithm is designed to approximate the solution to the Bellman optimality equation in each episode, rather than the solution for a fixed policy. They then need to argue that the bias in their estimate of the optimal Q function is not too large. Therefore, their regret analysis is based entirely on the Bellman optimality equation and martingale difference concentration rather than a performance difference lemma.

The author mentioned that [Gabbianelli et al., 2024] is the only paper that also discusses offline RL under average MDP. However, they only achieve O(\eps^{-4}). Could the author explain what the new idea or technique is adopted to improve the result to O(\eps^{-2})?

Thank you for the question. To explain our new technique, we first briefly summarize the main ideas of [Gabbianelli et al., 2024] using their notation. In their work, they assume a linear MDP structure, realized by $\Psi \in R^{d \times |S|}$, $w \in R^d$ satisfying $r(s, a) = \langle \phi(s, a), w \rangle, \qquad P(s'|s, a) = \langle \phi(s, a), \psi(s') \rangle,$ where $\psi(s')$ is the $s'$ column of $\Psi$. To adapt to average reward case, they also assume there is some $\varrho \in R^d$ with $\langle \phi(s, a), \varrho\rangle = 1$. Assuming that state action pairs are drawn i.i.d from some distribution $\mu_{B}$, the population covariance is $\Lambda = E_{\mu_B}[\phi(s, a)\phi(s, a)^{\top}]$. Under their linear MDP assumption, they observe the average reward (they call $\rho^{\pi}$) can be written as $\rho^{\pi} = \langle \lambda^{\pi}, w \rangle$ where $\lambda^{\pi} = E_{(s, a)\sim \mu^{\pi} \otimes \pi}[\phi(s, a)]$. They then formulate the problem as a linear program with Lagrangian $f(\beta, \pi; \rho, \theta) = \rho + \langle \beta, \Lambda[w + \Psi v_{\theta, \pi} - \theta - \rho \varrho]\rangle.$ Here $v_{\theta, \pi}(s) = \sum_{a}\pi(a|s)\langle \phi(s, a), \theta \rangle$ and $\beta = \Lambda^{-1}\lambda$ is a reparameterization used so that knowledge of $\Lambda$ is not necessary to run the algorithm. Their main idea is to solve the saddle point problem $\min_{\rho, \theta}\max_{\beta, \pi} f(\beta, \pi; \rho, \theta).$ Note, importantly, that this is a population level problem. Since they only assume access to i.i.d. samples from $\mu_{B}$, they solve the problem via stochastic gradient descent-ascent.

The main bottleneck in the analysis of [Gabbianelli et al., 2024] is their algorithm's double-loop structure that, in each outer loop $t$, uses an inner-loop to solve $\min_{\rho, \theta} f(\beta_t, \pi_t; \rho, \theta)$ nearly exactly using $O(\varepsilon^{-2})$ samples.Then as they need $O(\varepsilon^{-2})$ outer-loop total iterations, they could only get the $O(\varepsilon^{-4})$ upper bound.

The key reason for the sample complexity improvement in our paper is the efficient use of the entire dataset for the critic's optimization problem (6). Namely, the data to construct this problem only needs to be sampled once and is then re-used in each iteration of our algorithm. This avoids the inner-loop that uses additional samples. The only issue that needs to be dealt with is the statistical correlation between iterates of the algorithm arising from data re-use. This correlation is handled in the analysis using the $\epsilon$-net covering technique.

We would also like to point out that the linear program structure of [Gabianelli et al., 2024] crucially relies on the linear MDP structure. Our use of the constrained optimization approach is the key reason we can relax the assumptions to cover restricted closedness and approximate realizability.

The authors have clearly answered my question (1) the difference that I didn't notice compared with [Wei et al, 2021] (2) The reason why the proposed algorithm achieved much faster convergence rate than [Gabbianelli et al., 2024].

I highly encourage the authors to explain the advantage over [Gabbianelli et al., 2024] in their revised version. I will keep my positive score.

We thank the reviewer for their time and positive review. Please see our response below.

However, I would like to highlight that the algorithm is not fully computationally efficient. It is true that Algorithm 1 can be efficiently implemented with the help of an interior-point method to solve the optimization problem in step 5. While interior-point methods are highly precise, they are not exact. That said, the theoretical guarantees provided in the paper hold under the assumption that there exists an optimization solution which provides an exact solution to step 5. I strongly recommend that the authors add a comment to clarify this, especially in the abstract and other parts of the main text with explicit claims about computational efficiency.

Thank you for raising this point. While we assumed the use of exact solutions, inexact solutions still only incur an additive error proportional to the accuracy of the solution. Specifically, suppose that the optimization problem (6) is solved to within accuracy $\delta$, so that $J_k$ is within $\delta$ of the optimal solution and the constraint violation is also at most $\delta$. The computational complexity of many existing algorithms to obtain such a solution solution scale as $\log\delta^{-1}$ [Nesterov and Nemirovskii, 1994]. Then additional error will be incurred in Lemma 6.2 where we will instead have $J_k \leq J_k^* + \delta \le J^{\pi_k} + \kappa_{\mathcal{Q}(B_w), \Pi} + \delta,\quad \forall k \in [K].$ There will also be an error incurred in Lemma 6.3 due to the constraint violations, which can be updated to read $|E_{(s, a) \sim \mu} [q_k(s, a) - \phi(s, a)^{\top}\bar{w_k}]| \leq 2(\beta + \delta) \lVert\phi^{\mu}\rVert_{\hat{\Lambda}^{-1}} + \delta.$ Here the factor of $2\delta$ that is a coefficient of $\lVert\phi^{\mu}\rVert_{\hat{\Lambda}^{-1}}$ is due to the violation of the quadratic constraints and the additive factor is from a violation of the linear constraint. We will add a clarifying statement in the main text accordingly.

Do the authors assume that the initial state is fixed?

We do not need to assume that the initial state is fixed. It may be random. Under assumption 3.1, the average reward $J^{\pi}$ does not depend on the initial state. Therefore, the identity of the initial state, or distribution if it is random, does not change the problem or our results.

Thank you for pointing out the typos! We will correct these in the revision.

[Yurii Nestorov and Arkadii Nemirovskii, 'Interior-Point Polynomial Algorithms in Convex Programming', Society for Industrial and Applied Mathematics, 1994]

We thank the reviewer for the positive feedback.

The main obstacle in directly replacing the linear critic with a neural network is the method we use to enforce pessimism. Specifically, the variable $\xi$ in (6) is constrained to a confidence ellipsoid $\lVert \xi \rVert_{\hat{\Lambda}} \leq \beta$. The parameter $\beta$ quantifies the uncertainty in the dataset. This is determined precisely using martingale concentration results from the linear MDP literature. It is not yet clear how to extend this technique to more general models such as neural networks.

There are methods in the literature that implement pessimism for more general approximation. However, they need to sacrifice the precise uncertainty quantification that can be achieved in the linear setting. Two relevant papers mentioned in the related work section are those of [Xie et al., 2021] and [Cheng et al., 2022], which focus on the discounted setting. Algorithm 1 in [Xie et al., 2022] alternates between policy optimization steps and a critic which chooses a pessimistic approximation model from a general function class. When specialized to the case of linear function approximation, their algorithm can be viewed as solving a Lagrangian relaxation of a constrained problem similar to ours, but with only the linear constraint. Since they work with more general function approximation, they are not able to quantify the uncertainty in the dataset to the same level of precision as we are in our paper. Instead, the amount of pessimism is determined by the Lagrange multiplier, which is a tunable hyper-parameter. The algorithm is very general, but when specified to the linear setting, the lack of precise uncertainty quantification results in convergence rates of only $N^{-1/3}$ for $N$ data samples.

It would be interesting to see if similar ideas work for average reward MDPs, but it does not seem trivial due to the presence of the additional variable $J$ in the average reward setting. We leave these interesting questions for future work.