# Dynamic Selection of p-Norm in Linear Adaptive Filtering via
Online Kernel-Based Reinforcement Learning

###### Abstract

This study addresses the problem of selecting dynamically, at each time
instance, the “optimal” p-norm to combat outliers in linear adaptive filtering
without any knowledge on the potentially time-varying probability density
function of the outliers. To this end, an online and data-driven framework is
designed via kernel-based reinforcement learning (KBRL). Novel Bellman
mappings on reproducing kernel Hilbert spaces (RKHSs) are introduced that need
no knowledge on transition probabilities of Markov decision processes, and are
nonexpansive with respect to the underlying Hilbertian norm. An approximate
policy-iteration framework is finally offered via the introduction of a
finite-dimensional affine superset of the fixed-point set of the proposed
Bellman mappings. The well-known “curse of dimensionality” in RKHSs is
addressed by building a basis of vectors via an approximate linear dependency
criterion. Numerical tests on synthetic data demonstrate that the proposed
framework selects always the “optimal” p-norm for the outlier scenario at
hand, outperforming at the same time several non-RL and KBRL schemes.

##  1 Introduction

The least-squares (LS) error/loss (between an observed value and its predicted
one) plays a pivotal role in signal processing, e.g., adaptive filtering [1],
and machine learning [2]. For example, the least-mean squares (LMS) and
recursive (R)LS [1] are two celebrated algorithms in adaptive filtering and
stochastic approximation which are based on the LS-error criterion.
Notwithstanding, LS methods are notoriously sensitive to the presence of
outliers within data [3], where outliers are defined as (sparsely)
contaminating data that do not adhere to a nominal data generation model, and
are often modeled as random variables (RVs) with non-Gaussian heavy tailed
distributions, e.g., $\alpha$-stable ones [4, 5]. To combat the negative
effects of outliers, several non-LS criteria, such as least mean $p$-power
(LMP) [6, 7, 8, 9, 10, 11] and maximum correntropy (MC) [12], have been
studied. This work focuses on the LMP criterion, owing to the well-documented
robustness of LMP against outliers [13], while results on MC will be reported
elsewhere.

This study builds around the classical data-generation model
$y_{n}=\bm{\theta}_{*}^{\intercal}\mathbf{x}_{n}+o_{n}$, where
$n\in\mathbb{N}$ denotes discrete time ($\mathbb{N}$ is the set of all non-
negative integers), $\bm{\theta}_{*}$ is the $L\times 1$ vector whose entries
are the system parameters that need to be identified, $o_{n}$ is the RV which
models outliers/noise, $(\mathbf{x}_{n},y_{n})$ stands for the input-output
pair of available data, where $\mathbf{x}_{n}$ is an $L\times 1$ vector and
$y_{n}$ is real-valued, and $\intercal$ denotes vector/matrix transposition.
For an arbitrarily fixed $\bm{\theta}_{0}$, the LMP algorithm [6] generates
estimates $(\bm{\theta}_{n})_{n\in\mathbb{N}}$ of $\bm{\theta}_{*}$ according
to the following recursion:

$\bm{\theta}_{n+1}\coloneqq\bm{\theta}_{n}+\rho p\lvert
e_{n}\rvert^{p-2}e_{n}\mathbf{x}_{n}\,,$ (1)

where $e_{n}\coloneqq y_{n}-\mathbf{x}_{n}^{\intercal}\bm{\theta}_{n}$, $\rho$
is the learning rate (step size), and $p$ is a fixed user-defined real-valued
number within the interval $[1,2]$ to ensure that the $p$-norm loss $\lvert
y_{n}-\mathbf{x}_{n}^{\intercal}\bm{\theta}\rvert^{p}$ is a convex function of
$\bm{\theta}$ [6]. Notice that if $p=1$ and $2$, then (1) boils down to the
classical sign-LMS and LMS, respectively [1].

Intuition suggests that the choice of $p$ should be based on the distribution
of $o_{n}$ in the data-generation model. For example, if $o_{n}$ obeys a
Gaussian distribution, then $p=2$ should be chosen (recall the maximum-
likelihood criterion). To enhance robustness against outliers, combination of
adaptive filters with different forgetting factors, but with the same fixed
$p$-norm, have been also introduced [9]. However, it seems that an online and
data-driven solution to the problem of dynamically selecting $p$, without any
prior knowledge on the distribution of $o_{n}$, which may change with time, is
yet to be found.

This work offers a solution to the aforementioned open problem via
reinforcement learning (RL) [14]; a machine-learning paradigm where an “agent”
interacts with the surrounding environment to identify iteratively the policy
which minimizes the cost of its “actions.” More specifically, the well-known
policy-iteration (PI) framework [14] of RL is adopted, because of its well-
documented merits (e.g., [15, 16, 17]) over the alternative RL frameworks of
temporal-difference (TD) and Q-learning [14], especially for continuous and
high-dimensional state spaces such as the one considered here. PI comprises
two stages at every iteration $n$: policy evaluation and policy improvement.
At policy evaluation, the current policy is evaluated by a $Q$-function [14],
which represents, loosely speaking, the long-term cost that the agent would
suffer had the current policy been used to determine the next state, whereas
at the policy-improvement stage, the agent uses the $Q$-function values to
update the policy. The underlying state space is considered to be continuous
and high dimensional, due to the nature of the available data
$(\mathbf{x}_{n},y_{n})$, while the action space is considered to be discrete:
an action is a value of $p$ taken from a finite grid of the interval $[1,2]$.

Deep neural networks offer approximating spaces for Q-functions, e.g., [18],
but they may require processing of batch data (even re-training) during
online-mode operation, since they may face test data generated by probability
density functions (PDFs) different from those of the training ones (dynamic
environments). Such batch processing inflicts large computational times and
complexity, discouraging the application of deep neural networks to online
learning where a small complexity footprint is desired.

To meet the desired computational complexity requirements, this study builds
an approximate (A)PI framework for online RL along the lines of kernel-based
(KB)RL [15, 16, 17, 19, 20, 21, 22, 23, 24, 25, 26]. Central to the API design
is the construction of novel Bellman mappings [27, 14]. The proposed Bellman
mappings are defined on a reproducing kernel Hilbert space (RKHS)
$\mathscr{H}$ [28, 29], which serves as the approximating space for the
$Q$-functions. In contrast to the prevailing route in KBRL [15, 16, 17, 19,
20, 21, 22, 23, 24, 25, 26], which views Bellman mappings as contractions in
$\mathscr{L}_{\infty}$-norm Banach spaces (by definition, no inner product
available), this study introduces nonexpansive [30] Bellman operators on
$\mathscr{H}$ to capitalize on the reproducing property of the inner product
of $\mathscr{H}$ [28, 29], and to open the door to powerful Hilbertian tools
[30]. This path offers also flexibility to the user to choose any point from
the potentially non-singleton fixed-point set of the nonexpansive Bellman
mapping, as opposed to the case of a contraction mapping which is known to
have a unique fixed point. A superset of the fixed-point set of one of the
proposed Bellman mappings is designed onto which the proposed API framework is
based. To robustify the policy-improvement stage, the well-known methodology
of rollout [14] is employed. Moreover, to address the issue of the “curse of
dimensionality” that arises naturally in online learning in RKHSs, the
proposed framework is equipped with the approximate linear dependency (ALD)
criterion [31]. Note that [19, 26], being along the lines of TD and
Q-learning, do not include any discussion on Bellman mappings.

Unlike the classical Bellman operators, where information on transition
probabilities of a Markov decision process is needed [14], the proposed
Bellman mappings need neither such information nor any training/offline data
and past policies, but sample and average the state space on-the-fly, at each
$n$, to explore the surrounding environment. This suits the current adaptive-
filtering setting, where the presence of outliers, with a possibly time-
varying PDF, may render the information obtained offline or from past policies
outdated. As such, the proposed Bellman mappings fall closer to [16] than to
studies which use training data collected beforehand (offline), e.g., [32, 17,
33].

It is worth stressing here that the proposed framework, together with its
complementary study [34], appear to be the first attempts to apply RL
arguments to robust adaptive filtering. In contrast to [34], where the state
space is the low-dimensional $\mathbb{R}^{4}$, this study considers the high-
dimensional $\mathbb{R}^{2L+1}$ ($\mathbb{R}$ is the set of all real numbers).
This work constructs a finite-dimensional affine set as a superset of the
fixed-point set of one of the proposed Bellman mappings, as opposed to [34]
where a potentially infinite-dimensional hyperplane is designed. To reduce the
size of the computational footprint of the proposed framework, ALD is used
instead of the random Fourier features (RFF) [35] in [34]. Finally, rollout
[14] is employed for robustification, whereas experience replay [36] is
applied in [34].

Numerical tests on synthetic data demonstrate that the advocated framework
always spots the value of $p$ that leads to “optimal” performance, without any
knowledge on the PDF of the outliers, which is intentionally made to be time
varying. Due to space limitations, any proofs, results on convergence
analysis, and further numerical tests will be reported in the journal version
of this paper.

##  2 Nonexpansive Bellman Mappings on RKHSs

###  2.1 State-Action Space

Following (1), the state space is defined as the following continuous and
high-dimensional
$\mathfrak{S}\coloneqq\\{\mathbf{s}\coloneqq(\mathbf{x},y,\bm{\theta}){\mathbin{}\mid\mathbin{}}\mathbf{x}\in\mathbb{R}^{L},y\in\mathbb{R},\bm{\theta}\in\mathbb{R}^{L}\\}=\mathbb{R}^{2L+1}$,
where $\mathbb{R}$ stands for the set of all real numbers. The action space
$\mathfrak{A}$ is defined as any finite grid of the interval $[1,2]$, so that
an action $a\in\mathfrak{A}$ becomes any value of $p$ taken from that finite
grid. The state-action space is defined as
$\mathfrak{Z}\coloneqq\mathfrak{S}\times\mathfrak{A}$, and its element is
denoted as $\mathbf{z}=(\mathbf{s},a)$.

Along the lines of the general notation in [14], consider now the set of all
mappings
$\mathscr{M}\coloneqq\\{\mu(\cdot)\nonscript\,|\allowbreak\nonscript\,\mathopen{}\mu(\cdot):\mathfrak{S}\to\mathfrak{A}:\mathbf{s}\mapsto\mu(\mathbf{s})\\}$.
In other words, given a $\mu\in\mathscr{M}$, $\mu(\mathbf{s})$ denotes the
action that the “system” may take at state $\mathbf{s}$ to “move to” the state
$\mathbf{s}^{\prime}\in\mathfrak{S}$. The set $\Pi$ of policies is defined as
$\Pi\coloneqq\mathscr{M}^{\mathbb{N}}\coloneqq\\{(\mu_{0},\mu_{1},\ldots,\mu_{n},\ldots)\nonscript\,|\allowbreak\nonscript\,\mathopen{}\mu_{n}\in\mathscr{M},n\in\mathbb{N}\\}$.
A policy will be denoted by $\pi\in\Pi$. Given $\mu\in\mathscr{M}$, the
stationary policy $\pi_{\mu}\in\Pi$ is defined as
$\pi_{\mu}\coloneqq(\mu,\mu,\ldots,\mu,\ldots)$. It is customary for $\mu$ to
denote also the stationary policy $\pi_{\mu}$.

The one-step loss $g:\mathfrak{Z}\to\mathbb{R}:(\mathbf{s},a)\mapsto
g(\mathbf{s},a)$, which quantifies the cost of transition from the current
state $\mathbf{s}$ to the next one $\mathbf{s}^{\prime}$ under action $a$, is
defined for the current setting as:

$\displaystyle
g(\mathbf{s},a)\coloneqq|y-\bm{\theta}^{\intercal}\mathbf{x}|^{a}+|y-\bm{\theta}^{\prime}{}^{\intercal}\mathbf{x}|\,,$
(2)

where $\bm{\theta}^{\prime}\coloneqq\bm{\theta}+\rho a\lvert
y-\bm{\theta}^{\intercal}\mathbf{x}\rvert^{a-2}(y-\bm{\theta}^{\intercal}\mathbf{x})\mathbf{x}$,
according to (1). Recall that action $a$ is defined to take certain values of
$p$ from a user-defined grid in $[1,2]$. The loss in (2) is motivated by
classical adaptive filtering [1]; namely, the first term in (2) resembles the
prior loss, while the second one mimics the posterior loss. Only the prior
loss is chosen to be affected by action $a$ in (2), because $a$ applies to the
current state $\mathbf{s}$ and not $\mathbf{s}^{\prime}$. The long-term loss
$Q:\mathfrak{Z}\to\mathbb{R}:(\mathbf{s},a)\mapsto Q(\mathbf{s},a)$ quantifies
the long-term cost that the agent would suffer had the action $a$ been used to
determine the next state of $\mathbf{s}$.

###  2.2 Novel Bellman Mappings

Central to dynamic programming and RL [14] is the concept of Bellman mappings
which operate on $Q$-functions. Typical definitions are, e.g., [37],
$\forall(\mathbf{s},a)\in\mathfrak{Z}$,

$\displaystyle(T_{\mu}^{\diamond}Q)(\mathbf{s},a)$ $\displaystyle\coloneqq
g(\mathbf{s},a)+\alpha\mathbb{E}_{\mathbf{s}^{\prime}{\mathbin{}\mid\mathbin{}}(\mathbf{s},a)}\\{Q(\mathbf{s}^{\prime},\mu(\mathbf{s}^{\prime}))\\}\,,$
(3a) $\displaystyle(T^{\diamond}Q)(\mathbf{s},a)$ $\displaystyle\coloneqq
g(\mathbf{s},a)+\alpha\mathbb{E}_{\mathbf{s}^{\prime}{\mathbin{}\mid\mathbin{}}(\mathbf{s},a)}\\{\min_{a^{\prime}\in\mathfrak{A}}Q(\mathbf{s}^{\prime},a^{\prime})\\}\,,$
(3b)

where
$\mathbb{E}_{\mathbf{s}^{\prime}{\mathbin{}\mid\mathbin{}}(\mathbf{s},a)}\\{\cdot\\}$
stands for the conditional expectation operator with respect to
$\mathbf{s}^{\prime}$ conditioned on $(\mathbf{s},a)$, and $\alpha\in[0,1)$ is
the discount factor. In the case where $Q$ is considered an element of the
Banach space $\mathscr{L}_{\infty}$ of all (essentially) bounded functions
[38], equipped with the norm $\lVert\>\cdot\>\rVert_{\infty}$, then it can be
shown that the mappings in (3) are contractions [14], and according to the
Banach-Picard theorem [30], they possess unique fixed points
$Q^{\diamond}_{\mu},Q^{\diamond}$, i.e., points which solve the Bellman
equations $T_{\mu}^{\diamond}Q^{\diamond}_{\mu}=Q^{\diamond}_{\mu}$ and
$T^{\diamond}Q^{\diamond}=Q^{\diamond}$, and which characterize “optimal”
long-term losses [14].

Nevertheless, in most cases of practical interest, there is not sufficient
information on the conditional PDF to compute
$\mathbb{E}_{\mathbf{s}^{\prime}{\mathbin{}\mid\mathbin{}}(\mathbf{s},a)}\\{\cdot\\}$.
Motivated by this fact, this study proposes approximations of the Bellman
mappings in (3) by assuming that losses $g$ and $Q$ belong to an RKHS
$\mathscr{H}$, i.e., a Hilbert space with inner product
$\langle\cdot\mathop{}|\mathop{}\cdot\rangle_{\mathscr{H}}$, norm
$\lVert\cdot\rVert_{\mathscr{H}}\coloneqq\langle\cdot\mathop{}|\mathop{}\cdot\rangle_{\mathscr{H}}^{1/2}$,
and a reproducing kernel
$\kappa(\cdot,\cdot):\mathfrak{Z}\times\mathfrak{Z}\to\mathbb{R}$, such that
$\kappa(\mathbf{z},\cdot)\in\mathscr{H}$, $\forall\mathbf{z}\in\mathfrak{Z}$,
and the reproducing property holds true: $Q(\mathbf{z})=\langle
Q\mathop{}|\mathop{}\kappa(\mathbf{z},\cdot)\rangle_{\mathscr{H}}$, $\forall
Q\in\mathscr{H}$, $\forall\mathbf{z}\in\mathfrak{Z}$. Space $\mathscr{H}$ may
be infinite dimensional; e.g., $\kappa(\cdot,\cdot)$ is a Gaussian kernel [28,
29]. For compact notations, let
$\varphi(\mathbf{z})\coloneqq\kappa(\mathbf{z},\cdot)$, and
$Q^{\intercal}Q^{\prime}\coloneqq\langle
Q\mathop{}|\mathop{}Q^{\prime}\rangle_{\mathscr{H}}$.

Hereafter, losses $g,Q$ are assumed to belong to $\mathscr{H}$. The proposed
Bellman mappings $T_{\mu}:\mathscr{H}\to\mathscr{H}:Q\mapsto T_{\mu}Q$ and
$T:\mathscr{H}\to\mathscr{H}:Q\mapsto TQ$ are defined as:

$\displaystyle T_{\mu}Q$ $\displaystyle\coloneqq
g+\alpha\sum\nolimits_{j=1}^{N_{\textnormal{av}}}Q(\mathbf{s}^{\textnormal{av}}_{j},\mu(\mathbf{s}^{\textnormal{av}}_{j}))\cdot\psi_{j}\,,$
(4a) $\displaystyle TQ$ $\displaystyle\coloneqq
g+\alpha\sum\nolimits_{j=1}^{N_{\textnormal{av}}}\inf\nolimits_{a_{j}\in\mathfrak{A}}Q(\mathbf{s}^{\textnormal{av}}_{j},a_{j})\cdot\psi_{j}\,,$
(4b)

where $\\{\psi_{j}\\}_{j=1}^{N_{\textnormal{av}}}\subset\mathscr{H}$, for a
user-defined positive integer $N_{\textnormal{av}}$, and
$\\{\mathbf{s}^{\textnormal{av}}_{j}\\}_{j=1}^{N_{\textnormal{av}}}$ are state
vectors chosen by the user for the summations in (4) to approximate the
conditional expectations in (3). For example,
$\\{\mathbf{s}^{\textnormal{av}}_{j}\\}_{j=1}^{N_{\textnormal{av}}}$ may be
samples drawn from a Gaussian PDF centered at a state of interest (the current
state $\mathbf{s}_{n}$ in Section 4). For notational convenience, let
$\bm{\Psi}\coloneqq[\psi_{1},\ldots,\psi_{N_{\textnormal{av}}}]$, and its
$N_{\textnormal{av}}\times N_{\textnormal{av}}$ kernel matrix
$\mathbf{K}_{\Psi}\coloneqq\bm{\Psi}^{\intercal}\bm{\Psi}$ whose
$(j,j^{\prime})$ entry is equal to
$\langle\psi_{j}\mathop{}|\mathop{}\psi_{j^{\prime}}\rangle_{\mathscr{H}}$.

Consider a subset
$\mathfrak{B}\coloneqq\\{\mathbf{z}^{\textnormal{b}}_{i}\\}^{N_{\textnormal{b}}}_{i=1}\subset\mathfrak{Z}$,
for some $N_{\textnormal{b}}\in\mathbb{N}_{*}$. Define also
$\bm{\Phi}^{\textnormal{b}}\coloneqq[\varphi(\mathbf{z}^{\textnormal{b}}_{1}),\ldots,\varphi(\mathbf{z}^{\textnormal{b}}_{N_{\textnormal{b}}})]$,
$\bm{\Phi}_{\mu}^{\textnormal{av}}\coloneqq[\varphi^{\textnormal{av}}_{\mu,1},\ldots,\varphi^{\textnormal{av}}_{\mu,N_{\textnormal{av}}}]$,
where
$\varphi^{\textnormal{av}}_{\mu,j}\coloneqq\varphi(\mathbf{s}_{j}^{\textnormal{av}},\mu(\mathbf{s}_{j}^{\textnormal{av}}))$,
and let
$\mathbf{K}^{\textnormal{av,b}}_{\mu}\coloneqq{\bm{\Phi}}^{\textnormal{av}}_{\mu}{}^{\intercal}\bm{\Phi}^{\textnormal{b}}$.
Let also the kernel matrices
$\mathbf{K}_{\textnormal{b}}\coloneqq{\bm{\Phi}^{\textnormal{b}}}^{\intercal}\bm{\Phi}^{\textnormal{b}}$,
and
$\mathbf{K}^{\textnormal{av}}_{\mu}\coloneqq\bm{\Phi}_{\mu}^{\textnormal{av}}{}^{\intercal}\bm{\Phi}_{\mu}^{\textnormal{av}}$.
In the case vectors in $\mathfrak{B}$ are linearly independent,
$\mathbf{K}_{\textnormal{b}}$ is positive definite. Moreover, consider an
$N_{\textnormal{b}}\times N_{\textnormal{av}}$ matrix $\bm{\Upsilon}$. Define
then mappings
$T_{\mu}^{\sharp},T^{\sharp}:\mathbb{R}^{N_{\textnormal{b}}}\to\mathbb{R}^{N_{\textnormal{b}}}$
as follows: $\forall\bm{\xi}\in\mathbb{R}^{N_{\textnormal{b}}}$,

$\displaystyle T_{\mu}^{\sharp}\bm{\xi}$
$\displaystyle\coloneqq\bm{\eta}+\alpha\bm{\Upsilon}\mathbf{K}^{\textnormal{av,b}}_{\mu}\bm{\xi}\,,$
(5a) $\displaystyle T^{\sharp}\bm{\xi}$
$\displaystyle\coloneqq\bm{\eta}+\alpha\bm{\Upsilon}\inf\nolimits_{\mu\in\mathscr{M}}\mathbf{K}^{\textnormal{av,b}}_{\mu}\bm{\xi}\,,$
(5b)

It can be verified that the fixed-point set of $T_{\mu}^{\sharp}$ satisfies:

$\displaystyle\operatorname{Fix}T_{\mu}^{\sharp}$
$\displaystyle\coloneqq\\{\bm{\xi}\in\mathbb{R}^{N_{\textnormal{b}}}{\mathbin{}\mid\mathbin{}}T_{\mu}^{\sharp}\bm{\xi}=\bm{\xi}\\}$
$\displaystyle=\\{\bm{\xi}\in\mathbb{R}^{N_{\textnormal{b}}}{\mathbin{}\mid\mathbin{}}(\mathbf{I}_{N_{\textnormal{b}}}-\alpha\bm{\Upsilon}\mathbf{K}^{\textnormal{av,b}}_{\mu})\bm{\xi}=\bm{\eta}\\}\,,$
(6)

where $\mathbf{I}_{N_{\textnormal{b}}}$ stands for the
$N_{\textnormal{b}}\times N_{\textnormal{b}}$ identity matrix. Since
$\operatorname{Fix}T_{\mu}^{\sharp}$ may be empty, 9 of Algorithm 1 defines a
non-empty affine set as the superset of (6). More details on (6) are omitted
due to lack of space. Note here that [34] focuses on (4a) and follows a
different route by defining a potentially infinite-dimensional hyperplane as
the superset of the fixed-point of (4a).

1:Arbitrarily initialize $\mathfrak{B}_{0}\subset\mathfrak{Z}$, with
$N_{\textnormal{b}}[0]=\lvert\mathfrak{B}_{0}\rvert$,
$\bm{\xi}_{0}\in\mathbb{R}^{N_{\textnormal{b}}[0]}$,
$\bm{\Phi}^{\textnormal{b}}_{0}$, $Q_{0}$, $\mu_{0}\in\mathscr{M}$, and a
$\bm{\theta}_{0}\in\mathbb{R}^{L}$.

2:while $n\in\mathbb{N}$ do

3: Data $(\mathbf{x}_{n},y_{n})$ become available. Let
$\mathbf{s}_{n}=(\mathbf{x}_{n},y_{n},\bm{\theta}_{n})$.

4: Define
$\\{\mathbf{s}_{j}^{\textnormal{av}}[n]\\}_{j=1}^{N_{\textnormal{av}}[n]}$.

5: Policy improvement: Update $\mu_{n}(\mathbf{s}_{n})$ by (9).

6: Update $\bm{\theta}_{n+1}$ by (1), where
$p\coloneqq\mu_{n}(\mathbf{s}_{n})$.

7: Update $\mathfrak{B}_{n}$, $\bm{\Phi}^{\textnormal{b}}_{n}$, and identify
$\bm{\Upsilon}_{n}$ and $\mathbf{K}^{\textnormal{av,b}}_{\mu_{n}}$ as in
Section 3.

8: Update $\bm{\eta}_{n}$ by (8b).

9: Compute
$\bm{\xi}_{n+1}\coloneqq\arg\min_{\bm{\xi}}\left\lVert(\mathbf{I}_{N_{\textnormal{b}}[n]}-\alpha\bm{\Upsilon}_{n}\mathbf{K}^{\textnormal{av,b}}_{\mu_{n}})\bm{\xi}-\bm{\eta}_{n}\right\rVert^{2}\,.$

10: Policy evaluation: Update
$Q_{n+1}\coloneqq\bm{\Phi}^{\textnormal{b}}_{n}\bm{\xi}_{n+1}$.

11: Increase $n$ by one, and go to Line 2.

12:end while

Algorithm 1 Approximate policy-iteration framework.

###### Theorem 1

Let $\psi_{j}(\mathbf{z})\geq 0$, $\forall\mathbf{z}\in\mathfrak{Z}$, $\forall
j\in\\{1,\ldots,N_{\textnormal{av}}\\}$.

  1. (i)

If
$\alpha\leq\lVert\mathbf{K}_{\Psi}\rVert^{-1/2}(\sup_{\mu\in\mathscr{M}}\lVert\mathbf{K}^{\textnormal{av}}_{\mu}\rVert)^{-1/2}$,
then $\forall\mu\ \in\mathscr{M}$, the mapping $T_{\mu}$ in (4a) is affine
nonexpansive and $T$ in (4b) is nonexpansive within the Hilbert space
$(\mathscr{H},\langle\cdot\mathop{}|\mathop{}\cdot\rangle_{\mathscr{H}})$.
Norms
$\lVert\mathbf{K}_{\Psi}\rVert,\lVert\mathbf{K}^{\textnormal{av}}_{\mu}\rVert$
are the spectral norms of
$\mathbf{K}_{\Psi},\mathbf{K}^{\textnormal{av}}_{\mu}$.

  2. (ii)

Let
$\alpha\leq\lVert\mathbf{K}_{\Psi}\rVert^{-1/2}(\sup_{\mu\in\mathscr{M}}\lVert\mathbf{K}^{\textnormal{av}}_{\mu}\rVert)^{-1/2}$.
Consider also the case where vectors
$\\{\varphi(\mathbf{z}^{\textnormal{b}}_{i})\\}^{N_{\textnormal{b}}}_{i=1}$
are linearly independent, and
$\\{g,\\{\psi_{j}\\}_{j=1}^{N_{\textnormal{av}}}\\}\subset\mathscr{H}_{\textnormal{b}}\coloneqq\operatorname{span}(\\{\varphi(\mathbf{z}^{\textnormal{b}}_{i})\\}^{N_{\textnormal{b}}}_{i=1})$,
that is, there exist $\bm{\eta}\in\mathbb{R}^{N_{\textnormal{b}}}$ and
$\bm{\Upsilon}\in\mathbb{R}^{N_{\textnormal{b}}\times N_{\textnormal{av}}}$
such that $g=\bm{\Phi}^{\textnormal{b}}\bm{\eta}$ and
$\bm{\Psi}=\bm{\Phi}^{\textnormal{b}}\bm{\Upsilon}$. Then, for any
$Q\in\mathscr{H}_{\textnormal{b}}$,
$T_{\mu}Q=\bm{\Phi}^{\textnormal{b}}T_{\mu}^{\sharp}\bm{\xi}$,
$\forall\mu\in\mathscr{M}$, and
$TQ=\bm{\Phi}^{\textnormal{b}}T^{\sharp}\bm{\xi}$. Moreover, $\forall\mu\
\in\mathscr{M}$, the mapping $T_{\mu}^{\sharp}$ in (5a) is affine nonexpansive
and $T^{\sharp}$ in (5b) is nonexpansive within the Euclidean space
$(\mathbb{R}^{N_{\textnormal{b}}},\langle\cdot\mathop{}|\mathop{}\cdot\rangle_{\mathbf{K}_{\textnormal{b}}})$,
where
$\langle\bm{\xi}\mathop{}|\mathop{}\bm{\xi^{\prime}}\rangle_{\mathbf{K}_{\textnormal{b}}}\coloneqq\bm{\xi}^{\intercal}\mathbf{K}_{\textnormal{b}}\bm{\xi^{\prime}}$.

Nonexpansivity for a mapping $T$ in a (Euclidean) Hilbert space
$(\mathscr{H},\langle\cdot\mathop{}|\mathop{}\cdot\rangle_{\mathscr{H}})$
means $\lVert TQ-TQ^{\prime}\rVert_{\mathscr{H}}\leq\lVert
Q-Q^{\prime}\rVert_{\mathscr{H}}$, $\forall Q,Q^{\prime}\in\mathscr{H}$ [30].
Moreover, mapping $T_{\mu}:\mathscr{H}\to\mathscr{H}$ is affine iff
$T_{\mu}(\lambda Q+(1-\lambda)Q^{\prime})=\lambda TQ+(1-\lambda)TQ^{\prime}$,
$\forall Q,Q^{\prime}\in\mathscr{H}$, $\forall\lambda\in\mathbb{R}$.

Mappings (4) share similarities with the mappings in [15, 16, 20, 21, 22, 25].
However, in [15, 16, 20, 21, 22, 25] as well as in the classical context of
(3), Bellman mappings are viewed as contractions on the Banach space of
(essentially) bounded functions with the $\mathscr{L}_{\infty}$-norm [14],
while no discussion on RKHSs is reported. Recall that, by definition, Banach
spaces are not equipped with inner products. On the other hand, Theorem 1
opens the door to the rich toolbox of nonexpansive mappings in Hilbert spaces
[30], to the reproducing property of the inner product in RKHSs [28, 29], and
also to the familiar Euclidean spaces with their convenient norms by (5) and
Theorem 1(ii).

Due to lack of space, a detailed discussion on ways to choose/design
$\\{\psi_{j}\\}_{j=1}^{N_{\textnormal{av}}}$,
$\\{\mathbf{s}^{\textnormal{av}}_{j}\\}_{j=1}^{N_{\textnormal{av}}}$, and
$\\{\mathbf{z}^{\textnormal{b}}_{i}\\}^{N_{\textnormal{b}}}_{i=1}$ will be
reported in the journal version of the paper. A way to choose
$\\{\psi_{j}\\}_{j=1}^{N_{\textnormal{av}}}$ is reported in [34].

##  3 Approximate Policy-Iteration Framework

Due to the approximate Bellman mappings offered in (4) and (5), the proposed
approximate PI framework for the problem at hand is summarized in Algorithm 1.
The framework operates sequentially, with its iteration index $n$ coinciding
with the time index of the streaming data
$(\mathbf{x}_{n},y_{n})_{n\in\mathbb{N}}$ which are generated according to the
discussion surrounding (1). To this end, the arguments of Section 2 will be
adapted to include hereafter the extra time dimension $n$, which will be
indicated by the super-/sub-scripts $[n]$, $(n)$ or $n$ in notations.

With
$\mathfrak{B}[n-1]=\\{\mathbf{z}^{\textnormal{b}}_{i}[n-1]\\}^{N_{\textnormal{b}}[n-1]}_{i=1}$
available at time $n$, and with
$\mathscr{H}_{\textnormal{b}}[n-1]=\operatorname{span}(\\{\varphi(\mathbf{z}^{\textnormal{b}}_{i}[n-1])\\}^{N_{\textnormal{b}}[n-1]}_{i=1})$,
ALD [31] is used here to let Algorithm 1 decide whether $\mathfrak{B}[n-1]$ is
updated to $\mathfrak{B}[n]$ by including a new basis vector or not. To reduce
as much as possible the unpleasant effects of the “curse of dimensionality,”
$\mathbf{s}_{n}=(\mathbf{x}_{n},y_{n},\bm{\theta}_{n})$ is projected onto the
compact unit ball to define
$\mathbf{s}_{n}^{\|}\coloneqq\mathbf{s}_{n}/\lVert\mathbf{s}_{n}\rVert$, and
the potentially new basis vector
$\varphi(\mathbf{z}^{\textnormal{b}}_{N_{\textnormal{b}}[n]}[n])\coloneqq\kappa(\mathbf{z}^{\textnormal{b}}_{N_{\textnormal{b}}[n]}[n],\cdot)$
of $\mathscr{H}_{\textnormal{b}}[n]$, with
$\mathbf{z}^{\textnormal{b}}_{N_{\textnormal{b}}[n]}[n]\coloneqq(\mathbf{s}_{n}^{\|},\mu_{n}(\mathbf{s}_{n}))$,
is tested against $\mathscr{H}_{\textnormal{b}}[n-1]$ according to ALD [31].
Note also that ALD ensures that vectors in $\mathscr{H}_{\textnormal{b}}[n]$
are linearly independent, as needed in Theorem 1(ii). ALD is also utilized
upon the arrival of
$\\{\mathbf{s}^{\textnormal{av}}_{j}[n]\\}^{N_{\textnormal{av}}[n]}_{j=1}$ to
check whether any of these vectors enter the basis or not.

To utilize (5a), an $N_{\textnormal{b}}[n]\times 1$ vector $\bm{\eta}_{n}$
needs to be identified. According to Theorem 1(ii), such an $\bm{\eta}_{n}$
identifies a vector
$\tilde{g}_{n}\coloneqq\bm{\Phi}^{\textnormal{b}}_{n}\bm{\eta}_{n}$ in
$\mathscr{H}_{\textnormal{b}}[n]$. More specifically, $\bm{\eta}_{n}$ is
chosen such that $[\tilde{g}_{n}(\mathbf{z}_{n})-g(\mathbf{z}_{n})]^{2}$ is
minimized, where $\mathbf{z}_{n}=(\mathbf{s}_{n},\mu_{n}(\mathbf{s}_{n}))$,
and the value $g(\mathbf{z}_{n})$ is provided by (2). Here is where the
reproducing property of the inner product in $\mathscr{H}$ becomes handy:
$\tilde{g}_{n}(\mathbf{z}_{n})=\langle\tilde{g}_{n}\mathop{}|\mathop{}\varphi(\mathbf{z}_{n})\rangle_{\mathscr{H}}=\langle\bm{\Phi}^{\textnormal{b}}_{n}\bm{\eta}_{n}\mathop{}|\mathop{}\varphi(\mathbf{z}_{n})\rangle_{\mathscr{H}}=\mathbf{k}^{\intercal}(\mathbf{z}_{n})\bm{\eta}_{n}$,
where
$\mathbf{k}(\mathbf{z}_{n})\coloneqq\bm{\Phi}^{\textnormal{b}}_{n}{}^{\intercal}\varphi(\mathbf{z}_{n})=[\kappa(\mathbf{z}^{\textnormal{b}}_{1}[n],\mathbf{z}_{n}),\ldots,\kappa(\mathbf{z}^{\textnormal{b}}_{N_{\textnormal{b}}[n]}[n],\mathbf{z}_{n})]^{\intercal}$.
Hence, $\bm{\eta}_{n}$ can be chosen as:

$\displaystyle\bm{\eta}_{n}\in{}$
$\displaystyle{}\arg\min\nolimits_{\bm{\eta}\in\mathbb{R}^{N_{\textnormal{b}}[n]}}[\mathbf{k}^{\intercal}(\mathbf{z}_{n})\bm{\eta}-g(\mathbf{z}_{n})]^{2}$
$\displaystyle{}=\\{\bm{\eta}\in\mathbb{R}^{N_{\textnormal{b}}[n]}{\mathbin{}\mid\mathbin{}}\mathbf{k}^{\intercal}(\mathbf{z}_{n})\bm{\eta}=g(\mathbf{z}_{n})\\}\,.$
(7)

To this end, the following options are provided to choose from:

$\displaystyle\bm{\eta}_{n}\coloneqq$
$\displaystyle\frac{g(\mathbf{z}_{n})}{\mathbf{k}^{\intercal}(\mathbf{z}_{n})\mathbf{K}_{\textnormal{b}}[n]\mathbf{k}(\mathbf{z}_{n})}\mathbf{k}(\mathbf{z}_{n})\,,$
(8a) $\displaystyle\bm{\eta}_{n}\coloneqq$
$\displaystyle\bm{\eta}_{n-1}+\frac{g(\mathbf{z}_{n})-\bm{\eta}_{n-1}^{\intercal}\mathbf{K}_{\textnormal{b}}[n]\mathbf{k}(\mathbf{z}_{n})}{\mathbf{k}^{\intercal}(\mathbf{z}_{n})\mathbf{K}_{\textnormal{b}}[n]\mathbf{k}(\mathbf{z}_{n})}\mathbf{k}(\mathbf{z}_{n})\,,$
(8b)

where (8a) offers the minimum-norm solution of (7), while (8b) computes the
(metric) projection [30] of $\bm{\eta}_{n-1}$ onto the hyperplane in (7),
after zero-padding the $N_{\textnormal{b}}[n-1]\times 1$ vector
$\bm{\eta}_{n-1}$ to meet the length $N_{\textnormal{b}}[n]$, in case
$N_{\textnormal{b}}[n]\geq N_{\textnormal{b}}[n-1]+1$. Following Theorem
1(ii), the inner product
$\langle\cdot\mathop{}|\mathop{}\cdot\rangle_{\mathbf{K}_{\textnormal{b}}[n]}$
was utilized in (8b).

To identify $\bm{\Upsilon}_{n}$, needed for the computations in 9 of Algorithm
1, Theorem 1(ii) is followed where $\bm{\Upsilon}_{n}$ is computed via
$\bm{\Psi}_{n}$. Although $\bm{\Psi}_{n}$ can be chosen in many ways such that
$\\{\psi_{j}[n]\\}_{j=1}^{N_{\textnormal{av}}[n]}\subset\mathfrak{B}[n]$,
$\\{\psi_{j}[n]\\}_{j=1}^{N_{\textnormal{av}}[n]}$ are selected in Section 4
simply by taking the first $N_{\textnormal{av}}[n]$ vectors from
$\mathfrak{B}[n]$, under the assumption that $N_{\textnormal{av}}[n]\leq
N_{\textnormal{b}}[n]$. Hence, according to Theorem 1(ii),
$\bm{\Upsilon}_{n}=[\mathbf{I}_{N_{\textnormal{av}}[n]},\mathbf{0}_{N_{\textnormal{av}}[n]\times(N_{\textnormal{b}}[n]-N_{\textnormal{av}}[n])}]^{\intercal}$
in Section 4. With regards to
$\varphi^{\textnormal{av}}_{\mu,j}[n]\coloneqq\varphi(\mathbf{s}_{j}^{\textnormal{av}}[n],\mu_{n}(\mathbf{s}_{j}^{\textnormal{av}}[n]))$,
needed in the computation of $\mathbf{K}^{\textnormal{av,b}}_{\mu_{n}}$, since
$\mathbf{s}_{j}^{\textnormal{av}}[n]$ is drawn from a Gaussian PDF centered at
$\mathbf{s}_{n}$, $\mu_{n}(\mathbf{s}_{j}^{\textnormal{av}}[n])$ is set equal
to $\mu_{n}(\mathbf{s}_{n})$ in Section 4.

Typically, the greedy rule
$\mu_{n}\left(\mathbf{s}_{n}\right)=\arg\min_{a\in\mathfrak{A}}Q_{n}(\mathbf{s}_{n},a)$
is used for policy improvement [14, 32, 17]. However, it has been observed
that this greedy rule may lead into instabilities of the RL agent’s behavior
and hinder its cognition about the surrounding environment [14]. The
$\varepsilon$-greedy strategy [14], a variant of the greedy one, has not
performed well in the numerical tests of Section 4. To address this potential
drawback of greedy strategies, the popular rollout methodology [14] is
employed in Algorithm 1. With the user-defined $M\in\mathbb{N}_{*}$ denoting
the number of steps in rollout, the following rule is employed for policy
improvement:

$\displaystyle\mu_{n}(\mathbf{s}_{n})=\arg\min_{a\in\mathfrak{A}}\Big{[}$
$\displaystyle
g(\mathbf{s}_{n},a)+\sum\nolimits_{m=1}^{M-1}\alpha^{m}g(\mathbf{s}_{n+m}^{\triangleright},\mu^{\triangleright}(\mathbf{s}_{n+m}^{\triangleright}))$
$\displaystyle+\alpha^{M}Q_{n}(\mathbf{s}_{n+M}^{\triangleright},\mu^{\triangleright}(\mathbf{s}_{n+M}^{\triangleright}))\Big{]}\,,$
(9)

where $\mu^{\triangleright}\in\mathscr{M}$ is called the heuristic
(stationary) policy, which could be chosen to be random or $\mu_{n-1}$, and
$\\{\mathbf{s}_{n+m}^{\triangleright}\\}_{m=1}^{M}$, denote possible successor
states of $\mathbf{s}_{n}$ under action $a$ and policy $\mu^{\triangleright}$.
A large $M$ (long trajectory) may inflict long computational times and
diminish the effect of $Q_{n}$ on policy improvement through $\alpha^{M}$ in
(9); usually, $\alpha\leq 1$ [14]. Note here that experience replay [36] and
not rollout was employed in [34].

##  4 Numerical Tests

The proposed PI framework is compared numerically against  (i) (1), where $p$
is fixed throughout iterations with
$p\in\mathfrak{A}\coloneqq\\{1,1.25,1.5,1.75,2\\}$,  (ii) CAC-RL$p$ [9], which
uses a combination of $p$-norm adaptive filters with different forgetting
factors, and  (iii) [17], which is based on the kernel-based least-squares
temporal-difference methodology with no rollout employed.  Performance is
measured by the normalized deviation
$\lVert\bm{\theta}_{*}-\bm{\theta}_{n}\rVert/\lVert\bm{\theta}_{*}\rVert$ vs.
time index $n$. Multiple ($100$) independent tests were performed, with their
uniformly averaged results reported in Figures 1 and 2. The software code was
written in Julia [39], with $\alpha$-stable outliers generated by [5].

The length $L$ of the system $\bm{\theta}_{*}$ is set equal to $10$ in the
data-generation model of Section 1. The entries of $\bm{\theta}_{*}$ and
$\mathbf{x}_{n}$, $\forall n$, are generated by independent and identically
distributed (IID) normal RVs. In the data-generation model of Section 1, RV
$o_{n}$ describes both outlier and noise scenarios. In the case where no
outliers appear in the data-generation model, $o_{n}$ obeys the Gaussian PDF,
and whenever outliers appear, $o_{n}$ follows the $\alpha$-stable PDF [5]. In
all figures, no outliers appear for time instances $n<2\times 10^{4}$ and only
Gaussian noise, with $\textnormal{SNR}=20\textnormal{dB}$, corrupts the data.
Outliers appear for $n\geq 2\times 10^{4}$, following the $\alpha$-stable PDF
[4]. Note that when the “stability” parameter of the $\alpha$-stable PDF takes
the values of $2$, then the PDF boils down to the Gaussian one. Two types of
$\alpha$-stable outliers are examined:  (i) “Gaussian-like” ones, with
parameters $\textnormal{"stability"}=1.95$, $\textnormal{"skewness"}=0.5$,
$\textnormal{"location"}=0.5$, $\textnormal{"scale"}=10^{-2}$, which make the
tails of the PDF slightly heavier than those of the Gaussian PDF [5]; and
(ii) “Cauchy-like” ones, with parameters $\textnormal{"stability"}=1$,
$\textnormal{"skewness"}=0.5$, $\textnormal{"location"}=0.5$ and
$\textnormal{"scale"}=1$, which make the tails of the PDF rather heavy [5].

Number $N_{\textnormal{av}}[n]=10$, $\forall n$, $M=2$ in (9), while (8a) is
employed. The learning rate $\rho$ of LMP in (1) is set equal to $10^{-3}$. A
homogeneous polynomial kernel of degree $2$ is used for both Algorithm 1 and
KLSPI [17]. Note that this kernel satisfies the condition stated in the first
line of Theorem 1. Results on other kernels will be reported elsewhere. The
ALD criterion is employed with parameter
$\delta_{\textnormal{ALD}}=\sin(40\pi/180)$. The parameters of CAC-RL$p$ are
set as $p_{\textnormal{RLP}}=1.4$, $\beta_{\textnormal{RLP}}=0.9$ and
$\lambda_{\textnormal{RLP}}=0.99$.

All tests show that Algorithm 1 succeeds always in identifying the “optimal”
$p$-norm, regardless of the PDF of the outliers. Even if outliers are absent
from the data-generation model and only Gaussian noise corrupts the data,
Algorithm 1 chooses the $2$-norm, which is well-known to be optimal for
Gaussian noise. In contrast, KLSPI identifies fast the $2$-norm in the case
where only Gaussian noise appears, but performs poorly when Cauchy-like
outliers appear. CAC-RL$p$ [9] converges fast in cases where only Gaussian
noise appears, with sub-optimal performance in terms of the normalized
deviation, but diverges in the case where Cauchy-like outliers corrupt the
data. Due to space limitations, further numerical tests will be reported in
the journal version of this manuscript.

(a)

(b)

Fig. 1: Outliers appear only at $n\geq 2\times 10^{4}$. Curve markers:
Algorithm 1: , KLSPI [17]: , CAC-RL$p$ [9]: , optimal $p$-norm criterion,
$2$-norm:  for Figure 1(a) and $1$-norm:  for Figure 1(b).

(a)

(b)

Fig. 2: System $\bm{\theta}_{*}$ changes randomly at $n=2\times 10^{4}$, and
outliers appear only at $n\geq 2\times 10^{4}$. Curve markers follow those of
Figure 1.

## References

  * [1] A. H. Sayed, Adaptive Filters.  Wiley, 2011.
  * [2] S. Theodoridis, Machine Learning—A Bayesian and Optimization Perspective.  Elsevier, 2nd ed., 2020.
  * [3] P. J. Rousseeuw and A. Leroy, Robust Regression and Outlier Detection.  Wiley, 1987.
  * [4] M. Shao and C. L. Nikias, “Signal processing with fractional lower order moments: Stable processes and their applications,” Proceedings of the IEEE, vol. 81, no. 7, pp. 986–1010, 1993.
  * [5] J. M. Miotto, “Pylevy.” https://github.com/josemiotto/pylevy, 2020.
  * [6] S.-C. Pei and C.-C. Tseng, “Least mean p-power error criterion for adaptive FIR filter,” IEEE Journal on Selected Areas in Communications, vol. 12, no. 9, pp. 1540–1547, 1994.
  * [7] Y. Xiao, Y. Tadokoro, and K. Shida, “Adaptive algorithm based on least mean p-power error criterion for Fourier analysis in additive noise,” IEEE Transactions on Signal Processing, vol. 47, no. 4, pp. 1172–1181, 1999.
  * [8] E. E. Kuruoğlu, “Nonlinear least $\ell_{p}$-norm filters for nonlinear autoregressive $\alpha$-stable processes,” Digital Signal Processing, vol. 12, no. 1, pp. 119–142, 2002.
  * [9] A. Navia-Vazquez and J. Arenas-Garcia, “Combination of recursive least p-norm algorithms for robust adaptive filtering in alpha-stable noise,” IEEE Transactions on Signal Processing, vol. 60, no. 3, pp. 1478–1482, 2012.
  * [10] B. Chen, L. Xing, Z. Wu, J. Liang, J. C. Príncipe, and N. Zheng, “Smoothed least mean p-power error criterion for adaptive filtering,” Digital Signal Processing, vol. 40, pp. 154–163, May 2015.
  * [11] K. Slavakis and M. Yukawa, “Outlier-robust kernel hierarchical-optimization RLS on a budget with affine constraints,” in Proc. IEEE ICASSP, pp. 5335–5339, 2021.
  * [12] A. Singh and J. C. Príncipe, “Using correntropy as a cost function in linear adaptive filters,” in Proc. International Joint Conference on Neural Networks, pp. 2950–2955, 2009.
  * [13] C. Gentile, “The robustness of the p-norm algorithms,” Machine Learning, vol. 53, pp. 265–299, 2003.
  * [14] D. Bertsekas, Reinforcement Learning and Optimal Control.  Athena Scientific, 2019.
  * [15] D. Ormoneit and S. Sen, “Kernel-based reinforcement learning,” Machine Learning, vol. 49, pp. 161–178, 2002.
  * [16] D. Ormoneit and P. Glynn, “Kernel-based reinforcement learning in average-cost problems,” IEEE Transactions on Automatic Control, vol. 47, pp. 1624–1636, Oct. 2002.
  * [17] X. Xu, D. Hu, and X. Lu, “Kernel-based least squares policy iteration for reinforcement learning,” IEEE Transactions on Neural Networks, vol. 18, no. 4, pp. 973–992, 2007.
  * [18] H. Van Hasselt, A. Guez, and D. Silver, “Deep reinforcement learning with double Q-learning,” in Proc. AAAI conference on Artificial Intelligence, vol. 30, 2016.
  * [19] J. Bae, L. S. Giraldo, P. Chhatbar, J. Francis, J. Sanchez, and J. Príncipe, “Stochastic kernel temporal difference for reinforcement learning,” in Proc. IEEE MLSP, pp. 1–6, 2011.
  * [20] A. Barreto, D. Precup, and J. Pineau, “Reinforcement learning using kernel-based stochastic factorization,” in Proc. NIPS, vol. 24, 2011.
  * [21] A. Barreto, D. Precup, and J. Pineau, “On-line reinforcement learning using incremental kernel-based stochastic factorization,” in Proc. NIPS, vol. 25, 2012.
  * [22] B. Kveton and G. Theocharous, “Structured kernel-based reinforcement learning,” in Proc. AAAI Conference on Artificial Intelligence, vol. 27, pp. 569–575, June 2013.
  * [23] W. Sun and J. A. Bagnell, “Online Bellman residual and temporal difference algorithms with predictive error guarantees,” in Proc. International Joint Conference on Artificial Intelligence, pp. 4213–4217, 2016.
  * [24] A.-M. Farahmand, M. Ghavamzadeh, C. Szepesvári, and S. Mannor, “Regularized policy iteration with nonparametric function spaces,” J. Machine Learning Research, vol. 17, no. 1, pp. 4809–4874, 2016.
  * [25] B. Kveton and G. Theocharous, “Kernel-based reinforcement learning on representative states,” in Proc. AAAI Conference on Artificial Intelligence, vol. 26, pp. 977–983, Sept. 2021.
  * [26] Y. Wang and J. C. Príncipe, “Reinforcement learning in reproducing kernel Hilbert spaces,” IEEE Signal Processing Magazine, vol. 38, no. 4, pp. 34–45, 2021.
  * [27] R. E. Bellman, Dynamic Programming.  Dover Publications, 2003.
  * [28] N. Aronszajn, “Theory of reproducing kernels,” Transactions of the American Mathematical Society, vol. 68, pp. 337–404, 1950.
  * [29] B. Schölkopf and A. J. Smola, Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond.  Adaptive computation and machine learning, MIT Press, 2002.
  * [30] H. H. Bauschke and P. L. Combettes, Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  New York: Springer, 2011.
  * [31] Y. Engel, S. Mannor, and R. Meir, “The kernel recursive least-squares algorithm,” IEEE Transactions on Signal Processing, vol. 52, no. 8, pp. 2275–2285, 2004.
  * [32] M. G. Lagoudakis and R. Parr, “Least-squares policy iteration,” J. Mach. Learn. Res., vol. 4, pp. 1107–1149, Dec. 2003.
  * [33] K. Panaganti, Z. Xu, D. Kalathil, and M. Ghavamzadeh, “Robust reinforcement learning using offline data,” arXiv, 2022.  abs/2208.05129.
  * [34] Y. Akiyama, M. Vu, and K. Slavakis, “Online and lightweight kernel-based approximate policy iteration for dynamic p-norm linear adaptive filtering.” Submitted for publication to arXiv, Oct. 2022.
  * [35] A. Rahimi and B. Recht, “Random features for large-scale kernel machines,” in Proc. NIPS, vol. 20, 2007.
  * [36] T. Schaul, J. Quan, I. Antonoglou, and D. Silver, “Prioritized experience replay,” in Proc. International Conference on Learning Representations, 2016.
  * [37] M. G. Bellemare, G. Ostrovski, A. Guez, P. Thomas, and R. Munos, “Increasing the action gap: New operators for reinforcement learning,” in Proc. AAAI Conference on Artificial Intelligence, vol. 30, 2016.
  * [38] R. G. Bartle, The Elements of Integration and Lebesgue Measure.  John Wiley & Sons, 1995.
  * [39] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, “Julia: A fresh approach to numerical computing,” SIAM Review, vol. 59, no. 1, pp. 65–98, 2017.