Accelerated and instance-optimal policy evaluation
with linear function approximation

Tianjiao Li⋆ |  Guanghui Lan⋆ |  Ashwin Pananjady⋆,†
---|---|---
Schools of Industrial & Systems Engineering⋆ and Electrical & Computer
Engineering†
---
Georgia Institute of Technology
December 2021; Revised: August 2022
---

###### Abstract

We study the problem of policy evaluation with linear function approximation
and present efficient and practical algorithms that come with strong
optimality guarantees. We begin by proving lower bounds that establish
baselines on both the deterministic error and stochastic error in this
problem. In particular, we prove an oracle complexity lower bound on the
deterministic error in an instance-dependent norm associated with the
stationary distribution of the transition kernel, and use the local asymptotic
minimax machinery to prove an instance-dependent lower bound on the stochastic
error in the i.i.d. observation model. Existing algorithms fail to match at
least one of these lower bounds: To illustrate, we analyze a variance-reduced
variant of temporal difference learning, showing in particular that it fails
to achieve the oracle complexity lower bound. To remedy this issue, we develop
an accelerated, variance-reduced fast temporal difference algorithm (VRFTD)
that simultaneously matches both lower bounds and attains a strong notion of
instance-optimality. Finally, we extend the VRFTD algorithm to the setting
with Markovian observations, and provide instance-dependent convergence
results. Our theoretical guarantees of optimality are corroborated by
numerical experiments.

##  1 Introduction

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

A fundamental building block in RL is the problem of _policy evaluation_ , in
which we are interested in estimating the long-term (discounted) _value_ of
each state under a fixed policy with sample access to the transition and
reward functions. The literature considers three observation models for
transition and reward samples, namely the generative model, the so-called
“i.i.d.” model, and the Markovian noise model111Two of these observation
models are formally discussed in Section 2.. Furthermore, in modern
applications with large state spaces, it is common to seek an approximation to
the true value function within the span of a small number of basis functions,
a setting that is commonly known as _linear function approximation_. In the
canonical setting of the problem, one is interested in using random
observations to compute an approximate value function within the subspace,
with the distance between the true value function and its approximation being
measured according to an instance-dependent “weighted $\ell_{2}$-norm” that
depends on the stationary distribution of the transition kernel.

It is common to use stochastic approximation (SA) algorithms to solve the
policy evaluation problem in the setting described above. Given the iterative
nature of these algorithms, their convergence rates can be decomposed into two
types of error: a _deterministic error_ that measures how fast the algorithm
converges to its fixed point even in the absence of noise, and a _stochastic
error_ that measures the contribution of the noise. Either of these errors
could dominate in practice. While the stochastic error is typically larger in
noisy problems, the deterministic error of the algorithm can dominate, for
example, in settings with multiple processors. In particular, the collection
and use of multiple samples/trajectories in parallel can reduce the stochastic
error considerably.

Loosely speaking, the deterministic error is measured in terms of the _oracle
complexity_ of the algorithm and the stochastic error is measured in terms of
the _sample complexity_. The eventual goal of algorithm design is to develop
practical algorithms that have optimal oracle and sample complexities.
Ideally, these optimality guarantees should be _instance-specific_ , in that
they depend explicitly on the problem at hand and enable us to draw
distinctions between the performance profiles of different algorithms.

With this context in hand, let us briefly discuss the state of the art
relevant to characterizing the complexity of policy evaluation and related
problems. Classical work by Nemirovsky (1991, 1992) established oracle
complexity lower bounds for solving linear operator equations in
$\ell_{2}$-norm. However, these results do not extend to the more specific
policy evaluation setting under the weighted $\ell_{2}$-norm. On the other
hand, Khamaru et al. (2021a) recently provided an instance-specific analysis
of the sample complexity of policy evaluation under the $\ell_{\infty}$-norm,
focusing on the generative observation model without function approximation.
Mou et al. (2020) studied the broader problem of solving projected fixed point
equations, with a focus on characterizing the error incurred due to projection
onto a subspace. By virtue of studying the more general problem, the lower
bounds on the statistical error proved in this paper are not specific enough
to capture the policy evaluation setting with function approximation.
Concurrent work by Mou et al. (2021) studied SA methods for solving linear
fixed point equations with Markovian samples and established a non-asymptotic,
instance-dependent lower bound. Given this state of affairs, the central
question that motivates this paper is the following:

What are the optimal oracle and sample complexities of policy evaluation in
weighted $\ell_{2}$-norm with linear function approximation, and do existing
algorithms achieve these bounds?

###  1.1 Related work

There is a large literature on stochastic approximation for policy evaluation.
The most popular stochastic iterative algorithm used for policy evaluation is
temporal difference (TD) learning; see Dann et al. (2014) for a survey. The TD
learning algorithm was first introduced by Sutton (1988), and convergence
guarantees for TD have been proven in both asymptotic and non-asymptotic
settings. Asymptotic convergence of TD with linear function approximation was
established in Tsitsiklis and Van Roy (1997), and other classical asymptotic
guarantees include those due to Borkar and co-authors (Borkar and Meyn, 2000;
Borkar, 2009). The vanilla TD algorithm can also be combined with the iterate
averaging technique, and the asymptotic convergence of this algorithm was
shown by Tadic (2004), who extended the convergence results for solving noisy
linear systems (Polyak and Juditsky, 1992).

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

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

During the past decade, there has been a flurry of parallel work in stochastic
optimization on developing first-order methods with variance-reduction; early
examples include IAG (Blatt et al., 2007), SAG (Schmidt et al., 2017), SVRG
(Johnson and Zhang, 2013; Xiao and Zhang, 2014), and SAGA (Defazio et al.,
2014). There are several papers that apply variance reduction to reinforcement
learning, e.g., Du et al. (2017); Papini et al. (2018); Wai et al. (2018).
Recent work in policy evaluation has shown that variance reduction techniques
can also be applied to algorithms of the TD-type (Korda and La, 2015; Wai et
al., 2019; Xu et al., 2020; Khamaru et al., 2021a). Among these, the paper
Khamaru et al. (2021a) is motivated by the desire to draw distinctions between
RL algorithms with similar worst-case performance and follows a line of work
deriving instance-dependent bounds on the stochastic error in policy
evaluation (Pananjady and Wainwright, 2021; Li et al., 2020). Specifically,
the results of Khamaru et al. (2021a) capture the optimal instance-dependent
stochastic error in the $\ell_{\infty}$-norm. However, the optimal sample
complexity is not achieved in Khamaru et al. (2021a) since the algorithm
requires $\mathcal{O}\\{1/(1-\gamma)^{3}\\}$ samples in each epoch.

###  1.2 Contributions and organization

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

  * •

Lower bounds. We construct a worst-case instance that shows an oracle
complexity lower bound of order
$\Omega\\{(1-\gamma)^{-1}\cdot\log(1/\epsilon)\\}$ for any iterative method
whose iterates lie within the linear span of the initial point $v_{0}$ and
subsequent temporal differences, to converge to $\epsilon$-error in weighted
$\ell_{2}$-norm. We also prove a lower bound on sample complexity using the
classical local minimax theorem (Hájek, 1972; Le Cam, 1972; Le Cam and Yang,
2000) to provide an instance-specific baseline for algorithm design.

  * •

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

  * •

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

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

###  1.3 Notation

For a positive integer n, we define $[n]:=\\{1,2,...,n\\}$. We let
$\mathbf{1}$ denote the all-ones vector in $\mathbb{R}^{D}$. We let $e_{j}$
denote the $j$-th standard basis vector in $\mathbb{R}^{D}$. Let
$\mathbb{I}_{S}:X\rightarrow\\{0,1\\}$ denote the indicator function of the
subset $S\subseteq X$. Given a vector $x\in\mathbb{R}^{m}$, denote its $i$-th
entry by222In situations in which there is no ambiguity, we also use $x_{i}$
to denote the $i$-th coordinate of a vector $x$. $x_{(i)}$. Let
$\|x\|_{1}:=\textstyle{\sum}_{i=1}^{m}|x_{(i)}|$,
$\|x\|_{2}:=\sqrt{\textstyle{\sum}_{t=1}^{m}x_{(i)}^{2}}$ and
$\|x\|_{\infty}:=\max_{i\in[m]}|x_{(i)}|$ denote the $\ell_{1}$, $\ell_{2}$
and $\ell_{\infty}$-norms respectively. Given a matrix $A$, denote its
$(i,j)$-th entry by $P_{i,j}$. Let $\|A\|_{2}$ denote the spectral norm of
matrix $A$. We let $\lambda_{\min}(A)$ and $\lambda_{\max}(A)$ denote the
smallest and largest eigenvalue of a square matrix $A$, respectively. For a
symmetric positive definite matrix $A$, define the inner product $\langle
x,y\rangle_{A}:=x^{\top}Ay$ and the associated norm
$\|x\|_{A}:=\sqrt{x^{\top}Ax}$. We refer to $\|x\|_{A}$ as the $\ell_{A}$-norm
of $x$.

##  2 Background and problem setting

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

###  2.1 Markov reward process and policy evaluation

An MRP is described by a tuple $(\mathcal{S},\mathsf{P},R,\gamma)$, where
$\mathcal{S}=[D]$ denotes the state space, $\mathsf{P}$ is the transition
kernel, $R$ is the reward function and $\gamma\in(0,1)$ is the discount
factor. At each iteration, the system moves from the current state
$s\in\mathcal{S}$ to some state $s^{\prime}\in\mathcal{S}$ with probability
$\mathsf{P}(s^{\prime}|s)$, while the agent realizes the reward
$R(s,s^{\prime})$. We denote by
$r(s):=\textstyle{\sum}_{s^{\prime}\in\mathcal{S}}\mathsf{P}(s^{\prime}|s)R(s,s^{\prime})$
the expected instantaneous reward generated at state $s$. Let $P$ denote the
transition probability matrix having $(i,j)$-th entry
$P_{i,j}=\mathsf{P}(j|i)$. The reward $R$ can also be written in matrix form,
i.e., $R_{i,j}=R(i,j)$. The value function specifies the infinite-horizon
discounted reward as a function of the initial state:

$\displaystyle
v^{*}(s):={\mathbb{E}}\big{[}\textstyle{\sum}_{t=0}^{\infty}\gamma^{t}R(s_{t},s_{t+1})|s_{0}=s\big{]}.$

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

$\displaystyle v^{*}=\gamma Pv^{*}+r.$ (1)

Throughout this paper, we assume that the Markov chain is aperiodic, ergodic
and that there exists a unique stationary distribution
$\pi:=(\pi_{1},...,\pi_{D})$ with strictly positive entries, satisfying $\pi
P=\pi$. Let $\Pi:=\text{diag}(\pi_{1},...,\pi_{D})$ denote a $D\times D$
diagonal matrix whose non-zero elements are given by the entries of the
stationary distribution.

###  2.2 Linear function approximation

In modern applications with large state spaces, it is common to seek
approximate solutions to the Bellman equation (1), and the standard approach
is to choose a $d$-dimensional subspace $\mathbb{S}$ for the purposes of
approximation. In particular, one chooses
${\mathbb{S}}:=\text{span}\\{\psi_{1},...,\psi_{d}\\}$ for $d$ linearly
independent basis vectors $\psi_{1},\ldots,\psi_{d}$. For each state $s\in[D]$
we let $\psi(s):=[\psi_{1}(s),\psi_{2}(s),...,\psi_{d}(s)]^{\top}$ denote its
feature vector. Letting $\Pi_{\mathbb{S}}$ denote the projection onto the
subspace with respect to the $\|\cdot\|_{\Pi}$-norm, define $\bar{v}$ as the
solution to the _projected fixed point equation_

$\displaystyle\bar{v}=\Pi_{\mathbb{S}}(\gamma P\bar{v}+r).$ (2)

It is convenient to write this projection in matrix notation. Let
$\Psi:=[\psi_{1},\psi_{2},...,\psi_{d}]^{\top}$, and for $v^{\Diamond}$ in
$\mathbb{S}$, use $\theta^{\Diamond}$ to denote its corresponding
parameterization in $\mathbb{R}^{d}$, e.g.,
$\Psi^{\top}\theta^{\prime}=v^{\prime}$. With this shorthand, equation (2) can
be equivalently written as

$\displaystyle\Psi\Pi\Psi^{\top}\bar{\theta}=\Psi\Pi\gamma
P\Psi^{\top}\bar{\theta}+\Psi\Pi r.$ (3)

It is convenient in the analysis to have access to an orthonormal basis
spanning the projected space $\mathbb{S}$. Define the matrix
$B\in\mathbb{R}^{d\times d}$ by letting
$B_{i,j}:=\langle\psi_{i},\psi_{j}\rangle_{\Pi}$ for each $i,j$, and let

$\Phi:=[\phi_{1},\phi_{2},...,\phi_{d}]^{\top}=B^{-\frac{1}{2}}\Psi.$

By construction, the vectors $\phi_{1},\ldots,\phi_{d}$ satisfy
$\langle\phi_{i},\phi_{j}\rangle_{\Pi}=\mathbb{I}(i=j)$. Next, define the
scalars

$\displaystyle\beta:=\lambda_{\max}(B),~{}~{}\text{and}~{}\mu:=\lambda_{\min}(B),$
(4)

so that $\beta/\mu$ is the condition number of the covariance matrix of the
features. Finally, let

$\displaystyle M:=\gamma\Phi\Pi P\Phi^{\top}$

denote the $d$-dimensional matrix that describes the action of $\gamma P$ on
the projected space $\mathbb{S}$.

###  2.3 Observation models and problem statement

We start by introducing the i.i.d. observation model, in which we have access
to a black box or simulator that generates samples from the transition kernel
and reward functions. In particular, we observe independent tuples
$\xi_{i}=(s_{i},s_{i}^{\prime},R(s_{i},s_{i}^{\prime}))$, such that

$s_{i}\sim\omega,~{}~{}s_{i}^{\prime}\sim P(\cdot|s_{i}),$

where $\omega:=(\omega_{1},...,\omega_{D})$ is a distribution with strictly
positive entries, and we use the shorthand
$\Omega:=\operatorname{diag}([\omega_{1},...,\omega_{D}])$. A natural and
popular choice is $\omega=\pi$, in which case the i.i.d. model is meant to
approximate the stationary Markov chain.

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

Our goal in both cases is to use the observations to generate an estimator
$\widehat{v}_{n}$ of $v^{*}$ which satisfies an _oracle inequality_ of the
form

$\displaystyle{\mathbb{E}}\|\widehat{v}_{n}-v^{*}\|_{\Pi}^{2}\leq\mathcal{O}(1)\|\bar{v}-v^{*}\|_{\Pi}^{2}+\delta_{n}\|v_{0}-\bar{v}\|_{\Pi}^{2}+\epsilon_{n,\sigma},$
(5)

where $v_{0}$ is the initial iterate of the algorithm. The three terms
appearing on the RHS of inequality (5) all have concrete interpretations. The
first term $\|\bar{v}-v^{*}\|_{\Pi}^{2}$ characterizes the _approximation
error_ incurred by the linear function approximation. As a point of the
background, we recall the following instance-dependent upper bound on the
approximation error due to Mou et al. (2020):

$\displaystyle\|\bar{v}-v^{*}\|_{\Pi}^{2}\leq\mathcal{A}(M,\gamma)\inf_{v\in
S}\|v-v^{*}\|_{\Pi}^{2},$ (6)

where
$\mathcal{A}(M,\gamma)=1+\lambda_{\max}\left((I-M)^{-1}(\gamma^{2}I_{d}-MM^{\top})(I-M)^{-\top}\right)$.
See Mou et al. (2020) for a proof, alongside guarantees of information-
theoretic optimality.

This work focuses on sharply analyzing the last two terms on the RHS of
inequality (5), both of which have concrete operational interpretations. The
term $\delta_{n}\|v_{0}-\bar{v}\|_{\Pi}^{2}$ is the _deterministic error_ ,
which characterizes the convergence of the iterative algorithm in the purely
deterministic setting. Specifically, the term $\delta_{n}$, which should tend
to zero as the number of iterations (or oracle calls) $n$ goes to infinity,
quantifies how fast the discrepancy between the initialization $v_{0}$ and the
approximate solution $\bar{v}$ diminishes by running the iterative algorithm.
The third term $\epsilon_{n,\sigma}$ is the _stochastic error_ , which is
incurred due to the stochastic observation model. Here we use the notation
$\sigma$ as a placeholder for the “noise level” in the observed samples. One
should expect the stochastic error $\epsilon_{n,\sigma}$ to go to zero as $n$
goes to infinity or as $\sigma$ goes to zero. Several previous works mix the
deterministic error with stochastic error in their guarantees (see, e.g.,
Bhandari et al., 2018; Mou et al., 2020). However, the key benefit of
separating the deterministic error from the stochastic error is that it allows
a clean understanding of situations in which either the observations have low
noise or parallel implementation may be available. In these cases, the
deterministic error dominates the overall convergence rate of the algorithm,
and so having algorithms that attain the optimal deterministic error is a key
desideratum.

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

##  3 Lower bounds in weighted $\ell_{2}$-norm

We study the oracle complexity lower bound on deterministic error in Section
3.1 and the instance-specific stochastic error lower bound in Section 3.2.

###  3.1 Oracle complexity lower bound on deterministic error

It is well-known that a linear rate can be achieved for the deterministic
policy evaluation problem, and the convergence rate is highly dependent on the
effective horizon $(1-\gamma)^{-1}$ (Puterman, 2014). Accordingly, our goal in
this section is to prove an oracle complexity lower bound in terms of
$(1-\gamma)^{-1}$, which can be done even in the tabular setting in which the
subspace is all of $\mathbb{R}^{D}$. The following assumption on the oracle
captures algorithms in the temporal difference learning family.

######  Assumption 1 (Amenable iterative method).

An _amenable_ iterative method $\mathcal{M}$ generates a sequence of iterates
$v_{k}$ such that

$\displaystyle v_{k}\in
v_{0}+\mathrm{span}\\{G(v_{0}),G(v_{1}),...,G(v_{k-1})\\},\quad k\geq 1,$ (7)

where $G(v)=(I-\gamma P)v-r$.

Noting that $G(v)$ is precisely the temporal difference operator applied at
the point $v$, an amenable algorithm is one whose iterates are always in the
linear span of the initial point $v_{0}$ and subsequent temporal differences.
The linear span assumption is commonly used in proving oracle complexity lower
bounds (Nesterov, 2003; Ouyang and Xu, 2021), and as such, nearly all the
algorithms in the temporal difference family can be shown to be amenable. The
sole exceptions that we are aware of occur in cases where there are
projections involved in the algorithm, e.g., Bhandari et al. (2018). However,
in policy evaluation problems with unbounded feasible region $\mathbb{R}^{D}$,
projection steps are often unnatural and vanilla TD algorithms are able to
attain similar performance (see, e.g. Kotsalis et al., 2020b). The following
theorem provides an oracle complexity lower bound for policy evaluation
problem under the $\ell_{\Pi}$-norm for amenable algorithms.

######  Theorem 1.

Fix a constant $\gamma>\tfrac{1}{2}$ and an initialization $v_{0}$. There
exists a transition kernel $P$ and an expected reward vector $r$ such that any
iterative method $\mathcal{M}$ satisfying Assumption 1 produces iterates
$\\{v_{k}\\}_{k\geq 1}$ satisfying the following. If $(D,k)$ satisfies
$\tfrac{1-(2\gamma-1)^{2D-2k}}{1-(2\gamma-1)^{2D}}\geq\tfrac{1}{2}$, then

$\displaystyle\|v_{k}-v^{*}\|_{\Pi}^{2}\geq\tfrac{1}{2}(2\gamma-1)^{2k}\|v_{0}-v^{*}\|_{\Pi}^{2},$
(8)

where $v^{*}$ is the solution of equation (1).

See Section 7.2 for the proof of this theorem.

Noting that $2\gamma-1=1-2(1-\gamma)$, Theorem 1 shows an oracle complexity
lower bound
$\mathcal{O}\\{\tfrac{1}{1-\gamma}\log(\tfrac{\|v_{0}-v^{*}\|_{\Pi}^{2}}{\epsilon})\\}$
for finding a solution $\widehat{v}\in\mathbb{R}^{m}$ such that
$\|\widehat{v}-v^{*}\|_{\Pi}^{2}\leq\epsilon$. It should be noted that the
metric (i.e., the $\ell_{\Pi}$-norm) used in Theorem 1 depends on the problem
instance through the stationary distribution of the transition kernel $P$.
Such an instance-dependent metric makes the construction of our worst-case
instance non-standard and challenging.

On a related note, it is instructive to recall that classical oracle
complexity bounds for solving linear operator equations (Nemirovsky, 1991,
1992) allow the conjugate operator to be queried within the oracle, making the
method class wider than the class of amenable algorithms captured in
Assumption 1. On the one hand, the conjugate operator is not natural for
solving a policy evaluation problem under stochastic settings since the vector
$(I-\gamma P)^{\top}v$ is hard to estimate with transition and reward samples.
On the other hand, our construction used in proving Theorem 1 naturally
extends to this wider method class, and we provide an even stronger
deterministic error lower bound than Theorem 1 in Appendix A.

###  3.2 Instance-specific lower bound on stochastic error

We now turn our attention to proving lower bounds on the instance-specific
sample complexity under the i.i.d. observation model introduced in Section
2.3. We assume that the feature matrix $\Psi$ is fixed and known, and let
$\vartheta=(\omega,P,R)$ denote an individual problem instance parameterized
by the initial state distribution $\omega$, transition kernel $P$, and reward
function $R$. Note at this juncture that we do not require that $\omega=\pi$;
this is akin to the so-called _off-policy_ situation in which the sampling (or
behavior) policy may differ from the policy that we are interested in
evaluating. Our result will apply in this general case; but given that the
initial state is drawn from the distribution $\omega$, it is convenient to
consider solving the projected fixed point equation with respect to the
$\|\cdot\|_{\Omega}$-norm (cf. Eq. (3)), written as

$\displaystyle\Psi\Omega\Psi^{\top}\theta=\Psi\Omega\gamma
P\Psi^{\top}\theta+\Psi\Omega r.$

Use the function
$\bar{\theta}(\vartheta):=(\Psi\Omega\Psi^{\top}-\Psi\Omega\gamma
P\Psi^{\top})^{-1}\Psi\Omega r$ to denote the target of interest.

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

$\displaystyle\mathfrak{N}(\vartheta;\epsilon):=\\{\vartheta^{\prime}=(\omega^{\prime},P^{\prime},R^{\prime}):\|\omega-\omega^{\prime}\|_{2}+\|P-P^{\prime}\|_{F}+\|R-R^{\prime}\|_{F}\leq\epsilon\\}.$

Define the matrix $\widetilde{B}\in\mathbb{R}^{d\times d}$ by
$\widetilde{B}_{i,j}:=\langle\psi_{i},\psi_{j}\rangle_{\Omega}$ for
$i,j\in[d]$. Thus $\widetilde{B}$ satisfies

$\displaystyle\widetilde{B}^{-\frac{1}{2}}\Psi\Omega\Psi\widetilde{B}^{-\frac{1}{2}}=I_{d}.$
(9)

Adopting the $\ell_{\widetilde{B}}$-norm as our loss function, define the
following local asymptotic minimax risk (Hájek, 1972; Le Cam, 1972):

$\displaystyle\mathfrak{M}(\vartheta):=\lim_{c\rightarrow\infty}\lim_{N\rightarrow\infty}\inf_{\widehat{\theta}_{N}}\sup_{\vartheta^{\prime}\in\mathfrak{N}(\vartheta;c/\sqrt{N})}N\cdot{\mathbb{E}}_{\vartheta^{\prime}}\left[\left\|\widehat{\theta}_{N}-\bar{\theta}(\vartheta^{\prime})\big{)}\right\|_{\widetilde{B}}^{2}\right].$
(10)

The infimum in Eq. (10) is taken over all estimators $\widehat{\theta}_{N}$
that are measurable functions of $N$ observations drawn according to the
i.i.d. observation model. In contrast to the global minimax risk—which takes a
supremum of the risk over all the problem instances within a reasonable
class—the local minimax risk $\mathfrak{M}(\vartheta)$ looks for the hardest
alternative in a small neighborhood of the instance $\vartheta$ with diameter
$c/\sqrt{N}$. To capture the hardest local alternative (in an asymptotic
sense) it suffices to take the diameter of the neighborhood to be of the order
$1/\sqrt{N}$. Invoking Eq. (9) yields the equivalent definition

$\displaystyle\mathfrak{M}(\vartheta)=\lim_{c\rightarrow\infty}\lim_{N\rightarrow\infty}\inf_{\widehat{\theta}_{N}}\sup_{\vartheta^{\prime}\in\mathfrak{N}(\vartheta;c/\sqrt{N})}N\cdot{\mathbb{E}}_{\vartheta^{\prime}}\left[\left\|\Psi^{\top}\widehat{\theta}_{N}-\Psi^{\top}\bar{\theta}(\vartheta^{\prime})\big{)}\right\|_{\Omega}^{2}\right].$
(11)

The following proposition characterizes the local asymptotic risk
$\mathfrak{M}(\vartheta)$ explicitly.

######  Proposition 1.

Consider the i.i.d. observation model with the initial state drawn from the
distribution $\omega$. Let $Z\in\mathbb{R}^{d}$ be a multivariate Gaussian

$\displaystyle
Z\sim\mathcal{N}\big{(}0,(I_{d}-\widetilde{M})^{-1}\widetilde{\Sigma}(I_{d}-\widetilde{M})^{-T}\big{)},$

where
$\widetilde{\Sigma}:=\operatorname{cov}\big{[}\widetilde{B}^{-\frac{1}{2}}\big{(}\langle\psi(s)-\gamma\psi(s^{\prime}),\bar{\theta}\rangle-R(s,s^{\prime})\big{)}\psi(s)\big{]}$
and $\widetilde{M}:=\gamma\widetilde{B}^{-\frac{1}{2}}\Psi\Omega
P\Psi\widetilde{B}^{-\frac{1}{2}}$. Then we have

$\displaystyle\mathfrak{M}(\vartheta)={\mathbb{E}}[\|Z\|_{2}^{2}]=\mathrm{trace}\left\\{(I_{d}-\widetilde{M})^{-1}\widetilde{\Sigma}(I_{d}-\widetilde{M})^{-T}\big{)}\right\\}.$
(12)

See Section 7.3 for the proof of this theorem.

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

Let us now specialize Proposition 1 by taking $\omega=\pi$, where $\pi$ is the
stationary distribution of the transition kernel $P$. Denote by
$\vartheta_{\pi}:=(\pi,P,R)$ the instance of interest. Let

$\displaystyle\bar{\Sigma}_{\mathsf{iid}}:=\operatorname{cov}\big{[}B^{-\frac{1}{2}}\big{(}\langle\psi(s)-\gamma\psi(s^{\prime}),\bar{\theta}\rangle-R(s,s^{\prime})\big{)}\psi(s)\big{]}\quad\text{
for }s\sim\pi\text{ and }s^{\prime}\sim\mathsf{P}(\cdot|s).$ (13)

Applying Proposition 1, the local asymptotic minimax risk (11) under this
setting is then given by

$\displaystyle\lim_{c\rightarrow\infty}\lim_{N\rightarrow\infty}\inf_{\widehat{\theta}_{N}}\sup_{\vartheta^{\prime}\in\mathfrak{N}(\vartheta_{\pi};c/\sqrt{N})}N\cdot{\mathbb{E}}_{\vartheta^{\prime}}$
$\displaystyle\left[\left\|\Psi^{\top}\widehat{\theta}_{N}-\Psi^{\top}\bar{\theta}(\vartheta^{\prime})\big{)}\right\|_{\Pi}^{2}\right]$
$\displaystyle=\mathrm{trace}\left\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-T}\big{)}\right\\}.$
(14)

Taking stock, we have proved two lower bounds (8) and (3.2) on the
deterministic and stochastic errors in $\ell_{\Pi}$-norm under the i.i.d.
observation model $s\sim\pi$ and $s^{\prime}\sim\mathsf{P}(\cdot|s)$. Given
these baselines, it is natural to ask whether there is a practical iterative
algorithm in the TD family that can achieve both lower bounds, which is the
main focus of Section 4.

##  4 Algorithms for policy evaluation in the i.i.d. setting

Taking both lower bounds proved in Section 3 as our baseline, we now turn our
attention to the question of algorithm design. In this section, we assume the
i.i.d. observation model introduced in Section 2.3 with $\omega=\pi$. In order
to state the results clearly, we require some additional notation. For
$\theta\in\mathbb{R}^{d}$, we define the deterministic operator for solving
equation (3) as

$\displaystyle g(\theta)=\Psi\Pi(\Psi^{\top}\theta-r-\gamma
P\Psi^{\top}\theta);$ (15)

note that $\bar{\theta}$ is the solution to $g(\theta)=0$ The corresponding
stochastic operator calculated from sample $\xi_{i}$ is defined as

$\displaystyle\widetilde{g}(\theta,\xi_{i})=\left(\langle\psi(s_{i}),\theta\rangle-R(s_{i},s_{i}^{\prime})-\gamma\langle\psi(s_{i}^{\prime}),\theta\rangle\right)\psi(s_{i});$
(16)

note that
${\mathbb{E}}_{s_{i}\sim\pi,s_{i}^{\prime}\sim\mathsf{P}(\cdot|s_{i})}[\widetilde{g}(\theta,\xi_{i})]=g(\theta)$.
To characterize the “variance” of the stochastic operator under the i.i.d.
observation model, we make the following assumption:

######  Assumption 2.

There exists a constant $\varsigma\geq 0$ such that for every
$\theta,\theta^{\prime}\in\mathbb{R}^{d}$,

$\displaystyle{\mathbb{E}}\|\widetilde{g}(\theta,\xi)-\widetilde{g}(\theta^{\prime},\xi)-\big{(}g(\theta)-g(\theta^{\prime})\big{)}\|^{2}_{2}\leq\varsigma^{2}\|v-v^{\prime}\|_{\Pi}^{2},$
(17)

where $v=\Psi^{\top}\theta$ and $v^{\prime}=\Psi^{\top}\theta^{\prime}$.

In words, instead of bounding the “variance” of the stochastic operator
directly as in Kotsalis et al. (2020b), Assumption 2 guarantees that the
variance of the difference between stochastic operators with different
variables $\theta$ under the same data $\xi$ is upper bounded by the distance
between the variables. This assumption is critical for implementing the
variance-reduction techniques and capturing the instance-dependent stochastic
error at the approximate solution $\bar{\theta}$. Clearly, the parameter
$\varsigma^{2}$ is bounded provided the features $\psi(s)$ are bounded, and
provides a natural measure of “noise” in the problem. Accordingly, we make use
of Assumption 2 throughout Sections 4 and 5.

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

###  4.1 A warm-up algorithm: variance-reduced temporal difference learning

Variance-reduced temporal difference learning (VRTD) solves the policy
evaluation problem using epochs. With a slight ambiguity of notation, we let
$v_{t}$ and its corresponding parameterization $\theta_{t}$ denote the
iterates generated within each epoch, and let $v^{0}$ and its corresponding
parameterization $\theta^{0}$ denote the initialization of the algorithm. At
the beginning of each epoch $k$, the algorithm uses $N_{k}$ samples to compute
an averaged stochastic operator $\widehat{g}$ and evaluates it at a point
$\widetilde{\theta}$, where $\widetilde{\theta}$ should be understood as the
best current approximation of the optimal solution. The vector
$\widehat{g}(\widetilde{\theta})$ is used to recenter the updates in each
epoch.

Algorithm 1 Variance-reduced Temporal Difference Algorithm under i.i.d
observations

Input: $\theta^{0}=\widehat{\theta}_{0}\in\mathbb{R}^{d}$, $\eta>0$,
$\\{\zeta_{t}\\}_{t=1}^{T}\geq 0$ and
$\\{N_{k}\\}_{k=1}^{K}\subset\mathbb{Z}_{+}$.

for $k=1,\ldots,K$ do

Set $\theta_{1}=\widetilde{\theta}=\widehat{\theta}_{k-1}$. Collect $N_{k}$
samples $\xi_{i}^{k}=(s_{i},s_{i}^{\prime},R(s_{i},s_{i}^{\prime}))$ from the
i.i.d. model.

Calculate
$\widehat{g}(\widetilde{\theta})=\tfrac{1}{N_{k}}\textstyle{\sum}_{i=1}^{N_{k}}\widetilde{g}(\widetilde{\theta},\xi_{i}^{k}).$

for $t=1,\ldots,T$ do

Collect a sample $\xi_{t}=(s_{t},s_{t}^{\prime},R(s_{t},s_{t}^{\prime}))$ from
the i.i.d. observation model and compute
$\theta_{t+1}=\theta_{t}-\eta\left(\widetilde{g}(\theta_{t},\xi_{t})-\widetilde{g}(\widetilde{\theta},\xi_{t})+\widehat{g}(\widetilde{\theta})\right).$
(18)

end for

Output of the epoch:
$\widehat{\theta}_{k}=\frac{\sum_{t=1}^{T+1}\zeta_{t}\theta_{t}}{\sum_{t=1}^{T+1}\zeta_{t}}.$
(19)

end for

Note that this algorithm is distinct from previous instantiations of variance-
reduced temporal difference algorithms (Xu et al., 2020; Khamaru et al.,
2021a), since the output of each epoch (19) is a weighted average of the
iterates. The following theorem provides a convergence guarantee on the VRTD
algorithm.

######  Theorem 2.

Consider the i.i.d. observation model with the initial state drawn from the
distribution $\pi$. Fix the total number of epochs $K$ and a positive integer
$N$. Assume that for each epoch $k\in[K]$, the parameters $\eta$, $N_{k}$ and
$T$ satisfy

$\displaystyle\eta\leq\min\left\\{\tfrac{1-\gamma}{2\beta(1+\gamma)^{2}},\tfrac{1-\gamma}{32\varsigma^{2}}\right\\},~{}~{}T\geq\tfrac{32}{\mu(1-\gamma)\eta},~{}\text{and}~{}~{}N_{k}\geq\left\\{\tfrac{38\varsigma^{2}}{\mu(1-\gamma)^{2}},(\tfrac{3}{4})^{K-k}N\right\\}.$

Set the output of the epoch to be
$\widehat{v}_{k}:=\frac{\sum_{t=1}^{T}\eta(1-\gamma)v_{t}+(1/\beta)v_{T+1}}{T\eta(1-\gamma)+(1/\beta)}$.
Then for each $\delta>0$, we have

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{K}-v^{*}\|_{\Pi}^{2}]$
$\displaystyle\leq(1+\delta)\mathcal{A}(M,\gamma)\inf_{v\in
S}\|v-v^{*}\|_{\Pi}^{2}$
$\displaystyle\quad+(1+\tfrac{1}{\delta})\left[\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|_{\Pi}^{2}+\tfrac{15}{N}\mathrm{trace}\left((I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\right)\right].$
(20)

See Section 7.4 for the detailed proof of this theorem.

The first term in the bound (2) is the approximation error term alluded to
previously; let us extract the deterministic and stochastic errors from the
remaining terms. The number of epochs required by the VRTD method to find a
solution $\widehat{v}\in\mathbb{R}^{D}$, such that
${\mathbb{E}}[\|\widehat{v}-\bar{v}\|_{\Pi}^{2}]\leq\epsilon$ is bounded by
$\mathcal{O}\big{\\{}\log(\|v^{0}-\bar{v}\|_{\Pi}^{2}/\epsilon)\big{\\}}$. The
total number of samples used is $\textstyle{\sum}_{k=1}^{K}(T+N_{k})$, which
is of the order

$\displaystyle\underbrace{\tfrac{\beta}{(1-\gamma)^{2}\mu}\log(\tfrac{\|v^{0}-\bar{v}\|_{\Pi}^{2}}{\epsilon})}_{\text{deterministic
error}}+\underbrace{\tfrac{\varsigma^{2}}{(1-\gamma)^{2}\mu}\log(\tfrac{\|v^{0}-\bar{v}\|_{\Pi}^{2}}{\epsilon})+\tfrac{\mathrm{trace}\big{(}(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\big{)}}{\epsilon}}_{\text{stochastic
error}}.$ (21)

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

##### Comparing the upper and lower bounds

We first focus on the stochastic error in the bound (21). The VRTD algorithm
requires at least $\mathcal{O}\\{\tfrac{\varsigma^{2}}{(1-\gamma)^{2}\mu}\\}$
samples in each epoch, which accounts for the first term. Note that with noisy
observations, it is necessary to have
$\mathcal{O}\\{\tfrac{1}{(1-\gamma)^{2}}\\}$ samples in order to obtain an
estimate of the value function within $\mathcal{O}(1)$ error, so this higher-
order term is natural. The dominating stochastic error term is the last term,
and matches the lower bound in equation (3.2). Therefore the VRTD algorithm is
instance-optimal in terms of its stochastic error.

Next, we turn our attention to the deterministic error, noticing that the
dependence on $\tfrac{1}{1-\gamma}$ is quadratic. Comparing with the oracle
complexity lower bound proved in Theorem 1, this quadratic dependence is
suboptimal. This shortcoming motivates us to develop an accelerated algorithm
in the next subsection.

##### Comparing with related work

To our knowledge, the only work using variance reduction that captures the
correct instance-specific stochastic error is that of Khamaru et al. (2021a),
which showed that the VRTD algorithm can match the lower bound on stochastic
error in $\ell_{\infty}$-norm. However, their guarantees require
$\mathcal{O}\\{\tfrac{1}{(1-\gamma)^{3}}\\}$ samples in each epoch to compute
the recentered update, and this sample size is suboptimal. In addition, the
deterministic error proved in this paper is of the order
$\tfrac{1}{(1-\gamma)^{2}\epsilon}$ which is also suboptimal.333In more
detail, a family of such deterministic error guarantees is possible to extract
from the paper. The dependence on $\epsilon$ can be improved but the
dependence on $(1-\gamma)^{-1}$ is at least quadratic. The work of Mou et al.
(2020) provided an analysis for the Polyak–Ruppert averaged temporal
difference learning algorithm with linear function approximation in the
weighted $\ell_{2}$-norm, showing that the dominant stochastic error term
matches the stochastic lower bound proved in Proposition 1. However, the
sample complexity is suboptimal due to the presence of higher-order terms (see
Khamaru et al., 2021a, for simulations demonstrating this suboptimality), as
is the oracle complexity.

###  4.2 Variance-reduced fast temporal temporal difference algorithm

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

The VRFTD algorithm is formally presented in Algorithm 2, and we introduce the
basic idea of the algorithm below. First, it utilizes the idea of recentering
updates from VRTD with the operator $\widehat{g}(\widetilde{\theta})$ used in
each epoch. Second, in terms of iterate updating in the inner loop, it
involves an inner mini-batch that generates the averaged operator
$\widetilde{g}_{t}$, which allows the algorithm to be run with a much larger
stepsize. Finally, each iteration within an epoch involves an operator
extrapolation step (22). This is crucial to achieving the optimal
deterministic error (cf. the VRTD update (18)).

Algorithm 2 Variance-reduced Fast Temporal Difference Algorithm under i.i.d.
observations

Input: $\theta^{0}=\widehat{\theta}_{0}\in\mathbb{R}^{d}$, $\eta>0$,
$\lambda\geq 0$, $\\{\zeta_{t}\\}_{t=1}^{T}\geq 0$ and nonnegative integers
$m$, $\\{N_{k}\\}_{k=1}^{K}$.

for $k=1,\ldots,K$ do

Set $\theta_{0}=\theta_{1}=\widetilde{\theta}=\widehat{\theta}_{k-1}$. Collect
$N_{k}$ sample tuples
$\xi_{i}^{k}=(s_{i},s_{i}^{\prime},R(s_{i},s_{i}^{\prime}))$ from the i.i.d.
observation model.

Calculate
$\widehat{g}(\widetilde{\theta})=\tfrac{1}{N_{k}}\textstyle{\sum}_{i=1}^{N_{k}}\widetilde{g}(\widetilde{\theta},\xi_{i}^{k}).$

for $t=1,\ldots,T$ do

Collect $m$ sample tuples
$\xi_{j}^{t}=(s_{j},s_{j}^{\prime},R(s_{j},s_{j}^{\prime}))$ from the i.i.d.
observation model.

Calculate
$\widetilde{g}_{t}(\cdot)=\tfrac{1}{m}\textstyle{\sum}_{j=1}^{m}\widetilde{g}(\cdot,\xi_{j}^{t})$.

Denote
$\widetilde{F}_{t}(\theta_{t})=\widetilde{g}_{t}(\theta_{t})-\widetilde{g}_{t}(\widetilde{\theta})+\widehat{g}(\widetilde{\theta})$.
Set $\widetilde{F}_{0}(\theta_{0})=\widetilde{F}_{1}(\theta_{1})$. Let

$\theta_{t+1}=\theta_{t}-\eta\left[\widetilde{F}_{t}(\theta_{t})+\lambda\left(\widetilde{F}_{t}(\theta_{t})-\widetilde{F}_{t-1}(\theta_{t-1})\right)\right].$
(22)

end for

Output of the epoch:
$\widehat{\theta}_{k}=\frac{\sum_{t=1}^{T+1}\zeta_{t}\theta_{t}}{\sum_{t=1}^{T+1}\zeta_{t}}.$
(23)

end for

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

######  Theorem 3.

Fix the total number of epochs $K$ and a positive integer $N$. Assume that for
each epoch, the parameters $\eta$, $\lambda$, $m$, $N_{k}$, $T$ satisfy

$\eta\leq\tfrac{1}{4\beta(1+\gamma)},~{}\lambda=1,~{}T\geq\tfrac{32}{\mu(1-\gamma)\eta},~{}m\geq\max\left\\{1,\tfrac{256\eta\varsigma^{2}}{1-\gamma}\right\\}~{}\text{and}~{}N_{k}\geq\max\left\\{\tfrac{56\varsigma^{2}}{\mu(1-\gamma)^{2}},(\tfrac{3}{4})^{K-k}N\right\\}.$
(24)

Set the output of this epoch to be
$\widehat{v}_{k}:=\frac{\sum_{t=2}^{T+1}v_{t}}{T}$. Then for $\delta>0$,

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{K}-v^{*}\|_{\Pi}^{2}]$
$\displaystyle\leq(1+\delta)\mathcal{A}(M,\gamma)\inf_{v\in
S}\|v-v^{*}\|_{\Pi}^{2}$
$\displaystyle\quad+(1+\tfrac{1}{\delta})\left[\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|_{\Pi}^{2}+\tfrac{15}{N}\mathrm{trace}\left((I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\right)\right].$

See Section 7.5 for the detailed proof of this theorem.

In view of Theorem 3, the number of epochs required by the VRFTD method to
find a solution $\widehat{v}\in\mathbb{R}^{D}$, such that
${\mathbb{E}}[\|\widehat{v}-\bar{v}\|_{\Pi}^{2}]\leq\epsilon$ is bounded by
$\mathcal{O}\\{\log(\|v^{0}-\bar{v}\|^{2}_{\Pi}/\epsilon)\\}$. The total
number of samples used is $\textstyle{\sum}_{k=1}^{K}(mT+N_{k})$, which is
bounded on the order

$\displaystyle\underbrace{\tfrac{\beta}{(1-\gamma)\mu}\log(\tfrac{\|v^{0}-\bar{v}\|_{\Pi}^{2}}{\epsilon})}_{\text{deterministic
error}}+\underbrace{\tfrac{\varsigma^{2}}{(1-\gamma)^{2}\mu}\log(\tfrac{\|v^{0}-\bar{v}\|_{\Pi}^{2}}{\epsilon})+\tfrac{\mathrm{trace}\big{(}(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\big{)}}{\epsilon}}_{\text{stochastic
error}}.$ (25)

Similar to the VRTD algorithm, the VRFTD algorithm achieves optimal sample
complexity in terms of stochastic error. For the deterministic error, the
dependence on $1/(1-\gamma)$ matches the oracle complexity lower bound proved
in Theorem 1. Note that the term $\beta/\mu$ is the condition number of the
feature matrix in $\ell_{\Pi}$-norm. Therefore, for “well-conditioned” feature
matrices, the VRFTD algorithm achieves optimal oracle complexity444It is an
interesting open problem to prove an oracle complexity lower bound for policy
evaluation with linear function approximation having linear dependence on
$\beta/\mu$.. In summary, VRFTD is an accelerated and instance-optimal policy
evaluation algorithm, and answers the central question posed in this paper.

##  5 Algorithm for policy evaluation in the Markovian setting

Finally, we extend the VRFTD algorithm to the Markovian setting, noting in
passing that such an extension is also possible for the VRTD algorithm. The
challenge of Markovian noise stems from the presence of dependent data that
leads to biased samples. To control the bias caused by correlation, we need a
standard ergodicity assumption on the underlying Markov chain.

######  Assumption 3.

There exist constants $C_{P}>0$ and $\rho\in(0,1)$ such that

$\displaystyle\max_{s\in
S}\|\mathbb{P}(s_{t}=\cdot|s_{0}=s)-\pi\|_{\infty}\leq
C_{P}\cdot\rho^{t}\quad\text{ for all }\;t\in\mathbb{Z}_{+}.$ (26)

In other words, with the following definition of mixing time

$\displaystyle t_{\mathsf{mix}}:=\inf\\{t\in\mathbb{Z}_{+}\;\mid\;\max_{s\in
S}\|\mathbb{P}(s_{t}=\cdot|s_{0}=s)-\pi\|_{\infty}\leq 1/4\\},$

Assumption 3 guarantees that the mixing time is bounded as555Note that while
our choice of the constant $1/4$ in the definition is arbitrary, there is no
additional dependence on $\epsilon$ when accounting for the mixing time,
unlike in the assumptions made by (Bhandari et al., 2018, Eq. (21)).
$t_{\mathsf{mix}}\leq\frac{\log(4C_{P})}{\log(1/\rho)}$.

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

######  Lemma 1.

For every $t,\tau\in\mathbb{Z}_{+}$, with probability 1,

$\displaystyle\|{\mathbb{E}}[\widetilde{g}(\bar{\theta},\xi_{t+\tau})|\mathcal{F}_{t}]-g(\bar{\theta})\|_{2}\leq
C_{M}\cdot\rho^{\tau}\|\bar{v}-v^{*}\|_{\Pi}.$ (27)

where
$C_{M}:=\tfrac{C_{P}}{\sqrt{\min_{i\in[D]}\pi_{i}}}\|\Psi\|_{2}\|I-\gamma
P\|_{2}$ and $\mathcal{F}_{t}:=[\xi_{1},...,\xi_{t}]$.

See Appendix B.1 for a proof of this lemma. In words, Lemma 1 provides an
upper bound on the bias of the stochastic operator at the solution
$\bar{\theta}$ in terms of the approximation error, and the bound decays
exponentially with $\tau$.

######  Lemma 2.

For every $t,\tau\in\mathbb{Z}_{+}$ and
$\theta,\theta^{\prime}\in\mathbb{R}^{d}$, with probability 1,

$\displaystyle\|{\mathbb{E}}[\widetilde{g}(\theta,\xi_{t+\tau})|\mathcal{F}_{t}]-{\mathbb{E}}[\widetilde{g}(\theta^{\prime},\xi_{t+\tau})|\mathcal{F}_{t}]-[g(\theta)-g(\theta^{\prime})]\|_{2}\leq
C_{M}\cdot\rho^{\tau}\|v-v^{\prime}\|_{\Pi}.$ (28)

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

Algorithm 3 Variance-reduced Fast Temporal Difference Algorithm under
Markovian noise

Input: $\theta^{0}=\widehat{\theta}_{0}\in\mathbb{R}^{d}$, $\eta>0$,
$\lambda\geq 0$, $\\{\zeta_{t}\\}_{t=1}^{T}\geq 0$ and nonnegative integers
$m$, $m_{0}$, $n_{0}$, $\\{N_{k}\\}_{k=1}^{K}$.

for $k=1,\ldots,K$ do

Set $\theta_{1}=\widetilde{\theta}=\widehat{\theta}_{s-1}$. Collect $N_{k}$
successive samples $\xi_{i}^{k}:=(s_{i},s_{i+1},R(s_{i},s_{i+1}))$ from the
single Markov trajectory. Calculate
$\widehat{g}(\widetilde{\theta})=\tfrac{1}{N_{k}-n_{0}}\textstyle{\sum}_{i=n_{0}+1}^{N_{k}}\;\;\widetilde{g}(\widetilde{\theta},\xi_{i}^{k}).$
(29)

for $t=1,\ldots,T$ do

Collect $m$ successive samples
$\widehat{\xi}_{j}^{t}:=(s_{j},s_{j+1},R(s_{j},s_{j+1}))$ from the Markov
trajectory.

Calculate
$\widetilde{g}_{t}(\cdot)=\tfrac{1}{m-m_{0}}\textstyle{\sum}_{j=m_{0}+1}^{m}\widetilde{g}(\cdot,\widehat{\xi}_{j}^{t})$.

Let
$\widetilde{F}_{t}(\theta_{t})=\widetilde{g}_{t}(\theta_{t})-\widetilde{g}_{t}(\widetilde{\theta})+\widehat{g}(\widetilde{\theta})$
and set $\widetilde{F}_{0}(\theta_{0})=\widetilde{F}_{1}(\theta_{1})$. Let
$\theta_{t+1}=\theta_{t}-\eta\left[\widetilde{F}_{t}(\theta_{t})+\lambda\left(\widetilde{F}_{t}(\theta_{t})-\widetilde{F}_{t-1}(\theta_{t-1})\right)\right].$
(30)

end for

Output of the epoch:
$\widehat{\theta}_{k}=\frac{\sum_{t=1}^{T+1}\zeta_{t}\theta_{t}}{\sum_{t=1}^{T+1}\zeta_{t}}.$
(31)

end for

Before presenting our main convergence result for the VRFTD algorithm, we
first define the matrix $\bar{\Sigma}_{\mathsf{Mkv}}$, which is a covariance
matrix analog for the Markovian case (see Mou et al. (2021) and references
therein). Letting $\\{\widetilde{s}_{t}\\}_{t=-\infty}^{\infty}$ define a
sequence of samples obtained from a _stationary_ Markov trajectory, define

$\bar{\Sigma}_{\mathsf{Mkv}}:=\sum_{t=-\infty}^{\infty}B^{-\tfrac{1}{2}}~{}{\mathbb{E}}\left[\big{(}\widetilde{g}(\bar{\theta},\widetilde{\xi}_{t})-g(\bar{\theta})\big{)}\big{(}\widetilde{g}(\bar{\theta},\widetilde{\xi}_{0})-g(\bar{\theta})\big{)}^{\top}\right]B^{-\tfrac{1}{2}},$

where
$\widetilde{\xi}_{t}:=(\widetilde{s}_{t},\widetilde{s}_{t+1},R(\widetilde{s}_{t},\widetilde{s}_{t+1}))$.
This matrix is an infinite sum of matrices where one of the summands (when
$t=0$) is the matrix $\bar{\Sigma}_{\mathsf{iid}}$ defined in Eq. (13).

Similarly to the i.i.d. setting, the instance-dependent complexity of
Markovian linear stochastic approximation was shown in Mou et al. (2021) to be
governed by the trace of the matrix
$(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-T}$. To interpret this
functional, consider the special case in which the approximation error caused
by linear function approximation is $0$, i.e., $v^{*}=\bar{v}$. Let
$\widetilde{\mathcal{F}_{i}}$ denote the $\sigma$-field generated by samples
$\widetilde{\xi}_{0},...,\widetilde{\xi}_{i}$ and let
$\widetilde{\Pi}_{j}^{i}:=\operatorname{diag}\\{[\mathbb{P}(\widetilde{s}_{j}=1|\widetilde{s}_{i}),...,\mathbb{P}(\widetilde{s}_{j}=D|\widetilde{s}_{i})]\\}$
for $j\geq i$. Then we have

$\displaystyle{\mathbb{E}}\left\langle(I_{d}-M)^{-1}B^{-\frac{1}{2}}\big{(}\widetilde{g}(\bar{\theta},\widetilde{\xi}_{0})-g(\bar{\theta})\big{)},(I_{d}-M)^{-1}B^{-\frac{1}{2}}\big{(}\widetilde{g}(\bar{\theta},\widetilde{\xi}_{i})-g(\bar{\theta})\big{)}\right\rangle$
$\displaystyle={\mathbb{E}}\left\langle(I_{d}-M)^{-1}B^{-\frac{1}{2}}\big{(}\widetilde{g}(\bar{\theta},\widetilde{\xi}_{0})-g(\bar{\theta})\big{)},(I_{d}-M)^{-1}B^{-\frac{1}{2}}\big{(}{\mathbb{E}}[\widetilde{g}(\bar{\theta},\widetilde{\xi}_{i})|\widetilde{\mathcal{F}}_{0}]-g(\bar{\theta})\big{)}\right\rangle$
$\displaystyle={\mathbb{E}}\left\langle(I_{d}-M)^{-1}B^{-\frac{1}{2}}\big{(}\widetilde{g}(\bar{\theta},\widetilde{\xi}_{0})-g(\bar{\theta})\big{)},(I_{d}-M)^{-1}B^{-\frac{1}{2}}\Psi(\widetilde{\Pi}_{j}^{i}-\Pi)(\Psi^{\top}\bar{\theta}-\gamma
P\Psi^{\top}\bar{\theta}-r)\right\rangle=0,$

where the last equation follows from the fact that
$v^{*}=\bar{v}=\Psi^{\top}\bar{\theta}$ and the Bellman equation (1). Then

$\displaystyle\mathrm{trace}\left((I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-\top}\right)=\mathrm{trace}\left((I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\right).$

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

######  Theorem 4.

Fix the total number of epochs $K$ and a positive integer $N$. Consider an
integer $\tau$ satisfying
$\rho^{\tau}\leq\min\\{\frac{2(1-\rho)\varsigma}{3C_{M}},\frac{2(1-\rho)^{2}}{5C_{M}}\\}$.
Suppose the parameters $n_{0}$ and $m_{0}$ satisfy

$\displaystyle\rho^{n_{0}}\leq\tfrac{\min_{i\in[D]}\pi_{i}}{C_{P}},~{}\text{and}\quad\rho^{m_{0}}\leq\min\left\\{\tfrac{\min_{i\in[D]}\pi_{i}}{C_{P}},\tfrac{\sqrt{\mu}\eta\tau\varsigma^{2}(1-\rho)}{C_{M}}\right\\}.$
(32)

Assume that for each epoch $k\in[K]$, the parameters $\eta$, $\lambda$, $m$,
$N_{k}$, $T$ satisfy

$\displaystyle\eta$
$\displaystyle\leq\tfrac{1}{4\beta(1+\gamma)},~{}~{}~{}~{}\lambda=1,~{}~{}~{}~{}T\geq\tfrac{64}{\mu(1-\gamma)\eta},~{}~{}~{}~{}m-m_{0}\geq\max\left\\{1,\tfrac{792\eta(\tau+1)\varsigma^{2}}{1-\gamma}\right\\},$
$\displaystyle\rho^{N_{k}-n_{0}}\leq\tfrac{\tau(1-\rho)}{5C_{M}(N_{k}-n_{0})},~{}\text{and}~{}~{}~{}N_{k}-n_{0}\geq\left\\{\tfrac{206(\tau+1)\varsigma^{2}}{\mu(1-\gamma)^{2}},(\tfrac{3}{4})^{K-k}N\right\\}.$
(33)

Set the output of each epoch to be
$\widehat{v}_{k}:=\frac{\sum_{t=2}^{T+1}v_{t}}{T}$. Then the following results
hold.
(a) If $\bar{v}=v^{*}$, we have

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{K}-v^{*}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|_{\Pi}^{2}+\tfrac{30}{N}\cdot\mathrm{trace}\left((I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\right).$
(34)

(b) If $\bar{v}\neq v^{*}$, then for any $\delta>0$ we have

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{K}-v^{*}\|_{\Pi}^{2}]$
$\displaystyle\leq\left(1+\delta+\tfrac{18(\tau+1)(1+1/\delta)}{\mu(1-\gamma)^{2}N^{2}}\right)\cdot\mathcal{A}(M,\gamma)\cdot\inf_{v\in
S}\|v-v^{*}\|_{\Pi}^{2}$
$\displaystyle\quad+(1+\tfrac{1}{\delta})\left[\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|_{\Pi}^{2}+\tfrac{30}{N}\cdot\mathrm{trace}\left((I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-\top}\right)+\tfrac{\mathcal{H}}{N^{2}}\right],$
(35)

where
$\mathcal{H}:=(90\tau^{2}+18\tau+18)\cdot\mathrm{trace}\left((I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\right).$

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

##### Estimation of mixing time

From the conditions above, e.g., Ineq. (74), the parameters $\tau,n_{0},m_{0}$
scale linearly in the mixing time $t_{\mathsf{mix}}$ and logarithmically in
other problem parameters. As such, only some rough estimation of the mixing
time is sufficient, which has been the topic of active research. Nontrivial
confidence intervals for the reversible case can be found in Hsu et al.
(2019). There are also guarantees in the more challenging and prevalent case
when the underlying Markov chain is non-reversible (Wolfer and Kontorovich,
2019).

##### Sample complexities

We first consider the case when $\bar{v}=v^{*}$. In view of Ineq. (34) in
Theorem 4, the total number of samples required by the VRFTD method to find a
solution $\widehat{v}\in\mathbb{R}^{D}$, such that
${\mathbb{E}}[\|\widehat{v}-v^{*}\|_{\Pi}^{2}]\leq\epsilon$ is
$\textstyle{\sum}_{k=1}^{K}(mT+N_{k})$, which is bounded on the order666Note
that we omit the logarithmic dependence on the problem parameters, e.g.,
$\mu,\beta,\gamma$. Same for the following complexity.

$\displaystyle\underbrace{\tfrac{t_{\mathsf{mix}}\beta}{(1-\gamma)\mu}\log(\tfrac{\|v^{0}-\bar{v}\|_{\Pi}^{2}}{\epsilon})}_{\text{deterministic
error}}+\underbrace{\tfrac{t_{\mathsf{mix}}\varsigma^{2}}{(1-\gamma)^{2}\mu}\log(\tfrac{\|v^{0}-\bar{v}\|_{\Pi}^{2}}{\epsilon})+\tfrac{\mathrm{trace}\big{(}(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\big{)}}{\epsilon}}_{\text{stochastic
error}},$ (36)

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

When $\bar{v}\neq v^{*}$, the Markovian setting has a biased stochastic
operator at optimal solution $\bar{v}$, and a larger approximation error
$\|\bar{v}-v^{*}\|_{\Pi}^{2}$ caused by linear function approximation enlarges
the bias of the stochastic operator and consequently enlarges the dependence
on the approximation error in Eq. (4). Therefore, a natural stopping criterion
for the Markovian setting is to find a solution $\widehat{v}\in\mathbb{R}^{D}$
satisfying ${\mathbb{E}}[\|\widehat{v}-v^{*}\|_{\Pi}^{2}]\leq
c\|\bar{v}-v^{*}\|_{\Pi}^{2}+\epsilon$ for some absolute constant $c>0$. From
Ineq. (4), the total number of required samples
$\textstyle{\sum}_{k=1}^{K}(mT+N_{k})$ is bounded on the order

$\displaystyle\underbrace{\tfrac{t_{\mathsf{mix}}\beta}{(1-\gamma)\mu}\log(\tfrac{\|v^{0}-\bar{v}\|_{\Pi}^{2}}{\epsilon})}_{\text{deterministic
error}}+\underbrace{\tfrac{t_{\mathsf{mix}}\varsigma^{2}\log(\frac{\|v^{0}-\bar{v}\|_{\Pi}^{2}}{\epsilon})}{(1-\gamma)^{2}\mu}+\tfrac{\sqrt{{\mathcal{H}}}}{\sqrt{\epsilon}}+\tfrac{\mathrm{trace}\big{(}(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-\top}\big{)}}{\epsilon}}_{\text{stochastic
error}}.$ (37)

Note that, in this bound, the leading stochastic error matches the lower bound
proved in Mou et al. (2021), which can depend on the mixing time, but is
generally smaller than the product of the i.i.d. stochastic error and the
mixing time. These results show a delicate difference between how the mixing
time of the Markov chain enters the bound depending on whether the function
approximation is exact or not. It should be noted that, the stronger stopping
criterion, i.e., finding $\widehat{v}$ to satisfy
${\mathbb{E}}[\|\widehat{v}-\bar{v}\|_{\Pi}^{2}]\leq\epsilon$, can also be
applied in this setting. We can generate a similar sample complexity by
enlarging the constants $\tau$ and $N_{k}$ by an additive factor of
$\log(\|v^{*}-\bar{v}\|_{\Pi}^{2})$. However, given that the approximation
error is unavoidable and generally unknown, there is marginal benefit to using
this stronger stopping criterion.

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

##  6 Numerical experiments

In this section, we report numerical experiments for both VRTD and VRFTD,
comparing them against temporal difference learning (TD), conditional temporal
difference learning (CTD), and fast temporal difference learning (FTD)
(Kotsalis et al., 2020a, b). To generate a comprehensive performance profile,
we conduct experiments under both the i.i.d. and Markovian noise models.

###  6.1 The i.i.d. setting: A simple two-state construction

We consider a family of two-state MRPs inspired by the construction of Duan et
al. (2021). For a discount factor $\gamma\in(\tfrac{1}{2},1)$, the transition
kernel $P$ and reward vector $r$ are given by
$P=\begin{bmatrix}\tfrac{2\gamma-1}{\gamma}&\tfrac{1-\gamma}{\gamma}\\\
\tfrac{1-\gamma}{\gamma}&\tfrac{2\gamma-1}{\gamma}\end{bmatrix}$ and
$r=\begin{bmatrix}1\\\ -1\end{bmatrix}$. Clearly the transition kernel is
symmetric, thus the stationary distribution is $\pi=[0.5,0.5]$. For
simplicity, we choose the feature matrix
$\Psi=\operatorname{diag}([\sqrt{2},\sqrt{2}])$, which forms an orthonormal
basis under $\ell_{\Pi}$-norm. Assuming the i.i.d. model in which
$s_{i}\sim\pi$ and $s_{i}^{\prime}\sim P(\cdot|s_{i})$, it can be shown via
simple calculation that the stochastic error term is given by

$\displaystyle\mathrm{trace}\big{(}(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\big{)}=\tfrac{40}{81}\cdot\tfrac{2\gamma-1}{(1-\gamma)^{3}}.$
(38)

(a) TD algorithm

(b) FTD algorithm

(c) VRTD algorithm

(d) VRFTD algorithm

Figure 1:  Log-log plots of the squared $\ell_{\Pi}$-norm error versus
$1/(1-\gamma)$. Logarithms are to the natural base. The number of samples used
in each single experiment is $N=\lceil 5/(1-\gamma)^{2}\rceil$. Each point in
the plot is an average of $1000$ independent trials. The slope of the lower
bound is 1.

##### Instance-optimality

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

From subplots (a) and (b) of Figure 1, it is clear that the vanilla TD and FTD
algorithms with diminishing stepsizes (Kotsalis et al., 2020b) do not achieve
the lower bound calculated in equation (38). On the other hand, sub-plots (c)
and (d) show that the VRTD and VRFTD algorithms achieve the lower bound (38),
and that these behaviors are robust to the choice of stepsize parameters.
However, given their epoch-wise nature, the outputs of variance-reduced
algorithms are more volatile than TD and FTD. Another interesting observation
is that the accelerated algorithms—FTD and VRFTD—are less sensitive to
stepsize parameters. Our next set of experiments explores this further.

##### Ablation analysis of VRFTD

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

(a) Ablation study for OE (I)

(b) Ablation study for OE (II)

(c) Ablation study for mini-batching

Figure 2:  Log-log plots of the squared $\ell_{\Pi}$-norm error versus
$1/(1-\gamma)$. Logarithms are to the natural base. The number of samples used
in each single experiment is $N=\lceil 5/(1-\gamma)^{2}\rceil$. Each point in
the plot is an average of $1000$ independent trials.

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

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

###  6.2 The Markovian setting: 2D Grid World

Our experiments under the Markovian noise model are conducted on the 2D Grid
World environment. This is a classical problem in reinforcement learning with
finite state and action spaces. An agent realizes a positive reward when
reaching a predetermined goal and negative ones when going through “traps”.
The dimension of the state space is set to be $D=400$, among which we assign a
goal state (with reward $r=1$) and 30 traps (with reward $r=-0.2$). The
transition kernel is fixed as follows: With probability $0.95$, the agent
moves in a direction that points towards the goal and with probability $0.05$
in a random direction. Our goal is to compute the value function $v^{*}$—for
each possible initial state of the agent. We also incorporate linear function
approximation in these experiments. Specifically, we generate random features
with dimension $d=50$ to estimate the $D=400$ dimensional value function in
this problem.

Figure 3: Comparison of the algorithms for the 2D-Grid world example. From
left to right $\gamma$ is set to $0.99$, and $0.999$ respectively. In the
$y$-axis we report ratios in terms of the Euclidean norm $\|\cdot\|_{\Pi}$.

We test the performance of four algorithms, with the discount factor $\gamma$
set to $0.99$ and $0.999$. Figure 3 plots the normalized error in
$\ell_{\Pi}$-norm against the length of the trajectory. In both experiments,
the VRFTD algorithm exhibits the fastest convergence to the true value
function, thereby corroborating our theoretical results. Note that in both
experiments, the estimation errors do not converge to zero, because there is a
nontrivial error incurred by linear function approximation. Another salient
takeaway is the following: Closer to the optimal solution $v^{*}$, the
variance-reduced algorithms (VRTD/VRFTD) achieve faster convergence rate
compared to their counterparts that do not employ variance reduction.

##  7 Proofs

This section contains full proofs of our main results.

###  7.1 Preliminary supporting lemmas

We begin by establishing a few useful lemmas which are frequently used in the
proof of the main theorems. All the lemmas are simple to state and prove, and
provide the basis for several manipulations in our subsequent proofs.

######  Lemma 3.

For any $\theta\in\mathbb{R}^{d}$ and $v:=\Psi^{\top}\theta$, we have

$\displaystyle\|\theta\|_{2}=\|\Phi^{\top}\theta\|_{\Pi}\quad\text{and}\quad\|v\|_{\Pi}=\|B^{\frac{1}{2}}\theta\|_{2}.$
(39)

###### Proof.

The desired result follows immediately from $\Phi=B^{-\frac{1}{2}}\Psi$ and
$\Phi\Pi\Phi^{\top}=I_{d}$. ∎

######  Lemma 4.

For each vector $u\in\mathbb{R}^{D}$, we have

$\displaystyle\|Pu\|_{\Pi}\leq\|u\|_{\Pi}.$ (40)

###### Proof.

Based on the definition of $\|\cdot\|_{\Pi}$, we have

$\displaystyle\|Pu\|_{\Pi}^{2}=\textstyle{\sum}_{i=1}^{D}\pi_{i}(Pu)_{(i)}^{2}\overset{(i)}{\leq}\textstyle{\sum}_{i=1}^{D}\pi_{i}(\textstyle{\sum}_{j=1}^{D}P_{i,j}u_{(j)}^{2})=\textstyle{\sum}_{j=1}^{D}(\textstyle{\sum}_{i=1}^{D}\pi_{i}P_{i,j})u_{(j)}^{2}\overset{(ii)}{=}\|u\|_{\Pi}^{2},$

where step (i) follows from Jensen’s inequality and step (ii) follows from the
definition of the stationary distribution. ∎

######  Lemma 5.

For any $\theta,\theta^{\prime}\in\mathbb{R}^{d}$, we have

$\displaystyle\langle
g(\theta)-g(\theta^{\prime}),\theta-\theta^{\prime}\rangle\geq(1-\gamma)\|v-v^{\prime}\|_{\Pi}^{2}.$
(41)

###### Proof.

By the definition of $g(\cdot)$, we have

$\displaystyle\langle
g(\theta)-g(\theta^{\prime}),\theta-\theta^{\prime}\rangle$
$\displaystyle=\langle\Psi\Pi(I-\gamma
P)\Psi^{\top}(\theta-\theta^{\prime}),\theta-\theta^{\prime}\rangle$
$\displaystyle=\langle(I-\gamma
P)\Psi^{\top}(\theta-\theta^{\prime}),\Psi^{\top}(\theta-\theta^{\prime})\rangle_{\Pi}$
$\displaystyle\overset{(i)}{\geq}\|v-v^{\prime}\|_{\Pi}^{2}-\gamma\|P(v-v^{\prime})\|_{\Pi}\|v-v^{\prime}\|_{\Pi}$
$\displaystyle\overset{(ii)}{\geq}(1-\gamma)\|v-v^{\prime}\|_{\Pi}^{2},$

where step (i) follows from the Cauchy–Schwarz inequality, and step (ii)
follows from Lemma 4. ∎

######  Lemma 6.

For any $\theta,\theta^{\prime}\in\mathbb{R}^{d}$, we have

$\displaystyle\|g(\theta)-g(\theta^{\prime})\|_{2}\leq(1+\gamma)\sqrt{\beta}\|v-v^{\prime}\|_{\Pi}\leq(1+\gamma)\beta\|\theta-\theta^{\prime}\|_{2}.$
(42)

###### Proof.

Recall the definition of $\beta$ in Eq. (4). By the definition of $g(\cdot)$,
we have

$\displaystyle\|g(\theta)-g(\theta^{\prime})\|_{2}$
$\displaystyle=\|(\Psi\Pi\Psi^{\top}-\gamma\Psi\Pi
P\Psi^{\top})(\theta-\theta^{\prime})\|_{2}$
$\displaystyle\overset{(i)}{\leq}\sqrt{\beta}\|(\Phi\Pi\Psi^{\top}-\gamma\Phi\Pi
P\Psi^{\top})(\theta-\theta^{\prime})\|_{2}$
$\displaystyle\overset{(ii)}{=}\sqrt{\beta}\|(\Phi^{\top}\Phi\Pi-\gamma\Phi^{\top}\Phi\Pi
P)\Psi^{\top}(\theta-\theta^{\prime})\|_{\Pi}$ (43)

where step (i) follows from $\Psi=B^{\frac{1}{2}}\Phi$ and
$\|B^{\frac{1}{2}}\|_{2}\leq\sqrt{\beta}$, and step (ii) follows from Lemma 3.
Invoking the relationship $\Phi\Pi\Phi^{\top}=I_{d}$ gives us

$\displaystyle\|(\Phi^{\top}\Phi\Pi-\gamma\Phi^{\top}\Phi\Pi
P)\Psi^{\top}(\theta-\theta^{\prime})\|_{\Pi}=\|(I-\gamma\Pi_{\mathbb{S}}P)\Psi^{\top}(\theta-\theta^{\prime})\|_{\Pi}\overset{(i)}{\leq}(1+\gamma)\|v-v^{\prime}\|_{\Pi},$
(44)

where step (i) follows from
$\|\Pi_{\mathbb{S}}P\Psi^{\top}(\theta-\theta^{\prime})\|_{\Pi}\leq\|P\Psi^{\top}(\theta-\theta^{\prime})\|_{\Pi}\leq\|\Psi^{\top}(\theta-\theta^{\prime})\|_{\Pi}.$

Substituting Ineq. (44) into (7.1) yields the first inequality in the chain
(42). The second inequality in the chain (42) follows from the fact that

$\displaystyle\|\Psi(\theta-\theta^{\prime})\|_{\Pi}=\|B^{\frac{1}{2}}(\theta-\theta^{\prime})\|_{2}\leq\sqrt{\beta}\|\theta-\theta^{\prime}\|_{2},$

which completes the proof. ∎

######  Lemma 7.

For each vector $u\in\mathbb{R}^{D}$, we have

$\displaystyle\|(I-\gamma P)^{-1}u\|_{\Pi}\leq\tfrac{1}{1-\gamma}\|u\|_{\Pi}.$
(45)

###### Proof.

Applying the triangle inequality, we have

$\displaystyle\|u\|_{\Pi}$ $\displaystyle=\|(I-\gamma P)(I-\gamma
P)^{-1}u\|_{\Pi}$ $\displaystyle\geq\|(I-\gamma
P)^{-1}u\|_{\Pi}-\gamma\|P(I-\gamma P)^{-1}u\|_{\Pi}$
$\displaystyle\overset{(i)}{\geq}(1-\gamma)\|(I-\gamma P)^{-1}u\|_{\Pi}.$

Here, step (i) follows from Lemma 4. For an alternative proof, note that
$(1-\gamma)(I-\gamma P)^{-1}$ is a Markov transition matrix with stationary
distribution $\pi$, and the desired result follows from Lemma 4. ∎

######  Lemma 8.

For any vector $u\in\mathbb{R}^{D}$, we have

$\displaystyle(I-\Pi_{\mathbb{S}}\gamma
P)^{-1}\Pi_{\mathbb{S}}u=\Phi^{\top}(I_{d}-M)^{-1}\Phi\Pi u.$ (46)

###### Proof.

Note that $z:=(I-\Pi_{\mathbb{S}}\gamma
P)^{-1}\Pi_{\mathbb{S}}u\in{\mathbb{S}}$ since $z=\Pi_{\mathbb{S}}\gamma
Pz+\Pi_{\mathbb{S}}u$. Also, we have $z=\Pi_{\mathbb{S}}z=\Phi^{\top}\Phi\Pi
z$. Consequently,

$\Phi\Pi z=\gamma\Phi\Pi Pz+\Phi\Pi u=\Phi\gamma\Pi P\Phi^{\top}\Phi\Pi
z+\Phi\Pi u=M\Phi\Pi z+\Phi\Pi u.$

Rearranging yields $\Phi\Pi z=(I_{d}-M)^{-1}\Phi\Pi u$ and putting together
the pieces, we have $z=\Phi^{\top}\Phi\Pi z=\Phi^{\top}(I_{d}-M)^{-1}\Phi\Pi
u$. ∎

######  Lemma 9.

For each vector $\theta\in\mathbb{R}^{d}$, we have

$\displaystyle\|(I_{d}-M)^{-1}\theta\|_{2}\leq\tfrac{1}{1-\gamma}\|\theta\|_{2}.$
(47)

###### Proof.

Invoking Lemma 8, we have

$\displaystyle\|(I_{d}-M)^{-1}\theta\|_{2}$
$\displaystyle=\|\Phi^{\top}(I_{d}-M)^{-1}\theta\|_{\Pi}$
$\displaystyle=\|\Phi^{\top}(I_{d}-M)^{-1}\Phi\Pi\Phi^{\top}\theta\|_{\Pi}$
$\displaystyle\overset{(i)}{=}\|(I-\Pi_{\mathbb{S}}\gamma
P)^{-1}\Pi_{\mathbb{S}}\Phi^{\top}\theta\|_{\Pi}$
$\displaystyle=\|(I-\Pi_{\mathbb{S}}\gamma P)^{-1}\Phi^{\top}\theta\|_{\Pi},$

where step (i) follows from Lemma 8. At the same time, we have

$\displaystyle\|\theta\|_{2}=\|\Phi^{\top}\theta\|_{\Pi}$
$\displaystyle=\|(I-\Pi_{\mathbb{S}}\gamma P)(I-\Pi_{\mathbb{S}}\gamma
P)^{-1}\Phi^{\top}\theta\|_{\Pi}$
$\displaystyle\geq\|(I-\Pi_{\mathbb{S}}\gamma
P)^{-1}\Phi^{\top}\theta\|_{\Pi}-\|\Pi_{\mathbb{S}}\gamma
P(I-\Pi_{\mathbb{S}}\gamma P)^{-1}\Phi^{\top}\theta\|_{\Pi}$
$\displaystyle\overset{(i)}{\geq}\|(I-\Pi_{\mathbb{S}}\gamma
P)^{-1}\Phi^{\top}\theta\|_{\Pi}-\gamma\|P(I-\Pi_{\mathbb{S}}\gamma
P)^{-1}\Phi^{\top}\theta\|_{\Pi}$
$\displaystyle\overset{(ii)}{\geq}(1-\gamma)\|(I-\Pi_{\mathbb{S}}\gamma
P)^{-1}\Phi^{\top}\theta\|_{\Pi},$

where step (i) follows from the non-expansiveness of the projection, and step
(ii) follows from Lemma 4. Putting together the two pieces yields the result.
∎

###  7.2 Proof of Theorem 1

First, it is clear that the methods of this type are invariant to a
simultaneous shift of variables. The sequence of iterates for solving $G(v)=0$
starting from $v_{0}$ is just a shift of the sequence generated for solving
$G(v+v_{0})=0$ starting from the origin. Therefore, without loss of
generality, we assume $v_{0}=0$.

Now let us construct a specific instance $(P,r)$ to show the lower bound.
Consider the $D\times D$ matrix

$\displaystyle
P:=\begin{bmatrix}\tfrac{1}{2\gamma}&0&0&...&0&1-\tfrac{1}{2\gamma}\\\
1-\tfrac{1}{2\gamma}&\tfrac{1}{2\gamma}&0&...&0&0\\\
0&1-\tfrac{1}{2\gamma}&\tfrac{1}{2\gamma}&...&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&...&1-\tfrac{1}{2\gamma}&\tfrac{1}{2\gamma}\end{bmatrix},$ (48)

and note that

$\displaystyle I-\gamma
P=\begin{bmatrix}\tfrac{1}{2}&0&0&...&0&\tfrac{1}{2}-\gamma\\\
\tfrac{1}{2}-\gamma&\tfrac{1}{2}&0&...&0&0\\\
0&\tfrac{1}{2}-\gamma&\tfrac{1}{2}&...&0&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&...&\tfrac{1}{2}-\gamma&\tfrac{1}{2}\end{bmatrix}.$ (49)

Also define the $D$-dimensional vector

$\displaystyle
r:=[\gamma-\tfrac{1}{2}+(\tfrac{1}{2}-\gamma)(2\gamma-1)^{D},~{}0,~{}0,...,~{}0,~{}0]^{\top}$
(50)

The matrix $I-\gamma P$ is square and full rank, and it is straightforward to
verify that the unique solution of the linear equation $(I-\gamma P)v^{*}=r$
is

$\displaystyle(v^{*})_{(i)}=(2\gamma-1)^{i}\quad\text{for all }i\in[D].$ (51)

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

$\displaystyle\Pi=[\tfrac{1}{D},~{}\tfrac{1}{D},~{}\tfrac{1}{D},...,~{}\tfrac{1}{D}]^{\top}.$

Thus, we obtain

$\displaystyle\|v_{0}-v^{*}\|_{\Pi}^{2}=\tfrac{1}{D}\textstyle{\sum}_{i=1}^{D}(2\gamma-1)^{2i}=\tfrac{(2\gamma-1)^{2}[1-(2\gamma-1)^{2D}]}{D[1-(2\gamma-1)^{2}]}.$
(52)

Let $\mathbb{R}^{k,D}:=\\{v\in\mathbb{R}^{D}\;|\;v_{(i)}=0\text{ for all
}k+1\leq i\leq D\\}$ denote the set of all $D$-dimensional vectors lying in
the span of the first $k$ standard basis vectors. Since all entries of
$I-\gamma P$ below its subdiagonal are equal to $0$ and all entries of $r$
except for its first are equal to $0$, we conclude that
$v_{k}\in\mathbb{R}^{k,D}$. Therefore,

$\displaystyle\|v_{k}-v^{*}\|_{\Pi}^{2}\geq\tfrac{1}{D}\textstyle{\sum}_{i=k+1}^{D}(2\gamma-1)^{2i}=\tfrac{(2\gamma-1)^{2k+2}[1-(2\gamma-1)^{2D-2k}]}{D[1-(2\gamma-1)^{2}]}.$

If $D\gg k$ is such that
$\tfrac{1-(2\gamma-1)^{2D-2k}}{1-(2\gamma-1)^{2D}}\geq\tfrac{1}{2}$, then we
conclude that

$\displaystyle\|v_{k}-v^{*}\|_{\Pi}^{2}\geq\tfrac{1}{2}(2\gamma-1)^{2k}\|v_{0}-v^{*}\|_{\Pi}^{2},$

as desired.

###  7.3 Proof of Proposition 1

For the reader’s convenience, we begin by stating a version of the Hájek-Le
Cam local asymptotic minimax theorem.

######  Theorem 5.

Let $\\{\mathbb{P}_{\vartheta}\\}_{\vartheta\in\Theta}$ denote a family of
parametric models, quadratically mean differentiable with Fisher information
matrix $J_{\vartheta^{\prime}}$ . Fix some parameter $\vartheta\in\Theta$, and
consider a function $h:\Theta\rightarrow\mathbb{R}^{d}$ that is differentiable
at $\vartheta$. Then for any quasi-convex loss
$L:\mathbb{R}^{d}\rightarrow\mathbb{R}$, we have

$\displaystyle\lim_{c\rightarrow\infty}\lim_{N\rightarrow\infty}\inf_{\widehat{h}_{N}}\sup_{\vartheta^{\prime}:\|\vartheta^{\prime}-\vartheta\|_{2}\leq
c/\sqrt{N}}{\mathbb{E}}_{\vartheta^{\prime}}\left[L\big{(}\sqrt{N}(\widehat{h}_{N}-h(\vartheta^{\prime})\big{)}\right]={\mathbb{E}}[L(Z)],$
(53)

where the infimum is taken over all estimators $\widehat{h}_{N}$ that are
measurable functions of $N$ i.i.d. data points drawn from
$\mathbb{P}_{\vartheta^{\prime}}$ , and the expectation is taken over a
multivariate Gaussian

$Z\sim\mathcal{N}\left(0,\nabla
h(\vartheta)^{\top}J^{\dagger}_{\vartheta}\;\nabla h(\vartheta)\right).$

In our model, we have $\vartheta=(\omega,P,R)$ and
$h(\vartheta)=(\Psi\Omega\Psi^{\top}-\Psi\Omega\gamma
P\Psi^{\top})^{-1}\Psi\Omega r$. We set the loss function $L$ to be
$\|\cdot\|_{\widetilde{B}}^{2}$. Invoking Theorem 5 yields

$\displaystyle\mathfrak{M}(\vartheta)={\mathbb{E}}\left[\|Z\|_{\widetilde{B}}^{2}\right]\quad\text{where}\quad
Z\sim\mathcal{N}\left(0,\nabla
h(\vartheta)^{\top}J^{\dagger}_{\vartheta}\;\nabla h(\vartheta)\right).$ (54)

The covariance is explicitly computed in the following lemma.

######  Lemma 10.

We have

$\displaystyle\nabla h(\vartheta)^{\top}J^{\dagger}_{\vartheta}\;\nabla
h(\vartheta)=(\Psi\Omega\Psi^{\top}-\gamma\Psi\Omega
P\Psi^{\top})^{-1}\Sigma~{}(\Psi\Omega\Psi^{\top}-\gamma\Psi\Omega
P\Psi^{\top})^{-T},$ (55)

where
$\Sigma=\operatorname{cov}\big{[}\big{(}\langle\psi(s)-\gamma\psi(s^{\prime}),\bar{\theta}\rangle-R(s,s^{\prime})\big{)}\psi(s)\big{]}$
for $s\sim\omega,~{}s^{\prime}\sim\mathbb{P}(\cdot|s)$.

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

Recall the definition of $\widetilde{M}$ and $\widetilde{\Sigma}$ in the
statement of Proposition 1. By substituting equality (55) into (54) and
invoking the relation
$\widetilde{B}^{-\frac{1}{2}}\Psi\Omega\Psi^{\top}\widetilde{B}^{-\frac{1}{2}}=I_{d}$,
we obtain

$\widetilde{B}^{\frac{1}{2}}Z\sim\mathcal{N}\left(0,(I_{d}-\widetilde{M})^{-1}\widetilde{\Sigma}(I_{d}-\widetilde{M})^{-T}\big{)}\right),$

which completes the proof of the proposition.

###  7.4 Proof of Theorem 2

Let $\underline{\theta}$ satisfy
$g(\underline{\theta})-g(\widetilde{\theta})+\widehat{g}(\widetilde{\theta})=0$
and $\underline{v}=\Psi^{\top}\underline{\theta}$. Recalling the definition of
$\bar{v}$ from Eq. (2), the following lemma provides a bound for
$\|\underline{v}-\bar{v}\|_{\Pi}^{2}$. This bound is also valid for the VRFTD
algorithm in the i.i.d. setting.

######  Lemma 11.

Consider a single epoch with index $k\in[K]$. We have

$\displaystyle{\mathbb{E}}[\|\underline{v}-\bar{v}\|^{2}_{\Pi}]\leq\tfrac{2}{N_{k}}\mathrm{trace}\left((I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\right)+\tfrac{2\varsigma^{2}}{N_{k}(1-\gamma)^{2}\mu}{\mathbb{E}}\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}.$
(56)

See Appendix B.4 for the proof of this lemma.

Given Lemma 11, we can derive the following proposition which characterizes
the progress of the VRTD algorithm in a single epoch.

######  Proposition 2.

Consider a single epoch with index $k\in[K]$. Suppose that the parameters
$\eta$, $N_{k}$ and $T$ satisfy

$\eta\leq\min\\{\tfrac{(1-\gamma)}{2\beta(1+\gamma)^{2}},\tfrac{1-\gamma}{32\varsigma^{2}}\\},~{}~{}T\geq\tfrac{32}{\mu(1-\gamma)\eta},~{}~{}\text{and}~{}~{}N_{k}\geq\tfrac{38\varsigma^{2}}{\mu(1-\gamma)^{2}}.$
(57)

Set the output of this epoch to be
$\widehat{v}_{k}:=\frac{\sum_{t=1}^{T}\eta(1-\gamma)v_{t}+(1/\beta)v_{T+1}}{T\eta(1-\gamma)+(1/\beta)}$.
Then we have

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{2}{\mathbb{E}}[\|\widehat{v}_{k-1}-\bar{v}\|_{\Pi}^{2}]+\tfrac{5}{N_{k}}\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}.$
(58)

Taking Proposition 2 as given for the moment, let us complete the proof of the
theorem. The main idea is to bound the approximation error term
$\|\bar{v}-v^{*}\|_{\Pi}^{2}$ separately from the term
$\|\widehat{v}_{K}-\bar{v}\|_{\Pi}^{2}$. To bound
$\|\bar{v}-v^{*}\|_{\Pi}^{2}$, we use the instance-dependent upper bound in
Ineq. (6). We bound the term $\|\widehat{v}_{K}-\bar{v}\|_{\Pi}^{2}$ by using
Ineq. (58) as follows

$\displaystyle{\mathbb{E}}\|\widehat{v}_{K}-\bar{v}\|_{\Pi}^{2}$
$\displaystyle\leq\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|^{2}_{\Pi}+\textstyle{\sum}_{k=1}^{K}\tfrac{5}{2^{K-k}N_{k}}\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}$
$\displaystyle\overset{(i)}{=}\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|^{2}_{\Pi}+\textstyle{\sum}_{k=1}^{K}(\tfrac{2}{3})^{K-k}\tfrac{5}{N}\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}$
$\displaystyle\leq\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|^{2}_{\Pi}+\tfrac{15}{N}\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}.$
(59)

Here, step (i) from the condition that $N_{k}\geq(\tfrac{3}{4})^{K-k}N$ for
all $k\in[K]$. To conclude, we use Young’s inequality to obtain

$\displaystyle{\mathbb{E}}\|\widehat{v}_{K}-v^{*}\|_{\Pi}^{2}$
$\displaystyle\leq(1+\delta)\cdot{\mathbb{E}}\|\bar{v}-v^{*}\|_{\Pi}^{2}+(1+\tfrac{1}{\delta})\cdot{\mathbb{E}}\|\widehat{v}_{K}-\bar{v}\|_{\Pi}^{2}$
$\displaystyle\leq(1+\delta)\cdot\mathcal{A}(M,\gamma)\cdot\inf_{v\in
S}\|v-v^{*}\|_{\Pi}^{2}+(1+\tfrac{1}{\delta})\cdot{\mathbb{E}}\|\widehat{v}_{K}-\bar{v}\|_{\Pi}^{2}.$

####  7.4.1 Proof of Proposition 2

We first provide an upper bound on the term
${\mathbb{E}}\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}$. By the updates of
the VRTD algorithm, we have

$\displaystyle\theta_{t+1}-\underline{\theta}=\theta_{t}-\eta\big{(}\widetilde{g}(\theta_{t},\xi_{t})-\widetilde{g}(\widetilde{\theta},\xi_{t})+\widehat{g}(\widetilde{\theta})\big{)}-\underline{\theta}.$

Now write

$\displaystyle{\mathbb{E}}\|\theta_{t+1}-\underline{\theta}\|_{2}^{2}$
$\displaystyle\overset{(i)}{=}{\mathbb{E}}\|\theta_{t}-\underline{\theta}-\eta\big{(}g(\theta_{t})-g(\underline{\theta})\big{)}+\eta\big{(}g(\theta_{t})-\widetilde{g}(\theta_{t},\xi_{t})-g(\widetilde{\theta})+\widetilde{g}(\widetilde{\theta},\xi_{t})\big{)}\|_{2}^{2}$
$\displaystyle={\mathbb{E}}\|\theta_{t}-\underline{\theta}\|_{2}^{2}-2\eta{\mathbb{E}}\langle
g(\theta_{t})-g(\underline{\theta}),\theta_{t}-\underline{\theta}\rangle$
$\displaystyle\quad+\eta^{2}{\mathbb{E}}\|g(\theta_{t})-g(\underline{\theta})\|_{2}^{2}+\eta^{2}{\mathbb{E}}\|g(\theta_{t})-\widetilde{g}(\theta_{t},\xi_{t})-g(\widetilde{\theta})+\widetilde{g}(\widetilde{\theta},\xi_{t})\|_{2}^{2}$
$\displaystyle\overset{(ii)}{\leq}{\mathbb{E}}\|\theta_{t}-\underline{\theta}\|_{2}^{2}+\big{(}\eta^{2}\beta(1+\gamma)^{2}-2\eta(1-\gamma)\big{)}{\mathbb{E}}\|v_{t}-\underline{v}\|_{\Pi}^{2}+\eta^{2}\varsigma^{2}{\mathbb{E}}\|\widetilde{v}-v_{t}\|_{\Pi}^{2}$

where step (i) follows from the fact that
$\widehat{g}(\widetilde{\theta})=g(\widetilde{\theta})-g(\underline{\theta})$,
and step (ii) follows from using Lemma 5 to bound the second term, Lemma 6 to
bound the third term and Assumption 2 to bound the last term. Using Young’s
inequality then yields

$\displaystyle{\mathbb{E}}\|\theta_{t+1}$
$\displaystyle-\underline{\theta}\|_{2}^{2}$
$\displaystyle\leq{\mathbb{E}}\|\theta_{t}-\underline{\theta}\|_{2}^{2}+\big{(}\eta^{2}\beta(1+\gamma)^{2}+2\eta^{2}\varsigma^{2}-2\eta(1-\gamma)\big{)}{\mathbb{E}}\|v_{t}-\underline{v}\|_{\Pi}^{2}+2\eta^{2}\varsigma^{2}{\mathbb{E}}\|\widetilde{v}-\underline{v}\|_{\Pi}^{2}$
$\displaystyle\overset{(i)}{\leq}{\mathbb{E}}\|\theta_{t}-\underline{\theta}\|_{2}^{2}-\eta(1-\gamma){\mathbb{E}}\|v_{t}-\underline{v}\|_{\Pi}^{2}+4\eta^{2}\varsigma^{2}{\mathbb{E}}\|\bar{v}-\widetilde{v}\|_{\Pi}^{2}+4\eta^{2}\varsigma^{2}{\mathbb{E}}\|\underline{v}-\bar{v}\|_{\Pi}^{2},$
(60)

where step (i) follows from the assumption
$\eta\leq\min\\{\tfrac{1-\gamma}{2\beta(1+\gamma)},\tfrac{1-\gamma}{4\varsigma^{2}}\\}$
and a further application of Young’s inequality. By taking a telescopic sum of
Ineq. (7.4.1) from $t=1$ to $T$, we obtain

$\displaystyle\textstyle{\sum}_{t=1}^{T}\eta(1-\gamma){\mathbb{E}}[\|v_{t}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle+{\mathbb{E}}[\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}]$
$\displaystyle\leq{\mathbb{E}}[\|\theta_{1}-\underline{\theta}\|_{2}^{2}]+4T\eta^{2}\varsigma^{2}\big{(}{\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}]+{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]\big{)}.$

Utilizing the inequalities
$\|\theta\|_{2}^{2}\geq\tfrac{1}{\beta}\|B^{\frac{1}{2}}\theta\|_{2}^{2}=\tfrac{1}{\beta}\|v\|_{\Pi}^{2}$
and
$\|\theta\|_{2}^{2}\leq\tfrac{1}{\mu}\|B^{\frac{1}{2}}\theta\|_{2}^{2}=\tfrac{1}{\mu}\|v\|_{\Pi}^{2}$,
we then have

$\displaystyle\textstyle{\sum}_{t=1}^{T}\eta(1-\gamma){\mathbb{E}}[\|v_{t}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle+\tfrac{1}{\beta}{\mathbb{E}}[\|v_{T+1}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{\mu}{\mathbb{E}}[\|v_{1}-\underline{v}\|_{\Pi}^{2}]+4T\eta^{2}\varsigma^{2}({\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}]+{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]).$

By definition, we have
$\widehat{v}_{k}=\frac{\sum_{t=1}^{T}\eta(1-\gamma)v_{t}+(1/\beta)v_{T+1}}{T\eta(1-\gamma)+(1/\beta)}$;
noting in addition that $v_{1}=\widetilde{v}$ then yields

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{T\eta(1-\gamma)+(1/\beta)}\big{[}\textstyle{\sum}_{t=1}^{T}\eta(1-\gamma){\mathbb{E}}[\|v_{t}-\underline{v}\|_{\Pi}^{2}]+{\mathbb{E}}[\|v_{T+1}-\underline{v}\|_{\Pi}^{2}]\big{]}$
$\displaystyle\leq\tfrac{1}{T\eta(1-\gamma)+(1/\beta)}\big{[}\tfrac{1}{\mu}{\mathbb{E}}[\|\widetilde{v}-\underline{v}\|_{\Pi}^{2}]+4T\eta^{2}\varsigma^{2}({\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}]+{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}])\big{]}$
$\displaystyle\leq\tfrac{2/\mu+4T\eta^{2}\varsigma^{2}}{T\eta(1-\gamma)}({\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}]+{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]).$
(61)

Combining the bound on $\|\underline{v}-v^{*}\|_{\Pi}^{2}$ shown in Lemma 11
and the bound on $\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}$ from Ineq.
(7.4.1), we obtain

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}$
$\displaystyle-\bar{v}\|_{\Pi}^{2}]$ $\displaystyle\leq
2{\mathbb{E}}[\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}]+2{\mathbb{E}}[\|\bar{v}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{4/\mu+8T\eta^{2}\varsigma^{2}}{T\eta(1-\gamma)}{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]+\left(\tfrac{4/\mu+8T\eta^{2}\varsigma^{2}}{T\eta(1-\gamma)}+2\right){\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{4/\mu+8T\eta^{2}\varsigma^{2}}{T\eta(1-\gamma)}{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\quad+\left(\tfrac{4/\mu+8T\eta^{2}\varsigma^{2}}{T\eta(1-\gamma)}+2\right)\big{(}\tfrac{2\varsigma^{2}}{\mu(1-\gamma)^{2}N_{k}}{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]+\tfrac{2}{N_{k}}\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}\big{)}.$
(62)

The final term takes the desired form; we are left with bounding the first two
terms in the display above. First, recall by assumption that
$\eta\leq\tfrac{1-\gamma}{32\varsigma^{2}}$, and
$T\geq\tfrac{32}{\mu(1-\gamma)\eta}$, so that

$\displaystyle\tfrac{4/\mu+8T\eta^{2}\varsigma^{2}}{T\eta(1-\gamma)}\leq\tfrac{\frac{4}{\mu}+\frac{256\eta\varsigma^{2}}{\mu(1-\gamma)}}{32/\mu}\leq\tfrac{3}{8}.$
(63)

Second, utilizing the condition
$N_{k}\geq\tfrac{38\varsigma^{2}}{\mu(1-\gamma)^{2}}$, we have

$\displaystyle\left(\tfrac{4/\mu+8T\eta^{2}\varsigma^{2}}{T\eta(1-\gamma)}+2\right)\tfrac{2\varsigma^{2}}{\mu(1-\gamma)^{2}N_{k}}\leq\tfrac{1}{8}.$
(64)

Substituting the bounds (63) and (64) into inequality (7.4.1), we obtain

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{2}{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]+\tfrac{5}{N_{k}}\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}\big{)},$

which completes the proof of Proposition 2.

###  7.5 Proof of Theorem 3

The structure of the proof is similar to the analysis of VRTD. We first state
a proposition that characterizes the progress of the VRFTD algorithm in a
single epoch.

######  Proposition 3.

Assume that for each epoch $k\in[K]$, the parameter $\eta$, $\lambda$, $m$,
$N_{k}$ and $T$ satisfy

$\eta\leq\tfrac{1}{4\beta(1+\gamma)},~{}\lambda=1,~{}T\geq\tfrac{32}{\mu(1-\gamma)\eta},~{}m\geq\max\\{1,\tfrac{256\eta\varsigma^{2}}{1-\gamma}\\},~{}\text{and}~{}N_{k}\geq\tfrac{56\varsigma^{2}}{\mu(1-\gamma)^{2}}.$
(65)

Set the output of this epoch to be
$\widehat{v}_{k}:=\frac{\sum_{t=2}^{T+1}v_{t}}{T}$. Then we have

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{2}{\mathbb{E}}[\|\widehat{v}_{k-1}-\bar{v}\|_{\Pi}^{2}]+\tfrac{5}{N_{k}}\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}.$
(66)

Taking Proposition 3 as given, the proof of Theorem 3 follows exactly as the
proof of Theorem 2 in Section 7.4. Therefore, to complete the proof, it
remains to prove Proposition 3, which we do in the following subsection.

####  7.5.1 Proof of Proposition 3

The proof structure parallels that of Proposition 2: The main idea is to bound
the terms $\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}$ and
$\|\underline{v}-\bar{v}\|_{\Pi}^{2}$ separately and then use Young’s
inequality to bound $\|\widehat{v}_{k}-\bar{v}\|_{\Pi}^{2}$. As alluded to
before, the bound proved on $\|\underline{v}-\bar{v}\|_{\Pi}^{2}$ in Lemma 11
is still valid. It remains to show a bound on the term
$\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}$.

From step (22) of Algorithm 2, we have

$\displaystyle\theta_{t}-\underline{\theta}=\theta_{t+1}+\eta\left[\widetilde{F}_{t}(\theta_{t})+\lambda\left(\widetilde{F}_{t}(\theta_{t})-\widetilde{F}_{t-1}(\theta_{t-1})\right)\right]-\underline{\theta}.$

As before, we like to obtain a one-step inequality that will allow us to take
a telescopic sum. Toward that end, note that some manipulation yields

$\displaystyle\|\theta_{t}-\underline{\theta}\|_{2}^{2}$
$\displaystyle=\|\theta_{t+1}-\underline{\theta}\|_{2}^{2}+2\eta\langle\widetilde{F}_{t}(\theta_{t})+\lambda[\widetilde{F}_{t}(\theta_{t})-\widetilde{F}_{t-1}(\theta_{t-1})],\theta_{t+1}-\underline{\theta}\rangle+\|\theta_{t+1}-\theta_{t}\|_{2}^{2}$
$\displaystyle=\|\theta_{t+1}-\underline{\theta}\|_{2}^{2}+\|\theta_{t+1}-\theta_{t}\|_{2}^{2}+2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1}),\theta_{t+1}-\underline{\theta}\rangle$
$\displaystyle\quad-2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1})-\widetilde{F}_{t}(\theta_{t}),\theta_{t+1}-\underline{\theta}\rangle+2\eta\lambda\langle\widetilde{F}_{t}(\theta_{t})-\widetilde{F}_{t-1}(\theta_{t-1}),\theta_{t}-\underline{\theta}\rangle$
$\displaystyle\quad+2\eta\lambda\langle\widetilde{F}_{t}(\theta_{t})-\widetilde{F}_{t-1}(\theta_{t-1}),\theta_{t+1}-\theta_{t}\rangle.$

Defining $\widetilde{\Delta}_{t}(\theta):=\widetilde{g}_{t}(\theta)-g(\theta)$
for convenience and recalling the definition of $\widetilde{F}_{t}$ in
Algorithm 2, we can then write

$\displaystyle\|$ $\displaystyle\theta_{t}-\underline{\theta}\|_{2}^{2}$
$\displaystyle=\|\theta_{t+1}-\underline{\theta}\|_{2}^{2}+\|\theta_{t+1}-\theta_{t}\|_{2}^{2}+2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1}),\theta_{t+1}-\underline{\theta}\rangle-2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1})-\widetilde{F}_{t}(\theta_{t}),\theta_{t+1}-\underline{\theta}\rangle$
$\displaystyle\qquad\qquad+2\eta\lambda\langle\widetilde{F}_{t}(\theta_{t})-\widetilde{F}_{t-1}(\theta_{t-1}),\theta_{t}-\underline{\theta}\rangle+2\eta\lambda\langle
g(\theta_{t})-g(\theta_{t-1}),\theta_{t+1}-\theta_{t}\rangle$
$\displaystyle\qquad\qquad+2\eta\lambda\langle\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})-\widetilde{\Delta}_{t-1}(\theta_{t-1})+\widetilde{\Delta}_{t-1}(\widetilde{\theta}),\theta_{t+1}-\theta_{t}\rangle.$

Summing both sides of the equation above from $t=1$ to $T$ and noting that
$\lambda=1$, $\theta_{1}=\theta_{0}$, and
$\widetilde{F}_{0}(\theta_{0})=\widetilde{F}_{1}(\theta_{1})$, we have

$\displaystyle\|\theta_{1}-\underline{\theta}\|_{2}^{2}$
$\displaystyle=\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}-2\eta\langle\widetilde{F}_{T+1}(\theta_{T+1})-\widetilde{F}_{T}(\theta_{T}),\theta_{T+1}-\underline{\theta}\rangle+Q_{1}$
$\displaystyle\qquad\qquad\qquad+\textstyle{\sum}_{t=1}^{T}2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1}),\theta_{t+1}-\underline{\theta}\rangle,$
(67)

where we have defined

$\displaystyle Q_{1}$
$\displaystyle:=\textstyle{\sum}_{t=1}^{T}\|\theta_{t+1}-\theta_{t}\|_{2}^{2}+\textstyle{\sum}_{t=1}^{T}2\eta\langle
g(\theta_{t})-g(\theta_{t-1}),\theta_{t+1}-\theta_{t}\rangle$
$\displaystyle\qquad\qquad\qquad+\textstyle{\sum}_{t=1}^{T}2\eta\langle\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})-\widetilde{\Delta}_{t-1}(\theta_{t-1})+\widetilde{\Delta}_{t-1}(\widetilde{\theta}),\theta_{t+1}-\theta_{t}\rangle$

for convenience. Noting that it suffices to lower bound $Q$, write

$\displaystyle Q_{1}$
$\displaystyle\overset{(i)}{\geq}\textstyle{\sum}_{t=1}^{T}\|\theta_{t+1}-\theta_{t}\|_{2}^{2}-\textstyle{\sum}_{t=1}^{T}2\eta(1+\gamma)\beta\|\theta_{t}-\theta_{t-1}\|_{2}\|\theta_{t+1}-\theta_{t}\|_{2}$
$\displaystyle\qquad\qquad-\textstyle{\sum}_{t=1}^{T}2\eta\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})-\widetilde{\Delta}_{t-1}(\theta_{t-1})+\widetilde{\Delta}_{t-1}(\widetilde{\theta})\|_{2}\|\theta_{t+1}-\theta_{t}\|_{2}$
$\displaystyle=\tfrac{1}{4}\|\theta_{T+1}-\theta_{T}\|_{2}^{2}$
$\displaystyle\qquad+\textstyle{\sum}_{t=1}^{T}\left[\tfrac{1}{4}\|\theta_{t+1}-\theta_{t}\|_{2}^{2}+\tfrac{1}{4}\|\theta_{t}-\theta_{t-1}\|_{2}^{2}-2\eta(1+\gamma)\beta\|\theta_{t}-\theta_{t-1}\|_{2}\|\theta_{t+1}-\theta_{t}\|_{2}\right]$
$\displaystyle\qquad+\textstyle{\sum}_{t=1}^{T}\left[\tfrac{1}{2}\|\theta_{t+1}-\theta_{t}\|_{2}^{2}-2\eta\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})-\widetilde{\Delta}_{t-1}(\theta_{t-1})+\widetilde{\Delta}_{t-1}(\widetilde{\theta})\|_{2}\|\theta_{t+1}-\theta_{t}\|_{2}\right]$
$\displaystyle\overset{(ii)}{\geq}\tfrac{1}{4}\|\theta_{T+1}-\theta_{T}\|_{2}^{2}-\textstyle{\sum}_{t=1}^{T}2\eta^{2}\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})-\widetilde{\Delta}_{t-1}(\theta_{t-1})+\widetilde{\Delta}_{t-1}(\widetilde{\theta})\|_{2}^{2}$
$\displaystyle\overset{(iii)}{\geq}\tfrac{1}{4}\|\theta_{T+1}-\theta_{T}\|_{2}^{2}-\textstyle{\sum}_{t=1}^{T}4\eta^{2}(\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})\|_{2}^{2}+\|\widetilde{\Delta}_{t-1}(\widetilde{\theta})-\widetilde{\Delta}_{t-1}(\theta_{t-1})\|_{2}^{2}),$
(68)

where step (i) follows from the Cauchy–Schwarz inequality and Lemma 6, step
(ii) follows from the fact that $4\eta(1+\gamma)\beta\leq 1$ and Young’s
inequality, and step (iii) follows from Young’s inequality. Substituting Ineq.
(7.5.1) into Eq. (7.5.1), we obtain

$\displaystyle\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}+\tfrac{1}{4}\|\theta_{T+1}$
$\displaystyle-\theta_{T}\|_{2}^{2}+\textstyle{\sum}_{t=1}^{T}2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1}),\theta_{t+1}-\underline{\theta}\rangle$
$\displaystyle\leq\|\theta_{1}-\underline{\theta}\|_{2}^{2}+2\eta\langle\widetilde{F}_{T+1}(\theta_{T+1})-\widetilde{F}_{T}(\theta_{T}),\theta_{T+1}-\underline{\theta}\rangle$
$\displaystyle\qquad+\textstyle{\sum}_{t=1}^{T}4\eta^{2}(\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})\|_{2}^{2}+\|\widetilde{\Delta}_{t-1}(\widetilde{\theta})-\widetilde{\Delta}_{t-1}(\theta_{t-1})\|_{2}^{2}).$

Rearranging terms and invoking the definition of the function
$\widetilde{F}(\cdot)$, we have

$\displaystyle\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}+\textstyle{\sum}_{t=1}^{T}2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1}),\theta_{t+1}-\underline{\theta}\rangle$
$\displaystyle\quad\leq\|\theta_{1}-\underline{\theta}\|_{2}^{2}+4\eta^{2}\big{(}\textstyle{\sum}_{t=1}^{T}2\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})\|_{2}^{2}-\|\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T})\|_{2}^{2}\big{)}\underbrace{-\tfrac{1}{4}\|\theta_{T+1}-\theta_{T}\|_{2}^{2}}_{R_{1}}$
$\displaystyle\quad\qquad+\underbrace{2\eta\langle\widetilde{\Delta}_{T+1}(\theta_{T+1})-\widetilde{\Delta}_{T+1}(\widetilde{\theta})+\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T}),\theta_{T+1}-\underline{\theta}\rangle}_{R_{2}}$
$\displaystyle\quad\qquad\qquad+\underbrace{2\eta\langle
g(\theta_{T+1})-g(\theta_{T}),\theta_{T+1}-\underline{\theta}\rangle}_{R_{3}}.$
(69)

Using the shorthand

$Q_{2}:=2\eta\langle\widetilde{\Delta}_{T+1}(\theta_{T+1})-\widetilde{\Delta}_{T+1}(\widetilde{\theta}),\theta_{T+1}-\underline{\theta}\rangle+2\eta\langle\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T}),\theta_{T}-\underline{\theta}\rangle,$

note that the final three terms on the RHS of Eq. (69) can be written as

$\displaystyle R_{1}+R_{2}+R_{3}$ $\displaystyle\leq
Q_{2}+\big{(}2\eta\langle\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T}),\theta_{T+1}-\theta_{T}\rangle-\tfrac{1}{8}\|\theta_{T+1}-\theta_{T}\|_{2}^{2}\big{)}$
$\displaystyle\qquad\qquad+\big{(}2\eta\langle
g(\theta_{T+1})-g(\theta_{T}),\theta_{T+1}-\underline{\theta}\rangle-\tfrac{1}{8}\|\theta_{T+1}-\theta_{T}\|_{2}^{2}\big{)}$
$\displaystyle\overset{(i)}{\leq}Q_{2}+\big{(}2\eta\|\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T})\|_{2}\|\theta_{T+1}-\theta_{T}\|_{2}-\tfrac{1}{8}\|\theta_{T+1}-\theta_{T}\|_{2}^{2}\big{)}$
$\displaystyle\qquad\qquad+\big{(}2\eta(1+\gamma)\beta\|\theta_{T+1}-\theta_{T}\|_{2}\|\theta_{T+1}-\underline{\theta}\|_{2}-\tfrac{1}{8}\|\theta_{T+1}-\theta_{T}\|_{2}^{2}\big{)},$

where step (i) follows from the Cauchy–Schwarz inequality and Lemma 6.
Applying Young’s inequality and putting together the pieces, we have

$\displaystyle\big{(}1-8\eta^{2}\beta^{2}$
$\displaystyle(1+\gamma)^{2}\big{)}\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}+\textstyle{\sum}_{t=1}^{T}2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1}),\theta_{t+1}-\underline{\theta}\rangle$
$\displaystyle\leq$
$\displaystyle\|\theta_{1}-\underline{\theta}\|_{2}^{2}+\textstyle{\sum}_{t=1}^{T}8\eta^{2}\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})\|_{2}^{2}+4\eta^{2}\|\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T})\|_{2}^{2}+Q_{2}.$
(70)

Note that the Ineq. (7.5.1) holds pointwise; we now take expectation on both
sides, and note that ${\mathbb{E}}[Q_{2}]=0$. We now use Assumption 2 to upper
bound on the RHS of Ineq. (7.5.1) and Lemma 5 to lower bound on the LHS of
Ineq. (7.5.1). Putting these together with the fact that
$\|\theta_{1}-\underline{\theta}\|_{2}^{2}\leq\tfrac{1}{\mu}\|v_{1}-\underline{v}\|_{\Pi}^{2}$,
we obtain

$\displaystyle\big{(}1-8\eta^{2}\beta^{2}(1+\gamma)^{2}\big{)}{\mathbb{E}}[\|\theta_{T+1}$
$\displaystyle-\underline{\theta}\|_{2}^{2}]+\textstyle{\sum}_{t=2}^{T+1}2\eta(1-\gamma){\mathbb{E}}[\|v_{t}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{\mu}{\mathbb{E}}[\|v_{1}-\underline{v}\|_{\Pi}^{2}]+\tfrac{4\eta^{2}\varsigma^{2}}{m}{\mathbb{E}}[\|v_{T}-\widetilde{v}\|_{\Pi}^{2}]+\textstyle{\sum}_{t=1}^{T}\tfrac{8\eta^{2}\varsigma^{2}}{m}{\mathbb{E}}[\|v_{t}-\widetilde{v}\|_{\Pi}^{2}].$

We now utilize the relation $\|v_{t}-\widetilde{v}\|_{\Pi}^{2}\leq
2\|v_{t}-\underline{v}\|_{\Pi}^{2}+2\|\widetilde{v}-\underline{v}\|_{\Pi}^{2}$
in conjunction with

$\displaystyle v_{1}=\widetilde{v}\quad\text{ and }\quad
8\eta^{2}(1+\gamma)^{2}\beta^{2}\leq 1.$

Rearranging terms then yields

$\displaystyle\textstyle{\sum}_{t=2}^{T+1}\big{(}2\eta(1-\gamma)-\tfrac{24\eta^{2}\varsigma^{2}}{m}\big{)}{\mathbb{E}}[\|v_{t}-\underline{v}\|_{\Pi}^{2}]\leq\big{(}\tfrac{1}{\mu}+\tfrac{(16T+40)\eta^{2}\varsigma^{2}}{m}\big{)}{\mathbb{E}}[\|\widetilde{v}-\underline{v}\|_{\Pi}^{2}].$

Invoking the bounds $m\geq\tfrac{256\eta\varsigma^{2}}{1-\gamma}$, and
$\widehat{v}_{k}=\tfrac{\sum_{t=2}^{T+1}v_{t}}{T}$, we have

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{\sum_{t=2}^{T+1}\eta(1-\gamma){\mathbb{E}}[\|v_{t}-\underline{v}\|_{\Pi}^{2}]}{T\eta(1-\gamma)}\leq\tfrac{1/\mu+(16T+40)\eta^{2}\varsigma^{2}/m}{T\eta(1-\gamma)}{\mathbb{E}}[\|\widetilde{v}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1/\mu+(16T+40)\eta^{2}\varsigma^{2}/m}{T\eta(1-\gamma)}(2{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]+2{\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}]).$

Applying Young’s inequality and Lemma 11, we obtain

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\quad\leq
2{\mathbb{E}}[\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}]+2{\mathbb{E}}[\|\bar{v}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\quad\leq\tfrac{4/\mu+(64T+160)\eta^{2}\varsigma^{2}/m}{T\eta(1-\gamma)}{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]+\left(\tfrac{4/\mu+(64T+160)\eta^{2}\varsigma^{2}/m}{T\eta(1-\gamma)}+2\right){\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\quad\leq\tfrac{4/\mu+(64T+160)\eta^{2}\varsigma^{2}/m}{T\eta(1-\gamma)}{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\quad\quad+\left(\tfrac{4/\mu+(64T+160)\eta^{2}\varsigma^{2}/m}{T\eta(1-\gamma)}+2\right)\bigg{(}\tfrac{2\varsigma^{2}}{\mu(1-\gamma)^{2}N_{k}}{\mathbb{E}}[\|(\widetilde{v}-\bar{v})\|_{\Pi}^{2}]$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad+\tfrac{2}{N_{k}}\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}\bigg{)}.$
(71)

As before, we bound the first two terms by assumption. First, note that
$\eta\leq\tfrac{1}{4\beta(1+\gamma)}$, $T\geq\tfrac{32}{\mu(1-\gamma)\eta}$
and $m\geq\max\\{1,\tfrac{256\eta\varsigma^{2}}{1-\gamma}\\}$ whereby

$\displaystyle\tfrac{4/\mu+(64T+160)\eta^{2}\varsigma^{2}/m}{T\eta(1-\gamma)}\leq\tfrac{4/\mu+T\eta(1-\gamma)/4+5\eta(1-\gamma)/8}{T\eta(1-\gamma)}\leq\tfrac{4/\mu+8/\mu+5/(32\beta)}{32/\mu}\leq\tfrac{13}{32}.$
(72)

Next, owing to the assumption
$N_{k}\geq\tfrac{56\varsigma^{2}}{\mu(1-\gamma)^{2}}$, we have

$\displaystyle\left(\tfrac{4/\mu+(64T+160)\eta^{2}\varsigma^{2}/m}{T\eta(1-\gamma)}+2\right)\tfrac{2\varsigma^{2}}{\mu(1-\gamma)^{2}N_{k}}\leq\tfrac{3}{32}.$
(73)

Substituting Ineqs. (72) and (73) into Ineq. (7.5.1), we obtain

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-v^{*}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{2}{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]+\tfrac{5}{N_{k}}\text{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\},$

which completes the proof of Proposition 3.

###  7.6 Proof of Theorem 4

The structure of this proof is similar to proofs of Theorems 2 and 3: We first
derive a bound for a single epoch, and then apply it recursively to obtain the
eventual convergence result. The following proposition characterizes the
progress in a single epoch $k\in[K]$.

######  Proposition 4.

Consider a single epoch with index $k\in[K]$. Consider an integer $\tau$
satisfying
$\rho^{\tau}\leq\min\\{\frac{2(1-\rho)\varsigma}{3C_{M}},\frac{2(1-\rho)^{2}}{5C_{M}}\\}$.
Suppose the parameters $N_{k}$, $n_{0}$ and $m_{0}$ satisfy

$\displaystyle\rho^{N_{k}-n_{0}}\leq\tfrac{\tau(1-\rho)}{5C_{M}(N_{k}-n_{0})},~{}\rho^{n_{0}}\leq\tfrac{\min_{i\in[D]}\pi_{i}}{C_{P}},~{}\text{and}\quad\rho^{m_{0}}\leq\min\left\\{\tfrac{\min_{i\in[D]}\pi_{i}}{C_{P}},\tfrac{\sqrt{\mu}\eta\tau\varsigma^{2}(1-\rho)}{C_{M}}\right\\}.$
(74)

Suppose that the parameter $\eta$, $\lambda$, $m$, $N_{k}$ and $T$ satisfy

$\eta\leq\tfrac{1}{4\beta(1+\gamma)},~{}~{}\lambda=1,~{}~{}T\geq\tfrac{64}{\mu(1-\gamma)\eta},~{}~{}m-m_{0}\geq\max\\{1,\tfrac{792\eta(\tau+1)\varsigma^{2}}{1-\gamma}\\},~{}\text{and}~{}N_{k}-n_{0}\geq\tfrac{206(\tau+1)\varsigma^{2}}{\mu(1-\gamma)^{2}}.$
(75)

Set the output of this epoch to be
$\widehat{v}_{k}:=\frac{\sum_{t=2}^{T+1}v_{t}}{T}$. Then we have the following
results:
(a) If $v^{*}=\bar{v}$, we have

$\displaystyle{\mathbb{E}}\|\widehat{v}_{k}-v^{*}\|_{\Pi}^{2}\leq\tfrac{1}{2}{\mathbb{E}}\|\widehat{v}_{k-1}-\bar{v}\|_{\Pi}^{2}+\tfrac{10\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}}{N_{k}-N_{0}}.$
(76)

(b) If $v^{*}\neq\bar{v}$, we have

$\displaystyle{\mathbb{E}}\|\widehat{v}_{k}-v^{*}\|_{\Pi}^{2}\leq\tfrac{1}{2}{\mathbb{E}}\|\widehat{v}_{k-1}-\bar{v}\|_{\Pi}^{2}+\tfrac{10\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-\top}\\}}{N_{k}-N_{0}}+\tfrac{\widetilde{\mathcal{H}}}{(N_{k}-n_{0})^{2}},$
(77)

where
$\widetilde{\mathcal{H}}:=\tfrac{2(\tau+1)}{(1-\gamma)^{2}\mu}\|\bar{v}-v^{*}\|_{\Pi}^{2}+(10\tau^{2}+2\tau+2)\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\\}.$

Taking Proposition 4 as given for the moment, let us complete the proof of the
theorem. First, consider the case when $\bar{v}\neq v^{*}$. Recursively using
Ineq. (77) yields

$\displaystyle{\mathbb{E}}\|\widehat{v}_{K}-\bar{v}\|_{\Pi}^{2}$
$\displaystyle\leq\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|^{2}_{\Pi}+\textstyle{\sum}_{k=1}^{K}\left(\tfrac{10\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-\top}\\}}{2^{K-k}(N_{k}-n_{0})}+\tfrac{\widetilde{\mathcal{H}}}{2^{K-k}(N_{k}-n_{0})^{2}}\right)$
$\displaystyle\overset{(i)}{\leq}\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|^{2}_{\Pi}+\textstyle{\sum}_{k=1}^{K}(\tfrac{2}{3})^{K-k}\tfrac{10\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-\top}\\}}{N}+\textstyle{\sum}_{k=1}^{K}(\tfrac{8}{9})^{K-k}\tfrac{\widetilde{\mathcal{H}}}{N^{2}}$
$\displaystyle\leq\tfrac{1}{2^{K}}\|v^{0}-\bar{v}\|^{2}_{\Pi}+\tfrac{30\cdot\mathrm{trace}\\{(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-\top}\\}}{N}+\tfrac{9\widetilde{\mathcal{H}}}{N^{2}},$
(78)

where step (i) follows from the condition
$N_{k}-n_{0}\geq(\tfrac{3}{4})^{K-k}N$. The proof of the case when
$\bar{v}=v^{*}$ follows from the same derivation.

In order to complete the proof, it remains to prove Proposition 4, which we do
in the following subsection.

####  7.6.1 Proof of Proposition 4

The basic structure of the proof is similar to the proof of Proposition 3.
Recall that $\underline{\theta}$ satisfies
$g(\underline{\theta})-g(\widetilde{\theta})+\widehat{g}(\widetilde{\theta})=0$
and $\underline{v}=\Psi^{\top}\underline{\theta}$. The main idea of the proof
is to bound the terms $\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}$ and
$\|\underline{v}-\bar{v}\|_{\Pi}^{2}$ separately and use Young’s inequality to
bound $\|\widehat{v}_{k}-\bar{v}\|_{\Pi}^{2}$. Before we start, we require
some supporting lemmas that upper bound the variance and bias of the
stochastic operator $\widetilde{g}$ and averaged stochastic operators
$\widehat{g}$ and $\widetilde{g}_{t}$.

######  Lemma 12.

For every $t\in\mathbb{Z}_{+}$, $\theta,\theta^{\prime}\in\mathbb{R}^{d}$, and
$\tau\in\mathbb{Z}_{+}$ such that

$\displaystyle C_{P}\cdot\rho^{\tau}\leq\min_{i\in[D]}\pi_{i},$ (79)

we have

$\displaystyle{\mathbb{E}}[\|\widetilde{g}(\theta,\xi_{t+\tau})-\widetilde{g}(\theta^{\prime},\xi_{t+\tau})-g(\theta)+g(\theta^{\prime})\|^{2}_{2}|\mathcal{F}_{t}]\leq
2\varsigma^{2}\|v-v^{\prime}\|_{\Pi}^{2}$ (80)

with probability 1.

See Appendix B.5 for the proof of this lemma.

######  Lemma 13.

For every $\theta,\theta^{\prime}\in\mathbb{R}^{d}$,
$n_{0},\tau\in\mathbb{Z}_{+}$, if $n_{0}$ satisfies (79) and $\tau$ satisfies

$\displaystyle\rho^{\tau}\leq\tfrac{2(1-\rho)\varsigma}{3C_{M}},$ (81)

then we have with probability 1,

$\displaystyle{\mathbb{E}}\|\widehat{g}(\theta)-\widehat{g}(\theta^{\prime})-g(\theta)+g(\theta^{\prime})\|^{2}_{2}\leq\tfrac{4(\tau+1)\varsigma^{2}}{N_{k}-n_{0}}\|v-v^{\prime}\|_{\Pi}^{2}.$
(82)

See Appendix B.6 for the proof of this lemma.

Note that $\widetilde{g}_{t}(\cdot)$ has the same property (82) as
$\widehat{g}(\cdot)$. If $m_{0}$ satisfies Ineq. (79) and $\tau$ satisfies
Ineq. (81), then we have that for any $t\in[T]$,

$\displaystyle{\mathbb{E}}[\|\widetilde{g}_{t}(\theta)-\widetilde{g}_{t}(\theta^{\prime})-g(\theta)+g(\theta^{\prime})\|^{2}_{2}]\leq\tfrac{4(\tau+1)\varsigma^{2}}{m-m_{0}}\|v-v^{\prime}\|_{\Pi}^{2}.$
(83)

######  Lemma 14.

For every $\theta,\theta^{\prime}\in\mathbb{R}^{d}$, $m_{0}\in[m]$ and
$t\in[T]$, with probability 1,

$\displaystyle\|{\mathbb{E}}[\widetilde{g}_{t}(\theta)|\mathcal{F}_{t-1}]-{\mathbb{E}}[\widetilde{g}_{t}(\theta^{\prime})|\mathcal{F}_{t-1}]-g(\theta)+g(\theta^{\prime})\|_{2}\leq\tfrac{C_{M}\rho^{m_{0}}}{(1-\rho)(m-m_{0})}\|v-v^{\prime}\|_{\Pi}.$
(84)

See Appendix B.7 for the proof of this lemma.

With these supporting lemmas in hand, we are now ready to proceed to the meat
of the proof. We first focus on upper-bounding
$\|\bar{v}-\underline{v}\|_{\Pi}^{2}$ in the following lemma.

######  Lemma 15.

Consider a single epoch with index $k\in[K]$. Assume $n_{0}$ satisfies Ineq.
(79), $\tau$ satisfies Ineq. (81) and
$\rho^{\tau}\leq\frac{2(1-\rho)^{2}}{5C_{M}}$, and $N_{k}$ satisfies
$\rho^{N_{k}-n_{0}}\leq\tfrac{\tau(1-\rho)}{5C_{M}(N_{k}-n_{0})}$.
(a) If $\bar{v}=v^{*}$, then we have

$\displaystyle{\mathbb{E}}[\|\underline{v}-\bar{v}\|^{2}_{\Pi}]\leq$
$\displaystyle\tfrac{4\cdot\mathrm{trace}\big{\\{}(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\big{\\}}}{N_{k}-N_{0}}+\tfrac{8(\tau+1)\varsigma^{2}}{(1-\gamma)^{2}\mu(N_{k}-n_{0})}\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}.$
(85)

(b) If $\bar{v}\neq v^{*}$, then we have

$\displaystyle{\mathbb{E}}[\|\underline{v}-\bar{v}\|^{2}_{\Pi}]\leq$
$\displaystyle\tfrac{4\cdot\mathrm{trace}\big{\\{}(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{Mkv}}(I_{d}-M)^{-\top}\big{\\}}}{N_{k}-N_{0}}+\tfrac{4(\tau+1)}{5(1-\gamma)^{2}\mu(N_{k}-n_{0})^{2}}\|\bar{v}-v^{*}\|_{\Pi}^{2}$
$\displaystyle~{}~{}+\tfrac{(20\tau^{2}+4\tau+4)\cdot\mathrm{trace}\big{\\{}(I_{d}-M)^{-1}\bar{\Sigma}_{\mathsf{iid}}(I_{d}-M)^{-\top}\big{\\}}}{5(N_{k}-n_{0})^{2}}+\tfrac{8(\tau+1)\varsigma^{2}}{(1-\gamma)^{2}\mu(N_{k}-n_{0})}\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}.$
(86)

See Appendix B.8 for the proof of this lemma.

Next, we turn our attention to providing an upper bound on
$\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}$. At this juncture, it is
important to note that our notation is chosen to be consistent with the i.i.d.
setting. In particular, we keep the definition
$\widetilde{\Delta}_{t}(\theta):=\widetilde{g}_{t}(\theta)-g(\theta)$, and
note that $\widetilde{g}_{t}(\theta)$ here is defined in Algorithm 3 rather
than Algorithm 2. Let
$F(\theta):=g(\theta)-g(\widetilde{\theta})+\widehat{g}(\widetilde{\theta})$.

Now recall that in the proof of the corresponding proposition in the i.i.d.
setting (Proposition 3), our first few steps up to Ineq. (7.5.1) were obtained
purely via algebraic manipulation. Therefore, Ineq. (7.5.1) is still valid in
the Markovian noise setting, and rearranging terms in Ineq. (7.5.1) yields the
bound

$\displaystyle(1-8\eta^{2}(1+\gamma)^{2}\beta^{2})\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}+\textstyle{\sum}_{t=1}^{T}2\eta\langle
F(\theta_{t+1}),\theta_{t+1}-\underline{\theta}\rangle$
$\displaystyle\qquad\qquad\leq\|\theta_{1}-\underline{\theta}\|_{\Pi}^{2}+4\eta^{2}\|\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T})\|_{2}^{2}+Q_{3}+\textstyle{\sum}_{t=1}^{T}8\eta^{2}\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})\|_{2}^{2},$
(87)

where we have defined

$\displaystyle Q_{3}$
$\displaystyle:=2\eta\langle\widetilde{\Delta}_{T+1}(\theta_{T+1})+\widetilde{\Delta}_{T+1}(\widetilde{\theta}),\theta_{T+1}-\underline{\theta}\rangle-2\eta\langle\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T}),\theta_{T}-\underline{\theta}\rangle$
$\displaystyle\qquad\qquad-\textstyle{\sum}_{t=1}^{T}2\eta\langle\widetilde{F}_{t+1}(\theta_{t+1})-F(\theta_{t+1}),\theta_{t+1}-\underline{\theta}\rangle$
$\displaystyle~{}=2\eta\langle\widetilde{\Delta}_{T+1}(\theta_{T+1})-\widetilde{\Delta}_{T+1}(\widetilde{\theta}),\theta_{T+1}-\underline{\theta}\rangle+2\eta\langle\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T}),\theta_{T}-\underline{\theta}\rangle$
$\displaystyle\qquad\qquad-\textstyle{\sum}_{t=1}^{T}2\eta\langle\widetilde{\Delta}_{t+1}(\theta_{t+1})-\widetilde{\Delta}_{t+1}(\widetilde{\theta}),\theta_{t+1}-\underline{\theta}\rangle.$

Recall that $F(\underline{\theta})=0$. Invoking Lemma 5 yields

$\displaystyle\big{(}1-8\eta^{2}(1+\gamma)^{2}\beta^{2}\big{)}\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}+\textstyle{\sum}_{t=1}^{T}2\eta(1-\gamma)\|v_{t+1}-\underline{v}\|_{\Pi}^{2}$
$\displaystyle\qquad\qquad\qquad\leq\|\theta_{1}-\underline{\theta}\|_{\Pi}^{2}+4\eta^{2}\|\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T})\|_{2}^{2}+Q_{3}+\textstyle{\sum}_{t=1}^{T}8\eta^{2}\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})\|_{2}^{2}.$
(88)

As before, it remains to take an expectation on both sides of Ineq. (7.6.1).
Note that in the Markovian noise setting, the operator $\widetilde{g}_{t}$ is
not unbiased, and this necessitates a more careful analysis. We provide an
upper bound on ${\mathbb{E}}[Q_{3}]$ by writing

$\displaystyle{\mathbb{E}}[Q_{3}]$
$\displaystyle=2\eta{\mathbb{E}}\langle\widetilde{\Delta}_{T+1}(\theta_{T+1})-\widetilde{\Delta}_{T+1}(\widetilde{\theta}),\theta_{T+1}-\underline{\theta}\rangle+2\eta{\mathbb{E}}\langle\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T}),\theta_{T}-\underline{\theta}\rangle$
$\displaystyle\qquad\qquad-\textstyle{\sum}_{t=1}^{T}2\eta{\mathbb{E}}\langle\widetilde{\Delta}_{t+1}(\theta_{t+1})-\widetilde{\Delta}_{t+1}(\widetilde{\theta}),\theta_{t+1}-\underline{\theta}\rangle$
$\displaystyle\overset{(i)}{\leq}\tfrac{2\eta
C_{M}\rho^{m_{0}}\cdot{\mathbb{E}}[\|v_{T+1}-\widetilde{v}\|_{\Pi}\|\theta_{T+1}-\underline{\theta}\|_{2}]}{(1-\rho)(m-m_{0})}+\tfrac{2\eta
C_{M}\rho^{m_{0}}\cdot{\mathbb{E}}[\|v_{T}-\widetilde{v}\|_{\Pi}\|\theta_{T}-\underline{\theta}\|_{2}]}{(1-\rho)(m-m_{0})}$
$\displaystyle\qquad\qquad+\textstyle{\sum}_{t=1}^{T}\tfrac{2\eta
C_{M}\rho^{m_{0}}\cdot{\mathbb{E}}[\|v_{t+1}-\widetilde{v}\|_{\Pi}\|\theta_{t+1}-\underline{\theta}\|_{2}]}{(1-\rho)(m-m_{0})}$
$\displaystyle\overset{(ii)}{\leq}\tfrac{2\eta
C_{M}\rho^{m_{0}}\cdot{\mathbb{E}}[\|v_{T+1}-\widetilde{v}\|_{\Pi}\|v_{T+1}-\underline{v}\|_{\Pi}]}{\sqrt{\mu}(1-\rho)(m-m_{0})}+\tfrac{2\eta
C_{M}\rho^{m_{0}}}{\sqrt{\mu}(1-\rho)(m-m_{0})}\cdot{\mathbb{E}}[\|v_{T}-\widetilde{v}\|_{\Pi}\|v_{T}-\underline{v}\|_{\Pi}]$
$\displaystyle\qquad\qquad+\textstyle{\sum}_{t=1}^{T}\tfrac{2\eta
C_{M}\rho^{m_{0}}\cdot{\mathbb{E}}[\|v_{t+1}-\widetilde{v}\|_{\Pi}\|v_{t+1}-\underline{v}\|_{\Pi}]}{\sqrt{\mu}(1-\rho)(m-m_{0})},$
(89)

where step (i) follows from the Cauchy–Schwarz inequality and Lemma 14, and
step (ii) follows from the fact that
$\|\theta-\theta^{\prime}\|_{2}\leq\tfrac{1}{\sqrt{\mu}}\|v-v^{\prime}\|_{\Pi}$.
Next, note from Ineq. (83) that

$\displaystyle
4\eta^{2}{\mathbb{E}}\|\widetilde{\Delta}_{T}(\widetilde{\theta})-\widetilde{\Delta}_{T}(\theta_{T})\|_{2}^{2}$
$\displaystyle+\textstyle{\sum}_{t=1}^{T}8\eta^{2}{\mathbb{E}}\|\widetilde{\Delta}_{t}(\theta_{t})-\widetilde{\Delta}_{t}(\widetilde{\theta})\|_{2}^{2}$
$\displaystyle\leq\tfrac{16\eta^{2}(\tau+1)\varsigma^{2}}{m-m_{0}}{\mathbb{E}}[\|v_{T}-\widetilde{v}\|_{\Pi}^{2}]+\textstyle{\sum}_{t=1}^{T}\tfrac{32\eta^{2}(\tau+1)\varsigma^{2}}{m-m_{0}}{\mathbb{E}}[\|v_{t}-\widetilde{v}\|_{\Pi}^{2}].$

Putting together the pieces yields

$\displaystyle\big{(}1-8\eta^{2}$
$\displaystyle(1+\gamma)^{2}\beta^{2}\big{)}{\mathbb{E}}[\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}]+\textstyle{\sum}_{t=2}^{T+1}2\eta(1-\gamma){\mathbb{E}}[\|v_{t}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1}{\mu}{\mathbb{E}}[\|v_{1}-\underline{v}\|_{\Pi}^{2}]+\tfrac{2\eta
C_{M}\rho^{m_{0}}}{\sqrt{\mu}(1-\rho)(m-m_{0})}{\mathbb{E}}[\|v_{T}-\widetilde{v}\|_{\Pi}\|v_{T}-\underline{v}\|_{\Pi}]$
$\displaystyle\quad\quad+\tfrac{2\eta
C_{M}\rho^{m_{0}}}{\sqrt{\mu}(1-\rho)(m-m_{0})}{\mathbb{E}}[\|v_{T+1}-\widetilde{v}\|_{\Pi}\|v_{T+1}-\underline{v}\|_{\Pi}]+\tfrac{16\eta^{2}(\tau+1)\varsigma^{2}}{m-m_{0}}{\mathbb{E}}[\|v_{T}-\widetilde{v}\|_{\Pi}^{2}]$
$\displaystyle\quad\quad+\textstyle{\sum}_{t=1}^{T}\tfrac{32\eta^{2}(\tau+1)\varsigma^{2}}{m-m_{0}}{\mathbb{E}}[\|v_{t}-\widetilde{v}\|_{\Pi}^{2}]+\textstyle{\sum}_{t=2}^{T+1}\tfrac{2\eta
C_{M}\rho^{m_{0}}}{\sqrt{\mu}(1-\rho)(m-m_{0})}{\mathbb{E}}[\|v_{t}-\widetilde{v}\|_{\Pi}\|v_{t}-\underline{v}\|_{\Pi}].$

We now use the relations

$\displaystyle\|v_{t}-\widetilde{v}\|_{\Pi}^{2}\leq
2\|v_{t}-\underline{v}\|_{\Pi}^{2}+2\|\widetilde{v}-\underline{v}\|_{\Pi}^{2},$
$\displaystyle\quad\|v_{t}-\widetilde{v}\|_{\Pi}\|v_{t}-\underline{v}\|_{\Pi}\leq\tfrac{1}{2}(\|v_{t}-\widetilde{v}\|_{\Pi}^{2}+\|v_{t}-\underline{v}\|_{\Pi}^{2}),\text{
and }$ $\displaystyle v_{1}=v_{0}=\widetilde{v}.$

Rearranging terms above, we then obtain

$\displaystyle\textstyle{\sum}_{t=2}^{T+1}\big{(}2\eta(1-\gamma)-\tfrac{96\eta^{2}(\tau+1)\varsigma^{2}}{m-m_{0}}-\tfrac{3\eta
C_{M}\rho^{m_{0}}}{\sqrt{\mu}(1-\rho)(m-m_{0})}\big{)}{\mathbb{E}}\|v_{t}-\underline{v}\|_{\Pi}^{2}$
$\displaystyle\leq\big{(}\tfrac{1}{\mu}+\tfrac{(64T+32)\eta^{2}(\tau+1)\varsigma^{2}}{m-m_{0}}+\tfrac{2(T+2)\eta
C_{M}\rho^{m_{0}}}{\sqrt{\mu}(1-\rho)(m-m_{0})}\big{)}{\mathbb{E}}\|\widetilde{v}-\underline{v}\|_{\Pi}^{2}-\big{(}1-8\eta^{2}(1+\gamma)^{2}\beta^{2}\big{)}{\mathbb{E}}\|\theta_{T+1}-\underline{\theta}\|_{2}^{2}.$

Recall that $m-m_{0}\geq\tfrac{192\eta(\tau+1)\varsigma^{2}}{1-\gamma}$,
$1\geq 8\eta^{2}(1+\gamma)^{2}\beta^{2}$,
$\rho^{m_{0}}\leq\tfrac{\sqrt{\mu}\eta\tau\varsigma^{2}(1-\rho)}{C_{M}}$ and
$\widehat{v}_{k}=\tfrac{\sum_{t=2}^{T+1}v_{t}}{T}$. Combining these with
Jensen’s inequality yields

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{\sum_{t=2}^{T+1}\eta(1-\gamma){\mathbb{E}}[\|v_{t}-\underline{v}\|_{\Pi}^{2}]}{T\eta(1-\gamma)}\leq\tfrac{1/\mu+(66T+36)\eta^{2}(\tau+1)\varsigma^{2}/(m-m_{0})}{T\eta(1-\gamma)}{\mathbb{E}}[\|\widetilde{v}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq\tfrac{1/\mu+(66T+36)\eta^{2}(\tau+1)\varsigma^{2}/(m-m_{0})}{T\eta(1-\gamma)}(2{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]+2{\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}]).$

By Young’s inequality, the following chain of bounds holds:

$\displaystyle{\mathbb{E}}[\|\widehat{v}_{k}-\bar{v}\|_{\Pi}^{2}]$
$\displaystyle\leq
2{\mathbb{E}}[\|\widehat{v}_{k}-\underline{v}\|_{\Pi}^{2}]+2{\mathbb{E}}[\|\bar{v}-\underline{v}\|_{\Pi}^{2}]$
$\displaystyle\leq
D_{1}\cdot{\mathbb{E}}[\|\widetilde{v}-\bar{v}\|_{\Pi}^{2}]+\left(D_{1}+2\right){\mathbb{E}}[\|\underline{v}-\bar{v}\|_{\Pi}^{2}].$
(90)

where
$D_{1}:=\tfrac{4/\mu+(264T+144)\eta^{2}(\tau+1)\varsigma^{2}/(m-m_{0})}{T\eta(1-\gamma)}$.
By assumption, we have $\eta\leq\tfrac{1}{4\beta(1+\gamma)}$,
$T\geq\tfrac{64}{\mu(1-\gamma)\eta}$, and
$m-m_{0}\geq\max\\{1,\tfrac{792\eta(\tau+1)\varsigma^{2}}{1-\gamma}\\}$. Thus,

$\displaystyle
D_{1}\leq\tfrac{4/\mu+T\eta(1-\gamma)/3+2\eta(1-\gamma)/11}{T\eta(1-\gamma)}\leq\tfrac{4/\mu+64/3\mu+1/22\beta}{64/\mu}\leq\tfrac{26}{64}.$
(91)

Next, note that
$N_{k}-n_{0}\geq\tfrac{206(\tau+1)\varsigma^{2}}{\mu(1-\gamma)^{2}}$ to write

$\displaystyle\left(D_{1}+2\right)\cdot\tfrac{8(\tau+1)\varsigma^{2}}{(1-\gamma)^{2}\mu(N_{k}-n_{0})}\leq\tfrac{154}{64}\cdot\tfrac{8}{206}<\tfrac{6}{64}.$
(92)

Substituting Ineqs. (91) and (92) into Ineq. (7.6.1), and utilizing Lemma 15
to bound ${\mathbb{E}}[\|\bar{v}-\underline{v}\|_{\Pi}^{2}]$, we finally
obtain the desired results for both cases of $v^{*}=\bar{v}$ and
$v^{*}\neq\bar{v}$ as desired.

##  8 Discussion

In this paper, we investigated the problem of policy evaluation with linear
function approximation, making three contributions. First, we proved lower
bounds on both deterministic error and stochastic error. With these lower
bounds in hand, we presented an analysis of a variance-reduced variant of
temporal difference algorithm (VRTD) in the i.i.d. observation model and
showed that it fails to match the oracle complexity lower bound on the
deterministic error. In order to remedy this difficulty, we developed an
optimal variance-reduced fast temporal difference algorithm (VRFTD) that
nearly matches both lower bounds simultaneously. Finally, we extended the
VRFTD algorithm to the Markovian setting and provided instance-dependent
convergence results. The leading stochastic error matches the instance-
dependent lower bound for Markovian linear stochastic approximation (Mou et
al., 2021), and the deterministic error matches the i.i.d. setting up to a
multiplicative factor proportional to the mixing time of the chain. Our
theoretical guarantees were corroborated with numerical experiments in both
the i.i.d. and Markovian settings, showing that the VRFTD algorithm enjoys
several advantages over the prior state-of-the-art.

Our work leaves open severaal salient future directions; let us mention two.
First, our oracle complexity lower bound is proved in the tabular setting. On
the other hand, our upper bounds on the deterministic error indicate that with
linear function approximation, we pay a multiplicative factor depending on the
condition number of the feature matrix. It would be interesting to see if an
oracle complexity lower bound can be proved under linear function
approximation, and whether the linear dependence on the condition number in
our bounds is optimal. Second, and more broadly, note that our analysis relies
heavily on the linear structure of the problem. However, there are many
problems in the reinforcement learning literature that have nonlinear
structures, e.g., the policy optimization problem involving the Bellman
optimality operator. An interesting direction for future work is to understand
problems with nonlinear structure from an instance-specific point of view and
develop efficient algorithms to capture the optimal deterministic and
stochastic errors. For instance, variance reduction has been applied to the
policy optimization problem under the generative model (Sidford et al., 2018;
Wainwright, 2019) and some instance-dependent bounds are known (e.g., for
variants of $Q$-learning Khamaru et al., 2021b). It is an important open
question to develop acceleration schemes for such algorithms in a fashion
similar to our paper, while extending the results to the more realistic
Markovian setting.

### Acknowledgments

TL and GL were supported in part by Office of Naval Research grant
N00014-20-1-2089. TL and AP were supported in part by the National Science
Foundation grant CCF-2107455, and are thankful to the Simons Institute for the
Theory of Computing for their hospitality when part of this work was
performed.

## References

  * Bhandari et al. [2018] J. Bhandari, D. Russo, and R. Singal.  A finite time analysis of temporal difference learning with linear function approximation.  In _Conference on learning theory_ , pages 1691–1692. PMLR, 2018\.
  * Blatt et al. [2007] D. Blatt, A. O. Hero, and H. Gauchman.  A convergent incremental gradient method with a constant step size.  _SIAM Journal on Optimization_ , 18(1):29–51, 2007.
  * Borkar [2009] V. S. Borkar.  _Stochastic approximation: A dynamical systems viewpoint_ , volume 48.  Springer, 2009.
  * Borkar and Meyn [2000] V. S. Borkar and S. P. Meyn.  The ODE method for convergence of stochastic approximation and reinforcement learning.  _SIAM Journal on Control and Optimization_ , 38(2):447–469, 2000.
  * Chen et al. [2021] Z. Chen, S. T. Maguluri, S. Shakkottai, and K. Shanmugam.  A Lyapunov theory for finite-sample guarantees of asynchronous Q-learning and TD-learning variants.  _arXiv preprint arXiv:2102.01567_ , 2021.
  * Dann et al. [2014] C. Dann, G. Neumann, and J. Peters.  Policy evaluation with temporal differences: A survey and comparison.  _Journal of Machine Learning Research_ , 15:809–883, 2014\.
  * Defazio et al. [2014] A. Defazio, F. Bach, and S. Lacoste-Julien.  SAGA: A fast incremental gradient method with support for non-strongly convex composite objectives.  In _Advances in neural information processing systems_ , pages 1646–1654, 2014.
  * Du et al. [2017] S. S. Du, J. Chen, L. Li, L. Xiao, and D. Zhou.  Stochastic variance reduction methods for policy evaluation.  In _International Conference on Machine Learning_ , pages 1049–1058. PMLR, 2017.
  * Duan et al. [2021] Y. Duan, M. Wang, and M. J. Wainwright.  Optimal policy evaluation using kernel-based temporal difference methods.  _arXiv preprint arXiv:2109.12002_ , 2021.
  * Durmus et al. [2021] A. Durmus, E. Moulines, A. Naumov, S. Samsonov, and H. T. Wai.  On the stability of random matrix product with Markovian noise: Application to linear stochastic approximation and TD learning.  _arXiv preprint arXiv:2102.00185_ , 2021.
  * Hájek [1972] J. Hájek.  Local asymptotic minimax and admissibility in estimation.  In _Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 1: Theory of Statistics_ , pages 175–194. University of California Press, 1972.