Online Estimation and Inference for Robust Policy Evaluation in Reinforcement Learning

Reinforcement learning has offered immense success and remarkable breakthroughs in a variety of application domains, including autonomous driving, precision medicine, recommendation systems, and robotics (to name a few, e.g., Murphy, 2003; Kormushev et al., 2013; Mnih et al., 2015; Shi et al., 2018). From recommendation systems to mobile health (mHealth) intervention, reinforcement learning can be used to adaptively make personalized recommendations and optimize intervention strategies learned from retrospective behavioral and physiology data. While the achievements of reinforcement learning algorithms in applications are undisputed, the reproducibility of its results and reliability is still in many ways nascent.

Those recommendation and health applications enjoy great flexibility and affordability due to the development of reinforcement algorithms, despite calling for critical needs for a reliable and trustworthy uncertainty quantification for such implementation. The reliability of such implementations sometimes plays a life-threatening role in emerging applications. For example, in autonomous driving, it is critical to avoid deadly explorations based upon some uncertainty measures in the trial-and-error learning procedure. This substance also extends to other applications including precision medicine and autonomous robotics. From the statistical perspective, it is important to quantify the uncertainty of a point estimate with complementary hypothesis testing to reveal or justify the reliability of the learning procedure.

Policy evaluation plays a cornerstone role in typical reinforcement learning (RL) algorithms such as Temporal Difference (TD) learning. As one of the most commonly adopted algorithms for policy evaluation in RL, TD learning provides an estimator of the value function iteratively with regard to a given policy based on samples from a Markov chain. In large-scale RL tasks where the state space is infinitely expansive, a typical procedure to provide a scalable yet efficient estimation of the value function is via linear function approximation. This procedure can be formulated in a linear stochastic approximation problem (Sutton, 1988; Tsitsiklis and Van Roy, 1997; Sutton et al., 2009; Ramprasad et al., 2022), which is designed to sequentially solve a deterministic equation $A\theta=b$ by a matrix-vector pair of a sequence of unbiased random observations of $(A_{t},b_{t})$ governed by an ergodic Markov chain.

The earliest and most prototypical stochastic approximation algorithm is the Robbins-Monro algorithm introduced by Robbins and Monro (1951) for solving a root-finding problem, where the function is represented as an expected value, e.g., ${\mathbb{E}}[f(\theta)]=0$. The algorithm has generated profound interest in the field of stochastic optimization and machine learning to minimize a loss function using random samples. When referring to an optimization problem, its first-order condition can be represented as ${\mathbb{E}}[f(\theta)]=0$, and the corresponding Robbins-Monro algorithm is often referred to as first-order methods, or more widely known as stochastic gradient descent (SGD) in machine learning literature. It is well established in the literature that its averaged version (Ruppert, 1988; Polyak and Juditsky, 1992), as an online first-order method, achieves optimal statistical efficiency when estimating the model parameters in statistical models, which apparently kills the interest in developing second-order methods that use additional information to help the convergence. That being said, it is often observed in practice that first-order algorithms entail significant accuracy loss on non-asymptotic convergence as well as severe instability in the choice of hyperparameters, specifically, the stepsizes (a.k.a. learning rate). In addition, the stepsize tuning further complicates the quantification of uncertainty associated with the algorithm output. Despite the known drawbacks above, first-order stochastic methods are historically favored in machine learning tasks due to their computational efficiency while they primarily focus on estimation. On the other hand, when the emphasis of the task lies on statistical inference, certain computation of the second-order information is generally inevitable during the inferential procedure, which shakes the supremacy of first-order methods over second-order algorithms.

In light of that, we propose a second-order online algorithm which utilizes second-order information to perform the policy evaluation sequentially. Meanwhile, our algorithm can be used for conducting statistical inference in an online fashion, allowing for the characterization of uncertainty in the estimation of the value function. Such a procedure generates no extra per-unit computation and storage cost beyond $O(d^{2})$, which is at least the same as typical first-order stochastic methods featuring statistical inference (see, e.g., Chen et al., 2020 for SGD and Ramprasad et al., 2022 for TD). More importantly, we show theoretically that the proposed algorithm converges faster in terms of the remainder term compared with first-order stochastic approximation methods, and revealed significant discrepancies in numerical experiments. In addition, the proposed algorithm is free from tuning stepsizes, which has been well-established as a substantial criticism of first-order algorithms.

Another challenge to the reliability of reinforcement learning algorithms lies in the modeling assumptions. Most algorithms in RL have been in the “optimism in the face of uncertainty” paradigm where such procedures are vulnerable to manipulation (see some earlier exploration in e.g., Everitt et al., 2017; Wang et al., 2020). In practice, it is often unrealistic to believe that rewards on the entire trajectory follow exactly the same underlying model. Indeed, non-standard behavior of the rewards happens from time to time in practice. The model-misspecification and presence of outliers are indeed very common in an RL environment, especially that with a large time horizon $T$. It is of substantial interest to design a robust policy evaluation procedure. In pursuit of this, our proposed algorithm uses a smoothed Huber loss to replace the least-squares loss function used in classical TD learning, which is tailored to handle both outliers and heavy-tailed rewards. To model outlier observations of rewards in reinforcement learning, we bring the static $\alpha$-contamination model (Huber, 1992) to an online setting with dependent samples. In a static offline robust estimation problem, one aims to learn the distribution of interest, where a sample of size $n$ is drawn i.i.d. from a mixture distribution $(1-\alpha_{n})P+\alpha_{n}Q$, and $Q$ denotes an arbitrary outlier distribution. We adopt robust estimation in an online environment, where the observations are no longer independent, and the occurrence time of outliers is unknown. In contrast to the offline setting, future observations cannot be used in earlier periods in an online setting. Therefore, in the earlier periods, there is very limited information to help determine whether an observation is an outlier. In addition to this discrepancy between the online decision process and offline estimation, we further allow the outlier reward models to be potentially different for different time $t$ instead of being from a fixed distribution, and such rewards may be arbitrary and even adversarially adaptive to historical information. Besides the outlier model, our model also incorporates rewards with heavy-tailed distributions. This substantially relaxes the boundness condition that is commonly assumed in policy evaluation literature.

We summarize the challenges and contributions of this paper in the following facets.

• •

  We propose an online policy evaluation method with dependent samples and simultaneously conduct statistical inference for the model parameters in a fully online fashion. Furthermore, we build a Bahadur representation of the proposed estimator, which includes the main term corresponding to the asymptotic normal distribution and a higher-order remainder term. To the best of our knowledge, there is no literature establishing the Bahadur representation for online policy evaluation. Moreover, it shows that our algorithm matches the offline oracle and converges strictly faster to the asymptotic distribution than that of a prototypical first-order stochastic method such as TD learning.

• •

  Compared to existing reinforcement learning literature, our proposed algorithm features an online generalization of the $\alpha_{n}$-contamination model where the rewards contain outliers or arbitrary corruptions. Our proposed algorithm is robust to adversarial corruptions which can be adaptive to the trajectory, as well as heavy-tailed distribution of the rewards. Due to the existence of outliers, we use a smooth Huber loss where the thresholding parameter is carefully specified to change over time to accommodate the online streaming data. From a theoretical standpoint, a robust policy evaluation procedure forces the update step from $\theta_{t}$ to $\theta_{t+1}$ to be a non-linear function of $\theta_{t}$, which brings in additional technical challenges compared to the analysis of classical TD learning algorithms (see e.g., Ramprasad et al., 2022) based on linear stochastic approximation.

• •

  Our proposed algorithm is based on a dedicated averaged version of the second-order method, where in each iteration a surrogate Hessian is obtained and used in the update step. This second-order information enables the proposed algorithm to be free from stepsize tuning while still ensuring efficient implementation. Furthermore, our proposed algorithm stands out distinctly from conventional first-order stochastic approximation approaches which fall short of attaining the optimal offline remainder rate. On the other hand, while deterministic second-order methods do excel in offline scenarios, they lack the online adaptability crucial for real-time applications.

We first review the Least-squares temporal difference methods in RL. Consider a $4$-tuple $({\mathcal{S}},{\mathcal{U}},{\mathcal{P}},{\mathcal{R}})$. Here ${\mathcal{S}}=\{1,2,...,N\}$ is the global finite state space, ${\mathcal{U}}$ is the set of control (action), ${\mathcal{P}}$ is the transition kernel, and ${\mathcal{R}}$ is the reward function. One of the core steps in RL is to estimate the accumulative reward $J^{*}$ (which is also called the state value function) for a given policy:

where $\gamma\in(0,1]$ is a given discount factor and $s\in{\mathcal{S}}$ is any state. Here $\{s_{k}\}$ denote the environment states which are usually modeled by a Markov chain. In real-world RL applications, the state space is often very large such that one cannot directly compute value estimates for every state in the state space. A common approach in the modern RL is to approximate the value function $J^{*}(\cdot)$, i.e., let

be a linear approximation of $J^{*}(\cdot)$, where $\boldsymbol{\theta}=(\theta_{1},...,\theta_{d})^{\top}\in\mathbb{R}^{d}$ contains the model parameters, and $\boldsymbol{\phi}(s)=(\phi_{1}(s),...,\phi_{d}(s))^{\top}\in\mathbb{R}^{d}$ is a set of feature vectors that corresponds to the state $s\in{\mathcal{S}}$. Here we write $\boldsymbol{\phi}_{l}=(\phi_{l}(1),...,\phi_{l}(N))^{\top}$, $1\leq l\leq d$ are linearly independent vectors in ${\mathbb{R}}^{N}$ and $d\ll N$. That is, we use a low dimensional linear approximation (with the basis vectors $\{\boldsymbol{\phi}_{1},...,\boldsymbol{\phi}_{d}\}$) for $J^{*}(\cdot)$. Let the matrix $\boldsymbol{\Phi}=(\boldsymbol{\phi}_{1},...,\boldsymbol{\phi}_{d})$ and then $\widetilde{J}=\boldsymbol{\Phi}\boldsymbol{\theta}$.

The state value function $J^{*}$ satisfies the Bellman equation

where $s^{\prime}$ is the next state transferred from $s$. Let ${\mathcal{P}}=(p_{ss^{\prime}})_{N\times N}$ be the probability transition matrix and define the Bellman operator $\mathcal{T}$ by $\mathcal{T}(Q)={\mathcal{R}}+\gamma{\mathcal{P}}Q$. The state function $J^{*}$ is the unique fixed point of the Bellman operator $\mathcal{T}$. When $J^{*}$ is replaced by $\widetilde{J}=\boldsymbol{\Phi}\boldsymbol{\theta}$, the Bellman equation may not hold.

The classical temporal difference attempts to find $\boldsymbol{\theta}$ such that $\widetilde{J}(s,\boldsymbol{\theta})$ well approximates the value function $J^{*}(s)$. Particularly, the TD algorithm updates

where $\eta_{t}$ is the step size which often requires careful tuning. We next illustrate the key point that (2) leads to a good estimator such that $\widetilde{J}(s,\boldsymbol{\theta}_{t})$ is close to $J^{*}(s)$. It can be shown that $\boldsymbol{\theta}_{t}$ converges to an unknown population parameter $\boldsymbol{\theta}^{*}$ that minimizes the expected squared difference of the Bellman equation,

It is easy to see that such $\boldsymbol{\theta}^{*}$ exists and satisfies the following first-order condition,

Under the condition that $\textbf{H}:={\mathbb{E}}\boldsymbol{\phi}(s)(\boldsymbol{\phi}^{{\top}}(s)-\gamma\boldsymbol{\phi}^{{\top}}(s^{\prime}))$ is positive definite in the sense that $x^{{\top}}\textbf{H}x>0$ for all $x\in{\mathbb{R}}^{d}$; see Tsitsiklis and Van Roy (1997), the solution of (4) writes

The estimation equation (2) is a first-order stochastic algorithm that converges to the stochastic root-finding problem (4) using a sequence of observations $\{s_{t},r_{t}\}_{t\geq 1}$. In this paper, we refer to it as the Least-squares temporal difference estimator. See Kolter and Ng (2009) for details on the properties of the Least-squares TD estimator.

Least-squares-based methods are oftentimes criticized due to their sensitivity to outliers in data. When there may exist outliers in some observations of the reward ${\mathcal{R}}(s)$, it is a natural call for interest to design a robust estimator of $\boldsymbol{\theta}^{*}$. Following the widely-known classical literature (Huber, 1964; Charbonnier et al., 1994, 1997; Hastie et al., 2009) , we replace the square loss $|\cdot|^{2}$ by a smoothed (Pseudo) Huber loss $f_{\tau}(x)=\tau^{2}(\sqrt{1+(x/\tau)^{2}}-1)$, parametrized by a thresholding parameter $\tau$. We define a similar fixed point equation with the smoothed Huber loss by

In this section, when we motivate the algorithm, we assume that the fixed point $\boldsymbol{\theta}^{*}_{\tau}$ exists ¹¹1Here $\boldsymbol{\theta}^{*}_{\tau}$ is used for the motivation only. Our algorithm and theory do not depend on the existence of $\boldsymbol{\theta}^{*}_{\tau}$. . As the thresholding parameter $\tau$ goes to infinity, the objective equation in (5) is close to the least-squares loss in (3), and $\boldsymbol{\theta}_{\tau}^{*}$ should be close to $\boldsymbol{\theta}^{*}$. When $\tau$ tends to $0$, the problem (5) becomes $\boldsymbol{\theta}^{*}_{0}=\mathop{\rm{argmin}}_{\boldsymbol{u}\in{\mathbb{R}}^{d}}{\mathbb{E}}|\boldsymbol{\phi}^{{\top}}(s)\boldsymbol{u}-({\mathcal{R}}(s)+\gamma\boldsymbol{\phi}^{{\top}}(s^{\prime})\boldsymbol{\theta}^{*}_{0})|,$ with a nonsmooth least absolute deviation (LAD) loss, which is out of the scope of this paper. In this paper, we carefully specify $\tau$ to balance the statistical efficiency and the effect of potential outliers in an online fashion (see Theorem 1).

By the first-order condition of (5), we obtain a similar estimation function as (4),

where the function $g_{\tau}(x)=f^{\prime}_{\tau}(x)$ is the differential of the smoothed Huber loss $f_{\tau}(x)$. Instead of using the first-order iteration with the estimation equation (6) as in TD (2), we propose a Newton-type iterative estimator, which avoids the tuning of the learning rate. Newton-type estimators are often referred to as second-order methods when discussed in convex optimization. Nonetheless, the equation (5) is just for illustration purposes, which cannot be directly optimized since $\boldsymbol{\theta}_{\tau}^{*}$ appears in the objective function for minimization. On the contrary, our proposed method is a Newton-type method for solving the root-finding problem (6) using a sequence of the observations.

In this paper, we model two types of noise in observed rewards. The first is the classical Huber contamination model (Huber, 1992, 2004), where an $\alpha_{n}$-fraction of rewards comes from arbitrary distributions. The second is the heavy-tailed model, where the reward function may admit heavy-tailed distributions. In the following section, we show that our proposed estimator is robust to the aforementioned two noises.

In the following, we introduce an online Newton-type method for estimating the parameter $\boldsymbol{\theta}^{*}$ in the presence of outliers and heavy-tailed noise in the reward function ${\mathcal{R}}(s)$. For ease of presentation, we denote the observations by $\boldsymbol{X}_{i}=\boldsymbol{\phi}(s_{i})$, $\boldsymbol{Z}_{i}=\boldsymbol{\phi}(s_{i})-\gamma\boldsymbol{\phi}(s_{i+1})$, and $b_{i}={\mathcal{R}}(s_{i})$, where $s_{i}$ is the state at time $i$. Here $(\boldsymbol{X}_{i},\boldsymbol{Z}_{i},b_{i})\in\mathbb{R}^{d}\times\mathbb{R}^{d}\times\mathbb{R}$ for each observation in the sample. Our objective is to propose an online estimator by the estimation function (6), which can be rewritten as

based on a sequence of dependent observations $\{(\boldsymbol{X}_{i},\boldsymbol{Z}_{i},b_{i})\}_{i\geq 1}$.

At iteration $n+1$, our proposed Newton-type estimator updates $\widehat{\boldsymbol{\theta}}_{n+1}$ by

Here $\widehat{\boldsymbol{H}}_{n+1}$ is an empirical information matrix of the estimation equation (7), as

where $\tau_{i}$ is the thresholding parameter in the Huber loss. We let $\tau_{n}$ tend to infinity to eliminate the bias generated by the smoothed Huber loss.

It is noteworthy to mention that the matrix $\widehat{\boldsymbol{H}}_{n+1}$ is not the Hessian matrix of the objective function on the right-hand side of (5). As discussed in the previous section, (5) cannot be directly optimized, nor do they lead to M-estimation problems. Indeed, our proposed update (8) is a Newton-type method to find the root of (7) using a sequence of observations $\{(\boldsymbol{X}_{i},\boldsymbol{Z}_{i},b_{i})\}_{i\geq 1}$. In the following section, we will show the use of the matrix $\widehat{\boldsymbol{H}}_{n+1}$ helps for desirable convergence properties.

It should be noted that (8) can be implemented efficiently in a fully-online manner, i.e., without storage of the trajectory of historical information. Specifically, we write (8) as $\widehat{\boldsymbol{\theta}}_{n+1}=\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern 1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{\theta}$\kern-1.00006pt}}}\,_{n+1}-\widehat{\boldsymbol{H}}_{n+1}^{-1}\boldsymbol{G}_{n+1}$, with each of the item on the right-hand side being a running average. It is easy to see that the averaged estimator $\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern 1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{\theta}$\kern-1.00006pt}}}\,_{n+1}=\frac{1}{n+1}(n\,\hbox{\kern 0.09995pt\vbox{\hrule height=0.5pt\kern 1.42082pt\hbox{\kern-1.00006pt$\boldsymbol{\theta}$\kern-1.00006pt}}}\,_{n}+\widehat{\boldsymbol{\theta}}_{n})$ and the vector

can both be updated online. In addition, the inverse $\widehat{\boldsymbol{H}}_{n+1}^{-1}$ can be directly and efficiently computed online by the inverse recursion formulation. By the Sherman-Morrison formula, we have

Here we note that both terms $\frac{1}{n}\boldsymbol{Z}_{n+1}^{{\top}}\widehat{\boldsymbol{H}}_{n}^{-1}\boldsymbol{X}_{n+1}$ and $\{g_{\tau_{n+1}}^{\prime}(\boldsymbol{Z}_{n+1}^{{\top}}\widehat{\boldsymbol{\theta}}_{n}-b_{n+1})\}^{-1}$ inside the brackets on the right-hand side of (11) are scalars in $\mathbb{R}^{1}$, not matrices.

The complete algorithm is presented in Algorithm 1. We refer to it as the Robust Online Policy Evaluation (ROPE). Compared with existing stochastic approximation algorithms for TD learning (Durmus et al., 2021; Mou et al., 2021; Ramprasad et al., 2022), our ROPE algorithm does not need to tune the step size.

Since performing iterations of $\widehat{\boldsymbol{H}}_{n+1}^{-1}$ in (11) requires an initial invertible $\widehat{\boldsymbol{H}}_{n_{0}}$, we shall compute from the first $n_{0}$ samples that

which serves as the initial quantities of (10)–(11) in order to compute the forthcoming iterations. Here $\widehat{\boldsymbol{\theta}}_{0}$ is a given initial parameter, and $\tau_{0}$ is a specified initial threshold level.

In the above section, we provide an online Newton-type algorithm to estimate the parameter $\boldsymbol{\theta}^{*}$. As we can see in Theorem 4, the proposed estimator has an asymptotic variance with a sandwich structure $\boldsymbol{H}^{-1}\boldsymbol{\Sigma}(\boldsymbol{H}^{{\top}})^{-1}$. To conduct statistical inference simultaneously with ROPE method we proposed in Section 3, we propose an online plug-in estimator for this sandwich structure. Note that the online estimator $\widehat{\boldsymbol{H}}^{-1}_{n}$ of $\boldsymbol{H}^{-1}$ is proposed and utilized in (11) of the estimation algorithm. It remains to construct an online estimator for $\boldsymbol{\Sigma}$ in (17).

Developing an online estimator of $\boldsymbol{\Sigma}$ with dependent samples is intricate compared to the case of independent ones. As the samples are dependent, the covariance matrix $\boldsymbol{\Sigma}$ is an infinite sum of the series of time-lag covariance matrices. To estimate the above long-run covariance matrix $\boldsymbol{\Sigma}$, we first rewrite the definition of $\boldsymbol{\Sigma}$ in (17) into the following form

Next, we replace each term ${\mathbb{E}}\big{[}\boldsymbol{X}_{0}\boldsymbol{X}^{{\top}}_{k}(\boldsymbol{Z}^{{\top}}_{0}\boldsymbol{\theta}^{*}-b_{0})(\boldsymbol{Z}^{{\top}}_{k}\boldsymbol{\theta}^{*}-b_{k})\big{]}$ in (17) with its empirical counterpart using the $n$ samples. Meanwhile, to handle outliers, we replace $(\boldsymbol{Z}^{{\top}}_{i}\boldsymbol{\theta}^{*}-b_{i})$ by $g_{\tau_{i}}(\boldsymbol{Z}^{{\top}}_{i}\widehat{\boldsymbol{\theta}}_{i-1}-b_{i})$. In addition, the infinite sum in (17) is not feasible for direct estimation. Nonetheless, Condition (C1) on the mixing rate of $\{(\boldsymbol{X}_{t},\boldsymbol{Z}_{t})\}$ allows us to approximate it by effectively estimating the first $\lceil\lambda\log n\rceil$ terms only, where $\lambda$ is a pre-specified constant. In summary, we can construct the following estimator for the covariance matrix $\boldsymbol{\Sigma}$,

To perform a fully-online update for the estimator $\widehat{\boldsymbol{\Sigma}}_{n}$, we define

and the above equation (18) can be rewritten as,

In practice, at the $n$-th step, we keep only the terms $\{\boldsymbol{S}_{n},...,\boldsymbol{S}_{n-\lceil\lambda\log n\rceil\wedge(n-1)}\}$ in memory, and the covariance matrix can be updated by

The proposed estimator $\widehat{\boldsymbol{\Sigma}}_{n}$ complements the estimator $\widehat{\boldsymbol{H}}_{n}^{-1}$ in (11) to construct an estimator of the sandwich form $\boldsymbol{H}^{-1}\boldsymbol{\Sigma}(\boldsymbol{H}^{{\top}})^{-1}$ that appears in the asymptotic distribution in Theorem 4.