# Fundamental Parameters of $\sim$30,000 M dwarfs in LAMOST DR1
Using Data-driven Spectral Modeling

Brianna Galgano Fisk University, Department of Life and Physical Sciences,
W.E.B. DuBois Hall, 1717 Jackson St, Nashville, TN 37208, USA  Keivan Stassun
Vanderbilt University, Department of Physics & Astronomy, 6301 Stevenson
Center Lane, Nashville, TN 37235, USA  Bárbara Rojas-Ayala Instituto de Alta
Investigación, Universidad de Tarapacá, Casilla 7D, Arica, Chile

(Received September 12, 2019; Revised December 9, 2019; Accepted December 11,
2019)

###### Abstract

M dwarfs are the most common type of star in the Galaxy, and because of their
small size are favored targets for searches of Earth-sized transiting
exoplanets. Current and upcoming all-sky spectroscopic surveys, such as the
Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST), offer an
opportunity to systematically determine physical properties of many more M
dwarfs than has been previously possible. Here we present new effective
temperatures, radii, masses, and luminosities for 29,678 M dwarfs with
spectral types M0—M6 in the first data release (DR1) of LAMOST. We derived
these parameters from the supervised machine learning code, The Cannon,
trained with 1,388 M dwarfs in the Transiting Exoplanet Survey Satellite
(TESS) Cool Dwarf Catalog that were also present in LAMOST with high signal-
to-noise ratio ($>$250) spectra. Our validation tests show that the output
parameter uncertainties are strongly correlated with the signal-to-noise of
the LAMOST spectra, and we achieve typical uncertainties of 110 K in $T_{\rm
eff}$ ($\sim$3%), 0.065 $R_{\odot}$ ($\sim$14%) in radius, 0.054 $M_{\odot}$
($\sim$12%) in mass, and 0.012 $L_{\odot}$ ($\sim$20%) in luminosity. The
model presented here can be rapidly applied to future LAMOST data releases,
significantly extending the samples of well characterized M dwarfs across the
sky using new and exclusively data-based modeling methods.

stars: fundamental parameters — stars: low-mass — methods: statistical —
techniques: spectroscopic

††journal: AJ††software: The Cannon (Ness et al. 2015), corner (Foreman-Mackey
2016), matplotlib (Hunter 2007), numpy (Oliphant 2006)

tablenum SIXtablenum

##  1 Introduction

M dwarfs are preferred targets for exoplanet hunting due to their status as
the most common stellar type in the galaxy. (e.g., Mann et al. 2015), and for
their potential in detecting Earth and super-Earth sized exoplanets using both
ground and space-based telescopes via radial velocity measurements and transit
detection (e.g. Shields et al. 2016). The fundamental properties of any host
star such as radius and mass need to be known to a reasonable degree of
certainty in order to characterize the exoplanets they host.

For nearby, bright, well-known M dwarf samples identified via, e.g., very high
proper motions (e.g., Lépine et al. 2013), precise stellar properties can be
estimated from broadband colors and parallax alone. For example, Mann et al.
(e.g., 2019, 2015) have developed empirical calibrations that permit the
determination of radii and masses of M-dwarfs to better than 10% precision
from near-infrared colors and a near-infrared absolute magnitude. However, for
the multitude of M dwarfs that have not yet been identified or characterized
by such catalogs, there remains a need to develop robust spectroscopic methods
which can circumvent problems associated with reddening that can negatively
affect color-based methods, and additionally leverage current and upcoming
large-scale spectroscopic surveys.

However, automated M dwarf characterization via stellar spectroscopy with
synthetic models is difficult given the complicated nature of their stellar
atmospheres. For example, generating model synthetic spectra for stars cooler
than FGK-types, such as M dwarfs, is difficult because of the formation of
molecules in their photospheres e.g., TiO, VO, and CaH in the optical, and H2O
and CO in the near-IR ( e.g., Rojas-Ayala et al. 2012; Shields et al. 2016).
Very heavy absorption can be observed in an M dwarf spectrum by the presence
of molecular compounds which are allowed under their lower effective
temperatures ($T_{\rm eff}$ $<$ 3850 K), which lead to broad and overlapping
spectral lines that can be hard to parse ( e.g., Rojas-Ayala et al. 2012).
Modeling the absorption lines of these compounds is incomplete (e.g., Allard
et al. 2011; Husser et al. 2013), and the lack of a well-defined continuum is
also expected, with some wavelength regions being completely saturated with
absorption features, making normalization for standard spectroscopic analysis
also difficult. Additional complications due to the stars’ intrinsic
properties are deep convection cells and magnetic activity/rapid rotation,
effects that in general are not included in models of low-mass stellar
atmospheres.

The low luminosities of these stars also brings restrictions to both of these
methods; parallax of dwarfs can only be obtained for the very brightest, and
quality (high signal-to-noise ratio) spectra are difficult to obtain. A common
method of practice is to also parameterize the properties of dwarfs if they
are in binary with FGK-type companions, but these are limited to a select
number of dwarfs.

The result is that physics-based models are incomplete and computationally
expensive for M dwarfs which make the method of fitting to synthetic spectra
difficult. This presents a problem—there is a high demand for M dwarf basic
properties for exoplanet and stellar population characterization, but there is
much room for improvement in automated analysis in terms of both computational
speed and accuracy.

A promising alternative is using data-based modeling approaches that use
machine learning to predict what a spectrum would look like given its stellar
properties, and then infer those properties of uncharacterized spectra based
on the model. We discuss in this work how The Cannon
111https://github.com/annayqho/TheCannon (Ness et al. 2015) can successfully
model low resolution, low to moderate signal-to-noise M dwarf spectra, and we
use The Cannon model to find the properties of 29,678 previously
uncharacterized optical spectra of M dwarfs in the LAMOST DR1 catalog.

The Cannon has so far been successfully applied to characterize red giants
using APOGEE spectra and LAMOST spectra and using parameters derived via the
ASPCAP pipeline (Ho et al. 2017a, b). The Cannon has also been applied to much
higher resolution than LAMOST M dwarf spectra (Behmard et al. 2019). LAMOST is
an optical spectroscopic survey all-sky survey for the Northern celestial
hemisphere (Luo et al. 2015). There are $\sim$121,000 M dwarf spectra in the
first data release (DR1), the majority of which have not yet been
characterized beyond a simple spectral subtype determined from color
photometry. We cross-match the LAMOST targets with the TESS Cool Dwarf Catalog
(Muirhead et al. 2018), from which we create a training set to train a
spectral model using The Cannon.

In Section 2, we summarize the data selection and preparation in this work. In
Section 3, we discuss the methods by which we chose an optimal training set
for accuracy and the testing set for applicability. We also discuss the
methods by which we assessed the model’s accuracy in determining basic
parameters: $T_{\rm eff}$, radius, mass, and luminosity. In Section 4, we
report these properties and their uncertainties for the 29,678 M dwarfs in
LAMOST DR1 for which we were able to determine reliable properties. Finally,
Section 5 concludes with a summary of our conclusions.

##  2 Data

In this section, we describe the data that we use, their preparation for use
with The Cannon, and quality control steps taken to train, validate, and apply
The Cannon model to reliable data.

The data are drawn from two main catalogs: (1) the LAMOST first data release
(DR1) M dwarf catalog, is the spectral data set that we wish to classify (Luo
et al. 2015), and (2) the TESS Cool Dwarf Catalog (TCD) (Muirhead et al.
2018), whose parameters we utilize for training and validation of The Cannon
model. The parameters from the TCD (Muirhead et al. 2018) on which we have
chosen to train The Cannon model are effective temperature ($T_{\rm eff}$),
radius ($R$), mass ($M$), and luminosity ($L$); these are the parameters that
we will therefore be able to estimate for the LAMOST DR1 stars.

Our data are split into three subsets, following the approach described by
Ness et al. (2015). The subset of 30,152 TCD stars for which we also have
LAMOST spectra (hereafter, the cross-matched TCD) are divided into a training
set and a validation set. The training set contains the highest signal-to-
noise ratio (SNR) spectra of all TCD and LAMOST cross-matches, paired with
their known associated stellar characteristics or parameters (i.e., labels in
Ness et al. 2015); the training set is used to generate the data-driven model.
The validation set is the remaining cross-matched TCD stars that were not used
in the training set but whose parameters are known and therefore can be used
to assess the fidelity of the model. The LAMOST stars that are included in
neither the training set nor the validation set comprise the survey set. These
are spectra we wish to characterize using the model we created from the
training set, and vetted with the validation set.

###  2.1 LAMOST DR1 Spectra

The Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST) is a
low resolving power (R $\sim$ 1800), optical/near-infrared (3690–9100Å),
ground-based survey (Luo et al. 2015). One of the data products from LAMOST is
the LAMOST DR1 M dwarf catalog (Guo et al. 2015), a catalog of stars whose
spectra have been classified broadly as being of M spectral type but without
more detailed physical information. This catalog consists of 122,677 spectra
with simple spectral type classifications from M0 to M9; these spectra
represent 103,467 unique objects after excluding any duplicate observations of
the same target.

While the LAMOST catalog includes some stars with approximate classifications
as late as M9, we chose to include only stars with a LAMOST DR1 spectral
subtype of M6 and earlier, given that the training set includes only stars
with LAMOST subtypes from M0 to M6 (see Section 3.1). Figure 1 shows the
comparison of the LAMOST spectral subtypes to those from the TCD for the
30,152 spectra that were cross-matched, where the relatively small number of
objects with the latest spectral types is due to difficulty of observing very
faint objects with high signal.

Figure 1: The spectral subtype agreement between cross-matched M stars of
LAMOST DR1 and the TESS Cool Dwarf Catalog (TCD). Note the lack of cross-
matched stars with spectral types later than M6; this effectively sets the
cool limit of our analysis.

The LAMOST DR1 catalog provides the spectra for each of these stars. We
require that to be usable a LAMOST spectrum must include a wavelength solution
($\lambda$), flux ($f$), and flux error (or flux inverse variance;
$1/\sigma_{f}^{2}$), the last of which is used as a statistical weight for
training a model pixel-by-pixel.

####  2.1.1 Basic processing of LAMOST spectra

First, we applied a bad pixel mask using the code supplied with The Cannon
package for this purpose. The bad pixel mask flagged pixels in the spectra
from information provided in the LAMOST meta-data, including poor sky-
subtraction, bad CCD pixels, and infinite or negative flux inverse variance
values. The mask also manually flagged pixels with strong telluric lines and
discontinuities from the joining of LAMOST’s blue and red spectrographs at
5800–6000 Å. For pixels flagged in the bad pixel mask we manually set the
inverse variances to a very small value222We chose a value of $10^{-5}$, as
opposed to 0, which helps avoid errors arising from division by zero in future
operations..

In order for a spectrum to be input into The Cannon for analysis, the spectra
in the training, validation, and survey subsets also must be interpolated to a
common wavelength scale with uniform wavelength spacing. From inspection of
the spectra we found that the LAMOST data were most consistently free of bad
pixels or other problems at wavelengths shorter than $\sim$4500Å, and thus we
interpolated to a common wavelength range of 4500–7500Å, and used a uniform
spacing of 1 Å, as these parameters approximately represent the native format
of the LAMOST spectra for M stars. Each resulting spectrum is a $3000\times 3$
data array, with each of the 3000 pixels containing a standardized wavelength,
flux, and a flux inverse variance float value.

In order to compare fluxes from spectrum to spectrum to create The Cannon
model, each spectrum’s continuum must also be normalized. Traditional methods
of finding and then dividing an M dwarf’s spectrum by its continuum
(continuum-normalization) are challenging because of the complex nature of the
stellar spectra as mentioned in Section 1. Therefore, we adopted the so-called
“pseudo continuum-normalization” method used by Ho et al. (2017b); Behmard et
al. (2019), which first approximates the continuum at wavelength $\lambda_{0}$
via Gaussian smoothing over $L$ pixels and then normalizes the original
spectrum with this smoothed fit:

$f(\lambda_{0})=\frac{\sum_{n}f_{n}{\sigma_{n}}^{-2}w_{n}(\lambda_{0})}{\sum_{n}{\sigma_{n}}^{-2}w_{n}(\lambda_{0})}$
(1)

where index $\sum_{n}$ represents the sum over pixels, and the individual
pixel weights, $w_{n}$, are given by:

$w_{n}(\lambda_{0})=e^{-\frac{(\lambda_{0}-\lambda_{n})^{2}}{L^{2}}}$ (2)

where $\lambda_{n}$, $f_{n}$, $\sigma_{n}$, is the wavelength, flux, and flux
error at the $n$th pixel, respectively. We pseudo-normalized all spectra with
a Gaussian width of $L$ = 50 Å pixels, which is wider than most broader
absorption features, and used the flux inverse variance as weights as
recommended by Ho et al. (2017b). The procedure is visualized with an example
in Figure 2.

We note that this is not a standard robust normalization, in which the maximum
flux is well defined and set to unity, but it is a technique where the spectra
can be quantitatively and consistently compared because they are placed on a
common flux scale. This pseudo-normalization method sometimes over-fits too
deeply into absorption bands (e.g., Figure 2), but because the LAMOST spectra
have the same Gaussian function outlined in Ho et al. (2017b) systematically
applied, it is the preferable technique as it avoids human or physics-model
bias from other standard normalization methods.

Figure 2: An example of a pseudo continuum-normalized using Gaussian smoothing
LAMOST spectrum of an M dwarf from our training set. The effective Gaussian
smoothing width is 50 Å. All spectra from each data subset (training,
validation, and survey) are normalized with this continuum fitting method as
defined by Ness et al. (2015) (see Section 2.1).

Next, LAMOST DR1 objects with no redshift $Z$ parameter reported in the meta-
data cannot be shifted to a rest wavelength frame, a requirement for any
spectrum to be input into The Cannon. Approximately 10,000 stars were
eliminated for this reason.

Finally, through visual inspection of a random set of the LAMOST spectra for
quality checks, some were found to have unusually high flux at blue
wavelengths ($\sim$4500Å and bluer) which is atypical of M dwarf stars and
almost certainly not intrinsic. We chose to eliminate $\sim$1000 stars based
on having especially severe upward trends to the blue, which we defined on the
basis of the slope of the pseudo-continuum being more negative than $-3$ (in
units of normalized flux per unit wavelength). The cause of this high blueward
flux effect is unclear, but might be due to instrumental, calibration, or sky
subtraction error.

####  2.1.2 Spectrum filtering quality control

We require an objective measure of the quality of the LAMOST spectra, both to
select the highest quality subset for training, and to ensure that we do not
attempt to classify very poor spectra.

We chose as our quality metric a measure of signal-to-noise ratio (SNR), which
we calculated from the sum of the measured SNRs in the individual LAMOST
quasi-photometric bands ($ugriz$), hereafter SNRλ. This approach was preferred
over the more formal method of using flux error as it is computationally
faster and independent of the bad pixel masks, and it is readily obtainable
from the LAMOST meta-data. We determined a requirement of SNRλ $>$ 50 based on
validation testing of the model described in Section 4.3; consequently
$\sim$45,000 objects with SNRλ $<$ 50 were eliminated from any further
analysis.

For the training subset, we chose to use a more stringent quality threshold of
SNRλ $>$ 250 (see Figure 3). These choices of SNRλ thresholds are described
more fully in Section 3.2.

Figure 3: Signal to noise ratio measured in $ugriz$ bands (SNRλ) for both
survey and training sets. Cutoff minimums chosen were SNRλ$>$250 for the
training subset and $>$50 for the validation and survey subset (see Section
3.3).

###  2.2 TESS Cool Dwarf Parameters for Training and Validation

The training and validation of The Cannon model requires a set of stars whose
physical parameters (i.e., labels in Ness et al. 2015) are known. We adopt the
parameters provided in the TCD as the known labels with which to train and
validate our model. Those parameters were in turn determined by Muirhead et
al. (2018) on the basis of the empirical relations of Mann et al. (2015).

Note that while the TCD was created in order to enhance planet detection
around cool dwarfs by the TESS mission, the TCD was created more generally to
develop a systematically vetted, brightness and proper-motion limited sample
of cool dwarfs across the sky in a manner that is as unbiased as possible. We
cross-matched the LAMOST DR1 M-dwarf catalog with the entirety of the $\sim$1
million sources in the TCD catalog, not just the TCD targets that were
ultimately selected by the TESS mission for 2-minute cadence observations for
planet detection. Of course, our training set will necessarily inherit any
features inherent in the TCD sample, the most important of which are as
follows (see Muirhead et al. 2018): Unresolved binaries will appear to have
masses and radii that are too large for their color; metallicities are not
known for the majority of the TCD sample, hence we cannot train on metallicity
as a data label and we cannot incorporate any metallicity effects into our
trained model; the catalogs from which the TCD drew its proper motions were
highly incomplete in the Southern hemisphere; and some stars in the TCD relied
on inhomogeneous sources of photometry for colors.

Upon a quality check examination of $T_{\rm eff}$ vs. $M$ and $T_{\rm eff}$
vs. $R$ for these stars, we found a small number of stars whose TCD parameters
deviate from the Mann et al. (2015) empirical relations. These outliers are
most likely the result of inhomogeneous sources of photometry for some stars
and/or unresolved binaries, as noted above.

In order to ensure consistency in the parameters used for the training, we
eliminated these 18 outlier stars that deviated in $R$ by more than 0.01
$R_{\odot}$ and in $M$ by more than 0.01 $M_{\odot}$.

##  3 Methods

In this section we discuss the methods for ensuring the accuracy of the model
produced by The Cannon for M dwarfs and detail how we further optimized the
data-driven model by choosing appropriate SNRλ minimum cutoffs for training,
validation, and survey sets.

Parameter uncertainties are currently not incorporated in The Cannon’s
generative modeling; broadly, the model output uncertainties are formal
uncertainties based on the flux error and closeness of the model fit to a
spectrum. We therefore determine more robust uncertainties on the output
parameters using the scatter of the validation set about known parameters for
the 30,152 TCD/LAMOST cross-matches, to find a function that estimates
parameter error with SNRλ, as described in Section 4.3.

###  3.1 Model Training

In the first step, we trained The Cannon with the training set of 1,388 stars
for which we have high-quality LAMOST spectra as well as parameters ($T_{\rm
eff}$, $R$, $M$, and $L$) from the TCD (see Section 2.2), and that satisfied
the quality control cuts discussed in Section 2.

As outlined in previous work (e.g., Ness et al. 2015; Ho et al. 2017b; Behmard
et al. 2019), we can represent a given spectrum’s flux at wavelength
$\lambda$, $f_{\lambda}$, as the product of the model coefficients
$\theta_{\lambda}$ and the spectrum’s labels (i.e., physical parameters), $l$,
plus noise:

$f_{\lambda}=\theta_{\lambda}\cdot l+noise$ (3)

The noise is the quadrature sum of the measurement error in the observed flux,
$\sigma_{\lambda}$, and the uncertainty in the model coefficients,
$s_{\lambda}$:

$noise=[s_{\lambda}^{2}+\sigma_{\lambda}^{2}]\zeta_{\lambda}$ (4)

where $\zeta_{\lambda}$ is a random Gaussian deviate with zero mean and zero
unit variance (Ness et al. 2015). For our implementation with training labels
taken from the TCD, we define $l$ in terms of the physical parameters $T_{\rm
eff}$, $R$, $M$, and $L$, to second order and including cross-terms:

$\displaystyle{l}\equiv$ $\displaystyle[1,T_{\rm eff},R,M,L,T_{\rm eff}\cdot
R,T_{\rm eff}\cdot M,T_{\rm eff}\cdot L,R\cdot M,$ (5) $\displaystyle R\cdot
L,M\cdot L,T_{\rm eff}^{2},R^{2},M^{2},L^{2}]$

where the first term allows for a linear offset in the fitted flux values
(Behmard et al. 2019).

For a given spectrum, $n$, we can then get the single-pixel log likelihood
function ($\ln p$), which gives the most likely flux value at a specific
wavelength pixel, given model scatter, observational uncertainty, and set of
labels:

$\displaystyle\ln
p(f_{n\lambda}|\theta_{\lambda}^{T},l_{n},s_{\lambda}^{2})=\frac{1}{2}\frac{[f_{n\lambda}-\theta_{\lambda}^{T}\cdot
l_{n}]^{2}}{s^{2}_{\lambda}+\sigma^{2}_{n\lambda}}-\frac{1}{2}\ln(s^{2}_{\lambda}+\sigma^{2}_{n\lambda})$
(6)

where the superscript $T$ denotes the transpose of the matrix. The Cannon then
uses these pixel likelihood probabilities to derive model coefficients that
apply to the full set of Ness et al. (e.g., 2015):

$\theta_{\lambda},s_{\lambda}\leftarrow\operatorname*{argmax}_{\theta_{\lambda},s_{\lambda}}\sum_{n=1}^{N}\ln
p(f_{n\lambda}|\theta_{\lambda}^{T},l_{n},s_{\lambda}^{2})$ (7)

The output of this step is a set of model coefficients for each wavelength
pixel for each stellar parameter. The data-driven model can now be applied to
any other LAMOST spectrum, where the most probable value for a stellar
parameter can be calculated given flux and flux error across all pixels.

###  3.2 Model Validation

Next, to validate The Cannon model, we applied the trained data-driven model
from above to the 30,152 stars comprising the validation subset (see Section
2). We chose 18,844 LAMOST spectra after the same quality checks described in
2, except the SNRλ minimum is reduced to $>$50\. Figure 4 shows the results of
the validation, in which we compare the parameters output by The Cannon to the
input values as adopted from the TCD.

Figure 4: Comparison of output to input labels for the validation data set,
cross-matched between LAMOST DR1 and TCD. While there is a general one-to-one
relationship between the model-derived values and TCD, the model performs
worse at either end of the range of M subtypes (very early or late) as
discussed in Section 3.2. A total of 18,844 spectra, or 11,279 individual
LAMOST objects comprises the validation set.

Overall, we observe good one-to-one agreement between the TCD input and The
Cannon output labels for the validation set in all four of the training labels
($T_{\rm eff}$, $R$, $M$, and $L$). The scatter about the one-to-one line of
agreement is approximately Gaussian (see the insets in each panel of Figure
4), with Gaussian widths of $\sim$100 K, $\approx$0.06 $R_{\odot}$,
$\approx$0.05 $M_{\odot}$, and $\approx$0.01 $L_{\odot}$ in $T_{\rm eff}$,
$R$, $M$, and $L$, respectively. For reference, the uncertainties in the
training labels themselves as reported in the TCD are $\approx$100 K,
$\approx$0.06 $R_{\odot}$, $\approx$0.07 $M_{\odot}$, and $\approx$0.008
$L_{\odot}$, respectively in Muirhead et al. (2018).

We do observe some slight asymmetry in the residuals about the one-to-one
relation, or a “tail” of residuals where The Cannon output labels are a bit
larger than the analogous values of TCD input labels. We also observe some
evidence for breakdown of the model for the coolest stars, i.e., $T_{\rm
eff}$$\lesssim$3150 K (likely due to the small number of M5–M6 stars in the
training set), as well as for the hottest stars, i.e., $T_{\rm
eff}$$\gtrsim$4000 K, where the scatter in the residuals become slightly
discontinuous and asymmetric (Figure 4). The disagreement at the hotter end
($\gtrsim$4000 K) is not surprising, as the TCD’s upper $T_{\rm eff}$ limit is
$\sim$4000 K (Muirhead et al. 2018).

###  3.3 SNRλ Quality Cutoffs

Several factors such as wavelength range, quantity of objects in the training
set, type of continuum-normalization, label choice, and others, can affect the
accuracy of the labels derived by our trained model from The Cannon package.
In our tests with The Cannon model applied to the LAMOST spectra, we found
that SNR of spectra, both for training and validation, was the most important
factor for model accuracy.

In order to determine an appropriate requirement for SNRλ for both the final
training subset described above and for application of the trained model to
the final survey set described below, we examined The Cannon’s performance as
a function of SNRλ in Figure 5.

We found that the accuracy of the model worsened significantly at SNRλ $<$ 50,
and so have opted in the final results reported below to exclude such stars.
For example, at SNRλ$\lesssim 50$, the scatter in the inferred radii and
masses becomes $\gtrsim$15% (see Figure 5). In addition, we selected a
threshold of SNRλ $>$ 250 for the training set, as these represent the very
highest quality spectra while still possessing a significant number of objects
($\sim$2300) to include in the training set. Moreover, inspection of the
trends in Figure 5 reveal that the scatter in the inferred parameters does not
improve beyond SNRλ$\approx 250$, thus we assume that including all objects
with SNRλ$>250$ in the training set will maximize the performance of the
trained model. With this cutoff of SNRλ$<250$ and other quality checks, the
final training set consists of 1,388 spectra that are each unique objects
(with unique LAMOST designations) out all 24,546 unique objects available that
are cross-matches with TCD.

Figure 5: Relationship between SNRλ and label residuals between cross-matches
of LAMOST DR1 and TCD. Each point represents the bin average and its
horizontal bar represents that bin’s width. We note bin width increases at
high SNRλ so that each bin contains approximately the same number of objects
(1360$\pm$5; 860 in the bin for SNRλ $>$ 239). A small number of stars beyond
SNRλ $>$ 300 are not shown. The coefficients of the best-fit 2nd degree
polynomial are provided in Table 2. The solid line represents the adopted
error floor for SNRλ $>$ 225\.

##  4 Results

After application of the quality and other cuts described in the preceding
sections, we have been able to successfully apply The Cannon to 29,678 M
dwarfs, for which we report newly characterized $T_{\rm eff}$, $R$, $M$, and
$L$ labels. A sample of these final results are reported in Table 1.

Designation | $T_{\rm eff}$ |  $\sigma$$T_{\rm eff}$ | $R$ |  $\sigma$$R$ | $M$ |  $\sigma$$M$ | $L$ |  $\sigma$$L$ | SNRλ | SpType | $\chi^{2}$
---|---|---|---|---|---|---|---|---|---|---|---
J080546.84+202017.8 | 3469.6 | 128.5 | 0.376 | 0.064 | 0.384 | 0.077 | 0.026 | 0.013 | 58.4 | M4 | 1998.0
J064420.73+323838.6 | 3680.1 | 106.8 | 0.481 | 0.054 | 0.513 | 0.063 | 0.049 | 0.014 | 264.2 | M1 | 1899.0
J131151.74+565947.5 | 3837.4 | 99.0 | 0.560 | 0.049 | 0.607 | 0.059 | 0.070 | 0.012 | 119.9 | M1 | 1312.0
J095520.46+301012.0 | 3732.4 | 129.2 | 0.513 | 0.064 | 0.551 | 0.077 | 0.055 | 0.013 | 57.3 | M2 | 3592.7
J100744.45$-$020217.4 | 3995.9 | 124.2 | 0.626 | 0.062 | 0.677 | 0.074 | 0.095 | 0.013 | 65.2 | M0 | 2512.6
Table 1: Resulting stellar parameter values from applying The Cannon M dwarf
model to LAMOST spectra. The formal error (1$\sigma$) has been estimated from
the best fitted polynomial function of SNRλ (see Section 4.3, Figure 5). Units
for $T_{\rm eff}$, $R$, $M$, and $L$ are K, $R_{\odot}$, $M_{\odot}$, and
$L_{\odot}$, respectively. Spectral types are as reported by LAMOST DR1. The
full table is available in the electronic version of the Journal; a small
sample is provided here for guidance regarding its form and content.

###  4.1 Final Eliminations

We note that a small fraction of the stars analyzed and that otherwise
satisfied our quality cuts nonetheless failed to be characterized by The
Cannon. This was the case for 2,445 stars, for which at least one or more of
the output labels failed, identified by NaN or negative label values.

We visually inspected these failed cases, noting that in many instances the
LAMOST spectrum had a very large positive or negative outlier feature, perhaps
due to bad pixels that were not masked. However, there were also instances
where visually the spectrum did not appear unusual. There were no other
obvious trends in the characteristics of the failed spectra; they covered a
large SNRλ range and included both early and late subtypes. In conclusion, we
were not able to determine an obvious common cause for this failure. In any
case, these have been completely removed from our final characterized set in
Table 1 and what is shown in Figures 6, 7, and 8.

Finally, we also eliminated a small number of stars whose final goodness-of-
fit statistic ($\chi^{2}$) from The Cannon was very poor. This is shown in
Figure 6, where we adopted a cutoff of $\chi^{2}<$ 5,000, corresponding to a
$\chi^{2}$ of $\sim$3 per degree of freedom.

Figure 6: Distribution of ${\chi}^{2}$ values for The Cannon model fit for all
spectra in the survey set. Since the majority of stars had $\chi^{2}$ $<$
5,000 (corresponding to $\chi^{2}$ less than $\sim$3 per degree of freedom),
we eliminate any object that had a $\chi^{2}$ greater than this threshold (see
Section 4.1).

###  4.2 Model Parameters

The parameters of the M dwarfs The Cannon has calculated from LAMOST DR1
spectra are as expected for their spectral types. Their ranges are
2901$<$$T_{\rm eff}$$<$4113 K, 0.14$<R<$0.66 $R_{\odot}$, 0.10$<M<$0.71
$M_{\odot}$, and 0.002$<L<$0.115 $L_{\odot}$. Note that the relatively smaller
number of characterized stars at the coolest $T_{\rm eff}$ (latest spectral
subtypes) is at least in part due to our quality cuts that require SNRλ $>50$
and that our training set had fewer stars in this range.

While the LAMOST pipeline does not provide a $T_{\rm eff}$ estimate for most M
dwarfs in the catalog, the catalog does provide an estimated subtype, which we
used as a sanity check on the $T_{\rm eff}$ output by our trained model (see
Figure 7). We find that there is a reasonable progression from higher to lower
$T_{\rm eff}$ derived from The Cannon corresponding to increasingly later
subtypes reported by LAMOST. The overall distribution of model-derived
parameters of the newly characterized 33,095 LAMOST spectra can be seen in
Figure 8.

Figure 7: Comparison of subtype as reported in the LAMOST DR1 M catalog (Guo
et al. 2015) to the outputed model-derived $T_{\rm eff}$ by The Cannon.
Figure 8: The basic parameter space of our characterized survey set as
discussed in Section 4.2. Highlighted in gray is the $T_{\rm eff}$ for which
the The Cannon model we have trained is applicable; dashed lines represent the
other parameter ranges corresponding to this applicable $T_{\rm eff}$ range.

###  4.3 Error Estimation

As described above, we experimented by running The Cannon with many different
training subsets and assessing the impact of different input variables on the
scatter in the resulting output parameters. We found that SNRλ is by far the
most important predictor of the parameter scatter (see Figure 5) versus other
potential factors (e.g., continuum normalization method, different pixel
masks, etc). Therefore, we proceed to derive stellar parameter errors for a
particular spectrum as a function of its SNRλ. Figure 5 shows the empirical
relationship between the standard deviation of the residuals of each label
from the validation step and the SNRλ.

We have opted to approximate the relationship between parameter scatter and
SNRλ with a simple polynomial function. This is the simplest form that appears
to reasonably well represent the empirical relationship—it is not derived from
first principles—however, we believe it fully suffices for our purposes here.
Table 2 reports the coefficients of the fitted polynomial relationships, which
we use to report the final estimated uncertainties on the labels that we
report in Table 1.

$Y$ | units | $a$ | $b$ | $c$ | error floor
---|---|---|---|---|---
$\sigma$$T_{\rm eff}$ | K | $2.594\text{\times}{10}^{-3}\text{\,}$ | $-0.942$ | $174.639$ | $89.2$
$\sigma$$R$ | $R_{\odot}$ | $1.3\text{\times}{10}^{-6}$ | $-4.692\text{\times}{10}^{-4}$ | $8.688\text{\times}{10}^{-2}$ | $0.0432$
$\sigma$$M$ | $M_{\odot}$ | $1.6\text{\times}{10}^{-6}$ | $-5.841\text{\times}{10}^{-4}$ | $0.106$ | $0.0516$
$\sigma$$L$ | $L_{\odot}$ | $2.0\text{\times}{10}^{-7}$ | $-6.3\text{\times}{10}^{-5}$ | $1.623\text{\times}{10}^{-2}$ | $0.0097$

Table 2: Polynomial fit of parameter error versus SNRλ (see Figure 5). We
adopt a minimum error floor for stars with SNRλ $>$ 225, as indicated in the
final column (red solid lines in Figure 5).

Note. — $Y(X)=a+bX+cX^{2}$, where $X\equiv$ SNRλ.

We do note that the uncertainties we infer in this way for the high SNRλ
training sample are very similar to those reported by Muirhead et al. (2018),
even though The Cannon model does not incorporate label uncertainties in its
training. In other words, the model appears to be naturally reproducing the
expected errors for the highest quality sample, effectively setting the error
floor for the output parameters.

##  5 Summary

We were successfully able to extend the range of use for The Cannon in
determining properties of M dwarfs to low-resolution, low signal-to-noise,
optical spectra. We obtained spectra from the LAMOST DR1 M dwarf catalog and
properties from the TESS Cool Dwarf Catalog, and performed a cross-label
transfer similar to the procedure detailed in Ho et al. (2017b). First, we
performed quality checks on the spectra before implementing the model, such as
filtering out poor spectra and establishing minimum signal-to-noise ratio
(SNRλ) cutoffs (Section 2.1.2, Section 3.3). We selected a training subset
with the best available (SNRλ$>$250) LAMOST spectra that also had reasonable
TESS Cool Dwarf Catalog parameters to train the model (Section 2.2). We then
assessed the model’s validity and calculated parameter errors as a function of
SNRλ (Section 3.2, Section 4.3). We also eliminated objects from our final
characterized survey set according to poor $\chi^{2}$ and kept only objects
for which all four model parameters ($T_{\rm eff}$, $R$, $M$, $L$) could be
determined successfully.

In the end, were able to apply this method to determine the properties of
$\sim$30,000 M dwarfs observed by LAMOST that had not been otherwise
characterized. We achieve typical uncertainties of 110 K in $T_{\rm eff}$
($\sim$3%), 0.065 $R_{\odot}$ ($\sim$14%) in radius, 0.054 $M_{\odot}$
($\sim$12%) in mass, and 0.012 $L_{\odot}$ ($\sim$20%) in luminosity, driven
almost entirely by the precision and range of the training set. The model
presented here can be rapidly applied to future LAMOST data releases,
significantly extending the samples of well characterized M dwarfs across the
sky using new and exclusively data-based modeling methods.

There is room for future improvements. One particular area for improved
accuracy would be to incorporate parameter measurement uncertainties as as an
additional weight for the model. The Cannon approach can also be sensitive to
outliers, especially if the training set is small. We therefore recommend
careful scrutiny and removal of outliers, such as by the procedures we adopted
in Sections 2.1.2 and 2.2.

Even so, this work demonstrates that The Cannon model can be applied to M
dwarfs with a smaller training set than those used in previous works (e.g.,
Ness et al. 2015, 2016; Ho et al. 2017b, a), and can be used specifically to
determine basic stellar properties from LAMOST M dwarf spectra with high
accuracy and speed.

B.G. acknowledges partial funding support from NSF PAARE grant AST-1358862
through the Fisk-Vanderbilt Masters-to-PhD Bridge Program. B.R-A. acknowledges
funding support from FONDECYT through grant 11181295.

## References

  * Allard et al. (2011) Allard, F., Homeier, D., & Freytag, B. 2011, in Astronomical Society of the Pacific Conference Series, Vol. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. C. Johns-Krull, M. K. Browning, & A. A. West, 91
  * Behmard et al. (2019) Behmard, A., Petigura, E. A., & Howard, A. W. 2019, arXiv e-prints, arXiv:1904.00094
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, The Journal of Open Source Software, 24, doi:10.21105/joss.00024
  * Guo et al. (2015) Guo, Y.-X., Yi, Z.-P., Luo, A.-L., et al. 2015, Research in Astronomy and Astrophysics, 15, 1182
  * Ho et al. (2017a) Ho, A. Y. Q., Rix, H.-W., Ness, M. K., et al. 2017a, ApJ, 841, 40
  * Ho et al. (2017b) Ho, A. Y. Q., Ness, M. K., Hogg, D. W., et al. 2017b, ApJ, 836, 5
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6
  * Lépine et al. (2013) Lépine, S., Hilton, E. J., Mann, A. W., et al. 2013, AJ, 145, 102
  * Luo et al. (2015) Luo, A. L., Zhao, Y. H., Zhao, G., et al. 2015, arXiv e-prints, arXiv:1505.01570
  * Mann et al. (2015) Mann, A. W., Feiden, G. A., Gaidos, E., Boyajian, T., & von Braun, K. 2015, ApJ, 804, 64
  * Mann et al. (2019) Mann, A. W., Dupuy, T., Kraus, A. L., et al. 2019, ApJ, 871, 63
  * Muirhead et al. (2018) Muirhead, P. S., Dressing, C. D., Mann, A. W., et al. 2018, AJ, 155, 180
  * Ness et al. (2015) Ness, M., Hogg, D. W., Rix, H. W., Ho, A. Y. Q., & Zasowski, G. 2015, ApJ, 808, 16
  * Ness et al. (2016) Ness, M., Hogg, D. W., Rix, H. W., et al. 2016, ApJ, 823, 114
  * Oliphant (2006) Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (Trelgol Publishing USA)
  * Rojas-Ayala et al. (2012) Rojas-Ayala, B., Covey, K. R., Muirhead, P. S., & Lloyd, J. P. 2012, ApJ, 748, 93
  * Shields et al. (2016) Shields, A. L., Ballard, S., & Johnson, J. A. 2016, Phys. Rep., 663, 1


# Bag of Policies for Distributional Deep Exploration

Asen Nachkov  Department of Computing
Imperial College London
London, UK  Luchen Li  Department of Computing
Imperial College London
London, UK  Giulia Luise  Department of Computing
Imperial College London
London, UK  Filippo Valdettaro  Department of Computing
Imperial College London
London, UK  Aldo Faisal  Department of Computing
Imperial College London
London, UK  Department of Bioengineering
Imperial College London
London, UK  Lehrstuhlinhaber für Digital Health & Data Science
University of Bayreuth
Bayreuth, Bavaria, Germany

###### Abstract

Efficient exploration in complex environments remains a major challenge for
reinforcement learning (RL). Compared to previous Thompson sampling-inspired
mechanisms that enable temporally extended exploration, i.e., deep
exploration, we focus on deep exploration in distributional RL. We develop
here a general purpose approach, Bag of Policies (BoP), that can be built on
top of any return distribution estimator by maintaining a population of its
copies. BoP consists of an ensemble of multiple heads that are updated
independently. During training, each episode is controlled by only one of the
heads and the collected state-action pairs are used to update all heads off-
policy, leading to distinct learning signals for each head which diversify
learning and behaviour. To test whether optimistic ensemble method can improve
on distributional RL as did on scalar RL, by e.g. Bootstrapped DQN, we
implement the BoP approach with a population of distributional actor-critics
using Bayesian Distributional Policy Gradients (BDPG). The population thus
approximates a posterior distribution of return distributions along with a
posterior distribution of policies. Another benefit of building upon BDPG is
that it allows to analyze global posterior uncertainty along with local
curiosity bonus simultaneously for exploration. As BDPG is already an
optimistic method, this pairing helps to investigate if optimism is
accumulatable in distributional RL. Overall BoP results in greater robustness
and speed during learning as demonstrated by our experimental results on ALE
Atari games.

##  1 Introduction

Distributional RL (DiRL) has rapidly established its place among reinforcement
learning (RL) algorithms Bellemare et al. [2017] as a powerful improvement
over non-distributional value-based counterparts Lyle et al. [2019]. In DiRL,
the agent does not learn a single summary statistic of the return for each
state-action pair, but instead learns the whole return distribution. The
agent’s behaviour is being evaluated for multiple possible consequences which
in turn affect the policy update. While this does lead to more stable learning
and better performance Lyle et al. [2019], it does not itself change the way
actions are selected; as distributional extensions to value-based RL, in C51
Bellemare et al. [2017], QR-DQN Dabney et al. [2018b] the agent still takes
actions according to the mean of the estimated return distributions in each
state-action pair. Thus, estimating a return distribution, at its core,
provides performance advantages from better representation and evaluation
which are unrelated to the action-selection and the exploration behaviour of
the agent.

Figure 1: Schematic and data flow in BoP.

One of the major challenges in DiRL is to better leverage an estimated return
distribution for action selection Dabney et al. [2018a]. In a discrete action
setting, the agent could take the action maximizing the mean plus a time-
dependent scaled standard deviation of the return Mavrin et al. [2019]. This
approach, combined with computing the upper variance as an embodiment of the
“optimism in the face of uncertainty”, has yielded considerable improvements
in exploration. A DiRL solution based on Wasserstein-GAN Freirich et al.
[2019] captures uncertainty in the gradient magnitude of parameters. In the
recent Bayesian Distributional Policy Gradients (BDPG) Li and Faisal [2021], a
curiosity bonus in the form of an information gain is added to the reward
function, which motivates the agent to choose actions whose outcome has high
epistemic uncertainty in terms of the return distribution estimation. These
strategies all bias the data-collecting policy with a curiosity component
computed from the return distribution modelling, instead of using only the
mean returns as done in conventional RL.

In the meantime, one of the most effective exploration approaches centred
around the idea of optimism in the face of uncertainty in scalar RL, is those
based on posterior sampling or Thompson sampling Thompson [1935], as we shall
discuss in more detail in the Related Work section. These methods typically
maintain an ensemble of randomized copies of the same function, relying on
model diversity to select actions that are either relatively certainly optimal
or uncertain, and therefore efficiently expand the data space with maximal
information gathering.

In this work, we want to combine ensemble method with already optimistic DiRL
to see if extra optimism from diversified model estimation is beneficial. To
this end, we enhance BDPG Li and Faisal [2021] into multiple distributional
actor-critic models to achieve a combination of both Thompson sampling and
curiosity bonus induced optimism. Crucially, we are now able to leverage the
epistemic uncertainty of the agent in a two-pronged manner: First, through the
Thompson sampling enabled by an ensemble of policies and return distributions
(heads); Second, by using a local information gain inferred from the observed
inaccuracies in each head’s return distribution estimate.

Our contributions in this paper are as follows: we introduce the Bag of
Policies framework, analyze its properties and the off-policy training
aspects. We also explore variations of our algorithm and test its performance
on the Atari Arcade Learning Environment Bellemare et al. [2013]. Our results
suggest that Thompson sampling can indeed improve upon already optimistic
DiRL. This finding is promising, as we now know that optimism is accumulatable
in DiRL.

##  2 Preliminaries

A reinforcement learning (RL) task is modelled as a Markov decision process
(MDP) ($\mathcal{S},\mathcal{A},\mathcal{R},P,\gamma$) Puterman [1994] where
$\mathcal{S}$ is the set of possible states of the environment, $\mathcal{A}$
is the set of actions the agent can take, $\mathcal{R}$ is a reward function ,
$P$ is the transition probability density for transitioning to a new state
given the agent’s current state and their selected action, and $\gamma$ is the
discount factor. States and actions can be continuous in our framework. The
agent also learns or controls a possibly probabilistic policy
$\pi:\mathcal{S}\rightarrow P(\mathcal{A})$ mapping from states to a
distribution over actions.

In our implementation the policy $\pi$ is parameterized and updated with
policy gradient Sutton [1999], Mnih et al. [2016], Schulman et al. [2016],
maximizing an explicit objective function of the performance of the agent such
as $\mathbb{E}_{s\sim
d^{\pi},a\sim\pi(\cdot|s)}\big{[}A(s,a)\log\pi(a|s)\big{]}$, where $d^{\pi}$
is the stationary marginal state density induced by $\pi$ and $A(s,a)$ the
advantage function.

The return for a given state $s\in\mathcal{S}$ is defined as the sum of
discounted future rewards that the agent collects starting from the state $s$,
$G^{\pi}(s):=\sum_{t=0}^{\infty}\gamma^{t}r_{t},\ s_{0}=s$. We consider state-
dependent return is this works, while the idea is applicable also to action-
dependent return. The return is a random variable with potential aleatoric
uncertainty (e.g. from stochastic state transitions). DiRL methods learn not
only a single statistic of the return in a given state, but a representation
of the whole distribution. Similar to the Bellman operator defined in
conventional RL with respect to the mean of the return, one can use the
distributional Bellman operator Bellemare et al. [2017].

$\mathcal{T}^{\pi}G^{\pi}(s):\stackrel{{\scriptstyle D}}{{=}}R(s)+\gamma
G^{\pi}(s^{\prime}),$ (1)

where $R(s)$ is the random variable of the reward and $s^{\prime}\sim
P(\cdot|a,s),a\sim\pi(\cdot|s)$, to learn the state return distribution.

We test our idea of improving optimistic DiRL using (also optimistic) deep
exploration on the DiRL method BDPG Li and Faisal [2021], as it provides an
optimistic exploration bonus. BDPG interprets the distributional Bellman
operator as a variational auto-encoding (VAE) process. Meanwhile, this VAE
provides an estimate of information gain between samples from target and model
return distributions, giving rise to the quantification of the epistemic
uncertainty in return distribution modelling around the current state under
the current policy. This information gain is then used to augment reward
signals to encourage exploration. As an already (locally) optimistic DiRL
method, BDPG is a good baseline to pair with Thompson sampling for our
investigation of buildable optimism in DiRL.

We introduce the superscript $i$ in the notation $\pi^{i}$ to denote the
$i$-th policy in an ensemble of $K$ distributional estimator actor-critic
heads. Accordingly, $A_{t}^{i}$ is the advantage of the $i$-th policy at
timestep $t$. We adopt the notations used in BDPG to shorthand the model
output state-return as $g(s):=G^{\pi}(s)$, and its Bellman backup target as
$x(s):=\mathcal{T}^{\pi}G^{\pi}(s)$. They both belong to the return space.
BDPG matches $g(s)$ to $x(s)$ in Wasserstein metric (under which the
distributional Bellman operator is known to be a contraction Bellemare et al.
[2017]) with a deterministic decoder $p_{\theta}(x|z,s)$ that maps a latent
random variable $z$ to a return sample $g(s)$, a variational encoder
$q_{\phi}(z|x,s)$ that approximates the posterior over $z$ conditional on
observed data $x(s)$, and a discriminator $D_{\psi}(x,z,s)$ to enforce joint-
contrastive learning Donahue et al. [2017] now that the decoder is non-
differentiable. For algorithm details please refer to the original paper Li
and Faisal [2021]. Dependency on $s$ is omitted when no confusion is to be
incurred. We adapt per-head shorthand from BDPG notations for the Bellman
targets as $x_{t}^{i}$, the samples from the generators $g_{t}^{i}$, and
latent variables for the generative process $z_{t}^{i}$.

##  3 Related Work

Closest to our work is the multi-worker or multi-head approaches which
leverage Thompson sampling Thompson [1935] to diversify behaviour. Notice that
unlike parallel approaches such as A3C Mnih et al. [2016] and variants Liang
et al. [2018], Zhang et al. [2019] in which the multiple workers are updated
with the same gradient, Thompson sampling allows the workers to learn from
bootstrapped signals Efron and Tibshirani [1994] and to focus on where the
current understanding is insufficient rather than merely cover more state
space more quickly together than a single worker.

One early example of Thompson sampling in deep RL is the Bootstrapped DQN
Osband et al. [2016] which utilizes multiple $Q$-value “heads”. Each head is
initialized randomly and differently and then trained on their respective
learning signal (with slightly different data or even identical data but
distinct gradient). Thanks to this diversity, the heads’ evaluations for a
state-action pair tend to converge to the true value (which then allows the
optimal action to be chosen) as the pair is being sufficiently visited. Thus
said, an action is selected either because it is optimal or because its
evaluation is uncertain, thus optimistic. This uncertainty is referred to as
the epistemic uncertainty that would diminish as training data increases.
Bootstrapped DQN and variants Osband et al. [2019], Chen et al. [2017],
O’Donoghue et al. [2018] can be thought of as representing a posterior
distribution of the $Q$-function and thus gives a measure of epistemic
uncertainty in the $Q$-value estimation. Similarly, Tang and Agrawal [2018]
approximates Thompson sampling by drawing network parameters for an RL model.
To further encourage optimism in the face of uncertainty, gradient
diversification An et al. [2021] can effectively boost performance of an
ensemble method.

Ensemble method is important as it also allows for deep exploration by
engaging the same policy for the duration of the whole episode. However, all
ensemble-based exploration schemes for now are intended for conventional
scalar RL, and the purpose of this work is to investigate whether deep
exploration could improve performance in distributional RL, further still,
whether the global optimism incurred by the ensemble effect can further
enhance exploration efficiency when the agent is already optimistically biased
on per-state basis. Combining ensemble technique with DiRL is not a trivial
investigation and the benefit of diversity does not automatically propagate,
as now we would be learning a posterior distribution over _distributions_ ,
whilst admittedly the diversity over distributions is harder to maintain than
that over numbers.

We base our approach on BDPG Li and Faisal [2021] which is a DiRL approach
made up of a single distributional actor-critic and provides a local curiosity
bonus favouring states whose return distribution estimations suffer from high
epistemic uncertainty. BoP in contrast maintains a posterior distribution over
return distributions and over policies, keeping track of both aleatoric and
epistemic uncertainties about return distribution estimation simultaneously.

##  4 Bag of Policies (BoP)

In the following we lay out the structure, theory and variants of BoP. A
schematic can be found in Fig. 1 and for the pseudocode please refer to
Algorithm$\leavevmode\nobreak\ 1$.

Algorithm 1 Bag of Policies

1: Initialize prior $p_{\theta}(z|s)$, encoder $q_{\phi}(z|x,s)$,
discriminator $D_{\psi}(x,z,s)$, actor-critic population
$\\{G_{\theta}^{i}(z,s),\pi^{i}\\}_{i=1}^{K}$

2: while not converged do

3: // Roll-out stage

4: training batch $\mathcal{D}\leftarrow\emptyset$

5: Sample $k\sim\mathrm{Uniform}(1,...,K)$ // actor-critic to act

6: for $t=0$ to $T-1$ do

7: execute $a_{t}^{k}\sim\pi^{k}(\cdot|s_{t})$, get $r_{t},s_{t+1}$

8: sample return $z_{t}\sim p_{\theta}(\cdot|s_{t}),g_{t}^{i}\leftarrow
G_{\theta}^{i}(z_{t},s_{t}),\ \forall i$

9:
$\mathcal{D}\leftarrow\mathcal{D}\cup\big{(}s_{t},a_{t},r_{t},g_{t}^{i},\log\pi^{i}(a_{t}|s_{t})\big{)},\forall
i$

10: end for

11: sample last return
$z_{T}\sim p_{\theta}(\cdot|s_{T}),g_{T}^{i}\leftarrow
G_{\theta}^{i}(z_{T},s_{T}),\ \forall i$

12: // Data-driven estimation stage

13: for $t$ $\in$ $\mathcal{D},\leavevmode\nobreak\ \forall i$ do

14: estimate off-policy advantage $A_{t}^{i}$

15: compute Bellman target $x_{t}^{i}\leftarrow A_{t}^{i}+g_{t}^{i}$

16: compute mixed curiosity reward
$\hat{r}_{t}^{i}\propto\text{KL}\big{(}q_{\phi}(\cdot|x_{t}^{k},s_{t})\ ||\
q_{\phi}(\cdot|g_{t}^{i},s_{t})\big{)}$

17: augment $A_{t}^{i}$ by replacing $r_{t}$ with $r_{t}+\hat{r}_{t}^{i}$

18: end for

19: // Update stage

20: // Train with minibatch $B\subset\mathcal{D}$

21: for $t$ $\in$ $B$ do

22: encode $\tilde{z}_{t}^{i}\sim
q_{\phi}(\cdot|x_{t}^{i},s_{t}),\leavevmode\nobreak\ \forall i$

23: sample adversaries
$z_{t}^{i}\sim p_{\theta}(\cdot|s_{t}),\ \tilde{x}_{t}^{i}\leftarrow
G_{\bar{\theta}}^{i}(z_{t}^{i},s_{t}),\ \forall i$

24: take averages
$\bar{z}_{t}=\frac{1}{K}\sum_{i=1}^{K}\tilde{z}_{t}^{i},\bar{x}_{t}=\frac{1}{K}\sum_{i=1}^{K}\tilde{x}_{t}^{i}$

25: end for

26: update $D_{\psi}$ by ascending
$\frac{1}{|B|}\sum_{t\in B}\big{[}\log
D_{\psi}(\bar{x}_{t},z_{t},s_{t})+\newline \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\log\big{(}1-D_{\psi}(x_{t},\bar{z}_{t},s_{t})\big{)}\big{]}$

27: update encoder, prior by ascending
$\frac{1}{|B|}\sum_{t\in
B}\big{[}\log\big{(}1-D_{\psi}(\bar{x}_{t},z_{t},s_{t})\big{)}+\newline
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \log
D_{\psi}(x_{t},\bar{z}_{t},s_{t})\big{]}$

28: update $G_{\theta}^{i},\ \forall i$ by descending
$\frac{1}{|B|}\sum_{t\in
B}||x_{t}^{i}-G_{\theta}^{i}(\tilde{z}_{t}^{i},s_{t})||_{2}^{2}$

29: update $\pi^{i},\ \forall i$ by ascending
$\frac{1}{|B|}\sum_{t\in B}\log\pi^{i}(a_{t}|s_{t})A_{t}^{i}$

30: end while

###  4.1 Architecture

The Bag of Policies (BoP) framework can be applied to any DiRL estimator. In
this work we have chosen BDPG Li and Faisal [2021] as the implementation
framework because it exploits local epistemic uncertainty. We use the
architecture as that of BDPG except that we have multiple distributional
actor-critic heads. Specifically, we refer to each distributional actor-critic
pair (policy + return distribution) as a head. Each episode of trajectory is
rolled out by only one of the heads sampled uniformly randomly, and the
collected data is trained on by all of the heads for data efficiency. The head
that generated the trajectory learns on-policy, whilst all other heads learn
from the same batch of data off-policy. Throughout training, the heads are
updated by different gradients because they were initialized differently and
learn from their own Bellman target which entails random sampling from their
own return distributions at subsequent steps. This ensures a diversity in
policies, giving rise to the effect of Thompson sampling.

Figure 2: Comparison of BDPG and BoP on selected Atari environments.

The training of BoP entails an outer loop of three stages: a roll-out stage, a
data-driven estimation stage, and an update stage. The majority of the
algorithm is the same as that of BDPG, we only highlight the different or
extra aspects incurred by using multiple distributional actor-critic heads, as
well as the heed we paid to maintain diversity. During the roll-out stage, at
the start of each episode, the agent selects a head uniformly at random, which
we call the active head $k$, and executes actions according to the policy of
that head for the duration of the whole episode. Thus, the sequence of
collected transitions $(s_{t},a_{t},r_{t})$ is determined entirely by the
policy of the active head as a behaviour policy. Meanwhile, the other heads
still compute the log-probability of the selected action
$\log\pi^{i}(a_{t}|s_{t})$, and form an estimation of the return distribution
for each state in the episode by generating a sample $g_{t}^{i}$ from each
distributional critic. These local evaluations are used to provide them with a
unique sample estimate of the learning objective when updating their
parameters.

During the sample estimate stage, the agent computes the advantages
$A_{t}^{i}$ and the Bellman targets $x_{t}^{i}$, both being unique for each
respective head. These local estimates also engender a local curiosity bonus
$\hat{r}_{t}^{i}$ that is proportional to the entropy reduction in the latent
space conditioned on the return space due to the condition being a more
accurate estimate (the on-policy Bellman target $x_{t}^{k}$) as opposed to a
model prediction (a sample from the model $g_{t}^{i}$):

$\displaystyle\hat{r}_{t}^{i}$
$\displaystyle\propto\mathrm{Information\\_gain}(s_{t},i)$
$\displaystyle:=\mathrm{KL}\big{(}q_{\phi}(\cdot|x_{t}^{k},s_{t})\ ||\
q_{\phi}(\cdot|g_{t}^{i},s_{t})\big{)},$ (2)

in which the state $s_{t}$ is generated by the active head $k$. This entropy
reduction, or equivalently, information gain, quantifies the epistemic
uncertainty about how well each head is modelling the return distribution. The
posterior $x_{t}^{k}$ is shared across heads as it is the only Bellman backup
target among all $x_{t}^{i}$ that is estimated on-policy and thus suffers the
least degree of bias / variance. The curiosity bonus is then added to the
external reward for computing the augmented advantage function that would be
used to update the local policies.

During the update stage, all policies $\pi^{i}$ are updated on the data
generated by the active policy $\pi^{k}$. The active head is updated on-policy
as in BDPG. For the other heads we use an off-policy advantage function
estimated with V-trace Espeholt et al. [2018] at each timestep. Subsequently,
each head’s critic is updated by minimizing a Wasserstein distance between its
current return distribution model and the distribution of a Bellman target
estimated from its own predictions at subsequent timesteps. Like BDPG, this is
achieved by minimizing the squared distance between Bellman target $x_{t}^{i}$
and sample $g_{t}^{i}$ while in the meantime enforcing joint-contrastive
learning adversarially.

We used a single common set of hyperparameters for implementing the BoP on
Atari as shown in Table 1. The choices for the number of heads were laid out
in Table 2. The policy heads are updated by Proximal Policy Optimization (PPO)
Schulman et al. [2017].

Hyperparameter | Value | Hyperparameter | Value
---|---|---|---
Num. parallel envs | 16 | Evaluation frequency |  Every $2e4$ timesteps
Num. stacked frames | 4 | Num. heads | Variable
GAE $\lambda$ | 0.95 |  Discount $\gamma$ | 0.99
Num. episodes per roll-out | 4 | Minibatch size | 256
Roll-out timesteps | 128 | Value loss weight | 0.5
Entropy coefficient | 0.01 | Learning rate |  Linearly $2.5e-4\to 0$
PPO clip range |  Linearly $0.1\to 0$ | Curiosity | Mixed
Actor-critics updated | All & on all the data |  $x_{t}^{i}$ computed by |  Critic $i$
Actions during roll-out | Sampled | Actions during testing | Greedy
Table 1: Hyperparameters for implementing BoP model on Atari games.

###  4.2 Design considerations & implications

In this subsection we discuss the possible variants of BoP enabled by its
ensemble characteristic, and which of them are conducive to diversity. To
understand these and the design decisions leading to the final settings
described above, we first explain the need for diversification and the off-
policy side of our framework.

Note that the following considerations for diversification are not entailed in
ensemble methods applied to scalar RL, where behaviour and model diversity is
much easier to uphold, i.e., a surprising backup target is considered unseen
in a scalar model but merely less likely in a distribution model. Thus said,
the same level of diversity would provide lesser learning signal to the
latter. This means if the posterior over return distributions is not
sufficiently diverse, which would happen if ensemble techniques and DiRL were
naively combined, the return distribution models will only tweak the
likelihood for the current backup target rather than fundamentally change its
return prediction; and the ensemble members remain similar after light update,
forming a vicious cycle.

With multiple heads that are being maintained in BoP, each head $i$ can form a
return estimate $g_{t}^{i}$ at the current timestep and produce a Bellman
target $x_{t}^{i}$ by sampling from itself for each of the future timesteps.
Thus, there are choices about sharing these targets or not across heads in two
domains: 1. to update the critic and 2. to compute the local curiosity bonus.

First we considered sharing the Bellman target $x_{t}^{k}$ across critics,
i.e., using the target of the current active head $k$ to update all critics.
We found that in practice sharing Bellman targets causes the diversity between
policies to collapse, as measured by various metrics such as the sum of all
absolute differences between the assigned probabilities of any two policies
and the cosine similarity between two policies. This resulted in poorer
exploration outcomes than with individual Bellman targets. Therefore, in the
final BoP framework, the on-policy Bellman target is not shared and each
critic head is trained against its own target, which provides (more) unique
bootstrapped signals to each head. Formally, we use the following variation of
the distributional Bellman equation:

$\displaystyle\mathcal{T}^{\pi^{k}}G^{\pi^{i}}(s_{t})\stackrel{{\scriptstyle
D}}{{=}}r_{t}+\gamma G^{\pi^{i}}(s_{t+1})$ (3) $\displaystyle s_{t+1}\sim
P(s_{t+1}|s_{t},a_{t}),\leavevmode\nobreak\ a_{t}\sim\pi^{k}(s_{t}).$

Since the states, actions and rewards are collected by the active head $k$,
only this head is trained on-policy, while all other heads are trained off-
policy. Thus, the final version of BoP is an algorithm of a mix of on-policy
and off-policy learning, in contrast to other algorithms Bellemare et al.
[2017], Mnih et al. [2016], Li and Faisal [2021], Osband et al. [2016], Dabney
et al. [2018b], which are either one or the other.

Each distributional actor-critic member in BoP has a curiosity bonus
$\hat{r}_{t}^{i}$ applied to its personal policy. We have discussed the reason
for sharing the posterior return estimate $x^{k}_{t}$ in Eq. (4.1) which is
for its accuracy, we here consider the choices for the prior $g_{t}$. This
prior can either be computed from that head itself
$\mathrm{KL}\big{(}q(\cdot|x_{t}^{k},s_{t})\ ||\
q(\cdot|g_{t}^{i},s_{t})\big{)}$ as in Eq. (4.1) or shared
$\mathrm{KL}\big{(}q(\cdot|x_{t}^{k},s_{t})\ ||\
q(\cdot|g_{t}^{k},s_{t})\big{)}$, where $q(z|x,s)$ is the variational encoder
for the return distribution. We use local return samples $g_{t}^{i}$ as prior
to compute the information gain. This is because a sample surprising to one
return distribution model might not surprise another, and each model should
keep track of its accuracy to correctly bias local action selection, so that
the model agreement on policies informs also implicitly the agreement on the
level of surprise.

Next, we consider the schedule for updating heads: we can update all of them
on every batch of freshly collected data, or just some of them. It is indeed
feasible to update only one head - e.g. the active one or the most uncertain
one, as measured by the curiosity bonuses. By updating only the most uncertain
head at each iteration, the agent maximizes the reduction in epistemic
uncertainty while still keeping the heads sufficiently diverse (by updating
only one head in each update iteration). However, we found in practice that
this diminishes sample efficiency due to limiting the number of heads being
trained with a given amount of data. Similar empirical findings were made when
training with other setups that selectively update members like updating a
randomly sampled head according to its epistemic uncertainty, or training each
head on 50% of their most uncertain samples. For that reason in BoP we update
all heads on all the data, which provides the fastest learning at no cost of
sacrificing sufficient diversity of the heads and therefore maximizes their
ensemble benefits on balance.

During testing, on the other hand, the question arises as to how action
selection should be performed now that we have multiple policies. BoP selects
an action by averaging the action distributions of all policy heads and then
chooses the greedy action from that distribution Wiering and van Hasselt
[2008]. Since all policies in our implementation are Gaussian, this is easy to
realize. In practice we observed more stable performance and higher scores
when using this “average-then-argmax” approach compared to other options, such
as each head picking their greedy action and then performing a majority vote
selection of the most popular / similar action across heads (“argmax-then-
vote”), or even selecting actions only from one randomly selected head
(“sample-then-argmax”). The benefits of the “average-then-argmax” tactic over
the others result from relying on model agreement which is what has been truly
learned and understood rather than uncertainty-quenching exploration.

Finally, we considered how to choose the number of ensemble members. We found
that a small number of $3$ to $5$ strikes a good balance between effective
exploration and training speed. For calibration of this statement, we found
that each additional head adds about $15\%$ more FLOPS to the total amount of
computation. This increase in performance from having many heads saturates as
policies start to converge and thus become more and more correlated, which
justifies using a small number of heads to increase performance and
exploration efficiency without increasing training time too much.

Figure 3: Comparison of BoP with one $(K=1)$ vs many heads $(K=3)$ or $(K=5)$
on Breakout and Pong. The shading shows one standard error from the mean.
Exponential smoothing was applied for the scores on Breakout. The numbers of
heads in multi-head cases shown are the minimal numbers that enable
significant improvement over one-head baseline.

##  5 Experimental results

We tested the performance of BoP on a diverse set of Atari environments, and
compared to three baselines with which BoP shares similarity in different
aspects: to Bootstrapped DQN Osband et al. [2016] as an ensemble-based
exploration in scalar RL, to BDPG Li and Faisal [2021] as a locally optimistic
DiRL exploration approach and the algorithm we built upon, and to A3C Mnih et
al. [2016] which facilitates covering extensive data space mainly by deploying
multiple workers in parallel. The set of Atari environments were chosen to
represent an exemplar mix of hard exploration environments, mazes, and shooter
games. Pixel-based observations were pre-processed by standard Atari wrappers
as laid out in Dabney et al. [2018b], including cropping the image frame,
using only grayscaled frames, frame stacking, and taking the maximum values of
any two consecutive frames.

First of all, we compare performances between BDPG and BoP on
$\mathrm{Freeway}$, $\mathrm{Hero}$ and $\mathrm{MsPacman}$, with learning
curves shown in Fig. 2. Compared to BDPG, BoP improves significantly on sample
efficiency, albeit not necessarily on asymptotic performance. We note that the
only difference between BDPG and BoP is that BoP is a bootstrapped ensemble of
BDPG, so the improvement of learning speed is in principle attributed to the
ensemble effect. However, at this stage, it cannot be decided it is the deep
exploration enabled by Thompson sampling or the sheer multitude of workers
that has resulted the advantage of ensemble technique.

Atari env | # | BoP | # | A3C | # | Boot
---|---|---|---|---|---|---
Freeway | 3 | 34 | 16 | 0.1 | 10 | 34
Breakout | 3 | 686 | 16 | 682 | 10 | 855
Hero | 3 | 37,728 | 16 | 32,464 | 10 | 21,021
MsPacman | 3 | 9,057 | 16 | 654 | 10 | 2,983
Qbert | 3 | 20,583 | 16 | 15,149 | 10 | 15,092
Alien | 5 | 5,470 | 16 | 518 | 10 | 2,436
Asterix | 5 | 42,500 | 16 | 22,141 | 10 | 19,713
Frostbite | 3 | 5,940 | 16 | 191 | 10 | 2,181
Amidar | 5 | 1,123 | 16 | 264 | 10 | 1,272
BeamRider | 3 | 6,499 | 16 | 22,708 | 10 | 23,429
StarGunner | 3 | 56,200 | 16 | 138,218 | 10 | 55,725
Seaquest | 3 | 2,027 | 16 | 2,355 | 10 | 9,083
Table 2: The testing results achieved across 3 runs. The number of heads (“#”)
for BoP were set by us, the ones for A3C and Bootstrapped DQN (“Boot”) were
set in their original works Mnih et al. [2016], Osband et al. [2016].

To this end, we move on to investigate the impact of the number of ensemble
members, i.e. heads, on a multi-worker agent’s performance. In Table 2, we
compare BoP against multi-worker / multi-head methods Bootrapped DQN and A3C.
Although BoP does not beat baselines on all selected environments, we can see
that it indeed outperforms the others most frequently. Crucially, the number
of heads for A3C ($K=16$) and Bootstrapped DQN ($K=10$), as proposed by their
original publications, are substantially higher than those needed for BoP
($K=3$ or in two cases $K=5$). Moreover, for all games where deep exploration
is essential (such as $\mathrm{MsPacman}$ and $\mathrm{Frostbite}$), BoP
substantially outperforms baselines by $50\%-200\%$. These promising results
suggest that increasing the multitude of workers alone is not key to better
performance. Combining with the observations made in comparing to BDPG, one
can establish that the strengths of BoP are attributed to deep exploration. In
the meantime, we also observed that BoP is not as competitive in games with
moving targets but not exploration-demanding, such as $\mathrm{BeamRider}$ and
$\mathrm{Seaquest}$. We surmise that this is because BoP became overly
explorative seeing the moving subjects whereas the motion could not be
accounted for by being bold unlike in the case of $\mathrm{Freeway}$.

Figure 4: Testing scores for BoP on Atari. We test the agent periodically
while training. Shaded areas indicate variation in performance across multiple
runs.

Another experiment that we did to investigate the ensemble effect is to
compare a one-head BoP against a few-head BoP (Fig. 3). We noted a substantial
improvement when the number of heads increased from one to even a very small
value (e.g. $3$), while being a qualitative transformation from a flat
algorithm into an ensemble one. This finding reinforces the notion that the
benefit of BoP is not due to merely more workers, but a fundamentally distinct
way that an ensemble method works in – by posterior sampling and therefore
favouring uncertain areas for exploration.

This performance improvement outweighs the increase in compute time for a few-
head BoP. This is probably due to the better exploration capability of each
head operating with an independently updated policy learning from its own
estimation errors so that an action is selected either because it is agreed to
be relatively certainly optimal or is uncertain. However, when the number of
ensemble members becomes too large, BoP will have lost its advantage. As
mentioned earlier, we observed that on Atari each additional head adds about
$15\%$ more FLOPS for the full duration of the training compared to having
only one head. This linear relationship holds true also when having $10$ or
$15$ heads. The marginal cost in terms of FLOPS is constant. The marginal
benefit from adding another head decreases with the number of heads quickly,
resulting in our lower use of heads than other multi-worker algorithms. In
conclusion, BoP is most advantageous when using a handful of members.

The learning curves of all Atari environments that BoP was tested on can be
found in Fig. 4. For instance, $\mathrm{Frostbite}$ is a game that needs very
specific sequences of actions, similar to maze games like $\mathrm{MsPacman}$
or $\mathrm{QBert}$ that require deep exploration, as decisions which one
takes early on have long term impacts.

##  6 Discussion & future work

The Bag of Policies (BoP) algorithm presented here is a multi-head
distributional estimator applicable to a large class of algorithms and thus
extends deep exploration to distributional RL settings. Our most salient
finding is that the benefits of the optimism in the face of uncertainty _can
accumulate_ for DiRL exploration, as we were able to substantially improve
upon a DiRL approach that is optimistic on per-state basis simply by making an
ensemble of it and performing Thompson sampling. This suggests that at least
in DiRL, Thompson sampling and curiosity bonus, although both facilitate
optimism, can work in conjunction to enhance performance further.

In addtion, as our experimental results showed, the advantage of BoP is not
attributed to the mere multitude of workers as BoP can surpass other multi-
worker / multi-head algorithms with much fewer heads, and to boost performance
considerably with a few-head BoP from a one-head counterpart. These
experiments substantiate that the ensemble technique which empowered BoP
functions in a fundamentally different fashion than naively summing up the
efforts of parallel workers – by Thompson sampling, which enables optimistic
and deep exploration in the same time.

On average, BoP achieves better learning speed and asymptotic performance than
baselines. We have observed the biggest improvement from baselines in maze-
like games like $\mathrm{MsPacman}$ or $\mathrm{QBert}$ where the agent has to
choose a path in a labyrinth while collecting various items scattered
throughout it and avoiding the enemies many moves down the line. These kinds
of environments require deep exploration and the agent can experience vastly
different outcomes depending on which path it takes. Hence, these environments
provide a good example where the BoP exploration capabilities can improve the
agent’s performance.

In a nutshell, BoP has demonstrated that not only deep exploration is viable
in DiRL, offering extensions to any DiRL settings (e.g. Bellemare et al.
[2017], Dabney et al. [2018b, a], Freirich et al. [2019], Doan et al. [2018],
Martin et al. [2020], Barth-Maron et al. [2018], Singh et al. [2020],
Kuznetsov et al. [2020], Choi et al. [2019]), but also that deep exploration
can further improve learning from an already optimistic exploration strategy.

#### Acknowledgments

We are grateful for our funding support. At the time of this work, GL and AF
are sponsored by UKRI Turing AI Fellowship (EP/V025449/1), LL and FV by the
PhD sponsorship of the Department of Computing, Imperial College London.

## References

  * An et al. [2021] G. An, S. Moon, J. Kim, and H. O. Song.  Uncertainty-based offline reinforcement learning with diversified q-ensemble.  In _Advances in Neural Information Processing Systems_ , volume 34, pages 7436–7447, 2021.
  * Barth-Maron et al. [2018] G. Barth-Maron, M. W. Hoffman, D. Budden, W. Dabney, D. Horgan, D. TB, A. Muldal, N. Heess, and T. Lillicrap.  Distributed distributional deterministic policy gradients.  In _Proceedings of the 6th International Conference on Learning Representations (ICLR)_ , 2018.
  * Bellemare et al. [2013] M. G Bellemare, Y. Naddaf, J. Veness, and M. Bowling.  The Arcade Learning Environment: An evaluation platform for general agents.  _Journal of Artificial Intelligence Research_ , 47:253–279, 2013.
  * Bellemare et al. [2017] M. G. Bellemare, W. Dabney, and R. Munos.  A distributional perspective on reinforcement learning.  _Proceedings of the 34th International Conference on Machine Learning_ , 70:449–458, 2017.
  * Chen et al. [2017] R. Y. Chen, S. Sidor, P. Abbeel, and J. Schulman.  UCB exploration via q-ensembles, 2017.
  * Choi et al. [2019] Y. Choi, K. Lee, and S. Oh.  Distributional deep reinforcement learning with a mixture of Gaussians.  In _2019 International Conference on Robotics and Automation (ICRA)_ , pages 9791–9797, 2019.
  * Dabney et al. [2018a] W. Dabney, G. Ostrovski, D. Silver, and R. Munos.  Implicit quantile networks for distributional reinforcement learning.  _Proceedings of the 35th International Conference on Machine Learning_ , 80:1096–1105, 2018a.
  * Dabney et al. [2018b] W. Dabney, M. Rowland, M. G. Bellemare, and R. Munos.  Distributional reinforcement learning with quantile regression.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 2018b.
  * Doan et al. [2018] T. Doan, B. Mazoure, and C. Lyle.  GAN q-learning, 2018.
  * Donahue et al. [2017] J. Donahue, P. Krähenbühl, and T. Darrell.  Adversarial feature learning.  _Proceedings of the 5th International Conference on Learning Representations (ICLR)_ , 2017.
  * Efron and Tibshirani [1994] B. Efron and R. J. Tibshirani.  _An introduction to the bootstrap._ CRC press, 1994.
  * Espeholt et al. [2018] Lasse Espeholt, Hubert Soyer, Remi Munos, Karen Simonyan, Vlad Mnih, Tom Ward, Yotam Doron, Vlad Firoiu, Tim Harley, Iain Dunning, Shane Legg, and Koray Kavukcuoglu.  IMPALA: Scalable distributed deep-RL with importance weighted actor-learner architectures.  In _Proceedings of the 35th International Conference on Machine Learning_ , volume 80, pages 1407–1416, Stockholmsmässan, Stockholm, Sweden, 2018\.
  * Freirich et al. [2019] D. Freirich, T. Shimkin, R. Meir, and A. Tamar.  Distributional multivariate policy evaluation and exploration with the bellman gan.  In _Proc. the 36th Intl. Conf. on Machine Learning (ICML)_ , volume 97, pages 1983–1992, Long Beach, CA, USA, 2019.
  * Kuznetsov et al. [2020] A. Kuznetsov, P. Shvechikov, A. Grishin, and D. Vetrov.  Controlling overestimation bias with truncated mixture of continuous distributional quantile critics.  In _Proceedings of the 37th International Conference on Machine Learning_ , 2020.
  * Li and Faisal [2021] L. Li and A. Faisal.  Bayesian distributional policy gradients.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 35(1):8429–8437, 2021.
  * Liang et al. [2018] Jacky Liang, Viktor Makoviychuk, Ankur Handa, Nuttapong Chentanez, Miles Macklin, and Dieter Fox.  Gpu-accelerated robotic simulation for distributed reinforcement learning.  In _Conference on Robot Learning_ , pages 270–282. PMLR, 2018.
  * Lyle et al. [2019] C. Lyle, M. G. Bellemare, and P. S. Castro.  A comparative analysis of expected and distributional reinforcement learning.  _Proceedings of the AAAI Conference on Artificial Intelligence_ , 33:4504–4511, 2019.
  * Martin et al. [2020] J. Martin, M. Lyskawinski, X. Li, and B. Englot.  Stochastically dominant distributional reinforcement learning.  In _Proceedings of the 37th International Conference on Machine Learning_ , 2020.
  * Mavrin et al. [2019] B. Mavrin et al.  Distributional reinforcement learning for efficient exploration.  _Proceedings of the 36th International Conference on Machine Learning,_ , 97:4424–4434, 2019.
  * Mnih et al. [2016] V.; Mnih et al.  Asynchronous methods for deep reinforcement learning.  _Proceedings of The 33rd International Conference on Machine Learning_ , 48:1928–1937, 2016.
  * O’Donoghue et al. [2018] B. O’Donoghue, I. Osband, R. Munos, and V. Mnih.  The uncertainty Bellman equation and exploration.  In _Proceedings of the 35th International Conference on Machine Learning_ , volume 80, pages 3839–3848, Stockholmsmässan, Stockholm Sweden, 10–15 Jul 2018.
  * Osband et al. [2016] I. Osband, C. Blundell, A. Pritzel, and B. Van Roy.  Deep exploration via Bootstrapped dqn.  _Advances in Neural Information Processing Systems_ , 29:4026–4034, 2016.
  * Osband et al. [2019] I. Osband, B. Van Roy, D. J. Russo, and Z. Wen.  Deep exploration via randomized value functions.  _Journal of Machine Learning Research_ , 20:1–62, 2019\.
  * Puterman [1994] M. L. Puterman.  _Markov Decision Processes: Discrete stochastic dynamic programming._ John Wiley & Sons, 1994.
  * Schulman et al. [2016] J. Schulman, P. Moritz, S. Levine, M. Jordan, and P. Abbeel.  High-dimensional continuous control using generalized advantage estimation.  _Proceedings of the 4th International Conference on Learning Representations (ICLR)_ , 2016.
  * Schulman et al. [2017] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov.  Proximal policy optimization algorithms.  _arXiv preprint arXiv:1707.06347_ , 2017.
  * Singh et al. [2020] R. Singh, K. Lee, and Y. Chen.  Sample-based distributional policy gradient, 2020.
  * Sutton [1999] R. S. Sutton.  Policy gradient methods for reinforcement learning with function approximation.  _Advances in Neural Information Processing Systems_ , 12:1057–1063, 1999.
  * Tang and Agrawal [2018] Y. Tang and S. Agrawal.  Exploration by distributional reinforcement learning.  _Proceedings of the 27th International Joint Conference on Artificial Intelligence_ , pages 2710–2716, 2018.
  * Thompson [1935] W. R. Thompson.  On the theory of apportionment.  _American Journal of Mathematics_ , 57(2):450–456, 1935.
  * Wiering and van Hasselt [2008] M. A. Wiering and H. P. van Hasselt.  Ensemble algorithms in reinforcement learning.  _IEEE Transactions on Systems, Man, and Cybernetics, Part B_ , 38(4):930–936, 2008.
  * Zhang et al. [2019] Zhizheng Zhang, Jiale Chen, Zhibo Chen, and Weiping Li.  Asynchronous episodic deep deterministic policy gradient: Toward continuous control in computationally complex environments.  _IEEE transactions on cybernetics_ , 2019.