Generalized Gaussian Temporal Difference Error for Uncertainty-aware Reinforcement Learning

## Abstract

Abstract Conventional uncertainty-aware temporal difference (TD) learning methods often rely on simplistic assumptions, typically including a zero-mean Gaussian distribution for TD errors. Such oversimplification can lead to inaccurate error representations and compromised uncertainty estimation. In this paper, we introduce a novel framework for generalized Gaussian error modeling in deep reinforcement learning, applicable to both discrete and continuous control settings. Our framework enhances the flexibility of error distribution modeling by incorporating additional higher-order moment, particularly kurtosis, thereby improving the estimation and mitigation of data-dependent noise, i.e., aleatoric uncertainty. We examine the influence of the shape parameter of the generalized Gaussian distribution (GGD) on aleatoric uncertainty and provide a closed-form expression that demonstrates an inverse relationship between uncertainty and the shape parameter. Additionally, we propose a theoretically grounded weighting scheme to fully leverage the GGD. To address epistemic uncertainty, we enhance the batch inverse variance weighting by incorporating bias reduction and kurtosis considerations, resulting in improved robustness. Extensive experimental evaluations using policy gradient algorithms demonstrate the consistent efficacy of our method, showcasing significant performance improvements.

## 1 Introduction

Deep reinforcement learning (RL) has demonstrated promising potential across various real-world applications, e.g., finance (Moody & Saffell, ; Byun et al., ; Sun et al., ), and autonomous driving (Kahn et al., ; Emamifar & Ghoreishi, ; Hoel et al., ).
One critical avenue for improving the performance and robustness of RL agents in these complex, high-dimensional environments is the quantification and integration of *uncertainty* associated with the decisions made by agents or the environment (Lockwood & Si, ).
Effective management of uncertainty promotes the agents to make more informed decisions leading to enhanced sample efficiency in RL context, which is particularly beneficial in unseen or ambiguous situations.

Temporal difference (TD) learning is a fundamental component of many RL algorithms, facilitating value function estimation and policy derivation through iterative updates (Sutton, ).
Traditionally, these TD updates are typically grounded in $L_{2}$ loss, corresponding to maximum likelihood estimation (MLE) under the assumption of Gaussian error.
Such simplification may be overly restrictive, especially considering the noisy nature of TD errors, which are based on constantly changing estimates of value functions and policies.
This assumption compromises sample efficiency, necessitating the incorporation of additional distributional parameters for flexible but computationally efficient TD error modeling.

In statistics and probability theory, distributions are typically characterized by their central tendency, variability, and shape (DeCarlo, ; Milton et al., ).
Traditional deep RL methods, however, effectively exploit on the variance of the error distribution through the scale parameter, yet they often disregard its *shape*.
This oversight hinders these methods from fully capturing the true underlying uncertainty.
The kurtosis, the scale-independent moment, does significantly influence both inferential and descriptive statistics (Balanda & MacGillivray, ), and the reliability of uncertainty estimation.

Therefore, it is essential to incorporate the shape of the error distribution into TD learning to better reflect uncertainties present in RL environments, enabling more robust and reliable decision-making processes in dynamic and complex scenarios.

Figure: Figure 1: Generalized Gaussian distribution.
Refer to caption: extracted/5894656/figures/ggd.png

A notable enhancement to the normality hypothesis is the use of the *generalized Gaussian distribution* (GGD), also known as the generalized error distribution or exponential power distribution.
This flexible family of symmetric distributions, as depicted in Figure 1, encompasses a wide range of classical distributions, including Gaussian, Laplacian, and uniform distributions, all adjustable via a shape parameter (Box & Tiao, ).
This specific parameter allows for fine-tuning the distribution to match the characteristics of TD error distribution, providing a more reliable representation of uncertainty.

##### Contributions

Our primary contribution is the introduction of a novel framework of generalized Gaussian error modeling in deep RL, enabling robust training methodologies by incorporating the distribution’s shape.
This approach addresses both data-dependent noise, i.e., aleatoric uncertainty, and the variance of model estimates, i.e., epistemic uncertainty, ultimately enhancing model stability and performance.

The key contributions of our work are as follows:

- 1.
Empirical investigations (Section 3.1.1):
We conduct empirical investigations of TD error distributions, revealing substantial deviations from the Gaussian distribution, particularly in terms of tailedness.
These findings underscore the limitations of conventional Gaussian assumptions.
- 2.
Theoretical exploration (Section 3.1.2):
Building on empirical insights, we explore the theoretical suitability of modeling TD errors with a GGD.
Theorem 1 demonstrates the effectiveness and well-definedness of the proposed method under leptokurtic error distributions, characterized by $\beta\in(0,2]$.
Our experimental results suggest that the estimates of $\beta$ mostly converge within this range, aligning with the empirical findings.
- 3.
Aleatoric uncertainty mitigation (Section 3.1):
We investigate the implications of the distribution shape on the estimation and mitigation of aleatoric uncertainty.
GGD error modeling enables the quantification of aleatoric uncertainty in a closed form, indicating a negative relation to the shape parameter $\beta$ on an exponential scale, with a constant scale parameter $\alpha$.
We also leverage the second-order stochastic dominance of GGD to weight error terms proportional to $\beta$, enhancing the model’s robustness to heteroscedastic aleatoric noise by focusing on less spread-out samples.
- 4.
Epistemic uncertainty mitigation (Section 3.2):
We introduce the batch inverse error variance weighting scheme, adapted from the batch inverse variance scheme (Mai et al., ), to account for both variance and kurtosis of the estimation error distribution.
This scheme down-weights high-variance samples to prevent noisy data and improves model robustness by focusing on reliable error estimates.
- 5.
Experimental evaluations (Section 4):
We conduct extensive experimental evaluations using policy gradient algorithms, demonstrating the consistent efficacy of our method and significant performance enhancements.

## 2 Background

We consider a Markov decision process (MDP) governed by state transition probability ${\mathcal{P}}(s_{t+1}|s_{t},a_{t})$, with $s_{t}\in{\mathcal{S}}$ and $a_{t}\in{\mathcal{A}}$ represent the state and action at step $t$, respectively (Sutton & Barto, ).
Within this framework, an agent interacts with the environment via a policy $\pi(a_{t}|s_{t})$, leading to the acquisition of rewards $r(s_{t},a_{t})\sim{\mathcal{R}}(s_{t},a_{t})$.

Numerous model-free deep RL algorithms leverage TD updates for value function approximation (Mnih et al., , ; Schulman et al., ; Fujimoto et al., ; Haarnoja et al., ).
In these methods, neural networks parameterized by $\theta$ are trained to approximate the state-action value $Q(s_{t},a_{t})$ by minimizing the error between the target and predicted value:

$$ $\delta(s_{t},a_{t};\theta)=T(s_{t},a_{t};\theta)-Q(s_{t},a_{t};\theta),$ (1) $$

where the target is computed according to Bellman’s equation:

$$ $T(s_{t},a_{t};\theta)=r(s_{t},a_{t})+\gamma Q(s_{t+1},a_{t+1};\theta).$ (2) $$

Typically, TD updates involve minimizing the mean squared error (MSE) loss $\operatorname{MSE}_{\theta}=\mathbb{E}[(T(s_{t},a_{t};\theta)-Q(s_{t},a_{t};
\theta))^{2}]$ through stochastic gradient descent.
This minimization implicitly presumes that errors conform to a Gaussian distribution with zero-mean, consistent with the principles of MLE.

For clarity, we henceforth omit the learnable parameter and adopt subscript notation for function arguments, e.g., $\delta_{t}=T_{t}-Q_{t}$.

### 2.1 Uncertainty

Uncertainty in neural networks is commonly decomposed into two sources: aleatoric and epistemic (Der Kiureghian & Ditlevsen, ; Kendall & Gal, ; Depeweg et al., ; Valdenegro-Toro & Mori, ).
Epistemic uncertainty arises from limitations within the neural network and can potentially be reduced through further learning or model improvements.
In contrast, aleatoric uncertainty stems from the inherent stochasticity of the environment or the dynamics of the agent-environment interactions and is fundamentally irreducible.

This distinction is crucial in the context of RL, where areas with high epistemic uncertainty need to be explored, whereas exploring areas with high aleatoric uncertainty may lead to ineffective training, since the agent might have adequate knowledge but insufficient information for decisive actions.
Quantifying aleatoric uncertainty is known to facilitate learning dynamics of stochastic processes and enables risk-sensitive decision making (Dabney et al., ; Vlastelica et al., ; Seitzer et al., ).

To address aleatoric uncertainty, variance networks, denoted by $Q^{\sigma}$ to distinguish them from the value approximation head $Q^{\mu}$, are frequently employed (Bishop, ; Nix & Weigend, ; Kendall & Gal, ; Lakshminarayanan et al., ; Mai et al., ; Mavor-Parker et al., ).
The estimated variance from $Q^{\sigma}$ is utilized for heteroscedastic Gaussian error modeling, with $\delta_{t}\sim{\mathcal{N}}(Q^{\mu}_{t},{Q^{\sigma}_{t}}^{2})$, thereby quantifying aleatoric uncertainty (Seitzer et al., ).
This network is trained by maximizing the log likelihood function of the Gaussian distribution or minimizing the Gaussian negative log likelihood (NLL) loss.
The variance of estimated values is also used for batch inverse variance (BIV) weighting in inverse-variance RL (Mai et al., ).
This regularization effectively captures epistemic uncertainty with a regularizing temperature $\lambda$:

$$ $\displaystyle\mathcal{L}$ $\displaystyle=\mathcal{L}_{\text{GD-NLL}}+\lambda\mathcal{L}_{\text{BIV}}$ (3) $\displaystyle=\sum_{t}{\left((\delta_{t}/{Q^{\sigma}_{t}})^{2}+\log{{Q^{\sigma }_{t}}^{2}}\right)+\lambda\left(\frac{\omega^{\text{BIV}}_{t}}{\sum_{\tau} \omega^{\text{BIV}}_{\tau}}\delta_{t}^{2}\right)},$ $$

where the BIV weight $\omega^{\text{BIV}}_{t}=1/(\gamma^{2}{{\mathbb{V}}[Q^{\mu}_{t}]}+\xi)$ with empirical variance ${\mathbb{V}}$, and $\xi$ is either a hyperparameter or numerically computed to ensure a sufficient effective batch size.
In the official implementation, the sample variance with Bessel’s correction is used as a variance estimator, which is particularly apparent given a small sample size, i.e., an ensemble number of five.

### 2.2 Tailedness

While mainstream machine learning literature often prioritizes on capturing central tendencies, the significance of extreme events residing in the tails for enhancing performance and gaining a deeper understanding of learning dynamics cannot be overlooked.
This is especially relevant for MLE base on the normality assumption, which is commonly applied in variance network frameworks.
Focusing solely on averages or even deviations is proven to be inadequate in the presence of outlier samples (David, ; Gather & Kale, ).

For instance, consider the impact of non-normal samples on the estimate of the variance, as described in 1 with proof presented in Section B.1.

###### Proposition 1 (Biased variance estimator ( Yuan et al. , 2005 ) ) .

Let $X_{1},X_{2},...,X_{n}$ be a sequence of independent, non-normally distributed random variables from a population $X$ with mean $\mu$, variance $\sigma^{2}$, and kurtosis $\kappa$.
Then, with the MLE estimators under normality assumption, i.e., $\hat{\mu}=\sum_{i=1}^{n}X_{i}/n$ for mean and $\hat{\sigma}^{2}=\sum_{i=1}^{n}(X_{i}-\hat{\mu})^{2}/n$ for variance, the variance estimator $\hat{\sigma}^{2}$ exhibits bias.
Specifically, it will be negatively biased when $\kappa>0$ and positively biased when $\kappa<0$.

1 elucidates that the standard error of estimated TD error variance is a function of kurtosis $\kappa$.
With heavy-tail distributions, the standard error of variance estimates, through networks trained by Gaussian NLL (Nix & Weigend, ; Bishop, ), may also be underestimated, highlighting the influence of kurtosis on variance estimation.
Furthermore, it is shown that the likelihood ratio statistics for variance estimator depends on kurtosis even for large $n$ (Yuan et al., ).
This emphasizes the necessity of accounting for tailedness to derive robust variance estimates confidently.

###### Remark 1 (Varietal variance estimator ( Burch , 2014 ) ) .

The dependence of the standard error of the variance estimator on kurtosis impacts both the bias and variance of variance estimation.
Specifically, when the kurtosis of the underlying distribution exceeds zero, assuming normality tends to result in an underestimation of the confidence interval of variance estimation.

###### Proposition 1 (Biased variance estimator ( Yuan et al. , 2005 ) ) .

###### Remark 1 (Varietal variance estimator ( Burch , 2014 ) ) .

##### Gumbel error modeling

Recent applications of Gumbel distribution closely related to TD learning have emerged for estimating the maximum $Q$ value in the Bellman update process (Garg et al., ; Hui et al., ).
These approaches capitalize on the foundation of the extreme value theorem, which states that maximal values drawn from any exponential-tailed distribution follow the Gumbel distribution (Fisher & Tippett, ; Mood, ).
Compared to conventional distributional RL algorithms employing Gaussian mixture or quantile regression (Dabney et al., ; Shahriari et al., ), this approach showcases superior control performance.
However, it has been observed that while Gumbel modeling initially aligns with the error distribution propagated through the chain of $\max$ operations, its Gumbel-like attribute may diminish over the course of training (Garg et al., ).
We instead propose a novel approach utilizing GGD, which offers flexibility in expressing the tail behavior of diverse distributions.
This method is adaptable to wider range of MDPs, even those without $\max$ operators.

## 3 Methods

Our approach introduces enhancements to the loss function derived from equation 3, incorporateing tailedness into both loss attenuation and regularization terms:

$$ $\displaystyle\mathcal{L}$ $\displaystyle=\mathcal{L}^{\text{RA}}_{\text{GGD-NLL}}+\lambda\mathcal{L}_{ \text{BIEV}}$ (4) $\displaystyle=\sum_{t}{\frac{\omega^{\text{RA}}_{t}}{\sum_{\tau}\omega^{\text{ RA}}_{\tau}}\left(\left(\left|\delta_{t}\right|/Q^{\alpha}_{t}\right)^{Q^{ \beta}_{t}}-\log{Q^{\beta}_{t}/Q^{\alpha}_{t}}+\log{\Gamma(1/Q^{\beta}_{t})} \right)+\lambda\left(\frac{\omega^{\text{BIEV}}_{t}}{\sum_{\tau}\omega^{\text{ BIEV}}_{\tau}}\delta_{t}^{2}\right)},$ $$

where $Q^{\alpha}$ and $Q^{\beta}$ represent the alpha and beta networks, respectively.
Here, risk-averse weights $\omega^{\text{RA}}_{t}=Q^{\beta}_{t}$, and batch inverse error variance (BIEV) weights $\omega^{\text{BIEV}}_{t}=1/({\mathbb{V}}[\delta_{t}]+\xi)$.
These dual loss terms effectively capture aleatoric and epistemic uncertainties, respectively.
The subsequent subsections provide a detailed rationale for this modification.

### 3.1 Generalized Gaussian error modeling

One simple yet promising approach to address non-normal heteroscedastic error distributions involves modeling the per-error distribution using a zero-mean symmetric GGD (Zeckhauser & Thompson, ; Chai et al., ; Giacalone, ; Upadhyay et al., ):

$$ $\delta\sim\text{GGD}(0,\alpha,\beta)=\frac{\beta}{2\alpha\Gamma(1/\beta)}e^{- \left(|\delta|/\alpha\right)^{\beta}},$ (5) $$

where $\alpha$ and $\beta$ represent the scale and shape parameter, respectively.
This method not only allows for modeling each non-identical error by parameterizing the GGD with different $\alpha_{t}$ and $\beta_{t}$ at step $t$, but also offers a flexible parametric form that adapts across a spectrum of classical distributions from Gaussian to uniform as $\beta$ increases to infinity (Giller, ; Nadarajah, ; Novey et al., ; Dytso et al., ).

The shape parameter $\beta$ serves as a crucial structure characteristic, distinguishing underlying mechanisms.
The kurtosis $\kappa$, commonly used to discern different distribution shapes, is solely a function of $\beta$ and is defined as Pearson’s kurtosis minus three to emphasize the difference from Gaussian distribution (DeCarlo, ):

$$ $\kappa=\frac{\Gamma(5/\beta)\Gamma(1/\beta)}{\Gamma(3/\beta)^{2}}-3.$ (6) $$

This implies that distributions with $\beta<2$ are leptokurtic, i.e., $\kappa>0$, indicating a higher frequency of outlier errors compared to the normal error distribution.
With only one additional parameter to characterize the distribution, GGD effectively expresses differences in tail behavior, a capability distributions dependent solely on location or scale parameters lack.

###### Remark 2 .

Despite the GGD having three parameters, we only employ the $\beta$ estimation for GGD error modeling to minimize computational overhead, setting the scale parameter $\alpha$ to one.
While this may slightly limit the expressivity of the error model, it significantly enhances training stability.
The impact of omitting the $\alpha$ parameter is minimal, as moments influenced by $\alpha$ can also be represented by $\beta$.
Additionally, this simplification offers implementation advantages, requiring only a slight change in the loss function to migrate from variance networks, i.e., from Gaussian to GGD NLL.

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/tde-sac-ant.png

###### Remark 2 .

#### 3.1.1 Empirical evidence

Figure 2 presents empirical findings that reveal a significant deviation from Gaussian distribution in TD errors, evidenced by well-fitted GGDs and pronounced differences in the shape of distributions.
This non-normality becomes particularly apparent when contrasting initial and final evaluations, suggesting an increasing prominence of tailedness throughout the training process.

We hypothesize that such departure from normality stems from the exploratory nature of agent behavior.
During exploration, agents frequently encounter *unexpected* states and rewards, leading to a broader spectrum of TD errors than those seen in purely exploitative scenarios.
This increased variability likely contributes to the emergence of non-normal distributions, characterized by heavier tails.

Furthermore, the observed evolution in the tails of TD errors underscores the shifting interplay between aleatoric and epistemic uncertainties.
As training advances, aleatoric errors tend to become more influential while epistemic uncertainty diminishes, potentially resulting in non-normally distributed errors with heavier tails.
Notably, the diminishing property of Gumbel-like attribute, discussed in Section 2.2, is also reflected in the evolving TD error distributions.

#### 3.1.2 Theoretical analysis

Given the empirical suitability of GGD for modeling TD errors, we conduct an in-depth theoretical examination of GGD.
We begin by demonstrating the well-definedness of GGD regression, with a specific emphasis on the positive-definiteness of its PDF under certain conditions.

###### Theorem 1 (Positive-definiteness ( Bochner , 1937 ; Ushakov , 2011 ; Dytso et al. , 2018 ) ) .

The NLL of GGD is well-defined for $\beta\in(0,2]$.

This theorem guarantees the well-definedness of the NLL of the GGD.
It ensures that the PDF is guaranteed to be positive everywhere under highly probable conditions of $\beta$, thus affirming the suitability of the NLL as a loss function for a valid probability distribution.
In fact, it is known that parameter estimation for GGD can also be numerically accomplished by minimizing the NLL function, with asymptotic normality, consistency, and efficiency of the estimates ensured under suitable regularity conditions (Agro, ).
This positive-definiteness not only implies a theoretical property but also has practical implications in deep RL, where ensuring stability and convergence is crucial.
The positive-definiteness of the PDF of GGD shown in the proof, elaborated in Section B.2, also assures that the function integrates to one.

The shape parameter $\beta$ also introduces a desirable property of *risk-aversion* to GGD, which can be mathematically formulated by stochastic dominance (Levy, ; Martin et al., ).
Stochastic dominance, a concept assessing random variables via a stochastic order, reflects the shared preferences of rational decision-makers.

###### Theorem 2 (Second-order stochastic dominance ( Dytso et al. , 2018 ) ) .

Consider two random variables $X_{1}\sim\text{GGD}(0,\alpha,\beta_{1})$ and $X_{2}\sim\text{GGD}(0,\alpha,\beta_{2})$, where $\alpha>0$ and $\beta_{1},\beta_{2}>0$.
If $\beta_{1}\leq\beta_{2}$, then $X_{2}$ exhibits second-order stochastic dominance over $X_{1}$.
This dominance implies, for all $x$,

$$ $\int_{-\infty}^{x}\left[F_{X_{1}}(t)-F_{X_{2}}(t)\right]dt\geq 0,$ (7) $$

where $F$ denotes the cumulative distribution function.

The above theorem, with its proof detailed in Section B.3, suggests second-order stochastic dominance of TD errors.
This dominance signifies a preference for risk-aversion, wherein the dominant variable $X_{2}$ exhibits greater predictability and maintains expectations that are equal to or higher than those of $X_{1}$ for all concave and non-decreasing functions (Osband & Van Roy, ).
Formally, for such a function $u:{\mathbb{R}}\to{\mathbb{R}}$, we have $\mathbb{E}[u(X_{2})]\geq\mathbb{E}[u(X_{1})]$, with $\mathbb{E}$ denoting expectation.

The dominance relationship among GGD random variables, determined by the shape parameter $\beta$, reinforces the suitability of GGD for modeling errors in TD learning.
Leveraging this risk-aversion, the application of the risk-averse weighting $\omega^{\text{RA}}_{t}=Q^{\beta}_{t}$ in $\mathcal{L}^{\text{RA}}_{\text{GGD-NLL}}$ capitalizes on the tendency of GGD to learn from relatively less spread-out samples, thereby enhancing robustness to heteroscedastic noise.

###### Remark 3 .

The training of the critic is more directly influenced by aleatoric uncertainty, since only the critic loss is a direct function of the state, action, and reward, with the actor being downstream of the critic in uncertainty propagation.
GGD error modeling enables us to quantify aleatoric uncertainty as a closed form, i.e., $\sigma^{2}=\alpha^{2}\Gamma(3/\beta)/\Gamma(1/\beta)$ (Upadhyay et al., ).
Remarkably, aleatoric uncertainty exhibits a negative proportionality to the shape parameter $\beta$ on an exponential scale, with a constant scale parameter $\alpha=1$ adopted in our implementation employing a beta head exclusively.
Building on this, risk-averse weighting mitigates the negative impacts of noisy supervision by assigning higher weights to errors with lower aleatoric uncertainty for the loss attenuation term.

###### Theorem 1 (Positive-definiteness ( Bochner , 1937 ; Ushakov , 2011 ; Dytso et al. , 2018 ) ) .

###### Theorem 2 (Second-order stochastic dominance ( Dytso et al. , 2018 ) ) .

###### Remark 3 .

### 3.2 Batch inverse error variance regularization

When estimating the uncertainty of target estimates, as employed in BIV weighting, potential bias can arise (Janz et al., ; Liang et al., ).
Conversely, the bias of TD errors remains notably small with the assumption of constant value approximation bias (Flennerhag et al., ).
Motivated by this, we propose the BIEV weighting:

$$ $\omega^{\text{BIEV}}_{t}=\frac{1}{{\mathbb{V}}[\delta_{t}]},$ (8) $$

incorporating the concept of *error variance* explicitly into BIV weighting.

Recent investigations have explored advancements in variance estimation, particularly through a constant multiplier, i.e., $s_{\omega}^{2}=\omega(n-1)s^{2}$ (Kleffe, ; Wencheko & Chipoyera, ), where $s^{2}$ denotes sample variance, the MLE estimator of variance.
Although non-standard weights $\omega\neq 1/(n-1)$ may introduce bias in variance estimation, the estimation of inverse variance remains biased due to Jensen’s inequality, even with the use of the unbiased estimator $s^{2}$ (Walter et al., ).
Consequently, our focus shifts to relative efficiency (RE), where we derive the MSE-best biased estimator (MBBE) in 2.

###### Proposition 2 (MBBE of variance ( Searls & Intarapanich , 1990 ; Wencheko & Chipoyera , 2009 ) ) .

Let $s_{\omega}^{2}=\omega(n-1)s^{2}$ be the adjusted variance estimator, with the sample variance $s^{2}$ and the weight $\omega$ being a function of the sample size $n$ and the population kurtosis $\kappa$.
Then, the estimator with the least MSE is given by:

$$ $s_{\omega^{*}}^{2}=\left(\frac{\kappa}{n}+\frac{n+1}{n-1}\right)^{-1}s^{2}.$ (9) $$

Additionally, MBBE of variance $s_{\omega^{*}}^{2}$ has consistent superior efficiency over the sample variance $s^{2}$, i.e.:

$$ $\text{RE}_{n}=\frac{{\mathbb{V}}[s^{2}]}{\operatorname{MSE}(s_{\omega^{*}}^{2} )}=1+\frac{\kappa}{n}+\frac{2}{n-1}>1.$ (10) $$

The reciprocal relationship between sample size and the impact of kurtosis on variance estimation is intuitive, especially for small samples where kurtosis is much more significant.
Consequently, we advocate for the adoption of the MBBE in epistemic uncertainty estimation.
It’s worth noting that while the derivation of improved estimators presupposes known kurtosis, our method differs by utilizing estimated kurtosis.

As BIV weighting, the BIEV weighting plays a crucial role in enhancing the robustness and efficiency TD updates.
By normalizing the weight of each sample in a batch relative to the scale of its epistemic uncertainty compared to other samples, BIEV weighting ensures that the model appropriately accounts for the reliability of the estimate from each data point, resulting in robustness against inaccuracies in variance estimate calibration.

###### Proposition 2 (MBBE of variance ( Searls & Intarapanich , 1990 ; Wencheko & Chipoyera , 2009 ) ) .

## 4 Experiments

We conduct a comprehensive evaluation of our proposed method across well-established benchmarks, including MuJoCo (Todorov et al., ), and discrete control environments from Gymnasium (Towers et al., ).
Notably, we augment the discrete control environments through the introduction of supplementary uniform action noise to enhance environmental fidelity.

To underscore the versatility and robustness of our approach, we deliberately choose baseline algorithms that cover a wide spectrum of RL paradigms.
Specifically, we employ soft actor-critic (SAC) (Haarnoja et al., ), an off-policy $Q$-based policy gradient algorithm, and proximal policy optimization (PPO) (Schulman et al., ), an on-policy $V$-based method.
We focus on adversaries limited to variance networks due to the use of separate target networks in previously mentioned Gumbel error modeling methods.
This constraint is intended to focus on computationally efficient algorithms that only incorporate an additional layer, referred to as a *head*.

The algorithms are implemented using PyTorch (Paszke et al., ), within the Stable-Baselines3 framework (Raffin et al., ).
We use default configurations from Stable-Baselines3, with adaptations limited to newly introduced hyperparameters.
For both PPO and SAC, along with their variants, we employ five ensembled critics.
The parameter $\xi$ from equation 4 is computed with a minimum batch size of 16, and the regularizing temperature $\lambda$ is set to $0.1$.
Additional experimental details are provided in Appendix C.

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/result-sac-ant.png

The performance of SAC across different MuJoCo environments is presented in Figure 3.
While the variance head degrades performance in certain scenarios, SAC variants employing the beta head consistently lead to better sample efficiency and asymptotic performance.
Notable improvements are observed in the HalfCheetah-v4 and Hopper-v4 environments, where the variance head substantially reduces sample efficiency.
This suggests that the TD error distributions in these environments may exhibit heavy tails.
The impact of BIEV regularization varies by environment but generally performs at least as well as BIV regularization.

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/param-sac-ant.png

Figure 4 shows the coefficients of variation, defined as the ratio of the standard deviation to the mean, i.e., $\sqrt{{\mathbb{V}}[X]}/\mathbb{E}[X]$, for the variance and $\beta$ estimates.
This statistic demonstrates the scale-invariant volatility of the parameter estimation, given that the scale of the estimated variance is significantly larger than that of the $\beta$ estimates.
A lower coefficient of variation in $\beta$ estimation indicates greater stability compared to variance estimation.

The convergence of $\beta$ estimation is also more stable than variance estimation.
This observation aligns with 1, suggesting a potential underestimation of confidence intervals in variance estimation.
Considering the susceptibility of variance estimates to extreme values, such underestimation introduces considerable uncertainty in parameter estimation.
Our hypothesis regarding the escalating impact of extreme TD errors throughout training is consistent with this findings, as it exacerbates the challenge in variance estimation, leading to volatility or divergence of the variance head.

These findings support that utilizing the beta head results in lower and converging coefficients of variation in parameter estimation.

Figure: (a) CartPole-v1
Refer to caption: extracted/5894656/figures/result-ppo-cart-pole.png

Figure 5 demonstrates the results of training PPO on MuJoCo and discrete control environments with additional noise.
Remarkably, the incorporation of the beta head and BIEV regularization yields similar outcomes to those observed in SAC experiments.
This indicates the efficacy of GGD error modeling in state value $V$-based TD learning as well.

We present comprehensive ablation studies in Appendix D, examining the efficacy of key components, including risk-averse weighting, regularizing temperature $\lambda$, BIEV regularization applied to SAC and its Gaussian variant, and the integration of the alpha head.

## 5 Discussion

In this paper, we advocate for and substantiate the integration of GGD modeling for TD error analysis.
Our main contribution is the introduction of a novel framework that enables robust training methodologies by leveraging the distribution’s shape.
This approach accounts for both data-dependent noise, i.e., aleatoric uncertainty, and the uncertainty of value estimation, i.e., epistemic uncertainty, ultimately enhancing the model’s stability and accuracy.

##### Further investigation

An imperative avenue for further investigation is the application of GGD within the context of maximum entropy RL.
Similar to how the Gaussian distribution maximizes entropy with constraints up to the second moment, the GGD maximizes entropy subject to a constraint on the $p$-th absolute moment (MacKay, ).
Exploring higher moments of the distribution could provide new insights into maximum entropy RL frameworks.

The absence of a comprehensive regret analysis in our current study also presents an opportunity for future work.
Considering that aleatoric uncertainty in TD learning predominantly arises from reward dynamics, conducting a regret analysis on heavy-tailed TD error is warranted.
This is particularly relevant as previous research has extensively studied regret in the context of heavy-tailed rewards.

While our experiments focus on policy gradient methods, the implications of TD error tailedness extend to the $Q$-learning family of algorithms.
We provide empirical findings on this extension in Appendix E, which demonstrate the broader applicability of our approach.
Additionally, exploring the generalized extreme value distribution could further enrich our understanding of tailedness phenomena due to its close association with extreme value theory.

##### Relevant applications

The implications of tailedness in TD error distributions extend into various domains, notably robust RL and risk-sensitive RL.
The focus in robust RL lies on developing algorithms that are less sensitive to noise and outliers within the reward signal.
Recognizing potential deviations from normality in TD error distributions is critical for designing such algorithms.
Our work emphasizes the importance of considering non-normal error distributions, especially the tail behaviors, to enhance the robustness of RL algorithms.

Another significant direction is risk-sensitive RL, which seeks to assess and mitigate the risks associated with different action choices.
In noisy and outlier-prone environments, capturing the risk profile using a Gaussian assumption for TD errors might be inadequate.
By considering the GGD, which better models the tail behaviors of error distributions, we can develop more accurate and reliable risk-sensitive RL algorithms.

In summary, our exploration into GGD modeling of TD errors opens several promising research directions and applications, emphasizing the need to consider non-normal error distributions for enhancing the robustness and risk-sensitivity of RL algorithms.

## Acknowledgments and Disclosure of Funding

We are grateful for the helpful advice provided by Dr. Heather Battey in the Department of Mathematics, Imperial College London.
Her insights and feedback have been invaluable to the development and improvement of this work.

## Appendix A Extended results

### A.1 Temporal difference error plots

We present the distributions of TD errors sampled at the initial and final evaluation steps, depicted in Figures 2 and 6 for SAC, and Figure 7 for PPO, which highlights the heavy tailedness of TD errors and the tendency converge to heavy tail throughout training.
This finding also emphasizes how aleatoric uncertainty affects their distribution, as elaborated in Section 3.1.
Interestingly, both state-action values $Q$ and state values $V$ demonstrate similar characteristics in their TD error distributions.

Figure: (a) HalfCheetah-v4
Refer to caption: extracted/5894656/figures/tde-sac-halfcheetah.png

Figure: (a) HalfCheetah-v4
Refer to caption: extracted/5894656/figures/tde-ppo-halfcheetah.png

### A.2 Coefficients of variation of parameter estimation

Figure: (a) HalfCheetah-v4
Refer to caption: extracted/5894656/figures/param-sac-half-cheetah.png

### A.3 PPO on other environments

Figure: (a) Hopper-v4
Refer to caption: extracted/5894656/figures/result-ppo-hopper.png

### A.4 Parameter estimation

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/param2-sac-ant.png

The parameter estimates from the beta head reveal consistently low values of $\beta$ across all environments, indicating a leptokurtic TD error distribution, as depicted in Figure 10.
These findings align with the observations from the TD error plots, where TD errors exhibit a closer resemblance to the GGD rather than a Gaussian distribution.

## Appendix B Proofs

### B.1 Of 1

###### Proof.

Consider a finite sample $X_{1},X_{2},...,X_{n}$ of independent, normally distributed random variables with $X\sim{\mathbb{P}}_{\theta_{0}}$, where $\theta_{0}=(\mu,\sigma)\in\Theta$ represents the true generative parameters.
Under this assumption, both skewness $\gamma$ and kurtosis $\kappa$ are zero.
Consequently, the moments are given by $\mathbb{E}[X]=\mu$, $\mathbb{E}\left[(X-\mu)^{2}\right]=\sigma^{2}$, $\mathbb{E}\left[(X-\mu)^{3}\right]=\sigma^{3/2}\gamma=0$ and $\mathbb{E}\left[(X-\mu)^{4}\right]=\sigma^{4}(\kappa+3)=3\sigma^{4}$.

The MLE estimator of $\mu$ and $\sigma$ is $\hat{\theta}=(\hat{\mu},\hat{\sigma})$, given by

$$ $\hat{\mu}=\bar{X}=\frac{1}{n}\sum_{i=1}^{n}X_{i},\quad\text{and}\quad\hat{ \sigma}^{2}=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\bar{X})^{2}.$ $$

It is well known to be consistent (Casella & Berger, ).

Assuming $\hat{\theta}\overset{p}{\to}\theta_{0}$, where $\overset{p}{\to}$ denotes the convergence in probability, under appropriate regularity conditions, then the asymptotic normality theorem of Cramer leads to

$$ $\sqrt{n}\left(\hat{\theta}-\theta_{0}\right)\overset{p}{\to}{\mathcal{N}}(0,{ \mathcal{I}}(\theta_{0})^{-1}),$ $$

as $n\to\infty$ (Ferguson, ).
Here, ${\mathcal{I}}(\theta_{0})=\begin{pmatrix}{\mathcal{I}}_{11}&{\mathcal{I}}_{12}
\\
{\mathcal{I}}_{21}&{\mathcal{I}}_{22}\end{pmatrix}$ is the Fisher information matrix.

Through straightforward calculations involving the log likelihood function derivatives, we obtain

$$ ${\mathcal{I}}(\theta_{0})=\begin{pmatrix}{\mathcal{I}}_{11}&{\mathcal{I}}_{12} \\ {\mathcal{I}}_{21}&{\mathcal{I}}_{22}\end{pmatrix}=\begin{pmatrix}1/\sigma^{2} &0\\ 0&1/2\sigma^{4}\end{pmatrix}.$ $$

Thus, as $n\to\infty$, $\sqrt{n}(\hat{\theta}-\theta_{0})\overset{p}{\to}{\mathcal{N}}\left(0,\begin{
pmatrix}\sigma^{2}&0\\
0&2\sigma^{4}\end{pmatrix}\right)$.

Now, consider a finite sample $X_{1},X_{2},...,X_{n}$ of independent, non-normally distributed random variables, with MLE estimator $\hat{\theta}=(\hat{\mu},\hat{\sigma})$.
By applying analogous reasoning, we seek the asymptotic distribution of $\sqrt{n}(\hat{\theta}-\theta_{0})$.
Letting $\tilde{\sigma}^{2}=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\mu)^{2}$ gives us

$$ $\displaystyle\hat{\sigma}^{2}$ $\displaystyle=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\hat{\mu})^{2}$ $\displaystyle=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-\mu)^{2}-(\hat{\mu}-\mu)^{2}$ $\displaystyle=\tilde{\sigma}^{2}-\frac{1}{n}\sum_{i=1}^{n}(\hat{\mu}-\mu)^{2}.$ $$

Therefore, as $\sqrt{n}(\hat{\mu}-\mu)^{2}\overset{p}{\to}0$,

$$ $\displaystyle\sqrt{n}(\hat{\theta}-\theta_{0})$ $\displaystyle=\sqrt{n}\left(\hat{\mu}-\mu,\hat{\sigma}^{2}-\sigma^{2}\right)$ (11) $\displaystyle=\sqrt{n}\left(\hat{\mu}-\mu,\tilde{\sigma}^{2}-\sigma^{2}\right) -\left(0,\sqrt{n}(\hat{\mu}-\mu)^{2}\right)$ $\displaystyle\overset{p}{\to}\sqrt{n}\left(\hat{\mu}-\mu,\tilde{\sigma}^{2}- \sigma^{2}\right).$ $$

Denoting $\tilde{\theta}=(\hat{\mu},\tilde{\sigma}^{2})$, equation 11 states that $\hat{\theta}$ and $\tilde{\theta}$ are asymptotically equivalent, i.e.,

$$ $\left(\hat{\theta}-\tilde{\theta}\right)\overset{p}{\to}0.$ $$

By Slutsky’s theorem (Ferguson, ), $\tilde{\theta}$ has the same asymptotic normality as $\hat{\theta}$, i.e.,

$$ $\sqrt{n}\left(\tilde{\theta}-\theta_{0}\right)\overset{p}{\to}{\mathcal{N}}(0, {\mathcal{K}}),$ $$

where ${\mathcal{K}}=\begin{pmatrix}{\mathcal{K}}_{11}&{\mathcal{K}}_{12}\\
{\mathcal{K}}_{21}&{\mathcal{K}}_{22}\end{pmatrix}$ is a covariance matrix with

$$ $\displaystyle{\mathcal{K}}_{11}$ $\displaystyle=\mathbb{E}\left[(X_{i}-\mu)^{2}\right]=\sigma^{2},$ $\displaystyle{\mathcal{K}}_{12}$ $\displaystyle=\mathbb{E}\left[(X_{i}-\mu)((X_{i}-\mu)^{2}-\sigma^{2})\right]= \mathbb{E}\left[(X_{i}-\mu)^{3}\right]=\sigma^{\frac{3}{2}}\gamma,$ $\displaystyle{\mathcal{K}}_{21}$ $\displaystyle={\mathcal{K}}_{12},\text{ and}$ $\displaystyle{\mathcal{K}}_{22}$ $\displaystyle=\mathbb{E}\left[((X_{i}-\mu)^{2}-\sigma^{2})((X_{i}-\mu)^{2}- \sigma^{2})\right]$ $\displaystyle=\mathbb{E}\left[(X_{i}-\mu)^{4}\right]-\sigma^{4}=\sigma^{4}( \kappa+2).$ $$

From the equality between ${\mathcal{K}}$ and ${\mathcal{I}}(\theta_{0})^{-1}$, we find that

$$ $\sigma^{4}(\kappa+2)=2\sigma^{4}\iff\kappa=0.$ $$

Hence, when normality is assumed for non-normally distributed data, the bias of the standard error estimate $\hat{\sigma}^{2}$ based on ${\mathcal{I}}^{-1}$ depends on $\kappa$, i.e., negative when $\kappa>0$ and positive when $\kappa<0$.
∎

###### Proof.

### B.2 Of 1

To prove that the NLL function of GGD is well defined for $\beta\in(0,2]$, we show that the PDF of GGD is everywhere positive with it being a positive definite function.
An easy but effective proof can be done by demonstrating that the PDF of GGD for $\beta\in(0,2]$ is equivalent to the characteristic function of an $\alpha$-stable distribution, up to a normalizing constant (Dytso et al., ).
Given the positive definiteness of all characteristic functions (Ushakov, ; Dytso et al., ), this equivalence assures the positive definiteness of the GGD PDF.
However, we offer a proof rooted in the properties of the positive definite function class.

###### Proof.

To demonstrate the positive definiteness of the GGD PDF, it is sufficient to show the positivity of the function:

$$ $f_{\beta}(x)=e^{-|x|^{\beta}},$ $$

where $\beta\in(0,2]$.

For a class of positive definite functions ${\mathcal{F}}$, which can be interpreted as Fourier transforms of bounded non-negative distributions, the function $f\in{\mathcal{F}}$, i.e.,

$$ $f(x)=\int^{\infty}_{\infty}e^{i\chi x}dV(\chi),$ $$

satisfy the following properties:

- 1.
For any non-negative scalars $a_{1},a_{2}$ and functions $f_{1},f_{2}\in{\mathcal{F}}$, $a_{1}f_{1}+a_{2}f_{2}\in{\mathcal{F}}$.
- 2.
For $f_{1},f_{2}\in{\mathcal{F}}$, $f_{1}f_{2}\in{\mathcal{F}}$.
- 3.
If a sequence of functions $f_{n}\in{\mathcal{F}}$ converges uniformly in every finite interval, then the limit function $\lim_{n\to\infty}f_{n}\in{\mathcal{F}}$.

Now, let $\rho=\beta/2$, and we aim to prove that $f_{\rho}=\exp(-|x|^{2\rho})$ belongs to ${\mathcal{F}}$ for $\rho\in(0,1)$, excluding the trivial case $\rho=1$.
Since for $0<\rho<1$,

$$ $|x|^{2\rho}=c_{\rho}\int^{\infty}_{0}\frac{\chi^{2\rho-1}d\chi}{1+(\frac{\chi} {x})^{2}},\quad c_{\rho}>0,$ $$

$f_{\rho}$ can be expressed as a uniform limit of functions $f_{\rho}\approx\lim_{n\to\infty}f_{n}$, where each $f_{n}$ takes the form:

$$ $\displaystyle f_{n}$ $\displaystyle=\exp{\left(-\sum^{n}_{\nu=1}\frac{a_{\nu}^{2}}{1+\left(\frac{b_{ \nu}}{x}\right)^{2}}\right)}$ $\displaystyle=\exp{\left(-\sum^{n}_{\nu=1}\frac{a_{\nu}^{2}x^{2}}{x^{2}+b_{\nu }^{2}}\right)}$ $\displaystyle=\prod_{\nu=1}^{n}{\exp{\left(-\frac{a_{\nu}^{2}x^{2}}{x^{2}+b_{ \nu}^{2}}\right)}},$ $$

for some sequences $\{a_{\nu}\}$ and $\{b_{\nu}\}$ with $\nu\in\{1,...,n\}$.

By simplifying, we find that for some $a\in\{a_{\nu}\}$ and $b\in\{b_{\nu}\}$ and letting $c^{2}=a^{2}b^{2}$:

$$ $\displaystyle\exp{\left(-\frac{a^{2}x^{2}}{x^{2}+b^{2}}\right)}$ $\displaystyle=\exp{\left(\frac{-a^{2}x^{2}-a^{2}b^{2}+a^{2}b^{2}}{x^{2}+b^{2}} \right)}$ $\displaystyle=\exp{\left(-\frac{a^{2}(x^{2}+b^{2})}{x^{2}+b^{2}}+\frac{a^{2}b^ {2}}{x^{2}+b^{2}}\right)}$ $\displaystyle=\exp{\left(-a^{2}\right)}\exp{\left(\frac{c^{2}}{x^{2}+b^{2}} \right)}.$ $$

Note that $\exp{\left(c^{2}/(x^{2}+b^{2})\right)}=\sum_{n=0}^{\infty}c^{2n}/n!\times(x^{2
}+b^{2})^{-n}$ from Taylor expansions.

From the properties of a positive definite function class, it is now sufficient to show that $(x^{2}+b^{2})^{-1}$ belongs to ${\mathcal{F}}$ to prove $f_{\beta}\in{\mathcal{F}}$.
We demonstrate this by expressing it as a Fourier transform of a bounded non-negative distribution:

$$ $(x^{2}+b^{2})^{-1}=\frac{1}{2b}\int^{\infty}_{-\infty}e^{i\chi x}e^{-b|\chi|}d\chi.$ $$

∎

###### Proof.

### B.3 Of 2

###### Proof.

For random variable $X\sim\text{GGD}(0,\alpha,\beta)$, its cumulative distribution function $F_{X}(t)$ is defined as

$$ $F_{X}(t)=\frac{1}{2}+\operatorname{sign}(t)\frac{\gamma\left(1/\beta,\left(|t- \mu|/\alpha\right)^{\beta}\right)}{2\Gamma\left(1/\beta\right)},$ $$

where

$$ $\Gamma(x)=\int_{0}^{\infty}t^{x-1}e^{-t}dt,\quad\text{and}\quad\gamma(x,s)= \int_{0}^{s}t^{x-1}e^{-t}dt,$ $$

represent the gamma function and lower incomplete gamma function for a complex parameter $x$ with a positive real part.

Expanding the left-hand side of equation 7, we obtain

$$ $\displaystyle\int_{-\infty}^{x}F_{X_{1}}(t)-F_{X_{2}}(t)dt$ $\displaystyle=\int_{-\infty}^{x}\operatorname{sign}(t)\left(\frac{\gamma\left( 1/\beta_{1},\left(|t-\mu|/\alpha\right)^{\beta}_{1}\right)}{2\Gamma\left(1/ \beta_{1}\right)}-\frac{\gamma\left(1/\beta_{2},\left(|t-\mu|/\alpha\right)^{ \beta}_{2}\right)}{2\Gamma\left(1/\beta_{2}\right)}\right)dt$ (12) $\displaystyle=\int_{x}^{\infty}\frac{\gamma\left(1/\beta_{2},\left(|t-\mu|/ \alpha\right)^{\beta}_{2}\right)}{2\Gamma\left(1/\beta_{2}\right)}-\frac{ \gamma\left(1/\beta_{1},\left(|t-\mu|/\alpha\right)^{\beta}_{1}\right)}{2 \Gamma\left(1/\beta_{1}\right)}dt.$ $$

Defining $f(\beta,t)=\gamma(1/\beta,\left(|t-\mu|/\alpha\right)^{\beta})/2\Gamma(1/\beta)$, we aim to demonstrate the monotonicity of $f$ to conclude the proof, as monotonically increasing $f$ leads to $f(\beta_{2},x)-f(\beta_{1},x)\geq 0$ for $\beta_{2}\geq\beta_{1}$, i.e., the integral in equation 12 is greater than or equal to zero.

With the definition of the gamma and lower incomplete gamma function,

$$ $f(\beta,t)=\frac{\gamma\left(1/\beta,\left(|t-\mu|/\alpha\right)^{\beta}\right )}{2\Gamma\left(1/\beta\right)}=\frac{\int_{0}^{(x/\alpha)^{\beta}}t^{1/\beta- 1}e^{-t}dt}{\int_{0}^{\infty}t^{1/\beta-1}e^{-t}dt}.$ (13) $$

Employing integration by substitution with $u=\left(\alpha^{\beta}t\right)^{\frac{1}{\beta}}$, equation 13 transforms to

$$ $\displaystyle f(\beta,t)$ $\displaystyle=\frac{\int_{0}^{x}\beta/\alpha\left((u/\alpha)^{\beta}\right)^{1 /\beta-1}e^{-(u/\alpha)^{\beta}}\left(\left((u/\alpha)^{\beta}\right)^{1/\beta -1}\right)^{-1}du}{\int_{0}^{\infty}\beta/\alpha\left((u/\alpha)^{\beta}\right )^{1/\beta-1}e^{-(u/\alpha)^{\beta}}\left(\left((u/\alpha)^{\beta}\right)^{1/ \beta-1}\right)^{-1}du}$ $\displaystyle=\frac{\int^{x}_{0}e^{-(u/\alpha)^{\beta}}du}{\int^{\infty}_{0}e^ {-(u/\alpha)^{\beta}}du}.$ $$

The function $f$ is thus increasing if, for $\beta_{1}\leq\beta_{1}$ and any $\alpha>0$,

$$ $\displaystyle\frac{\int^{x}_{0}e^{-(u/\alpha)^{\beta_{1}}}du}{\int^{\infty}_{0 }e^{-(u/\alpha)^{\beta_{1}}}du}\leq\frac{\int^{x}_{0}e^{-(u/\alpha)^{\beta_{2} }}du}{\int^{\infty}_{0}e^{-(u/\alpha)^{\beta_{2}}}du}$ $\displaystyle\iff\int^{x}_{0}\int^{\infty}_{0}e^{-\left((u/\alpha)^{\beta_{1}} +(v/\alpha)^{\beta_{2}}\right)}dvdu\leq\int^{x}_{0}\int^{\infty}_{0}e^{-\left( (u/\alpha)^{\beta_{1}}+(v/\alpha)^{\beta_{2}}\right)}dudv,$ $$

which follows from the monotonicity of the exponential function (Dytso et al., ).

Consequently, $f(\beta_{2},t)-f(\beta_{1},t)\geq 0$, implying that equation 12 is greater than or equal to zero.
∎

###### Proof.

### B.4 Of 2

###### Proof.

It is well known that the variance of the unbiased estimator $s^{2}$ is given by

$$ ${\mathbb{V}}[s^{2}]=\frac{1}{n}\left[\kappa-\frac{n-3}{n-1}\sigma^{2}\right].$ $$

And the MSE of a biased estimator $s_{\omega}^{2}=\omega(n-1)s^{2}$ is

$$ $\operatorname{MSE}(s_{\omega}^{2})=\omega^{2}(n-1)^{2}{\mathbb{V}}[s^{2}]+ \left[(n-1)\omega-1\right]^{2}\sigma^{4}.$ (14) $$

By differentiating equation 14 with respect to $\omega$, we can calculate the optimal $\omega$ value with minimal $\operatorname{MSE}(s_{\omega}^{2})$ as

$$ $\left.\frac{d\operatorname{MSE}(s_{\omega}^{2})}{d\omega}\right|_{\omega= \omega^{*}}=2(n-1)^{2}\omega^{*}{\mathbb{V}}[s^{2}]+2\left[(n-1)\omega^{*}-1 \right]\sigma^{4}=0.$ $$

Therefore,

$$ $\omega^{*}=\left[\frac{\sigma^{4}}{(n-1)\left({\mathbb{V}}[s^{2}]+\sigma^{4} \right)}\right]=\left[\frac{n-1}{n}\kappa+(n+1)\right]^{-1}.$ $$

It is easy to show that this is the optimal value, given that the second derivative is positive, i.e.,

$$ $\frac{d^{2}\operatorname{MSE}(s_{\omega}^{2})}{d\omega^{2}}=2(n-1)^{2}{\mathbb {V}}[s^{2}]+2(n-1)^{2}\sigma^{4}>0.$ $$

∎

###### Proof.

## Appendix C Experimental details

All experiments are conducted on a computational infrastructure consisting of 8 NVIDIA A100 80GB PCIe GPUs and 256 AMD EPYC™ 7742 processors.
Further details on the software setup can be accessed on GitHub(^1^11[https://github.com/ait-lab/ggtde](https://github.com/ait-lab/ggtde)).

##### Additional noise

To extend our study to discrete control scenarios, we enhance the baseline environments with stochastic perturbations.
In CartPole-v1 and MountainCar-v0, we introduce uniform noise to manipulate forces or torques.
Furthermore, in LunarLander-v2, wind dynamics are activated.
For a comprehensive analysis of how wind impacts the LunarLander-v2 dynamics, please consult the official documentation(^2^22[https://gymnasium.farama.org/environments/box2d/lunar_lander](https://gymnasium.farama.org/environments/box2d/lunar_lander)).

##### Implementation

To bolster numerical stability and enforce positivity constraints, we apply the softplus function, a smooth approximation to the ReLU function (Dugas et al., ), to the outputs from the variance or beta head.
Furthermore, we modify the NLL loss for the GGD by employing $Q^{\beta}$ as a multiplier rather than as an exponent:

$$ $\mathcal{L}_{\text{GGD-NLL}}\approx\sum_{t}{\left(\left|\delta_{t}\right|/Q^{ \alpha}_{t}\right)\times Q^{\beta}_{t}-\log{Q^{\beta}_{t}/Q^{\alpha}_{t}}+\log {\Gamma(1/Q^{\beta}_{t})}},$ $$

Although this adaptation deviates from the exact formulation of the NLL loss of the GGD, the difference in computed loss between the original and modified forms remains negligible across practical ranges of the TD error and $\beta$.
This adjustment offers a practical advantage by addressing the issue of flat regions in the original loss function, thus yielding more informative gradients for model updates.
Importantly, despite this modification, the positive-definiteness ensured by 1 is preserved, as the modified loss remains proportionally related to the original when $\beta$ is set to 1.
This property is crucial for maintaining stability and convergence in error modeling based on GGD.

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/loss-ant.png

Additionally, introducing scale invariance into the loss framework can be achieved by a simple adjustment: replacing the squared TD error term with the absolute error.
As the BIEV regularization framework is agnostic to the specific choice of loss function, the impact of the loss selection on the overall performance of BIEV regularization is minimal, as demonstrated in Figure 11.
With MAE loss, the resulting loss function depends only on the absolute TD errors, reducing the influence of the regularization temperature $\lambda$.
In contrast, for other loss functions, the performance tends to be more sensitive to the choice of temperature, motivating us to perform an ablation study on the effect of the regularization temperature, detailed in Section D.3.

## Appendix D Ablation studies

We conduct a series of ablation studies using the SAC algorithm across selected MuJoCo environments.

### D.1 On risk-averse weighting

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/raw-ant.png

The theoretical foundation of the risk-averse weighting $\omega^{\text{RA}}_{t}=Q^{\beta}_{t}$ is provided by 2.
Its empirical effectiveness, in comparison to the original GGD NLL and risk-seeking weighting $\omega^{\text{RA}}_{t}=1/Q^{\beta}_{t}$, is demonstrated in Figure 12.
It is evident that risk-aversion fosters sample-efficient training but does not necessarily translate to improved asymptotic performance.
Notably, the adoption of risk-seeking behavior results in periodic performance jumps, indicating active exploration during side-stepping.

### D.2 On BIEV regularization

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/biev-ant.png

The BIEV loss term is inherently independent of the parameter head, as it relies solely on the variance of TD errors.
This encourages us to investigate its efficacy independently of the GGD error modeling scheme.
As depicted in Figure 13, the impact of BIEV on sample efficiency is nearly equivalent to that of BIV regularization, but slightly more advantageous in terms of asymptotic performance.

### D.3 On regularizing temperature

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/temp-ant.png

The loss function presented in equation 4 is influenced by the scale of the rewards, requiring the regularizing temperature to be adjusted accordingly for different environments.
As shown in Figure 14, the regularizing temperature $\lambda$ affects sample efficiency, training stability, and asymptotic performance.
However, this sensitivity is notable only with large differences in scale, such as between 0.01 and 100.
Furthermore, as the regularizing temperature approaches zero, the model’s performance converges to that of the GGD error modeling scheme without BIEV regularization, labeled as ‘0’, as expected.

### D.4 On alpha head

Figure: (a) Ant-v4
Refer to caption: extracted/5894656/figures/alpha-ant.png

We exclusively utilize the beta head for computational efficiency and training stability, as elaborated in 2.
Integrating the alpha head, as anticipated, diminishes sample efficiency and even leads to decreased asymptotic performance, particularly evident in Ant-v4, as depicted in Figure 15.

## Appendix E Extension

Figure: (a) CartPole-v1
Refer to caption: extracted/5894656/figures/tde-dqn-cartpole.png

As discussed in Section 5, we expand our empirical analysis to include $Q$-learning, with particular attention to DQN Mnih et al. ().
Figure 16 demonstrates that TD errors generated by DQN exhibit deviations from a Gaussian distribution and closely follow the GGD.

Figure: (a) CartPole-v1
Refer to caption: extracted/5894656/figures/result-dqn-cart-pole.png

The performance of DQN in noisy discrete control environments, as depicted in Figure 17, highlights that simply applying our method to $Q$-learning is not sufficiently effective.
Interestingly, regularization appears to have minimal impact, as evidenced by the similar performance observed among DQN variants equipped with GGD error modeling.

Since our implementation of DQN utilizes a greedy action policy, the approximated value function inherently represents the policy.
Consequently, employing risk-averse weighting might impede exploration.
Conversely, in policy gradient algorithms featuring separate policy networks, such risk-aversion serves to mitigate noisy supervision, aligning with its intended purpose.
The adverse effects of risk-aversion are evident in Figure 18, where the model trained with the original NLL exhibits slightly superior sample efficiency compared to the proposed method.

Figure: (a) CartPole-v1
Refer to caption: extracted/5894656/figures/raw-cart-pole.png

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: - It seems like the distribution associated with environment transition stochasticity could easily lead to multi modal distributions of the cumulative return, which might be a much more interesting/important aspect than estimating the shape of a unimodal distribution better. Please address whether you observed any evidence of multimodality in your empirical results, and if so, how your method handles or could be extended to handle such cases.

- Rewards are typically bounded, so I would expect any estimation of the tail behavior to be quite sensitive to various practical assumptions. Please describe any preprocessing steps applied to the rewards, and how the performance varies with different reward scales or bounds.

- Given that your method requires predicting an extra head for the beta/kurtosis parameter, and training this (see Equation (4)) requires gradient descent through the gamma function of the output, how stable is the learning and/or optimization?  Please consider providing empirical evidence of the optimization stability, such as plots of the beta parameter estimates over time, or  discuss any specific techniques used to ensure stable training.
soundness: 3
strengths: - Effective uncertainty estimation is an important problem for RL and the authors study a valid relaxation of commonly held simplifying assumptions around Gaussianity.
- The authors conduct an experimental study that shows plausible benefits of the proposed generalizations in simple environments
summary: The authors consider uncertainty estimation in TD learning and seek to generalize the conceptually simple and implicit Gaussian assumption behind the MSE loss. To do so, they consider  the generalized gaussian distribution, which has an additional shape parameter that can modulate heavy/light tailedness for the errors. The authors extend a previous approach to estimating uncertainty in TD learning to include an additional network to predict the kurtosis/shape parameter.

Edit: updated my score after reviewing the author replies.
weaknesses: - There appear to be several ways to generalize the Normal distribution to have different tail behaviors, for example see q-Gaussian for yet another alternative. In the tradeoff between simplicity and more expressive modeling, it is unclear what is the right axis to explore. Even more generally, all of these are unimodal, symmetric (i.e. zero skew) distributions and arguably one could consider a full return distribution as well. In fact, prior work on distributional RL has explored this (https://arxiv.org/abs/1707.06887). 

- In terms of the learning algorithm, the paper modifies the loss function proposed in [Mai et al 2022] to the GGD case aided with an additional predictor for the kurtosis (the beta term). 

- The empirical results are mostly within error bars of the baselines. This seems especially so, when considering the marginal benefits in going from a basic uncertainty modeling with the gaussian assumption to the extra parameter estimation. Considering that this involves a whole extra network (not just an extra hyper-parameter) to predict the beta, this seems like a much more complex method for relatively little gain.

- The tailedness of the distribution seems like a very sensitive parameter to estimate, and sensitive to simple reward transformations and/or clipping so I would be surprised if these observations are robust to such changes.

---

confidence: 3
contribution: 3
initial_rating: 6
presentation: 2
questions: What specific tasks or environments in the real world do the authors envision as most beneficial for applying their method?
soundness: 3
strengths: 1. The paper offers a closed-form expression demonstrating the relationship between uncertainty and the generalized Gaussian distribution's shape parameter, adding depth to the theoretical framework.
2. The framework's applicability to both discrete and continuous control settings makes it relevant across various reinforcement learning contexts.
3. The emphasis on data-dependent noise and aleatoric uncertainty is timely and important for improving the robustness of reinforcement learning algorithms.
summary: The paper presents a novel framework for generalized Gaussian error modeling in uncertainty-aware temporal difference (TD) learning. It critiques conventional methods that assume a zero-mean Gaussian distribution for TD errors, leading to inaccurate uncertainty estimations. The proposed framework incorporates higher-order moments, specifically kurtosis, to enhance error modeling in reinforcement learning.
weaknesses: 1. The introduction of higher-order moments may complicate the implementation in practical scenarios, which could deter application by practitioners.
2. In Figure 5, the return curves of Ant, HalfCheetah, and Humanoid are still increasing at the end of training step. The training steps can be increased to compare the final performances when all algorithms are converged.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 1
questions: - Q1: Related to W1 and W2, what is the intuition behind Equation 3? In the original paper (Mai et al., 2022), this is Equation (10), which is preceded by a detailed explanation of many complex terms in the equation. As it stands, I do not think the main text is self-contained, as it requires substantial prior knowledge of that specific reference.

- Q2: How are ensembles used in this method? Ensembles are mentioned, but there is no explanation of how they are applied. Is this related to the variance in epistemic uncertainty estimation in the regularization term?

- Q3: How does the regularization term capture epistemic uncertainty (as the uncertainty that can be reduced through learning)?

- Q4: In Equation 4, how are risk-averse weights introduced? Was this derived from the inclusion of GGD in Mai et al. (2022), or was it manually designed to incorporate risk sensitivity?

- Q5: Why is estimating only beta sufficient to account for uncertainty compared to variance? Since the GGD has three parameters, with the first two set to 0 and 1 respectively, how can variance and kurtosis be represented simultaneously by a single number, beta? If I understand correctly, variance and kurtosis can be derived from beta (and alpha, in the case of variance). How do both quantities vary for a fixed alpha and a changing beta? There may be a range of TD-error distributions that cannot be fully described by just beta, which could potentially decrease the agent's performance.

- Q6: If Equation 4 holds, why does optimizing for both alpha and beta jointly decrease the agent’s performance, as shown in Figure 15?

- Q7: If considering higher moments of the TD-error distribution improves performance, does this mean that greater improvements are seen in tasks where TD-error tails are larger or smaller? How does performance compare to the evolution of the moments in the observed TD-errors? Does this method show a greater improvement over the baseline when higher-order moments are more or less pronounced?

- Q8: Typo in lines 198 and 987, and the font size of labels and ticks in figures could be increased significantly to facilitate reading.
soundness: 3
strengths: - (S1): The authors provide strong motivation for improving uncertainty-aware reinforcement learning (RL), it is an active area with relevance in applications like risk-averse decision-making, managing the exploration-exploitation trade-off, and enhancing sample efficiency.

- (S2): They conduct extensive experimentation, including ablation studies and parameter sweeps, to demonstrate that their proposed method outperforms the baseline (the original uncertainty-aware method).

- (S3): The proposed method is supported by theoretical validations and guarantees. Additionally, the proofs and mathematical analyses are well-documented in the appendix, making them relatively straightforward to follow.
summary: In this work, the authors present a method to estimate uncertainty based on the Generalized Gaussian Distribution (GGD) to better characterize temporal difference error (TD-error) distributions. Unlike previous work, which assumes Gaussianity, the GGD captures additional parameters, specifically kurtosis, to account for higher-order moments, thus enhancing the description of the TD-error distribution. The authors then modify an existing algorithm to operate under GGD assumptions, demonstrating improved performance across various settings. They also provide theoretical guarantees for well-behaved probability density function optimization under GGD assumptions and offer insights into why their proposed method is effective.
weaknesses: - (W1): The writing of the paper could be significantly improved. Many concepts and terms are used in the main text without proper definitions, making it difficult to follow.

- (W2): Some important aspects of the proposed method are not well-developed and are instead left to references from the baseline method. While this approach is generally acceptable, in this case, understanding key elements like the inclusion of GGD and the consideration of uncertainty during training requires substantial familiarity with the referenced material. As a result, the text is not self-contained, and combined with the clarity issues, makes the main content challenging to understand.

- (W3): Although extensive experimentation demonstrates the effectiveness of the proposal, some results are counterintuitive. There is also insufficient exposition of the proposed method’s inner workings and its differences from the baseline, for which performance plots alone are insufficient.

---

confidence: 4
contribution: 2
initial_rating: 3
presentation: 1
questions: 1. Please explain further why the Gaussian fits in Fig.2 are very off, and why looking at the TD errors of all examples can justify the underlying GGD assumption as discussed above.

2. How can Thm.2 lead to the specific weighting in (4)?

3. What is the main theoretical contribution of this work?
soundness: 3
strengths: - Discussion on some advantageous properties of GGD modelling
- Empirical success with the proposed objective
summary: Instead of Gaussian distribution, this paper applies generalized Gaussian distribution (GGD) to model temporal difference (TD) error, and proposes an uncertainty-aware objective function to minimize the TD error. It discussed several properties of the GGD distributions and empirically the proposed objective can work better than some baselines on common benchmark RL problems.
weaknesses: The main problem with the current paper is its presentation, especially the lack of motivation for the objective function.

1. Fig.2 requires further clarification. It is surprising to see that, for example, the Gaussian fitted PDF is so different from the empirical histogram for Ant-v4 (see also the MountainCar-v0 in Fig.7(d) for a more obvious mismatch). For Hopper-v4, the standard deviations of the two Gaussians are very close to each other despite the empirical histograms having very different shapes. More importantly, this figure is about the distribution of ALL TD errors across state-action pairs, while the NLL objective is about the underlying distribution of ONE single state-action pair. For this figure to make sense, one needs to assume that all TD errors come from the same underlying distribution, which doesn’t seem to be the case as Eq.(4) uses different betas for different step t. 

2. L273 offers a hypothesis, which is unverified. How do we define “unexpected states and rewards” and could you provide empirical support for this? Similarly, the explanation of interplay is unsatisfactory for the paragraph in L278. No evidence to support claims like “As training advances, aleatoric errors tend to become more influential while epistemic uncertainty diminishes, potentially resulting in non-normally distributed errors with heavier tails.”

3. It remains unclear how Thm.2 can lead to the specific weighting in Eq.(4). A detailed derivation would be helpful. Similarly, it is unclear why Eq.(8) can address the bias issue mentioned in L332. The current paper has too many hand-wavy arguments without sufficient rigorous discussions.

4. Finally, it seems that all theoretical properties in the current paper are stemmed from prior work, as shown by the reference in every theorem/proposition. If this is indeed the case, then what would be the main theoretical contribution of the current paper, other than an application of an existing tool (GGD) to TD learning?

Meta-review: This paper addresses the limitations of conventional uncertainty-aware Temporal Difference (TD) learning methods, and proposes a more sophisticated and flexible way to model uncertainty in TD learning by using the GGD and explicitly considering higher-order moments, leading to improved performance in policy gradient algorithms.

One of the major contribution from the paper, i.e., the loss function design, is modified from [Mai et al 2022] with the GGD case aided with an additional predictor for the kurtosis (the beta term). Mwanehile, all the theoretial properties are stemmed from prior work. Moreover, there are several alternative ways to generalize the Normal distribution to have different tail behaviors, especially the distributional RL. However, there is no empirical comparison conducted, which makes the significance of the paper difficult to be justified. 

In sum, this paper is promising, however, the current version is not ready to be published. I encourage the authors to take the reviewers' suggestions into account to improve the paper.

Decision: Reject