Stochastic Variance Reduction Methods for Policy Evaluation
Simon S. Du Jianshu Chen Lihong Li Lin Xiao Dengyong Zhou



§ INTRODUCTION



Reinforcement learning (RL) is a powerful learning paradigm for 
sequential decision making 
<cit.>.
An RL agent interacts with the environment by repeatedly
observing the current state, taking an action according to a certain policy,
receiving a reward signal and transitioning to a next state. 
A policy specifies which action to take given the current state. 
Policy evaluation estimates a value function that predicts expected
cumulative reward the agent would receive by following a fixed policy
starting at a certain state.
In addition to quantifying long-term values of states, which can be of interest on its own,

value functions also provide important information for the agent to optimize its policy.
For example, policy-iteration algorithms iterate between policy-evaluation steps 

and policy-improvement steps, until a (near-)optimal policy is found <cit.>.

Therefore, estimating the value function efficiently and accurately is essential in RL.

There has been substantial work on policy evaluation, with temporal-difference (TD) methods being perhaps the most popular.


These methods use the Bellman equation to bootstrap the estimation process. 
Different cost functions are formulated to exploit this idea, 
leading to different policy evaluation algorithms; 






see <cit.> for a comprehensive survey.
In this paper, we study policy evaluation by minimizing 
the mean squared projected Bellman error (MSPBE) 
with linear approximation of the value function.
We focus on the batch setting where a fixed, finite dataset is given.




This fixed-data setting is not only important in itself <cit.>, but also an important component in other
RL methods such as experience replay <cit.>.









The finite-data regime makes it possible to solve policy evaluation more efficiently with recently developed fast optimization methods based on 
stochastic variance reduction, such as 
SVRG <cit.> and SAGA <cit.>.
For minimizing strongly convex functions with a finite-sum structure, such methods enjoy the same low computational cost per iteration as the classical stochastic gradient method, but also achieve fast, linear convergence rates (i.e., exponential decay of the optimality gap in the objective).  However, they cannot be applied directly to minimize the MSPBE, whose objective does not have the finite-sum structure.
In this paper,
we overcome this obstacle by transforming the empirical MSPBE problem 

to an equivalent convex-concave saddle-point problem that possesses the desired finite-sum structure.

In the saddle-point problem, we consider the model parameters as the primal
variables, which are coupled with the dual variables through a bilinear term. 
Moreover, without an ℓ_2-regularization on the model parameters, the objective 

is only strongly concave in the dual variables,
but not 

in the primal variables. 

We propose a primal-dual batch gradient method,
as well as two stochastic variance-reduction methods based on
SVRG and SAGA, respectively.
Surprisingly, we show that when the coupling matrix is full rank, these algorithms 
achieve linear convergence in both the primal and dual spaces, despite the lack of strong convexity of the objective in the primal variables.

Our results also extend to off-policy learning and TD with
eligibility traces <cit.>.


We note that <cit.> have extended 
both SVRG and SAGA to solve convex-concave saddle-point problems with
linear-convergence guarantees. The main difference between our results and theirs
are

  
  * Linear convergence in <cit.> relies on the assumption that the objective is strongly convex in the primal variables and strongly concave in the dual.



    Our results show, somewhat surprisingly, that only one of them is necessary if the primal-dual
    coupling is bilinear and the coupling matrix is full rank. 
    In fact, we are not aware of similar previous results even for the 
    primal-dual batch gradient method, which we show in this paper.
  
  * Even if a strongly convex regularization on the primal variables is introduced to the MSPBE objective,
    the algorithms in <cit.> 
    cannot be applied efficiently.


    Their algorithms require that the proximal mappings of 
    the strongly convex and concave regularization functions 
    
    be computed efficiently.
    In our saddle-point formulation, the strong concavity of the dual variables
    comes from a quadratic function defined by the feature covariance matrix,
    which cannot be inverted efficiently and makes the proximal mapping costly
    to compute.
    Instead, our algorithms only use its (stochastic) gradients and hence are much more efficient.

We compare various gradient based algorithms on a Random MDP and Mountain Car 
data sets. 
The experiments demonstrate the effectiveness of our proposed methods.






§ PRELIMINARIES






We consider a Markov Decision Process 
(MDP) <cit.> described by
(𝒮,𝒜,𝒫_ss'^a,ℛ,γ), 
where 𝒮 is the set of states, 
𝒜 the set of actions,
𝒫_ss'^a the transition probability from state 
s to state s' after taking action 
a, ℛ(s,a) the reward received after taking action a in state s, and γ∈ [0,1) a discount factor.

The goal of an agent is to find an action-selection policy π, 

so that the long-term reward under this policy is maximized.  For ease of exposition, we assume 𝒮 is finite, but none of our results relies on this assumption.

A key step in many algorithms in RL is to estimate the value function of 
a given policy π, defined as
V^π(s)	 [
						∑_t=0^∞γ^t ℛ(s_t,a_t) 
						| s_0 = s, π
					].


Let V^π denote a vector constructed by stacking the values of V^π(1),…, V^π(|S|) on top of each other. Then V^π is the unique fixed point of the Bellman operator T^π:
	
    V^π = T^π V^π≜ R^π + γ P^π V^π ,
 
where R^π is the expected reward vector under policy π,
defined elementwise as R^π(s)=_π(a|s)ℛ(s,a);
and P^π is the transition matrix induced by the policy applying π,
defined entrywise as P^π(s,s')=_π(a|s)𝒫^a_ss'.



 §.§ Mean squared projected Bellman error (MSPBE)


One approach to scale up when the state space size 𝒮 is large or infinite is to use a linear approximation for V^π.
Formally, we use a feature map ϕ: 𝒮→ℝ^d and
approximate the value function by 
V^π(s) = ϕ(s)^T θ,
where θ∈ℝ^d is the model parameter to be estimated. 
Here, we want to find θ that minimizes 


the mean squared projected Bellman error, or MSPBE:


    MSPBE(θ) 
    			   1/2V^π - Π T^πV^π_Ξ^2 ,

where Ξ is a diagonal matrix with diagonal elements being the stationary
distribution over 𝒮 induced by the policy π, and Π is the
weighted projection matrix onto the linear space spanned by
ϕ(1),…, ϕ(|S|), that is, 

    Π   =
    				Φ (Φ^T ΞΦ)^-1Φ^T Ξ
 
where Φ [ ϕ^T(1),…,ϕ^T(|S|)] is the matrix obtained
by stacking the feature vectors row by row.
Substituting (<ref>) and (<ref>) 
into (<ref>),

we obtain <cit.>

    MSPBE(θ)
    			   =
    				1/2Φ^T Ξ (V^π - T^πV^π ) 
    				^2_(Φ^T ΞΦ)^-1 .

We can further rewrite the above expression for MSPBE as a standard weighted 
least-squares problem:

    MSPBE(θ)
    			   =
    				1/2 
    					A θ - b
    				_C^-1^2,

with properly defined A, b and C, described as follows.
Suppose the MDP under policy π settles at its stationary distribution 
and generates an infinite transition sequence
{(s_t, a_t,r_t,s_t+1)}_t=1^∞,
where s_t is the current state, a_t is the action, 
r_t is the reward, and s_t+1 is the next state.
Then with the definitions
ϕ_t ϕ(s_t) and ϕ_t' ϕ(s_t+1),
we have

    A = [ ϕ_t (ϕ_t - γϕ_t')^T ], 
       b = [ ϕ_t r_t ], 
       C =  [ϕ_t ϕ_t^T ],

where [·] are with respect to the stationary distribution.  Many TD  solutions converge to a minimizer of MSPBE in the limit <cit.>.



 §.§ Empirical MSPBE



In practice, quantities in (<ref>) are often unknown, and we only have access to a finite dataset with n transitions 𝒟 = {(s_t, a_t,r_t,s_t+1)}_t=1^n.



By replacing the unknown statistics with their finite-sample estimates, we obtain the Empirical MSPBE, or EM-MSPBE.  Specifically, let

    A   ≜1/n∑_t=1^n A_t,   b≜1/n∑_t=1^n b_t,   C≜1/n∑_t=1^nC_t,

where for t=1,…,n, 

    A_t ≜ϕ_t(ϕ_t-γϕ'_t)^T,   
      b_t ≜ r_tϕ_t,    
      C_t ≜ϕ_tϕ_t^T.

EM-MSPBE with an optional ℓ_2-regularization is given by:

    EM-MSPBE(θ) 
    		   = 
    			1/2Aθ - b_C^-1^2
    			+
    			ρ/2θ^2 ,

where ρ≥ 0 is a regularization factor.

Observe that (<ref>) is a (regularized) weighted least squares problem. 
Assuming C is invertible, its optimal solution is

    θ^⋆ = (A^⊤C^-1A + ρ I)^-1A^⊤C^-1b.

Computing θ^⋆ directly requires O(nd^2) operations to form the 
matrices A, b and C, and then 
O(d^3) operations to complete the calculation.
This method, known as least-squares temporal difference or LSTD <cit.>,
can be very expensive when n and d are large.
One can also skip forming the matrices explicitly and
compute θ^⋆ using n recusive rank-one updates 
<cit.>.
Since each rank-one update costs O(d^2), the total cost is O(nd^2).


In the sequel, we develop efficient algorithms to minimize
EM-MSPBE by using stochastic variance reduction methods, which samples one
(ϕ_t, ϕ_t') per update without pre-computing A,
b and C.  These algorithms not only maintain a low O(d) per-iteration computation cost, but also attain fast linear convergence rates with a log(1/ϵ) dependence 
on the desired accuracy ϵ.






§ SADDLE-POINT FORMULATION OF EM-MSPBE



Our algorithms (in Section <ref>) are based on the stochastic variance
reduction techniques developed for minimizing a finite sum of convex functions, more specifically, SVRG <cit.> and 
SAGA <cit.>.
They deal with problems of the form

    min_x∈^d { f(x) ≜1/n∑_i=1^n f_i(x) },

where each f_i is convex.
We immediately notice that the EM-MSPBE in (<ref>) cannot 
be put into such a form, even though the matrices A, b
and C have the finite-sum structure given in (<ref>).
Thus, extending variance reduction techniques to EM-MSPBE minimization
is not straightforward.

Nevertheless, we will show that the minimizing the EM-MSPBE is equivalent to
solving a convex-concave saddle-point problem which actually possesses the
desired finite-sum structure.
To proceed, we resort to the machinery of conjugate functions
<cit.>. 
For a function f: R^d →R, its
conjugate function f^⋆: R^d →R is defined as
f^⋆(y) sup_x (y^T x - f(x)). 
Note that the conjugate function of 1/2 x _C^2 is
1/2 y _C^-1^2, i.e.,  

    1/2 y _C^-1^2
    			   = max_x(y^T x - 1/2 x _C^2 ).

With this relation, we can rewrite EM-MSPBE in (<ref>) as

    max_w( w^T(b-Aθ)
    					- 1/2 w _C^2
    				)
    				+
    				ρ/2θ^2 ,

so that minimizing EM-MSPBE is equivalent to solving

    min_θ∈R^dmax_w ∈R^d {ℒ(θ, w) 
        = 1/n∑_t=1^n ℒ_t (θ, w) },

where the Lagrangian, defined as

    ℒ(θ, w) ≜ρ/2θ^2
    				-
    				w^T Aθ
    				-
    				(
    					1/2 w _C^2
    					-
    					w^T b) ,


may be decomposed using  (<ref>), with

    ℒ_t(θ, w) ρ/2θ^2 - w^T A_t θ - (1/2w_C_t^2 - w^T b_t) .

Therefore, minimizing the EM-MSPBE is equivalent to solving the saddle-point
problem (<ref>), which is convex in the primal 
variable θ and concave in the dual variable w.
Moreover, it has a finite-sum structure similar to (<ref>).

<cit.> and <cit.> independently 
showed that the GTD2 algorithm <cit.>
is indeed a stochastic gradient method for solving the 
saddle-point problem (<ref>), although
they obtained the saddle-point formulation with different derivations.
More recently, <cit.> used the conjugate function approach
to obtain saddle-point formulations for a more general class of problems
and derived primal-dual stochastic gradient algorithms for solving them.
However, 
these algorithms have sublinear convergence rates, which leaves much room
to improve when applied to problems with finite datasets. Recently, <cit.>
developed SVRG methods for a general finite-sum composition optimization that achieve
linear convergence rate. Different from our methods, their stochastic gradients are biased and they have
worse dependency on the condition numbers (κ^3 and κ^4).


The fast linear convergence of our algorithms presented in
Sections <ref> and <ref>
requires the following assumption:

A has full rank, C is strictly positive definite, and the feature vector ϕ_t is uniformly bounded.





Under mild regularity conditions <cit.>, 
we have A and C converge in probability
to A and C defined in (<ref>), respectively.
Thus, if the true statistics A is non-singular and
C is positive definite, 
and we have enough training samples, these assumptions are
usually satisfied. 
They have been widely used in previous works on gradient-based algorithms
<cit.>.

A direct consequence of Assumption <ref> is that 
θ^⋆ in (<ref>) is the unique minimizer of
the EM-MSPBE in (<ref>), even without any strongly convex
regularization on θ (i.e., even if ρ=0).
However, if ρ=0, then the Lagrangian ℒ(θ,w) is 
only strongly concave in w, but not strongly convex in θ.
In this case, we will show that non-singularity of the coupling 
matrix A can “pass” an implicit strong convexity 
on θ, which is exploited by our algorithms to obtain 
linear convergence in both the primal and dual spaces.







§ A PRIMAL-DUAL BATCH GRADIENT METHOD



Before diving into the stochastic variance reduction algorithms, we first
present Algorithm <ref>, which is a primal-dual 
batch gradient (PDBG) algorithm for solving the
saddle-point problem (<ref>).
In Step 2, the vector B(θ,w) is obtained by
stacking the primal and negative dual gradients:

    B(θ,w)  ≜[  ∇_θ L(θ,w); -∇_w L(θ,w) ]
                        =
    					[   ρθ - A^T w; Aθ - b + C w ] .

Some notation is needed in order to characterize the convergence rate of 
Algorithm <ref>.
For any symmetric and positive definite matrix S, 
let λ_max(S) and λ_min(S) denote its 
maximum and minimum eigenvalues respectively, and define its condition number to be κ(S)λ_max(S)/λ_min(S).
We also define L_ρ and μ_ρ for any ρ≥0:

    L_ρ   λ_max(ρ I + A^T C^-1A), 
    μ_ρ   λ_min(ρ I + A^T C^-1A).

By Assumption <ref>, we have L_ρ≥μ_ρ>0.
The following theorem is proved in Appendix <ref>.

Suppose Assumption <ref> holds and let 
(θ_⋆, w_⋆) be the (unique) solution 
of (<ref>).
If the step sizes are chosen as
σ_θ = 1/ 9 L_ρκ(C) and
σ_w = 8/9λ_max(C),
then the number of iterations of Algorithm <ref> 
to achieve θ-θ_⋆^2 + w - w_⋆^2 ≤ϵ^2 
is upper bounded by

    O
    	(		
    		κ(
    			ρ I + A^T C^-1A) 
    		·κ(C)
    		·log(1/ϵ)
    	) .


We assigned specific values to the step sizes σ_θ and σ_w
for clarity.
In general, we can use similar step sizes while keeping their ratio
roughly constant as
σ_w/σ_θ≈8 L_ρ/λ_min(C); see Appendices <ref> and <ref> for more details.  In practice, one can use a parameter search on a small subset of data to find reasonable step sizes.  It is an interesting open problem how to automatically select and adjust step sizes.





Note that the linear rate is determined by two parts: (i) the strongly convex
regularization parameter ρ, and (ii) the positive definiteness of
A^T C^-1A.
The second part could be
interpreted as transferring  strong concavity in dual variables via the
full-rank bi-linear coupling matrix A. For this reason, even if
the saddle-point problem (<ref>) has only strong
concavity in dual variables (when ρ=0), 
the algorithm still enjoys a linear convergence rate. 

Moreover, even if ρ>0, it will be inefficient to solve 
problem (<ref>) using primal-dual algorithms based on 
proximal mappings of the strongly convex and concave terms
<cit.>.
The reason is that, in (<ref>), the strong concavity of the 
Lagrangian with respect to the dual lies in the quadratic
function (1/2)w_C, whose proximal mapping cannot be computed
efficiently.
In contrast, the PDBG algorithm only needs its gradients.

If we pre-compute and store A, b and
C, which costs O(nd^2) operations, then computing
the gradient operator B(θ,w) in (<ref>)
during each iteration of PDBG costs O(d^2) operations.
Alternatively, if we do not want to store these d× d matrices
(especially if d is large),
then we can compute B(θ,w) as finite sums on the fly. 
More specifically,
B(θ,w) = 1/n∑_t=1^n B_t(θ, w),
where for each t=1,…,n,

    B_t(θ,w) = [[         ρθ - A_t w; A_tθ - b_t + C_t w ]].

Since A_t, b_t and C_t are all rank-one matrices,
as given in (<ref>),
computing each B_t(θ,w) only requires O(d) operations. 
Therefore, computing B(θ,w) costs O(nd) 
operations as it averages B_t(θ,w) over n samples. 










§ STOCHASTIC VARIANCE REDUCTION METHODS



If we replace B(θ,w) in Algorithm <ref> (line 2) 
by the stochastic gradient B_t(θ,w) in (<ref>), 
then we recover the GTD2 algorithm of <cit.>, 
applied to a fixed dataset, possibly with multiple passes. 
It has a low per-iteration cost but a slow, sublinear convergence rate.
In this section, we provide two stochastic variance reduction methods 
and show they achieve fast linear convergence.





 §.§ SVRG for policy evaluation

Algorithm <ref> is adapted from the stochastic variance
reduction gradient (SVRG) method <cit.>.
It uses two layers of loops and maintains two sets of parameters
(θ̃,w̃) and (θ,w). In the outer loop, the algorithm
computes a full gradient B(θ̃,w̃) using (θ̃,
w̃), which takes O(nd) operations. Afterwards, the
algorithm executes the inner loop, which randomly samples an index t_j and 
updates (θ,w) using variance-reduced stochastic gradient: 
	
    B_t_j(θ,w,θ̃,w̃) = B_t_j(θ,w) + B(θ̃,w̃) - B_t_j(θ̃,w̃).

Here, B_t_j(θ,w) contains the stochastic gradients at (θ,w) 
computed using the random sample with index t_j, and
B(θ̃,w̃) - B_t_j(θ̃,w̃) is a term
used to reduce the variance in B_t_j(θ,w) while keeping
B_t_j(θ,w,θ̃,w̃) 
an unbiased estimate of B(θ,w). 

Since B(θ̃,w̃) is computed once during each iteration of 
the outer loop with cost O(nd) 
(as explained at the end of Section <ref>),
and each of the N iterations of the inner loop cost O(d) operations,
the total computational cost of for each outer loop is O(nd+Nd). 
We will present the overall complexity analysis of 
Algorithm <ref> in Section <ref>.



 §.§ SAGA for policy evaluation





The second stochastic variance reduction method for policy evaluation
is adapted from SAGA <cit.>;
see Algorithm <ref>.
It uses a single loop, and maintains a single set of parameters (θ,w). 
Algorithm <ref> starts by first computing each component
gradients g_t=B_t(θ,w) at the initial point, and also form their average
B=∑_t^n g_t.
At each iteration, the algorithm randomly picks an index t_m∈{1,…,n} 
and computes the stochastic gradient h_t_m=B_t_m(θ,w).
Then, it updates (θ,w) using a variance reduced stochastic gradient:
B + h_t_m - g_t_m,
where g_t_m is the previously computed stochastic gradient using the 
t_m-th sample (associated with certain past values of θ and w).
Afterwards, it updates the batch gradient estimate B as 
B+1/n(h_t_m-g_t_m) and replaces g_t_m with h_t_m. 

As Algorithm <ref> proceeds, different vectors g_t are 
computed using different values of θ and w
(depending on when the index t was sampled).
So in general we need to store all vectors g_t, for t=1,…,n, 
to facilitate individual updates, which will cost additional O(nd) storage.
However, by exploiting the rank-one structure in (<ref>),
we only need to store three scalars
(ϕ_t-γ_ϕ')^Tθ, (ϕ_t-γ_ϕ')^T w, and ϕ_t^Tw, 
and form g_t_m on the fly using O(d) computation.
Overall, each iteration of SAGA costs O(d) operations.




 §.§ Theoretical analyses of SVRG and SAGA



In order to study the convergence properties of SVRG and SAGA for 
policy evaluation, we introduce a smoothness parameter L_G
based on the stochastic gradients B_t(θ,w). 
Let β = σ_w/σ_θ be the ratio between the
primal and dual step-sizes, and define a pair of weighted Euclidean norms 

    Ω(θ,w)     (θ^2 + β^-1w^2)^1/2, 
    Ω^*(θ,w)     (θ^2 + βw^2)^1/2.

Note that Ω(·,·) upper bounds the error in optimizing θ: Ω(θ-θ_⋆,w-w_⋆) ≥θ-θ_⋆.  Therefore, any bound on Ω(θ-θ_⋆,w-w_⋆) applies automatically to θ-θ_⋆.

Next, we define the parameter L_G through its square:
	
    L_G^2
    			   sup_θ_1,w_1, θ_2, w_21/n∑_t=1^n
    					Ω^*(B_t(θ_1,w_1)-B_t(θ_2,w_2))^2
    				/Ω(θ_1-θ_2, w_1-w_2)^2
    				.

This definition is similar to the smoothness constant L̅ used in
<cit.> except that we used the step-size ratio
β rather than the strong convexity and concavity parameters
of the Lagrangian to define
Ω and Ω^*.[Since our saddle-point problem is not
  necessarily strongly convex in θ (when ρ=0), 
  we could not define Ω and Ω^* in the same way 
  as <cit.>.]
Substituting the definition of B_t(θ,w) 
in (<ref>), we have 

    L_G^2 = 1/n∑_t=1^n  G_t^T G_t,
        
    		G_t
    			[         ρ I -√(β) A_t^T;    √(β) A_t       β C_t ].
















With the above definitions, 
we characterize the convergence of
Ω(θ_m-θ_⋆, w_m-w_⋆),
where (θ_⋆, w_⋆) is the solution 
of (<ref>), and (θ_m,w_m) is the output of 
the algorithms after the m-th iteration.
For SVRG, it is the m-th outer iteration
in Algorithm <ref>.
The following two theorems are proved in 
Appendices <ref> and <ref>, 
respectively.



Suppose Assumption <ref> holds.
If we choose
σ_θ = μ_ρ/48κ(C) L_G^2,  
σ_w=8L_ρ/λ_min(C)σ_θ,
N=51 κ^2(C) L_G^2/μ^2_ρ, 

where L_ρ and μ_ρ are defined in (<ref>) 
and (<ref>),
then 

    [Ω(θ_m-θ_⋆, w_m-w_⋆)^2]
     
      ≤(4/5)^m
      Ω(θ_0-θ_⋆, w_0-w_⋆)^2 .



The overall computational cost for reaching
[Ω(θ_m-θ_⋆, w_m-w_⋆)]≤ϵ
is upper bounded by

    O(
    (n+κ(C) L_G^2/λ_min^2(ρ I + A^TC^-1A))d log(1/ϵ)
    ).






Suppose Assumption <ref> holds.
If we choose 
σ_θ = μ_ρ/3(8κ^2(C)L_G^2 + n μ_ρ^2)
and σ_w=8L_ρ/λ_min(C)σ_θ
in Algorithm <ref>,
then

    [Ω(θ_m-θ_⋆, w_m-w_⋆)^2]
      ≤ 2 (1-ρ)^m
      Ω(θ_0-θ_⋆, w_0-w_⋆)^2 ,

where 
ρ≥μ_ρ^2/9(8κ^2(C)L_G^2 + nμ_ρ^2).
The total cost to achieve 
[Ω(θ_m-θ_⋆, w_m-w_⋆)]≤ϵ
has the same bound in (<ref>).


Similar to our PDBG results in (<ref>), 
both the SVRG and SAGA algorithms for policy evaluation enjoy
linear convergence even if there is no strong convexity in the 
saddle-point problem (<ref>) (i.e., when ρ = 0).
This is mainly due to the
positive definiteness of A^T C^-1A when
C is positive-definite and A is full-rank. 
In contrast, the linear convergence of SVRG and SAGA in
<cit.> requires 
the Lagrangian to be both strongly convex in θ and
strongly concave in w.


Moreover, in the policy evaluation problem,
the strong concavity with respect to the dual variable w
comes from a weighted quadratic norm (1/2)w_C, which
does not admit an efficient proximal mapping as required by
the proximal versions of SVRG and SAGA in
<cit.>.
Our algorithms only require computing the stochastic gradients of 
this function, which is easy to do due to its finite sum structure.

<cit.> also proposed accelerated variants of 
SVRG and SAGA using the “catalyst” framework of 
<cit.>.
Such extensions can be done similarly for the three algorithms presented
in this paper, and we omit the details due to space limit.









§ COMPARISON OF DIFFERENT ALGORITHMS


This section compares the computation complexities of several representative policy-evaluation algorithms that minimize EM-MSPBE, as summarized in
Table <ref>. 

The upper part of the table lists algorithms
whose complexity is linear in feature dimension d, including the two new algorithms presented in the previous section.

We can also apply GTD2
to a finite dataset with samples drawn uniformly at random with replacement.

It costs O(d) per iteration, but has a sublinear 
convergence rate regarding ϵ.  
In practice, people may choose ϵ=Ω(1/n) for generalization reasons (see, e.g., <cit.>), leading to an O(κ'nd) overall complexity for GTD2, where κ' is a condition number related to the algorithm.  However, as verified by our experiments, the bounds in the table show that our SVRG/SAGA-based algorithms are much faster as their effective condition numbers vanish when n becomes large.  TDC has a similar complexity to GTD2.

In the table, we list two different implementations of PDBG. PDBG-(I)
computes the gradients by averaging the stochastic gradients over the entire
dataset at each iteration, which costs O(nd) operations; see discussions at the end of Section <ref>.
PDBG-(II) first pre-computes the matrices A, b and 
C using O(nd^2) operations, then computes the batch
gradient at each iteration with O(d^2) operations. 
If d is very large (e.g., when d ≫ n), then PDBG-(I) would have an advantage over PDBG-(II). 
The lower part of the table also includes LSTD, which has O(nd^2) complexity if rank-one updates are used.

SVRG and SAGA are more efficient than the other algorithms, when either d or n is very large. 
In particular, they have a lower complexity than LSTD when 

 d > (1 + κ(C)κ_G^2/n ) log(1/ϵ),
This condition is easy to satisfy, when n is very large.

On the other hand, SVRG and SAGA algorithms are more efficient than PDBG-(I) if n is large, say 
n > κ(C)κ_G^2/
(κ(C)κ - 1),
where κ and κ_G are described in the caption of
Table <ref>.



There are other algorithms whose complexity scales linearly with n and d, including iLSTD <cit.>, and TDC <cit.>, fLSTD-SA <cit.>,
and the more recent algorithms of <cit.> and 
<cit.>.

However, their convergence is slow: 
the number of iterations required to reach a desired accuracy ϵ 
grows as 1/ϵ or worse.

The CTD algorithm <cit.> uses a similar idea as SVRG to reduce variance in TD updates.  This algorithm is shown to have a similar linear convergence rate in an online setting where the data stream is generated by a Markov process with finite states and exponential mixing.  The method solves for a fixed-point solution by stochastic approximation.  As a result, they can be non-convergent in off-policy learning, while our algorithms remain stable (c.f., Section <ref>).


























§ EXTENSIONS



Other objectives than MSPBE can be used for
policy evaluation. We choose MSPBE in this paper
because it has both theoretical and empirical advantages over other objective
functions. Nevertheless, we believe that the saddle-point formulation in
Section <ref> and the algorithms in 
Sections<ref> and <ref>) could be
generalized to other similar objective functions such as mean squared Bellman
error (MSBE) and the norm of the expected TD update (NEU). 
We refer the reader to
the survey by <cit.> for a detailed account
of these objective functions.


It is possible to extend our approach to accelerate optimization of other objectives such as MSBE and NEU <cit.>.  In this section, we briefly describe two extensions of the algorithms developed earlier.






 §.§ Off-policy learning



In some cases, we may want to estimate the value function of a
policy π from a set of data 𝒟 generated 
by a different “behavior” policy π_b. 
This is called off-policy learning <cit.>.

In the off-policy case, samples are generated from the distribution induced by the behavior policy π_b, not the the target policy π.  While such a mismatch often causes stochastic-approximation-based methods to diverge <cit.>, our gradient-based algorithms remain convergent with the same (fast) convergence rate.


Consider the RL framework outlined in Section <ref>.  For each state-action pair (s_t,a_t) such that π_b(a_t|s_t) > 0, we define the importance ratio, ρ_t π(a_t|s_t)/π_b(a_t|s_t).  The EM-MSPBE for off-policy learning has the same expression as in (<ref>) except that A_t, b_t and C_t are modified by the weight factor ρ_t, as listed in Table <ref>; see also <cit.> for a related discussion.)  Algorithms <ref>–<ref> remain the same for the off-policy case after A_t, b_t and C_t are modified correspondingly.








 §.§ Learning with eligibility traces


Eligibility traces are a useful technique to trade off bias and variance in TD learning <cit.>.
When they are used, we can pre-compute z_t in Table <ref> before running our new algorithms. 
Note that EM-MSPBE with eligibility traces has the same form of
(<ref>), with A_t, b_t and C_t defined differently according to the last row of Table <ref>. 
At the m-th step of the learning process, the algorithm randomly samples z_t_m, ϕ_t_m, ϕ_t_m' and r_t_m from the fixed dataset and computes the corresponding stochastic gradients, where the index t_m is uniformly distributed over {1,…,n} and are independent for different values of m. 

Algorithms <ref>–<ref>
immediately work for this case, enjoying a similar linear convergence rate and a computation complexity linear in n and d.
We need  additional O(nd) operations to pre-compute z_t recursively and an additional O(nd) storage for z_t. 
However, it does not change the order of the total complexity for SVRG/SAGA. 

 














§ EXPERIMENTS































In this section, we compare the following algorithms on two benchmark problems: 
(i) PDBG (Algorithm <ref>); 
(ii) GTD2 with samples drawn randomly with replacement from a dataset;


(iii) TD: the fLSTD-SA algorithm of <cit.>;



(iv) SVRG (Algorithm <ref>);


and (v) SAGA (Algorithm <ref>).

Note that when ρ > 0, the TD solution and EM-MSPBE minimizer

differ, so we do not include TD.


For step size tuning, σ_θ is chosen from {10^-1,10^-2,…,10^-6}1/L_ρκ(C) and σ_w is chosen from  {1,10^-1,10^-2}1/λ_max(C). 
We only report the results of each algorithm which correspond to the best-tuned step sizes; for SVRG we choose N=2n.
















In the first task, we consider a randomly generated MDP with 400 states and 10 actions <cit.>. 
The transition probabilities are defined as P(s'|a,s) ∝ p_ss'^a + 10^-5, where p_ss'^a ∼ U[0,1].

The data-generating policy and start distribution were generated in a similar way. 
Each state is represented by a 201-dimensional feature vector, where 200 of the features were sampled from a uniform distribution, and the last feature was constant one.

We chose γ=0.95.

Fig. <ref> shows the performance of various algorithms for n=20000.
First, notice that the stochastic variance methods converge much faster than others.  
In fact, our proposed methods achieve linear convergence. 
Second, as we increase ρ, the performances of PDBG, SVRG and SAGA improve significantly due to better conditioning, as predicted by our theoretical results.





Next, we test these algorithms on Mountain Car <cit.>.  To collect the dataset, we first ran Sarsa with d=300 CMAC features to obtain a good policy.  Then, we ran this policy to collect trajectories that comprise the dataset.


Figs. <ref> and <ref> show our proposed stochastic variance reduction methods dominate other first-order methods.  Moreover, with better conditioning (through a larger ρ), PDBG, SVRG and SAGA achieve faster convergence rate.
Finally, as we increase sample size n, SVRG and SAGA converge faster.
This simulation verifies our theoretical finding in Table <ref> that SVRG/SAGA need fewer epochs for large n.





























§ CONCLUSIONS



In this paper, we reformulated the EM-MSPBE minimization problem in policy evaluation into an empirical saddle-point problem, and developed and analyzed a batch gradient method and two first-order stochastic variance reduction methods to solve the problem. 

An important result we obtained is that even when the reformulated saddle-point problem lacks strong convexity in primal variables and has only strong concavity in dual variables, the proposed algorithms are still able to achieve a linear convergence rate.  We are not aware of any similar results for primal-dual batch gradient methods or stochastic variance reduction methods. Furthermore, we showed that when both the feature dimension d and the number of samples n are large, the developed stochastic variance reduction methods are more efficient than any other gradient-based methods which are convergent in off-policy settings. 


This work leads to several interesting directions for research.  First, we believe it is important to extend the stochastic variance reduction methods to nonlinear approximation paradigms <cit.>, especially with deep neural networks.

Moreover, it remains an important open problem how to apply stochastic variance reduction techniques to policy optimization.
