Stochastic Variance Reduction Methods for Policy Evaluation

Simon S. Du 1 Jianshu Chen 2 Lihong Li 2 Lin Xiao 2 Dengyong Zhou 2

arXiv:1702.07944v2 [cs.LG] 9 Jun 2017

Abstract
Policy evaluation is concerned with estimating
the value function that predicts long-term values of states under a given policy. It is a crucial step in many reinforcement-learning algorithms. In this paper, we focus on policy evaluation with linear function approximation over
a fixed dataset. We first transform the empirical policy evaluation problem into a (quadratic)
convex-concave saddle-point problem, and then
present a primal-dual batch gradient method, as
well as two stochastic variance reduction methods for solving the problem. These algorithms
scale linearly in both sample size and feature dimension. Moreover, they achieve linear convergence even when the saddle-point problem has
only strong concavity in the dual variables but no
strong convexity in the primal variables. Numerical experiments on benchmark problems demonstrate the effectiveness of our methods.

1

Introduction

Reinforcement learning (RL) is a powerful learning
paradigm for sequential decision making (see, e.g., Bertsekas & Tsitsiklis, 1995; Sutton & Barto, 1998). 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
1

Machine Learning Department, Carnegie Mellon University, Pittsburgh, Pennsylvania 15213, USA. 2 Microsoft Research, Redmond, Washington 98052, USA.. Correspondence
to: Simon S. Du <ssdu@cs.cmu.edu>, Jianshu Chen <jianshuc@microsoft.com>, Lihong Li <lihongli@microsoft.com>,
Lin Xiao <lin.xiao@microsoft.com>, Dengyong Zhou <denzho@microsoft.com>.
Proceedings of the 34 th International Conference on Machine
Learning, Sydney, Australia, PMLR 70, 2017. Copyright 2017
by the author(s).

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 (Bertsekas & Tsitsiklis,
1995; Lagoudakis & Parr, 2003). 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 Dann et al. (2014) 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 (Lange et al., 2011), but also an important component
in other RL methods such as experience replay (Lin, 1992).
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 (Johnson & Zhang, 2013) and SAGA (Defazio et al., 2014). 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

Stochastic Variance Reduction Methods for Policy Evaluation

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 (Sutton &
Barto, 1998; Precup et al., 2001).

where Rπ is the expected reward vector under policy π,
defined elementwise as Rπ (s) = Eπ(a|s) R(s, a); and P π
is the transition matrix induced by the policy applying π,
a
defined entrywise as P π (s, s0 ) = Eπ(a|s) Pss
0.

We note that Balamurugan & Bach (2016) have extended
both SVRG and SAGA to solve convex-concave saddlepoint problems with linear-convergence guarantees. The
main difference between our results and theirs are

2.1

• Linear convergence in Balamurugan & Bach (2016)
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 Balamurugan & Bach (2016) 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.

2

Preliminaries

We consider a Markov Decision Process (MDP) (Puterman,
a
2005) described by (S, A, Pss
0 , R, γ), where S is the set of
a
states, A the set of actions, Pss
0 the transition probability
0
from state s to state s after taking action a, R (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 actionselection policy π, so that the long-term reward under this
policy is maximized. For ease of exposition, we assume S
is finite, but none of our results relies on this assumption.
A key step in many algorithms in RL is to estimate the
value
of a given policy π, defined as V π (s) ,
P∞function
t
E[ t=0 γ R(st , at )|s0 = 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 π :

Mean squared projected Bellman error (MSPBE)

One approach to scale up when the state space size |S| is
large or infinite is to use a linear approximation for V π .
Formally, we use a feature map φ : S → Rd and approximate the value function by Vb π (s) = φ(s)T θ, where
θ ∈ Rd is the model parameter to be estimated. Here, we
want to find θ that minimizes the mean squared projected
Bellman error, or MSPBE:
MSPBE (θ) ,

V

π

=T V

π

π

π

π

, R + γP V ,

(1)

(2)

where Ξ is a diagonal matrix with diagonal elements being the stationary distribution over S induced by the policy π, and Π is the weighted projection matrix onto the
linear space spanned by φ(1), . . . , φ(|S|), that is,
Π = Φ(ΦT ΞΦ)−1 ΦT Ξ

(3)

where Φ , [φT (1), . . . , φT (|S|)] is the matrix obtained by
stacking the feature vectors row by row. Substituting (3)
and (1) into (2), we obtain (see, e.g., Dann et al., 2014)
MSPBE(θ) =

1 T bπ
kΦ Ξ(V − T π Vb π )k2(ΦT ΞΦ)−1 .
2

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

1
kAθ − bk2C −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
∞
{(st , at , rt , st+1 )}t=1 , where st is the current state, at is
the action, rt is the reward, and st+1 is the next state. Then
with the definitions φt , φ(st ) and φ0t , φ(st+1 ), we have
A = E[φt (φt − γφ0t )T ], b = E[φt rt ], C = E[φt φTt ], (4)
where E[·] are with respect to the stationary distribution.
Many TD solutions converge to a minimizer of MSPBE in
the limit (Tsitsiklis & Van Roy, 1997; Dann et al., 2014).
2.2

Empirical MSPBE

In practice, quantities in (4) are often unknown, and we
only have access to a finite dataset with n transitions D =
n
{(st , at , rt , st+1 )}t=1 . By replacing the unknown statistics
with their finite-sample estimates, we obtain the Empirical
MSPBE, or EM-MSPBE. Specifically, let
n

π

1 bπ
kV − ΠT π Vb π k2Ξ ,
2

X
b, 1
A
At ,
n t=1

n

X
bb , 1
bt ,
n t=1

n

X
b, 1
C
Ct ,
n t=1

(5)

Stochastic Variance Reduction Methods for Policy Evaluation

where for t = 1, . . . , n,
At , φt (φt − γφ0t )T ,

bt , rt φt ,

Ct , φt φTt .

(6)

EM-MSPBE with an optional `2 -regularization is given by:
EM-MSPBE (θ) =

ρ
1 b
kAθ − bbk2Cb−1 + kθk2 ,
2
2

(7)

where ρ ≥ 0 is a regularization factor.
Observe that (7) is a (regularized) weighted least squares
b is invertible, its optimal solution is
problem. Assuming C
b> C
b −1 A
b + ρI)−1 A
b> C
b −1bb.
θ? = (A

(8)

Computing θ? directly requires O(nd2 ) operations to form
b bb and C,
b and then O(d3 ) operations to
the matrices A,
complete the calculation. This method, known as leastsquares temporal difference or LSTD (Bradtke & Barto,
1996; Boyan, 2002), 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
(Nedić & Bertsekas, 2003). Since each rank-one update
costs O(d2 ), the total cost is O(nd2 ).
In the sequel, we develop efficient algorithms to minimize
EM-MSPBE by using stochastic variance reduction methods, which samples one (φt , φ0t ) per update without preb bb and C.
b These algorithms not only maintain
computing A,
a low O(d) per-iteration computation cost, but also attain
fast linear convergence rates with a log(1/) dependence
on the desired accuracy .

3

Saddle-Point Formulation of EM-MSPBE

Our algorithms (in Section 5) are based on the stochastic
variance reduction techniques developed for minimizing a
finite sum of convex functions, more specifically, SVRG
(Johnson & Zhang, 2013) and SAGA (Defazio et al., 2014).
They deal with problems of the form

min

x∈Rd

f (x) ,


n
1X
fi (x) ,
n i=1

(9)

where each fi is convex. We immediately notice that the
EM-MSPBE in (7) cannot be put into such a form, even
b bb and C
b have the finite-sum structhough the matrices A,
ture given in (5). Thus, extending variance reduction techniques to EM-MSPBE minimization is not straightforward.
Nevertheless, we will show that the minimizing the EMMSPBE is equivalent to solving a convex-concave saddlepoint problem which actually possesses the desired finitesum structure. To proceed, we resort to the machinery of
conjugate functions (e.g. Rockafellar, 1970, Section 12).
For a function f : Rd → R, its conjugate function f ? :

Rd → R is defined as f ? (y) , supx (y T x − f (x)). Note
that the conjugate function of 12 kxk2Cb is 12 kyk2Cb−1 , i.e.,


1
1
kyk2Cb−1 = max y T x − kxk2Cb .
x
2
2
With this relation, we can rewrite EM-MSPBE in (7) as


b − 1 kwk2 + ρ kθk2 ,
max wT (bb − Aθ)
b
C
w
2
2
so that minimizing EM-MSPBE is equivalent to solving


n
1X
min max L(θ, w) =
Lt (θ, w) ,
(10)
n t=1
θ∈Rd w∈Rd
where the Lagrangian, defined as


ρ
b − 1 kwk2 − wT bb , (11)
L(θ, w) , kθk2 − wT Aθ
b
C
2
2
may be decomposed using (5), with
1

ρ
kwk2Ct − wT bt .
Lt (θ, w) , kθk2 − wT At θ −
2
2
Therefore, minimizing the EM-MSPBE is equivalent to
solving the saddle-point problem (10), which is convex in
the primal variable θ and concave in the dual variable w.
Moreover, it has a finite-sum structure similar to (9).
Liu et al. (2015) and Valcarcel Macua et al. (2015) independently showed that the GTD2 algorithm (Sutton et al.,
2009b) is indeed a stochastic gradient method for solving
the saddle-point problem (10), although they obtained the
saddle-point formulation with different derivations. More
recently, Dai et al. (2016) 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, Lian et al. (2017) 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 4 and 5 requires the following assumption:
b has full rank, C
b is strictly positive defiAssumption 1. A
nite, and the feature vector φt is uniformly bounded.
Under mild regularity conditions (e.g., Wasserman, 2013,
b and C
b converge in probability to A
Chapter 5), we have A
and C defined in (4), 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 (e.g., Sutton et al., 2009a;b).

Stochastic Variance Reduction Methods for Policy Evaluation

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

4

A Primal-Dual Batch Gradient Method

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

 "
bT w
ρθ − A
∇θ L(θ, w)
. (12)
B (θ, w) ,
= b
b
−∇w L(θ, w)
Aθ − bb + Cw
Some notation is needed in order to characterize the convergence rate of Algorithm 1. 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:
bT C
b −1 A),
b
Lρ , λmax (ρI + A
bT C
b −1 A).
b
µρ , λmin (ρI + A

(13)
(14)

By Assumption 1, we have Lρ ≥ µρ > 0. The following
theorem is proved in Appendix B.
Theorem 1. Suppose Assumption 1 holds and let (θ? , w? )
be the (unique) solution of (10). If the step sizes are chosen
8
as σθ = 9L 1κ(C)
b and σw = 9λ
b , then the number
(C)
ρ

max

of iterations of Algorithm 1 to achieve kθ − θ? k2 + kw −
w? k2 ≤ 2 is upper bounded by
 

 1 
T b −1 b
b
b
O κ ρI + A C A · κ(C) · log
.
(15)

We assigned specific values to the step sizes σθ and σw
for clarity. In general, we can use similar step sizes while
8L
keeping their ratio roughly constant as σσwθ ≈ λ (ρC)
b ; see
min
Appendices A and B 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
bT C
b −1 A.
b The second part could
positive definiteness of A
be interpreted as transferring strong concavity in dual varib For
ables via the full-rank bi-linear coupling matrix A.

Algorithm 1 PDBG for Policy Evaluation
Inputs: initial point (θ, w), step sizes σθ and σw , and
number of epochs M .
1: for i = 1 to M do
 
  

θ
θ
σ
0
2:
←
− θ
B(θ, w)
w
w
0 σw
where B(θ, w) is computed according to (12).
3: end for

this reason, even if the saddle-point problem (10) 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 (10) using primal-dual algorithms based on proximal
mappings of the strongly convex and concave terms (e.g.,
Chambolle & Pock, 2011; Balamurugan & Bach, 2016).
The reason is that, in (10), the strong concavity of the Lagrangian with respect to the dual lies in the quadratic function (1/2)kwkCb , whose proximal mapping cannot be computed efficiently. In contrast, the PDBG algorithm only
needs its gradients.
b bb and C,
b which costs
If we pre-compute and store A,
O(nd2 ) operations, then computing the gradient operator
B(θ, w) in (12) during each iteration of PDBG costs O(d2 )
operations. Alternatively, if we do not want to store these
d × d matrices (especially if d is large), then we can compute B(θ, w) P
as finite sums on the fly. More specifically,
n
B(θ, w) = n1 t=1 Bt (θ, w), where for each t = 1, . . . , n,


ρθ − At w
Bt (θ, w) =
.
(16)
At θ − bt + C t w
Since At , bt and Ct are all rank-one matrices, as given
in (6), computing each Bt (θ, w) only requires O(d) operations. Therefore, computing B(θ, w) costs O(nd) operations as it averages Bt (θ, w) over n samples.

5

Stochastic Variance Reduction Methods

If we replace B(θ, w) in Algorithm 1 (line 2) by the
stochastic gradient Bt (θ, w) in (16), then we recover the
GTD2 algorithm of Sutton et al. (2009b), applied to a fixed
dataset, possibly with multiple passes. It has a low periteration cost but a slow, sublinear convergence rate. In
this section, we provide two stochastic variance reduction
methods and show they achieve fast linear convergence.
5.1

SVRG for policy evaluation

Algorithm 2 is adapted from the stochastic variance reduction gradient (SVRG) method (Johnson & Zhang, 2013).
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

Stochastic Variance Reduction Methods for Policy Evaluation

Algorithm 2 SVRG for Policy Evaluation
Inputs: initial point (θ, w), step sizes {σθ , σw }, number of
outer iterations M , and number of inner iterations N .
1: for m = 1 to M do
2:
Initialize (θ̃, w̃) = (θ, w) and compute B(θ̃, w̃).
3:
for j = 1 to N do
4:
Sample an index tj from {1, · · · , n} and do
5:
Compute Btj (θ, w) and Btj (θ̃, w̃).
 
  

θ
θ
σθ 0
6:
←
−
Btj (θ, w, θ̃, w̃)
w
w
0 σw
where Btj (θ, w, θ̃, w̃) is given in (17).
7:
end for
8: end for
Algorithm 3 SAGA for Policy Evaluation
Inputs: initial point (θ, w), step sizes σθ and σw , and
number of iterations M .
1: Compute each gt = Bt (θ, w)P
for t = 1, . . . , n.
n
2: Compute B = B(θ, w) = n1
t=1 gt .
3: for m = 1 to M do
4:
Sample an index tm from {1, · · · , n}.
5:
Compute htm = Btm (θ, w).
 
  

θ
θ
σθ 0
6:
←
−
(B + htm − gtm ).
w
w
0 σw
7:
B ← B + n1 (htm − gtm )
8:
gtm ← htm .
9: end for

O(nd) operations. Afterwards, the algorithm executes the
inner loop, which randomly samples an index tj and updates (θ, w) using variance-reduced stochastic gradient:
Btj(θ, w, θ̃, w̃) = Btj(θ, w) + B(θ̃, w̃) − Btj(θ̃, w̃). (17)
Here, Btj (θ, w) contains the stochastic gradients at (θ, w)
computed using the random sample with index tj , and
B(θ̃, w̃) − Btj (θ̃, w̃) is a term used to reduce the variance
in Btj (θ, w) while keeping Btj(θ, 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 4), 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 + N d). We will present the overall
complexity analysis of Algorithm 2 in Section 5.3.
5.2

SAGA for policy evaluation

The second stochastic variance reduction method for policy
evaluation is adapted from SAGA (Defazio et al., 2014);
see Algorithm 3. It uses a single loop, and maintains a single set of parameters (θ, w). Algorithm 3 starts by first

computing each component gradients gt = Bt (θ,Pw) at
n
the initial point, and also form their average B = t gt .
At each iteration, the algorithm randomly picks an index
tm ∈ {1, . . . , n} and computes the stochastic gradient
htm = Btm (θ, w). Then, it updates (θ, w) using a variance reduced stochastic gradient: B + htm − gtm , where
gtm is the previously computed stochastic gradient using
the tm -th sample (associated with certain past values of θ
and w). Afterwards, it updates the batch gradient estimate
B as B + n1 (htm − gtm ) and replaces gtm with htm .
As Algorithm 3 proceeds, different vectors gt 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 gt , for t = 1, . . . , n, to facilitate individual updates,
which will cost additional O(nd) storage. However, by exploiting the rank-one structure in (6), we only need to store
three scalars (φt − γφ0 )T θ, (φt − γφ0 )T w, and φTt w, and
form gtm on the fly using O(d) computation. Overall, each
iteration of SAGA costs O(d) operations.
5.3

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 LG based on the stochastic gradients Bt (θ, w).
Let β = σw /σθ be the ratio between the primal and dual
step-sizes, and define a pair of weighted Euclidean norms
Ω(θ, w) , (kθk2 + β −1 kwk2 )1/2 ,
Ω∗ (θ, w) , (kθk2 + βkwk2 )1/2 .
Note that Ω(·, ·) upper bounds the error in optimizing θ:
Ω(θ − θ? , w − w? ) ≥ kθ − θ? k. Therefore, any bound on
Ω(θ − θ? , w − w? ) applies automatically to kθ − θ? k.
Next, we define the parameter LG through its square:
2
Pn
1
∗
t=1 Ω Bt (θ1 , w1 ) − Bt (θ2 , w2 )
2
n
.
LG ,
sup
Ω(θ1 − θ2 , w1 − w2 )2
θ1 ,w1 ,θ2 ,w2
This definition is similar to the smoothness constant L̄ used
in Balamurugan & Bach (2016) except that we used the
step-size ratio β rather than the strong convexity and concavity parameters of the Lagrangian to define Ω and Ω∗ .1
Substituting the definition of Bt (θ, w) in (16), we have


√
n
1X T
ρI − βATt
L2G =
.
Gt Gt , where Gt , √
βAt
βCt
n t=1
(18)
With the above definitions, we characterize the convergence of Ω(θm − θ? , wm − w? ), where (θ? , w? ) is the solution of (10), and (θm , wm ) is the output of the algorithms
1

Since our saddle-point problem is not necessarily strongly
convex in θ (when ρ = 0), we could not define Ω and Ω∗ in the
same way as Balamurugan & Bach (2016).

Stochastic Variance Reduction Methods for Policy Evaluation

after the m-th iteration. For SVRG, it is the m-th outer
iteration in Algorithm 2. The following two theorems are
proved in Appendices C and D, respectively.
Theorem 2 (Convergence rate of SVRG). Suppose Asµρ
sumption 1 holds. If we choose σθ = 48κ(C)L
b 2 , σw =
G

b 2
51κ2(C)L
G
,
µ2ρ

8Lρ
b σθ ,
λmin (C)

N =
in (13) and (14), then

where Lρ and µρ are defined


  4 m
E Ω(θm −θ? , wm −w? )2 ≤
Ω(θ0 −θ? , w0 −w? )2 .
5

The overall computational
cost for reaching E Ω(θm −

θ? , wm − w? ) ≤  is upper bounded by


 1 
b 2
κ(C)L
G
. (19)
O
n+
d log
bT C
b −1 A)
b

λ2min (ρI + A
Theorem 3 (Convergence rate of SAGA). Suppose Asµρ
sumption 1 holds. If we choose σθ = 3 8κ2 (C)L
b 2 +nµ2 )
(
ρ
G
8L
and σw = λ (ρC)
σ
in
Algorithm
3,
then
b θ
min



E Ω(θm −θ? , wm −w? )2 ≤ 2(1−ρ)m Ω(θ0 −θ? , w0 −w? )2 ,
µ2

ρ
where ρ ≥ 9 8κ2 (C)L
b 2 +nµ2 ) . The total cost to achieve
(
ρ
G 

E Ω(θm − θ? , wm − w? ) ≤  has the same bound in (19).

Similar to our PDBG results in (15), both the SVRG and
SAGA algorithms for policy evaluation enjoy linear convergence even if there is no strong convexity in the saddlepoint problem (10) (i.e., when ρ = 0). This is mainly due to
b when C
b is positiveb −1 A
bT C
the positive definiteness of A
b is full-rank. In contrast, the linear converdefinite and A
gence of SVRG and SAGA in Balamurugan & Bach (2016)
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)kwkCb , which does not admit an efficient proximal mapping as required by the proximal versions of SVRG and SAGA in Balamurugan &
Bach (2016). Our algorithms only require computing the
stochastic gradients of this function, which is easy to do
due to its finite sum structure.
Balamurugan & Bach (2016) also proposed accelerated
variants of SVRG and SAGA using the “catalyst” framework of Lin et al. (2015). Such extensions can be done
similarly for the three algorithms presented in this paper,
and we omit the details due to space limit.

6

Comparison of Different Algorithms

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

Table 1. Complexity of different policy evaluation algorithms. In
the table, d is feature dimension, n is dataset size, κ , κ(ρI +
bT C
b −1 A);
b κG , LG /λmin (ρI + A
bT C
b −1 A);
b and κ0 is a conA
dition number related to GTD2.
Algorithm
SVRG / SAGA
GTD2
PDBG-(I)
PDBG-(II)
LSTD

Total Complexity




b 2
κ(C)κ
G
O nd · 1 +
· log 1/
n
O (d · κ0 /)


b · log(1/)
O nd · κ(C)κ


b · log(1/)
O nd2 + d2 κ(C)κ


O nd2 or O nd2 + d3

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., Lazaric et al. (2010)), leading to an
O(κ0 nd) overall complexity for GTD2, where κ0 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
b bb
Section 4. PDBG-(II) first pre-computes the matrices A,
2
b
and C using O(nd ) operations, then computes the batch
gradient at each iteration with O(d2 ) 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(nd2 ) complexity if rankone 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 >
b 2
κ(C)κ

G
(1 +
) log 1 , This condition is easy to satisfy,
n
when n is very large. On the other hand, SVRG and SAGA
algorithms are more
 efficient than
 PDBG-(I) if n is large,
b 2 κ(C)κ
b − 1 , where κ and κG are desay n > κ(C)κ
G
scribed in the caption of Table 1.

There are other algorithms whose complexity scales
linearly with n and d, including iLSTD (Geramifard
et al., 2007), and TDC (Sutton et al., 2009b), fLSTDSA (Prashanth et al., 2014), and the more recent algorithms
of Wang et al. (2016) and Dai et al. (2016). However, their

Stochastic Variance Reduction Methods for Policy Evaluation

convergence is slow: the number of iterations required to
reach a desired accuracy  grows as 1/ or worse. The
CTD algorithm (Korda & Prashanth, 2015) 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 7.1).

7

Extensions

It is possible to extend our approach to accelerate optimization of other objectives such as MSBE and NEU (Dann
et al., 2014). In this section, we briefly describe two extensions of the algorithms developed earlier.
7.1

Off-policy learning

In some cases, we may want to estimate the value function
of a policy π from a set of data D generated by a different “behavior” policy πb . This is called off-policy learning (Sutton & Barto, 1998, Chapter 8).
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 (Tsitsiklis & Van Roy, 1997), our gradient-based algorithms remain convergent with the same (fast) convergence rate.
Consider the RL framework outlined in Section 2. For each
state-action pair (st , at ) such that πb (at |st ) > 0, we define
the importance ratio, ρt , π(at |st )/πb (at |st ). The EMMSPBE for off-policy learning has the same expression as
in (7) except that At , bt and Ct are modified by the weight
factor ρt , as listed in Table 2; see also Liu et al. (2015,
Eqn 6) for a related discussion.) Algorithms 1–3 remain the
same for the off-policy case after At , bt and Ct are modified
correspondingly.
7.2

Learning with eligibility traces

Eligibility traces are a useful technique to trade off bias and
variance in TD learning (Singh & Sutton, 1996; Kearns &
Singh, 2000). When they are used, we can pre-compute zt
in Table 2 before running our new algorithms. Note that
EM-MSPBE with eligibility traces has the same form of
(7), with At , bt and Ct defined differently according to the
last row of Table 2. At the m-th step of the learning process, the algorithm randomly samples ztm , φtm , φ0tm and
rtm from the fixed dataset and computes the corresponding
stochastic gradients, where the index tm is uniformly distributed over {1, . . . , n} and are independent for different
values of m. Algorithms 1–3 immediately work for this
case, enjoying a similar linear convergence rate and a com-

Table 2. Expressions of At , bt and Ct for different cases of policy
Here, ρt , π(at |st )/πb (at |st ); and zt ,
Pt evaluation.
t−i
φi , where λ ≥ 0 is a given parameter.
i=1 (λγ)
On-policy
Off-policy
Eligibility trace

At
φt (φt − γφ0t )>
ρt φt (φt − γφ0t )>
zt (φt − γφ0t )>

bt
rt φt
ρt r t φ t
r t zt

Ct
φt φ>
t
φt φ>
t
φt φ>
t

putation complexity linear in n and d. We need additional
O(nd) operations to pre-compute zt recursively and an additional O(nd) storage for zt . However, it does not change
the order of the total complexity for SVRG/SAGA.

8

Experiments

In this section, we compare the following algorithms on
two benchmark problems: (i) PDBG (Algorithm 1); (ii)
GTD2 with samples drawn randomly with replacement
from a dataset; (iii) TD: the fLSTD-SA algorithm of
Prashanth et al. (2014); (iv) SVRG (Algorithm 2); and
(v) SAGA (Algorithm 3). 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
 −1
1
10 , 10−2 , . . . , 10−6 L κ(
b and σw is chosen from
C)
ρ

1
−1
−2
1, 10 , 10
b . We only report the results of
λmax (C)
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 (Dann et al., 2014). The transition probabilities are defined as P (s0 |a, s) ∝ pass0 +10−5 ,
where pass0 ∼ 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. 1 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 (Sutton &
Barto, 1998, Chapter 8). 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. 2 and 3 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

0

0

−2

−2

−4
−6
PDBG
GTD2
TD
SVRG
SAGA

−8
−10
−12
0

1000

2000
Epochs

3000

−4
−6
−8

−12
0

4000

PDBG
GTD2
SVRG
SAGA

−10

1000

2000
Epochs

3000

log10(Objective Value)

0
−2

log10(Objective Value)

log10(Objective Value)

Stochastic Variance Reduction Methods for Policy Evaluation

−4
−6
−8

−12
0

4000

1000

2000
Epochs

3000

4000



b> C
b −1 A
b
(c) ρ = λmax A

q
b> C
b −1 A)
b
(b) ρ = λmax (A

(a) ρ = 0

PDBG
GTD2
SVRG
SAGA

−10

0

0

−2

−2

−2

−4
−6
PDBG
GTD2
TD
SVRG
SAGA

−8
−10
−12
0

500

1000

1500

−4
−6
−8
−10
−12
0

2000

PDBG
GTD2
SVRG
SAGA
500

Epochs

1000

1500

log10(Objective Value)

0

log10(Objective Value)

log10(Objective Value)

Figure 1. Random MDP with s = 400, d = 200, and n = 20000.

−4
−6
−8
PDBG
GTD2
SVRG
SAGA

−10
−12
0

2000

500

Epochs

b
b −1 A)
b> C
(b) ρ = 0.01λmax (A

(a) ρ = 0

1000

1500

2000

Epochs

b
b −1 A)
b> C
(c) ρ = λmax (A

0

0

−2

−2

−2

−4
−6
PDBG
GTD2
TD
SVRG
SAGA

−8
−10
−12
0

500

1000

1500

2000

−4
−6
−8
PDBG
GTD2
SVRG
SAGA

−10
−12
0

500

Epochs

(a) ρ = 0

1000

1500

2000

log10(Objective Value)

0

log10(Objective Value)

log10(Objective Value)

Figure 2. Mountain Car Data Set with d = 300 and n = 5000.

−4
−6
−8
PDBG
GTD2
SVRG
SAGA

−10
−12
0

500

Epochs

b> C
b −1 A)
b
(b) ρ = 0.01λmax (A

1000
Epochs

1500

2000

b> C
b −1 A)
b
(c) ρ = λmax (A

Figure 3. Mountain Car Data Set with d = 300 and n = 20000.

theoretical finding in Table 1 that SVRG/SAGA need fewer
epochs for large n.

9

Conclusions

In this paper, we reformulated the EM-MSPBE minimization problem in policy evaluation into an empirical saddlepoint 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 (Bhatnagar et al., 2009), especially with
deep neural networks. Moreover, it remains an important
open problem how to apply stochastic variance reduction
techniques to policy optimization.

Stochastic Variance Reduction Methods for Policy Evaluation

References
Balamurugan, P and Bach, Francis. Stochastic variance reduction methods for saddle-point problems. In Advances
in Neural Information Processing Systems 29, pp. 1416–
1424, 2016.
Bertsekas, Dimitri P and Tsitsiklis, John N. Neurodynamic programming: An overview. In Decision and
Control, 1995., Proceedings of the 34th IEEE Conference on, volume 1, pp. 560–564. IEEE, 1995.
Bhatnagar, Shalabh, Precup, Doina, Silver, David, Sutton, Richard S, Maei, Hamid R, and Szepesvári, Csaba.
Convergent temporal-difference learning with arbitrary
smooth function approximation. In Advances in Neural
Information Processing Systems, pp. 1204–1212, 2009.
Boyan, Justin A. Technical update: Least-squares temporal difference learning. Machine Learning, 49:233–246,
2002.
Bradtke, Steven J and Barto, Andrew G. Linear leastsquares algorithms for temporal difference learning. Machine Learning, 22:33–57, 1996.
Chambolle, Antonin and Pock, Thomas. A first-order
primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging
and Vision, 40(1):120–145, 2011.
Dai, Bo, He, Niao, Pan, Yunpeng, Boots, Byron, and Song,
Le. Learning from conditional distributions via dual embeddings. arXiv:1607.04579, 2016.
Dann, Christoph, Neumann, Gerhard, and Peters, Jan. Policy evaluation with temporal differences: a survey and
comparison. Journal of Machine Learning Research, 15
(1):809–883, 2014.
Defazio, Aaron, Bach, Francis, and Lacoste-Julien, Simon.
SAGA: A fast incremental gradient method with support
for non-strongly convex composite objectives. In Advances in Neural Information Processing Systems, pp.
1646–1654, 2014.
Geramifard, Alborz, Bowling, Michael H., Zinkevich,
Martin, and Sutton, Richard S. iLSTD: Eligibility traces
and convergence analysis. In Advances in Neural Information Processing Systems 19, pp. 441–448, 2007.
Gohberg, Israel, Lancaster, Peter, and Rodman, Leiba. Indefinite linear algebra and applications. Springer Science & Business Media, 2006.
Johnson, Rie and Zhang, Tong. Accelerating stochastic
gradient descent using predictive variance reduction. In
Advances in Neural Information Processing Systems, pp.
315–323, 2013.

Kearns, Michael J. and Singh, Satinder P. “Bias-variance”
error bounds for temporal difference updates. In Proceedings of the Thirteenth Annual Conference on Computational Learning Theory (COLT-00), pp. 142–147,
2000.
Korda, Nathaniel and Prashanth, L.A. On TD(0) with function approximation: Concentration bounds and a centered variant with exponential convergence. In Proceedings of the Thirty-Second International Conference on
Machine Learning (ICML-15), pp. 626–634, 2015.
Lagoudakis, Michail G and Parr, Ronald. Least-squares
policy iteration. Journal of Machine Learning Research,
4(Dec):1107–1149, 2003.
Lange, Sascha, Gabel, Thomas, and Riedmiller, Martin.
Batch reinforcement learning. In Wiering, Marco and
van Otterlo, Martijn (eds.), Reinforcement Learning:
State of the Art, pp. 45–73. Springer Verlag, 2011.
Lazaric, Alessandro, Ghavamzadeh, Mohammad, and
Munos, Rémi. Finite-sample analysis of LSTD. In
Proceedings of the Twenty-Seventh International Conference on Machine Learning, pp. 615–622, 2010.
Lian, Xiangru, Wang, Mengdi, and Liu, Ji. Finite-sum
composition optimization via variance reduced gradient
descent. In Proceedings of Artificial Intelligence and
Statistics Conference (AISTATS), pp. 1159–1167, 2017.
Liesen, Jörg and Parlett, Beresford N. On nonsymmetric
saddle point matrices that allow conjugate gradient iterations. Numerische Mathematik, 108(4):605–624, 2008.
Lin, Hongzhou, Mairal, Julien, and Harchaoui, Zaid.
A universal catalyst for first-order optimization. In
Advances in Neural Information Processing Systems
(NIPS) 28, pp. 3384–3392, 2015.
Lin, Long-Ji. Self-improving reactive agents based on reinforcement learning, planning and teaching. Machine
Learning, 8(3–4):293–321, 1992.
Liu, Bo, Liu, Ji, Ghavamzadeh, Mohammad, Mahadevan,
Sridhar, and Petrik, Marek. Finite-sample analysis of
proximal gradient TD algorithms. In Proc. The 31st
Conf. Uncertainty in Artificial Intelligence, Amsterdam,
Netherlands, 2015.
Nedić, A. and Bertsekas, Dimitri P. Least squares policy evaluation algorithms with linear function approximation. Discrete Event Dynamics Systems: Theory and
Applications, 13(1):79–110, 2003.
Prashanth, LA, Korda, Nathaniel, and Munos, Rémi. Fast
LSTD using stochastic approximation: Finite time analysis and application to traffic control. In Joint European

Stochastic Variance Reduction Methods for Policy Evaluation

Conference on Machine Learning and Knowledge Discovery in Databases, pp. 66–81. Springer, 2014.
Precup, Doina, Sutton, Richard S., and Dasgupta, Sanjoy.
Off-policy temporal-difference learning with funtion approximation. In Proceedings of the Eighteenth Conference on Machine Learning (ICML-01), pp. 417–424,
2001.
Puterman, Martin L. Markov Decision Processes: Discrete
Stochastic Dynamic Programming. John Wiley & Sons,
2005.
Rockafellar, R. Tyrrell. Convex Analysis. Princeton University Press, 1970.
Shen, Shu-Qian, Huang, Ting-Zhu, and Cheng, GuangHui. A condition for the nonsymmetric saddle point
matrix being diagonalizable and having real and positive eigenvalues. Journal of Computational and Applied
Mathematics, 220(1):8–12, 2008.
Singh, Satinder P. and Sutton, Richard S. Reinforcement learning with replacing eligibility traces. Machine
Learning, 22(1–3):123–158, 1996.
Sutton, Richard S and Barto, Andrew G. Reinforcement
Learning: An Introduction. MIT Press, Cambridge, MA,
1998.
Sutton, Richard S, Maei, Hamid R, and Szepesvári, Csaba.
A convergent o(n) temporal-difference algorithm for
off-policy learning with linear function approximation.
In Advances in neural information processing systems,
pp. 1609–1616, 2009a.
Sutton, Richard S, Maei, Hamid Reza, Precup, Doina,
Bhatnagar, Shalabh, Silver, David, Szepesvári, Csaba,
and Wiewiora, Eric. Fast gradient-descent methods
for temporal-difference learning with linear function approximation. In Proceedings of the 26th Annual International Conference on Machine Learning, pp. 993–1000.
ACM, 2009b.
Tsitsiklis, John N. and Van Roy, Benjamin. An analysis
of temporal-difference learning with function approximation. IEEE Transactions on Automatic Control, 42:
674–690, 1997.
Valcarcel Macua, Sergio, Chen, Jianshu, Zazo, Santiago,
and Sayed, Ali H. Distributed policy evaluation under
multiple behavior strategies. Automatic Control, IEEE
Transactions on, 60(5):1260–1274, 2015.
Wang, Mengdi, Liu, Ji, and Fang, Ethan. Accelerating
stochastic composition optimization. In Advances in
Neural Information Processing Systems (NIPS) 29, pp.
1714–1722, 2016.

Wasserman, Larry. All of Statistics: A Concise Course in
Statistical Inference. Springer Science & Business Media, 2013.

Stochastic Variance Reduction Methods for Policy Evaluation

A

Eigen-analysis of G

A.2

In this section, we give a thorough analysis of the spectral
properties of the matrix
"
#
bT
ρI
−β 1/2 A
G = 1/2 b
,
(20)
b
β A
βC
which is critical in analyzing the convergence of the PDBG,
SAGA and SVRG algorithms for policy evaluation. Here
β = σw /σθ is the ratio between the dual and primal step
sizes in these algorithms. For convenience, we use the following notation:
bT C
b −1 A),
b
L , λmax (A
bT C
b −1 A).
b
µ , λmin (A

Analysis of eigenvectors

If the matrix G is diagonalizable, then it can be written as
G = QΛQ−1 ,
where Λ is a diagonal matrix whose diagonal entries are the
eigenvalues of G, and Q consists of it eigenvectors (each
with unit norm) as columns. Our goal here is to bound
κ(Q), the condition number of the matrix Q. Our analysis is inspired by Liesen & Parlett (2008). The core is the
following fundamental result from linear algebra.
Theorem 4 (Theorem 5.1.1 of Gohberg et al. (2006)). Suppose G is diagonalizable. If H is a symmetric positive definite matrix and HG is symmetric, then there exist a complete set of eigenvectors of G, such that they are orthonormal with respect to the inner product induced by H:

Under Assumption 1, they are well defined and we have
L ≥ µ > 0.
A.1

Diagonalizability of G

First, we examine the condition of β that ensures the diagonalizability of the matrix G. We cite the following result
from (Shen et al., 2008).
Lemma 1. Consider the matrix A defined as


A −B >
A=
,
B
C

(21)

where A  0, C  0, and B is full rank. Let τ = λmin√
(C),
δ = λmax (A) and σ = λmax (B > C −1 B). If τ > δ +2 τ σ
holds, then A is diagonalizable with all its eigenvalues real
and positive.
Applying this lemma to the matrix G in (20), we have
b = βλmin (C),
b
τ = λmin (β C)
δ = λmax (ρI) = ρ,

b> (β C)
b −1 β 1/2 A
b = λmax (A
b> C
b −1 A).
b
σ = λmax β 1/2 A
√

The condition τ > δ + 2 τ σ translates into
q
b > ρ + 2 βλmin (C)λ
b max (A
b> C
b −1 A),
b
βλmin (C)
which can be solved as
q
q
b> C
b −1 A)+
b> C
b −1 A)
b
b
p
λmax (A
ρ+λmax (A
q
β>
.
b
λmin (C)
In the rest of our discussion, we choose β to be


b> C
b −1 A
b
8 ρ + λmax A
8(ρ + L)
β=
=
,
b
b
λmin (C)
λmin (C)
which satisfies the inequality above.

(22)

Q> HQ = I.

(23)

If H satisfies the conditions in Theorem 4, then we have
H = Q−> Q−1 , which implies κ(H) = κ2 (Q). Therefore,
in order to bound κ(Q), we only need to find such an H
and analyze its conditioning. To this end, we consider the
matrix of the following form:
"
√ > #
b
βA
(δ − ρ)I
√
(24)
H=
b
b − δI .
βA
βC
It is straightforward to check that HG is a symmetric matrix. The following lemma states the conditions for H being positive definite.
b − δI − β A
b b>  0,
Lemma 2. If δ − ρ > 0 and β C
δ−ρ A
then H is positive definite.
Proof. The matrix
decomposition:
"
I
H = √β b
A
δ−ρ

H in (24) admits the following Schur
0
I

#
(δ − ρ)I

"
I
S 0

#
√
β b>
δ−ρ A
,
I

b − δI − β A
b b>
where S = β C
δ−ρ A . Thus H is congruence to
the block diagonal matrix in the middle, which is positive
definite under the specified conditions. Therefore, the matrix H is positive definite under the same conditions.
In addition to the choice of β in (22), we choose δ to be
δ = 4(ρ + L).

(25)

It is not hard to verify that this choice ensures δ−ρ > 0 and
b − δI − β A
b b>
βC
δ−ρ A  0 so that H is positive definite. We
now derive an upper bound on the condition number of H.
Let λ be an eigenvalue of H and [xT y T ]T be its associated
eigenvector, where kxk2 + kyk2 > 0. Then it holds that
p T
b y = λx,
(δ − ρ)x + β A
(26)

Stochastic Variance Reduction Methods for Policy Evaluation

p
b + (β C
b − δI)y = λy.
β Ax

(b)

(27)

≥ δ−ρ−L−

From (26), we have



3ρ + 4L
= (ρ + L) 3 −
2(ρ + L)

(c)

√

x=

δ(δ − ρ)
b
βλmin (C)

β
bT y.
A
λ−δ+ρ

(28)

Note that λ − δ + ρ 6= 0 because if λ − δ + ρ = 0 we have
bT y = 0 so that y = 0 since A
b is full rank. With y = 0 in
A
b = 0 so that x = 0, which contradicts
(27), we will have Ax
the assumption that kxk2 + kyk2 > 0.
Substituting (28) into (27) and multiplying both sides with
y T , we obtain the following equation after some algebra
2

λ − pλ + q = 0,

(29)

≥ ρ + L,

(32)

where step (a) uses the fact that both the numerator and
denominator are positive, step (b) uses the fact


T b bT
bT C
b −1 A
b ≥ y AA y ,
L , λmax A
b
y T Cy
and step (c) substitutes the expressions of δ and β. Therefore, we can upper bound the condition number of H, and
thus that of Q, as follows:

where
κ2 (Q) = κ(H) ≤

b − δI)y
y T (β C
,
p,δ−ρ+
kyk2
b − δI)y
bA
bT y
y T (β C
yT A
q , (δ − ρ)
−β
.
2
kyk
kyk2

A.3

We can verify that both p and q are positive with our choice
of δ and β. The roots of the quadratic equation in (29) are
given by
p
p ± p2 − 4q
λ=
.
(30)
2
Therefore, we can upper bound the largest eigenvalue as
p
p + p2 − 4q
λmax (H) ≤
2
b
y T Cy
≤p=δ−ρ−δ+β
kyk2
b
≤ −ρ + βλmax (C)
8(ρ + L)
b
λ
(C)
= −ρ +
b max
λmin (C)

(33)

Analysis of eigenvalues

>
Suppose λ is an eigenvalue of G and let ξ > , η > be its
corresponding eigenvector. By definition, we have
 
 
ξ
ξ
G
=λ
,
η
η
which is equivalent to the following two equations:
p >
b η = λξ,
ρξ − β A
p
b + β Cη
b = λη.
β Aξ
Solve ξ in the first equation in terms of η, then plug into the
second equation, we obtain:
b + β(A
bA
b> η + ρCη)
b = 0.
λ2 η − λ(ρη + β Cη)
2

Now left multiply η > , then divide by the kηk2 , we have:
λ2 − pλ + q = 0.

b
≤ 8(ρ + L)κ(C).

(31)

Likewise, we can lower bound the smallest eigenvalue:
p
p − p2 − 4q
p − p + 2q/p
q
λmin (H) ≥
≥
=
2
2
p


b
bA
bT y
y T Cy
yT A
β (δ − ρ) kyk2 − kyk2
− δ(δ − ρ)
=
T
b
−ρ + β ykykCy
2


T b
T b bT
y Cy
y AA y
− δ(δ − ρ)
(a) β (δ − ρ) kyk2 −
kyk2
≥
T b
β ykykCy
2
bA
bT y δ(δ − ρ)
yT A
=δ−ρ−
−
·
b
β
y T Cy

b
8(ρ + L)κ(C)
b
= 8κ(C).
ρ+L

1
b
y T Cy
kyk2

where p and q are defined as
b
η > Cη
,
kηk2
bA
b> η
ηT A

p,ρ+β
q,β

kηk2

b
η T Cη
+ρ
kηk2

Therefore the eigenvalues of G satisfy:
p
p ± p2 − 4q
.
λ=
2

!
.

(34)

(35)

Recall that our choice of β ensures that G is diagonalizable and has positive real eigenvalues. Indeed, we can verify that the diagonalization condition guarantees p2 ≥ 4q

Stochastic Variance Reduction Methods for Policy Evaluation

so that all eigenvalues are real and positive. Now we can
obtain upper and lower bounds based on (35). For upper
bound, notice that
b
λmax (G) ≤ p ≤ ρ + βλmax (C)
=ρ+

The gradients of the Lagrangian with respect to θ and w,
respectively, are
b> w
∇θ L (θ, w) = ρθ − A
b − Cw
b + bb.
∇w L (θ, w) = −Aθ

8(ρ + L)
b
λ
(C)
b max
λmin (C

b
= ρ + 8(ρ + L)κ(C)

b ρ+L
≤ 9κ(C)

b max ρI + A
bT C
b −1 A
b .
= 9κ(C)λ

(36)

For lower bound, notice that
p
p − p2 − 4q
p − p + 2q/p
λmin (G) ≥
≥
= q/p
2

 T b2bT
A η
+ρ
β η ηA
T Cη
b
=
kηk2
ρ ηT Cη
b +β
(a)

b
β(ρ + µ)
βλmin (C)(ρ
+ µ)
=
b +β
b
ρ/λmin (C)
ρ + βλmin (C)
(b) 8(ρ + L)(ρ + µ)
=
ρ + 8(ρ + L)
8
≥ (ρ + µ)
9
8
bT C
b −1 A))
b
= (ρ + λmin (A
9
8
bT C
b −1 A),
b
= λmin (ρI + A
(37)
9
≥

where the second inequality is by the concavity property of
the square root function, step (a) used the fact


T b bT
bT C
b −1 A
b ≤ y AA y ,
µ , λmin A
b
y T Cy
and step (b) substitutes the expressions of β.
Since G is not a normal matrix, we cannot use their eigenvalue bounds to bound its condition number κ(G).

B



b> θ? .
b −1 bb − A
w? = C

Linear convergence of PDBG

The first-order optimality condition is obtained by setting
them to zero, which is satisfied by (θ? , w? ):
"
#   
b> θ?
0
ρI −A
= b .
(39)
b
b
w?
b
A
C
The PDBG method in Algorithm 1 takes the following iteration:

   

θm+1
θm
σθ 0
=
−
B(θm , wm ),
wm+1
wm
0 σw
where
 "
∇θ L(θ, w)
ρI
B(θ, w) =
= b
−∇w L(θ, w)
A


b>
−A
b
C

#   
0
θ
− b .
w
b

Letting β = σw /σθ , we have
"
#    !

  
b> θm
0
θm+1
θm
ρI −A
=
− σθ
− b .
b
b
wm+1
wm
w
βb
βA βC
m
Subtracting both sides of the above recursion by (θ? , w? )
and using (39), we obtain
"
#

 


bT θm −θ?
θm+1 − θ?
θm − θ?
ρI −A
=
−σθ
b βC
b wm −w? .
wm+1 −w?
wm −w?
βA
We analyze the convergence of the algorithms by examining the differences between the current parameters to the
optimal solution. More specifically, we define a scaled
residue vector


θm − θ?
,
(40)
∆m , √1
(wm − w? )
β

Recall the saddle-point problem we need to solve:
which obeys the following iteration:
min max L(θ, w),
θ

w

∆m+1 = (I − σθ G) ∆m ,

where the Lagrangian is defined as
ρ
b − 1 w> Cw
b + bb> w.
L(θ, w) = kθk2 − w> Aθ
2
2

(38)

b is positive definite and A
b has full
Our assumption is that C
rank. The optimal solution can be expressed as

−1
b> C
b −1 A
b + ρI
b> C
b −1bb,
θ? = A
A

(41)

where G is exactly the matrix defined in (20). As analyzed in Section A.1, if we choose β sufficiently large,
such as in (22), then G is diagonalizable with all its eigenvalues real and positive. In this case, we let Q be the
matrix of eigenvectors in the eigenvalue decomposition
G = QΛQ−1 , and use the potential function
Pm , Q−1 ∆m

2
2

Stochastic Variance Reduction Methods for Policy Evaluation

in our convergence analysis. We can bound the usual Euclidean distance by Pm as
2

2

kθm − θ? k + kwm − w? k ≤ (1 +

2
β)σmax
(Q)Pm .

If we have linear convergence in Pm , then the extra factor
2
(1 + β)σmax
(Q) will appear inside a logarithmic term.
Remark: This potential function has an intrinsic geometric
interpretation. We can view column vectors of Q−1 a basis
for the vector space, which is not orthogonal. Our goal is to
show that in this coordinate system, the distance to optimal
solution shrinks at every iteration.

C

Analysis of SVRG

Here we establish the linear convergence of the SVRG algorithm for policy evaluation described in Algorithm 2.
b bb and C:
b
Recall the finite sum structure in A,
n

X
b= 1
A
At ,
n t=1

n

= Q

−1

QQ

−1

− σθ QΛQ

= (I − σθ Λ) Q−1 ∆m





2
∆m 2

Bt (θ, w) =

2
2

2

≤ kI − σθ Λk2 Q−1 ∆m

2
2

2

= kI − σθ Λk2 Pm

(42)

The inequality above uses sub-multiplicity of spectral
norm. We choose σθ to be
σθ =

1
λmax (Λ)

=

1
λmax (G)

,

(43)

Since all eigenvalues of G are real and positive, we have


2

kI − σθ Λk =
≤

λmin (G)
1−
λmax (G)

2

1
8
·
1−
b
bT C
b −1 A)
b
81 κ(C)κ(ρI
+A

!2
,

where we used the bounds on the eigenvalues λmax (G) and
λmin (G) in (36) and (37) respectively. Therefore, we can
achieve an -close solution with

 
b
bT C
b −1 A)
b log P0
m = O κ(C)κ(ρI
+A

iterations of the PDBG algorithm.
In order to minimize kI − σθ Λk, we can choose
σθ =

1X
Bt (θ, w),
n t=1

where

2
2

= Q−1 (I − σθ G) ∆m
−1

B(θ, w) =

2
2

n

X
b= 1
C
Ct .
n t=1

This structure carries over to the Lagrangian L(θ, w) as
well as the gradient operator B(θ, w), so we have

We proceed to bound the growth of Pm :
Pm+1 = Q−1 ∆m+1

n

X
bb = 1
bt ,
n t=1

2
,
λmax (G) + λmin (G)

which results in kI − σθ Λk = 1 − 2/(1 + κ(Λ)) instead of
1−1/κ(Λ). The resulting complexity stays the same order.
The step sizes stated in Theorem 1 is obtained by replacing
λmax in (43) with its upper bound in (36) and setting σw
through the ratio β = σw /σθ as in (22).

ρI
At

−A>
t
Ct

   
θ
0
−
.
w
bt

(44)

Algorithm 2 has both an outer loop and an inner loop. We
use the index m for the outer iteration and j for the inner
iteration. Fixing the outer loop index m, we look at the
inner loop of Algorithm 2. Similar to full gradient method,
we first simplify the dynamics of SVRG.

 
 
 
θm,j+1
θ
σ
= m,j − θ
× B(θm−1 , wm−1 )
wm,j+1
wm,j
σw

+ Btj (θm,j , wm,j ) − Bt (θm−1 , wm−1 )

 

θ
σ
= m,j − θ
wm,j
σw
"
#
  
>
b
0
ρI −A
θm−1
×
− b
b
b
w
b
A
C
m−1






θm,j
0
ρI −A>
t
+
−
wm,j
bt
At
Ct





!
θm−1
0
ρI −A>
t
−
+
.
wm−1
bt
At
Ct
Subtracting (θ? , w? ) from both sides and using the optimality condition (39), we have

 
 

θm,j+1 − θ?
θm,j − θ?
σθ
=
−
wm,j+1 − w?
wm,j − w?
σw
"
#

b>
ρI −A
θm−1 − θ?
×
b
b
wm−1 − w?
A
C




θm,j − θ?
ρI −A>
t
+
wm,j − w?
At
Ct



!
θm−1 − θ?
ρI −A>
t
−
.
wm−1 − w?
At
Ct

Stochastic Variance Reduction Methods for Policy Evaluation

Multiplying
√ both sides of the above recursion by
diag(I, 1/ βI), and using a residue vector ∆m,j defined
similarly as in (40), we obtain
∆m,j+1 = ∆m,j − σθ (G∆m−1 + Gtj ∆m,j − Gtj ∆m−1 )
= (I − σθ G) ∆m,j

+ σθ G − Gtj (∆m,j − ∆m−1 ) ,

(45)

where in the second inequality we used the definition of L2G
in (18), i.e., L2G = kE[GTtj Gtj ]k. In addition, we have

where Gtj is defined in (18).
For SVRG, we use the following potential functions to facilitate our analysis:
h
i
2
Pm , E Q−1 ∆m
,
(46)
h
i
2
Pm,j , E Q−1 ∆m,j
.
(47)
Unlike the analysis for the batch gradient methods, the nonorthogonality of the eigenvectors will lead to additional dependency of the iteration complexity on the condition number of Q, for which we give a bound in (33).
Multiplying both sides of Eqn. (45) by Q−1 , taking squared
2-norm and taking expectation, we obtain
h

Pm,j+1 = E Q−1 (I − σθ G) ∆m,j

 2i
+ σθ G − Gtj (∆m,j − ∆m−1 )
h
i
(a)
2
= E (I − σθ Λ) Q−1 ∆m,j
h
i

2
+ σθ2 E Q−1 G−Gtj (∆m,j −∆m−1 )
h
i
(b)
2
2
≤ kI − σθ Λk E Q−1 ∆m,j
h
i
2
+ σθ2 E Q−1 Gtj (∆m,j − ∆m−1 )
(c)

2

= kI − σθ Λk Pm,j
h
+ σθ2 E Q−1 Gtj (∆m,j − ∆m−1 )

2

i
. (48)

where step (a) used the facts that Gtj is independent of
∆m,j and ∆m−1 and E[Gtj ] = G so the cross terms are
zero, step (b) used again the same independence and that
the variance of a random variable is less than its second
moment, and step (c) used the definition of Pm,j in (47).
To bound the last term in the above inequality, we use the
simple notation δ = ∆m,j − ∆m−1 and have
Q−1 Gtj δ

2

= δ T GTtj Q−T Q−1 Gtj δ
≤ λmax (Q−T Q−1 )δ T GTtj Gtj δ.

Therefore, we can bound the expectation as

E Q−1 Gtj δ

2



≤λmax (Q−T Q−1 )E δ T GTtj Gtj δ



=λmax (Q−T Q−1 )E δ T E[GTtj Gtj ]δ


≤λmax (Q−T Q−1 )L2G E δ T δ


=λmax (Q−T Q−1 )L2G E δ T Q−T QT QQ−1 δ


=λmax (Q−T Q−1 )λmax (QT Q)L2G E δ T Q−T Q−1 δ


≤κ(Q)2 L2G E kQ−1 δk2 ,
(49)




2
E kQ−1 δk2 =E Q−1 (∆m,j − ∆m−1 )


2
≤2 E Q−1 ∆m,j
+ 2 E Q−1 ∆m−1

2

= 2Pm,j + 2Pm−1 .
Then it follows from (48) that
Pm,j+1 ≤kI − σθ Λk2 Pm,j
+ 2σθ2 κ2 (Q)L2G (Pm,j + Pm−1 ).
Next, let λmax and λmin denote the largest and smallest
diagonal elements of Λ (eigenvalues of G), respectively.
Then we have

2
kI − σθ Λk = max (1 − σθ λmin )2 , (1 − σθ λmin )2
≤ 1 − 2σθ λmin + σθ2 λ2max
≤ 1 − 2σθ λmin + σθ2 κ2 (Q)L2G ,
where the last inequality uses the relation
λ2max ≤ kGk2 = kEGt k2 ≤ kEGTt Gt k = L2G ≤ κ2 (Q)L2G .
It follows that

Pm,j+1 ≤ 1 − 2σθ λmin + σθ2 κ2 (Q) L2G Pm,j
+ 2σθ2 κ2 (Q) L2G (Pm,j + Pm−1 )


= 1 − 2σθ λmin + 3σθ2 κ2 (Q) L2G Pm,j
+ 2σθ2 κ2 (Q) L2G Pm−1
If we choose σθ to satisfy
0 < σθ ≤

λmin
,
(Q) L2G

3κ2

(50)

then 3σθ2 κ2 (Q) L2G < σθ λmin , which implies
Pm,j+1 ≤ (1 − σθ λmin ) Pm,j + 2σθ2 κ2 (Q) L2G Pm−1 .
Iterating the above inequality over j = 1, · · · , N − 1 and
using Pm,0 = Pm−1 and Pm,N = Pm , we obtain
Pm = Pm,N


N
−1
X
N
j
2 2
2
≤ 1−σθ λmin +2σθ κ (Q)LG
1−σθ λmin Pm−1
j=0

Stochastic Variance Reduction Methods for Policy Evaluation


N
1−(1−σθ λmin )N
1−σθ λmin +2σθ2 κ2 (Q)L2G
Pm−1
1−(1−σθ λmin )


N 2σ 2 κ2 (Q) L2G
≤ 1 − σθ λmin + θ
Pm−1
σθ λmin


N 2σθ κ2 (Q) L2G
= 1 − σθ λmin +
Pm−1 .
(51)
λmin


=

at m-th iteration, gt is computed using θφm
and wφm
. With
t
t
this definition, φm
t has the following dynamics:
(
φm
if tm 6= t,
t
m+1
φt
=
(53)
m if tm = t.
We can write the m-th iteration’s full gradient as

We can choose
σθ =

n

λmin
,
5κ2 (Q)L2G

N=

1
=
σθ λmin

5κ (Q)L2G
,
λ2min

B=

2

(52)

which satisfies the condition in (50) and results in
Pm ≤ (e−1 + 2/5)Pm−1 ≤ (4/5)Pm−1 .
There are many other similar choices, for example,
σθ =

λmin
,
2
3κ (Q)L2G

N=

3
=
σθ λmin

2

9κ (Q)L2G
,
λ2min

which results in
Pm ≤ (e

−3

+ 2/3)Pm−1 ≤ (3/4)Pm−1 .

These results imply that the number of outer iterations
needed to have E[Pm ] ≤ ] is log(P0 /). For each outer
iteration, the SVRG algorithm need O(nd) operations to
compute the full gradient operator B(θ, w), and then N =
O(κ2 (Q)L2G /λ2min ) inner iterations with each costing O(d)
operations. Therefore the overall computational cost is


 
κ2 (Q) L2G
P0
O
n+
d
log
.
λ2min

Substituting (33) and (37) in the above bound, we get the
overall cost estimate
!
 !
b 2
κ(C)L
P0
G
.
O
n+
d log
2
T
−1
b
b
b

λ (ρI + A C A)
min

Finally, substituting the bounds in (33) and (37) into (52),
we obtain the σθ and N stated in Theorem 2:
σθ =

bT C
b −1 A)
b
λmin (ρI + A
,
b 2
48κ(C)L
G

2

N=

λ2min

b 2
51κ (C)L
G
,
bT C
b −1 A)
b
(ρI + A

which achieves the same complexity.

D

Analysis of SAGA

SAGA in Algorithm 3 maintains a table of previously computed gradients. Notation wise, we use φm
t to denote that


1X
Bt θφm
, wφm
.
t
t
n t=1

For convergence analysis, we define the following quantity:


θφm
− θ?
t
m
.
(54)
∆φt , √1
(wφm
− w? )
t
β
Similar to (53), it satisfies the following iterative relation:
(
∆φ m
if tm 6= t,
t
∆φm+1 =
t
∆m
if tm = t.
With these notations, we can express the vectors used in
SAGA as


n 
n  
1 X ρI −ATt
1X 0
θφm
t
Bm =
,
−
wφm
Ct
n t=1 At
n t=1 bt
t
   

0
ρI −ATtm θm
−
,
htm =
wm
btm
At m
Ctm






0
ρI −ATtm θφm
t
gtm =
−
.
wφm
b
At m
Ctm
tm
t
The dynamics of SAGA can be written as

   

θm+1
θm
σθ
=
−
(Bm + htm − gtm )
wm+1
wm
σw
  

θ
σ
= m − θ
wm
σw
( n 


n  
1X 0
1 X ρI −ATt
θφm
t
+
wφm
Ct
n t=1 At
n t=1 bt
t





)
ρI −ATtm θm
ρI −ATtm θφm
tm
+
−
wm
wφm
Atm Ctm
Atm Ctm
tm
Subtracting (θ? , w? ) from both sides, and using the optimality condition in (39), we obtain

 
 

θm+1 − θ?
θm − θ?
σθ
=
−
wm+1 − w?
wm − w?
σw
( n 


1 X ρI −ATt
θφm
− θ?
t
wφm
− w?
Ct
n t=1 At
t



ρI −ATtm θm − θ?
+
wm − w?
Atm Ctm

Stochastic Variance Reduction Methods for Policy Evaluation

ρI −ATtm
−
Atm Ctm


θφm
− θ?
tm
wφm
−
w?
tm



)
.

√
Multiplying both sides by diag(I, 1/ βI), we get
!
n
σθ X
∆m+1 = ∆m −
G t ∆φ m
t
n t=1


.
− σθ Gtm ∆m − ∆φm
tm

(55)

where Gtm is defined in (18).
For SAGA, we use the following two potential functions:
2

Pm = E Q−1 ∆m 2 ,
 X
n
1
Q−1 Gt ∆φm
Qm = E
t
nt=1

2
2





2
.
= E Q−1 Gtm ∆φm
tm
2

The last equality holds because we use uniform sampling.
We first look at how Pm evolves. To simplify notation, let
!
n


σθ X
vm =
Gt ∆φm
+ σθ Gtm ∆m − ∆φm
,
t
tm
n t=1
so that (55) becomes ∆m+1 = ∆m − vm . We have
h
i
2
Pm+1 = E Q−1 ∆m+1 2
i
h
2
= E Q−1 (∆m − vm )
i
h
2
2
−> −1
Q vm + Q−1 vm 2
= E Q−1∆m 2 −2∆>
mQ
h
i


2
−1
−> −1
Q
v
= Pm − E 2∆>
Q
Q
v
+
E
.
m
m
m
2

where the inequality used λmin , λmin (Λ) = λmin (G) > 0,
which is true under our choice of β = σw /σθ in Section A.1. Next, we bound the last term of Eqn. (56):
h
i
2
E Q−1 vm 2

n
  2

 X
−1 σθ
Gt ∆φm
+σθ Gtm ∆m −∆φm
=E Q
t
tm
n t=1
h
i
2
≤ 2σθ2 E Q−1 Gtm ∆m 2

n
 X
 2
−1 1
2
m
m
+ 2σθ E Q
Gt ∆φt − Gtm ∆φtm
n t=1
i
h
i
h
2
2
k
≤ 2σθ2 E Q−1Gtm ∆m 2 + 2σθ2 E kQ−1 Gtm∆φm
tm
h
i
2
2
−1
2
= 2σθ E Q Gtm ∆m 2 + 2σθ Qm ,
2

2

2

where the first inequality uses ka + bk2 ≤ 2 kak2 + 2 kbk2 ,
and the second inequality holds because for any random
2
2
2
2
variable ξ, E kξ − E [ξ]k2 = E kξk2 − kEξk2 ≤ E kξk2 .
Using similar arguments as in (49), we have
h
i
2
E Q−1Gtm ∆m 2 ≤ κ2(Q)L2G Pm ,
(57)
Therefore, we have

Pm+1 ≤ 1 − 2σθ λmin + 2σθ2 κ2 (Q) L2G Pm
+ 2σθ2 Qm .

(58)

The inequality (58) shows that the dynamics of Pm depends
on both Pm itself and Qm . So we need to find another
iterative relation for Pm and Qm . To this end, we have
" n
#
2
1X
−1
Qm+1 = E
Q Gt ∆φm+1
t
n t=1
2

Since ∆m is independent of tm , we have
1
= E kQ−1 Gtm ∆φm+1 k2
h
i
h
i
tm
n
> −> −1
> −> −1
E 2∆m Q Q vm = E 2∆m Q Q Etm [vm ] ,

X
1
+
kQ−1 Gt ∆φm+1 k2
t
n
where the inner expectation is with respect to tm condit6=tm

tioned on all previous random variables. Notice that
1
(a)
= E kQ−1 Gtm ∆m k2
n
n


1X

Etm Gtm ∆φm
=
Gt ∆φm
,
1 X
tm
t
2
n t=1
+
kQ−1 Gt ∆φm
k
t
n
t6=tm

which implies Etm [vm ] = σθ Etm [Gtm ]∆m = σθ G∆m .
1
1
= E kQ−1 Gtm ∆m k2 − kQ−1 Gtm ∆φm
k2
Therefore, we have
tm
n
n
h
i
h
i

n
2
1 X −1
2
Pm+1 = Pm − E 2σθ ∆Tm Q−T Q−1 G∆m + E Q−1 vm 2
+
kQ Gt ∆φm
k
t
n t=1
h
i
h
i
2
T
−T
−1
−1
= Pm − E2σθ ∆m Q ΛQ ∆m + E Q vm 2
1
1
h
i
h
i
k2 ]
= E[kQ−1 Gtm ∆m k2 ] − E[kQ−1 Gtm ∆φm
tm
2
2
−1
−1
n
n
≤ Pm − 2σθ λmin E Q ∆m
+ E Q vm 2
 X

n
1
h
i
−1
2
2
m
+
E
kQ
G
∆
k
t φt
= (1 − 2σθ λmin )Pm + E Q−1 vm 2 ,
(56)
n t=1

Stochastic Variance Reduction Methods for Policy Evaluation

=

1
1
E[kQ−1 Gtm ∆m k2 ] − E[kQ−1 Gtm ∆φm
k2 ]
tm
n
n


+ E kQ−1 Gtm ∆φm
k2
tm

n−1
1
Qm
= E[kQ−1 Gtm ∆m k2 ] +
n
n
2
2
(b) κ (Q)L
n−1
G
≤
Pm +
Qm .
n
n

λ2min
λ2min
≤
ρ
≤
,
2
2
9 (κ2 (Q) LG +nλmin )
3 (κ2 (Q) L2G +nλ2min )
(59)

where step (a) uses (53) and step (b) uses (57).

which implies
nρ − 1 ≤

To facilitate our convergence analysis on Pm , we construct
a new Lyapunov function which is a linear combination of
Eqn. (58) and Eqn. (59). Specifically, consider
Tm

(or smaller), the second term on the right-hand side of (60)
is non-positive. To see this, we first notice that with this
choice of σθ , we have

nσθ λmin (1 − σθ λmin )
Qm .
= Pm +
κ2 (Q)L2G

Now consider the dynamics of Tm . We have

3 (κ2

1
nλ2min
2
−1≤ −1=− .
(Q) L2G +nλ2min )
3
3

Then, it holds that
2σθ κ2 (Q)L2G + (nρ − 1)λmin (1 − σθ λmin )
2
≤2σθ κ2 (Q)L2G − λmin (1 − σθ λmin )
3
(6n − 2)λ3min
=−
< 0.
9 (κ2 (Q)L2G + nλ2min )

nσθ λmin (1 − σθ λmin )
Qm+1
Therefore (60) implies
κ2 (Q)L2G

≤ 1 − 2σθ λmin + 2σθ2 κ2 (Q) L2G Pm + 2σθ2 Qm
Tm+1 ≤ (1 − ρ)Tm .


n−1
nσθ λmin (1−σθ λmin ) κ2 (Q)L2G
Pm +
Qm Notice that P ≤ T and Q = P . Therefore we have
+
m
m
0
0
κ2 (Q)L2G
n
n

T0 ≤ 2P0 and
2 2
2
2 2
= 1 − σθ λmin + 2σθ κ (Q) LG − σθ λmin Pm

Tm+1 = Pm+1 +

+

2σθ2 κ2 (Q) L2G +(n−1) σθ λmin (1−σθ λmin )
Qm .
κ2 (Q)L2G

Pm ≤ 2(1 − ρ)m P0 .
Using (61), we have

Let’s define
ρ = σθ λmin (G) − 2σθ2 κ2 (Q)L2G ≥

ρ = σθ λmin − 2σθ2 κ2 (Q) L2G .
The coefficient for Pm in the previous inequality can be
upper bounded by 1 − ρ because 1 − ρ − σθ2 λ2min ≤ 1 − ρ.
Then we have
Tm+1
2σθ2 κ2 (Q) L2G +(n−1) σθ λmin (1−σθ λmin )
Qm
κ2 (Q)L2G


nσθ λmin (1 − σθ λmin )
= (1 − ρ) Pm +
Q
m
κ2 (Q)L2G

= (1 − ρ) Tm
+ σθ

2σθ κ

(Q)L2G

To achieve Pm ≤ , we need at most

 

κ2 (Q) L2G
P0
m=O
n+
log
2
λmin


Finally, using the bounds in (33) and (37), we can replace
the step size in (61) by
µρ

σθ =
+ (nρ − 1)λmin (1 − σθ λmin )
Qm .
κ2 (Q)L2G
(60)

Next we show that with the step size
σθ =

3 (κ2

λmin
(Q) L2G + nλ2min )

.

min

2σθ κ2 (Q)L2G + (nρ − 1)λmin (1 − σθ λmin )
Qm
κ2 (Q)L2G
2

nλ2min

iterations. Substituting (37) and (33) in the above bound,
we get the desired iteration complexity
!
 !
b 2
κ(C)L
P0
G
O
n+
.
log
2
T
−1
b
b
b

λ (ρI + A C A)

≤ (1 − ρ) Pm +

+ σθ

9

λ2min
2
κ (Q)L2G +

(61)

3



b 2
8κ2 (C)L
G

+ nµ2ρ

,

bT C
b −1 A)
b as defined in (14).
where µρ = λ2min (ρI + A

