11institutetext:  Instituto de Astrofísica de Canarias, Vía Láctea s/n,
E-38205 La Laguna, Tenerife, Spain 22institutetext: Departamento de
Astrofísica, Universidad de La Laguna, E-38200 La Laguna, Tenerife, Spain
33institutetext: Astrophysics Research Institute, Liverpool John Moores
University, 146 Brownlow Hill, Liverpool L3 5RF, UK

# BAYES-LOSVD: a bayesian framework for non-parametric extraction of the line-
of-sight velocity distribution of galaxies

J. Falcón-Barroso Email<EMAIL_ADDRESS>a bayesian framework for
non-parametric extraction of the line-of-sight velocity distribution of
galaxiesBAYES-LOSVD: a bayesian framework for non-parametric extraction of the
line-of-sight velocity distribution of galaxiesBAYES-LOSVD: a bayesian
framework for non-parametric extraction of the line-of-sight velocity
distribution of galaxiesBAYES-LOSVD: a bayesian framework for non-parametric
extraction of the line-of-sight velocity distribution of galaxies M. Martig
Email<EMAIL_ADDRESS>a bayesian framework for non-parametric
extraction of the line-of-sight velocity distribution of galaxiesBAYES-LOSVD:
a bayesian framework for non-parametric extraction of the line-of-sight
velocity distribution of galaxies

(Received October 9, 2020; accepted)

We introduce BAYES-LOSVD, a novel implementation of the non-parametric
extraction of line-of-sight velocity distributions (LOSVDs) in galaxies. We
employ bayesian inference to obtain robust LOSVDs and associated
uncertainties. Our method relies on principal component analysis to reduce the
dimensionality of the base of templates required for the extraction and thus
increase the performance of the code. In addition, we implement several
options to regularise the output solutions. Our tests, conducted on mock
spectra, confirm the ability of our approach to model a wide range of LOSVD
shapes, overcoming limitations of the most widely used parametric methods
(e.g. Gauss-Hermite expansion). We present examples of LOSVD extractions for
real galaxies with known peculiar LOSVD shapes, i.e. NGC 4371, IC 0719 and NGC
4550, using MUSE and SAURON integral-field unit (IFU) data. Our implementation
can also handle data from other popular IFU surveys (e.g. ATLAS3D, CALIFA,
MaNGA, SAMI). Details of the code and relevant documentation are freely
available to the community in the dedicated repositories.

###### Key Words.:

Methods: data analysis – Techniques: spectroscopic – Galaxies: general –
Galaxies: kinematics and dynamics – Galaxies: elliptical and lenticular, cD –
Galaxies: spiral

##  1 Introduction

Galaxies are made of stars that move on orbits with different degrees of
coherence around their nuclei. Each orbital group contains detailed
information about the assembly history of each component and thus retains some
memory of the different accretion events suffered by the galaxy over time.
This information is encoded in their line-of-sight velocity distribution
(LOSVDs), and thus by extracting this property we have access to vital clues
to unravel the formation and evolution of galaxies we see today. The analysis
of the LOSVD can be done directly in the Milky Way and in the Local Universe
by tracing the motions of stars with common stellar populations along a given
line-of-sight (e.g. Norris 1986; Tolstoy et al. 2004; Deason et al. 2011;
Kunder et al. 2012; Ness et al. 2013; Debattista et al. 2015; Zoccali et al.
2017; Du et al. 2020). This is a much more difficult task in nearby galaxies
beyond the Local Group, where stars are unresolved and thus the LOSVD at a
given position in a galaxy represents multiple populations along that line-of-
sight.

The extraction of the LOSVD of galaxies has been an active field of research
for many decades. Both parametric and non-parametric approaches have been
developed over the years to address this issue. Interestingly original
implementations prioritised non-parametric over parametric recoveries,
something that has radically changed in the last few decades. Most popular
approaches include: Fourier correlation quotient (FCQ) (e.g. Simkin 1974;
Sargent et al. 1977; Franx & Illingworth 1988; Bender 1990). Cross-correlation
(XC or CCF) (e.g. Tonry & Davis 1979; Statler 1995). Maximum penalised
likelihood (MPL) (e.g. Saha & Williams 1994; Merritt 1997; Pinkney et al.
2003). Direct fitting in pixel space (e.g. Rix et al. 1992; Kuijken &
Merrifield 1993; van der Marel & Franx 1993; Gebhardt et al. 2000; Kelson et
al. 2000; Cappellari & Emsellem 2004; Ocvirk et al. 2006; Chilingarian et al.
2007).

The extraction of the LOSVD is a degenerate problem, as there are an
indefinite number of combinations of stellar populations and LOSVD shapes that
can explain a particular spectroscopic observation. Breaking the degeneracies
implies a perfect knowledge of the underlying stellar populations that
contribute to a particular line-of-sight. In the last few decades, this issue
has been mitigated with the advent of a large number of intermediate-
resolution stellar libraries and stellar population models (see e.g. Bruzual &
Charlot 2003; Valdes et al. 2004; Coelho et al. 2005; Sánchez-Blázquez et al.
2006; Prugniel et al. 2007; Vazdekis et al. 2010; Gonneau et al. 2020;
Maraston et al. 2020), that have helped to greatly reduce the so-called effect
of template mismatch (e.g. Falcón-Barroso et al. 2003). Another important
aspect is the uncertainty of the recovered LOSVD. The methods cited above
handle this in different ways. Most of them deal with it with the aid of Monte
Carlo simulations by perturbing the input spectrum several times with some
known observed uncertainty (but see de Bruyne et al. 2003 for a more realistic
approach). While in recent years the use of parametric forms of the LOSVD has
prevailed in the literature (e.g. Emsellem et al. 2011; Falcón-Barroso et al.
2017; van de Sande et al. 2017), it is becoming increasingly clear that non-
parametric approaches are needed to describe the complexity in the LOSVD
shapes observed in real data (e.g. Jore et al. 1996; Kuijken et al. 1996;
Halliday et al. 2001; González-García et al. 2006; Katkov et al. 2011; Coccato
et al. 2013; Fabricius et al. 2014; Pizzella et al. 2018) and numerical
simulations (e.g. Jesseit et al. 2007; Martig et al. 2014; Schulze et al.
2017).

Bayesian inference methods (e.g. Hoffman & Gelman 2014) offer a natural way of
treating both: (1) the uncertainties in the fitting process, with the
advantage that they allow the inclusion of our knowledge of the problem during
the fitting process via priors on the input parameters, and (2) the true
complex, non-parametric nature of the LOSVDs. Saha & Williams (1994)
(hereafter SW94) already thought about this particular way of framing the
problem, but limitations in computer performance did not allow them to perform
a fully general optimisation of both the LOSVDs and templates. In this paper,
we revise the SW94 approach and use the latest developments on Bayesian
inference and dimensionality reduction techniques to develop a
Python111https://www.python.org/ implementation that can efficiently handle
simultaneously template optimisation and robust LOSVD extraction.

The paper is organised as follows. We describe the problem of LOSVD
extraction, the techniques for dimensionality reduction of the templates, and
LOSVD regularisation in section 2. We present our tests on mock spectra in §3,
and apply our extraction methods to real data in §5. We provide all the
necessary technical details of our python implementation in section 4, and
provide a summary of the paper and outlook for future applications of this
methodology in §6.

##  2 The LOSVD extraction

A LOSVD represents the distribution of the number of stars as a function of
velocity along a particular line-of-sight in a galaxy. In spectroscopic data
of nearby galaxies, the LOSVD is the broadening function to be applied to the
spectrum of the underlying stellar populations. In essence, the extraction of
the LOSVD is a deconvolution problem.

###  2.1 The equation

In its simplest form, the equation that describes the model to be fit to the
data can be expressed in mathematical terms as:

$G_{model}\,(\lambda)=\sum_{k=1}^{K}[w_{k}\cdot T_{k}(\lambda)]\star B+C(n),$
(1)

where $w_{k}$ are the weights for each stellar population template ($T_{k}$),
$B$ is the broadening function (i.e. the LOSVD), the $\star$ operator a
convolution, and $C(n)$ an additive polynomial of order $n$. The polynomial
term is convenient to reduce the impact of template mismatch (which can happen
even with the most complete template libraries), and other imperfections
during data reduction (e.g. non-perfect sky subtraction and/or scattered
light). The equation can become more complicated if the effects of dust
attenuation and/or calibration issues between data and templates are to be
taken into account (see equation 11 in Cappellari 2017a for a more complete
example). The implementation presented in this paper uses the prescription
shown in Eq. 1, as we have checked that it works for a wide range of datasets,
but it can be easily extended to include other correction terms if required.

The optimisation procedure involves the minimisation of the residuals:

$r_{p}=\frac{G_{data}(x_{p})-G_{model}(x_{p})}{\Delta G_{data}(x_{p})},$ (2)

where $x_{p}$ is the value of the data or model at a given pixel, $G_{data}$
is the observed spectrum, $\Delta G_{data}$ are the observed uncertainties,
and $G_{model}$ is the model presented in eq. 1. The most commonly used method
for the minimisation of eq. 2 is Least Squares, as there are plenty of
computer implementations (e.g. Lawson & Hanson 1974; Moré et al. 2001; Jones
et al. 2001) to perform the fits very efficiently.

Inspired by the work of SW94, we have opted for revisiting their bayesian non-
parametric approach for the LOSVD extraction. There were three main aspects of
SW94 work that were not possible to explore given the computer capabilities
and/or mathematical methods available at the time: (1) the Markov Chain Monte
Carlo (MCMC) sampling strategy, (2) template optimisation, (3) different forms
of regularisation for the LOSVD. We describe the new improvements in the
following sections.

Figure 1: Comparison of the quality of spectral fitting and LOSVD recovery for
different numbers of PCA components. Each column presents input spectra with
different stellar populations (young on the left, and intermediate/old
populations on the right). The input spectra have a S/N of 50 per pixel. All
panels are plotted on the same scale. Black line on main panels are the input
test data, while the red line shows the best fitting model. Residuals are
indicated in green. The spectral fits are carried out with 2, 5 and 10 PCA
templates (from top to bottom as indicated). In the insets, the input LOSVD is
a Gaussian centred at zero and a velocity dispersion of 150 km s-1 (indicated
in red). Units of the abscissae in the insets are km s-1. The recovered median
values of the LOSVDs are indicated with a thick black line. 16%$-$84% and
1%$-$99% confidence limits at each point are indicated in dark and light blue,
respectively. An order 2 auto-regressive prior was used to perform the fitting
(see §2.4 for details).

###  2.2 Markov Chain Monte Carlo sampling

There are multiple possible strategies for the sampling of parameter space in
bayesian inference frameworks. Classical approaches include Metropolis-
Hastings (Metropolis et al. 1953; Hastings 1970), or Gibbs (Geman & Geman
1984) samplers. With more complex models, the field has experienced a spur of
new methods to efficiently probe large parameter spaces: Nested sampling (e.g.
Buchner 2014), Hamiltonian Monte Carlo (e.g. Duane et al. 1987), or Stein
Variational Gradient Descent (e.g. Liu & Wang 2016) to cite a few. We refer
the interested reader to Chi Feng’s Github webpage222http://chi-
feng.github.io/mcmc-demo/ for a demo on the performance of different samplers.

The SW94 approach relied on the Metropolis algorithm for the exploration of
parameter space. Our implementation is based on the No-U-Turn-Sampler (NUTS)
introduced by Hoffman & Gelman (2014), a much more effective sampler. This is
part of the Stan333https://mc-stan.org/ package (Carpenter et al. 2017), our
software of choice for the MCMC sampling. Stan is a probabilistic programming
language for statistical modelling and data analysis used in many fields, from
physics or engineering to social sciences. We refer the interested reader to
Stenning et al. (2016); Asensio Ramos et al. (2017); Parviainen (2018);
Lamperti et al. (2019); Dullo et al. (2020) for a few examples of Stan
applications in astronomy.

Besides the sampling scheme, one of the main advantages of our implementation
in Stan is the use of the simplex (i.e. a vector of positive values whose sum
is equal to one) to describe the LOSVD. This type of parametrisation provides,
naturally, physical and normalisation constraints, and allows for a very
efficient exploration of parameter space during minimisation (Betancourt
2012).

###  2.3 Template optimisation

A crucial element in the recovery of the LOSVD is the basis of stellar
templates used to fit the observed spectra. Ideally, such basis should contain
stellar spectra covering the widest range possible of stellar parameters (i.e.
Teff, [Fe/H], log(g), and stellar abundances) or stellar population models
with a good sampling of, e.g., ages, metallicities, IMF shapes and slopes, and
possibly different chemical abundance ratios too. As already mentioned in the
introduction, this is now possible with the advent of latest the stellar
libraries and models (see §1 for details).

Using the typical set of $\sim$500-1000 templates, the minimisation involves
finding the weights ($w_{k}$) for each of them, and performing the convolution
of the LOSVD on $\sim$1000 spectral pixels per fitting iteration. This is a
very time consuming task even for the most efficient bayesian samplers. Since
reducing the number of pixels to fit may not be an option, we are thus left
with the only alternative of decreasing the number of templates to be used in
the optimisation process. This will have the advantage of not only reducing
the number of parameters (i.e. $w_{k}$), but also boosting minimisation
performance by very large factors (e.g. from days to minutes).

There are a number of well known techniques for dimensionality reduction that
one could use to whittle down the number of templates, e.g., Independent
Component Analysis (Lu et al. 2006), Factor analysis (Nolan et al. 2006), Non-
Negative Matrix factorization (Blanton & Roweis 2007), Diffusion Maps
(Richards et al. 2009), or $K$-means (Sánchez Almeida & Allende Prieto 2013).
Among all alternatives, we favoured the well-known principal component
analysis (PCA) method. PCA has been extensively used in astrophysics for
problems very much related to spectral fitting (e.g. Ronen et al. 1999; Li et
al. 2005; Chen et al. 2012) and has proven to be both robust and very
effective in speeding up calculations. In our particular problem the use of
PCA has improved performance from several days to a few minutes (i.e. a
500$\times$ fold decrease of computing time).

In practice the use of PCA means capturing well over 99% of the variance
present in a template library of $\sim$1000 spectra with less than 10 PCA
components. Figure 1 shows the effect of using different number of PCA
components during the fit for two mock spectra444For reference, for this
particular case, the continuum-to-line ratio (defined as 100
$\times\,\sigma_{residuals}/\sigma_{spectrum}$) is 13% and 23% for the young
and old populations respectively, almost independent of the number of PCA
components used in the fitting procedure. with signal-to-noise ratio per pixel
(hereafter S/N) of 50. It is remarkable how a number as low as 2 PCA
components can already reproduce with great accuracy the observed spectrum.
This is, in part, due to the use of additive polynomial that helps to reduce
potential template mismatch issues. We have checked the recovery for different
types of stellar populations (as shown in the two columns) and the extracted
LOSVD is less accurate for young stellar populations with a low number of PCA
components. This not unexpected given that spectra of younger populations show
greater variation, which is difficult to capture with just 2 PCA components.
We noticed that the closer the input spectrum is to the average spectra of the
templates (i.e. intermediate stellar populations in our particular case), the
better the recovery is for low number of PCA components. In general, our tests
suggest that 5 PCA components suffice to recover the input LOSVD with great
accuracy, with a clear improvement as S/N increases.

The optimal selection of the number of PCA components is not a well defined
quantity though. Choosing a number of PCA components that would explain a
certain variance of the spectra (e.g. compatible with the noise of the
spectrum) seems a reasonable approach. However, in practice, for real data,
this means selecting a different number of components for each spectrum of the
dataset, which is very impractical. In our experience, for typical S/N ratios
between 50 and 100, it is sufficient to choose a number between 5 and 10 PCA
components555This is using the MILES models (Vazdekis et al. 2010). In our
implementation, this number is a free parameter that the user has to decide on
before runtime.

###  2.4 LOSVD regularisation

The extraction of the LOSVD is a degenerate problem with a large number of
LOSVD shapes and combinations of templates that can reproduce the observed
spectrum with great accuracy. Nevertheless the range of allowed LOSVD
solutions can be constrained for reasonably high S/N data.

There are different ways to impose some level of regularisation onto the
output LOSVD. In our bayesian approach they are expressed in the form of
priors over the LOSVD. We have explored four different types of priors:

  1. 1.

No regularisation

$\mathrm{LOSVD}_{i}\sim\mathcal{N}(0,\sigma^{2}),$ (3)

where $\mathrm{LOSVD}_{i}$ is the i-th velocity element of the LOSVD, and
$\mathcal{N}(0,\sigma^{2})$ represents a normal distribution with mean zero
and variance $\sigma^{2}$. This is a fairly uninformative prior that assumes
the same level of uncertainty on all LOSVD elements.

  2. 2.

Random-Walk prior

$\mathrm{LOSVD}_{i}\sim\mathcal{N}(\mathrm{LOSVD}_{i-1},\sigma^{2}),$ (4)

where the i-th element of the LOSVD is linked to the previous one. This type
of prior was the one proposed by SW94. We set the prior for the first element
to: $\mathrm{LOSVD}_{0}\sim\mathcal{N}(0,\sigma^{2})$.

  3. 3.

Auto-regressive prior

$\mathrm{LOSVD}_{i}\sim\mathcal{N}\left(\alpha+\sum_{k=1}^{N}\beta_{k}\times\mathrm{LOSVD}_{i-k},\sigma^{2}\right).$
(5)

This is a more general form of prior that becomes equal to the Random-Walk for
$\alpha=0$ and $\beta=1$ for $k=1$. Stronger regularisation is obtained by
increasing the order $k$, which links more consecutive LOSVD elements. We
impose the same prior used for the Random-Walk approach to the first element
of the LOSVD, and use weakly informative normal priors for $\alpha$ and
$\beta$.

  4. 4.

Penalised B-splines

B-splines are a special type of piecewise polynomials controlled by knots that
are often used for interpolation (e.g. Press et al. 2003). One of the major
difficulties on the definition of B-splines is the choice of the number of
knots to define the polynomial function. In our Stan implementation, we
followed the example provided by Milad Kharratzadeh666https://mc-
stan.org/users/documentation/case-studies/splines_in_stan.html, and rather
than establishing a relation between different elements of the LOSVD, we apply
Random-Walk priors on the B-splines coefficients such that:

$a_{1}\sim\mathcal{N}(0,1),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
a_{i}\sim\mathcal{N}(a_{i-1},\tau),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \tau\sim\mathcal{N}(0,1).$ (6)

This approach promotes smoothness in the resulting LOSVD, preventing excessive
wiggling in the solution. Another important parameter that controls the level
of flexibility of the B-splines is the B-spline order $k$.

Figure 2: LOSVD recovery for no regularisation and an auto-regressive (order
2) prior. Colors as in Fig. 1

An interesting feature of our procedure is that we do not impose any
predefined value to $\sigma^{2}$ or $\tau$ during the fitting process. In
fact, they are considered nuisance parameters and we let the quality of the
data establish their optimal distributions. Far from being unconstrained, it
turns out that both parameters are well behaved and display fairly tight
distributions.

We illustrate the effect of the choice of prior for the LOSVD and its
uncertainties in Figure 2. We have created a test spectrum for this purpose
with a Gaussian LOSVD and a S/N=100 per spectral pixel. The figure shows the
difference between no regularisation and an auto-regressive prior of order 2.
Both approaches deliver an indistinguishable fit to the input spectrum, and
capture the Gaussian nature of the input LOSVD. However, the level of
uncertainty displayed by the case without regularisation is far larger than
that of the auto-regressive prior. These results are very much in agreement
with Saha & Williams (1994) findings. As we will show in Section 3, this
difference persists even at the highest S/N levels and is related to the
number of degrees of freedom of one method versus the other.

Figure 3: LOSVD recovery for different S/N ratios and types of regularisation.
We use a double Gaussian LOSVD profile as an example and it is indicated in
red. Left column displays solutions for extraction with no regularisation,
while right column shows results for an order 2 auto-regressive prior. Each
row represents two cases for a different S/N per pixel as indicated. All
panels are plotted on the same scale. Colors as in Fig. 2.

###  2.5 Likelihood

Besides the priors, for the minimisation process, we need to define the form
of the likelihood of our data given a model. In our problem we assume that our
spectroscopic observations can be explained by a normal distribution such
that:

$G_{data}\,(\lambda)\sim\mathcal{N}(G_{model}\,(\lambda),\leavevmode\nobreak\
\sigma_{G_{data}}^{2}),$ (7)

where $G_{data}\,(\lambda)$ is the observed input spectrum,
$G_{model}\,(\lambda)$ is our model, and $\sigma_{G_{data}}^{2}$ is the
variance of the observed spectrum.

Equation 1 is adequate to define our model when a full template library is to
be used during the minimisation process. In our case, however, the use of PCA
components transforms that equation to:

$G_{model}\,(\lambda)=\left[\widetilde{T}+\sum_{k=1}^{K}w_{k}\cdot
PCA_{k}(\lambda)\right]\star B+C(n),$ (8)

where now $w_{k}$ are the weights for each PCA component (PCAk),
$\widetilde{T}$ is the mean template of the input library, $B$ is the
broadening function (i.e. the LOSVD), the $\star$ operator a convolution, and
$C(n)$ an additive polynomial of order $n$. Our implementation in Stan
supports both forms of equations, but it is significantly more efficient with
Eq. 8.

##  3 Tests on simulated data

We have checked the performance of our implementation under different
circumstances by creating mock spectra for a wide range of S/N ratios and
input LOSVD shapes. While some of the LOSVDs were arbitrarily defined by
combining Gaussian functions, others come from numerical simulations and are
thus more realistic. We choose to illustrate here a case with an input
spectrum with an intermediate-age stellar population, but results are
consistent when other populations were used. Our fits were performed using 5
PCA components.

Figure 4: LOSVD recovery for different input LOSVDs and types of
regularisation. Colors as in Fig. 2. Each row represents a particular LOSVD
shape for different types of regularisation. These are solutions for input
spectra with S/N=50 per pixel. All panels are plotted on the same scale.

Figure 3 shows an example of the recovery of an extreme case of LOSVD shape
(i.e. a double Gaussian profile) for a range of S/N ratios and two prior
assumptions. The effect of regularisation is evident from the lowest to the
highest S/N ratios. With no regularisation there is a large range of possible
solutions, as displayed by the confidence intervals, even at S/N=200. It is
interesting to note that the level of uncertainty decreases drastically with
S/N when no regularisation is applied, while the improvement is not so large
when priors are used. Regularised solutions come with the price of introducing
some bias (see §5 for an example). This is most acute for the lowest S/N
ratios, where the regularised solution fails to capture the double peaked
nature of the input LOSVD. It is clear though that at S/N=10, even non-
regularised solutions cannot reproduce the input LOSVD. It seems that at S/N
of 25 the non-regularised solution clearly reveals already the double peak
feature, while the regularised one does not. This is an important result, as
it shows that non-regularised solutions are more accurate at low S/N regimes,
at the expense of larger uncertainties.

In Figure 4 we present how the different types of regularisation methods
listed in §2.4 influence the recovery of the range of input LOSVDs we prepared
for these tests. The S/N ratio of the input data is 50 per spectral pixel. The
first thing to notice is that the overall level of accuracy in the LOSVD
recovery is rather good (i.e. the input LOSVD is contained within the
confidence levels) for all LOSVD shapes. Nevertheless, as expected from
previous figures, the use of different priors result in different confidence
intervals. An interesting feature is the difference in the confidence
intervals between Random Walk and Auto-Regressive (order 1) priors. This is
entirely driven by the reduced number of degrees of freedom of the former over
the latter (see §2.4). Perhaps one of the most important messages from the
figure is that subtle details can be extracted at S/N=50. This is particularly
striking for the case displayed in the bottom row. This LOSVD is made of two
Gaussians: a prominent one centred at $-200$ km s-1, and a faint component
located at 500 km s-1. The goal of this test was to check the ability of the
code to recover the presence of potential small satellites being accreted onto
a galaxy. This appears to have been achieved for any choice of regularisation,
with slightly better recovery of the faint component with the Random-Walk and
Auto-Regressive priors. These results are quite encouraging, as it means that
there is not need for very high S/N ratios to detect this kind of features. We
provide versions of Fig. 4 for different S/N rations in Appendix A.

##  4 The implementation with pyStan

There are many alternatives one could choose from to implement all the ideas
proposed in section 2. Specifically in Python some of the most popular
packages include emcee777https://emcee.readthedocs.io/ (Foreman-Mackey et al.
2013) and pyMC3888https://docs.pymc.io/ (Salvatier et al. 2016), but see
Gabriel Perren’s webpage999https://gabriel-p.github.io/pythonMCMC/ for a
comprehensive list of options. We picked
pyStan101010https://pystan.readthedocs.io as our package of choice to develop
BAYES-LOSVD, as it offers a convenient python interface to Stan.

BAYES-LOSVD is built in a modular fashion to make it very simple for the end
user to extend its capabilities. This includes the addition of: (1) read-in
routines for data of a new instrument, (2) new Stan models with different
kinds of regularisation, (3) new template libraries. Our current
implementation supports datacubes from the MUSE-WFM, SAURON, ATLAS3D, CALIFA,
MaNGA, and SAMI surveys, besides the possibility to read standard 2D FITS
format files with spectra along rows. The linearly-sampled input spectra is
preprocessed before execution allowing the user to pick the level of Voronoi
binning (using Cappellari & Copin 2003 implementation), and velocity range and
sampling of the output LOSVD. In this process, the designated templates will
be prepared accordingly, with the possibility of switching off the reduction
of the template basis with the PCA scheme described in Sect. 2.3.

On execution, the user can decide whether to analyse the entire set of Voronoi
bins or a just a selection of them. Distributed computing is implemented
natively, so that multiple spectra can be executed in parallel on multi-CPU
machines. On output, by default, only summary statistics are stored on disk.
MCMC posterior distributions are described with highest density interval
estimators (e.g. Kruschke 2014), which are more accurate to describe highly
skewed distributions (as opposed to the standard percentiles approach). They
are stored on disk in HDF5 format111111https://www.h5py.org/. Diagnostic plots
are also created if requested. Users wanting to delve into the details can
chose to save the entire posterior distribution values, which can be then
easily analysed using the Arviz121212https://arviz-devs.github.io/arviz/
package.

Performance and execution times depend very much on the data, and parameters
used for the LOSVD extraction. Thanks to Stan though, convergence is usually
achieved with a very small number of iterations. In all the tests presented
here 3 chains with 500 iterations (i.e. warm-up + sampling) sufficed to obtain
well-behaved posterior distributions. A typical spectrum with $\sim$500 pixels
(e.g. 4800$-$5500 Å region) and S/N=50 per pixel, a velocity range of $\pm$
700 km s-1with a sampling of 50 km s-1, and 5 PCA components would require
$\sim$10 min on a cluster with Intel Xeon E5-2630 (v4) CPUs.

There is one major bottleneck in our implementation: convolution is performed
in direct space given Stan’s current inability to handle complex numbers. This
has a very strong impact when wide spectral ranges are to be fitted. Based on
discussions on the Stan forum131313https://discourse.mc-stan.org/ we are aware
that Fast Fourier Transforms will be possible in the not-so-distant future.
Another aspect where we could already optimise performance is the likelihood
evaluation. The latest version of Stan has introduced new features that can
speed this process by large factors by distributing its computation over many
CPUs. This is unfortunately not available in the pyStan version (v2.19) we are
using, but it will be available with the release of pyStan (v3). Therefore
there is still room for performance improvements in the midterm. We will be
alert and update the code to keep up with any new development.

The code can be accessed via the Github repository
https://github.com/jfalconbarroso/BAYES-LOSVD. Detailed documentation can be
found on the same page.

Figure 5: BAYES-LOSVD extraction for NGC 4371. (Top panel) Hubble Space
Telescope colour image based on the F475W and F850LP filters. North is up and
East to the left. (Bottom panels) LOSVDs extracted at different locations (as
indicated with black crosses on the image). On the left, without
regularisation, and on the right with an order 2 auto-regressive prior. Colors
as in Fig. 2. Red lines correspond to the best Gauss-Hermite LOSVD extracted
with pPXF (see §5.1 for details).

##  5 Application to real data

We now turn the attention to the recovery of LOSVD using real data from
different instruments. We show here examples of LOSVDs for three galaxies
presenting LOSVDs with varied degrees of complexity. We chose data from MUSE-
WFM and SAURON IFUs, but note that the code has been benchmarked with data
from some of the most popular IFU surveys (e.g. ATLAS3D, CALIFA, MaNGA, SAMI)
too.

###  5.1 NGC 4371

The first case we present is NGC 4371. This galaxy was studied in great detail
by Gadotti et al. (2015) using the MUSE IFU, and was part of the pilot
programme for the TIMER survey (Gadotti et al. 2019). NGC 4371 is interesting
for being a fairly old system (i.e. $\geq$ 7 Gyr throughout) with evidence for
a fossil nuclear stellar ring void of star formation (e.g. Erwin et al. 2001).
Following Gadotti et al. (2015) study, we extracted spectra with a
3.0″aperture in three positions of the galaxy: at the centre, the nuclear
stellar ring, and a location along the bar. The S/N of the spectra in the
three apertures is well above 100 per pixel. We impose a velocity sampling of
50 km s-1.

Figure 5 shows the results of our analysis. The three LOSVDs are different
from each other in different aspects. While the LOSVDs at the centre and bar
regions display fairly Gaussian profiles, the ring-dominated region is clearly
asymmetric and skewed. The panels on the left hand side show solutions with no
regularisation, and the ones on the right used an order 2 auto-regressive
prior. In essence, we see the same behaviour observed with the test data in
previous sections with non-regularised solutions giving larger confidence
intervals than the regularised ones. The two sets of solutions are very much
consistent with each other, as also observed in our experiments. These results
presented in the figure have been obtained applying our method to spectra in
the wavelength range 4800$-$5300Å. We also computed solutions based on spectra
around the Calcium Triplet region (8450$-$8700Å) and obtained identical
solutions (and thus not shown here). This not totally unexpected since there
is no evidence for multiple stellar populations in this galaxy.

In addition, for comparison, we plot the best Gauss-Hermite LOSVD extracted
with pPXF (Cappellari 2017b). In this particular case, the agreement between
the recovered LOSVDs with our method and pPXF is very good. This is specially
true for the almost Gaussian LOSVDs at the centre and bar locations of the
galaxy. At the ring, the regularised solution does not match the pPXF result
as closely as the non-regularised one, but differences are still within a
1%$-$99% percentiles of our non parametric extraction. While for cases like
this one, the advantage of the non-parametric approach may not seem evident,
it is important to note that the Gauss-Hermite parametrization allows for
negative values on the wings of the LOSVD. This situation occurs for h3 and h4
values of 0.1, $-0.1$ respectively, which are not uncommon in the kinematic
maps presented in many IFU surveys. Since we constrain the LOSVD during the
fit to admit only positive values, our method overcomes this limitation of and
provides naturally physically meaningful LOSVDs.

Figure 6: BAYES-LOSVD extraction for IC 0719. (Top panel) Sloan Digital Sky
Survey colour image based on the $g,r,i$ filters. The image has been rotated
so that the major axis is parallel to the abscissae. (Bottom panels) LOSVDs
extracted at different locations along the major axis of the galaxy (as
indicated with black crosses on the image). On the left column, without
regularisation, and on the right with an order 2 auto-regressive prior. Colors
as in Fig. 2. Red lines correspond to the best Gauss-Hermite LOSVD extracted
with pPXF (see §5.2 for details). Figure 7: BAYES-LOSVD extraction for NGC
4550. (Left panel) Hubble Space Telescope colour image based on the F555W and
F814W filters. North is up and East to the left. (Other panels) LOSVDs
extracted at several locations along the major axis of the galaxy (as marked
by the black crosses). Each row corresponds to a different location from the
centre of the galaxy (in arcsec) as shown on the first column. Each column
uses a different kind of regularisation, as indicated on the top row. Colors
as in Fig. 2. Red lines show Rix et al. (1992) LOSVDs at those locations.

###  5.2 IC 0719

The second case we study is IC 0719, a spectacular case displaying multiple
kinematic components (Pizzella et al. 2018). The galaxy is made of stars in
two counter-rotating large scale discs with distinct stellar populations. In
addition it has an ionised-gas component showing the same sense of rotation of
the secondary, lower-mass, younger stellar disc. We used MUSE observations
around the 4800$-$5300Å region to extract three apertures along the major axis
of the galaxy. As for the case of NGC 4371, the S/N of the spectra is well
above 100 per pixel. We impose a velocity sampling of 50 km s-1.

Figure 6 shows an almost perfect Gaussian LOSVD profile at the centre of the
galaxy, while LOSVDs along the major axis display clear double-peaked shapes.
This is in perfect agreement to the analysis performed by Pizzella et al.
(2018). For comparison, we also plot the best fit Gauss-Hermite LOSVD
parametrisation obtained with pPXF (in red). Here it becomes obvious that the
Gauss-Hermite expansion cannot reproduce well such complicated shapes and have
LOSVD wings that have slightly negative values. Although not obvious due to
the normalisation, for the aperture at 17″the pPXF extraction creates a third
peak in the LOSVD at velocities $\sim$500 km s-1, where the non-parametric
approach goes effectively to zero. It is worth highlighting the level of
complexity of the recovered LOSVDs despite the smooth morphological
appearance. The same applies to NGC 4550, as we discuss in Sec. 5.3, and
emphasizes the need for the non-parametric LOSVD extraction in galaxies. This
topic is gaining the attention in the literature and non-parametric LOSVDs are
now being routinely included in the dynamical modelling of early-type galaxies
(e.g. Mehrgan et al. 2019; Neureiter et al. 2020).

Another interesting aspect to explore in this galaxy is the non-parametric
extraction of the LOSVD in wavelength regions sensitive to different stellar
populations. This is actually possible with MUSE data and it will be a matter
of analysis in Rubino et al. (in preparation).

###  5.3 NGC 4550

The last case we analysed is NGC 4550, another classical showcase galaxy with
prominent double-peaked LOSVD profiles. We have extracted LOSVDs from SAURON
spectra (Emsellem et al. 2004). We performed our calculations in the
wavelength range between 4800$-$5300Å and a S/N=150 per spectral pixel. We
analysed the results along the major axis of the galaxy at three of the
locations presented in Rix et al. (1992). Since their data was not in
electronic form, we digitised it using the
WebPlotDigitiser141414https://automeris.io/WebPlotDigitizer tool (Marin et al.
2017) and then fitted the LOSVDs with double Gaussian profiles for best
reproduction. Figure 7 (top panel) shows the F555W/F814W color image from the
Hubble Space Telescope. Remaining panels show our recovered LOSVDs at the
three positions along the major axis, as indicated. Each row corresponds to a
location, while each column uses different priors for the LOSVD extraction: no
regularisation, auto-regressive (order 2), auto-regressive (order 1). The red
lines show the Rix et al. (1992) results.

The first thing to notice is the excellent agreement between our LOSVDs and
those of Rix et al. (1992) when no regularisation is used. The complexity of
the LOSVDs increases as we move away from the centre of the galaxy, becoming
double-peaked from 7.6″. The Rix et al. (1992) results are well within our
16%$-$84% confidence intervals (dark blue shaded region) despite the different
sampling in velocity.

As opposed to IC 079, the situation is drastically different when
regularisation is applied though. Our auto-regressive (order 2) solutions are
not capable of capturing the double-peaked nature of the LOSVDs at larger
distances from the centre. We have investigated the source for this
discrepancy and concluded that it is related to the intrinsic difference in
velocity between the two peaks and the velocity sampling used to extract the
LOSVDs. In other words, the level of correlation between velocity bins imposed
by this prior is too strong and smooths the solution too much. In order to
check this we have extracted the LOSVDs with an auto-regressive (order 1)
prior, but sampling the LOSVD in steps of 30 km s-1 instead of the 60 km s-1
used in all the other extractions. This is shown in the right-most column. It
clearly shows that the two peaks can be recovered with a less stringent prior
and finer sampling in velocity.

Based on these findings, we warn the reader to understand the implications of
using regularisation in their analysis. We therefore recommend potential users
to perform non-regularised fits to their data, and consider carefully the
velocity sampling to be used in the LOSVD extraction.

##  6 Summary & Conclusions

The advent of very high quality data from many integral-field spectrographs
and surveys has opened the possibility of efficiently extracting LOSVDs in
galaxies. At the same time, great progress in computer performance, algorithms
and mathematical methods, makes it possible now to handle large datasets.
Inspired by the work of SW94, we have developed a bayesian inference approach
to the LOSVD extraction from spectra. The code improves on SW94 work on 3 main
aspects: (1) the Markov Chain Monte Carlo sampling strategy, (2) the
possibility of different forms of regularization for the LOSVD, and (3)
template optimisation based on PCA components. Our tests on mock data indicate
that LOSVD recovery is accurate for spectra with S/N$>$50 with as few as 5 PCA
templates. Regularised solutions provide less uncertain LOSVDs, but it is at
the expense of biased solutions for low S/N ratios. We have also successfully
applied our approach to MUSE and SAURON data, displaying many interesting
features and warning the users to be careful with regularised solutions in
some situations (see §5.3 for an example). The use of non-regularised
solutions should therefore be preferred as it provides non-biased solutions.

On the technical side, our implementation is very versatile and allows the
possibility of extending its capabilities in different fronts (i.e. inclusion
of read-in routines for data of new instruments, new Stan models with
different kinds of regularisation, and/or addition of new template libraries.
The code and documentation can be downloaded from the repository indicated in
§4.

The complexity in the kinematics observed in IFU surveys (e.g. Krajnović et
al. 2011), but also numerical simulations (e.g. Martig et al. 2014; Schulze et
al. 2017; Walo-Martín et al. 2020), clearly indicates that a non-parametric
approach is necessary to capture the great level of detail that current data
offers. This has been shown for decades in early-type galaxies, but the
potential is much greater in late-type spiral systems which display much more
complex structures. In this respect, recent and upcoming large-scale IFU
facilities (e.g. VIRUS-W, Fabricius et al. 2008; MEGARA, Gil de Paz et al.
2018 and WEAVE-LIFU, Dalton et al. 2018) operating at spectral resolutions
above $R\geq 5\,000$ open the door to explore the details of the LOSVDs in low
velocity dispersion regimes (e.g. galaxy disks, dwarf galaxies), where it has
been very hard to operate with current instrumentation. The non-parametric
description of the LOSVDs will also have a big impact in the decomposition of
galaxies into their kinematic/dynamical components. Current efforts rely
heavily on the smooth LOSVDs provided by Gauss-Hermite parametrisations (e.g.
Tabor et al. 2017; Coccato et al. 2018; Oh et al. 2020). There is thus a great
potential to go beyond those (necessary) efforts to explore galaxy mass
assembly.

###### Acknowledgements.

We thank the referee Prasenjit Saha for very useful comments that have helped
improving the manuscript. We are grateful to Prashin Jethwa, Ignacio Martín-
Navarro, Marc Sarzi, Glenn van de Ven, and Eugene Vasiliev for many inspiring
discussions over the course of this project. We are also indebted to Michela
Rubino for her feedback and suggestions while testing earlier versions of the
code, and Alireza Molaeinezhad for assisting in the preparation of the
software package for Github. We thank Alessandro Pizzella for letting us use
the IC 0719 MUSE datacube in our analysis. J. F-B thanks Andrés Asensio Ramos
for introducing him to the world of bayesian inference methods and Stan in
particular. We thank Michele Cappellari for letting us include some of his
python functions in our software package. NGC 4371 results are based on
observations collected at the European Southern Observatory under ESO
programme 060.A-9313(A).
Software acknowledgements: Our code uses Astropy151515https://www.astropy.org/
a community-developed core Python package for Astronomy (Astropy Collaboration
et al. 2013; Price-Whelan et al. 2018), as well as
NumPy161616https://numpy.org/ (Oliphant 2006), SciPy171717https://scipy.org/
(Jones et al. 2001) and Matplotlib181818https://matplotlib.org/ (Hunter 2007)
Funding and financial support acknowledgements: J. F-B acknowledges support
through the RAVET project by the grant PID2019-107427GB-C32 from the Spanish
Ministry of Science, Innovation and Universities (MCIU), and through the IAC
project TRACES which is partially supported through the state budget and the
regional budget of the Consejería de Economía, Industria, Comercio y
Conocimiento of the Canary Islands Autonomous Community.

## References

  * Asensio Ramos et al. (2017) Asensio Ramos, A., de la Cruz Rodríguez, J., Martínez González, M. J., & Socas-Navarro, H. 2017, A&A, 599, A133
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Bender (1990) Bender, R. 1990, A&A, 229, 441
  * Betancourt (2012) Betancourt, M. 2012, in American Institute of Physics Conference Series, Vol. 1443, American Institute of Physics Conference Series, ed. P. Goyal, A. Giffin, K. H. Knuth, & E. Vrscay, 157–164
  * Blanton & Roweis (2007) Blanton, M. R. & Roweis, S. 2007, AJ, 133, 734
  * Bruzual & Charlot (2003) Bruzual, G. & Charlot, S. 2003, MNRAS, 344, 1000
  * Buchner (2014) Buchner, J. 2014, arXiv e-prints, arXiv:1407.5459
  * Cappellari (2017a) Cappellari, M. 2017a, MNRAS, 466, 798
  * Cappellari (2017b) Cappellari, M. 2017b, MNRAS, 466, 798
  * Cappellari & Copin (2003) Cappellari, M. & Copin, Y. 2003, MNRAS, 342, 345
  * Cappellari & Emsellem (2004) Cappellari, M. & Emsellem, E. 2004, PASP, 116, 138
  * Carpenter et al. (2017) Carpenter, B., Gelman, A., Hoffman, M., et al. 2017, Journal of Statistical Software, Articles, 76, 1
  * Chen et al. (2012) Chen, Y.-M., Kauffmann, G., Tremonti, C. A., et al. 2012, MNRAS, 421, 314
  * Chilingarian et al. (2007) Chilingarian, I., Prugniel, P., Sil’Chenko, O., & Koleva, M. 2007, in IAU Symposium, Vol. 241, Stellar Populations as Building Blocks of Galaxies, ed. A. Vazdekis & R. Peletier, 175–176
  * Coccato et al. (2018) Coccato, L., Fabricius, M. H., Saglia, R. P., et al. 2018, MNRAS, 477, 1958
  * Coccato et al. (2013) Coccato, L., Morelli, L., Pizzella, A., et al. 2013, A&A, 549, A3
  * Coelho et al. (2005) Coelho, P., Barbuy, B., Meléndez, J., Schiavon, R. P., & Castilho, B. V. 2005, A&A, 443, 735
  * Dalton et al. (2018) Dalton, G., Trager, S., Abrams, D. C., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, 107021B
  * de Bruyne et al. (2003) de Bruyne, V., Vauterin, P., de Rijcke, S., & Dejonghe, H. 2003, MNRAS, 339, 215
  * Deason et al. (2011) Deason, A. J., Belokurov, V., & Evans, N. W. 2011, MNRAS, 411, 1480
  * Debattista et al. (2015) Debattista, V. P., Ness, M., Earp, S. W. F., & Cole, D. R. 2015, ApJ, 812, L16
  * Du et al. (2020) Du, H., Mao, S., Athanassoula, E., Shen, J., & Pietrukowicz, P. 2020, arXiv e-prints, arXiv:2007.01102
  * Duane et al. (1987) Duane, S., Kennedy, A. D., Pendleton, B. J., & Roweth, D. 1987, Physics Letters B, 195, 216
  * Dullo et al. (2020) Dullo, B. T., Bouquin, A. Y. K., Gil de Paz, A., Knapen, J. H., & Gorgas, J. 2020, ApJ, 898, 83
  * Emsellem et al. (2011) Emsellem, E., Cappellari, M., Krajnović, D., et al. 2011, MNRAS, 414, 888
  * Emsellem et al. (2004) Emsellem, E., Cappellari, M., Peletier, R. F., et al. 2004, MNRAS, 352, 721
  * Erwin et al. (2001) Erwin, P., Vega Beltrán, J. C., & Beckman, J. E. 2001, in Astronomical Society of the Pacific Conference Series, Vol. 249, The Central Kiloparsec of Starbursts and AGN: The La Palma Connection, ed. J. H. Knapen, J. E. Beckman, I. Shlosman, & T. J. Mahoney, 171
  * Fabricius et al. (2008) Fabricius, M. H., Barnes, S., Bender, R., et al. 2008, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 7014, Ground-based and Airborne Instrumentation for Astronomy II, 701473
  * Fabricius et al. (2014) Fabricius, M. H., Coccato, L., Bender, R., et al. 2014, MNRAS, 441, 2212
  * Falcón-Barroso et al. (2003) Falcón-Barroso, J., Balcells, M., Peletier, R. F., & Vazdekis, A. 2003, A&A, 405, 455
  * Falcón-Barroso et al. (2017) Falcón-Barroso, J., Lyubenova, M., van de Ven, G., et al. 2017, A&A, 597, A48
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Franx & Illingworth (1988) Franx, M. & Illingworth, G. D. 1988, ApJ, 327, L55
  * Gadotti et al. (2019) Gadotti, D. A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2019, MNRAS, 482, 506
  * Gadotti et al. (2015) Gadotti, D. A., Seidel, M. K., Sánchez-Blázquez, P., et al. 2015, A&A, 584, A90
  * Gebhardt et al. (2000) Gebhardt, K., Richstone, D., Kormendy, J., et al. 2000, AJ, 119, 1157
  * Geman & Geman (1984) Geman, S. & Geman, D. 1984, IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMI-6, 721
  * Gil de Paz et al. (2018) Gil de Paz, A., Carrasco, E., Gallego, J., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, 1070217
  * Gonneau et al. (2020) Gonneau, A., Lyubenova, M., Lançon, A., et al. 2020, A&A, 634, A133
  * González-García et al. (2006) González-García, A. C., Balcells, M., & Olshevsky, V. S. 2006, MNRAS, 372, L78
  * Halliday et al. (2001) Halliday, C., Davies, R. L., Kuntschner, H., et al. 2001, MNRAS, 326, 473
  * Hastings (1970) Hastings, W. K. 1970, Biometrika, 57, 97
  * Hoffman & Gelman (2014) Hoffman, M. D. & Gelman, A. 2014, Journal of Machine Learning Research, 15, 1593
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Jesseit et al. (2007) Jesseit, R., Naab, T., Peletier, R. F., & Burkert, A. 2007, MNRAS, 376, 997
  * Jones et al. (2001) Jones, E., Oliphant, T., Peterson, P., & et al. 2001, SciPy: Open source scientific tools for Python (http://www.scipy.org/)
  * Jore et al. (1996) Jore, K. P., Broeils, A. H., & Haynes, M. P. 1996, AJ, 112, 438
  * Katkov et al. (2011) Katkov, I., Chilingarian, I., Sil’chenko, O., Zasov, A., & Afanasiev, V. 2011, Baltic Astronomy, 20, 453
  * Kelson et al. (2000) Kelson, D. D., Illingworth, G. D., van Dokkum, P. G., & Franx, M. 2000, ApJ, 531, 159
  * Krajnović et al. (2011) Krajnović, D., Emsellem, E., Cappellari, M., et al. 2011, MNRAS, 414, 2923
  * Kruschke (2014) Kruschke, J. 2014, Doing Bayesian Data Analysis: A Tutorial with R, JAGS, and Stan (Academic Press)
  * Kuijken et al. (1996) Kuijken, K., Fisher, D., & Merrifield, M. R. 1996, MNRAS, 283, 543
  * Kuijken & Merrifield (1993) Kuijken, K. & Merrifield, M. R. 1993, MNRAS, 264, 712
  * Kunder et al. (2012) Kunder, A., Koch, A., Rich, R. M., et al. 2012, AJ, 143, 57
  * Lamperti et al. (2019) Lamperti, I., Saintonge, A., De Looze, I., et al. 2019, MNRAS, 489, 4389
  * Lawson & Hanson (1974) Lawson, C. L. & Hanson, R. J. 1974, Solving least squares problems
  * Li et al. (2005) Li, C., Wang, T.-G., Zhou, H.-Y., Dong, X.-B., & Cheng, F.-Z. 2005, AJ, 129, 669
  * Liu & Wang (2016) Liu, Q. & Wang, D. 2016, arXiv e-prints [arXiv:1608.04471]
  * Lu et al. (2006) Lu, H., Zhou, H., Wang, J., et al. 2006, AJ, 131, 790
  * Maraston et al. (2020) Maraston, C., Hill, L., Thomas, D., et al. 2020, MNRAS, 496, 2962
  * Marin et al. (2017) Marin, F., Rohatgi, A., & Charlot, S. 2017, in SF2A-2017: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, ed. C. Reylé, P. Di Matteo, F. Herpin, E. Lagadec, A. Lançon, Z. Meliani, & F. Royer, Di
  * Martig et al. (2014) Martig, M., Minchev, I., & Flynn, C. 2014, MNRAS, 443, 2452
  * Mehrgan et al. (2019) Mehrgan, K., Thomas, J., Saglia, R., et al. 2019, ApJ, 887, 195
  * Merritt (1997) Merritt, D. 1997, AJ, 114, 228
  * Metropolis et al. (1953) Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., & Teller, E. 1953, J. Chem. Phys., 21, 1087
  * Moré et al. (2001) Moré, J., Garbow, B., & Hillstrom, K. 2001, User guide for MINPACK-1. Argonne National Laboratory Argonne, IL (http://cds.cern.ch/record/126569)
  * Ness et al. (2013) Ness, M., Freeman, K., Athanassoula, E., et al. 2013, MNRAS, 432, 2092
  * Neureiter et al. (2020) Neureiter, B., Thomas, J., Saglia, R., et al. 2020, arXiv e-prints, arXiv:2009.08979
  * Nolan et al. (2006) Nolan, L. A., Harva, M. O., Kabán, A., & Raychaudhury, S. 2006, MNRAS, 366, 321
  * Norris (1986) Norris, J. 1986, ApJS, 61, 667
  * Ocvirk et al. (2006) Ocvirk, P., Pichon, C., Lançon, A., & Thiébaut, E. 2006, MNRAS, 365, 74
  * Oh et al. (2020) Oh, S., Colless, M., Barsanti, S., et al. 2020, MNRAS, 495, 4638
  * Oliphant (2006) Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (Trelgol Publishing USA)
  * Parviainen (2018) Parviainen, H. 2018, Bayesian Methods for Exoplanet Science, 149
  * Pinkney et al. (2003) Pinkney, J., Gebhardt, K., Bender, R., et al. 2003, ApJ, 596, 903
  * Pizzella et al. (2018) Pizzella, A., Morelli, L., Coccato, L., et al. 2018, A&A, 616, A22
  * Press et al. (2003) Press, W. H., Teukolsky, S. A., Vettering, W. T., & Flannery, B. P. 2003, European Journal of Physics, 24, 329
  * Price-Whelan et al. (2018) Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018, AJ, 156, 123
  * Prugniel et al. (2007) Prugniel, P., Soubiran, C., Koleva, M., & Le Borgne, D. 2007, arXiv e-prints, astro
  * Richards et al. (2009) Richards, J. W., Freeman, P. E., Lee, A. B., & Schafer, C. M. 2009, ApJ, 691, 32
  * Rix et al. (1992) Rix, H.-W., Franx, M., Fisher, D., & Illingworth, G. 1992, ApJ, 400, L5
  * Ronen et al. (1999) Ronen, S., Aragon-Salamanca, A., & Lahav, O. 1999, MNRAS, 303, 284
  * Saha & Williams (1994) Saha, P. & Williams, T. B. 1994, AJ, 107, 1295
  * Salvatier et al. (2016) Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. 2016, PeerJ Comput. Sci., 2, 55
  * Sánchez Almeida & Allende Prieto (2013) Sánchez Almeida, J. & Allende Prieto, C. 2013, ApJ, 763, 50
  * Sánchez-Blázquez et al. (2006) Sánchez-Blázquez, P., Peletier, R. F., Jiménez-Vicente, J., et al. 2006, MNRAS, 371, 703
  * Sargent et al. (1977) Sargent, W. L. W., Schechter, P. L., Boksenberg, A., & Shortridge, K. 1977, ApJ, 212, 326
  * Schulze et al. (2017) Schulze, F., Remus, R.-S., & Dolag, K. 2017, Galaxies, 5, 41
  * Simkin (1974) Simkin, S. M. 1974, A&A, 31, 129
  * Statler (1995) Statler, T. 1995, AJ, 109, 1371
  * Stenning et al. (2016) Stenning, D. C., Wagner-Kaiser, R., Robinson, E., et al. 2016, ApJ, 826, 41
  * Tabor et al. (2017) Tabor, M., Merrifield, M., Aragón-Salamanca, A., et al. 2017, MNRAS, 466, 2024
  * Tolstoy et al. (2004) Tolstoy, E., Irwin, M. J., Helmi, A., et al. 2004, ApJ, 617, L119
  * Tonry & Davis (1979) Tonry, J. & Davis, M. 1979, AJ, 84, 1511
  * Valdes et al. (2004) Valdes, F., Gupta, R., Rose, J. A., Singh, H. P., & Bell, D. J. 2004, ApJS, 152, 251
  * van de Sande et al. (2017) van de Sande, J., Bland-Hawthorn, J., Fogarty, L. M. R., et al. 2017, ApJ, 835, 104
  * van der Marel & Franx (1993) van der Marel, R. P. & Franx, M. 1993, ApJ, 407, 525
  * Vazdekis et al. (2010) Vazdekis, A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2010, MNRAS, 404, 1639
  * Walo-Martín et al. (2020) Walo-Martín, D., Falcón-Barroso, J., Dalla Vecchia, C., Pérez, I., & Negri, A. 2020, MNRAS, 494, 5652
  * Zoccali et al. (2017) Zoccali, M., Vasquez, S., Gonzalez, O. A., et al. 2017, A&A, 599, A12

##  Appendix A Trends with S/N ratio

In this appendix we present the equivalent figures to Fig. 4 for different S/N
ratios. It is evident that the LOSVD recovery worsens as S/N decreases. It is
interesting to see that not all the features are captured well even at S/N=100
(e.g. small bump on the positive wing of the case in the bottom row in all
figures).

Figure 8: Same as Fig. 4 for S/N=10. Figure 9: Same as Fig. 4 for S/N=25.
Figure 10: Same as Fig. 4 for S/N=100. Figure 11: Same as Fig. 4 for S/N=200.


# Humanoid-Gym: Reinforcement Learning for Humanoid Robot with Zero-Shot
Sim2Real Transfer

Xinyang Gu2∗, Yen-Jen Wang13∗, Jianyu Chen123 ∗Equal contribution. Listed
alphabetically.1Shanghai Qi Zhi Institute, Shanghai, China.2RobotEra
TECHNOLOGY CO., LTD., Beijing, China3Tsinghua University, Beijing, China.

###### Abstract

Humanoid-Gym is an easy-to-use reinforcement learning (RL) framework based on
Nvidia Isaac Gym, designed to train locomotion skills for humanoid robots,
emphasizing zero-shot transfer from simulation to the real-world environment.
Humanoid-Gym also integrates a sim-to-sim framework from Isaac Gym to Mujoco
that allows users to verify the trained policies in different physical
simulations to ensure the robustness and generalization of the policies. This
framework is verified by RobotEra’s XBot-S (1.2-meter tall humanoid robot) and
XBot-L (1.65-meter tall humanoid robot) in a real-world environment with zero-
shot sim-to-real transfer. The project website and source code can be found
at: sites.google.com/view/humanoid-gym.

##  I INTRODUCTION

Modern environments are primarily designed for humans. Therefore, humanoid
robots, with their human-like skeletal structure, are especially suited for
tasks in human-centric environments, offering unique advantages over other
types of robots. Recently, massively parallel deep reinforcement learning (RL)
in simulation has become popular [1, 2, 3]. However, due to the complex
structure of humanoid robots, the sim-to-real gap [4, 5, 6, 7] exists and is
larger than that of quadrupedal robots. Therefore, we release Humanoid-Gym, an
easy-to-use RL framework based on Nvidia Isaac Gym [8], designed to train
locomotion skills for humanoid robots, emphasizing zero-shot transfer from
simulation to the real-world environment. Humanoid-Gym features specialized
rewards and domain randomization techniques for humanoid robots, simplifying
the difficulty of sim-to-real transfer. Furthermore, it also integrates a sim-
to-sim framework from Isaac Gym [8] to MuJoCo [9] that allows users to verify
the trained policies in different physical simulations to ensure the
robustness and generalization of the policies, shown in Fig. 1. Currently,
Humanoid-Gym is verified by multiple humanoid robots with different sizes in a
real-world environment with zero-shot sim-to-real transfer, including
RobotEra’s XBot-S (1.2-meter tall humanoid robot) and XBot-L (1.65-meter tall
humanoid robot) [10]. The contribution of Humanoid-Gym can be summarized as
follows:

Figure 1: Humanoid-Gym enables users to train their policies within Nvidia
Isaac Gym and validate them in MuJoCo. Additionally, we have successfully
tested the complete pipeline with two humanoid robots. They were trained in
Humanoid-Gym and transferred to real-world environments in a zero-shot manner.

  * •

We launch an open-source reinforcement learning (RL) framework with meticulous
system design.

  * •

Our framework enables zero-shot transfer from simulation to the real world,
which has been rigorously tested across humanoid robots of various sizes.

  * •

Our open-source library features a sim-to-sim validation tool, enabling users
to test their policies across diverse environmental dynamics rigorously.

##  II Related Works

###  II-A Robot Learning on Locomotion Tasks

Reinforcement learning (RL) has shown promise in enabling robots to achieve
stable locomotion [6, 11, 12]. Compared to prior RL efforts with quadrupedal
robots [1, 13] and bipedal robots like Cassie [14, 15], our work with humanoid
robots introduces a more challenging scenario for robot control. Recent
studies [16, 17] have applied transformer architecture to improve the walking
performance of humanoid robots on flat surfaces. Beyond lower-body control,
some works [18, 19] have also explored more complex upper-body skills for
humanoid robot control. However, the sim-to-real transition for humanoid
locomotion remains a significant challenge, with a notable lack of open-source
resources in the robot learning community. To contribute to this area, we have
developed Humanoid-Gym, an accessible framework with full codebase.

##  III Method

The workflow of Humanoid-Gym is illustrated in Fig. 2. In this section, we
will introduce the problem setting, system design, and reward design of our
Humanoid-Gym.

Figure 2: Pipeline of Humanoid-Gym. Initially, we employ massively parallel
deep reinforcement learning (RL) within Nvidia Isaac Gym, incorporating
diverse terrains and dynamics randomization. Subsequently, we undertake sim-
to-sim transfer to test policies. Due to our meticulous calibration, the
performance in both MuJoCo and real-world settings aligns closely.

###  III-A Reinforcement Learning For Robot Control

Our approach employs a reinforcement learning model
$\mathcal{M}=\langle\mathcal{S},\mathcal{A},T,\mathcal{O},R,\gamma\rangle$,
with $\mathcal{S}$ and $\mathcal{A}$ denoting state and action spaces,
$T(\mathbf{s}^{\prime}|\mathbf{s},\mathbf{a})$ the transition dynamics,
$R(\mathbf{s},\mathbf{a})$ the reward function, $\gamma\in[0,1]$ the discount
factor, and $\mathcal{O}$ the observation space. The framework is designed for
both simulated and real-world settings, transitioning from full observability
in simulations $(\mathbf{s}\in\mathcal{S})$ to partial observability in the
real world ($\mathbf{o}\in\mathcal{O}$). This necessitates operating within a
Partially Observable Markov Decision Process (POMDP) [20], with the policy
$\pi(\mathbf{a}|\mathbf{o}_{\leq t})$ mapping observations to action
distributions to maximize the expected return
$J=\mathbb{E}[R_{t}]=\mathbb{E}\left[\sum_{t}\gamma^{t}r_{t}\right]$.

We leverage Proximal Policy Optimization (PPO) [21] loss, supplemented by the
Asymmetric Actor Critic [22] method and the integration of privileged
information during training, shifting to partial observations during
deployment. The policy loss is defined as:

$\displaystyle\mathcal{L}_{\pi}$ $\displaystyle=\min\left[\frac{\pi(a_{t}\mid
o_{\leq t})}{\pi_{b}(a_{t}\mid o_{\leq t})}A^{\pi_{b}}(o_{\leq
t},a_{t}),\right.$ (1)
$\displaystyle\quad\left.\text{clip}\left(\frac{\pi(a_{t}\mid o_{\leq
t})}{\pi_{b}(a_{t}\mid o_{\leq t})},c_{1},c_{2}\right)A^{\pi_{b}}(o_{\leq
t},a_{t})\right]$

Advantage estimation utilizes Generalized Advantage Estimation (GAE) [23],
requiring an updated value function:

$\mathcal{L}_{v}=\|R_{t}-V(s_{t})\|_{2},$ (2)

###  III-B System design

The base poses of the robot, denoted as $P^{b}$, are six-dimensional vectors
$[x,y,z,\alpha,\beta,\gamma]$, representing both the position coordinates
$x,y,z$ and the orientation angles $\alpha,\beta,\gamma$ in Euler notation.
The joint position for each motor is represented by $\theta$, and the
corresponding joint velocity by $\dot{\theta}$. Furthermore, we define a gait
phase [24, 25], which comprises two double support phases (DS) and two single
support phases (SS) within each gait cycle. The cycle time, denoted as
$C_{T}$, is the duration of one full gait cycle. A sinusoidal wave is employed
to generate reference motion, reflecting the repetitive nature of the gait
cycle involving pitch, knee, and ankle movements. Notably, we also designed a
periodic stance mask $I_{p}(t)$ (Fig 6) that indicates foot contact patterns
in synchronization with the reference motion. For instance, if the reference
motion lifts the left foot, the right foot should be in the single support
phase, with the foot contact mask indicated as $[0,1]$; during DS phases, it
would be $[1,1]$.

The chosen action is the target joint position for the Proportional-Derivative
(PD) controller. The policy network integrates proprioceptive sensor data, a
periodic clock signal $[\sin(2\pi t/C_{T}),\cos(2\pi t/C_{T})]$, and velocity
commands $\dot{P}_{x,y,\gamma}$. A single frame of input are elaborated in
Table I. Additionally, the state frame includes feet contact detect $I_{d}(t)$
and other privileged observations.

Table I: Summary of Observation Space. The table categorizes the components of the observation space into observation and state. The table also details their dimensions. Components | Dims | Observation | State
---|---|---|---
Clock Input $(\sin(t),\cos(t))$ | 2 | ✓ | ✓
Commands $(\dot{P}_{x,y,\gamma})$ | 3 | ✓ | ✓
Joint Position $(\theta)$ | 12 | ✓ | ✓
Joint Velocity $(\dot{\theta})$ | 12 | ✓ | ✓
Angular Velocity $(\dot{P}^{b}_{\alpha\beta\gamma})$ | 3 | ✓ | ✓
Euler Angle $(P^{b}_{\alpha\beta\gamma})$ | 3 | ✓ | ✓
Last Actions $(a_{t-1})$ | 12 | ✓ | ✓
Frictions | 1 |  | ✓
Body Mass | 1 |  | ✓
Base Linear Velocity | 3 |  | ✓
Push Force | 2 |  | ✓
Push Torques | 3 |  | ✓
Tracking Difference | 12 |  | ✓
Periodic Stance Mask | 2 |  | ✓
Feet Contact detection | 2 |  | ✓

Our control policy operates at a high frequency of $100$Hz, providing enhanced
granularity and precision beyond standard RL locomotion approaches. The
internal PD controller runs at an even higher frequency of $1000$Hz. For
training simulations, Isaac Gym is utilized [8], while MuJoCo, known for its
accurate physical dynamics, is chosen for sim2sim validation. This approach
combines the benefits of high-speed GPU-based parallel simulation, albeit with
less accuracy, and the high accuracy but slower CPU-based simulation.

The detailed settings for both algorithms and the environment designed are
shown in Appendix TABLE II. We use multi-frames of observations and privilege
observation, which is crucial for locomotion tasks on uneven terrain.

###  III-C Reward Design

Our reward function directs the robot to adhere to velocity commands, sustain
a stable gait, and achieve smooth contact. The reward function is structured
into four key components: (1) velocity tracking, (2) gait reward, and (3)
regularization terms.

The reward function is summarized in Appendix Table IV. It is important to
note that the commands $\text{CMD}_{z,\gamma,\beta}$ (velocity mismatch term)
are intentionally set to zero. This is because we do not control them; rather,
we aim to maintain their values at zero to ensure stable and smooth walking.
In addition, the reward (contact pattern) encourages feet to align with their
contact masks, denoting swing, and stance phases, as illustrated in Appendix
Fig. 6. Therefore, the total reward at any time step $t$ is computed as the
weighted sum of individual reward components, expressed as
$r_{t}=\sum_{i}r_{i}\cdot\mu_{i}$, where $\mu_{i}$ represents the weighting
factor for each reward component $r_{i}$.

##  IV Experiments

In this section, we will illustrate the result of zero-shot transfer for both
sim-to-sim and sim-to-real scenarios. Additionally, we also provide
visualization of the calibration for MuJoCo to verify the effectiveness of
sim-to-sim. For the validation, we utilized Robot Era’s humanoid robots,
XBot-S and XBot-L, measuring 1.2 meters and 1.65 meters in height, shown in
Appendix Fig. 5, respectively.

###  IV-A Zero-shot Transfer

We carefully design domain randomization terms, as detailed in Appendix TABLE
III, to minimize the sim2real gap, following the approach outlined in [26].
Our agents are capable of transitioning to real-world environments via zero-
shot sim-to-real transfer, which is illustrated in Fig. 1. The standard
procedure involves training agents on GPUs, followed by policy analysis in
MuJoCo. For a comprehensive evaluation, we developed two types of terrains:
flat and uneven, as depicted in Appendix Fig. 7. The flat terrain replicates
the environment encountered during training in Isaac Gym, while the uneven
terrain offers a substantially more challenging landscape, differing
significantly from our initial training scenarios. Remarkably, our trained
policies enable the robots to traverse both types of terrain successfully.

###  IV-B Calibration for MuJoCo

We meticulously calibrated the MuJoCo environment to align its dynamics and
performance more closely with that of the real world. By comparing the leg
swing sine waves in both the MuJoCo and the real-world environments, which
represent joint positions over time that track a specified sine wave, we
observed nearly identical trajectories. These are depicted in Fig. 3.
Furthermore, we also compare the resulting phase portrait of the left knee
joint and left ankle pitch joint[17] within 5-second trajectories, 0.5 n/s
heading velocity, as shown in Fig. 4. It is clear to see that the dynamics in
MuJoCo are closer to the real environment than Isaac Gym.

Figure 3: Sine wave in Both MuJoCo and real-world environment. It can be found
that the trajectories of the two are very close after calibration. Figure 4:
Phase Portrait for MuJoCo, Real-World Environment, and Isaac Gym.

##  V CONCLUSIONS

Humanoid-Gym facilitates zero-shot transfer for humanoid robots of two
distinct sizes, from sim-to-sim and sim-to-real, via a specialized reward
function tailored for humanoid robotics. Our experimental outcomes reveal that
the adjusted MuJoCo simulation closely mirrors the dynamics and performance of
the real-world environment. This congruence enables researchers lacking
physical robots to validate training policies through sim-to-sim,
significantly enhancing the potential for successful sim-to-real transfers.

## ACKNOWLEDGMENT

The implementation of Humanoid-Gym relies on resources from legged_gym and
rsl_rl projects[1] created by the Robotic Systems Lab. We specifically utilize
the ‘LeggedRobot‘ implementation from their research to enhance our codebase.

## References

  * [1] N. Rudin, D. Hoeller, P. Reist, and M. Hutter, “Learning to walk in minutes using massively parallel deep reinforcement learning,” in _Conference on Robot Learning_. PMLR, 2022, pp. 91–100.
  * [2] A. Kumar, Z. Fu, D. Pathak, and J. Malik, “Rma: Rapid motor adaptation for legged robots,” _arXiv preprint arXiv:2107.04034_ , 2021.
  * [3] Y. Guo, Z. Jiang, Y.-J. Wang, J. Gao, and J. Chen, “Decentralized motor skill learning for complex robotic systems,” _IEEE Robotics and Automation Letters_ , 2023.
  * [4] X. B. Peng, M. Andrychowicz, W. Zaremba, and P. Abbeel, “Sim-to-real transfer of robotic control with dynamics randomization,” in _2018 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2018, pp. 3803–3810.
  * [5] W. Zhao, J. P. Queralta, and T. Westerlund, “Sim-to-real transfer in deep reinforcement learning for robotics: a survey,” in _2020 IEEE symposium series on computational intelligence (SSCI)_. IEEE, 2020, pp. 737–744.
  * [6] J. Tan, T. Zhang, E. Coumans, A. Iscen, Y. Bai, D. Hafner, S. Bohez, and V. Vanhoucke, “Sim-to-real: Learning agile locomotion for quadruped robots,” _arXiv preprint arXiv:1804.10332_ , 2018.
  * [7] A. Kadian, J. Truong, A. Gokaslan, A. Clegg, E. Wijmans, S. Lee, M. Savva, S. Chernova, and D. Batra, “Sim2real predictivity: Does evaluation in simulation predict real-world performance?” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 6670–6677, 2020.
  * [8] V. Makoviychuk, L. Wawrzyniak, Y. Guo, M. Lu, K. Storey, M. Macklin, D. Hoeller, N. Rudin, A. Allshire, A. Handa, _et al._ , “Isaac gym: High performance gpu-based physics simulation for robot learning,” _arXiv preprint arXiv:2108.10470_ , 2021.
  * [9] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in _2012 IEEE/RSJ international conference on intelligent robots and systems_. IEEE, 2012, pp. 5026–5033.
  * [10] RobotEra, “Robotera technology co.,ltd.” [Online]. Available: https://www.robotera.com/
  * [11] J. Hwangbo, J. Lee, A. Dosovitskiy, D. Bellicoso, V. Tsounis, V. Koltun, and M. Hutter, “Learning agile and dynamic motor skills for legged robots,” _Science Robotics_ , vol. 4, no. 26, p. eaau5872, 2019.
  * [12] J. Lee, J. Hwangbo, L. Wellhausen, V. Koltun, and M. Hutter, “Learning quadrupedal locomotion over challenging terrain,” _Science robotics_ , vol. 5, no. 47, p. eabc5986, 2020.
  * [13] G. B. Margolis, G. Yang, K. Paigwar, T. Chen, and P. Agrawal, “Rapid locomotion via reinforcement learning,” _The International Journal of Robotics Research_ , p. 02783649231224053, 2022.
  * [14] Z. Li, X. Cheng, X. B. Peng, P. Abbeel, S. Levine, G. Berseth, and K. Sreenath, “Reinforcement learning for robust parameterized locomotion control of bipedal robots,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 2811–2817.
  * [15] A. Kumar, Z. Li, J. Zeng, D. Pathak, K. Sreenath, and J. Malik, “Adapting rapid motor adaptation for bipedal robots,” in _2022 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2022, pp. 1161–1168.
  * [16] I. Radosavovic, T. Xiao, B. Zhang, T. Darrell, J. Malik, and K. Sreenath, “Learning humanoid locomotion with transformers,” _arXiv preprint arXiv:2303.03381_ , 2023.
  * [17] I. Radosavovic, B. Zhang, B. Shi, J. Rajasegaran, S. Kamat, T. Darrell, K. Sreenath, and J. Malik, “Humanoid locomotion as next token prediction,” _arXiv preprint arXiv:2402.19469_ , 2024.
  * [18] T. He, Z. Luo, W. Xiao, C. Zhang, K. Kitani, C. Liu, and G. Shi, “Learning human-to-humanoid real-time whole-body teleoperation,” _arXiv preprint arXiv:2403.04436_ , 2024.
  * [19] X. Cheng, Y. Ji, J. Chen, R. Yang, G. Yang, and X. Wang, “Expressive whole-body control for humanoid robots,” _arXiv preprint arXiv:2402.16796_ , 2024.
  * [20] M. T. Spaan, “Partially observable markov decision processes,” in _Reinforcement learning: State-of-the-art_. Springer, 2012, pp. 387–414.
  * [21] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal policy optimization algorithms,” _arXiv preprint arXiv:1707.06347_ , 2017.
  * [22] L. Pinto, M. Andrychowicz, P. Welinder, W. Zaremba, and P. Abbeel, “Asymmetric actor critic for image-based robot learning,” _arXiv preprint arXiv:1710.06542_ , 2017.
  * [23] J. Schulman, P. Moritz, S. Levine, M. Jordan, and P. Abbeel, “High-dimensional continuous control using generalized advantage estimation,” _arXiv preprint arXiv:1506.02438_ , 2015.
  * [24] J. Siekmann, Y. Godse, A. Fern, and J. Hurst, “Sim-to-real learning of all common bipedal gaits via periodic reward composition,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 7309–7315.
  * [25] Y. Yang, T. Zhang, E. Coumans, J. Tan, and B. Boots, “Fast and efficient locomotion via learned gait transitions,” in _Conference on robot learning_. PMLR, 2022, pp. 773–783.
  * [26] J. Tobin, R. Fong, A. Ray, J. Schneider, W. Zaremba, and P. Abbeel, “Domain randomization for transferring deep neural networks from simulation to the real world,” in _2017 IEEE/RSJ international conference on intelligent robots and systems (IROS)_. IEEE, 2017, pp. 23–30.

Figure 5: Hardware Platform. Our Humanoid-Gym framework is tested on two distinct sizes of humanoid robots, XBot-S and XBot-L, provided by Robot Era. Figure 6: The stance mask is the contact planning for the left (L) and right (R) feet, where 0 indicates the swing phase and 1 indicates the stance phase is expected. Figure 7: Terrains in MuJoCo. Humanoid-Gym provides two types of terrains utilized for sim-to-sim validation: flat planes and uneven terrains. Table II: Hyperparameters. Parameter | Value
---|---
Number of Environments | 8192
Number Training Epochs | 2
Batch size | $8192\times 24$
Episode Length | 2400 steps
Discount Factor | 0.994
GAE discount factor | 0.95
Entropy Regularization Coefficient | 0.001
$c1$ | 0.8
$c2$ | 1.2
Learning rate | 1e-5
Frame Stack of Single Observation | 15
Frame Stack of Single Privileged Observation | 3
Number of Single Observation | 47
Number of Single Privileged Observation | 73
Table III: Overview of Domain Randomization. Presented are the domain randomization terms and the associated parameter ranges. Additive randomization increments the parameter by a value within the specified range while scaling randomization adjusts it by a multiplicative factor from the same range. Parameter | Unit | Range | Operator | Type
---|---|---|---|---
Joint Position | rad | [-0.05, 0.05] | additive | Gaussian ($1\sigma$)
Joint Velocity | rad/s | [-0.5, 0.5] | additive | Gaussian ($1\sigma$)
Angular Velocity | rad/s | [-0.1, 0.1] | additive | Gaussian ($1\sigma$)
Euler Angle | rad | [-0.03, 0.03] | additive | Gaussian ($1\sigma$)
System Delay | ms | [0, 10] | - | Uniform
Friction | - | [0.1, 2.0] | - | Uniform
Motor Strength | % | [95, 105] | scaling | Gaussian ($1\sigma$)
Payload | kg | [-5, 5] | additive | Gaussian ($1\sigma$)
Table IV: In defining the reward function, we use a tracking error metric denoted by $\phi(e,w)$. This metric is expressed as $\phi(e,w):=\exp(-w\cdot\|e\|^{2}),$ where $e$ represents the tracking error, and $w$ is the associated weight. The target base height is set to $0.7\,\text{m}$.  Reward | Equation ($r_{i}$) | reward scale($\mu_{i}$)
---|---|---
Lin. velocity tracking | $\phi(\dot{P}^{b}_{xyz}-\text{CMD}_{xyz},5)$ | 1.2
Ang. velocity tracking | $\phi(\dot{P}^{b}_{\alpha\beta\gamma}-\text{CMD}_{\alpha\beta\gamma},5)$ | 1.0
Orientation tracking | $\phi(P^{b}_{\alpha\beta},5)$ | 1.0
Base height tracking | $\phi(P^{b}_{z}-0.7,100)$ | 0.5
Velocity mismatch | $\phi(\dot{P}^{b}_{z,\gamma,\beta}-\text{CMD}_{z,\gamma,\beta},5)$ | 0.5
Contact Pattern | $\phi(I_{p}(t)-I_{d}(t),\infty)$ | 1.0
Joint Position Tracking | $\phi(\theta-\theta_{\text{target}},2)$ | 1.5
Default Joint | $\phi(\theta_{t}-\theta_{0},2)$ | 0.2
Energy Cost | $|\tau||\dot{\theta}|$ | -0.0001
Action Smoothness | $\|a_{t}-2a_{t-1}+a_{t-2}\|_{2}$ | -0.01
Large contact | $\text{max}(F_{L,R}-400,0,100)$ | -0.01