1. Introduction. Reinforcement learning (RL) problems are generally formulated in terms of Markov decision processes (MDPs). At each time step, the agent observes the current state and subsequently takes an action, which leads to the realization of some reward as well as a transition to the next state according to the underlying, but unknown, stochastic transition function. The eventual goal of the agent is to learn a policy, i.e., a mapping from states to actions, to optimize the reward accrued over time. The setting is a very general one, with applications ranging from engineering to the natural and social sciences; see, e.g., [14,17] for surveys of RL applications.

A fundamental building block in RL is the problem of policy evaluation, in which we are interested in estimating the long-term (discounted) value of each state under a fixed policy with sample access to the transition and reward functions. The literature considers three observation models for transition and reward samples, namely the generative model, the socalled "i.i.d." model, and the Markovian noise model 1 . Furthermore, in modern applications with large state spaces, it is common to seek an approximation to the true value function within the span of a small number of basis functions, a setting that is commonly known as linear function approximation. In the canonical setting of the problem, one is interested in using convergence of this algorithm was shown by Tadic [42], who extended the convergence results for solving noisy linear systems [35].

While asymptotic convergence results offer a proof-of-concept, the algorithm is often run in large-scale applications with relatively small sample sizes. The first results proving finite time convergence under i.i.d. setting were proposed by Sutton et al. [41] and later extended by Lakshminarayanan and Szepesvári [21]. Finite-time analysis of TD learning under Markovian noise was carried out by Bhandari et al. [1], where the authors employed nonsmooth analysis to a variant of TD learning, by requiring projections at each iteration onto a pre-specified ball. A consequence of the nonsmooth approach is that there is no obvious way of benefiting from the variance reduction effect of parallel computing. In recent work, a subset of the current authors [20] provided an improved analysis of vanilla TD algorithm that overcomes this hurdle. There are also several other notable finite sample analyses of policy evaluation in various settings, e.g., [39,5,10,24,26], and statistical lower bounds have also been shown for offline reinforcement learning with linear function approximation, e.g., [47,51].

While some of these analyses are sharp, to our knowledge, vanilla TD learning is not known to attain the optimal oracle complexity and instance-dependent sample complexity. Li et al. [24] recently proved that TD learning achieves the minimax lower bound on stochastic error up to logarithmic factors, but their analysis is not instance dependent. Recent work by Kotsalis et al. [20] presented two new algorithms, the conditional temporal difference (CTD) and the fast temporal difference (FTD) learning, where FTD exhibits an accelerated rate in deterministic error. However, these algorithms fail to capture the correct stochastic error in the policy evaluation problem, and the bounds can be shown to be suboptimal for policy evaluation both in an instance-dependent sense, and in the worst-case over natural problem classes.

During the past decade, there has been a flurry of parallel work in stochastic optimization on developing first-order methods with variance-reduction; early examples include IAG [2], SAG [37], SVRG [13,49], and SAGA [7]. There are several papers that apply variance reduction to reinforcement learning, e.g., [8,34,45]. Recent work in policy evaluation has shown that variance reduction techniques can also be applied to algorithms of the TD-type [18,44,50,15]. Among these, the paper [15] is motivated by the desire to draw distinctions between RL algorithms with similar worst-case performance and follows a line of work deriving instance-dependent bounds on the stochastic error in policy evaluation [33,25]. Specifically, the results of [15] capture the optimal instance-dependent stochastic error in the ℓ ∞ -norm.

However, the optimal sample complexity is not achieved in [15] since the algorithm requires O{1/(1 -γ) 3 } samples in each epoch.

Towards answering the question posed at the end of Section 1, we make three distinct contributions:

• Lower bounds. We construct a worst-case instance that shows an oracle complexity lower bound of order Ω{(1 -γ) -1 • log(1/ϵ)} for any iterative method whose iterates lie within the linear span of the initial point v 0 and subsequent temporal differences, to converge to ϵ-error in weighted ℓ 2 -norm. We also prove a lower bound on sample complexity using the classical local minimax theorem [11,22,23] to provide an instance-specific baseline for algorithm design.

This manuscript is for review purposes only.

• Algorithm design in the i.i.d. setting. We start by applying the variance reduction technique to the classical TD algorithm, showing that the resulting variancereduced temporal difference (VRTD) algorithm nearly matches the optimal stochastic error, but the analysis suggests a suboptimal deterministic error. This motivates us to further improve the VRTD algorithm with the stochastic operator extrapolation (SOE) device [19]. We provide a sharp analysis of our new algorithm-termed variancereduced fast temporal difference (VRFTD)-showing that it achieves a convergence rate nearly matching both the deterministic error lower bound (for well-conditioned feature matrices) and the stochastic error lower bound.

• Extension to the Markovian setting. We extend the VRFTD algorithm to the Markovian setting by introducing a burn-in period during sample collection. We show that the resulting algorithm also achieves similarly fast convergence, with a dominating stochastic error term that matches the instance-dependent lower bound proved in [28] and a deterministic error that matches that of the i.i.d. setting up to a multiplicative factor of the mixing time. In particular, in the so-called realizable case when the approximation error caused by linear function approximation is 0 (e.g., in the tabular setting), the leading Markovian stochastic error term is equal to the i.i.d. stochastic error term, indicating that the additional dependence on mixing time only appears in terms whose dependence on the final tolerance ϵ is weak.

The rest of this paper is organized as follows. In Section 2, we formally present the problem setting. The three aforementioned main contributions are presented in Sections 3-5. In Section 6, we provide numerical experiments that corroborate our optimality guarantees. The proofs of our main results are postponed to Section 7, and auxiliary results are collected in the supplementary materials.

1.3. Notation. For a positive integer n, we define [n] := {1, 2, ..., n}. We let 1 denote the all-ones vector in R D . We let e j denote the j-th standard basis vector in R D . Let I S : X → {0, 1} denote the indicator function of the subset S ⊆ X. Given a vector x ∈ R m , denote its i-th entry byfoot_0 x (i) . Let

denote the ℓ 1 , ℓ 2 and ℓ ∞ -norms respectively. Given a matrix A, denote its (i, j)-th entry by P i,j . Let ∥A∥ 2 denote the spectral norm of matrix A. We let λ min (A) and λ max (A) denote the smallest and largest eigenvalue of a square matrix A, respectively. For a symmetric positive definite matrix A, define the inner product ⟨x, y⟩ A := x ⊤ Ay and the associated norm

We refer to ∥x∥ A as the ℓ A -norm of x.

2. Background and problem setting. In this section, we formally introduce Markov reward processes (MRPs) and the (discounted) policy evaluation problem. We also define linear function approximation of the value function, and present the concrete observation models that we study.

2.1. Markov reward process and policy evaluation. An MRP is described by a tuple (S, P, R, γ), where S = [D] denotes the state space, P is the transition kernel, R is the reward function and γ ∈ (0, 1) is the discount factor. At each iteration, the system moves from the current state s ∈ S to some state s ′ ∈ S with probability P(s ′ |s), while the agent realizes the reward R(s, s ′ ). We denote by r(s) := s ′ ∈S P(s ′ |s)R(s, s ′ ) the expected instantaneous reward generated at state s. Let P denote the transition probability matrix having (i, j)-th entry P i,j = P(j|i). The reward R can also be written in matrix form, i.e., R i,j = R(i, j).

The value function specifies the infinite-horizon discounted reward as a function of the initial state:

In the case where the number of states is finite and equal to D, both the expected reward function r and the value function v * are D-dimensional vectors of reals. The value function is given by the solution to the Bellman equation

Throughout this paper, we assume that the Markov chain is aperiodic, ergodic and that there exists a unique stationary distribution π := (π 1 , ..., π D ) with strictly positive entries, satisfying πP = π. Let Π := diag(π 1 , ..., π D ) denote a D × D diagonal matrix whose non-zero elements are given by the entries of the stationary distribution.

In modern applications with large state spaces, it is common to seek approximate solutions to the Bellman equation (2.1), and the standard approach is to choose a d-dimensional subspace S for the purposes of approximation. In particular, one chooses S := span{ψ 1 , ..., ψ d } for d linearly independent basis vectors ψ 1 , . . . , ψ d . For each state s ∈ [D] we let ψ(s) := [ψ 1 (s), ψ 2 (s), ..., ψ d (s)] ⊤ denote its feature vector. Letting Π S denote the projection onto the subspace with respect to the ∥ • ∥ Π -norm, define v as the solution to the projected fixed point equation

It is convenient to write this projection in matrix notation. Let Ψ := [ψ 1 , ψ 2 , ..., ψ d ] ⊤ , and for v ♢ in S, use θ ♢ to denote its corresponding parameterization in R d , e.g., Ψ ⊤ θ ′ = v ′ . With this shorthand, equation (2.2) can be equivalently written as

It is convenient in the analysis to have access to an orthonormal basis spanning the projected space S. Define the matrix B ∈ R d×d by letting B i,j := ⟨ψ i , ψ j ⟩ Π for each i, j, and let

By construction, the vectors ϕ 1 , . . . , ϕ d satisfy ⟨ϕ i , ϕ j ⟩ Π = I(i = j). Next, define the scalars β := λ max (B), and µ := λ min (B), (2.4) so that β/µ is the condition number of the covariance matrix of the features. Finally, let M := γΦΠP Φ ⊤ denote the d-dimensional matrix that describes the action of γP on the projected space S.

This manuscript is for review purposes only.

2.3. Observation models and problem statement. We start by introducing the i.i.d. observation model, in which we have access to a black box or simulator that generates samples from the transition kernel and reward functions. In particular, we observe independent tuples

where ω := (ω 1 , ..., ω D ) is a distribution with strictly positive entries, and we use the shorthand Ω := diag([ω 1 , ..., ω D ]). A natural and popular choice is ω = π, in which case the i.i.d. model is meant to approximate the stationary Markov chain.

In the Markovian noise model, we assume that all of our observations come from a single trajectory of a Markov chain. Precisely, the sequence of states {s 0 , s 1 ...., } generated by the MRP is a time-homogeneous Markov chain, with s 0 ∼ π. The tuple ξ t = (s t , s t+1 , R(s t , s t+1 )) is observed at each time t. The highly correlated nature of these observations renders algorithm design and analysis in the Markovian setting more challenging than in the i.i.d. setting.

Our goal in both cases is to use the observations to generate an estimator v n of v * which satisfies an oracle inequality of the form

where v 0 is the initial iterate of the algorithm. The three terms appearing on the RHS of inequality (2.5) all have concrete interpretations. The first term ∥v -v * ∥ 2 Π characterizes the approximation error incurred by the linear function approximation. As a point of the background, we recall the following instance-dependent upper bound on the approximation error due to Mou et al. [27]:

where A(M, γ) = 1 + λ max (I -M ) -1 (γ 2 I d -M M ⊤ )(I -M ) -⊤ . See Mou et al. [27] for a proof, alongside guarantees of information-theoretic optimality. This work focuses on sharply analyzing the last two terms on the RHS of inequality (2.5), both of which have concrete operational interpretations. The term δ n ∥v 0 -v∥ 2 Π is the deterministic error, which characterizes the convergence of the iterative algorithm in the purely deterministic setting. Specifically, the term δ n , which should tend to zero as the number of iterations (or oracle calls) n goes to infinity, quantifies how fast the discrepancy between the initialization v 0 and the approximate solution v diminishes by running the iterative algorithm.

The third term ϵ n,σ is the stochastic error, which is incurred due to the stochastic observation model. Here we use the notation σ as a placeholder for the "noise level" in the observed samples. One should expect the stochastic error ϵ n,σ to go to zero as n goes to infinity or as σ goes to zero. Several previous works mix the deterministic error with stochastic error in their guarantees (see, e.g., [1,27]). However, the key benefit of separating the deterministic error from the stochastic error is that it allows a clean understanding of situations in which either the observations have low noise or parallel implementation may be available. In these cases, the deterministic error dominates the overall convergence rate of the algorithm, and so having algorithms that attain the optimal deterministic error is a key desideratum.

Having precisely defined the deterministic and stochastic errors, we are now in a position to present our first set of results on lower bounds for both of these terms.

This manuscript is for review purposes only.

3. Lower bounds in weighted ℓ 2 -norm. We study the oracle complexity lower bound on deterministic error in Section 3.1 and the instance-specific stochastic error lower bound in Section 3.2.

3.1. Oracle complexity lower bound on deterministic error. It is well-known that a linear rate can be achieved for the deterministic policy evaluation problem, and the convergence rate is highly dependent on the effective horizon (1-γ) -1 [36]. Accordingly, our goal in this section is to prove an oracle complexity lower bound in terms of (1 -γ) -1 , which can be done even in the tabular setting in which the subspace is all of R D . The following assumption on the oracle captures algorithms in the temporal difference learning family.

Assumption 1 (Amenable iterative method). An amenable iterative method M generates a sequence of iterates v k such that

where

Noting that G(v) is precisely the temporal difference operator applied at the point v, an amenable algorithm is one whose iterates are always in the linear span of the initial point v 0 and subsequent temporal differences. The linear span assumption is commonly used in proving oracle complexity lower bounds [31,32], and as such, nearly all the algorithms in the temporal difference family can be shown to be amenable. The sole exceptions that we are aware of occur in cases where there are projections involved in the algorithm, e.g., [1].

However, in policy evaluation problems with unbounded feasible region R D , projection steps are often unnatural and vanilla TD algorithms are able to attain similar performance (see, e.g., [20]). The following theorem provides an oracle complexity lower bound for policy evaluation problem under the ℓ Π -norm for amenable algorithms.

Theorem 3.1. Fix a constant γ > 1 2 and an initialization v 0 . There exists a transition kernel P and an expected reward vector r such that any iterative method M satisfying Assumption 1 produces iterates {v k } k≥1 satisfying the following.

where v * is the solution of equation (2.1). See Section 7.1 for the proof of this theorem.

Noting that 2γ -1 = 1 -2(1 -γ), Theorem 3.1 shows an oracle complexity lower bound

It should be noted that the metric (i.e., the ℓ Π -norm) used in Theorem 3.1 depends on the problem instance through the stationary distribution of the transition kernel P . Such an instance-dependent metric makes the construction of our worst-case instance non-standard and challenging.

On a related note, it is instructive to recall that classical oracle complexity bounds for solving linear operator equations [29,30] allow the conjugate operator to be queried within the oracle, making the method class wider than the class of amenable algorithms captured in Assumption 1. On the one hand, the conjugate operator is not natural for solving a policy This manuscript is for review purposes only.

evaluation problem under stochastic settings since the vector (I -γP ) ⊤ v is hard to estimate with transition and reward samples. On the other hand, our construction used in proving Theorem 3.1 naturally extends to this wider method class, and we provide an even stronger deterministic error lower bound than Theorem 3.1 in Appendix SM1.

Instance-specific lower bound on stochastic error. We now turn our attention to proving lower bounds on the instance-specific sample complexity under the i.i.d. observation model introduced in Section 2.3. We assume that the feature matrix Ψ is fixed and known, and let ϑ = (ω, P, R) denote an individual problem instance parameterized by the initial state distribution ω, transition kernel P , and reward function R. Note at this juncture that we do not require that ω = π; this is akin to the so-called off-policy situation in which the sampling (or behavior) policy may differ from the policy that we are interested in evaluating. Our result will apply in this general case; but given that the initial state is drawn from the distribution ω, it is convenient to consider solving the projected fixed point equation with respect to the

3)), written as

Use the function θ(ϑ) := (ΨΩΨ ⊤ -ΨΩγP Ψ ⊤ ) -1 ΨΩr to denote the target of interest.

In order to state our result, we require some additional notation. Fix an instance ϑ = (ω, P, R), and for any ϵ > 0, define an ϵ-neighborhood of problem instances by

Adopting the ℓ B -norm as our loss function, define the following local asymptotic minimax risk [11,22]:

The infimum in Eq. (3.4) is taken over all estimators θ N that are measurable functions of N observations drawn according to the i.i.d. observation model. In contrast to the global minimax risk-which takes a supremum of the risk over all the problem instances within a reasonable class-the local minimax risk M(ϑ) looks for the hardest alternative in a small neighborhood of the instance ϑ with diameter c/ √ N . To capture the hardest local alternative (in an asymptotic sense) it suffices to take the diameter of the neighborhood to be of the order 1/ √ N . Invoking Eq. (3.3) yields the equivalent definition

The following proposition characterizes the local asymptotic risk M(ϑ) explicitly.

This manuscript is for review purposes only. Proposition 3.2. Consider the i.i.d. observation model with the initial state drawn from the distribution ω. Let Z ∈ R d be a multivariate Gaussian

Then we have

See Section 7.2 for the proof of this theorem.

A few comments are in order. First, it should be noted that this lower bound is distinct from the asymptotic minimax lower bound shown in Khamaru et al. [15], in which a generative observation model is assumed (where we observe transitions from all D initial states) and there is no function approximation. Consequently, our choice of a problem instance of interest is ϑ = (ω, P, R) rather than (P, R) in [15]. Second, and on a related note, it is important that ω be unknown and included in the set of parameters ϑ; if in contrast ω is known a priori, then the local asymptotic minimax risk differs from the characterization (3.6). Finally, we note that Mou et al. [27] provide non-asymptotic, instance-dependent lower bounds on stochastic error for solving projected fixed-point equations using the Bayesian Cramér-Rao bound. However, these lower bounds do not directly apply here, since the family of hardest local alternatives constructed in [27] may not be valid instances in the policy evaluation setting.

Let us now specialize Proposition 3.2 by taking ω = π, where π is the stationary distribution of the transition kernel P . Denote by ϑ π := (π, P, R) the instance of interest. Let Σiid := cov B -1 2 ⟨ψ(s) -γψ(s ′ ), θ⟩ -R(s, s ′ ) ψ(s) for s ∼ π and s ′ ∼ P(•|s).

Applying Proposition 3.2, the local asymptotic minimax risk (3.5) under this setting is then given by lim

Taking stock, we have proved two lower bounds (3.2) and (3.7) on the deterministic and stochastic errors in ℓ Π -norm under the i.i.d. observation model s ∼ π and s ′ ∼ P(•|s). Given these baselines, it is natural to ask whether there is a practical iterative algorithm in the TD family that can achieve both lower bounds, which is the main focus of Section 4.

Algorithms for policy evaluation in the i.i.d. setting. Taking both lower bounds proved in Section 3 as our baseline, we now turn our attention to the question of algorithm design. In this section, we assume the i.i.d. observation model introduced in Section 2.3 with ω = π. In order to state the results clearly, we require some additional notation. For θ ∈ R d , we define the deterministic operator for solving equation (2.3) as

note that θ is the solution to g(θ) = 0 The corresponding stochastic operator calculated from sample ξ i is defined as

. To characterize the "variance" of the stochastic operator under the i.i.d. observation model, we make the following assumption:

Assumption 2. There exists a constant ς ≥ 0 such that for every θ, θ

In words, instead of bounding the "variance" of the stochastic operator directly as in [20], Assumption 2 guarantees that the variance of the difference between stochastic operators with different variables θ under the same data ξ is upper bounded by the distance between the variables. This assumption is critical for implementing the variance-reduction techniques and capturing the instance-dependent stochastic error at the approximate solution θ. Clearly, the parameter ς 2 is bounded provided the features ψ(s) are bounded, and provides a natural measure of "noise" in the problem. Accordingly, we make use of Assumption 2 throughout Sections 4 and 5.

We are now ready to present our algorithms. We start with a variance-reduced version of the TD algorithm that captures the instance-specific stochastic error lower bound but fails to achieve the oracle complexity lower bound on deterministic error. To remedy this issue, we develop an accelerated variance-reduced TD algorithm that matches both lower bounds proved in Section 3.

A warm-up algorithm: variance-reduced temporal difference learning. Variancereduced temporal difference learning (VRTD) solves the policy evaluation problem using epochs. With a slight ambiguity of notation, we let v t and its corresponding parameterization θ t denote the iterates generated within each epoch, and let v 0 and its corresponding parameterization θ 0 denote the initialization of the algorithm. At the beginning of each epoch k, the algorithm uses N k samples to compute an averaged stochastic operator g and evaluates it at a point θ, where θ should be understood as the best current approximation of the optimal solution. The vector g( θ) is used to recenter the updates in each epoch.

This manuscript is for review purposes only.

end for Output of the epoch:

Note that this algorithm is distinct from previous instantiations of variance-reduced temporal difference algorithms [50,15], since the output of each epoch (4.5) is a weighted average of the iterates. The following theorem provides a convergence guarantee on the VRTD algorithm.

Set the output of the epoch to be

. Then for each δ > 0, we have

which is of the order

A few comments on the upper bound provided in Theorem 4.1 are in order.

This manuscript is for review purposes only.

Comparing the upper and lower bounds. We first focus on the stochastic error in the bound (4.7). The VRTD algorithm requires at least O{ ς 2

(1-γ) 2 µ } samples in each epoch, which accounts for the first term. Note that with noisy observations, it is necessary to have O{ 1

samples in order to obtain an estimate of the value function within O(1) error, so this higherorder term is natural. The dominating stochastic error term is the last term, and matches the lower bound in equation (3.7). Therefore the VRTD algorithm is instance-optimal in terms of its stochastic error.

Next, we turn our attention to the deterministic error, noticing that the dependence on Comparing with related work. To our knowledge, the only work using variance reduction that captures the correct instance-specific stochastic error is that of Khamaru et al. [15], which

showed that the VRTD algorithm can match the lower bound on stochastic error in ℓ ∞ -norm.

However, their guarantees require O{ 1 (1-γ) 3 } samples in each epoch to compute the recentered update, and this sample size is suboptimal. In addition, the deterministic error proved in this paper is of the order

which is also suboptimal. 3 The work of Mou et al. [27] provided an analysis for the Polyak-Ruppert averaged temporal difference learning algorithm with linear function approximation in the weighted ℓ 2 -norm, showing that the dominant stochastic error term matches the stochastic lower bound proved in Proposition 3.2. However, the sample complexity is suboptimal due to the presence of higher-order terms (see [15] for simulations demonstrating this suboptimality), as is the oracle complexity.

Variance-reduced fast temporal temporal difference algorithm. Motivated by the suboptimality of VRTD in its oracle complexity, we now present a variance-reduced "fast" temporal difference (VRFTD) algorithm, which incorporates the idea of operator extrapolation introduced in [19]. This serves to accelerate the algorithm, and our analysis of VRFTD shows a convergence rate matching both the deterministic and stochastic error lower bounds.

The VRFTD algorithm is formally presented in Algorithm 4.2, and we introduce the basic idea of the algorithm below. First, it utilizes the idea of recentering updates from VRTD with the operator g( θ) used in each epoch. Second, in terms of iterate updating in the inner loop, it involves an inner mini-batch that generates the averaged operator g t , which

allows the algorithm to be run with a much larger stepsize. Finally, each iteration within an epoch involves an operator extrapolation step (4.8). This is crucial to achieving the optimal deterministic error (cf. the VRTD update (4.4)).

The following theorem establishes a convergence rate for the VRFTD algorithm.

Theorem 4.2. Fix the total number of epochs K and a positive integer N . Assume that for each epoch, the parameters η, λ, m, N k , T satisfy

Collect m sample tuples ξ t j = (s j , s ′ j , R(s j , s ′ j )) from the i.i.d. observation model. Calculate

end for Output of the epoch:

Set the output of this epoch to be

. Then for δ > 0,

See Section 7.4 for the detailed proof of this theorem.

In view of Theorem 4.2, the number of epochs required by the VRFTD method to find a

The total number of samples used is K k=1 (mT + N k ), which is bounded on the order

Similar to the VRTD algorithm, the VRFTD algorithm achieves optimal sample complexity in terms of stochastic error. For the deterministic error, the dependence on 1/(1 -γ) matches the oracle complexity lower bound proved in Theorem 3.1. Note that the term β/µ is the condition number of the feature matrix in ℓ Π -norm. Therefore, for "well-conditioned" feature matrices, the VRFTD algorithm achieves optimal oracle complexity 4 . In summary, VRFTD is an accelerated and instance-optimal policy evaluation algorithm, and answers the central question posed in this paper.

5. Algorithm for policy evaluation in the Markovian setting. Finally, we extend the VRFTD algorithm to the Markovian setting, noting in passing that such an extension is also possible for the VRTD algorithm. The challenge of Markovian noise stems from the presence of dependent data that leads to biased samples. To control the bias caused by correlation, we need a standard ergodicity assumption on the underlying Markov chain. Assumption 3. There exist constants C P > 0 and ρ ∈ (0, 1) such that

In other words, with the following definition of mixing time

Assumption 3 guarantees that the mixing time is bounded asfoot_5 t mix ≤ log(4C P ) log(1/ρ) .

In order to overcome the difficulty caused by highly-correlated data, we introduce a burn-in period for sample collection. For instance, to compute the operator g defined in Algorithm 4.2, we collect N k successive samples and only use the last N k -n 0 of them. With this method, we are able to reduce the bias induced by Markovian samples and achieve the desired variance reduction properties. The following two lemmas make this quantitative.

Lemma 5.1. For every t, τ ∈ Z + , with probability 1,

where C M :=

See Appendix SM3.1 for a proof of this lemma. In words, Lemma 5.1 provides an upper bound on the bias of the stochastic operator at the solution θ in terms of the approximation error, and the bound decays exponentially with τ . Lemma 5.2. For every t, τ ∈ Z + and θ, θ ′ ∈ R d , with probability 1,

See Appendix SM3.2 for a proof of this lemma. In contrast to Lemma 5.1, Lemma 5.2 provides an upper bound on the bias of the difference of stochastic operators, which allows us to get rid of any dependence on the approximation error. We are now ready to formally state the VRFTD algorithm in the Markovian noise setting in Algorithm 5.1.

for t = 1, . . . , T do Collect m successive samples ξ t j := (s j , s j+1 , R(s j , s j+1 )) from the Markov trajectory.

end for Output of the epoch:

end for Before presenting our main convergence result for the VRFTD algorithm, we first define the matrix ΣMkv , which is a covariance matrix analog for the Markovian case (see Mou et al. [28] and references therein). Letting { s t } ∞ t=-∞ define a sequence of samples obtained from a stationary Markov trajectory, define

where ξ t := ( s t , s t+1 , R( s t , s t+1 )). This matrix is an infinite sum of matrices where one of the summands (when t = 0) is the matrix Σiid defined in Eq. (3.2).

Similarly to the i.i.d. setting, the instance-dependent complexity of Markovian linear stochastic approximation was shown in [28] to be governed by the trace of the matrix (I d -

To interpret this functional, consider the special case in which the approximation error caused by linear function approximation is 0, i.e., v * = v. Let F i denote the σ-field generated by samples ξ 0 , ..., ξ i and let Π i j := diag{[P( s j = 1| s i ), ..., P(

This manuscript is for review purposes only.

for j ≥ i. Then we have

where the last equation follows from the fact that v * = v = Ψ ⊤ θ and the Bellman equation

Armed with this intuition, we are now ready to establish the main convergence result for the VRFTD algorithm under Markovian noise. Given the calculation above, we discuss the cases v = v * and v ̸ = v * separately for clarity.

Theorem 5.3. Fix the total number of epochs K and a positive integer N . Consider an

Suppose the parameters n 0 and m 0 satisfy

, and ρ m 0 ≤ min

where

See Section 7.5 for detailed proofs of Theorem 5.3. Let us now discuss a few aspects of the theorem.

Estimation of mixing time. From the conditions above, e.g., Ineq. (5.7), the parameters τ, n 0 , m 0 scale linearly in the mixing time t mix and logarithmically in other problem parameters.

As such, only some rough estimation of the mixing time is sufficient, which has been the topic of active research. Nontrivial confidence intervals for the reversible case can be found in Hsu et al. [12]. There are also guarantees in the more challenging and prevalent case when the underlying Markov chain is non-reversible [48].

This manuscript is for review purposes only.

Sample complexities. We first consider the case when v = v * . In view of Ineq. (5.9) in Theorem 5.3, the total number of samples required by the VRFTD method to find a solution

, which is bounded on the orderfoot_6

where the mixing time only enters along with terms that scale logarithmically in 1/ϵ. The phenomenon that the mixing time does not enter multiplicatively with the leading-order stochastic error term was also noticed by Li et al. [24] for vanilla TD learning, but as mentioned before, this algorithm does not attain the correct instance-dependent stochastic error.

When v ̸ = v * , the Markovian setting has a biased stochastic operator at optimal solution v, and a larger approximation error ∥v -v * ∥ 2 Π caused by linear function approximation enlarges the bias of the stochastic operator and consequently enlarges the dependence on the approximation error in Eq. (5.10). Therefore, a natural stopping criterion for the Markovian

Π + ϵ for some absolute constant c > 0. From Ineq. (5.10), the total number of required samples K k=1 (mT + N k ) is bounded on the order

Note that, in this bound, the leading stochastic error matches the lower bound proved in [28], which can depend on the mixing time, but is generally smaller than the product of the i.i.d. stochastic error and the mixing time. These results show a delicate difference between how the mixing time of the Markov chain enters the bound depending on whether the function approximation is exact or not. It should be noted that, the stronger stopping criterion, i.e., finding v to satisfy E[∥ v -v∥ 2 Π ] ≤ ϵ, can also be applied in this setting. We can generate a similar sample complexity by enlarging the constants τ and N k by an additive factor of log(∥v * -v∥ 2 Π ). However, given that the approximation error is unavoidable and generally unknown, there is marginal benefit to using this stronger stopping criterion.

Mou et al. [28] established convergence guarantees for TD with averaging in the Markovian noise setting, showing a similar leading order stochastic error term but without accelerating the deterministic error. Besides improving on the deterministic error, Theorem 5.3 also guarantees that the higher-order terms on stochastic error are smaller than those proved in [28].

6. Numerical experiments. In this section, we report numerical experiments for both VRTD and VRFTD, comparing them against temporal difference learning (TD), conditional temporal difference learning (CTD), and fast temporal difference learning (FTD) [19,20]. To generate a comprehensive performance profile, we conduct experiments under both the i.i.d.

and Markovian noise models.

6.1. The i.i.d. setting: A simple two-state construction. We consider a family of twostate MRPs inspired by the construction of Duan et al. [9]. For a discount factor γ ∈ ( 1 2 , 1), the transition kernel P and reward vector r are given by P = 2γ-1 γ

Clearly the transition kernel is symmetric, thus the stationary distribution is π = [0.5, 0.5].

For simplicity, we choose the feature matrix

), which forms an orthonormal basis under ℓ Π -norm. Assuming the i.i.d. model in which s i ∼ π and s ′ i ∼ P (•|s i ), it can be shown via simple calculation that the stochastic error term is given by trace

(b) FTD algorithm Instance-optimality. We generate a range of MRPs with different values of discount factor γ and run the four aforementioned algorithms on each MRP. In order to test the robustness of our results, we simulate various step-sizes for each algorithm. To be fair in our comparison, we also include a simulation of the best-tuned stepsize for each algorithm. We plot the prediction from the lower bound (6.1) as well.

From subplots (a) and (b) of Figure 1, it is clear that the vanilla TD and FTD algorithms with diminishing stepsizes [20] do not achieve the lower bound calculated in equation (6.1).

This manuscript is for review purposes only.

On the other hand, sub-plots (c) and (d) show that the VRTD and VRFTD algorithms achieve the lower bound (6.1), and that these behaviors are robust to the choice of stepsize parameters. However, given their epoch-wise nature, the outputs of variance-reduced algorithms are more volatile than TD and FTD. Another interesting observation is that the accelerated algorithms-FTD and VRFTD-are less sensitive to stepsize parameters. Our next set of experiments explores this further.

Ablation analysis of VRFTD. Notice that VRFTD includes two new ingredients when compared with VRTD: mini-batching and operator extrapolation (OE). We now perform an ablation analysis to disentangle the contribution of both ingredients. We generate a range of MRPs with different values of γ and run the experimental and control groups on each MRP.

(a) Ablation study for OE (I)

(b) Ablation study for OE (II)

(c) Ablation study for minibatching In the first experiments, we ran the experimental group with OE steps and the control group without OE steps. We first ran both groups with the stepsize policies suggested by the theoretical analysis (see subplot (a) of Figure 2). The results indicate that the experimental group significantly outperforms the control group. However, this performance difference can largely be attributed to the conservative stepsize for the control group, as prescribed by the theory. To be more fair to both algorithms, we further fine-tuned the stepsize parameters of both algorithms and obtained subplot (b), where the two algorithms exhibit similar convergence rates. Taking stock, the first set of experiments shows that the analysis and stepsize policy of the VRFTD (with OE steps) serves as a better theoretical guideline for practical applications.

To demonstrate the benefits of mini-batching in the inner loop, we ran a second experiment with two groups with and without mini-batching. Note that we keep the stepsizes and the total number of samples the same for both groups (which means that the control group without mini-batching has larger epoch lengths). From subplot (c), we can see that the performance difference is significant, showing that without mini-batching, the algorithm exhibits instability when run with aggressive stepsize policies.

This manuscript is for review purposes only.

6.2. The Markovian setting: 2D Grid World. Our experiments under the Markovian noise model are conducted on the 2D Grid World environment. This is a classical problem in reinforcement learning with finite state and action spaces. An agent realizes a positive reward when reaching a predetermined goal and negative ones when going through "traps". The dimension of the state space is set to be D = 400, among which we assign a goal state (with reward r = 1) and 30 traps (with reward r = -0.2). The transition kernel is fixed as follows:

With probability 0.95, the agent moves in a direction that points towards the goal and with probability 0.05 in a random direction. Our goal is to compute the value function v * -for each possible initial state of the agent. We also incorporate linear function approximation in these experiments. Specifically, we generate random features with dimension d = 50 to estimate the D = 400 dimensional value function in this problem. We test the performance of four algorithms, with the discount factor γ set to 0.99 and 0.999. Figure 3 plots the normalized error in ℓ Π -norm against the length of the trajectory. In both experiments, the VRFTD algorithm exhibits the fastest convergence to the true value function, thereby corroborating our theoretical results. Note that in both experiments, the estimation errors do not converge to zero, because there is a nontrivial error incurred by linear function approximation. Another salient takeaway is the following: Closer to the optimal solution v * , the variance-reduced algorithms (VRTD/VRFTD) achieve faster convergence rate compared to their counterparts that do not employ variance reduction.

7. Proofs. In this section, we provide the proofs of our main results. The proof of other auxiliary results are collected in the supplementary material. The matrix I -γP is square and full rank, and it is straightforward to verify that the unique solution of the linear equation (I -γP

It is also easy to see that the Markov chain induced by the transition kernel P is irreducible and aperiodic. From the cyclical behavior of the Markov chain, we have that the stationary distribution is uniform, i.e.,

Thus, we obtain

vectors lying in the span of the first k standard basis vectors. Since all entries of I -γP below its subdiagonal are equal to 0 and all entries of r except for its first are equal to 0, we conclude that v k ∈ R k,D . Therefore,

as desired.

where the infimum is taken over all estimators h N that are measurable functions of N i.i.d. data points drawn from P ϑ ′ , and the expectation is taken over a multivariate Gaussian Z ∼ N 0, ∇h(ϑ) ⊤ J † ϑ ∇h(ϑ) .

In our model, we have ϑ = (ω, P, R) and h(ϑ) = (ΨΩΨ ⊤ -ΨΩγP Ψ ⊤ ) -1 ΨΩr. We set the loss function L to be ∥ • ∥ 2

. Invoking Theorem 7.1 yields

The covariance is explicitly computed in the following lemma. Lemma 7.2. We have

The proof of this lemma is straightforward but involves some lengthy calculations, which we defer to Appendix SM3.3.

Recall the definition of M and Σ in the statement of Proposition 3.2. By substituting equality (7.7) into (7.6) and invoking the relation

which completes the proof of Proposition 3.2. We have

Set the output of this epoch to be

. Then we have

This manuscript is for review purposes only.

See Appendix SM2.2 for the proof of this proposition. The basic idea of the proof is first providing an upper bound on the term E∥ v k -v∥ 2 Π and then combining it with Lemma 7.3.

Taking Proposition 7.4 as given for the moment, let us complete the proof of the theorem.

The main idea is to bound the approximation error term ∥v -v * ∥ 2 Π separately from the term

To bound ∥v -v * ∥ 2 Π , we use the instance-dependent upper bound in Ineq. (2.6).

We bound the term ∥ v K -v∥ 2 Π by using Ineq. (7.10) as follows

To conclude, we use Young's inequality to obtain

which completes the proof. See Appendix SM2.4 for the proof of this proposition.

Taking Proposition 7.6 as given for the moment, let us complete the proof of the theorem.

First, consider the case when v ̸ = v * . Recursively using Ineq. (7.17) yields

where step (i) follows from the condition N k -n 0 ≥ ( 3 4 ) K-k N . The proof of the case when v = v * follows from the same derivation.

8. Discussion. In this paper, we investigated the problem of policy evaluation with linear function approximation, making three contributions. First, we proved lower bounds on both deterministic error and stochastic error. With these lower bounds in hand, we presented an analysis of a variance-reduced variant of temporal difference algorithm (VRTD) in the i.i.d.

observation model and showed that it fails to match the oracle complexity lower bound on the deterministic error. In order to remedy this difficulty, we developed an optimal variancereduced fast temporal difference algorithm (VRFTD) that nearly matches both lower bounds simultaneously. Finally, we extended the VRFTD algorithm to the Markovian setting and provided instance-dependent convergence results. The leading stochastic error matches the instance-dependent lower bound for Markovian linear stochastic approximation [28], and the deterministic error matches the i.i.d. setting up to a multiplicative factor proportional to the mixing time of the chain. Our theoretical guarantees were corroborated with numerical experiments in both the i.i.d. and Markovian settings, showing that the VRFTD algorithm enjoys several advantages over the prior state-of-the-art.

Our work leaves open severaal salient future directions; let us mention two. First, our oracle complexity lower bound is proved in the tabular setting. On the other hand, our upper bounds on the deterministic error indicate that with linear function approximation, we pay

This manuscript is for review purposes only.

a multiplicative factor depending on the condition number of the feature matrix. It would be interesting to see if an oracle complexity lower bound can be proved under linear function approximation, and whether the linear dependence on the condition number in our bounds is optimal. Second, and more broadly, note that our analysis relies heavily on the linear structure of the problem. However, there are many problems in the reinforcement learning literature that have nonlinear structures, e.g., the policy optimization problem involving the Bellman optimality operator. An interesting direction for future work is to understand problems with nonlinear structure from an instance-specific point of view and develop efficient algorithms to capture the optimal deterministic and stochastic errors. For instance, variance reduction has been applied to the policy optimization problem under the generative model [38,46] and some instance-dependent bounds are known (e.g., for variants of Q-learning [16]). It is an important open question to develop acceleration schemes for such algorithms in a fashion similar to our paper, while extending the results to the more realistic Markovian setting.