# Bayesian Nonparametrics for Offline Skill Discovery

Valentin Villecroze  Harry J. Braviner  Panteha Naderian  Chris J. Maddison
Gabriel Loaiza-Ganem

###### Abstract

Skills or low-level policies in reinforcement learning are temporally extended
actions that can speed up learning and enable complex behaviours. Recent work
in offline reinforcement learning and imitation learning has proposed several
techniques for skill discovery from a set of expert trajectories. While these
methods are promising, the number $K$ of skills to discover is always a fixed
hyperparameter, which requires either prior knowledge about the environment or
an additional parameter search to tune it. We first propose a method for
offline learning of options (a particular skill framework) exploiting advances
in variational inference and continuous relaxations. We then highlight an
unexplored connection between Bayesian nonparametrics and offline skill
discovery, and show how to obtain a nonparametric version of our model. This
version is tractable thanks to a carefully structured approximate posterior
with a dynamically-changing number of options, removing the need to specify
$K$. We also show how our nonparametric extension can be applied in other
skill frameworks, and empirically demonstrate that our method can outperform
state-of-the-art offline skill learning algorithms across a variety of
environments. Our code is available at https://github.com/layer6ai-labs/BNPO.

Imitation learning, Bayesian nonparametrics, Reinforcement learning,
Variational inference

##  1 Introduction

Hierarchical policies have been studied in reinforcement learning for decades
(Dietterich et al., 1998; Vezhnevets et al., 2017). Low-level policies in
hierarchical frameworks can be interpreted as skills endowing agents with
temporal abstraction. These skills also offer a divide-and-conquer approach to
complex reinforcement learning environments, and learning them is therefore a
relevant problem.

One of the most commonly used skills frameworks is options (Sutton et al.,
1999). Methods for online learning of options have been proposed (Bacon et
al., 2017; Khetarpal et al., 2020), but cannot leverage data from expert
trajectories. To the best of our knowledge, the only method allowing offline
option discovery is that of Fox et al. (2017) (deep discovery of options,
DDO), which forgoes the use of highly successful variational inference
advances (Kingma & Welling, 2013; Rezende et al., 2014) for discrete latent
variables (Maddison et al., 2016; Jang et al., 2016) which in this case
correspond to the options to be inferred. Our first contribution is proposing
a method for offline learning of options combining these previously neglected
advances along with a judiciously constructed approximate posterior, which we
show empirically outperforms not only DDO, but also other offline skill
discovery algorithms.

Additionally, all (discrete) skill learning approaches we are aware of –
including options – require specifying the number of skills $K$ as a
hyperparameter, rather than learning it. This is a significant practical
limitation, as $K$ cannot realistically be known in advance for complex
environments. Our second contributions is highlighting a similarity between
mixture models and inferring skills from expert trajectories: both are a form
of soft clustering where observed data (state-action pairs) are
probabilistically assigned to their corresponding cluster (skill). Bayesian
nonparametric approaches such as Dirichlet process mixture models (Neal, 2000)
are not only mathematically elegant, but also circumvent the practical need to
prespecify the number of clusters. We therefore argue that Bayesian
nonparametrics should be more heavily used for option and skill discovery.

Our third contribution is proposing a scheme allowing for our option discovery
method to accommodate a nonparametric prior over options (thus modelling $K$
as infinite), which can also be applied in other skill learning frameworks. By
further adding structure to the variational posterior, allowing it to
dynamically change the number of options/skills being assigned positive
probability, we recover a method that is not only tractable, but also retains
the nonparametric aspect of the model and does not require treating $K$ as a
hyperparameter. We show that our nonparametric versions of skill learning
algorithms match the performance of their parametric counterparts with tuned
$K$, thus successfully “learning a sufficient number of options” in practice.
Finally, we hope that our nonparametric variational inference scheme will find
uses outside of offline skill recovery.

##  2 Related Work

Imitation learning (Hussein et al., 2017), also called learning from
demonstrations (Schaal et al., 1997; Argall et al., 2009), consists of
learning to solve a task from a set of expert demonstrations. This can be
achieved by methods such as behavioural cloning (Esmaili et al., 1995; Ross et
al., 2011) or inverse reinforcement learning (Ho & Ermon, 2016; Sharma et al.,
2018; Krishnan et al., 2019). In this paper we focus on the former. Many
recent works on imitation learning attempt to divide the demonstrations into
several smaller segments, before fitting models on each of them (Niekum et
al., 2013; Murali et al., 2016; Krishnan et al., 2018; Shiarlis et al., 2018).
These works usually consider these two steps as two distinct stages, and
unlike our proposed approach, do not learn skills end-to-end. Niekum et al.
(2013) and Krishnan et al. (2018) are of particular interest as they are the
only ones, to our knowledge, to make use of a Bayesian nonparametric model to
segment the trajectories. However, they significantly differ from our method
in that they use handcrafted movement primitives and linear dynamical models
to fit the resulting segments.

Our work is closer to the field of option and skill discovery, which leverages
work on online hierarchical reinforcement learning (Sutton et al., 1999;
Gregor et al., 2016; Bacon et al., 2017; Achiam et al., 2018; Eysenbach et
al., 2018; Sharma et al., 2019; Florensa et al., 2017) and aims to learn
adequate options from a set of expert trajectories. DDO (Fox et al., 2017)
uses an EM-type algorithm (Dempster et al., 1977) that allows multi-level
hierarchies. As previously mentioned, this method forgoes advances in
variational inference and continuous relaxations, and we will later show that
option learning can be significantly improved upon through the use of these
techniques. While not following the option framework, Kipf et al. (2019)
propose CompILE, which uses a variational approach and continuous relaxations
to both segment the trajectories and encode the resulting slices as discrete
skills. We will show that with these advances, options also outperform
CompILE.

Both DDO and CompILE need the number of options/skills, $K$, to be specified
beforehand. We borrow from the Bayesian nonparametrics literature (Teh et al.,
2006; Caron et al., 2007; Dunson & Xing, 2009) – in particular from models
involving Dirichlet processes (Ferguson, 1973) – in order to model options. We
follow Nalisnick & Smyth (2016) and use a variational-autoencoder-inspired
(Kingma & Welling, 2013; Rezende et al., 2014) inference scheme that avoids
the need for the expensive MCMC samplers commonly associated with these types
of models (Neal, 2000). As a result, the nonparametric versions of the models
we consider do not need to prescpecify nor tune $K$.

Finally, Shankar & Gupta (2020) and Ajay et al. (2021) also use variational
inference to discover skills but choose to represent them as continuous latent
variables instead of categorical variables. While this choice also technically
leads to infinitely many skills, as we will see in our experiments, discrete
skills learned from expert trajectories are a useful way to enhance online
agents. The use of continuous skills results not only in less interpretable
skills, but also loses the ability to perform these enhancements.

##  3 Preliminaries

We now review the main concepts needed for our model: behavioural cloning, the
options framework, continuous relaxations, and the stick-breaking process. We
also review CompILE, a method that we shall later make nonparametric.

###  3.1 Behavioural Cloning

The goal of behavioural cloning is to imitate an expert who demonstrates how
to accomplish a task. More formally, we consider a Markov Decision Process
without rewards (MDP\R). An MDP\R is a tuple
$\mathcal{M}:\langle\mathcal{S},\mathcal{A},P,\rho\rangle$, where
$\mathcal{S}$ is the state space, $\mathcal{A}$ the action space,
$P\left(s_{t+1}\mid s_{t},a_{t}\right)$ a transition distribution, and $\rho$
the starting state distribution. We assume access to a dataset of expert
trajectories
$\xi:=\\{\xi^{(i)}:=(s_{0}^{(i)},a_{0}^{(i)},s_{1}^{(i)},a_{1}^{(i)},\dots,s_{T}^{(i)})\\}^{N}_{i=1}$
of states $s_{t}^{(i)}\in\mathcal{S}$ and actions $a_{t}^{(i)}\in\mathcal{A}$.
The trajectory length $T$ need not be identical across trajectories and could
be changed to $T_{i}$, but we keep $T$ to avoid further complicating notation.
We want to find a policy $\pi_{\theta}:\mathcal{S}\mapsto\Delta(\mathcal{A})$
(where $\Delta(\mathcal{A})$ denotes the set of distributions over
$\mathcal{A}$) parameterized by $\theta$ that maximizes the logarithm of
following expression:

$\displaystyle p_{\theta}(\xi)$
$\displaystyle\propto\prod_{i=1}^{N}\prod_{t=0}^{T-1}\pi_{\theta}(a_{t}^{(i)}|s_{t}^{(i)}),$
(1)

i.e. the policy maximizing the likelihood assigned to the expert trajectories.
We have dropped terms that are constant with respect to $\theta$ in Equation
1.

###  3.2 The Options Framework

As previously mentioned, the option framework introduces temporally extended
actions, allowing for temporal abstraction and more complex behaviour.
Formally, an option $h\in\Omega$ (Sutton et al., 1999) is a tuple
$\langle\mathcal{I}_{h},\pi_{h},\psi_{h}\rangle$, where
$\mathcal{I}_{h}\subseteq\mathcal{S}$ is the initiation set, $\pi_{h}(a|s)$
the control or low-level policy, and $\psi_{h}:\mathcal{S}\rightarrow[0,1]$
the termination function. An option $h$ can be invoked in any state
$s\in\mathcal{I}_{h}$, in which case actions are drawn according to $\pi_{h}$
until the option terminates, which happens with probability
$\psi_{h}(s^{\prime})$ at each subsequent state $s^{\prime}$. In our setting
we assume that $\mathcal{I}_{h}=\mathcal{S}$, which is a common assumption
(Bacon et al., 2017; Fox et al., 2017; Shankar & Gupta, 2020). In order to use
these options to solve a task, a high-level policy (or policy over options)
$\eta:\mathcal{S}\mapsto\Delta(\Omega)$ is used to select a new option $h$
after the previous one terminates. This two-level structure is called a
hierarchical policy. The resulting generative process is described in
Algorithm 1.

Learning a hierarchical policy from expert demonstrations using behavioural
cloning is much more challenging than with a single policy. Indeed, not only
do multiple policies need to be learned concurrently (one for each option),
but, as the options are unobserved, they must be treated as latent variables.
Fox et al. (2017) propose an EM-based approach to do this, which we improve
upon.

Algorithm 1 Trajectory generation with options.

1: $s_{0}\sim\rho(\cdot)$, $b_{0}\leftarrow 1$

2: for $t\in[0,\dots,T-1]$ do

3: if $b_{t}=1$ then

4: $h_{t}\sim\eta(\cdot|s_{t})$ \\\Draw option from high-level policy

5: else

6: $h_{t}\leftarrow h_{t-1}$

7: end if

8: $a_{t}\sim\pi_{h_{t}}(\cdot|s_{t})$ \\\Draw action from low-level policy

9: $s_{t+1}\sim P(\cdot|s_{t},a_{t})$ \\\Draw the next state

10: $b_{t+1}\sim\operatorname{Bernoulli}(\cdot|\psi_{h_{t}}(s_{t+1}))$
\\\Terminate option?

11: end for

###  3.3 CompILE

CompILE (Kipf et al., 2019) is used for hierarchical imitation learning but
does not rely on the option framework. Binary termination variables are not
sampled at every timestep to determine if the current skill should terminate.
Rather, whenever a new skill is drawn from the high-level policy, CompILE also
samples the number of steps for which the skill will be active from a Poisson
distribution with parameter $\lambda$. Also, rather than using a state-
dependent high-level policy, $\eta$ is assumed to be uniform over skills. We
refer to any time interval between subsequent skill selections as a segment.
This procedure is summarized in Algorithm 2, where we have abused notation and
still use $b$ to denote termination variables in order to highlight the
similarity with the options framework, even though these are not the same
variables as in options. Similarly, $h_{j}$ does not denote an option, but
rather the skill being used during segment $j$.

Doing behavioural cloning in ComPILE is harder than just maximizing Equation
1, since computing $\log p_{\theta}(\xi)$, where $\theta$ now parameterizes
all the low-level policies, requires an intractable marginalization over the
unobserved variables ($b$’s and $h$’s). In order to circumvent this issue, an
approximate posterior $q_{\phi}(\zeta|\xi)$ is introduced, where $\zeta$
denotes all the unobserved variables for all trajectories and $\phi$ the
variational parameters. Instead of directly maximizing $\log p_{\theta}(\xi)$,
the following lower bound, which is called the ELBO, is maximized over
$(\theta,\phi)$:

$\displaystyle L(\theta,\phi)$
$\displaystyle:=\mathbb{E}_{q_{\phi}(\zeta|\xi)}[\log
p_{\theta}(\zeta,\xi)-\log q_{\phi}(\zeta|\xi)]$ (2) $\displaystyle\leq\log
p_{\theta}(\xi).$ (3)

We omit the details on how $q_{\phi}(\zeta|\xi)$ is structured, but highlight
that the number of segments needs to be specified as a hyperparameter of this
variational approximation rather than being properly treated as random. In
order to maximize Equation 2, the authors use the reparameterization trick of
variational autoencoders (Kingma & Welling, 2013; Rezende et al., 2014) along
with continuous relaxations. In particular, they use the Concrete or Gumbel-
Softmax (GS) distribution (Maddison et al., 2016; Jang et al., 2016), which we
briefly review in the next section, in order to efficiently backpropagate
through Equation 2.

Algorithm 2 Trajectory generation with CompILE.

1: $s_{0}\sim\rho(\cdot)$, $b_{0}\leftarrow 0$, $j\leftarrow 0$

2: for $t\in[0,\dots,T-1]$ do

3: while $t=b_{j}$ do

4: $j\leftarrow j+1$ \\\Consider the next segment

5: $h_{j}\sim\eta(\cdot)=\operatorname{Uniform}(\cdot|K)$ \\\Draw skill

6: $b_{j}\sim\operatorname{Poisson}(\cdot|\lambda)$ \\\Draw segment length

7: $b_{j}\leftarrow b_{j}+b_{j-1}$ \\\Next segment boundary

8: end while

9: $a_{t}\sim\pi_{h_{j}}(\cdot|s_{t})$ \\\Draw action from low-level policy

10: $s_{t+1}\sim P(\cdot|s_{t},a_{t})$ \\\Draw the next state

11: end for

###  3.4 Continuous Relaxations

The reparameterization trick is used in variational autoencoders to
backpropagate through expressions of the form
$\mathbb{E}_{q_{\phi}(\zeta)}[f_{\phi}(\zeta)]$, with respect to $\phi$, where
$f_{\phi}$ is a real-valued function. Since the distribution with respect to
which the expectation is taken depends on $\phi$, one cannot simply bring the
gradient inside of the expectation. Gradient estimators such as REINFORCE
(Glynn, 1990; Williams, 1992) typically exhibit high variance, which the
reparameterization trick empirically reduces. When $\zeta$ is a continuous
random variable, it is often (but not always) the case that one can easily
find a continuously differentiable function $g$ such that $\zeta\sim
q_{\phi}(\cdot)\iff\zeta=g(\epsilon,\phi)$ where $\epsilon$ follows some
continuous distribution which does not depend on $\phi$. In this case, the
gradient is given by:

$\nabla_{\phi}\mathbb{E}_{q_{\phi}(\zeta)}[f_{\phi}(\zeta)]=\mathbb{E}_{\epsilon}[\nabla_{\phi}f_{\phi}(g(\epsilon,\phi))],$
(4)

and a Monte Carlo estimate can be easily obtained.

When $\zeta$ is categorical, $g$ has to be piece-wise constant and Equation 4
no longer holds. Continuous relaxations approximate
$\mathbb{E}_{q_{\phi}(\zeta)}[f_{\phi}(\zeta)]$ with an expectation over a
continuous random variable, so that the reparameterization trick can be used.
The Gumbel-Softmax (GS) or Concrete distribution is a distribution on the
$K$-simplex,
$\Delta(K):=\\{x\in\mathbb{R}^{K}:x_{k}>0,\sum_{k=0}^{K-1}x_{k}=1\\}$,
parameterized by $q\in\Delta(K)$ and a temperature hyperparameter $\tau>0$,
designed to address this issue. It is reparameterized as follows:

$\zeta\sim
GS_{\tau}(\cdot|q)\iff\zeta=\operatorname{softmax}\left(\dfrac{\epsilon+\log
q}{\tau}\right),$ (5)

where $\epsilon$ is a $K$-dimensional vector with independent Gumbel$(0,1)$
entries, and the $\log$ is taken elementwise. As $\tau\rightarrow 0$, the
$GS_{\tau}(\cdot|q)$ distribution converges to the discrete distribution $q$.
By thinking of categorical $\zeta$’s as one-hot vectors of length $K$, the GS
thus provides a continuous relaxation of $\zeta$, and
$\mathbb{E}_{q_{\phi}(\zeta)}[f_{\phi}(\zeta)]$ can be approximated by:

$\mathbb{E}_{q_{\phi}(\zeta)}[f_{\phi}(\zeta)]\approx\mathbb{E}_{GS_{\tau}(\zeta|q_{\phi}(\cdot))}[\tilde{f}(\zeta)],$
(6)

where we think of the discrete distribution $q_{\phi}(\cdot)$ as a
$K$-dimensional vector, and $\tilde{f}$ is a relaxation of $f$ mapping all of
$\Delta(K)$ (rather than just the vertices, i.e. one-hot vectors) to
$\mathbb{R}$. Equation 6 admits the gradient estimator of Equation 4, and Kipf
et al. (2019) use it to optimize Equation 2.

###  3.5 The Stick-Breaking Process

The stick-breaking process (Sethuraman, 1994; Ishwaran & James, 2001), which
is deeply connected to Dirichlet processes, places a distribution over
probability vectors of infinite length. Equivalently, it is a distribution
over distributions on a countably infinite set. This process will allow us to
both assume that there are infinitely many options, and to place a proper
prior on the high-level policy $\eta$ itself (i.e. a distribution over
distributions of options). When performing posterior inference some options
will have extremely small probabilities; this enables us to “learn the number
of options” from the data. More formally, the Griffiths-Engen-McCloskey
distribution, denoted $\operatorname{GEM}(\alpha)$ and parameterized by
$\alpha>0$, produces an infinitely long vector $(\beta_{0},\beta_{1},\dots)$
such that $\beta_{k}>0$ and $\sum_{k=0}^{\infty}\beta_{k}=1$ almost surely.
Samples are obtained by first sampling
$\beta^{\prime}_{k}\sim\operatorname{Beta}(1,\alpha)$ for $k=0,1,\dots$ and
then setting:

$\beta_{k}=\beta_{k}^{\prime}\prod_{l=0}^{k-1}(1-\beta_{l}^{\prime}).$ (7)

Intuitively, we start off with a stick of length $1$, and at every step we
remove $\beta_{k}$ from the stick. $\beta^{\prime}_{k}$ is the fraction of the
remaining stick that is assigned to $\beta_{k}$.

We shall later use the $\operatorname{GEM}$ distribution as a prior, and will
want to perform variational inference. Applying the stick-breaking procedure
to Beta random variables with different parameters might seem like the most
natural approximate posterior. However, the Beta distribution is not easily
reparameterized, and thus inference with such a posterior becomes
computationally challenging. It is therefore common to instead use the
Kumaraswamy distribution in the approximate posterior (Nalisnick & Smyth,
2016; Stirn et al., 2019). Like the Beta distribution, this is a
$(0,1)$-valued distribution with two parameters $a_{1},a_{2}>0$, whose density
is given by:

$p(x|a_{1},a_{2})=a_{1}a_{2}x^{a_{1}-1}(1-x^{a_{1}})^{a_{2}-1},$ (8)

and which can easily be reparameterized.

##  4 Our Model

We first describe a parametric version of our model, where the number of
options $K$ is a fixed hyperparameter. We will detail in Sections 5.1 and 5.2
how we use a nonparametric prior to circumvent the need to specify this
number. Throughout this section and until Section 5.3 (exclusive), our
notation refers to options and not to CompILE.

###  4.1 Overview

Now that we have covered the required concepts, we introduce our method in
more detail. We assume that the expert trajectories $\xi$ are generated by a
two-level option hierarchy, as described in Algorithm 1, with the high-level
policy $\eta$ being shared across trajectories. We denote the corresponding
trajectories of hidden options $h$ and binary termination variables $b$ as
$\zeta:=\\{\zeta^{(i)}:=(b_{0}^{(i)},h_{0}^{(i)},b_{1}^{(i)},h_{1}^{(i)},\dots,h_{T-1}^{(i)})\\}^{N}_{i=1}$,
and $\eta$ the high-level policy. Fox et al. (2017) and Kipf et al. (2019)
found that taking $\eta$ as a uniform policy rather than learning it results
in equally useful learned skills. We simplify $\eta$ with the less restrictive
assumption that it does not depend on the current state $s$, a choice that
will later simplify our nonparametric component. That is, in Algorithm 1, we
draw $h_{t}$ according to $\eta(\cdot)$ rather than $\eta(\cdot|s_{t})$. We
then treat $\eta$ in a Bayesian way, as a latent variable to be inferred from
observed trajectories, and assume a $K$-dimensional prior over it, obtained by
truncating the stick-breaking process of Beta$(1,\alpha)$ variables after
$K-1$ steps (and having the last entry be such that the resulting vector adds
up to $1$). The resulting graphical model can be seen in Figure 1.

$s_{0}^{(i)}$$a_{0}^{(i)}$$b_{0}^{(i)}$$h_{0}^{(i)}$$s_{1}^{(i)}$$a_{1}^{(i)}$$b_{1}^{(i)}$$h_{1}^{(i)}$$\ldots$$\ldots$$\ldots$$\ldots$$s_{T-1}^{(i)}$$a_{T-1}^{(i)}$$b_{T-1}^{(i)}$$h_{T-1}^{(i)}$$s_{T}^{(i)}$$i=1,\dots,N$$\eta$

Figure 1: Graphical model for options. Shaded nodes correspond to observed
variables ($\xi^{(i)}$) while the blank ones are latent ($\zeta^{(i)}$).

We denote as $\theta$ all the parameters from the prior (i.e. $\alpha$), the
low-level policies, and the terminations functions. We implement the
termination functions with a single neural network taking $s$ as input and
outputting the termination probabilities for every option. We keep a separate
neural network for every low-level policy, each of which takes $s$ as input
and outputs the parameters of a categorical over actions (we assume a discrete
action space, but our method can easily be modified for continuous action
spaces). Note that, other than the Bayesian treatment of $\eta$ and the
inclusion of $\alpha$ in $\theta$, our model is identical to DDO. We will
detail differences in how learning is performed in this section, and later
show how our model can be made nonparametric.

As with CompILE, naïvely doing behavioural cloning is intractable, and we
therefore introduce an approximate posterior $q_{\phi}(\eta,\zeta|\xi)$. As we
will cover in detail later, we endow this posterior with a structure
respecting the conditional independences implied by Figure 1, and unlike
CompILE, we do not have to arbitrarily specify a number of segments.

###  4.2 Objective

We jointly train $\theta$ and $\phi$ by using the ELBO as the maximization
objective, which is now given by:

$\displaystyle L(\theta,\phi):=$
$\displaystyle\;\mathbb{E}_{q_{\phi}(\eta,\zeta|\xi)}[\log
p_{\theta}(\zeta,\xi|\eta)-\log q_{\phi}(\zeta|\eta,\xi)]$
$\displaystyle-\mathbb{KL}(q_{\phi}(\eta)||p(\eta))\leq\log p_{\theta}(\xi).$
(9)

We use the GS distribution in order to perform the reparameterization trick
and backpropagate through the ELBO. If we write the first term in more detail,
we get:

$\displaystyle\log p_{\theta}(\zeta,\xi|\eta)=\displaystyle\sum_{i=1}^{N}\log
p_{\theta}(\zeta^{(i)},\xi^{(i)}|\eta),$ (10)

where each term is given by:

$\displaystyle\medmath{\log p_{\theta}(\zeta^{(i)},\xi^{(i)}|\eta)}$
$\displaystyle=\medmath{\log\rho(s_{0}^{(i)})+\log\delta_{b_{0}^{(i)}=1}}$
$\displaystyle+\displaystyle\medmath{\log\eta(h_{0}^{(i)})+\sum_{t=1}^{T-1}\log
p_{\theta}(b_{t}^{(i)},h_{t}^{(i)}|h_{t-1}^{(i)},s_{t}^{(i)},\eta)}$ (11)
$\displaystyle+\medmath{\sum_{t=0}^{T-1}\log\pi_{h_{t}^{(i)}}(a_{t}^{(i)}|s_{t}^{(i)})+\log
P(s_{t+1}^{(i)}|s_{t}^{(i)},a_{t}^{(i)})}.$

Note that the initial state distribution and environment dynamics terms cannot
be evaluated but, as they do not depend on the parameters being optimized,
they can be ignored. We can further decompose the term for $b_{t}^{(i)}$ and
$h_{t}^{(i)}$:

$\displaystyle
p_{\theta}(b_{t}^{(i)}=1,h_{t}^{(i)}|h_{t-1}^{(i)},s_{t}^{(i)},\eta)=$
$\displaystyle\;\psi_{h_{t-1}^{(i)}}(s_{t}^{(i)})\eta(h_{t}^{(i)})$ (12)
$\displaystyle
p_{\theta}(b_{t}^{(i)}=0,h_{t}^{(i)}|h_{t-1}^{(i)},s_{t}^{(i)},\eta)=$
$\displaystyle\;(1-\psi_{h_{t-1}^{(i)}}(s_{t}^{(i)}))$
$\displaystyle\cdot\delta_{h_{t}^{(i)}=h_{t-1}^{(i)}},$ (13)

where $\\{\psi_{h}\\}_{h}$ are the termination functions. Note that equations
4.2, 12 and 13 assume that the $b$’s are binary and the $h$’s are categorical
(one-hot) through the Kronecker delta terms and evaluating $\eta$ at different
$h$ terms. Since we use a Gumbel-Softmax relaxation to optimize the ELBO, the
sampled values of $b$’s and $h$’s will be $(0,1)$ and simplex-valued,
respectively. We therefore relax the corresponding terms so that they can be
evaluated for such values and gradients can be propagated. We detail the
relaxations in Appendix A.

###  4.3 Variational Posterior

We now describe the structure of our approximate posterior
$q_{\phi}(\eta,\zeta|\xi)$. First, we observe from Figure 1 through the rules
of $d$-separation (Koller & Friedman, 2009) that $\zeta^{(i)}$ is independent
of $\zeta^{(j)}$ given $\eta$ and $\xi$, provided $i\neq j$. Additionally,
given $\eta$, $\zeta^{(i)}$ depends on $\xi$ only through $\xi^{(i)}$. We thus
take $q_{\phi}(\eta,\zeta|\xi)$ to obey the following conditional independence
relationship:

$\displaystyle
q_{\phi}(\eta,\zeta|\xi)=q_{\phi}(\eta)\displaystyle\prod_{i=1}^{N}q_{\phi}(\zeta^{(i)}|\eta,\xi^{(i)}),$
(14)

which holds for the true posterior as well.

Since $\eta$ is a global variable (i.e. it does not have components for each
trajectory $i$) while $\zeta=\\{\zeta^{(i)}\\}_{i}$ is composed of local
variables assigned to the corresponding $\xi^{(i)}$, we treat $\eta$ in a non-
amortized (Gershman & Goodman, 2014) way (i.e. the parameters of
$q_{\phi}(\eta)$ are directly optimized, there is no neural network), while we
treat $\zeta^{(i)}$ in an amortized way (i.e. a neural network takes $\eta$
and $\xi^{(i)}$ as inputs to produce the parameters of a distribution).

As previously hinted, we parameterize $q_{\phi}(\eta)$ as a sequence of
Kumaraswamy distributions to which the stick-breaking procedure is applied,
which allows us to straightforwardly use the reparameterization trick.

As mentioned above, $q_{\phi}(\zeta^{(i)}|\eta,\xi^{(i)})$ is treated in an
amortized way, so that a neural network takes $\eta$ and $\xi^{(i)}$ as input
and produces the parameters for the distribution of $\zeta^{(i)}$, i.e. the
Bernoulli parameters of each $b_{t}^{(i)}$ and categorical parameters of
$h_{t}^{(i)}$ for $t=0,\dots,T-1$. We use an autoregressive structure for
$q_{\phi}(\zeta^{(i)}|\eta,\xi^{(i)})$, and assume conditional independence of
options and terminations at every time step:

$\displaystyle\medmath{q_{\phi}(\zeta^{(i)}|\eta,\xi^{(i)})}=$
$\displaystyle\;\medmath{q_{\phi}(b_{0}^{(i)}|\eta,\xi^{(i)})q_{\phi}(h_{0}^{(i)}|\eta,\xi^{(i)})}$
$\displaystyle\medmath{\cdot\prod_{t=1}^{T-1}q_{\phi}(b_{t}^{(i)}|\zeta_{t-1}^{(i)},\eta,\xi^{(i)}_{t:T-1})q_{\phi}(h_{t}^{(i)}|\zeta_{t-1}^{(i)},\eta,\xi^{(i)}_{t:T-1})},$
(15)

where $\xi_{t}^{(i)}=(s_{t}^{(i)},a_{t}^{(i)})$ and
$\zeta_{t}^{(i)}=(b_{t}^{(i)},h_{t}^{(i)})$. Note that conditioned on
$(\zeta_{t-1}^{(i)},\eta,\xi^{(i)})$, $(b_{t}^{(i)},h_{t}^{(i)})$ is
independent of $\xi_{0:t-1}$, which can again be verified through Figure 1. In
other words, states and actions after time $t$ can be useful to infer the
option and termination at time $t$; but previous states and actions are not
(assuming the previous option and termination are known). This is why we only
condition on $\xi^{(i)}_{t:T-1}$ in Equation 4.3, rather than on $\xi^{(i)}$.
In practice, we use an LSTM (Hochreiter & Schmidhuber, 1997) to parse
$\xi^{(i)}$ in reverse order, then sequentially sample $\zeta^{(i)}$ using
multi-layer perceptrons (MLPs) to output the distributions’ parameters as
shown below (the MLPs’ inputs are concatenated):

$\displaystyle\medmath{b_{t}^{(i)}|\zeta_{t-1}^{(i)},\eta,\xi^{(i)}}$
$\displaystyle\medmath{\sim\text{\scriptsize{GS}}_{\tau}(\cdot|\text{\scriptsize{MLP}}_{b}(\text{\scriptsize{LSTM}}(\xi^{(i)}_{t:T-1}),\eta,\zeta_{t-1}^{(i)}))},$
(16) $\displaystyle\medmath{h_{t}^{(i)}|\zeta_{t-1}^{(i)},\eta,\xi^{(i)}}$
$\displaystyle\medmath{\sim\text{\scriptsize{GS}}_{\tau}(\cdot|\text{\scriptsize{MLP}}_{h}(\text{\scriptsize{LSTM}}(\xi^{(i)}_{t:T-1}),\eta,\zeta_{t-1}^{(i)}))}.$
(17)

The $b$ and $h$ indices reflect the fact that we use two separate output
heads. The MLPs share their parameters except those of their last linear
layer. The $\text{\small{LSTM}}(\xi^{(i)}_{t:T-1})$ term denotes the hidden
state of the LSTM layer taken at time step $t$. Also note that we abuse
notation and use $b_{t}^{(i)}$, $h_{t}^{(i)}$ and $\zeta_{t}^{(i)}$ for both
the discrete variables and their relaxed counterparts. We refer to the LSTM
and the MLP heads as the encoder.

###  4.4 Entropy Regularizer

In some preliminary experiments, we observed that our model could get stuck in
some local optima where too few options were used. To address this issue, we
add a regularizing term to the ELBO in Equation 4.2, defined as the entropy
over the average of the sampled options $\\{h_{t}^{(i)}\\}$. More formally, we
define $h^{(i)}_{avg}:=\frac{1}{T}\sum_{t=0}^{T-1}h_{t}^{(i)}$, which is a
probability vector of size $K$ and consider the following regularizing term
that we want to maximize:

$l_{ent}:=-\mathbb{E}_{q_{\phi}(\zeta|\xi)}\left[\sum_{i=1}^{N}\sum_{k=0}^{K-1}h^{(i)}_{avg,k}\log\left(h^{(i)}_{avg,k}\right)\right],$
(18)

where $h^{(i)}_{avg,k}$ is the $k$-th coordinate of $h^{(i)}_{avg}$. We also
use the reparameterization trick to backpropagate through Equation 18. This
term encourages the model to use all the available options in equal amounts.
Our final maximization objective is given by the following expression:

$L(\theta,\phi)+\lambda_{ent}l_{ent},$ (19)

where we anneal $\lambda_{ent}\geq 0$ throughout training with a fixed decay
rate. We highlight that we only anneal $\lambda_{ent}$ since the ELBO is a
highly principled objective, and annealing causes the objective being
optimized to be closer and closer to the ELBO. That being said, we did not
observe significant empirical differences between annealing and fixing
$\lambda_{ent}$ throughout training (see Appendix D).

##  5 Nonparametric Models

We now describe how to use a nonparametric prior to remove the need for a
prespecified number of options.

###  5.1 GEM Prior

We now put a $\operatorname{GEM}(\alpha)$ prior over $\eta$, replacing the
prior obtained from truncating the stick-breaking process at $K$ steps. The
result is a nonparametric model which assumes a countably infinite number of
options. Doing posterior inference over $\eta$ then enables learning the
number of options present in the observed trajectories.

Note that $\theta$ should have an infinite number of parameters, since we
consider an infinite number of options, but we detail below how we manage this
in practice.

###  5.2 Truncation

Recall that $\eta$ is now an infinite dimensional vector whose entries add up
to one. A naïve attempt to optimize the objective from Equation 4.2 would thus
involve not only sampling infinitely long $\eta$’s, but also neural networks
with infinitely many parameters (for terminations), or infinitely many neural
networks (for low-level policies); this is clearly infeasible. Nalisnick &
Smyth (2016) truncate $\eta$ and consider only the $K$ first coordinates,
treating $K$ as a hyperparameter. We highlight that fixing $K$ in the
approximate posterior yields the same objective as reverting the nonparametric
$\operatorname{GEM}$ prior back to a $K$-dimensional prior; thus discarding
the nonparametric aspect of the model. Even in this setting, our fixed-$K$
model remains a novel way of learning options offline. Nonetheless, in order
to properly retain the nonparametric aspect of the model, we allow $K$ to
increase throughout training. We check at regular intervals (every $n_{K}$
training steps) whether we should increment $K$ by looking at the “usage”,
$U(h)$, for each option $h$:

$\displaystyle\medmath{U(h):=\frac{1}{NT}\sum_{i=1}^{N}\sum_{t=0}^{T-1}\mathds{1}\left(h=\underset{h_{k}:0\leq
k<K}{\operatorname{argmax}}\;\pi_{h_{k}}(a_{t}^{(i)}|s_{t}^{(i)})\right)},$
(20)

i.e. $U(h)$ is the fraction of steps for which we inferred $h$ to be the
option most likely to have been in use. Our rule for adding skills (clusters)
is inspired by the small variance asymptotics Bayesian nonparametrics
literature (Kulis & Jordan, 2011; Broderick et al., 2013). We increase $K$ by
$1$ if there is no option $h_{k}$ for $k\in\\{0,...,K-1\\}$ such that
$U(h_{k})<\delta/K$ for a fixed hyperparameter $\delta<1$. When increasing
$K$, we add a new low-level policy, and expand the termination function by
adding a new row to the last linear layer, so that the output changes from
size $K$ to size $K+1$. We handle this new option in the encoder’s MLPs (see
Equations 16 and 17) by adding two new columns to their first layer (to take
into account the new input size, as both $\eta$ and $\zeta_{t-1}^{(i)}$ have
increased in length by 1) and a new row to the last layer of $\text{MLP}_{h}$
(to output the correct number of parameters for $\zeta_{t}^{(i)}$). We also
add parameters for another Kumaraswamy distribution to $\phi$. All new
parameters are initialized randomly. The logic behind this heuristic is that
there is no need to add a new option (i.e. increment $K$) if there is an
existing option which is rarely used.

###  5.3 Nonparametric CompILE

Finally, we show that CompILE can now be easily turned into a nonparametric
model thanks to the machinery we have developed. When placing a
$\operatorname{GEM}(\alpha)$ prior over $\eta$, the CompILE ELBO from Equation
2 changes to accommodate the fact that $\eta$ is now being treated in a
Bayesian way rather than being fixed and uniform. The ELBO then becomes:

$\displaystyle L(\theta,\phi):=$
$\displaystyle\;\mathbb{E}_{q_{\phi}(\eta,\zeta|\xi)}[\log
p_{\theta}(\zeta,\xi|\eta)-\log q_{\phi}(\zeta|\eta,\xi)]$
$\displaystyle-\mathbb{KL}(q_{\phi}(\eta)||p(\eta)).$ (21)

While this objective looks identical to ours, we highlight once again that we
have switched back to CompILE notation, and that the terms are not identical
to those of Equation 4.2. The only practical changes to CompILE are then:
$(1)$ that we also use a non-amortized posterior $q_{\phi}(\eta)$ obtained by
stick-breaking Kumaraswamy random variables, and so the $\mathbb{KL}$ term in
the ELBO is identical to that of Equation 4.2; $(2)$ the approximate posterior
becomes $q_{\phi}(\eta,\zeta|\xi)=q_{\phi}(\eta)q_{\phi}(\zeta|\eta,\xi)$,
which implies that, as with our own model, the neural network modeling
$q_{\phi}(\zeta|\eta,\xi)$ now needs to take $\eta$ as an extra input,
although we leave the rest of the structure used in CompILE intact; and $(3)$
when evaluating $\eta(h_{j}^{(i)})$ for $p_{\theta}(\zeta,\xi|\eta)$, where
$h_{j}^{(i)}$ is the skill used for the $j$-th segment in the $i$-th
trajectory, we actually evaluate $\eta$ (as we do in options), rather than
using $1/K$ as implied by the uniform distribution. Additionally, for the
nonparametric version of CompILE we use the same truncation and entropy
regularizer as we used for our options model.

##  6 Experiments

Figure 2: Results for our proof-of-concept environment. We consider 10 random
seeds per vocabulary size $n_{V}$, and a run is considered a success if one of
the learned options selects the correct action at $t=4$ with probability at
least $0.95$. The red line shows the number of options that should be
recovered ($n_{V}+1$), while the blue line is the linear regression of the
actual recovered numbers and closely mimics the red line.

Figure 3: Results for Atari environments showing mean reward per episode
(after training) averaged across at least $10$ seeds for different numbers of
skills, with corresponding standard errors. Dashed lines correspond to
nonparametric models where $K$ is not specified. The label ‘ours’ denotes the
fixed-$K$ version of our options model, and the ‘-NP’ suffix denotes
nonparametric versions, both for ours and CompILE.

The goal of our experiments is twofold: to show that our options framework
learns more useful skills than DDO and CompILE, and also that the
nonparametric extensions of our own model and CompILE (which circumvent the
need to specify $K$) match the performance of their respective parametric
versions with $K$ tuned as a hyperparameter. The former goal highlights the
usefulness of incorporating variational inference advances to offline option
learning, and the latter highlights the benefits of using Bayesian
nonparametrics for skill discovery. All experimental details are given in
Appendix C.

###  6.1 Proof-of-Concept Environment

We evaluate our model in an environment we designed, where we know the number
of skills required to fully reconstruct the expert trajectories. Our
nonparametric model recovers sufficient options for reconstruction, without
this number being specified beforehand. In this environment, an agent receives
a message $m$ from a vocabulary $\\{0,...,n_{V}-1\\}$ at time $t=0$ and has to
emit the same message $m$ at time $t=4$. The agent’s observation is $(t,m)$ at
$t=0$ and $(t,-1)$ for all subsequent timesteps, and they can take a discrete
action $a_{t}\in\\{0,...,n_{V}-1\\}$ at each timestep $t$. The task is
considered a success if $a_{4}=m$. Note that this task cannot be solved by a
Markovian policy, since at time $t=4$ the agent needs to remember the message
they received at time $t=0$. However, it is possible to solve the environment
with a hierarchical policy with at least $n_{V}$ options: if every possible
message has a corresponding option whose low-level policy emits that message
at time $t=4$, then an optimal agent needs just select the appropriate option
at time $t=0$ when observing $m$ and not terminate before time $t=4$. We
generate expert trajectories for various vocabulary sizes. Considering the
heuristic we use to update $K$ (see Section 5.2), we expect our method to
recover at least $K=n_{V}+1$ ($K$ options equally used, and one unused
option). We show the results in Figure 2: indeed, we recover enough options
without overestimating $K$ by more than $1$. We also perform ablations in
Appendix D, where we show that removing the entropy regularizer from Section
4.4 significantly hurts performance.

###  6.2 Atari Environments

We further test our model on several games from the Atari learning environment
(Bellemare et al., 2013). For each game, we use expert trajectories generated
by a trained Ape-X agent (Horgan et al., 2018; Such et al., 2019). To evaluate
the quality of the discovered options, we train a PPO agent (Schulman et al.,
2017) on an environment for which the action space has been ‘augmented’ using
the skills, similarly to Fox et al. (2017). In this augmented environment,
alongside the default game actions, the agent can choose to perform ‘enhanced
actions’, which consist of acting according to one of the learned skills until
it terminates. We modify the PPO algorithm to take these enhanced actions into
account (see Appendix B). We compare our method to DDO and CompILE. We compare
both the parametric version of our model where $K$ is treated as a
hyperparameter, and the nonparametric versions of our model and CompILE.
Results are summarized in Figures 3 and 4. We can see that $(1)$ our model
learns more useful skills than DDO and CompILE (green lines consistently above
purple and orange ones, except for the Asteroids environment), and that $(2)$
the nonparametric version of our model matches the performance of tuning $K$
(dashed green lines match the peak of the green lines for all environments,
except Asterix and Space Invaders), and similarly for CompILE’s nonparametric
version (analogously for orange lines, except for Alien and Chopper Command
this time). Note that the performance of the agents in the original (i.e. not
augmented) versions of the environments correspond to the ‘$K=0$’ entries in
the figures: the significant improvements with $K>0$ highlight the relevance
of offline skill discovery. The crosses in Figures 3 and 4 show the recovered
values of $K$ (averaged across runs) for the nonparametric models, and we can
see that, except for a few cases, the number of recovered options closely
matches the smallest optimal number that would be obtained by tuning $K$
through many runs, underlining the benefits of our nonparametric approaches.
We include visualizations of our learned options in the github repository
containing our code.111https://github.com/layer6ai-
labs/BNPO/blob/main/visualization/visualization.md

We highlight that we did not tune $\lambda_{ent}$ nor its decay schedule and
used the same settings across all experiments, so that not tuning $K$ does not
come at the cost of having to tune other parameters. We present ablations with
respect to the entropy regularizer in Appendix D, where we found that while it
did not hurt performance, this regularizer was not as fundamental as in the
proof-of-concept environment. We also include in Appendix D ablations over
trajectory length, and find that while results vary for all methods, ours
remains the strongest regardless of trajectory length.

Figure 4: Additional Atari environment, same setup as in Figure 3.

##  7 Conclusions and Future Work

We introduced a novel approach for offline option discovery taking advantage
of variational inference and the Gumbel-Softmax distribution. We leave an
exploration of continuous relaxation (Kool et al., 2020; Potapczynski et al.,
2020; Paulus et al., 2020) and variational inference advances (Burda et al.,
2015; Roeder et al., 2017) within our framework for future research.

We also highlighted an unexplored connection between skill discovery and
Bayesian nonparametrics, and showed that specifying $K$ as a hyperparameter
can be avoided by placing a $\operatorname{GEM}$ prior over skills, both in
our model and other skill discovery frameworks. We hope our work will motivate
further research combining reinforcement learning and Bayesian nonparametrics,
for example through small-variance asymptotics (Kulis & Jordan, 2011; Jiang et
al., 2012; Broderick et al., 2013; Roychowdhury et al., 2013) for hard rather
than soft clustering of skills, or with hierarchical models (Teh et al.,
2006). Finally, we also hope that the technique we developed to increase $K$
throughout training will find uses in other Bayesian nonparametric models,
outside of skill discovery.

## Acknowledgements

We thank the anonymous reviewers for their feedback, which helped improve our
paper. We also thank Junfeng Wen for a useful suggestion on an ablation
experiment. Our code was written in Python (vanRossum, 1995; Oliphant, 2007)
and particularly relied on the following packages: Matplotlib (Hunter, 2007),
TensorFlow (Abadi et al., 2015) (in particular for TensorBoard), gym (Brockman
et al., 2016), Jupyter Notebook (Kluyver et al., 2016), PyTorch (Paszke et
al., 2019), NumPy (Harris et al., 2020), and Stable-Baselines3 (Raffin et al.,
2021).

## References

  * Abadi et al. (2015) Abadi, M., Agarwal, A., Barham, P., Brevdo, E., Chen, Z., Citro, C., Corrado, G. S., Davis, A., Dean, J., Devin, M., Ghemawat, S., Goodfellow, I., Harp, A., Irving, G., Isard, M., Jia, Y., Jozefowicz, R., Kaiser, L., Kudlur, M., Levenberg, J., Mané, D., Monga, R., Moore, S., Murray, D., Olah, C., Schuster, M., Shlens, J., Steiner, B., Sutskever, I., Talwar, K., Tucker, P., Vanhoucke, V., Vasudevan, V., Viégas, F., Vinyals, O., Warden, P., Wattenberg, M., Wicke, M., Yu, Y., and Zheng, X.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Achiam et al. (2018) Achiam, J., Edwards, H., Amodei, D., and Abbeel, P.  Variational option discovery algorithms.  _arXiv preprint arXiv:1807.10299_ , 2018.
  * Ajay et al. (2021) Ajay, A., Kumar, A., Agrawal, P., Levine, S., and Nachum, O.  {OPAL}: Offline primitive discovery for accelerating offline reinforcement learning.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=V69LGwJ0lIN.
  * Argall et al. (2009) Argall, B. D., Chernova, S., Veloso, M., and Browning, B.  A survey of robot learning from demonstration.  _Robotics and autonomous systems_ , 57(5):469–483, 2009.
  * Bacon et al. (2017) Bacon, P.-L., Harb, J., and Precup, D.  The option-critic architecture.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 31, 2017.
  * Bellemare et al. (2013) Bellemare, M. G., Naddaf, Y., Veness, J., and Bowling, M.  The arcade learning environment: An evaluation platform for general agents.  _Journal of Artificial Intelligence Research_ , 47:253–279, 2013.
  * Brockman et al. (2016) Brockman, G., Cheung, V., Pettersson, L., Schneider, J., Schulman, J., Tang, J., and Zaremba, W.  Openai gym.  _arXiv preprint arXiv:1606.01540_ , 2016.
  * Broderick et al. (2013) Broderick, T., Kulis, B., and Jordan, M.  Mad-bayes: Map-based asymptotic derivations from bayes.  In _International Conference on Machine Learning_ , pp. 226–234. PMLR, 2013.
  * Burda et al. (2015) Burda, Y., Grosse, R., and Salakhutdinov, R.  Importance weighted autoencoders.  _arXiv preprint arXiv:1509.00519_ , 2015.
  * Caron et al. (2007) Caron, F., Davy, M., and Doucet, A.  Generalized polya urn for time-varying dirichlet process mixtures.  In _Proceedings of the Twenty-Third Conference on Uncertainty in Artificial Intelligence_ , pp. 33–40, 2007.
  * Dempster et al. (1977) Dempster, A. P., Laird, N. M., and Rubin, D. B.  Maximum likelihood from incomplete data via the em algorithm.  _Journal of the Royal Statistical Society: Series B (Methodological)_ , 39(1):1–22, 1977.
  * Dietterich et al. (1998) Dietterich, T. G. et al.  The maxq method for hierarchical reinforcement learning.  In _ICML_ , volume 98, pp. 118–126. Citeseer, 1998.
  * Dunson & Xing (2009) Dunson, D. B. and Xing, C.  Nonparametric bayes modeling of multivariate categorical data.  _Journal of the American Statistical Association_ , 104(487):1042–1051, 2009.
  * Esmaili et al. (1995) Esmaili, N., Sammut, C., and Shirazi, G.  Behavioural cloning in control of a dynamic system.  In _1995 IEEE International Conference on Systems, Man and Cybernetics. Intelligent Systems for the 21st Century_ , volume 3, pp. 2904–2909. IEEE, 1995.
  * Eysenbach et al. (2018) Eysenbach, B., Gupta, A., Ibarz, J., and Levine, S.  Diversity is all you need: Learning skills without a reward function.  _arXiv preprint arXiv:1802.06070_ , 2018.
  * Ferguson (1973) Ferguson, T. S.  A bayesian analysis of some nonparametric problems.  _The annals of statistics_ , pp. 209–230, 1973.
  * Florensa et al. (2017) Florensa, C., Duan, Y., and Abbeel, P.  Stochastic neural networks for hierarchical reinforcement learning.  _arXiv preprint arXiv:1704.03012_ , 2017.
  * Fox et al. (2017) Fox, R., Krishnan, S., Stoica, I., and Goldberg, K.  Multi-level discovery of deep options.  _arXiv preprint arXiv:1703.08294_ , 2017.
  * Gershman & Goodman (2014) Gershman, S. and Goodman, N.  Amortized inference in probabilistic reasoning.  In _Proceedings of the annual meeting of the cognitive science society_ , volume 36, 2014.
  * Glynn (1990) Glynn, P. W.  Likelihood ratio gradient estimation for stochastic systems.  _Communications of the ACM_ , 33(10):75–84, 1990\.
  * Gregor et al. (2016) Gregor, K., Rezende, D. J., and Wierstra, D.  Variational intrinsic control.  _arXiv preprint arXiv:1611.07507_ , 2016.
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Cournapeau, D., Wieser, E., Taylor, J., Berg, S., Smith, N. J., Kern, R., Picus, M., Hoyer, S., van Kerkwijk, M. H., Brett, M., Haldane, A., del Río, J. F., Wiebe, M., Peterson, P., Gérard-Marchant, P., Sheppard, K., Reddy, T., Weckesser, W., Abbasi, H., Gohlke, C., and Oliphant, T. E.  Array programming with NumPy.  _Nature_ , 585(7825):357–362, September 2020\.  doi: 10.1038/s41586-020-2649-2.  URL https://doi.org/10.1038/s41586-020-2649-2.
  * Ho & Ermon (2016) Ho, J. and Ermon, S.  Generative adversarial imitation learning.  _Advances in neural information processing systems_ , 29:4565–4573, 2016.
  * Hochreiter & Schmidhuber (1997) Hochreiter, S. and Schmidhuber, J.  Long short-term memory.  _Neural computation_ , 9:1735–1780, 1997.
  * Horgan et al. (2018) Horgan, D., Quan, J., Budden, D., Barth-Maron, G., Hessel, M., Van Hasselt, H., and Silver, D.  Distributed prioritized experience replay.  _arXiv preprint arXiv:1803.00933_ , 2018.
  * Hunter (2007) Hunter, J. D.  Matplotlib: A 2d graphics environment.  _Computing in Science & Engineering_, 9(3):90–95, 2007.  doi: 10.1109/MCSE.2007.55.
  * Hussein et al. (2017) Hussein, A., Gaber, M. M., Elyan, E., and Jayne, C.  Imitation learning: A survey of learning methods.  _ACM Computing Surveys (CSUR)_ , 50(2):1–35, 2017\.
  * Ishwaran & James (2001) Ishwaran, H. and James, L. F.  Gibbs sampling methods for stick-breaking priors.  _Journal of the American Statistical Association_ , 96(453):161–173, 2001.
  * Jang et al. (2016) Jang, E., Gu, S., and Poole, B.  Categorical reparameterization with gumbel-softmax.  _arXiv preprint arXiv:1611.01144_ , 2016.
  * Jiang et al. (2012) Jiang, K., Kulis, B., and Jordan, M.  Small-variance asymptotics for exponential family dirichlet process mixture models.  _Advances in Neural Information Processing Systems_ , 25:3158–3166, 2012.
  * Khetarpal et al. (2020) Khetarpal, K., Klissarov, M., Chevalier-Boisvert, M., Bacon, P.-L., and Precup, D.  Options of interest: Temporal abstraction with interest functions.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 34, pp. 4444–4451, 2020.
  * Kingma & Ba (2014) Kingma, D. P. and Ba, J.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Kingma & Welling (2013) Kingma, D. P. and Welling, M.  Auto-encoding variational bayes.  _arXiv preprint arXiv:1312.6114_ , 2013.
  * Kipf et al. (2019) Kipf, T., Li, Y., Dai, H., Zambaldi, V., Sanchez-Gonzalez, A., Grefenstette, E., Kohli, P., and Battaglia, P.  Compile: Compositional imitation learning and execution.  In _International Conference on Machine Learning_ , pp. 3418–3428. PMLR, 2019.
  * Kluyver et al. (2016) Kluyver, T., Ragan-Kelley, B., Pérez, F., Granger, B., Bussonnier, M., Frederic, J., Kelley, K., Hamrick, J., Grout, J., Corlay, S., Ivanov, P., Avila, D., Abdalla, S., and Willing, C.  Jupyter notebooks – a publishing format for reproducible computational workflows.  In Loizides, F. and Schmidt, B. (eds.), _Positioning and Power in Academic Publishing: Players, Agents and Agendas_ , pp. 87 – 90. IOS Press, 2016.
  * Koller & Friedman (2009) Koller, D. and Friedman, N.  _Probabilistic graphical models: principles and techniques_.  MIT press, 2009.
  * Kool et al. (2020) Kool, W., van Hoof, H., and Welling, M.  Estimating gradients for discrete random variables by sampling without replacement.  _arXiv preprint arXiv:2002.06043_ , 2020.
  * Krishnan et al. (2018) Krishnan, S., Garg, A., Patil, S., Lea, C., Hager, G., Abbeel, P., and Goldberg, K.  Transition state clustering: Unsupervised surgical trajectory segmentation for robot learning.  In _Robotics Research_ , pp. 91–110. Springer, 2018.
  * Krishnan et al. (2019) Krishnan, S., Garg, A., Liaw, R., Thananjeyan, B., Miller, L., Pokorny, F. T., and Goldberg, K.  Swirl: A sequential windowed inverse reinforcement learning algorithm for robot tasks with delayed rewards.  _The International Journal of Robotics Research_ , 38(2-3):126–145, 2019.
  * Kulis & Jordan (2011) Kulis, B. and Jordan, M. I.  Revisiting k-means: New algorithms via bayesian nonparametrics.  _arXiv preprint arXiv:1111.0352_ , 2011.
  * Maddison et al. (2016) Maddison, C. J., Mnih, A., and Teh, Y. W.  The concrete distribution: A continuous relaxation of discrete random variables.  _arXiv preprint arXiv:1611.00712_ , 2016.
  * Mnih et al. (2015) Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A. A., Veness, J., Bellemare, M. G., Graves, A., Riedmiller, M., Fidjeland, A. K., Ostrovski, G., et al.  Human-level control through deep reinforcement learning.  _nature_ , 518(7540):529–533, 2015.
  * Murali et al. (2016) Murali, A., Garg, A., Krishnan, S., Pokorny, F. T., Abbeel, P., Darrell, T., and Goldberg, K.  Tsc-dl: Unsupervised trajectory segmentation of multi-modal surgical demonstrations with deep learning.  In _2016 IEEE International Conference on Robotics and Automation (ICRA)_ , pp. 4150–4157. IEEE, 2016.
  * Nalisnick & Smyth (2016) Nalisnick, E. and Smyth, P.  Stick-breaking variational autoencoders.  _arXiv preprint arXiv:1605.06197_ , 2016.
  * Neal (2000) Neal, R. M.  Markov chain sampling methods for dirichlet process mixture models.  _Journal of computational and graphical statistics_ , 9(2):249–265, 2000.
  * Niekum et al. (2013) Niekum, S., Chitta, S., Barto, A. G., Marthi, B., and Osentoski, S.  Incremental semantically grounded learning from demonstration.  In _Robotics: Science and Systems_ , volume 9, pp. 10–15607. Berlin, Germany, 2013.
  * Oliphant (2007) Oliphant, T. E.  Python for scientific computing.  _Computing in science & engineering_, 9(3):10–20, 2007.
  * Paszke et al. (2019) Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., et al.  Pytorch: An imperative style, high-performance deep learning library.  _Advances in neural information processing systems_ , 32:8026–8037, 2019.
  * Paulus et al. (2020) Paulus, M. B., Maddison, C. J., and Krause, A.  Rao-blackwellizing the straight-through gumbel-softmax gradient estimator.  _arXiv preprint arXiv:2010.04838_ , 2020.
  * Potapczynski et al. (2020) Potapczynski, A., Loaiza-Ganem, G., and Cunningham, J. P.  Invertible gaussian reparameterization: Revisiting the gumbel-softmax.  _Advances in Neural Information Processing Systems_ , 33, 2020.
  * Raffin et al. (2021) Raffin, A., Hill, A., Gleave, A., Kanervisto, A., Ernestus, M., and Dormann, N.  Stable-baselines3: Reliable reinforcement learning implementations.  _Journal of Machine Learning Research_ , 22(268):1–8, 2021.  URL http://jmlr.org/papers/v22/20-1364.html.
  * Rezende et al. (2014) Rezende, D. J., Mohamed, S., and Wierstra, D.  Stochastic backpropagation and approximate inference in deep generative models.  In _International conference on machine learning_ , pp. 1278–1286. PMLR, 2014.
  * Roeder et al. (2017) Roeder, G., Wu, Y., and Duvenaud, D. K.  Sticking the landing: Simple, lower-variance gradient estimators for variational inference.  _Advances in Neural Information Processing Systems_ , 30:6925–6934, 2017.
  * Ross et al. (2011) Ross, S., Gordon, G., and Bagnell, D.  A reduction of imitation learning and structured prediction to no-regret online learning.  In _Proceedings of the fourteenth international conference on artificial intelligence and statistics_ , pp. 627–635. JMLR Workshop and Conference Proceedings, 2011.
  * Roychowdhury et al. (2013) Roychowdhury, A., Jiang, K., and Kulis, B.  Small-variance asymptotics for hidden markov models.  In _Advances in Neural Information Processing Systems_ , pp. 2103–2111, 2013.
  * Schaal et al. (1997) Schaal, S. et al.  Learning from demonstration.  _Advances in neural information processing systems_ , pp. 1040–1046, 1997.
  * Schulman et al. (2017) Schulman, J., Wolski, F., Dhariwal, P., Radford, A., and Klimov, O.  Proximal policy optimization algorithms.  _arXiv preprint arXiv:1707.06347_ , 2017.
  * Sethuraman (1994) Sethuraman, J.  A constructive definition of dirichlet priors.  _Statistica sinica_ , pp. 639–650, 1994.
  * Shankar & Gupta (2020) Shankar, T. and Gupta, A.  Learning robot skills with temporal variational inference.  In _International Conference on Machine Learning_ , pp. 8624–8633. PMLR, 2020.
  * Sharma et al. (2018) Sharma, A., Sharma, M., Rhinehart, N., and Kitani, K. M.  Directed-info gail: Learning hierarchical policies from unsegmented demonstrations using directed information.  _arXiv preprint arXiv:1810.01266_ , 2018.
  * Sharma et al. (2019) Sharma, A., Gu, S., Levine, S., Kumar, V., and Hausman, K.  Dynamics-aware unsupervised discovery of skills.  _arXiv preprint arXiv:1907.01657_ , 2019.
  * Shiarlis et al. (2018) Shiarlis, K., Wulfmeier, M., Salter, S., Whiteson, S., and Posner, I.  TACO: Learning task decomposition via temporal alignment for control.  In Dy, J. and Krause, A. (eds.), _Proceedings of the 35th International Conference on Machine Learning_ , volume 80 of _Proceedings of Machine Learning Research_ , pp. 4654–4663. PMLR, 10–15 Jul 2018.  URL https://proceedings.mlr.press/v80/shiarlis18a.html.
  * Stirn et al. (2019) Stirn, A., Jebara, T., and Knowles, D.  A new distribution on the simplex with auto-encoding applications.  _Advances in Neural Information Processing Systems_ , 32:13670–13680, 2019.
  * Such et al. (2019) Such, F. P., Madhavan, V., Liu, R., Wang, R., Castro, P. S., Li, Y., Zhi, J., Schubert, L., Bellemare, M. G., Clune, J., et al.  An atari model zoo for analyzing, visualizing, and comparing deep reinforcement learning agents.  _Proceedings of IJCAI 2019_ , 2019.  URL https://github.com/uber-research/atari-model-zoo.
  * Sutton et al. (1999) Sutton, R. S., Precup, D., and Singh, S.  Between mdps and semi-mdps: A framework for temporal abstraction in reinforcement learning.  _Artificial intelligence_ , 112(1-2):181–211, 1999.
  * Teh et al. (2006) Teh, Y. W., Jordan, M. I., Beal, M. J., and Blei, D. M.  Hierarchical dirichlet processes.  _Journal of the american statistical association_ , 101(476):1566–1581, 2006.
  * vanRossum (1995) vanRossum, G.  Python reference manual.  _Department of Computer Science [CS]_ , (R 9525), 1995.
  * Vezhnevets et al. (2017) Vezhnevets, A. S., Osindero, S., Schaul, T., Heess, N., Jaderberg, M., Silver, D., and Kavukcuoglu, K.  Feudal networks for hierarchical reinforcement learning.  In _International Conference on Machine Learning_ , pp. 3540–3549. PMLR, 2017.
  * Williams (1992) Williams, R. J.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  _Machine learning_ , 8(3):229–256, 1992.

##  Appendix A Relaxations

We choose the following relaxations:

$\displaystyle\delta_{b_{0}^{(i)}=1}\approx b_{0}^{(i)}$ (22)
$\displaystyle\delta_{h_{t}^{(i)}=h_{t-1}^{(i)}}\approx
1-||h_{t}^{(i)}-h_{t-1}^{(i)}||_{1}/2$ (23)
$\displaystyle\eta(h_{t}^{(i)})\approx\eta^{\top}h_{t}^{(i)}$ (24)
$\displaystyle
p_{\theta}(b_{t}^{(i)},h_{t}^{(i)}|h_{t-1}^{(i)},s_{t}^{(i)})\approx
b_{t}^{(i)}\psi_{h_{t-1}^{(i)}}(s_{t}^{(i)})\eta^{\top}h_{t}^{(i)}+(1-b_{t}^{(i)})(1-\psi_{h_{t-1}^{(i)}}(s_{t}^{(i)}))\cdot(1-||h_{t}^{(i)}-h_{t-1}^{(i)}||_{1}/2)$
(25)

Note that these relaxations $(1)$ match the objective they are relaxing when
$b$’s are binary and $h$’s are one-hot; and $(2)$ that they are not unique:
there are many sensible choices that could be used.

##  Appendix B Learning in Augmented Environments

In the augmented environment, taking an action corresponding to one of the
skills means that the agent may actually interact with the environment for
several timesteps. Naïvely, one might think that simply using the usual
Bellman backup update considering these enhanced actions as actions would be
sufficient. However, this neglects the fact that the reward credited to the
enhanced action arises from multiple environment interactions, and that the
discounting of the $Q$-value of the next state depends on the length of the
enhanced action. Recall that the usual Bellman equation is:

$\displaystyle Q(s,a)$
$\displaystyle=\mathbb{E}_{\pi}\left[\displaystyle\sum_{t=1}^{\infty}\gamma^{t-1}R_{t}|S_{0}=s,A_{0}=a\right]$
(26) $\displaystyle=\mathbb{E}_{\pi}\left[R_{1}+\gamma
Q(S_{1},A_{1})|S_{0}=s,A_{0}=a\right]$ (27)

where the first term in the last equality ($R_{1}$) corresponds to the reward
from a single timestep, and the second term ($\gamma Q(S_{1},A_{1})$) to the
discounted return of all the following states.

For an augmented environment, let action $a$ result in $\tau$ interactions
with the (unaugmented) environment. Note that $\tau$ is a random variable,
since skill terminations are stochastic, though of course $\tau=1$
deterministically if $a$ is a primitive action. We then have:

$Q(s,a)=\mathbb{E}_{\pi,\tau}\left[\displaystyle\sum_{t=1}^{\tau}\gamma^{t-1}R_{t}+\gamma^{\tau}Q(S_{\tau},A_{\tau})|S_{0}=s,A_{0}=a\right]$
(28)

As before, $t$ indexes interactions with the unaugmented environment during
execution of a single enhanced action. The first term
($\sum_{t=1}^{\tau}\gamma^{t-1}R_{t}$) corresponds to the reward directly
credited to the (possibly enhanced) action, and the second term
($\gamma^{\tau}Q(S_{\tau},A_{\tau})$) is the discounted expected return from
future (possibly enhanced) actions. As a sanity check, note that in the case
where $\tau$ is always $1$, both equations match.

This modified Bellman equation implies that we must keep a slightly modified
replay buffer. The typical buffer $\\{(s_{t},a_{t},s_{t+1},r_{t})\\}_{t}$,
storing every interaction with the unaugmented environment, is insufficient.
We would be unable to determine when enhanced actions were taken or when they
terminated, and therefore unable to credit them accordingly. Further, such a
buffer is wasteful, assuming that enhanced actions are common. Instead we
should keep a buffer
$\\{(s_{t},a_{t},s_{t+\tau_{t}},\sum_{t^{\prime}=t}^{t+\tau_{t}}\gamma^{t^{\prime}-t}r_{t^{\prime}},\tau_{t})\\}$.
Here $a_{t}$ denotes the action taken in the augmented action space at time
$t$ (possibly a primitive action, possibly an enhanced action). $\tau_{t}$ is
the (random) duration of that action. We must keep $\tau_{t}$ to be able to
discount $Q$ on the right hand side of the Bellman equation appropriately. We
need only put entries into the buffer at timesteps at which we take an action
in the augmented environment. We do not need to record every interaction with
the unaugmented environment.

##  Appendix C Experimental Details

#### Proof-of-concept environment

For this experiment, we generate $1000$ expert trajectories using a manually
designed policy. We then train our model for $500$ epochs. The options’ sub-
policies and termination functions consist of MLPs with two hidden layers of
16 units separated by a ReLU activation and followed by a Softmax activation.
The parameters of all layers except the last one are shared across options.
For the encoder, we use an LSTM layer with 32 hidden units, and MLPs with two
hidden layers of 32 units and ReLU activations for both heads, with weights
shared except for the last layer. We use a learning rate of 0.005 with the
Adam optimizer (Kingma & Ba, 2014) and a batch size of 128. The GS temperature
parameter is initialized at 1 and annealed by a factor of 0.995 each epoch.
$\lambda_{ent}$ is initialized at 5 and also annealed by a factor of 0.995
each epoch. We check the option usage every 10 epochs ($n_{K}=10^{4}$) and use
$\delta=0.5$ for our rule of increasing $K$.

#### Atari environments

For these experiments, we use $1000$ expert trajectories of length $300$ and
train for $1500$ epochs. We use the RAM state of the Atari environments, where
each state is a 128 bytes array, that we unpack into a 1024 bit array. The
options’ sub-policies and termination functions consist of a single linear
layer with a Softmax activation per option. This choice is motivated by the
fact that we do not want a single sub-policy to be able to fully reconstruct
the expert trajectories. We use the same encoder as for the proof-of-concept
environment, and the same optimizer with a learning rate of 0.001. The GS
temperature is now annealed by a factor of 0.999. We keep the same values for
$n_{K}$ and $\delta$. We use similar architectures for sub-policies and
termination in both CompILE and DDO, as well as for the encoder in CompILE. We
fix the number of segments in CompILE to 7 for all runs.

To learn in the augmented environment, we use the PPO agent implemented by
Raffin et al. (2021), with the default ’CnnPolicy’ that takes as input the
image state of the environment with the same preprocessing as done by Mnih et
al. (2015). We use the implementation default parameters except for the
’n_steps’ variable (the number of environment steps used per update) that we
set to 512. We also modify the replay buffer used during training to take into
account the specific aspects of learning in an augmented environment mentioned
in Appendix B. This agent is trained for $3\cdot 10^{6}$ (augmented) steps and
the mean reward across the last $5\cdot 10^{5}$ steps is used in Fig. 3. For
CompILE, we automatically terminate each enhanced action after 15 time steps
for all augmented environments, as doing so was preferable to following the
Poisson-sampled termination.

##  Appendix D Additional Experiments

###  D.1 Proof-of-Concept Environment

Figure 5: Results on our proof-of-concept environment, without entropy
regularization (left panel), and with fixed entropy regularization (right
panel).

We show in Figure 5 results analogous to those of Figure 2, except we do not
use the entropy regularizer from Section 4.4 (left panel), or we simply do not
anneal it (right panel). We can see that, as mentioned in the main manuscript,
not using the regularizer significantly degrades performance, although not
using annealing (and keeping the regularizing coefficient fixed throughout
training) does not have much of an impact.

###  D.2 Atari Environments

As mentioned in the main manuscript, we ablate some of our choices. Figure 6
shows our ablation results of using the entropy regularizer, green curves show
results with the regularizer, and blue ones without. Across environments,
performance is comparable, so that the entropy regularizer is not truly needed
here. We highlight this result is opposite to what we found in our proof-of-
concept environment, where the regularizer was fundamental to getting our
method to work. Since adding the regularizer does not hurt performance in
Atari environments and greatly helps in our proof-of-concept environment, we
nonetheless recommend to use it as a default.

Tables 1, 2, and 3 show the results of ablations where the length of expert
trajectories is changed for a subset of the Atari environments that we
considered. We highlight that we did not cherry pick these environments, and
the fact that we do not show analogous results for the missing environments
was merely a matter of computational costs. While results do change
significantly by varying trajectory length, both for our methods and the
baselines, we can see that: $(1)$ our fixed-$K$ method consistently
outperforms or remains competitive with both CompILE and DDO, the only
exception being Montezuma’s revenge with trajectories of length $50$, and
$(2)$ our nonparametric method outperforms or remains competitive with
nonparametric CompILE across all settings. We also highlight that DDO ran out
of memory when using trajectories of length $1000$. We can thus see that our
empirical superiority shown in the main manuscript was not due to a lucky
choice of expert trajectory length.

Figure 6: Evolution of the mean reward per episode during training with and without the entropy regularizer. Results are averaged across 5 runs. Table 1: Reward on Montezuma’s revenge, varying the length of expert trajectories. Columns requiring $K$, i.e. “Ours”,“CompILE”, and “DDO”, use $K=7$. Length | Ours | CompILE | DDO | Ours-NP | CompILE-NP
---|---|---|---|---|---
50 | 53.8 $\pm$ 7.6 | 113.1 $\pm$ 65.6 | 41.8 $\pm$ 20.9 | 274.0 $\pm$ 51.6 | 149.9 $\pm$ 73.7
150 | 350.9 $\pm$ 44.8 | 0.3 $\pm$ 0.3 | 24.3 $\pm$ 19.2 | 344.9 $\pm$ 44.6 | 102.7 $\pm$ 91.7
300 | 317.4 $\pm$ 48.6 | 0.0 $\pm$ 0.0 | 0.3 $\pm$ 0.2 | 261.0 $\pm$ 38.8 | 92.6 $\pm$ 80.0
500 | 350.9 $\pm$ 42.5 | 0.0 $\pm$ 0.0 | 27.0 $\pm$ 0.0 | 391.2 $\pm$ 31.8 | 55.1 $\pm$ 48.4
1000 | 328.3 $\pm$ 86.2 | 0.0 $\pm$ 0.0 | NA | 610.9 $\pm$ 74.4 | 66.6 $\pm$ 45.0
Table 2: Reward on Breakout, varying the length of expert trajectories. Columns requiring $K$ use $K=7$ for “Ours” and “CompILE” and $K=11$ for “DDO”. Length | Ours | CompILE | DDO | Ours-NP | CompILE-NP
---|---|---|---|---|---
50 | 32.8 $\pm$ 1.0 | 16.0 $\pm$ 2.9 | 27.4 $\pm$ 1.3 | 31.6 $\pm$ 1.5 | 25.9 $\pm$ 2.6
150 | 27.1 $\pm$ 1.3 | 20.0 $\pm$ 1.2 | 25.3 $\pm$ 1.9 | 22.9 $\pm$ 1.1 | 21.5 $\pm$ 1.1
300 | 36.6 $\pm$ 3.0 | 18.6 $\pm$ 2.1 | 26.5 $\pm$ 1.7 | 31.6 $\pm$ 2.9 | 27.4 $\pm$ 2.9
500 | 28.0 $\pm$ 2.0 | 18.4 $\pm$ 1.6 | 20.6 $\pm$ 0.0 | 22.6 $\pm$ 1.1 | 24.0 $\pm$ 3.4
1000 | 23.0 $\pm$ 1.7 | 16.4 $\pm$ 3.2 | NA | 19.8 $\pm$ 1.1 | 17.9 $\pm$ 0.9
Table 3: Reward on Space Invaders, varying the length of expert trajectories. Columns requiring $K$, i.e. “Ours”,“CompILE”, and “DDO”, use $K=7$. Length | Ours | CompILE | DDO | Ours-NP | CompILE-NP
---|---|---|---|---|---
50 | 488.1 $\pm$ 10.5 | 367.7 $\pm$ 18.1 | 414.1 $\pm$ 15.6 | 418.2 $\pm$ 26.6 | 440.2 $\pm$ 17.8
150 | 562.0 $\pm$ 21.6 | 401.2 $\pm$ 14.8 | 475.0 $\pm$ 17.3 | 492.0 $\pm$ 16.1 | 447.3 $\pm$ 20.3
300 | 583.0 $\pm$ 16.2 | 438.8 $\pm$ 45.6 | 479.3 $\pm$ 19.9 | 531.0 $\pm$ 17.7 | 469.9 $\pm$ 15.9
500 | 562.9 $\pm$ 18.9 | 432.6 $\pm$ 29.2 | 468.8 $\pm$ 27.1 | 478.4 $\pm$ 12.9 | 469.2 $\pm$ 17.2
1000 | 552.5 $\pm$ 19.4 | 420.8 $\pm$ 24.4 | NA | 490.8 $\pm$ 16.3 | 502.3 $\pm$ 25.5


# Constant-roll inflation from a fermionic field

Mehdi Shokri<EMAIL_ADDRESS>Department of Physics, Campus of
Bijar, University of Kurdistan, Bijar, Iran  Jafar Sadeghi<EMAIL_ADDRESS>Department of Physics, University of Mazandaran, P. O. Box 47416-95447,
Babolsar, Iran  Mohammad Reza Setare<EMAIL_ADDRESS>Department of Physics,
Campus of Bijar, University of Kurdistan, Bijar, Iran

###### Abstract

We study the inflationary period driven by a fermionic field which is non-
minimally coupled to gravity in the context of the constant-roll approach. We
consider the model for a specific form of coupling and perform the
corresponding inflationary analysis. By comparing the result with the Planck
observations coming from CMB anisotropies, we find the observational
constraints on the parameters space of the model and also the predictions the
model. We find that the values of $r$ and $n_{s}$ for $-1.5<\beta\leq-0.9$ are
in good agreement with the observations when $|\xi|=0.1$ and $N=60$.

PACS: 98.80.Cq; 98.80.$-$K.
Keywords: Fermionic field; Constant-roll inflation; Cosmic Microwave
Background.

##  I Introduction

Cosmic inflation has been identified as the most appealing theory removing
shortcomings of the hot big bang model and also the main reliable approach for
structure formation of the universe. Furthermore, primordial gravitational
waves are produced through primordial perturbations and they could be
observed, in principle, in the stochastic background Guth ; Sato:1980yn ;
Kazanas:1980tx ; Linde:1981my ; Albrecht:1982wi ; Lyth:1998xn . Based on the
standard formalism, a single scalar field (inflaton) can drive inflation, and
then by decaying at the last stage of the period, inflation will be terminated
throughout a reheating process Kofman2 ; Shtanov . Several inflationary models
are equipped with the slow-roll approximation which says inflaton rolls slowly
down from the top of potential to the minimum point with a canonical kinetic
term. For a long time, the single field inflationary models have been
investigated in literature in detail so that by comparing with the Planck
data, some of them are restricted or even ruled out martin . However, some
models are still viable with the recent high precision observations staro ;
barrow ; bezrukov ; kallosh1 . One of the most significant properties of these
models is that they do not present any non-Gaussianity in their primordial
spectrum because of uncorrelated modes of the spectrum Chen . In such a case,
if the future observations predict non-Gaussianities in the perturbations
spectrum, then this class of models will be felt in serious uncertainty. In
order to remove the expected obstacles, a new class of inflationary models has
been suggested in which inflaton is rolling down with a constant rate during
inflation martin2 ; Motohashi1 ; Motohashi2 . This means that $\ddot{\varphi}$
is non-negligible and can be expressed as

$\ddot{\varphi}=\beta H\dot{\varphi}$ (1)

where $\beta=-(3+\alpha)$ and $\alpha$ is a non-zero parameter. For
$\alpha=-3$, the model is reduced to the standard slow-roll. The need to go
beyond of the standard slow-roll approximation goes back to the ultra slow-
roll regime where the term of $\ddot{\varphi}$ is non-negligible in the Klein-
Gordon equation as $\ddot{\varphi}=3H\dot{\varphi}$. This class of
inflationary models shows a finite value for the non-decaying mode of
curvature perturbations and can be expressed with the familiar form of
standard inflation Inoue . Although the ultra slow-roll model predicts a large
$\eta$ because of the non-negligible term $\ddot{\varphi}$, it is not able to
solve the $\eta$ problem introduced in supergravity for the hybrid
inflationary models Kinney . Moreover, the inflationary solutions of the ultra
model are situated in the non-attractor phase of inflation but show a scale-
invariance for the scalar perturbation spectrum. Besides, the main problem of
ultra models is that the non-Gaussinaity consistency relation of single field
models is violated in the presence of ultra condition through super-Hubble
evolution of the scalar perturbation Namjoo . Despite the above problems, the
ultra slow-roll inflationary model is really interesting since it reveals the
plenty of unknown and unexplained issues around the inflationary dynamics.
Fast-roll models are known as another class of inflationary models introduced
in the beyond of the slow-roll approximation in which a fast-rolling stage is
considered at the start of inflation and will be connected to the standard
slow-roll only after a few e-folds Contaldi ; Lello ; Hazra . This approach
can be considered as a correction to slow-roll models in order to improve the
consistency with observations.

Constant-roll approach has opened a new window to analyze cosmic inflation due
to the consideration of a constant rate of rolling for inflaton. The model is
becoming very popular since a wide range of inflationary models has been
investigated in the context of constant-roll inflation Odintsov ; Nojiri ;
Motohashi9 ; Cicciarella ; Awad ; Anguelova ; Ito ; Yi ; Morse ; karam1 ;
Ghersi ; Lin ; Micu ; Oliveros ; Motohashi3 ; Kamali ; diego ; setare ; new2 .
In the present manuscript, we study a fermionic inflationary model in which a
fermionic field which is non-minimally coupled to gravity is the main
responsible to drive inflation instead of inflaton introduced in the standard
inflationary paradigm. In fact, Noether symmetry arguments of the model show
an exponentially expanding universe if a non-minimal coupling between the
fermionic field and the Ricci scalar is considered kramers1 . Fermionic fields
also are known as gravitational sources of dark energy as the late-time
accelerating of the universe in addition to playing the role of inflaton for
the early time acceleration kramers2 . Moreover, the interaction of a
fermionic field with a scalar field can be recognized as the source of
interacting dark matter-dark energy lepe . Such theories that deal with
fermionic fields instead of scalar fields are named the Fermions Tensor
Theories (FTT) q1 ; q2 ; q3 ; q4 .

The main aim of this paper is studying the constant-roll inflation for a
fermionic inflationary model in which a fermionic field is non-minimally
coupled to gravity. We investigate the model for a specific form of coupling
and then perform the corresponding inflationary analysis. Finally, we compare
the obtained results with the Planck data coming from CMB anisotopies in order
to find the observational constraints on the parameters space of the model.
The above discussion motivates us to arrange the paper as follows. In §II, we
introduce the properties of the FTT with the corresponding dynamical
equations. In §III, we present the fermionic inflationary model in the context
of the constant-roll inflation for the coupling form $f\sim 1+\xi\phi$. In
§IV, we analyze the obtained results of the model by comparing with the
observational datasets from the Planck and the BICEP2/Keck array satellites.
In §V, we conclude the analysis of the models and draw the possible outlooks.

##  II Fermions Tensor Theory

Let us start with the inflationary action for a fermionic field $\psi$ non-
minimally coupled to to gravity as follows Bojowald

$S=\int{d^{4}x\sqrt{-g}\bigg{(}\frac{1}{2}f(\phi)R+\frac{i}{2}\big{(}\bar{\psi}\tilde{\gamma}^{\mu}D_{\mu}\psi-(D_{\mu}\bar{\psi})\tilde{\gamma}^{\mu}\psi\big{)}}-V(\phi)\bigg{)}$
(2)

where $g$ is the determinant of the metric $g_{\mu\nu}$,
$R=g^{\mu\nu}R_{\mu\nu}$ is the Ricci scalar Ricci scalar, $\psi$ and
$\bar{\psi}=\psi^{\dagger}\gamma^{0}$ are the fermionic field and its adjoint
field, respectively. $\phi\equiv|\bar{\psi}\psi|$ is the scalar bilinear and
$f$ is a function that couples the fermionic field to gravity. Also,
$\tilde{\gamma}^{\mu}=e^{\mu}_{\nu}\gamma^{\nu}$ are the generalized gamma
matrices where $e^{\mu}_{\nu}$ are the tetrad fields and $D_{\mu}$ is the
covariant derivative expressed by

$D_{\mu}\psi=\partial_{\mu}\psi-\Omega_{\mu}\psi,\hskip
28.45274ptD_{\mu}\bar{\psi}=\partial_{\mu}\bar{\psi}+\bar{\psi}\Omega_{\mu}\hskip
28.45274ptwith\hskip
28.45274pt\Omega_{\mu}=-\frac{1}{4}g_{\rho\sigma}\big{(}\Gamma^{\rho}_{\mu\delta}-e^{\rho}_{b}(\partial_{\mu}e^{b}_{\delta})\big{)}\tilde{\gamma}^{\delta}\tilde{\gamma}^{\sigma}$
(3)

where $\Omega_{\mu}$ is the spin connection and $\Gamma^{\mu}_{\nu\delta}$ is
the Christoffel symbol. Moreover, $V$ as the fermionic potential can be
written in terms of the fermion mass $m$ as follows

$V=m\bar{\psi}\psi.$ (4)

In the following, we consider a spatially flat universe described by a
Friedmann-Robertson-Walker (FRW) metric as
$ds^{2}=-dt^{2}+a(t)^{2}(dx^{2}+dy^{2}+dz^{2})$ where $t$ and $a$ depict to
cosmic time and scale factor, respectively. Therefore, we can state the
effective Lagrangian of the model by

$\mathcal{L}=3a\dot{a}\big{(}\dot{a}f+a\dot{\phi}f^{\prime}\big{)}+\frac{i}{2}a^{3}\big{(}\dot{\bar{\psi}}\gamma^{0}\psi-\bar{\psi}\gamma^{0}\dot{\psi}\big{)}+a^{3}V$
(5)

where dot and prime refer to the derivative with respect to the cosmic time
$t$ and the bilinear field $\phi$, respectively. Using the Euler-Lagrange
equation of $\bar{\psi}$ and $\psi$ for the Lagrangian (5), we obtain the
following equations of motion

$\dot{\bar{\psi}}+\frac{3}{2}H\bar{\psi}-iV^{\prime}_{0}\bar{\psi}\gamma^{0}+6if^{\prime}\bar{\psi}\gamma^{0}(\dot{H}+2H^{2})=0$
(6) $\dot{\psi}+\frac{3}{2}H\psi-
iV^{\prime}_{0}\gamma^{0}\psi+6if^{\prime}\gamma^{0}\psi(\dot{H}+2H^{2})=0.$
(7)

Multiplying the Eqs. (6) and (7) by $\psi$ and $\bar{\psi}$, respectively and
then adding two equations, the equation of motion of the bilinear field $\phi$
takes the following form

$\frac{d\phi}{dt}+3H\phi=0$ (8)

where $H=\dot{a}/a$ denotes the Hubble parameter. Also, by using the Euler-
Lagrange equation of $a$ and setting the energy function of the Lagrangian (5)
as zero, we obtain the acceleration and Friedmann equations, respectively as

$\frac{\ddot{a}}{a}=-\frac{\rho_{f}+3p_{f}}{6f},\hskip
28.45274ptH^{2}=\frac{\rho_{f}}{3f}$ (9)

where the energy density $\rho_{f}$ and the pressure $p$ of the fermionic
field are given by

$\rho_{f}=V-3Hf^{\prime}\dot{\phi},\hskip
28.45274ptp_{f}=\phi\big{(}V^{\prime}-3f^{\prime}(\dot{H}+2H^{2})\big{)}-V+\ddot{\phi}f^{\prime}+2H\dot{\phi}f^{\prime}+f^{\prime\prime}\dot{\phi}^{2}.$
(10)

Plugging the above expressions to the Eq. (9) and using the the definition of
the Hubble constant $H=\dot{a}/a$, the dynamical equations are obtained as

$3fH^{2}=V-3Hf^{\prime}\dot{\phi}$ (11) $-2f\dot{H}=\phi
V^{\prime}-Hf^{\prime}\dot{\phi}-3f^{\prime}\phi(\dot{H}+2H^{2})+\ddot{\phi}f^{\prime}+f^{\prime\prime}\dot{\varphi}^{2}$
(12)

where again dot and prime denotes the derivative with respect to $t$ and
$\phi$.

##  III Fermionic Constant-roll Inflation

Now that we know the properties of the FTT, let’s study a fermionic
inflationary model in the context of the constant-roll idea. By setting the
constant-roll condition (1) and then using $\dot{H}=\dot{\phi}H^{\prime}$, the
Friedmann equation (11) and the equation of motion of the bilinear field (8),
from (12) we have

$H(\phi)=C\phi^{-\frac{\beta+3}{3}}$ (13)

where $C$ is the integration constant. Also, from the Eq. (8), the evolution
of the blinear field is given by

$\phi(t)=\big{(}-(\beta+3)Ct\big{)}^{\frac{3}{3+\beta}}.$ (14)

Since Noether symmetry discussions of the model show an exponential expansion
for the universe if $f^{\prime}\neq 0$ kramers1 ; kramers2 , we consider the
coupling form

$f=1+\xi\phi$ (15)

where $\xi$ is the coupling constant. From the Eq. (11), we obtain the
corresponding potential as

$V(\phi)=3C^{2}(1-2\xi\phi)\phi^{-\frac{2(\beta+3)}{3}}$ (16)

which $\phi_{critical}=\frac{1}{2\xi}$ is the barrier that the bilinear field
can not cross it. In Fig. 1, we present the behaviour of the above potential
for different values of $\beta$ when $|\xi|=0.1$ in which the bilinear field
$\phi$ sharply rolls down from top of the potential to the minimum point. Now,
by using the conformal transformation $\hat{g}_{\mu\nu}=\Omega^{2}g_{\mu\nu}$
where the conformal factor $\Omega=\Omega(\phi(x))$ is a non-null and
differentiable function, we can move from the Jordan frame to the Einstein
frame with the conformal factor $\Omega^{2}=f$. Therefore, the action in the
Einstein frame is given by

$S_{E}=\int
d^{4}x\sqrt{-\hat{g}}\bigg{(}\frac{\hat{R}}{2}-\frac{i}{2}F^{2}\big{(}\hat{\bar{\psi}}\tilde{\gamma}^{\mu}\hat{D}_{\mu}\hat{\psi}-(\hat{D}_{\mu}\hat{\bar{\psi}})\tilde{\gamma}^{\mu}\hat{\psi}\big{)}-\hat{V}(\hat{\phi})\bigg{)}.$
(17)

In the new frame, we deal with a redefined bilinear field $\hat{\phi}$ and its
potential

$F^{2}\equiv\bigg{(}\frac{d\hat{\phi}}{d\phi}\bigg{)}^{2}=\frac{2f+3f^{\prime
2}}{2f^{2}},\quad\quad\hat{V}(\hat{\phi})\equiv\frac{V(\phi)}{\Omega^{4}}=\frac{V(\phi)}{f^{2}}.$
(18)

Now, we can calculate the slow-roll parameters of the model

$\hat{\epsilon}\equiv\frac{1}{2}\bigg{(}\frac{\hat{V}^{\prime}}{\hat{V}}\bigg{)}^{2}=\frac{4\Big{(}\phi^{2}(2\beta+9)\xi^{2}+\phi(\beta-3)\xi-(\beta+3)\Big{)}^{2}}{9\phi^{2}\big{(}2\phi\xi+3\xi^{2}+2\big{)}(2\phi\xi-1)^{2}},$
(19)

Figure 1: The potential (16) plotted versus the bilinear field $\phi$ for
different values of $\beta$ when $|\xi|=0.1$.

$\displaystyle\\!\hat{\eta}\equiv\frac{\hat{V}^{\prime\prime}}{\hat{V}}=\frac{1}{9\phi^{2}\big{(}2\phi\xi+3\xi^{2}+2\big{)}^{2}(2\phi\xi-1)}\bigg{\\{}12(2\beta+9)^{2}\phi^{3}\xi^{5}+4\Big{(}(8\beta^{2}+78\beta+189)\phi^{2}+18\beta^{2}+90\beta-27\Big{)}\phi^{2}\xi^{4}+$
$\displaystyle\\!+\Big{(}(80\beta^{2}+564\beta+648)\phi^{3}-(252\beta+756)\phi\Big{)}\xi^{3}+12(2\beta+9)\Big{(}2(\beta-3)\phi^{2}-(\beta+3)\Big{)}\xi^{2}-4(\beta+3)\phi(4\beta+63)\xi-$
$\displaystyle\\!-(16\beta^{2}+120\beta+216)\bigg{\\}},$ (20)
$\displaystyle\\!\hat{\zeta}^{2}\equiv\frac{\hat{V}^{\prime}\hat{V}^{\prime\prime\prime}}{\hat{V}^{2}}=\frac{16\Big{(}\phi^{2}(2\beta+9)\xi^{2}+\phi(\beta-3)\xi-(\beta+3)\Big{)}}{567\phi^{4}\big{(}2\phi\xi+3\xi^{2}+2\big{)}^{4}(2\phi\xi-1)^{2}}\Bigg{\\{}63(2\beta+9)^{3}\phi^{4}\xi^{8}+21\Big{(}2(8\beta+33)(\beta+6)(2\beta+9)\phi^{2}+$
$\displaystyle\\!+60\beta^{3}+648\beta^{2}+1836\beta+405\Big{)}\phi^{3}\xi^{7}+\Big{(}28(\beta+6)(2\beta+9)(4\beta+21)\phi^{4}+21(\beta+6)\big{(}112\beta^{2}+714\beta+675\big{)}\phi^{2}+$
$\displaystyle\\!+189(4\beta^{2}+6\beta-111)(\beta+3)\Big{)}\phi^{2}\xi^{6}+63\Big{(}2(\beta+6)\big{(}8\beta^{2}+62\beta+87\big{)}\phi^{4}+\big{(}\frac{128}{3}\beta^{3}+378\beta^{2}+471\beta-1350\big{)}\phi^{2}-2$
$\displaystyle\\!\times(2\beta+9)(2\beta+39)(\beta+3)\Big{)}\phi\xi^{5}+7\Big{(}8\big{(}333\beta^{2}+891\beta+30\beta^{3}\big{)}\phi^{4}+(32\beta^{2}-330\beta-2997)(\beta+3)\phi^{2}-18(\beta+6)$
$\displaystyle\\!\times(2\beta+9)(\beta+3)\Big{)}\xi^{4}+7\Big{(}2(80\beta^{3}+468\beta^{2}-1422\beta-7371)\phi^{2}-3(2\beta+9)(\beta+3)(16\beta+213)\Big{)}\phi\xi^{3}-168$
$\displaystyle\\!\times\Big{(}\frac{3}{2}(26\beta+147)\phi^{2}+2\beta^{2}+21\beta+54\Big{)}(\beta+3)\xi^{2}-42(4\beta+45)(2\beta+9)(\beta+3)\phi\xi-56(\beta+6)(2\beta+9)(\beta+3)\Bigg{\\}}$
(21)

where prime implies to the derivative with respect to the bilinear field
$\phi$ and inflation ends when the condition $\hat{\epsilon}=1$ or
$\hat{\eta}=1$ is fulfilled. The number of e-folds is given by

$\hat{N}\equiv\int^{\hat{\phi}_{i}}_{\hat{\phi}_{f}}{\frac{1}{\sqrt{2\hat{\epsilon}}}}d\hat{\phi}$
(22)

and also the spectral parameters, i.e. the spectral index, its running and the
tensor-to-scalar ratio are defined by

$\hat{n}_{s}=1-6\hat{\epsilon}+2\hat{\eta},\quad\quad\quad\hat{\alpha}_{s}=16\hat{\epsilon}{\hat{\eta}}-24\hat{\epsilon}^{2}-2\hat{\zeta}^{2},\quad\quad\quad\hat{r}=16\hat{\epsilon}.$
(23)

Figure 2: The marginalized joint 68% and 95% CL regions for $n_{s}$ and $r$ at
$k=0.002$ Mpc-1 from Planck alone and in combination with BK15 or BK15+BAO
data cmb  and the $n_{s}-r$ constraints on the model with the coupling (15).
The dashed and solid lines represent $N=50$ and $N=60$, respectively. The
results are obtained for different values of $\beta$ when $|\xi|=0.1$.

##  IV Constraints from CMB Anisotropies

Let’s compare the obtained results with the observational datasets coming from
CMB anisotropies. In Figure 2, we present the $n_{s}-r$ constraints coming
from the marginalized joint 68% and 95% CL regions of the Planck 2018 alone
and in combination with BK15 or BK15+BAO data on the fermionic inflationary
model with the coupling (15) introduced in the context of the constant-roll
approach. The figure is drawn for different values of $\beta$ in the cases
$N=50$ (dashed line) and $N=60$ (solid line) when $|\xi|=0.1$. By considering
the CMB anisotropies datasets from the Planck alone, we find that all values
of $\beta$ in the case of $N=50$ show a reasonable tensor-to-scalar ratio $r$
below the upper limit while the corresponding values of the spectral index
$n_{s}$ are not in good agreement with observations. The situation in the case
of $N=60$ is much better than $N=50$. As we can see, the cases of $\beta>-0.9$
and $\beta\leq-1.5$ are excluded since the obtained values of $r$ and $n_{s}$
are out of the observational region. The cases of $-1.15\leq\beta\leq-0.9$ and
$-1.4<\beta<-1.15$ show observationally acceptable values of $r$ and $n_{s}$
at the 68% and 95% CL, respectively. By combination of the BK15 data and the
Planck data, the situation of $N=50$ is the similar to the previous
observational case and all values of $\beta$ provide disfavored values of $r$
and $n_{s}$. For $N=60$ the cases of $\beta\geq-1.1$ and $\beta\leq-1.5$ show
undesirable values of $r$ and $n_{s}$ while the cases of $-1.25\leq\beta<-1.1$
and $-1.4<\beta<-1.25$ present the values compatible with the observations at
the 68% and $95\%$ CL, respectively. As a full consideration, we compare the
obtained results with the observations coming from the Planck in combination
with BK15+BAO data. In such a case, still $N=50$ is fully excluded by the
observations. Also, $\beta>-1.1$ and $\beta<-1.3$ are not in good agreement
with data while the cases of $-1.25\leq\beta\leq-1.1$ and
$-1.3\leq\beta<-1.25$ show desirable values of $r$ and $n_{s}$ at the 68% and
95% CL, respectively.

##  V Conclusion

In this paper, we have studied the constant-roll inflation driven by a
fermionic field instead of inflaton which is non-minimally coupled to gravity.
First, we have introduced the Fermions Tensors Theories (FTT) briefly and then
have investigated the constant-roll condition for the model with the coupling
form $f\sim 1+\xi\varphi$ which has been approved by the Noether symmetry. At
the next step, we calculated the spectral parameters of the model and finally
have compared the results with the CMB observations in order to find the
observational constraints on the parameters space, in particular, the
constant-roll parameter $\beta$. We have found that the obtained values of $r$
and $n_{s}$ for $-1.5<\beta\leq-0.9$ are fully compatible with the CMB
anisotriopes observations when $|\xi|=0.1$ and $N=60$.

## References

  * (1) A. H. Guth, “The inflationary universe: A possible solution to the horizon and flatness problems,” Phys. Rev. D, vol. 23, p. 347, (1981).
  * (2) K. Sato, “First order phase transition of a vacuum and expansion of the universe,” Mon. Not. Roy. Astron. Soc., vol. 195, p. 467, (1981).
  * (3) D. Kazanas, “Dynamics of the universe and spontaneous symmetry breaking,” Astrophys. J., vol. 241, p. L59, (1980).
  * (4) A. D. Linde, “A new inflationary universe scenario: A possible solution of the horizon, flatness, homogeneity, isotropy problems,” Phys. Lett. B, vol. 108, p. 389, (1982).
  * (5) A. Albrecht and P. J. Steinhardt, “Cosmology for grand unified theories with radiatively induced symmetry breaking,” Phys. Rev. Lett., vol. 48, p. 1220, (1982).
  * (6) D. H. Lyth and A. Riotto, “Particle physics models of inflation and the perturbation,” Phys. Rept., vol. 314, p. 1, (1999).
  * (7) L. Kofman, A. D. Linde, and A. A. Starobinsky, “Reheating after inflation,” Phys. Rev. Lett., vol. 73, p. 3195, (1994).
  * (8) Y. Shtanov, J. H. Traschen, and R. H. Brandenberger, “Universe reheating after inflation,” Phys. Rev. D, vol. 51, p. 5438, (1995).
  * (9) J. Martin, “What have the planck data taught us about inflation?,” Class. Quant. Grav., vol. 33, p. 034001, (2016).
  * (10) A. A. Starobinsky, “A new type of isotropic cosmological models without singularity,” Phys. Lett. B, vol. 91, p. 99, (1980).
  * (11) J. D. Barrow and S. Cotsakis, “Inflation and the conformal structure of higher order gravity theories,” Phys. Lett. B, vol. 214, p. 515, (1988).
  * (12) F. L. Bezrukov and M. Shaposhnikov, “The standard model higgs boson as the inflaton,” Phys. Lett. B, vol. 659, p. 703, (2008).
  * (13) R. Kallosh and A. Linde, “Universality class in conformal inflation,” JCAP, vol. 07, p. 002, (2013).
  * (14) X. Chen, “Primordial non-gaussianities from inflation models,” Adv. Astron., vol. 2010, p. 638979, (2010).
  * (15) J. Martin, H. Motohashi, and T. Suyama, “Primordial non-gaussianities from inflation models,” Phys. Rev. D, vol. 87, p. 023514, (2013).
  * (16) H. Motohashi, A. A. Starobinsky, and J. Yokoyama, “Inflation with a constant rate of roll,” JCAP, vol. 09, p. 018, (2015).
  * (17) H. Motohashi and A. A. Starobinsky, “Constant-roll inflation: confrontation with recent observational data,” Europhys. Lett., vol. 117, p. 39001, (2017).
  * (18) S. Inoue and J. Yokoyama, “Curvature perturbation at the local extremum of the inflaton’s potential,” Phys. Lett. B, vol. 524, p. 15, (2002).
  * (19) W. H. Kinney, “Horizon crossing and inflation with large $\eta$,” Phys. Rev. D, vol. 72, p. 023515, (2005).
  * (20) M. H. Namjoo, H. Firouzjahi, and M. Sasaki, “Violation of non-gaussianity consistency relation in a single field inflationary model,” Europhys. Lett., vol. 101, p. 39001, (2013).
  * (21) C. R. Contaldi, L. K. M. Peloso, and A. D. Linde, “Suppressing the lower multipoles in the cmb anisotropies,” JCAP, vol. 0307, p. 002, (2003).
  * (22) L. Lello and D. Boyanovsky, “Tensor to scalar ratio and large scale power suppression from pre-slow roll initial conditions,” JCAP, vol. 1405, p. 029, (2014).
  * (23) D. K. Hazra, A. Shafieloo, G. F. Smoot, , and A. A. Starobinsky, “Whipped inflation,” Phys. Rev. Lett., vol. 113, p. 071301, (2014).
  * (24) S. Odintsov and V. Oikonomou, “Inflationary dynamics with a smooth slow-roll to constant-roll era transition,” Phys. Rev. D, vol. 96, p. 024029, (2017).
  * (25) S. Nojiri, D. Odintsov, and V. K. Oikonomou, “Constant-roll inflation in f(R) gravity,” Class. Quant. Grav., vol. 34, p. 245012, (2017).
  * (26) H. Motohashi and A. A. Starobinsky, “f(R) constant-roll inflation,” Eur. Phys. J. C, vol. 77, p. 538, (2017).
  * (27) F. Cicciarella, J. Mabillard, and M. Pieroni, “New perspectives on constant-roll inflation,” JCAP, vol. 01, p. 024, (2018).
  * (28) A. Awad, W. E. Hanafy, G. Nashed, S. Odintsov, and V. Oikonomou, “Constant-roll inflation in f(T) teleparallel gravity,” JCAP, vol. 07, p. 026, (2017).
  * (29) L. Anguelova, P. Suranyi, and L. Wijewardhana, “Systematics of constant roll inflation,” JCAP, vol. 02, p. 004, (2018).
  * (30) A. Ito and J. Soda, “Anisotropic constant-roll inflation,” Eur. Phys. J. C, vol. 78, p. 55, (2018).
  * (31) Z. Yi and Y. Gong, “On the constant-roll inflation,” JCAP, vol. 1803, p. 052, (2018).
  * (32) M. J. P. Morse and W. H. Kinney, “Large $\eta$ constant-roll inflation is never an attractor,” Phys. Rev. D, vol. 97, p. 123519, (2018).
  * (33) A. Karam, L. Marzola, T. Pappas, A. Racioppi, and K. Tamvakis, “Constant-roll (quasi-)linear inflation,” JCAP, vol. 05, p. 011, (2018).
  * (34) J. T. G. Ghersi, A. Zucca, y, and A. V. Frolov, “Observational constraints on constant roll inflation,” JCAP, vol. 05, p. 030, (2019).
  * (35) W.-C. Lin and M. J. P. Morsey, “Dynamical analysis of attractor behavior in constant roll inflation,” JCAP, vol. 09, p. 063, (2019).
  * (36) A. Micu, “Two-field constant roll inflation,” JCAP, vol. 11, p. 003, (2019).
  * (37) A. Oliveros and H. E. Noriega, “Constant-roll inflation driven by a scalar field with nonminimal derivative coupling,” Int. J. Mod. Phys. D, vol. 28, p. 1950159, (2019).
  * (38) H. Motohashi and A. A. Starobinsky, “Constant-roll inflation in scalar-tensor gravity,” JCAP, vol. 11, p. 025, (2019).
  * (39) V. Kamali, M. Artymowski, and M. R. Setare, “Constant roll warm inflation in high dissipative regime,” JCAP, vol. 07, p. 002, (2020).
  * (40) M. Guerrero, D. Rubiera-Garcia, and D. S.-C. Gomez, “Constant roll inflation in multifield models,” Phys. Rev. D, vol. 102, p. 123528, (2020).
  * (41) M. Shokri, J. Sadeghi, M. R. Setare, and S. Capozziello, “Nonminimal coupling inflation with constant slow roll,” To be publish in Int. J. Mod. Phys. D [gr-qc:2104.00596].
  * (42) M. Shokri, J. Sadeghi, and M. R. Setare, “The generalized $sl(2,\textit{R})$ and $su(1,1)$ in non-minimal constant-roll inflation,” Ann. Phys., vol. 429, p. 168487, (2021).
  * (43) R. C. de Souza and G. M. Kremer, “Noether symmetry for non-minimally coupled fermion fields,” Class. Quant. Grav., vol. 25, p. 225006, (2008).
  * (44) G. Grams, R. C. de Souza, and G. M. Kremer, “Fermion field as inflaton, dark energy and dark matter,” Class. Quant. Grav., vol. 31, p. 185008, (2014).
  * (45) S. Lepe, J. Lorca, F. Pena, and Y. Vasquez, “Fermionic and scalar fields as sources of interacting dark matter-dark energy,” Int. J. Mod. Phys. D, vol. 20, p. 2543, (2011).
  * (46) A. Kumar, “Inflation and reheating with a fermionic field,” [gr-qc:1811.12237].
  * (47) D. Benisty, E. I. Guendelman, E. N. Saridakis, H. Stoecker, J. Struckmeier, and D. Vasak, “Inflation from fermions with curvature-dependent mass,” Phys. Rev. D, vol. 100, p. 043523, (2019).
  * (48) D. Benisty, “Decaying coupled fermions to curvature and the $\textit{H}_{0}$ tension,” [gr-qc: 1912.11124].
  * (49) B. Saha, “Non-minimally coupled nonlinear spinor field in FRW cosmology,” Astrophys.Space Sci., vol. 365, p. 68, (2020).
  * (50) M. Bojowald and R. Das, “Canonical gravity with fermions,” Phys. Rev. D, vol. 78, p. 064009, (2008).
  * (51) Y. Akrami et al., “Planck 2018 results. constraints on inflation,” [astro-ph: 1807.06211].