The detection of moving objects is a trivial task when performed by vertebrate retinas, yet a complex computer vision task. This PhD research programme has made three key contributions, namely: 1) a mhsnn architecture for detecting horizontal and vertical movements, 2) a hsmd algorithm for detecting object motion and 3) the neurohsmd , a real-time neuromorphic implementation of the hsmd  algorithm.

The mhsnn is a customised 4 layers snn architecture designed to reflect the basic connectivity, similar to canonical behaviours found in the majority of vertebrate retinas (including human retinas). The architecture, was trained using images from a custom dataset generated in laboratory settings. Simulation results revealed that each cell model is sensitive to vertical and horizontal movements, with a detection error of 6.75% contrasted against the teaching signals (expected output signals) used to train the mhsnn. The experimental evaluation of the methodology shows that the mhsnn was not scalable because of the overall number of neurons and synapses which lead to the development of the hsmd.

The hsmd  algorithm enhanced an existing dbs algorithm using a customised 3-layer snn. The customised 3-layer snn was used to stabilise the foreground information of moving objects in the scene, which improves the object motion detection. The algorithm was compared against existing background subtraction approaches, available on the opencv library, specifically on the cdnet2012 and the cdnet2014 benchmark datasets. The accuracy results show that the hsmd was ranked overall first and performed better than all the other benchmarked algorithms on four of the categories, across all eight test metrics. Furthermore, the hsmd is the first to use an snn to enhance the existing dynamic background subtraction algorithm without a substantial degradation of the frame rate, being capable of processing images $720\times480$ at 13.82 fps (cdnet2014) and $720\times480$ at 13.92 fps (cdnet2012) on a High Performance computer (96 cores and 756 GB of RAM). Although the hsmd analysis shows good pcc on the cdnet2012 and cdnet2014, it was identified that the 3-layer customised snn was the bottleneck, in terms of speed, and could be improved using dedicated hardware.

The neurohsmd is thus an adaptation of the hsmd algorithm whereby the SNN component has been fully implemented on dedicated hardware [Terasic DE10-pro fpga board]. opencl was used to simplify the fpga design flow and allow the code portability to other devices such as fpga and gpu. The neurohsmd was also tested against the cdnet2012 and cdnet2014 datasets with an acceleration of 82% over the hsmd algorithm, being capable of processing $720 \times 480$ images at 28.06 fps (cdnet2012) and 28.71 fps (cdnet2014).

The copyright in this work is held by the author. You may copy up to 5% of this work for private study, or personal, non-commercial research. Any re-use of the information contained within this document should be fully referenced, quoting the author, title, university, degree level and pagination. Queries or requests for any other use, or if a more substantial copy is required, should be directed to the author.

[title=List of Acronyms,type=]

CHAPTER: APPENDIX

The mhsnn architecture presented in Chapter <ref> was published on the Conference proceedings of the International Joint Conferences on Neural Networks 2018 <cit.>.

P. Machado, A. Oikonomou, G. Cosma and T. M. McGinnity, "Bio-Inspired Ganglion Cell Models for Detecting Horizontal and Vertical Movements" 2018 ijcnn, 2018, pp. 1-8, doi: 10.1109/IJCNN.2018.8489439.

The hsmd algorithm presented in Chapter <ref> was published on the ieee Access <cit.>.
P. Machado, A. Oikonomou, J. F. Ferreira and T. M. McGinnity, "HSMD: An object motion detection algorithm using a Hybrid Spiking Neural Network Architecture," in ieee Access, doi: 10.1109/ACCESS.2021.3111005.

The neurohsmd algorithm presented in Chapter <ref> was submitted to the ACM Transactions on Reconfigurable Technology and Systems.


# Shift manipulation of intrinsic localized mode
in ac driven Klein Gordon lattice

Hirotaka Araki<EMAIL_ADDRESS>Takashi Hikihara
<EMAIL_ADDRESS>Department of Electrical Engineering, Kyoto
University, Katsura, Nishikyo, Kyoto 615-8510, Japan.

###### Abstract

Shift manipulation of intrinsic localized mode (ILM) is numerically discussed
in an ac driven Klein Gordon lattice. Before the manipulation, we introduce
the 2-degree of freedom nonlinear system, which is obtained by reducing the
lattice. In the reduced system, two localized modes correspond to ILMs in the
original system. These two localized modes are switched by changing a coupling
constant adiabatically. Based on the result obtained in the reduced system, we
numerically show that the shift manipulation of ILM is achieved, like a
particle, by using an adiabatic change of coupling constant between the
targeted sites in the lattice.

###### keywords:

Klein Gordon lattice , intrinsic localized mode , discrete breather ,
localized mode

##  1 Introduction

Spatially localized and periodic vibration modes existing in nonlinear
discrete systems are called Intrinsic localized modes(ILMs) or discrete
breathers(DBs). ILMs were theoretically discovered by Sievers and Takeno in
1988 [1]. After the discovery, numerous analytical, numerical, and
experimental studies have been reported [2]. Many researchers have
analytically and numerically investigated ILMs from various perspectives such
as the existence [3, 4], the stability [5, 6] and the movability [7, 8]. ILMs,
on the other hand, have been experimentally observed in various physical
systems, for example, Josephson-junction array [9], waveguide arrays [10],
coupled cantilever array [11], and nonlinear transmission line [12, 13]. Also,
molecular dynamics have revealed that ILMs can appear in nanoscale systems
[14, 15, 16, 17, 18]. These experimental and numerical studies suggest the
phenomenological universality of ILMs and the possibility of applications in
both fundamental science and practical engineering. Indeed, recently, studies
for ILM applications have been reported in energy management in device [19],
energy harvesting [20], and sensors, actuators [21, 22, 23], and so on. ILMs
also have the potential to be used for phonon engineering in nonlinear lattice
because ILMs can scatter phonons [24, 25] and affect thermal conductivity
[26].

Control of ILMs is smoothly considered as a key to achieve ILM applications as
mentioned above. It includes generation, annihilation, and spatial position
manipulation of ILMs. These manipulation has been investigated numerically and
experimentally [19, 27]. Based on the results, we are going to focus on the
spatial position manipulation. Sato et al. experimentally and numerically
demonstrated a shift of the position of ILMs using a local impurity in
nonlinear transsmision line and coupled cuntilever array [19, 28]. The
manipulation is based on the observations that ILMs is attracted or repulsed
by the impurity. Kimura at el. proposed another method for the spatial
position manipulation of ILMs in coupled cantilever array, which is called the
capture and release manipulation of the traveling ILMs [23]. In the method, a
nonlinear coupling coefficient is varied at the appropriate time to switching
between the pinned stable ILM and the traveling ILM. It is based on the
stability change due to the saddle-node bifurcation of ILMs.

Spatial position manipulation of ILMs is also investigated for engineering
applications in nonlinear optics. In this field, ILMs appear as discrete
optical solitons in periodic structures describing arrays of coupled optical
waveguides. Light propagation in the waveguide arrays is described by
nonlinear Schrödinger equation in tight binding approximation. In the
equation, the propagation distance along the waveguide plays the role of time
in the dynamic lattice equations. Mobility of discrete optical solitons is
related to an application for multiport optical switching [29]. So, for the
unambiguously controlled switching using an appropriate perturbation, there
are various proposed methods. Rodrigo et al. numerically demonstrate digitized
switching [30, 31]. They used a steplike variation of the waveguide coupling
for adjusting Peierls-Nabarro(PN) potential appeared from the lattice
discreteness [32]. To enhance the mobility of discrete optical solitons, other
studies investigate bifurcations of solitons and behavior of PN potential when
varying the power of solitons [33, 34, 35].

The main purpose of this letter is to propose a new method for shift
manipulation of ILMs in an ac driven Klein Gordon lattice. Here, shift
manipulation of ILM corresponds to the moving of an ILM from the site on which
the ILM exists to its neighboring site in the specified direction without
decaying. It is the achievement of spatial position manipulation of ILMs. In
the following, at first, we reduce the lattice to the 2-degree of freedom
nonlinear system so as to avoid the multi-degree of freedom of the lattice.
Two coexisting localized modes in the reduced system is numerically
investigated and it is revealed that switching of these localized modes are
achieved by adiabatic change of the coupling constant. Based on this result,
secondly, the manipulation method is prepared and achieved by adiabatic change
of the coupling constant between targeted sites.

##  2 The model and ILMs

In this letter, we consider shift manipulation of ILM in ac driven Klein
Gordon lattice The Hamiltonian of Klein Gordon lattice is written as

$H=\sum_{n=1}^{N}\left[\frac{\dot{u}_{n}^{2}}{2m}+V(u_{n})+\frac{C}{2}(u_{n+1}-u_{n})^{2}\right],$
(1)

where $u_{n}$ is displacement of the $n$-th site particle from its equilibrium
point, $\dot{u}_{n}$ is its velocity and $C$ is a coupling constant. We take

$V(z)=\frac{z^{2}}{2}+\frac{z^{4}}{4},$ (2)

for the on-site potential with hard-type anharmonicity. The investigation of
this manipulation is summarized for the case of soft-type anharmonicity in the
appendix.

The equation of motion of Klein Gordon lattice with dissipation and periodic
external driver is written as

$\displaystyle m\ddot{u}_{n}+\gamma\dot{u}_{n}+u_{n}+u_{n}^{3}-C(u_{n-1}$
$\displaystyle+u_{n}-u_{n+1})$ $\displaystyle=(-1)^{n}F\cos\omega{t},$ (3)

where $\gamma$ is dissipation coefficient, $F$ is the driving amplitude, and
$\omega$ is driver frequency. Each particles are periodically driven in
reverse phase.

This model is derived from the coupled cantilever arrays devised by Sato and
Sievers. For simplicity, however, we neglect the influence of longer than the
nearest-neighbor interaction range and the higher-order nonlinearity in the
interaction potential. We also use the reverse-phase driver instead of the in-
phase driver. Each parameters are set as follows based on previous researches:
$m=1,\ C=0.07953,$ and $\gamma=1.534\times 10^{-4}$. We also set the driving
amplitude and frequency: $F=0.008$ and $\omega=1.5$ respectively.

There are two kinds of ILM with different spatial symmetry in this lattice.
One is Sievers-Takeno mode which is spatially symmetric and centered on sites
[1]. The other is Page mode which is spatially antisymmetric and centered
between sites [36]. Both of them show in phase to the driver. Fig. 1(a) shows
the amplitude distribution of the Sievers-Takeno mode. Fig. 1(b) shows the
amplitude distribution of Page mode. By using the linear stability analysis
method based on Floquet theory, it is confirmed that Sievers-Takeno mode is
stable and Page mode is unstable. In this letter, we focus Sievers-Takeno mode
as the target of shift manipulation.

Figure 1: Amplitude distribution of Sievers-Takeno mode (a) and Page mode (b).
Each points indicate the position of particles when they have largest
displacements. Each particles oscillate with the same period.

##  3 The reduced system and NNMs

At this stage, it seems to be difficult to directly devise a method of the
manipulation in this lattice due to its nonlinearity and multi-degree-of-
freedom nature. So, in the following, we introduce a reduced system of the
lattice to consider the shift manipulation of ILM. Here, we focus on two
sites. One is a site where the ILM exists and the other is its neighboring
site at the direction to move the ILM. The manipulation is a transition of ILM
between two adjacent sites. The setting makes the influence from outside those
two sites neglect. This is allowed by the localization of ILM and the
sufficiently small coupling constant. As a result, we reduce the ac driven
Klein Gordon lattice to 2-degrees of freedom system. The equation of motion of
this 2-degree of freedom nonlinear system is written as

$\displaystyle\ddot{u}_{1}+\gamma\dot{u}_{1}+u_{1}+u_{1}^{3}-C(u_{2}-u_{1})=-F\cos\omega{t},$
(4)
$\displaystyle\ddot{u}_{2}+\gamma\dot{u}_{2}+u_{2}+u_{2}^{3}-C(u_{1}-u_{2})=F\cos\omega{t},$
(5)

where subscripts are replaced by 1 and 2 for simplicity. The reduction allows
us to avoid the disturbanceto the mode by the structure of the original
lattice.

We investigate the system described by Eq. (4) and (5) using the concept of
nonlinear normalized modes(NNMs) [37, 38, 39]. NNMs is proposed by extending
the concept of linear normal modes to examine nonlinear vibration systems in
which superposition principle and no-interference between modes no longer
hold. The most simplest definition of NNMs is a vibration in unison of the
system.

We obtain three coexisting NNMs by the shooting method using Newton-Raphson
method [39]. Fig. 2 shows the behavior of each NNM. Fig. 2(a) shows the anti-
phase mode which is unstable. Fig. 2(b) shows the localized mode which are the
NNM localized at oscillator 1. Fig. 2(c) shows the localized mode which are
the NNM localized at oscillator 2. These two localized modes are specific to
NNMs and stable. The anti-phase mode corresponds to unstable Page-mode in the
original system. Two localized modes correspond to stable Sievers-Takeno
modes. So, shift manipulation of ILM in the ac driven Klein Gordon lattice
corresponds to switching these two coexisting localized modes in the reduced
system.

Figure 2: Time development of three NNMs in the reduced system. Solid lines is
displacement of oscillator 1. Dotted lines is displacement of oscillator 2.
(a) An anti-phase mode. (b) and (c) Two localized modes.

##  4 Switching of two coexisting localized modes

In this section, we consider switching of two coexisting localized modes in
the reduced system.

Here, we need to investigate the coupling constant dependency of coexisting
NNMs in the system, because the coupling constant is the only parameter that
we can tune to achieve the switching. We vary the coupling constant and
numerically obtain coexisting NNMs at each coupling constant value. The NNMs
are characterized by the amplitude ratio $k$ defined by $k=r_{1}/r_{2}$, where
$r_{1}$ denotes the amplitude of oscillator 1, and $r_{2}$ the amplitude of
oscillator 2.

Fig. 3 shows amplitude ratio $k$ in each of the coexisting NNMs with respect
to the coupling constant $C$. The solid line indicates stable NNMs and the
dashed line unstable NNMs. The straight line at k=1 corresponds to the anti-
phase mode. The curve above the line corresponds to NNMs localized at
oscillator 2. The curve below the line corresponds to NNMs localized at
oscillator 1. When the coupling constant is large, only stable anti-phase mode
exists. However, the anti-phase mode loses its stability at $C\simeq 0.425$
and two coexisting localized modes appear. This is a pitchfork bifurcation of
NNMs.

We consider the switching of two coexisting localized modes on the basis of
the bifurcation diagram of NNMs. In the diagram, two coexisting localized
modes are smoothly connected via the pitchfork bifurcation point. So, we
implement the switching two coexisting localized modes by tuning the coupling
constant. This manipulation needs to be achieved adiabatically because a rapid
change of the coupling constant destroys NNMs. Here, adiabatic change implies
sufficiently slow change compared to the period of driver $T=2\pi/\omega$.
Fig. 4 explains how we change the coupling constant. First, we increase the
coupling constant linearly from the original value $C=0.07953$ to $C=0.45$
which is greater than the bifurcating point at $C=0.425$. This process take
$T_{\mathrm{I}}$. Next, we keep the coupling constant at the value $C=0.45$
during period $T_{\mathrm{C}}$. Lastly, we reduce the coupling constant to the
original value during period $T_{\mathrm{R}}$.

Here, we switch the mode from NNM localized at oscillator 1 to one at
oscillator 2. Fig. 5(a) and (b) show simulation results when $T_{\mathrm{I}},\
T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set at $250T,\ 20T$, and $250T$,
respectively. Points indicate Poincare maps of orbits of each oscillator
during the switching. The solid line is the trajectory for a period of each
oscillator before the switching and the dotted line is during period
$T_{\mathrm{C}}$, and the dashed line is after the switching. In the
following, we explain the development of the state during the switching by
points. First, during the period $T_{\mathrm{I}}$, the initial state(solid
line) gradually approaches the anti-phase mode(dotted line). The anti-phase
mode appears during the period $T_{\mathrm{C}}$. Finally, during the period
$T_{\mathrm{R}}$, the anti-phase mode approaches NNM localized at oscillator
2(dashed line). Thus, under this condition, the switching is achieved between
the coexisting localized modes.

Fig. 5(c) and (d) show simulation results when $T_{\mathrm{I}},\
T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set at $250T,\ 30T$, and $250T$,
respectively. The initial state once becomes the anti-phase mode, however, it
finally returns to the original state. Therefore, this condition will not
achieve the switching of coexisting localized modes is not achieved. From
these results, the result of this switching depends on $T_{\mathrm{C}}$.
Further numerical investigations reveal that result of the manipulation
depends not only $T_{\mathrm{C}}$, but also other factors present in the
process such as $T_{\mathrm{I}}$ and $T_{\mathrm{R}}$. To be more precise, the
result of the switching depends on how to change the coupling constant.

Figure 3: Relationship between coupling constant and amplitude ratio $k$ in
each of the coexisting NNMs. The solid line corresponds to stable NNMs. The
dashed line corresponds to unstable NNMs.  Figure 4: Time development of
stiffness during the manipulation. The unit time on the horizontal axis is $T$
which is the period of driver. Figure 5: (a) and (b) Simulation results of the
switching when $T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set
at $250T,\ 20T$, and $250T$ respectively. (c) and (d) Simulation results of
the switching when $T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are
$250T,\ 30T$, and $250T$ respectively. Points indicates Poincare map of orbits
of each oscillators during the manipulation. Solid line is trajectory for a
period of each oscillators before the switching and dotted line is during
period $T_{\mathrm{C}}$, and dashed line is after the switching.

##  5 Shift manipulation of ILM

At last, we consider the shift manipulation of ILM in the ac driven Klein
Gordon lattice. Based on the results obtained in the reduced system, we expect
that the manipulation is achieved by the adiabatic change of a coupling
constant between two sites at which ILM originally exists and the neighbor
site of the direction to move the ILM. Here, we handle an ILM at the 5th site
to the 6th site. So, we change a coupling constant between the 5th site and
the 6th site adiabatically as in the case of switching coexisting localized
modes.

Fig. 6(a) shows the behavior of the ILM during the manipulation when
$T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set at $250T,\
20T$, and $250T$, respectively. The heat map indicates the time development of
energy at each site. The dark region indicates a high energy state. Before the
manipulation, most of the energy is concentracted at the 5th site, which
corresponds to the ILM at the 5th site. We can prepare this ILM by the
shooting method using Newton-Raphson method [40]. As we linearly increase the
coupling constant to get closer to 0.45 during the period $T_{\mathrm{I}}$,
the energy of the 5th site gradually excludes into the 6th site. When keeping
the coupling constant 0.45 during the period $T_{\mathrm{C}}$, most of the
energy is almost equally distributed at both of the 5th and the 6th sites.
During this period, the amplitude distribution looks like the Page mode.
However, this is not the Page mode because the lattice is not homogeneous
during this period. As we linearly reduce the coupling constant to return to
the original value during the period $T_{\mathrm{R}}$, the energy starts to
gather in the 6th site. After the manipulation, most of the energy is
concentracted at the 6th site. As a result, the manipulation succeeds.

Fig. 6(b) shows the behavior of the ILM during the manipulation when
$T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set at $250T,\
30T$, and $250T$ respectively. Before the period $T_{\mathrm{R}}$, this figure
shows the same behavior when $T_{\mathrm{C}}$ is set at $20T$. However, after
the manipulation, the ILM returns to the 5th site. It shows the failure of the
mamnipulation in this case. Thus, the shift manipulation of ILM depends on
$T_{\mathrm{C}}$ as in the case of the switching of two coexisting localized
modes. It would be more accurate to say that it depends on how to change the
coupling constant during the manipulation. Then, we investigate
$T_{\mathrm{C}}$ dependency of the manipulation.

To investigate $T_{\mathrm{C}}$ dependency of the shift manipulation, we
implement the same simulation at each $T$ when $T_{\mathrm{C}}$ is varied from
$T$ to $100T$ every $T$. Fig. 7(a) plots the result of the manipulation on the
displacement of the 5th site. Circle-points indicates the success of the
manipulation and cross-points the failure. We found some regions that the
manipulation succeeds and regions that the manipulation fails. However, no
rules are found on the result of manipulation. We focus on a boundary between
regions of the success and the failure. So, we implement the same simulation
at each $T$ when $T_{\mathrm{C}}$ is varied from $22T$ to $25T$ every
$3T/100$. Fig. 7(b) shows the result. We found nested structures regarding the
success and the failure of the manipulation. These structures are found in the
same simulation near $T_{\mathrm{C}}=96$ in Fig. 7(b). Fig. 7(c) shows the
result. There exist the finer nested structures. So, the uncertainty of the
manipulation seems to depend on regions for $T_{\mathrm{C}}$ where it is
implemented.

The shift manipulation of ILM depends on the process itself. Especially, we
cannot shift an ILM definitely near the regions in which the manipulation is
effected by the ambigious structure of the boundary. It implies that we need
to repeat the manipulation until it is achieved.

Figure 6: Time development of energy of each sites during the shift
manipulation of ILM. The unit time on the horizontal axis is $T$ which is the
period of driver. (a) Simulation results when $T_{\mathrm{I}},\
T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set at $250T,\ 20T$, and $250T$
respectively. (b) Simulation results when $T_{\mathrm{I}},\ T_{\mathrm{C}}$,
and $T_{\mathrm{R}}$ are set at $250T,\ 30T$, and $250T$ respectively.  Figure
7: $T_{\mathrm{C}}$ dependency of the shift manipulation of ILM. Results of
the manipulation are plotted on the displacement of 5th site. Circle-points
indicates success of the manipulation and cross-points fails. (a) Result when
$T_{\mathrm{C}}$ is varied from $T$ to $100T$ every $T$. The unit time is $T$.
(b) Result when $T_{\mathrm{C}}$ is varied from $22T$ to $25T$ every $3T/100$.
(c) Result of the same simulation near $T_{\mathrm{C}}=96$

##  6 Conclusion

In this letter, we have examined shift manipulation of ILM in ac driven Klein
Gordon lattice.

At first, we introduced a reduced system of the lattice, which is obtained by
neglecting the influence from outside two sites which are the subject of the
manipulation. We investigated this reduced system using the concept of NNM,
and obtained three NNMs; two stable localized modes corresponding to the
stable ILM and an unstable anti-phase mode corresponding to the unstable Page
mode. We investigated these NNMs’ dependency on the coupling constant and
found that pitchfork bifurcation of NNMs occurs. Based on the bifurcation
diagram, we numerically implement the switching of two coexisting localized
modes by changing the coupling constant adiabatically.

Lastly, we confirmed that the shift manipulation of ILM is achieved by an
adiabatic change of coupling constant between two targeted sites as is the
case in the reduced system. However, it was also revealed that the result of
the manipulation depends on how to change the coupling constant during the
manipulation. Under certain conditions, there are fine structures regard to
whether the manipulation is achieved or not. This is an interesting phenomenon
for the manipulation of ILMs.

In the appendix, we show that this manipulation is also achieved in Klein
Gordon lattice with soft-type anharmonicity. So, it is seems that this
manipulation method is practicable in a wide range of Klein Gordon lattices.

This manipulation have uncertainly described above. Thus, we need to repeat
this manipulation in oder to move ILM definitely. We expect that this
manipulation is able to be used in other nonlinear systems not included in
this paper by following the same schemes as in this paper.

## Appendix A. Shift manipulation of ILM in Klein Gordon lattice with soft-
type potential

Shift manipulation of ILM in Klein Gordon lattice with soft-type potential is
implemented with the same procedure as the manipulation in the lattice with
hard-type anharmonicity. We take Morse potential expressed by the following
equation

$V(z)=\frac{1}{2}(1-\exp(-z))^{2},$ (6)

for the on-site potential with soft-type anharmonicity. Thus, the equation of
motion of the Klein Gordon lattice with dissipation and the periodic external
driver is written as

$\displaystyle
m\ddot{u}_{n}+\gamma\dot{u}_{n}+\exp(-u_{n})-\exp(-2u_{n})-C(u_{n-1}$
$\displaystyle+u_{n}-u_{n+1})$ $\displaystyle=F\cos\omega{t},$ (7)

where the values of each parameter except $\omega$ are the same as in the
lattice with hard-type anharmonicity. Here, we use the in-phase driver to
sustain ILMs. The driving frequency $\omega$ is set at $0.8$ because ILMs have
frequencies lying below the lowest band frequency which is 1 in this case.

There are two kinds of symmetrical ILM which are reverse-phase with the
driver. One is a stable site centered mode, the other is an unstable bond
centered mode. Fig. 8(a) shows the amplitude distribution of the site centered
mode. Fig. 8(b) shows the amplitude distribution of the bond centered mode.
There are additionally the symmetrical site centered and the symmetrical bond
centered mode which are in-phase with the driver. However, we don’t consider
these modes since both of them are unstable.

Figure 8: Amplitude distribution of site centered mode (a) and bond centered
mode (b). Each points indicate the position of particles when they have
largest displacements. Each particles oscillate with the same period.

By reducing the ac driven Klein Gordon lattice to a two-degrees of freedom
system, we obtain the equation of motion of the two-degree of freedom
nonlinear system described by

$\displaystyle\ddot{u}_{1}+\gamma\dot{u}_{1}+\exp(-u_{1})-\exp(-2u_{1})-C(u_{2}-u_{1})=F\cos\omega{t},$
(8)
$\displaystyle\ddot{u}_{2}+\gamma\dot{u}_{2}+\exp(-u_{2})-\exp(-2u_{2})-C(u_{1}-u_{2})=F\cos\omega{t},$
(9)

where subscripts are replaced by 1 and 2 for simplicity.

By using the Newton-Raphson method, we obtain three coexisting NNMs which have
reverse phase with the driver in this system. One is the in-phase mode with
equal amplitude, which is unstable. This mode corresponds to unstable bond
centered mode. The other two modes are localized modes which are the NNMs
localized on each oscillator. These modes correspond to stable site centered
mode. Fig. 9 shows a pitchfork bifurcation of these NNMs when the coupling
constant is varied. When the coupling constant is large, only stable anti-
phase mode exists. However, the in-phase mode loses its stability at $C\simeq
0.132$ and two stable coexisting localized modes appear.

Figure 9: Relationship between coupling constant and amplitude ratio $k$ in
each of the coexisting NNMs. The solid line corresponds to stable NNMs. The
dashed line corresponds to unstable NNMs.

Based on this bifurcation diagram, we implement the switching of two
coexisting localized mode by appropriately tuning the coupling constant.
First, we increase the coupling constant linearly from the initial value
$C=0.07953$ to $C=0.15$, which is greater than the bifurcation point at
$C=0.132$. This process takes $T_{\mathrm{I}}$. Next, we keep the coupling
constant at the value $C=0.15$ during period $T_{\mathrm{C}}$. Lastly, we
reduce the coupling constant to the original value during period
$T_{\mathrm{D}}$.

We switch the mode from NNM localized at oscillator 1 to one localized at
oscillator 2. Figures. 10(a) and (b) show simulation results when
$T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{D}}$ are set at $250T,\
20T$, and $250T$, respectively. These figures shows that switching between
coexisting localized modes is achieved. Figures. 10(c) and (d) show the
simulation results when $T_{\mathrm{I}},\ T_{\mathrm{C}}$, and
$T_{\mathrm{R}}$ are set at $250T,\ 40T$, and $250T$, respectively. The
initial state once becomes the in-phase mode. However, it finally returns to
the original state. These showed that the switching depends on the way of
tuning the coupling constant as in the case of hard-type anharmonicity.

Figure 10: (a) and (b) Simulation results of the switching when
$T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set at $250T,\
20T$, and $250T$ respectively. (c) and (d) Simulation results of the switching
when $T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set at
$250T,\ 40T$, and $250T$ respectively. Points indicates Poincare map of orbits
of each oscillators during the manipulation. Solid line is trajectory for a
period of each oscillators before the switching and dotted line is during
period $T_{\mathrm{C}}$, and dashed line is after the switching.

Based on the results of the reduced system, we move an ILM from the 5th site
to the 6th site. Thus, we change a coupling constant between the 5th and the
6th site adiabatically as in the case of switching between coexisting
localized modes. Figure. 11(a) shows the behavior of ILM during the
manipulation when $T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{D}}$ are
set at $250T,\ 30T$, and $250T$, respectively. In this case, the manipulation
is achieved. Figure. 11(b) shows the behavior of ILM during the manipulation
when $T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{D}}$ are set at
$250T,\ 60T$, and $250T$, respectively. This figure shows the failure of the
manipulation. These results indicate that the manipulation depends on the way
of tuning the coupling constant. Further investigation of $T_{\mathrm{C}}$
dependency reveals that there are the fine structures in some regarding to the
success or the failure of this manipulation. So, the manipulation similarly
determines the result uncertainly in some regions for $T_{\mathrm{C}}$.

Figure 11: Time development of energy of each sites during the shift
manipulation of ILM in the lattice with soft-type anharmonicity. The unit time
on the horizontal axis is $T$ which is the period of driver. (a) Simulation
results when $T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set
at $250T,\ 30T$, and $250T$ respectively. (b) Simulation results when
$T_{\mathrm{I}},\ T_{\mathrm{C}}$, and $T_{\mathrm{R}}$ are set at $250T,\
60T$, and $250T$ respectively.

## References

  * [1] A. J. Sievers, S. Takeno, Intrinsic localized modes in anharmonic crystals, Phys. Rev. Lett. 61 (8) (1988) 970–973.  doi:10.1103/PhysRevLett.61.970.
URL https://link.aps.org/doi/10.1103/PhysRevLett.61.970

  * [2] S. Flach, A. V. Gorbach, Discrete breathers—advances in theory and applications, Physics Reports 467 (1-3) (2008) 1–116.  doi:10.1016/j.physrep.2008.05.002.
URL https://doi.org/10.1016/j.physrep.2008.05.002

  * [3] R. S. MacKay, S. Aubry, Proof of existence of breathers for time-reversible or hamiltonian networks of weakly coupled oscillators, Nonlinearity 7 (6) (1994) 1623–1643.  doi:10.1088/0951-7715/7/6/006.
URL https://doi.org/10.1088/0951-7715/7/6/006

  * [4] S. Flach, Existence of localized excitations in nonlinear hamiltonian lattices, Phys. Rev. E 51 (1995) 1503–1507.  doi:10.1103/PhysRevE.51.1503.
URL https://link.aps.org/doi/10.1103/PhysRevE.51.1503

  * [5] S. Aubry, Breathers in nonlinear lattices: Existence, linear stability and quantization, Physica D: Nonlinear Phenomena 103 (1-4) (1997) 201–250.  doi:10.1016/S0167-2789(96)00261-8.
URL https://doi.org/10.1016/S0167-2789(96)00261-8

  * [6] R. S. MacKay, J. A. Sepulchre, Stability of discrete breathers, Physica D: Nonlinear Phenomena 119 (1-2) (1998) 148–162.  doi:10.1016/S0167-2789(98)00073-6.
URL https://doi.org/10.1016/S0167-2789(98)00073-6

  * [7] D. Chen, S. Aubry, G. P. Tsironis, Breather mobility in discrete $\varphi^{4}$ nonlinear lattices, Physical Review Letters 77 (23) (1996) 4776–4779.  doi:10.1103/PhysRevLett.77.4776.
URL https://link.aps.org/doi/10.1103/PhysRevLett.77.4776

  * [8] S. Aubry, T. Cretegny, Mobility and reactivity of discrete breathers, Physica D: Nonlinear Phenomena 119 (1-2) (1998) 34–46.  doi:10.1016/S0167-2789(98)00062-1.
URL https://doi.org/10.1016/S0167-2789(98)00062-1

  * [9] E. Trias, J. J. Mazo, T. P. Orlando, Discrete breathers in nonlinear lattices: Experimental detection in a josephson array, Physical Review Letters 84 (4) (2000) 741–744.  doi:10.1103/PhysRevLett.84.741.
URL https://link.aps.org/doi/10.1103/PhysRevLett.84.741

  * [10] H. S. Eisenberg, Y. Silberberg, R. Morandotti, A. R. Boyd, J. S. Aitchison, Discrete spatial optical solitons in waveguide arrays, Physical Review Letters 81 (16) (1998) 3383–3386.  doi:10.1103/PhysRevLett.81.3383.
URL https://link.aps.org/doi/10.1103/PhysRevLett.81.3383

  * [11] M. Sato, B. E. Hubbard, A. J. Sievers, B. Ilic, D. A. Czaplewski, H. G. Craighead, Observation of locked intrinsic localized vibrational modes in a micromechanical oscillator array, Physical Review Letters 90 (4) (2003) 044102.  doi:10.1103/PhysRevLett.90.044102.
URL https://link.aps.org/doi/10.1103/PhysRevLett.90.044102

  * [12] R. Stearrett, L. Q. English, Experimental generation of intrinsic localized modes in a discrete electrical transmission line, Journal of Physics D: Applied Physics 40 (17) (2007) 5394.  doi:10.1088/0022-3727/40/17/058.
URL https://doi.org/10.1088/0022-3727/40/17/058

  * [13] K. Narahara, Dissipative discrete breathers in series-connected tunnel diode oscillator lattice, Journal of the Physical Society of Japan 89 (7) (2020) 074005.  doi:10.7566/JPSJ.89.074005.
URL https://doi.org/10.7566/JPSJ.89.074005

  * [14] A. V. Savin, Y. S. Kivshar, Discrete breathers in carbon nanotubes, EPL (Europhysics Letters) 82 (6) (2008) 66002.  doi:10.1209/0295-5075/82/66002.
URL https://doi.org/10.1209/0295-5075/82/66002

  * [15] L. Z. Khadeeva, S. V. Dmitriev, Discrete breathers in crystals with nacl structure, Phys. Rev. B 81 (2010) 214306.  doi:10.1103/PhysRevB.81.214306.
URL https://link.aps.org/doi/10.1103/PhysRevB.81.214306

  * [16] M. Haas, V. Hizhnyakov, A. Shelkan, M. Klopov, A. J. Sievers, Prediction of high-frequency intrinsic localized modes in ni and nb, Phys. Rev. B 84 (2011) 144303.  doi:10.1103/PhysRevB.84.144303.
URL https://link.aps.org/doi/10.1103/PhysRevB.84.144303

  * [17] L. Z. Khadeeva, S. V. Dmitriev, Y. S. Kivshar, Discrete breathers in deformed graphene, JETP letters 94 (7) (2011) 539–543.  doi:10.1134/S0021364011190106.
URL https://doi.org/10.1134/S0021364011190106

  * [18] R. T. Murzaev, D. V. Bachurin, E. A. Korznikova, S. V. Dmitriev, Localized vibrational modes in diamond, Physics Letters A 381 (11) (2017) 1003–1008.  doi:10.1016/j.physleta.2017.01.014.
URL https://doi.org/10.1016/j.physleta.2017.01.014

  * [19] M. Sato, S. Yasui, M. Kimura, T. Hikihara, A. J. Sievers, Management of localized energy in discrete nonlinear transmission lines, EPL (Europhysics Letters) 80 (3) (2007) 30002.  doi:10.1209/0295-5075/80/30002.
URL https://doi.org/10.1209/0295-5075/80/30002

  * [20] L. Jin, J. Mei, L. Li, Analysis of intrinsic localised mode for a new energy harvesting cantilever array, The European Physical Journal Applied Physics 66 (1) (2014) 10902.  doi:10.1051/epjap/2014130565.
URL https://doi.org/10.1051/epjap/2014130565

  * [21] M. Kimura, T. Hikihara, Stability change of intrinsic localized mode in finite nonlinear coupled oscillators, Physics Letters A 372 (25) (2008) 4592–4595.  doi:10.1016/j.physleta.2008.04.054.
URL https://doi.org/10.1016/j.physleta.2008.04.054

  * [22] M. Kimura, T. Hikihara, Coupled cantilever array with tunable on-site nonlinearity and observation of localized oscillations, Physics Letters A 373 (14) (2009) 1257–1260.  doi:10.1016/j.physleta.2009.02.005.
URL https://doi.org/10.1016/j.physleta.2009.02.005

  * [23] M. Kimura, T. Hikihara, Capture and release of traveling intrinsic localized mode in coupled cantilever array, Chaos: An Interdisciplinary Journal of Nonlinear Science 19 (1) (2009) 013138\.  doi:10.1063/1.3097068.
URL https://doi.org/10.1063/1.3097068

  * [24] T. Cretegny, S. Aubry, S. Flach, 1d phonon scattering by discrete breathers, Physica D: Nonlinear Phenomena 119 (1-2) (1998) 73–87.  doi:10.1016/S0167-2789(98)00066-9,.
  * [25] F. Hadipour, D. Saadatmand, M. Ashhadi, A. M. Marjaneh, I. Evazzade, A. Askari, S. V. Dmitriev, Interaction of phonons with discrete breathers in one-dimensional chain with tunable type of anharmonicity, Physics Letters A 384 (4) (2020) 126100.  doi:10.1016/j.physleta.2019.126100.
URL https://doi.org/10.1016/j.physleta.2019.126100

  * [26] G. P. Tsironis, A. R. Bishop, A. V. Savin, A. V. Zolotaryuk, Dependence of thermal conductivity on discrete breathers in lattices, Physical Review E 60 (6) (1999) 6610–6613.  doi:10.1103/PhysRevE.60.6610.
URL https://link.aps.org/doi/10.1103/PhysRevE.60.6610

  * [27] M. Sato, B. E. Hubbard, A. J. Sievers, Colloquium: Nonlinear energy localization and its manipulation in micromechanical oscillator arrays, Reviews of Modern Physics 78 (1) (2006) 137–157.  doi:10.1103/RevModPhys.78.137.
URL https://link.aps.org/doi/10.1103/RevModPhys.78.137

  * [28] M. Sato, B. E. Hubbard, A. J. Sievers, B. Ilic, H. G. Craighead, Optical manipulation of intrinsic localized vibrational energy in cantilever arrays, EPL (Europhysics Letters) 66 (3) (2004) 318–323.  doi:10.1209/epl/i2003-10224-x.
URL https://doi.org/10.1209/epl/i2003-10224-x

  * [29] O. Bang, P. D. Miller, Exploiting discreteness for switching in waveguide arrays, Optics letters 21 (15) (1996) 1105–1107.  doi:10.1364/OL.21.001105.
URL https://doi.org/10.1364/OL.21.001105

  * [30] R. A. Vicencio, M. I. Molina, Y. S. Kivshar, Controlled switching of discrete solitons in waveguide arrays, Optics letters 28 (20) (2003) 1942–1944.  doi:10.1364/OL.28.001942.
URL https://doi.org/10.1364/OL.28.001942

  * [31] R. A. Vicencio, M. I. Molina, Y. S. Kivshar, Switching of discrete optical solitons in engineered waveguide arrays, Physical Review E 70 (2) (2004) 026602.  doi:10.1103/PhysRevE.70.026602.
URL https://doi.org/10.1103/PhysRevE.70.026602

  * [32] Y. S. Kivshar, D. K. Campbell, Peierls-nabarro potential barrier for highly localized nonlinear modes, Physical Review E 48 (4) (1993) 3077.  doi:10.1103/PhysRevE.48.3077.
URL https://doi.org/10.1103/PhysRevE.48.3077

  * [33] M. J. M. Öster, A. Eriksson, Enhanced mobility of strongly localized modes in waveguide arrays by inversion of stability, Physical Review E 67 (5) (2003) 056606.  doi:10.1103/PhysRevE.67.056606.
URL https://doi.org/10.1103/PhysRevE.67.056606

  * [34] L. Hadžievski, A. Maluckov, M. Stepić, D. Kip, Power controlled soliton stability and steering in lattices with saturable nonlinearity, Physical review letters 93 (3) (2004) 033901.  doi:10.1103/PhysRevLett.93.033901.
URL https://doi.org/10.1103/PhysRevLett.93.033901

  * [35] R. A. Vicencio, M. Johansson, Discrete soliton mobility in two-dimensional waveguide arrays with saturable nonlinearity, Physical Review E 73 (4) (2006) 046602.  doi:10.1103/PhysRevE.73.046602.
URL https://doi.org/10.1103/PhysRevE.73.046602

  * [36] J. B. Page, Asymptotic solutions for localized vibrational modes in strongly anharmonic periodic systems, Phys. Rev. B 41 (11) (1990) 7835.  doi:10.1103/PhysRevB.41.7835.
URL https://link.aps.org/doi/10.1103/PhysRevB.41.7835

  * [37] R. M. Rosenberg, Normal modes of nonlinear dual-mode systems, J. Appl. Mech. 27 (1960) 263–268.  doi:10.1115/1.3643948.
URL https://doi.org/10.1115/1.3643948

  * [38] G. Kerschen, M. Peeters, J. C. Golinval, A. F. Vakakis, Nonlinear normal modes, part I: A useful framework for the structural dynamicist, Mechanical Systems and Signal Processing 23 (1) (2009) 170–194.  doi:10.1016/j.ymssp.2008.04.002.
URL https://doi.org/10.1016/j.ymssp.2008.04.002

  * [39] M. Peeters, R. Viguié, G. Sérandour, G. Kerschen, J. C. Golinval, Nonlinear normal modes, part II: toward a practical computation using numerical continuation techniques, Mechanical systems and signal processing 23 (1) (2009) 195–216.  doi:10.1016/j.ymssp.2008.04.003.
URL https://doi.org/10.1016/j.ymssp.2008.04.003

  * [40] J. L. Marin, S. Aubry, Breathers in nonlinear lattices: Numerical calculation from the anticontinuous limit, Nonlinearity 9 (6) (1996) 1501–1528.  doi:10.1088/0951-7715/9/6/007.
URL https://doi.org/10.1088/0951-7715/9/6/007


# Generalized Capsule Networks with
Trainable Routing Procedure

Zhenhua Chen
School of School of Informatics,
Computing, and Engineering
Indiana University
Bloomington, IN 47408, USA
<EMAIL_ADDRESS>
&David Crandall
School of School of Informatics,
Computing, and Engineering
Indiana University
Bloomington, IN 47408, USA
<EMAIL_ADDRESS>

###### Abstract

CapsNet (Capsule Network) was first proposed by Sabour et al. (2017) and later
another version of CapsNet was proposed by Hinton et al. (2018). CapsNet has
been proved effective in modeling spatial features with much fewer parameters.
However, the routing procedures in both papers are not well incorporated into
the whole training process. The optimal number of routing procedure is misery
which has to be found manually. To overcome this disadvantages of current
routing procedures in CapsNet, we embed the routing procedure into the
optimization procedure with all other parameters in neural networks, namely,
make coupling coefficients in the routing procedure become completely
trainable. We call it Generalized CapsNet (G-CapsNet). We implement both
“full-connected” version of G-CapsNet and “convolutional” version of
G-CapsNet. G-CapsNet achieves a similar performance in the dataset MNIST as in
the original papers. We also test two capsule packing method (cross feature
maps or with feature maps) from previous convolutional layers and see no
evident difference. Besides, we also explored possibility of stacking multiple
capsule layers. The code is shared on CAFFE-CapsNet.

##  1 Introduction

According to Sabour et al. (2017), a capsule is a group of neurons whose
activity vector represents the instantiation parameters of a specific type of
entity such as an object or an object part. Intuitively, CapsNet can better
model spatial relationship by using much more fewer parameters. The
experiments in Sabour et al. (2017) proved this argument in some small
datasets like MNIST and smallNORB . However, the routing procedure is very
computationally expensive, and the number of routing iterations has to be
manually set by testing. To overcome this issue, we propose Generalized
CapseNet. The key idea of CapsNet is incorporating the routing procedure to
the overall optimization procedure. In other words, it makes the coupling
coefficients trainable instead of being calculated by dynamic routing Sabour
et al. (2017) or EM routing Hinton et al. (2018).

Another interesting question yet to answer is that how to package a capsule
from activations of previous convolutional layers. We can select the elements
if each capsule at the same position across different feature maps (as Figure
2 Left shows). We can also select the elements of each capsule within each
feature map (for example, we can choose each row or column as a capsule as
Figure 2 Right shows). Two methods seem different, and the previous one makes
more sense since capsules are supposed to capture different spatial features.
However, according to our experiment, both ways have similar performance. We
do not know the exact reason for the question but will try to answer it in the
experiment session.

CapsNet seems more general than standard neural networks (neurons can be
considered as a unit length of capsules), so we also explore the scalability
of CapsNet by building a two-layer G-CapsNet. Unfortunately, the network is
inclined to get saturated even after a capsule version of ReLU layer is added.
The same thing happens when we try to ‘capsulize’ the whole network, namely
capsules other than neurons are the atomic units. For example, a colored input
image can be considered as $h\times w$, 3-dimension capsules, so are the
following layers.

##  2 Generalized CapsNet

To better illustrate the idea of G-Capsnet, let’s check what the normal neural
networks’ loss function looks like. Assume we have a dataset
$\\{(x^{(1)},y^{(1)}),\ldots,(x^{(m)},y^{(m)})\\}$. If we can use neural
networks to train a model, the loss function can be defined as Equation 1:

$J(W,b)=\left[\frac{1}{m}\sum_{i=1}^{m}\left(\frac{1}{2}\left\|h_{W,b}(x^{(i)})-y^{(i)}\right\|^{2}\right)\right]+\frac{\lambda}{2}\sum_{l=1}^{n_{l}-1}\;\sum_{i=1}^{s_{l}}\;\sum_{j=1}^{s_{l+1}}\left(W^{(l)}_{ji}\right)^{2}$
(1)

The loss function for G-CapsNet is similar. The only difference is we have
extra coupling coefficients (through routing procedure) to train. Note that
$\left\|h_{W,b}(x^{(i)})\right\|$ is the mathematical form of the neural
network that we assume would fit the dataset. $n_{l}$, $s_{l}$, $s_{l+1}$ are
the number of layers, the number of neurons of layer $l$ and the number of
neurons of layer $l+1$. $W$ and $b$ are the parameters we are supposed to get
by training the neural network.

For the generalized CapsNet, $W^{(l)}_{ji}$ is the transformation matrix that
maps one type of capsule in one layer to another type of capsule on top of it.
$c^{(l)}_{ji}$ is the coupling coefficient between adjacent layers.

$J(W,b,c)=\left[\frac{1}{m}\sum_{i=1}^{m}\left(\frac{1}{2}\left\|h_{W,b,c}(x^{(i)})-y^{(i)}\right\|^{2}\right)\right]+\frac{\lambda}{2}\sum_{l=1}^{n_{l}-1}\;\left(\sum_{i=1}^{s_{l}}\;\sum_{j=1}^{s_{l+1}}\left(W^{(l)}_{ji}\right)^{2}+\sum_{i=1}^{p_{l}}\;\sum_{j=1}^{p_{l+1}}\left(c^{(l)}_{ji}\right)^{2}\right)$
(2)

###  2.1 Structure of Generalized CapsNet

Like the structure of the previous CapsNets, the generalized CapsNet also
includes two phases, capsule transformation, and capsule routing, as Figure 1
shows. Capsule transformation converse one type of capsules into another type
of capsules. For example, Sabour et al. (2017) transforms 8-dimension capsules
into 16-dimension capsules while Hinton et al. (2018) transforms $4\times
4$capsules into $4\times 4$ capsules. In theory, we can transform type of
capsules into any other type of capsules. The capsules can be of any shape
(vector, matrix, cube or even hyper-cube). Assume the shape of capsules
$\mathbf{u}$ in the lower layer is $(m_{1},m_{2},\dots,m_{k})$, the
transformation matrix $\mathbf{W}$ is $m_{k},n_{1},n_{2},\dots,n_{k}$, then
the shape of output capsules $\mathbf{v}$ in the higher layer is
$m_{1},m_{2},\dots,m_{k},n_{1},n_{2},\dots,n_{k}$. Capsule routing ensures
capsules in lower layers are scaled and sent to their parent capsules in
higher layers, as Equation 3 shows.

$\mathbf{v_{j}}=\sum_{i}c_{ij}\mathbf{u_{j|i}},\quad\quad\mathbf{u_{j|i}}=\mathbf{W_{ij}}\mathbf{u_{i}}$
(3)

Figure 1: The structure of generalized CapsNet.

###  2.2 Squash Function

We adopt two versions of squash functions which are proposed by Sabour et al.
(2017) and Edgar et al. (2017). As Equation 4 shows. Squash functions here
equal to a normalization operation which also bring non linearity.

$\mathbf{v_{j}^{{}^{\prime}}}=\frac{\mathbf{\|v_{j}\|}^{2}}{1+\mathbf{\|v_{j}\|}^{2}}\frac{\mathbf{v_{j}}}{\|\mathbf{v_{j}}\|}\quad\quad\mathbf{v_{j}^{{}^{\prime}}}=\left(1-\frac{1}{\mathbf{e}^{\|\mathbf{v_{j}}\|}}\right)\frac{\mathbf{v_{j}}}{\|\mathbf{v_{j}}\|}$
(4)

###  2.3 Loss Function

We adopt the same margin loss function in Sabour et al. (2017), as Equation 5
shows.

$L_{k}=T_{k}\max(0,m^{+}-\|\mathbf{v_{k}}\|)^{2}+\lambda(1-T_{k})\max(0,\|\mathbf{v_{k}}\|-m^{-})^{2}$
(5)

##  3 Experiments

###  3.1 full connected G-CapsNet on MNIST

We adopt the same baseline as in paper . The first convolutional layer outputs
256 feature maps. The second convolutional layer outputs 256 feature maps or
32$\times$6$\times$6 8D capsules. The last two layers are fully connected
layers. Please check our released code or the original paper written by Sabour
et al. (2017) for more details.

We call the Capsule structure that Sabour et al. (2017) proposed full
connected CapsNet since each capsule in the higher layer connects to every
capsule in the lower layer. As Table 1 shows, no matter whether the
reconstruction procedure involved, G-CapsNet can always achieve better
performance by using much less number of parameters. Note that the performance
of the baseline and the G-CapsNet reported here is a little lower than in
Sabour et al. (2017), we consider the difference is caused by different deep
learning frame (We use Caffe other than TensorFlow).

Figure 2: Packaging capsules across feature maps VS within feature maps

Another interesting thing we found is that the way of stacking the activation
of neurons does not matter a lot (0.68 VS 0.66). It makes more sense to
package capsules across different feature maps since the capsules are supposed
to capture a couple of different types of features at the same position, as
Figure 2 left shows. However, what we found is that packaging capsules within
each feature map (as Figure 2 right shows) as opposed to across feature maps
can achieve similar performance. The potential reason for this is that no
matter how we package capsules, once the organization of these capsules is
fixed, the network will finally learn potential spatial relationship.

###  3.2 Convolutional G-CapsNet on MNIST

Similarly, we call the structure in Hinton et al. (2018) as convolutional
CapsNet since the same type of capsules of different positions shares the same
transformation matrices. The convolutional G-CapsNet we use here is similar to
the previously mentioned full connected G-CapsNet except we use 6 by 6 kernel
for calculating the and 4 by 4 “matrix” capsule instead of 1 by 16 “vector”
capsule for the capsule layer. Please check our project for more details. As
Table 1 shows, convolutional G-CapsNet can achieve better performance compared
to the baseline by using fewer parameters.

Above G-CapsNets all have a single capsule layer, then a natural question
arises, can we build a Multi-layer G-CapsNet? After all, neural networks can
be considered as a special type of CapsNet (the length of each capsule is 1).
One straightforward way is stacking multiple capsule layers on top of each
other. However, we found this type of CapsNet is easy to get saturated. To
solve this issue, we design a capsule version of “ReLU” layer which makes the
performance better but still far from satisfying. How to make a CapsNet more
scalable will be our next work.

Algorithm | error rate(%) | number of parameters
---|---|---
baseline | 0.83 | 35.4M
full connected G-CapsNet | 0.66 | 8.2M
full connected G-CapsNet* | 0.66 | 6.8M
full connected G-CapsNet** | 0.76 | 8.2M
Convolutional G-CapsNet | 0.75 | 6.9M
Convolutional G-CapsNet* | 0.70 | 5.5M
Table 1: Error rate VS number of parameters on MNIST. Note that “*” means no
reconstruction, “**” means the alternate way of stacking capsules.

##  4 Conclusion

G-CapsNet incorporates the routing procedure of capsules into the whole
optimization process which gets rid of the set of routing times and guarantees
the convergence. Also, we tried two different ways of packaging capsules and
conclude that how to packaging capsules does not matter a lot. Finally, we
evaluated muli-layer CapsNets and found that multi-layer CapsNets are easy to
get saturated. How to make CapsNet scalable is still an open question.

## References

  * Edgar et al. (2017) Xi Edgar, Bing Selina, and Jin Yang.  Capsule network performance on complex data.  2017\.  URL https://arxiv.org/abs/1712.03480.
  * Hinton et al. (2018) Geoffrey E Hinton, Sara Sabour, and Nicholas Frosst.  Matrix capsules with EM routing.  In _International Conference on Learning Representations_ , 2018.  URL https://openreview.net/forum?id=HJWLfGWRb.
  * Sabour et al. (2017) Sara Sabour, Nicholas Frosst, and Geoffrey E. Hinton.  Dynamic routing between capsules.  _CoRR_ , abs/1710.09829, 2017.  URL http://arxiv.org/abs/1710.09829.


# A Brief Account of
Klein’s Icosahedral Extensions

Leonardo Solanilla, Erick S. Barreto and Viviana Morales Departamento de
Matemáticas y Estadística, Universidad del Tolima, Barrio Santa Elena, Ibagué,
Tolima, Colombia<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(February 2021)

We present an alternative relatively easy way to understand and determine the
zeros of a quintic whose Galois group is isomorphic to the group of rotational
symmetries of a regular icosahedron. The extensive algebraic procedures of
Klein in his famous Vorlesungen über das Ikosaeder und die Auflösung der
Gleichungen vom fünften Grade are here shortened via Heymann’s theory of
transformations. Also, we give a complete explanation of the so-called
icosahedral equation and its solution in terms of Gaussian hypergeometric
functions. As an innovative element, we construct this solution by using
algebraic transformations of hypergeometric series. Within this framework, we
develop a practical algorithm to compute the zeros of the quintic.

Keywords (MSC 2020): Computational methods for problems pertaining to field
theory (12-08); Polynomials in real and complex fields: location of zeros
(algebraic theorems) (12D10); Real polynomials: location of zeros (26C10);
Zeros of polynomials, rational functions, and other analytic functions of one
complex variable (30C15); Numerical computation of solutions to single
equations (65H05).

##  1 Introduction

Let $f\in K[X]$ be an irreducible quintic with coefficients in a field
$K\leq\mathbb{C}$. Let $L$ be the splitting field of $f$, $K\leq
L\leq\mathbb{C}$, and let $\textrm{Gal}(K\leq L)$ denote the Galois group of
$L$ over $K$. It is a well-known fact that $f$ is solvable by radicals if and
only if $\textrm{Gal}(K\leq L)$ is a subgroup of the Frobenius group $F_{20}$,
i.e. if and only if $\textrm{Gal}(K\leq L)$ is isomorphic (up to conjugacy) to
$F_{20}$ of order 20, to the dihedral group $D_{10}$ of order 10, or to the
cyclic group $\mathbb{Z}/5\mathbb{Z}$. Since the Galois group of certain
quintics is $A_{5}$ (isomorphic to the group of symmetries of a regular
icosahedron) and this group is not solvable, the solutions to this type of
quintic cannot be expressed by field operations and radicals from the
polynomial coefficients. We will say that $L$ is an icosahedral extension of
$K$ if $\textrm{Gal}(K\leq L)=A_{5}$. We assume the reader already knows the
basics of Galois Theory and quintics, as they are given e.g. in Cox [5] and
the earliest works of Abel [1] and Ruffini [12].

In this paper we give a method to solve any quintic whose splitting field is
icosahedral. Furthermore, we also provide a numerical implementation of the
method. Our approach is based on Klein’s [10] foundational paper, as well as
on the more recent renditions of Heymann [8], Shurman [13] and Nash [11]. For
the numerical part, we have been inspired by Cox, Little & O’Shea [6] and
Trott [15]. Some incomplete sketchy elementary ideas on the subject can be
found in [2] and [3].

Tschirnhaus transformations [16] are key ingredients of Klein’s method. As
usual, given a general quintic
$x^{5}+c_{4}x^{4}+c_{3}x^{3}+c_{2}x^{2}+c_{1}x+c_{0}$, the first degree
transformation $z=x+(c_{4}/5)$ eliminates the term of 4th degree and reduces
the problem to a quintic $z^{5}+pz^{3}+qz^{2}+rz+s$. After a further quadratic
transformation $y=z^{2}-az-b$, we obtain the so-called principal quintic

$y^{5}+5\alpha y^{2}+5\beta y+\gamma=0.$ (1)

The quadratic transformation extends the coefficient field $K$ to $K(\delta)$,
where
$\delta^{2}=\displaystyle\frac{9q^{2}}{4p^{2}}+\displaystyle\frac{3p}{5}\displaystyle\displaystyle-\frac{2r}{p}$.
For simplicity’s sake, we will denote this new field by the same letter $K$.

From a purely field-theoretical approach, our main task consists of proving
the following theorem, that we will call Klein’s theorem. C. f. Slodowy [14].

######  Theorem 1.

Let $K$ be a subfield of the complex numbers containing $\delta$ and the fifth
roots of unity. We suppose that $L$ is a Galois extension of $K$ with
$\textup{Gal}(K\leq L)=A_{5}$, the alternating group of five letters. Then,
there exists a $J\in K$ such that the solution $Y$ of the icosahedral equation
$q(Y)=J$ (Section 2 below) yields $L=K(Y)$, that is, $L$ is obtained from $K$
by adjoining $Y$.

Notably, along the proof we will develop a practical numerical method for
finding the zeros of our principal quintic.

Since Klein’s theorem essentially asserts that the solution of the quintic
reduces to the solution of the icosahedral equation, in Section 2 we introduce
this equation together with a method of solution. Gaussian hypergeometric
functions provide the means to determine such a solution. Our approach is
comprehensive and we provide full details. Although this solution has been
described barely in Nash [11], we use a simpler procedure based on the
algebraic transformations of the hypergeometric series. We devote Section 3 to
the algebraic treatment of the principal quintic. We have closely followed the
theory of Heymann’s [8] resolvents. They allow us to show the way the solution
of the quintic effectively amounts to the solution of the icosahedral
equation. With this, Klein’s theorem is proved. In Section 4 we rewrite the
previous sections as a useful algorithm that produces as outputs the five
solutions of an icosahedral quintic. In the end, we draw some conclusions
regarding the main features of our method.

##  2 Solution of the icosahedral equation

The solution of a Klein’s quintic relies heavily on the geometry of the
icosahedron.

### Icosahedral symmetry

The icosahedron is the surface of the regular polyhedron or Platonic solid
having 12 vertices, 30 edges and 20 triangular faces (the ancient Greek number
$\epsilon\iota\kappa o\sigma\iota$, eíkosi, means twenty). Once it is
inscribed in a sphere, we project radially from the center to obtain the image
of the icosahedron on the sphere. Then, we project the circumsphere onto the
complex plane via the stereographic projection to get a plane image of the
icosahedral surface. As it is usual, the extended complex plane or Riemann
sphere $\mathbb{S}$ is conveniently identified with the complex projective
space $\mathbb{P}^{1}$.

The icosahedral rotation group $\mathcal{I}$ has order 60 and is isomorphic to
$A_{5}$, the group of the even permutations of five letters. It comprises all
spherical transformations or rotations under which the icosahedron remains
invariant. By virtue of the identification $\mathbb{S}\simeq\mathbb{P}^{1}$,
the elements of $\mathcal{I}$ can be understood as Möbius transformations or,
alternatively, as homogeneous forms. As well, $A_{5}$ is isomorphic to the
projective special linear group $\textrm{PSL}(2,5)$. The action of the
icosahedral group on the sphere yields three abnormal orbits corresponding to
the vertices (12 points), edge midpoints (30 points) and face centers (20
points). The rest of the orbits are normal in the sense that the group acts
transitively on them and so, they all consist of 60 points. The icosahedral
equation arises from the problem of inverting the branched covering
$\mathbb{P}^{1}\rightarrow\mathbb{P}^{1}/\mathcal{I}$.

### Invariant forms of 60th degree

Klein (1884) describes an inverse of this covering in relation to polynomials
or homogeneous forms. For instance, the vertices of the icosahedron are
precisely the zeros of the invariant polynomial

$f=f(z,w)=zw(z^{10}+11z^{5}w^{5}-w^{10}).$ (2)

Similarly, the Hessian determinant of $f$ produces the invariant form

$H=H(z,w)=-(z^{20}+w^{20})+228(z^{15}w^{5}-z^{5}w^{15})-494z^{10}w^{10},$ (3)

which vanishes at the centers of the icosahedron faces. Also, the Jacobian
determinant of $f,H$ provides the 30th degree invariant

$T=T(z,w)=(z^{30}+w^{30})+522(z^{25}w^{5}-z^{5}w^{25})-10005(z^{20}w^{10}+z^{10}w^{20}),$
(4)

vanishing at the midpoints of the edges. We should standardize these invariant
forms in the 60th degree in order to deal with the spherical orbits under the
action of $\mathcal{I}$. Once we do so, the set $\mathbb{O}$ of the invariant
forms of 60th degree is easily furnished with a structure of vector space of
dimension two over $\mathbb{C}$. Among other choices, the forms $H^{3}$ and
$f^{5}$ form a basis for $\mathbb{O}$. Remarkably enough,
$T^{2}=1728f^{5}-H^{3}$ in this basis.

It is not hard to establish an identification between the homogeneous
polynomials in $\mathbb{O}$ with the orbits in $\mathbb{P}^{1}/\mathcal{I}$,
where they vanish.

### The isomorphism $\mathbb{P}^{1}\simeq\mathbb{P}^{1}/\mathcal{I}$

The canonical projection $\mathbb{P}^{1}\rightarrow\mathbb{P}^{1}/\mathcal{I}$
takes a point $(z,w)$ to its orbit $[z,w]$. It corresponds to a map
$\mathbb{P}^{1}\rightarrow\mathbb{O}$, taking $(z,w)$ to the form $P(z,w)$,
which equals zero at $(z,w)$. An arbitrary $P\in\mathbb{O}$ is written $P=\mu
H^{3}-\lambda f^{5}$, in the basis $H^{3},f^{5}$. Then, $P=0$ implies

$\frac{H^{3}}{f^{5}}=\frac{\lambda}{\mu}.$

So, the covering can be written $(z,w)\mapsto(\lambda,\mu)$. In non-
homogeneous coordinates we have
$\left(Y=\frac{z}{w},1\right)\mapsto\left(\frac{\lambda}{\mu}=J,1\right)$, or
simply $Y\mapsto J$.

The inversion of the covering reveals interesting facts. First, we must
restrict ourselves to one of the possible branches or restricted domains.
Then, we construct the inverse map $[z,w]\mapsto(z,w)$ by means of
$\left(H^{3}(z,w),f^{5}(z,w)\right)\mapsto(z,w)$ and thus $J\mapsto Y$. In the
subset of the normal orbits $[z,w]\in\mathbb{P}^{1}/\mathcal{I}$, it is
possible to choose a unique pair of holomorphic functions $(z,w)$. However,
the abnormal orbits introduce discontinuities corresponding to the vertices,
edge midpoints and face midpoints. The nature of this inverse map is next
elucidated with the aid of Complex Analysis.

### Differential resolvents

Even more, we can construct an explicit inverse map $J\mapsto Y$ of the
covering. For sure, $H^{3}$ and $f^{5}$ are connected through variable $J$ and
the constant rank theorem reduces the problem to the equivalent form

$\displaystyle f(z,w)$ $\displaystyle=$ $\displaystyle k,$ $\displaystyle
H(z,w)$ $\displaystyle=$ $\displaystyle u,$

where $k$ is constant and $u$ is a complex variable. We recall that the
jacobian determinant of the map $\mathbb{P}^{1}\rightarrow\mathbb{P}^{1}$,
$(z,w)\mapsto(k,u)$, equals $-20T$, where $T$ is the icosahedral form (4).
Also, the derivatives of $z,w$ with respect to $u$ can be computed implicitly
through

$\left(\begin{array}[]{l}\frac{\displaystyle dz}{\displaystyle du}\\\
\frac{\displaystyle dw}{\displaystyle
du}\end{array}\right)=-\frac{1}{20T}\left(\begin{array}[]{rr}\frac{\displaystyle\partial
H}{\displaystyle\partial w}&-\frac{\displaystyle\partial
f}{\displaystyle\partial w}\\\ -\frac{\displaystyle\partial
H}{\displaystyle\partial z}&\frac{\displaystyle\partial
f}{\displaystyle\partial
z}\end{array}\right)\left(\begin{array}[]{l}\displaystyle 0\\\ \displaystyle
1\end{array}\right).$

Since

$\frac{d}{du}\left(\frac{\partial f}{\partial
z}\right)=\frac{\partial^{2}f}{\partial
z^{2}}\cdot\frac{dz}{du}+\frac{\partial^{2}f}{\partial z\partial
w}\cdot\frac{dw}{du}\quad\textrm{and}\quad\frac{d}{du}\left(\frac{\partial
f}{\partial w}\right)=\frac{\partial^{2}f}{\partial
w^{2}}\cdot\frac{dw}{du}+\frac{\partial^{2}f}{\partial z\partial
w}\cdot\frac{dz}{du},$

this yields the linear ordinary homogeneous equations

$T^{2}\frac{d^{2}z}{du^{2}}+T\frac{dT}{du}\cdot\frac{dz}{du}+\frac{11Hz}{400}=T^{2}\frac{d^{2}w}{du^{2}}+T\frac{dT}{du}\cdot\frac{dw}{du}+\frac{11Hw}{400}=0.$
(5)

By using $T^{2}=12^{3}f^{5}-H^{3}$ and setting $u=H$, we obtain

$(12^{3}f^{5}-H^{3})\frac{d^{2}z}{dH^{2}}-\frac{3}{2}H^{2}\frac{dz}{dH}+\frac{11Hz}{400}=0$
(6)

and a similar equation for $w$. Then, we suitably recover
$J=H^{3}/12^{3}f^{5}$ and perform the change of variables $H\mapsto J$ to get

$\frac{dz}{dH}=\frac{3J}{H}\cdot\frac{dz}{dJ}\quad\textrm{and}\quad\frac{d^{2}z}{dH^{2}}=\frac{9J^{2}}{H^{2}}\cdot\frac{d^{2}z}{dJ^{2}}+\frac{2}{H}\cdot\frac{dz}{dH}.$

The substitution of these derivatives in (6) produces the well-known
hypergeometric equation

$J(1-J)\frac{d^{2}z}{dJ^{2}}+[c-(a+b+1)J]\frac{dz}{dJ}-abz=0,\quad\textrm{with}\quad
a=\frac{11}{60},\ b=-\frac{1}{60},\ c=\frac{2}{3},$ (7)

and a corresponding equation for $w$. This equation has three regular
singularities located respectively at $J=0,1,\infty$. They correspond to the
branch points of the covering: $J=0$ means $H=0$, $J=1$ means $T=0$ and
$J=\infty$ means $f=0$. A general exposition of the nature of the solutions to
general hypergeometric equations similar to (7) near ordinary points and
regular singularities can be found in Copson [4, Chapter X] and Hille [9,
Section 6.1].

The solutions to (7) can be expressed in terms of the hypergeometric series

$F(a,b,c;J)=1+\sum_{n=1}^{\infty}\frac{(a)_{n}(b)_{n}}{(1)_{n}(c)_{n}}J^{n},$

where $(p)_{n}=p(p+1)\cdots(p+n-1)$. For the values of $a,b,c$ in (7), there
are fundamental systems $z,w$ of solutions at the singular points
$J=0,1,\infty$. In particular, at $J=\infty$, since $a-b$ is not an integer,

$z(J)=\frac{1}{J^{\frac{11}{60}}}F\left(\frac{11}{60},\frac{31}{60},\frac{6}{5},\frac{1}{J}\right),\quad
w(J)=J^{\frac{1}{60}}F\left(-\frac{1}{60},\frac{19}{60},\frac{4}{5},\frac{1}{J}\right)$

are linearly independent for $|J|>1$. It is a known fact that the functions
$z,w$ have an analytical continuation in the whole plane, save for the
singularities.

### Solving for $Y$

Inasmuch as our main interest lies in the ratio of functions $z,w$ and not in
the actual solutions to the differential equation, we consider the function

$s=s(J)=\frac{z(J)}{w(J)}=\frac{F\left(\frac{11}{60},\frac{31}{60},\frac{6}{5},\frac{1}{J}\right)}{J^{\frac{1}{5}}~{}F\left(-\frac{1}{60},\frac{19}{60},\frac{4}{5},\frac{1}{J}\right)}.$
(8)

With it, the solution to the icosahedral equation is accomplished by composing
$s$ with a Möbius transformation

$\mathfrak{M}(s)=\mathfrak{M}(s(J))=Y(J)=Y.$

The spherical symmetry $\mathfrak{M}$ rotates $s$ to the correct value of $Y$,
which satisfies the correct “initial or boundary conditions” at three
convenient points. We can determine $\mathfrak{M}$ assisted by Table 1.

$J$ | $0$ | $1$ |  $\infty$
---|---|---|---
$s=s(J)$ | $\infty$ | $-(-1)^{\frac{1}{10}}\sqrt[5]{1728}$ |  $0$
$Y=\mathfrak{M}(s)$ | $\infty$ | $-(-1)^{\frac{1}{10}}$ |  $0$
Table 1: Determining $Y$.

The first row contains three different values of $J$. The second row is
calculated from the first row by using the formula defining $s(J)$ and some
hypergeometric identities that we will state shortly. The values of $Y$ in the
third row are obtained from $J=H^{3}/1728f^{5}$. For instance, $J=\infty$
implies $f=0$ and this is achieved through $f(Y_{\infty}=0,1)=0$. For the case
$J=0$ we realize that, at the infinity,

$\frac{1}{J}=\frac{H^{3}}{1728f^{5}}(1,Y_{0}=\infty)=\infty.$

Finally, an easy calculation proves that

$J=\frac{H^{3}}{1728f^{5}}\left({\displaystyle-(-1)^{\frac{1}{10}}},1\right)=1.$

The values of $s(\infty),s(0)$ are easily determined. $s(1)$ is found with the
aid of the following algebraic transformations of the Gaussian hypergeometric
function, cf. Vidunas [17, Section 6.3]:

$\displaystyle
F\left(\frac{11}{60},\frac{31}{60},\frac{6}{5},\varphi_{1}(x)\right)$
$\displaystyle=\frac{(1-228x+494x^{2}+228x^{3}+x^{4})^{\frac{11}{20}}}{1+11x-x^{2}},$
(9) $\displaystyle
F\left(-\frac{1}{60},\frac{19}{60},\frac{4}{5},\varphi_{1}(x)\right)$
$\displaystyle=(1-228x+494x^{2}+228x^{3}+x^{4})^{-\frac{1}{20}},$ (10)

where

$\varphi_{1}(x)=\frac{1728x(x^{2}-11x-1)^{5}}{(1-228x+494x^{2}+228x^{3}+x^{4})^{3}}.$
(11)

The quotient of (9) into (10) provides a feasible expression for (8):

$s(J)=\frac{F\left(\frac{11}{60},\frac{31}{60},\frac{6}{5},\varphi_{1}(x)\right)}{J^{\frac{1}{5}}~{}F\left(-\frac{1}{60},\frac{19}{60},\frac{4}{5},\varphi_{1}(x)\right)}=\frac{(1-228x+494x^{2}+228x^{3}+x^{4})^{\frac{3}{5}}}{J^{\frac{1}{5}}(1+11x-x^{2})}.$
(12)

Thus, we may take at convenience $x=i$ to obtain $\varphi_{1}(i)=1=1/J$, that
is $J=1$. Therefore,

$\displaystyle
s(1)=\frac{F\left(\frac{11}{60},\frac{31}{60},\frac{6}{5},\varphi_{1}(i)\right)}{(1)^{\frac{1}{5}}~{}F\left(-\frac{1}{60},\frac{19}{60},\frac{4}{5},\varphi_{1}(i)\right)}$
$\displaystyle=$
$\displaystyle\frac{(1-228(i)+494(i)^{2}+228(i)^{3}+(i)^{4})^{\frac{3}{5}}}{(1)^{\frac{1}{5}}(1+11(i)-(i)^{2})}$
$\displaystyle=$ $\displaystyle-(-1)^{\frac{1}{10}}\sqrt[5]{1728}.$

Finally, we write

$Y=\mathfrak{M}(s)=\frac{\mathfrak{a}s+\mathfrak{b}}{\mathfrak{c}s+\mathfrak{d}},$

for some $\mathfrak{a},\mathfrak{b},\mathfrak{c},\mathfrak{d}\in\mathbb{C}$.
Since $s(\infty)=0$ yields $Y=0$, we must have $\mathfrak{b}=0$ and
$\mathfrak{d}\neq 0$. Without loss of generality, we set $\mathfrak{a}=1$.
Then, $s(0)=\infty$ and $Y=\infty$ imply $1/\mathfrak{c}=\infty$ and so,
$\mathfrak{c}=0$. This leaves us with $s(1)/\mathfrak{d}=-(-1)^{\frac{1}{10}}$
and an elementary calculation gives $\mathfrak{d}=\sqrt[5]{1728}$. Hence, $Y$
is given by

$Y=\mathfrak{M}(s)=\frac{1}{\sqrt[5]{1728}}\frac{z(J)}{w(J)}=\frac{F\left(\frac{11}{60},\frac{31}{60},\frac{6}{5},\frac{1}{J}\right)}{\sqrt[5]{1728\cdot
J}\times F\left(-\frac{1}{60},\frac{19}{60},\frac{4}{5},\frac{1}{J}\right)}.$
(13)

This non-homogeneous coordinate solves the problem.

##  3 Proof of Klein’s Theorem

### Heymann’s resolvents

The coefficients $\alpha,\beta,\gamma$ of the principal quintic $y^{5}+5\alpha
y^{2}+5\beta y+\gamma=0$ belong to a quadratic Tschirnhaus extension $K$ of
the original coefficient field. We remind that

$\Delta=3125\cdot(108\alpha^{5}\gamma-135\alpha^{4}\beta^{2}+90\alpha^{2}\beta\gamma^{2}-320\alpha\beta^{3}\gamma+256\beta^{5}+\gamma^{4})$

is the so-called discriminant of the quintic. By hypothesis, $K$ also contains
the fifth roots of unity and $\sqrt{\Delta}$, because $\textrm{Gal}(K\leq
L)=A_{5}$.

First and foremost, the immemorial symmetric relation $yz=y+z$ transforms this
quintic into

$(1+5\alpha+5\beta+\gamma)z^{5}-5(3\alpha+4\beta+\gamma)z^{4}+5(3\alpha+6\beta+2\gamma)z^{3}-5(\alpha+4\beta+2\gamma)z^{2}+5(\beta+\gamma)z-\gamma.$

The fourth and third degree terms vanish whenever $\alpha=\gamma/3$ and
$\beta=-\gamma/2$ (the actual solution does not have to fulfill these strident
conditions). By renaming $y=\eta_{1}$, $z=\eta_{2}$, $h_{1}=-6\gamma^{-1}$ and
$h_{2}=1-h_{1}$, we obtain Heymann’s $\eta$-resolvents

$\displaystyle h_{1}\eta_{1}^{5}-10\eta_{1}^{2}+15\eta_{1}-6=0,$
$\displaystyle h_{2}\eta_{2}^{5}-10\eta_{2}^{2}+15\eta_{2}-6=0.$

We observe that $h_{1}+h_{2}=1$ and $\eta_{1}\eta_{2}=\eta_{1}+\eta_{2}$. Many
higher degree polynomials in $\eta_{1},\eta_{2}$ are easily reduced to
quadratic or linear polynomials by these identities. For example,
$h_{1}\eta_{1}^{4}\eta_{2}=10\eta_{1}+\eta_{2}-6$. Such polynomials are called
Heymann’s simultaneous resolvents.

### The first form of the solution

We are interested in solutions of the form

$y=p\eta_{1}+q\eta_{2},\quad p,q\in\mathbb{C}.$ (14)

The substitution of this expression in (1) yields the simultaneous resolvent

$M\eta_{1}^{2}+N\eta_{2}^{2}+P\eta_{1}+Q\eta_{2}+R,$

where $M,N,P,Q$ and $R$ are given by

$\displaystyle M$ $\displaystyle=$ $\displaystyle
5p^{2}(2h_{2}p^{3}+2h_{1}q^{3}+\alpha h_{1}h_{2}),$ $\displaystyle N$
$\displaystyle=$ $\displaystyle 5q^{2}(2h_{1}q^{3}+2h_{2}p^{3}+\alpha
h_{1}h_{2}),$ $\displaystyle P$ $\displaystyle=$ $\displaystyle
5p(-3h_{2}p^{4}+10h_{2}p^{3}q+6h_{1}pq^{3}+h_{1}q^{4}+2\alpha
h_{1}h_{2}q+\beta h_{1}h_{2}),$ $\displaystyle Q$ $\displaystyle=$
$\displaystyle 5q(-3h_{1}q^{4}+10h_{1}q^{3}p+6h_{2}qp^{3}+h_{2}p^{4}+2\alpha
h_{1}h_{2}p+\beta h_{1}h_{2}),$ $\displaystyle R$ $\displaystyle=$
$\displaystyle
6(h_{2}p^{5}-5h_{2}p^{4}q+10h_{2}q^{2}p^{3}+10h_{1}p^{2}q^{3}-5h_{1}pq^{4}+h_{1}q^{5})+\gamma
h_{1}h_{2}.$

In order to find values of the coefficients $\alpha,\beta,\gamma$, we require
$M=N=R=0$. With this,

$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle-2(p^{3}h_{1}^{-1}+q^{3}h_{2}^{-1}),$ $\displaystyle\beta$
$\displaystyle=$ $\displaystyle
3(p^{3}h_{1}^{-1}(p-2q)-q^{3}h_{2}^{-1}(2p-q)),$ $\displaystyle\gamma$
$\displaystyle=$
$\displaystyle-6(p^{3}h_{1}^{-1}(p^{2}-5pq+10q^{2})+q^{3}h_{2}^{-1}(10p^{2}-5pq+q^{2})).$

These formulas provide a map $(p,q,h_{1})\mapsto(\alpha,\beta,\gamma)$. We
also need certain map $(\alpha,\beta,\gamma)\mapsto(p,q,h_{1})$. To do so, the
previous equations for $\alpha$ and $\beta$ produce

$h_{1}=-\frac{18p^{3}(p-q)}{3\alpha(2p-q)-2\beta},\quad
h_{2}=-\frac{18q^{3}(p-q)}{3\alpha(p-2q)+2\beta}.$

By replacing these expressions in the equation for $\gamma$,

$\displaystyle 12\alpha r+6\beta s-\gamma=0,$ (15)

where we have introduced the new variables $r=(p-q)^{2}$, $s=p+q$ (this
variable $s$ is different from variable $s(J)$ in Section 2). Clearly, $r$ and
$s$ depend on each other. Now we replace $h_{1},h_{2}$ in the fundamental
relation $h_{1}+h_{2}=1$ and use (15) to obtain the central quadratic relation

$(\alpha^{4}+\alpha\beta\gamma-\beta^{3})(12r)^{2}-(2\alpha^{3}\gamma+11\alpha^{2}\beta^{2}+\beta\gamma^{2})(12r)+(\alpha\gamma-8\beta^{2})^{2}=0.$
(16)

We notice that the quadratic discriminant needed to solve for $r$ is

$108\alpha^{5}\gamma-135\alpha^{4}\beta^{2}+90\alpha^{2}\beta\gamma^{2}-320\alpha\beta^{3}\gamma+256\beta^{5}+\gamma^{4}=\frac{\Delta}{3125}.$

Therefore, the solution of (16) does not extend our field $K$. Each value of
$r$ provides a value of $s$ and so, values of $p$ and $q$. Heymann’s
parameters $h_{1},h_{2}$ are the solutions of the product-sum or trinomial-
factoring relations $h_{1}+h_{2}=1$ and

$h_{1}h_{2}=\frac{432\beta
p^{3}q^{3}}{12(\alpha\gamma-\beta^{2})r-\gamma^{2}}.$ (17)

### Eulerian resolvents and icosahedral invariants forms

Heymann’s ideas recall some results of Euler (1764) on solvable quintics of
the form (1). His solutions have the form

$y=\epsilon z-\epsilon^{2}w,$

where $\epsilon$ denotes a fifth root of unity. In particular, $\epsilon=1$
yields $y=z-w$. Also, $z,w$ are related to the polynomial coefficients through

$\alpha=-zw^{2},\quad\beta=z^{3}w,\quad\gamma=-(z^{5}-w^{5}).$

Euler’s solution is not our solution. However, we can use the map
$(\alpha,\beta,\gamma)\mapsto(p,q,h_{1})$ and give it a shot. First, we notice
that the second degree coefficient of equation (16) vanishes with these
expressions of $\alpha,\beta,\gamma$. The remaining linear equation gives

$r=\frac{w^{4}(7z^{5}+w^{5})^{2}}{12f},$

where $f$ is the well-known icosahedral invariant form (2). From (15) it is
easy to find the corresponding value

$s=\frac{z^{3}(-z^{10}+39z^{5}w^{5}+26w^{10})}{6f}.$

With $r$ and $s$ we determine $p$ and $q$. Then, from (17) we get

$4h_{1}h_{2}=\frac{H^{3}}{12^{3}f^{5}},$ (18)

where $H$ is the expression in (3). In other words, $J=4h_{1}h_{2}$ is the
icosahedral parameter.

Now, from $h_{1}+h_{2}=1$ and (18) we obtain

$h_{1}-h_{2}=\frac{\sqrt{12^{3}f^{5}-H^{3}}}{24f^{2}\sqrt{3f}}=\frac{T}{24f^{2}\sqrt{3f}},$
(19)

where $T$ is just but the invariant (4).

The smooth way to compute the parameters $h_{1},h_{2}$ is through $r$ and $s$:

$\displaystyle
h_{1}=-\frac{9(s+\sqrt{r})^{3}\sqrt{r}}{6\alpha(s+3\sqrt{r})-8\beta},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}h_{2}=\frac{9(s-\sqrt{r})^{3}\sqrt{r}}{6\alpha(s-3\sqrt{r})-8\beta}.$

Henceforth, by virtue of (19),

$h_{1}-h_{2}=\kappa(r)\sqrt{r}=\frac{T}{24f^{2}\sqrt{3f}},$ (20)

where $\kappa=\kappa(r)$ is a rational function of $r$ and so,

$\sqrt{3rf}=\frac{T}{24\kappa f}.$

Since $h_{1}-h_{2}\in K$, neither $\sqrt{r}$ nor $\sqrt{3f}$ introduces extra
elements into field $K$.

### The second form of the solution

The main idea to find the solutions of the quintic runs as follows. From
coefficients $\alpha,\beta,\gamma$ of (1), we find $r,s$. With them, we
determine $p,q$ together with $h_{1},h_{2}$. Then, we pose the problem by
forcing the solution to have also Euler’s form but not its original meaning.
That is,

$y=p\eta_{1}+q\eta_{2}=z-w,\quad\eta_{1}\eta_{2}=\eta_{1}+\eta_{2}.$ (21)

The order pair $(z,w)$ represents homogeneous coordinates of $\mathbb{P}^{1}$
and should be found. The resolution for Heymann’s resolvents gives

$\eta_{1}=-\frac{2\sqrt{3f}}{t-\sqrt{3f}}\quad\textrm{and}\quad\eta_{2}=\frac{2\sqrt{3f}}{t+\sqrt{3f}},$

where $t$ is the well-konwn octahedral form

$t(z,w)=z^{6}+2z^{5}w-5z^{4}w^{2}-5z^{2}w^{4}-2zw^{5}+w^{6}.$

In this way, the quadratic term $\eta_{1}\eta_{2}$ in (21) does not introduce
new elements to field $K$. Further, the solution $y$ to the reduced quintic is
obtained by substituting $\eta_{1},\eta_{2}$ in system (21), thus

$y=-\displaystyle\frac{6sf+2t\sqrt{3rf}}{t^{2}-3f}.$

We observe that $r,s$ are fixed by the coefficients of the quintic and $f$ is
invariant under the action of the icosahedral group $\mathcal{I}$. The only
quantity that changes is $t$. Certainly, the action of $\mathcal{I}$ yields
the five values

$t_{\nu}=\epsilon^{3\nu}z^{6}+2\epsilon^{2\nu}z^{5}w-5\epsilon^{\nu}z^{4}w^{2}-5\epsilon^{4\nu}z^{2}w^{4}-2\epsilon^{3\nu}zw^{5}+\epsilon^{2\nu}w^{6},$
(22)

in which $\epsilon=e^{2\pi i/5}\in K$ and $\nu=0,1,2,3,4$. In brief, the
solution set is

$\left\\{y_{\nu}=-\displaystyle\frac{6sf+2t_{\nu}\sqrt{3rf}}{t_{\nu}^{2}-3f}:\nu=0,1,2,3,4\right\\}.$
(23)

Now, from $\alpha,\beta,\gamma$ and the values $p,q$, we determine the
icosahedral parameter

$J=4h_{1}h_{2}=\frac{H^{3}(z,w)}{12^{3}f^{5}(z,w)}.$ (24)

With this, the problem is reduced to determine a pair $(z,w)$ from an orbit
$[z,w]\in\mathbb{P}^{1}/\mathcal{I}$, i.e. a solution of the icosahedral
equation. The corresponding $Y=Y(J)$ is determined as in Section 2. Then, we
replace $f(Y,1),H(Y,1)$ and $t_{\nu}(Y,1)$ in (23) to find the zeros of the
principal quintic.

In such a way, Klein’s theorem has been proved.

##  4 Algorithm

The foregoing procedure forms a finite sequence of computer-implementable
instructions.

Read : coefficients $\alpha,\beta,\gamma$ of a reduced quintic with Galois
group $A_{5}$;

Function  _f(_$z$_)_:

Data: $z$, a complex number

Result: inhomogeneous icosahedral form $f(z)=f(z,1)$

$f(z)\leftarrow z(z^{10}+11z^{5}-1)$;

return ;

Function  _t(_$z$_)_:

input: $z$, a complex number

output: inhomogeneous octahedral forms $t_{\nu}(z)=t_{\nu}(z,1)$, a vector of
five components

$\epsilon\leftarrow\exp\left(\frac{2\pi i}{5}\right)$;

for  _$\nu\leftarrow 0$ to $4$_ do

$t_{\nu}(z)\leftarrow\epsilon^{3\nu}z^{6}+2\epsilon^{2\nu}z^{5}-5\epsilon^{\nu}z^{4}-5\epsilon^{4\nu}z^{2}-2\epsilon^{3\nu}z+\epsilon^{2\nu}$;

end for

return ;

Function  _F(_$a,b,c;z$_)_:

Data: argument $z$; parameters $a,b,c$

Result: Gaussian hypergeometric function $F(a,b,c;z)=_{2}F_{1}(a,b,c;z)$

$F(a,b,c;z)\leftarrow\displaystyle
1+\sum_{n=1}^{\infty}\frac{(a)_{n}(b)_{n}}{(1)_{n}(c)_{n}}z^{n}$;

return ;

/* initialization */

$r\leftarrow$roots of
$(\alpha^{4}+\alpha\beta\gamma-\beta^{3})(12r)^{2}-(2\alpha^{3}\gamma+11\alpha^{2}\beta^{2}+\beta\gamma^{2})(12r)+(\alpha\gamma-8\beta^{2})^{2}=0$;

$s\leftarrow$solutions of $12\alpha r+6\beta s-\gamma=0$;

$p\leftarrow\frac{1}{2}\left(\sqrt{r}+s\right)$;
$q\leftarrow\frac{1}{2}\left(s-\sqrt{r}\right)$;

$J\leftarrow\displaystyle 4h_{1}h_{2}=4\frac{432\beta
p^{3}q^{3}}{12(\alpha\gamma-\beta^{2})r-\gamma^{2}}$;

$Y\leftarrow\displaystyle\frac{F\left(\frac{11}{60},\frac{31}{60},\frac{6}{5},\frac{1}{J}\right)}{\sqrt[5]{1728\cdot
J}\times F\left(-\frac{1}{60},\frac{19}{60},\frac{4}{5},\frac{1}{J}\right)}$;

/* computation of zeros */

for  _$\nu\leftarrow 0$ to $4$_ do

$y_{\nu}\leftarrow\displaystyle-\frac{6sf(Y)+2t_{\nu}(Y)\sqrt{3rf(Y)}}{t_{\nu}^{2}(Y)-3f(Y)}$;

end for

Write : zeros $y_{\nu}$, $\nu=0,1,2,3,4$, of the quintic.

### Example

We consider the following quintic with coefficients in
$\mathbb{Q}[i,\sqrt{\Delta}]$:

$y^{5}+5iy^{2}-12y+(1-i)=0.$

We know its Galois group is $A_{5}$. We have
$\alpha=i,\beta=-\frac{12}{5},\gamma=1-i$. Firstly, we conveniently pick the
solution $r=-0.140712-1.06363i$ of

$(12r)^{2}\left(\alpha^{4}+\alpha\beta\gamma-\beta^{3}\right)-(12r)\left(2\alpha^{3}\gamma+11\alpha^{2}\beta^{2}+\beta\gamma^{2}\right)+\left(\alpha\gamma-8\beta^{2}\right)^{2}=0.$

The corresponding value $s=0.816914-0.0478157i$ is obtained from $12\alpha
r+6\beta s-\gamma=0$. Next, $p=\frac{1}{2}\left(\sqrt{r}+s\right)$ and
$q=\frac{1}{2}\left(s-\sqrt{r}\right)$ provide

$~{}~{}~{}~{}~{}p=0.749812\,-0.413396i~{}~{}~{}~{}~{}~{}\textrm{and}~{}~{}~{}~{}~{}~{}q=0.0671022\,+0.365581i.~{}~{}~{}~{}$

With all this, we obtain immediately

$J=4h_{1}h_{2}=\frac{H^{3}}{1728f^{5}}=-0.324158-2.04659i.$

A solution of this icosahedral equation is

$\displaystyle Y=Y(-0.324158-2.04659i)$ $\displaystyle=$
$\displaystyle\frac{{}_{2}F_{1}(\frac{11}{60},\frac{31}{60};\frac{6}{5};J^{-1})}{\sqrt[5]{1728J}~{}_{2}F_{1}(-\frac{1}{60},\frac{19}{60};\frac{4}{5};J^{-1})}$
$\displaystyle=$ $\displaystyle 0.178352\,+0.0718131i.$

Afterwards, we calculate the values $f(Y)$ and $t_{\nu}(Y)$:

$\displaystyle f$ $\displaystyle=-0.178721-0.0713975i,$ $\displaystyle t_{0}$
$\displaystyle=0.509555\,-0.278001i,$ $\displaystyle t_{1}$
$\displaystyle=-0.761539+0.997924i,$ $\displaystyle t_{2}$
$\displaystyle=0.372515\,-1.14707i,$ $\displaystyle t_{3}$
$\displaystyle=0.240993\,+1.27692i,$ $\displaystyle t_{4}$
$\displaystyle=-0.361523-0.849771i.$

Finally, we substitute these values in

$y_{\nu}=-\frac{6sf(Y)+2t_{\nu}(Y)\sqrt{3rf(Y)}}{t_{\nu}^{2}(Y)-3f(Y)}$

to get the desired solutions

$\displaystyle y_{0}$ $\displaystyle=0.0895118\,-0.0828539i,$ $\displaystyle
y_{1}$ $\displaystyle=-0.0120031+2.20094i,$ $\displaystyle y_{2}$
$\displaystyle=-0.0430531-1.43083i,$ $\displaystyle y_{3}$
$\displaystyle=-1.90456-0.333135i,$ $\displaystyle y_{4}$
$\displaystyle=1.87011\,-0.354121i.$

This example has been developed with the routines of the program Wolfram
Mathematica, version 12.0.0.0.

##  5 Concluding remarks

We have constructed a simple algorithm for the evaluation of the zeros of any
quintic polynomial whose splitting field is icosahedral. The mathematical
machinery used to construct the method is not that arduous and more abridged
than the usual approaches employed to compute such zeros. In particular, the
theory of Heymann [8] permits to rapidly establish a central relation between
the solutions to the equation and the icosahedral invariants. Also, the
algebraic transformations of hypergeometric series furnish an alternative
simplifying way to solve the icosahedral equation.

Some other interesting questions arise from these results. They will be
reported at a later date.

## Acknowledgments

This research was partially funded by the Comité Central de Investigaciones,
Universidad del Tolima, Ibagué, Colombia, grant number 60120. We also thank
the Facultad de Ciencias, Universidad del Tolima, for its logistic support to
craft the manuscript.

## References

## References

  * [1] N. H. Abel (1824) Mémoire sur les équations algébriques où on démontre l’impossibilité de la résolution de l’équation générale du cinquième dégré. Christiania: Groendahl.
  * [2] C. Almario, J. Rodríguez (2017) Algunas propiedades geométricas del icosaedro. Ibagué, Colombia: Universidad del Tolima, Facultad de Ciencias, report.
  * [3] E. S. Barreto, V. Morales (2020) El grupo del icosaedro y la solubilidad de la quíntica. Ibagué, Colombia: Universidad del Tolima, Facultad de Ciencias de la Educación, report.
  * [4] E. T. Copson (1972) An Introduction the Theory of Functions of a Complex Variable. Oxford at the Clarendon Press. The first edition dates back to 1935.
  * [5] A. Cox (2012) Galois Theory. Hoboken, New Jersey: John Wiley & Sons.
  * [6] A. Cox, J. Little & D. O’Shea (1997) Ideals, Varieties and Algorithms. New York: Springer.
  * [7] L. Euler (1764) De resolutione aequationum cuiusvis gradus. Novi commentarii academiae scientiarum Petropolitanae 9, 70-98. English translation: H. J. Stevens & T. C. Stevens (2012) On the solution of equations of any degree. The Euler Archive http://eulerarchive.maa.org/, retrieved on 02.17.2020.
  * [8] W. Heymann (1897) Die Transformation und Auflösung der Gleichung fünften Grades in elementarer Darstellung. Zeitschrift für Mathematik und Physik 42. Jahrg. 1897. 2. Heft. Spanish translation: L. Solanilla (2019) Transformación y resolución de la ecuación de quinto grado: una presentación elemental. Ibagué, Colombia: Universidad del Tolima, Departamento de Matemáticas y Estadística.
  * [9] E. Hille (1976) Ordinary Differential Equations in the Complex Domain. New York: John Wiley & Sons.
  * [10] F. Klein (1884) Vorlesungen über das Ikosaeder und die Auflösung der Gleichungen vom fünften Grade. Leipzig: Teubner. English translation: G. G. Morrice (1956) Lectures on the Icosahedron and the solution of equations of fifth degree. New York: Dover Publications.
  * [11] O. Nash (2014) On Klein’s Icosahedral Solution of the Quintic. Expositiones Mathematicae, 32, 99–120.
  * [12] P. Ruffini (1799) Teoria generale delle equazioni: in cui si dimostra impossibile la soluzione algebraicha delle equazioni generali di grado superiori al quarto. Bologna: Nella Stamperia di S. Tommaso d’Aquino.
  * [13] J. Shurman (1997) Geometry of the Quintic. New York: John Wiley & Sons.
  * [14] P. Slodowy (1986) Das Ikosaeder und die Gleichungen fünften Grades. In Arithmetik und Geometrie. Vier Vorlesungen, 71-112. Basel: Springer.
  * [15] M. Trott (2006) The Mathematica GuideBook for Symbolics. New York: Springer.
  * [16] E. Tschirnhaus (1683) Methodus auferendi omnes terminos intermedios ex data aequatione. Acta Eruditorum, 204-207. English translation: R. F. Green (2003) A method for removing all intermediate terms from a given aquation. ACM SIGSAM Bulletin, 37 1.
  * [17] R. Vidunas (2018) Transformations of algebraic Gauss hypergeometric functions. http://arxiv.org/abs/0807.4808v1, [math.CA].

Conflicts of interest: The authors declare no conflict of interest. The
founding sponsors had no role in the design of the study, in the analyses, in
the writing of the manuscript and in the decision to publish the results.