# Robust and Adaptive Temporal-Difference Learning Using An Ensemble of
Gaussian Processes

Qin Lu1 Georgios B. Giannakis1
1University of Minnesota
{qlu<EMAIL_ADDRESS>

###### Abstract

Value function approximation is a crucial module for policy evaluation in
reinforcement learning when the state space is large or continuous. The
present paper takes a generative perspective on policy evaluation via
temporal-difference (TD) learning, where a Gaussian process (GP) prior is
presumed on the sought value function, and instantaneous rewards are
probabilistically generated based on value function evaluations at two
consecutive states. Capitalizing on a random feature-based approximant of the
GP prior, an online scalable (OS) approach, termed OS-GPTD, is developed to
estimate the value function for a given policy by observing a sequence of
state-reward pairs. To benchmark the performance of OS-GPTD even in an
adversarial setting, where the modeling assumptions are violated,
complementary worst-case analyses are performed by upper-bounding the
cumulative Bellman error as well as the long-term reward prediction error,
relative to their counterparts from a fixed value function estimator with the
entire state-reward trajectory in hindsight. Moreover, to alleviate the
limited expressiveness associated with a single fixed kernel, a weighted
ensemble (E) of GP priors is employed to yield an alternative scheme, termed
OS-EGPTD, that can jointly infer the value function, and select interactively
the EGP kernel on-the-fly. Finally, performances of the novel OS-(E)GPTD
schemes are evaluated on two benchmark problems.

##  1 Introduction

In reinforcement learning (RL), an agent learns a policy function yielding as
output a sequence of actions with input the observed states of an environment,
so that a cumulative reward criterion is maximized. As a dynamic system, the
environment can be _unknown, nonlinear, and complex_ , while the agent has to
learn how to act directly from state-reward observations. RL has documented
success in several disciplines, including e.g., control, information theory,
operations research, machine learning, and statistics, to list a few [28].

Many algorithms for obtaining an optimal policy achieving the highest
cumulative reward from each state are based on the policy iteration method
[20], in which the so-termed value function must be estimated for a sequence
of fixed policies, making _value function estimation_ (policy evaluation) a
critical algorithmic component. Temporal-difference (TD) learning [28], is
arguably the best known method for value function estimation. It relies on
forming the so-called TDs to update an estimate of the value function on-the-
fly. Although intuitively simple, the classic TD learning algorithm can be
intractable when the state space of the environment is large or continuous,
often the case in contemporary control and artificial intelligence tasks [19].
A scalable practice to overcome this hurdle leverages approximants of the
value function, including e.g., linear ones [28], deep neural networks [19],
or Gaussian processes (GPs) [7, 12].

Speaking of GPs, cross-fertilizing the merits of kernel learning and Bayesian
methods, they have been widely employed in several learning tasks, including
regression, classification, clustering, and dimensionality reduction, in
addition to RL [23]. Postulating a probabilistic generative model, GP-based
methods not only allow learning over a rich class of nonlinear functions, but
also quantify the associated uncertainty. Leveraging GPs for RL has led to
efforts in model-based approaches including [12, 5, 6, 11], as well as model-
free ones, as in [7, 8, 9, 24]. Specifically in [7], the developed GPTD
algorithm estimates the value function from a Bayesian perspective, where the
value function is assumed drawn from a GP prior, and a probabilistic
observation model is postulated to generate the instantaneous reward from
value function realizations at two consecutive states. Further, a revised
observation model has been proposed in [8], along with a GP-based SARSA
algorithm.

Concerning _performance analysis_ of TD methods, several attempts have been
made. Under some probabilistic assumptions about the environment to effect a
sequence of Markovian observations, or of samples from a static Markov
decision process (MDP), both asymptotic and non-asymptotic analyses of TD
learning with function approximation have been widely studied; see e.g., [29,
3, 4, 2, 26]. Without any statistical assumptions about the environment,
worst-case analysis on the long-term reward predictive error has been reported
for TD learning and its variants in an online setting [25, 14]. Further along
this line, connections between online long-term reward prediction error and
Bellman error were drawn in [27]. However, their results do not directly apply
to GPTD-type algorithms. Among other GP-based RL frameworks, sample complexity
analysis has been conducted for both model-based and model-free methods in
[9].

Contributions. Building on the GPTD framework for value estimation, the
present work develops _online scalable_ (OS) schemes whose performance can be
also quantified using worst-case analyses. Relative to past works, the
contributions can be summarized as follows.111Part of this work is published
in[15].

  1. c1)

Leveraging the random feature (RF) approximant of the GP prior, we develop an
OS-GPTD scheme, which readily accommodates streaming reward-state observations
to learn the value function on-the-fly.

  2. c2)

To benchmark the performance of OS-GPTD even in the adversarial setting where
state-reward pairs can be chosen adversarially to violate the generative
modeling assumptions, two complementary worst-case analyses are provided. The
first analysis establishes the _no-regret_ property of the online function
estimator given by OS-GPTD in terms of Bellman error, while the second one
further bounds the accumulated future reward prediction error incurred by OS-
GPTD compared with that corresponding to the best fixed function estimator
with data in hindsight.

  3. c3)

To enrich the expressiveness of the sought value function, a weighted ensemble
(E-) of OS-GPTDs from multiple GP experts with different priors, are employed
to not only infer the value function, but also select the appropriate
kernel(s) adaptively.

##  2 Preliminaries

To elicit the main contributions of this paper, some preliminaries on RL and
GPTD learning in particular, are outlined in this section.

Characterizing the interactions between an agent and the environment, an RL
problem is generally formulated as a discounted MDP specified by the quintuple
$\mathcal{M}:=\\{\mathcal{S},\mathcal{A},\mathcal{P},{R},\gamma\\}$, where
$\mathcal{S}$ is the state space of the environment, $\mathcal{A}$ is the
action space of the agent, $\mathcal{P}$ defines the state dynamics, ${R}$ is
the reward function, and $\gamma$ is the discount factor belonging to $[0,1)$.

Capitalizing on simulated or real-world experiences, the objective of the
agent is to obtain an optimal policy that maximizes the cumulative rewards
from each state. To this end, policy- and value-based approaches have been
devised; see, e.g., [20]. The state value function is defined as

$v(\mathbf{s}):=\mathbb{E}\\!\left[\sum_{\tau=1}^{\infty}\gamma^{\tau-1}r_{\tau}|\mathbf{s}_{1}=\mathbf{s}\right]$
(1)

where the expectation is taken over the entire MDP trajectory conditioned on
the initial state $\mathbf{s}_{1}=\mathbf{s}\in\mathcal{S}$, while following
some fixed policy to take future actions. Value-based methods proceed by
alternating between value function estimation for a given policy, and policy
improvement using the predicted value function. Value function estimation
relies on the popular TD-based iteration indexed by superscript $i$ with
scalar $\alpha>0$ [28]

$v^{(i+1)}(\mathbf{s}_{t})=v^{(i)}(\mathbf{s}_{t})+\alpha\\!\left[r_{t}+\gamma
v^{(i)}(\mathbf{s}_{t+1})-v^{(i)}(\mathbf{s}_{t})\right].$ (2)

When the state space becomes large or continuous, the computational
requirements of value function estimation in (2) are overwhelming as all
solvers of Bellman’s equation. This has motivated scalable approaches
leveraging low-dimensional approximants of $v$, including linear [28], deep
neural networks [19], and Bayesian ones [7].

In Bayesian approximators, $v$ is assumed drawn from a GP prior denoted as
$v\sim\mathcal{GP}(0,\kappa(\mathbf{s},\mathbf{s}^{\prime}))$, where $\kappa$
is a kernel function measuring the pairwise similarity between the function
values at $d$-dimensional states $\mathbf{s}$ and $\mathbf{s}^{\prime}$. The
GP prior on $v$ implies that the joint probability density function (pdf) for
any number of function values
$\mathbf{v}_{t}:=[v(\mathbf{s}_{1})~{}v(\mathbf{s}_{2})~{}\cdots~{}v(\mathbf{s}_{t})]^{\top}$
(⊤ stands for transposition) evaluated at
$\mathbf{S}_{t}:=[\mathbf{s}_{1}~{}\mathbf{s}_{2}~{}\cdots~{}\mathbf{s}_{t}]$
is multivariate Gaussian [23]

$p(\mathbf{v}_{t}|\mathbf{S}_{t})=\mathcal{N}(\mathbf{v}_{t};\mathbf{0}_{t},\mathbf{K}_{t})$
(3)

where $\mathbf{0}_{t}$ is the $t$-dimensional all-zero vector, and the
$(i,j)$-th entry of $\mathbf{K}_{t}\in\mathbb{R}^{t\times t}$ denotes the
covariance between $v(\mathbf{s}_{i})$ and $v(\mathbf{s}_{j})$; that is,
$[\mathbf{K}_{t}]_{i,j}:={\rm
cov}(v(\mathbf{s}_{i}),v(\mathbf{s}_{j}))=\kappa(\mathbf{s}_{i},\mathbf{s}_{j})$.

Upon convergence, the update in (2) implies that the term in square brackets
must vanish. This in turn suggests a generative model relating the reward
$r_{t}$ with the value function $v(\mathbf{s})$ as [7]

$r_{t}=v(\mathbf{s}_{t})-\gamma v(\mathbf{s}_{t+1})+n_{t}$ (4)

where the noise $n_{t}$ is assumed white, uncorrelated with any
$v(\mathbf{s})$, and zero-mean Gaussian with variance $\sigma_{n}^{2}$.

Starting from some initial state $\mathbf{s}_{1}\in\mathcal{S}$, taking
subsequent actions by following a given policy generates a trajectory of
states and rewards, namely $\\{\mathbf{s}_{1},r_{1},\mathbf{s}_{2},r_{2},$
$\mathbf{s}_{3},\ldots\\}$. Upon collecting the rewards into vector
$\mathbf{r}_{t}:=[r_{1}~{}r_{2}~{}\cdots~{}r_{t}]^{\top}$, the conditional
likelihood is

$p(\mathbf{r}_{t}|\mathbf{v}_{t+1},\mathbf{S}_{t+1})=\mathcal{N}(\mathbf{r}_{t};\mathbf{H}_{t}\mathbf{v}_{t+1},\sigma^{2}_{n}\mathbf{I}_{t})$
(5)

where $\mathbf{I}_{t}$ is the $t\times t$ identity matrix and the
$t\times(t+1)$ observation matrix is

$\mathbf{H}_{t}:=\begin{bmatrix}1&-\gamma&0&\cdots&0\\\
0&1&-\gamma&\cdots&0\\\ \vdots&&\ddots&\ddots&\vdots\\\
0&0&\cdots&1&-\gamma\end{bmatrix}.$

Given the GP prior in (3) and the conditional likelihood in (5) for
$\mathbf{r}_{t}=\mathbf{H}_{t}\mathbf{v}_{t+1}+\mathbf{n}_{t}$, the goal is to
predict the pdf of $v(\cdot)$ at some unseen state $\mathbf{s}$. With this
goal in mind, let us start off by writing the joint pdf of $\mathbf{r}_{t}$
and $v(\mathbf{s})$ as

$\displaystyle\begin{bmatrix}\mathbf{r}_{t}\\\
v(\mathbf{s})\end{bmatrix}\sim\mathcal{N}\\!\left(\\!\mathbf{0}_{t+1},\begin{bmatrix}\mathbf{Q}_{t}&\mathbf{H}_{t}\mathbf{k}_{t+1}(\mathbf{s})\\\
\mathbf{k}_{t+1}^{\top}(\mathbf{s})\mathbf{H}_{t}^{\top}&\kappa(\mathbf{s},\mathbf{s})\end{bmatrix}\\!\right)$
(6)

where
$\mathbf{Q}_{t}=\mathbf{H}_{t}\mathbf{K}_{t+1}\mathbf{H}_{t}^{\top}+\sigma^{2}_{n}\mathbf{I}_{t}$
is the covariance matrix of the marginal pdf

$\displaystyle p(\mathbf{r}_{t}|\mathbf{S}_{t+1})$ $\displaystyle=\int
p(\mathbf{r}_{t}|\mathbf{v}_{t+1})p(\mathbf{v}_{t+1}|\mathbf{S}_{t+1})d\mathbf{v}_{t+1}$
$\displaystyle=\mathcal{N}(\mathbf{r}_{t};\mathbf{0}_{t},\mathbf{Q}_{t})$

and $\mathbf{H}_{t}\mathbf{k}_{t+1}(\mathbf{s})$ denotes the cross-covariance

$\mathbb{E}[\mathbf{r}_{t}v(\mathbf{s})|\mathbf{S}_{t+1}]=\mathbb{E}[\mathbf{H}_{t}\mathbf{v}_{t+1}v(\mathbf{s})|\mathbf{S}_{t+1}]=\mathbf{H}_{t}\mathbf{k}_{t+1}(\mathbf{s})$

with
$[\mathbf{k}_{t+1}(\mathbf{s})]_{\tau}:=\kappa(\mathbf{s}_{\tau},\mathbf{s})$
for $\tau=1,2,\ldots,t+1$.

Leveraging the pdf (6), the conditional pdf of $v(\mathbf{s})$ is

$p(v(\mathbf{s})|\mathbf{r}_{t},\mathbf{S}_{t+1})=\mathcal{N}(v(\mathbf{s});\hat{v}(\mathbf{s}),\sigma_{v}^{2}(\mathbf{s}))$
(7)

with mean and variance given by

$\displaystyle\hat{v}(\mathbf{s})$
$\displaystyle=\mathbf{k}_{t+1}^{\top}(\mathbf{s})\mathbf{H}_{t}^{\top}\mathbf{Q}_{t}^{-1}\mathbf{r}_{t}$
(8a) $\displaystyle\sigma_{v}^{2}(\mathbf{s})$
$\displaystyle=\kappa(\mathbf{s},\mathbf{s})\\!-\\!\mathbf{k}_{t+1}^{\top}(\mathbf{s})\mathbf{H}_{t}^{\top}\mathbf{Q}_{t}^{-1}\mathbf{H}_{t}\mathbf{k}_{t+1}(\mathbf{s}).$
(8b)

At negligible extra cost, the accompanying variance (8b) not only provides a
self-assessment of the accuracy of estimate (8a), but can also guide
exploration for policy improvement. Clearly, evaluating the function (7)
incurs computational complexity $\mathcal{O}(t^{3})$ and storage complexity
$\mathcal{O}(t^{2})$, which can become prohibitive as $t$ grows. Although a
scalable solution has been devised for the online GPTD algorithm in [7], the
estimation performance is compromised by the information loss incurred by the
sparsification employed.

Our pursuit in the rest of this paper aims at online scalable solvers and
their associated performance analyses.

##  3 Scalable TD learning with a single GP

Several attempts have been made to effect scalability in GP-based learning;
see, e.g., [21, 13]. For GPTD in particular, a sparsification scheme was
advocated in [7]; and a scalable value function estimator was reported in [18]
using an approximant [21], which summarizes the training data through a
reduced number of pseudo data with inducing inputs. Unfortunately, performance
of both low-rank approximants of $\mathbf{K}_{t}$ has not been analyzed.

Aiming at a low-rank approximant with quantifiable performance, we rely here
on a shift-invariant, real, and even $\bar{\kappa}(\cdot)$, which is given by
the inverse Fourier transform of a power spectral density
$\pi_{\bar{\kappa}}$; that is,

$\bar{\kappa}(\mathbf{s},\mathbf{s}^{\prime})=\bar{\kappa}(\mathbf{s}-\mathbf{s}^{\prime})=\int\pi_{\bar{\kappa}}({\mbox{\boldmath$\zeta$}})e^{j{\mbox{\boldmath$\zeta$}}^{\top}(\mathbf{s}-\mathbf{s}^{\prime})}d{\mbox{\boldmath$\zeta$}}.$
(9)

If $\bar{\kappa}(\cdot)$ is standardized so that
$\bar{\kappa}(\mathbf{0})=\int\pi_{\bar{\kappa}}({\mbox{\boldmath$\zeta$}})d{\mbox{\boldmath$\zeta$}}=1$,
then $\pi_{\bar{\kappa}}$ boils down to a pdf. Thus, the integral in (9) can
be expressed as the expectation

$\bar{\kappa}(\mathbf{s},\mathbf{s}^{\prime})\\!=\\!\mathbb{E}_{\pi_{\bar{\kappa}}}\\!\left[e^{j{\mbox{\boldmath$\zeta$}}^{\top}(\mathbf{s}-\mathbf{s}^{\prime})}\right]=\mathbb{E}_{\pi_{\bar{\kappa}}}\\!\left[\cos\\!\left({\mbox{\boldmath$\zeta$}}^{\top}(\mathbf{s}-\mathbf{s}^{\prime})\right)\right].$

Upon drawing a sufficiently large number $D$ ($\gg d$), of independent and
identically distributed (i.i.d.) samples (a.k.a. features)
$\\{{\mbox{\boldmath$\zeta$}}_{i}\\}_{i=1}^{D}$ from $\pi_{\bar{\kappa}}$, the
last expectation over $\zeta$ can be approximated by the sample average,
222Quantities with $\check{(\cdot)}$ involve the RF approximation. which can
be further rewritten as

$\check{\bar{\kappa}}(\mathbf{s},\mathbf{s}^{\prime})={\mbox{\boldmath$\phi$}}_{\zeta}^{\top}(\mathbf{s}){\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}^{\prime})$
(10)

where ${\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s})$ is the real $2D\times 1$
random feature (RF) vector given by [22]

$\displaystyle{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}):=$ (11)
$\displaystyle\
\frac{1}{\sqrt{D}}\\!\left[\sin({\mbox{\boldmath$\zeta$}}_{1}^{\top}\mathbf{s})~{}\cos({\mbox{\boldmath$\zeta$}}_{1}^{\top}\mathbf{s})~{}\cdots~{}\sin({\mbox{\boldmath$\zeta$}}_{D}^{\top}\mathbf{s})~{}\cos({\mbox{\boldmath$\zeta$}}_{D}^{\top}\mathbf{s})\right]^{\top}\\!\\!.$

Consider now the RF-based linear function approximant

${\check{v}}(\mathbf{s})={\mbox{\boldmath$\phi$}}_{\zeta}^{\top}(\mathbf{s}){\mbox{\boldmath$\theta$}},~{}{\rm
with}~{}p({\mbox{\boldmath$\theta$}})=\mathcal{N}({\mbox{\boldmath$\theta$}};\mathbf{0}_{2D},\sigma_{\theta}^{2}\mathbf{I}_{2D})$
(12)

where $\sigma_{\theta}^{2}:=\kappa(0)=\sigma_{\theta}^{2}\bar{\kappa}(0)$. The
prior pdf of
$\check{\mathbf{v}}_{t}:=[\check{v}(\mathbf{s}_{1})~{}\check{v}(\mathbf{s}_{2})~{}\cdots$
$\check{v}(\mathbf{s}_{t})]^{\top}$ then becomes

$p(\check{\mathbf{v}}_{t}|\mathbf{S}_{t})=\mathcal{N}(\check{\mathbf{v}}_{t};\mathbf{0}_{t},\check{\mathbf{K}}_{t}),~{}{\rm
with}~{}\,\check{\mathbf{K}}_{t}=\sigma_{\theta}^{2}{\mbox{\boldmath$\Phi$}}_{t}{\mbox{\boldmath$\Phi$}}_{t}^{\top}$
(13)

where
${\mbox{\boldmath$\Phi$}}_{t}:=[{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{1})~{}{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{2})~{}\cdots~{}{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t})]^{\top}$.
Evidently,
$\check{\mathbf{K}}_{t}=\sigma_{\theta}^{2}{\mbox{\boldmath$\Phi$}}_{t}{\mbox{\boldmath$\Phi$}}_{t}^{\top}$
is a low-rank approximant of $\mathbf{K}_{t}$ in (3) when $t>2D$. With this
approximate GP prior, the parametric function ${\check{v}}$ enjoys the
nonparametric GP model expressiveness, with the simplicity of linear models.

Specifically for GPTD, this RF-based GP prior (13) allows us to approximate
the predictor in (7) through

$\displaystyle\hat{\check{v}}(\mathbf{s})\\!=\\!{\mbox{\boldmath$\phi$}}_{\zeta}^{\top}\\!(\mathbf{s})\bigg{(}\\!{\mbox{\boldmath$\Phi$}}_{t+1}^{\top}\mathbf{H}_{t}^{\top}\mathbf{H}_{t}{\mbox{\boldmath$\Phi$}}_{t+1}\\!+\\!\frac{\sigma_{n}^{2}}{\sigma_{\theta}^{2}}\mathbf{I}_{2D}\bigg{)}^{-1}\\!\\!\\!\\!\\!\\!{\mbox{\boldmath$\Phi$}}_{t+1}^{\top}\mathbf{H}_{t}^{\top}\mathbf{r}_{t}\vspace{-.4cm}$
(14a)
$\displaystyle\sigma_{\check{v}}^{2}(\mathbf{s})\\!=\\!{\mbox{\boldmath$\phi$}}_{\zeta}^{\top}\\!(\mathbf{s})\bigg{(}\\!\frac{{\mbox{\boldmath$\Phi$}}_{t+1}^{\top}\mathbf{H}_{t}^{\top}\\!\mathbf{H}_{t}{\mbox{\boldmath$\Phi$}}_{t+1}}{\sigma_{n}^{2}}\\!+\\!\frac{\mathbf{I}_{2D}}{\sigma_{\theta}^{2}}\bigg{)}^{-1}\\!\\!\\!\\!\\!\\!\\!{\mbox{\boldmath$\phi$}}_{\zeta}\\!(\mathbf{s})$
(14b)

which incurs computational complexity $\mathcal{O}(2Dt(t+1)+(2D)^{3})$, that
is certainly more affordable than the plain-vanilla GP predictor in (8a) and
(8b), for $t\geq 2D$.

Besides being scalable, the parametric function approximant $\check{v}$
readily accommodates streaming data $\mathbf{s}_{t+1}$ and $r_{t}$ in a
recursive Bayesian fashion, leading to our online scalable (OS-) GPTD
approach. Upon replacing $v$ by $\check{v}$ in the observation model (4), the
per datum likelihood is

$p(r_{t}|{\mbox{\boldmath$\theta$}},\mathbf{S}_{t+1})=\mathcal{N}(r_{t};\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\theta$}},\sigma_{n}^{2})$
(15)

where
$\mathbf{h}_{t}:={\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t})-\gamma{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t+1})$.
This linear Gaussian likelihood of (15) then propagates the Gaussianity from
the prior (cf. (12)) to the posterior across time.

Capitalizing on the past rewards $\mathbf{r}_{t-1}$ and states
$\mathbf{S}_{t}$ up until slot $t-1$, the agent keeps track of the posterior
pdf of $\theta$ as
$p({\mbox{\boldmath$\theta$}}|\mathbf{r}_{t-1},\mathbf{S}_{t})=\mathcal{N}({\mbox{\boldmath$\theta$}};\hat{{\mbox{\boldmath$\theta$}}}_{t-1},{\mbox{\boldmath$\Sigma$}}_{t-1})$,
where $\hat{{\mbox{\boldmath$\theta$}}}_{t-1}$ and
${\mbox{\boldmath$\Sigma$}}_{t-1}$ are the posterior mean and the posterior
covariance matrix up until slot $t-1$, respectively; and with initial
$\bm{\theta}_{0}=\mathbf{0}_{2D}$ and
$\bm{\Sigma}_{0}=\sigma_{\theta}^{2}\mathbf{I}_{2D}$ based on (12). Before
receiving future rewards $\\{r_{t},r_{t+1},\ldots\\}$ that determine the value
function $\check{v}(\cdot)$ at $\mathbf{s}_{t}$, the agent predicts the pdf of
$\check{v}_{t}:=\check{v}(\mathbf{s}_{t})$ using

$p(\check{v}_{t}|\mathbf{r}_{t-1},\mathbf{S}_{t})=\mathcal{N}(\check{v}_{t};\hat{\check{v}}_{t|t-1},\sigma_{t|t-1}^{2})$

where the predictive mean and variance are

$\displaystyle\hat{\check{v}}_{t|t-1}$
$\displaystyle={\mbox{\boldmath$\phi$}}_{\zeta}^{\top}(\mathbf{s}_{t})\hat{{\mbox{\boldmath$\theta$}}}_{t-1}$
(16a) $\displaystyle\sigma_{t|t-1}^{2}$
$\displaystyle={\mbox{\boldmath$\phi$}}_{\zeta}^{\top}(\mathbf{s}_{t}){\mbox{\boldmath$\Sigma$}}_{t-1}{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t}).$
(16b)

Upon taking action $a_{t}\in\mathcal{A}$ by following given policy $\mu$ at
the beginning of time slot $t$, the agent receives reward $r_{t}$ and observes
the new state of the environment $\mathbf{s}_{t+1}$, based on which the
posterior pdf of $\theta$ can be found by invoking the Bayes’ rule as

$\displaystyle p({\mbox{\boldmath$\theta$}}|\mathbf{r}_{t},\mathbf{S}_{t+1})$
$\displaystyle=\frac{p({\mbox{\boldmath$\theta$}}|\mathbf{r}_{t-1},\mathbf{S}_{t})p(r_{t}|{\mbox{\boldmath$\theta$}},\mathbf{S}_{t+1})}{p(r_{t}|\mathbf{r}_{t-1})}$
$\displaystyle=\mathcal{N}({\mbox{\boldmath$\theta$}};\hat{{\mbox{\boldmath$\theta$}}}_{t},{\mbox{\boldmath$\Sigma$}}_{t})$

where the moments are recursively updated as follows

$\displaystyle\hat{{\mbox{\boldmath$\theta$}}}_{t}$
$\displaystyle=\hat{{\mbox{\boldmath$\theta$}}}_{t-1}+\frac{{\mbox{\boldmath$\Sigma$}}_{t-1}\mathbf{h}_{t}(r_{t}-\mathbf{h}_{t}^{\top}\hat{{\mbox{\boldmath$\theta$}}}_{t-1})}{\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t-1}\mathbf{h}_{t}+\sigma_{n}^{2}}$
(17a) $\displaystyle{\mbox{\boldmath$\Sigma$}}_{t}$
$\displaystyle={\mbox{\boldmath$\Sigma$}}_{t-1}-\frac{{\mbox{\boldmath$\Sigma$}}_{t-1}\mathbf{h}_{t}\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t-1}}{\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t-1}\mathbf{h}_{t}+\sigma_{n}^{2}}.$
(17b)

Besides offering a means of quantifying the uncertainty of the estimate in
(17a), the covariance matrix (17b) also adjusts the learning rate adaptively
from a stochastic optimization point of view. Moreover, it is evident that the
computational overhead over $t$ slots is $\mathcal{O}(t(2D)^{2})$ while the
storage requirement is $\mathcal{O}((2D)^{2})$. Hence, scalability is also
guaranteed for our OS-GPTD scheme.

##  4 Worst-case performance analysis

The value function estimator along with its variance in (16) were derived
based on the assumed generative model (4) and a GP prior on $v(\cdot)$. One
would be interested in benchmarking the performance of the sequence of
estimates (17a) even if such modeling assumptions are violated. Considering
the worst-case scenario in which reward-state pairs are chosen by an
adversary, performance is commonly evaluated in terms of the so-called regret.
Regret analysis has well-documented merits in online learning and optimization
[10]. For a given loss, it benchmarks the performance of an online learner
(a.k.a., algorithm) against the best fixed strategy that can utilize data in
hindsight. In our online TD learning context, the loss can be simply selected
as Bellman’s error. Beyond establishing logarithmic regret on this error, we
can further bound the long-term reward prediction error of OS-GPTD relative to
the benchmark value function.

Worst-case Bellman’s error and the long-term reward prediction error will be
analyzed in the rest of this section. It is worth stressing that for a given
trajectory of state-reward pairs, our analyses pertain to the sequence of
parameter estimates produced by (17), without any probabilistic assumptions
stated in the previous section.

###  4.1 Bellman error regret analysis

Regret analysis of our OS-GPTD algorithm will be quantified first in terms of
Bellman’s error, which demonstrates the self-consistency of the value function
evaluated at two successive states. For the RF-based OS-GPTD, the
instantaneous Bellman error at slot $t$ for any $\theta$ is simply
$l_{t}({\mbox{\boldmath$\theta$}}):=(r_{t}-\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\theta$}})^{2}$
[28, 27]. With batch data in hindsight, the benchmark function $v^{*}$ belongs
to the $\kappa$-induced reproducing kernel Hilbert space (RKHS), as we have
seen in (8a). The Bellman error functional of $v^{\ast}$ is thus
$l_{t}(v^{*}):=(r_{t}-v^{*}(\mathbf{s}_{t})+\gamma
v^{*}(\mathbf{s}_{t+1}))^{2}$. Hence, the cumulative regret of OS-GPTD over
$T$ slots is

$\mathcal{R}(T):=\sum_{t=1}^{T}l_{t}(\hat{{\mbox{\boldmath$\theta$}}}_{t-1})-\sum_{t=1}^{T}l_{t}(v^{*}).$
(18)

Let ${\mbox{\boldmath$\theta$}}^{*}$ denote the best fixed RF-based parameter
vector given the entire state-reward trajectory. The regret in (18) can then
be readily decomposed as

$\mathcal{R}(T)=\underbrace{\sum_{t=1}^{T}\\!l_{t}(\hat{{\mbox{\boldmath$\theta$}}}_{t-1})\\!-\\!\\!\sum_{t=1}^{T}\\!l_{t}({\mbox{\boldmath$\theta$}}^{*})\\!}_{:=\mathcal{R}_{1}(T)}+\\!\underbrace{\sum_{t=1}^{T}\\!l_{t}({\mbox{\boldmath$\theta$}}^{*})\\!-\\!\\!\sum_{t=1}^{T}\\!l_{t}(v^{*})}_{:=\mathcal{R}_{2}(T)}$

where $\mathcal{R}_{1}(T)$ corresponds to the regret incurred by confining
ourselves to learn within the class of RF-induced parametric functions, and
$\mathcal{R}_{2}(T)$ captures the accumulated RF-based approximation error of
the kernel. To upper-bound $\mathcal{R}_{1}(T)$ and $\mathcal{R}_{2}(T)$, we
will start with two standard assumptions.

  * as1)

The Bellman error is bounded by $B_{e}$; and,

  * as2)

The kernel $\bar{\kappa}$ is shift-invariant, standardized, and bounded such
that $\bar{\kappa}(\mathbf{s},\mathbf{s}^{\prime})\leq 1$,
$\forall\mathbf{s},\mathbf{s}^{\prime}$.

Boundedness under as1) is always satisfied for bounded rewards and state-
conditioned features, when the RF-based parameter vector or the RKHS-based
function is bounded as well. On the other hand, conditions in as2) hold for a
broad family of kernels including e.g., Gaussian, Laplace, and Cauchy ones
[22].

The ensuing lemma upper bounds the regret $\mathcal{R}_{1}(T)$, whose proof is
relegated to Appendix A of the supplementary material.

###### Lemma 1

Under as1), the following bound holds for the cumulative Bellman error
incurred by the OS-GPTD relative to its counterpart by the best fixed
parameter vector ${\mbox{\boldmath$\theta$}}^{*}$ in hindsight

$\mathcal{R}_{1}(T)\leq
2DB_{e}^{2}\log\\!\left(\frac{(1+\gamma)^{2}\sigma_{\theta}^{2}}{\sigma_{n}^{2}}T+1\right)+\frac{\|{\mbox{\boldmath$\theta$}}^{*}\|_{2}^{2}}{2\sigma_{\theta}^{2}}.$
(19)

Lemma 1 establishes the regret bound concerning the cumulative Bellman error
from the online RF-based function estimator relative to that from the best
fixed one with data in hindsight. In the next theorem, we further account for
the difference of the Bellman errors arising from approximating the RKHS-based
function estimator with the RF-based linear one. This will lead us to a bound
for the overall regret $\mathcal{R}(T)$ in (18).

###### Theorem 1

Under as1)-as2) and with $v^{*}$ belonging to the RKHS $\mathcal{H}$ induced
by $\kappa$, for any fixed $\epsilon>0$, the following bound holds with
probability at least
$1-\frac{2^{8}\sigma_{\zeta}^{2}}{\epsilon^{2}}\,{e}^{-\frac{D\epsilon^{2}}{4d+8}}$

$\displaystyle\sum_{t=1}^{T}l_{t}(\hat{{\mbox{\boldmath$\theta$}}}_{t-1})\\!-\\!\sum_{t=1}^{T}l_{t}(v^{*})\leq
2(1+\gamma)\epsilon B_{e}CT$
$\displaystyle~{}+\frac{(1+\epsilon)C^{2}}{2\sigma_{\theta}^{2}}+2DB_{e}^{2}\log\\!\left(\frac{(1+\gamma)^{2}\sigma_{\theta}^{2}}{\sigma_{n}^{2}}T+1\right)$
(20)

where $C>0$ is some constant, and
$\sigma_{\zeta}^{2}:=\mathbb{E}_{\pi_{\bar{\kappa}}}[\|{\mbox{\boldmath$\zeta$}}\|^{2}]$
is the second-order moment of the RF vector $\zeta$. Setting
$\epsilon=\mathcal{O}(\log T/T)$, the regret in (20) boils down to

$\mathcal{R}(T)=\mathcal{O}(\log T).$ (21)

With the total regret over $T$ slots being logarithmic in $T$, our proposed
OS-GPTD algorithm achieves _no regret_ on average in terms of the Bellman
error.

###  4.2 Long-term reward prediction error bound

For a given state-reward trajectory, the empirical value function or the
discounted sum of future rewards at slot $t$ is
$\bar{v}_{t}:=\sum_{\tau=t}^{\infty}\gamma^{\tau-t}r_{\tau}$. OS-GPTD then
incurs online prediction error of the long-term reward given by
$e_{t}:=\hat{\check{v}}_{t|t-1}-\bar{v}_{t}$, where $\hat{\check{v}}_{t|t-1}$
is the discounted future reward predicted as in (16a) based on the historical
rewards up to slot $t-1$. Although the no-regret property of OS-GPTD in the
Bellman error demonstrates that the estimated value function is self-
consistent at any two consecutive states, it does not necessarily imply small
long-term reward prediction errors, as it is also noted by [27]. To further
benchmark the performance of OS-GPTD regarding its long-term reward prediction
error, we will compare the cumulative online prediction error
$\sum_{t=1}^{T}e_{t}$ with its counterpart from the best fixed function
$v^{*}$, namely, $\sum_{t=1}^{T}e_{t}^{*}$ with
$e_{t}^{*}:=v^{*}(\mathbf{s}_{t})-\bar{v}_{t}$. While being not necessary, the
following performance bound conveniently builds upon the no-regret property
established in the previous subsection.

###### Theorem 2

Let as1)-as2) hold, and consider the sequence of RF-based linear value
function estimators produced recursively by (17a). The next bound holds w.h.p.
on the average online prediction error of the long-term rewards relative to
its counterpart from the best fixed estimator $v^{*}$, which belongs to the
$\kappa$-induced RKHS $\mathcal{H}$

$\lim_{T\rightarrow\infty}\frac{\sum_{t=1}^{T}e_{t}^{2}}{T}\leq\frac{2(1+\gamma)^{2}}{(1-\gamma)^{2}}\frac{\sum_{t=1}^{T}e_{t}^{*2}}{T}.$
(22)

##  5 Scalable TD learning with ensemble GPs

While a scalable approach has been devised in Section 3 to estimate the state
value function on-the-fly, its performance hinges largely on the _preselected_
kernel of the GP prior, which confines expressiveness of the function space.
Aiming at constructing a richer function space[17, 16], we advocate in this
section an _ensemble_ of $M$ GPs (experts) for the agent, each postulating a
unique GP prior on $v$, denoted as
$v|m\sim\mathcal{GP}(0,\kappa^{m}(\mathbf{s},\mathbf{s}^{\prime}))$, where
$m\in\mathcal{M}:=\\{1,\ldots,M\\}$ is the expert index and $\kappa^{m}$ is a
shift-invariant kernel selected from a prescribed dictionary
$\mathcal{K}:=\\{\kappa^{1},\ldots,\kappa^{M}\\}$. Per expert $m$, the prior
pdf of the value function $v$ evaluated at $\mathbf{S}_{t}$ is

$p(\mathbf{v}_{t}|z=m,\mathbf{S}_{t})\\!=\\!\mathcal{N}(\mathbf{v}_{t};{\bf
0}_{t},{\bf K}^{m}_{t})$

where the hidden association variable $z\in\mathcal{M}$ is introduced to
indicate the existence of each expert in the ensemble (E) GP, and ${\bf
K}_{t}^{m}$ is the covariance matrix induced by $\kappa^{m}$. Considering the
$M$ GP experts jointly, the EGP priors on $\mathbf{v}_{t}$ from the Gaussian
mixture (GM) given by

$p(\mathbf{v}_{t}|\mathbf{S}_{t})=\sum_{m=1}^{M}w^{m}\mathcal{N}(\mathbf{v}_{t};{\bf
0}_{t},\mathbf{K}^{m}_{t})$ (23)

where the weights $\\{w^{m}\\}_{m=1}^{M}$ are viewed as the prior
probabilities of each GP expert to be present in the EGPs.

As with a single GP prior or expert, the computational complexity of exact
function estimation with the EGP prior (23) will be intractable when the
number of the state-rewards grows large. Again, we will leverage RF-based
kernel approximants for the value function to effect scalability in the EGP
prior, contributing to our novel OS-EGPTD approach.

Seeking a tractable function predictor, each GP expert $m$ resorts to the
$\kappa^{m}$ induced RF-based approximant, where
$\kappa^{m}=\sigma_{\theta^{m}}^{2}\bar{\kappa}^{m}$. Upon drawing vectors
$\\{{\mbox{\boldmath$\zeta$}}_{i}^{m}\\}_{i=1}^{D}$ i.i.d. from
$\pi_{\bar{\kappa}}^{m}({\mbox{\boldmath$\zeta$}})$, expert $m$ relies on the
RF vector (11) and the per-expert parameter vector
${\mbox{\boldmath$\theta$}}^{m}$ to yield the generative model for the linear
approximant $\check{v}$

$\displaystyle{p}({\mbox{\boldmath$\theta$}}^{m}|z=m)$
$\displaystyle=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\mathbf{0}_{2D},\sigma_{\theta^{m}}^{2}\mathbf{I}_{2D})$
(24a) $\displaystyle
p(\check{v}(\mathbf{s})|{\mbox{\boldmath$\theta$}}^{m},z=m)$
$\displaystyle=\delta\big{(}\check{v}(\mathbf{s})-{\mbox{\boldmath$\phi$}}_{\zeta}^{m\top}(\mathbf{s}){\mbox{\boldmath$\theta$}}^{m}\big{)}$
(24b)

where $\delta$ is Kronecker’s delta. This linear approximant implies that the
data likelihood at slot $t$ for expert $m$ is

$p(r_{t}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{S}_{t+1},z=m)=\mathcal{N}\big{(}r_{t};\mathbf{h}_{t}^{m\top}{\mbox{\boldmath$\theta$}}^{m},\sigma_{n}^{2}\big{)}$
(25)

where the per-expert $m$ observation vector is

$\mathbf{h}_{t}^{m}:={\mbox{\boldmath$\phi$}}^{m}_{\zeta}(\mathbf{s}_{t})-\gamma{\mbox{\boldmath$\phi$}}^{m}_{\zeta}(\mathbf{s}_{t+1}).$
(26)

Aided by the Gaussian prior for ${\mbox{\boldmath$\theta$}}^{m}$ (cf. (24a)),
this linear Gaussian likelihood (25) will lead to a Gaussian posterior
$p({\mbox{\boldmath$\theta$}}^{m}|\mathbf{r}_{t},\mathbf{S}_{t+1},z=m)=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m},{\mbox{\boldmath$\Sigma$}}^{m}_{t})$
characterized by its mean $\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m}$ and
covariance matrix ${\mbox{\boldmath$\Sigma$}}^{m}_{t}$. The agent will rely on
$w_{t}^{m}:={\rm Pr}(z=m|\mathbf{r}_{t},\mathbf{S}_{t+1})$ to account for the
contribution of expert $m$. With sequential pairs $(r_{t},\mathbf{s}_{t+1})$
collected by the agent, the $M$ triplets
$\\{w_{t}^{m},\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m},{\mbox{\boldmath$\Sigma$}}^{m}_{t}\\}_{m=1}^{M}$
summarizing the past information about the state value function, will be
updated by alternating between prediction and correction steps as follows.

###  5.1 Prediction

Upon collecting
$\\{w_{t-1}^{m},\hat{{\mbox{\boldmath$\theta$}}}_{t-1}^{m},{\mbox{\boldmath$\Sigma$}}^{m}_{t-1}\\}_{m=1}^{M}$
from the $M$ GP experts at the end of slot $t-1$, the agent predicts the per-
expert pdf of $v(\mathbf{s}_{t})$ through

$\displaystyle
p(\check{v}_{t}|z\\!=\\!m,\mathbf{r}_{t-1},\mathbf{S}_{t})\\!=\\!\\!\int\\!\\!p(\check{v}_{t},{\mbox{\boldmath$\theta$}}^{m}|z\\!=\\!m,\mathbf{r}_{t-1},\mathbf{S}_{t})d{\mbox{\boldmath$\theta$}}^{m}$
$\displaystyle=\\!\int\\!\delta\big{(}\check{v}_{t}-{\mbox{\boldmath$\phi$}}_{\zeta}^{m\top}(\mathbf{s}_{t}){\mbox{\boldmath$\theta$}}^{m}\big{)}p({\mbox{\boldmath$\theta$}}^{m}|\mathbf{r}_{t-1},\mathbf{S}_{t})d{\mbox{\boldmath$\theta$}}^{m}$
$\displaystyle=\mathcal{N}\big{(}\check{v}_{t};\hat{\check{v}}_{t|t-1}^{m},(\sigma^{m}_{t|t-1})^{2}\big{)}$
(27)

where the predictive mean and variance are

$\displaystyle\hat{\check{v}}_{t|t-1}^{m}$
$\displaystyle={\mbox{\boldmath$\phi$}}_{\zeta}^{m\top}(\mathbf{s}_{t})\hat{{\mbox{\boldmath$\theta$}}}_{t-1}^{m}$
(28a) $\displaystyle(\sigma^{m}_{t|t-1})^{2}$
$\displaystyle={\mbox{\boldmath$\phi$}}_{\zeta}^{m\top}(\mathbf{s}_{t}){\mbox{\boldmath$\Sigma$}}^{m}_{t-1}{\mbox{\boldmath$\phi$}}_{\zeta}^{m}(\mathbf{s}_{t}).$
(28b)

The sum and product probability rules further enable the agent to obtain the
ensemble predictive pdf of $v(\mathbf{s}_{t})$ as

$\displaystyle
p(\check{v}(\mathbf{s}_{t})|\mathbf{r}_{t-1},\mathbf{S}_{t})=\sum_{m=1}^{M}p(\check{v}(\mathbf{s}_{t}),z=m|\mathbf{r}_{t-1},\mathbf{S}_{t})$
$\displaystyle=\sum_{m=1}^{M}w_{t-1}^{m}\mathcal{N}\big{(}\check{v}(\mathbf{s}_{t});\hat{\check{v}}_{t|t-1}^{m},(\sigma^{m}_{t|t-1})^{2}\big{)}$
(29)

which is a weighted superposition of the predictive pdfs from all experts. The
minimum mean-square error predictor of $v(\mathbf{s}_{t})$ can be obtained
along with its variance as

$\displaystyle\bar{\check{v}}_{t|t-1}\\!$
$\displaystyle=\\!\sum_{m=1}^{M}w_{t-1}^{m}\hat{\check{v}}_{t|t-1}^{m}$ (30a)
$\displaystyle\bar{\sigma}_{t|t-1}^{2}\\!$
$\displaystyle=\\!\sum_{s=1}^{m}w_{t-1}^{m}\big{[}(\sigma^{m}_{t|t-1})^{2}\\!+\\!(\bar{\check{v}}_{t|t-1}\\!-\\!\hat{\check{v}}_{t|t-1}^{m})^{2}\big{]}.$
(30b)

Algorithm 1 OS-EGPTD

Input: $\kappa^{m}$, $\sigma_{\theta^{m}}^{2}$, $m=1,\ldots,M$,
$\sigma_{n}^{2}$, $D$, $\gamma$ and $\mathbf{s}_{1}$.

Initialization:

for m = 1, 2, …, M do

Draw $D$ RF vectors $\\{{\mbox{\boldmath$\zeta$}}_{i}^{m}\\}_{i=1}^{D}$;

$w_{0}^{m}=1/M$; $\hat{{\mbox{\boldmath$\theta$}}}_{0}^{m}=\mathbf{0}_{2D}$;
${\mbox{\boldmath$\Sigma$}}_{0}^{m}=\sigma_{\theta^{m}}^{2}\mathbf{I}_{2D}$;

end for

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

Find predictive mean/variance of $\check{v}(\mathbf{s}_{t})$ via (30);

Receive $(r_{t},\mathbf{s}_{t+1})$;

for $m=1,2,\ldots,M$ do

Construct RF ${\mbox{\boldmath$\phi$}}_{\zeta}^{m}(\mathbf{s}_{t+1})$ via
(11);

Construct the observation vector $\mathbf{h}^{m}_{t}$ via (26);

Update per-expert pdf of ${\mbox{\boldmath$\theta$}}^{m}$ via (31);

Update $w_{t}^{m}$ via (33);

end for

end for

###  5.2 Correction

After receiving $r_{t}$ and $\mathbf{s}_{t+1}$ at slot $t$, GP expert $m$
employs the data likelihood (25) to propagate the posterior pdf as

$\displaystyle{p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{r}_{t},\mathbf{S}_{t+1},z=m)$
(31)
$\displaystyle=\frac{{p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{r}_{t-1},\mathbf{S}_{t},z=m){p}(r_{t}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{S}_{t},z=m)}{{p}(r_{t}|\mathbf{r}_{t-1},\mathbf{S}_{t+1},z=m)}$
$\displaystyle=\mathcal{N}({\mbox{\boldmath$\theta$}}^{m};\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m},{\mbox{\boldmath$\Sigma$}}^{m}_{t})$

where the mean and covariance matrix are updated using

$\displaystyle\hat{{\mbox{\boldmath$\theta$}}}_{t}^{m}$
$\displaystyle=\hat{{\mbox{\boldmath$\theta$}}}_{t-1}^{m}+\frac{{\mbox{\boldmath$\Sigma$}}_{t-1}^{m}\mathbf{h}_{t}^{m}\big{(}r_{t}-\mathbf{h}_{t}^{m\top}\hat{{\mbox{\boldmath$\theta$}}}^{m}_{t-1}\big{)}}{\mathbf{h}_{t}^{m\top}{\mbox{\boldmath$\Sigma$}}_{t-1}^{m}\mathbf{h}_{t}^{m}+\sigma_{n}^{2}}$
(32a) $\displaystyle{\mbox{\boldmath$\Sigma$}}_{t}^{m}$
$\displaystyle={\mbox{\boldmath$\Sigma$}}_{t-1}^{m}-\frac{{\mbox{\boldmath$\Sigma$}}_{t-1}^{m}\mathbf{h}_{t}^{m}\mathbf{h}_{t}^{m\top}{\mbox{\boldmath$\Sigma$}}_{t-1}^{m}}{\mathbf{h}_{t}^{m\top}{\mbox{\boldmath$\Sigma$}}_{t-1}^{m}\mathbf{h}_{t}^{m}+\sigma_{n}^{2}}.$
(32b)

Meanwhile, the per-expert posterior weight is corrected using Bayes’ rule as

$\displaystyle w_{t}^{m}={\rm Pr}(z=m|\mathbf{r}_{t},\mathbf{S}_{t+1})$
$\displaystyle=\frac{{\rm
Pr}(z=m|\mathbf{r}_{t-1},\mathbf{S}_{t}){p}(r_{t}|\mathbf{r}_{t-1},\mathbf{S}_{t+1},z=m)}{\sum_{m=1}^{M}{\rm
Pr}(z=m|\mathbf{r}_{t-1},\mathbf{S}_{t}){p}(r_{t}|\mathbf{r}_{t-1},\mathbf{S}_{t+1},z=m)}$

which, upon utilizing

$\displaystyle{p}(r_{t}|\mathbf{r}_{t-1},\mathbf{S}_{t+1},z=m)$
$\displaystyle=\int\\!{p}(r_{t}|{\mbox{\boldmath$\theta$}}^{m},\mathbf{S}_{t+1},z=m){p}({\mbox{\boldmath$\theta$}}^{m}|\mathbf{r}_{t-1},\mathbf{S}_{t},z=m)d{\mbox{\boldmath$\theta$}}^{m}$
$\displaystyle=\mathcal{N}\big{(}r_{t};\mathbf{h}_{t}^{m\top}\hat{{\mbox{\boldmath$\theta$}}}^{m}_{t-1},\mathbf{h}_{t}^{m\top}{\mbox{\boldmath$\Sigma$}}^{m}_{t-1}\mathbf{h}_{t}^{m}+\sigma_{n}^{2}\big{)}$

boils down to

$\displaystyle w_{t}^{m}=$ (33)
$\displaystyle\frac{w_{t-1}^{m}\mathcal{N}\big{(}r_{t};\mathbf{h}_{t}^{m\top}\hat{{\mbox{\boldmath$\theta$}}}^{m}_{t-1},\mathbf{h}_{t}^{m\top}{\mbox{\boldmath$\Sigma$}}^{m}_{t-1}\mathbf{h}_{t}^{m}+\sigma_{n}^{2}\big{)}}{\sum_{m^{\prime}=1}^{M}w_{t-1}^{m^{\prime}}\mathcal{N}\big{(}r_{t};\mathbf{h}_{t}^{m^{\prime}\top}\hat{{\mbox{\boldmath$\theta$}}}^{m^{\prime}}_{t-1},\mathbf{h}_{t}^{m^{\prime}\top}{\mbox{\boldmath$\Sigma$}}^{m^{\prime}}_{t-1}\mathbf{h}_{t}^{m^{\prime}}\\!+\\!\sigma_{n}^{2}\big{)}}.$

Accounting for all $M$-expert updates, the novel OS-EGPTD scheme, summarized
in Algorithm 1, runs $M$ OS-GPTDs in parallel, thus incurring a per-slot
complexity of $\mathcal{O}(M(2D)^{2})$. Hence, scalability is not compromised
by the ensemble approach that also endows the value function $v$ with richer
expressiveness.

###### Remark 1

Although the focus so far has been on the state value function estimation, our
OS-(E)GPTD methods can be readily employed to estimate state-action value
functions, thus accommodating policy improvement without any additional
knowledge on the environment or MDP.

$200$$400$$600$$800$$1{,}000$$2.8$$3$$3.2$$3.4$Slot $t$Average prediction
errorSPA-GPTDOS-GPTDOS-EGPTD Figure 1: Average prediction error on Random
Walk. $200$$400$$600$$800$$1{,}000$$1.5$$2$$2.5$$3$Slot $t$Average Bellman
errorSPA-GPTDOS-GPTDOS-EGPTD Figure 2: Average Bellman error on Random Walk.

##  6 Experiments

In this section, our novel OS-(E)GPTD schemes are tested using two simulated
policy evaluation benchmarks. To better assess their performance, we compare
them with the online sparsification approach based GPTD devised in [7],
henceforth abbreviated as SPA-GPTD. Performance metrics include the average
cumulative long-term reward prediction error
$\sqrt{\sum_{\tau=1}^{t}e_{\tau}^{2}/t}$ and Bellman error
$\sqrt{\sum_{\tau=1}^{t}l_{\tau}(\hat{{\mbox{\boldmath$\theta$}}}_{\tau-1})/t}$
up to slot $t$, as well as the average runtime. The kernel dictionary of OS-
EGPTD comprises three standardized Gaussian kernels, whose characteristic
lengthscales ($\sigma$’s) are $\mathcal{K}:=\\{0.1,1,10\\}$, while both OS-
GPTD and SPA-GPTD employ a single Gaussian kernel with $\sigma\in\mathcal{K}$
specified below. The reported results are averaged over $100$ trajectories.

$500$$1{,}000$$1{,}500$$2{,}000$$10$$20$$30$$40$Slot $t$Average prediction
errorSPA-GPTDOS-GPTDOS-EGPTD Figure 3: Average prediction error on Puddle
World. $0$$500$$1{,}000$$1{,}500$$2{,}000$$6$$8$$10$$12$Slot $t$Average
Bellman errorSPA-GPTDOS-GPTDOS-EGPTD Figure 4: Average Bellman error on Puddle
World .

Random Walk. This problem is a variant of the Hall problem introduced in [1],
and also used in [14, 27]. The state space consists of $N=50$ states, which
are linked as a ring without terminal states. Per slot $t$, the current state
transitions clockwise to the next state with probability $1$. Each state is
associated with a vector $\mathbf{s}_{i}\in\mathbb{R}^{10}$, $i=1,\ldots,50$,
randomly generated from the standard normal distribution. Instantaneous
rewards $r_{t}$ are uniformly sampled from $[-3,3]$, and each trajectory
contains at most $T=1,000$ slots. For OS-GPTD and SPA-GPTD, the single GP
prior uses a standardized Gaussian kernel with $\sigma=0.1$. The discount
factor $\gamma$ is set to $0.75$, and $\sigma_{n}^{2}=0.01$.

RandomWalkPuddleWorld$10^{-1}$$10^{0}$$10^{1}$Normalized RuntimeSPA-GPTDOS-
EGPTDOS-GPTD Figure 5: Average CPU times.

PuddleWorld was adapted from [28] with its state space being a unit square
with ‘puddles.’ The agent’s state is represented by its $x$ and $y$
coordinates. In each episode, the starting state is uniformly sampled in the
region $[0,0.2]\times[0,0.2]$. The policy per step selects ‘go north or east’
with probability $0.5$. The terminal region is a circle with radius $0.1$
centered at $(1,1)$. The reward $r_{t}$ is $-1$ if the agent does not step in
a puddle, and decreases to $-2$ as the agent steps into the puddles or
boundaries. Both OS-GPTD and SPA-GPTD employ a Gaussian kernel with
$\sigma=1$, while $\sigma_{n}^{2}=0.001$ and $\gamma=0.7$.

As evident from Figs. 1-4, our OS-EGPTD enjoys lower prediction and Bellman
errors than single GP-based alternatives. Although achieving similar
estimation performances, the RF-based OS-GPTD scheme is much more
computationally efficient than SPA-GPTD as shown in Fig. 5. Additional tests
to validate the regret bound in Theorem 1 can be found in the supplementary
file.

##  7 Conclusions

Building on the GPTD framework for value estimation in RL, this paper
developed OS-(E)GPTD schemes, leveraging RF-based kernel approximants using an
(ensemble of) GP prior(s) to effect online function estimation along with
performance quantification. To account for violations of the GPTD modeling
assumptions in practice, worst-case analyses have been performed in terms of
both the Bellman error as well as the long-term reward prediction error, that
benchmark the performance of online function estimators produced by OS-GPTD
relative to the best fixed one with batch data in hindsight. Numerical tests
corroborated the merits of the novel OS-(E)GPTD schemes.

## References

  * [1] L. Baird, “Residual algorithms: Reinforcement learning with function approximation,” in _ICML_ , 1995, pp. 30–37.
  * [2] J. Bhandari, D. Russo, and R. Singal, “A finite time analysis of temporal difference learning with linear function approximation,” in _COLT_ , 2019, pp. 1691–1692.
  * [3] S. Bhatnagar, D. Precup, D. Silver, R. S. Sutton, H. R. Maei, and C. Szepesvári, “Convergent temporal-difference learning with arbitrary smooth function approximation,” in _NeurIPS_ , 2009, pp. 1204–1212.
  * [4] G. Dalal, B. Szörényi, G. Thoppe, and S. Mannor, “Finite sample analyses for TD(0) with function approximation,” in _AAAI_ , 2018, pp. 6144–6152.
  * [5] M. Deisenroth and C. E. Rasmussen, “PILCO: A model-based and data-efficient approach to policy search,” in _ICML_ , 2011, pp. 465–472.
  * [6] M. P. Deisenroth, D. Fox, and C. E. Rasmussen, “Gaussian processes for data-efficient learning in robotics and control,” _IEEE TPAMI_ , vol. 37, no. 2, pp. 408–423, 2017.
  * [7] Y. Engel, S. Mannor, and R. Meir, “Bayes meets Bellman: The Gaussian process approach to temporal difference learning,” in _ICML_ , 2003, pp. 154–161.
  * [8] ——, “Reinforcement learning with Gaussian processes,” in _ICML_ , 2005, pp. 201–208.
  * [9] R. Grande, T. Walsh, and J. How, “Sample efficient reinforcement learning with Gaussian processes,” in _ICML_ , 2014, pp. 1332–1340.
  * [10] E. Hazan, A. Agarwal, and S. Kale, “Logarithmic regret algorithms for online convex optimization,” _Mach. Learn._ , vol. 69, no. 2-3, pp. 169–192, 2007\.
  * [11] M. Jin, A. Damianou, P. Abbeel, and C. Spanos, “Inverse reinforcement learning via deep Gaussian process,” _UAI_ , 2017.
  * [12] M. Kuss and C. E. Rasmussen, “Gaussian processes in reinforcement learning,” in _NeurIPS_ , 2004, pp. 751–758.
  * [13] M. Lázaro-Gredilla, J. Quiñonero Candela, C. E. Rasmussen, and A. Figueiras-Vidal, “Sparse spectrum Gaussian process regression,” _JMLR_ , vol. 11, no. Jun, pp. 1865–1881, 2010.
  * [14] L. Li, “A worst-case comparison between temporal difference and residual gradient with linear function approximation,” in _ICML_ , 2008, pp. 560–567.
  * [15] Q. Lu and G. B. Giannakis, “Gaussian process temporal-difference learning with scalability and worst-case performance guarantees,” in _IEEE ICASSP_ , 2021, pp. 3485–3489.
  * [16] Q. Lu, G. Karanikolas, and G. B. Giannakis, “Incremental ensemble Gaussian processes,” _arXiv preprint arXiv:2110.06777_ , 2021.
  * [17] Q. Lu, G. Karanikolas, Y. Shen, and G. B. Giannakis, “Ensemble Gaussian processes with spectral features for online interactive learning with scalability,” in _AISTATS_ , Palermo, Italy, June 2020.
  * [18] J. Martin, J. Wang, and B. Englot, “Sparse gptd learning for marine robot navigation,” in _ICRL_ , 2018, pp. 179–189.
  * [19] V. Mnih, K. Kavukcuoglu, D. Silver _et al._ , “Human-level control through deep reinforcement learning,” _Nature_ , vol. 518, no. 7540, p. 529, 2015\.
  * [20] M. L. Puterman, _Markov Decision Processes: Discrete Stochastic Dynamic Programming_. John Wiley & Sons, 2014.
  * [21] J. Quiñonero-Candela and C. E. Rasmussen, “A unifying view of sparse approximate Gaussian process regression,” _JMLR_ , vol. 6, no. Dec, pp. 1939–1959, 2005.
  * [22] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” in _NeurIPS_ , 2008, pp. 1177–1184.
  * [23] C. E. Rasmussen and C. K. Williams, _Gaussian Processes for Machine Learning_. MIT press Cambridge, MA, 2006\.
  * [24] S. Sæmundsson, K. Hofmann, and M. Deisenroth, “Meta reinforcement learning with latent variable Gaussian processes,” in _UAI_ , vol. 34, 2018, pp. 642–652.
  * [25] R. E. Schapire and M. K. Warmuth, “On the worst-case analysis of temporal-difference learning algorithms,” _Mach. Learn._ , vol. 22, no. 1-3, pp. 95–121, 1996.
  * [26] R. Srikant and L. Ying, “Finite-time error bounds for linear stochastic approximation and TD learning,” in _COLT_ , 2019, pp. 1–11.
  * [27] W. Sun and J. A. Bagnell, “Online Bellman residual algorithms with predictive error guarantees,” in _UAI_ , 2015.
  * [28] R. S. Sutton and A. G. Barto, _Reinforcement Learning: An Introduction_. MIT press, 2018.
  * [29] J. Tsitsiklis and B. Van Roy, “An analysis of temporal-difference learning with function approximation,” _IEEE TAC_ , vol. 42, no. 5, 1997.

##  Appendix A Proof of Lemma 1

It will be useful to rewrite the update in (17) as

$\hat{{\mbox{\boldmath$\theta$}}}_{t}=\hat{{\mbox{\boldmath$\theta$}}}_{t-1}-\frac{1}{2\sigma_{n}^{2}}\mathbf{J}_{t}^{-1}\nabla_{t|t-1}$
(34)

where the gradient of the Bellman error at slot $t$ evaluated at
$\hat{{\mbox{\boldmath$\theta$}}}_{t-1}$ is

$\nabla_{t|t-1}:=2\mathbf{h}_{t}(\mathbf{h}_{t}^{\top}\hat{{\mbox{\boldmath$\theta$}}}_{t-1}-r_{t})$
(35)

and $\mathbf{J}_{t}:={\mbox{\boldmath$\Sigma$}}_{t}^{-1}$ denotes the
information matrix, which starting from
$\mathbf{J}_{0}:=\sigma_{\theta}^{-2}\mathbf{I}_{2D}$, can be propagated using

$\mathbf{J}_{t}=\mathbf{J}_{t-1}+\sigma_{n}^{-2}\mathbf{h}_{t}\mathbf{h}_{t}^{\top}.$
(36)

Invoking Taylor’s expansion, we can write

$\displaystyle
l_{t}(\hat{{\mbox{\boldmath$\theta$}}}_{t-1})-l_{t}({\mbox{\boldmath$\theta$}}^{*})$
$\displaystyle=$
$\displaystyle-\\!\nabla_{t|t-1}^{\top}({\mbox{\boldmath$\theta$}}^{*}\\!-\hat{{\mbox{\boldmath$\theta$}}}_{t-1}\\!)\\!-\\!({\mbox{\boldmath$\theta$}}^{*}\\!-\hat{{\mbox{\boldmath$\theta$}}}_{t-1}\\!)^{\top}\mathbf{h}_{t}\mathbf{h}_{t}^{\top}\\!({\mbox{\boldmath$\theta$}}^{*}\\!-\hat{{\mbox{\boldmath$\theta$}}}_{t-1}\\!)$
$\displaystyle\overset{(a)}{=}$
$\displaystyle-\\!\nabla_{t|t-1}^{\top}({\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t-1})\\!-\\!\sigma_{n}^{2}({\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t-1}\\!)^{\top}\mathbf{J}_{t}({\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t-1})$
$\displaystyle+\sigma_{n}^{2}({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{t-1}\\!)^{\top}\mathbf{J}_{t-1}({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{t-1})$

where $(a)$ utilizes
$\mathbf{h}_{t}\mathbf{h}_{t}^{\top}=\sigma_{n}^{2}(\mathbf{J}_{t}-\mathbf{J}_{t-1})$
from (36). The cumulative regret $\mathcal{R}_{1}(T)$ now becomes

$\displaystyle\mathcal{R}_{1}$
$\displaystyle(T)=\sum_{t=1}^{T}l_{t}(\hat{{\mbox{\boldmath$\theta$}}}_{t-1})-\sum_{t=1}^{T}l_{t}({\mbox{\boldmath$\theta$}}^{*})$
$\displaystyle=$
$\displaystyle\sum_{t=1}^{T}\Big{[}\\!-\\!\nabla_{t|t-1}^{\top}({\mbox{\boldmath$\theta$}}^{*}\\!-\hat{{\mbox{\boldmath$\theta$}}}_{t-1})\\!+\\!\sigma_{n}^{2}({\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t})^{\top}\\!\mathbf{J}_{t}({\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t})$
$\displaystyle\quad-\sigma_{n}^{2}({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{t-1})^{\top}\mathbf{J}_{t}({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{t-1})\Big{]}+{\mbox{\boldmath$\theta$}}^{*\top}\mathbf{J}_{0}{\mbox{\boldmath$\theta$}}^{*}$
$\displaystyle\quad-({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{T})^{\top}\mathbf{J}_{T}({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{T}).$
(37)

On the other hand, the update in (34) confirms that

$\displaystyle({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{t})^{\top}\mathbf{J}_{t}({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{t})$
$\displaystyle=\\!\bigg{(}\\!{\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t-1}\\!+\\!\frac{\mathbf{J}_{t}^{-1}\nabla_{t|t-1}}{2\sigma_{n}^{2}}\\!\\!\bigg{)}^{\\!\top}\\!\\!\\!\mathbf{J}_{t}\\!\bigg{(}\\!{\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t-1}\\!+\\!\frac{\mathbf{J}_{t}^{-1}\nabla_{t|t-1}}{2\sigma_{n}^{2}}\\!\bigg{)}$
$\displaystyle=({\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t-1})^{\top}\mathbf{J}_{t}({\mbox{\boldmath$\theta$}}^{*}\\!-\\!\hat{{\mbox{\boldmath$\theta$}}}_{t-1})+\frac{1}{\sigma_{n}^{2}}\nabla_{t|t-1}^{\top}({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{t-1})$
$\displaystyle\quad+\frac{1}{4\sigma_{n}^{4}}\nabla_{t|t-1}^{\top}\mathbf{J}_{t}^{-1}\nabla_{t|t-1}$

substituting which into (37), yields

$\displaystyle\mathcal{R}_{1}(T)$
$\displaystyle=\sum_{t=1}^{T}\frac{1}{4\sigma_{n}^{2}}\nabla_{t|t-1}^{\top}\mathbf{J}_{t}^{-1}\nabla_{t|t-1}+{\mbox{\boldmath$\theta$}}^{*\top}\mathbf{J}_{0}{\mbox{\boldmath$\theta$}}^{*}$
$\displaystyle\quad-({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{T})^{\top}\mathbf{J}_{T}({\mbox{\boldmath$\theta$}}^{*}-\hat{{\mbox{\boldmath$\theta$}}}_{T})$
$\displaystyle\overset{(a)}{\leq}\sum_{t=1}^{T}\frac{1}{4\sigma_{n}^{2}}\nabla_{t|t-1}^{\top}\mathbf{J}_{t}^{-1}\nabla_{t|t-1}+\frac{\|{\mbox{\boldmath$\theta$}}^{*}\|_{2}^{2}}{2\sigma_{\theta}^{2}}$
$\displaystyle\overset{(b)}{\leq}\sum_{t=1}^{T}\frac{4B_{e}^{2}}{4\sigma_{n}^{2}}\mathbf{h}_{t}^{\top}\mathbf{J}_{t}^{-1}\mathbf{h}_{t}+\frac{\|{\mbox{\boldmath$\theta$}}^{*}\|_{2}^{2}}{2\sigma_{\theta}^{2}}$
$\displaystyle\overset{(c)}{\leq}2DB_{e}^{2}\log\\!\left(\frac{(1+\gamma)^{2}\sigma_{\theta}^{2}}{\sigma_{n}^{2}}T+1\right)+\frac{\|{\mbox{\boldmath$\theta$}}^{*}\|_{2}^{2}}{2\sigma_{\theta}^{2}}$

where $(a)$ follows from
$({\mbox{\boldmath$\theta$}}^{*}-{\mbox{\boldmath$\theta$}}_{T})^{\top}\mathbf{J}_{T}({\mbox{\boldmath$\theta$}}^{*}-{\mbox{\boldmath$\theta$}}_{T})\geq
0$ and $\mathbf{J}_{0}=\sigma_{\theta}^{-2}\mathbf{I}_{2D}$, $(b)$ leverages
the gradient evaluation in (35), and $(c)$ relies on a standard result in [10,
Appendix 2] with $\|\mathbf{h}_{t}\|\leq 1+\gamma$ and
$|r_{t}-\mathbf{h}_{t}^{\top}\hat{{\mbox{\boldmath$\theta$}}}_{t-1}|\leq
B_{e}$ under our assumptions. This concludes the proof of Lemma 1.

##  Appendix B Proof of Theorem 1

For a given shift-invariant, standardized kernel $\bar{\kappa}^{s}$, the
maximum point-wise error of the RF kernel approximant can be uniformly bounded
with probability at least
$1-\frac{2^{8}\sigma_{\zeta}^{2}}{\epsilon^{2}}\,{e}^{-\frac{D\epsilon^{2}}{4d+8}}$
by [22]

$\sup_{\mathbf{s}_{i},\mathbf{s}_{j}\in\mathcal{S}}\;\left|{\mbox{\boldmath$\phi$}}^{\top}_{\zeta}(\mathbf{s}_{i}){\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{j})-\bar{\kappa}(\mathbf{s}_{i},\mathbf{s}_{j})\right|<\epsilon$
(38)

where $\epsilon$ is a given constant, $D$ is the number of RF vectors, $d$ is
the dimension of $\mathbf{s}$, and
$\sigma_{\zeta}^{2}=\mathbb{E}_{\pi_{\bar{\kappa}}}[\|{\mbox{\boldmath$\zeta$}}\|^{2}]$
is the second-order moment of the RF vector $\zeta$.

The optimal function estimator in $\mathcal{H}$ induced by $\kappa$ is given
by
$v^{*}(\mathbf{s}):=\sum_{t=1}^{T}\hat{\alpha}_{t}\kappa(\mathbf{s},\mathbf{s}_{t})=\sigma_{\theta}^{2}\sum_{t=1}^{T}\hat{\alpha}_{t}\bar{\kappa}(\mathbf{s},\mathbf{s}_{t})$;
and its RF-based approximant is
$\check{v}^{*}(\mathbf{s}):={\mbox{\boldmath$\phi$}}^{\top}_{\zeta}(\mathbf{s}){\mbox{\boldmath$\theta$}}^{*}$
with
${\mbox{\boldmath$\theta$}}^{*}:=\sigma_{\theta}^{2}\sum_{t=1}^{T}\hat{\alpha}_{t}{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t})$.
We then have that

$\displaystyle\big{|}l_{t}(\check{v}^{*})-l_{t}(v^{*})\big{|}$
$\displaystyle\overset{(a)}{\leq}$
$\displaystyle\,2B_{e}\big{|}\check{v}^{*}(\mathbf{s}_{t-1})-\gamma\check{v}^{*}(\mathbf{s}_{t})-[v^{*}(\mathbf{s}_{t-1})-\gamma
v^{*}(\mathbf{s}_{t})]\big{|}$ $\displaystyle\overset{}{\leq}$ $\displaystyle
2B_{e}\big{[}\big{|}\check{v}^{*}(\mathbf{s}_{t-1})-v^{*}(\mathbf{s}_{t-1})\big{|}+\gamma\big{|}\check{v}^{*}(\mathbf{s}_{t})-v^{*}(\mathbf{s}_{t})\big{|}\big{]}$
$\displaystyle\overset{}{=}$
$\displaystyle\,2B_{e}\\!\sigma_{\theta}^{2}\Big{[}\Big{|}\sum_{t^{\prime}=1}^{T}\hat{\alpha}_{t^{\prime}}\big{(}{\mbox{\boldmath$\phi$}}_{\zeta}^{\top}\\!(\mathbf{s}_{t-1}){\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t^{\prime}})-\bar{\kappa}(\mathbf{s}_{t-1},\mathbf{s}_{t^{\prime}})\big{)}\Big{|}$
$\displaystyle+\gamma\Big{|}\sum_{t^{\prime}=1}^{T}\hat{\alpha}_{t^{\prime}}\big{(}{\mbox{\boldmath$\phi$}}_{\zeta}^{\top}(\mathbf{s}_{t}){\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t^{\prime}})-\bar{\kappa}(\mathbf{s}_{t},\mathbf{s}_{t^{\prime}})\big{)}\Big{|}\Big{]}$
$\displaystyle\overset{}{\leq}$ $\displaystyle
2B_{e}\sigma_{\theta}^{2}\sum_{t^{\prime}=1}^{T}|\hat{\alpha}_{t^{\prime}}|\Big{(}\big{|}{\mbox{\boldmath$\phi$}}_{\zeta}^{\top}(\mathbf{s}_{t-1}){\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t^{\prime}})-\bar{\kappa}(\mathbf{s}_{t-1},\mathbf{s}_{t^{\prime}})\big{|}$
$\displaystyle+\gamma\big{|}{\mbox{\boldmath$\phi$}}_{\zeta}^{\top}(\mathbf{s}_{t}){\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t^{\prime}})-\bar{\kappa}(\mathbf{s}_{t},\mathbf{s}_{t^{\prime}})\big{|}\Big{)}$
(39)

where $(a)$ follows because $l_{t}(v)$ is convex wrt
$v({\mathbf{s}_{t}})-\gamma v({\mathbf{s}_{t+1}})$, with gradient
$2(v({\mathbf{s}_{t}})-\gamma v({\mathbf{s}_{t+1}})+r_{t})$, bounded by
$2B_{e}$. Combining with (38), we deduce that

$|l_{t}(\check{v}^{*})-l_{t}(v^{*})|\leq 2(1+\gamma)B_{e}C\epsilon,~{}~{}~{}\
{\rm w.h.p.}$

where $C:=\sum_{t=1}^{T}\sigma_{\theta}^{2}|\hat{\alpha}_{t}|$. It thus holds
that

$\sum_{t=1}^{T}l_{t}(\check{v}^{*})\\!-\\!\sum_{t=1}^{T}l_{t}(v^{*})\leq
2(1+\gamma)\epsilon B_{e}CT,\,{\rm w.h.p.}$ (40)

On the other hand, the uniform convergence bound in (38) along with as2) imply
that

$\underset{\mathbf{s}_{t},\mathbf{s}_{t^{\prime}}\in\mathcal{S}}{\sup}~{}~{}{\mbox{\boldmath$\phi$}}^{\top}_{\zeta}(\mathbf{s}_{t}){\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t^{\prime}})\leq
1+\epsilon,\quad{\rm w.h.p.}$ (41)

which leads to

$\displaystyle\|{\mbox{\boldmath$\theta$}}^{*}\|^{2}$
$\displaystyle\\!:=\Big{\|}\sigma_{\theta}^{2}\sum_{t=1}^{T}\hat{\alpha}_{t}{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t})\Big{\|}^{2}$
(42)
$\displaystyle=\sigma_{\theta}^{4}\sum_{t=1}^{T}\sum_{t^{\prime}=1}^{T}\hat{\alpha}_{t}\hat{\alpha}_{t^{\prime}}{\mbox{\boldmath$\phi$}}^{\top}_{\zeta}(\mathbf{s}_{t}){\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t^{\prime}})$
$\displaystyle\leq(1+\epsilon)C^{2}.$

Putting together the bounds in (42) and (40), along with Lemma 1, we arrive at
(20), which completes the proof of Theorem 1.

##  Appendix C Proof of Theorem 2

Recall first that [27, Thm. 3.3]

$\displaystyle(1-\gamma)^{2}\sum_{t=1}^{T}e_{t}^{2}\leq
2(1+\gamma)^{2}\sum_{t=1}^{T}e_{t}^{*2}+C^{\prime}$ (43)
$\displaystyle+\\!2\gamma^{2}\sum_{t=1}^{T}\\!(\hat{\check{v}}_{t}\big{(}\mathbf{s}_{t})-\hat{\check{v}}_{t-1}(\mathbf{s}_{t})\big{)}^{2}\\!+\\!2\sum_{t=1}^{T}\big{(}l_{t}(\hat{{\mbox{\boldmath$\theta$}}}_{t-1})\\!-\\!l_{t}(v^{*})\big{)}$

where $C^{\prime}$ is a constant. Using Theorem 1, it follows that

$\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{t=1}^{T}\big{(}l_{t}(\hat{{\mbox{\boldmath$\theta$}}}_{t-1})-l_{t}(v^{*})\big{)}=0.$
(44)

To deduce (22) from (43), it suffices to show that

$\lim_{T\rightarrow\infty}\frac{1}{T}\sum_{t=1}^{T}\big{(}\hat{\check{v}}_{t}(\mathbf{s}_{t})\\!-\\!\hat{\check{v}}_{t-1}(\mathbf{s}_{t})\big{)}^{2}=0.$
(45)

which is designated as the _online stability_ condition in e.g., [27]. To
establish (45), we need the following inequality

$\displaystyle(\hat{\check{v}}_{t}(\mathbf{s}_{t})-\hat{\check{v}}_{t-1}(\mathbf{s}_{t}))^{2}=({\mbox{\boldmath$\phi$}}_{\zeta}^{\top}(\mathbf{s}_{t})(\hat{{\mbox{\boldmath$\theta$}}}_{t}-\hat{{\mbox{\boldmath$\theta$}}}_{t-1}))^{2}$
$\displaystyle\overset{(a)}{\leq}\\!\|{\mbox{\boldmath$\phi$}}_{\zeta}(\mathbf{s}_{t})\|^{2}\|\hat{{\mbox{\boldmath$\theta$}}}_{t}-\hat{{\mbox{\boldmath$\theta$}}}_{t-1}\|^{2}\overset{}{=}\\!\left\|\frac{{\mbox{\boldmath$\Sigma$}}_{t-1}\mathbf{h}_{t}(r_{t}\\!-\\!\mathbf{h}_{t}^{\top}\hat{{\mbox{\boldmath$\theta$}}}_{t-1})}{\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t}\mathbf{h}_{t}\\!+\\!\sigma_{n}^{2}}\right\|^{2}$
$\displaystyle\overset{(b)}{\leq}\frac{B_{e}^{2}}{\sigma_{n}^{2}}\frac{\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t-1}{\mbox{\boldmath$\Sigma$}}_{t-1}\mathbf{h}_{t}}{\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t}\mathbf{h}_{t}\\!+\\!\sigma_{n}^{2}}=\frac{B_{e}^{2}}{\sigma_{n}^{2}}{\rm
Tr}\\!\left(\frac{{\mbox{\boldmath$\Sigma$}}_{t-1}\mathbf{h}_{t}\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t-1}}{\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t-1}\mathbf{h}_{t}\\!+\\!\sigma_{n}^{2}}\right)$
$\displaystyle\overset{(c)}{=}\frac{B_{e}^{2}}{\sigma_{n}^{2}}{\rm
Tr}({\mbox{\boldmath$\Sigma$}}_{t-1}-{\mbox{\boldmath$\Sigma$}}_{t})$ (46)

where $(a)$ arises from the Cauchy-Schwarz inequality; $(b)$ is due to as1)
and $\mathbf{h}_{t}^{\top}{\mbox{\boldmath$\Sigma$}}_{t}\mathbf{h}_{t}>0$; and
the identity $(c)$ employs the update of the covariance matrix in (17b).

Summing (46) from $t=1$ to $T$, yields

$\vspace{-.2cm}\sum_{t=1}^{T}\big{(}\hat{\check{v}}_{t}(\mathbf{s}_{t})-\hat{\check{v}}_{t-1}(\mathbf{s}_{t})\big{)}^{2}\leq\frac{B_{e}^{2}}{\sigma_{n}^{2}}{\rm
Tr}({\mbox{\boldmath$\Sigma$}}_{0}-{\mbox{\boldmath$\Sigma$}}_{T})$ (47)

which, in conjunction with (44), validates (22), thus completing the proof of
Theorem 2.

##  Appendix D Additional tests

$0$$200$$400$$600$$800$$1{,}000$$0$$200$$400$$600$$800$$1{,}000$$T$Regret
Figure 6: Regret of OS-GPTD vs. $T$ on the Random Walk problem.