# A note on reducing spurious pressure oscillations in fully conservative
discontinuous Galerkin simulations of multicomponent flows

Eric J. Ching, Ryan F. Johnson, and Andrew D. Kercher  Laboratories for
Computational Physics and Fluid Dynamics, U.S. Naval Research Laboratory, 4555
Overlook Ave SW, Washington, DC 20375

###### keywords:

Discontinuous Galerkin method; pressure equilibrium; spurious pressure
oscillations; contact discontinuity; contact interface

††footnotetext:
DISTRIBUTION STATEMENT A. Approved for public release. Distribution is
unlimited.

##  1 Introduction

The compressible, inviscid, multicomponent Euler equations in $d$ spatial
dimensions are given by

$\frac{\partial y}{\partial t}+\nabla\cdot\mathcal{F}\left(y\right)=0$ (1.1)

where $t$ is the time,
$y(x,t):\mathbb{R}^{d}\times\mathbb{R}^{+}\rightarrow\mathbb{R}^{m}$ is the
vector of $m$ conservative variables, expanded as

$y=\left(\rho v_{1},\ldots,\rho v_{d},\rho
e_{t},C_{1},\ldots,C_{n_{s}}\right)^{T},$ (1.2)

and $\mathcal{F}(y):\mathbb{R}^{m}\rightarrow\mathbb{R}^{m\times d}$ is the
convective flux, the $k$th spatial component of which is written as

$\mathcal{F}_{k}\left(y\right)=\left(\rho v_{k}v_{1}+P\delta_{k1},\ldots,\rho
v_{k}v_{d}+P\delta_{kd},v_{k}\left(\rho
e_{t}+P\right),v_{k}C_{1},\ldots,v_{k}C_{n_{s}}\right)^{T}.$ (1.3)

$x=\left(x_{1},\ldots,x_{d}\right)$ is the vector of physical coordinates,
$\rho$ is the density, $v=\left(v_{1},\ldots,v_{d}\right)$ is the velocity
vector, $e_{t}$ is the mass-specific total energy,
$C=\left(C_{1},\ldots,C_{n_{s}}\right)$ is the vector of $n_{s}$ species
concentrations, and $P$ is the pressure, which is computed using the ideal gas
law as

$P=R^{0}T\sum_{i=1}^{n_{s}}C_{i},$ (1.4)

where $R^{0}$ is the universal gas constant and $T$ is the temperature. The
mass fraction of the $i$th species is computed as

$Y_{i}=\frac{\rho_{i}}{\rho},$

where $\rho_{i}=W_{i}C_{i}$ is the partial density, with $W_{i}$ denoting the
molecular weight of the $i$th species, and $\rho=\sum_{i=1}^{n_{s}}W_{i}C_{i}$
is the total density. In this work, we assume $d=1$ and mixtures of thermally
perfect gases, such that the specific heat capacities vary with temperature.

It is well-known that spurious pressure oscillations are typically generated
at contact interfaces when using fully conservative numerical schemes to solve
Equation (1.1) [1, 2, 3], an issue not observed when considering a
monocomponent calorically perfect gas. These pressure oscillations may lead to
solver failure and can still occur even if discrete entropy stability is
satisfied [4]. Quasi-conservative methods are usually employed to prevent the
occurrence of such oscillations. For example, the double-flux approach [5]
employs elementwise-constant, thermodynamically frozen auxiliary variables to
mimic a calorically perfect gas and mathematically guarantee preservation of
pressure equilibrium at contact interfaces.

A popular numerical method that has successfully been applied to
multicomponent flows (both nonreacting and reacting) in recent years is the
discontinuous Galerkin (DG) method [6, 7, 8]. This family of numerical schemes
boasts a number of desirable properties, such as arbitrarily high order of
accuracy on unstructured meshes, compatibility with local grid and polynomial
adaptation, and suitability for modern compute hardware. Initial efforts to
extend DG schemes to multicomponent flows employed the double-flux approach
[9, 10, 11]. This work, however, focuses on fully conservative DG schemes so
that any issues associated with not maintaining conservation are circumvented.
A handful of fully conservative DG schemes have been introduced that can
maintain pressure equilibrium at contact interfaces (while preserving order of
accuracy in smooth regions of the flow and without artificial dissipation) in
an approximate sense, i.e., pressure oscillations can occur but remain small
over long periods of time and do not cause the solver to diverge if the
solution is adequately resolved. For example, Johnson and Kercher [12]
computed a canonical test case involving the advection of a constant-pressure,
hydrogen/oxygen thermal bubble. They found that a colocated scheme
successfully maintained pressure equilibrium, but overintegration, which is
important for reducing aliasing errors, rapidly led to solver divergence. As
such, they developed an overintegration approach in which the flux is
evaluated with a modified state based on an approximate pressure that is
projected onto the selected finite element space via interpolation. This
approach was used to successfully compute complex multicomponent flows,
including a moving detonation wave and a reacting shear layer configuration
[12]. Additionally, it can maintain pressure equilibrium when using
multidimensional, curved elements [13] and can be easily incorporated into a
positivity-preserving, entropy-bounded framework [14, 13]. Bando [15] further
investigated this flux-evaluation approach and compared it to another approach
in which the pressure is projected onto the finite element space via
$L^{2}$-projection (as opposed to interpolation). He found that in the same
thermal-bubble test case, the interpolation-based approach is simpler and
yielded smaller deviations from pressure equilibrium although both approaches
(unlike standard flux evaluation) maintained stability over long times. He
also investigated the influences of initialization strategy and choice of
quadrature rule. A somewhat similar approach was proposed by Franchina et al.
[16], in which the primitive variables are treated as the unknowns and an
$L^{2}$-projection of the pressure is employed. Implicit time integration was
used in [16]; though compatible with explicit time stepping, a solution-
dependent Jacobian term must be inverted at each iteration, making it arguably
less appropriate for explicit time integration (which was used in [12] and
[15]).

In this short note, we revisit the flux-evaluation approaches considered in
[12] and [15], which successfully maintained pressure equilibrium in the
canonical hydrogen/oxygen thermal-bubble advection test. We consider high-
pressure, nitrogen/n-dodecane thermal-bubble advection at low and high
velocities. This type of test case has been previously used to assess
numerical schemes designed for simulating transcritical/supercritical, real-
fluid flows with, for example, cubic equations of state and more complicated
thermodynamic relations [17, 18]. Here, although we restrict ourselves to the
thermally perfect gas model, the increased nonlinearity of the thermodynamics
of the nitrogen/n-dodecane mixture is nevertheless perhaps more representative
of certain complicated mixtures at realistic conditions than the simpler
hydrogen/oxygen case. Furthermore, it makes the considered case more effective
at revealing deficiencies in numerical techniques. Also of interest is the
influence of projecting additional variables (other than solely pressure) to
the finite element trial space.

##  2 Mathematical formulation

###  2.1 Discontinuous Galerkin discretization

The semi-discrete form of Equation (1.1) is given as: find $y\in V_{h}^{p}$
such that

$\displaystyle\sum_{\kappa\in\mathcal{T}}\left(\frac{\partial y}{\partial
t},\mathfrak{v}\right)_{\kappa}-\sum_{\kappa\in\mathcal{T}}\left(\mathcal{F}\left(y\right),\nabla\mathfrak{v}\right)_{\kappa}+\sum_{\epsilon\in\mathcal{E}}\left(\mathcal{F}^{\dagger}\left(y^{+},y^{-},n\right),\left\llbracket\mathfrak{v}\right\rrbracket\right)_{\mathcal{E}}=0\qquad\forall\>\mathfrak{v}\in
V_{h}^{p},$ (2.1)

where $\left(\cdot,\cdot\right)$ denotes the inner product, $\mathcal{T}$ is
the set of cells $\kappa$, $\mathcal{E}$ is the set of interfaces $\epsilon$,
$n$ is the normal vector at a given interface, $y^{+}$ and $y^{-}$ are the
states at both sides of a given interface,
$\mathcal{F}^{\dagger}\left(y^{+},y^{-},n\right)$ is the numerical flux, and
$V_{h}^{p}$ is the space of basis and test functions

$\displaystyle V_{h}^{p}$ $\displaystyle=$
$\displaystyle\left\\{\mathfrak{v}\in\left[L^{2}\left(\Omega\right)\right]^{m}\,\middle|\,\forall\kappa\in\mathcal{T},\left.\mathfrak{v}\right|_{\kappa}\in\left[\mathsf{P}_{p}(\kappa)\right]^{m}\right\\},$
(2.2)

with $\mathsf{P}_{p}\left(\kappa\right)$ denoting the space of polynomial
functions of degree no greater than $p$ in $\kappa$. Only periodic boundary
conditions are considered in this work, such that all interfaces are interior
interfaces and the jump operator, $\left\llbracket\cdot\right\rrbracket$, is
defined as
$\left\llbracket\cdot\right\rrbracket=\left(\cdot\right)^{+}-\left(\cdot\right)^{-}$.
We employ a nodal basis using Gauss-Lobatto points, such that the element-
local solution approximation is given by

$y_{\kappa}=\sum_{j=1}^{n_{b}}y_{\kappa}(x_{j})\phi_{j},$ (2.3)

where $n_{b}$ is the number of basis functions,
$\left\\{\phi_{1},\ldots,\phi_{n_{b}}\right\\}$ is the set of basis functions,
and $\left\\{x_{1},\ldots,x_{n_{b}}\right\\}$ is the set of node coordinates.
The flux in Equation (2.1) is approximated as

$\mathcal{F_{\kappa}}\approx\sum_{k=1}^{n_{c}}\mathcal{F}\left(\left.\text{$\mathcal{P}$}\left(z\left(y_{\kappa}\right)\right)\right|_{x_{k}}\right)\varphi_{k},$
(2.4)

where $n_{c}\geq n_{b}$ (with $n_{c}>n_{b}$ corresponding to overintegration),
$\left\\{\varphi_{1},\ldots,\varphi_{n_{c}}\right\\}$ is the corresponding set
of basis functions (also based on Gauss-Lobatto points), $\mathcal{P}$ is a
projection operator, and
$z\left(y\right):\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ is a vector of
intermediate state variables. A major focus of this study is assessing the
effects of various choices of $\mathcal{P}$ and $z$, which will be detailed in
the next subsection. The integrals in Equation (2.1) are computed with a
quadrature-free approach [19, 20] that is also employed in [12]. The overall
trends observed here are expected to also apply to a more conventional
quadrature-based approach, which was employed in [15]; indeed, in the context
of maintaining pressure equilibrium in the hydrogen/oxygen thermal-bubble
configuration, both the quadrature-free and quadrature-based approaches
yielded very similar results.

###  2.2 Flux evaluation

This subsection discusses the choices of $\mathcal{P}$ and $z$ considered in
this study.

####  2.2.1 Projection operators

As in [15], three projection operators are considered:

  * 1.

$\mathcal{P}_{1}\left(\cdot\right)=\mathrm{id}\left(\cdot\right)$, the
identity function. This corresponds to a standard flux evaluation, such that
the RHS of Equation (2.4) reduces to
$\sum_{k=1}^{n_{c}}\mathcal{F}\left(y_{\kappa}\left(x_{k}\right)\right)\varphi_{k}$,
regardless of the choice of $z$.

  * 2.

$\mathcal{P}_{2}\left(\cdot\right)=\mathcal{I}\left(\cdot\right)$,
interpolatory projection onto $V_{h}^{p}$. Note that using a colocated scheme,
this is equivalent to the identity function. This approach was first
introduced by Johnson and Kercher [12].

  * 3.

$\mathcal{P}_{3}\left(\cdot\right)=\Pi\left(\cdot\right)$, $L^{2}$-projection
onto $V_{h}^{p}$.

In this study, non-identity projection is applied only when overintegration is
employed. Optimal convergence was previously observed when using interpolatory
projection, $\mathcal{P}_{2}$, [12, 15] and $L^{2}$-projection,
$\mathcal{P}_{3}$, [15]; therefore, we do not assess order of convergence
here.

####  2.2.2 Intermediate state variables

Three choices of $z$ are considered:

  * 1.

$z_{1}=\left(\rho v_{1},\ldots,\rho
v_{d},P,C_{1},\ldots,C_{n_{s}}\right)^{T}$, where total energy is replaced
with pressure, which was proposed by Johnson and Kercher [12]. Note that this
was the only choice of intermediate state variables considered in [12] and
[15].

  * 2.

$z_{2}=\left(v_{1},\ldots,v_{d},P,C_{1},\ldots,C_{n_{s}}\right)^{T}$, where
total energy is replaced with pressure and momentum is replaced with velocity.

  * 3.

$z_{3}=\left(v_{1},\ldots,v_{d},P,T,Y_{1},\ldots,Y_{n_{s}-1}\right)^{T}$,
which represents a full set of primitive variables.

####  2.2.3 Integration

We compute $p=2$ solutions using both colocated integration and
overintegration. The former was found to maintain stability over long times in
the hydrogen-oxygen thermal-bubble advection test [12, 15], whereas the
latter, in conjunction with $\mathcal{P}_{1}$ (i.e., standard flux
evaluation), rapidly led to solver divergence. For colocated integration, only
$\mathcal{P}_{1}$ is applied, which corresponds to standard colocated
integration. For overintegration, the 5-point Gauss-Lobatto nodal set is used
for the integration points. Note that Bando [15] assessed the effects of
different quadrature rules. Such an investigation is not repeated here.

##  3 Results

This section presents results for the advection of a high-pressure,
nitrogen/n-dodecane thermal bubble [17, 18]. Although this test case was
presented in [17] and [18] with discontinuous initial conditions,
discontinuity-capturing schemes were applied in those studies. Here, we are
specifically interested in the performance of DG schemes without any
discontinuity-capturing techniques (e.g., dissipative limiters or artificial
viscosity). Similar to the hydrogen/oxygen thermal-bubble configuration
considered in [10, 12, 15], we modify the initial condition to be smooth (but
still with high gradients) as

$\displaystyle Y_{n\text{-}\mathrm{C_{12}H_{26}}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1-\tanh\left(25|x|-0.2\right)\right],$
$\displaystyle Y_{\mathrm{N_{2}}}$ $\displaystyle=$ $\displaystyle
1-Y_{n\text{-}\mathrm{C_{12}H_{26}}},$ (3.1) $\displaystyle T$
$\displaystyle=$
$\displaystyle\frac{T_{\min}+T_{\max}}{2}-\frac{T_{\max}-T_{\min}}{2}\tanh\left(25|x|-0.2\right)\textrm{
K},$ $\displaystyle P$ $\displaystyle=$ $\displaystyle 6\textrm{ MPa},$

where $T_{\min}=363\text{ K}$ and $T_{\max}=900\text{ K}$. In the
hydrogen/oxygen case, despite the smooth initial condition, underresolution
leads to spurious pressure oscillations that, with standard overintegration,
grow rapidly and lead to solver divergence [12, 15]. We consider two advection
velocities, $v=1\text{ m/s}$, as in [17], and $v=600\text{ m/s}$, as in [18],
since we have found that changing the velocity can yield markedly different
results. The computational domain is $\left[-0.5,0.5\right]\text{ m}$,
partitioned into $50$ line cells. Both sides of the domain are periodic. The
HLLC [21] flux function is employed. All solutions are initialized using
interpolation and then integrated forward in time for ten advection periods
using the third-order strong-stability-preserving Runge-Kutta scheme [22] with
$\mathrm{CFL}=0.8$ based on the order-dependent linear-stability constraint.
To assess deviations from pressure equilibrium, we compute the following
global quantity as a function of time [15]:

$\Delta P(t)=\max_{x}P(t,x)-\min_{x}P(t,x),$ (3.2)

where the maximum and minimum are evaluated over all solution nodes. All
simulations are performed using a modified version of the JENRE® Multiphysics
Framework [23, 12] that incorporates the techniques described in this work. In
the following, $\tau$ denotes the time to advect the solution one period, and
$P_{0}$ is the initial pressure of $6\text{ MPa}$.

###  3.1 Low velocity, $v=1\text{ m/s}$

Figure 3.1a presents the temporal variation of $\Delta P$ for the following
approaches:

  * 1.

Colocated integration.

  * 2.

Standard overintegration, i.e., $\mathcal{P}_{1}$.

  * 3.

Modified overintegration using all combinations of
$\left\\{\mathcal{P}_{2},\mathcal{P}_{3}\right\\}$ and
$\left\\{z_{1},z_{2},z_{3}\right\\}$.

$\Delta P$ is sampled every $\tau/10$ seconds until the solution either
diverges or is advected for ten periods (i.e., $t=10\tau$). The solution
diverges rapidly in the case of standard overintegration or modified
integration with $\left\\{\mathcal{P}_{2},z_{3}\right\\}$. $\mathcal{P}_{2}$
in conjunction with $z_{2}$ or $z_{3}$ maintains solution stability for longer
times, although the solution nevertheless blows up well before $t=10\tau$.
Only colocated integration and $\mathcal{P}_{3}$, regardless of choice of $z$,
maintain stability through ten advection periods. Deviations from pressure
equilibrium are overall smallest with colocated integration. In the context of
$\mathcal{P}_{3}$, $z_{3}$ preserves pressure equilibrium slightly better than
$z_{1}$ and $z_{2}$, both of which give nearly identical results. At early
times, however, $\mathcal{P}_{2}$ yields the smallest pressure oscillations
(which then grow rapidly before leading to solver failure). Figure 3.1b
presents the temperature profiles for the stable solutions at $t=10\tau$.
Despite most effectively maintaining pressure equilibrium, colocated
integration leads to appreciable oscillations in the temperature profile that
are not observed in the $\mathcal{P}_{3}$ solutions, all of which exhibit good
agreement with the exact solution.

(a) Temporal variation of $\Delta P$. (b) Temperature profiles at $t=10\tau$.

Figure 3.1: Solution to the advection of a nitrogen/n-dodecane thermal bubble
at $v=1\text{ m/s}$. Figure 3.1a shows the temporal variation of $\Delta P$,
sampled at intervals of $\tau/10$ (for a total of 100 samples if the solution
does not diverge). For the $\mathcal{P}_{3}$ results, only every fourth marker
is shown. Figure 3.1b displays the temperature profiles at $t=10\tau$ for the
stable solutions.

###  3.2 High velocity, $v=600\text{ m/s}$

As in the previous subsection, Figure 3.2a displays the temporal variation of
$\Delta P$ for the considered approaches. The $\mathcal{P}_{2}$ solutions
again initially exhibit very small pressure oscillations before diverging,
with the $z_{3}$ case diverging most rapidly. Furthermore, when using
$\mathcal{P}_{3}$, $z_{3}$ again gives smaller deviations from pressure
equilibrium than $z_{1}$ and $z_{2}$. However, key differences with the low-
velocity case are observed. First, the magnitude of pressure oscillations is
generally greater in the high-velocity case, although the $\mathcal{P}_{3}$
solutions yield deviations of similar magnitude between both velocities.
Second, standard overintegration maintains stability throughout the
simulation. Additionally, of the stable solutions, the colocated case exhibits
the largest pressure deviations, which would likely continue to grow if the
simulation is run for longer times. In contrast, the maximum pressure
deviations for the other stable solutions seem to have either plateaued or
begun to plateau. Pressure equilibrium is the most well-preserved in the
$\mathcal{P}_{3}$ solutions, particularly in conjunction with $z_{3}$.

Figure 3.2b presents the temperature profiles of the stable solutions at
$t=10\tau$. Temperature oscillations are observed in the cases of colocated
integration and standard overintegration. In contrast, the $\mathcal{P}_{3}$
solutions exhibit better agreement with the exact temperature.

(a) Temporal variation of $\Delta P$. (b) Temperature profiles at $t=10\tau$.

Figure 3.2: Solution to the advection of a nitrogen/n-dodecane thermal bubble
at $v=600\text{ m/s}$. Figure 3.2a shows the temporal variation of $\Delta P$,
sampled at intervals of $\tau/10$ (for a total of 100 samples if the solution
does not diverge). For the $\mathcal{P}_{3}$ results, only every fourth marker
is shown. Figure 3.2b displays the temperature profiles at $t=10\tau$ for the
stable solutions.

To help explain why overintegration with $\mathcal{P}_{3}$ can lead to higher
solution stability than both standard overintegration (i.e.,
$\mathcal{P}_{1}$) and overintegration with $\mathcal{P}_{2}$, we take the
$\mathcal{P}_{1}$ solution at $t=10\tau$ and then evaluate the pressure in the
cell $\kappa=\left[-0.22,-0.2\right]\text{ m}$ with the $\mathcal{P}_{1}$,
$\mathcal{P}_{2}$, and $\mathcal{P}_{3}$ projection operators. The resulting
pressure-deviation profiles in the given cell are displayed in Figure 3.3.
Noticeable pressure deviations are located at the faces in the
$\mathcal{P}_{1}$ case (i.e., the pressure is evaluated normally); since
$\mathcal{P}_{2}$ (at least using Gauss-Lobatto points) does not modify
pressure at the faces, overall deviation from pressure equilibrium is not
appreciably reduced. In contrast, $\mathcal{P}_{3}$ markedly improves
preservation of pressure equilibrium throughout the cell, indicating that
$L^{2}$-projection can act as a mechanism to reduce spurious oscillations that
initially emerge via, for instance, underresolution or inexact evaluation of
the flux. In the case of $z_{2}$ or $z_{3}$, oscillations of other variables
can also be mitigated. However, it should be noted that $\mathcal{P}_{2}$ can
reduce pressure oscillations if deviations are large _inside_ the cell but not
at the faces. Furthermore, there is no guarantee that $\mathcal{P}_{3}$ will
always reduce spurious oscillations.

Figure 3.3: Pressure-deviation profiles in the cell
$\kappa=\left[-0.22,-0.2\right]\text{ m}$ evaluated with the
$\mathcal{P}_{1}$, $\mathcal{P}_{2}$, and $\mathcal{P}_{3}$ projection
operators for the $\mathcal{P}_{1}$ solution at $t=10\tau$. The dotted blue
line corresponds to exact pressure equilibrium.

##  4 Concluding remarks

In this short note, we reevaluated previously introduced techniques designed
to reduce spurious pressure oscillations at contact interfaces in
multicomponent flows in the context of overintegrated DG discretizations [12,
15]. Specifically, we focused on strategies that do not (a) introduce
conservation error, (b) rely on artificial viscosity or limiting, or (c)
degrade order of accuracy in smooth regions of the flow. The considered
techniques, which employ a projection of the pressure (and potentially
additional variables) onto the finite element trial space via either
interpolation or $L^{2}$-projection, were previously shown to maintain
approximate pressure equilibrium and thus stable solutions over long times (in
contrast to a standard DG scheme with overintegration) in the advection of a
hydrogen/oxygen thermal bubble. In [12], interpolation-based projection was
considered the preferred approach due to better preservation of pressure
equilibrium.

In this work, we considered a more challenging test case: advection of a high-
pressure nitrogen/n-dodecane thermal bubble at both low and high velocities.
All simulations were performed on a 50-cell grid using $p=2$. Key observations
from this test case are as follows:

  * 1.

Interpolation-based projection always led to solver divergence. Although it
overall maintained stability for longer times than standard overintegration in
the low-velocity case, it was outperformed by standard overintegration in the
high-velocity case.

  * 2.

$L^{2}$-projection always maintained solution stability. It produced larger
pressure deviations than colocated integration in the low-velocity case but
more effectively maintained pressure equilibrium than all other approaches in
the high-velocity case.

  * 3.

Performing $L^{2}$-projection of a full set of primitive variables
(specifically, mass fractions, pressure, temperature, and velocity) more
effectively preserved pressure equilibrium than if only pressure or if both
pressure and velocity were projected. However, projecting a smaller set of
variables may be sufficient for stability.

  * 4.

$L^{2}$-projection led to better predictions of temperature than colocated
integration and standard overintegration.

  * 5.

Colocated integration always maintained solution stability, but led to
inferior temperature predictions, indicating that higher resolution (than in
the case of overintegration) may be needed to offset the greater integration
error.

These findings suggest that this configuration is more effective at testing
the ability of a numerical scheme to preserve pressure equilibrium, even
without consideration of a cubic equation of state or more complicated
thermodynamic relations commonly associated with it [17, 18]. Furthermore, it
is also valuable to consider different advection velocities. It should be also
be noted that although the $L^{2}$-projection-based approach is the only
overintegration strategy that prevented solver divergence across all
considered conditions, the interpolation-based approach is simpler and was
shown to maintain robustness across a variety of challenging test cases
involving more realistic flow conditions and geometries, including moving
detonation waves and a chemically reacting shear layer [12], suggesting that
it may still be a reliable choice.

Future work will involve consideration of real-fluid effects, in particular a
cubic equation of state and thermodynamic departure functions [17], in
multiple dimensions. The ability of the considered techniques, specifically
the $L^{2}$-projection-based strategy, will be further assessed in this more
challenging context of transcritical and supercritical flows. Furthermore, a
fully conservative finite volume scheme that mathematically guarantees
preservation of pressure equilibrium was recently developed by Fujiwara et al.
[24]; an extension to DG schemes may indeed be worth pursuing.

## Acknowledgments

This work is sponsored by the Office of Naval Research through the Naval
Research Laboratory 6.1 Computational Physics Task Area.

## References

  * [1] R. Abgrall, Generalisation of the Roe scheme for the computation of mixture of perfect gases, La Recherche Aérospatiale 6 (1988) 31–43.
  * [2] S. Karni, Multicomponent flow calculations by a consistent primitive algorithm, Journal of Computational Physics 112 (1) (1994) 31 – 43.  doi:https://doi.org/10.1006/jcph.1994.1080.
  * [3] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations: A quasi conservative approach, Journal of Computational Physics 125 (1) (1996) 150 – 160.  doi:https://doi.org/10.1006/jcph.1996.0085.
  * [4] A. Gouasmi, K. Duraisamy, S. M. Murman, Formulation of entropy-stable schemes for the multicomponent compressible Euler equations, Computer Methods in Applied Mechanics and Engineering 363 (2020) 112912.
  * [5] R. Abgrall, S. Karni, Computations of compressible multifluids, Journal of Computational Physics 169 (2) (2001) 594 – 623.  doi:https://doi.org/10.1006/jcph.2000.6685.
  * [6] W. H. Reed, T. Hill, Triangular mesh methods for the neutron transport equation, Tech. rep., Los Alamos Scientific Lab., N. Mex.(USA) (1973).
  * [7] F. Bassi, S. Rebay, High-order accurate discontinuous finite element solution of the 2D Euler equations, Journal of Computational Physics 138 (2) (1997) 251–285.
  * [8] B. Cockburn, G. Karniadakis, C.-W. Shu, The development of discontinuous Galerkin methods, in: Discontinuous Galerkin Methods, Springer, 2000, pp. 3–50.
  * [9] G. Billet, J. Ryan, A Runge–Kutta discontinuous Galerkin approach to solve reactive flows: The hyperbolic operator, Journal of Computational Physics 230 (4) (2011) 1064 – 1083.  doi:https://doi.org/10.1016/j.jcp.2010.10.025.
  * [10] Y. Lv, M. Ihme, Discontinuous Galerkin method for multicomponent chemically reacting flows and combustion, Journal of Computational Physics 270 (2014) 105 – 137.  doi:https://doi.org/10.1016/j.jcp.2014.03.029.
  * [11] K. Bando, M. Sekachev, M. Ihme, Comparison of algorithms for simulating multi-component reacting flows using high-order discontinuous Galerkin methods (2020).  doi:10.2514/6.2020-1751.
  * [12] R. F. Johnson, A. D. Kercher, A conservative discontinuous Galerkin discretization for the chemically reacting Navier–Stokes equations, Journal of Computational Physics 423 (2020) 109826.  doi:10.1016/j.jcp.2020.109826.
  * [13] E. J. Ching, R. F. Johnson, A. D. Kercher, Positivity-preserving and entropy-bounded discontinuous Galerkin method for the chemically reacting, compressible Euler equations. Part II: The multidimensional case, arXiv preprint arXiv:2211.16297 https://arxiv.org/abs/2211.16297 (2022).
  * [14] E. J. Ching, R. F. Johnson, A. D. Kercher, Positivity-preserving and entropy-bounded discontinuous Galerkin method for the chemically reacting, compressible Euler equations. Part I: The one-dimensional case, arXiv preprint arXiv:2211.16254 https://arxiv.org/abs/2211.16254 (2022).
  * [15] K. Bando, Towards high-performance discontinuous Galerkin simulations of reacting flows using Legion, Ph.D. thesis, Stanford University (2023).
  * [16] N. Franchina, M. Savini, F. Bassi, Multicomponent gas flow computations by a discontinuous Galerkin scheme using L2-projection of perfect gas EOS, Journal of Computational Physics 315 (2016) 302–322.
  * [17] P. C. Ma, Y. Lv, M. Ihme, An entropy-stable hybrid scheme for simulations of transcritical real-fluid flows, Journal of Computational Physics 340 (2017) 330–357.
  * [18] B. Boyd, D. Jarrahbashi, A diffuse-interface method for reducing spurious pressure oscillations in multicomponent transcritical flow simulations, Computers & Fluids 222 (2021) 104924.
  * [19] H. Atkins, C. Shu, Quadrature-free implementation of discontinuous Galerkin methods for hyperbolic equations, ICASE Report 96-51, 1996, Tech. rep., NASA Langley Research Center, nASA-CR-201594 (August 1996).
  * [20] H. L. Atkins, C.-W. Shu, Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations, AIAA Journal 36 (5) (1998) 775–782.
  * [21] E. Toro, Riemann solvers and numerical methods for fluid dynamics: A practical introduction, Springer Science & Business Media, 2013.
  * [22] S. Gottlieb, C. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM review 43 (1) (2001) 89–112.
  * [23] A. Corrigan, A. Kercher, J. Liu, K. Kailasanath, Jet noise simulation using a higher-order discontinuous Galerkin method, in: 2018 AIAA SciTech Forum, 2018, AIAA-2018-1247.
  * [24] Y. Fujiwara, Y. Tamaki, S. Kawai, Fully conservative and pressure-equilibrium preserving scheme for compressible multi-component flows, Journal of Computational Physics 478 (2023) 111973.


* (22) A. Gadde and V. Krishna, _Comments on the bulk replica symmetry of multi-partite entanglement measures_ , Work in progress.
  * (23) X. Dong and A. Lewkowycz, _Entropy, Extremality, Euclidean Variations, and the Equations of Motion_ , _JHEP_ 01 (2018) 081 [1705.08453].
  * (24) B. Schumacher and M.A. Nielsen, _Quantum data processing and error correction_ , _Phys. Rev. A_ 54 (1996) 2629.
  * (25) X. Dong, D. Harlow and A.C. Wall, _Reconstruction of Bulk Operators within the Entanglement Wedge in Gauge-Gravity Duality_ , _Phys. Rev. Lett._ 117 (2016) 021601 [1601.05416].
  * (26) A. Almheiri, X. Dong and D. Harlow, _Bulk Locality and Quantum Error Correction in AdS/CFT_ , _JHEP_ 04 (2015) 163 [1411.7041].
  * (27) D. Harlow, _The Ryu–Takayanagi Formula from Quantum Error Correction_ , _Commun. Math. Phys._ 354 (2017) 865 [1607.03901].


# Deep Reinforcement Learning with Weighted Q-Learning

Andrea Cini  Carlo D’Eramo  Jan Peters  Cesare Alippi

###### Abstract

Overestimation of the maximum action-value is a well-known problem that
hinders Q-Learning performance, leading to suboptimal policies and unstable
learning. Among several Q-Learning variants proposed to address this issue,
Weighted Q-Learning (WQL) effectively reduces the bias and shows remarkable
results in stochastic environments. WQL uses a weighted sum of the estimated
action-values, where the weights correspond to the probability of each action-
value being the maximum; however, the computation of these probabilities is
only practical in the tabular setting. In this work, we provide the
methodological advances to benefit from the WQL properties in Deep
Reinforcement Learning (DRL), by using neural networks with Dropout
Variational Inference as an effective approximation of deep Gaussian
processes. In particular, we adopt the Concrete Dropout variant to obtain
calibrated estimates of epistemic uncertainty in DRL. We show that model
uncertainty in DRL can be useful not only for action selection, but also
action evaluation. We analyze how our novel Deep Weighted Q-Learning algorithm
reduces the bias w.r.t. relevant baselines and provide empirical evidence of
its advantages on several representative benchmarks.

Reinforcement Learning, Q-Learning, Deep Reinforcement Learning, DQN, ICML

##  1 Introduction

Reinforcement Learning (RL) aims at learning how to take optimal decisions in
unknown environments by solving credit assignment problems that extend in
time. In order to be sample efficient learners, agents are required to
constantly update their own beliefs about the world, about which actions are
good and which are not. Temporal difference (TD) (Sutton & Barto, 1998) and
off-policy learning are the constitutional elements of this kind of behavior.
TD allows agents to bootstrap their current knowledge to learn from a new
observation as soon as it is available. Off-policy learning gives the means
for exploration and enables experience replay (Lin, 1991). Q-Learning
(Watkins, 1989) implements both paradigms.

Algorithms based on Q-learning are, in fact, driving Deep Reinforcement
Learning (DRL) research towards solving complex problems and achieving super-
human performance on many of them (Mnih et al., 2015; Hessel et al., 2018).
Nonetheless, Q-Learning is known to be positively biased (Van Hasselt, 2010)
since it learns by using the maximum over the - noisy - bootstrapped TD
estimates. This overoptimism can be particularly harmful in stochastic
environments and when using function approximation (Thrun & Schwartz, 1993),
notably also in the case where the approximators are deep neural networks (Van
Hasselt et al., 2016). Systematic overestimation of the action-values coupled
with the inherently high variance of DRL methods can lead to incrementally
accumulate errors, causing the learning algorithm to diverge.

Among the possible solutions, the Double Q-Learning algorithm (Van Hasselt,
2010) and its DRL variant Double DQN (Van Hasselt et al., 2016) tackle the
overestimation problem by disentangling the choice of the target action and
its evaluation. The resulting estimator, while achieving superior performance
in many problems, is negatively biased (Van Hasselt, 2013). Underestimation,
in fact, can lead in some environments to lower performance and slower
convergence rates compared to standard Q-Learning (D’Eramo et al., 2016; Lan
et al., 2020). Overoptimism, in general, is not uniform over the state space
and may induce to overestimate the value of arbitrary bad actions, throwing
the agent completely off. The same holds true, symmetrically, for overly
pessimistic estimates that might undervalue a good course of action. Ideally,
we would like DRL agents to be aware of their own uncertainty about the
optimality of each action, and be able to exploit it to make more informed
estimations of the expected return. This is exactly what we achieve in this
work.

We exploit recent developments in Bayesian Deep Learning to model the
uncertainty of DRL agents using neural networks trained with dropout
variational inference (Kingma et al., 2015; Gal & Ghahramani, 2016). We
combine, in a novel way, the dropout uncertainty estimates with the Weighted
Q-Learning algorithm (D’Eramo et al., 2016), extending it to the DRL settings.
The proposed Deep Weighted Q-Learning algorithm, or Weighted DQN (WDQN),
leverages an approximated posterior distribution on Q-networks to reduce the
bias of deep Q-learning. WDQN bias is neither always positive, neither
negative, but depends on the state and the problem at hand. WDQN only requires
minor modifications to the baseline algorithm and its computational overhead
is negligible on specialized hardware.

The paper is organized as follows. In Section 2 we define the problem
settings, introducing key aspects of value-based RL. In Section 3 we analyze
in depth the problem of estimation biases in Q-Learning and sequential
decision making problems. Then, in Section 4, we first discuss how neural
networks trained with dropout can be used for Bayesian inference in RL and,
from that, we derive the WDQN algorithm. In Section 5 we empirically evaluate
the proposed method against relevant baselines on several benchmarks. Finally,
we provide an overview of related works in Section 6, and we draw our
conclusions and discuss future works in Section 7.

##  2 Preliminaries

A Markov Decision Process (MDP) is a tuple
$\langle\mathcal{S},\mathcal{A},\mathcal{P},\mathcal{R},\gamma\rangle$ where
$\mathcal{S}$ is a state space, $\mathcal{A}$ is an action space,
$\mathcal{P}:\mathcal{S}\times\mathcal{A}\rightarrow\mathcal{S}$ is a
Markovian transition function,
$\mathcal{R}:\mathcal{S}\times\mathcal{A}\rightarrow\mathbb{R}$ is a reward
function, and $\gamma\in[0,1]$ is a discount factor. A sequential decision
maker ought to estimate, for each state $s$, the optimal value $Q^{*}(s,a)$ of
each action $a$, i.e., the expected return obtained by taking action $a$ in
$s$ and following the optimal policy $\pi^{*}$ afterwards. We can write
$Q^{*}$ using the Bellman optimality equation (Bellman, 1954)

$Q^{*}(s,a)={}\\\
\mathbb{E}\left[r_{t+1}+\gamma\max_{a^{\prime}}Q^{*}(s_{t+1},a^{\prime})|s_{t}=s,a_{t}=a\right].$
(1)

#### (Deep) Q-Learning

A classical approach to solve finite MDPs is the Q-Learning algorithm
(Watkins, 1989), an off-policy value-based RL algorithm, based on TD. A
Q-Learning agent learns the optimal value function using the following update
rule:

$Q(s_{t},a_{t})\leftarrow
Q(s_{t},a_{t})+\alpha\left(y_{t}^{QL}-Q(s_{t},a_{t})\right),$ (2)

where $\alpha$ is the learning rate and, following the notation introduced by
(Van Hasselt et al., 2016),

$y_{t}^{QL}=r_{t}+\gamma\max_{a}Q(s_{t+1},a).$ (3)

The popular Deep Q-Network algorithm (DQN) (Mnih et al., 2015) is a variant of
Q-Learning designed to stabilize off-policy learning with deep neural networks
in highly dimensional state spaces. The two most relevant architectural
changes to standard Q-Learning introduced by DQN are the adoption of a replay
memory, to learn offline from past experience, and the use of a target
network, to reduce correlation between the current model estimate and the
bootstrapped target value.

In practice, DQN learns the Q-values online, using a neural network with
parameters $\bm{\theta}$, sampling the replay memory, and with a target
network whose parameters $\bm{\theta}^{-}$ are updated to match those of the
online model every $C$ steps. The model is trained to minimize the loss

$L(\bm{\theta})=\underset{{\langle s_{i},a_{i},r_{i},s_{i}^{\prime}\rangle\sim
m}}{\mathbb{E}}\left[\left(y_{i}^{DQN}-Q(s_{i},a_{i};\bm{\theta})\right)^{2}\right],$
(4)

where $m$ is a uniform distribution over the transitions stored in the replay
buffer and $y^{DQN}$ is defined as

$y^{DQN}_{i}=r_{i}+\gamma\max_{a}Q(s_{i}^{\prime},a;\bm{\theta}^{-}).$ (5)

#### Double DQN

Among the many studied improvements and extensions of the baseline DQN
algorithm (Wang et al., 2016; Schaul et al., 2016; Bellemare et al., 2017;
Hessel et al., 2018), Double DQN (DDQN) (Van Hasselt et al., 2016) reduces the
overestimation bias of DQN with a simple modification of the update rule. In
particular, DDQN uses the target network to decouple action selection and
evaluation, and estimates the target value as

$y^{DDQN}_{i}=r_{i}+\gamma
Q(s_{i}^{\prime},\operatorname*{argmax}_{a}Q\left(s_{i}^{\prime},a;\bm{\theta});\bm{\theta}^{-}\right).$
(6)

DDQN improves on DQN converging to a more accurate approximation of the value
function, while maintaining the same model complexity and adding a minimal
computational overhead.

##  3 Estimation biases in Q-Learning

Choosing a target value for the Q-Learning update rule can be seen as an
instance of the Maximum Expected Value (MEV) estimation problem for a set of
random variables, here the action-values $Q(s_{t+1},{}\cdot{})$. Q-Learning
uses the Maximum Estimator (ME) 111Details about the estimators considered in
this section are provided in the appendix. to estimate the maximum expected
return and exploits it for policy improvement. It is well-known that ME is a
positively biased estimator of MEV (Smith & Winkler, 2006). The divergent
behaviour that may occur in Q-Learning, then, may be explained by the
amplification over time effect on the action-value estimates caused by the
overestimation bias, which introduces a positive error at each update (Van
Hasselt, 2010). Double Q-Learning (Van Hasselt, 2010), on the other hand,
learns two value functions in parallel and uses an update scheme based on the
Double Estimator (DE). It is shown that DE is a negatively biased estimator of
MEV, which helps to avoid catastrophic overestimates of the Q-values. The DDQN
algorithm, introduced in Section 2, is the extension of Double Q-Learning to
the DRL settings. In DDQN, the target network is used as a proxy of the second
value function of Double Q-Learning to preserve sample and model complexity.

In practice, as shown by D’Eramo et al. (2016) and Lan et al. (2020), the
overestimation bias of Q-Learning is not always harmful, and may also be
convenient when the action-values are significantly different among each other
(e.g., deterministic environments with a short time horizon, or small action
spaces). Conversely, the underestimation of Double Q-Learning is effective
when all the action-values are very similar (e.g., highly stochastic
environments with a long or infinite time horizon, or large action spaces). In
fact, depending on the problem, both algorithms have properties that can be
detrimental for learning. Unfortunately, a prior knowledge about the
environment is not always available and, whenever it is, the problem may be
too complex to decide which estimator should be preferred. Given the above, it
is desirable to use a methodology which can robustly deal with heterogeneous
problems.

###  3.1 Weighted Q-Learning

D’Eramo et al. (2016) proposes the Weighted Q-Learning (WQL) algorithm, a
variant of Q-Learning based on the therein introduced Weighted Estimator (WE).
The WE estimates MEV as the weighted sum of the random variables sample means,
weighted according to their probability of corresponding to the maximum.
Intuitively, the amount of uncertainty, i.e., the entropy of the WE weights,
will depend on the nature of the problem, the number of samples and the
variance of the mean estimator (critical when using function approximation).
WE bias is bounded by the biases of ME and DE (D’Eramo et al., 2016).

The target value of WQL can be computed as

$y_{t}^{WQL}=r_{t}+\gamma\sum_{a\in\mathcal{A}}w_{a}^{s_{t+1}}Q(s_{t+1},a),$
(7)

where $w_{a}^{s_{t+1}}$ are the weights of the WE and correspond to the
probability of each action-value being the maximum:

$w^{s}_{a}=P\left(a=\operatorname*{argmax}_{a^{\prime}}Q(s,a^{\prime})\right).$
(8)

The update rule of WQL can be obtained replacing $y_{t}^{QL}$ with
$y_{t}^{WQL}$ in Equation 2. The weights of WQL are estimated in the tabular
setting assuming the sample means to be normally distributed.

WE has been studied also in the Batch RL settings, with continuous state and
action spaces, by using Gaussian Process regression (D’Eramo et al., 2017).

##  4 Deep Weighted Q-Learning

Algorithm 1 Weighted DQN

Input: Q-network parameters $\bm{\theta}$, dropout rates $p_{1},\dots,p_{L}$,
a policy $\pi$, replay memory $\mathcal{D}$

$\bm{\theta}^{-}\leftarrow\bm{\theta}$

Initialize memory $\mathcal{D}$.

for step $t=0,\dots$ do

Select an action $a_{i}$ according to some policy $\pi$ given the distribution
over action-value functions $Q(s,{}\cdot{};\bm{\theta},\bm{\omega})$

Execute $a_{t}$ and add $\langle s_{t},a_{t},r_{t},s_{t+i}\rangle$ to
$\mathcal{D}$

Sample a mini-batch of transitions $\\{\langle
s_{i},a_{i},r_{i},s_{i}^{\prime}\rangle,i=1,\dots,M\\}$ from $\mathcal{D}$

for $i=1,\dots,M$ do $\rhd$ COMMENTcan be done in parallel

Take $K$ samples from $Q(s_{i},{}\cdot{};\bm{\theta}^{-},\bm{\omega})$ by
performing $K$ stochastic forward passes

Use the samples to compute the WDQN weights (Eq.4.2) and targets (Eq.15)

end for

Perform a SGD step on the WDQN loss (Eq. 4.2)

Eventually update $\bm{\theta}^{-}$

end for

A natural way to extend the WQL algorithm to the DRL settings is to consider
the uncertainty over the model parameters using a Bayesian approach. Among the
possible solutions to estimate uncertainty, bootstrapping has been the most
successful in RL problems, with BootstrappedDQN (BDQN) (Osband et al., 2016,
2018) achieving impressive results in environments where exploration is
critical. On the other hand, using bootstrapping necessitates significant
modifications to the baseline DQN architecture and requires to train a model
for each sample of the approximate posterior distribution. This limits the
number of samples available considerably and is a major drawback in using BDQN
to approximate the WE weights. Using dropout, conversely, does not impact
model complexity and allows to compute the weights of the WE by using
infinitely many samples.

In the following we first introduce how neural networks trained with dropout
can be used for approximated Bayesian inference and discuss how this approach
has been used with success in RL problems. Then, we propose a novel approach
to exploit the uncertainty over the model parameters for action evaluation,
adapting the WE to the DRL settings. Finally we analyze a possible shortcoming
of the proposed method and identify a solution from the literature to address
it.

###  4.1 Bayesian inference with dropout

Dropout (Srivastava et al., 2014) is regularization technique used to train
large neural networks by randomly dropping units during learning. In recent
years, dropout has been analyzed from a Bayesian perspective (Kingma et al.,
2015; Gal & Ghahramani, 2016), and interpreted as a variational approximation
of a posterior distribution over the neural network parameters. In particular,
Gal & Ghahramani (2016) show how a neural network trained with dropout and
weight decay can bee seen as an approximation of a deep Gaussian process
(Damianou & Lawrence, 2013). The result is a theoretically grounded
interpretation of dropout and a class of Bayesian neural networks (BNNs) that
are cheap to train and can be queried to obtain uncertainty estimates. In
fact, a single stochastic forward pass through the BNN can be interpreted as
taking a sample from the model’s predictive distribution, while the predictive
mean can be computed as the average of multiple samples. This inference
technique is known as Monte Carlo (MC) dropout and can be efficiently
parallelized in modern GPUs.

A straightforward application of Bayesian models to RL is Thompson Sampling
(TS) (Thompson, 1933). TS is an exploration technique that aims at improving
the sample complexity of RL algorithms by selecting actions according to their
probability of being optimal given the current agent’s beliefs. A practical
way to use TS in Deep Reinforcement Learning is to take a single sample from a
Q-network trained with dropout and select the action that corresponds to the
maximum sampled action-value (Gal & Ghahramani, 2016). TS based on dropout
achieves superior data efficiency compared against naïve exploration
strategies, such as $\varepsilon$-greedy, both in sequential decision making
problems (Gal & Ghahramani, 2016; Stadie et al., 2015) and contextual bandits
(Riquelme et al., 2018; Collier & Llorens, 2018). Furthermore, dropout has
been successfully used in model-based RL, to estimate the agent’s uncertainty
over the environment dynamics (Gal et al., 2016; Kahn et al., 2017; Malik et
al., 2019).

Here we focus on the problem of action evaluation. We show how to use
approximate Bayesian inference to evaluate the WE estimator by introducing a
novel approach to exploit uncertainty estimates in DRL. Our method empirically
reduces Q-Learning bias, is grounded in theory and simple to implement.

###  4.2 Weighted DQN

Let $Q({}\cdot{},{}\cdot{};\bm{\theta},\bm{\omega})$ be a BNN with weights
$\bm{\theta}$ trained with a Gaussian prior and Dropout Variational Inference
to learn the optimal action-value function of a certain MDP. We indicate with
$\bm{\omega}$ the set of random variables that represents the dropout masks,
with $\bm{\omega}_{i}$ the $i$-th realization of the random variables and with
$\bm{\Omega}$ their joint distribution:

$\displaystyle\bm{\omega}=\\{\omega_{lk}:l=1,\dots,L,k=1,\dots,K_{l}\\},$ (9)
$\displaystyle\begin{aligned}
&\omega_{lk}\sim\text{Bernoulli}(p_{l}),&\bm{\omega}_{i}\sim\bm{\Omega}(p_{1},\dots,p_{L}),\end{aligned}$
(10)

where $L$ is the number of weight layers of the network and $K_{l}$ is the
number of units in layer $l$.

Consider a sample $q(s,a)$ the MDP return, obtained taking action $a$ in $s$
and following the optimal policy afterwards. Following the GP interpretation
of dropout of Gal & Ghahramani (2016), we can approximate the likelihood of
this observation as a Gaussian such that

$q(s,a)\sim\mathcal{N}\left(Q(s,a;\bm{\theta},\bm{\omega});\tau^{-1}\right),$
(11)

where $\tau$ is the model precision.

We can approximate the predictive mean of the process, and the expectation
over the posterior distribution of the Q-value estimates, as the average of
$T$ stochastic forward passes through the network:

$\mathbb{E}\left[Q(s,a;\bm{\theta},\bm{\omega})\right]\approx{}\\\
\hat{Q}_{T}(s,a;\bm{\theta})=\frac{1}{T}\sum_{t=1}^{T}Q(s,a;\bm{\theta},\bm{\omega}_{t}).$
(12)

$\hat{Q}_{T}(s,a;\bm{\theta})$ is the BNN prediction of the action-values
associated to state $s$ and action $a$. A similar, more computationally
efficient, approximation can be obtained through weight averaging, which
consists in scaling the output of the neurons in layer $l$ by $1-p_{l}$ during
training and leaving them unchanged at inference time. We indicate this
estimate as $\tilde{Q}(s,a;\bm{\theta})$ and we use it for action selection
during training.

The model epistemic uncertainty, i.e., the model uncertainty over its
parameters, can be measured similarly as the sample variance across $T$
realizations of the dropout random variables:

${Var}\left[Q(s,a;\bm{\theta},\bm{\omega})\right]\approx{}\\\
\frac{1}{T}\sum_{t=1}^{T}Q(s,a;\bm{\theta},\bm{\omega}_{t})^{2}-\left(\frac{1}{T}\sum_{t=1}^{T}Q(s,a;\bm{\theta},\bm{\omega}_{t})\right)^{2}.$
(13)

As shown in Gal & Ghahramani (2016) the predictive variance can be
approximated with the variance of the estimator in Eq. 13 plus the model
inverse precision $\tau^{-1}$.

We can estimate the probability required to calculate the WE in a similar way.
Given an action $a$, the probability that $a$ corresponds the maximum expected
action-value can be approximated as the number of times in which, given $T$
samples, the sampled action-value of $a$ is the maximum over the number of
samples

$\displaystyle w^{s}_{a}(\bm{\theta})$
$\displaystyle=P\left(a=\operatorname*{argmax}_{a^{\prime}}Q(s,a^{\prime};\bm{\theta},\bm{\omega})\right)$
$\displaystyle\approx\frac{1}{T}\sum_{t=1}^{T}\left\llbracket
a=\operatorname*{argmax}_{a^{\prime}}Q(s,a^{\prime};\bm{\theta},\bm{\omega}_{t})\right\rrbracket,$
(14)

where $\llbracket\dots\rrbracket$ are the Iverson brackets ($\llbracket
P\rrbracket$ is $1$ if $P$ is true, $0$ otherwise). The weights can be
efficiently inferred in parallel with no impact in computational time.

We can define the WE given the Bayesian target Q-network estimates using the
obtained weights as:

$y^{WDQN}_{i}=r_{i}+\gamma\sum_{a\in\mathcal{A}}w^{s^{\prime}_{i}}_{a}(\bm{\theta}^{-})\hat{Q}_{T}(s_{i}^{\prime},a;\bm{\theta}^{-}).$
(15)

Finally we report for completeness the loss minimized by WDQN, where the
parameter updates are backpropagated using the dropout masks:

$\displaystyle L(\bm{\theta})$ $\displaystyle={}$
$\displaystyle\quad\underset{{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\hfil\cr\langle
s_{i},a_{i},r_{i},s_{i}^{\prime}\rangle&\sim m\\\
\bm{\omega}_{i}&\sim\bm{\Omega}\crcr}}}}{\mathbb{E}}\left[\left(y_{i}^{WDQN}-Q(s_{i},a_{i};\bm{\theta},\bm{\omega_{i}})\right)^{2}\right]$
(18)
$\displaystyle\quad+\lambda\sum_{l=1}^{L}\left\|\bm{\theta}_{l}\right\|^{2}_{2},$
(19)

where $\bm{\theta}_{l}$ are the weights of layer $l$ and $\lambda$ the weight
decay coefficient. Using weight decay is necessary for the the variational
approximation being valid. The complete WDQN algorithm is reported in
Algorithm 1.

###  4.3 Concrete Dropout

The dropout probabilities are variational parameters and influence the quality
of the approximation. Ideally, they should be tuned to maximize the log-
likelihood of the observations using a validation method. This is clearly not
possible in RL where the available samples, and the underlying distribution
generating them, change as the policy improves. In fact, using dropout with a
fixed probability might lead to a poor uncertainty estimates (Osband et al.,
2016, 2018; Gal et al., 2017). Concrete Dropout (Gal et al., 2017) mitigates
this problem by using a differentiable continuous relaxation of the Bernoulli
distribution and learning the dropout rate from data.

In practice, this means that the distribution of the dropout random variables
$\omega_{lk}$ becomes

$\displaystyle\omega_{lk}=\sigma\left({\beta}\left(\log
p_{l}-\log(1-p_{l})+\log u-\log(1-p_{l})\right)\right),$ (20)

where $\beta$ is a temperature parameter (fixed at $\beta=10$), $u$ is a
uniform random variable $u\sim\mathcal{U}(0,1)$ and $\sigma({}\cdot{})$ is the
sigmoid function. With this formulation the sampling procedure becomes
differentiable and the loss in Eq. 4.2 can be rewritten as:

$\displaystyle L(\bm{\theta},\bm{p})$ $\displaystyle={}$
$\displaystyle\quad\underset{{\vbox{\Let@\restore@math@cr\default@tag\halign{\hfil$\m@th\scriptstyle#$&$\m@th\scriptstyle{}#$\hfil\cr\langle
s_{i},a_{i},r_{i},s_{i}^{\prime}\rangle&\sim m\\\
\bm{\omega}_{i}&\sim\bm{\Omega}\crcr}}}}{\mathbb{E}}\left[\left(y_{i}^{WDQN}-Q(s_{i},a_{i};\bm{\theta},\bm{\omega_{i}})\right)^{2}\right]$
(23)
$\displaystyle\quad+\sum_{l=1}^{L}\left(\lambda\left\|\bm{\theta}_{l}\right\|^{2}_{2}-\zeta
K_{l}\mathcal{H}(p_{l})\right),$ (24)

where $\zeta$ is a dropout regularization coefficient, $\mathcal{H}(p)$ is the
entropy of a Bernoulli distribution with parameter $p$ and $K_{l}$ the number
of neurons in layer $l$.

##  5 Experiments

Figure 1: Learning curve of the analyzed DQN variants on Asterix. The curves
on the left show the evaluation scores. The figure in the middle reports, for
each agent, the estimate of the expected return w.r.t. the starting screen of
the game; the dashed lines indicate the real discounted return for each agent.
The figure on the right shows the moving average of the entropy of the WDQN
weights. Each epoch corresponds to 1M frames. The results shown here are the
average of $3$ independent runs and the shaded areas represent $95\%$
confidence intervals. The curves are smoothed using a moving average of 10
epochs to improve readability.

In this section we compare WDQN against the standard DQN algorithm and its
DDQN variant to measure the effect of the different estimators on the quality
of the learned value functions and policies. First we run a proof of concept
experiment on the Lunar Lander environment (Brockman et al., 2016). Then we
perform an in depth analysis on Asterix, one of the Atari games in the Arcade
Learning Environment (ALE) (Bellemare et al., 2013) where DQN is known to
overestimate the action-values (Van Hasselt et al., 2016). Finally we test
WDQN in three environments of the MinAtar benchmark (Young & Tian, 2019). WDQN
achieves a more accurate estimation of the expected return which results, on
average, in better policies.

All the agents are evaluated after each epoch of training using the greedy
policy, for Asterix we follow the evaluation protocol of Mnih et al. (2015).
For each experiment we report both the average cumulative reward at each
evaluation step and prediction of the expected return, compared against the
real discounted return obtained by the agents. Even if TS is the natural
choice for WDQN, all the algorithms, unless explicitly stated, use an
$\varepsilon$-greedy policy to guarantee the fairness of the comparison. For
WDQN the greedy action is selected during training as the action corresponding
the maximum Q-value estimated with weight averaging, while during evaluation
we take the action with the highest probability of being optimal computed as
in Eq. 4.2. We found WDQN to be robust to the number of dropout samples used
to compute the WE (e.g., $T\geq 30$). We used low values for the weight decay
term, as in (Farebrother et al., 2018), and tuned the dropout regularization
coefficient for the problem at hand. The algorithms and the experimental setup
have been developed using the open-source RL libraries MushroomRL (D’Eramo et
al., 2020) and OpenAI Gym (Brockman et al., 2016). Full description of the
complete experimental setup and further results for each experiment can be
found in the appendix.

###  5.1 Lunar Lander

Figure 2: Learning curves on Lunar Lander. On the left the evaluation scores
of each epoch ($20$k steps). On the right the maximum action-value estimated
by each agent at the initial state of the MDP; the dashed lines are the real
obtained discounted return. The results shown are the average of 20
independent runs, shaded areas are 95% confidence intervals. The curves are
smoothed using a moving average of 10 epochs.

Lunar Lander is an MDP from the Gym collection. In order to solve the
environment, the agent has to control the thrusters of a spacecraft to safely
land at a specific location. In order to make prediction and control more
challenging, we increased the stochasticity of the environment adding a 10%
probability of repeating the last executed action, instead of the one selected
by the agent. The three agents use the same exact network architecture and
hyperparameters, the only difference is that WDQN uses Concrete Dropout in
each hidden layer. Figure 2 shows the learning curves of the three agents,
with WDQN achieving a significantly higher average reward and prediction
accuracy.

Figure 3: Comparison of WDQN against WDQN with TS and vanilla DQN with
Concrete Dropout (DropDQN). On the left the evaluation scores. On the right
the maximum estimated action-value at the starting screen of the game; dashed
lines are the real obtained discounted return. The results shown are the
average across 3 different random seeds, the shaded areas are 95% confidence
intervals. The curves are smoothed using a moving average of 10 epochs.

###  5.2 Asterix

Figure 4: Learning curves of the analyzed DQN variants on three Minatar games.
In the leftmost column, we show the evaluation scores after each of the 200
training epochs, where $1$ epoch corresponds to $25$k steps. In the middle
column we report for each game and agent the estimate of the expected return
w.r.t. the initial state of the environment; the dashed lines indicate the
real discounted return obtained by the agents. In the rightmost column we show
the moving average of the entropy of the WDQN weights. The results shown here
are averaged over 20 independent runs and the shaded areas represent $95\%$
confidence intervals. The curves are smoothed using a moving average of 5
epochs.

We use the same neural network and hyperparameters of (Mnih et al., 2015),
except for the optimizer that we replace with Adam (Kingma & Ba, 2015),
following more recent best practices (Hessel et al., 2018; Dabney et al.,
2018). For WDQN we use Concrete Dropout only in the fully connected layer
after the convolutional block. We use sticky-actions, with a probability of
repeating last action of $25\%$, as proposed by (Machado et al., 2018). Figure
1 shows the result of the comparison in terms of average reward and prediction
accuracy. Differently from what observed by (Van Hasselt et al., 2016), DQN
does not diverge in our setting. In fact, while widely overestimating the
discounted return, DQN achieves better sample-complexity and higher cumulative
reward w.r.t. DDQN.

WDQN, on the other hand, is clearly the less biased option and despite the
slow start - that can be explained by the regularization introduced by dropout
- reaches the highest average reward. After an initial phase, the average
entropy of the WE weights remains almost constant with two possible
interpretations: 1) in many states there is not a single action that dominates
the others, 2) the agent is not able to completely resolve its uncertainty
over the action-value functions.

In order to gather more insights on WDQN and on how each component (namely the
policy and the MEV estimator) influence the observed behavior, in Figure 3 we
compare it against a baseline DQN agent trained with Concrete Dropout (DropDQN
in the figure) and a version of WDQN that uses a TS policy as in (Gal &
Ghahramani, 2016). In this setting TS is not beneficial and underperforms.
DropDQN is clearly the worst performer of the tested algorithms, indicating
that the use of WE is indeed beneficial and suggesting that WDQN may improve
with ad hoc hyperparameter tuning.

###  5.3 MinAtar

MinAtar (Young & Tian, 2019) is a RL testbed with environments mimicking the
dynamics of games from ALE, but with a simplified state representation.
MinAtar implements also sticky actions and difficulty ramping (Machado et al.,
2018). The results of the experiment are shown in Figure 4. We use the same
convolutional neural network and hyperparameters of (Young & Tian, 2019), but
we replace RMSProp with the Adam optimizer 222We found Adam to provide
significantly more stable learning across all the agents. with learning rate
$1\mathrm{e}{-4}$. For WDQN we choose the dropout regularization coefficient
with random search and use the same value in the three games.

WDQN achieves higher average reward in Breakout, performs on par with the
other two algorithms in Seaquest and is more sample efficient in Freeway. For
what concerns the estimation of the expected return, WDQN shows a much lower
bias in the first epochs of the learning procedure and converge to more
accurate estimates in Breakout and Freeway, while performing similarly to the
other algorithms in Seaquest. Since the number of actions and the
regularization coefficient are the same among the three environments, it is
clear that the entropy of the WE weights heavily depends on the dynamics and
input space of the problem.

##  6 Related works

The Weighted Q-Learning algorithm (D’Eramo et al., 2016) joins previous works
in the attempt to address the bias of Q-Learning. Double Q-Learning (Van
Hasselt, 2010) firstly showed how the use of DE, to underestimate the action-
values, can be helpful to stabilize learning in highly stochastic
environments. Bias-corrected Q-Learning (Lee et al., 2013) improves the
learning stability of Q-Learning by subtracting a quantity to the target value
that depends on the standard deviation of the reward. The Double Weighted
Q-Learning algorithm (Zhang et al., 2017) (which, despite the name, is not
related to WQL) uses instead a combination of ME and DE to balance between the
overestimation and underestimation of the two estimators. DE and WE have also
been introduced in Batch RL as alternatives to the ME used in Fitted
Q-Iteration (Ernst et al., 2005), respectively in the Double Fitted
Q-Iteration and Weighted Fitted Q-Iteration algorithms (D’Eramo et al., 2017).
Among the three, Weighted Fitted Q-Iteration is the only algorithm able to
handle continuous action spaces.

Overestimation of the action-values can be even more problematic in the DQN
algorithm (Mnih et al., 2015), due to the high variance typical of DRL
approaches. The Double DQN algorithm (Van Hasselt et al., 2016) introduces the
use of DE in DQN, and shows better estimate of action-values and superior
performance w.r.t. vanilla DQN. Other variants of DQN addressing the
overestimation problem are based on exploiting multiple estimates of the
action-values. For instance, Averaged DQN (Anschel et al., 2017) controls the
variance of the estimation by averaging the target of the update over an
arbitrary number of previous checkpoints of the target network. Then, the
recent Maxmin DQN (Lan et al., 2020) reduces the bias keeping several
estimates of the action-values in parallel, and computing the target values
using the maximum of the minimums of each action-value estimate.
Overestimation is also detrimental in actor-critic DRL algorithms, e.g. Deep
Deterministic Policy Gradient (DDPG) (Lillicrap et al., 2015). Notably, the
Twin Delayed DDPG (TD3) algorithm (Fujimoto et al., 2018) is a variant of DDPG
exploiting several tricks, e.g. a clipped version of Double Q-Learning update,
that significantly improves the approximation of the value function and the
overall performance.

##  7 Conclusion and future works

We present WDQN, a new value-based Deep Reinforcement Learning algorithm, that
extends the Weighted Q-Learning algorithm to work in environments with an
highly dimensional state representation. WDQN is a principled and robust
method to exploit uncertainty in DRL to accurately estimate the maximum
action-value in a given state. We empirically support our claims by showing
that WDQN consistently reduces the bias of DQN across different tasks. Our
results corroborate the findings of previous works (Gal & Ghahramani, 2016;
Gal et al., 2017), confirming that dropout can be used successfully for
approximate Bayesian inference in DRL. Future works may explore the
combination of WDQN with the other, orthogonal, DQN extensions and may attempt
to adapt the WDQN approach to other techniques modelling uncertainty in deep
neural networks.

## References

  * Anschel et al. (2017) Anschel, O., Baram, N., and Shimkin, N.  Averaged-dqn: Variance reduction and stabilization for deep reinforcement learning.  In _Proceedings of the 34th International Conference on Machine Learning-Volume 70_ , pp. 176–185. JMLR. org, 2017.
  * Bellemare et al. (2013) Bellemare, M. G., Naddaf, Y., Veness, J., and Bowling, M.  The arcade learning environment: An evaluation platform for general agents.  _Journal of Artificial Intelligence Research_ , 47:253–279, jun 2013.
  * Bellemare et al. (2017) Bellemare, M. G., Dabney, W., and Munos, R.  A distributional perspective on reinforcement learning.  In _Proceedings of the 34th International Conference on Machine Learning - Volume 70_ , ICML’17, pp. 449–458. JMLR.org, 2017.
  * Bellman (1954) Bellman, R.  The theory of dynamic programming.  _Bulletin of the American Mathematical Society_ , 60(6):503–515, 1954.
  * Blumenthal & Cohen (1968) Blumenthal, S. and Cohen, A.  Estimation of the larger of two normal means.  _Journal of the American Statistical Association_ , 63(323):861–876, 1968.
  * Brockman et al. (2016) Brockman, G., Cheung, V., Pettersson, L., Schneider, J., Schulman, J., Tang, J., and Zaremba, W.  Openai gym, 2016.
  * Collier & Llorens (2018) Collier, M. and Llorens, H. U.  Deep contextual multi-armed bandits.  _CoRR_ , abs/1807.09809, 2018.
  * Dabney et al. (2018) Dabney, W., Rowland, M., Bellemare, M., and Munos, R.  Distributional reinforcement learning with quantile regression.  In _AAAI Conference on Artificial Intelligence_ , 2018.
  * Damianou & Lawrence (2013) Damianou, A. and Lawrence, N.  Deep gaussian processes.  In _Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics_ , volume 31 of _Proceedings of Machine Learning Research_ , pp. 207–215, Scottsdale, Arizona, USA, 29 Apr–01 May 2013. PMLR.
  * D’Eramo et al. (2017) D’Eramo, C., Nuara, A., Pirotta, M., and Restelli, M.  Estimating the maximum expected value in continuous reinforcement learning problems.  In _Thirty-First AAAI Conference on Artificial Intelligence_ , 2017\.
  * D’Eramo et al. (2020) D’Eramo, C., Tateo, D., Bonarini, A., Restelli, M., and Peters, J.  Mushroomrl: Simplifying reinforcement learning research.  _arXiv preprint arXiv:2001.01102_ , 2020.
  * D’Eramo et al. (2016) D’Eramo, C., Restelli, M., and Nuara, A.  Estimating maximum expected value through gaussian approximation.  In _International Conference on Machine Learning_ , pp. 1032–1040, 2016.
  * Ernst et al. (2005) Ernst, D., Geurts, P., and Wehenkel, L.  Tree-based batch mode reinforcement learning.  _Journal of Machine Learning Research_ , 6(Apr):503–556, 2005.
  * Farebrother et al. (2018) Farebrother, J., Machado, M. C., and Bowling, M.  Generalization and regularization in dqn, 2018.
  * Fujimoto et al. (2018) Fujimoto, S., Hoof, H., and Meger, D.  Addressing function approximation error in actor-critic methods.  In _International Conference on Machine Learning_ , pp. 1587–1596, 2018.
  * Gal & Ghahramani (2016) Gal, Y. and Ghahramani, Z.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning.  In _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48 of _Proceedings of Machine Learning Research_ , pp. 1050–1059, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * Gal et al. (2016) Gal, Y., McAllister, R. T., and Rasmussen, C. E.  Improving pilco with bayesian neural network dynamics models.  In _Data-Efficient Machine Learning workshop, ICML 2016_ , 2016.
  * Gal et al. (2017) Gal, Y., Hron, J., and Kendall, A.  Concrete dropout.  In _Advances in Neural Information Processing Systems 30_ , pp. 3581–3590. Curran Associates, Inc., 2017.
  * Hessel et al. (2018) Hessel, M., Modayil, J., van Hasselt, H., Schaul, T., Ostrovski, G., Dabney, W., Horgan, D., Piot, B., Azar, M. G., and Silver, D.  Rainbow: Combining improvements in deep reinforcement learning.  In _The Thirty-Second AAAI Conference on Artificial Intelligence (AAAI-18)_ , 2018.
  * Kahn et al. (2017) Kahn, G., Villaflor, A., Pong, V., Abbeel, P., and Levine, S.  Uncertainty-aware reinforcement learning for collision avoidance.  _CoRR_ , abs/1702.01182, 2017.
  * Kingma & Ba (2015) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  In _Proceedings of the 3rd International Conference on Learning Representations (ICLR)_ , 2015.
  * Kingma et al. (2015) Kingma, D. P., Salimans, T., and Welling, M.  Variational dropout and the local reparameterization trick.  In _Advances in Neural Information Processing Systems 28_ , pp. 2575–2583. Curran Associates, Inc., 2015.
  * Lan et al. (2020) Lan, Q., Pan, Y., Fyshe, A., and White, M.  Maxmin q-learning: Controlling the estimation bias of q-learning.  In _International Conference on Learning Representations_ , 2020.
  * Lee et al. (2013) Lee, D., Defourny, B., and Powell, W. B.  Bias-corrected q-learning to control max-operator bias in q-learning.  In _2013 IEEE Symposium on Adaptive Dynamic Programming and Reinforcement Learning (ADPRL)_ , pp. 93–99. IEEE, 2013.
  * Lillicrap et al. (2015) Lillicrap, T. P., Hunt, J. J., Pritzel, A., Heess, N., Erez, T., Tassa, Y., Silver, D., and Wierstra, D.  Continuous control with deep reinforcement learning.  In _International Conference on Learning Representations_ , 2015.
  * Lin (1991) Lin, L.-J.  Programming robots using reinforcement learning and teaching.  In _Proceedings of the Ninth National Conference on Artificial Intelligence - Volume 2_ , AAAI’91, pp. 781–786. AAAI Press, 1991.  ISBN 0262510596.
  * Machado et al. (2018) Machado, M. C., Bellemare, M. G., Talvitie, E., Veness, J., Hausknecht, M. J., and Bowling, M.  Revisiting the arcade learning environment: Evaluation protocols and open problems for general agents.  _Journal of Artificial Intelligence Research_ , 61:523–562, 2018.
  * Malik et al. (2019) Malik, A., Kuleshov, V., Song, J., Nemer, D., Seymour, H., and Ermon, S.  Calibrated model-based deep reinforcement learning.  In _Proceedings of the 36th International Conference on Machine Learning_ , volume 97 of _Proceedings of Machine Learning Research_ , pp. 4314–4323, Long Beach, California, USA, 09–15 Jun 2019. PMLR.
  * Mnih et al. (2015) Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A. A., Veness, J., Bellemare, M. G., Graves, A., Riedmiller, M., Fidjeland, A. K., and Ostrovski, G.  Human-level control through deep reinforcement learning.  _Nature_ , 518(7540):529, 2015.
  * Osband et al. (2016) Osband, I., Blundell, C., Pritzel, A., and Van Roy, B.  Deep exploration via bootstrapped dqn.  In _Advances in Neural Information Processing Systems_ , pp. 4026–4034, 2016.
  * Osband et al. (2018) Osband, I., Aslanides, J., and Cassirer, A.  Randomized prior functions for deep reinforcement learning.  In _Advances in Neural Information Processing Systems 31_ , pp. 8617–8629. Curran Associates, Inc., 2018.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S.  Pytorch: An imperative style, high-performance deep learning library.  In _Advances in Neural Information Processing Systems 32_ , pp. 8024–8035. Curran Associates, Inc., 2019.
  * Riquelme et al. (2018) Riquelme, C., Tucker, G., and Snoek, J.  Deep bayesian bandits showdown: An empirical comparison of bayesian deep networks for thompson sampling.  In _International Conference on Learning Representations_ , 2018.
  * Schaul et al. (2016) Schaul, T., Quan, J., Antonoglou, I., and Silver, D.  Prioritized experience replay.  In _International Conference on Learning Representations_ , Puerto Rico, 2016.
  * Smith & Winkler (2006) Smith, J. E. and Winkler, R. L.  The optimizer’s curse: Skepticism and postdecision surprise in decision analysis.  _Management Science_ , 52(3):311–322, 2006.
  * Srivastava et al. (2014) Srivastava, N., Hinton, G., Krizhevsky, A., Sutskever, I., and Salakhutdinov, R.  Dropout: A simple way to prevent neural networks from overfitting.  _Journal of Machine Learning Research_ , 15:1929–1958, 2014\.
  * Stadie et al. (2015) Stadie, B. C., Levine, S., and Abbeel, P.  Incentivizing exploration in reinforcement learning with deep predictive models.  _CoRR_ , abs/1507.00814, 2015.
  * Stone (1974) Stone, M.  Cross-validatory choice and assessment of statistical predictions.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 36(2):111–133, 1974.
  * Sutton & Barto (1998) Sutton, R. S. and Barto, A. G.  _Reinforcement learning: An introduction_ , volume 1.  MIT press Cambridge, 1998.
  * Thompson (1933) Thompson, W. R.  On the likelihood that one unknown probability exceeds another in view of the evidence of two samples.  _Biometrika_ , 25(3/4):285–294, 1933.
  * Thrun & Schwartz (1993) Thrun, S. and Schwartz, A.  Issues in using function approximation for reinforcement learning.  In _Proceedings of the 1993 Connectionist Models Summer School_ , pp. 255–263. Lawrence Erlbaum, 1993.
  * Van Hasselt (2010) Van Hasselt, H.  Double q-learning.  In _Advances in Neural Information Processing Systems_ , 2010.
  * Van Hasselt (2013) Van Hasselt, H.  Estimating the maximum expected value: an analysis of (nested) cross-validation and the maximum sample average.  _arXiv preprint arXiv:1302.7175_ , 2013.
  * Van Hasselt et al. (2016) Van Hasselt, H., Guez, A., and Silver, D.  Deep reinforcement learning with double q-learning.  In _AAAI_ , volume 16, pp. 2094–2100, 2016.
  * Wang et al. (2016) Wang, Z., Schaul, T., Hessel, M., Hasselt, H., Lanctot, M., and Freitas, N.  Dueling network architectures for deep reinforcement learning.  In _Proceedings of The 33rd International Conference on Machine Learning_ , volume 48 of _Proceedings of Machine Learning Research_ , pp. 1995–2003, New York, New York, USA, 20–22 Jun 2016. PMLR.
  * Watkins (1989) Watkins, C. J. C. H.  _Learning from Delayed Rewards_.  PhD thesis, King’s College, Cambridge, UK, May 1989.
  * Young & Tian (2019) Young, K. and Tian, T.  Minatar: An atari-inspired testbed for thorough and reproducible reinforcement learning experiments.  _arXiv preprint arXiv:1903.03176_ , 2019.
  * Zhang et al. (2017) Zhang, Z., Pan, Z., and Kochenderfer, M. J.  Weighted double q-learning.  In _Proceedings of the Twenty-Sixth International Joint Conference on Artificial Intelligence, IJCAI-17_ , pp. 3455–3461, 2017.  doi: 10.24963/ijcai.2017/483.

## Appendix

The appendix is organized as follows:

  * •

in Section A we further discuss the problem of estimating the maximum expected
value (MEV) of a set of random variables;

  * •

in Section B we provide, together with some additional empirical results,
details on the hyperparameters and experimental setup used to evaluate the
agents.

##  Appendix A Estimating the Maximum Expected Value

Let $X=\\{X_{1},\dots,X_{M}\\}$ be a set of $M\geq 2$ independent random
variables with unknown means $\mu_{i}=\mathbb{E}[X_{i}]$, and
$S=\bigcup_{i=1}^{M}S_{i}$ be a set of noisy samples, where $S_{i}$
corresponds to the subset of samples drawn from the unknown distribution of
the random variable $X_{i}$. The problem of estimating the MEV consists in
using the sample means $\hat{\mu}_{i}(S)$ to approximate the largest true mean
$\max_{i}\mu_{i}=\mu_{*}$. Unfortunately, for many distributions (e.g.
Gaussian), no unbiased estimator of the MEV exists (Blumenthal & Cohen, 1968;
Van Hasselt, 2013).

#### Maximum Estimator

The Maximum Estimator (ME) can be defined as

$\hat{\mu}_{*}^{\text{ME}}(S)=\max_{i}\hat{\mu}_{i}\approx\mu_{*}.$ (25)

The ME is the most straightforward way of estimating the MEV and is positively
biased (Smith & Winkler, 2006).

#### Double Estimator

The Double Estimator (DE) is a negatively biased estimator for the MEV (Stone,
1974; Van Hasselt, 2013). DE splits the set of samples $S$ into two
independent sets $S^{A}$ and $S^{B}$, with sample means $\hat{\mu}_{i}(S^{A})$
and $\hat{\mu}_{i}(S^{B})$. Using an cross-validation approach, the DE is
computed by selecting the random variable with the highest sample mean in one
set, and computing its sample mean from samples in the other one. The
resulting values are
$\hat{\mu}_{b^{*}}(S^{A})=\hat{\mu}_{\operatorname*{argmax}_{i}\hat{\mu}_{i}(S^{B})}(S^{A})$
and
$\hat{\mu}_{a^{*}}(S^{B})=\hat{\mu}_{\operatorname*{argmax}_{i}\hat{\mu}_{i}(S^{A})}(S^{B})$.
Finally, the DE can be obtained as

$\hat{\mu}_{*}^{\text{DE}}(S)=\frac{\hat{\mu}_{b^{*}}(S^{A})+\hat{\mu}_{a^{*}}(S^{B})}{2}.$
(26)

#### Weighted Estimator

D’Eramo et al. (2016) introduced the Weighted Estimator (WE), defined as

$\hat{\mu}_{*}^{\text{WE}}(S)=\sum_{i=1}^{M}w_{i}^{S}\hat{\mu}_{i}(S),$ (27)

where $w_{i}^{S}$ is the probability of $\hat{\mu}_{i}$ being the maximum
sample mean. Since the distribution of the sample mean of a random variable
with an unknown distribution is itself unknown, D’Eramo et al. (2016)
approximates the sample means as normally distributed random variables by the
central limit theorem. WE can be both positively and negatively biased, and
its bias is bounded by the one of ME and DE:

$-\frac{1}{2}\left(\sqrt{\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{S_{i}^{A}}}+\sqrt{\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{S_{i}^{B}}}\right)\leq\text{Bias}(\hat{\mu}_{*}^{\text{DE}}(S))\leq\text{Bias}(\hat{\mu}_{*}^{\text{WE}}(S))\leq\text{Bias}(\hat{\mu}_{*}^{\text{ME}}(S))\leq\sqrt{\frac{M-1}{M}\sum_{i=1}^{M}\frac{\sigma_{i}^{2}}{S_{i}}},$
(28)

where $\sigma_{i}$ is the variance of the $i$-th random variable.

Furthermore, WE has often a variance that is empirically smaller to the the
variance of the other two estimators (D’Eramo et al., 2016). A variant of WE
that extends to infinite random variables has been proposed in D’Eramo et al.
(2017).

##  Appendix B Experiments details

In this section we provide details on the experimental setup used for the
empirical evaluation of the proposed method.

For the implementation of the algorithms and the simulation environments we
rely on the following open-source libraries:

  * •

MushroomRL (D’Eramo et al., 2020);

  * •

Gym (Brockman et al., 2016);

  * •

ALE (Bellemare et al., 2013);

  * •

MinAtar (Young & Tian, 2019);

  * •

PyTorch (Paszke et al., 2019).

###  B.1 Asterix

For the experiment on the Asterix Atari game we use the settings of Mnih et
al. (2015), but we add sticky actions (Machado et al., 2018), i.e., a
probability of $p_{\text{repeat}}=25\%$ of repeating the action executed at
the previous frame instead of the one selected by the agent. For the agents we
use the same neural network of Mnih et al. (2015) with, as already mentioned,
the Adam (Kingma & Ba, 2015) optimizer. The additional hyperparameters
introduced by WDQN are tuned with a - small - random search. Concrete Dropout
is used only on the neurons of the last hidden layer. Agents are evaluated
during training every $1$M frames, with the 30 no-op starting condition.
During evaluation the episode length is capped at 30 minutes. Table 1 reports
a list of additional relevant hyperparameters.

Table 1: DQN hyperparameters on Asterix. Hyperparameters marked with a $*$ are used only for WDQN. Hyperparameter | Value
---|---
Optimizer | Adam
Learning rate | $5\mathrm{e}{-5}$
Adam epsilon ($\epsilon_{\text{Adam}}$) | $0.01/32$
Batch size | $32$
Loss function | Huber
Training frequency | $4$
Target network update frequency | $10000$
Min. memory size | $50000$
Max. memory size | $1\text{M}$
Discount factor ($\gamma$) | $0.99$
Initial exploration rate ($\varepsilon_{start}$) | $1.0$
Final exploration rate ($\varepsilon_{end}$) | $0.1$
Exploration steps | $1\text{M}$
Evaluation exploration rate ($\varepsilon_{test}$) | $0.001$
$\text{MC dropout samples}^{*}$ | $100$
Weight decay coefficient $(\lambda)^{*}$ | $1\mathrm{e}{-6}$
Dropout regularization coefficient $(\zeta)^{*}$ | $5\mathrm{e}{-4}$
Initial dropout rate $(p)^{*}$ | $0.5$

###  B.2 Lunar Lander

Table 2: DQN hyperparameters on Lunar Lander. Hyperparameters marked with a $*$ are used only for WDQN. Hyperparameter | Value
---|---
Units per layer | $[100,100]$
Activation | $relu$
Optimizer | Adam
Learning rate | $3\mathrm{e}{-4}$
Batch size | $32$
Loss function | MSE
Training frequency | $1$
Target network update frequency | $300$
Min. memory size | $250$
Max. memory size | $10000$
Discount factor ($\gamma$) | $0.99$
Initial exploration rate ($\varepsilon_{start}$) | $1.0$
Final exploration rate ($\varepsilon_{end}$) | $0.01$
Exploration steps | $1000$
Evaluation exploration rate ($\varepsilon_{test}$) | $0.0$
$\text{MC dropout samples}^{*}$ | $50$
Weight decay coefficient $(\lambda)^{*}$ | $1\mathrm{e}{-6}$
Dropout regularization coefficient $(\zeta)^{*}$ | $2.5\mathrm{e}{-3}$
Initial dropout rate $(p)^{*}$ | $0.2$

For Lunar Lander (Brockman et al., 2016) we limit the length of an episode
(during both evaluation and training) at 1000 steps. To increase the
complexity of the problem, we make the environment stochastic using sticky
actions (see previous subsection) with a $p_{\text{repeat}}=10\%$. We use a
small neural network with only two hidden layers. For WDQN, we use Concrete
Dropout on each hidden layer. Table 2 reports the hyperparameters used to
train the agents.

###  B.3 MinAtar

MinAtar (Young & Tian, 2019) offers a collection of environments resembling
games from the Atari Learning Environment (Bellemare et al., 2013). The state
representation of MinAtar environments is a matrix with multiple channels,
where each channel gives specific information about some aspects of the
environment (e.g., position and speed of a moving object). MinAtar implements
the ALE modifications suggested by Machado et al. (2018), i.e., sticky actions
and difficulty ramping.

We use the same hyperparameters and convolutional neural network of Young &
Tian (2019), but we use Adam for training. The WDQN dropout regularization
coefficient is tuned with random search: we select the value providing
qualitatively better learning curves and keep it fixed across the three games.
Table 3 shows the relevant hyperparameters.

Table 3: DQN hyperparameters on Minatar environments. Hyperparameters marked with a $*$ are used only for WDQN. Hyperparameter | Value
---|---
Optimizer | Adam
Learning rate | $1\mathrm{e}{-4}$
Batch size | $32$
Loss function | Huber
Training frequency | $1$
Target network update frequency | $1000$
Min. memory size | $5000$
Max. memory size | $100000$
Discount factor ($\gamma$) | $0.99$
Initial exploration rate ($\varepsilon_{start}$) | $1.0$
Final exploration rate ($\varepsilon_{end}$) | $0.1$
Exploration steps | $100000$
Evaluation exploration rate ($\varepsilon_{test}$) | $0.0$
$\text{MC dropout samples}^{*}$ | $100$
Weight decay coefficient $(\lambda)^{*}$ | $1\mathrm{e}{-6}$
Dropout regularization coefficient $(\zeta)^{*}$ | $1\mathrm{e}{-4}$
Initial dropout rate $(p)^{*}$ | $0.1$

For WDQN, we run an additional experiment to asses the impact of the dropout
regularization coefficient. We set the initial dropout rate at $p=0.5$ (which
corresponds to the maximum entropy) and test the agents using regularization
coefficients of different magnitude. The results, reported in Figure 5, show
how higher levels of regularization generally correspond to higher entropy of
the weights used to compute the WE.

Figure 5: Learning curves for $3$ different regularization values on three
Minatar games. In the leftmost column, we show the evaluation scores after
each of the 200 training epochs, where $1$ epoch corresponds to $25$k steps.
In the middle column we report for each game and agent the estimate of the
expected return w.r.t. the initial state of the environment; the dashed lines
indicate the real discounted return obtained by the agents. In the rightmost
column we show the moving average of the entropy of the WDQN weights. The
results shown here are averaged over 20 independent runs and the shaded areas
represent $95\%$ confidence intervals. The curves are smoothed using a moving
average of 5 epochs.