Cover Page
On the Variance of Temporal Difference Learning and its Reduction Using Control Variates Hsiao-Ru Pan, Bernhard Schölkopf Keywords: temporal difference learning, Monte Carlo estimation, advantage function, control variate
arXiv:2606.20357v1 [cs.LG] 18 Jun 2026
Summary We analyze the variance of temporal difference (TD) learning using the phased setting with tabular representation, and show that one of the mechanisms behind its ability to reduce variance is by effectively aggregating over a larger number of independent trajectories. Based on this insight, we demonstrate that (1) the variance of TD is asymptotically bounded from above by Monte Carlo (MC) estimators, and (2) shorter horizon updates incurs less variance for a fixed number of samples. Beyond TD, we show that Direct Advantage Estimation (DAE), a method for estimating the advantage function, can be seen as a type of regression-adjusted control variate, which achieves a tighter bound on the variance compared to TD in the largesample limit. Finally, we numerically illustrate the behaviors of these estimators with carefully designed environments.
Contribution(s) 1. We prove that, in the phased (synchronous) and the IID (asynchronous) settings, the asymptotic variance of temporal difference (TD) learning is upper-bounded by that of Monte Carlo methods, and that one mechanism underlying TD’s variance reduction is its effective aggregation over a larger pool of trajectories. Context: TD learning is widely assumed to reduce variance via bootstrapping (Sutton et al., 1998; Szepesvari, 2010; Dann et al., 2014). While previous works have shown that bootstrapping can be seen as a form of bias-variance tradeoff for multi-step TD (Kearns & Singh, 2000), and that batch TD is more efficient than MC (Grunewalder et al., 2007; Cheikhi & Russo, 2023), the variance comparison is less clear in the online regime. We address this gap by analyzing TD’s variance reduction in the phased (synchronous) setting. 2. We show that (1) the advantage function can be used as control variates for policy evaluation, and (2) Direct Advantage Estimation (DAE) (Pan et al., 2022) is a form of control variate regression that simultaneously estimates both the value and advantage functions. Context: While Pan et al. (2022) demonstrated empirically that DAE improves deep RL performance, its theoretical justification has remained limited to convergence of the estimator. We provide a complementary perspective by showing that DAE also reduces variance by analyzing it through the lens of control variates in policy evaluation, thereby offering a theoretical explanation for its empirical performance gains.
On the Variance of Temporal Difference Learning and its Reduction Using Control Variates Hsiao-Ru Pan1 , Bernhard Schölkopf1,2,3 {hpan,bs}@tuebingen.mpg.de 1
Max Planck Institute for Intelligent Systems, Tübingen ELLIS Institute Tübingen 3 ETH Zürich 2
Abstract We analyze the variance of temporal difference (TD) learning using the phased setting with tabular representation, and show that one of the mechanisms behind its ability to reduce variance is by effectively aggregating over a larger number of independent trajectories. Based on this insight, we demonstrate that (1) the variance of TD is asymptotically bounded from above by Monte Carlo (MC) estimators, and (2) shorter horizon updates incurs less variance for a fixed number of samples. Beyond TD, we show that Direct Advantage Estimation (DAE), a method for estimating the advantage function, can be seen as a type of regression-adjusted control variate, which achieves a tighter bound on the variance compared to TD in the large-sample limit. Finally, we numerically illustrate the behaviors of these estimators with carefully designed environments.
1
Introduction
Policy evaluation is a central problem in reinforcement learning (RL) (Sutton et al., 1998). Among various methods, temporal difference (TD) learning (Sutton, 1988) stands out as a cornerstone technique for this class of problems. Traditional methods like Monte Carlo (MC) methods estimate the expected return by averaging returns of sample trajectories, typically resulting in unbiased but high variance estimates. In contrast, TD learning updates value estimates iteratively through bootstrapping (i.e., estimate based on previous estimates), thereby avoiding the need for full trajectories and tends to exhibit lower variance. Kearns & Singh (2000) showed that bootstrapping can be seen as a form of bias-variance tradeoff, where high variance estimates from sample trajectories are replaced with low variance biased bootstrapped values. While this intuition largely holds true, it is not difficult to see that the full story is more complicated. For example, if we initialize the value estimates with the ground truth values and update them using TD, then we essentially inject variance (e.g., stochastic rewards) into the estimates, and it is not immediately clear how this variance affects bootstrapping asymptotically. This shows that the way bootstrapping reduces variance may be more nuanced than simply "replace a high variance trajectory with a biased estimate". In the present work, we analyze the variance of multi-step TD in the phased (synchronous) setting (Kearns & Singh, 1998), which abstracts away some of the complexities due to stochastic approximations. Beyond TD, we also draw connections to Direct Advantage Estimation (DAE) (Pan et al., 2022), a method that simultaneously estimates the value function and the advantage function, and show that DAE can be seen as a type of control variate regression that enjoys a tighter bound on the variance compared to TD. To summarize, we show that: • The asymptotic variance of multi-step TD is bounded above by that of MC, and one way TD reduces variance is by effectively averaging over a larger pool of trajectories. 1
• The advantage function can be seen as control variates for the MC estimator, which reduces its variance, and DAE is a type of regression-adjusted control variate for TD. Finally, we construct examples to illustrate the asymptotic behaviors of these estimators.
2
Background
We consider a discounted Markov Decision Process (MDP) (Puterman, 2014) (S, A, p, r, γ) with finite state space S, finite action space A, transition probability p(s′ |s, a), reward function r(s, a), and discount factor γ ∈ [0, 1). For simplicity, we assume the reward function is deterministic unless otherwise stated and denote r(st , at ) by rt when the context is clear. A policy π(·|s) is a function that maps states to distributions over A, and we focus P∞ on the case where the policy π isP fixed. The value functions are defined by V π (s) = Eπ [ t=0 γ t rt |s0 =s], and Qπ (s, a) = ∞ Eπ [ t=0 γ t rt |s0 =s, a0 =a], respectively (Eπ indicates that actions are sampled from π). The advantage function is defined by Aπ (s, a) = Qπ (s, a) − V π (s). For the present work, we focus on tabular representations, and use Vπ ∈ R|S| , Aπ ∈ R|S|×|A| to denote the value and the advantage functions, respectively, and use their functional forms for indexing. Learning value functions is at the core of various policy optimization algorithms, and one classical method is TD learning (Sutton, 1988), which updates the estimates through bootstrapping. Specifically, one-step TD estimates the value function by sampling transition tuples (st , at , st+1 ) and updating the values by V (st ) ← V (st ) + α(rt + γV (st+1 ) − V (st )), where α ∈ R is the learning rate. An extension of this is multi-step TD, which samples multiple timesteps before updating the value via V (st ) ← V (st ) + α(rt + · · · + γ k V (st+k ) − V (st )). We denote this update by TD(k)1 . Phased TD We presently consider the phased setting (Kearns & Singh, 1998), which removes some of the complexities of TD leaning due to asynchronous updates and stochastic approximations. In phased TD, value estimates are updated in phases, where each phase consists of two steps: (1) sampling n k-step trajectories for each state: D = {τs,i }s∈S,i∈[1,··· ,n] , where τs,i = (s, ai0 , r0i , si1 , · · · , sik ), and (2) updating the value of every state synchronously by n
V T +1 (s) ←
1X i i r0 + · · · + γ k−1 rk−1 + γ k V T sik , n i=1
(1)
where T denotes the phase. It was shown that phased TD is analogous to TD learning with a fixed learning rate under mild assumptions. Kearns & Singh (2000) used this setting to analyze the biasvariance tradeoff of multi-step TD, and showed that the estimation error satisfies: ! " k−1 ! # n k−1 X 1X X t i T +1 π t max |V (s) − V (s)| := ∆T +1 ≤ γ rt − E γ rt |s0 =s +γ k ∆T , s n i=1 t=0 t=0 |{z} | {z } variance
bias
where the variance term accounts for the stochasticity from the sample rewards, while the bias term accounts for the error from bootstrapping. The authors then derived the following probabilistic bound (assume ∆0 = 1 and |r(·)| ≤ 1): r 1 − γ kT 3 log(k/δ) ∆T ≤ + γ kT , (2) 1−γ n with probability 1 − δ. Intuitively, increasing k reduces the bias from bootstrapping at the cost of increased variances from the rewards. However, this bias-variance decomposition fails to consider the variance from finite samples of sik , and the variance from previous estimates V T . Furthermore, as k → ∞ (MC estimation), this bound becomes vacuous even though the rewards are bounded. 1 This is not to be confused with TD(λ).
2
Control Variate Control variate (Hammersley & Handscomb, 1964; Asmussen & Glynn, 2007) is a technique for reducing variances of MC methods. We briefly review the basics of control variate in the one-dimensional setting, and refer the reader to Owen (2013) for a more general treatment. Pn iid Recall that MC methods estimate µX = E[X] by µ̂X = n1 i=1 Xi , where Xi ∼ X. Given another random variable Y , the control variate, that satisfies Cov(X, Y ) ̸=P0 and has known E[Y ] n (assume E[Y ] = 0, otherwise use Z = Y − E[Y ]). Let µ̂X,λY = n1 i=1 Xi + λYi , where iid
(Xi , Yi ) ∼ (X, Y ) and λ is a tunable constant. This new estimator remains unbiased, but has a 2 ∗ different variance Var(µ̂X,λY ) = Var(µ̂X ) + λn Var(Y ) + 2λ n Cov(X, Y ). If we choose λ = λ = 2
) Cov(X,Y ) ∗ − Cov(X,Y Var(Y ) , then Var(µ̂X,λ Y ) = Var(µ̂X )− nVar(Y ) ≤ Var(µ̂X ), which is never worse than the ) original estimator. In practice, however, λ∗ is usually unknown, and estimated by λ̂ = − Cov(X,Y d ) . Var(Y Interestingly, this estimator arises naturally from solving the following least squares: d
n X (µ̂X,λ̂Y , λ̂) = arg min (θ − xi − λyi )2 . θ,λ
(3)
i=1
This approach is also known as regression-adjusted control variate, and can be readily generalized to more complex settings (e.g., multiple control variates). The error of this estimator is: ∗
µ̂X,λ̂Y − µX = (λ̂ − λ )
n X Yi i=1
|
{z
O(1/n)
n }
+ µ̂X,λ∗ Y − µX . |
√ O(1/ n)
{z
(4)
}
Pn √ As both λ̂ − λ∗ and i=1 Yni approach zero with rates O(1/ n), their product converges to 0 with rate O(1/n). Consequently, the second term dominates the error asymptotically, and µ̂X,λ∗ Y can be seen as a large-sample approximation of µ̂X,λ̂Y . Finally, we note that the estimator is, in general, biased since Yi can be correlated with λ̂ (i.e., E[Yi λ̂] ̸= 0); however, the estimator remains consistent as the errors approach zero for n → ∞. Direct Advantage Estimation Direct Advantage Estimation (DAE) (Pan et al., 2022) is a method for estimating the advantage function directly from sample trajectories. Similar to multi-step TD, DAE can update values by bootstrapping previous estimates. Specifically, DAE estimates the values by iteratively minimizing the following constrained least-squares: L(Â, V̂ ) = Eπ
k−1 X
!2 γ t (rt − Ât ) + γ k Vtar (sk ) − V̂ (s0 )
s.t.
X
π(a|s)Â(s, a) = 0, (5)
a
t=0
where Ât = Â(st , at ), and Vtar is the bootstrapping target. It was shown that the minimizer of this objective can be viewed as multi-step estimates of the advantage function and the value function. In practice, the expectation is replaced with an average over sample trajectories. We note that, if we force  ≡ 0 and only optimize with respect to V̂ , then DAE reduces to multi-step TD. DAE demonstrated strong empirical performance in the deep RL setting; however, it remains unclear whether estimating the value function this way is beneficial compared to classical approaches.
3
Variance of Monte Carlo (MC) and TD Learning
We begin by showing that the variance of MC can be broken down into segments. This will become useful when compared to multi-step TD, whose variance follows a similar form. For clarity, we summarize the notations used in subsequent analyses in Table 1. 3
Table 1: Variables names and their definitions. i denotes the ith trajectory. k and n are constants unless otherwise stated. Note that P is a (row) vector of dimension |S|. Variable G D R P
Def. Pn P∞ t i 1 i=1 t=0 γ rt n (s, ai0 , r0i , si1 , · · · , sik ) s∈S,i=1,...,n Pn Pk−1 t i 1 i=1 n Pn t=0 γ rt Ps′ = n1 i=1 I(sik = s′ )
Description Average of sample returns k-step partial trajectories Average of sample k-step rewards Empirical k-step transition distribution
Variance of MC Recall that MC estimates values via VMC (s) = G. The following lemma shows that its variance is equal to the sum of the variances of partial trajectories. P∞ Lemma 1. Var(VMC (s)) = Var(G) = m=0 γ 2km E Var(R + γ k PVπ |s0 =skm ) s0 =s , This is a multi-step (k ≥ 1) multi-sample (n ≥ 1) generalization of the variance recursion derived by Sobel (1982), see Appendix 8.1 for a proof. As we will see, the variance of TD shares a similar structure, except Vπ is replaced by bootstrapping targets. Variance of Multi-Step TD
We first rewrite the phased update (Equation 1) into: T −1 T VTD(k) (s) = R + γ k PVTD(k) .
(6)
0 T Without loss of generality, we set VTD(k) ≡ 0. Note that VTD(k) is now a random vector for T > 0 T as the update involves random variables R and P, and its expectation, denoted V̄TD(k) , is equal to: " k−1 # h i X T −1 T k t k T −1 V̄TD(k) (s) = E [R] + γ E [P] E VTD(k) = E γ rt + γ V̄TD(k) (sk ) s0 =s , t=0
and we recover the k-step Bellman update. We emphasize that, while this expectation converges to T Vπ in the limit (T → ∞), VTD(k) in general does not. This is due to the finite-sample variance, which causes the estimates to oscillate around Vπ . We now analyze this variance and show that, similar to Lemma 1, the upper bound of the variance of TD also satisfies a recurrence relation. Lemma 2. i h T −1 T −1 T s0 =s + γ 2k Esk Var VTD(k) Var(VTD(k) (s)) ≤ Var R + γ k PV̄TD(k) (sk ) s0 =s . (7) Proof. T −1 T Var(VTD(k) (s)) = Var R + γ k PVTD(k) s0 =s i h i h T −1 T −1 =Var E R + γ k PVTD(k) D s0 =s + E Var R + γ k PVTD(k) D s0 =s h i T −1 T −1 =Var R + γ k PV̄TD(k) s0 =s + γ 2k E Var PVTD(k) D s0 =s h i T −1 T −1 ≤Var R + γ k PV̄TD(k) s0 =s + γ 2k E PVar VTD(k) s0 =s h i T −1 T −1 =Var R + γ k PV̄TD(k) s0 =s + γ 2k Esk Var VTD(k) (sk ) s0 =s .
(8) (9) (10) (11) (12)
If we expand this recursion, then T Var(VTD(k) (s)) ≤
T −1 X
h i T −1−m γ 2km E Var R + γ k PV̄TD(k) s0 =skm s0 =s ,
(13)
m=0
and we arrive at a similar expression to the MC estimator (Lemma 1). In fact, we get, in the limit: 4
s1k
s2k s0
sik
s3k
s0
(a) Independent bootstrapping states
(b) Correlated bootstrapping states
Figure 1: Illustration of how TD can reduce variance. Solid and dashed arrows represent trajectories collected in the current phase and the previous phase, respectively. (a) sik are independent, and averaging their values effectively aggregates over a larger pool of trajectories. (b) sik are all the same, and averaging their values provides no variance reduction, as the effective number of trajectories remains the same.
T Theorem 1. limT →∞ Var(VTD(k) (s)) ≤ Var(VMC (s)).
See Appendix 8.2 for a proof and a discussion of how this can be extended to the asynchronous IID setting. This shows that, asymptotically, TD is no worse than MC, independent of k; however, it also suggests that TD, in the worst case, can suffer from the same variance as MC. When hold? Let usexamine the only inequality used in the derivation, namely, does the equality T −1 T −1 Var PVTD(k) D ≤ PVar VTD(k) . Recall that the left-hand side is simply the variance of the P T −1 (sik ) D , and this inequality becomes average of the bootstrap values, namely Var n1 i VTD(k) an equality when the value estimates have correlation 1 (e.g., sik are the same for all i). This suggests that one way bootstrapping reduces variance is by effectively aggregating over diverse states with weak dependencies between their value estimates, as illustrated in Figure 1. In Section 5, we also construct toy examples to illustrate this effect. A remark on sample efficiency It is not immediately clear whether the comparison in Theorem 1 is fair in terms of sample efficiency as the result is asymptotic. Here, we compare two settings with the same number of samples: (1) T phases with k-step updates, and (2) 1 phase with kT -step update. Both settings require nkT |S| number of samples, and share the same expected value given the same initialization. However, their variances satisfy the following ordering: τ +T τ +1 τ τ Corollary 1. For τ ≥ 0, we have Var(VTD(k) (s)|VTD(k) =V) ≤ Var(VTD(kT ) (s)|VTD(kT ) =V).
See Appendix 8.3 for a proof. This shows that short horizon updates incur less variance than a single long horizon update given the same number of interactions and bootstrapping targets.
4
Control Variate and the Advantage Function
We now shift our focus to the second main contribution of the present work, namely, how to use the advantage function to reduce the variance of value estimations. We begin by considering the π-centered functionP class Fπ = {f |Eπ [f (s, a)|s] =P0 ∀s}, which ∞ ∞ has the following property (given f ∈ Fπ ): Eπ [ t=0 γ t (rt − f (st , at ))] = Eπ [ t=0 γ t rt ]. In other words, introducing f does not bias the MC estimate, and f can be seen as a control variate. A natural question is, then, what would be the optimal choice of f ∗ that minimizes the variance, that P∞ ∗ is, f = arg minf ∈Fπ Var ( t=0 γ t (rt − f (st , at ))). Pan et al. (2022) proved that the advantage function Aπ is the unique minimizer under mild assumptions, and used it as a way to estimate Aπ . 5
Now, if Aπ is known, then we can construct a new estimator via: ! n ∞ 1X X t i π i i VMC-A (s) = γ rt − A (st , at ) . n i=1 t=0
(14)
One may wonder to what extent can this control variate reduce the variance. To answer this, we use the return decomposition by Pan & Schölkopf (2024): ∞ X
γ t rt = V π (s0 ) +
t=0
∞ X
γ t (Aπ (st , at ) + B π (st , at , st+1 )),
(15)
t=0
where B π (st , at , rt , st+1 ) = rt + γV π (st+1 ) − E[r + γV π (s′ )|st , at ] quantifies how much of the return is caused by stochastic transitions2 . If B π ≡ 0 (e.g., the environment is deterministic), then: ∞ X
γ t rt = V π (s0 ) +
t=0
∞ X
γ t Aπ (st , at ).
(16)
t=0
Combine this with Equation 14, we have Var(VMC-A (s)) = 0, meaning that Aπ can fully explain away the variance of MC in this case. For more general cases, B π is required to account for the variance caused by stochastic transitions, and it remains open whether B π can be easily estimated in model-free settings. As such, we focus on the advantage function in the present work. In practice, Aπ is rarely known a priori, so we have to estimate both V π and Aπ simultaneously. Next, we show that DAE (Pan et al., 2022) can be seen as a type of regression-adjusted control variate, which achieves this. 4.1
Direct Advantage Estimation and Control Variate Regression
DAE estimates V π and Aπ by solving a constrained least-square problem (Equation 5). For the present work, we focus only on the value estimates and treat  as nuisance parameters. This allows us to remove the constraint and reformulate the objective into: !2 k−1 X T LT (Â, V̂ ) = Eπ γ t rt − Ăt + γ k VDAE(k) (sk ) − V̂ (s) , t=0
P P where Ăt = Â(st , at )− a π(a|st )Â(st , a) (this parametrization ensures that a π(a|s)Ă(s, a) = 0). Under this formulation, the minimizer of  may no longer be unique, but the minimizer of V̂ remains unchanged. Let us now consider DAE in the phased setting, where the objective becomes: n k−1 X X i=1
γ
t
rti − Ăit
!2 +γ
k
T VDAE(k) (sik ) − V̂ (s)
.
t=0
T +1 Let (VDAE(k) , AT +1 ) be a minimizer of this objective, and M ∈ R|S|×|A| be the following:
M(s,a) =
n k−1 1 XX t γ I(sit =s, ait =a) − π(a|s)I(sit =s) , n i=1 t=0
which compares the empirical occupancy measure to its expectation over the actions, then T +1 T VDAE(k) (s) = R + γ k PVDAE(k) − MAT +1 .
(17)
Note that, this update rule differs from the TD update (cf. Equation 6) only by the term MAT +1 . Comparing this to Equation 3, we see that DAE is a case of regression-adjusted control variate, 2 We use a slightly more general definition to incorporate stochastic rewards
6
1
.. .
2
.. .
···
8
Table 2: Parameters of the experiment.
.. .
Param. |A| n k pr ps
Figure 2: Chain MDP with S = {1, 2, ..., 8}, a (indepenA = {1, ..., |A|}, and r(s, a) = (−1) 4 dent of s). Agents return to state 1 after state 8.
Description action space size number of sample trajectories backup length probability of reward masking probability of sticky transition
where M is the control variate with E[M] = 0, and AT is the corresponding coefficients. AddiT tionally, recall that this regression estimator (VDAE(k) ) behaves similarly to the one with optimal T control variate coefficients (VDAE∗ (k) ) under the large-sample approximation (cf. Section 2). Here, T we focus on VDAE ∗ and leave it for future work to analyze the small-sample bias due to control 3 variates. The following result shows the asymptotic behavior of DAE: T Theorem 2. lim supT →∞ Var(VDAE ∗ (k) (s)) ≤ Var(VMC-A (s))
See Appendix 8.4 for a proof. Compared to Theorem 1, this shows that DAE enjoys a tighter upper bound on the variance by using control variates (recall that Var(VMC-A (s)) ≤ Var(VMC (s))). Finally, we note that since TD is a special case of Equation 17 with  ≡ 0, it follows that T +1 T +1 T T Var(VDAE ∗ (k) (s)|VDAE∗ (k) =V) ≤ Var(VTD(k) (s)|VTD(k) =V).
5
Empirical Illustration
In this section, we illustrate the behaviors of different estimators through experiments based on variants of the chain environment shown in Figure 2. Despite its simplicity, the environment is sufficiently expressive to illustrate various properties of the estimators analyzed in the present study, such as how the variances depend on the bootstrapping states. All experiments are based on the phased setting, where values are updated synchronously at the end of each phase. We fix the number of phases at 2500, which we found sufficient convergence, the policy π to be uniform, the discount factor at γ = 0.99, and consider two types of stochasticity: 1. Reward masking: Rewards are masked out with probability pr (i.e., r = 0 with probability pr ). 2. Sticky transition: With probability ps , the agent stays in the current state instead of advancing to the next state (i.e., st+1 = st with probability ps ). For simplicity, we consider cases where |A| is even, such that V π ≡ 0. The variances of the MC and the MC-A estimators are equal to Var(VMC (s)) =
(1 − pr ) , and Var(VMC-A (s)) = pr Var(VMC (s)), 16(1 − γ 2 )n
(18)
respectively, and will be used as baselines for comparing the estimators. We note that the variances do not depend on ps due to the symmetric nature of the states. Each run (configuration) is repeated for 1000 different random seeds to ensure statistical significance. We compare the mean squared error between the true value function V π and the estimated value (MSE(V̂ , V π ), MSEπ for short) , which reflects how quickly the estimator converges to V π . We also measure the error between the estimated value and the expected Bellman update at each phase (MSE(V̂ , VBE ), MSEBE for short). For TD, MSEBE corresponds to the variance of the upT T date, since the expected TD update follows the Bellman update (i.e., E[VTD ] = VBE ). For DAE, it additionally captures the small-sample bias introduced by the regression estimator. 3 The large-sample approximation is common in the literature (Lin, 2013; Owen, 2013; Davidson & MacKinnon, 1992), since the small-sample bias from imperfect control variate coefficients is often much smaller than the variance.
7
1.0
0.4
0.5
0.2
0.0
100
101
102
MSE (V, VBE)
MSE (V, V )
1.0 0.5 0.0
100
101
102
103
1.0 0.5 0.0
0.0
103
100
101
102 Phase
103
0.3 0.2 0.1 0.0 0.3 0.2 0.1 0.0
100
101
102
103
100
101
102
103
100
101
102 Phase
103
k-step 64 16 4 Method TD DAE Bellman MC MC-A
Figure 3: The Deterministic case (top), the Sticky Transition case (middle), and the Reward Masking case (bottom). Lines and shadings represent (mean ± 3 standard error). The Deterministic Case (|A| = 2, n = 8, k ∈ {4, 16, 64}, ps = 0, pr = 0) We first examine the effect of k in the simplest setting with a deterministic environment. Since the transitions are fully deterministic and independent of actions, the bootstrapping states are identical for a given starting state. Our analysis (Section 3) suggests that bootstrapping loses its ability to reduce variance when the value estimates of the bootstrapping states are fully correlated, and TD would behave like MC asymptotically in this case. Indeed, Figure 3 (top) shows that TD learning, independent of k, approaches the same MSEπ as MC. Note that this cannot be explained by the bound given by Equation 2, which predicts that the asymptotic error would grow as k increases. From MSEBE , we also see the effect of variance injection, where the variances increase with respect to the number of updates. On the other hand, we see DAE achieving almost zero MSEBE , demonstrating the effectiveness of the variance reduction. The Sticky Transition Case (|A| = 2, n = 8, k ∈ {4, 16, 64}, ps = 0.25, pr = 0) We now examine how stochastic transitions affect the estimators. Figure 3 (middle) shows that, counterintuitively, stochasticity reduces the variance of TD. Furthermore, the learning curves now follow the common bias-variance tradeoff intuition of multi-step learning (i.e., larger k learns faster but leads to higher variance and vice versa). We emphasize that this setting effectively differs from the deterministic case only in how the bootstrapping states are sampled. As such, the variance reduction can only be explained by TD learning’s ability to effectively aggregate value estimates over more diverse states (cf. Figure 1). DAE converges with rates similar to Bellman iterations again, since B π ≡ 0 does not depend on ps . So far, both cases have almost full coverage of the state-action space from the data in each phase, and have variances that can be fully explained by the advantage function (cf. Equation 16). This allowed DAE to converge almost as fast as Bellman iterations. To see when this breaks down, we next consider the stochastic reward (B π ̸≡ 0) and the large |A| settings. The Reward Masking Case (|A| = 2, n = 8, k ∈ {4, 16, 64}, ps = 0, pr = 0.2) Similar to the deterministic case, Figure 3 (bottom) shows that, asymptotically, TD performs similar to MC irrespective of k, as the bootstrapping states are deterministic. On the other hand, the variance of 8
MC-A is no longer zero since the variance of the rewards cannot be explained by the actions. This causes the MSEπ of DAE to converge slightly above Bellman iterations asymptotically. The Low Coverage Case (|A| ∈ {4, 16, 64}, n = 8, k = 16, ps = 0, pr = 0) Finally, we consider the low coverage case, where the least squares objective become underdetermined for DAE. In the context of RL, this means that most of the actions are not sampled, and the advantage estimates become unreliable. We follow the common practice of choosing the minimum Euclidean-norm solution. Figure 4 shows how the small-sample bias affects learning in early phases (phase < 10) and its relationship to the size of the action space. Interestingly, we find that even when |A| ≥ 16, where DAE converges to suboptimal solutions compared to the full coverage cases, DAE still converges to a lower MSEπ than TD, suggesting that the advantage estimates may help reduce variance even when they are poorly estimated. 4
0.4 MSE (V, VBE)
MSE (V, V )
1.0 0.5 0.0
100
101
102 Phase
4
| | 16
64
0.2 0.0
103 TD
DAE
100
101 Method Bellman
102 Phase MC
103 MC-A
Figure 4: The low coverage case. Lines and shadings represent (mean ± 3 standard error). Note that TD (black, solid) does not depend on |A|.
6
Related Work
The bias-variance tradeoff between MC and TD has been discussed in classical texts and surveys (Sutton et al., 1998; Szepesvari, 2010; Dann et al., 2014); however, a rigorous analysis remained desirable. Kearns & Singh (2000) analyzed the error bounds of multi-step TD learning in the phased setting (Kearns & Singh, 1998), which also serves as the basis of the present work. In the batch setting, Grunewalder et al. (2007) showed that LSTD (Bradtke & Barto, 1996) is statistically more efficient than MC when the Markovian structure of the environment can provide additional information. More recently, Cheikhi & Russo (2023) derived a more precise statistical relationship between batch TD and MC by analyzing trajectory pooling. However, the batch setting abstracts away the iterative nature of TD, and their implications in the online setting remain unclear. Similar problems have also been studied in function approximation settings (Tsitsiklis & Van Roy, 1996; Dalal et al., 2018; Bhandari et al., 2018), where convergence and finite-sample (time) error bounds were explored. The advantage function (Baird, 1995) is commonly used as control variates for policy gradient methods (Sutton et al., 1999; Greensmith et al., 2004). The present work demonstrates that the advantage function can also be used as control variates for policy evaluation, and shows that DAE (Pan et al., 2022) can be seen as regression-adjusted control variates.
7
Discussion
We analyzed the behaviors of MC, TD, and DAE, and revealed one mechanism behind the variance reduction property of bootstrapping, namely, the ability to aggregate over a larger number of independent trajectories. Furthermore, we established a connection between DAE and control variate regression, demonstrating how it can further reduce the variance of TD learning. At its core, DAE 4 We find DAE converging to the true value function again using iterative solvers (see Appendix 9.1).
9
exploits our knowledge of the policy to reduce variances and an interesting future direction would be to explore other types of control variates for policy evaluation. Finally, we note some limitations: (1) In the present work, we considered only the phased and the IID settings, and it remains unclear to what extent the results generalize to more realistic sampling schemes (e.g., Markovian sampling). (2) We focused mainly on the variance reduction property of DAE, but it should be noted that DAE also incurs additional space (store M) and time (solve least-squares) complexities. (3) The theoretical results for DAE only hold up to the large-sample approximation case where the control variate coefficients are optimal, and while empirical results seem to suggest that the variance reduction is beneficial even in the small-sample regime, a more rigorous analysis remains desirable. Acknowledgments Hsiao-Ru Pan thanks Claire Vernade and Onno Eberhard for the fruitful discussions. The authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Hsiao-Ru Pan.
References E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, and D. Sorensen. LAPACK Users’ Guide. Society for Industrial and Applied Mathematics, Philadelphia, PA, third edition, 1999. ISBN 0-89871-447-8. Søren Asmussen and Peter W Glynn. Stochastic simulation: algorithms and analysis, volume 57. Springer, 2007. Leemon Baird. Residual algorithms: Reinforcement learning with function approximation. In Machine Learning Proceedings 1995, pp. 30–37. Elsevier, 1995. Jalaj Bhandari, Daniel Russo, and Raghav Singal. A finite time analysis of temporal difference learning with linear function approximation. In Conference on learning theory, pp. 1691–1692. PMLR, 2018. James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: composable transformations of Python+NumPy programs, 2018. URL http: //github.com/jax-ml/jax. Steven J Bradtke and Andrew G Barto. Linear least-squares algorithms for temporal difference learning. Machine learning, 22(1):33–57, 1996. Robert M Burton and Uwe Rösler. An l2 convergence theorem for random affine mappings. Journal of applied probability, 32(1):183–192, 1995. David Cheikhi and Daniel Russo. On the statistical benefits of temporal difference learning. In International Conference on Machine Learning, pp. 4269–4293. PMLR, 2023. Gal Dalal, Balázs Szörényi, Gugan Thoppe, and Shie Mannor. Finite sample analyses for td (0) with function approximation. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 32, 2018. Christoph Dann, Gerhard Neumann, and Jan Peters. Policy evaluation with temporal differences: A survey and comparison. The Journal of Machine Learning Research, 15(1):809–883, 2014. Russell Davidson and James G MacKinnon. Regression-based methods for using control variates in monte carlo experiments. Journal of Econometrics, 54(1-3):203–222, 1992. 10
Evan Greensmith, Peter L Bartlett, and Jonathan Baxter. Variance reduction techniques for gradient estimates in reinforcement learning. Journal of Machine Learning Research, 5(9), 2004. Steffen Grunewalder, Sepp Hochreiter, and Klaus Obermayer. Optimality of lstd and its relation to mc. In 2007 International Joint Conference on Neural Networks, pp. 338–343. IEEE, 2007. J. M. Hammersley and D. C. Handscomb. Monte Carlo Methods. Springer Netherlands, 1964. ISBN 9789400958197. DOI: 10.1007/978-94-009-5819-7. URL http://dx.doi.org/10. 1007/978-94-009-5819-7. Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, Robert Kern, Matti Picus, Stephan Hoyer, Marten H. van Kerkwijk, Matthew Brett, Allan Haldane, Jaime Fernández del Río, Mark Wiebe, Pearu Peterson, Pierre Gérard-Marchant, Kevin Sheppard, Tyler Reddy, Warren Weckesser, Hameer Abbasi, Christoph Gohlke, and Travis E. Oliphant. Array programming with NumPy. Nature, 585(7825):357–362, September 2020. DOI: 10.1038/ s41586-020-2649-2. URL https://doi.org/10.1038/s41586-020-2649-2. Michael Kearns and Satinder Singh. Finite-sample convergence rates for q-learning and indirect algorithms. Advances in neural information processing systems, 11, 1998. Michael J Kearns and Satinder Singh. Bias-variance error bounds for temporal difference updates. In COLT, pp. 142–147, 2000. Winston Lin. Agnostic notes on regression adjustments to experimental data: Reexamining freedman’s critique. The Annals of Applied Statistics, pp. 295–318, 2013. Art B. Owen. Monte Carlo theory, methods and examples. https://artowen.su.domains/ mc/, 2013. Christopher C Paige and Michael A Saunders. Lsqr: An algorithm for sparse linear equations and sparse least squares. ACM Transactions on Mathematical Software (TOMS), 8(1):43–71, 1982. Hsiao-Ru Pan and Bernhard Schölkopf. Skill or luck? return decomposition via advantage functions. arXiv preprint arXiv:2402.12874, 2024. Hsiao-Ru Pan, Nico Gürtler, Alexander Neitz, and Bernhard Schölkopf. Direct advantage estimation. Advances in Neural Information Processing Systems, 35:11869–11880, 2022. Martin L Puterman. Markov decision processes: discrete stochastic dynamic programming. John Wiley & Sons, 2014. Matthew J Sobel. The variance of discounted markov decision processes. Journal of Applied Probability, 19(4):794–802, 1982. Richard S Sutton. Learning to predict by the methods of temporal differences. Machine learning, 3 (1):9–44, 1988. Richard S Sutton, Andrew G Barto, et al. Introduction to reinforcement learning. 1998. Richard S Sutton, David McAllester, Satinder Singh, and Yishay Mansour. Policy gradient methods for reinforcement learning with function approximation. Advances in neural information processing systems, 12, 1999. Csaba Szepesvari. Algorithms for Reinforcement Learning. Morgan & Claypool Publishers, 2010. John Tsitsiklis and Benjamin Van Roy. Analysis of temporal-diffference learning with function approximation. Advances in neural information processing systems, 9, 1996. 11
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern, Eric Larson, C J Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261–272, 2020. DOI: 10.1038/s41592-019-0686-2.
12
Supplementary Materials The following content was not necessarily subject to peer review.
8
Proofs
8.1
Proof of Lemma 1
Lemma 1 Var(VMC (s)) = Var(G) =
∞ X
γ 2km E Var(R + γ k PVπ |s0 =skm ) s0 =s .
m=0
Proof. Note that, when conditioned on D, both R and P become constants. Consequently, we have: Var(VMC (s)) = Var ( G| s0 =s) = Var ( E [ G| D, s0 =s]| s0 =s) + E [ Var ( G| D, s0 =s)| s0 =s] # ! " ∞ n X X 1 γ t rti D s0 =s = Var R + γ k PVπ s0 =s + E Var n i=1 t=k # " n X 1 k π k i = Var R + γ PV s0 =s + E Var γ VMC (sk ) s0 =s n i=1 = Var R + γ k PVπ s0 =s + γ 2k E [Var (VMC (sk )) |s0 =s] . Since Var(VMC (·)) is bounded, Lemma 1 follows from expanding this recursion. 8.2
Proof of Theorem 1
Theorem 1 T lim Var(VTD(k) (s)) ≤ Var(VMC (s)).
T →∞
T Proof. First, we show that limT →∞ Var(VTD(k) (s)) converges. Note that T −1 T VTD(k) (s) = R + γ k PVTD(k) ,
is a special case of random affine iterated system of the form: Xt = Mt Xt−1 + Nt ,
(19)
where (Mt , Nt ) are IID random variables. Furthermore, since γ k P is a contraction respect to || · ||∞ T and R is bounded, we know that the random vector VTD(k) converges in distribution with respect to the Wassterstein distance W∞ (Burton & Rösler, 1995). Consequently, all finite moments of T VTD(k) also converges as T → ∞. By Lemma 2, we have: T Var(VTD(k) (s)) ≤
T −1 X
h i T −1−m γ 2km E Var R + γ k PV̄TD(k) s0 =skm s0 =s .
m=0
It is enough to show that the hsummation on the right hand side converges i to Var(VMC (s)) as T −1−m T → ∞. Let xT −m,m = E Var R + γ k PV̄TD(k) s0 =skm the following limit T −1 X lim γ 2km xT −m,m . T →∞
m=0
13
s0 =s , we are interested in
Since x∞,m := limT →∞ xT −m,m = E Var R + γ k PVπ s0 =skm s0 =s , there exists N ∈ N such that if T − m > N then |xT −m,m − x∞,m | < ϵ. In addition, since R, P and Vπ are all bounded, there exists M ∈ R such that |x∞,m | < M and |xT −m,m | < M for all m, T − m ∈ Z+ . Consequently, T −1 X
γ 2km xT −m,m −
m=0
≤
T −1 X
∞ X
γ 2km x∞,m
m=0
γ 2km (xT −m,m − x∞,m ) +
m=0
≤
T −n−1 X
∞ X m=T
γ 2km (xT −m,m − x∞,m ) +
T −1 X
γ 2km (xT −m,m − x∞,m ) +
m=T −n
m=0 2k(T −n)
≤
γ 2km x∞,m
ϵ 2M γ + 1 − γ 2k 1 − γ 2k
M γ 2kT 1 − γ 2k
2kT
+
Mγ 1 − γ 2k
is arbitrarily small as T → ∞, and T lim Var(VTD(k) (s)) ≤ lim
T →∞
T →∞
T −1 X
γ 2km xT −m,m =
∞ X
γ 2km x∞,m = Var(VMC (s)).
m=0
m=0
Extension to IID Sampling In the IID setting, we sample states independently from a fixed distribution µ (e.g., the occupancy measure) and update their values accordingly. Assuming that µ has full support over the state space S, the above analysis can be extended to the asynchronous setting with minor modifications. Specifically, at each iteration T , we sample a starting state s ∼ µ and update its value using the k-step TD learning rule. Accordingly, D now denotes a collection of k-step trajectories originating from s, rather than samples covering the entire state space S. We are now ready to proof the following theorem: Theorem 3. Given distribution µ over S. We define the IID phased TD iteration as a two-step process: (1) sample a state s ∼ µ, and (2) update its value estimate by ( T −1 R + γ k PVTD(k) s̃ = s T VTD(k) (s̃) = . T −1 VTD (s̃) otherwise T If µ covers S, then limT →∞ Var(VTD(k) (s)) ≤ Var(VMC (s)).
Proof. Under this sampling scheme, the expected value vector evolves according to a slightly difT −1 T −1 T ferent equation V̄TD(k) (s) = µ(s)(E[R] + γ k E[P]V̄TD(k) ) + (1 − µ(s))V̄TD(k) (s). If µ covers S T π (i.e., µ(s) > 0 for all s ∈ S), then V̄TD(k) (s) converges to V . Note that Lemma 2 does not depend on the sampling distribution µ, and continues to hold for the sampled state. However, we now need to account for the variance from whether a state gets updated or not. More precisely, we have T +1 T T Var(VTD(k) (s)) =µ(s)Var(R + γ k PVTD(k) ) + (1 − µ(s))Var(VTD(k) (s)) 2 T T + µ(s)(1 − µ(s)) E[R + γ k PVTD(k) ] − E[VTD(k) (s)] .
We now denote 2 T T ∆T = E[R + γ k PVTD(k) ] − E[VTD(k) (s)] . 14
T Since the expectation of VTD(k) converges to Vπ , we have ∆T → 0. Expanding this recursion gives T +1 Var(VTD(k) (s)) = T T X X T −t µ(s) (1 − µ(s))t Var(R + γ k PVTD(k) )+ (1 − µ(s))t+1 ∆T −t t=0
! .
t=0
T In the limit T → ∞, the first sum converges to limT →∞ Var(R + γ k PVTD(k) ), and the second T −t sum converges to limT →∞ ∆ = 0. The results then follow from Lemma 2.
Remark: The key requirement for this result is that trajectory sampling be independent of the value estimates. This requirement may be violated in the Markovian sampling setting, where the value T −1 estimates can be correlated with the next-state distribution (i.e., Such P and V are dependent). T −1 T −1 dependence invalidates the conditional variance bound Var PVTD(k) D ≤ PVar VTD(k) , which is crucial for the analysis. 8.3
Proof of Corollary 1
τ τ Corollary 1 If VTD(k) =VTD(kT ) =V, then τ +T τ +1 τ τ Var(VTD(k) (s)|VTD(k) =V) ≤ Var(VTD(kT ) (s)|VTD(kT ) =V).
Proof. Without loss of generality, we let τ = 0, since V0 is arbitrary. We proof by induction on T . For T = 1, equality holds. For T > 1, assume the inequality holds for T − 1. By Lemma 2, we have: h i T −1 T −1 T s0 =s + γ 2k Esk Var VTD(k) (sk ) s0 =s . Var(VTD(k) (s)) ≤ Var R + γ k PV̄TD(k) Similar to Lemma 1, we can break down the variance of long trajectories into segments, which leads to: ! ! n kT −1 1X X t i 1 kT 0 i Var(VTD(kT ) (s)) = Var γ rt + γ VTD(kT ) (skT ) s0 =s n i=1 t=0 h i 1 2k 1 = Var R + γ k PV̄TD(k(T s =s + γ E Var(V (s ))|s =s 0 s k 0 k −1)) TD(k(T −1)) T −1 1 Recall that the expected value estimates follow the Bellman update, so V̄TD(k) = V̄TD(k(T −1)) . Since the inequality holds for T − 1, we have: h i h i T −1 1 (s ) s =s , Esk Var VTD(k) (sk ) s0 =s ≤ Esk Var VTD(k(T 0 −1)) k
which concludes the proof. 8.4
Proof of Theorem 2
Theorem 2 T lim sup Var(VDAE ∗ (k) (s)) ≤ Var(VMC-A (s)). T →∞
T ∗T Proof. Note that VDAE . Conse∗ (k) (s) is DAE with the optimal control variate coefficients A quently, we must have: T −1 T Var(VDAE R + γ k PVDAE(k) − MA∗T s0 =s (20) ∗ (k) (s)) = Var T −1 ≤ Var R + γ k PVDAE(k) − MAπ s0 =s , (21)
15
where Aπ is the true advantage function. Since Aπ is a constant vector, and M is uniformly bounded, we can rewrite the right hand side of this inequality by T −1 s0 = s , (22) Var R′ + γ k PVDAE(k) where R′ = R−MAπ . By Corollay 1, we know that this variance is bounded above by the variance of the MC estimator with this new reward function, which is precisely the variance of VMC-A . Finally, we make a remark about this corollary. Since the advantage estimate AT now also depends on the bootstrapping value function, it is not clear whether the update remains a contraction, or whether higher moments (e.g., variance) also converge. As such, we only prove the supremum limit is upper bounded by VMC-A .
9
Experimental Details
Algorithm 1 shows the pseudocode. All experiments are based on Python with least-square solvers implemented by NumPy (Harris et al., 2020), SciPy (Virtanen et al., 2020) or JAX (Bradbury et al., 2018). A single run (1 seed, 2500 phases) takes less than a minute on commercial CPUs, except for the large |A| experiment, where we leveraged GPUs (Nvidia A100) to parallelize the least-square solver. We use LSQR (Paige & Saunders, 1982) as the default least-square solver as we found it to be slightly faster. The only exception is the large |A| experiment, where we used the SVD-based minimum norm solver (Anderson et al., 1999) to ensure reproducibility. Algorithm 1 Phased TD/DAE Require: n, k, alg∈{TD, DAE}, LSTSQ_SOLVER 1: Initialize V ≡ 0 2: for T = 1, 2, . . . do 3: D = {} 4: for s ∈ S do 5: for i = 1, . . . , n do 6: Sample k-step trajectory τ from environment 7: D ← D ∪ {τ } 8: end for 9: end for 10: if alg == TD then 11: Compute R, P from D 12: V ← R + γ k PV 13: else 14: Compute M, R, P from D 15: V, A ← LSTSQ_SOLVER(||R + γ k PV − V̂ − MÂ||2 ) 16: end if 17: end for 9.1
Additional Experiments
The Low Coverage Case With an Iterative Solver (|A| ∈ {4, 16, 64}, n = 8, k = 16, ps = 0, pr = 0) In Section 5, we showed that increasing the size of the action space results in DAE converging to suboptimal solutions when regularized with minimum norm solutions. In Figure 5, we rerun the same experiment but with an iterative least-sqaures solver (LSQR (Paige & Saunders, 1982) in this case), where the optimum in the previous phase is used as the initialization for the current phase. We find that DAE converges again to the true value function, although at a slower rate as |A| increases. This might partially explain the success of DAE in the deep RL setting (Pan et al., 2022), where gradient-based optimization is used. 16
0.4 MSE (V, VBE)
MSE (V, V )
1.0 0.5 0.0
100
101
102 Phase
4
| | 16
64
0.2 0.0
103 TD
DAE
100
101 Method Bellman
102 Phase MC
103 MC-A
Figure 5: The low coverage case with an iterative solver. Lines and shadings represent (mean ± 3 standard error). Note that TD does not depend on |A|.
17