REANALYSIS OF VARIANCE REDUCED TEMPORAL DIF-
Year: 2020

Abstract:
Temporal difference (TD) learning is a popular algorithm for policy evaluation in reinforcement learning, but the vanilla TD can substantially suffer from the inherent optimization variance. A variance reduced TD (VRTD) algorithm was proposed by Korda and La (2015), which applies the variance reduction technique directly to the online TD learning with Markovian samples. In this work, we first point out the technical errors in the analysis of VRTD in Korda and La (2015), and then provide a mathematically solid analysis of the non-asymptotic convergence of VRTD and its variance reduction performance. We show that VRTD is guaranteed to converge to a neighborhood of the fixed-point solution of TD at a linear convergence rate. Furthermore, the variance error (for both i.i.d. and Markovian sampling) and the bias error (for Markovian sampling) of VRTD are significantly reduced by the batch size of variance reduction in comparison to those of vanilla TD. As a result, the overall computational complexity of VRTD to attain a given accurate solution outperforms that of TD under Markov sampling and outperforms that of TD under i.i.d. sampling for a sufficiently small conditional number.

heading: 1 INTRODUCTION
text: In reinforcement learning (RL), policy evaluation aims to obtain the expected long-term reward of a given policy and plays an important role in identifying the optimal policy that achieves the maximal cumulative reward over time Bertsekas and Tsitsiklis (1995); Dayan and Watkins (1992); Rummery and Niranjan (1994). The temporal difference (TD) learning algorithm, originally proposed by Sutton (1988), is one of the most widely used policy evaluation methods, which uses the Bellman equation to iteratively bootstrap the estimation process and continually update the value function in an incremental way. In practice, if the state space is large or infinite, function approximation is often used to find an approximate value function efficiently. Theoretically, TD with linear function approximation has been shown to converge to the fixed point solution with i.i.d. samples and Markovian samples in Sutton (1988); Tsitsiklis and Van Roy (1997). The finite sample analysis of TD has also been studied in Bhandari et al. (2018); Dalal et al. (2018a); Cai et al. (2019); Srikant and Ying (2019).
Since each iteration of TD uses one or a mini-batch of samples to estimate the mean of the pseudogradient 1, TD learning usually suffers from the inherent variance, which substantially degrades the convergence accuracy. Although a diminishing stepsize or very small constant stepsize can reduce the variance Bhandari et al. (2018); Srikant and Ying (2019), they also slow down the convergence significantly.
Two approaches have been proposed to reduce the variance. The first approach is the so-called batch TD, which takes a fixed sample set and transforms the empirical mean square projected Bellman error (MSPBE) into an equivalent convex-concave saddle-point problem Du et al. (2017). Due to the finite-sample nature of such a problem, stochastic variance reduction techniques for conventional optimization can be directly applied here to reduce the variance. In particular, Du et al. (2017) showed that SVRG Johnson and Zhang (2013) and SAGA Defazio et al. (2014) can be applied to improve
1We call the increment in each iteration of TD as "pseudo-gradient" because such a quantity is not a gradient of any objective function, but the role that it serves in the TD algorithm is analogous to the role of the gradient in the gradient descent algorithm.
the performance of batch TD algorithms, and Peng et al. (2019) proposed two variants of SVRG to further save the computation cost. However, the analysis of batch TD does not take into account the statistical nature of the training samples, which are generated by a MDP. Hence, there is no guarantee of such obtained solutions to be close to the fixed point of TD learning.
The second approach is the so-called TD with centering (CTD) algorithm proposed in Korda and La (2015), which introduces the variance reduction idea to the original TD learning algorithm. For the sake of better reflecting its major feature, we refer to CTD as Variance Reduced TD (VRTD) throughout this paper. Similarly to the SVRG in Johnson and Zhang (2013), VRTD has outer and inner loops. The beginning of each inner-loop (i.e. each epoch) computes a batch of sample pseudogradients so that each subsequent inner loop iteration modifies only one sample pseudo-gradient in the batch pseudo-gradients to reduce the variance. The main difference between VRTD and batch TD is that VRTD applies the variance reduction directly to TD learning rather than to a transformed optimization problem in batch TD. Though Korda and La (2015) empirically verified that VRTD has better convergence accuracy than vanilla TD learning, some technical errors in the analysis in Korda and La (2015) have been pointed out in follow up studies Dalal et al. (2018a); Narayanan and Szepesvári (2017). Furthermore, as we discuss in Section 3, the technical proof in Korda and La (2015) regarding the convergence of VRTD also has technical errors so that their results do not correctly characterize the impact of variance reduction on TD learning. Given the recent surge of interest in the finite time analysis of the vanilla TD Dalal et al. (2018a); Bhandari et al. (2018); Dalal et al. (2018b); Srikant and Ying (2019), it becomes imperative to reanalyze the VRTD and accurately understand whether and how variance reduction can help to improve the convergence accuracy over vanilla TD. Towards this end, this paper specifically addresses the following central questions. • For i.i.d. sampling, it has been shown in Dalal et al. (2018a); Bhandari et al. (2018) that vanilla
TD converges only to a neighborhood of the fixed point for a constant stepsize and suffers from a constant error term caused by the variance of the stochastic pseudo-gradient at each iteration. For VRTD, does the variance reduction help to reduce such an error and improve the accuracy of convergence? How does the error depend on the variance reduction parameter, i.e., the batch size for variance reduction? • For Markovian sampling, it has been shown in Bhandari et al. (2018) that the convergence of vanilla TD further suffers from a bias error due to the correlation among samples in addition to the variance error as in i.i.d. sampling. Does VRTD, which was designed to have reduced variance, also enjoy reduced bias error? If so, how does the bias error depend on the batch size for variance reduction? • For both i.i.d. and Markovian sampling, to attain an -accurate solution, what is the overall computational complexity (the total number of computations of pseudo-gradients) of VRTD, and does VRTD have a reduced overall computational complexity compared to TD?

heading: 1.1 OUR CONTRIBUTIONS
text: Our main contributions are summarized in Table 1 and are described as follows.
For i.i.d. sampling, we show that a slightly modified version of VRTD (for avoiding bias error) converges linearly to a neighborhood of the fixed point solution for a constant stepsize α, with the variance error at the order ofO(α/M), where M is the batch size for variance reduction. To attain an -accurate solution, the overall computational complexity (i.e., the total number of pseudo-gradient computations) of VRTD outperforms the vanilla TD algorithm Bhandari et al. (2018) if the condition number is small.
For Markovian sampling, we show that VRTD has the same linear convergence and the same variance error reduction over the vanilla TD Dalal et al. (2018a); Bhandari et al. (2018) as i.i.d. sampling. More importantly, the variance reduction in VRTD also attains a substantially reduced bias error at the order of O(1/M) over the vanilla TD Bhandari et al. (2018), where the bias error is at the order of O(α). As a result, VRTD outperforms vanilla TD in terms of the total computational complexity by a factor of log 1 .
At the technical level, our analysis of bias error for Markovian sampling takes a different path from how the existing analysis of TD handles Markovian samples in Bhandari et al. (2018); Wang et al. (2017); Srikant and Ying (2019). Due to the batch average of stochastic pseudo-gradients adopted by VRTD to reduce the variance, the correlation among samples in different epochs is eliminated. Such an analysis explicitly explains why the variance reduction step helps to further reduce the bias error.

heading: 1.2 RELATED WORK
text: On-policy TD learning and variance reduction. On-policy TD learning aims to minimize the Mean Squared Bellman Error (MSBE) Sutton (1988) when samples are drawn independently from the stationary distribution of the corresponding MDP. The non-asymptotic convergence under i.i.d. sampling has been established in Dalal et al. (2018a) for TD with linear function approximation and for TD with overparameterized neural network approximation Cai et al. (2019). The convergence of averaged linear SA with constant stepsize has been studied in Lakshminarayanan and Szepesvari (2018). In the Markovian setting, the non-asymptotic convergence has been studied for on-policy TD in Bhandari et al. (2018); Karmakar and Bhatnagar (2016); Wang et al. (2019); Srikant and Ying (2019). Korda and La (2015) proposed a variance reduced CTD algorithm (called VRTD in this paper), which directly applies variance reduction technique to the TD algorithm. The analysis of VRTD provided in Korda and La (2015) has technical errors. The aim of this paper is to provide a technically solid analysis for VRTD to characterize the advantage of variance reduction.
Variance reduced batch TD learning. Batch TD Lange et al. (2012) algorithms are generally designed for policy evaluation by solving an optimization problem on a fixed dataset. In Du et al. (2017), the empirical MSPBE is first transformed into a quadratic convex-concave saddle-point optimization problem and variance reduction methods of SVRG Johnson and Zhang (2013) and SAGA Defazio et al. (2014) were then incorporated into a primal-dual batch gradient method. Furthermore, Peng et al. (2019) applied two variants of variance reduction methods to solve the same saddle point problems, and showed that those two methods can save pseudo-gradient computation cost.
We note that due to the extensive research in TD learning, we include here only studies that are highly related to our work, and cannot cover many other interesting topics on TD learning such as asymptotic convergence of TD learning Tadić (2001); Hu and Syed (2019), off-policy TD learning Sutton et al. (2008; 2009); Liu et al. (2015); Wang et al. (2017); Karmakar and Bhatnagar (2017), two time-scale TD algorithms Xu et al. (2019); Dalal et al. (2018b); Yu (2017), fitted TD algorithms Lee and He (2019), SARSA Zou et al. (2019) etc. The idea of the variance reduction algorithm proposed in Korda and La (2015) as well as the analysis techniques that we develop in this paper can potentially be useful for these algorithms.

heading: 2 PROBLEM FORMULATION AND PRELIMINARIES

heading: 2.1 ON-POLICY VALUE FUNCTION EVALUATION
text: We describe the problem of value function evaluation over a Markov decision process (MDP) (S,A,P, r, γ), where each component is explained in the sequel. Suppose S ⊂ Rd is a compact state space, and A is a finite action set. Consider a stationary policy π, which maps a state s ∈ S to the actions in A via a probability distribution π(·|s). At time-step t, suppose the process is in some state st ∈ S, and an action at ∈ A is taken based on the policy π(·|st). Then the transition kernel P = P(st+1|st, at) determines the probability of being at state st+1 ∈ S in the next time-step, and the reward rt = r(st, at, st+1) is received, which is assumed to be bounded by rmax. We denote the associated Markov chain by p(s′|s) = ∑ a∈A p(s
′|s, a)π(a|s), and assume that it is ergodic. Let µπ be the induced stationary distribution, i.e., ∑ s p(s
′|s)µπ(s) = µπ(s′). We define the value function for a policy π as vπ (s) = E[ ∑∞ t=0 γ tr(st, at, st+1)|s0 = s, π], where γ ∈ (0, 1) is the discount
factor. Define the Bellman operator Tπ for any function ξ(s) as Tπξ(s) := rπ(s) + γEs′|sξ(s′), where rπ(s) = Ea,s′|sr(s, a, s′) is the expected reward of the Markov chain induced by the policy π. It is known that vπ(s) is the unique fixed point of the Bellman operator Tπ, i.e., vπ(s) = Tπvπ(s). In practice, since the MDP is unknown, the value function vπ(s) cannot be directly obtained. The goal of policy evaluation is to find the value function vπ(s) via sampling the MDP.

heading: 2.2 TD LEARNING WITH LINEAR FUNCTION APPROXIMATION
text: In order to find the value function efficiently particularly for large or infinite state space S, we take the standard linear function approximation v̂(s, θ) = φ(s)>θ of the value function, where φ(s)> = [φ1(s), · · · , φd(s)] with φi(s) for i = 1, 2, · · · d denoting the fixed basis feature functions of state s, and θ ∈ Rd is a parameter vector. Let Φ be the |S| × d feature matrix (with rows indexed by the state and columns corresponding to components of θ). The linear function approximation can be written in the vector form as v̂(θ) = Φθ. Our goal is to find the fixed-point parameter θ∗ ∈ Rd that satisfies Eµπ v̂(s, θ∗) = EµπTπ v̂(s, θ∗). The TD learning algorithm performs the following fixed-point iterative update to find such θ∗.
θt+1 = θt + αtgxt(θt) = θt + αt(Axtθt + bxt), (1)
where αt > 0 is the stepsize, and Axt and bxt are specified below. For i.i.d. samples generated from the distribution µπ, we denote the sample as xt = (st, rt, s′t), and Axt = φ(st)(γφ(s ′ t) − φ(st)) > and bxt = r(st)φ(st). For Markovian samples generated sequentially from a trajectory, we denote the sample as xt = (st, rt, st+1), and in this case Axt = φ(st)(γφ(st+1) − φ(st))> and bxt = r(st)φ(st). We further define the mean pseudo-gradient g(θ) = Aθ + b where A = Eµπ [φ(s)(γφ(s′)−φ(s))>] and b = Eµπ [r(s)φ(s)]. We call g(θ) as pseudo-gradient for convenience due to its analogous role as in the gradient descent algorithm. It has been shown that the iteration in eq. (1) converges to the fix point θ∗ = −A−1b at a sublinear rate O(1/t) with diminishing stepsize αt = O(1/t) using both Markovian and i.i.d. samples Bhandari et al. (2018); Dalal et al. (2018a). Throughout the paper, we make the following standard assumptions Wang et al. (2017); Korda and La (2015); Tsitsiklis and Van Roy (1997); Bhandari et al. (2018). Assumption 1 (Problem solvability). The matrix A is non-singular. Assumption 2 (Bounded feature). ‖φ(s)‖2 ≤ 1 for all s ∈ S. Assumption 3 (Geometric ergodicity). The considered MDP is irreducible and aperiodic, and there exist constants κ > 0 and ρ ∈ (0, 1) such that
sup s∈S
dTV (P(st ∈ ·|s0 = s), µπ(s)) ≤ κρt, ∀t ≥ 0,
where dTV (P,Q) denotes the total-variation distance between the probability measures P and Q.
Assumption 1 requires the matrix A to be non-singular so that the optimal parameter θ∗ = −A−1b is well defined. Assumption 2 can be ensured by normalizing the basis functions {φi}di=1. Assumption 3 holds for any time-homogeneous Markov chain with finite state-space and any uniformly ergodic Markov chains with general state space.

heading: 3 THE VARIANCE REDUCED TD ALGORITHM
text: In this section, we first introduce the variance-reduced TD (VRTD) algorithm proposed in Korda and La (2015) for Markovian sampling and then discuss the technical errors in the analysis of VRTD in Korda and La (2015).

heading: 3.1 VRTD ALGORITHM KORDA AND LA (2015)
text: Since the standard TD learning takes only one sample in each update as can be seen in eq. (1), it typically suffers from a large variance. This motivates the development of the VRTD algorithm in Korda and La (2015) (named as CTD in Korda and La (2015)). VRTD is formally presented in Algorithm 2, and we briefly introduce the idea below. The algorithm runs in a nested fashion with each inner-loop (i.e., each epoch) consists of M updates. At the beginning of the m-th epoch, a batch of M samples are acquired and a batch pseudo-gradient gm(θ̃m−1) is computed based on these samples as an estimator of the mean pseudo-gradient. Then, each inner-loop update randomly takes
Algorithm 1 Variance Reduced TD with iid samples
Input: batch size M , learning rate α and initialization θ̃0
1: for m = 1, 2, ..., S do 2: θm,0 = θ̃m−1 3: Sample a set Bm with M samples indepedently
from the distribution µπ 4: gm(θ̃m−1) = 1M ∑ xi∈Bm gxi(θ̃m−1) 5: for t = 0, 1, ...,M − 1 do 6: Sample xjm,t indepedently from the distribu-
tion µπ 7: θm,t+1 = θm,t + α ( gxjm,t (θm,t)
8: −gxjm,t (θ̃m−1) + gm(θ̃m−1) )
9: end for 10: set θ̃m = θm,t for randomly chosen t ∈ {1, 2, ...,M} 11: end for Output: θ̃S
Algorithm 2 Variance Reduced TD with Markovian samples Korda and La (2015)
Input: batch size M , learning rate α and initialization θ̃0
1: for m = 1, 2, ..., S do 2: θm,0 = θ̃m−1 3: gm(θ̃m−1) = 1M ∑mM−1 i=(m−1)M gxi(θ̃m−1) 4: for t = 0, 1, ...,M − 1 do 5: Sample jm,t uniformly at random in {(m −
1)M, ...,mM − 1} from trajetory 6: θm,t+1 = ΠRθ ( θm,t + α ( gxjm,t (θm,t)
7: −gxjm,t (θ̃m−1) + gm(θ̃m−1) )) 8: end for 9: set θ̃m = θm,t for randomly chosen t ∈ {1, 2, ...,M}
10: end for Output: θ̃S
one sample from the batch, and updates the corresponding component in gm(θ̃m−1). Here, ΠRθ in Algorithm 2 denotes the projection operator onto a norm ball with the radius Rθ. The idea is similar to the SVRG algorithm proposed in Johnson and Zhang (2013) for conventional optimization. Since a batch pseudo-gradient is used at each inner-loop update, the variance of the pseudo-gradient is expected to be reduced.

heading: 3.2 TECHNICAL ERRORS IN KORDA AND LA (2015)
text: In this subsection, we point out the technical errors in the analysis of VRTD in Korda and La (2015), which thus fails to provide the correct variance reduction performance for VRTD.
At the high level, the batch pseudo-gradient gm(θ̃m−1) computed at the beginning of each epoch m should necessarily introduce a non-vanishing variance error for a fixed stepsize, because it cannot exactly equal the mean (i.e. population) pseudo-gradient g(θ̃m−1). Furthermore, due to the correlation among samples, the pseudo-gradient estimator in expectation (with regard to the randomness of the sample trajectory) does not equal to the mean pseudo-gradient, which should further cause a non-vanishing bias error in the convergence bound. Unfortunately, the convergence bound in Korda and La (2015) indicates an exact convergence to the fixed point, which contradicts the aforementioned general understanding. More specifically, if the batch size M = 1 (with properly chosen λA defined as λA := 2|λmax(A + A>)|), VRTD reduces to the vanilla TD. However, the exact convergence result in Theorem 3 in Korda and La (2015) does not agree with that of vanilla TD characterized in the recent studies Bhandari et al. (2018), which has variance and bias errors.
In Appendix B, we further provide a counter-example to show that one major technical step for characterizing the convergence bound in Korda and La (2015) does not hold. The goal of this paper is to provide a rigorous analysis of VRTD to characterize its variance reduction performance.

heading: 4 MAIN RESULTS
text: As aforementioned, the convergence of VRTD consists of two types of errors: the variance error due to inexact estimation of the mean pseudo-gradient and the bias error due to Markovian sampling. In this section, we first focus on the first type of error and study the convergence of VRTD under i.i.d. sampling. We then study the Markovian case to further analyze the bias. In both cases, we compare the performance of VRTD to that of the vanilla TD described in eq. (1) to demonstrate its advantage.

heading: 4.1 CONVERGENCE ANALYSIS OF VRTD WITH I.I.D. SAMPLES
text: For i.i.d. samples, it is expected that the bias error due to the time correlation among samples does not exist. However, if we directly apply VRTD (Algorithm 2) originally designed for Markovian samples, there would be a bias term due to the correlation between the batch pseudo-gradient estimate and every inner-loop updates. Thus, we slightly modify Algorithm 2 to Algorithm 1 to avoid the bias error in the convergence analysis with i.i.d. samples. Namely, at each inner-loop iteration, we draw a new sample from the stationary distribution µπ for the update rather than randomly selecting one from the batch of samples drawn at the beginning of the epoch as in Algorithm 2. In this way, the new independent samples avoid the correlation with the batch pseudo-gradient evaluated at the beginning of the epoch. Hence, Algorithm 1 does not suffer from an extra bias error.
To understand the convergence of Algorithm 1 at the high level, we first note that the sample batch pseudo-gradient cannot estimate the mean pseudo-gradient g(θ̃m−1) exactly due to its population nature. Then, we define em(θ̃m) = gm(θ̃m−1)− g(θ̃m−1) as such a pseudo-gradient estimation error, our analysis (see Appendix D) shows that after each epoch update, we have
E [∥∥∥θ̃m − θ∗∥∥∥2
2 ∣∣∣Fm,0] ≤ 1/M + 4α 2(1 + γ)2
αλA − 4α2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 +
2α λA − 4α(1 + γ)2 E [∥∥∥em(θ̃m−1)∥∥∥2 2 ∣∣∣Fm,0] , (2) where Fm,0 denotes the σ-field that includes all the randomness in sampling and updates before the m-th epoch. The first term in the right-hand side of eq. (2) captures the contraction property of Algorithm 1 and the second term corresponds to the variance of the pseudo-gradient estimation error. It can be seen that due to such an error term, Algorithm 1 is expected to have guaranteed convergence only to a neighborhood of θ∗, when applying eq. (2) iteratively. Our further analysis shows that such an error term can still be well controlled (to be small) by choosing an appropriate value for the batch size M , which captures the advantage of the variance reduction. The following theorem precisely characterizes the non-asymptotic convergence of Algorithm 1. Theorem 1. Consider the VRTD algorithm in Algorithm 1. Suppose Assumptions 1–3 hold. Set a constant stepsize α < λA8(1+γ)2 and the batch size M > 4(1+γ)2α2+1 α[λA−8α(1+γ)2] . Then, for all m ∈ N,
E [∥∥∥θ̃m − θ∗∥∥∥2
2
] ≤ Cm1 ∥∥∥θ̃0 − θ∗∥∥∥2 2 + 2D2α (1− C1)(λA − 4α(1 + γ)2)M , (3)
where C1 := ( 4α(1 + γ)2 + 4(1+γ) 2α2+1
αM
) 1
λA−4α(1+γ)2 (with C1 < 1 due to the choices of α and
M ), and D2 = 4((1 + γ)2R2θ + r 2 max).
We note that the convergence rate in eq. (3) can be written in a simpler form as E[||θ̃m − θ∗||2] ≤ Cm1 ||θ̃0 − θ∗||2 +O(α/M). Theorem 1 shows that Algorithm 1 converges linearly (under a properly chosen constant stepsize) to a neighborhood of the fixed point solution, and the size of the neighborhood (i.e., the error term) has the order of O( αM ), which can be made as small as possible by properly increasing the batch size M . This is in contrast to the convergence result of the vanilla TD, which suffers from the constant error term with order O(α) Bhandari et al. (2018) for a fixed stepsize. Thus, a small stepsize α is required in vanilla TD to reduce the variance error, which, however, slows down the practical convergence significantly. In contrast, this is not a problem for VRTD, which can attain a high accuracy solution while still maintaining fast convergence at a desirable stepsize.
We further note that if we have access to the mean pseudo-gradient g(θ̃m−1) in each epoch m, then the error term ||θm − θ∗||22 becomes zero, and Algorithm 1 converges linearly to the exact fixed point solution, as the iteration number m goes to infinity with respect to the conditional number C1, which is a positive constant and less than 1. This is similar to the conventional convergence of SVRG for strongly convex optimization Johnson and Zhang (2013). However, the proof here is very different. In Johnson and Zhang (2013), the convergence proof relies on the relationship between the gradient and the value of the objective function, but there is not such an objective function in the TD learning problem. Thus, the convergence of the parameter θ needs to be developed by exploiting the structure of the Bellman operator.
Based on the convergence rate of VRTD characterized in Theorem 1 under i.i.d. sampling, we obtain the following bound on the overall computational complexity.
Corollary 1. Suppose Assumptions 1-3 hold. Let α = λA16(1+γ)2 , M = dmax{ D2 3(1−C1) 1 ,
33(1+γ)2
λ2A }e.
Then, for any > 0, an -accuracy solution (i.e., E||θ̃m − θ∗||2 ≤ ) can be attained with at most m = dlog 2‖θ̃0−θ∗‖2 2
/ log 1 C1 e iterations. Correspondingly, the total number of pseudo-gradient
computations required by VRTD (i.e., Algorithm 1) under i.i.d. sampling to attain such an -accuracy solution is at most
O ( max {1 , 1
λ2A
} log (1 )) .
Proof. Given the values of α and M in the theorem, it can be easily checked that E||θ̃m − θ∗||2 ≤ for m = dlog 2‖θ̃0−θ∗‖2 2
/ log 1 C1 e. Then the total number of pseudo-gradient computations is given
by 2mM that yields the desired order given in the theorem.
As a comparison, consider the vanilla TD algorithm studied in Bhandari et al. (2018) with the constant stepsize α = O( ). If the samples are i.i.d. generated, it can be shown (see Appendix F.1) that the vanilla TD requires O ( 1 λ2A log ( 1 )) pseudo-gradient computations in total to obtain an -accuracy solution. Clearly, VRTD has lower computational complexity than vanilla TD if λA is small. Such a comparison is similar in nature to the comparison between SVRG Johnson and Zhang (2013) and SGD in traditional optimization, where SVRG achieves better computational complexity than SGD for strongly convex objectives if the conditional number of the loss is small.

heading: 4.2 CONVERGENCE ANALYSIS OF VRTD WITH MAKOVIAN SAMPLES
text: In this section, we study the VRTD algorithm (i.e., Algorithm 2) with Markovian samples, in which samples are generated from one single MDP path. In such a case, we expect that the convergence of VRTD to have both the variance error due to the pseudo-gradient estimation (similar to the case with i.i.d. samples) and the bias error due to the correlation among samples. To understand this at the high level, we define the bias at each iteration as ξm(θ) = (θ − θ∗)>(gm(θ)− g(θ)). Then our analysis (see Appendix E) shows that after the update of each epoch, we have
E [∥∥∥θ̃m − θ∗∥∥∥2
2 ∣∣∣Fm,0] ≤ 1/M + 3α 2(1 + γ)2
αλA − 3α2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 +
3α λA − 3α(1 + γ)2 E [ ‖gm(θ∗)‖22 ∣∣∣Fm,0] + 2
[λA − 3α(1 + γ)2]M M−1∑ i=0 E [ ξm(θm,i) ∣∣∣Fm,0] (4) The first term on the right-hand side of eq. (4) captures the epochwise contraction property of Algorithm 2. The second term is due to the variance of the pseudo-gradient estimation, which captures how well the batch pseudo-gradient gm(θ∗) approximates the mean pseudo-gradient g(θ∗) (note that g(θ∗) = 0). Such a variance term can be shown to decay to zero as the batch size gets large similarly to the i.i.d. case. The third term captures the bias introduced by the correlation among samples in the m-th epoch. To quantitatively understand this error term, we provide the following lemma that characterizes how the bias error is controlled by the batch size M .
Lemma 1. For any m > 0 and any θ ∈ Bθ, which is a ball with the radius Rθ, we have
E[ξm(θ)] ≤ λA 4 E[‖θ − θ∗‖22 |Fn,0] + 8[1 + (κ− 1)ρ] λA(1− ρ)M [R2θ(1 + γ) 2 + r2max],
where the expectation is over the random trajectory, θ is treated as a fixed variable, and 0 < C0 <∞ is a constant depending only on the MDP.
Lemma 1 shows that the bias error diminishes as the batch size M increases and the algorithm approaches to the fixed point θ∗. To explain why this happens, the definition of ξm(θ) immediately yields the following bound:
ξm(θ) ≤ 1
λA ‖gn(θ)− g(θ)‖22 + λA 4 ‖θ − θ∗‖22 . (5)
The first term on the right-hand-side of eq. (5) can be bounded by the concentration property for the ergodic process as gm(θ) = 1M ∑mM−1 i=(m−1)M gxi(θ)
a.s.→ g(θ). As M increases, the randomness due to the pseudo-gradient estimation is essentially averaged out due to the variance reduction step in VRTD, which implicitly eliminates its correlation from samples in the previous epochs.
As a comparison, the bias error in vanilla TD has been shown to be bounded by E[ξn(θ)] = O(α log(1/α)) Bhandari et al. (2018). In order to reduce the bias and achieve a high convergence accuracy, the stepsize α is required to be small, which causes the algorithm to run very slowly. The advantage of VRTD is that the bias can be reduced by choosing a sufficiently large batch size M so that the stepsize can still be kept at a desirable constant to guarantee fast convergence. Theorem 2. Consider the VRTD algorithm in Algorithm 2. Suppose Assumptions 1–3 hold. Set the constant stepsize α < λA12(1+γ)2 and the batch size M > 1 0.5αλA−6α2(1+γ)2 . Then, we have
E [∥∥∥θ̃m − θ∗∥∥∥2
2
] ≤ Cm1 ∥∥∥θ̃0 − θ∗∥∥∥2 2 + 3C4α+ C2/λA (1− C1)[0.5λA − 3α(1 + γ)2]M , (6)
where C1 = 1/M+3α2(1+γ)2
0.5αλA−3α2(1+γ)2 (with C1 < 1 due to the choices for α and M ), C2 = 16[1+(κ−1)ρ][R2θ(1+γ)
2+r2max] 1−ρ and C4 = [(1 + γ)Rθ + rmax] 2 + 2ρκG[(1+γ)Rθ+rmax]1−ρ .
We note that the convergence rate in eq. (6) can be written in a simpler form as E[||θ̃m − θ∗||2] ≤ Cm1 ||θ̃0 − θ∗||2 +O(1/M). Theorem 2 shows that VRTD (i.e., Algorithm 2) with Markovian samples converges to a neighborhood of θ∗ at a linear rate, and the size of the neighborhood (i.e., the convergence error) decays sublinearly with the batch size M . More specifically, the first term in the right-hand side of eq. (6) captures the linear convergence of the algorithm, the second term corresponds to the sum of the cumulative pseudogradient estimation error and the cumulative bias error. For the fixed stepsize, the total convergence error is dominated by the sum of those two error terms with the order O(1/M). Therefore, the variance reduction in Algorithm 2 reduces both the variance and the bias of the pseudo-gradient estimator.
Based on the convergence rate of VRTD characterized in Theorem 2 under Markovian sampling, we obtain the following bound on the corresponding computational complexity.
Corollary 2. Suppose Assumptions 1-3 hold. Let α = λA24(1+γ)2 and M = d( 32C2+4C4 3(1−C1) + 100(1 + γ)2) max{ 1 , 1 λ2A }e. Then, for any > 0, an -accuracy solution (i.e., E||θ̃m − θ∗||2 ≤ ) can be attained with at most m = dlog 2‖θ̃0−θ∗‖2 2
/ log 1 C1 e iterations. Correspondingly, the total number of
pseudo-gradient computations required by VRTD (i.e., Algorithm 2) under Markovian sampling to attain such an -accuracy solution is at most
O ( max {1 , 1
λ2A
} log 1 ) .
Proof. Given the values of α and M in the theorem, it can be easily checked that E||θ̃m − θ∗||2 ≤ for m = dlog 2‖θ̃0−θ∗‖2 2
/ log 1 C1 e. Then the total number of pseudo-gradient computations is given
by 2mM that yields the desired order given in the theorem.
As a comparison, consider the vanilla TD algorithm studied in Bhandari et al. (2018) with the constant stepsize α = O( / log(1/ )). Under Markovian sampling, it can be shown (see Appendix F.2) that vanilla TD requires O ( 1 λ2A log2 ( 1 )) pseudo-gradient computations in total to obtain an -accuracy
solution. Hence, in the Markovian setting, VRTD outperforms vanilla TD in terms of the total computational complexity by a factor of log 1 . To intuitively explain, we first note that the correlation among data samples in the Markovian case also causes a bias error in addition to the variance error. For VRTD, due to the variance reduction scheme, the bias and variance errors are kept at the same level (with respect to the batch size) so that the bias error does not cause order-level increase in the computational complexity for VRTD. However, for vanilla TD, the bias error dominates the variance error, which turns out to require more iterations to attain an -accurate solution, and yields an additional log 1 factor in the total complexity compared to VRTD.

heading: 5 EXPERIMENTS
text: In this section, we provide numerical results to verify our theoretical results. Note that in Appendix A, we provide further experiments on two problems in OpenAI Gym Brockman et al. (2016) and one experiment to demonstrate that VRTD is more sample-efficient than vanilla TD.
We consider an MDP with γ = 0.95 and |S| = 50. Each transition probability are randomly sampled from [0,1] and the transitions were normalized to one. The expected reward for each transition is also generated randomly in [0,1] and the reward on each transition was sampled without noise. Each component of the feature matrix Φ ∈ R50×4 is randomly and uniformly sampled between 0 and 1. The baseline for comparison is the vanilla TD algorithm, which corresponds to the case with M = 1 in our figure. We conduct two experiments to investigate how the batch size M for variance reduction affects the performance of VRTD with i.i.d. and Markovian samples. In the Markovian setting, we sample the data from a MDP trajectory. In the i.i.d. setting, we sample the data independently from the corresponding stationary distribution. In both experiments, we set the constant stepsize to be α = 0.1 and we run the experiments for five different batch sizes: M = 1, 50, 500, 1000, 2000. Our results are reported in Figure 1 and 2. All the plots report the square error over 1000 independent runs. In each case, the left figure illustrates the convergence process over the number of pseudo-gradient computations and the right figure shows the convergence errors averaged over the last 10000 iterations for different batch size values. It can be seen that in both i.i.d. and Markovian settings, the averaged error decreases as the batch size increases, which corroborates both Theorem 1 and Theorem 2. We also observe that increased batch size substantially reduces the error without much slowing down the convergence, demonstrating the desired advantage of variance reduction. Moreover, we observe that the error of VRTD with i.i.d samples is smaller than that of VRTD with Markovian samples under all batch size settings, which indicates that the correlation among Markovian samples introduces additional errors.

heading: 6 CONCLUSION
text: In this paper, we provided the convergence analysis for VRTD with both i.i.d. and Markovian samples. We developed a novel technique to bound the bias of the VRTD pseudo-gradient estimator. Our result demonstrate the advantage of VRTD over vanilla TD on the reduced variance and bias errors by the batch size. We anticipate that such a variance reduction technique and our analysis tools can be further applied to other RL algorithms.

heading: ACKNOWLEDGMENTS
text: The work was supported in part by US National Science Foundation under the grants CCF-1801855, ECCS-1818904, and CCF-1909291. The authors would like to thank Bowen Weng at the Ohio State University for the helpful discussions on the experiments. The authors would also like to thank a few anonymous reviewers for their suggestions on the analysis of the overall computational complexity as well as additional experiments, which significantly help to improve the quality of the paper.

heading: A ADDITIONAL EXPERIMENTS
text: In this section, we assess the practical performance of VRTD (Algorithm 2) on two problems in OpenAI Gym Brockman et al. (2016), which are Frozen Lake (4× 4) and Mountain Car. We also provide an additional experiment to demonstrate that VRTD is more sample-efficient than vanilla TD.
A.1 FROZEN LAKE
Frozen Lake is a game in OpenAI Gym, which is designed by an MDP with a finite state and action space. An agent starts from the starting point at t = 0 and can only transport to the neighbor blocks. It returns to the start-point every time it reaches a "hole" or the "goal". The agent receives a reward 1 only when it reaches the goal and 0 otherwise. Each transition probability is randomly sampled from [0, 1] and normalized to one, and each component of the feature matrix Φ ∈ R16×4 is also randomly sampled from [0, 1]. Given the feature matrix and the transition probability, the ground truth value of θ∗ can be calculated, which is used to evaluate the error in the experiments. We set the stepsize to be α = 0.1 and run vanilla TD (M = 1) and VRTD with the batch sizes M = 50, 500, 1000, 2000. Note that M = 1 corresponds to the base line vanilla TD. We compute the squared error over 1000 independent runs. The left plot in Figure 3 shows the convergence process over the number of pseudo-gradient computations and the right plot in Figure 3 shows the convergence error averaged over the last 10000 iterations. It can be observed that VRTD achieves much smaller error than TD, and increasing the batch size for VRTD substantially reduces the error without much slowing down the convergence.
A.2 MOUNTAIN CAR
Mountain Car is a game in OpenAI Gym, which is driven by an MDP with an infinite state space and a finite action space. At each time step, an agent randomly chooses an action ∈ {push left, push right, no push}. In this problem, the ground truth value of θ∗ is not known. In order to quantify the performance of VRTD, we apply the error metric known as the norm of the expected TD update given by NEU= ‖E[δφ]‖22, where δ is the temporal difference Sutton et al. (2009); Maei (2011). The state sample is transformed into a feature vector with the dimension 20 using an approximation of a RBF kernel. The agent follows a random policy in our experiment and we initialize θ0 = 0. At t = 0, the agent starts from the lowest point, receives a reward of −1 at each time step, and returns to the starting point every time it reaches the goal. We set the stepsize to be α = 0.2 and run vanilla TD (M = 1) and VRTD with batch size M = 1000. After every 10000 pseudo-gradient computations, learning is paused and the NEU is computed by averaging over 1000 test samples. We conduct 1000 independent runs and the results are reported by averaging over these
runs. Figure 4 shows the convergence process of the NEU versus the number of pseudo-gradient computations. It can been seen that VRTD achieves smaller NEU than vanilla TD.
A.3 COMPARISON BETWEEN VRTD AND TD WITH A CHANGING STEPSIZE
In this subsection, we provide an additional experiment to compare the performance of VRTD given in Algorithm 2 (under constant stepsize) with the TD algorithm (under a changing stepsize as suggested by the reviewer). We adopt the same setting of Frozen Lake as in Appendix A.1. Let VRTD take a batch size M = 5000 and stepsize α = 0.1. For a fair comparison, we start TD with the same constant stepsize α = 0.1 and then reduce the stepsize by half whenever the error stops decrease. The comparison is reported in Figure 5, where both curves are averaged over 1000 independent runs. The two algorithms are compared in terms of the squared error versus the total number of pseudo-gradient computations (equivalently, the total number of samples being used). It can be seen that VRTD reaches the required accuracy much faster than TD.

heading: B A COUNTER EXAMPLE
text: In this section, we use a counter-example to show that one major technical step for characterizing the convergence bound in Korda and La (2015) does not hold. Consider Step 4 in the proof of Theorem 3 in Korda and La (2015). For the following defined (θ)
(θ) = (θ − θ∗)>[E(v>v|Fn)− EΨ,θn(v>v)](θ − θ∗), (7)
where Ψ denotes the stationary distribution of the corresponding Markov chain, Korda and La (2015) claimed that the following inequality holds
‖ (θ)‖2 ≤ 2H ‖E(v|Fn)− EΨ,θn(v)‖2 . (8) This is not correct. Consider the following counter-example. Let the batch size M = 3 and the dimension of the feature vector be one, i.e., Φ ∈ R|S|×1. Hence, all variables in eq. (8) and eq. (7) are scalars. Since the steps for proving eq. (8) in Korda and La (2015) do not have specific requirements for the transition kernel, eq. (8) should hold for any distribution of v. Thus, suppose v follows the uniform distribution over [−3, 3]. Further assume that in the n-th epoch, the samples of v are given by {1, 2,−3}. Recall that E(·|Fn) is the average over the batch samples in the n-th epoch. We have:
EΨ,θn(v) = 0, EΨ,θn(v2) = 3, E(v|Fn) = 0, E(v2|Fn) = 14
3 .
Substituting the above values into eq. (8) yields ‖ (θ)‖2 = (14 3 − 3 ) (θ − θ∗)2 ≤ 2H × 0 = 0, (9) which obviously does not hold in general when θ 6= θ∗. Consequently the second statement in Theorem 3 of Korda and La (2015), which is critically based on the above erroneous steps, does not hold. Hence, the first statement in the same theorem whose proof is based on the second statement cannot hold either.

heading: C USEFUL LEMMAS
text: In the rest of the paper, for any matrix W ∈ Rd×d, we denote ‖W‖2 as the spectral norm of W and ‖W‖F as the Frobenius norm of W . Lemma 2. For any xi = (si, ri, s′i) (i.i.d. sample) or xi = (si, ri, si+1) (Markovian sample), we have ‖Axi‖2 ≤ 1 + γ and ‖bxi‖2 ≤ rmax.
Proof. First consider the case when samples are i.i.d. Due to the definition of Axi , we have ‖Axi‖2 = ∥∥φ(si)(γφ(s′i)− φ(si))>∥∥2
≤ ∥∥φ(si)(γφ(s′i)− φ(si))>∥∥F
≤ γ ∥∥φ(si)φ(s′i)>∥∥F + ∥∥φ(si)φ(si)>∥∥F
≤ 1 + γ. Then, consider bxi :
‖bxi‖2 = ‖rxiφ(si)‖2 ≤ rmax ‖φ(si)‖2 ≤ rmax. Following similar steps, we can obtain the same upper bounds for the case with Markovian samples.
Lemma 3. Let G = (1 + γ)Rθ + rmax. Consider Algorithm 2. For any m > 0 and 0 ≤ t ≤M − 1, we have ∥∥gxjm,t(θm,t)∥∥2 , ∥∥∥gxjm,t(θ̃m−1)∥∥∥2 , ∥∥∥gm(θ̃m−1)∥∥∥2 ≤ G. Proof. First, we bound
∥∥gxjm,t(θm,t)∥∥2 as follows.∥∥gxjm,t(θm,t)∥∥2 = ∥∥Axjm,tθm,t + bθm,t∥∥2 ≤ ∥∥Axjm,t∥∥2 ‖θm,t‖2 + ∥∥bθm,t∥∥2
≤ (1 + γ)Rθ + rmax. Following the steps similar to the above, we have ∥∥∥gxjm,t(θ̃m−1)∥∥∥
2 ≤ G. Finally for∥∥∥gxjm,t(θ̃m−1)∥∥∥
2 , we have
∥∥∥gxjm,t(θ̃m−1)∥∥∥ 2 = ∥∥∥∥∥∥ 1M mM−1∑
i=(m−1)M
gxi(θ̃m−1) ∥∥∥∥∥∥ 2
≤ 1 M mM−1∑ i=(m−1)M ∥∥∥gxi(θ̃m−1)∥∥∥ 2
≤ G, (10) where eq. (10) follows from the last fact ∥∥∥gxjm,t(θ̃m−1)∥∥∥
2 ≤ G.
Lemma 4. Define D1 = 2(1 + γ)2 and D2 = 4((1 + γ)2R2θ + r2max). For any θ ∈ Rd, we have ‖gxi(θ)‖ 2 2 ≤ D1 ‖θ − θ∗‖ 2 2 +D2.
Proof. Recalling the definition of gxi , and applying Lemma 2, we have
‖gxi(θ)‖ 2 2 = ‖Axiθ + bxi‖ 2 2
= ‖Axi(θ − θ∗) + (Axiθ∗ + bxi)‖ 2 2 ≤ 2 ‖Axi(θ − θ∗)‖ 2 2 + 2 ‖Axiθ ∗ + bxi‖ 2 2 ≤ 2 ‖Axi‖ 2 2 ‖θ − θ ∗‖22 + 4(‖Axi‖ 2 2 ‖θ ∗‖22 + ‖bxi‖ 2 2) ≤ 2(1 + γ)2 ‖θ − θ∗‖22 + 4((1 + γ) 2R2θ + r 2 max) = D1 ‖θ − θ∗‖22 +D2.
Lemma 5. Considering Algorithm 2 with Markovian samples. We have ‖E[Aj |Pi]−A‖F ≤ (1 + γ)κρj−i and ‖E[bj |Pi]− b‖2 ≤ rmaxκρ j−i for 0 < i < j.
Proof. We first derive ‖E[Aj |Pi]−A‖F = ∥∥∥∥∫ AxidP (xi|Pj)− ∫ Axidµπ∥∥∥∥
F ≤ ∫ ‖AxidP (xi|Pj)−Axidµπ‖F
≤ ∫ ‖Axi‖F |dP (xi|Pj)− dµπ|
≤ (1 + γ) ‖P (xi|Pj), µπ‖TV ≤ (1 + γ)κρj−i.
Following the steps similar to the above, we can derive ‖E[bj |Pi]− b‖2 ≤ 2rmaxκρ j−i.

heading: D PROOF OF THEOREM 1: CONVERGENCE OF VRTD WITH I.I.D. SAMPLES
text: Recall that Bm is the sample batch drawn at the beginning of each m-th epoch and xi,j denotes the sample picked at the j-th iteration in the i-th epoch in Algorithm 1. We denote σ(θ̃0) as a trivial σ-field when θ̃0 is a deterministic vector. Let σ(A ∪ B) indicate the smallest σ-field that contains both A and B. Then, we construct a set of σ-fields in the following incremental way.
F1,0 = σ(θ̃0), F1,1 = σ(F1,0 ∪ σ(B1) ∪ σ(x1,1)), ..., F1,M = σ(F1,(M−1) ∪ σ(x1,M )), F2,0 = σ(F1,M ∪ σ(θ̃1)), F21 = σ(F2,0 ∪ σ(B2) ∪ σ(x2,1)), ..., F2,m = σ(F2,(M−1) ∪ σ(x2,M )), ...
Fm,0 = σ(F(m−1),M ∪ σ(θ̃m−1)), Fm1 = σ(Fm,0 ∪ σ(Bm) ∪ σ(xm,1)), ..., Fm,M = σ(Fm,(M−1) ∪ σ(xm,M )).
The proof of Theorem 1 proceeds along the following steps.
Step 1: Iteration within the m-th epoch
For the m-th epoch, we consider the last update (i.e., the M -th iteration in the epoch), and decompose its error into the following form.
‖θm,M − θ∗‖22 = ∥∥∥θm,M−1 + α(gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1))− θ∗∥∥∥22
= ‖θm,M−1 − θ∗‖22 + 2α(θm,M−1 − θ ∗)> ( gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1) ) + α2
∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)∥∥∥22 . (11) First, consider the third term in the right-hand side of eq. (11), we have∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)∥∥∥22 ≤ 2
∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + g(θ̃m−1)∥∥∥22 + 2 ∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥22 = 2
∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ∗)− [(gxjm,M (θ̃m−1)− gxjm,M (θ∗))− (g(θ̃m−1)− g(θ∗))]∥∥∥22 + 2 ∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2 2
≤ 4 ∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ∗)∥∥∥22 + 4 ∥∥∥(gxjm,M (θ̃m−1)− gxjm,M (θ∗))− (g(θ̃m−1)− g(θ∗))∥∥∥22
+ 2 ∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2
2 . (12)
Then, by taking the expectation conditioned on Fm,M−1 on both sides of eq. (12), we have E [∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)∥∥∥22 ∣∣∣Fm,M−1 ] (i)
≤ 4E [∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ∗)∥∥∥22 ∣∣∣Fm,M−1 ] + 4E [∥∥∥(gxjm,M (θ̃m−1)− gxjm,M (θ∗))− E[gxjm,M (θ̃m−1)− gxjm,M (θ∗)∣∣Fm,M−1]∥∥∥22 ∣∣∣Fm,M−1 ]
+ 2E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2
2 ∣∣∣Fm,M−1] (ii)
≤ 4(1 + γ)2E [ ‖θm,M−1 − θ∗‖22 |Fm,M−1 ] + 4(1 + γ)2E [∥∥∥θ̃m−1 − θ∗∥∥∥2 2 |Fm,M−1 ] + 2E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2 2
∣∣∣Fm,M−1] where (i) follows from the fact that E[ ( gxjm,M (θ̃m−1)−gxjm,M (θ ∗) ) |Fm,M−1] = g(θ̃m−1)−g(θ∗), and (ii) follows from the inequality E[(X − EX)2] ≤ EX2 and Lemma 2. Then, taking the expectation conditioned on Fm,M−1 on both sides of eq. (11) yields
E [ ‖θm,M − θ∗‖22 ∣∣∣Fm,M−1] = ‖θm,M−1 − θ∗‖22 + 2α(θm,M−1 − θ ∗)>E [ gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)
∣∣∣Fm,M−1] + α2E [∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)∥∥∥22 ∣∣∣Fm,M−1 ]
(i) ≤ ‖θm,M−1 − θ∗‖22 + 2α(θm,M−1 − θ ∗)>g(θm,M−1) + 2α(θm,M−1 − θ∗)> ( E [ gm(θ̃m−1) ∣∣∣Fm,M−1]− g(θ̃m−1)) + 4α2(1 + γ)2 ‖θm,M−1 − θ∗‖22 + 4α 2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2
2
+ 2α2E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2
2 ∣∣∣Fm,M−1] (ii)
≤ ‖θm,M−1 − θ∗‖22 − αλA ‖θm,M−1 − θ ∗‖22 + 2αE [ ξm(θ̃m−1) ∣∣∣Fm,M−1] + 4α2(1 + γ)2 ‖θm,M−1 − θ∗‖22 + 4α 2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2
2 + 2α2E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2
2 ∣∣∣Fm,M−1] (iii)
≤ ‖θm,M−1 − θ∗‖22 − [αλA − 4α 2(1 + γ)2] ‖θm,M−1 − θ∗‖22 + 4α 2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2
2 + 2αE [ ξm(θ̃m−1) ∣∣∣Fm,M−1]+ 2α2E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2 2 ∣∣∣Fm,M−1] , (13) where (i) follows from the fact that E [ gxjm,M (θ̃m−1)
∣∣∣Fm,M−1] = g(θ̃m−1). In (ii) we define λA as the absolute value of the largest eigenvalue of matrix (AT + A), which is negative definite according to Tsitsiklis and Van Roy (1997). In (iii) we define ξm(θ) = (θ − θ∗)>(gm(θ)− g(θ)) for θ ∈ Rd. Then, by applying eq. (13) iteratively, we have
E [ ‖θm,1 − θ∗‖22 ∣∣∣Fm,0] ≤ ‖θm,0 − θ∗‖22 − [αλA − 4α 2(1 + γ)2]
M−1∑ i=0 E [ ‖θm,i − θ∗‖22 ∣∣∣Fm,0]+ 4Mα2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2
+ 2αME [ ξm(θ̃m−1) ∣∣∣Fm,0]+ 2Mα2E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2 2 ∣∣∣Fm,0] . (14) For all 1 ≤ i ≤M , we have
E [ ξm(θ̃m−1) ∣∣∣Fm,0] = E [gm(θ̃m−1)∣∣∣Fm,0]− g(θ̃m−1) = 1
M ∑ i∈Bm E[Axi θ̃m + bxi |Fm,0]− (Aθ̃m + b)
= [( 1 M ∑ i∈Bm E[Axi |Fm,0] ) −A ] θ̃m + [( 1 M ∑ i∈Bm E[bxi |Fm,0] ) − b ]
= 0.
Then, arranging terms in eq. (14) and using the above fact yield
[αλA − 4α2(1 + γ)2] M−1∑ i=0 E [ ‖θm,i − θ∗‖22 ∣∣∣Fm,0] ≤ [1 + 4Mα2(1 + γ)2] ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 + 2Mα2E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2 2
∣∣∣Fm,0] . (15) Finally, dividing eq. (15) by [αλA − 4α2(1 + γ)2]M on both sides yields
E [∥∥∥θ̃m − θ∗∥∥∥2
2 ∣∣∣Fm,0] ≤ 1/M + 4α 2(1 + γ)2
αλA − 4α2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 +
2α λA − 4α(1 + γ)2 E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2 2 ∣∣∣Fm,0] . (16)
Step 2: Bounding the variance error
For any 0 ≤ k ≤ m− 1, we have E [∥∥∥gm(θ̃m−1)− g(θ̃m−1)∥∥∥2
2
∣∣∣Fm,0] (17)
= E ∥∥∥∥∥ 1M ∑ i∈Bm gxi(θ̃m−1)− g(θ̃m−1) ∥∥∥∥∥ 2
2
∣∣∣Fm,0  = 1
M2 E ∥∥∥∥∥ ∑ i∈Bm gxi(θ̃m−1)− g(θ̃m−1) ∥∥∥∥∥ 2
2
∣∣∣Fm,0 
= 1
M2 E ( ∑ i∈Bm gxi(θ̃m−1)− g(θ̃m−1) )>( ∑ j∈Bm gxj (θ̃m−1)− g(θ̃m−1) )∣∣∣Fm,0  = 1
M2 ∑ i∈Bm ∑ j∈Bm E [〈 gxi(θ̃m−1)− g(θ̃m−1), gxj (θ̃m−1)− g(θ̃m−1) 〉∣∣∣Fm,0] = 1
M2 ∑ i=j E [∥∥∥gxi(θ̃m−1)− g(θ̃m−1)∥∥∥2 2 ∣∣∣Fm,0] = 1
M2 ∑ i=j E [∥∥∥gxi(θ̃m−1)− E[gxi(θ̃m−1)∣∣∣Fm,0]∥∥∥ 2 ∣∣∣Fm,0] ≤ 1 M2 ∑ i=j E [∥∥∥gxi(θ̃m−1)∥∥∥2 2 ∣∣∣Fm,0] ≤ 1 M ( D1 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 +D2 ) , (18)
where eq. (18) follows from Lemma 4.
Step 3: Iteration over m epoches
First, we substitute eq. (18) into eq. (16) to obtain E [∥∥∥θ̃m − θ∗∥∥∥2
2 ∣∣∣Fm,0] ≤ C1 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 + 2D2α (λA − 4α(1 + γ)2)M , (19)
where we define C1 = ( 4α(1 + γ)2 + 2D1α 2+1
αM
) 1
λA−4α(1+γ)2 .
Taking the expectation of eq. (19) conditioned on Fm−1,0 and following the steps similar to those in step 1 to upper bound E [∥∥∥θ̃m−1 − θ∗∥∥∥2
2 ∣∣∣Fm−1,0], we obtain E [∥∥∥θ̃m − θ∗∥∥∥2
2 ∣∣∣Fm−1,0] ≤ C1E [∥∥∥θ̃m−1 − θ∗∥∥∥2 2 ∣∣∣Fm−1,0]+ 2D2α (λA − 4α(1 + γ)2)M
≤ C21 ∥∥∥θ̃m−2 − θ∗∥∥∥2
2 +
2D2α
(λA − 4α(1 + γ)2)M 1∑ k=0 Ck1 .
Then, by following the above steps for (m− 1) times, we have E [∥∥∥θ̃m − θ∗∥∥∥2
2
] ≤ Cm1 ∥∥∥θ̃0 − θ∗∥∥∥2 2 + 2D2α (λA − 4α(1 + γ)2)M m−1∑ k=0 Ck1
≤ Cm1 ∥∥∥θ̃0 − θ∗∥∥∥2
2 +
2D2α
(1− C1)(λA − 4α(1 + γ)2)M ,
which yields the desirable result.

heading: E PROOF OF THEOREM 2: CONVERGENCE OF VRTD WITH MARKOVIAN SAMPLES
text: We define σ(Si) to be the σ-field of all samples up to the i-th epoch and recall that jm,t is the index of the sample picked at the t-th iteration in the m-th epoch in Algorithm 2. Then we define a set of σ-fields in the following incremental way: F1,0 = σ(S0), F1,1 = σ(F1,0 ∪ σ(j1,1)), ..., F1,M = σ(F1,(M−1) ∪ σ(j1,M )), F2,0 = σ(σ(S1) ∪ F1,M ∪ σ(θ̃1)), F21 = σ(F2,0 ∪ σ(j2,1)), ..., F2,m = σ(F2,(M−1) ∪ σ(j2,M )), ...
Fm,0 = σ(σ(Sm−1) ∪ F(m−1),M ∪ σ(θ̃m−1)), Fm,1 = σ(Fm,0 ∪ σ(jm,1)), ..., Fm,M = σ(Fm,(M−1) ∪ σ(jm,M )).
E.1 PROOF OF LEMMA 1
We first prove Lemma 1, which is useful for step 4 in the main proof in Theorem 2 provided in Section E.2.
Proof. Recall the definition of the bias term: ξn(θ) = (θ − θ∗)>(gn(θ)− g(θ)). We have
E [ξn(θ)| Fn,0] = E [ (θ − θ∗)>(gn(θ)− g(θ))|Fn,0 ] = E (θ − θ∗)>[( 1 M nM−1∑ i=(n−1)M Axi −A ) θ + ( 1 M nM−1∑ i=(n−1)M bxi − b )]∣∣∣∣∣Fn,0
 ≤ λA
4 E[‖θ − θ∗‖22 |Fn,0] +
1
λA E  ∥∥∥∥∥∥ ( 1 M nM−1∑ i=(n−1)M Axi −A ) θ + ( 1 M nM−1∑ i=(n−1)M bxi − b )∥∥∥∥∥∥ 2
2
∣∣∣∣∣Fn,0 
≤ λA 4 E[‖θ − θ∗‖22 |Fn,0] + 2 λA E  ∥∥∥∥∥∥ ( 1 M nM−1∑ i=(n−1)M Axi −A ) θ ∥∥∥∥∥∥ 2
2
∣∣∣∣∣Fn,0 
+ 2
λA E  ∥∥∥∥∥∥ 1M nM−1∑ i=(n−1)M bxi − b ∥∥∥∥∥∥ 2
2
∣∣∣∣∣Fn,0 
≤ λA 4 E[‖θ − θ∗‖22 |Fn,0] + 2R2θ λA E  ∥∥∥∥∥∥ 1M nM−1∑ i=(n−1)M Axi −A ∥∥∥∥∥∥ 2
2
∣∣∣∣∣Fn,0 
+ 2
λA E  ∥∥∥∥∥∥ 1M nM−1∑ i=(n−1)M bxi − b ∥∥∥∥∥∥ 2
2
∣∣∣∣∣Fn,0 
(i) ≤ λA 4 E[‖θ − θ∗‖22 |Fn,0] + 2R2θ λA E  ∥∥∥∥∥∥ 1M nM−1∑ i=(n−1)M Axi −A ∥∥∥∥∥∥ 2
F
∣∣∣∣∣Fn,0 
+ 2
λA E  ∥∥∥∥∥∥ 1M nM−1∑ i=(n−1)M bxi − b ∥∥∥∥∥∥ 2
2
∣∣∣∣∣Fn,0  , (20)
where (i) follows from the fact that ‖W‖2 ≤ ‖W‖F for all matrices W ∈ Rd×d. We define the interproduct between two maticesW,V ∈ Rd×d as 〈W,V 〉 = ∑ ij ∑ ijWijVij . Consider the second
term in eq. (20): E [∥∥∥ 1M ∑nM−1i=(n−1)M Axi −A∥∥∥2
F ∣∣∣Fn,0], we have E  ∥∥∥∥∥∥ 1M nM−1∑ i=(n−1)M Axi −A ∥∥∥∥∥∥ 2
F
∣∣∣∣∣Fn,0 
= 1
M2 nM−1∑ i=(n−1)M nM−1∑ j=(n−1)M E[〈Axi −A,Axj −A〉|Fn,0]
= 1
M2 [∑ i=j E[‖Axi −A‖ 2 F |Fn,0] + ∑ i 6=j E[〈Axi −A,Axj −A〉|Fn,0] ]
≤ 1 M2 [∑ i=j E[(‖Axi‖F + ‖A‖F ) 2|Fn,0] + ∑ i6=j E[〈Axi −A,Axj −A〉|Fn,0] ]
≤ 1 M2
[ 4(1 + γ)2M + ∑ i 6=j E[〈Axi −A,Axj −A〉|Fn,0]. ]
(21)
Now consider upper bound the term E[〈Axi − A,Axj − A〉|Fn,0] in eq. (21). Without loss of generality, we consider the case when i > j as follows:
E[〈Axi −A,Axj −A〉|Fn,0] = E [ E [ 〈Axi −A,Axj −A〉|xj ]∣∣∣Fn,0] = E [ 〈E [ Axi |xj ] −A,Axj −A〉
∣∣∣Fn,0] ≤ E
[ ∥∥E[Axi |xj]−A∥∥F ∥∥Axj −A∥∥F ∣∣∣Fn,0] ≤ E
[ ∥∥E[Axi |xj]−A∥∥F ( ∥∥Axj∥∥F + ‖A‖F )∣∣∣Fn,0] ≤ 2(1 + γ)E
[ ∥∥E[Axi |xj]−A∥∥F ∣∣∣Fn,0] ≤ 2κ(1 + γ)2ρi−j .
We can further obtain∑ i6=j E[〈Axi −A,Axj −A〉|Fn,0] ≤ 2κ(1 + γ)2 ∑ i 6=j ρ|i−j| ≤ 2κ(1 + γ)2 M−1∑ k=1 k∑ l=1 ρl
≤ 2κ(1 + γ)2 2ρ 1− ρ M−1∑ k=1 (1− ρk) ≤ 4(1 + γ) 2Mκρ 1− ρ . (22)
Then substituting eq. (22) into eq. (21) yields
E  ∥∥∥∥∥∥ 1M nM−1∑ i=(n−1)M Axi −A ∥∥∥∥∥∥ 2
F
∣∣∣∣∣Fn,0  ≤ 1 M [ 4(1 + γ)2 + 4(1 + γ)2κρ 1− ρ ] .
Then consider the third term in eq. (20): E [∥∥∥ 1M ∑nM−1i=(n−1)M bxi − b∥∥∥2 2 ∣∣∣∣∣Fn,0 ] , similarly we have
E  ∥∥∥∥∥∥ 1M nM−1∑ i=(n−1)M bxi − b ∥∥∥∥∥∥ 2
2
∣∣∣∣∣Fn,0 
= 1
M2 nM−1∑ i=(n−1)M nM−1∑ j=(n−1)M E[(bxi − b)>(bxj − b)|Fn,0]
= 1
M2 [∑ i=j E[‖bxi − b‖ 2 2 |Fn,0] + ∑ i6=j E[(bxi − b)>(bxj − b)|Fn,0] ]
≤ 1 M2 [∑ i=j E[(‖bxi‖2 + ‖b‖2) 2|Fn,0] + ∑ i 6=j E[(bxi − b)>(bxj − b)|Fn,0] ]
≤ 1 M2
[ 4r2maxM + ∑ i6=j E[(bxi − b)>(bxj − b)|Fn,0] (23)
Now to upper-bound the term E[(bxi − b)>(bxj − b)|Fn,0], without loss of generality, we assume i > j
E[(bxi − b)>(bxj − b)|Fn,0]
= E [ E [ (bxi − b)>(bxj − b)|xj ]∣∣∣Fn,0] = E [ (E [ bxi |xj ] − b)>(bxj − b)
∣∣∣Fn,0] ≤ E
[ ∥∥E[bxi |xj]− b∥∥2 ∥∥bxj − b∥∥2 ∣∣∣Fn,0] ≤ E
[ ∥∥E[bxi |xj]− b∥∥2 ( ∥∥bxj∥∥2 − ‖b‖2 )∣∣∣Fn,0] ≤ 2rmaxE
[ ∥∥E[bxi |xj]− b∥∥2 ∣∣∣Fn,0] ≤ 2κr2maxρi−j .
Thus, we have∑ i 6=j E[(bxi − b)>(bxj − b)|Fn,0] ≤ 2κr2max ∑ i6=j ρ|i−j| ≤ 2κr2max M−1∑ k=1 k∑ l=1 ρl
≤ 2κr2max 2ρ
1− ρ M−1∑ k=1 (1− ρk) ≤ 4r 2 maxMκρ 1− ρ . (24)
Combing all of the above pieces together yields the following final bound
E [ξn(θ)| Fn,0] ≤ λA 4 E[‖θ − θ∗‖22 |Fn,0] + 8[1 + (κ− 1)ρ] λA(1− ρ)M [R2θ(1 + γ) 2 + r2max]. (25)
E.2 PROOF OF THEOREM 2
Step 1: Iteration within the m-th inner loop
For the m-th inner loop, we consider the last update (i.e., the M -th iteration in the epoch), and decompose its error into the following form.
‖θm,M − θ∗‖22 = ∥∥∥ΠRθ (θm,M−1 + α(gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)))− θ∗∥∥∥22
≤ ∥∥∥θm,M−1 + α(gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1))− θ∗∥∥∥22
= ‖θm,M−1 − θ∗‖22 + 2α(θm,M−1 − θ ∗)> ( gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1) ) + α2
∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)∥∥∥22 . (26) First, consider the third term in the right-hand side of eq. (26).∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)∥∥∥22 = ∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ∗)− [(gxjm,M (θ̃m−1)− gxjm,M (θ∗))− (gm(θ̃m−1)− gm(θ∗))]+ gm(θ∗)∥∥∥22
≤ 3 ∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ∗)∥∥∥22 + 3 ∥∥∥(gxjm,M (θ̃m−1)− gxjm,M (θ∗))− (gm(θ̃m−1)− gm(θ∗))∥∥∥22
+ 3 ‖gm(θ∗)‖22 . (27) Then, by taking the expectation conditioned on Fm,(M−1) on both sides of eq. (27), we have
E [∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)∥∥∥22 ∣∣∣Fm,M−1 ] (i)
≤ 3E [∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ∗)∥∥∥22 ∣∣∣Fm,M−1 ] + 3E [∥∥∥(gxjm,M (θ̃m−1)− gxjm,M (θ∗))− E[gxjm,M (θ̃m−1)− gxjm,M (θ∗)|Fm,M−1]∥∥∥22 ∣∣∣Fm,M−1 ]
+ 3E [ ‖gm(θ∗)‖22 ∣∣∣Fm,M−1] ≤ 3E [∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ∗)∥∥∥22 ∣∣∣Fm,M−1 ] + 3E [∥∥∥gxjm,M (θ̃m−1)− gxjm,M (θ∗)∥∥∥22 ∣∣∣Fm,M−1
] (ii)
≤ 3E [ ‖Am,M‖22 ‖θm,M−1 − θ ∗‖22 ∣∣∣Fm,M−1]+ 3E [‖Am,M‖22 ∥∥∥θ̃m−1 − θ∗∥∥∥2
2 ∣∣∣Fm,M−1] + 3E [ ‖gm(θ∗)‖22
∣∣∣F1,0] (iii)
≤ 3(1 + γ)2 ‖θm,M−1 − θ∗‖22 + 3(1 + γ) 2 ∥∥∥θ̃m−1 − θ∗∥∥∥2
2 + 3E
[ ‖gm(θ∗)‖22 ∣∣∣F1,0] (28) where (i) follows from the fact that E[ ( gxjm,M (θ̃m−1) − gxjm,M (θ ∗) ) |Fm,M−1] = gm(θ̃m−1) − gm(θ ∗), (ii) follows from the inequality E[(X − EX)2] ≤ EX2, and (iii) follows from Lemma 2. We further consider the last term in eq. (28):
E [ ‖gm(θ∗)‖22 ∣∣∣F1,0] = ∥∥∥∥∥∥ ( 1 M mM−1∑ i=(m−1)M Ai ) θ∗ + ( 1 M mM−1∑ i=(m−1)M bi )∥∥∥∥∥∥ 2
2
.
Then, taking the expectation conditioned on Fm,M−1 on both sides of eq. (26) yields E [ ‖θm,M − θ∗‖22 ∣∣∣Fm,M−1] ≤ ‖θm,M−1 − θ∗‖22 + 2α(θm,M−1 − θ ∗)>E [ gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + gm(θ̃m−1)
∣∣∣Fm,M−1] + α2E [∥∥∥gxjm,M (θm,M−1)− gxjm,M (θ̃m−1) + g̃m∥∥∥22 ∣∣∣Fm,M−1 ]
(i) ≤ ‖θm,M−1 − θ∗‖22 + 2α(θm,M−1 − θ ∗)>E [ gxjm,M (θm,M−1) ∣∣∣Fm,M−1]+ 3α2(1 + γ)2 ‖θm,M−1 − θ∗‖22 + 3α2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 + 3α2E [ ‖gm(θ∗)‖22
∣∣∣F1,0] (ii) = ‖θm,M−1 − θ∗‖22 + 2α(θm,M−1 − θ ∗)>g(θm,M−1) + 2αE [ ξm(θm,M−1)
∣∣∣Fm,M−1] + 3α2(1 + γ)2 ‖θm,M−1 − θ∗‖22 + 3α 2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2
2 + 3α2E
[ ‖gm(θ∗)‖22 ∣∣∣F1,0] ≤ ‖θm,M−1 − θ∗‖22 − [αλA − 3α 2(1 + γ)2] ‖θm,M−1 − θ∗‖22 + 3α 2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2
+ 2αE [ ξm(θm,M−1) ∣∣∣Fm,M−1]+ 3α2E [‖gm(θ∗)‖22 ∣∣∣F1,0] , (29) where (i) follows by plugging eq. (28) into its preceding step and from the fact that E [ gxjm,M (θ̃m−1)− gm(θ̃m−1)
∣∣∣Fm,M−1] = 0. In (ii) we define ξm(θ) = (θ−θ∗)>(gm(θ)−g(θ)) for θ ∈ Rd. Then, by applying eq. (29) iteratively, we have
E [ ‖θm,1 − θ∗‖22 ∣∣∣Fm,0] ≤ ‖θm,0 − θ∗‖22 − [αλA − 3α 2(1 + γ)2]
M−1∑ i=0 E [ ‖θm,i − θ∗‖22 ∣∣∣Fm,0]+ 3Mα2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2
+ 2α M−1∑ i=0 E [ ξm(θm,i) ∣∣∣Fm,0]+ 3Mα2E [‖gm(θ∗)‖22 ∣∣∣F1,0] . (30) Arranging the terms in eq. (30) yields [αλA − 3α2(1 + γ)2] M−1∑ i=0 E [ ‖θm,i − θ∗‖22 ∣∣∣Fm,0]
≤ [1 + 3Mα2(1 + γ)2] ∥∥∥θ̃m−1 − θ∗∥∥∥2
2 + 2α M−1∑ i=0 E [ ξm(θm,i) ∣∣∣Fm,0]+ 3Mα2E [‖gm(θ∗)‖22 ∣∣∣F1,0] . (31)
Then, substituting eq. (25) into eq. (31), we obtain
[αλA − 3α2(1 + γ)2] M−1∑ i=0 E [ ‖θm,i − θ∗‖22 ∣∣∣Fm,0] ≤ [1 + 3Mα2(1 + γ)2] ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 + λAα 2 M−1∑ i=0 E [ ‖θm,i − θ∗‖22 |Fm,0
] +
16[1 + (κ− 1)ρ]α λA(1− ρ) [R2θ(1 + γ) 2 + r2max] + 3Mα
2E [ ‖gm(θ∗)‖22 ∣∣∣F1,0] . (32) Subtracting 0.5λAα ∑M−1 i=0 E[‖θm,i − θ∗‖ 2 2 |Fm,0] on both sides of eq. (32) yields
[0.5αλA − 3α2(1 + γ)2] M−1∑ i=0 E [ ‖θm,i − θ∗‖22 ∣∣∣Fm,0] ≤ [1 + 3Mα2(1 + γ)2] ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 + 16[1 + (κ− 1)ρ]α λA(1− ρ) [R2θ(1 + γ) 2 + r2max]
+ 3Mα2E [ ‖gm(θ∗)‖22 ∣∣∣F1,0] . (33) Then, dividing eq. (31) by [0.5αλA − 3α2(1 + γ)2]M on both sides, we obtain
E [∥∥∥θ̃m − θ∗∥∥∥2
2 ∣∣∣Fm,0] ≤ 1/M + 3α 2(1 + γ)2
0.5αλA − 3α2(1 + γ)2 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 + 16[1 + (κ− 1)ρ][R2θ(1 + γ)2 + r2max]α λA(1− ρ)[0.5αλA − 3α2(1 + γ)2]M
+ 3α 0.5λA − 3α(1 + γ)2 E [ ‖gm(θ∗)‖22 ∣∣∣F1,0] . (34) For simplicity, let C1 = 1/M+3α2(1+γ)2 0.5αλA−3α2(1+γ)2 , C2 = 16[1+(κ−1)ρ][R2θ(1+γ) 2+r2max] 1−ρ and C3 = 3α 0.5λA−3α(1+γ)2 . Then we rewrite eq. (34):
E [∥∥∥θ̃m − θ∗∥∥∥2
2 ∣∣∣Fm,0] ≤ C1 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 + C2 [0.5λA − 3α(1 + γ)2]λAM + C3E [ ‖gm(θ∗)‖22 ∣∣∣F1,0] . (35)
Step 2: Iteration over m epochs Taking the expectation of eq. (35) conditioned on Fm−1,0 and upper-bounding E [∥∥∥θ̃m−1 − θ∗∥∥∥2
2 ] by following similar steps in the previous steps, we obtained
E [∥∥∥θ̃m − θ∗∥∥∥2
2 ∣∣∣Fm−1,0] ≤ C1 ∥∥∥θ̃m−1 − θ∗∥∥∥2 2 + C2 [0.5λA − 3α(1 + γ)2]λAM + C3E [ ‖gm(θ∗)‖22
∣∣∣F1,0] ≤ C21 ∥∥∥θ̃m−2 − θ∗∥∥∥2 2 + C2 [0.5λA − 3α(1 + γ)2]λAM 1∑ k=0 Ck1 + C3 1∑ k=0 Ck1E [ ‖gm−k(θ∗)‖22
∣∣∣F1,0] . By following the above steps for (m− 1) times, we have
E [∥∥∥θ̃m − θ∗∥∥∥2
2
∣∣∣F1,0]
≤ Cm1 ∥∥∥θ̃0 − θ∗∥∥∥2
2 + C2 [0.5λA − 3α(1 + γ)2]λAM m−1∑ k=0 Ck1 + C3 m−1∑ k=0 Ck1E [ ‖gm−k(θ∗)‖22 ∣∣∣F1,0] . (36)
Then taking the expectation of σ(S) (which contains the randomness of the entire sample trajectory) on both sides of eq. (36) yields
E [∥∥∥θ̃m − θ∗∥∥∥2
2 ] ≤ Cm1 ∥∥∥θ̃0 − θ∗∥∥∥2 2 + C2 [0.5λA − 3α(1 + γ)2]λAM m−1∑ k=0 Ck1 + C3 m−1∑ k=0 Ck1E [ ‖gm−k(θ∗)‖22 ] ,
(37)
where the second term in the right hand side of eq. (37) corresponds to the bias error and the third term corresponds to the variance error.
Step 3: Bounding the variance error
For any 0 ≤ k ≤ m− 1, we have
‖gm−k(θ∗)‖22 = ∥∥∥∥∥∥ 1M (m−k)M−1∑ i=(m−k−1)M gxi(θ ∗) ∥∥∥∥∥∥ 2
2
= 1
M2 ( (m−k)M−1∑ i=(m−k−1)M g>xi(θ ∗) )( (m−k)M−1∑ j=(m−k−1)M gxj (θ ∗) )
= 1
M2 (m−k)M−1∑ i=(m−k−1)M (m−k)M−1∑ j=(m−k−1)M g>xi(θ ∗)gxj (θ ∗)
= 1
M2 ∑ i=j ‖gxi(θ∗)‖ 2 2 + 1 M2 ∑ i 6=j g>xi(θ ∗)gxj (θ ∗)
(i) ≤ G 2
M +
1
M2 ∑ i 6=j g>xi(θ ∗)gxj (θ ∗), (38)
where (i) follows from Lemma 3. Consider the expectation of the second term in eq. (38), which is given by
1
M2 ∑ i 6=j E[g>xi(θ ∗)gxj (θ ∗)]. (39)
Without loss of generality, we consider the case when j > i as follows:
E[g>xi(θ ∗)gxj (θ ∗)] = E[E[gxj (θ∗)|Pi]>gxi(θ∗)] ≤ E[ ∥∥E[gxj (θ∗)|Pi]∥∥2 ‖gxi(θ∗)‖2] ≤ GE[
∥∥E[gxj (θ∗)|Pi]∥∥2] = GE[‖E[(Ajθ∗ + bj)|Pi]‖2] ≤ GE[‖E[Aj |Pi]θ∗ + E[bj |Pi]‖2] = GE[‖(E[Aj |Pi]−A)θ∗ + (E[bj |Pi]− b)‖2] ≤ GE[‖(E[Aj |Pi]−A)θ∗‖2 + ‖E[bj |Pi]− b‖2] ≤ GE[‖E[Aj |Pi]−A‖2 ‖θ
∗‖2 + ‖E[bj |Pi]− b‖2] ≤ κG[(1 + γ)Rθ + rmax]ρj−i. (40)
Substituting eq. (40) into eq. (39), we obtain 1
M2 ∑ i 6=j E[g>xi(θ ∗)gxj (θ ∗)] ≤ κG[(1 + γ)Rθ + rmax] M2 ∑ i6=j ρ|i−j|
≤ κG[(1 + γ)Rθ + rmax] M2 (2M dM2 e∑ k=1 ρk)
≤ 2ρκG[(1 + γ)Rθ + rmax] (1− ρ)M . (41)
Then substituting eq. (41) into eq. (38) yields
E[‖gm−k(θ∗)‖22] ≤ 1
M
( G2 + 2ρκG[(1 + γ)Rθ + rmax]
(1− ρ)
) ≤ C4 M , (42)
where C4 = G2 + 2ρκG[(1+γ)Rθ+rmax]
(1−ρ) . Finally, substituting eq. (42) into the accumulated residual variance term in eq. (37), we have
C3 m−1∑ k=0 Ck1E [ ‖gm−k(θ∗)‖22 ] ≤ C3C4 M m−1∑ k=0 Ck1 ≤ C3C4 (1− C1)M . (43)
Step 4: Combining all error terms
Finally, substituting eq. (43) and substituting the values of C2 and C3 into eq. (37), we have E [∥∥∥θ̃m − θ∗∥∥∥2
2
] ≤ Cm1 ∥∥∥θ̃0 − θ∗∥∥∥2 2 + 3C4α+ C2/λA (1− C1)[0.5λA − 3α(1 + γ)2]M ,
which yields the desired result.

heading: F SAMPLE COMPLEXITY OF TD
text: The finite-time convergence rate of vanilla TD under i.i.d. and Markovian sampling has been characterized in Bhandari et al. (2018); Srikant and Ying (2019). However, these studies did not provide the overall computational complexity, i.e., the total number of pseudo-gradient computations to achieve an -accuracy solution. This section provides such an analysis based on their convergence results for completeness.
F.1 TD WITH I.I.D. SAMPLES
Consider the vanilla TD update in Bhandari et al. (2018). Following the steps similar to those in Bhandari et al. (2018) for proving Theorem 2, and letting the constant stepsize α ≤ min{ λA4(1+γ)2 , 2 λA }, we have
E ‖θt − θ∗‖22 ≤ (1− 1
2 λAα)
t ‖θ0 − θ∗‖22 + 2C5α
λA
≤ e− 12λAαt ‖θ0 − θ∗‖22 + 2C5α
λA ,
where 0 < C5 < ∞ is a constant. Let α = min{ λA4(1+γ)2 , 2 λA , λA4C5 }. It can be checked easily that with the total number of iterations at most
t = d 2 λAα
log( 2 ‖θ0 − θ∗‖22 )e = d2 max{4(1 + γ) 2
λ2A , 1, 4C5 λ2A } log(
2 ‖θ0 − θ∗‖22
)e
= O (( 1
λ2A
) log (1 )) , (44)
an -accurate solution can be attained, i.e., E ‖θt − θ∗‖22 ≤ . Since each iteration requires one pseudo-gradient computation, the total number of pseudo-gradient computations is also given by eq. (44).
F.2 TD WITH MARKOVIAN SAMPLES
Consider the vanilla TD update in Bhandari et al. (2018). Following the steps similar to those in Bhandari et al. (2018) for proving Theorem 3, and letting the constant stepsize α ≤ 1λA , we have
E ‖θt − θ∗‖22 ≤ (1− λAα) t ‖θ0 − θ∗‖22 +
C6α λA + C7α log(
1 α )
λA
≤ e−λAαt ‖θ0 − θ∗‖22 + C6α λA + C7α log(
1 α )
λA
where 0 < C6 < ∞ and 0 < C7 < ∞ are constants. Now let α = min{ C8 log(1/C8 ) , 1 λA } where C8 = λA min{ 1C6 , 1 6C7 }, it can be checked easily that with the total number of iterations at most
t = d 2 λAα
log( 3 ‖θ0 − θ∗‖22 )e
= dmax{ 2 min{ 1C6 , 1 6C7 }λ2A
, 2} log ( 1
C8
) log( 3 ‖θ0 − θ∗‖22 )e
= O (( 1
λ2A
) log2 (1 )) , (45)
an -accurate solution can be attained, i.e., E ‖θt − θ∗‖22 ≤ . Since each iteration requires one pseudo-gradient computation, the total number of pseudo-gradient computations is also given by eq. (45).

Review summaries:
The paper is on temporal difference learning, specifically variance reduction of it. As per the claims of the paper, a previous method from (Korda and La, 2015) had technical errors, which the paper corrects and provides a better analysis of variance reduction. In the end, the paper focuses on the variance of the gradient estimator in temporal difference learning, and on analysis of the bias error. The final method is validated on two experiments with iid and Markovian sampling.

In this paper, the authors study the variance reduced TD (VRTD) algorithm, by Korda and Prashanth (2015)  (KP15), for policy evaluation in RL. They first highlight technical errors in the analysis of KP15, and then provide new convergence analysis for this algorithm. The new analysis is based on a new technique to bound the bias of the VRTD gradient estimator, and shows the advantage of VRTD over vanilla TD (the analyses by Bhandari et al. 2018 and Srikant and Ying 2019), both in terms of variance and bias, that are reduced by increasing the batch size. The authors show that while the variance and bias of vanilla TD are both of O(\alpha) (where \alpha is the step-size), they are of O(\alpha/M) and O(1/\sqrt{M}) (where M is the batch size) in VRTD. This shows that a good convergence is possible for VRTD without reducing the step-size \alpha, that causes slow convergence, by increasing the batch size M. In the middle of their analysis, the authors propose a slight modification of VRTD for the case that the samples are obtained iid from the stationary distribution of the evaluating policy, and provide its analysis. Finally, the authors provide simple experiments to support their theoretical findings.