# Passive self-propulsion and enhanced annihilation of $p-$atic defects

Dimitrios Krommydas  Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506,
2300 RA Leiden, The Netherlands  Livio Carenza  Instituut-Lorentz,
Universiteit Leiden, P.O. Box 9506, 2300 RA Leiden, The Netherlands  Luca
Giomi<EMAIL_ADDRESS>Instituut-Lorentz, Universiteit Leiden,
P.O. Box 9506, 2300 RA Leiden, The Netherlands

###### Abstract

We investigate the effects of hydrodynamics on the motion of topological
defects in $p-$atic liquid crystals: i.e. two-dimensional liquid crystals
characterised by $p-$fold rotational symmetry. Using numerical simulations and
analytical work, we demonstrate the existence of a generic passive self-
propulsion mechanism for defects of winding number $s=(p-1)/p$ and arbitrary
$p$ values. Remarkably, whereas this mechanism is not unique to nematics (i.e.
$p=2$), nematics are the only type of liquid crystals where passive self-
propulsion is thermodynamically stable. Furthermore, we find that
hydrodynamics can accelerate the annihilation dynamics of pairs of $\pm 1/p$
defects and that this effect increases with $p$ before approaching a plateau.

The physics of topological defects in liquid crystals have experienced, in the
last decade, a tremendous revival, thanks to a wealth of exciting discoveries
at the interface between soft condensed matter and biological physics Sanchez
_et al._ (2012); Keber  _et al._ (2014); Giomi  _et al._ (2014a); Giomi
(2015); Doostmohammadi  _et al._ (2017); Guillamat  _et al._ (2018);
Doostmohammadi  _et al._ (2018); Carenza  _et al._ (2019a, b). The most common
class of liquid crystals defects, known as disclinations, consists of point or
line singularities around which the average orientation of the anisotropic
building blocks undergoes one or more complete revolutions, thereby disrupting
the local orientational order Mermin (1979); Chaikin and Lubensky (1995);
Giomi  _et al._ (2017); Pollard  _et al._ (2019); Carenza  _et al._ (2022);
Lopez-Leon and Fernandez-Nieves (2011). In $p-$atic liquid crystals $-$ i.e.
two-dimensional liquid crystals with $p-$fold rotational symmetry, among which
nematics (i.e. $p=2$) and hexatics (i.e. $p=6$) are the best known examples
$-$ defects can be classified in terms of the winding number or strength $s$,
defined as the number of revolutions of the orientation field along an
arbitrary loop enclosing the core of the defect: i.e. $s=\pm 1/p,\,\pm 2/p$
etc.

Whereas the equilibrium physics of liquid crystal defects represents a mature
topic across several areas of physics $-$ from cosmology Kibble (1976), down
to condensed matter Chaikin and Lubensky (1995); Halperin and Nelson (1978)
and particle physics Dirac (1931); Hooft (1974); Polyakov (1974) $-$ our
understanding of their dynamics is still in a phase of rapid expansion,
especially in the realm of biological matter, where defects have been
suggested to accomplish various vital functions. These include driving the
extrusion of apoptotic cells in epithelial layers Saw  _et al._ (2017);
Monfared  _et al._ (2021), coordinating large scale cellular flows during
wound healing and morphogenetic events Brugués  _et al._ (2014); Kawaguchi
_et al._ (2017), and seeding the development of non-planar features, such as
tentacles and protrusion in simple organisms, such as Hydra Livshits  _et al._
(2017); Braun and Keren (2018); Livshits  _et al._ (2021); Maroudas-Sacks  _et
al._ (2021).

Whilst the biochemical and chemotactic apport to these processes is mostly
understood, much less is known about the physical interactions. Nevertheless,
their origin can be single-handedly ascribed to the existence of a
hydrodynamic phenomenon known as backflow. This effect $-$ whose name directly
refers to the “bounce” in the optical transmission of a twisted nematic cell
between polarizers after switching off the applied field Berreman (1975);
Brochard (1973); Kos Ž (2020); Tsuji and Chono (2021); Svenšek and Žumer
(2001); Missaoui  _et al._ (2021) $-$ consists of the hydrodynamic flow
resulting from spatial variations of the average microscopic orientation. In
passive liquid crystals, departure from the uniformly oriented equilibrium
configuration is generally transient and often originates from a sudden change
in the environmental conditions, such as the abrupt variation of an external
electric or magnetic field in optical devices Brochard  _et al._ (1972);
Tiribocchi  _et al._ (2014). Conversely, in active systems, distortions occur
spontaneously as a consequence of the internal stresses collectively exerted
by the active subunits Sanchez  _et al._ (2012); Hatwalne  _et al._ (2004);
Giomi (2015).

Whether passive or active, topological defects play a primary role with
respect to backflow phenomena, by virtue of the extended and yet persistence
distortion resulting from their appearance. In passive nematic liquid
crystals, for instance, this effect is known to affect the annihilation
dynamics of neutral pairs of elementary disclinations Tóth  _et al._ (2003),
while in active nematics $-$ such as in vitro mixture of cytoskeletal
filaments and motor proteins Sanchez  _et al._ (2012); Keber  _et al._ (2014);
Martínez-Prat  _et al._ (2021); Colen  _et al._ (2021) or certain kind of
prokaryotic Dunkel  _et al._ (2013); Wensink  _et al._ (2012) or eukaryotic
Saw  _et al._ (2017); Balasubramaniam  _et al._ (2021) cells $-$ backflow
drives the propulsion of $s=1/2$ defects and influences the hydrodynamic
stability of active layers with respect to non-planar deformations Keber  _et
al._ (2014); Metselaar  _et al._ (2019); Hoffmann  _et al._ (2022); Alert
(2022). Importantly, a hydrodynamic theory capable of capturing backflows
effects in liquid crystals with generic $p-$atic symmetry has been developed
very recently Giomi  _et al._ (2021a, b). Hence, the role of hydrodynamics in
dynamics of defects in liquid crystalline systems with generic $p-$fold
symmetry is mostly unknown.

In this article we take a step in bridging this gap. Leveraging on recent
progress toward generalizing the classic hydrodynamic theory of hexatic liquid
crystals Zippelius (1980); Zippelius  _et al._ (1980) to account for arbitrary
$p-$fold symmetry symmetry, we construct the velocity field of arbitrary
$p-$atic disclinations and demonstrate the existence of a generic passive
self-propulsion mechanism for defects of winding number $s=(p-1)/p$ and
arbitrary $p$ values. Remarkably, whereas this mechanism is not unique to
nematics, nematics are the only type of liquid crystals where passive self-
propulsion is thermodynamically stable. Furthermore, we analyze the effect of
hydrodynamics in the annihilation of pairs of elementary $\pm 1/p$
disclinations and find that backflow always accelerates the annihilation
dynamics and, contrary to expectations, becomes increasingly more relevant as
$p$ increases.

Let us consider an incompressible $p-$atic liquid crystal, whose microscopic
orientation is characterized by the unit vector
$\bm{\nu}=\cos\vartheta\,\bm{e}_{x}+\sin\vartheta\,\bm{e}_{y}$ and whose
physical properties are invariant under rotations by $2\pi/p$. For $p=2$, this
is the direction of the rod-like building blocks comprising a nematic liquid
crystals, whereas for $p=3$, $\bm{\nu}$ is either one of the three equivalent
directions of $3-$legged star, etc (see e.g. Fig. 1). At length scales larger
than that of the individual constituents and yet infinitesimal compared to the
system size, $p-$atic order can be conveniently described in terms of the
tensor order parameter
$\bm{Q}_{p}=Q_{i_{1}i_{2}\ldots\,i_{p}}\bm{e}_{i_{1}}\otimes\bm{e}_{i_{2}}\otimes\ldots\otimes\bm{e}_{i_{p}}$,
where $i_{n}=\\{x,y\\}$ and $n=1,\,2\,\ldots\,p$, constructed upon averaging
the $p-$fold tensorial power of the local orientation $\bm{\nu}$. That is:
$\bm{Q}_{p}=\sqrt{2^{p-2}}\left\llbracket\bm{\nu}^{\otimes
p}\right\rrbracket=|\Psi_{p}|\left\llbracket\bm{n}^{\otimes
p}\right\rrbracket$ , where
$\bm{n}=\cos\theta\,\bm{e}_{x}+\sin\theta\,\bm{e}_{y}$ and
$|\Psi_{p}|=|\Psi_{p}|(\bm{r})$ and $\theta=\theta(\bm{r})$ are respectively
the magnitude and the phase of the coarse-grained complex order parameter
$\langle e^{ip\vartheta}\rangle=|\Psi_{p}|e^{ip\theta}$, while the operator
$\left\llbracket\cdots\right\rrbracket$ renders its argument traceless (i.e.
$Q_{jji_{3}\ldots\,i_{p}}=0$) and symmetric with respect to the exchange of
any two indices Giomi  _et al._ (2021a, b). For nematics (i.e. $p=2$), this
construction gives the standard rank$-2$ tensor
$\bm{Q}_{2}=|\Psi_{2}|(\bm{n}\otimes\bm{n}-\mathbb{1}/2)$, with $\mathbb{1}$
the identity tensor and $|\Psi_{2}|=\langle
2|\bm{\nu}\cdot\bm{n}|^{2}-1\rangle$ the scalar order parameter. The dynamics
of the fluid is then governed by the following set hydrodynamic equations for
momentum density $\rho\bm{v}$ and the tensor order parameter $\bm{Q}_{p}$:

$\displaystyle\rho\,\frac{D\bm{v}}{Dt}=\nabla\cdot\bm{\sigma}\,,$ (1a)
$\displaystyle\frac{D\bm{Q}_{p}}{Dt}=\Gamma_{p}\bm{H}_{p}+p\big{\llbracket}\bm{Q}_{p}\cdot\bm{\omega}\big{\rrbracket}$
$\displaystyle+\lambda_{p}\big{\llbracket}\nabla^{\otimes
p-2}\bm{u}\big{\rrbracket}+\nu_{p}\big{\llbracket}\nabla^{\otimes p\,{\rm
mod}\,2}\bm{u}^{\otimes\lfloor p/2\rfloor}\big{\rrbracket}\,,$ (1b)

where $D/Dt=\partial_{t}+\bm{v}\cdot\nabla$ is the material derivative, $\rho$
is a constant density, $\bm{v}$ is the incompressible velocity field
($\nabla\cdot\bm{v}=0$) and $\bm{\sigma}$ is the total stress tensor, to be
defined later. The tensors $\bm{u}=[\nabla\bm{v}+(\nabla\bm{v})^{\rm T}]/2$
and $\bm{\omega}=[\nabla\bm{v}-(\nabla\bm{v})^{\rm T}]/2$, with ${\rm T}$
indicating transposition, are respectively the strain rate and vorticity
fields and entail the coupling between $p-$atic order and flow, with
$\lambda_{p}$ and $\nu_{p}$ material constants. The dot product indicates the
contraction between the last index of $\bm{Q}_{p}$ and the first index of
$\bm{\omega}$ $-$ i.e.
$(\bm{Q}_{p}\cdot\bm{\omega})_{i_{1}i_{2}\ldots\,i_{p}}=Q_{i_{1}i_{2}\ldots\,j}\omega_{ji_{p}}$
$-$ and $\left(\nabla^{\otimes
n}\right)_{i_{1}i_{2}...i_{n}}=\partial_{i_{1}}\partial_{i_{2}}\ldots\,\partial_{i_{n}}$,
whereas $\lfloor\cdots\rfloor$ denotes the floor function and $p\,{\rm
mod}\,2=p-2\lfloor p/2\rfloor$ is zero for even $p$ values and one for odd $p$
values. The quantity $\bm{H}_{p}=-\delta F/\delta\bm{Q}_{p}$ is the $p-$atic
analog of the molecular tensor, dictating the relaxation dynamics of the order
parameter tensor toward the minimum of the orientational free energy
$F=\int{\rm
d}A\,\left(L_{p}/2\,|\nabla\bm{Q}_{p}|^{2}+A_{p}/2\,|\bm{Q}_{p}|^{2}+B_{p}/4\,|\bm{Q}_{p}|^{4}\right)$,
where $|\cdots|^{2}$ is the Euclidean norm and is such that
$|\bm{Q}_{p}|^{2}=|\Psi_{p}|^{2}/2$. The constant $L_{p}$ is the order
paramter stiffness and $A_{p}$ and $B_{p}$ are phenomenological constants
setting the magnitude of the scalar order parameter at equilibrium: i.e.
$|\Psi_{p}|=|\Psi_{p}^{(0)}|=\sqrt{-2A_{p}/B_{p}}$, when
$H_{i_{1}i_{2}\cdots\,i_{p}}=0$.

Figure 1: Examples of elementary $s=1/p$ (red) and $s=-1/p$ (blue) defects in
$p-$atics and their associated backflow (black and white) for
$p=2,\,3\ldots\,6$. The streamlines are obtained from Eq. (4).

The stress tensor on the left-hand side of Eq. (1a) is customarily decomposed
into a static and a dynamic contribution: i.e. $\bm{\sigma}=\bm{\sigma}^{({\rm
s})}+\bm{\sigma}^{({\rm d})}$. The static stress tensor is given by
$\bm{\sigma}^{({\rm s})}=-P\mathbb{1}+\bm{\sigma}^{({\rm e})}$, where $P$ is
the pressure and $\sigma_{ij}^{({\rm
e})}=-L_{p}\partial_{i}\bm{Q}_{p}\odot\partial_{j}\bm{Q}_{p}$ is the elastic
stress resulting from a static distortion of the $p-$atic orientation, with
the symbol $\odot$ indicating a contraction of all matching indices of the two
operands yielding a tensor whose rank equates the number of unmatched indices
(two in this case). The dynamic stress, on the other hand, can be further
decomposed into a viscous or energy dissipating part and a reactive or energy
preserving part: i.e. $\bm{\sigma}^{({\rm d})}=\bm{\sigma}^{({\rm
v})}+\bm{\sigma}^{({\rm r})}$. The former given by $\bm{\sigma}^{({\rm
v})}=2\eta\left\llbracket\bm{u}\right\rrbracket$, with $\eta$ the shear
viscosity; while the latter takes the form

$\displaystyle\bm{\sigma}^{({\rm r})}=\lambda_{p}(-1)^{p-1}\nabla^{\otimes
p-2}\odot\bm{H}_{p}+\frac{p}{2}\left(\bm{Q}_{p}\cdot\bm{H}_{p}-\bm{H}_{p}\cdot\bm{Q}_{p}\right).$
(2)

Both this terms result from a departure from the lowest free energy state and,
together with the elastic stress $\bm{\sigma}^{({\rm e})}$, can drive backflow
(see e.g. Ref. Oswald and Pieranski (2005)).

To investigate the role of backflow in the dynamics of topological defects, we
next consider the case of an isolated disclination of strength $s=\pm
1/p,\,\pm 2/p\ldots$ at the origin of an unbounded domain and take
$\theta=s\phi+\theta_{0}$, with $\phi=\arctan(y/x)$. The constant angle
$\theta_{0}$ determines the overall orientation of the defect Vromans and
Giomi (2016) and can be set to zero without loss of generality: i.e.
$\theta_{0}=0$. The order parameter $|\Psi_{p}|$ is assumed uniformly equal to
its equilibrium value $|\Psi_{p}^{(0)}|$ outside the core of the defects $-$
i.e. for $|\bm{r}|>a$ with $a\sim\mathcal{O}(\sqrt{-A_{p}/L_{p}})$ the defect
core radius $-$ and vanishing inside. Following Refs. Giomi  _et al._ (2013,
2014b), we assume flow alignment effects to be negligible and compute the
stationary solution of Eq. (1a). Remarkably, of the three backflow-driving
terms entering in the total stress tensor, only one contributes to the flow in
the surrounding of a defect $-$ i.e. the second term in Eq. (2) $-$ resulting
from the interplay between $p-$atic order and flow. By contrast, the second
term in Eq. (2), originating from the corotational derivative of the order
parameter tensor, is proportional to $\nabla^{2}\theta$ when $|\Psi_{p}|={\rm
const}$ (see e.g. Ref. Giomi  _et al._ (2021b)), thus it vanishes identically
in this case. The elastic stress, on the other hand, yields the isotropic
force density $\nabla\cdot\bm{\sigma}^{({\rm
e})}=(ps|\Psi_{p}|)^{2}L_{p}\,\bm{r}/(2|\bm{r}|^{4})$, which, in turn, leads
to a local pressure variation $-$ i.e. $P\to
P-(ps|\Psi_{p}|)^{2}L_{p}/(2|\bm{r}|^{2})$ $-$ without influencing the flow.
The backflow sourced by isolated defects, can then be found by convoluting the
two-dimensional Oseen-Green tensor SI  with the body force
$\bm{f}=\nabla\cdot\bm{\sigma}^{({\rm
r})}=c_{p}/|\bm{r}|^{p+1}[\cos(n\phi)\,\bm{e}_{x}+\sin(n\phi)\,\bm{e}_{y}]$,
where $c_{p}$ and $n$ are given by

$\displaystyle
c_{p}=\frac{(-1)^{p}(ps)^{2}\lambda_{p}L_{p}}{2^{p/2}}~{}\prod^{p-1}_{k=1}~{}\left[ps-2(p-k)\right]\;,$
(3a) $\displaystyle n=p(s-1)+1\;.$ (3b)

We stress that Eqs. (3) holds for arbitrary $p$ values, thus also for nematics
(i.e. $p=2$), where, despite the hydrodynamics of topological defects having
being thoroughly investigated Tóth  _et al._ (2003); Kats  _et al._ (2002);
Svenšek and Žumer (2002); Missaoui  _et al._ (2020); Kos and Ravnik (2020),
the origin of defect-induced backflow had, to the best of our knowledge, never
been clarified in these terms.

Eqs. (3) allows one to approximatively compute the velocity field in the
surrounding of isolated defects. This is given by

$\displaystyle\frac{\bm{v}}{c_{p}/(4\eta)}=\frac{\delta_{n0}}{(p-1)^{2}}\left[\left(\frac{|\bm{r}|}{R}\right)^{1-p}+(p-1)\left(\frac{a}{R}\right)^{1-p}\left(\log\frac{|\bm{r}|}{R}+\frac{3}{2}\right)-(p+1)\right]\,R^{1-p}\bm{e}_{x}+\frac{\delta_{n1}}{p-2}\,\left[\left(\frac{a}{|\bm{r}|}\right)^{2-p}-1\right]\frac{\bm{r}}{|\bm{r}|^{p}}$
$\displaystyle-\frac{1-\delta_{n0}}{n}\,\left\\{\frac{1}{n+p-1}\left(\frac{|\bm{r}|}{a}\right)^{n+p-1}+\frac{1}{n-p+1}\,\left[\left(\frac{|\bm{r}|}{R}\right)^{n+1-p}-1\right]\right\\}|\bm{r}|^{1-p}(\cos
n\phi\,\bm{e}_{x}+\sin n\phi\,\bm{e}_{y})$
$\displaystyle-(1-\delta_{n1})\cos[(n-1)\phi]\left\\{\frac{1}{n+p-1}\left(\frac{|\bm{r}|}{a}\right)^{n+p-1}+\frac{1}{n+p-3}\left[\frac{2}{n+p-1}-\left(\frac{|\bm{r}|}{a}\right)^{n+p-3}\right]\right\\}\frac{\bm{r}}{|\bm{r}|^{p}}\;,$
(4a)

where $R$ is a large-distance cut-off, analogous to the system size. A plot of
the flow associated with this velocity field is given in Fig. 1 for
$p=1,\,2\ldots\,6$ and $s=\pm 1/p$. As intuitive, these flows are
characterized by the same $(|n|+1)-$fold rotational symmetry of the driving
force $\bm{f}$, with $n$ as given in Eq. (3b). Thus, a $s=1/2$ disclination in
nematics, sources a typical Stokeslet-like flow consisting of two counter-
rotating vortices meeting along the defect’s longitudinal direction (i.e. the
$x-$direction in this case), whereas a $s=-1/2$ disclinations gives rise to a
$3-$fold symmetric flow consisting of six vortices with alternating positive
and negative vorticity. Similarly, a $s=1/3$ ($s=-1/3$) disclinations in
triatics drives a $2-$fold ($4-$fold) symmetric flow, etc. Because of this
rotational symmetry, and with only exception for $p=2$, these flows have the
effect of stirring the fluid around a defect and stablize the core by trapping
this at the central stagnation point. By contrast, the Stokeslet-like flow
sourced by a $s=1/2$ disclination yields a net momentum current, whose effect
is to propel the defect forward.

Such an interesting example of passive self-propulsion is not unique of
nematics (Fig. 2a-c), but nematics are the only type of liquid crystal where
passive self-propulsion is thermodynamically stable. To substantiate the
latter statement we express the defects strength in terms of the so called
topological charge $q=ps$ and note that the constant $c$ in Eq. (3a), hence
the flow speed, vanishes for any even $q$ value in the range $2\leq q\leq
2p-2$ (Fig. 2d). For odd $q$ values, the flow speed is instead always finite,
but the $(|n|+1)-$fold symmetry of the force $\bm{f}$ is inherited by the
flow, thereby preventing the occurrence of the Stokeslet-like structure
required for self-propulsion, unless $n=0$, or, equivalently, $s=(p-1)/p$. In
this case, $\bm{f}=c_{p}/|\bm{r}|^{p+1}\,\bm{e}_{x}$ and the symmetry of the
flow is explicitly broken, ultimately leading to the propulsion of the defect
in the $x-$direction at a speed

$v_{\rm d}=v_{0}~{}R^{1-p}\\\
\left\\{p+1+\left(\frac{a}{R}\right)^{1-p}\left[\frac{3p-1}{2}+(p-1)\log\frac{a}{R}\right]\right\\}\;,$
(5)

where
$v_{0}=(-1)^{p+1}\pi\lambda_{p}L_{p}/(2^{p/2+1}\eta)\prod_{k=1}^{p-1}(k-p-1)$.
With exception of $p=2$, however, none of these self-propelled defects
features an elementary winding number $s=\pm 1/p$ and when let free to evolve
it rapidly splits into $p-1$ elementary defects. Finally, the proportionality
between the speed of isolated $s=1/2$ disclinations and the flow alignment
parameter, brings to light an exciting opportunity for estimating the flow
alignment parameter $\lambda_{2}$ $-$ a notoriously elusive material parameter
in liquid crystals (see e.g. Ref. de Andrade Lima and Rey (2004)) $-$ from
measurements of the self-propulsion speed of elementary nematic defects.

Figure 2: (a) Example motile defects with $s=1/2$, $3/4$ and $5/6$. (b) Speed
of the defect core for $0\leq s\leq 1$. Notice that speed is zero for all
defects of strength $s=1/p$ apart from the case of nematics, where $s=1/2$.

Next we focus on the annihilation dynamics of pairs of elementary $s=\pm 1/p$
defect pairs. In the absence of hydrodynamics effects, two-dimensional
disclinations of opposite winding numbers are known to attract each other via
a Coulomb-like force and eventually annihilate. In nematics, Tóth et al.
showed that hydrodynamic affects this process in a two-fold way Tóth  _et al._
(2003): on the one hand, advection by the backflow causes the defects to move
faster, thereby speeding up their annihilation dynamics; on the other hand,
the different configuration of the velocity field in the surrounding of
positive and negative defects introduces an asymmetry in the annihilation
trajectory (Fig. 3a), which is then no longer symmetric about the mid-plane
separating the defects at $t=0$. Although tempting to explain this
phenomenology in the light of the aforementioned passive-self propulsion of
$s=1/2$ defects, in the following we demonstrate that both effects arise
instead from the antisymmetric part of the dynamic stress tensor: i.e. the
second term on the right-hand side of Eq. (2).

Figure 3: Power $\mathcal{P}$ delivered by each and every contribution to the
total stress, for $p=2,\,3\ldots\,6$. The power of dominant contribution to
the backflow (yellow tones) is sourced by the antisymmetric component of the
stress and dissipated by viscous stresses (red tones). Inset: magnification of
the power delivered by all components of the stress in proximity of the
annihilation of elementary nematic defects (i.e. $p=2$, $s=\pm 1/2$). In this
case only the stresses originating from the flow alignment behavior of
nematics (blue) and that elicit the self-propulsion of isolated $s=1/2$
defects, are comparable in magnitude with those arising from the antisymmetric
stresses. The power has been rescaled by $\mathcal{P}_{0}=L_{p}/\tau$, with
$\tau=|\Delta x(0)|^{2}/(\Gamma_{p}L_{p})$, with $\Delta x(0)$ the initial
distance between the defects.

To demonstrate the latter statement we show in Fig. 3 the power
$\mathcal{P}=\int{\rm d}A\,[\nabla\cdot\bm{\sigma}]\cdot\bm{v}$ delivered by
each and every contribution to the total stres, before, during and after an
annihilation event Carenza  _et al._ (2020a, b). Data are generated by a
numerical integration of Eqs. (1) on a square domain with periodic boundary
and whose initial configuration consists of a neutral pair of elementary
$p-$atic defects SI . As evident from our numerical data, for all $p$ values
the annihilation dynamics is dominated by the antisymmetric part of the stress
tensor (yellow tones). This converts the energy stored in the distorted
configuration of the $p-$atic director into kinetic energy, which is in turn
dissipated by viscous stresses (red tones). By contrast the stresses
originating from flow alignment, which source the propulsion of isolated
defects, contributes to the annihilation dynamics only for $p=2$ $-$ i.e. when
the second and third terms in Eq. (2) have the same differential order $-$ and
only in close proximity of annihilation (Fig. 3 inset).

Figure 4:  (a) Example of flow field in the surrounding of a $\pm 1/p$ pair
obtained from a numerical integration of Eqs. (1). (b) Trajectories of the
annihilating defects versus time. The annihilation time decreases
monotonically with $p$ and eventually plateaus for $p\to\infty$.

In addition to highlighting the origin of the hydrodynamic enhancement of pair
annihilation, Fig. 3 illustrates another surprising result: annihilation
occurs more rapidly as $p$ increases. To further clarify this phenomenon we
focus on the trajectories of the annihilating defects, with and without
backflow (Fig. 4a). In the absence of hydrodynamics, when annihilation is
solely driven by elastic forces, the trajectories are symmetric about the
$y-$axis and the defects annihilate at the origin (Fig. 4b). Conversely,
switching on hydrodynamics causes the positive defect to move faster toward
the positive $x-$direction and annihilation occurs in the half-plane that was
initially occupied by the negative defect. This effect is visibly more
pronounced in hexatics, thus confirming our previous observation that
annihilation occurs more rapidly as $p$ increases. Comparing the annihilation
trajectories of all $p$ values in the range $p=2,\,3\ldots\,6$, however, we
find that, even though the annihilation time decreases with $p$, it does so at
a decreasing rate and eventually plateaus as $p\to\infty$ (4b). This behavior
can be rationalized by requesting all forces experienced by each of the defect
to balance each other along the $x-$direction. Calling $x_{\pm}$ the
$x-$coordinate of the positive and negative defect and $v_{\pm}$ the
contribution to their velocity due to backflow, this gives

$\frac{{\rm d}x_{\pm}}{{\rm d}t}=v_{\pm}-\frac{ks^{2}}{x_{\pm}-x_{\mp}}\;.$
(6)

The first term on the right-hand side embodies the effective drag force acting
on the defects and originating from the dissipation of the defects elastic
energy Denniston (1996), whereas the last term arises from the elastic Coulomb
attraction between defects, with $k\sim\Gamma_{p}L_{p}$ a constant Giomi  _et
al._ (2013). From Eq. (6), taking $\Delta v=v_{+}-v_{-}\geq 0$, solving Eq.
(6) and imposing $x_{+}(t_{\rm a})=x_{-}(t_{\rm a})$ allows one to compute the
annihilation time $t_{\rm a}$ in the form

$t_{\rm a}=\frac{|\Delta x(0)|}{\Delta v}-\frac{2ks^{2}}{\Delta
v^{2}}\log\left[1+\frac{|\Delta x(0)|}{2ks^{2}}\,\Delta v\right]\;,$ (7)

where $\Delta x(0)=x_{+}(0)-x_{-}(0)$. Now, in the absence of backflow,
$\Delta v=0$ and Eq. (7) reduces to $t_{\rm a}=|\Delta x(0)|^{2}/(4ks^{2})$.
For finite $\Delta v$ values, on the other hand, $t_{\rm a}$ decreases
monotonically with $\Delta v$ and approaches $t_{\rm a}\approx|\Delta
x(0)|/\Delta v$ for large $\Delta v$ values. In turn, $\Delta v$ increase with
$p$ for small $p$ values SI , but vanishes for $p\to\infty$ when isotropy is
restored at the macroscopic scale and the defects themselves disappear: i.e.
$s=\pm 1/p\to 0$.

In conclusion, we demonstrated that hydrodynamic backflow profoundly affects
the dynamics of $p-$atic defects. In the case of isolated defects,
hydrodynamics can either stabilize defects, by trapping them at the stagnation
point of the resulting backflow, or propel them forward, by means of a
Stokeslet-type backflow consisting of two counterrotating vortices meeting
along the defect’s longitudinal direction. The latter phenomenon, however,
requires $s=(p-1)/p$, this is thermodynamically stable only in nematics (i.e.
$p=2$), where this coincides with the elementary winding number $s=1/2$. In
the case of annihilating defect pairs, we have established that hydrodynamic
effects become progressively more relevant as $p$ increases and we have
identified the precise origin of this effect. Our work paves the way for a
deeper understanding of cell intercalation and other remodelling events in
epithelial layers, where small-scale hexatic order (i.e. $p=6$) has been
recently discovered Armengol-Collado  _et al._ (2022a, b).

###### Acknowledgements.

This work is supported by the ERC-CoG grant HexaTissue and by Netherlands
Organization for Scientific Research (NWO/OCW). Part of this work was carried
out on the Dutch national e-infrastructure with the support of SURF through
the Grant 2021.028 for computational time. The authors acknowledge Ludwig
Hoffmann for fruitful discussions.

## References

  * Sanchez  _et al._ (2012) T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Nature 491, 431 (2012).
  * Keber  _et al._ (2014) F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp, L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, and A. R. Bausch, Science 345, 1135 (2014).
  * Giomi  _et al._ (2014a) L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek, and M. Cristina Marchetti, Philos. Trans. Royal Soc. A 372, 20130365 (2014a).
  * Giomi (2015) L. Giomi, Phys. Rev. X 5, 031003 (2015).
  * Doostmohammadi  _et al._ (2017) A. Doostmohammadi, T. N. Shendruk, K. Thijssen, and J. M. Yeomans, Nat. Commun. 8, 1 (2017).
  * Guillamat  _et al._ (2018) P. Guillamat, Ž. Kos, J. Hardoüin, J. Ignés-Mullol, M. Ravnik, and F. Sagués, Sci. Adv. 4, eaao1470 (2018).
  * Doostmohammadi  _et al._ (2018) A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Nat. Commun. 9, 1 (2018).
  * Carenza  _et al._ (2019a) L. N. Carenza, G. Gonnella, A. Lamura, G. Negro, and A. Tiribocchi, Eur. Phys. J. E 42, 1 (2019a).
  * Carenza  _et al._ (2019b) L. N. Carenza, G. Gonnella, D. Marenduzzo, and G. Negro, Proc. Nat. Acad. Sci. 116, 22065 (2019b).
  * Mermin (1979) N. D. Mermin, Rev. Mod. Phys. 51, 591 (1979).
  * Chaikin and Lubensky (1995) P. M. Chaikin and T. C. Lubensky,  _Principles of Condensed Matter Physics_ (Cambridge University Press, 1995).
  * Giomi  _et al._ (2017) L. Giomi, Žiga Kos, M. Ravnik, and A. Sengupta, Proc. Nat. Acad. Sci. 114, E5771 (2017).
  * Pollard  _et al._ (2019) J. Pollard, G. Posnjak, S. Čopar, I. Muševič, and G. P. Alexander, Phys. Rev. X 9, 021004 (2019).
  * Carenza  _et al._ (2022) L. N. Carenza, G. Gonnella, D. Marenduzzo, G. Negro, and E. Orlandini, Phys. Rev. Lett. 128, 027801 (2022).
  * Lopez-Leon and Fernandez-Nieves (2011) T. Lopez-Leon and A. Fernandez-Nieves, Colloid Polym Sci 289, 345 (2011).
  * Kibble (1976) T. W. B. Kibble, J. Phys. A: Math. Gen. 9, 1387 (1976).
  * Halperin and Nelson (1978) B. I. Halperin and D. R. Nelson, Phys. Rev. Lett. 41, 121 (1978).
  * Dirac (1931) P. A. M. Dirac, Proc. Math. Phys. 133, 60 (1931).
  * Hooft (1974) G. Hooft, Nucl. Phys. B. 79, 276 (1974).
  * Polyakov (1974) A. M. Polyakov, JETP Lett. 20, 194 (1974).
  * Saw  _et al._ (2017) T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu, S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeomans, and B. Ladoux, Nature 544, 212 (2017).
  * Monfared  _et al._ (2021) S. Monfared, G. Ravichandran, J. E. Andrade, and A. Doostmohammadi, arXiv:2108.07657 (2021).
  * Brugués  _et al._ (2014) A. Brugués, E. Anon, V. Conte, J. H. Veldhuis, M. Gupta, J. Colombelli, J. J. Muñoz, G. W. Brodland, B. Ladoux, and X. Trepat, Nat. Phys. 10, 683 (2014).
  * Kawaguchi  _et al._ (2017) K. Kawaguchi, R. Kageyama, and M. Sano, Nature 545, 327 (2017).
  * Livshits  _et al._ (2017) A. Livshits, L. Shani-Zerbib, Y. Maroudas-Sacks, E. Braun, and K. Keren, Cell Rep. 18, 1410 (2017).
  * Braun and Keren (2018) E. Braun and K. Keren, BioEssays 40, 1700204 (2018).
  * Livshits  _et al._ (2021) A. Livshits, L. Garion, Y. Maroudas-Sacks, L. Shani-Zerbib, K. Keren, and E. Braun, bioRxiv 2021.02.04.429818 (2021).
  * Maroudas-Sacks  _et al._ (2021) Y. Maroudas-Sacks, L. Garion, L. Shani-Zerbib, A. Livshits, E. Braun, and K. Keren, Nat. Phys. 17, 251 (2021).
  * Berreman (1975) D. W. Berreman, J. Appl. Phys. 46, 3746 (1975).
  * Brochard (1973) F. Brochard, Mol. Cryst. Liq. Cryst. 23, 51 (1973).
  * Kos Ž (2020) R. M. Kos Ž, Sci. Rep. 10, 10(1):1446 (2020).
  * Tsuji and Chono (2021) T. Tsuji and S. Chono, Sens. Actuator A Phys. 318, 112386 (2021).
  * Svenšek and Žumer (2001) D. Svenšek and S. Žumer, Liq. Cryst. 28, 1389 (2001).
  * Missaoui  _et al._ (2021) A. Missaoui, E. Lacaze, A. Eremin, and R. Stannarius, Crystals 11 (2021).
  * Brochard  _et al._ (1972) F. Brochard, P. Pieranski, and E. Guyon, Phys. Rev. Lett. 29, 320 (1972).
  * Tiribocchi  _et al._ (2014) A. Tiribocchi, O. Henrich, J. S. Lintuvuori, and D. Marenduzzo, Soft Matter 10, 4580 (2014).
  * Hatwalne  _et al._ (2004) Y. Hatwalne, S. Ramaswamy, M. Rao, and R. Simha, Phys. Rev. Lett. 92, 118101 (2004).
  * Tóth  _et al._ (2003) G. Tóth, C. Denniston, and J. Yeomans, Phys. Rev. E 67, 051705 (2003).
  * Martínez-Prat  _et al._ (2021) B. Martínez-Prat, R. Alert, F. Meng, J. Ignés-Mullol, J.-F. m. c. Joanny, J. Casademunt, R. Golestanian, and F. Sagués, Phys. Rev. X 11, 031065 (2021).
  * Colen  _et al._ (2021) J. Colen, M. Han, R. Zhang, S. A. Redford, L. M. Lemma, L. Morgan, P. V. Ruijgrok, R. Adkins, Z. Bryant, Z. Dogic, M. L. Gardel, J. J. de Pablo, and V. Vitelli, Proc. Nat. Acad. Sci. 118, e2016708118 (2021).
  * Dunkel  _et al._ (2013) J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär, and R. E. Goldstein, Phys. Rev. Lett. 110, 228102 (2013).
  * Wensink  _et al._ (2012) H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Proc. Natl. Acad. Sci. 109, 14308 (2012).
  * Balasubramaniam  _et al._ (2021) L. Balasubramaniam, A. Doostmohammadi, T. Saw, G. Narayana, R. Mueller, T. Dang, M. Thomas, S. Gupta, S. Sonam, A. Yap, Y. Toyama, R.-M. Mége, J. Yeomans, and B. Ladoux, Nat. Mater 20, 1156 (2021).
  * Metselaar  _et al._ (2019) L. Metselaar, J. M. Yeomans, and A. Doostmohammadi, Phys. Rev. Lett. 123, 208001 (2019).
  * Hoffmann  _et al._ (2022) L. A. Hoffmann, L. N. Carenza, J. Eckert, and L. Giomi, Sci. Adv. in press (2022).
  * Alert (2022) R. Alert, arXiv:2202.14027 (2022).
  * Giomi  _et al._ (2021a) L. Giomi, J. Toner, and N. Sarkar, arXiv:2111.04720 (2021a).
  * Giomi  _et al._ (2021b) L. Giomi, J. Toner, and N. Sarkar, arXiv:2106.11957 (2021b).
  * Zippelius (1980) A. Zippelius, Phys. Rev. A 22, 732 (1980).
  * Zippelius  _et al._ (1980) A. Zippelius, B. Halperin, and D. R. Nelson, Phys. Rev. B 22, 2514 (1980).
  * Oswald and Pieranski (2005) P. Oswald and P. Pieranski,  _Nematic and cholesteric liquid crystals: concepts and physical properties illustrated by experiments_ (CRC press, 2005).
  * Vromans and Giomi (2016) A. J. Vromans and L. Giomi, Soft Matter 12, 6490 (2016).
  * Giomi  _et al._ (2013) L. Giomi, M. J. Bowick, X. Ma, and M. C. Marchetti, Phys. Rev. Lett. 110, 228101 (2013).
  * Giomi  _et al._ (2014b) L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek, and M. Cristina Marchetti, Philos. Trans. Royal Soc. A 372, 20130365 (2014b).
  * (55) Supplementary information.
  * Kats  _et al._ (2002) E. Kats, V. Lebedev, and S. Malinin, J. Exp. Theor. Phys. 95, 714 (2002).
  * Svenšek and Žumer (2002) D. Svenšek and S. Žumer, Phys. Rev. E 66, 021712 (2002).
  * Missaoui  _et al._ (2020) A. Missaoui, K. Harth, P. Salamon, and R. Stannarius, Phys. Rev. Res. 2, 013080 (2020).
  * Kos and Ravnik (2020) Ž. Kos and M. Ravnik, Sci. Rep. 10, 1 (2020).
  * de Andrade Lima and Rey (2004) L. R. P. de Andrade Lima and A. D. Rey, Phys. Rev. E 70, 011701 (2004).
  * Carenza  _et al._ (2020a) L. N. Carenza, L. Biferale, and G. Gonnella, Eur. Phys. Lett. 132, 44003 (2020a).
  * Carenza  _et al._ (2020b) L. N. Carenza, L. Biferale, and G. Gonnella, Phys. Rev. Fluids 5, 011302 (2020b).
  * Denniston (1996) C. Denniston, Phys. Rev. B 54, 6272 (1996).
  * Armengol-Collado  _et al._ (2022a) J.-M. Armengol-Collado, Carenza, and L. Giomi, arXiv:2202.00651 (2022a).
  * Armengol-Collado  _et al._ (2022b) J.-M. Armengol-Collado, L. N. Carenza, J. Eckert, D. Krommydas, and L. Giomi, arXiv:2202.00668 (2022b).


11institutetext: Ronald R. Coifman 22institutetext: Department of Mathematics,
Program in Applied Mathematics, Yale University, New Haven, CT 06510, USA,
22email<EMAIL_ADDRESS>33institutetext: Jacques Peyrière
44institutetext: Institut de Mathématiques d’Orsay, CNRS, Université Paris-
Saclay, 91405 Orsay, France, 44email: jacques.peyriere@universite-paris-
saclay.fr

# Multiscale decompositions of Hardy spaces

Ronald R. Coifman and Jacques Peyrière

##  1 Introduction

We would like to elaborate on a program of analysis pursued by Alex Grossmann
and his collaborators on the analytic utilization of the phase of Hardy
functions, as a multiscale signal processing tool.

An inspiration at the origin of ”wavelet” analysis (when Grossmann, Morlet,
Meyer and collaborators were interacting and exploring versions of multiscale
representations) was provided, by the analysis of holomorphic signals, for
which, the images of the phase of Cauchy wavelets were remarkable in their
ability to reveal intricate singularities or dynamic structures, such as
instantaneous frequency jumps, in musical recordings. This work which was
pursued by Grossmann, Kronland Martinet et al sound exploiting phase and
amplitude variability of holomorphic signals was challenged by computational
complexity as well as by the lack of simple, efficient, mathematical
processing, and generalizations to higher dimensional signals. It was mostly
bypassed by the orthogonal wavelet transforms. We aim to show that these ideas
are powerful nonlinear subtle tools.

Our goal here is to follow their seminal work and introduce recent
developments in nonlinear analysis. In particular we will sketch methods
extending conventional Fourier analysis, exploiting both phase and amplitudes
of holomorphic functions.

The miracles of nonlinear complex analysis, such as factorization and
composition of functions lead to new versions of holomorphic wavelets, and
relate them to multiscale dynamical systems.

Our story interlaces the role of the phase of signals with their
analytic/geometric properties. The Blaschke factors are a key ingredient, in
building analytic tools, starting with the Malmquist Takenaka orthonormal
bases of the Hardy space $\mathsf{H}^{2}({\mathbb{T}})$, continuing with
”best” adapted bases obtained through phase unwinding, and concluding with
relations to composition of Blaschke products and their dynamics (on the disc,
and on invariant subspaces of $\mathsf{H}^{2}({\mathbb{T}})$). Specifically we
discuss multiscale orthonormal holomorphic wavelet bases, related to
Grossmann’s and Morlet’s program wavelet, and associated generalized scaled
holomorphic orthogonal bases, to dynamical systems, obtained by composing
Blaschke factors.

We also, remark, that the phase of a Blaschke product is a one layer neural
net with ($\arctan$ as an activation sigmoid) and that the composition is a
”Deep Neural Net” whose depth is the number of compositions, our results
provide a wealth of related libraries of orthogonal bases .

We sketch these ideas in various ”vignette” subsections and refer for more
details on analytic methods CP, related to the Blaschke based nonlinear phase
unwinding decompositions coifman; CSW; nahon, we also consider orthogonal
decompositions of invariant subspaces of Hardy spaces. In particular we
constructed a multiscale decomposition, described below, of the Hardy space of
the upper half-plane.

Such a decomposition can be carried in the unit disk by conformal mapping. A
somewhat different multiscale decomposition of the space
$\mathsf{H}^{2}({\mathbb{T}})$ has been constructed by using Malmquist-
Takenaka bases associated with Blaschke products whose zéros are
$\displaystyle(1-2^{-n})\mathrm{e}^{2\mathrm{i}\pi j/2^{n}}$ where $n\geq 1$
and $0\leq j<2^{n}$ feichtinger. Here we provide a variety of multiscale
decompositions by considering iterations of Blaschke products.

##  2 Preliminaries and notation

For $p\geq 1$, $\mathsf{H}^{p}({\mathbb{T}})$ stands for the space of analytic
functions $f$ on the unit disk ${\mathbb{D}}$ such that

$\sup_{0<r<1}\int_{0}^{2\pi}|f(r\mathrm{e}^{\mathrm{i}\theta})|^{p}\frac{\mathrm{d}\theta}{2\pi}<+\infty.$

Such functions have boundary values almost everywhere, and the Hardy space
$\mathsf{H}^{p}({\mathbb{T}})$ can be identified with the set of $L^{p}$
functions on the torus ${\mathbb{T}}=\partial{\mathbb{D}}$ whose Fourier
coefficients of negative order vanish.

A subspace of $\mathsf{H}^{p}({\mathbb{T}})$ is _invariant_ if it is invariant
under multiplication by $\mathrm{e}^{\mathrm{i}\theta}$ (or by $z$), depending
whether these functions are considered as functions on ${\mathbb{T}}$ or
${\mathbb{D}}$. An inner function is a bounded analytic function on the unit
disk whose boundary values have modulus 1 almost everywhere. It is known that
the invariant subspaces are of the form $u\mathsf{H}^{p}({\mathbb{T}})$ where
$u$ is an inner function helson; hoffman. The inner function $u$ is determined
by the invariant subspace up to multiplication by a constant of modulus 1.

If $f$ and $g$ are two functions on ${\mathbb{T}}$ (in $L^{p}$ and
$L^{p/(p-1)}$ for some $p\in[1,+\infty)$), let

$\left\langle
f,g\right\rangle=\frac{1}{2\pi}\int_{0}^{2\pi}f(\mathrm{e}^{\mathrm{i}\theta})\overline{g(\mathrm{e}^{\mathrm{i}\theta})}\,\mathrm{d}\theta.$

Let $\mathscr{H}$ be the operator of orthogonal projection of
$L^{2}({\mathbb{T}})$ onto $\mathsf{H}^{2}({\mathbb{T}})$. It results from the
properties of the Hilbert transform that this operator extends as a bounded
operator from $L^{p}({\mathbb{T}})$ to $\mathsf{H}^{p}({\mathbb{T}})$ for
$1<p<+\infty$.

If $u$ is an inner function, let $\chi_{u}$ be the operator of multiplication
by $u$ (which is an isometry of all the $L^{P}$). Then the operator
${\mathscr{H}}_{u}=\chi_{u}{\mathscr{H}}\chi_{u}^{-1}$ is the operator of
orthogonal projection of $L^{2}$ onto $u\mathsf{H}^{2}({\mathbb{T}})$. It
results that this operator extends as a bounded operator from
$L^{p}({\mathbb{T}})$ to $\mathsf{H}^{p}({\mathbb{T}})$ for all
$p\in(1,+\infty)$ with a norm independent of $u$. In other terms, for all
$p>1$, there exists $C_{p}$ such that, for all $u$ and all $f\in
L^{p}({\mathbb{T}})$,

$\|{\mathscr{H}}_{u}f\|_{p}\leq C_{p}\|f\|_{p}.$ (1)

There is a parallel theory for analytic functions on the upper half plane
${\mathbb{H}}=\\{x+\mathrm{i}y\ :\ y>0\\}$. The space of analytic functions
$f$ on ${\mathbb{H}}$ such that

$\sup_{y>0}\|f(\cdot+\mathrm{i}y)\|_{L^{p}({\mathbb{R}})}<+\infty$

is denoted by $\mathsf{H}^{p}({\mathbb{R}})$. These functions have boundary
values in $L^{p}({\mathbb{R}})$ when $p\geq 1$. The space
$\mathsf{H}^{p}({\mathbb{R}})$ is identified to the space of $L^{p}$ functions
whose Fourier transform vanishes on the negative half line $(-\infty,0)$.

A subspace of $\mathsf{H}^{2}({\mathbb{R}})$ is said to be invariant if it is
stable by multiplication by the functions $\mathrm{e}^{2\mathrm{i}\pi\xi x}$
for all $\xi>0$. As previously, the invariant subspaces are of the form
$u\,\mathsf{H}^{2}$ where $u$ is an inner function, i.e., a bounded analytic
function on ${\mathbb{H}}$ whose boundary values are of modulus 1 almost
everywhere.

As previously, the operators of orthogonal projections on invariant subspaces
extend, for any $p\in(1,+\infty)$, as continuous operators on
$\mathsf{H}^{p}({\mathbb{R}})$ with a uniform bound for their norms.

##  3 Malmquist-Takenaka bases on the torus

###### Lemma 1

Let $a$ be a complex number of modulus less than 1 and $u$ be an inner
function. Then $(z-a)u\mathsf{H}^{2}$ has codimension 1 in $u\mathsf{H}^{2}$
and $\displaystyle\frac{\sqrt{1-|a|^{2}}}{1-\overline{a}z}u$ is a unit vector
in the orthogonal complement of $(z-a)u\mathsf{H}^{2}$ in $u\mathsf{H}^{2}$.

###### Proof

Since $f\mapsto uf$ is an isometry of $\mathsf{H}^{2}$ onto $u\mathsf{H}^{2}$,
it is enough to consider the case $u=1$. One has

$\displaystyle\left\langle(z-a)f(z),\frac{\sqrt{1-|a|^{2}}}{1-\overline{a}z}\right\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}(\mathrm{e}^{\mathrm{i}\theta}-a)f(\mathrm{e}^{\mathrm{i}\theta})\frac{\sqrt{1-|a|^{2}}}{1-a\mathrm{e}^{-\mathrm{i}\theta}}\,\mathrm{d}\theta$
$\displaystyle=$
$\displaystyle\frac{\sqrt{1-|a|^{2}}}{2\pi}\int_{-\pi}^{\pi}\mathrm{e}^{\mathrm{i}\theta}f(\mathrm{e}^{\mathrm{i}\theta})\,\mathrm{d}\theta=0.$

Also, if $f$ is orthogonal to $(1-\overline{a}z)^{-1}$ one has

$0=\frac{1}{2\pi}\int_{-\pi}^{\pi}\frac{f(\mathrm{e}^{\mathrm{i}\theta})}{1-a\mathrm{e}^{-\mathrm{i}\theta}}\,\mathrm{d}\theta\\\
=\frac{1}{2\mathrm{i}\pi}\oint\frac{f(z)}{z-a}\,\mathrm{d}z=f(a),$

so $f\in(z-a)\mathsf{H}^{2}$.

Now let $(a_{n})_{n>0}$ be a sequence of complex numbers of modulus less than
1. For $n\geq 0$, let

${\mathsf{B}}_{n}(z)=\prod_{0\leq
j<n}\frac{z-a_{j}}{1-\overline{a}_{j}z}\text{\quad
and\quad}\phi_{n}(z)={\mathsf{B}}_{n}(z)\frac{\sqrt{1-|a_{n}|^{2}}}{1-\overline{a}_{n}z}.$

It results from Lemma 1 that, if

$\sum_{n\geq 1}(1-|a_{j}|^{2})=+\infty,$

the functions $\phi_{n}$ form an orthonormal basis of $\mathsf{H}^{2}$.

If $\displaystyle\sum_{n\geq 1}(1-|a_{j}|^{2})<+\infty$ the functions
$\phi_{n}$ form an orthonormal basis of
$\mathsf{H}^{2}\ominus{\mathsf{B}}\mathsf{H}^{2}$, where ${\mathsf{B}}$ is the
convergent Blaschke product

${\mathsf{B}}(z)=\prod_{j>0}\frac{\overline{a}_{j}}{|a_{j}|}\frac{z-a_{j}}{1-\overline{a}_{j}z}.$

Consider a sequence $({\mathsf{B}}_{m})_{m\geq 1}$:

${\mathsf{B}}_{m}(z)=\prod_{j>0}\frac{\overline{a}_{m,j}}{|a_{m,j}|}\frac{z-a_{m,j}}{1-\overline{a}_{m,j}z}$

of convergent Blaschke products such that
$\displaystyle\sum_{m,j>0}(1-|a_{m,j}|^{2})=+\infty$.

Let ${\mathsf{B}}_{0}=1$ and, for $n\geq 0$ and $m\geq 1$,

$\displaystyle{\mathsf{B}}_{m,n}(z)$ $\displaystyle=$
$\displaystyle\prod_{0\leq j<n}\frac{z-a_{m,j}}{1-\overline{a}_{m,j}z}$
$\displaystyle\phi_{m,n}(z)$ $\displaystyle=$
$\displaystyle{\mathsf{B}}_{m-1}(z){\mathsf{B}}_{m,n}(z)\frac{\sqrt{1-|a_{m,n}|^{2}}}{1-\overline{a}_{m,n}z}.$

Then, $(\phi_{m,n})_{n\geq 1}$ is an orthonormal basis of
${\mathsf{B}}_{m-1}\mathsf{H}^{2}\ominus{\mathsf{B}}_{m-1}{\mathsf{B}}_{m}\mathsf{H}^{2}$,
and $(\phi_{m,n})_{m\geq 1,n\geq 1}$ is an orthonormal basis of
$\mathsf{H}^{2}$.

The bases so obtained are the _Malmquist-Takenaka bases_ takenaka.

##  4 The upper half plane

We present some prior results CP, without proof. In this section one simply
writes $\mathsf{H}^{2}$ instead of $\mathsf{H}^{2}({\mathbb{R}})$.

###  4.1 Malmquist-Takenaka bases

Let $(a_{j})_{1\leq j}$ be a sequence (finite or not)) of complex numbers with
positive imaginary parts and such that

$\sum_{j\geq 0}\frac{\Im a_{j}}{1+|a_{j}|^{2}}<+\infty.$ (2)

The corresponding Blaschke product is

${\mathsf{B}}(x)=\prod_{j\geq
0}\frac{\left|1+a_{j}^{2}\right|}{1+a_{j}^{2}}\,\frac{x-a_{j}}{x-\overline{a}_{j}},$

where, $0/0$, which appears if $a_{j}=\mathrm{i}$, should be understood as 1.
The factors $\displaystyle\frac{\left|1+a_{j}^{2}\right|}{1+a_{j}^{2}}$ insure
the convergence of this product when there are infinitely many zeroes. But, in
some situations, it is more convenient to use other convergence factors as we
shall see below.

Whether the series (2) is convergent or not, one defines (for $n\geq 0$) the
functions

$\phi_{n}(x)=\frac{1}{\sqrt{\pi}}\left(\prod_{0\leq
j<n}\frac{x-a_{j}}{x-\overline{a}_{j}}\right)\,\frac{1}{x-\overline{a}_{n}}.$

Then these functions form a orthonormal system in $\mathsf{H}^{2}$. If the
series (2) diverges, it is a basis of $\mathsf{H}^{2}$, otherwise it is a
basis of the orthogonal complement of ${\mathsf{B}}\,\mathsf{H}^{2}$ in
$\mathsf{H}^{2}$.

###  4.2 A multiscale Wavelet decomposition

The infinite products

$G_{n}(x)=\prod_{j\leq
n}\frac{j-\mathrm{i}}{j+\mathrm{i}}\,\frac{x-j-\mathrm{i}}{x-j+\mathrm{i}}\text{\quad
and\quad}G(x)=\prod_{j\in{\mathbb{Z}}}\frac{j-\mathrm{i}}{j+\mathrm{i}}\,\frac{x-j-\mathrm{i}}{x-j+\mathrm{i}}$

can be expressed in terms of known functions:

$G_{n}(x)=\frac{\Gamma(-\mathrm{i}-n)}{\Gamma(\mathrm{i}-n)}\,\frac{\Gamma(x-n+\mathrm{i})}{\Gamma(x-n-\mathrm{i})}\text{\quad
and\quad}G(x)=\frac{\sin\pi(\mathrm{i}-x)}{\sin\pi(\mathrm{i}+x)}.$

###  4.3 An orthonormal system

Consider the function

$\phi(x)=\frac{\Gamma(x-1+\mathrm{i})}{\sqrt{\pi}\Gamma(x-\mathrm{i})}.$

It is easily checked that

$\phi(x-n)=\frac{\Gamma(\mathrm{i}-n)}{\Gamma(-\mathrm{i}-n)}\,\frac{G_{n}(x)}{\sqrt{\pi}\bigl{(}x-(n+1)+\mathrm{i}\bigr{)}}.$

Set $\phi_{n}(x)=\phi(x-n)$. For fixed $m$, the functions $\phi_{n}/G_{m}$,
for $n\geq m$, form a Malmquist-Takenaka basis of $(G/G_{m})\mathsf{H}^{2}$.
In other terms, the functions $\phi_{n}$, for $n\geq m$, form an orthonormal
basis of $G_{m}\mathsf{H}^{2}\ominus G\mathsf{H}^{2}$. This means that the
functions $\phi_{n}$ (for $n\in{\mathbb{Z}}$) form a Malmquist-Takenaka basis
of the orthogonal complement of $G\mathsf{H}^{2}$ in $\mathsf{H}^{2}$.

#### Multiscale decomposition

As $|1-G(2^{n}x)|\leq C2^{n}$ all the products

${\mathscr{B}}_{n}(x)=\prod_{j<n}G(2^{j}x)$ (3)

are convergent and $\displaystyle\lim_{n\to-\infty}{\mathscr{B}}_{n}=1$
uniformly.

Let ${\mathscr{B}}={\mathscr{B}}_{0}$. Obviously,
${\mathscr{B}}_{n}(x)={\mathscr{B}}(2^{n}x)$.

Consider the following subspaces of $\mathsf{H}^{2}$:

${\mathsf{E}}_{n}={\mathscr{B}}_{n}\mathsf{H}^{2}.$

This is a decreasing sequence. The space
$\displaystyle{\mathsf{E}}_{+\infty}=\bigcap_{n\in{\mathbb{Z}}}{\mathsf{E}}_{n}$
is equal to $\\{0\\}$ since a function orthogonal to this space would have too
many zeros, and the space
$\displaystyle{\mathsf{E}}_{-\infty}=\mathrm{closure~{}of}\bigcup_{n\in{\mathbb{Z}}}{\mathsf{E}}_{n}$
is equal to $\mathsf{H}^{2}$ since ${\mathscr{B}}_{n}$ converges uniformly to
1 when $n$ goes to $-\infty$.

For all $n$ and $j$, let

$\phi_{n,j}(x)=2^{n/2}\phi(2^{n}x-j){\mathscr{B}}(2^{n}x).$

Then, for all $n$, $(\phi_{n,j})_{j\in{\mathbb{Z}}}$ is an orthonormal basis
of ${\mathsf{E}}_{n}\ominus{\mathsf{E}}_{n+1}$. At last
$(\phi_{n,j})_{n,j\in{\mathbb{Z}}}$ is an orthonormal basis of
$\mathsf{H}^{2}$.

##  5  Adapted MT bases, ”phase unwinding”

Our goal is to find a ”best” adapted Malmquist Takenaka basis to analyze a
given function the idea is to peel off the oscillation of a function by
dividing by its Blaschke product, this procedure is iterated to yield an
expansion in an orthogonal collection of functions or Blaschke products which
of course are naturally embedded in a MT basis, once the zeroes are ordered.

###  5.1 The unwinding series.

There is a natural way to iterate the Blaschke factorization, it is inspired
by the power series expansion of a holomorphic function. If $G$ has no zeroes
inside $\mathbb{D}$, its Blaschke factorization is the trivial one $G=1\cdot
G$, however, the function $G(z)-G(0)$ certainly has at least one root inside
the unit disk $\mathbb{D}$ and will therefore yield some nontrivial Blaschke
factorization $G(z)-G(0)={\mathsf{B}}_{1}G_{1}$. We write

$\displaystyle F$ $\displaystyle={\mathsf{B}}\cdot G$
$\displaystyle={\mathsf{B}}\cdot(G(0)+(G(z)-G(0)))$
$\displaystyle={\mathsf{B}}\cdot(G(0)+{\mathsf{B}}_{1}G_{1})$
$\displaystyle=G(0){\mathsf{B}}+{\mathsf{B}}{\mathsf{B}}_{1}G_{1}.$

An iterative application gives rise to the unwinding series

$F=a_{1}{\mathsf{B}}_{1}+a_{2}{\mathsf{B}}_{1}{\mathsf{B}}_{2}+a_{3}{\mathsf{B}}_{1}{\mathsf{B}}_{2}{\mathsf{B}}_{3}+a_{4}{\mathsf{B}}_{1}{\mathsf{B}}_{2}{\mathsf{B}}_{3}{\mathsf{B}}_{4}+\dots$

This orthogonal expansion first appeared in the PhD thesis of Michel Nahon
nahon and independently by T. Qian in qtao; qw. Given a general function $F$
it is not numerically feasible to actually compute the roots of the function;
a crucial insight in nahon is that this is not necessary – one can numerically
obtain the Blaschke product in a stable way by using a method that was first
mentioned in a paper of Guido and Mary Weiss ww and has been investigated with
respect to stability by Nahon nahon Using the boundedness of the Hilbert
transform one can prove easily convergence in $L^{p},1<p<\infty$ .

###  5.2  The fast algorithm of Guido and Mary Weiss ww

Our starting point is the theorem that any Hardy function can be decomposed as

$F={\mathsf{B}}\cdot G,$

where $B$ is a Blaschke product, that is a function of the form

${\mathsf{B}}(z)=z^{m}\prod_{i\in
I}{\frac{\overline{a_{i}}}{|a_{i}|}\frac{z-a_{i}}{1-\overline{a_{i}}z}},$

where $m\in\mathbb{N}_{0}$ and $a_{1},a_{2},\dots\in\mathbb{D}$ are zeroes
inside the unit disk $\mathbb{D}$ and $G$ has no roots in $\mathbb{D}$. For
$|z|=1$ we have $|{\mathsf{B}}(z)|=1$ which motivates the analogy

${\mathsf{B}}\sim\mbox{frequency and}~{}G\sim\mbox{amplitude}$

for the function restricted to the boundary. However, the function $G$ need
not be constant: it can be any function that never vanishes inside the unit
disk. If $F$ has roots inside the unit disk, then the Blaschke factorization
$F={\mathsf{B}}\cdot G$ is going to be nontrivial (meaning
${\mathsf{B}}\not\equiv 1$ and $G\not\equiv F$). $G$ should be ’simpler’ than
$F$ because the winding number around the origin decreases.

In fact since $|F|=|G|$ and $\ln(G)$ is analytic in the disk we have formally
that

$G=\exp(\ln|F|+\mathrm{i}(\ln|F|)^{\sim}))=\exp(\mathscr{H}(\ln|F|))$

and ${\mathsf{B}}=F/G$

G can be computed easily using the FFT nahon.

##  6 Iteration of Blaschke products.

We are interested in the case we iterate finite Blaschke products:

${\mathsf{B}}(z)=\mathrm{e}^{\mathrm{i}\theta}z^{\nu}\prod_{j=1}^{\mu}\frac{z+a_{j}}{1+\overline{a}_{j}z},$

where $\mu$ and $\nu$ are nonnegative integers and the $a_{j}$ are complex
numbers of modulus less than 1.

It is well known that ${\mathbb{T}}$ and ${\mathbb{D}}$ are globally invariant
under ${\mathsf{B}}$, as well as the complement of $\overline{\mathbb{D}}$ in
the Riemann sphere.

We have

$\mathrm{e}^{-\mathrm{i}\varphi}{\mathsf{B}}(\mathrm{e}^{\mathrm{i}\varphi}z)=\mathrm{e}^{\mathrm{i}\bigl{(}\theta+(\nu+\mu-1)\varphi\bigr{)}}z^{\nu}\prod_{j=1}^{\mu}\frac{z+b_{j}}{1+\overline{b}_{j}z},$

where $b_{j}=\mathrm{e}^{-\mathrm{i}\varphi}a_{j}$. This means that for
iteration purpose we may assume $\theta=0$.

When $\nu\geq 1$, it results from the Schwarz lemma that $0$ is the unique
fixed point of ${\mathsf{B}}$ in ${\mathbb{D}}$ and that this fixed point is
attracting. The basin of attraction of $0$ contains a disk centered at $0$.
The sequence of iterates ${\mathsf{B}}_{n}$ is a normal family, since it
converges towards 0 on some neighborhood of $0$, it has a unique limit point.
This means that ${\mathsf{B}}_{n}$ converges to $0$ uniformly on any compact
subset of ${\mathbb{D}}$.

Now, suppose $\nu=0$. We have

${\mathsf{B}}(z)\overline{{\mathsf{B}}(1/\overline{z})}=1,$

so, if $\alpha$ is a fixed point of ${\mathsf{B}}$, so is
$1/\overline{\alpha}$.

we have two possibilities:

  1. 1.

There exists a fixed point $\alpha\in{\mathbb{D}}$ of ${\mathsf{B}}$. Then if
$\phi(z)=\frac{z+\alpha}{1+\overline{\alpha}z}$, then $0$ is a fixed point of
the Blaschke product $\phi^{-1}\circ{\mathsf{B}}\circ\phi$. It results from
the preceding discussion that $\alpha$ is the unique fixed point of
${\mathsf{B}}$ in ${\mathbb{D}}$, that it is attracting, and that the sequence
of iterates of ${\mathsf{B}}$ converges to the constant $\alpha$.

  2. 2.

All the fixed points lie in ${\mathbb{T}}$. Then according to the Denjoy-Wolff
theorem carleson, again the sequence of iterates of ${\mathsf{B}}$ converges
to some constant.

Moreover, the fixed points of ${\mathsf{B}}$ are the roots of an equation of
degree $\mu+1$, so if ${\mathsf{B}}$ has a fixed point in ${\mathbb{D}}$, it
has $\mu-1$ fixed points in $\mathbb{T}$ (none of them being attracting). In
this case, the dynamics on $\mathbb{T}$ is that of a cookie cutter.

When ${\mathsf{B}}_{1}$ and ${\mathsf{B}}_{2}$ are finite Blaschke products
with $n_{1}$ and $n_{2}$ zeros, then ${\mathsf{B}}_{2}\circ{\mathsf{B}}_{1}$
is a finite Blaschke products with $n_{1}n_{2}$ zeros (of course zeros are
counted with their multiplicities). This is obvious by considering the
variation of arguments on the boundary of the disk.

Let ${\mathsf{B}}$ be a finite Blaschke product having at least two zeros.
Then, one may consider the dynamical system which it generates. Let
${\mathsf{B}}_{0}(z)=z$, and
${\mathsf{B}}_{n+1}={\mathsf{B}}_{n}\circ{\mathsf{B}}$.

###### Lemma 2

Consider a Blaschke product $F$ of the form $F(z)=z{\mathsf{B}}(z)$, where
${\mathsf{B}}$ is nonconstant. Then there exists a sequence
$0,a_{1},a_{2},\dots,a_{j},\dots$ of complex numbers in the unit disk and an
increasing sequence $(\nu_{j})_{j\geq 1}$ of positive integers such that
$a_{1},a_{2},\dots,a_{\nu_{n}}$ are the zeros, counted according to their
multiplicity, of $F_{n}$. Moreover $\displaystyle\sum_{j\geq
1}(1-|a_{j}|)=+\infty$.


# Improving adaptability to new environments and removing catastrophic
forgetting in Reinforcement Learning by using an eco-system of agents

Olivier Moulin1 Vincent Francois-Lavet1 Paul Elbers2 Mark Hoogendoorn1 1Vrije
Universiteit Amsterdam
2Amsterdam UMC

###### Abstract

Adapting a Reinforcement Learning (RL) agent to an unseen environment is a
difficult task due to typical over-fitting on the training environment. RL
agents are often capable of solving environments very close to the trained
environment, but when environments become substantially different, their
performance quickly drops. When agents are retrained on new environments, a
second issue arises: there is a risk of catastrophic forgetting, where the
performance on previously seen environments is seriously hampered. This paper
proposes a novel approach that exploits an ecosystem of agents to address both
concerns. Hereby, the (limited) adaptive power of individual agents is
harvested to build a highly adaptive ecosystem. This allows to transfer part
of the workload from learning to inference. An evaluation of the approach on
two distinct distributions of environments shows that our approach outperforms
state-of-the-art techniques in terms of adaptability/generalization as well as
avoids catastrophic forgetting.

##  1 Introduction

In Reinforcement Learning, over-fitting to the environment on which the agent
has been trained is a common problem (cf. (KCobbe2019, , Cobbe et al., 2019)),
resulting in poor performance in unseen environments that are substantially
different from the one trained upon.

EnvironmentTask learnedAdaptability Figure 1: Our approach: eco-system
adaptability concept

This is particularly challenging in cases where agents are regularly
confronted with new environments. Adaptability to never seen or modified
environments is commonly referred to as the generalization capability of a
given approach (cf. (CPacker2018, , Packer et al., 2018)). A variety of
approaches have been introduced to address this challenge. The approaches
proposed in the literature often focus on a single agent being trained on
several environments and trying to learn a higher-level strategy to solve this
set of environments (see e.g. (asonar2020, , Sonar et al., 2018), (Migl2019, ,
Igl et al., 2019)). These approaches are typically shown to have better
performance on multiple environments. While such existing approaches certainly
improve performance on generalizability, they have several important
downsides. First of all, they typically suffer from so-called catastrophic
forgetting, which means that the agent looses the capability to perform well
in environments it has explored in details / solved in the past. This issue
has been detailed and solutions have been proposed by continuing to sample
experiences from past environments when training on new ones (cf. in
(catkinson2020, , Atkinson et al., 2020)), but still no guarantees can be
given on not forgetting anything which might be crucial in critical settings.
The second downside is that retraining on new environments can be very
expensive, especially in cases where limited computational resources are
available.

In this paper, we propose an approach based on the observation that agents
trained on an environment can generalize to environments that are close to the
original training environment (Figure 1). We combine a set of such agents in
an ecosystem, where each agent is trained on one specific environment. This
means that we have a whole set of agents with limited capacity to generalize,
which through combining them grows into an ecosystem that can generalize to a
rich set of environments. The ecosystem expands by adding new agents when new
environments are encountered that cannot be solved by the current ecosystem,
which means that each agent in the pool is never retrained after being added
to the pool. We evaluate the approach in two different environments, one being
the well-known minigrid framework and the other based on a newly defined
setting. We also benchmark our approach against an existing approach (cf.
(Migl2019, , Igl et al., 2019)) for trying to increase adaptability /
generalization and compare to a simple voting ensemble technique (cf.
(MAWiering2008, , Wiering et al., 2008)).

This paper is organized as follows. Section 2 presents the related work.
Section 3 follows by an explanation of our approach. Next, Section 4 presents
the experimental setup used to evaluate the approach. Section 5 presents the
results of the experiments and compares our approach with the chosen
benchmarks. We end with a discussion in Section 6.

##  2 Related work

###  2.1 Generalization

A significant number of papers have looked at generalization or adaptability
of a reinforcement learning agent to never seen environments. Agent
adaptability can be seen as a specific generalization capability. It has been
defined by Packer et al. CPacker2018  and assessed on a set of environments
and RL algorithms by Cobbe et al. KCobbe2019 . Two sets of approaches are
essentially proposed to improve generalizability which are described below.

The first category of approaches are focused on trying to create a
representation of the environment to learn a good policy with a more limited
risk of over-fitting. An example of this work can be found in the paper from
Sonar et al. asonar2020  where the researchers propose a new algorithm in
which the agent learns a representation on which the action predictor is
optimal on all the training environments. These approaches, although they
improve the adaptability of the agent by reducing over-fitting to an
environment, still suffer from catastrophic forgetting.

The second category of approaches are focused on adding noise and/or
information bottlenecks in the learner (e.g. neural network) to avoid over-
fitting to the environment. This approach is for instance used by Chen
jzchen2020  by implementing a Surprise Minimization algorithm, which is an
algorithm that provides an estimated reward based on probability according to
the history of rewards by state. A similar approach is used by Lu et al.
Xlu2020  and Igl et al. Migl2019  where the researchers propose to add an
information bottleneck between the environment and the agent as well as
including noise to avoid over-fitting. These types of approaches, have proven
to be successful on reducing over-fitting, but do not solve the catastrophic
forgetting problem.

###  2.2 Catastrophic forgetting

To address catastrophic forgetting, most of the literature focuses on three
alternatives: (1) making sure the agent retrains on previous tasks (e.g.
(catkinson2020, , Atkinson et al., 2020)); (2) modifying weights used for a
past task only slowly (cf. (JKickpatrick2016, , Kickpatrick et al., 2016)), or
(3) keeping several versions of the neural network with a progressive network
architecture (e.g. (AARusu2016, , Rusu et al., 2016)). These approaches
significantly reduce catastrophic forgetting but cannot completely avoid it,
which is the case for our proposed approach.

##  3 Approach

Our eco-system approach is based on the assumption that each agent trained on
an environment embeds some adaptability, which allows it to perform properly
on a limited number of other environments. The eco-system is composed of
multiple agents (a pool of agents). Each agent part of the eco-system is
trained as a standard Reinforcement Learning agent.

###  3.1 Reinforcement Learning formulation

We consider an agent interacting with its environment over discrete time
steps. The environment is formalized as an MDP defined by (i) a state space
$\mathcal{S}$ that is possibly continuous, (ii) an action space
$\mathcal{A}=\\{1,\ldots,N_{\mathcal{A}}\\}$, (iii) the transition function
$T:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\to[0,1]$, (iv) the reward
function $R:\mathcal{S}\times\mathcal{A}\times\mathcal{S}\to\mathcal{R}$ where
$\mathcal{R}$ is a continuous set of possible rewards in a range
$R_{\text{max}}\in\mathbb{R}^{+}$ ($[0,R_{\text{max}}]$). An MDP $M$ starts in
a distribution of initial states $b_{0}(s)$. At time step $t$, the agent
chooses an action based on the state of the system $s_{t}\in\mathcal{S}$
according to a policy $\pi:\mathcal{S}\times\mathcal{A}\rightarrow[0,1]$.
After taking action $a_{t}\sim\pi(s_{t},\cdot)$, the agent then observes a new
state $s_{t+1}\in\mathcal{S}$ as well as a reward signal
$r_{t}\in\mathcal{R}$.

####  3.1.1 The meta-learning setting

In this paper, we do not consider the task of learning in a single environment
but a distribution of different (yet related) environments $M\sim\mathcal{M}$
that differ in the reachable states from the state space $\mathcal{S}$ while
the action space $\mathcal{A}$, the reward function $R$ and the transition
function $T$ are the same (see Section 4 for two examples). The agent aims at
finding a policy $\pi(s,a)$ with the objective of maximizing its expected
return, defined (in the discounted setting) as

$V^{\pi}(b_{0})=\underset{M_{i}\sim\mathcal{M}}{\mathbb{E}}\left[\sum\nolimits_{k=0}^{\infty}\gamma^{k}r_{t+k}\mid
s_{t}=s,\pi\right].$

where

$r_{t}=\underset{a\sim\pi(s_{t},\cdot)}{\mathbb{E}}R\big{(}s_{t},a,s_{t+1}\big{)}$

and

$\mathbb{P}\big{(}s_{t+1}|s_{t},a_{t}\big{)}=T(s_{t},a_{t},s_{t+1})$

. In our approach, we consider that an agent has ”solved an environment” when
the total reward collected $\mathcal{R}$ is bigger than a predefined threshold
$l$, even if it is not the optimal total reward. This threshold is defined to
be what is considered a good enough policy to solve an environment, this
changes based on the environment presented to the system.
The Double Deep Q Network (DDQN) algorithm and the Proximal Policy
Optimization algorithm (PPO) have been used for running the experiments
(detailed in 4). Using these two different algorithms showed that the
improvements brought by the eco-system approach is independent from the
optimization algorithm used.

####  3.1.2 Q-learning

In Reinforcement Learning, the agent learns a policy (sequence of actions)
which allows it to get the best reward in its environment. In order to
leverage all past experiences, we use an off-policy learning algorithm that
samples transition tuples $(s,a,r,\gamma,s^{\prime})$ from a replay buffer.
The Q-function is learned using the Double Deep Q-Network (DDQN) algorithm
(cf. (HVanHasselt2019, , Van Hasselt et al., 2019)), which uses the target:

$Y=r+\gamma
Q(s^{\prime},\operatorname*{argmax}_{a^{\prime}\in\mathcal{A}}Q(s^{\prime},a^{\prime};\theta_{Q});\theta_{Q^{-}})$

where $\theta_{Q^{-}}$ are parameters of an earlier buffered Q-function (or
our target Q-function). The agent then minimizes the following loss:

$L_{Q}(\theta_{Q})=(Q(s,a;\theta_{Q})-Y)^{2}$

From this Q-value function we can derive the policy $\pi$ learned by the agent
as

$\pi(s)=\underset{a\in\mathcal{A}}{\operatorname*{argmax}}\,Q(s,a;\theta_{Q}),\forall
s\in\mathcal{S}$

The specific version of DDQN used in this paper is based on DeeR
implementation (VFrancoisLavet2016, , Francois-Lavet, 2016).

####  3.1.3 PPO

The Proximal Policy Optimization (PPO) method (cf. (jschulman2017, , Schulman
et al., 2017)) is an extension of the actor-critic method where the parameters
$w$ of the policy $\pi_{w}(s,a)$ are updated to optimize
$A^{\pi_{w}}(s,a)=Q^{\pi_{w}}(s,a)-V^{\pi_{w}}(s)$. The PPO algorithm also
enforces a limit on the policy changes, which results in maximizing the
following objective in expectation over
$s~{}\sim~{}\rho^{\pi_{w}},a~{}\sim~{}\pi_{w}$:

$\min\Big{(}r_{t}(w)A^{\pi_{w}}(s,a),\text{clip}\big{(}r_{t}(w),1-\epsilon,1+\epsilon\big{)}A^{\pi_{w}}(s,a)\Big{)}$

  * •

$r_{t}(w)=\frac{\pi_{w+\bigtriangleup w}(s,a)}{\pi_{w}(s,a)}$,

  * •

$\rho^{\pi_{w}}$ is the discounted state distribution, which is defined as
$\rho^{\pi_{w}}(s)=\sum_{t=0}^{\infty}\gamma^{t}Pr\\{s_{t}=s|s_{0},{\pi_{w}}\\}$,
and

  * •

$\epsilon\in\mathbb{R}$ is a hyper-parameter.

The specific version of PPO used in this paper is based on stable baselines 3
(cf. (sb3PPO, , OpenAI)).

###  3.2 Agent eco-system

Our approach, based on an eco-system of agents, makes use of the adaptability
(ability to perform properly to a certain degree on other environments without
additional training, also referred to as generalization) of each agent to
cover a higher number of environments in a given distribution of environments
$\mathcal{M}$ (Figure 1).

An agent is considered trained on a given environment when it can reach a
reward threshold $l$. The policy learned is then considered as the
satisfactory policy to solve the environment.
Our algorithm works as follow: new environments are drawn sequentially and
randomly from the distribution $M_{i}\sim\mathcal{M}$, $i~{}\in~{}\mathbb{N}$.
When the first environment is drawn, there is no agent in the pool. When a new
environment $M_{i}$ is presented to the algorithm, it first looks whether one
agent $e_{j}$ from the pool $\mathcal{E}=\\{e_{1},\ldots,e_{n}\\}$ can solve
$M_{i}$ without additional training. When a new environment is presented to
the eco-system, the policy $\pi$ of a given agent $e_{j}$, noted as
$\pi^{e_{j}}$, will be tested on the new environment. In the worst case, each
agent $e_{j}$ being part of the pool of agents $\mathcal{E}$ need to be
checked. In order to optimize the search of a good enough policy $\pi^{e*}$ in
the pool of agents $\mathcal{E}$, the agents in the pool are ordered by number
of environments they have already solved, thus checking the policy
$\pi^{e_{j}}$ that solved the highest number of environments first. If $e_{j}$
can solve $M_{i}$ then the new environment is added to the list of
environments the agent $e_{j}$ can solve. Otherwise, a new agent $e_{n+1}$ is
created, trained on $M_{i}$ until it can reach the reward threshold $l$ and is
added to the pool $\mathcal{E}$.

This new agent $e_{n+1}$ is then tested on the list of environments previously
solved by the pool $\mathcal{E}$, and each environment solved by $e_{n+1}$ is
added to its list of environments. If $e_{n+1}$ can solve all the environment
of an existing $e_{j}$ agent in the pool $\mathcal{E}$ then $e_{j}$ is removed
from the pool. This part is optional as adaptability will work without it, but
this ensure that the pool of agents $\mathcal{E}$ is kept to a reasonable
number of agents $e$ at the expense of having to test new agents on many
environments. This setting is needed in a system with limited resources
(memory, computational capabilities) in order to keep the pool of agents
manageable. The proposed solution is defined as an ecosystem of agents, which
means that agents are created as part of the pool, they are used for some time
and removed when a superior agent is found. Testing the newly created agents
on all past environment to see if they can replace an older agent is the way
chosen to implement the ”removal” part of the ecosystem. This approach can be
discussed, as it creates an additional inference workload each time a new
agent is added to the pool. However, by replacing less effective agents by a
smaller number of better agents, we also reduce the test time of the agents in
the pool each time an unknown environment is presented, which we have observed
happens more often than adding a new agent to the pool. Also this optimization
part is as indicated above not mandatory to ensure good performance of our
approach and it is only composed of inference and not training actions, which
are quite low on the computational side.

The pesudo-code of our approach is detailed in Algorithm 1 with $e_{n}$ being
the agent $n$ of the pool $\mathcal{E}$, $\mathcal{R}^{\pi_{e_{n}}}_{M_{i}}$
being the total reward gathered by the agent $e_{n}$ on the environment
$M_{i}$, $l$ representing the threshold to consider the policy $\pi^{e_{n}}$
successful, $e^{*}$ identifying the agent which solved successfully the new
environment, and $\delta^{e_{n}}$ being the list of environments solved
successfully by the agent $e_{n}$.

Algorithm 1 Eco-system - learn($M_{i}$

$e^{*}$ $\leftarrow$ $\emptyset$ #best agent found

$n$ $\leftarrow$ 0  #loop var.

while $e^{*}$ = $\emptyset$ and $\bigcup$ $e_{0...n}$ $\neq$ $\mathcal{E}$ do

#while good policy not found

#and not all agents reviewed

$\mathcal{R}^{\pi_{e_{n}}}_{M_{i}}=$ test_agent($e_{n}$,$M_{i}$)

#Total reward from $e_{n}$ on $M_{i}$

if $\mathcal{R}^{\pi_{e_{n}}}_{M_{i}}$ $\geq$ $l$ then #if $e_{n}$ solve
$M_{i}$

$e^{*}$ $\leftarrow$ $e_{n}$ #best agent found = $e_{n}$

else

$n$ $\leftarrow$ $n+1$

if $e^{*}$ = $\emptyset$ then #if $e^{*}$ not found

$e$ $\leftarrow$ new_agent()

while $\mathcal{R}^{\pi_{e}}_{M_{i}}$ $\leq$ $l$ do #while $e$ cannot solve
$M_{i}$

learn-epoch($e$,$M_{i}$)

$\mathcal{R}^{\pi_{e}}_{M_{i}}=$ test_agent($e$,$M_{i}$)

$\mathcal{E}$ $\leftarrow$ $\mathcal{E}$ \+ $e$ #add $e$ to the pool

#Following For statement is optional

#Only needed if optimization of the pool is needed

for $f$ $\in$ $\mathcal{E}$ do #for all agent $f$ in the pool

for $w$ $\in$ $\delta^{f}$ do #for all env. $w$ solved by $f$

$\mathcal{R}^{\pi_{e^{*}}}_{w}=$ test_agent($e^{*}$,$w$)

if $\mathcal{R}^{\pi_{e^{*}}}_{w}$ $\geq$ $l$ then #if $e$ can solve $w$

$\delta^{e}$ $\leftarrow$ $\delta^{e}$+$w$ #add $w$ to $e$ list

if $\delta^{f}$ $\in$ $\delta^{e}$ then

#if $e$ can solve all env. of $f$

$\mathcal{E}$ $\leftarrow$ $\mathcal{E}$ \- $f$ #remove $f$ from pool

Sort $\mathcal{E}$ by size $\delta$ descending order

else

$\delta^{e^{*}}$ $\leftarrow$ $\delta^{e^{*}}$+$M_{i}$ #add $M_{i}$ to $e^{*}$
list

Each agent of the pool is only trained on one environment and never retrained
on anything else during its life in the pool, ensuring that catastrophic
forgetting is not an issue. The distribution $\mathcal{D}$ of tasks
successfully solved by the eco-system, given the distribution of environments
$d(e)$ successfully solved by the agent $e$ in the pool $\mathcal{E}$ (seen or
not already seen), is defined as:

$\mathcal{D}=\bigcup_{e=0}^{n}d(e),e\in\mathcal{E}$

Essentially, the eco-system behaves as a knowledge accumulator, meaning that
it continues expanding its knowledge and optimizes its storage, while ensuring
it maintains agents to solve at least all the environments it has been
previously trained on.

##  4 Experimental setup

###  4.1 Environments

The experiments have been conducted on two different environments (a newly
developed Submarine environment and Minigrid), each having two different
settings.

##### Submarine

First, a submarine environment (Figure 2) is introduced for this paper as a
simple environment (with only 3 actions) while still being a challenging
environment (all obstacles are randomly placed at each level). The
observations can be easily customized to check how our approach performs on
different settings.

Figure 2: Submarine environment

In this environment, the agent controls a submarine which is moving from left
to right by one column after each step. The agent can choose between 3 actions
(up, stay, and down) and should reach the right side without hitting any
block. The blocks are randomly placed at each level using a seed (Figure 2).
Hitting a block or reaching the right side ends the game, with a -100 reward
in the first case and +100 in the second one. Two versions of the submarine
environment were created. The first one, submarine-easy, has a state
representation that allows for a better adaptability which includes the
current y position of the agent as well as a view of the next 5 columns
$[x:x+5]$ and all rows in these columns $[0:11]$. This state is a pattern
which can be easily found in multiple environments, which helps reducing over-
fitting and increase the probability that the output given by the neural
network in an environment for this state will also work on another
environment. The second version, submarine-hard, makes it more difficult to
learn a general strategy, which leads to smaller adaptability, as the state is
composed of the x and y position of the agent as well as the next 15 columns
$[x:x+15]$ and all rows in these columns $[0:11]$.

##### Minigrid

On top, we use Minigrid with two environments, namely the FourRoom and the
MultiRoom environments (cf. (gym_minigrid, , Chevalier-Boisvert et al.,
2018)). We have chosen this Minigrid setting to match the approach made by
Cobbe et al. KCobbe2019  using procedurally generated environments, which
means in our case that the map, start position, goal position, and obstacles
are positioned randomly according to the seed of each level, while the other
components of the experiments like reward given to the agent are kept the same
for each level. In order to benchmark our proposal we have chosen to compare
to the paper from Igl et al. Migl2019  as it performs very well on the
Minigrid environment. We run experiments on the FourRoom environment (Figure
3) and the MultiRoom environment (Figure 3). The MultiRoom environment was
used in the experiment of the benchmarking paper.

Figure 3: Minigrid FourRoom and Multiroom environment

In these environments the basic pixel view of Minigrid is used as state and
the agent gets a reward for reaching the goal in the maze. The reward is
defined as :

$1-0.9*(stepsUsed/maxStepsAllowed)$

###  4.2 Baselines

The baselines (single agent, IBAC-SNI introduced by Igl et al. Migl2019 ) are
chosen to match the approaches used in papers listed in the related work
section: a single agent trying to adapt to new environments. Another baseline
has been added to compare the performance of the eco-system to the
Reinforcement Learning voting ensemble technique described by Wiering et al.
MAWiering2008 .

##### Submarine setup

For the submarine environment the performance analysis was made between a
single agent on one side and the eco-system on the other side trying to solve
1000 never seen environments after having been trained on 1000 other
environments using the DeeR VFrancoisLavet2016  framework. The submarine
environment experiments have also been run using stable-baselines3 and PPO to
ensure that the results were not linked to a particular algorithm. The choice
has been made to use the exact same code, used for training each agent in the
eco-system, to train the single agent approach in order the performance
increase of the eco-system can be easily identified. We have decided not to
compare to IBAC-SNI on this environment has we wanted to compare to the
benchmarking paper only on environments used in it (Minigrid), to ensure a
fair comparison.

##### Minigrid setup

For the Minigrid FourRoom and MultiRoom environments, the performance analysis
was made between the code associated with the paper from Igl et al. Migl2019
and the eco-system after being trained for $10^{8}$ steps on different
environments randomly chosen, using the stable-baselines3 framework and the
PPO learning algorithm. On top, a comparison has been made for the Minigrid
FourRoom environment, as it was the most challenging one from both Minigrid
environments, between the eco-system and our implementation of the vorting
ensemble technique described by Wiering et al. MAWiering2008  and indicated as
the most efficient ensemble technique they have tested. For this experiment,
we used 200 environments and trained one agent on each of them and added the
agents to a pool. During the training, every 20 environments, we have tested
the pool of agents (composed of 20 agents, then 40, then 60 …) on 1000 never
seen environments and computed the average of total reward gathered. For
testing the pool of agents on each environment, the following approach was
taken: at each step on each environment, each agent in the pool was voting for
which action to take. The majority won and the corresponding action was taken.

We have also tested, in the Minigrid FourRoom Environment, if changing the
order for testing the agents in the pool has an effect on the results, by
ordering the agent in the descending order of solved environments and then on
the ascending order.

###  4.3 Training algorithms

In the experiments conducted, different learning algorithms (Double DQN, PPO)
and frameworks (DeeR, stable-baselines3) have been used. The threshold is
noted $l$, in the following formulas, but it will differ based on the type of
environment used, for example $l_{s}$ is the threshold for the submarine
environment. The threshold used is the same for all environments of a given
type. The threshold $l_{s}$ has been set to 100 for the submarine game as any
total reward $\geq$ 100 means that the agent has reached the goal. The
threshold $l_{m}$ has been set to 0.8 for the Minigrid environments (1 for
reaching the goal minus a penalty taken for each step). In order to define the
Minigrid threshold at 0.8, we have assessed if changing this threshold has a
direct impact on the generalization result of the experiment. We have tested 3
different threshold (0.8, 0.85 and 0.9). The results are indicated in section
5. More details about the technical implementation, and hyper-parameters are
indicated in the appendix of this paper.

###  4.4 Performance metrics

In order to asses the performance of each approach, 3 main metrics are
considered : the adaptability index, the catastrophic avoidance index and the
amount of access to the environments.

##### The adaptability index

This metric is based on testing on a number of unseen environments
$M_{i}\in\mathcal{M}$. It expressed whether an approach can perform over a
certain threshold $l$ or what average of total reward $\mathcal{R}$ is
gathered by the approach. The adaptability index based on a threshold, noted
as $\omega$, is indicated as a percentage representing the percentage of times
the approach has performed over the threshold $l$ on all the new environments
$M_{i}\in\mathcal{M}$ on which it was tested. It can be noted as:

$\omega=\frac{\sum_{i=0}^{n}\begin{Bmatrix}R_{M_{i}}>=l\rightarrow 1\\\
R_{M_{i}}<l\rightarrow 0\end{Bmatrix}*100}{n}$

The adaptability index based on the average reward gathered, noted as $\zeta$,
is indicated as a float value, showing the average of total rewards
$\mathcal{R}$ gathered over all the new environments $M_{i}\in\mathcal{M}$ on
which the approach was tested. It is formalized as follows:

$\zeta=\frac{\sum_{i=0}^{n}R_{M_{i}}}{n}$

##### The catastrophic forgetting avoidance index

This metric, noted $\xi$, is based on testing periodically during the initial
training how each approach performs on the environment $M_{i}\in\mathcal{M}$
on which it has already been trained. It checks that the approach does not
forget the knowledge already accumulated. The metric is indicated as a
percentage representing the percentage of time the agent $e$ total reward
$\mathcal{R}$ was over the threshold $l$ on previously learned environments.
It can be noted as:

$\xi=\frac{\sum_{i=0}^{n}\begin{Bmatrix}R_{M_{i}}>=l\rightarrow 1\\\
R_{M_{i}}<l\rightarrow 0\end{Bmatrix}*100}{n}$

##### The amount of access to the environments

: this metric is based on testing periodically how much access (number of
states observed) has been made to the environment by each approach. This is a
fair comparison as it accounts for the number of access made during learning
on a given environment but also the number of access used during inference
when testing if a given agent can solve an environment. The metrics are
calculated periodically after each approach has been presented to 50
additional environments (and completed the associated training if necessary).

##  5 Results

###  5.1 Submarine environment

##### Submarine-easy

The first set of results has been run on the submarine-easy environment easier
to generalize. For the adaptability index based on threshold $\omega$,
considering a distribution of 1000 new environment for testing (Figure 4), the
eco-system is far more stable than the one agent approach and also gives a
better result at the end of the training.

Figure 4: Adaptability and catastrophic forgetting DDQN

Not forgetting any good policy learned before is a key concept of the eco-
system approach, which is confirmed by the catastrophic forgetting avoidance
index $\xi$ (Figure 4). The single agent approach forgots environments
previously learned while learning new environments and the eco-system is
always able to solve any environment previously solved. The PPO algorithm with
stable-baselines3 framework gives similar results, even if the performance of
the single agent is lower than when using DDQN. (Figure 5).

Figure 5: Adaptability and catastrophic forgetting PPO

##### Submarine-hard

The second experiment was conducted on the harder to generalize submarine
environment (submarine-hard). The eco-system achieves similar performance on
the adaptability index based on threshold $\omega$ (Figure 6) on new
environments while the performance of the single agent decreases.

Figure 6: Adaptability and catastrophic forgetting DDQN

The catastrophic forgetting avoidance index $\xi$ (Figure 6) explains the
results seen with the other chart, the catastrophic forgetting effect creates
a lot of instability for the single agent. The agent is forgetting a
significant part of what it has learned before. This is typical, as the agent
aims to solve the current environment, it adapts his neural network to solve
specific states, which means over-fitting to this environment. The likelihood
to find similar states on another environment leading to the same action is
really low, which means that the single agent has to modify its neural network
to adapt to the new states and actions, creating what is called catastrophic
forgetting.

Figure 7: Adaptability and catastrophic forgetting PPO

We also run the same experiments using the PPO algorithm with stable-
baselines3 framework, and got results confirming the superiority of the Eco-
system approach (Figure 7).

On both experiments we can see that the better performance in adaptability and
the removal of catastrophic forgetting comes at the expense of the number of
accesses made to the environments by the eco-system, see Figure 8.

Figure 8: Number of accesses to the environment (easy and hard)

We can however see that when the environment is more challenging for an agent
to generalize on, the difference between the single agent and the eco-system
in term of accesses is becoming less and less (see Figure 8). If we look at
the duration to train the single agent and the eco-system on the initial
environments for the hard to generalize version (Figure 9), it is quite
similar. This indicates that even if the eco-system is making more access to
the environment, most of them are made as inference and not as learning, which
has a lower impact on the overall duration, where the single-agent focus most
of the access on learning and Neural Network back-propagation.

Figure 9: Training duration comparison

###  5.2 Minigrid

##### Minigrid FourRoom

The next experiment was run on the Minigrid FourRoom environment to compare
how the eco-system performs in comparison with the proposed approach by Igl et
al. Migl2019 . The generalization assessment for the FourRoom and MultiRoom
environments is based on the average of total return $\zeta$ to match the
metrics used in the benchmarking paper, which is independent of the threshold.
The threshold is only used to end the training of the new agent on a given
environment before it gets added to the pool. The FourRoom environment is hard
to generalize on, as a lot of parameters are randomly chosen: the start
position of the agent, the position of the goal to reach, and the opening in
the walls). The eco-system performs better than IBAC-SNI algorithm used in the
benchmarking paper ((Migl2019, , Igl et al., 2019)) on new randomly chosen
environments (see Figure 10) after being trained on $10^{8}$ steps on other
randomly chosen environments. This confirms the results gathered on the
submarine environment, when the generalization / adaptability is difficult for
a single agent, the eco-systems performs significantly better when compared to
a strong baseline, also reaching quicker an higher average total reward. The
standard deviation of the eco-system results is also smaller than the one of
the benchmarking approach, showing it is far more stable. The eco-system also
outperforms the ensemble voting technique described by Wiering et al.
MAWiering2008  (cf. Figure 10) We can also notice that the ensemble voting
technique seems to see a decrease in its performance when more agents are
added. Having the choice of the action to be taken linked with the answer of
the majority of agents can create cases where the agents in the pool vote for
a bad action if only a small number of agents were able to properly solve the
environment.

Figure 10: Benchmark on Minigrid FourRoom and Multiroom

We also study the impact of the ordering function applied to the pool of
agents. In the standard version, the agent that has solved the highest number
of environments is tested first when looking at a new environment, based on
the fact that an agent scoring high on generalization may have a better chance
to solve the new environment. Another approach is to consider that the more
environments an agent can solve, the more we may reach the limit of the neural
network to solve new environments and that ordering by the reverse order may
bring more performance. The results show that this has no impact on the
performance of the eco-system approach (see Figure 10). In order to find the
proper threshold for the Minigrid environment, the following experiment has
been run using 3 options : 0.8, 0.85 and 0.9. (Figure 11) The results of the
experiment indicate that changing the threshold does not have a significant
impact on the result of the experiments.

Figure 11: Impact on the threshold modification for Minigrid FourRoom
environment

The modification of the threshold has an impact on the number of agents needed
in the pool to reach the same result. We have then made the choice to use 0.8
as it was the one using the minimum number of agents and offering a better
generalizability by agent. The following results show the number of agents in
the pool after training the eco-system on 150 environments.

Threshold | Number of agents
---|---
0.8 | 79
0.85 | 90
0.9 | 107

This shows that our approach is very robust and not impacted by modifications
/ adjustments of some of its core parameters.

##### Minigrid MultiRoom

The last experiment was run on the Minigrid MultiRoom environment used by Igl
et al. Migl2019 . This experiment allows to benchmark our approach against a
state-of-the-art published approach by using the exact same environment and
performance indicator. Our approach and the IBAC-SNI approach are quite
different, the eco-system focus on solving a given environment before moving
to another one while IBAC-SNI trains on multiple environments without trying
to fully solve one, which explains why after a similar number of training
steps, the eco-system offers better generalization performance. As seen on
Figure 10, the eco-system consistently outperforms the IBAC-SNI approach, even
if IBAC-SNI gets closer to the eco-system results by the end of the training.

##  6 Discussion / conclusion

We have introduced a novel approach using an eco-system of agents to improve
generalizability of Reinforcement Learning. The eco-system approach stores a
limited number of agents in its pool and harvest their adaptability /
generalization capabilities to produce a system which can generalize to a
higher number of environments, without suffering from catastrophic forgetting.
Our approach of an eco-system of agents adapts better (i.e. a higher total
reward $\mathcal{R}$) to never seen environments and keeps performing better
on already learned environments, in all configurations explored in this paper,
than a single agent or ensemble techniques.

The experiments also show that our approach suppress fully catastrophic
forgetting when adapting to new environments. This means that when a given
environment is learned, it will never be forgotten. The agent solving it from
the pool may change, but there will always be at least one agent from the pool
that can solve it. This is a key point of our solution, which makes it a
perfect fit for setup where forgetting a previously trained environment is not
an option. We have also showed that our proposed solution performs at least as
well as the other solutions explored in this paper, and that it is less prone
to adaptability level drop or training time increase than the other solutions
as the generalizability of the environment becomes harder. Our approach, by
reducing the amount of retrain needs can also be used successfully on setups
where the computational cost of multiple retrainings is too high (IOT), by
pre-training the ecosystem and then deploying it on a low computational setup
where it will leverage the adaptability accumulated by the eco-system.

Future work will be focused on trying to optimize the initialization of the
agents in the pool and ensuring a better generalization capabilities to the
newly created agents. We will also focus on trying to find a way to reduce the
amount of environment access being done by the agents, especially when trying
the agents in the pool to find one which can solve a never seen environment.

## References

  * (1) [Cobbe et al., 2019] Cobbe, K., Klimov, O., Hesse, C., Kim, T. & Schulman, J.. (2019). Quantifying Generalization in Reinforcement Learning. Proceedings of the 36th International Conference on Machine Learning, in Proceedings of Machine Learning Research 97:1282-1289 Avai lable from https://proceedings.mlr.press/v97/cobbe19a.html
  * (2) [Farebrother et al., 2018] Farebrother, J., Machado, M. C., & Bowling, M. (2018). Generalization and regularization in DQN. arXiv preprint arXiv:1810.00123.
  * (3) [Justesen et al., 2018] Justesen, N., Torrado, R. R., Bontrager, P., Khalifa, A., Togelius, J., & Risi, S. (2018). Illuminating generalization in deep reinforcement learning through procedural level generation. arXiv preprint arXiv:1806.10729.
  * (4) [Ghiassian et al., 2019] Ghiassian, S., Rafiee, B., Lo, Y. L., & White, A. (2020). Improving performance in reinforcement learning by breaking generalization in neural networks. arXiv preprint arXiv:2003.07417.
  * (5) [Van Hasselt et al., 2019] Van Hasselt, H., Guez, A., & Silver, D. (2016, March). Deep reinforcement learning with double q-learning. In Proceedings of the AAAI conference on artificial intelligence (Vol. 30, No. 1).
  * (6) [Nichol et al., 2018] Nichol, A., Pfau, V., Hesse, C., Klimov, O., & Schulman, J. (2018). Gotta learn fast: A new benchmark for generalization in rl. arXiv preprint arXiv:1804.03720.
  * (7) [Packer et al., 2018] Packer, C., Gao, K., Kos, J., Krähenbühl, P., Koltun, V., & Song, D. (2018). Assessing generalization in deep reinforcement learning. arXiv preprint arXiv:1810.12282.
  * (8) [François-Lavet et al., 2016] V. François-Lavet et al. (2016).DeeR. Available from https://deer.readthedocs.io/.
  * (9) [Raffin et al., 2019] A.Raffin, A.Hill, M.Ernestus, A.Gleave, A.Kanervisto, Anssi, & N.Dormann. (2019). Stable Baselines3. Available from https://github.com/DLR-RM/stable-baselines3
  * (10) [François-Lavet et al., 2018] François-Lavet, V., Henderson, P., Islam, R., Bellemare, M. G., & Pineau, J. (2018). An introduction to deep reinforcement learning. arXiv preprint arXiv:1811.12560.
  * (11) [Chevalier-Boisvert et al., 2018] M. Chevalier-Boisvert, L. Willems, & S. Pal. (2018). Minimalistic Gridworld Environment for OpenAI Gym. Available from https://github.com/maximecb/gym-minigrid
  * (12) [Igl et al., 2019] Igl, M., Ciosek, K., Li, Y., Tschiatschek, S., Zhang, C., Devlin, S., & Hofmann, K. (2019). Generalization in reinforcement learning with selective noise injection and information bottleneck. arXiv preprint arXiv:1910.12911.
  * (13) [Sonar et al., 2021] Sonar, A., Pacelli, V., & Majumdar, A. (2021, May). Invariant policy optimization: Towards stronger generalization in reinforcement learning. In Learning for Dynamics and Control (pp. 21-33). PMLR.
  * (14) [Alver et al., 2020] Alver, S., & Precup, D. (2020). A brief look at generalization in visual meta-reinforcement learning. arXiv preprint arXiv:2006.07262.
  * (15) [Duan et al.,2016] Duan, Y., Schulman, J., Chen, X., Bartlett, P. L., Sutskever, I., & Abbeel, P. (2016). Rl2: Fast reinforcement learning via slow reinforcement learning. arXiv preprint arXiv:1611.02779.
  * (16) [Chen, 2020] Chen, J. Z. (2020). Reinforcement Learning Generalization with Surprise Minimization. arXiv preprint arXiv:2004.12399.
  * (17) [Lu et al., 2020] Lu, X., Lee, K., Abbeel, P., & Tiomkin, S. (2020). Dynamics Generalization via Information Bottleneck in Deep Reinforcement Learning. arXiv preprint arXiv:2008.00614
  * (18) [Sutton and Barto, 1998] Sutton, R. S., & Barto, A. G. (1998). Introduction to reinforcement learning (Vol. 135) Cambridge: MIT press.
  * (19) [Schulman et al., 2017] Schulman, J., Wolski, F., Dhariwal, P., Radford, A., & Klimov, O. (2017). Proximal policy optimization algorithms. arXiv preprint arXiv:1707.06347.
  * (20) [Wiering et Van Otterlo, 2012] Wiering, M. A., & Van Otterlo, M. (2012). Reinforcement learning. Adaptation, learning, and optimization.
  * (21) [Atkinson et al., 2021] Atkinson, C., McCane, B., Szymanski, L., & Robins, A. (2021). Pseudo-rehearsal: Achieving deep reinforcement learning without catastrophic forgetting. Neurocomputing 428, 291-307.
  * (22) [Kirkpatrick et al., 2016] Kirkpatrick, J., Pascanu, R., Rabinowitz, N., Veness, J., Desjardins, G., Rusu, A. A., … & Hadsell, R. (2017). Overcoming catastrophic forgetting in neural networks. Proceedings of the national academy of sciences 114(13), 3521-3526.
  * (23) [Rusu et al., 2016] Rusu, A. A., Rabinowitz, N. C., Desjardins, G., Soyer, H., Kirkpatrick, J., Kavukcuoglu, K., … & Hadsell, R. (2016). Progressive neural networks. arXiv preprint arXiv:1606.04671.
  * (24) [OpenAI] Stable baselines 3. Details on stable baselines 3 PPO implementation. Available from https://spinningup.openai.com/en/latest/algorithms/ppo.html.
  * (25) [Wiering et al., 2008] Marco A Wiering and H. V. Hasselt Ensemble Algorithms in Reinforcement Learning IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics)