# A nearly Blackwell-optimal policy gradient method

Vektor Dewanto, Marcus Gallagher
School of Information Technology and Electrical Engineering
University of Queensland, Australia
<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Abstract

For continuing environments, reinforcement learning methods commonly maximize
a discounted reward criterion with discount factor close to 1 in order to
approximate the steady-state reward (the gain). However, such a criterion only
considers the long-run performance, ignoring the transient behaviour. In this
work, we develop a policy gradient method that optimizes the gain, then the
bias (which indicates the transient performance and is important to capably
select from policies with equal gain). We derive expressions that enable
sampling for the gradient of the bias, and its preconditioning Fisher matrix.
We further propose an algorithm that solves the corresponding bi-level
optimization using a logarithmic barrier. Experimental results provide
insights into the fundamental mechanisms of our proposal.

##  1 Introduction

Consider an environment in Fig 1, for which any stationary policy
$\pi\in\Pi_{\mathrm{S}}$ has equal gain values (the expected long-run average
reward in steady-state), namely $v_{g}(\pi)=-1$. Therefore, all stationary
policies are gain optimal, satisfying $v_{g}(\pi)\geq v_{g}(\pi^{\prime}),\
\forall\pi,\pi^{\prime}\in\Pi_{\mathrm{S}}$. In this case, the gain optimality
(which is equivalent to $(n=-1)$-discount optimality) is underselective
because there is no way to discriminate policies based on their gain values
(to carry out policy selection). Nonetheless, there exists a more selective
optimality criterion that can be obtained by increasing the value of $n$ from
$-1$ to $0$ (or greater, if needed) in the family of $n$-discount optimality
(puterman_1994_mdp, Ch 10).

$s^{0}$$s^{1}$$(a^{0},5,0.5)$$(a^{0},5,0.5)$$(a^{1},r=10,p=1)$$(a^{2},-1,1)$
Figure 1: A two-state environment with two policies whose gains are equally
$-1$ (puterman_1994_mdp, p500). Edge labels tell the action $a$, reward $r$,
and transition probability $p$.

It turns out that for the environment in Fig 1, $(n=0)$-discount optimality is
the most selective because its optimal policy, i.e. selecting action $a^{0}$,
is also optimal with respect to all other $n$-discount optimality criteria for
all $n>0$. This selectiveness is due to consideration of the transient reward,
in addition to the gain (steady-state reward). It can be achieved by
$\gamma$-discounted optimality _only if_ the discount factor
$\gamma>10/11\approx 0.909$. Such a critical lower threshold for $\gamma$ is
difficult to estimate for environments with unknown dynamics. Moreover, there
is an implementation-dependent upper threshold such that the discounted value
$v_{\gamma}(\pi)$ is numerically not too close to the gain $v_{g}(\pi)$.
Otherwise, it becomes underselective as well because $\lim_{\gamma\to
1}(1-\gamma)v_{\gamma}(\pi)=v_{g}(\pi)$ (puterman_1994_mdp, Cor 8.2.5).

We believe that one always desires the most selective criterion. This can be
obtained by solving ($n=|\mathcal{S}|-2$, at most)-discount optimality for
some type of environments with a finite state set $\mathcal{S}$, as shown by
puterman_1994_mdp. As a first step in leveraging this fact in reinforcement
learning (RL), we focus on $(n=0)$-discount optimality, for which we have the
notion of nearly-Blackwell optimality blackwell_1962_ddp; feinberg_2002_hmdp,
Def 3.1 as follows.

######  Definition 1.1.

A stationary policy $\pi_{b}^{*}\in\Pi_{\mathrm{S}}$ is nearly-Blackwell
optimal (nBw-optimal, also termed bias optimal) if it is gain optimal, and in
addition

$v_{b}(\pi_{b}^{*},s)\geq v_{b}(\pi_{g}^{*},s),\quad\forall
s\in\mathcal{S},\forall\pi_{g}^{*}\in\Pi_{g}^{*},\quad\text{with}\
\Pi_{g}^{*}\coloneqq\operatorname*{\arg\\!\max}_{\pi\in\Pi_{\mathrm{S}}}v_{g}(\pi),\
\text{and}$ $\text{the bias value, i.e.~{}}\
v_{b}(\pi,s)\coloneqq\lim_{t_{\mathrm{max}}\to\infty}\mathbb{E}\mathopen{}\mathclose{{}\left[\sum_{t=0}^{t_{\mathrm{max}}-1}\big{(}r(S_{t},A_{t})-v_{g}(\pi)\big{)}\Big{|}S_{0}=s,\pi}\right],$
(1)

where $v_{g}$ is the gain value and $r$ is a reward function. The above limit
is assumed to exist.

To our knowledge, mahadevan_1996_sensitivedo proposed the first (and the only)
tabular Q-learning that can obtain nBw-optimal policies through optimizing the
family of $n$-discount optimality. It relies on stochastic approximation (SA)
to estimate the optimal values of gain $v_{g}^{*}(=v_{-1}^{*})$, bias
$v_{b}^{*}(=v_{0}^{*})$, and bias-offset $v_{1}^{*}$. Since those values are
in two _nested_ equations, standard SA is not guaranteed to converge, as noted
by the author. Moreover, it is likely to become unstable if used in conjuction
with function approximators, similar to other value-iteration approaches that
generally operate off-policy (sutton_2018_irl, Ch 11.3). For other related
works in the broader scope, refer to Appendix LABEL:app:rwork.

We take the alternative policy-iteration approach: a policy gradient method.
It works by maximizing the expected value of a policy with respect to some
initial state distribution, i.e.
$\mathbb{E}_{p_{s_{0}}}\mathopen{}\mathclose{{}\left[v(\pi,S_{0})}\right]$.
This requires _stochastic_ optimization on a _non-concave_ landscape with
gradient ascent update given by

$\boldsymbol{\theta}\leftarrow\boldsymbol{\theta}+\alpha\ \boldsymbol{C}^{-1}\
\nabla v(\boldsymbol{\theta},s_{0}),\quad\text{ for a parameterized policy
$\pi(\boldsymbol{\theta}),\forall\boldsymbol{\theta}\in\Theta=\real{\dim(\boldsymbol{\theta})}$,
}$ (2)

where $\alpha$ is a positive step length,
$\boldsymbol{C}\in\real{\dim(\boldsymbol{\theta})\times\dim(\boldsymbol{\theta})}$
is a positive-definite preconditioning matrix, $s_{0}$ is a deterministic
initial state (without loss of generality),
$\nabla\coloneqq\partial/\partial\boldsymbol{\theta}$, and
$v(\boldsymbol{\theta},s_{0})\coloneqq v(\pi(\boldsymbol{\theta}),s_{0})$. We
assume that $\pi(\boldsymbol{\theta})$ is twice continuously differentiable
with bounded first and second derivatives.

In essense, our nBw-optimal policy gradient method switches the optimization
objective from gain $v_{g}$ to bias $v_{b}$ once the gain optimization is
deemed to have converged. This needs three novel ingredients, to which we
contribute. First is a sampling-based approximator for the (standard, vanilla)
gradient of bias (Sec 2). Second is a suitable preconditioning matrix
$\boldsymbol{C}$, for which we define the Fisher (information) matrix of bias,
yielding natural gradients of bias (Sec 3). Third is an algorithm that
utilizes the gradient and Fisher estimates to carry out gain then bias (hence,
bi-level) optimization (Sec 4). To this end, we formulate a logarithmic
barrier function that ensures the optimization iterates stay in the gain-
optimal region (as per Def 1.1). The experiment results in Sec 5 provide
insights into the fundamental mechanisms of our proposal and comparison with
the discounted reward method.

#### Preliminaries and notations:

We model the agent-environment interaction as a Markov decision process (MDP)
with a finite state set $\mathcal{S}$, a finite action set $\mathcal{A}$, and
an infinite decision-epoch set $\mathcal{T}$. At each decision-epoch (discrete
timestep) $t\in\mathcal{T}$ for $t=0,1,\ldots,t_{\mathrm{max}}=\infty$, the
agent is at a state $s_{t}\in\mathcal{S}$ and chooses to then execute an
action $a_{t}\in\mathcal{A}$. Consequently, it arrives at a next state
$s_{t+1}\in\mathcal{S}$ and earns a finite reward $r_{t+1}\in\real{}$. Here,
$S_{0}\sim p_{s_{0}}$, $S_{t+1}\sim p(\cdot|s_{t},a_{t})$, and
$r_{t+1}=r(s_{t},a_{t})=\mathbb{E}\mathopen{}\mathclose{{}\left[r(s_{t},a_{t},S_{t+1})}\right]<\infty$,
where $p_{s_{0}}$ denotes the initial state distribution, $p$ the transition
distribution, and $r$ the reward function (these three entities specify the
time-homogenous environment dynamics). Each action comes from
$a_{t}\sim\pi(\cdot|s_{t})$, where
$\pi:\mathcal{S}\times\mathcal{A}\mapsto[0,1]$ is called a policy. In this
work, any $\pi$ belongs to a stationary policy set $\Pi_{\mathrm{S}}$ which
consists of the stationary randomized set $\Pi_{\mathrm{SR}}$ and its
degenerate counterpart, the stationary deterministic set $\Pi_{\mathrm{SD}}$.

An MDP induces Markov chains (MCs), one for each stationary policy $\pi$. Each
induced MC has a stochastic transition matrix under $\pi$, denoted as
$\boldsymbol{P}_{\pi}$ and of size $|\mathcal{S}|$-by-$|\mathcal{S}|$. Raising
it to the $t$-th power gives $\boldsymbol{P}_{\pi}^{t}$, whose
$[s_{0},s_{t}]$-th component indicates the $t$-step transition probability
$p_{\pi}^{t}(s_{t}|s_{0})\coloneqq\mathrm{Pr}\\{S_{t}=s_{t}|s_{0};\pi\\}$. The
limiting matrix is defined as
$\boldsymbol{P}_{\pi}^{\star}\coloneqq\lim_{t\to\infty}\boldsymbol{P}_{\pi}^{t}$,
which exists whenever the MC is aperiodic. Its component
$p_{\pi}^{\star}(s|s_{0})$ can be interpreted as the steady-state (stationary,
equilibrium) probability of being at state $s$ when starting at initial state
$s_{0}$. In some MCs, stationarity is achieved in finite time. A mixing time
measures the time required by an MC such that the distance to the stationarity
is small, i.e.
$t_{\mathrm{mix}}^{\pi}(\varepsilon)\coloneqq\min\\{t:\max_{s_{0}\in\mathcal{S}}\\{\frac{1}{2}\sum_{s\in\mathcal{S}}|\boldsymbol{P}_{\pi}^{t}[s_{0},s]-\boldsymbol{P}_{\pi}^{\star}[s_{0},s]|\\}\leq\varepsilon\\}$.

In this work, the modeling MDP is unichain, meaning that every induced MC
consists of a single recurrent class _plus_ a (possible empty) set of
transient states. A recurrent class is a _closed_ set of states, where
$\boldsymbol{P}_{\pi}^{t}(s^{\prime}|s)>0$ for some $t\geq 0$, for all
$s,s^{\prime}\in\mathcal{S}$. A transient state is a non-recurrent state since
after some finite time, such a state will never be visited again. In
particular, $t_{\mathrm{abs}}^{\pi|s_{0}}$ indicates the time when the chain
is absorbed into the recurrent class given that it starts at $s_{0}$ and
follows $\pi$.

##  2 Gradients of the bias

The bias (1) from every initial state $s_{0}\in\mathcal{S}$ is equivalent to

$\displaystyle v_{b}(\pi,s_{0})$
$\displaystyle=\sum_{s\in\mathcal{S}}\sum_{a\in\mathcal{A}}z_{\pi}(s|s_{0})\pi(a|s)r(s,a)-\sum_{s\in\mathcal{S}}z_{\pi}(s|s_{0})v_{g}(\pi),\quad\text{hence}$
$\displaystyle\nabla v_{b}(\pi,s_{0})$
$\displaystyle=\sum_{s\in\mathcal{S}}\sum_{a\in\mathcal{A}}z(s|s_{0})\pi(a|s)r(s,a)\big{\\{}\nabla\log
z_{\pi}(s|s_{0})+\nabla\log\pi(a|s)\big{\\}}$
$\displaystyle\quad-\sum_{s\in\mathcal{S}}z_{\pi}(s|s_{0})\big{\\{}v_{g}(\pi)\nabla\log
z_{\pi}(s|s_{0})+\nabla v_{g}(\pi)\big{\\}},$ (cf. that of gain in
(LABEL:equ:gain_grad_naive))

where $z_{\pi}(s|s_{0})$ is the $[s_{0},s]$-component of a non-stochastic
matrix
$(\boldsymbol{I}-\boldsymbol{P}_{\pi}+\boldsymbol{P}_{\pi}^{\star})^{-1}$.
Because $z_{\pi}$ is _not_ a probability distribution, the above expression
does not enable sampling-based approximation. Therefore, we derive the
following theorem that gives us a score-function gradient estimator for bias.

######  Theorem 2.1.

The gradient of the bias of a randomized policy
$\pi(\boldsymbol{\theta})\in\Pi_{\mathrm{SR}}$ is given by

$\displaystyle\nabla v_{b}(\boldsymbol{\theta},s_{0})$
$\displaystyle={\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Big{\\{}}\underbrace{\sum_{t=0}^{t_{\mathrm{mix}}^{\pi}-1}\Big{\\{}\mathbb{E}_{S_{t}\sim
p_{\pi}^{t}(\cdot|s_{0}),A_{t}\sim\pi(\cdot|S_{t};\boldsymbol{\theta})}\Big{[}q_{b}^{\pi}(S_{t},A_{t})\nabla\log\pi(A_{t}|S_{t};\boldsymbol{\theta})\Big{]}-\nabla
v_{g}(\boldsymbol{\theta})\Big{\\}}}_{\text{The pre-mixing
part}}{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\Big{\\}}}$
$\displaystyle\quad+{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\Big{\\{}}\underbrace{\mathbb{E}_{S_{t_{\mathrm{mix}}^{\pi}}\sim
p_{\pi}^{\star},A_{t_{\mathrm{mix}}^{\pi}}\sim\pi(\cdot|S_{t_{\mathrm{mix}}^{\pi}};\boldsymbol{\theta})}\Big{[}\nabla
q_{b}^{\pi}(S_{t_{\mathrm{mix}}^{\pi}},A_{t_{\mathrm{mix}}^{\pi}})\Big{]}+\nabla
v_{g}(\boldsymbol{\theta})}_{\text{The post-mixing
part}}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\Big{\\}}},$
(3)

for all $s_{0}\in\mathcal{S}$, and for all $\boldsymbol{\theta}\in\Theta$,
where $t_{\mathrm{mix}}^{\pi}\coloneqq t_{\mathrm{mix}}^{\pi}(\varepsilon=0)$,
and $q_{b}^{\pi}$ denotes the bias state-action value of $\pi$ (which is
defined correspondingly to the bias state value $v_{b}^{\pi}$ in (1)).

###### Proof.

See Appendix LABEL:app:biasgrad. ∎

The bias gradient expression in (3) has pre- and post-mixing parts. Both
involve the gain gradient $\nabla v_{g}(\boldsymbol{\theta})$, which is
$\boldsymbol{0}$ at every point in the gain-optimal region (in which the bias
optimization should be carried out as per Def 1.1). In Appendix
LABEL:app:gradbias_decompose, we provide empirical contributions of each term
from $t=0$ to $t_{\mathrm{mix}}^{\pi}$ by plotting the error between the
partial summation so-far and the exact $\nabla v_{b}(\boldsymbol{\theta})$. In
general, the error drops significantly due to the first few pre-mixing terms,
then levels-off until (just before) mixing. This agrees with the limiting
behaviour of the deviation of $p_{\pi}^{t}$ from $p_{\pi}^{\star}$. Finally,
another error drop occurs (notably in some cases) due to a single post-mixing
term. This indicates that truncation on the pre-mixing part is possible with
relatively small errors to lower the variance in sampling-based approximation.
For such a purpose, we can also substract some action-independent baseline,
e.g. $v_{b}^{\pi}$, from $q_{b}^{\pi}$ (leading to the use of the so-called
bias state-action advantages).

In the following, we explain each part of (3) further, beginning with the
post-mixing part as it shares a similarity with the existing gain gradient
expression (sutton_2000_pgfnapprox, Thm 1): both require state samples from
the stationary state distribution $p_{\pi}^{\star}$ (and the same policy
evaluation quantity, i.e. $q_{b}^{\pi}$).

#### The post-mixing part:

A single unbiased state sample from $p_{\pi}^{\star}$ can be obtained at
$t_{\mathrm{mix}}^{\pi}$. It can also be sampled at any time after mixing, but
multiple samples will not be independent (they are Markovian and identically
distributed, yielding _biased_ sample-mean gradient estimates).

In practice, $S_{t_{\mathrm{mix}}^{\pi}}\sim p_{\pi}^{\star}$ occurs at the
last timestep of an experiment-episode, assuming that the length of an
experiment-episode is equal or larger than $t_{\mathrm{mix}}^{\pi}+1$. Thus,
an unbiased sampling-based estimate of $\nabla v_{b}(\boldsymbol{\theta})$ is
available after every experiment-episode, yielding experiment-episode-wise
policy parameter updates (2), as shown in Algo
LABEL:alg:gainbiasbarrier_optim. We use the term “experiment-episode
($\mathrm{xep}$)” because the environments of interest do not have a natural
notion of episodes; they are non-episodic (continuing). Such environments are
modeled as infinite-horizon MDPs with $t_{\mathrm{max}}=\infty$, which can be
approximated by a finite length experiment-episode with maximum timestep
$\hat{t}_{\mathrm{max}}^{\mathrm{xep}}$.

As observed in (3), the post-mixing part requires that the estimator for
$q_{b}^{\pi}$ (termed the critic) depends on the policy parameter
$\boldsymbol{\theta}$. This can be achieved, for example, by using the
following state-action features for the critic,
$\boldsymbol{f}_{\boldsymbol{\theta}}(s,a)\coloneqq\nabla\log\pi(a|s;\boldsymbol{\theta}),\forall(s,a)\in\mathcal{S}\times\mathcal{A}$.
There also exists a class of critic that takes as input the policy parameter
(faccio_2021_pvf; harb_2020_pen).

Alternatively, whenever it is preferable to use _disjoint_ policy and critic
parameter sets (where there is no shared parameters), one can utilize our Thm
2.2 below. It presents an equivalent form of the post-mixing part that is free
from $\nabla
q_{b}^{\pi}(S_{t_{\mathrm{mix}}^{\pi}},A_{t_{\mathrm{mix}}^{\pi}})$.

######  Theorem 2.2.

The post-mixing part of (3), denoted as $\nabla
v_{b,post}(\boldsymbol{\theta},s_{0})$, has the following identity,

$\nabla
v_{b,post}(\boldsymbol{\theta},s_{0})=\mathbb{E}_{S_{t_{\mathrm{mix}}^{\pi}}\sim
p_{\pi}^{\star},A_{t_{\mathrm{mix}}^{\pi}}\sim\pi(\cdot|S_{t_{\mathrm{mix}}^{\pi}};\boldsymbol{\theta})}\Big{[}q_{1}^{\pi}(S_{t_{\mathrm{mix}}^{\pi}},A_{t_{\mathrm{mix}}^{\pi}})\nabla\log\pi(A_{t_{\mathrm{mix}}^{\pi}}|S_{t_{\mathrm{mix}}^{\pi}};\boldsymbol{\theta})\Big{]},$

where $q_{1}^{\pi}$ denotes the state-action value of a policy $\pi$ in terms
of the ($n=1$)-discount optimality.

###### Proof.

See Appendix LABEL:app:postmix_alt. ∎

#### The pre-mixing part:

The state samples for this part of (3) are all but the last state of an
experiment-episode. That is, $S_{0}\sim p_{\pi}^{0}$, $S_{1}\sim p_{\pi}^{1}$,
$\ldots$, $S_{\hat{t}_{\mathrm{max}}^{\mathrm{xep}}-1}\sim
p_{\pi}^{\hat{t}_{\mathrm{max}}^{\mathrm{xep}}-1}$, where
$\hat{t}_{\mathrm{max}}^{\mathrm{xep}}=t_{\mathrm{mix}}^{\pi}$. This is
because the state distributions from which states are sampled at every
timestep, i.e. $p_{\pi}^{t}(\cdot|s_{0})$, are the correct distributions of
each corresponding term in the pre-mixing part.

Algo LABEL:alg:gradfisher_sampling shows a way to obtain pre-mixing samples,
which are unbiased whenever $\hat{t}_{\mathrm{max}}^{\mathrm{xep}}\geq
t_{\mathrm{mix}}^{\pi}$. It postpones the substraction of $\nabla
v_{g}(\boldsymbol{\theta})$ in each pre-mixing term until the gain gradient
estimate is available, i.e. at the last step of an experiment-episode. Note
that if an experiment-episode runs longer than $t_{\mathrm{mix}}^{\pi}$, the
sum of pre-mixing terms from $t=t_{\mathrm{mix}}^{\pi}$ to
$\hat{t}_{\mathrm{max}}^{\mathrm{xep}}$ is equal to $\boldsymbol{0}$ (in exact
cases).

##  3 Natural gradients of the bias

The standard gradient ascent of the form (2) without preconditioner is _not_
invariant under policy parameterization. To overcome this, we need natural
(covariant) gradients of the bias, i.e. $\boldsymbol{F}_{b}^{-1}\nabla
v_{b}(\boldsymbol{\theta})$, where $\boldsymbol{F}_{b}$ is a bias Fisher
matrix derived in this Section. This $\boldsymbol{F}_{b}$ also increases the
convergence rate.

kakade_2002_npg firstly proposed a gain Fisher $\boldsymbol{F}_{g,a}$ of a
parameterized policy $\pi(\boldsymbol{\theta})$. That is,

$\boldsymbol{F}_{g,a}(\boldsymbol{\theta})\coloneqq\sum_{s\in\mathcal{S}}p_{\pi}^{\star}(s)\boldsymbol{F}_{a}(\boldsymbol{\theta},s),\quad\text{with}\
\boldsymbol{F}_{a}(\boldsymbol{\theta},s)\coloneqq\underbrace{\sum_{a\in\mathcal{A}}\pi(a|s;\boldsymbol{\theta})\nabla\log\pi(a|s;\boldsymbol{\theta})\nabla^{\intercal}\log\pi(a|s;\boldsymbol{\theta}),}_{\text{The
action Fisher, which is independent of optimality criteria}}$ (4)

where
$\boldsymbol{F}_{a}(\boldsymbol{\theta},s)\in\real{\dim(\boldsymbol{\theta})\times\dim(\boldsymbol{\theta})}$
is a positive semidefinite Fisher information matrix that defines a semi-
Riemannian metric (i.e. an inner product on the tangent space) on the action
probability manifold (“surface”) at a state $s\in\mathcal{S}$ (the term metric
is _not_ in its typical sense as a distance function). Note that
$\frac{1}{2}\delta_{\boldsymbol{\theta}}^{\intercal}\boldsymbol{F}_{a}(\boldsymbol{\theta},s)\delta_{\boldsymbol{\theta}}+\mathcal{O}\mathopen{}\mathclose{{}\left(\delta_{\boldsymbol{\theta}}^{3}}\right)$
is the second-order Taylor expansion of the Kullback-Leibler divergence at a
state $s$ between two policies parameterized by $\boldsymbol{\theta}$ and
$\boldsymbol{\theta}+\delta_{\boldsymbol{\theta}}$ for an infinitesimal
$\delta_{\boldsymbol{\theta}}$.

The gain Fisher (4) is defined as the weighted sum of the action Fisher
$\boldsymbol{F}_{a}$ whose weights are the components of
$\boldsymbol{P}_{\pi}^{\star}$, by which the gain is formulated as
$v_{g}(\pi)\cdot\boldsymbol{1}=\boldsymbol{P}_{\pi}^{\star}\boldsymbol{r}_{\pi}$,
where $\boldsymbol{r}_{\pi}\in\real{|\mathcal{S}|}$ is a reward vector with
elements
$\boldsymbol{r}_{\pi}[s]=\mathbb{E}_{A\sim\pi(\cdot|s)}\mathopen{}\mathclose{{}\left[r(s,A)}\right],\forall
s\in\mathcal{S}$, and $\boldsymbol{P}_{\pi}^{\star}$ has identical rows in
unichain MDPs. Now, the bias is given by
$\boldsymbol{v}_{b}(\pi)=\boldsymbol{D}_{\pi}\boldsymbol{r}_{\pi}$, where
$\boldsymbol{D}_{\pi}\in\real{|\mathcal{S}|\times|\mathcal{S}|}$ is a non-
stochastic deviation matrix of a policy $\pi$. Therefore, we expect that the
absolute value of the $[s_{0},s]$-component of $\boldsymbol{D}_{\pi}$ can be
used to weight $\boldsymbol{F}_{a}$ to yield a bias Fisher. We need to capture
only the deviation magnitude in order to maintain the positive
semidefiniteness of $\boldsymbol{F}_{a}$.

From puterman_1994_mdp, we know the following about the deviation matrix,

$\boldsymbol{D}_{\pi}=(\boldsymbol{I}-\boldsymbol{P}_{\pi}+\boldsymbol{P}_{\pi}^{\star})^{-1}(\boldsymbol{I}-\boldsymbol{P}_{\pi}^{\star})=\lim_{t_{\mathrm{max}}\to\infty}\sum_{t=0}^{t_{\mathrm{max}}-1}(\boldsymbol{P}_{\pi}^{t}-\boldsymbol{P}_{\pi}^{\star}).$
(5)

Since the induced MC is aperiodic, the limit (5) exists. Hence, we have

$|d_{\pi}(s|s_{0})|=\Big{|}\sum_{t=0}^{\infty}\big{(}p_{\pi}^{t}(s|s_{0})-p_{\pi}^{\star}(s|\cancel{s_{0}})\big{)}\Big{|}\leq\sum_{t=0}^{\infty}|p_{\pi}^{t}(s|s_{0})-p_{\pi}^{\star}(s)|,\qquad\forall
s,s_{0}\in\mathcal{S},$ (6)

where $d_{\pi}(s|s_{0})$, $p_{\pi}^{t}(s|s_{0})$, and
$p_{\pi}^{\star}(s|\cancel{s_{0}})$ denote the $[s_{0},s]$-entry of
$\boldsymbol{D}_{\pi}$, $\boldsymbol{P}_{\pi}^{t}$, and
$\boldsymbol{P}_{\pi}^{\star}$, respectively. Taking the absolute value per
timestep is beneficial for trajectory-sampling-based approximation in model-
free RL, as explained in the remainder of this Section.

We are now ready to define the bias Fisher $\boldsymbol{F}_{b,a}$ as follows.
For all $\pi(\boldsymbol{\theta})\in\Pi_{\mathrm{SR}}$ and for all
$s_{0}\in\mathcal{S}$,

$\displaystyle\boldsymbol{F}_{b,a}(\boldsymbol{\theta}|s_{0})\coloneqq\sum_{s\in\mathcal{S}}\Big{\\{}\underbrace{\sum_{t=0}^{\infty}|p_{\pi}^{t}(s|s_{0})-p_{\pi}^{\star}(s)|}_{\text{$\geq|d_{\pi}(s|s_{0})|$
in
\eqref{equ:deviation_abs}}}\Big{\\}}\boldsymbol{F}_{a}(\boldsymbol{\theta},s)=\sum_{t=0}^{\infty}\Big{\\{}\sum_{s\in\mathcal{S}}\underbrace{|p_{\pi}^{t}(s|s_{0})-p_{\pi}^{\star}(s)|}_{\text{a
weight (not a
probability)}}\boldsymbol{F}_{a}(\boldsymbol{\theta},s)\Big{\\}}$
$\displaystyle=\sum_{t=0}^{t_{\mathrm{abs}}^{\pi|s_{0}}-1}\Big{\\{}\sum_{s_{\mathrm{tr}}^{\pi}\in\mathcal{S}_{\mathrm{tr}}^{\pi}}|p_{\pi}^{t}(s_{\mathrm{tr}}^{\pi}|s_{0})-\cancelto{0}{p_{\pi}^{\star}(s_{\mathrm{tr}}^{\pi})}|\quad\boldsymbol{F}_{a}(\boldsymbol{\theta},s_{\mathrm{tr}}^{\pi})+\sum_{s_{\mathrm{rc}}^{\pi}\in\mathcal{S}_{\mathrm{rc}}^{\pi}}|\cancelto{0}{p_{\pi}^{t}(s_{\mathrm{rc}}^{\pi}|s_{0})}-p_{\pi}^{\star}(s_{\mathrm{rc}}^{\pi})|\
\boldsymbol{F}_{a}(\boldsymbol{\theta},s_{\mathrm{rc}}^{\pi})\Big{\\}}$
$\displaystyle\quad+\sum_{t=t_{\mathrm{abs}}^{\pi|s_{0}}}^{t_{\mathrm{mix}}^{\pi}-1}\Big{\\{}\sum_{s_{\mathrm{tr}}^{\pi}\in\mathcal{S}_{\mathrm{tr}}^{\pi}}\cancelto{0}{|p_{\pi}^{t}(s_{\mathrm{tr}}^{\pi}|s_{0})-p_{\pi}^{\star}(s_{\mathrm{tr}}^{\pi})|}\
\boldsymbol{F}_{a}(\boldsymbol{\theta},s_{\mathrm{tr}}^{\pi})+\underbrace{\sum_{s_{\mathrm{rc}}^{\pi}\in\mathcal{S}_{\mathrm{rc}}^{\pi}}|p_{\pi}^{t}(s_{\mathrm{rc}}^{\pi}|s_{0})-p_{\pi}^{\star}(s_{\mathrm{rc}}^{\pi})|\
\boldsymbol{F}_{a}(\boldsymbol{\theta},s_{\mathrm{rc}}^{\pi})}_{\text{ignored
in \eqref{equ:fisher_b_sampling}}}\Big{\\}}$
$\displaystyle\quad+\sum_{t=t_{\mathrm{mix}}^{\pi}}^{\infty}\Big{\\{}\sum_{s_{\mathrm{tr}}^{\pi}\in\mathcal{S}_{\mathrm{tr}}^{\pi}}\cancelto{0}{|p_{\pi}^{t}(s_{\mathrm{tr}}^{\pi}|s_{0})-p_{\pi}^{\star}(s_{\mathrm{tr}}^{\pi})|}\
\boldsymbol{F}_{a}(\boldsymbol{\theta},s_{\mathrm{tr}}^{\pi})+\sum_{s_{\mathrm{rc}}^{\pi}\in\mathcal{S}_{\mathrm{rc}}^{\pi}}\cancelto{0}{|p_{\pi}^{t}(s_{\mathrm{rc}}^{\pi}|s_{0})-p_{\pi}^{\star}(s_{\mathrm{rc}}^{\pi})|}\
\boldsymbol{F}_{a}(\boldsymbol{\theta},s_{\mathrm{rc}}^{\pi})\Big{\\}},$ (7)

where we decompose the time summation into pre-absorption, absorption-to-
mixing, and post-mixing, as well as the state summation into transient and
recurrent state sets under $\pi(\boldsymbol{\theta})$, i.e.
$\mathcal{S}=\mathcal{S}_{\mathrm{tr}}^{\pi}\cup\mathcal{S}_{\mathrm{rc}}^{\pi}$.
Note that since $\boldsymbol{F}_{a}(s)$ may have negative off-diagonal
entries, it can _not_ be brought inside the absolute operator to obtain
$\sum_{t=0}^{\infty}\sum_{s\in\mathcal{S}}|p_{\pi}^{t}(s|s_{0})\boldsymbol{F}_{a}(\boldsymbol{\theta},s)-p_{\pi}^{\star}(s)\boldsymbol{F}_{a}(\boldsymbol{\theta},s)|\neq\boldsymbol{F}_{b,a}(\boldsymbol{\theta}|s_{0})$.

After absorption but before mixing (i.e. from $t=t_{\mathrm{abs}}^{\pi|s_{0}}$
to $t_{\mathrm{mix}}^{\pi}-1$), the factor of
$|p_{\pi}^{t}(s_{\mathrm{rc}}^{\pi}|s_{0})-p_{\pi}^{\star}(s_{\mathrm{rc}}^{\pi})|$
is generally non-zero (although, it decreases to $0$ as $t$ approaches
$t_{\mathrm{mix}}^{\pi}$). In general, such an absolute deviation per timestep
is also not a state probability (inheriting the non-stochastic characteristic
of $\boldsymbol{D}_{\pi}$). Therefore, there is no choice but to ignore all
non-zero terms at the aforementioned timesteps in (7), in order to obtain an
expression that enables sampling. That is,

$\tilde{\boldsymbol{F}}_{b,a}^{t_{\mathrm{abs}}^{\pi|s_{0}}}(\boldsymbol{\theta}|s_{0})\coloneqq\Big{\\{}\underbrace{\sum_{t=0}^{t_{\mathrm{abs}}^{\pi|s_{0}}-1}\sum_{s\in\mathcal{S}}p_{\pi}^{t}(s|s_{0})\boldsymbol{F}_{a}(\boldsymbol{\theta},s)}_{\text{$\boldsymbol{F}_{\tau(t_{\mathrm{abs}}^{\pi|s_{0}})}(\boldsymbol{\theta}|s_{0})$
in
\eqref{equ:finite_traj_fisher}}}\Big{\\}}+t_{\mathrm{abs}}^{\pi|s_{0}}\cdot\underbrace{\sum_{s\in\mathcal{S}}p_{\pi}^{\star}(s)\boldsymbol{F}_{a}(\boldsymbol{\theta},s)}_{\text{$\boldsymbol{F}_{\tau(\infty)}(\boldsymbol{\theta})$
in
\eqref{equ:infinite_traj_fisher}}}\approx\underbrace{\boldsymbol{F}_{b,a}(\boldsymbol{\theta}|s_{0})}_{\text{See
\eqref{equ:fisher_b}}}.$ (8)

In both state summations, we do not distinguish transient and recurrent states
notationally, and write the sums over all states $s\in\mathcal{S}$ because
$p_{\pi}^{t}(s_{\mathrm{rc}})=0$ before absorption, and
$p_{\pi}^{\star}(s_{\mathrm{tr}})=0$ at all times.

###  3.1 Alternative interpretations

#### Interpretation 1:

The sampling expression for bias Fisher (8) can also be interpreted as the sum
of Fishers of finite and infinite trajectories (of states and actions). These
Fisher matrices were introduced by bagnell_2003_covps; peters_2003_rlhum who
identified the gain Fisher (4) as the Fisher matrix on the manifold of the
infinite trajectory distribution. That is,
$\boldsymbol{F}_{g,a}(\boldsymbol{\theta})$ is equal to

$\boldsymbol{F}_{\tau(\infty)}(\boldsymbol{\theta})\coloneqq\lim_{\ell\to\infty}\frac{1}{\ell}\underbrace{\sum_{t=0}^{\ell-1}\sum_{s\in\mathcal{S}}p_{\pi}^{t}(s|s_{0};\boldsymbol{\theta})\boldsymbol{F}_{a}(\boldsymbol{\theta},s)}_{\text{$\boldsymbol{F}_{\tau(\ell)}(\boldsymbol{\theta}|s_{0})$
in \eqref{equ:finite_traj_fisher}
below}}=\sum_{s\in\mathcal{S}}\Big{\\{}\underbrace{\lim_{\ell\to\infty}\frac{1}{\ell}\sum_{t=0}^{\ell-1}p_{\pi}^{t}(s|s_{0};\boldsymbol{\theta})}_{p_{\pi}^{\star}(s)}\Big{\\}}\boldsymbol{F}_{a}(\boldsymbol{\theta},s),$
(9)

which involves the trajectory Fisher $\boldsymbol{F}_{\tau(\ell)}$ defined as

$\displaystyle\boldsymbol{F}_{\tau(\ell)}(\boldsymbol{\theta}|s_{0})$
$\displaystyle\coloneqq\mathbb{E}_{\mathrm{T}\sim
p_{\tau(\ell)}}\mathopen{}\mathclose{{}\left[\big{\\{}\nabla\log
p(\mathrm{T};\boldsymbol{\theta})-\mathbb{E}\mathopen{}\mathclose{{}\left[\nabla\log
p(\mathrm{T};\boldsymbol{\theta})}\right]\big{\\}}^{2}}\right]$
$\displaystyle=\mathbb{E}_{\mathrm{T}\sim
p_{\tau(\ell)}}\mathopen{}\mathclose{{}\left[\sum_{t=0}^{\ell-1}\nabla\log\pi(a_{t}|s_{t};\boldsymbol{\theta})\nabla^{\intercal}\log\pi(a_{t}|s_{t};\boldsymbol{\theta})}\right]$
(See below text)
$\displaystyle=\sum_{t=0}^{\ell-1}\sum_{s_{t}\in\mathcal{S}}\sum_{a_{t}\in\mathcal{A}}\pi(a_{t}|s_{t};\boldsymbol{\theta})p_{\pi}^{t}(s_{t}|s_{0})\Big{\\{}\nabla\log\pi(a_{t}|s_{t};\boldsymbol{\theta})\nabla^{\intercal}\log\pi(a_{t}|s_{t};\boldsymbol{\theta})\Big{\\}}.$
(10)

Here,
$\tau_{\ell}(\pi)\coloneqq(s_{0},a_{0},s_{1},a_{1},\ldots,s_{\ell-1},a_{\ell-1}|\pi)$
denotes a trajectory $\tau$ of length $\ell$ under a policy
$\pi(\boldsymbol{\theta})$; note that the last action $a_{\ell-1}$ (at the
last state $s_{\ell-1}$) is included. There exists a distribution, from which
a trajectory is sampled, i.e. $\mathrm{T}\sim p_{\tau(\ell)}$. The penultimate
equality in (10) is because the mean of the score is
$\mathbb{E}\mathopen{}\mathclose{{}\left[\nabla\log
p(\mathrm{T};\boldsymbol{\theta})}\right]=\boldsymbol{0}$, and $\nabla\log
p(\tau_{\ell}(\boldsymbol{\theta}))=\sum_{t=0}^{\ell-1}\nabla\log\pi(a_{t}|s_{t};\boldsymbol{\theta})$
since $\nabla p(s_{0})=\boldsymbol{0}$ obviously, and $\nabla
p(s_{t+1}|s_{t},a_{t})=\boldsymbol{0}$ for randomized policies
$\pi(\boldsymbol{\theta})\in\Pi_{\mathrm{SR}}$.

#### Interpretation 2:

The bias Fisher (8) uses the Fisher of a policy $\pi(\boldsymbol{\theta})$,
that is $\boldsymbol{F}_{a}$ in (4). However, there is another statistical
model, i.e. the state distribution $p_{\pi}^{t}$, that also changes due to the
changes of the policy parameter $\boldsymbol{\theta}$. The bias Fisher,
therefore, can be defined using the Fisher of the joint state-action
distribution
$p_{\pi}^{t}(s,a|s_{0})=p_{\pi}^{t}(s|s_{0})\pi(a|s,\cancel{s_{0}};\boldsymbol{\theta})$.
This follows morimura_2008_nnpg who defined the state-action gain Fisher as
$\boldsymbol{F}_{g,sa}(\boldsymbol{\theta})\coloneqq\boldsymbol{F}_{s}^{\star}(\boldsymbol{\theta})+\boldsymbol{F}_{g,a}(\boldsymbol{\theta})$,
where the stationary state Fisher
$\boldsymbol{F}_{s}^{\star}(\boldsymbol{\theta})\coloneqq\sum_{s\in\mathcal{S}}p_{\pi}^{\star}(s)\nabla\log
p_{\pi}^{\star}(s)\nabla^{\intercal}\log p_{\pi}^{\star}(s)$.

Similarly to the (action) bias Fisher $\tilde{\boldsymbol{F}}_{b,a}$ (8), we
define the state-action bias Fisher as follows,

$\displaystyle\tilde{\boldsymbol{F}}_{b,sa}^{\hat{t}_{\mathrm{abs}}^{\pi|s_{0}}}(\boldsymbol{\theta}|s_{0})\coloneqq\sum_{t=0}^{\hat{t}_{\mathrm{abs}}^{\pi|s_{0}}-1}\sum_{s\in\mathcal{S}}\sum_{a\in\mathcal{A}}p_{\pi}^{t}(s,a|s_{0})\nabla\log
p_{\pi}^{t}(s,a|s_{0})\nabla^{\intercal}\log
p_{\pi}^{t}(s,a|s_{0})+\hat{t}_{\mathrm{abs}}^{\pi|s_{0}}\cdot\boldsymbol{F}_{g,sa}(\boldsymbol{\theta})$
$\displaystyle=\underbrace{\sum_{t=0}^{\hat{t}_{\mathrm{abs}}^{\pi|s_{0}}-1}\boldsymbol{F}_{s}^{t}(\boldsymbol{\theta})+\hat{t}_{\mathrm{abs}}^{\pi|s_{0}}\cdot\boldsymbol{F}_{s}^{\star}(\boldsymbol{\theta})}_{\text{neglected
(set to $\boldsymbol{0}$) in
\eqref{equ:fisher_b_sampling}}}+\underbrace{\sum_{t=0}^{\hat{t}_{\mathrm{abs}}^{\pi|s_{0}}-1}\sum_{s\in\mathcal{S}}p_{\pi}^{t}(s|s_{0})\boldsymbol{F}_{a}(s,\boldsymbol{\theta})+\hat{t}_{\mathrm{abs}}^{\pi|s_{0}}\cdot\boldsymbol{F}_{g,a}(\boldsymbol{\theta})}_{\text{$\tilde{\boldsymbol{F}}_{b,a}^{\hat{t}_{\mathrm{abs}}^{\pi|s_{0}}}(\boldsymbol{\theta}|s_{0})$
in \eqref{equ:fisher_b_sampling}}},$

where
$\boldsymbol{F}_{s}^{t}(\boldsymbol{\theta})\coloneqq\sum_{s\in\mathcal{S}}p_{\pi}^{t}(s|s_{0})\nabla\log
p_{\pi}^{t}(s|s_{0})\nabla^{\intercal}\log p_{\pi}^{t}(s|s_{0})$. As can be
seen, the bias Fisher (8) does not take into account the changes of the state
distribution under a policy $\pi(\boldsymbol{\theta})$ due to perturbations in
$\boldsymbol{\theta}$, i.e. setting $\nabla
p_{\pi}^{t}(s|s_{0})\leftarrow\boldsymbol{0}$ and $\nabla
p_{\pi}^{\star}(s)\leftarrow\boldsymbol{0}$.

###  3.2 Practical considerations

In model-free RL, we estimate the absorption time
$t_{\mathrm{abs}}^{\pi|s_{0}}\approx\hat{t}_{\mathrm{abs}}$, which is held
constant for all policies and all initial states for simplicity. Larger
$\hat{t}_{\mathrm{abs}}$ does not necessarily mean lower approximation error
in (8) because the process becomes more likely to already enter (and will
never escape) the recurrent class. Although larger $\hat{t}_{\mathrm{abs}}$
may compensate for the ignored terms of (7), the compensation may be
excessive. Moreover, it yields larger variance in sampling-based
approximation. In contrast, smaller $\hat{t}_{\mathrm{abs}}$ lowers the sample
variance with the cost of missing a number (i.e.
$t_{\mathrm{abs}}^{\pi|s_{0}}-\hat{t}_{\mathrm{abs}}$ ) of expectation terms
$\mathbb{E}_{S\sim
p_{\pi}^{t}}\mathopen{}\mathclose{{}\left[\boldsymbol{F}_{a}(\boldsymbol{\theta},S)}\right]$
and $\mathbb{E}_{S\sim
p_{\pi}^{\star}}\mathopen{}\mathclose{{}\left[\boldsymbol{F}_{a}(\boldsymbol{\theta},S)}\right]$
in (8).

Assuming that the process begins at a transient state (i.e.
$s_{0}\in\mathcal{S}_{\mathrm{tr}}^{\pi}$ under any policy $\pi$), we
conjecture that setting $\hat{t}_{\mathrm{abs}}\leftarrow 2$ is reasonable
(for conditioning the bias gradient in (2)) in terms of the aforementioned
trade-off across most RL environments (see Appendix LABEL:app:biasopt for
empirical evidence). This corresponds to the bias Fisher (8) that consists of
a 1-step trajectory Fisher and twice the infinite trajectory Fisher. Note that
setting $\hat{t}_{\mathrm{abs}}\leftarrow 1$ (the lowest because $s_{0}$ is a
transient state) diminishes the involvement of the transition probability as
no finite-step trajectory is taken into account in (8).

The additional experiment results (Appendix LABEL:app:biasopt) suggest that
the proposed bias Fisher (8) is beneficial in bias optimization. As
anticipated, the natural bias gradient
($\boldsymbol{C}\leftarrow\tilde{\boldsymbol{F}}_{b,a}^{(2)}$) outperforms the
standard bias gradient ($\boldsymbol{C}\leftarrow\boldsymbol{I}$) in both
exact and stochastic settings. It also often surpasses a Newton method that
uses the modified Hessian of the bias.

##  4 A nearly Blackwell optimal policy gradient algorithm

In this Section, we proposed an nBw-optimal algorithm via formulation as bi-
level optimization:

$\pi_{b}^{*}\in\operatorname*{\arg\\!\max}_{\pi_{g}^{*}}b(\pi_{g}^{*},s_{0}),\quad\text{subject
to}\
\pi_{g}^{*}\in\Pi_{g}^{*}\coloneqq\operatorname*{\arg\\!\max}_{\pi\in\Pi_{\mathrm{S}}}g(\pi),$
(11)

where $b(\pi)\coloneqq v_{b}(\pi)$, $g(\pi)\coloneqq v_{g}(\pi)$, and the
policy parameterization $\pi(\boldsymbol{\theta})$ is notationally hidden.
Here, the lower-level gain optimization yields the feasible set $\Pi_{g}^{*}$
for the upper-level bias optimization. puterman_1994_mdp showed that there
exist a stationary deterministic nBw-optimal policy $\pi_{b}^{*}$, so does a
stationary deterministic gain optimal policy $\pi_{g}^{*}$, for finite MDPs
(as long as the policy parameterization contains such a policy).

Numerically in practice, the lower-level gain optimization in (11) is deemed
converged once the norm of the gain gradient falls below a small positive
$\epsilon$, without considering the negative definiteness condition, i.e.
$\nabla^{2}g(\boldsymbol{\theta})\prec 0$, for a maximum. In addition, such
gain optimization is carried out using a well-established local optimizer so
there is no guarantee of finding a global maximum. These yield a set of
_approximately_ gain optimal policies, i.e.
$\Pi_{g}^{*_{\epsilon}}=\\{\pi(\boldsymbol{\theta}):\|\nabla
g(\boldsymbol{\theta})\|\leq\epsilon\\}$. Thus computationally, we approximate
(11) by imposing a set constraint $\pi_{g}^{*}\in\Pi_{g}^{*_{\epsilon}}$. Note
that one can perform multiple (parallelized) runs of local optimization with
randomized initial points for approximate global optimization in non-concave
objectives (zhigljavsky_1991_tgrs, Ch 2.1).

We solve (11) using a barrier method (with the above approximation), assuming
that $\Pi_{g}^{*_{\epsilon}}$ has an interior and any boundary point can be
attained by approaching it from the interior (for comparison with penalty
methods, see Appendix LABEL:app:barrier_vs_penalty). As an interior-point
technique, the barrier method must be initiated at a strictly feasible point
$\boldsymbol{\theta}$ satisfying $\|\nabla g(\boldsymbol{\theta})\|<\epsilon$.
This is fulfilled by performing gain optimization. Let
$\pi_{g}^{*_{\epsilon}}$ be the resulting gain optimal policy, whose gain is
denoted by $g^{*}_{\epsilon}$. Then, we can construct another set
$\Pi_{g+}^{*_{\epsilon}}=\\{\pi\ |\ (g^{*}_{\epsilon}-\zeta)-g(\pi)\leq 0\\}$,
which can serve the same role as $\Pi_{g}^{*_{\epsilon}}$ (as a constraint
set) and is beneficial because it does not involve the gradients of gain.
Here, $\zeta$ is a positive scalar specifying slackness in $g^{*}_{\epsilon}$
such that $\pi_{g}^{*_{\epsilon}}$ becomes an interior point of
$\Pi_{g+}^{*_{\epsilon}}$. It functions similarly as $\epsilon$ (but recall
$\epsilon$ is used while searching for $\pi_{g}^{*_{\epsilon}}$, whereas
$\zeta$ is used after $\pi_{g}^{*_{\epsilon}}$ is found).

Given a constraint set $\Pi_{g+}^{*_{\epsilon}}$ and its interior point
$\pi_{g}^{*_{\epsilon}}$, we have the following approximation to (11),

$\operatorname*{\arg\\!\max}_{\pi}\big{[}\underbrace{\mathbbm{b}(\pi,s_{0},g^{*}_{\epsilon},\beta,\zeta)}_{\text{Bias-
barrier
function}}\coloneqq\underbrace{b(\pi,s_{0})}_{\text{Bias}}+\underbrace{\beta\log(g(\pi)-g^{*}_{\epsilon}+\zeta)}_{\text{Barrier
function}}\big{]},\quad\underbrace{\text{subject to}\
\pi\in\mathrm{interior}(\Pi_{g+}^{*_{\epsilon}})}_{\text{No need to consider
in practice}},$ (12)

where $\beta$ is a non-negative barrier parameter. Its first and second
derivatives are given by

$\displaystyle\nabla\mathbbm{b}(\boldsymbol{\theta},s_{0})=\nabla
b(\boldsymbol{\theta},s_{0})+\tilde{\beta}\nabla
g(\boldsymbol{\theta}),\quad\text{and}$
$\displaystyle\nabla^{2}\mathbbm{b}(\boldsymbol{\theta},s_{0})=\nabla^{2}b(\boldsymbol{\theta},s_{0})+\tilde{\beta}\nabla^{2}g(\boldsymbol{\theta})-(\tilde{\beta}^{2}/\beta)\nabla
g(\boldsymbol{\theta})\nabla^{\intercal}g(\boldsymbol{\theta}),$ (13)

where $\tilde{\beta}=\beta/(g(\pi)-g^{*}_{\epsilon}+\zeta)$. A preconditioning
matrix for gradient-ascent optimizers (2) can be obtained by substituting the
gain and bias Hessians with their corresponding Fisher matrices, i.e. (4) and
(8), and multiplying the last term of RHS in (13) by $-1$ to make it (always)
positive semidefinite.

We argue setting $\zeta\leftarrow 1$ as follows. For environments where all
policies are gain optimal, we want to suppress the contribution of the barrier
to $\mathbbm{b}$ so that (12) becomes bias-only optimization. This can be
achieved by $\log(g(\pi)-g^{*}_{\epsilon}+(\zeta=1))=0$ since
$g^{*}_{\epsilon}=g(\pi),\forall\pi$ for such environments. For other
environments, a distance of one unit gain (between the interior
$\pi_{g}^{*_{\epsilon}}$ and the boundary) can be justified by observing the
behaviour of a logarithmic function, where
$\log(g(\pi)-g^{*}_{\epsilon}+(\zeta=1))$ drops (without bound) quickly as
$(g(\pi)-g^{*}_{\epsilon})$ approaches $-1$ from its initial value of $0$.
This prevents the optimization iterate from going to a policy parameter that
induces lower gain $g(\pi)<g^{*}_{\epsilon}$.

The barrier method proceeds by solving a sequence of (12) for $k=0,1,\ldots$
with monotonically decreasing barrier parameter $\beta\to 0$ as $k\to\infty$,
where the solution of the previous $k$-th iteration is used as the initial
point for the next. It is expected that as $\beta$ decreases, the solution of
(12) becomes close to that of (11), see boyd_2004_cvxopt. However, the Hessian
(13) tends to be increasingly ill-conditioned as $\beta$ shrinks, demanding a
proper preconditioning matrix (Sec 3).

Algo LABEL:alg:gainbiasbarrier_optim implements the above-mentioned barrier
method to attain nBw optimality. It employs Algo LABEL:alg:gradfisher_sampling
to estimate the gradients and Fisher matrices of gain and bias. These
algorithms altogether can be regarded as approximation to the $n$-discount
optimality policy iteration (puterman_1994_mdp, Ch 10.3).

##  5 Experimental results

(a) $\gamma<\gamma_{\mathrm{Bw}}$

(b) $\gamma\approx\gamma_{\mathrm{Bw}}$

(c) $\gamma_{\mathrm{Bw}}<\gamma\approx 1$

(d) Gain landscape

(e) Bias landscape

(f) Ours

(g) Ours (sampling)

Figure 2: Evaluations of our Algo LABEL:alg:gainbiasbarrier_optim with
$\boldsymbol{\theta}\in\real{2}$ (horizontal and vertical axes) on 6
environments (rows). The landscape plots map low=dark-blue and high=yellow.
The other plots show the absolute difference (low=blue to high=magenta)
between the bias of the nBw-optimal deterministic policy and of the resulting
randomized policy, where the 2D-coordinate indicates the initial value of
$\boldsymbol{\theta}$. The numbers in the middle are the mean and standard
deviation across $41^{2}$ grid values. Here, the lowest absolute bias
difference (blue color) is desirable. For experiment setups, see Appendix
LABEL:app:xprmt_setup.

Our primary experimental aim is to validate and gain insights into the
fundamental mechanisms of our algorithms, rather than demonstrating
performance on a large RL problem. To this end, we consider a small but
sufficient setup of benchmark environments and a policy parameterized as a
single neuron with two weights, i.e.
$\boldsymbol{\theta}=[\theta_{0},\theta_{1}]^{\intercal}\in\real{2}$ (see
Appendix LABEL:app:xprmt_setup for details). This allows us to visualize the
algorithm performance across a large portion of the policy space via
enumeration.

Fig 2 presents our experimental results. First, in exact settings, our
proposed Algo LABEL:alg:gainbiasbarrier_optim compares favourably with the
discounted reward method whose performance is sensitive to the discount factor
$\gamma$. Our proposal is superior to those with poor choices of $\gamma$
(either too low or too high). It is better or at least on par with its
discounted courterpart whose $\gamma$ is close to its critical value
$\gamma_{\mathrm{Bw}}$, which is unknown in practice. The hyperparameter
$\beta_{0}$ of Algo LABEL:alg:gainbiasbarrier_optim is much less sensitive
than $\gamma$, as empirically shown in Appendix LABEL:app:barrierparamtuning.
Second, we compare the exact variant of Algo LABEL:alg:gainbiasbarrier_optim
with its approximation, where we used sampling-based estimates of gradients
and Fisher matrices of gain and bias-barrier function, as in Algo
LABEL:alg:gradfisher_sampling. It can be seen that the approximation results
(with few samples) are reasonably close to the exact’s in most part of the
policy parameter space, following the law of large numbers.

##  6 Conclusions, limitations and future works

We propose a policy gradient method for obtaining nearly Blackwell optimal
policies. It uses an RL-specific barrier, and a novel expression that enables
sampling for natural gradients of the bias.

We have not incorporated any value function approximator to the proposed
algorithm. For larger environments, we need a nearly-Blackwell-optimal actor-
critic RL method. It will accomodate the bias-variance trade-off of natural
gradient estimates of the bias, and enable more efficient 2-timescale step-
wise policy parameter updates. We also envision analysis on the sample
complexity. It will be interesting to investigate for obtaining Blackwell
optimal policies, whether techniques based on $n$-discount optimality require
less samples and computation than those targeting $\gamma$-discounted
optimality (which demand an extra work for finding a proper discount factor).


# Explosive death in coupled oscillators with higher-order interactions

Richita Ghosh  Department of Physics, Central University of Rajasthan,
Rajasthan, Ajmer-305 817, India  Umesh Kumar Verma  Complex Systems Lab,
Department of Physics, Indian Institute of Technology Indore, Khandwa Road,
Simrol, Indore-453 552, India  Sarika Jalan  Complex Systems Lab, Department
of Physics, Indian Institute of Technology Indore, Khandwa Road, Simrol,
Indore-453 552, India  Manish Dev Shrimali<EMAIL_ADDRESS>Department of
Physics, Central University of Rajasthan, Rajasthan, Ajmer-305 817, India

###### Abstract

We investigate the dynamical evolution of globally connected Stuart-Landau
oscillators coupled through conjugate or dis-similar variables on simplicial
complexes. We report a first-order explosive phase transition from oscillatory
state to death state, with $2$-simplex (triadic) interactions, as opposed to
the second-order transition with only $1$-simplex (dyadic) interactions.
Moreover, the system displays four distinct homogeneous steady states in the
presence of triadic interactions, in contrast to the two homogeneous steady
states observed with dyadic interactions. We calculate the backward transition
point analytically, confirming the numerical results and providing the origin
of the dynamical states in the transition region. The study will be useful in
understanding complex systems, such as ecological and epidemiological, having
higher-order interactions and coupling through conjugate variables.

Suggested keywords

###### pacs:

Valid PACS appear here

## Introduction

The dynamics of interacting systems are usually modeled using the concept of
complex networks [1]. The units of each system are denoted as nodes and the
interactions by links joining the nodes form networks. Most of the studies on
networks deal with pairwise interactions between the nodes. However, such a
simplified version of pairwise interactions is not justified while dealing
with systems such as those formed by coauthors of research articles, multiple
species competing for food in an ecosystem, or phenomena of diseases and
rumors in society, etc. In nature, systems such as those governed by social
[2, 3], biological [4, 5], neuroscience [6, 7] or ecological dynamics [8, 9]
are reigned by group interactions. The mathematical framework of pairwise
links in networks is insufficient to represent such systems. Simplicial
complexes are one of the simplest candidates for encoding higher-order
interactions mathematically. A simplex of $N$ dimensions can be conceptualized
as a set of $(n+1)$ nodes. The node set $\mathscr{D}$ may be written as
$\mathscr{D}=[x_{0},x_{1},....x_{n}]$ where $[x_{0},x_{1},....x_{n}]$ is the
set of the vertices (nodes) of the simplex. A one-dimensional simplex
represents a link between two nodes, while a two-dimensional simplex
represents a triangle (denoting interaction between three nodes). Hence,
simplicial complexes may be considered topological constructions constituting
simplices of different dimensions [10].

Further, dynamical processes evolving on complex networks exhibit collective
emergent phenomena such as synchronization [11], self-organization [12],
oscillation quenching [13] and chimera states [14]. Among these, oscillation
suppression in dynamical systems is a well-documented and interesting behavior
arising due to interactions in complex systems. The onset of suppression of
oscillations has been observed through various coupling schemes such as mean-
field diffusive [15], through conjugate [16] and environmental variables [17],
and with attractive-repulsive links [18] in networks of both limit-cycle and
chaotic oscillators.

Suppression of oscillations can occur through two routes; amplitude death (AD)
[19] and oscillation death (OD) [20]. For AD, the coupled oscillators
stabilize to a common steady state which is also a stable global fixed point
of the system whereas in OD, they can settle to different fixed points
referred to as in-homogeneous (IHSS) steady state via symmetry-breaking
bifurcation or to a same fixed point referred as homogeneous (HSS) steady
states. The suppression of oscillations is relevant at various places, for
instance, in medical science in controlling diseases like Parkinson’s and
Alzheimer’s. Other examples include neuron models [21], lasers [22], and
climate systems [23].

Further, the transition from oscillations to the quenching of oscillations may
be either continuous and smooth or discontinuous and explosive. There have
been extensive studies on first-order transition to synchronization referred
to as explosive synchronization [24], explosive percolation [25], and
explosive death [26] on networks with dyadic interactions. Later, explosive
death (ED) [27, 28, 29, 30, 31, 32, 33] and semi-explosive death [34] were
observed in a network of conjugate-coupled limit-cycle oscillators. The
phenomenon of ED brings with it a hysteresis or bi-stable region, i.e., a
region in the parameter space where the oscillatory and the steady states
coexist. Such bi-stability has been observed in a number of physical [35] and
chemical systems [36]. With the advent of studies on higher-order
interactions, collective dynamics, such as synchronization [37, 38] and
chimera [39] have been reported for coupled dynamics on simplicial complexes.
In explosive transitions, a system suddenly experiences a transition to
synchronization from an incoherence state and vice-versa [37]. Further,
multilayer networks of simplicial complexes have exhibited multiple routes to
explosive synchronization [40, 41]. Recently, a sudden transition to anti-
phase synchronization on adaptive simplicial complexes has also been explored
[42].

In this study, we extend the mathematical framework of triadic interactions
[43] to induct conjugate coupling, or coupling through dis-similar variables
[44]. We investigate the effect of such a coupling scheme on limit-cycle
oscillators and discover a first-order phase transition from oscillatory to
steady state in the presence of triadic interactions instead of the continuous
second-order phase transition that we observe with only dyadic interactions.
We perform linear stability analysis for all the fixed points and derive the
condition for the transition point from oscillation to steady state to further
confirm our numerical results.

## Model:

Our model system with $1$-simplex (pairwise or dyadic) and $2$-simplex
(triadic) interactions with conjugate or dis-similar variables has been
mathematically formulated in the following way:

$\displaystyle\dot{\mathbf{X_{i}}}$ $\displaystyle=$
$\displaystyle\mathbf{F}(\mathbf{X_{i}})+\frac{\epsilon_{p}}{N}\sum_{j=1}^{N}C_{ij}\mathbf{f_{p}}(\mathbf{X_{i}},\mathbf{X^{*}_{j}})$
(1)
$\displaystyle+\frac{\epsilon_{h}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}D_{ijk}\mathbf{f_{h}}(\mathbf{X_{i}},\mathbf{X^{*}_{j}},\mathbf{X_{k}})$

where $\mathbf{X^{*}_{i}}$ is the conjugate variable of the state variable
$\mathbf{X_{i}}$. The $w$-dimensional state vector $\mathbf{X_{i}}$ denotes
the dynamics of each unit $i$ where $i=1,2,....N$. The dynamics of each
isolated node are assumed to be identical and described by
$F:\mathbb{R}^{w}\rightarrow\mathbb{R}^{w}$. $f_{p}:\mathbb{R}^{2\times
w}\rightarrow\mathbb{R}^{w}$ and $f_{h}:\mathbb{R}^{3\times
w}\rightarrow\mathbb{R}^{w}$ denote the conjugate coupling functions for the
$1$-simplex and $2$-simplex interactions respectively. $\epsilon_{p}$
illustrates the coupling strength for the pairwise interactions and
$\epsilon_{h}$ for the triadic ones. $C$ is the $N\times N$ adjacency matrix
associated to the $1$-simplex, where $C_{ij}=1$, if $i^{th}$ and $j^{th}$ node
are connected and $C_{ij}=0$, otherwise. Similarly, $D$ is the $N\times
N\times N$ adjacency tensor associated to the $2$-simplex where $D_{ijk}=1$ if
the nodes $i^{th}$, $j^{th}$ and $k^{th}$ are connected and $D_{ijk}=0$,
otherwise. Here the adjacency tensor can be defined as
$D_{ijk}=C_{ij}C_{jk}C_{ki}$.

We consider a globally-connected network of $N$ Stuart-Landau ($SL$) limit
cycle oscillators, interacting with each other through $1$-simplex and
$2$-simplex interactions. The equations governing the dynamics of the system
are:

$\displaystyle\dot{x_{i}}$ $\displaystyle=$
$\displaystyle(1-x_{i}^{2}-y_{i}^{2})x_{i}-wy_{i}+\frac{\epsilon_{p}}{N}\sum_{j=1}^{N}C_{ij}(y_{j}-x_{i})$
$\displaystyle+\frac{\epsilon_{h}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}D_{ijk}(x_{j}y_{k}^{2}-x_{i}^{3})$
$\displaystyle\dot{y_{i}}$ $\displaystyle=$
$\displaystyle(1-x_{i}^{2}-y_{i}^{2})y_{i}+wx_{i}+\frac{\epsilon_{p}}{N}\sum_{j=1}^{N}C_{ij}(x_{j}-y_{i})$
(2)
$\displaystyle+\frac{\epsilon_{h}}{N^{2}}\sum_{j=1}^{N}\sum_{k=1}^{N}D_{ijk}(y_{j}x_{k}^{2}-y_{i}^{3})$

with $i=1,2,...N$. $\omega$ is the intrinsic frequency of the oscillators. The
$N$ nodes in the network are coupled to each other by ‘conjugate’ variables.

The dynamical equations (Eq. 2) have been solved numerically using the $RK4$
method with a step-size of $0.01$ after removing transients of the order of
$10^{5}$. The initial conditions were randomly chosen within the interval
$[-1,1]$.

To characterize the nature of the transition from oscillatory state to death
state in our model system, we adopt an order parameter in terms of average
amplitude [45]

$A(\epsilon)=\frac{a(\epsilon)}{a(0)}$ (3)

where $a(\epsilon)=\sum_{i=1}^{N}\frac{1}{N}(\langle
x_{i,max}\rangle_{t}-\langle x_{i,min}\rangle_{t})$ with
$\langle...\rangle_{t}$ representing the average over long time. $a(\epsilon)$
is a measure of the difference of the global maxima and minima of the time
series data of all the oscillators, computed over a long time after discarding
a sufficient amount of transients at a particular coupling strength
$\epsilon$. $a(0)$ is the average amplitude at $\epsilon=0$. The normalized
order parameter $A(\epsilon)$ is greater than $0$ when the coupled system is
in the oscillatory state. However, when the system attains a steady state,
i.e., $A(\epsilon)=0$. This order parameter has been calculated adiabatic by
varying the value of $\epsilon$ ($\epsilon_{p}$ or $\epsilon_{h}$) in forward
and backward continuation. This process is termed ‘adiabatic’ since the
system’s state in the previous run is used as the initial conditions of the
next run of the simulations [26]. During the forward continuation process, the
coupling strength $\epsilon$ increases slowly in steps of $\delta\epsilon$.
When a specified value of $\epsilon$ is reached, the backward continuation
process starts where $\epsilon$ is now decreased slowly in a similar fashion
and continued until the minimum value of $\epsilon$ is reached.

Figure 1:  (Color online) Amplitude order parameter $A$ (a) for $1$-simplex
interactions and (b) for $2$-simplex interactions. The fraction of oscillators
$N_{d}/N$ is calculated for (c) $1$-simplex interactions and (d) $2$-simplex
interactions. Here $\omega=1$ and $N=100$.

## Numerical Results for 1-simplex (dyadic) interactions:

To investigate the dynamics with only $1$-simplex or pairwise interactions, we
set $\epsilon_{h}$ = 0 in Eq. (2). The coupling strength $\epsilon_{p}$ varies
adiabatically in steps of $0.02$ in both forward and backward continuation.
Fig. 1(a) exhibits a second-order transition from oscillation to death of
oscillations. The order-parameter $A(\epsilon_{p})$ gradually decreases to
zero, indicating suppression of oscillations in the system with a gradual
increase in coupling strength. The forward and backward continuations show no
change in the transition point, i.e., in both cases, the critical value of
$\epsilon_{p}$ remains the same, as observed in Fig. 1(a). Hence, in the
presence of only conjugate dyadic links, the system shows a continuous and
reversible second-order phase transition from oscillation to the death state,
without any hysteresis via supercritical Hopf bifurcation.

Next, we discuss the dynamics of the system at the transition points. The
bifurcation diagram for the $SL$ oscillators with $1$-simplex interactions has
been illustrated in Fig. 2(a). Here, a supercritical Hopf bifurcation destroys
the oscillations from the stable limit cycle, displayed in green ($\gamma$). A
subsequent supercritical pitch-fork bifurcation also occurs immediately,
giving rise to the two stable solutions which have been displayed in red
($\rho$) in Fig. 2(a)).

## Numerical results for 2-simplex (triadic) interactions:

Next, we consider the case of $\epsilon_{p}=0$ in Eq. 2, i.e., only the
triadic interactions are prevalent in the model. $\epsilon_{h}$ is varied in
an adiabatic manner, from $0$ to $\epsilon_{h,max}$ in steps of
$\delta\epsilon_{h}=0.02$ for the forward continuation. During this process,
the magnitude of the order parameter $A(\epsilon_{h})$ suddenly drops to $0$,
indicating an explosive transition from oscillatory to steady state, as shown
in Fig. 1(b). The transition point from oscillatory to death state occurs for
a critical $\epsilon_{h}=5.2$. For the backward continuation process, we
decrease $\epsilon_{h}$ from $\epsilon_{h,max}$ to $0$ similarly and observe
that the amplitude order-parameter jumps from $0$ to non-zero values at around
$\epsilon_{h}=4.82$. Hence, it is clear that the forward and the backward
transition points do not coincide. This indicates a hysteresis region in the
parameter space, i.e., a bi-stable region where both oscillatory and HSS
states co-exist. Such an abrupt and sudden transition with a hysteresis region
can be categorized as a first-order transition.

To investigate whether all the oscillators are going to death state
simultaneously or some oscillators may attain death state faster than others,
we define a parameter $N_{d}/N$ where $N_{d}$ is the number of oscillators
going to steady state, and $N$ is the total number of oscillators. This
parameter is also computed in a similar adiabatic way of forward continuation,
as described previously. Fig. 1(c) and Fig. 1(d) depict the abrupt jump from
$0$ to $1$ for $1$-simplex and $2$-simplex interactions respectively. This
indicates that the oscillators attain a death state simultaneously after
crossing a threshold value of $\epsilon_{h}$. We note that this value of
$\epsilon_{h}$ agrees well with our results for forward continuation.

The dynamics at the transition point for the $2$-simplex interactions differ
from that of the $1$-simplex ones. The oscillations of the system remain in
complete synchronization till a critical value of $\epsilon_{h}$ is reached.
At this point, the system stabilizes to newly created fixed points through
saddle-node bifurcations. Due to the bifurcations, two pairs of stable and
unstable fixed points appear. The oscillations from the stable limit cycle
continue and collide with the unstable fixed points. This leads to the
disappearance of the limit cycle, and the system settles down to any of the
stable fixed points, depending on the initial conditions. This scenario of a
stable limit cycle being destroyed by an unstable fixed point is the mechanism
of a sub-critical Hopf bifurcation. The entire dynamics of the model system at
the transition point are illustrated in Fig. 2(b).

Fig. 3 demonstrates the transition of the model system from oscillation to
oscillation-quenching. The system shows periodic and completely synchronous
oscillations (Fig. 3(a)) for coupling strengths below the critical value of
$\epsilon_{h}$ for which the transition to a steady state takes place. The
system goes to a fixed point as $\epsilon_{h}$ increases beyond this point.
Depending on the choice of initial conditions, the oscillators may settle down
to any of the four homogeneous steady states shown in Fig. 3(b). Fig. 3(b)
exhibits the fixed points $x_{i}^{*}\approx\pm 0.55$, displayed in dark purple
($\Psi$) and $x_{i}^{*}\approx\pm 0.68$, illustrated in dark green ($\Gamma$).
It is interesting to note that the system has a kind of mirror-image symmetry.
That is to say, the dynamics occurring in the phase space above the origin are
mirrored in the phase space below the origin. Fig. 4(a) displays the basin of
attraction at the hysteresis region. The numerical results have been
calculated for $\epsilon_{h}=5.0$ since the hysteresis region stretches from
$\epsilon_{h}\approx 4.82$ to $\epsilon_{h}\approx 5.25$. Here, we confirm the
presence of bi-stability through the two regions illustrated in the figure.
The pink ($\phi$) region corresponds to the initial conditions triggering the
oscillatory state, and the light blue ($\Lambda$) region denotes the ones for
the death state. The initial conditions $x_{i}$ and $y_{i}$ have been varied
from $-1$ to $1$. In Fig. 4(b), the basin of attraction for the four fixed
point attractors in the system has been displayed. The value of the coupling
strength has been fixed at around $5.6$ since we are exploring the basin after
the system achieves a steady state. In this figure, each of the four
differently shaded (colored) regions is populated by the initial conditions
for the four different fixed points that the oscillations settle down to
during the steady state. In Fig. 4, we consider $x=x_{i}$ and $y=y_{i}$ since
the system is in either complete synchronization or HSS in the hysteresis
region.

Figure 2:  (Color online) The bifurcation diagram of the coupled $SL$
oscillators is plotted. (a) for $1$-simplex interactions and (b) for
$2$-simplex interactions. The red ($\rho$) solid lines correspond to the
stable fixed points, the dashed black lines correspond to the unstable ones,
and the green ($\gamma$) solid lines indicate the oscillatory solutions from
the stable limit cycle. $i=1,2.....N$. Figure 3:  (Color online) Time series
of the coupled $SL$ oscillators with $2$-simplex interactions is plotted. (a)
at $\epsilon_{h}=2.45$, i.e., in the oscillatory region (b) at
$\epsilon_{h}=6.11$, i.e., at steady state region, one set of fixed points
($x_{i}\approx\pm 0.55)$ which are displayed in dark purple ($\Psi$) and the
other set of fixed points ($x_{i}\approx\pm 0.68)$ which are displayed in dark
green ($\Gamma$) also appear depending on the choice of initial conditions for
the same value of $\epsilon_{h}=6.11$. The other parameters are $\omega=1$,
$N=100$ and $i=1,2.....N$.

We further explore the impact of the parameter $\omega$ on the dynamics of the
$SL$ oscillators with $2$-simplex interactions. The parameter-space
($\epsilon_{h}-\omega$) diagram Fig. 5(a) illustrates that the transition from
oscillatory to steady state occurs for $\omega\leq 1$ via first-order
transition. It is important to note that here in HSS, the system may stabilize
to any of the four fixed points, as depicted in Fig. 3(b).

## Both 1-simplex and 2-simplex interactions:

We also study a case where $SL$ oscillators interact through $1$-and
$2$-simplex interactions. For a weak $1$-simplex interaction
($\epsilon_{p}\leq 0.8$), the system manifests a first-order transition from
the oscillatory to the steady state with an increase in $\epsilon_{h}$. For
intermittent values, $1<\epsilon_{p}<0.8$, a region of bi-stability or
hysteresis is observed regardless of $\epsilon_{h}$ values. The system with
stronger $\epsilon_{p}\geq 1.0$ always attains a steady state and has no
additional effect of $2$-simplex interaction. Hence, we can conclude that the
2-simplex conjugate interactions lead to more complex dynamics with the first-
order transition from oscillatory to steady state in the globally coupled
Stuart-Landau oscillators as compared to second-order transition in the
presence of $1$-simplex conjugate interaction.

## Analytical Calculations:

Analysing Eq. 2 and putting $\epsilon_{h}=0$ along with
$\epsilon_{p}=\epsilon$ we derive the following stable fixed point:

$\displaystyle x^{*}$
$\displaystyle=\sqrt{\frac{{1-\epsilon+\omega-\frac{\omega}{\epsilon}+\eta}}{2}}$
$\displaystyle y^{*}$
$\displaystyle=\frac{\epsilon(\epsilon-\omega)(\epsilon+\omega)\sqrt{1+\omega+\frac{\epsilon^{2}\eta-\epsilon(\epsilon^{2}+\omega)}{\epsilon^{2}}}}{\sqrt{2}\eta\epsilon^{2}}$
(4)

where
$\eta=\frac{\sqrt{\epsilon^{2}(\epsilon-\omega)^{3}(\epsilon+\omega)}}{\epsilon^{2}}$.
To find the condition of the transition from oscillatory to death state, the
Jacobian $\mathscr{W}$ is written in the form of block circulant matrices
[46], i.e., $\bf{\mathcal{J}}=circ(\zeta_{1},\zeta_{2},...\zeta_{2})$.

$\zeta_{1}=$
$\begin{bmatrix}1-\epsilon-3x^{2}-y^{2}&\frac{\epsilon}{2}-\omega-2xy\\\
\frac{\epsilon}{2}+\omega-2xy&1-\epsilon-x^{2}-3y^{2}\\\ \end{bmatrix}$ and
$\zeta_{2}$ = $\begin{bmatrix}0&\frac{\epsilon}{2}\\\
\frac{\epsilon}{2}&0\end{bmatrix}$ A $2(N-1)$ degeneracy exists for $2(N-1)$
eigenvalues of $\mathscr{W}$. These eigenvalues are equal to the eigenvalues
of the matrix $(\zeta_{1}-\zeta_{2})$. The other eigenvalues of $\mathscr{W}$
will equal to those of the matrix $\zeta_{1}+(N-1)\zeta_{2}$. The
characteristic equation corresponding to the matrix $\zeta_{1}+(N-1)\zeta_{2}$
is given by

$\displaystyle\lambda^{2}+d_{1}\lambda+d_{0}=0$ (5)

where $d_{1}=-2+2\epsilon+4x^{*^{2}}+4y^{*^{2}}$ and
$d_{0}=1-2\epsilon+\omega^{2}-4x^{*^{2}}+4\epsilon
x^{*^{2}}+3x^{*^{4}}+4\epsilon
x^{*}y^{*}-4y^{*^{2}}+6x^{*^{2}}y^{*^{2}}+3y^{*^{4}}$. The condition for the
Hopf bifurcation point is derived according to Routh-Hurwitz criterion [47] by
putting $d_{1}=0$. Hence, the critical value of $\epsilon$ for the transition
of the system from the oscillatory to OD state through Hopf bifurcation is
given by,

$\displaystyle\epsilon_{HB}=\frac{1}{3}\Big{(}-1+2\sqrt{1+3\omega^{2}}\Big{)}$
(6)

The condition for pitchfork bifurcation may also be derived by putting
$d_{0}=0$.

$\displaystyle\epsilon_{PB}=\frac{1}{2}(1+\omega^{2})$ (7)

Putting $\omega=1$ in Eq. 6 and Eq. 7, $\epsilon_{HB}$ and $\epsilon_{PB}$ are
both found to be $1$, which is in exact agreement with our numerical results.

We analyze Eq. 2 again, and this time, by considering only the $2$-simplex
interactions, i.e., $\epsilon_{p}=0$ and $\epsilon_{h}$ is generalized to
$\epsilon$. Numerical simulations reflect that the coupled system stabilizes
at the HSS. Hence the backward transition point for the explosive transition
can be computed by the stability analysis of HSS. The stable fixed points for
this system are, $x_{i}^{*}=x^{*},y_{i}^{*}=y^{*},\forall i=1,...,N$ , where,

$\displaystyle x^{*}$
$\displaystyle=\pm\frac{1}{2}\sqrt{\frac{{2+\epsilon+\alpha\pm\sqrt{2}(1+\epsilon)\beta}}{1+\epsilon}}$
(8) $\displaystyle y^{*}$
$\displaystyle=\pm\frac{1}{4(1+\epsilon)w}(\pm\epsilon\pm\alpha+\sqrt{2}\epsilon\beta+\sqrt{2}\epsilon^{2}\beta)x*$

where $\alpha=\sqrt{\epsilon^{2}-4(1+\epsilon)w^{2}}$ and
$\beta=\sqrt{\frac{{\epsilon+2w^{2}+2\epsilon
w^{2}-\sqrt{\epsilon^{2}-4(1+\epsilon)w^{2}}}}{\epsilon(1+\epsilon)^{2}}}$.
$2N\times 2N$ Jacobian matrices can be written in the form of block circulant
matrices
$\bf{\mathscr{J}=circ(\mathscr{M}_{1},\mathscr{M}_{2},...\mathscr{M}_{2})}$.

$\mathscr{M}_{1}=$ $\begin{bmatrix}1-\Omega
x^{2}+\frac{y^{2}}{3}\kappa&-\omega+\frac{2}{3}\kappa xy\\\
\omega+\frac{2}{3}\kappa xy&1-\Omega y^{2}+\frac{x^{2}}{3}\kappa\\\
\end{bmatrix}$ and $\mathscr{M}_{2}$ = $\begin{bmatrix}\frac{\epsilon
y^{2}}{3}&\frac{2\epsilon xy}{3}\\\ \frac{2\epsilon xy}{3}&\frac{\epsilon
x^{2}}{3}\end{bmatrix}$ with $\Omega=3(1+\epsilon)$ and $\kappa=(-3+\epsilon)$

Figure 4: (Color online) Basin of attraction illustrating the stable states in
the system. (a) At the hysteresis region $\epsilon_{h}=5.0$. The pink $(\phi)$
color represents the initial conditions for the oscillatory state and light
blue ($\Lambda$) for the homogeneous steady state (b) at the death region
$\epsilon_{h}=5.6$. The four distinct colors represent the four different
choices of fixed points for the system at a steady state. Here, $N=100$,
$\omega=1$. $i=1,2.....N$

As explained earlier, there will be $2(N-1)$ eigenvalues of $\mathscr{J}$ with
$2(N-1)$ degeneracy. These eigenvalues are equal to the eigenvalues of the
matrix $(\mathscr{M}_{1}-\mathscr{M}_{2})$. The rest of the eigenvalues will
be equal to the eigenvalues of the matrix
$\mathscr{M}_{1}+(N-1)\mathscr{M}_{2}$ and can be analyzed to know the
stability of the HSS state. The characteristic equation corresponding to this
matrix is given by

$\displaystyle\lambda^{2}+c_{1}\lambda+c_{0}=0$ (9)

where $c_{1}=-2+4x^{*^{2}}+2\epsilon x^{*^{2}}+4y^{*^{2}}+2\epsilon y^{*^{2}}$
and $c_{0}=1+\omega^{2}-4x^{*^{2}}-2\epsilon
x^{*^{2}}+3x^{*}{{}^{4}}-3\epsilon^{2}x^{*^{4}}-4y^{*^{2}}-2\epsilon
y^{*^{2}}+6x^{*^{2}}y^{*^{2}}+24\epsilon
x^{*^{2}}y^{*^{2}}+6\epsilon^{2}x^{*^{2}}y^{*^{2}}+3y^{*^{4}}-3\epsilon^{2}y^{*^{4}}$.
Solving this Eq. 9, we have two eigenvalues:

$\displaystyle\lambda_{1,2}=\frac{1}{2(1+\epsilon)}\Big{(}-2-2\epsilon-\epsilon^{2}-2\alpha-\epsilon\alpha\pm\gamma\Big{)}$
(10)

where
$\gamma=\sqrt{2(2-4\alpha+8\omega^{2}+\epsilon(4+16\omega^{2}-6\alpha)+\epsilon^{2}\xi)}$,
where $\xi=2+6\omega^{2}+\epsilon(\alpha+\epsilon-2\omega^{2})$ and
$\alpha=\sqrt{\epsilon^{2}-4(1+\epsilon)\omega^{2}}$ Equating the real parts
of the eigenvalues to zero, we obtain the expression for the backward
transition point through saddle-node bifurcation,

$\displaystyle\epsilon_{SN}=2(\omega^{2}+\sqrt{\omega^{2}+\omega^{4}})$ (11)

Hence, we confirm that this model system stabilizes to a homogeneous steady
state through a saddle-node bifurcation. The calculated backward transition
point from Eq. 11 after considering $\omega=1$ is found to be 4.82, which
agrees with our numerical results.

Figure 5:  (Color online) Different dynamical domains of $N$ coupled $SL$
oscillators are plotted in the parameter plane (a) $(\epsilon_{h}-\omega)$ for
$\epsilon_{p}=0$ and (b) $(\epsilon_{h}-\epsilon_{p})$ for $\omega=1$. Here
OS, HA, and HSS denote oscillatory state, hysteresis area, and homogeneous
steady state.

## Conclusion

Here, we have investigated the $SL$ oscillators coupled through conjugate
variables and having higher-order interactions. In a number of real-world
systems modeling, coupling through a dissimilar variable is desired. Examples
of such situations include certain electronic circuits [48] and
epidemiological models [49]. Conjugate coupling is also found in laser
experiments, where a photodiode detector monitors the light emitted by each
laser diode through feedback [50] and also in ecological settings, where
predators consume prey from external systems and push the food web to new
steady states [51]. In ecological systems, it is difficult to map the co-
existence in communities with multiple species without resorting to the
concept of higher-order interactions. Cross-predation is a common phenomenon
in ecosystems. Hence, one potential application of conjugate coupling with
higher-order interactions is to study the stabilization of the food web in
ecosystems characterized by frequent cross-predation among multiple species.

We report the first-order transition, an irreversible transition from
oscillation to steady state by considering $2$-simplex conjugate-coupled
interactions instead of the reversible and continuous transition for
$1$-simplex conjugate interactions only. As expected, triadic connections
yield more complex dynamics with the coexistence of oscillatory and steady
states in the transition region. We have observed a thin hysteresis region in
the parameter space. This region contains two stable states; oscillatory and
steady. This bi-stability in the hysteresis region has important applications
in the natural systems where oscillatory and steady states have been found to
coexist, such as fruit flies in eclosion rhythm [52] and injection of stimulus
pulse to restart the activity in the sino-arterial node [53] among others.
After the transition to the steady state, the system stabilizes at four
different fixed points, depending on the choice of the initial conditions. We
analyzed the basins of attraction in the hysteresis and steady-state regions
and confirmed our numerical results analytically by performing linear
stability analysis. The analytical calculations have further allowed us to
derive the backward transition point that matches the numerical results. We
have reported the occurrence of consecutive saddle-node bifurcations at the
transition region, leading the system to settle down to HSS.

It should be interesting to extend the present study for other network
topologies, e.g., Erdös-Réyni random networks and small-world networks, etc.
Further, this study is restricted to only dyadic and triadic interactions. The
inclusion of higher-order simplicial structures is a straightforward extension
of the present work. To conclude, the study may be used to understand the
origin of hysteresis and oscillation suppression phenomena in biological,
ecological, and chemical systems, which have multiple agents interacting with
each other and need to be mapped with higher-order interactions.

## Acknowledgements

RG acknowledges financial support from UGC New Delhi for SJSGC fellowship. SJ
gratefully acknowledges SERB Power grant SPF/2021/000136. MDS acknowledges
financial support from SERB, Department of Science and Technology (DST), India
(Grant No. CRG/2021/003301).

## References

  * Pastor-Satorras  _et al._ [2015] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Epidemic processes in complex networks, Reviews of modern physics 87, 925 (2015).
  * Benson  _et al._ [2016] A. R. Benson, D. F. Gleich, and J. Leskovec, Higher-order organization of complex networks, Science 353, 163 (2016).
  * Centola  _et al._ [2018] D. Centola, J. Becker, D. Brackbill, and A. Baronchelli, Experimental evidence for tipping points in social convention, Science 360, 1116 (2018).
  * Sanchez-Gorostiaga  _et al._ [2019] A. Sanchez-Gorostiaga, D. Bajić, M. L. Osborne, J. F. Poyatos, and A. Sanchez, High-order interactions distort the functional landscape of microbial consortia, PLoS Biology 17, e3000550 (2019).
  * Ritz  _et al._ [2014] A. Ritz, A. N. Tegge, H. Kim, C. L. Poirel, and T. Murali, Signaling hypergraphs, Trends in biotechnology 32, 356 (2014).
  * Petri  _et al._ [2014] G. Petri, P. Expert, F. Turkheimer, R. Carhart-Harris, D. Nutt, P. J. Hellyer, and F. Vaccarino, Homological scaffolds of brain functional networks, Journal of The Royal Society Interface 11, 20140873 (2014).
  * Yu  _et al._ [2011] S. Yu, H. Yang, H. Nakahara, G. S. Santos, D. Nikolić, and D. Plenz, Higher-order interactions characterized in cortical activity, Journal of neuroscience 31, 17514 (2011).
  * Grilli  _et al._ [2017] J. Grilli, G. Barabás, M. J. Michalska-Smith, and S. Allesina, Higher-order interactions stabilize dynamics in competitive network models, Nature 548, 210 (2017).
  * Bairey  _et al._ [2016] E. Bairey, E. D. Kelsic, and R. Kishony, High-order species interactions shape ecosystem diversity, Nature communications 7, 12285 (2016).
  * Boccaletti  _et al._ [2023] S. Boccaletti, P. De Lellis, C. del Genio, K. Alfaro-Bittner, R. Criado, S. Jalan, and M. Romance, The structure and dynamics of networks with higher order interactions, Physics Reports 1018, 1 (2023).
  * Pecora and Carroll [1990] L. M. Pecora and T. L. Carroll, Synchronization in chaotic systems, Physical review letters 64, 821 (1990).
  * Lehn [2002] J.-M. Lehn, Toward self-organization and complex matter, Science 295, 2400 (2002).
  * Koseska  _et al._ [2013] A. Koseska, E. Volkov, and J. Kurths, Oscillation quenching mechanisms: Amplitude vs. oscillation death, Physics Reports 531, 173 (2013).
  * Abrams and Strogatz [2004] D. M. Abrams and S. H. Strogatz, Chimera states for coupled oscillators, Physical review letters 93, 174102 (2004).
  * Verma  _et al._ [2018] U. K. Verma, A. Sharma, N. K. Kamal, and M. D. Shrimali, First order transition to oscillation death through an environment, Physics Letters A 382, 2122 (2018).
  * Karnatak  _et al._ [2009] R. Karnatak, R. Ramaswamy, and A. Prasad, Synchronization regimes in conjugate coupled chaotic oscillators, Chaos: An Interdisciplinary Journal of Nonlinear Science 19, 033143 (2009).
  * Arumugam  _et al._ [2016] R. Arumugam, P. S. Dutta, and T. Banerjee, Environmental coupling in ecosystems: From oscillation quenching to rhythmogenesis, Physical Review E 94, 022206 (2016).
  * Dixit and Shrimali [2020] S. Dixit and M. D. Shrimali, Static and dynamic attractive–repulsive interactions in two coupled nonlinear oscillators, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 033114 (2020).
  * Saxena  _et al._ [2012] G. Saxena, A. Prasad, and R. Ramaswamy, Amplitude death: The emergence of stationarity in coupled nonlinear systems, Physics Reports 521, 205 (2012).
  * Koseska  _et al._ [2010] A. Koseska, E. Volkov, and J. Kurths, Parameter mismatches and oscillation death in coupled oscillators, Chaos: An Interdisciplinary Journal of Nonlinear Science 20, 023132 (2010).
  * Ermentrout [1990] G. B. Ermentrout, Oscillator death in populations of “all to all” coupled nonlinear oscillators, Physica D: Nonlinear Phenomena 41, 219 (1990).
  * Wei and Lun [2007] M.-D. Wei and J.-C. Lun, Amplitude death in coupled chaotic solid-state lasers with cavity-configuration-dependent instabilities, Applied Physics Letters 91, 061121 (2007).
  * Gallego and Cessi [2001] B. Gallego and P. Cessi, Decadal variability of two oceans and an atmosphere, Journal of Climate 14, 2815 (2001).
  * Gómez-Gardenes  _et al._ [2011] J. Gómez-Gardenes, S. Gómez, A. Arenas, and Y. Moreno, Explosive synchronization transitions in scale-free networks, Physical review letters 106, 128701 (2011).
  * Bastas  _et al._ [2014] N. Bastas, P. Giazitzidis, M. Maragakis, and K. Kosmidis, Explosive percolation: Unusual transitions of a simple model, Physica A: Statistical Mechanics and its Applications 407, 54 (2014).
  * Verma  _et al._ [2017] U. K. Verma, A. Sharma, N. K. Kamal, J. Kurths, and M. D. Shrimali, Explosive death induced by mean–field diffusion in identical oscillators, Scientific reports 7, 7936 (2017).
  * Zhao  _et al._ [2017] N. Zhao, Z. Sun, X. Yang, and W. Xu, Restoration of oscillation from conjugate-coupling–induced amplitude death, Europhysics Letters 118, 30005 (2017).
  * Zhao  _et al._ [2018a] N. Zhao, Z. Sun, X. Yang, and W. Xu, Explosive death of conjugate coupled van der pol oscillators on networks, Physical Review E 97, 062203 (2018a).
  * Han  _et al._ [2016] W. Han, H. Cheng, Q. Dai, H. Li, P. Ju, and J. Yang, Amplitude death, oscillation death, wave, and multistability in identical stuart–landau oscillators with conjugate coupling, Communications in Nonlinear Science and Numerical Simulation 39, 73 (2016).
  * Sharma [2019] A. Sharma, Time delay induced partial death patterns with conjugate coupling in relay oscillators, Physics Letters A 383, 1865 (2019).
  * Zhao  _et al._ [2018b] N. Zhao, Z. Sun, and W. Xu, Enhancing coherence via tuning coupling range in nonlocally coupled stuart–landau oscillators, Scientific Reports 8, 8721 (2018b).
  * Ponrasu  _et al._ [2018] K. Ponrasu, K. Sathiyadevi, V. Chandrasekar, and M. Lakshmanan, Conjugate coupling-induced symmetry breaking and quenched oscillations, Europhysics Letters 124, 20007 (2018).
  * Dutta  _et al._ [2023] S. Dutta, O. Alamoudi, Y. S. Vakilna, S. Pati, and S. Jalan, Oscillation quenching in stuart-landau oscillators via dissimilar repulsive coupling, Physical Review Research 5, 013074 (2023).
  * Sun  _et al._ [2021] Z. Sun, S. Liu, and N. Zhao, Explosive and semi-explosive death in coupled oscillators, Chaos, Solitons & Fractals 142, 110514 (2021).
  * Herrero  _et al._ [2000] R. Herrero, M. Figueras, J. Rius, F. Pi, and G. Orriols, Experimental observation of the amplitude death effect in two coupled nonlinear oscillators, Physical review letters 84, 5312 (2000).
  * Bar-Eli and Reuveni [1985] K. Bar-Eli and S. Reuveni, Stable stationary states of coupled chemical oscillators. experimental evidence, The Journal of Physical Chemistry 89, 1329 (1985).
  * Skardal  _et al._ [2021] P. S. Skardal, L. Arola-Fernández, D. Taylor, and A. Arenas, Higher-order interactions can better optimize network synchronization, Physical Review Research 3, 043193 (2021).
  * Parastesh  _et al._ [2022] F. Parastesh, M. Mehrabbeik, K. Rajagopal, S. Jafari, and M. Perc, Synchronization in hindmarsh–rose neurons subject to higher-order interactions, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 013125 (2022).
  * Kundu and Ghosh [2022] S. Kundu and D. Ghosh, Higher-order interactions promote chimera states, Physical Review E 105, L042202 (2022).
  * Millán  _et al._ [2020] A. P. Millán, J. J. Torres, and G. Bianconi, Explosive higher-order kuramoto dynamics on simplicial complexes, Physical Review Letters 124, 218301 (2020).
  * Jalan and Suman [2022] S. Jalan and A. Suman, Multiple first-order transitions in simplicial complexes on multilayer systems, Physical Review E 106, 044304 (2022).
  * Kachhvah and Jalan [2022] A. D. Kachhvah and S. Jalan, First-order route to antiphase clustering in adaptive simplicial complexes, Physical Review E 105, L062203 (2022).
  * Gambuzza  _et al._ [2021] L. V. Gambuzza, F. Di Patti, L. Gallo, S. Lepri, M. Romance, R. Criado, M. Frasca, V. Latora, and S. Boccaletti, Stability of synchronization in simplicial complexes, Nature communications 12, 1255 (2021).
  * Karnatak  _et al._ [2007] R. Karnatak, R. Ramaswamy, and A. Prasad, Amplitude death in the absence of time delays in identical coupled oscillators, Physical Review E 76, 035201 (2007).
  * Resmi  _et al._ [2012] V. Resmi, G. Ambika, R. Amritkar, and G. Rangarajan, Amplitude death in complex networks induced by environment, Physical Review E 85, 046211 (2012).
  * Davis [1979] P. J. Davis,  _Circulant matrices_ , Vol. 120 (Wiley New York, 1979).
  * Liu [1994] W.-M. Liu, Criterion of hopf bifurcations without using eigenvalues, Journal of Mathematical Analysis and Applications 182, 250 (1994).
  * Singla  _et al._ [2011] T. Singla, N. Pawar, and P. Parmananda, Exploring the dynamics of conjugate coupled chua circuits: Simulations and experiments, Physical Review E 83, 026210 (2011).
  * Goleniewski [1996] G. Goleniewski, Modelling cultivar mixtures using seir compartmental models, Biometrical journal 38, 281 (1996).
  * Kumar  _et al._ [2008] P. Kumar, A. Prasad, and R. Ghosh, Stable phase-locking of an external-cavity diode laser subjected to external optical injection, Journal of Physics B: Atomic, Molecular and Optical Physics 41, 135402 (2008).
  * Karnatak  _et al._ [2014] R. Karnatak, R. Ramaswamy, and U. Feudel, Conjugate coupling in ecosystems: Cross-predation stabilizes food webs, Chaos, Solitons & Fractals 68, 48 (2014).
  * Beuter  _et al._ [2003] A. Beuter, L. Glass, M. C. Mackey, and M. S. Titcombe,  _Nonlinear dynamics in physiology and medicine_ (Springer, 2003).
  * Guevara and Jongsma [1992] M. R. Guevara and H. J. Jongsma, Three ways of abolishing automaticity in sinoatrial node: ionic modeling and nonlinear dynamics, American Journal of Physiology-Heart and Circulatory Physiology 262, H1268 (1992).