# Impact of phase lag on synchronization in frustrated Kuramoto model with
higher-order interactions

Sangita Dutta1<EMAIL_ADDRESS>Abhijit Mondal1 Prosenjit Kundu2
<EMAIL_ADDRESS>Pitambar Khanra3<EMAIL_ADDRESS>Pinaki
Pal1<EMAIL_ADDRESS>Chittaranjan Hens4 1Department of
Mathematics, National Institute of Technology, Durgapur 713209, India
2Dhirubhai Ambani Institute of Information and Communication Technology,
Gandhinagar, Gujarat, 382007, India  3Department of Mathematics, State
University of New York at Buffalo, Buffalo, USA  4Center for Computational
Natural Science and Bioinformatics, International Institute of Informational
Technology, Gachibowli, Hyderabad 500032, India

###### Abstract

The study of first order transition (explosive synchronization) in an ensemble
(network) of coupled oscillators has been the topic of paramount interest
among the researchers for more than one decade. Several frameworks have been
proposed to induce explosive synchronization in a network and it has been
reported that phase frustration in a network usually suppresses first order
transition in the presence of pairwise interactions among the oscillators.
However, on the contrary, by considering networks of phase frustrated coupled
oscillators in the presence of higher order interactions (upto $2$-simplexes)
we show here under certain conditions, phase frustration can promote explosive
synchronization in a network. A reduced order model of the network in the
thermodynamic limit is derived using the Ott-Antonsen ansatz to explain this
surprising result. Analytical treatment of the reduced order model including
bifurcation analysis explains the apparent counter intuitive result quite
clearly.

##  I Introduction

Synchronization Strogatz (2004); Pikovsky  _et al._ (2003); Kuramoto (1984) is
a captivating phenomenon observed in natural and artificial systems Wang  _et
al._ (2017); Kumar  _et al._ (2015); Lee  _et al._ (2018); Blasius and Stone
(2000); Pantaleone (2002); Grzybowski  _et al._ (2016); Kuramoto (2003); Buck
(1988); Yamaguchi  _et al._ (2003); Wiesenfeld  _et al._ (1998); Dörfler and
Bullo (2014). A large set of coupled oscillators undergoes a continuous
synchronization transition if the interaction strength among themselves is
gradually increased. To describe such transition, a simple and analytically
tractable phase model was developed by Yoshiki Kuramoto Kuramoto (2003), in
which each oscillator has its own intrinsic frequency ($\omega_{i}$), and
cross-talks to each other via a periodic coupling function representing
pairwise interactions. The equation of motion of each oscillator under all-to-
all coupling configuration is captured by

$\displaystyle\dot{\theta_{i}}=\omega_{i}+K\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}),\hskip
11.38092pti=1,2,\dots,N,$ (1)

where $K$ is the coupling strength. To measure the level of synchronization,
the order parameter $r$ is defined as

$\displaystyle re^{i\psi}=\sum_{j=1}^{N}e^{i\theta_{j}}.$ (2)

The classical Kuramoto model (KM) described above shows second order
(continuous) transition to global synchronization ($r\rightarrow 1$) state as
the coupling strength is gradually increased from a small value. All the
oscillators rotate with the mean of the natural frequencies in a higher
coupling strength. However, different scenario may appear if a phase-lag or
frustration parameter ($\beta$) Sakaguchi and Kuramoto (1986); Omel’Chenko and
Wolfrum (2012); Kundu  _et al._ (2017); Kundu and Pal (2019); Kundu  _et al._
(2018); Brede and Kalloniatis (2016); Lohe (2015); English  _et al._ (2016);
Wiesenfeld  _et al._ (1998); Lohe (2010); Dorfler and Bullo (2012) is
introduced in the pairwise coupling function
($\sin(\theta_{j}-\theta_{i}-\beta)$). This phase-lag parameter shifts the
system’s synchronous frequency from the mean natural frequency. Thus, $\beta$
can be used as a control parameter that can tune the mean frequency to a
desired one Lohe (2015). The model was first introduced by H. Sakaguchi
jointly with Y. Kuramoto and presently known as Sakaguchi-Kuramoto (SK) model.
In SK model, the phase-lag $\beta$ generally opposes the coupled system to
reach to the global synchronous state and consequently, higher coupling
strength compared to KM is required for achieving a desired level of
synchronization. Also, for a particular choice of frequencies, the
synchronization may deviate from the universal continuous transition route
Omel’Chenko and Wolfrum (2012); Omel’chenko and Wolfrum (2013, 2016). Further
it has been shown, that in a complex network setup, the SK model cannot reach
a global synchronous state (called as erosion of synchronization) Skardal  _et
al._ (2015). Even over a critical $\beta$, where the natural frequencies are
correlated with degree, the hysteresis or explosive width can be entirely
vanished, and the level of synchronization would be significantly reduced
Kundu  _et al._ (2017). To overcome such situations, the design of suitable
frequencies Brede and Kalloniatis (2016); Kundu  _et al._ (2018); Skardal  _et
al._ (2014), usage of multiple layers of networks Khanra  _et al._ (2018);
Kundu  _et al._ (2020); Khanra and Pal (2021), and time-dependent coupling
functions Khanra  _et al._ (2020) have been proposed. Note that, models of
Josephson junction arrays Wiesenfeld  _et al._ (1998), power network systems
Dorfler and Bullo (2012), dynamics of mechanical rotors Mertens and Weaver
(2011), quantum networks Lohe (2010) all can be captured by such phase-lag
oscillators under the paradigm of pairwise interactions.

Figure 1: Schematic diagram of a simplicial complex consisting of $0$-simplex
(circles), $1$-simplex (edges) and $2$-simplex (triangles). Figure 2:
Synchronization profile for the all-to-all connected network with different
phase-lag values as obtained from the numerical simulation. Forward and
backward transitions are indicated by up and down arrows respectively. (a)
Numerically computed $r$ is shown as a function of pairwise coupling $K_{1}$
for fixed triadic coupling $K_{2}=8$. As $\beta$ changes from $0$ to higher
values, the transition gradually becomes continuous from discontinuous. (b)
For fixed $K_{1}=2.5$, $r$ is plotted as a function of $K_{2}$ for various
$\beta$. Larger values of lag induces bistable states in the system.

However, in recent studies the importance of incorporating higher order
interactions (HOI) along with its pairwise counterpart has been emphasized
Majhi  _et al._ (2022); Battiston  _et al._ (2020); Alvarez-Rodriguez  _et
al._ (2021); Iacopini  _et al._ (2019); Martignon  _et al._ (1995); Yu  _et
al._ (2011); Giusti  _et al._ (2016); Reimann  _et al._ (2017); Tanaka and
Aoyagi (2011). Specially the works on neuroscience Yu  _et al._ (2011); Giusti
_et al._ (2016); Reimann  _et al._ (2017); Sizemore  _et al._ (2018), ecology
Chatterjee  _et al._ (2022); Grilli  _et al._ (2017); Bairey  _et al._ (2016);
Kareiva (1994) and physics Ghasemi and Kantz (2022) have highlighted the
crucial role of such higher order interactions ( e.g. in a collaboration
network multi authors collaborate Vasilyeva  _et al._ (2021)) in addition to
pairwise one. Simplicial complexes Salnikov  _et al._ (2018); Costa and Farber
(2016) can successfully encode these higher order structures. A $n$-simplex is
formed by $(n+1)$ interacting units, consisting all the $d$ ($<n$) simplexes.
For example, $1$-simplex denotes the pairwise interactions, $2$-simplex
denotes the three way interactions including pairwise one and so on. These
simplices adhere to one another along their sides, and form a simplicial
complex. FIG. 1 represents schematic diagram of a simplicial complex with
higher order interactions. Here nodes (circles) are connected with edges and
triangles (blue). In case of triangles, pairwise connections between the three
nodes are also taken into account. Till now some of the dynamical behaviors
induced by the HOI have been explored including multistability Skardal and
Arenas (2019); Xu  _et al._ (2020), chimera states Zhang  _et al._ (2021),
chaos Bick  _et al._ (2016) etc. Also it has been shown that the presence of
HOI along with pairwise interactions can produce high level of synchronization
in weaker coupling strength Skardal  _et al._ (2021) and may lead to the first
order or discontinuous transition to synchronization Skardal and Arenas
(2019).

Recently, Skardal et al. has reported in Skardal and Arenas (2020) that the
HOI in a simplicial complex can exhibit abrupt transition to synchrony forming
hysteresis, without any degree frequency correlation. These synchronized
states are stable even in negative pairwise coupling for the larger values of
HOI coupling strengths. They analytically derived the evolution equation for
the order parameter to observe the macroscopic dynamics in all-to-all network
and analyze the bifurcations happen in the considered parameter space. In
Rajwani  _et al._ (2023), it has been demonstrated that adaptation to higher
order coupling leads to tired synchronization transition in addition with
second order transition.

In this letter, we consider SK model with higher order interactions (upto
$2$-simplex) under the all-to-all coupling configuration and observe a
significantly different effect of phase frustration compared to the ones
observed in the presence of pairwise interactions only. It is observed that in
the presence of HOI, the increment of phase frustration eliminates first order
transition in the system as the pairwise coupling is varied for a fixed HOI
coupling strength. On the other hand, for a fixed pairwise coupling strength
as the HOI coupling strength is varied, first-order transition is induced in
the system even for very high values of $\beta$.

Along with numerical simulations, here we present an analytical framework
using Ott-Antonsen anastz Ott and Antonsen (2008) to derive a reduced order
model of the system. The reduced order model captures the dynamics of the
entire system quite accurately in the parameter regime of our interest and
explains the origin of different kinds of transitions in the system. It is
observed that the first order transition is associated with the subcritical
pitchfork bifurcation, while, second order transition is associated with
supercritical pitchfork bifurcation. In other words, the bifurcation analysis
of the reduced order equation provides comprehensive understanding of the
incoherent, coherent and bistable states observed in numerical simulations.

##  II Model Description

The dynamics of a phase frustrated undirected simplicial complex of $N$ number
of nodes with global connectivity is governed by the following system of
equations

$\displaystyle\dot{\theta_{i}}=\omega_{i}$ $\displaystyle+$
$\displaystyle\frac{K_{1}}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}-\beta)$
$\displaystyle+$
$\displaystyle\frac{K_{2}}{N^{2}}\sum_{j=1}^{N}\sum_{l=1}^{N}\sin(2\theta_{j}-\theta_{l}-\theta_{i}-\beta),$
$\displaystyle\hskip 99.58464pti=1,2,\dots,N.$

This is the generalization of classical Sakaguchi-Kuramoto model to a model
consisting of higher order interactions in addition to pairwise ones. Here
$\theta_{i}$ represents the phase and $\omega_{i}$ is the natural frequency of
the $i$th oscillator, $\omega$ is drawn from a Lorentzian distribution with
density function
$g(\omega)=\frac{\Delta}{\pi[\Delta^{2}+(\omega-\omega_{0}^{2})]}$ with mean
$\omega_{0}$. $\beta$ denotes the uniform phase lag of the system. $K_{1}$ and
$K_{2}$ are the coupling strengths of pairwise interactions and HOI
respectively.

Two complex order parameters associated with $1$-simplex and $2$-simplex are
defined by

$\displaystyle
z=re^{i\psi_{1}}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}},\hskip
14.22636ptz_{2}=r_{2}e^{i\psi_{2}}=\frac{1}{N}\sum_{j=1}^{N}e^{2i\theta_{j}}$
(4)

where the amplitudes $r=|z|$ and $r_{2}=|z_{2}|$ of the order parameters $z$
and $z_{2}$ respectively measures the level of synchronization, $\psi_{1}$ and
$\psi_{2}$ are the respective arguments, representing the average phases of
the oscillators.

##  III Numerical observation

To find the effect of the phase frustration on SK model with HOI, we simulate
the equation (II) with network size $N=10^{3}$. For the simulation, $4$-th
order Runge-Kutta method with time step $\delta t=0.01$ is used. First we
check the effect of phase lag $\beta$ with fixed $K_{2}$ ($=8$). Here the
values of $\beta$ are varied in the range of $[0,\frac{\pi}{2})$. Initially
the phases of the oscillators are spread uniformly at random around a circle
which means that the initial condition are chosen randomly from $[-\pi,\pi]$.
The natural frequencies are drawn from a Lorenzian distribution with mean $0$
and half width $\Delta=1$. To determine the nature of transitions to
synchronization we continue the simulation of the system (II) by varying the
relevant parameter both in forward and backward directions. We start the
simulation with $K_{1}=-1$ and gradually increase the pairwise coupling
strength $K_{1}$ to achieve the synchronized state for fixed $K_{2}=8$. For
backward transition we start from the synchronized state achieved with the
simulation for forward transition and gradually decrease $K_{1}$ until we
reach to $K_{1}=-1$. Note that, during forward and backward continuation, the
final state of the previous simulation is taken as the initial condition for
the present simulation. FIG. 2(a) shows the variation of $r$ with $K_{1}$ for
different values of $\beta$ namely, $\beta=0$ (blue), $\beta=1$ (magenta),
$\beta=1.2$ (cyan), $\beta=1.318$ (green), and $\beta=1.35$ (red) for both
forward and backward continuations. The figure shows that the system exhibits
explosive synchronization transition for smaller values of $\beta$ and the
transition becomes continuous as $\beta$ is gradually increased from $0$ to
$\frac{\pi}{2}$. Also, it is evident from FIG. 2(a) that the hysteresis width
for the synchronization transition decreases with the increase of $\beta$ and
the transition become second order or continuous for $\beta>1.318$.

Next, to understand the impact of HOI coupling strength $K_{2}$ on the
synchronization transition, we simulate the system (II) by varying $K_{2}$
with fixed $K_{1}=2.5$ for different $\beta$. The variation of the
synchronization order parameter $r$ with respect to $K_{2}$ (both for forward
and backward continuation, and for different $\beta$) is presented in the FIG.
2(b). From the figure one can observe for a low $\beta$ value ( $\beta=0$
(blue) and $0.6435$ (bright green)), $r$ follows a continuous transition, i.e.
forward and backward continuation data follow the same path. On the other
hand, for $\beta>0.6435$, say for $\beta=0.8$ (brown), $1$ (magenta) and $1.2$
(cyan) discontinuous (first order) transition to synchronization is observed.
Therefore, from the simulation result, it is clearly observed that in this
case with the increase of $\beta$ the backward transition points move forward
but explosive synchronization is not suppressed by higher phase frustration of
the system. This unexpected result seems to be a unique feature of HOI. Here
the effect is opposite to the one that is observed in presence of pairwise
interactions only, where phase frustration is known to suppress the
discontinuous transition to synchronization Kundu  _et al._ (2017). To develop
a deeper understanding about the above surprising result we proceed to derive
a reduced order model of the system (II) using Ott-Antonsen ansatz Ott and
Antonsen (2008). The details are described in the next section.

##  IV A reduced order model

As mentioned earlier, in this section we derive a reduced order model to
explain the numerically observed transition to synchronization in the coupled
system under consideration. First to simplify the onward calculation we write
the equation (II) with the help of equation (4) as follows

$\dot{\theta_{i}}=\omega_{i}+\frac{1}{2i}\left[e^{-i(\theta_{i}+\beta)}H-e^{i(\theta_{i}+\beta)}\bar{H}\right],$
(5)

where $H=K_{1}z+K_{2}z_{2}\bar{z}$ and overbar denotes complex conjugate. In
the thermodynamic limit, we assume that the density of the oscillators with
phase $\theta$, frequency $\omega$ at time $t$ is given by
$f(\theta,\omega,t)$. Since the natural frequency is drawn from a distribution
$g(\omega)$ the function $f$ can be expanded in a Fourier series as

$f=\frac{g(\omega)}{2\pi}\left[1+\sum_{n=1}^{\infty}[f_{n}(\omega,t)e^{in\theta}+c.c.]\right],$
(6)

here $c.c.$ denotes the complex conjugate and $f_{n}(\omega,t)$ is the co-
efficient of the $n$-th term of the series. Now following Ott-Antonsen ansatz
Ott and Antonsen (2008) we take $f_{n}=\alpha^{n}$ for some analytic function
$\alpha$. Since the oscillators are conserved in a network, the density
function must satisfy the continuity equation

$\frac{\partial f}{\partial t}+\frac{\partial}{\partial\theta}{(fv)}=0,$ (7)

where $v=\frac{d\theta}{dt}$ is the velocity given by equation (5). Replacing
$f$ by its Fourier series representation (6) in the continuity equation we
obtain

$\dot{\alpha}+i\alpha\omega+\frac{1}{2}\left[H\alpha^{2}e^{-i\beta}-\bar{H}e^{i\beta}\right]=0.$
(8)

Now, in continuous form, the complex order parameter can be written as

$\displaystyle\bar{z}$ $\displaystyle=$ $\displaystyle\int
d\omega\int\bar{f}(\theta,\omega,t)e^{-i\theta}d\theta$ (9) $\displaystyle=$
$\displaystyle\int
d\omega\int_{0}^{2\pi}\frac{g(\omega)}{2\pi}[1+\sum_{n=1}^{\infty}[\bar{\alpha}^{n}(\omega,t)e^{-in\theta}$
$\displaystyle\hskip
71.13188pt+\alpha^{n}(\omega,t)e^{in\theta}]]e^{-i\theta}d\theta$
$\displaystyle=$ $\displaystyle\int\alpha(\omega,t)g(\omega)d\omega.$

Integrating this in the lower half of the complex plane using the Cauchy’s
integral theorem we get, $\bar{z}=\alpha(\omega_{0}-i\Delta,t)$. Similarly we
can derive $\bar{z}_{2}=\alpha^{2}(\omega_{0}-i\Delta,t)={\bar{z}}^{2}$.

At $\omega=\omega_{0}-i\Delta$, equation (8) leads to

$\displaystyle\dot{\bar{z}}=-i\bar{z}(\omega_{0}-i\Delta)+\frac{1}{2}[(K_{1}\bar{z}+K_{2}\bar{z_{2}}z)e^{i\beta}$
$\displaystyle-{\bar{z}}^{2}(K_{1}z+K_{2}z_{2}\bar{z})e^{-i\beta}].$

Taking complex conjugate on both sides we get

$\displaystyle\dot{z}=iz\omega_{0}-z\Delta+\frac{1}{2}[(K_{1}z+K_{2}z^{2}\bar{z})e^{-i\beta}$
$\displaystyle-z^{2}(K_{1}\bar{z}+K_{2}{\bar{z}}^{2}z)e^{i\beta}].$ (10)

Putting $z=re^{i\psi}$ and comparing real and imaginary parts on both sides of
equation (IV) we obtain the following reduced order model for the system (II)

$\displaystyle\dot{r}$ $\displaystyle=$ $\displaystyle-\Delta
r+\frac{1}{2}[(K_{1}r+K_{2}r^{3})\cos\beta$ (11) $\displaystyle-$
$\displaystyle(K_{1}r+K_{2}r^{3})r^{2}\cos\beta],$
$\displaystyle\dot{\psi}=\omega_{0}+\frac{1}{2}[-(K_{1}+K_{2}r^{2})\sin\beta$
$\displaystyle-(K_{1}+K_{2}r^{2})r^{2}\sin\beta].$ (12)

Now using the above reduced order model we proceed to explain the numerically
observed collective dynamics of the considered system.

##  V Analytical VS. Numerical results

We perform bifurcation analysis of the reduced order model (11) using the
software MATCONT Dhooge  _et al._ (2003) to understand the origin of the
numerically observed solutions obtained from the system (II). First, we take
$K_{2}=8$ and construct the bifurcation diagrams for $\beta=1$ and $1.318$ by
varying $K_{1}$. FIG. 3(a) and (b) shows the bifurcation diagrams where
variation of $r$ with $K_{1}$ have been presented along with the stability
information. For $\beta=1$, the trivial $r=0$ solution (solid and dashed black
line) becomes unstable via a subcrcitical pitchfork bifurcation at
$K_{1}=3.6$. An unstable branch originated from there moves backward (dashed
magenta curve) and eventually becomes stable via a saddle-node (SN)
bifurcation at $K_{1}=2.9$, and move forward (solid magenta curve) as a stable
branch. Thus, in the range $2.9\leq K_{1}\leq 3.6$ bistability induced by the
subcritical pitchfork bifurcation is observed and it is responsible for the
appearance of first order transition at $\beta=1$. Further, on top of the
bifurcation diagram, some data points obtained from the numerical simulation
(magenta and black circles) are also plotted and a very good agreement is
observed. On the other hand, for $\beta=1.318$, $r=0$ branch (solid and black
line) becomes unstable via supercritical pitchfork bifurcation $K_{1}=8$ and a
stable branch is originated (solid green branch). As a result, a continuous
transition occurs for $\beta=1.318$. In this case also, numerical simulation
data (green and black circles) closely matches with the analytical ones.

Next with fixed $K_{1}=2.5$ we construct two bifurcation diagrams for
$\beta=0$ and $0.8$ by varying $K_{2}$ (FIG. 3(c) and (d)). For $\beta=0$
(solid blue curve), we observe continuous transition and the trivial $r=0$
branch (dashed black line) is always unstable in this case. On the other hand,
for $\beta=0.8$, zero solution is always stable (solid black line). In
addition to that, another pair of non-zero fixed points of the reduced system
originated from a SN bifurcation point at $K_{2}=5.3$, of which one is stable
(solid brown curve) and other one is unstable (dashed brown curve). As a
result, bistability occurs for all $K_{2}\geq 5.3.$ Here, the SN bifurcation
is responsible for the discontinuous transition to synchronization. In this
case also along with the analytical curve we plot the data obtained from the
numerical simulation (brown and black circles) and we observe a very close
agreement. It is now apparent that the reduced order model nicely mimics the
dynamics of the whole system. Therefore we proceed for more detailed analysis
using the model.

It is important to note here that for fixed $K_{2}$ as $\beta$ is increased
starting from $0$, the system moves from discontinuous to continuous
transition states with the variation of $K_{1}$. The discontinuous transition
here is associated with the existence of the SN bifurcation point and it
ceased to exist as the SN point meets the trivial ($r=0$) fixed point. This
observation helps to determine the starting point of continuous transition.

The equilibrium points of the reduced order model (11) are given by

$r=0~{}\mathrm{and}~{}r=\sqrt{\frac{K_{2}\cos\beta-
K_{1}\cos\beta\pm\sqrt{(K_{1}\cos\beta+K_{2}\cos\beta)^{2}-8K_{2}\cos\beta}}{2K_{2}\cos\beta}}.$
(13)

Figure 3: Bifurcation diagrams prepared using the reduced order model. First
column [(a)-(b)]: $r$ as a function of $K_{1}$ for $K_{2}=8$ with $\beta=1$
and $1.318$; Second column [(c)-(d)]: the variation of $r$ with $K_{2}$ for
$K_{1}=2.5$ with $\beta=0$ and $0.8$ respectively. Solid and dashed lines are
obtained from the ROM, and respectively represent the stable and unstable
equilibria. The filled black circles show the saddle-node (SN) bifurcation
points. The colored circles show the values of $r$ corresponding to the stable
solutions as obtained in the numerical simulation.

Out of two nonzero equilibrium points given in (13), the one obtained with the
positive sign (say $r^{+}$) is stable and the other one (say $r^{-}$) is
unstable in the entire ranges of their existence in the parameter space. Now
the SN bifurcation will occur for the parameter values for which $r^{+}=r^{-}$
and it leads to the relation

$(K_{1}+K_{2})^{2}\cos\beta=8K_{2}$ (14)

involving the parameters $K_{1}$, $K_{2}$ and $\beta$. Thus, for any two given
parameters, a feasible solution of the above equation for the third one will
provide the SN bifurcation point. Now when the saddle node point merges with
the fixed point $r=0$, the hysteresis will vanish and the system will follow a
continuous transition path. Imposing the condition (14) and $r=0$ in the
equation (13) we can find the point marking the onset of continuous
transition, which gives

$K_{2}\cos\beta-K_{1}\cos\beta=0.$ (15)

Solving the equation (14) and (15), the condition for the onset of continuous
transition is obtained as

$K_{1}=K_{2}=\frac{2}{\cos\beta}.$ (16)

Now in particular, for fixed $K_{2}=8$, the onset of continuous transition
occurs for $K_{1}=8$ and $\beta=\cos^{-1}(1/4)=1.3181.$ In the numerical
simulation also, the continuous transition start to occurs exactly at this
analytically predicted value.

On the other hand, for fixed $K_{1}$, as $\beta$ is increased, completely
opposite scenario occurs, i.e. in this case, the the system moves from
continuous to discontinuous transition states with the variation of $K_{2}$.
The trivial fixed point $r=0$ is unstable or stable according as
$\beta<~{}\mathrm{or}>~{}\cos^{-1}\left(\frac{2\Delta}{K_{1}}\right)$. In
particular, for $\Delta=1$ and $K_{1}=2.5$,
$\cos^{-1}\left(\frac{2\Delta}{K_{1}}\right)=0.6435$. Thus, for
$\beta<0.6435$, the trivial equilibrium point is unstable and it is stable for
$\beta>0.6435$. It can be easily checked that for $\Delta=1$, $K_{1}=2.5$, and
$\beta<0.6435$ only one nontrivial equilibrium point i.e. $r^{+}$ will exist
and stable. As a result, only second order transition to synchronization is
expected to occur in this parameter regime and it is clearly visible from the
numerically obtained blue and bright green curves presented in the FIG. 2(b).
Whenever, $\beta>0.6435$ along with the trivial stable fixed point, both the
nontrivial fixed points $r^{+}$ (stable) and $r^{-}$ (unstable) will exist
which marks the existence of bistability in the system that is responsible for
first order transition. The analytical result matches very closely with the
simulation results presented in the FIG. 2(b).

Figure 4: Two parameter stability diagrams prepared from the ROM showing the
incoherent (blue), coherent (cyan to maroon) and bistable (grey) regimes.
Stability diagrams: (a) The solid magenta (backward transition or SN points)
and black curves (forward transition point) delimit different synchronization
regimes on the $K_{1}-\beta$ plane for $K_{2}=8.$ The magenta and the black
curves intersect at the point $(8,1.318)$ marking the onset of continuous
transition. (b) Shows the synchronization regimes on the $K_{2}-\beta$ plane
for $K_{1}=2.5$. Below the critical lag (black line) synchronization level is
high and above this value incoherent and bistable regime are separated by the
saddle node curve (solid magenta curve). (c) Different regimes on the
$K_{1}-K_{2}$ plane for $\beta=0.8$. The bistable regime starts from the
intersection point $(2.87,2.87)$ of saddle node and pitchfork bifurcation
curves (supercritical for $K_{2}<2.87$, subcritical for $K_{2}>2.87$) and
expands with increasing $K_{2}$. Black and magenta dots respectively represent
the corresponding transition points obtained from the numerical simulations in
all the cases. The black filled squares in (b) and (c) respectively represent
the numerically obtained critical $\beta$ and $K_{1}$ for the onset of
discontinuous transitions.

Inspired by the close agreement of the reduced order model (ROM) and the
numerical simulations results, using the ROM we construct stability diagrams
on different projections of the parameter space and present in the FIG. 4.
FIG. 4(a) shows the bifurcations in the $K_{1}$ \- $\beta$ plane for fixed
$K_{2}=8$. The magenta curve indicates the backward transition points or the
points of SN bifurcation (given by equation (14))and the black curve indicates
the forward transition points (given by equation (16)). The region bounded by
these two curves depicts the bistable region. It is observed that the area of
bistability reduces with the increase of $K_{1}$. Furthermore, when these two
curves meet, the bistability vanishes at the critical point $(8,1.318)$ and
the system exhibits synchronization via continuous paths (see FIG. 2(a)) only.
At these forward transition points for different lag, the system shows
pitchfork bifurcation. When $\beta$ is greater than (less than) the critical
value $1.318$, the supercritical (subcritical) pitchfork bifurcation is
observed. Here we use colorbar to indicate the values of $r$ calculated using
the equation(13) considering the stable one. The blue side is for incoherent
state, as the color changes from blue to maroon the value of the order
parameter $r$ increases from $0$ to $1$. The bistable region is shaded by grey
colour. Next on top of the two parameter diagram (FIG. 4(a)) we put the
forward (black dots) and backward transition points (magenta dots) obtained
from numerical simulations for a comparison. For numerical simulations we use
$N=1000$ node and use the similar setup described in section III. The
numerically computed points are found to lie precisely on the analytical
curve. These findings demonstrate that the analytical theory correctly
captures the system’s dynamics.

Moving on to the stability diagram on the $K_{2}$ \- $\beta$ plane for
$K_{1}=2.5$ (FIG. 4(b)), it is found that the region is divided into three
different parts, namely, incoherent (blue), coherent (maroon) and bistable
(grey). Note that for $K_{1}=2.5$, the onset of SN bifurcation point as
calculated from the condition (16) takes place at the critical phase lag
$\beta=0.643501$ and $K_{2}=2.5$. For $\beta$ less than this critical value,
the system follows continuous path and the synchronization level is also high
(seen in FIG. 2(b)). As $\beta$ is increased a little more, the system
exhibits bistability. The numerically computed transition points (magenta
dots) shown in the figure show very close match with the analytical one.

Finally, in FIG. 4(c) we present the stability diagram on the $K_{1}$ \-
$K_{2}$ plane for $\beta=0.8$. Here also the region is parted into coherent,
incoherent and bistable states as obtained in the previous case. The bistable
regime starts at the intersection point of the pitchfork line and the saddle
node curve and the area of this regime broadened with increasing $K_{2}$. From
FIG. 4(b) and (c) we can conclude that $K_{2}$ contributes significantly to
promote bistability along with providing high level of synchronization in the
system. It has been verified that for other choices of the fixed values of one
of the parameters, similar stability diagrams are obtained. Therefore, in
summary, it is observed that $K_{2}$ promotes bistability in the system.
Moreover, using the results of the ROM, we can appropriately tune the
parameters in such a way that $\beta$ will promote first order transition in
the system and it depends on the existence of two stable fixed points (one
trivial and another nontrivial) separated by a nontrivial unstable fixed
point.

##  VI Conclusion

In this letter, the results of our investigation on the effect of phase-lag on
the transition to synchronization in all-to-all connected networks of phase
oscillators in the presence of higher order interactions (upto $2$-simplex)
have been presented. For the investigation, along with the detailed
performance of numerical simulations of the full system, a reduced order model
has been derived based on the Ott-Antonsen anastz to understand the origin of
different solutions in numerical simulations.

The numerical simulations with Lorenzian natural frequency reveals two
distinct effects of phase lag: (1) For fixed $K_{2}$, the system exhibits
discontinuous transition to synchronization with accompanying hysteresis loop
in presence of low $\beta$ including $0$ value. As $\beta$ is gradually
increased, the width of the hysteresis loop gradually decreases and
eventually, the system exhibits continuous transition to synchronization when
$\beta$ crosses a critical value. The appearance of discontinuous transition
in globally coupled phase oscillators in presence of HOI is consistent with
the earlier report Skardal and Arenas (2020). Further, the migration of the
system from the discontinuous to continuous transition states with the
increment of $\beta$ in this case occurs due the enhancement of the importance
of pairwise coupling in the network which is typically observed in a network
in the presence of only pairwise couplings Kundu  _et al._ (2017); Kundu and
Pal (2019). (2) Surprisingly, completely opposite scenario is observed in
numerical simulations as the higher order coupling parameter ($K_{2}$) is
varied for a fixed $K_{1}$. For lower $\beta$, transition is found to be
continuous and system exhibits discontinuous transition as $\beta$ is
increased beyond a critical value. Therefore, contrary to negative role of
phase frustration on the discontinuous transition observed in networks with
pairwise coupling, using the proposed framework with HOI, discontinuous
transitions can be achieved in a network with very high phase frustration.

To understand the origin of the counter intuitive numerical simulation results
mentioned above, a reduced order model is derived in the thermodynamic limit
which faithfully mimic the dynamics of the original system. The bifurcation
analysis of the ROM reveals a complex dependence of the transition scenario on
the parameters of the system. It is observed that the first order transition
is always associated with a SN bifurcation. In the first case mentioned above,
the SN responsible for discontinuous transition is associated with a
subcritical pitchfork bifurcation, while in the second case, SN bifurcation
independently occur in a region of the parameter space. Further, all the
analytically derived transition points determined from the ROM show a very
good agreement with the ones determined from the numerical simulation.

## References

  * Strogatz (2004) S. Strogatz,  _Sync: The emerging science of spontaneous order_ (Penguin UK, 2004).
  * Pikovsky  _et al._ (2003) A. Pikovsky, M. Rosenblum, J. Kurths, and J. Kurths,  _Synchronization: a universal concept in nonlinear sciences_ , Vol. 12 (Cambridge University Press, 2003).
  * Kuramoto (1984) Y. Kuramoto, in  _Chemical oscillations, waves, and turbulence_ (Springer, 1984) pp. 111–140.
  * Wang  _et al._ (2017) C.-Q. Wang, A. Pumir, N. B. Garnier, and Z.-H. Liu, Frontiers of Physics 12, 1 (2017).
  * Kumar  _et al._ (2015) P. Kumar, D. K. Verma, P. Parmananda, and S. Boccaletti, Physical Review E 91, 062909 (2015).
  * Lee  _et al._ (2018) U. Lee, M. Kim, K. Lee, C. M. Kaplan, D. J. Clauw, S. Kim, G. A. Mashour, and R. E. Harris, Scientific reports 8, 243 (2018).
  * Blasius and Stone (2000) B. Blasius and L. Stone, International Journal of Bifurcation and Chaos 10, 2361 (2000).
  * Pantaleone (2002) J. Pantaleone, American Journal of Physics 70, 992 (2002).
  * Grzybowski  _et al._ (2016) J. Grzybowski, E. Macau, and T. Yoneyama, Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 113113 (2016).
  * Kuramoto (2003) Y. Kuramoto,  _Chemical oscillations, waves, and turbulence_ (Courier Corporation, 2003).
  * Buck (1988) J. Buck, The Quarterly review of biology 63, 265 (1988).
  * Yamaguchi  _et al._ (2003) S. Yamaguchi, H. Isejima, T. Matsuo, R. Okura, K. Yagita, M. Kobayashi, and H. Okamura, Science 302, 1408 (2003).
  * Wiesenfeld  _et al._ (1998) K. Wiesenfeld, P. Colet, and S. H. Strogatz, Physical Review E 57, 1563 (1998).
  * Dörfler and Bullo (2014) F. Dörfler and F. Bullo, Automatica 50, 1539 (2014).
  * Sakaguchi and Kuramoto (1986) H. Sakaguchi and Y. Kuramoto, Progress of Theoretical Physics 76, 576 (1986).
  * Omel’Chenko and Wolfrum (2012) E. Omel’Chenko and M. Wolfrum, Physical Review Letters 109, 164101 (2012).
  * Kundu  _et al._ (2017) P. Kundu, P. Khanra, C. Hens, and P. Pal, Physical Review E 96, 052216 (2017).
  * Kundu and Pal (2019) P. Kundu and P. Pal, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 013123 (2019).
  * Kundu  _et al._ (2018) P. Kundu, C. Hens, B. Barzel, and P. Pal, EPL (Europhysics Letters) 120, 40002 (2018).
  * Brede and Kalloniatis (2016) M. Brede and A. C. Kalloniatis, Physical Review E 93, 062315 (2016).
  * Lohe (2015) M. Lohe, Automatica 54, 114 (2015).
  * English  _et al._ (2016) L. Q. English, D. Mertens, S. Abdoulkary, C. B. Fritz, K. Skowronski, and P. Kevrekidis, Physical Review E 94, 062212 (2016).
  * Lohe (2010) M. Lohe, Journal of Physics A: Mathematical and Theoretical 43, 465301 (2010).
  * Dorfler and Bullo (2012) F. Dorfler and F. Bullo, SIAM Journal on Control and Optimization 50, 1616 (2012).
  * Omel’chenko and Wolfrum (2013) E. Omel’chenko and M. Wolfrum, Physica D: Nonlinear Phenomena 263, 74 (2013).
  * Omel’chenko and Wolfrum (2016) O. E. Omel’chenko and M. Wolfrum, Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 094806 (2016).
  * Skardal  _et al._ (2015) P. S. Skardal, D. Taylor, J. Sun, and A. Arenas, Physical Review E 91, 010802 (2015).
  * Skardal  _et al._ (2014) P. S. Skardal, D. Taylor, and J. Sun, Physical Review Letters 113, 144101 (2014).
  * Khanra  _et al._ (2018) P. Khanra, P. Kundu, C. Hens, and P. Pal, Physical Review E 98, 052315 (2018).
  * Kundu  _et al._ (2020) P. Kundu, P. Khanra, C. Hens, and P. Pal, EPL (Europhysics Letters) 129, 30004 (2020).
  * Khanra and Pal (2021) P. Khanra and P. Pal, Chaos, Solitons & Fractals 143, 110621 (2021).
  * Khanra  _et al._ (2020) P. Khanra, P. Kundu, P. Pal, P. Ji, and C. Hens, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 031101 (2020).
  * Mertens and Weaver (2011) D. Mertens and R. Weaver, Physical Review E 83, 046221 (2011).
  * Majhi  _et al._ (2022) S. Majhi, M. Perc, and D. Ghosh, Journal of the Royal Society Interface 19, 20220043 (2022).
  * Battiston  _et al._ (2020) F. Battiston, G. Cencetti, I. Iacopini, V. Latora, M. Lucas, A. Patania, J.-G. Young, and G. Petri, Physics Reports 874, 1 (2020).
  * Alvarez-Rodriguez  _et al._ (2021) U. Alvarez-Rodriguez, F. Battiston, G. F. de Arruda, Y. Moreno, M. Perc, and V. Latora, Nature Human Behaviour 5, 586 (2021).
  * Iacopini  _et al._ (2019) I. Iacopini, G. Petri, A. Barrat, and V. Latora, Nature Communications 10, 1 (2019).
  * Martignon  _et al._ (1995) L. Martignon, H. Von Hassein, S. Grün, A. Aertsen, and G. Palm, Biological Cybernetics 73, 69 (1995).
  * Yu  _et al._ (2011) S. Yu, H. Yang, H. Nakahara, G. S. Santos, D. Nikolić, and D. Plenz, Journal of Neuroscience 31, 17514 (2011).
  * Giusti  _et al._ (2016) C. Giusti, R. Ghrist, and D. S. Bassett, Journal of Computational Neuroscience 41, 1 (2016).
  * Reimann  _et al._ (2017) M. W. Reimann, M. Nolte, M. Scolamiero, K. Turner, R. Perin, G. Chindemi, P. Dłotko, R. Levi, K. Hess, and H. Markram, Frontiers in Computational Neuroscience 11, 48 (2017).
  * Tanaka and Aoyagi (2011) T. Tanaka and T. Aoyagi, Physical Review Letters 106, 224101 (2011).
  * Sizemore  _et al._ (2018) A. E. Sizemore, C. Giusti, A. Kahn, J. M. Vettel, R. F. Betzel, and D. S. Bassett, Journal of Computational Neuroscience 44, 115 (2018).
  * Chatterjee  _et al._ (2022) S. Chatterjee, S. Nag Chowdhury, D. Ghosh, and C. Hens, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 103122 (2022).
  * Grilli  _et al._ (2017) J. Grilli, G. Barabás, M. J. Michalska-Smith, and S. Allesina, Nature 548, 210 (2017).
  * Bairey  _et al._ (2016) E. Bairey, E. D. Kelsic, and R. Kishony, Nature Communications 7, 1 (2016).
  * Kareiva (1994) P. Kareiva, Ecology 75, 1527–28 (1994).
  * Ghasemi and Kantz (2022) A. Ghasemi and H. Kantz, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 073101 (2022).
  * Vasilyeva  _et al._ (2021) E. Vasilyeva, A. Kozlov, K. Alfaro-Bittner, D. Musatov, A. Raigorodskii, M. Perc, and S. Boccaletti, Scientific Reports 11, 5666 (2021).
  * Salnikov  _et al._ (2018) V. Salnikov, D. Cassese, and R. Lambiotte, European Journal of Physics 40, 014001 (2018).
  * Costa and Farber (2016) A. Costa and M. Farber, in  _Configuration Spaces: Geometry, Topology and Representation Theory_ (Springer, 2016) pp. 129–153.
  * Skardal and Arenas (2019) P. S. Skardal and A. Arenas, Physical Review Letters 122, 248301 (2019).
  * Xu  _et al._ (2020) C. Xu, X. Wang, and P. S. Skardal, Physical Review Research 2, 023281 (2020).
  * Zhang  _et al._ (2021) Y. Zhang, V. Latora, and A. E. Motter, Communications Physics 4, 195 (2021).
  * Bick  _et al._ (2016) C. Bick, P. Ashwin, and A. Rodrigues, Chaos: An Interdisciplinary Journal of Nonlinear Science 26, 094814 (2016).
  * Skardal  _et al._ (2021) P. S. Skardal, L. Arola-Fernández, D. Taylor, and A. Arenas, Physical Review Research 3, 043193 (2021).
  * Skardal and Arenas (2020) P. S. Skardal and A. Arenas, Communications Physics 3, 1 (2020).
  * Rajwani  _et al._ (2023) P. Rajwani, A. Suman, and S. Jalan, arXiv preprint arXiv:2302.12076 (2023).
  * Ott and Antonsen (2008) E. Ott and T. M. Antonsen, Chaos: An Interdisciplinary Journal of Nonlinear Science 18, 037113 (2008).
  * Dhooge  _et al._ (2003) A. Dhooge, W. Govaerts, and Y. A. Kuznetsov, ACM Transactions on Mathematical Software (TOMS) 29, 141 (2003).


11institutetext: Institute for Neural Computation, Ruhr-University Bochum,
Bochum, Germany
11email<EMAIL_ADDRESS>
https://www.ini.rub.de/

# Modular Networks Prevent Catastrophic Interference in Model-Based Multi-Task
Reinforcement Learning

Robin Schiewer 11 0000-0002-4021-2476 Laurenz Wiskott 11 0000-0001-6237-740X

###### Abstract

In a multi-task reinforcement learning setting, the learner commonly benefits
from training on multiple related tasks by exploiting similarities among them.
At the same time, the trained agent is able to solve a wider range of
different problems. While this effect is well documented for model-free multi-
task methods, we demonstrate a detrimental effect when using a single learned
dynamics model for multiple tasks. Thus, we address the fundamental question
of whether model-based multi-task reinforcement learning benefits from shared
dynamics models in a similar way model-free methods do from shared policy
networks. Using a single dynamics model, we see clear evidence of task
confusion and reduced performance. As a remedy, enforcing an internal
structure for the learned dynamics model by training isolated sub-networks for
each task notably improves performance while using the same amount of
parameters. We illustrate our findings by comparing both methods on a simple
gridworld and a more complex vizdoom multi-task experiment.

###### Keywords:

model-based reinforcement learning multi-task reinforcement learning latent
space models catastrophic interference task confusion.

##  1 Introduction

In recent years, deep reinforcement learning (RL) has shown impressive results
in problem domains such as robotics and game playing [20, 19, 15, 27, 12].
However, sample inefficiency is still a major shortcoming of many of the
methods. To achieve superhuman performance e.g. in video games, the required
number of interactions for complex tasks lies in the tenths of millions.
Model-based approaches mitigate this problem by integrating the collected
sample information into a coherent model of the environmental dynamics [30].
Those learned models are used for direct policy learning [11], planning [13,
7] or to augment existing model-free approaches [8, 24, 21]. In the presence
of a readily available analytical model of the environment, that does not have
to be learned, perfect planning can be used to learn a policy with great
success [31, 29].

Despite the success of deep RL across a wide range of problem domains, most
state of the art approaches have to be re-trained for every new task. Since
similar problems might very well share a common underlying structure, this is
a potential waste of resources [3]. Furthermore, an agent can also encounter a
different mixture of tasks in a single environment. If the agent is created
without the concept of distinct tasks in mind, an alteration of the task
mixture may negatively impact the agent’s performance. While greater sample
efficiency is one reason to engage in multi-task learning, it also produces a
more capable agent that can solve a wider range of problems. Whereas there
exists a vast amount of research around multi-task RL in general [26, 9, 10,
6], the combination of model-based latent space RL and multiple tasks is still
largely unexplored. Since the former has brought impressive improvements for
sample efficiency and the latter promises greater flexibility and reusability,
we combine latent space models with multi-task RL in this paper and
investigate its usefulness.

We focus on model-based multi-task deep RL for complex observation domains
like images. The main question we address is if training the same dynamics
model on multiple similar tasks helps performance through knowledge transfer
or if catastrophic interference outweighs the benefits. If so, how can a
model-based multi-task agent be structured to avoid catastrophic interference?
Our contributions in answering these questions are the following:

  * •

In our experiments we show a detrimental effect on performance through
catastrophic interference when training a single dynamics model on multiple
tasks simultaneously. To mitigate the effect, we propose a world model that
uses multiple distinct latent space dynamics models which are activated
through a context-based task classification network.

  * •

By strict separation of the dynamics networks, the probability that dynamics
of different tasks interfere with each other is minimized. As a side effect,
retraining individual tasks and recombining learned tasks is straightforward.

  * •

We demonstrate the performance difference between a single dynamics model and
our method in a 3-task gridworld and a more complex 2-task 3D environment. To
evaluate both approaches, we perform planning on the learned dynamics models
and measure the obtained reward.

##  2 Related Work

####  2.0.1 Latent Space Models:

Previous work on learning latent space models from high dimensional inputs has
produced impressive results especially in the atari learning environment and
continuous control benchmarks [22, 14, 17]. The architectures share a similar
structure where a convolutional encoder embeds inputs to a latent
representation which is fed into a recurrent prediction network for states and
rewards. By stacking predictions on previous predictions, imaginary rollout
trajectories based on the learned dynamics are produced. Those are used for
planning [13], training an actor-critic agent [13, 14] or to drive exploration
[28]. To mitigate accumulating errors caused by imperfect dynamics models,
probabilistic prediction models and long rollout trajectories during training
are used. While the above mentioned methods demonstrate impressive performance
of latent space models in RL, we specifically focus on resource efficiency in
a multi-task setting.

####  2.0.2 Multi-Task RL:

To overcome catastrophic interference and pave the way for lifelong learning,
in [26] a policy network from multiple modules called columns is built. Each
column is a mini network with individual weights while the same architecture
is shared among all columns. For each task, a new column is added and layer-
wise laterally connected to all previous columns. Since the connections are
learned, this enables new columns to selectively benefit from previously
acquired behavior. The authors train their method using A3C [19] and
demonstrate performance improvements by comparing it to various versions of a
single-column baseline network.

By learning task-specific pathways of active subsections within one large
neural network, [9] use a similar concept than [26]. A subsection corresponds
to a small, localized region in the larger network. However, instead of
hardwiring the connections between subsections, they are learned by an
evolutionary strategy. Given an active path, training of the involved
subsections is done via gradient descent. At the same time, the remaining
weights are frozen so that change only occurs in the currently active modules.
The authors describe their approach as an evolutionary version of dropout,
where thinned out networks are evolved instead of randomly generated (as in
regular dropout).

A significant difference between the above mentioned approaches and our method
is that our algorithm learns modular dynamics models instead of modular
policies. Furthermore, in our proposed method lateral connections between the
models are absent since they would introduce dependencies which prevent
arbitrary recombination of the learned task models.

##  3 Method Description

The foundation of our approach lies in the combination of multiple recurrent
dynamics models (RDMs), a vector-quantized variational autoencoder (VQ-VAE)
and a task classification network (TCN) for online task detection. The VQ-VAE
is used to encode the image observation stream from the environment to a
latent space representation which we simply refer to as embedding from now on.
The TCN receives the stream of embeddings and predicts the probability of each
RDM to be responsible for the task at hand. The most probable RDM is then used
to generate predictions in latent space. A graphical overview of the method is
presented in Figure 1. In the following sections, we will explain each of the
components in more detail. While we discuss the hyperparameters we use for
some parts of the architecture, we do not mention all of them and refer to the
link to our code repository in the conclusion for an exhaustive list.

Figure 1: Overview of the multi RDM architecture. In the first timestep,
observation $o_{0}$ is embedded by the encoder of the VQ-VAE. The resulting
embedding is combined with the action and presented to the RDMs and the TCN.
The latter chooses one of the RDMs (visualized by the argmax node), which is
consequently used to predict the embedding for the next timestep. For
simplicity, reward and terminal state probability outputs are not visualized
here. Dashed lines indicate recurrent components that maintain information
persistent between the individual time steps.

###  3.1 Vector-Quantized Variational Autoencoder

To embed an image observation into a latent space, we use the Vector-Quantized
Variational Autoencoder (VQ-VAE) from [23]. It consists of a convolutional
encoder block, a quantization layer and a deconvolutional decoder block. In
the following, we summarize the parameters of all VQ-VAE components with the
parameter vector $\phi$. The encoder transforms an image observation
$\mathbf{o}$ into a 3D tensor $E_{\phi}(\mathbf{o})$ of size $w\times h\times
c$. It can be interpreted as a more compact image of size $w$ times $h$ with
$c$ channels. This image tensor is passed to the quantization layer that
replaces each pixel along the channel dimension with the closest from a set of
codebook vectors. The resulting tensor $\mathbf{z}$ is a per-pixel quantized
version of $E_{\phi}(\mathbf{o})$. Conceptually equivalent to a regular VAE,
the deconvolutional decoder $D$ transforms $\mathbf{z}$ back to the
(reconstructed) original input image $\mathbf{\hat{o}}=D_{\phi}(\mathbf{z})$.
The VQ-VAE is trained by minimizing the following loss function:

$\mathcal{L}_{\phi}(\mathbf{o},D_{\phi}(\mathbf{z}))=\underbrace{||\mathbf{o}-D_{\phi}(\mathbf{z})||^{2}_{2}}_{\text{decoder}}+\underbrace{||sg[E_{\phi}(\mathbf{o})]-\mathbf{z}||^{2}_{2}}_{\text{codebook}}+\beta\underbrace{||sg[\mathbf{z}]-E_{\phi}(\mathbf{o})||^{2}_{2}}_{\text{encoder}}$
(1)

Whereas the parameter $\beta$ trades off encoder vs. decoder loss and the
$sg[\cdot]$ operator indicates the stop of gradient backpropagation. We used
$\beta=0.5$ in all our experiments, since it performed best among a collection
of tested values between $0.25$ and $0.75$. According to the recommendations
of [23], we swap the codebook loss in equation 1 for an exponentially moving
average update of the codebook vectors $e_{i}$:

$N_{i}^{(t)}:=N_{i}^{(t-1)}\gamma+n_{i}^{(t)}(1-\gamma),\hskip
10.00002ptm_{i}^{(t)}:=m_{i}^{(t-1)}\gamma+\sum_{j}^{n_{i}^{(t)}}E(o)_{i,j}^{(t)}(1-\gamma),\hskip
10.00002pte_{i}^{(t)}:=\frac{m_{i}^{(t)}}{N_{i}^{(t)}}$

where $n_{i}^{(t)}$ is for the current batch the number of pixels/vectors in
$E(\mathbf{o})$ that will be replaced with codebook vector $e_{i}$ and
$\gamma$ is a decay parameter between 0 and 1. We used the default
$\gamma=0.99$ for all experiments and a Keras [2, 4] port of the VQ-VAE
implementation in the Sonnet library [1]. The autoencoder architecture of our
choice naturally forms clusters in embedding space, which has been reported to
increase prediction performance of dynamics models like ours [23] and
potentially helps the TCN in differentiating between tasks.

###  3.2 Recurrent Dynamics Models

For the design of our RDMs, we follow ideas from [13] and combine
probabilistic and deterministic model components to increase predictive
performance and prevent overfitting [5]. As previously explained, following
the work in [22] we additionally embed environment observations into latent
space representations $z$ and operate entirely in that latent space from there
on. Since $z$ is a 3D tensor, the architecture of our RDM contains mainly
convolutional layers (conv) to carry through any spatial information contained
in the image observations. An RDM is comprised of four submodules: The belief-
state network (conv), the embedding network (conv), the reward network (dense)
and the terminal state probability network (dense).

As already briefly mentioned, to make the prediction of $z_{t}$ more robust to
uncertainty, we combine a stochastic and a deterministic prediction path. The
deterministic path is realized through a convolutional LSTM network which
computes the belief state $h_{t}$ from the previous belief state $h_{t-1}$,
action $a_{t-1}$ and embedding $z_{t-1}$. Given $h_{t}$, the stochastic
prediction path is a categorical distribution from which $z_{t}$ is sampled.
This way, information about past observations, actions and belief states can
be accumulated in $h_{t}$ and passed on to future time steps. At the same
time, uncertainty regarding the next $z_{t}$ can be expressed through the
sampling process. Because the embeddings are defined completely by the indices
of their codebook vectors, we use a categorical distribution to sample those
indices, namely the Gumbel Softmax [16, 18], which is a reparameterized
approximation to a categorical one-hot distribution. It can be used with
backpropagation and at the same time offers almost discrete sampling of the
required one-hot vectors. The rewards are sampled from a diagonal Gaussian and
the terminal transition probabilities are sampled from a Bernoulli
distribution. The structure of the RDM is summarized as follows:

belief state: $\displaystyle h_{t}=f_{\theta}(h_{t-1},z_{t-1},a_{t-1})$ (2)
embedding: $\displaystyle z_{t}\sim q_{\theta}(z_{t}|h_{t})$ (3) reward:
$\displaystyle r_{t}\sim q_{\theta}(r_{t}|h_{t})$ (4) terminal transition:
$\displaystyle\gamma_{t}\sim q_{\theta}(\gamma_{t}|h_{t})$ (5)

Note that we summarize all RDM component’s parameters into one vector $\theta$
for simplicity. To train the RDM, the following loss function is minimized per
time step:

$\mathcal{L}_{\theta}=\mathcal{L}_{\theta}^{e}+\mathcal{L}_{\theta}^{r}+\mathcal{L}_{\theta}^{\gamma}$
(6)

whereas the individual prediction loss terms are the parameter maximum
likelihood solutions given the transition data:

$\displaystyle\text{embedding}:\ $
$\displaystyle\mathcal{L}_{\theta}^{e}=\operatorname{\mathbb{E}}_{q_{\phi}(z_{t}|o_{t})}\left[-\ln
q_{\theta}(z_{t}|h_{t})\right]$ (7) $\displaystyle\text{reward}:\ $
$\displaystyle\mathcal{L}_{\theta}^{r}=\operatorname{\mathbb{E}}_{q_{\theta}(z_{t}|h_{t}),p(r_{t}|o_{t},a_{t},o_{t+1})}\left[-\ln
q_{\theta}(r_{t}|h_{t},z_{t})\right]$ (8) $\displaystyle\text{terminal}:\ $
$\displaystyle\mathcal{L}_{\theta}^{\gamma}=\operatorname{\mathbb{E}}_{q_{\theta}(z_{t}|h_{t}),p(\gamma_{t}|o_{t},a_{t},o_{t+1})}\left[-\ln
q_{\theta}(\gamma_{t}|h_{t},z_{t})\right]$ (9)

while we denote all learned distributions by $q_{\theta}(\cdot)$ or
$q_{\phi}(\cdot)$ and the true environment dynamics by $p(\cdot)$.

###  3.3 Context Detection

In contrast to the algorithms presented in [13, 11, 22], we use distinct
dynamics models to combat catastrophic interference when training on multiple
tasks. In theory, a single (monolithic) dynamics model can exploit
similarities in the tasks it learns, which positively influences learning
speed and predictive performance. But since weight updates for different tasks
may interfere with each other negatively, a monolithic dynamics model could
instead suffer from reduced performance. This effect has been shown by [33]
for policy networks in multi-task reinforcement learning. To prevent this
catastrophic interference, we isolate knowledge about different task dynamics
in separate RDMs and use the TCN to orchestrate them. This essentially makes
the TCN a context detector that learns to classify an embedding-action-stream
and choose the correct RDM for the task at hand. In every prediction step, the
TCN receives the current observation embedding and action as input. While it
is theoretically able to choose a different RDM in every time step, we did not
observe this behavior in practice. Because different tasks can have locally
similar or equivalent embeddings, the TCN contains convolutional LSTM layers
and is thus able to remember past embeddings and actions. It is thereby
capable to disambiguate temporarily similar embedding streams. The TCN is
trained in a supervised manner via categorical crossentropy loss.

###  3.4 Planning

We directly use our learned RDMs for planning in the environments. Thereby, we
follow the crossentropy method that has originally been introduced by [25] and
is also employed by [11] and [5]. A planning procedure consists of $n$
iterations. At the beginning of the first iteration, the current observation
is embedded into latent space by the VQ-VAE and serves as a starting point for
$k$ rollouts of length $T$. To generate the action sequences required for the
rollouts, we initialize one categorical distribution per time step at random
from which we sample $k$ actions. This results in a $k\times T$ matrix of
actions where each row represents the action sequence for one rollout
trajectory. The resulting trajectories produced by the monolithic or multi RDM
are ranked w.r.t. their discounted return and the top $\rho$ percent are
chosen as winners. In the next iteration, we use the maximum likelihood
parameters of the previous winners for the action distributions and again
sample $k\times T$ actions to generate new rollouts. For the gridworld
setting, we use $n=20$ iterations, $k=500$ rollout trajectories and $T=100$
planning steps, because especially in the third task the reward is far away
from the starting location. We take the top $\rho=10/500$ of all rollout
trajectories in each iteration. For the 3D environment, we use $n=5$
iterations, $k=400$ rollout trajectories and $T=60$ time steps due to
technical limitations regarding the capacity and of the used GPUs (the
algorithm is running on a single Tesla v100 with 16 GB of VRAM). We take the
top $\rho=10/400$ of all rollout trajectories in each iteration. For both the
gridworld and the 3D environment, we add noise of $0.05$ to the maximum
likelihood parameters during planning because this effectively improved the
results.

###  3.5 Training

Presented with a batch of encoded trajectories from $m$ tasks, the TCN outputs
the probabilities for each of the $m$ RDMs per sample (i.e. per trajectory per
time step). Since the individual prediction errors of the RDMs are weighted
with their probability of being chosen, only those with a high probability
receive a strong learning signal and consequently adapt their weights. To
discourage the TCN from switching around RDMs during a task, we provide the
task ID as a supervised learning signal during training. The TCN produces a
less concentrated probability distribution in the beginning of the training
when the prediction error of the RDMs dominates the loss. At this time,
multiple RDMs can be chosen for one task without hurting the loss
significantly, which means that most RDMs are exposed to more than one task
during the course of the training. This increases variety in training samples
for each of the RDMs and encourages them to learn characteristics shared
across tasks. However, this situation only lasts for the first few training
epochs. As soon as the prediction error decreases, the TCN has to allocate the
correct models to their respective tasks in order to further minimize the
loss. The resulting system is then able to activate a task-specific RDM for
each task it was trained on. As a consequence of using distinct networks for
each task, individual RDMs may be re-trained or trained longer on tasks where
the performance is not yet as desired without the risk of sacrificing
performance in other tasks.

We compare the multi RDM model with a larger, monolithic RDM in two different
experiments: A simple 3-task gridworld and a more demanding 2-task setting
comprised of two vizdoom [32] environments. The latter setting can be
considered more demanding since it features a dynamic 3D environment, more
diverse observations and a more complex reward function. We want to emphasize
that per experiment, both the multi RDM and the monolithic RDM have the same
amount of trainable parameters. Since the monolithic RDM does not have a TCN
and to maintain a fair comparison, the architecture of the monolithic RDM is
slightly modified. Another network head is added and trained to output the
current task ID similarly to the TCN in the multi RDM architecture. This way,
the monolithic architecture gets exactly the same training information as the
multi RDM architecture. Training is summarized as follows:

  1. 1.

Collect around 30 000 (gridworld) / 200 000 (vizdoom) transitions from each
task by random action sampling. Since only complete trajectories are
collected, the final number of collected transitions might end up slightly
higher.

  2. 2.

Train the VQ-VAE on the collected samples for 200 (both) epochs.

  3. 3.

Train either the monolithic or the multi RDM predictors on the trajectories
using the VQ-VAE to embed observations. Training is done for 150 (gridworld) /
500 (vizdoom) epochs.

  4. 4.

Perform planning utilizing the trained predictors.

When fully trained, we use the VQ-VAE and RDMs to perform imaginary rollouts
as described in Section 3.4. Since we specifically want to assess the quality
of the learned dynamics for both the monolithic and the multi RDM, we do not
re-plan actions after each step. This way, groundtruth observation data would
be injected into the system continuously, which potentially helps a suboptimal
model. Instead, in the gridworld experiment our algorithm uses only the
starting observation of the agent to perform a full planning routine until the
end of the episode. This has the beneficial side effect of being considerably
less time consuming than re-planning after every step. Note that the starting
observations of the gridworld tasks are unique, so they should in theory
provide sufficient information for the algorithm. In the vizdoom experiment,
one of the two tasks features a dynamic environment that can change without
interaction of the agent. Thus, periodic re-planning is necessary to
incorporate new information from the environment. In an effort to find a good
balance between computational cost and planning accuracy, we use the first 30
actions from the planning procedure and then re-plans from there on.

To obtain the gridworld experiment data discussed in the next section, we
repeated the above mentioned four-step process three times with different
random seeds and averaged the results appropriately. To obtain the vizdoom
experiment data, we did a single iteration of the above mentioned four-step
process. We want to clarify that we trained one pair of dynamics models (a
monolithic and a multi RDM) per experiment. This means, the dynamics models
trained on the gridworld experiment did not receive any data from the vizdoom
experiment and vice versa.

##  4 Experiments

Figure 2: The three training tasks of the gridworld experiment. The agent
(green square in the lower left corner) has to find the reward (red square
marked with ”R”). Every action moves the agent one grid cell into the
respective direction if no obstacles are met. Walls (black squares) block
movement without reward penalty. The shaded blue area in the first environment
(left panel) shows an example of the top-down viewport of the agent. If part
of the viewport is outside of the environment at any time, the respective part
is filled with zeros (i.e. it is colored black).

To answer the question of whether dynamics models can benefit from training on
multiple tasks at once, we use a simple gridworld experiment (cf. Figure 2)
and a more complex vizdoom experiment (cf. Figure 3). For both experiments,
observations from different tasks can easily be confused with each other. Yet,
the system has to learn to identify the tasks and memorize task dynamics
without confusing them. This is especially difficult for the gridworld
experiment, since although having very different global structures, the
individual image observations look similar across all three tasks. On the one
hand, this makes it easier to learn shared properties across them, e.g. that
all objects in the viewport are shifted down one cell if an “up” move is
performed. We expect this to be advantageous for the monolithic RDM. On the
other hand, in the absence of striking visual cues it becomes harder to
distinguish between different tasks, which should put the multi RDM at an
advantage because it uses the specialized TCN for that. In case of the vizdoom
experiment, distinguishing the tasks is easier because there are salient
optical cues in every observation (e.g. the gun and HUD indicating the
“VidzoomBasic” task). We expect this to benefit the monolithic RDM. Yet,
although the environment dynamics are partially similar for both tasks (agent
movement), we see quite different observations and reward functions. We expect
this reduced degree of common structure among the vizdoom tasks to put the
multi RDM at an advantage.

Figure 3: The two training tasks of the vizdoom experiment. In both cases, the
agent is situated in a rectangular room with its back against one wall, facing
the opposite wall. The agent can move sideways left, sideways right and attack
in both tasks. In the ‘VizdoomBasic” task (left image), the agent has to shoot
the monster on the opposite wall. In the “TakeCover” task (right image), the
agent has to dodge incoming fireballs shot by enemies that randomly spawn on
the opposite wall.

####  4.0.1 Gridworld Experiment

All mazes have a size of 10 by 10 cells, a fixed layout per maze for
reproducibility and a maximum episode length of 100 steps. Per environment,
there exists one reward of value $1.0$ at a fixed position. Additionally, the
agent receives a $-0.01$ penalty per step, resulting in a positive reward only
if the agent is able to find the goal state fast enough. During data
collection, the agent starts at a random, unobstructed cell somewhere in the
maze and performs random actions. This takes exploration strategies, which
exceed the scope of this work, out of the equation. To assess the performance
of the learned dynamics models in a comparable way, the agent’s starting
position is fixed to the lower left corner of every maze during control. The
action space is discrete and provides four choices: Up, down, left and right.
The observation emitted by the environment is an image of a 5x5 cells top-down
view centered on the agent position, visualized in Figure 2 in the leftmost
panel. Although we conduct our experiments on a gridworld, we emphasize that
the input space consists of proper image observations.

####  4.0.2 Vizdoom Experiment

In both tasks the agent is positioned in a rectangular 3D room viewed from the
ego perspective (c.f. Figure 3). Standing with its back against one wall, the
agent faces the opposite wall and can choose from three actions: Move sideways
left, move sideways right and attack. Since the agent does not hold a weapon
in the second task, the attack action here simply advances the environment one
time step. In the first task (“VizdoomBasic”), the agent has to shoot the
monster at the opposite wall for a reward of $100$. Every move action results
in a reward penalty of $-1$ and every fired shot in a reward penalty of $-5$.
The maximum number of allowed time steps is 300. In the second task
(“TakeCover”) the agent has to evade fireballs shot by enemies randomly
spawning at the opposite side of the room. Since every time step the agent is
alive results in a reward of $1$, the goal is to survive as long as possible.
Although there is formally no time limit for the second task, more and more
monsters will spawn over time. This makes evading all fireballs increasingly
difficult and ultimately impossible. For both tasks, the agent and monster
starting positions are randomized during training data collection and fixed
for the control experiments. Both, monolithic and multi RDM get the same
sequence of starting positions.

Figure 4: Control results for the gridworld and the vizdoom experiments. Each
task is depicted as a column on the x-axis, in turn separated into two columns
representing the monolithic (blue) or multi (orange) RDM. The y-axis shows the
average episode rewards. Left: Results for the gridworld experiment obtained
via planning rollouts averaged over 60 trials per task. Note that the 60
trials originate from 3 full training runs as explained in Section 3.5 with 20
trials per run. Right: Results for the vizdoom experiment obtained via
planning rollouts averaged over 35 trials per task. Note that the 35 trials
originate from 1 full training run as explained in Section 3.5.

###  4.1 Evaluation

To evaluate the quality of the trained monolithic as well as multi RDM, we use
them to directly generate agent behavior via planning. This way, the
performance of the approaches can be assessed with minimal additional
complexity. The results obtained with our trained RDMs for the gridworld
experiment are shown in Figure 4 (left panel). The overall difficulty increase
from task one to three is reflected in the reward decrease for both
approaches. While the first task can be completed with a comparable
performance by both architectures, the second and third task show significant
differences between the two approaches, the split architecture outperforming
the monolithic network ($p<0.0001$, Mann-Whitney-U-test). The third task is
generally the hardest, since the way to the reward is most obstructed and
longer than for the other tasks. The results for the vizdoom experiment are
shown in Figure 4 (right panel) and follow a similar pattern as those of the
gridworld experiment. The multi RDM architecture manages to deliver
significantly higher average rewards than the monolithic architecture in both
tasks ($p<0.002$, Mann-Whitney-U-test). Note that the maximum reward that can
be obtained in the first task is capped at around $85$ (depending on the agent
and monster starting positions) while it is potentially infinite for the
second task.

Figure 5: For each plot, 64 rollouts of length 50 with randomly chosen initial
observations and random actions were performed. On the y-axis all samples of a
single time step are listed, color-coded depending on their most probable
class. Samples belonging to task (1, 2, 3) are colored in (blue, orange,
green) and samples that could not be assigned to any task are colored in red.
On the x-axis, the distinct time steps are shown. For visual clarity, data
obtained exclusively from one of the three training runs is shown. Top row:
Results of the monolithic architecture. Bottom row: Results of the multi RDM.

With these findings in mind, the question is why the monolithic RDM performs
worse than the multi RDM in planning. To shed light on possible reasons for
the performance difference, we carefully inspected the generated rollout
samples of both approaches in the gridworld experiment setting. Per task and
architecture, we generated 64 trajectories of length 50 from 64 arbitrary
starting observations. The actions were randomly sampled from the tasks’
action spaces. Figure 5 shows the most likely environment of the produced
latent space embedding $z_{i,t}$ for every trajectory $i$ and time step $t$.
To find the most likely environment, each $z_{i,t}$ was decoded into an image
using the VQ-VAE decoder and compared to a list of representative images. This
list contains one image for every possible agent location in the three tasks,
which means that if the decoded $z_{i,t}$ resembles any valid observation,
this was detected with absolute certainty. If the average per-pixel distance
to the closest match was above a certain threshold (here $0.01$), the
$z_{i,t}$ was classified as “undefined” (red in the plots). Otherwise, it was
classified to belong to the environment of the closest match (blue for task 1,
orange for task 2, green for task 3). Note that this analysis is feasible only
for environments with a relatively small amount of states, which is another
reason why we decided to use the gridworld tasks for this analysis.

It is clearly visible that the monolithic RDM produces trajectories of lower
overall quality (cf. Figure 5, top row). First of all, compared to the multi
RDM the fraction of undefined observations is notably higher in tasks 2 and 3.
Even more importantly, the monolithic RDM has a higher chance to generate $z$
that belong to a different task, most evidently reflected in task 2 (cf.
Figure 5, top row, middle panel). This can lead to malformed rollout
trajectories which miss rewards or report rewards where there are none,
harming the planning process in general. In the worst case, the monolithic RDM
switches tasks in the middle of a rollout trajectory. By visual inspection of
the rollout trajectories, we observed this to happen regularly for the
monolithic RDM in both the gridworld and the vizdoom experiments. Through
separation of the task dynamics in the multi RDM, considerably less confusion
among the tasks arises and the fraction of undefined samples is reduced as
well. This stands in line with measurably better planning performance for the
multi RDM architecture. We conclude from our analysis that catastrophic
interference outweighs possible advantages of transfer learning in the
monolithic RDM. In our experiments, the simple use of isolated dynamics models
mitigates the harmful interference to a large degree.

##  5 Conclusion

In this work, we demonstrate an unintuitive effect in model-based multi-task
RL. Contrary to expectations fueled by model-free multi-task RL, using a
single monolithic RDM for all tasks can harm performance in model-based multi-
task RL instead of improving it. We conclusively assume that the positive
effect of transfer learning between similar tasks is either absent for
dynamics models or outweighed by catastrophic interference. Moreover, while
using the same parameter budget for both approaches we show that imposing an
internal structure to the dynamics model can lead to notably improved learned
dynamics. We show this by using both approaches for control through planning,
where the multi RDM scores measurably higher rewards than the monolithic RDM.
By analyzing the rollout trajectories of both the monolithic and the multi RDM
in the gridworld setting, we find cleaner, more task-specific trajectories and
less overall wrong predictions in case of the multi RDM. We conclude that
separating task recognition from task dynamics effectively prevents
catastrophic interference (or task confusion) to a large degree, leading to
the measured performance improvement. The code for reproducing our
experiments, including all hyperparameters, can be found at
https://github.com/rschiewer/lrdm.

Acknowledgements: We thank Jan Bollenbacher, Dr. Anand Subramoney and Prof.
Dr. Tobias Glasmachers for their feedback and help, which greatly influenced
this work.

## References

  * [1] Sonnet (2017), https://github.com/deepmind/sonnet
  * [2] Abadi, M., et al.: TensorFlow: Large-scale machine learning on heterogeneous systems (2015), https://www.tensorflow.org/
  * [3] Caruana, R.: Multitask learning. Machine Learning 28 (07 1997). https://doi.org/10.1023/A:1007379606734
  * [4] Chollet, F., et al.: Keras. https://keras.io (2015)
  * [5] Chua, K., Calandra, R., McAllister, R., Levine, S.: Deep Reinforcement Learning in a Handful of Trials using Probabilistic Dynamics Models. Advances in Neural Information Processing Systems 2018-Decem(Nips), 4754–4765 (may 2018), http://arxiv.org/abs/1805.12114
  * [6] Duan, Y., Schulman, J., Chen, X., Bartlett, P.L., Sutskever, I., Abbeel, P.: Rl2: Fast reinforcement learning via slow reinforcement learning. arXiv preprint arXiv:1611.02779 abs/1611.02779 (2016), http://arxiv.org/abs/1611.02779
  * [7] Ebert, F., Finn, C., Dasari, S., Xie, A., Lee, A., Levine, S.: Visual foresight: Model-based deep reinforcement learning for vision-based robotic control. arXiv preprint arXiv:1812.00568 (2018)
  * [8] Feinberg, V., Wan, A., Stoica, I., Jordan, M.I., Gonzalez, J.E., Levine, S.: Model-based value expansion for efficient model-free reinforcement learning. In: Proceedings of the 35th International Conference on Machine Learning (ICML 2018) (2018)
  * [9] Fernando, C., Banarse, D., Blundell, C., Zwols, Y., Ha, D., Rusu, A.A., Pritzel, A., Wierstra, D.: PathNet: Evolution channels gradient descent in super neural networks. arXiv preprint arXiv:1701.08734 (2017)
  * [10] Finn, C., Abbeel, P., Levine, S.: Model-agnostic meta-learning for fast adaptation of deep networks. In: International Conference on Machine Learning. pp. 1126–1135. PMLR (2017)
  * [11] Ha, D., Schmidhuber, J.: Recurrent world models facilitate policy evolution. arXiv preprint arXiv:1809.01999 abs/1809.01999 (2018), http://arxiv.org/abs/1809.01999
  * [12] Haarnoja, T., Zhou, A., Abbeel, P., Levine, S.: Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor. In: International Conference on Machine Learning. pp. 1861–1870. PMLR (2018)
  * [13] Hafner, D., Lillicrap, T., Fischer, I., Villegas, R., Ha, D., Lee, H., Davidson, J.: Learning latent dynamics for planning from pixels. In: International Conference on Machine Learning. pp. 2555–2565. PMLR (2019)
  * [14] Hafner, D., Lillicrap, T.P., Norouzi, M., Ba, J.: Mastering atari with discrete world models. arXiv preprint arXiv:2010.02193 abs/2010.02193 (2020), https://arxiv.org/abs/2010.02193
  * [15] Hessel, M., Modayil, J., van Hasselt, H., Schaul, T., Ostrovski, G., Dabney, W., Horgan, D., Piot, B., Azar, M.G., Silver, D.: Rainbow: Combining improvements in deep reinforcement learning. In: AAAI. pp. 3215–3222 (2018), https://www.aaai.org/ocs/index.php/AAAI/AAAI18/paper/view/17204
  * [16] Jang, E., Gu, S., Poole, B.: Categorical reparameterization with gumbel-softmax. In: 5th International Conference on Learning ReprFesentations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net (2017)
  * [17] Kaiser, L., Babaeizadeh, M., Milos, P., Osinski, B., Campbell, R.H., Czechowski, K., Erhan, D., Finn, C., Kozakowski, P., Levine, S., Mohiuddin, A., Sepassi, R., Tucker, G., Michalewski, H.: Model-Based Reinforcement Learning for Atari. arXiv preprint arXiv:1903.00374 (mar 2019), http://arxiv.org/abs/1903.00374
  * [18] Maddison, C.J., Mnih, A., Teh, Y.W.: The concrete distribution: A continuous relaxation of discrete random variables. In: 5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings. OpenReview.net (2017)
  * [19] Mnih, V., Badia, A.P., Mirza, M., Graves, A., Lillicrap, T., Harley, T., Silver, D., Kavukcuoglu, K.: Asynchronous methods for deep reinforcement learning. In: International conference on machine learning. pp. 1928–1937. PMLR (2016)
  * [20] Mnih, V., Kavukcuoglu, K., Silver, D., Graves, A., Antonoglou, I., Wierstra, D., Riedmiller, M.: Playing atari with deep reinforcement learning. arXiv preprint arXiv:1312.5602 (2013)
  * [21] Nagabandi, A., Kahn, G., Fearing, R.S., Levine, S.: Neural network dynamics for model-based deep reinforcement learning with model-free fine-tuning. In: 2018 IEEE International Conference on Robotics and Automation (ICRA). pp. 7559–7566 (2018). https://doi.org/10.1109/ICRA.2018.8463189
  * [22] Oh, J., Guo, X., Lee, H., Lewis, R., Singh, S.: Action-Conditional Video Prediction using Deep Networks in Atari Games. CoRR 10(1), 52–55 (jul 2015), http://papers.nips.cc/paper/5859-action-conditional-video-prediction-using-deep-networks-in-atari-games.pdf http://arxiv.org/abs/1507.08750
  * [23] van den Oord, A., Vinyals, O., kavukcuoglu, k.: Neural discrete representation learning. In: Guyon, I., Luxburg, U.V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., Garnett, R. (eds.) Advances in Neural Information Processing Systems. vol. 30. Curran Associates, Inc. (2017), https://proceedings.neurips.cc/paper/2017/file/7a98af17e63a0ac09ce2e96d03992fbc-Paper.pdf
  * [24] Racanière, S., Weber, T., Reichert, D.P., Buesing, L., Guez, A., Rezende, D., Badia, A.P., Vinyals, O., Heess, N., Li, Y., et al.: Imagination-augmented agents for deep reinforcement learning. In: Proceedings of the 31st International Conference on Neural Information Processing Systems. pp. 5694–5705 (2017)
  * [25] Rubinstein, R.Y.: Optimization of computer simulation models with rare events. European Journal of Operational Research 99(1), 89–112 (1997). https://doi.org/https://doi.org/10.1016/S0377-2217(96)00385-2, https://www.sciencedirect.com/science/article/pii/S0377221796003852
  * [26] Rusu, A.A., Rabinowitz, N.C., Desjardins, G., Soyer, H., Kirkpatrick, J., Kavukcuoglu, K., Pascanu, R., Hadsell, R.: Progressive neural networks. arXiv preprint arXiv:1606.04671 abs/1606.04671 (2016), http://arxiv.org/abs/1606.04671
  * [27] Schulman, J., Wolski, F., Dhariwal, P., Radford, A., Klimov, O.: Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347 (2017)
  * [28] Sekar, R., Rybkin, O., Daniilidis, K., Abbeel, P., Hafner, D., Pathak, D.: Planning to explore via self-supervised world models. In: International Conference on Machine Learning. pp. 8583–8592. PMLR (2020)
  * [29] Silver, D., Hubert, T., Schrittwieser, J., Antonoglou, I., Lai, M., Guez, A., Lanctot, M., Sifre, L., Kumaran, D., Graepel, T., Lillicrap, T., Simonyan, K., Hassabis, D.: A general reinforcement learning algorithm that masters chess, shogi, and Go through self-play. Science 362(6419), 1140–1144 (dec 2018). https://doi.org/10.1126/science.aar6404, http://science.sciencemag.org/
  * [30] Sutton, R.S.: Dyna, an integrated architecture for learning, planning, and reacting. SIGART Bull. 2(4), 160–163 (Jul 1991). https://doi.org/10.1145/122344.122377, https://doi.org/10.1145/122344.122377
  * [31] Tesauro, G.: Temporal difference learning and TD-Gammon. Communications of the ACM 38(3), 58–68 (mar 1995). https://doi.org/10.1145/203330.203343, https://dl.acm.org/doi/abs/10.1145/203330.203343
  * [32] Wydmuch, M., Kempka, M., Jaśkowski, W.: Vizdoom competitions: Playing doom from pixels. IEEE Transactions on Games (2018)
  * [33] Yu, T., Kumar, S., Gupta, A., Levine, S., Hausman, K., Finn, C.: Multi-Task Reinforcement Learning without Interference. NeurIPS pp. 1–9 (2019)