# A Hyperbolic-to-Hyperbolic Graph Convolutional Network

Jindou Dai  Yuwei Wu Corresponding author Zhi Gao  and Yunde Jia
Laboratory of Intelligent Information Technology  School of Computer Science
Beijing Institute of Technology (BIT)  Beijing  100081  China.
{daijindou, wuyuwei, gaozhi_2017<EMAIL_ADDRESS>
###### Abstract

Hyperbolic graph convolutional networks (GCNs) demonstrate powerful
representation ability to model graphs with hierarchical structure. Existing
hyperbolic GCNs resort to tangent spaces to realize graph convolution on
hyperbolic manifolds, which is inferior because tangent space is only a local
approximation of a manifold. In this paper, we propose a hyperbolic-to-
hyperbolic graph convolutional network (H2H-GCN) that directly works on
hyperbolic manifolds. Specifically, we developed a manifold-preserving graph
convolution that consists of a hyperbolic feature transformation and a
hyperbolic neighborhood aggregation. The hyperbolic feature transformation
works as linear transformation on hyperbolic manifolds. It ensures the
transformed node representations still lie on the hyperbolic manifold by
imposing the orthogonal constraint on the transformation sub-matrix. The
hyperbolic neighborhood aggregation updates each node representation via the
Einstein midpoint. The H2H-GCN avoids the distortion caused by tangent space
approximations and keeps the global hyperbolic structure. Extensive
experiments show that the H2H-GCN achieves substantial improvements on the
link prediction, node classification, and graph classification tasks.

(a) Existing hyperbolic GCNs.

(b) The proposed H2H-GCN.

Figure 1:  Comparisons of exisiting hyperbolic GCNs and the proposed H2H-GCN.
At the $\ell$-th layer, (a) existing hyperbolic GCNs performs Euclidean graph
convolutional opeartions, _e.g_., feature transformation and neighborhood
aggregation, in the tangent space $\mathcal{T}_{{\bm{x}}}\mathcal{L}$ that is
a local approximation of the hyperbolic manifold $\mathcal{L}$; (b) H2H-GCN
directly performs a hyperbolic feature transformation and a hyperbolic
neighborhood aggregation on the hyperbolic manifold to learn node
representations, keeping the global hyperbolic structure.

##  1 Introduction

Graph convolutional networks (GCNs) have attracted increasing attention for
graph representation learning, where nodes in a graph are typically embedded
into Euclidean spaces [23, 46, 44, 40, 18, 19]. Several works reveal that many
graphs, such as social networks and citation networks, exhibit a highly
hierarchical structure [11, 24, 30]. Recent studies have shown that hyperbolic
spaces can well capture such hierarchical structure compared to Euclidean
spaces [28, 29, 12]. Different from Euclidean spaces with zero curvature,
hyperbolic spaces possess a constant negative curvature, which allows for an
exponential growth of space volume with radius. This property of hyperbolic
spaces pretty meets the requirements of hierarchical data (_e.g_., trees) that
need an exponential amount of space for branching, and encourages the
development of GCNs in hyperbolic spaces to capture the hierarchical structure
underlying graphs.

Existing hyperbolic GCNs [25, 10, 47] resort to tangent spaces to realize
graph convolution in hyperbolic spaces. Since the hyperbolic space is a
Riemannian manifold rather than a vector space, basic operations (such as
matrix-vector multiplication and vector addition) well defined in Euclidean
spaces are not applicable in hyperbolic space. To generalize graph convolution
to the hyperbolic space, the works in [25, 10, 47] first flatten a hyperbolic
manifold, and then apply Euclidean graph convolutional operations in the
tangent space. The results are projected back to the hyperbolic manifold. The
procedures follow a manifold-tangent-manifold scheme, as shown in Figure 1(a).
These methods has promoted the development of GCNs in hyperbolic spaces and
achieved good performance. However, the mapping between the manifold and the
tangent space is only locally diffeomorphic, which may distort the global
structure of the hyperbolic manifold, especially frequently using tangent
space approximations [21, 37].

In this paper, we propose to design a hyperbolic GCN that directly works on
the hyperbolic manifold to keep global hyperbolic structure, rather than
relying on the tangent space. This requires that each step of graph
convolution, _e.g_., feature transformation and neighborhood aggregation,
satisfies a manifold-to-manifold principle. To this end, we present a
hyperbolic-to-hyperbolic graph convolutional network (H2H-GCN), where graph
convolutional operations are directly conducted on the hyperbolic manifold.

Specifically, we developed a manifold-preserving graph convolution consisting
of a hyperbolic feature transformation and a hyperbolic neighborhood
aggregation. The hyperbolic feature transformation plays the role of linear
transformation on hyperbolic manifolds, which requires multiplication of node
representations by a transformation matrix. We constrain the transformation
matrix to be a block diagonal matrix composed of a scalar 1 and an orthogonal
matrix to ensure the transformed node representations still reside on the
hyperbolic manifold. For hyperbolic neighborhood aggregation, we adopt the
Einstein midpoint as the weighted message of neighbor nodes to update a node
representation. Figure 1(b) depicts that H2H-GCN directly carries out the two
steps on hyperbolic manifolds. In constrast to existing hyperbolic GCNs, the
proposed H2H-GCN can avoid the distortion caused by tangent space
approximations and keep the global hyperbolic structure underlying graphs. We
summarize the contributions of this paper as follows.

  * •

We propose a hyperbolic-to-hyperbolic graph convolutional network that
directly performs graph convolution on hyperbolic manifolds, keeping the
global hyperbolic structure underlying graphs. To the best of our knowledge,
this is the first hyperbolic GCN without relying on tangent spaces.

  * •

We developed a hyperbolic feature transformation that is a linear
transformation on hyperbolic manifolds. The manifold constraint on the
transformed hyperbolic representations is ensured by imposing the orthogonal
constraint on the transformation sub-matrix.

##  2 Related Work

GCNs generalize classical convolutional neural networks to graph domains. To
realize the convolution on graphs, there are two types of GCNs. Spectral-based
GCNs [7, 13, 23, 20] are based on the convolution theorem to perform
convolution by transforming graph signals into the spectral domain via the
graph Fourier transform. Spatial-based GCNs [18, 19, 39, 31, 40, 44] update
node representations by aggregating the message from its neighbor nodes, just
like applying convolutional kernel on a local image patch. Despite a solid
theoretical foundation of spectral-based GCNs, spatial-based GCNs have shown
more superiorities due to efficiency, generality and flexibility.

Researchers discoverd that many graphs, _e.g_., social networks and biological
networks, usually exhibit a highly hierarchical structure [24, 30]. Krioukov
et al. [24] pointed that the properties of strong clustering and power-law
degree distribution in such graphs can be explained as a hidden hierarchy.
Recent works have demonstrated powerful representation ability of hyperbolic
spaces to model hierarchies that underlie taxonomies [28, 29], knowledge
graphs [36, 2], images [22], semantic classes [26], actions [27], etc [35, 9,
42], achieving promising performance. Liu et al. [25] and Chami et al. [10]
proposed hyperbolic GCNs that extend GCNs to hyperbolic spaces to capture the
hierarchy underlying graphs. The main difference with our work is that they
perform Euclidean graph convolutional operations in the tangent space,
following a manifold-tangent-manifold scheme. The proposed H2H-GCN developed a
hyperbolic graph convolution in the hyperbolic space without relying on
tangent spaces. We designed a Lorentz linear transformation for feature
transformation on the Lorentz model, and adopted Einstein midpoint to
calculate manifold statistics [fletcher2013geodesic] as aggregation function.
We claim that such a manifold-to-manifold learning principle can avoid the
distortion caused by tangent space approximations and keep the global
hyperbolic structure, that is beneficial to graph representation learning.

##  3 Preliminaries

###  3.1 Hyperbolic Spaces

A Riemannian manifold $(\mathcal{M},g)$ is a differentiable manifold
$\mathcal{M}$ equipped with a metric tensor $g$. It can be locally
approximated to a linear Euclidean space at an arbitrary point
${\bm{x}}\in\mathcal{M}$, and the approximated space is termed as a tangent
space $\mathcal{T}_{{\bm{x}}}\mathcal{M}$. Hyperbolic spaces are smooth
Riemannian manifolds with a constant negative curvature [4]. There are five
isometric models of hyperbolic spaces: the Lorentz model (a.k.a the hyperbolid
model), the Klein model, the Jemisphere model, the Poincar$\mathrm{\acute{e}}$
ball model, and the Poincar$\mathrm{\acute{e}}$ half-space model [8]. In this
paper, we choose the Lorentz model due to its numerical stablity [29].

Formally, the Lorentz model of an $n$-dimensional hyperbolic space is defined
by the manifold
$\mathcal{L}=\\{{\bm{x}}=[x_{0},x_{1},\cdots,x_{n}]\in\mathbb{R}^{n+1}:\langle{\bm{x}},{\bm{x}}\rangle_{\mathcal{L}}=-1,{\bm{x}}_{0}>0\\}$
endowed with the metric tensor $g=\mathrm{diag}([-1,{\bm{1}}_{n}^{\top}])$
where $\mathrm{diag}(\cdot)$ function transforms a vector to a diagonal
matrix. The Lorentz inner product induced by $g$ is defined as

$\langle{\bm{x}},{\bm{y}}\rangle_{\mathcal{L}}={\bm{x}}^{\top}g{\bm{y}}=-{x}_{0}{y}_{0}+\sum_{i=1}^{n}{x}_{i}{y}_{i}.$
(1)

In the following, we describe necessary operations.

Distance. The distance on a manifold is termed as a geodesic that is commonly
a curve representing the shortest path between two nodes. For
$\forall{\bm{x}},{\bm{y}}\in\mathcal{L}$, the distance between them is given
by

$d_{\mathcal{L}}({\bm{x}},{\bm{y}})=\mathrm{arcosh}(-\langle{\bm{x}},{\bm{y}}\rangle_{\mathcal{L}}),\vspace{-2pt}$
(2)

where $\mathrm{arcosh}(\cdot)$ is the inverse hyperbolic cosine function.

Exponential and logarithmic maps. An exponential map
$\mathrm{exp}_{{\bm{x}}}({\bm{v}})$ is the function projecting a tangent
vector ${\bm{v}}\in\mathcal{T}_{{\bm{x}}}\mathcal{M}$ onto $\mathcal{M}$. A
logarithmic map projects vectors on the manifold back to the tangent space
satisfying
$\mathrm{log}_{{\bm{x}}}(\mathrm{exp}_{{\bm{x}}}({\bm{v}}))={\bm{v}}$. For
${\bm{x}},{\bm{y}}\in\mathcal{L}$, and
${\bm{v}}\in\mathcal{T}_{{\bm{x}}}\mathcal{L}$ the exponential map
$\mathrm{exp}_{{\bm{x}}}:\mathcal{T}_{{\bm{x}}}\mathcal{L}\to\mathcal{L}$ and
the logarithmic map
$\mathrm{log}_{{\bm{x}}}:\mathcal{L}\to\mathcal{T}_{{\bm{x}}}\mathcal{L}$ are
given by

$\displaystyle\mathrm{exp}_{{\bm{x}}}({\bm{v}})$
$\displaystyle=\mathrm{cosh}(\|{\bm{v}}\|_{\mathcal{L}}){\bm{x}}+\mathrm{sinh}(\|{\bm{v}}\|_{\mathcal{L}})\frac{{\bm{v}}}{\|{\bm{v}}\|_{\mathcal{L}}},$
(3) $\displaystyle\mathrm{log}_{{\bm{x}}}({\bm{y}})$
$\displaystyle=\frac{\mathrm{arcosh}(-\langle{\bm{x}},{\bm{y}}\rangle_{\mathcal{L}})}{\sqrt{\langle{\bm{x}},{\bm{y}}\rangle_{\mathcal{L}}^{2}-1}}({\bm{y}}+\langle{\bm{x}},{\bm{y}}\rangle_{\mathcal{L}}{\bm{x}}),$
(4)

where
$\|{\bm{v}}\|_{\mathcal{L}}=\sqrt{\langle{\bm{v}},{\bm{v}}\rangle_{\mathcal{L}}}$
is the norm of ${\bm{v}}$.

Isometric isomorphism. The Poincar$\mathrm{\acute{e}}$ ball model
$\mathcal{B}$ and the Klein model $\mathcal{K}$ are two other models of
hyperbolic spaces. The bijections between a node
${\bm{x}}=[{x}_{0},{x}_{1},\cdots,{x}_{n}]\in\mathcal{L}$ and its unique
corresponding node ${\bm{b}}=[b_{0},b_{1},\cdots,b_{n-1}]\in\mathcal{B}$ are
given by

$\displaystyle
p_{\mathcal{L}\to\mathcal{B}}({\bm{x}})=\frac{[{x}_{1},\cdots,{x}_{n}]}{{x}_{0}+1},\
p_{\mathcal{B}\to\mathcal{L}}({\bm{b}})=\frac{[1+\|{\bm{b}}\|^{2},2{\bm{b}}]}{1-\|{\bm{b}}\|^{2}}.$
(5)

The bijections between
${\bm{x}}=[{x}_{0},{x}_{1},\cdots,{x}_{n}]\in\mathcal{L}$ and its unique
corresponding node ${\bm{k}}=[k_{0},k_{1},\cdots,k_{n-1}]\in\mathcal{K}$ are
given by

$\displaystyle
p_{\mathcal{L}\to\mathcal{K}}({\bm{x}})=\frac{[{x}_{1},\cdots,{x}_{n}]}{{x}_{0}},\
p_{\mathcal{K}\to\mathcal{L}}({\bm{k}})=\frac{1}{\sqrt{1-\|{\bm{k}}\|^{2}}}[1,{\bm{k}}].$
(6)

Geometric relationships among the Lorentz model $\mathcal{L}$, the
Poincar$\mathrm{\acute{e}}$ ball model $\mathcal{B}$ and the Klein model
$\mathcal{K}$ are presented in Figure 2.

###  3.2 Graph Convolutional Networks

Figure 2: Geometric relationships among $\mathcal{L}$, $\mathcal{B}$ and
$\mathcal{K}$.

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a graph with a vertex set
$\mathcal{V}$ and an edge set $\mathcal{E}$, and
$\\{{\bm{x}}_{i}^{E}\\}_{i\in\mathcal{V}}$ be node features where E denotes
Euclidean representations. At the $\ell$-th layer of GCNs, the graph
convolution can be formulated into two steps.

Feature Transformation:

${\bm{\bar{h}}}_{i}^{\ell,E}={\bm{W}}^{\ell}{\bm{h}}_{i}^{\ell-1,E}.$ (7)

where ${\bm{h}}_{i}^{\ell-1,E}$ denotes the $i$-th node’s representation at
the ($\ell$-$1$)-th layer and ${\bm{h}}_{i}^{0,E}={\bm{x}}_{i}^{E}$.
${\bm{W}}^{\ell}$ is the learnable transformation matrix at the $\ell$-th
layer. ${\bm{\bar{h}}}_{i}^{\ell,E}$ is the intermediate representation of the
$i$-th node, ready for the next step.

Neighborhood Aggregation:

$\left\\{\begin{aligned}
{\bm{m}}_{i}^{\ell,E}&=({\bm{\bar{h}}}_{i}^{\ell,E}+\sum_{j\in\mathcal{N}(i)}w_{ij}{\bm{\bar{h}}}_{j}^{\ell,E})\\\
{\bm{h}}_{i}^{\ell,E}&=\sigma({\bm{m}}_{i}^{\ell,E})\\\ \end{aligned}\right.,$
(8)

where $\mathcal{N}(i)$ denotes the set of neighbor nodes of the $i$-th node,
and $w_{ij}$ is the aggregation weight. ${\bm{m}}_{i}^{\ell,E}$ is the
aggregated message, that is sent to a non-linear activation function
$\sigma(\cdot)$ to output the node representation ${\bm{h}}_{i}^{\ell,E}$ at
the $\ell$-th layer.

By stacking multiple graph convolutional layers, the feature transformation
enables GCNs to learn desirable node embeddings for a target task, _e.g_.,
more discriminative for classifications. The neighborhood aggregation enables
GCNs to exploit graph topology structures.

##  4 Hyperbolic-to-Hyperbolic GCN

We present H2H-GCN that directly performs graph convolution on hyperbolic
manifolds to keep global hyperbolic structure. First, we explain how to
generate hyperbolic node representations as input node features are usually
Euclidean. Then, we elaborate the developed hyperbolic feature transformation
and hyperbolic neighborhood aggregation. Next, we construct the H2H-GCN
architecture used for link prediction, node classification and graph
classification. Finally, we describe how to optimize parameters in the
H2H-GCN.

###  4.1 Hyperbolic Node Representations

Let $\\{{\bm{x}}_{i}^{E}\\}_{i\in\mathcal{V}}$ be input Euclidean node
features, and ${\bm{o}}\coloneqq[1,0,\cdots,0]$ denote the origin on the
manifold $\mathcal{L}$ of the Lorentz model. There is
$\langle{\bm{o}},[0,{\bm{x}}_{i}^{E}]\rangle_{\mathcal{L}}=0$, where
$\langle\cdot,\cdot\rangle_{\mathcal{L}}$ denotes the Lorentz inner product
defined in Eq. (1). We can reasonably regard $[0,{\bm{x}}_{i}^{E}]$ as a node
on the tangent space at the origin ${\bm{o}}$. H2H-GCN uses the exponential
map defined in Eq. (3) to generate hyperbolic node representations on the
Lorentz model:

$\displaystyle{\bm{x}}_{i}^{\mathcal{L}}$
$\displaystyle=\mathrm{exp}_{{\bm{o}}}\big{(}[0,{\bm{x}}_{i}^{E}]\big{)}$ (9)
$\displaystyle=\Big{[}\mathrm{cosh}\big{(}\|{\bm{x}}_{i}^{E}\|_{2}\big{)},\mathrm{sinh}\big{(}\|{\bm{x}}_{i}^{E}\|_{2}\big{)}\frac{{\bm{x}}_{i}^{E}}{\|{\bm{x}}_{i}^{E}\|_{2}}\Big{]}.$

###  4.2 Hyperbolic Feature Transformation

The feature transformation in (Euclidean) GCNs defined in Eq. (7) is a linear
transformation realized via a matrix-vector multiplication. Nevertheless, it
will break the hyperbolic manifold constraint while applying matrix-vector
multiplication to hyperbolic node representations, making the transformed
nodes not lie on hyperbolic manifolds. We developed a Lorentz linear
transformation to tackle this problem.

###### Definition 1

(The Lorentz linear transformation). For any ${\bm{x}}\in\mathcal{L}$, the
Lorentz linear transformation is defined as

$\displaystyle{\bm{y}}$ $\displaystyle={\bm{W}}{\bm{x}}$ (10)
$\displaystyle\mathrm{s.t.}~{}{\bm{W}}=$
$\displaystyle\begin{bmatrix}1&{\bm{0}}^{\top}\\\
{\bm{0}}&\widehat{{\bm{W}}}\end{bmatrix},\widehat{{\bm{W}}}^{\top}\widehat{{\bm{W}}}={\bm{I}},$

where ${\bm{W}}$ is a transformation matrix, and $\widehat{{\bm{W}}}$ is
called a transformation sub-matrix. ${\bm{0}}$ is a column vector of zeros,
and ${\bm{I}}$ is an identity matrix.

###### Proposition 1

The Lorentz linear transformation defined in Definition. 1 is manifold-
preserving. It ensures that the output ${\bm{y}}$ still lies on the manifold
$\mathcal{L}$ of the Lorentz model.

######  Proof 4.1.

For any ${\bm{x}}=[{x}_{0},{x}_{1},\cdots,{x}_{n}]\in\mathcal{L}$, we have

$\displaystyle-{x}_{0}^{2}+{\bm{x}}_{1:n}^{\top}{\bm{x}}_{1:n}=-1,~{}\textrm{and
}{x}_{0}>0,$

where ${\bm{x}}_{1:n}=[{x}_{1},{x}_{2},\cdots,{x}_{n}]$. After applying the
Lorentz linear transformation to ${\bm{x}}$, we have

${\bm{y}}={\bm{W}}{\bm{x}}=\Big{[}{x}_{0},\widehat{{\bm{W}}}{\bm{x}}_{1:n}\Big{]},$

satisfying

$y_{0}=x_{0}>0,\\\ $

and

$\displaystyle\langle{\bm{y}},{\bm{y}}\rangle_{\mathcal{L}}=-{x}_{0}^{2}+(\widehat{{\bm{W}}}{\bm{x}}_{1:n})^{\top}\widehat{{\bm{W}}}{\bm{x}}_{1:n}$
$\displaystyle=-{x}_{0}^{2}+{\bm{x}}_{1:n}^{\top}(\widehat{{\bm{W}}}^{\top}\widehat{{\bm{W}}}){\bm{x}}_{1:n}$
$\displaystyle=-{x}_{0}^{2}+{\bm{x}}_{1:n}^{\top}{\bm{x}}_{1:n}$
$\displaystyle=-1.$

Thus, ${\bm{y}}$ lies on the manifold $\mathcal{L}$ of the Lorentz model.

We utilize the Lorentz linear transformation as the hyperbolic feature
transformation in H2H-GCN. At the $\ell$-th layer, we take the node
representation from the previous layer ${\bm{h}}_{i}^{\ell-1,\mathcal{L}}$ and
the transformation matrix ${\bm{W}}^{\ell}$ as input. The $i$-th node’s
intermediate representation is calculated by

$\displaystyle{\bm{\bar{h}}}_{i}^{\ell,\mathcal{L}}$
$\displaystyle={\bm{W}}^{\ell}{\bm{h}}_{i}^{\ell-1,\mathcal{L}}$ (11)
$\displaystyle\mathrm{s.t.}~{}{\bm{W}}^{\ell}$
$\displaystyle=\begin{bmatrix}1&{\bm{0}}^{\top}\\\
{\bm{0}}&\widehat{{\bm{W}}}^{\ell}\end{bmatrix},\widehat{{\bm{W}}}^{\ell^{\top}}\widehat{{\bm{W}}}={\bm{I}},$

where ${\bm{h}}_{i}^{0,\mathcal{L}}={\bm{x}}_{i}^{\mathcal{L}}$. The
intermediate node representation ${\bm{\bar{h}}}_{i}^{\ell,\mathcal{L}}$ is
ready for hyperbolic neighborhood aggregation in Section 4.3. We describe an
effective way to learn transformation matrix ${\bm{W}}^{\ell}$, a constrained
parameter, via optimization on a Stiefel manifold in Section 4.5.

Figure 3: Einstein midpoint on the Klein model $\mathcal{K}$, taking three
nodes ${\bm{u}},{\bm{p}},{\bm{q}}\in\mathcal{K}$ for example.

###  4.3 Hyperbolic Neighborhood Aggregation

The neighborhood aggregation in GCNs defined in Eq.(8) updates a node
representation by aggregating the message from its neighbor node set, enabling
GCNs to capture graph topological structure. A generalization of Euclidean
mean aggregation in hyperbolic spaces is Fréchet mean [frechet1948elements].
However, Fréchet mean is difficult to apply because it does not have a closed-
form solution. We adopt the Einstein midpoint [38] as the hyperbolic
neighborhood aggregation in H2H-GCN. In this case, our hyperbolic neighborhood
aggregation possesses two desirable properties: translation invariance and
rotation invariance. The aggregated hyperbolic average is invariant to
translating the input node set by a same distance in a common direction, and
invariant to rotating the input node set by a same angle around the origin.

The Einstein midpoint takes the form in the Klein model, illustrated in Figure
3. We first project the intermediate node representations from the Lorentz
model to the Klein model, and then calculate the hyperbolic average via the
Einstein midpoint. The aggregated hyperbolic average on the Klein model is
projected back to the Lorentz model. Formally, given the intermediate
representation of a node ${\bm{\bar{h}}}_{i}^{\ell,\mathcal{L}}$ and the
intermediate representations of its neighbor nodes
$\\{{\bm{\bar{h}}}_{j}^{\ell,\mathcal{L}}\\}_{j\in\mathcal{N}(i)}$, the
hyperbolic neighborhood aggregation on the Lorentz model is given by

$\left\\{\begin{aligned}
{\bm{\bar{h}}}_{j}^{\ell,\mathcal{K}}&=p_{\mathcal{L}\to\mathcal{K}}({\bm{\bar{h}}}_{j}^{\ell,\mathcal{L}})\\\
{\bm{m}}_{i}^{\ell,\mathcal{K}}&=\sum_{j\in\widehat{\mathcal{N}}(i)}\gamma_{j}{\bm{\bar{h}}}_{j}^{\ell,\mathcal{K}}/\sum_{j\in\widehat{\mathcal{N}}(i)}\gamma_{j}\\\
{\bm{m}}_{i}^{\ell,\mathcal{L}}&=p_{\mathcal{K}\to\mathcal{L}}({\bm{m}}_{i}^{\ell,\mathcal{K}})\\\
\end{aligned}\right.,$ (12)

where ${\bm{\bar{h}}}_{j}^{\ell,\mathcal{K}}$ is the $j$-th node’s
intermediate representation on the Klein model.
$p_{\mathcal{L}\to\mathcal{K}}(\cdot)$ and
$p_{\mathcal{K}\to\mathcal{L}}(\cdot)$ are the isometric and isomorphic
bijections between the Lorentz model and the Klein model as defined in Eq.
(6). $\widehat{\mathcal{N}}(i)$ is a node set consisting of the $i$-th node
and its neighbor nodes.
$\gamma_{j}=\frac{1}{\sqrt{1-\|\bar{{\bm{h}}}_{j}^{\mathcal{K}}\|^{2}}}$
denotes the Lorentz factor. ${\bm{m}}_{i}^{\ell,\mathcal{K}}$ is the
hyperbolic average on the Klein model that aggregates the message from
$\widehat{\mathcal{N}}(i)$ via the Einstein midpoint. We get the hyperbolic
average on the Lorentz model ${\bm{m}}_{i}^{\ell,\mathcal{L}}$ by projecting
${\bm{m}}_{i}^{\ell,\mathcal{K}}$ to $\mathcal{L}$.

The non-linear activation plays an important role in GCNs, which prevents a
multi-layer network from collapsing into a single layer network. However,
applying commonly-used non-linear activation functions (_e.g_., ReLU) on the
Lorentz representation will break the manifold constraint of the Lorentz
model. We notice that the non-linear activation on the
Poincar$\mathrm{\acute{e}}$ ball model $\mathcal{B}$ is manifold-preserving:
for any ${\bm{b}}\in\mathcal{B}$, we have $\sigma({\bm{b}})\in\mathcal{B}$.
Inspired by this, we project hyperbolic average
${\bm{m}}_{i}^{\ell,\mathcal{L}}$ to the Poincar$\mathrm{\acute{e}}$ ball
model to apply non-linear activation, and then project the result back to the
Lorentz model, given by

${\bm{h}}_{i}^{\ell,\mathcal{L}}=p_{\mathcal{B}\to\mathcal{L}}\Big{(}\sigma\big{(}p_{\mathcal{L}\to\mathcal{B}}({\bm{m}}_{i}^{\ell,\mathcal{L}})\big{)}\Big{)},$
(13)

where $p_{\mathcal{B}\to\mathcal{L}}(\cdot)$ and
$p_{\mathcal{L}\to\mathcal{B}}(\cdot)$ are the isometric and isomorphic
bijections between the Lorentz model and the Poincar$\mathrm{\acute{e}}$ ball
model as defined in Eq. (5). After Eq.(13), H2H-GCN obtains the output of the
$\ell$-th layer: the $i$-th node’s representation
${\bm{h}}_{i}^{\ell,\mathcal{L}}$ on the Lorentz model.

| Datasets | Disease | Airport | PubMed | Cora
---|---|---|---|---|---
| Methods | LP | NC | LP | NC | LP | NC | LP | NC
Shallow |  Euc [10] | 59.8 $\pm$ 2.0 | 32.5 $\pm$ 1.1 | 92.0 $\pm$ 0.0 | 60.9 $\pm$ 3.4 | 83.3 $\pm$ 0.1 | 48.2 $\pm$ 0.7 | 82.5 $\pm$ 0.3 | 23.8 $\pm$ 0.7
Hyp [28] | 63.5 $\pm$ 0.6 | 45.5 $\pm$ 3.3 | 94.5 $\pm$ 0.0 | 70.2 $\pm$ 0.1 | 87.5 $\pm$ 0.1 | 68.5 $\pm$ 0.3 | 87.6 $\pm$ 0.2 | 22.0 $\pm$ 1.5
Euc-Mixed [10] | 49.6 $\pm$ 1.1 | 35.2 $\pm$ 3.4 | 91.5 $\pm$ 0.1 | 68.3 $\pm$ 2.3 | 86.0 $\pm$ 1.3 | 63.0 $\pm$ 0.3 | 84.4 $\pm$ 0.2 | 46.1 $\pm$ 0.4
Hyp-Mixed [10] | 55.1 $\pm$ 1.3 | 56.9 $\pm$ 1.5 | 93.3 $\pm$ 0.0 | 69.6 $\pm$ 0.1 | 83.8 $\pm$ 0.3 | 73.9 $\pm$ 0.2 | 85.6 $\pm$ 0.5 | 45.9 $\pm$ 0.3
NNs |  MLP [10] | 72.6 $\pm$ 0.6 | 28.8 $\pm$ 2.5 | 89.8 $\pm$ 0.5 | 68.6 $\pm$ 0.6 | 84.1 $\pm$ 0.9 | 72.4 $\pm$ 0.2 | 83.1 $\pm$ 0.5 | 51.5 $\pm$ 1.0
HNN [17] | 75.1 $\pm$ 0.3 | 41.0 $\pm$ 1.8 | 90.8 $\pm$ 0.2 | 80.5 $\pm$ 0.5 | 94.9 $\pm$ 0.1 | 69.8 $\pm$ 0.4 | 89.0 $\pm$ 0.1 | 54.6 $\pm$ 0.4
GCNs |  GCN [23] | 64.7 $\pm$ 0.5 | 69.7 $\pm$ 0.4 | 89.3 $\pm$ 0.4 | 81.4 $\pm$ 0.6 | 91.1 $\pm$ 0.5 | 78.1 $\pm$ 0.2 | 90.4 $\pm$ 0.2 | 81.3 $\pm$ 0.3
GAT [39] | 69.8 $\pm$ 0.3 | 70.4 $\pm$ 0.4 | 90.5 $\pm$ 0.3 | 81.5 $\pm$ 0.3 | 91.2 $\pm$ 0.1 | 79.0 $\pm$ 0.3 | 93.7 $\pm$ 0.1 | 83.0 $\pm$ 0.7
GraphSAGE [19] | 65.9 $\pm$ 0.3 | 69.1 $\pm$ 0.6 | 90.4 $\pm$ 0.5 | 82.1 $\pm$ 0.5 | 86.2 $\pm$ 1.0 | 77.4 $\pm$ 2.2 | 85.5 $\pm$ 0.6 | 77.9 $\pm$ 2.4
Hyp GCNs |  SGC [43] | 65.1 $\pm$ 0.2 | 69.5 $\pm$ 0.2 | 89.8 $\pm$ 0.3 | 80.6 $\pm$ 0.1 | 94.1 $\pm$ 0.0 | 78.9 $\pm$ 0.0 | 91.5 $\pm$ 0.1 | 81.0 $\pm$ $0.1$
HGCN [10] | 90.8 $\pm$ 0.3 | 74.5 $\pm$ 0.9 | 96.4 $\pm$ 0.1 | 90.6 $\pm$ 0.2 | 96.3 $\pm$ 0.0 | 80.3 $\pm$ 0.3 | 92.9 $\pm$ 0.1 | 79.9 $\pm$ 0.2
H2H-GCN (Ours) |  97.8 $\pm$ 0.0 |  90.3 $\pm$ 0.8 |  96.7 $\pm$ 0.0 |  91.0 $\pm$ 0.3 |  97.1 $\pm$ 0.0 |  82.3 $\pm$ 0.4 |  95.4 $\pm$ 0.0 |  83.6 $\pm$ 0.8

Table 1: ROC AUC for Link Prediction (LP), and F1 score (Disease, binary
class) and accuracy (the others, multi-class) for Node Classification (NC)
tasks.

###  4.4 H2H-GCN Architecture

Input: Graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$; node features
$\\{{\bm{x}}_{i}^{E}\\}_{i\in\mathcal{V}}$; number of layers $L$;
transformation matrices $\\{{\bm{W}}^{\ell}\\}_{\ell=1}^{L}$; non-linearity
activation function $\sigma(\cdot)$.

Output: H2H-GCN node embeddings
$\\{{\bm{h}}_{i}^{L,\mathcal{L}}\\}_{i\in\mathcal{V}}$.

1 Orthogonally initialize transformation sub-matrices
$\\{\widehat{{\bm{W}}}^{\ell}\\}_{\ell=1}^{L}$;

2 Construct ${\bm{W}}^{\ell}=\begin{bmatrix}1&{\bm{0}}^{\top}\\\
{\bm{0}}&\widehat{{\bm{W}}}^{\ell}\end{bmatrix},\forall\ell\in\\{1,\cdots,L\\}$;

3 Generate hyperbolic node representations
$\\{{\bm{x}}_{i}^{\mathcal{L}}\\}_{i\in\mathcal{V}}$ via Eq. (9);

4 ${\bm{h}}_{i}^{0,\mathcal{L}}={\bm{x}}_{i}^{\mathcal{L}},\forall
i\in\mathcal{V}$, ;

5 for  _$\ell=1$ to $L$_ do

6 Generate intermediate node representations
$\\{{\bm{\bar{h}}}_{i}^{\ell,\mathcal{L}}\\}_{i\in\mathcal{V}}$ via the
hyperbolic feature transformation in Eq. (11);

7 for  _$i\in\mathcal{V}$_ do

8 Generate hyperbolic average ${\bm{m}}_{i}^{\ell,\mathcal{L}}$ by aggregating
message from
$\\{{\bm{\bar{h}}}_{j}^{\ell,\mathcal{L}}\\}_{j\in\hat{\mathcal{N}}(i)}$ via
the hyperbolic heighborhood aggregation in Eq. (12);

9 Generate the node representation ${\bm{{h}}}_{i}^{\ell,\mathcal{L}}$ via the
non-linear activation on ${\bm{m}}_{i}^{\ell,\mathcal{L}}$ in Eq. (13);

10

11 end for

12

13 end for

return H2H-GCN node embeddings
$\\{{\bm{h}}_{i}^{L,\mathcal{L}}\\}_{i\in\mathcal{V}}$.

Algorithm 1 H2H-GCN embedding generation.

We summarize the H2H-GCN embedding generation algorithm as shown in Algorithm
1. Given a graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ with a vertex set
$\mathcal{V}$ and an edge set $\mathcal{E}$, H2H-GCN first maps input
Euclidean node features $\\{{\bm{x}}_{i}^{E}\\}_{i\in\mathcal{V}}$ to
hyperbolic space via Eq. (9). The obtained hyperbolic node representations
$\\{{\bm{h}}_{i}^{0,\mathcal{L}}\\}_{i\in\mathcal{V}}$ are sent to a multi-
layer H2H-GCN. At the $\ell$-th layer, the input node representation
${\bm{h}}_{i}^{\ell-1,\mathcal{L}}$ from previous layer is passed through the
hyperbolic feature transformation in Eq. (11), and is updated via the
hyperbolic neighborhood aggregation in Eq. (12) and the non-linear activation
in Eq.(13). After $L$ layers, we obtain the H2H-GCN node embeddings
$\\{{\bm{h}}_{i}^{L,\mathcal{L}}\\}_{i\in\mathcal{V}}$.

For link prediction, we use the Fermi-Dirac decoder [24, 28] to calculate
probability scores for the edge between the $i$-th node and the $j$-th node,
given by

${p}\big{(}(i,j)\in\mathcal{E}|{\bm{h}}_{i}^{L,\mathcal{L}},{\bm{h}}_{j}^{L,\mathcal{L}}\big{)}=\big{[}e^{(d_{\mathcal{L}}{({\bm{h}}_{i}^{L,\mathcal{L}},{\bm{h}}_{j}^{L,\mathcal{L}})^{2}}-r)/t}+1\big{]}^{-1},$
(14)

where $d_{\mathcal{L}}(\cdot,\cdot)$ is the hyperbolic distance function
defined in Eq. (2), and $r$ and $t$ are hyper-parameters. Following [10], we
use the negative sampling strategy and the cross entropy loss for training.

For node classification and graph classification, we exploit a centroid-based
classification method studied in [25]. Specifically, we introduce a set of
centroids
$C=\\{{\bm{c}}_{1}^{\mathcal{L}},{\bm{c}}_{2}^{\mathcal{L}},\cdots,{\bm{c}}_{|C|}^{\mathcal{L}}\\}$
lying on the Lorentz model, then calculate a distance matrix
${\bm{D}}\in\mathbb{R}^{|\mathcal{V}|\times|C|}$ whose element
${\bm{D}}_{i,j}=d_{\mathcal{L}}({\bm{h}}_{i}^{L,\mathcal{L}},{\bm{c}}_{j}^{\mathcal{L}})$
represents the distance between the $i$-th node embedding
${\bm{h}}_{i}^{L,\mathcal{L}}$ and the $j$-th centroid
${\bm{c}}_{j}^{\mathcal{L}}$. For node classification, we send ${\bm{D}}_{i}$,
the $i$-th row of distance matrix that contains the distance information
between the $i$-th node and all centroids, to a classifier to predict the
category of the $i$-th node. For graph classification, we apply average
pooling to $\\{{\bm{D}}_{i}\\}_{i=1}^{|\mathcal{V}|}$ as a readout operation
to yield a graph embedding
$\frac{1}{|\mathcal{V}|}\sum_{i=1}^{|\mathcal{V}|}{\bm{D}}_{i}$, followed a
classifier for prediction. On both classification tasks, we use softmax
classifiers and cross entropy loss functions.

###  4.5 Optimization

We explain how to learn parameters in H2H-GCN, especially the transformation
matrix ${\bm{W}}$ (omitting the layer number ℓ for convenience) in the
hyperbolic feature transformation Eq. (11), that is an optimization problem
with the orthogonal constraint. Other parameters can be learned by a standard
gradient descent optimizer straightforwardly.

The transformation matrix ${\bm{W}}=\begin{bmatrix}1&{\bm{0}}^{\top}\\\
{\bm{0}}&\widehat{{\bm{W}}}\end{bmatrix}$ is a block diagonal matrix that
consists of a scalar $1$ and an orthogonal matrix
$\widehat{{\bm{W}}}\in\mathrm{St}(n^{\prime},n)$ which resides on the Stiefel
manifold.

######  Definition 4.2.

(The Stiefel manifold). The set of ($n^{\prime}\times n$)-dimensional
matrices, $n\leq n^{\prime}$, with orthonormal columns forms a compact
Riemannian manifold called the Stiefel manifold $\mathrm{St(n^{\prime},n)}$
[5].

$\mathrm{St(n^{\prime},n)}\triangleq\\{{\bm{M}}\in\mathbb{R}^{n^{\prime}\times
n}:{\bm{M}}^{\top}{\bm{M}}={\bm{I}}\\}.$ (15)

While updating ${\bm{W}}$, we keep $1$ unchanged and introduce a Riemannian
stochastic gradient descent optimizer to update $\widehat{{\bm{W}}}$.
Formally, let $J$ be the loss function, _e.g_., cross entropy loss for
classifications. $\widehat{{\bm{W}}}$ is updated by

$\vspace{-5pt}\left\\{\begin{aligned}
&{\bm{P}}^{(t)}=\eta\pi_{\widehat{{\bm{W}}}^{(t)}}(\nabla_{\widehat{{\bm{W}}}}^{(t)})\\\
&\widehat{{\bm{W}}}^{(t+1)}=r_{\widehat{{\bm{W}}}^{(t)}}(-{\bm{P}}^{(t)})\\\
\end{aligned}\right.,\vspace{0pt}$ (16)

where $\eta$ is the learning rate.
$\nabla_{\widehat{{\bm{W}}}}^{(t)}=\mathrm{d}J/\widehat{{\bm{W}}}^{(t)}$
denotes the Eulcidean gradient of the loss function $J$ with respect to
$\widehat{{\bm{W}}}^{(t)}$ calculated at time $t$.
$\pi_{{\widehat{{\bm{W}}}^{(t)}}}(\cdot)$ is an orthogonal projection that
transforms the Euclidean gradient to the Riemannian gradient

$\pi_{{\widehat{{\bm{W}}}^{(t)}}}\big{(}{\nabla_{\widehat{{\bm{W}}}}^{(t)}}\big{)}={\nabla_{\widehat{{\bm{W}}}}^{(t)}}-\frac{1}{2}{\widehat{{\bm{W}}}^{(t)}}\Big{(}{\widehat{{\bm{W}}}^{(t)\top}}{\nabla_{\widehat{{\bm{W}}}}^{(t)}}+{\nabla_{\widehat{{\bm{W}}}}^{(t)}}^{\top}{\widehat{{\bm{W}}}^{(t)}}\Big{)}.$
(17)

$r_{{\widehat{{\bm{W}}}^{(t)}}}(\cdot)$ is a retraction operation, defined as

$\displaystyle
r_{{\widehat{{\bm{W}}}^{(t)}}}(-{\bm{P}}^{(t)})=\mathrm{{qf}}({\widehat{{\bm{W}}}^{(t)}}-{\bm{P}}^{(t)}),$
(18)

where $\mathrm{qf}(\cdot)$ extracts the orthogonal factor in the QR
decomposition. The retraction operation prevents the updated
$\widehat{{\bm{W}}}^{(t+1)}$ from falling off the Stiefel manifold.

##  5 Experiments

We evaluate the proposed H2H-GCN on the link prediction, node classification
and graph classification tasks, and comprehensively compare H2H-GCN with a
variety of state-of-the-art Euclidean GCNs and hyperbolic GCNs.

| Dimensionality
---|---
Methods | 3 | 5 | 10 | 20 | 256
Euclidean [25] | 77.2 $\pm$ 0.12 | 90.0 $\pm$ 0.21 | 90.6 $\pm$ 0.17 | 94.8 $\pm$ 0.25 | 95.3 $\pm$ 0.17
HGNN [25] | 94.1 $\pm$ 0.03 | 95.6 $\pm$ 0.14 | 96.4 $\pm$ 0.23 | 96.6 $\pm$ 0.22 | 95.3 $\pm$ 0.28
H2H-GCN (Ours) |  95.4 $\pm$ 0.26 |  96.7 $\pm$ 0.12 |  96.8 $\pm$ 0.04 |  97.0 $\pm$ 0.05 |  97.2 $\pm$ 0.03
Table 2: Results on synthetic graph classification where F1 (macro) score and
standard deviation are reported.

###  5.1 Link Prediction and Node Classification

The link prediction (LP) task is to predict the existence of links among nodes
in a graph, and the node classification (NC) task is to predict labels of
nodes in a graph. They have many applications such as predicting friendships
among users in a social network and predicting research directions of papers
in a citation network. For link prediction, we report area under the ROC curve
(AUC), and for node classification, we report F1 score for binary-class
datasets and accuracy for multi-class datasets. Following HGCN [10], a link
prediction regularization objective is added in the node classification task.

##### Datasets.

Disease [10] is constructed by simulating the SIR disease spreading model [1],
where the label of a node indicates whether the node was infected or not, and
the feature of a node indicates the susceptibility of the node to the disease.
Cora and Pubmed [33] are citation network datasets where nodes are scientific
papers and edges represent citation links. Cora contains $7$ classes of
machine learning papers, and there are $2,708$ nodes, $5,429$ edges and
$1,433$ features per node. Pubmed contains $3$ classes of medicine
publications, and there are $19,717$ nodes, $44,338$ edges and $500$ features
per node. Airport [10] is a flight network dataset where nodes are airports
and edges represent the airline routes. There are $2,236$ nodes in total and
the label of a node is the population of the country where the airport (node)
belongs to.

##### Baselines.

We consider four types of baseline methods: shallow methods, neural networks
(NNs) methods, GCNs, and hyperbolic GCNs (Hyp GCNs). Shallow methods optimize
to minimize a reconstruction loss, and the parameters in models act as an
embedding look-up table. We consider Euclidean embeddings (Euc) [10] and its
hyperbolic extension (Hyp) [28]. As the two embeddings fail to leverage node
features, Euc-Mixed [10] and Hyp-Mixed [10] concatenate the shallow embeddings
with node features for a fair comparison with other methods using node
features. NNs methods only utilize the node features but does not consider
graph structures. Compared NNs methods include Euclidean multi-layer
percerptron (MLP) and its hyperbolic extension: hyperbolic neural networks
(HNN) [17]. For GCNs, we compare H2H-GCN with several Euclidean state-of-the-
art GCNs models: GCN [23], GraphSAGE [19], graph attention networks (GAT)
[39], and simplified graph convolution (SGC) [43]. For Hyp GCNs, we consider
HGCN [10] that performs graph convolutional operations in tangent spaces.

The comparisons are presented in Table 1. We notice that HNN, a generalization
of MLP to hyperbolic spaces, outperforms MLP on most tasks. It indicates that
hyperbolic spaces are more suitable for modeling graphs compared than
Euclidean spaces. A similar conclusion can be drawn while comparing Euclidean
GNNs with hyperbolic GCNs. HGCN works better than state-of-the-art Euclidean
GCNs on most datasets.

The proposed H2H-GCN consistently achieve the best performance on both tasks
on all datasets. We take the Disease dataset as an example to analyze the
effectiveness of H2H-GCN. The Disease dataset is a tree network that possesses
a strong hierarchical structure. As hyperbolic spaces can be viewed as smooth
versions of trees, hyperbolic GCNs should work better than Euclidean GCNs. The
results are in line with expectation that both HGCN and H2H-GCN show
significant improvement on the LP and NC tasks compared with Euclidean
methods. In particular, H2H-GCN achieves an average of $22.7\%$ (LP) and
$19.9\%$ (NC) performance gains than state-of-the-art Euclidean GCNs, and
$7.0\%$ (LP) and $15.8\%$ (NC) performance gains than HGCN. It demonstrates
that our H2H-GCN is superior to hyperbolic GCNs which rely on tangent spaces.
We owe it to our developed hyperbolic graph convolution that directly works on
the hyperbolic manifold. Both the hyperbolic feature transformation and the
hyperbolic neighborhood aggregation are manifold-preserving. It can avoid the
distortion caused by tangent space approximations and keep the global
hyperbolic geometry underlying graphs.

###  5.2 Graph Classification

Graph classification is important to many real-world applications associated
with graph-structured data, _e.g_., chemical drug analysis. We first evaluate
our method on a synthetic graph dataset, and then on several commonly-used
molecular graph datasets.

####  5.2.1 Synthetic Graphs

Following [25], we take three graph generation algorithms: Erdős-Rényi [15],
Barabási-Albert [3] and Watts-Strogatz [41] to construct a synthetic graph
dataset. For each graph generation algorithm, we generate $6,000$ graphs and
divide them into 3 equal parts for the training, validation, and test (see
[25] for more generation details). Typical properties, such as small-world
property of graphs generated by Watts-Strogatz algorithm and scale-free
property of graphs generated by Barabási-Albert algorithm, can be explained by
an underlying hyperbolic geometry [24], thus it is more suitable for modeling
such graphs in hyperbolic spaces than in Euclidean spaces.

We compare the proposed H2H-GCN with Euclidean embeddings [25] and HGNN [25].
The F1 scores of different embedding dimensionalities are presented in Table
2. The performance of HGNN and our method surpasses Euclidean embeddings by a
large margin when embedding dimensionality is low. It is because hyperbolic
spaces can well capture the hyperbolic geometry underlying these synthetic
graphs. As the embedding dimensionality increases (_e.g_., $256$), HGNN tends
to be comparable with Euclidean alternative $95.3\%$, while our method
achieves $97.2\%$, $1.9\%$ higher than HGNN. The reason is that the distortion
caused by tangent space approximations in HGNN becomes signficant with the
increase of dimensionality of embedding spaces, leading to an inferior
performance. H2H-GCN tackles this problem by directly learning node embeddings
in the hyperbolic space. It shows the best performance from $3$-dimensional
embeddings to $256$-dimensional embeddings.

####  5.2.2 Molecular Graphs

We evaluate our method on several chemical datasets to predict the function of
molecular graphs. D&D [14] consists of X-ray crystal structures. It has
$1,178$ graphs in total, and $2$ classes indicting the molecular is an enzyme
or not. Proteins [6] has $1,113$ graphs, and $3$ classes of graphs
representing helix, sheet or turn. Enzymes [32] contains $600$ graphs and $6$
classes in total.

We compare our method with several state-of-the-art Euclidean GCNs, including
DGCNN [46], DiffPool [45], ECC [34], GIN [44] and GraphSAGE [19], and a
hyperbolic GCN, _i.e_., HGNN [25] that performs Euclidean graph convolutional
operations in tangent spaces. In general, researchers adopt tenfold cross
validation for model evaluation. However, as pointed out by [16], the data
splits are different and the experimental procedures are often ambiguous of
different works, which results in unfair comparisons. To this end, the work in
[16] provides a uniform and rigorous benchmarking of state-of-the-art models.
In this part, we follow the same experimental procedures and use the same data
splits as [16] for fair comparisons.

Methods | D&D | Proteins | Enzymes
---|---|---|---
DGCNN [46] | 76.6 $\pm$ 4.3 | 72.9 $\pm$ 3.5 | 38.9 $\pm$ 5.7
DiffPool [45] | 75.0 $\pm$ 3.5 | 73.7 $\pm$ 3.5 | 59.5 $\pm$ 5.6
ECC [34] | 72.6 $\pm$ 4.1 | 72.3 $\pm$ 3.4 | 29.5 $\pm$ 8.2
GIN [44] | 75.3 $\pm$ 2.9 | 73.3 $\pm$ 4.0 | 59.6 $\pm$ 4.5
GraphSAGE [19] | 72.9 $\pm$ 2.0 | 73.0 $\pm$ 4.5 | 58.2 $\pm$ 6.0
HGNN [25] | 75.8 $\pm$ 3.3 | 73.7 $\pm$ 2.3 | $51.3\pm 6.1$
H2H-GCN (Ours) |  78.2 $\pm$ 3.3 |  74.4 $\pm$ 3.0 |  61.3 $\pm$ 4.9

Table 3: Results on chemical graph classification where mean accuracy and
standard deviation are reported.

The mean accuracy and standard deviation are reported in Table 3. We observe
that HGNN are comparable with Euclidean methods, which illustrates two
possible reasons: either hyperbolic GCNs are not suitable for the three
datasets, or some factors in HGNN limit representation ability of hyperbolic
GCNs. The performance of H2H-GCN may give the answer. It achieves the best
performance on all datasets: $1.6\%$ higher than DGCNN on D&D, $0.7\%$ higher
than DiffPool on Proteins, and $1.7\%$ higher than GIN on Enzymes. It
indicates that there are some limitations in HGNN and H2H-GCN overcomes it.
Compared with HGNN that does graph convolutional operations in the tangent
space, the key difference lies in the propoed H2H-GCN performs a hyperbolic
graph convolution in the hyperbolic space. In this way, H2H-GCN preservs the
global hyperbolic geometry, leading to a superior performance.

| Dimensionality
---|---
Methods | 3 | 5 | 10 | 20 | 128
Euclidean [25] |  64.2 $\pm$ 4.9 |  71.2 $\pm$ 3.4 |  76.2 $\pm$ 1.5 |  78.1 $\pm$ 2.1 |  80.4 $\pm$ 0.9
HGNN [25] |  65.3 $\pm$ 3.6 |  71.0 $\pm$ 3.4 |  76.1 $\pm$ 1.5 |  79.2 $\pm$ 1.6 |  80.1 $\pm$ 0.9
HGCN [10] |  70.8 $\pm$ 1.6 |  75.4 $\pm$ 1.7 |  78.1 $\pm$ 0.8 |  79.7 $\pm$ 1.4 |  81.7 $\pm$ 0.7
H2H-GCN (Ours) |  73.1 $\pm$ 2.5 |  77.8 $\pm$ 0.6 |  79.9 $\pm$ 0.9 |  81.2 $\pm$ 0.9 |  83.6 $\pm$ 0.8
Table 4: Comparisons of embedding dimensionality for node classification on
Cora where accuracy and standard deviation are reported.

###  5.3 Dimensionality Comparisons

We test the effect of embedding dimensionality from $3$ to $128$ for node
classification on Cora, and report the performance of Euclidean embedding
[25], HGNN [25], HGCN [10] and the proposed H2H-GCN in Table 4. HGNN gets
comparable results with Euclidean embeddings, while HGCN shows pretty
improvments. It outperforms HGCN, $2.3\%$ and $1.9\%$ higher than HGCN when
embedding dimensionalities are $3$ and $256$, respectively. We claim that the
distortion caused by tangent space approximations exists in both low and high
embedding dimensionalities. Although increasing the dimensionality can improve
performance, it cannot solve this problem. H2H-GCN tackles it by proposing a
hyperbolic graph convolution to directly work on the hyperbolic manifold. Such
a manifold-to-manifold method achieves remarkable improvments. H2H-GCN works
the best in all embedding dimensionalities.

##  6 Conclusion

In this paper, we have presented a hyperbolic-to-hyperbolic graph
convolutional network (H2H-GCN) for embedding graph with hierarchical
structure into hyperbolic spaces. The developed hyperbolic graph convolution
which consists of a hyperbolic feature transformation and a hyperbolic
neighborhood aggregation, can be directly conducted on hyperbolic manifolds.
The both operations can ensure that the output still lies on the hyperbolic
manifold. In constrast to existing hyperbolic GCNs relying on tangent spaces,
H2H-GCN can avoid the distortion caused by tangent space approximations and
keep the global hyperbolic geometry underlying graphs. Extensive experiments
on link prediction, node classification and graph classification have showed
that H2H-GCN achieves competitive results compared with state-of-the-art
Euclidean GCNs and existing hyperbolic GCNs.

##### Acknowledgments.

This work was supported by the Natural Science Foundation of China (NSFC)
under Grants No. 61773062 and No. 62072041.

## References

  * [1] Roy M Anderson, B Anderson, and Robert M May.  Infectious diseases of humans: dynamics and control.  Oxford university press, 1992.
  * [2] Ivana Balazevic, Carl Allen, and Timothy Hospedales.  Multi-relational poincaré graph embeddings.  In Advances in Neural Information Processing Systems (NeurIPS), pages 4463–4473, 2019.
  * [3] Albert-László Barabási and Réka Albert.  Emergence of scaling in random networks.  Science, 286(5439):509–512, 1999.
  * [4] Riccardo Benedetti and Carlo Petronio.  Lectures on hyperbolic geometry.  Springer Science & Business Media, 2012.
  * [5] William M Boothby.  An introduction to differentiable manifolds and riemannian geometry.  In Academic press, 1986.
  * [6] Karsten M Borgwardt, Cheng Soon Ong, Stefan Schönauer, SVN Vishwanathan, Alex J Smola, and Hans-Peter Kriegel.  Protein function prediction via graph kernels.  Bioinformatics, 21(suppl_1):i47–i56, 2005.
  * [7] Joan Bruna, Wojciech Zaremba, Arthur Szlam, and Yann LeCun.  Spectral networks and locally connected networks on graphs.  In International Conference on Learning Representations (ICLR), pages 1–14, 2014.
  * [8] James W Cannon, William J Floyd, Richard Kenyon, Walter R Parry, et al.  Hyperbolic geometry.  Flavors of geometry, 31:59–115, 1997.
  * [9] Ines Chami, Albert Gu, Vaggos Chatziafratis, and Christopher Ré.  From trees to continuous embeddings and back: Hyperbolic hierarchical clustering.  In Advances in Neural Information Processing Systems (NeurIPS), volume 33, 2020.
  * [10] Ines Chami, Zhitao Ying, Christopher Ré, and Jure Leskovec.  Hyperbolic graph convolutional neural networks.  In Advances in Neural Information Processing Systems (NeurIPS), pages 4868–4879, 2019.
  * [11] Aaron Clauset, Cristopher Moore, and Mark EJ Newman.  Hierarchical structure and the prediction of missing links in networks.  Nature, 453(7191):98–101, 2008.
  * [12] Christopher De Sa, Albert Gu, Christopher Ré, and Frederic Sala.  Representation tradeoffs for hyperbolic embeddings.  Proceedings of Machine Learning Research, 80:4460, 2018.
  * [13] Michaël Defferrard, Xavier Bresson, and Pierre Vandergheynst.  Convolutional neural networks on graphs with fast localized spectral filtering.  In Advances in Neural Information Processing Systems (NeurIPS), pages 3844–3852, 2016.
  * [14] Paul D Dobson and Andrew J Doig.  Distinguishing enzyme structures from non-enzymes without alignments.  Journal of molecular biology, 330(4):771–783, 2003.
  * [15] Paul Erdős and Alfréd Rényi.  On random graphs, i.  Publicationes Mathematicae (Debrecen), 6:290–297, 1959.
  * [16] Federico Errica, Marco Podda, Davide Bacciu, and Alessio Micheli.  A fair comparison of graph neural networks for graph classification.  In International Conference on Learning Representations (ICLR), 2020\.
  * [17] Octavian Ganea, Gary Bécigneul, and Thomas Hofmann.  Hyperbolic neural networks.  In Advances in Neural Information Processing Systems (NeurIPS), pages 5345–5355, 2018.
  * [18] Justin Gilmer, Samuel S Schoenholz, Patrick F Riley, Oriol Vinyals, and George E Dahl.  Neural message passing for quantum chemistry.  In International Conference on Machine Learning (ICML), pages 1263–1272, 2017.
  * [19] Will Hamilton, Zhitao Ying, and Jure Leskovec.  Inductive representation learning on large graphs.  In Advances in Neural Information Processing Systems (NeurIPS), pages 1024–1034, 2017.
  * [20] Ruiqi Hu, Shirui Pan, Guodong Long, Qinghua Lu, Liming Zhu, and Jing Jiang.  Going deep: Graph convolutional ladder-shape networks.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), 2020.
  * [21] Zhiwu Huang, Ruiping Wang, Shiguang Shan, Luc Van Gool, and Xilin Chen.  Cross euclidean-to-riemannian metric learning with application to face recognition from video.  IEEE Transactions on Pattern Analysis and Machine Intelligence (T-PAMI), 40(12):2827–2840, 2017.
  * [22] Valentin Khrulkov, Leyla Mirvakhabova, Evgeniya Ustinova, Ivan Oseledets, and Victor Lempitsky.  Hyperbolic image embeddings.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 6418–6428, 2020.
  * [23] Thomas N Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In International Conference on Learning Representations (ICLR), 2017\.
  * [24] Dmitri Krioukov, Fragkiskos Papadopoulos, Maksim Kitsak, Amin Vahdat, and Marián Boguná.  Hyperbolic geometry of complex networks.  Physical Review E, 82(3):036106, 2010.
  * [25] Qi Liu, Maximilian Nickel, and Douwe Kiela.  Hyperbolic graph neural networks.  In Advances in Neural Information Processing Systems (NeurIPS), pages 8230–8241, 2019.
  * [26] Shaoteng Liu, Jingjing Chen, Liangming Pan, Chong-Wah Ngo, Tat-Seng Chua, and Yu-Gang Jiang.  Hyperbolic visual embedding learning for zero-shot recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 9273–9281, 2020.
  * [27] Teng Long, Pascal Mettes, Heng Tao Shen, and Cees GM Snoek.  Searching for actions on the hyperbole.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 1141–1150, 2020.
  * [28] Maximillian Nickel and Douwe Kiela.  Poincaré embeddings for learning hierarchical representations.  In Advances in Neural Information Processing Systems (NeurIPS), pages 6338–6347, 2017.
  * [29] Maximilian Nickel and Douwe Kiela.  Learning continuous hierarchies in the lorentz model of hyperbolic geometry.  In International Conference on Machine Learning (ICML), pages 3776–3785, 2018.
  * [30] Fragkiskos Papadopoulos, Maksim Kitsak, M Ángeles Serrano, Marián Boguná, and Dmitri Krioukov.  Popularity versus similarity in growing networks.  Nature, 489(7417):537–540, 2012.
  * [31] Hongbin Pei, Bingzhe Wei, Kevin Chen-Chuan Chang, Yu Lei, and Bo Yang.  Geom-gcn: Geometric graph convolutional networks.  In In International Conference on Learning Representations (ICLR), 2020.
  * [32] Ida Schomburg, Antje Chang, Christian Ebeling, Marion Gremse, Christian Heldt, Gregor Huhn, and Dietmar Schomburg.  Brenda, the enzyme database: updates and major new developments.  Nucleic acids research, 32(suppl_1):D431–D433, 2004.
  * [33] Prithviraj Sen, Galileo Namata, Mustafa Bilgic, Lise Getoor, Brian Galligher, and Tina Eliassi-Rad.  Collective classification in network data.  AI magazine, 29(3):93–93, 2008.
  * [34] Martin Simonovsky and Nikos Komodakis.  Dynamic edge-conditioned filters in convolutional neural networks on graphs.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), pages 3693–3702, 2017.
  * [35] Rishi Sonthalia and Anna C Gilbert.  Tree! i am no tree! i am a low dimensional hyperbolic embedding.  In Advances in Neural Information Processing Systems (NeurIPS), 2020\.
  * [36] Zequn Sun, Muhao Chen, Wei Hu, Chengming Wang, Jian Dai, and Wei Zhang.  Knowledge association with hyperbolic knowledge graph embeddings.  In Proceedings of the Conference on Empirical Methods in Natural Language Processing (EMNLP), 2020.
  * [37] Oncel Tuzel, Fatih Porikli, and Peter Meer.  Pedestrian detection via classification on riemannian manifolds.  IEEE Transactions on Pattern Analysis and Machine Intelligence (T-PAMI), 30(10):1713–1727, 2008.
  * [38] Abraham A Ungar.  Analytic hyperbolic geometry: Mathematical foundations and applications.  World Scientific, 2005.
  * [39] Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio.  Graph attention networks.  In International Conference on Learning Representations (ICLR), pages 1–12, 2017.
  * [40] Petar Velickovic, William Fedus, William L Hamilton, Pietro Liò, Yoshua Bengio, and R Devon Hjelm.  Deep graph infomax.  In International Conference on Learning Representations (ICLR), 2019\.
  * [41] Duncan J Watts and Steven H Strogatz.  Collective dynamics of ‘small-world’ networks.  Nature, 393(6684):440, 1998.
  * [42] Melanie Weber, Manzil Zaheer, Ankit Singh Rawat, Aditya Menon, and Sanjiv Kumar.  Robust large-margin learning in hyperbolic space.  In Advances in Neural Information Processing Systems (NeurIPS), 2020\.
  * [43] Felix Wu, Tianyi Zhang, Amauri Holanda de Souza Jr, Christopher Fifty, Tao Yu, and Kilian Q Weinberger.  Simplifying graph convolutional networks.  In International Conference on Machine Learning (ICML), 2019.
  * [44] Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka.  How powerful are graph neural networks?  In International Conference on Learning Representations (ICLR), pages 1–17, 2019.
  * [45] Zhitao Ying, Jiaxuan You, Christopher Morris, Xiang Ren, Will Hamilton, and Jure Leskovec.  Hierarchical graph representation learning with differentiable pooling.  In Advances in Neural Information Processing Systems (NeurIPS), pages 4800–4810, 2018.
  * [46] Muhan Zhang, Zhicheng Cui, Marion Neumann, and Yixin Chen.  An end-to-end deep learning architecture for graph classification.  In Proceedings of the AAAI Conference on Artificial Intelligence (AAAI), 2018.
  * [47] Yiding Zhang, Xiao Wang, Xunqiang Jiang, Chuan Shi, and Yanfang Ye.  Hyperbolic graph attention network.  arXiv preprint arXiv:1912.03046, 2019.


# Computational Hierarchy of Elementary Cellular Automata

Barbora Hudcová1,2 Tomáš Mikolov2

1Charles University, Prague
2Czech Institute of Informatics, Robotics and Cybernetics, CTU, Prague

###### Abstract

The complexity of cellular automata is traditionally measured by their
computational capacity. However, it is difficult to choose a challenging set
of computational tasks suitable for the parallel nature of such systems. We
study the ability of automata to emulate one another, and we use this notion
to define such a set of naturally emerging tasks. We present the results for
elementary cellular automata, although the core ideas can be extended to other
computational systems. We compute a graph showing which elementary cellular
automata can be emulated by which and show that certain chaotic automata are
the only ones that cannot emulate any automata non-trivially. Finally, we use
the emulation notion to suggest a novel definition of chaos that we believe is
suitable for discrete computational systems. We believe our work can help
design parallel computational systems that are Turing-complete and also
computationally efficient.

## Introduction

Discrete systems exhibit a wide variety of dynamical behavior ranging from
ordered and easily predictable to a very complex, disordered one. In this
paper, we study cellular automata (CA) that have intriguing visualizations of
their space-time dynamics. Observing them helps us build intuition for
distinguishing the different types of dynamics, as studied by Wolfram, (2002).
Informally, complex CA produce higher-order structures, whereas the space-time
diagrams of chaotic CA are seemingly random. Though CA dynamics has been
studied extensively (Wolfram, (1984), Kůrka, (2009), Wuensche and Lesser,
(2001), Gutowitz, (1990), Zenil, (2009)), it is still a very difficult problem
to formally define the notions of complexity and chaos in CA.

Traditionally, complexity of CA is studied through their computational
capacity (Toffoli, (1977), Cook, (2004)), whereas chaos is studied via
topological dynamics (Devaney, (1989)). As the two approaches are not
connected in any obvious way, it is an interesting open problem of whether
chaotic CA can compute non-trivial tasks (Martinez et al., (2013)).

In this paper, we study the complexity of CA through their computational
capacity; namely, we study their ability to emulate one another; this notion
was first introduced by Mazoyer and Rapaport, (1999). We present the emulation
relation between a pair of CA, and we demonstrate the results on a toy class
of elementary CA; namely, we present their computational hierarchy. The
results inspired us to define a new notion of chaos for discrete dynamical
systems connected to their computational capacity.

## Introducing Cellular Automata

A cellular automaton can be perceived as a $k$-dimensional grid consisting of
identical finite state automata with the same local neighborhood. They are all
updated synchronously in discrete time steps according to a fixed local update
rule which is a function of the neighbors’ states. A formal definition can be
found in Kari, (2005).

### Basic Notions

We say that $\mathbb{Z}$ is a one-dimensional cellular grid and we call its
elements the cells. Let $S$ be a finite set. An $S$-configuration of the grid
is a mapping $c:\mathbb{Z}\rightarrow S$, we write $c_{i}=c(i)$ for each $i$.
We define the  nearest-neighbors relative neighborhood of each cell
$i\in\mathbb{Z}$ to be the triple $(i-1,i,i+1)$. In this paper, we will study
1-dimensional CA with nearest neighbors. Each such CA is characterized by a
tuple $(S,f)$ where $S$ is a finite set of states and $f:S^{3}\rightarrow S$
is a local transition rule of the CA. The global rule of the CA $(S,f)$
operating on an infinite grid is a mapping $F:S^{\mathbb{Z}}\rightarrow
S^{\mathbb{Z}}$ defined as:

$F(c)_{i}=f(c_{i-1},c_{i},c_{i+1}).$

For practical purposes, when observing the CA simulations, we consider the
grid to be of finite size with a periodic boundary condition. In such a case,
we compute the cells in the relative neighborhood modulo the size of the grid.

#### Elementary CA

Elementary cellular automata (ECA) are 1D nearest-neighbors CA with states
$S=\\{0,1\\}$. We identify each local rule $f$ determining an ECA with the
Wolfram number of $f$ defined as:

$2^{0}f(0,0,0)+2^{1}f(0,0,1)+2^{2}f(0,1,0)+\ldots+2^{7}f(1,1,1).$

We will refer to each ECA as a “rule $k$” where $k$ is the corresponding
Wolfram number of its underlying local rule. The class of ECA is a frequently
used toy model for studying different CA properties due to its relatively
small size; there are only 256 of them.

For an ECA operating on a cyclic grid of size $n$ with global rule $F$ we
define the trajectory of a configuration $u\in\\{0,1\\}^{n}$ to be
$(u,F(u),F^{2}(u),\ldots)$. The space-time diagram of such a simulation is
obtained by plotting the configurations as horizontal rows of black and white
squares (corresponding to states 1 and 0) with time progressing downwards.

### CA Complexity via Computational Capacity

Classically, the complexity of a CA is demonstrated by its computational
capacity. Intuitively, we believe that CA capable of computing non-trivial
tasks should be more complex than those that are not. In the past, many
different computational problems were considered, such as the majority
computation task (Mitchell et al., (2000)) or the firing squad synchronization
(Mazoyer, (1986)), a detailed overview was written by Mitchell, (1998).
Nevertheless, the most classical task is the simulation of a computationally
universal system (a Turing machine, a tag system, etc.). Over the years, many
different CA were designed or showed to be Turing complete. However, it seems
unnatural to demonstrate the complexity of an inherently parallel system by
embedding a sequential computational model into it. In our opinion, an ideal
set of benchmark tasks helping us determine the computational capacity of CA
should

  * •

consist of tasks suitable for the parallel computational environment

  * •

be challenging enough

  * •

for a given CA and a task $T$, it should be effectively verifiable whether the
CA can compute $T$.

For a fixed class $\mathcal{C}$ of CA, a natural task is the following:

$\textit{Given a CA in }\mathcal{C},\textit{ how many other CA in
}\mathcal{C}\textit{ can it simulate?}$

The key problem is finding a suitable definition of CA ”simulating” one-
another. Various approaches have been suggested — simulations can be
interpreted via CA coarse-grainings (Israeli and Goldenfeld, (2006)), or
through embedding a local rule to larger size cell-blocks (Mazoyer and
Rapaport, (1999), Ollinger, (2001)). Many interesting theoretical results stem
from such definitions. For instance, Ollinger, (2001) defined a CA simulation
notion which admits a universal CA — one that is able to simulate all other CA
with the same dimensionality; it was designed by Culik II and Albert, (1987).
In this paper, we consider a notion very similar to the one suggested by
Mazoyer and Rapaport, (1999). Compared to their definition, ours is stricter
and, thus, slightly easier to verify. We call it CA emulation and introduce it
in the subsequent section.

## CA Emulations

### Basic Notions

######  Definition 1 (Subautomaton).

Let $\mathrm{ca}_{1}=(S,f)$, $\mathrm{ca}_{2}=(T,g)$ be two nearest-neighbor
1D CA. We say that $\mathrm{ca}_{1}$ is a subautomaton of $\mathrm{ca}_{2}$ if
there exists a one-to-one encoding $e:S\rightarrow T$ such that

$e(f(s_{1},s_{2},s_{3}))=g(e(s_{1}),e(s_{2}),e(s_{3}))$

for all $(s_{1},s_{2},s_{3})\in S^{3}$.

This simply means that we can embed the rule table of $\mathrm{ca}_{1}$ into
the rule table of $\mathrm{ca}_{2}$; the definition is equivalent to saying
that $(S,f)$ is a subalgebra of $(T,g)$.

######  Example 2.

For an ECA $\mathrm{ca}=(\\{0,1\\},f)$, we define the dual ECA
$\mathrm{ca}^{\prime}=(\\{0,1\\},f^{\prime})$ as

$f^{\prime}(b_{1},b_{2},b_{3})=1-f(1-b_{1},1-b_{2},1-b_{3}).$

$f^{\prime}$ simply changes the role of 0 and 1 states. The encoding
$e(b)=1-b$, $b\in\\{0,1\\}$ witnesses that $\mathrm{ca}$ is a subautomaton of
$\mathrm{ca}^{\prime}$ and vice versa.

Figure 1: Space-time diagram of ECA rule 110 on the left and its dual rule 137
on the right.

Each ECA contains at most two subautomata from the ECA class: itself and its
dual ECA. Hence, the subautomaton relation would not produce a rich hierarchy.
It is, however, a key concept for defining the CA emulation. Below, we
restrict the terminology to ECA, but the theory can be generalized to any 1D
CA in a straightforward way.

We will write ${\mathbf{2}}=\\{0,1\\}$ and define
${\mathbf{2}}^{+}=\bigcup_{k=1}^{\infty}{\mathbf{2}}^{k}$ to be the set of all
finite nonempty binary sequences.

######  Definition 3 (Unravelling a Boolean function).

Let $f:{\mathbf{2}}^{3}\rightarrow{\mathbf{2}}$. We define
$\widetilde{f}:{\mathbf{2}}^{+}\rightarrow{\mathbf{2}}^{+}$ as

$\widetilde{f}(b_{1},b_{2},\ldots,b_{n})=f(b_{1},b_{2},b_{3})\ldots
f(b_{n-2},b_{n-1},b_{n})$

for any binary sequence $b_{1}b_{2}\ldots b_{n}$, $n\geq 3$.

By $\widetilde{f}^{k}$ we simply mean the composition

$\widetilde{f}^{k}=\underbrace{\widetilde{f}\circ\widetilde{f}\circ\ldots\circ\widetilde{f}}_{k-\text{times}}.$

Each iteration of $\widetilde{f}$ shortens the input size by 2, therefore we
can notice that when we restrict the domain of $\widetilde{f}^{k}$ we get

${\left.\kern-1.2pt\widetilde{f}^{k}\vphantom{\big{|}}\right|_{{\mathbf{2}}^{3k}}}:{\mathbf{2}}^{k}\times{\mathbf{2}}^{k}\times{\mathbf{2}}^{k}\rightarrow{\mathbf{2}}^{k}.$

Hence, $\widetilde{f}^{k}$ can be interpreted as a ternary function operating
on supercells of $k$ bits and we obtain a CA
$({\mathbf{2}}^{k},\widetilde{f}^{k})$ whose dynamics is completely governed
by the simple local rule $f$. Therefore, each ECA $({\mathbf{2}},f)$ gives
rise to a series of CA

$({\mathbf{2}},f),\,({\mathbf{2}}^{2},\widetilde{f}^{2}),\,({\mathbf{2}}^{3},\widetilde{f}^{3}),\,({\mathbf{2}}^{4},\widetilde{f}^{4}),\,\ldots$

######  Example 4.

Suppose we have $f:{\mathbf{2}}^{3}\rightarrow{\mathbf{2}}$,
$f(b_{1},b_{2},b_{3})=b_{1}\oplus b_{2}\oplus b_{3}$, where $\oplus$ is the
XOR operation. In Figure 2 we show the diagram of $\widetilde{f}^{2}$
computation.

$\widetilde{f}^{2}$ Figure 2:  Diagram of $\widetilde{f}^{2}$ computing on
supercells of size 2.

Let $e:S\rightarrow T$ be a mapping between some finite sets. We define
$\overline{e}:S^{+}\rightarrow T^{+}$ simply as

$\overline{e}(s_{1},s_{2},\ldots,s_{n})=e(s_{1}),e(s_{2}),\ldots,e(s_{n})$

for each $s_{1},\ldots,s_{n}\in S$, $n\in\mathbb{N}$.

######  Definition 5 (ECA Emulation).

We say that $\mathrm{ca}_{1}=({\mathbf{2}},f)$ can be emulated by
$\mathrm{ca}_{2}=({\mathbf{2}},g)$ with a supercell size $k$, if
$\mathrm{ca}_{1}$ is a subautomaton of $({\mathbf{2}}^{k},\widetilde{g}^{k})$.
In such a case, we write $\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$.

Hence, $({\mathbf{2}},f)\leq_{k}({\mathbf{2}},g)$ holds if and only if there
exists a one-to-one encoding

$\mathrm{enc}:{\mathbf{2}}\rightarrow{\mathbf{2}}^{k}$

such that for any initial configuration $c\in{\mathbf{2}}^{3}$ it holds that

$\mathrm{enc}(f(c))=\widetilde{g}^{k}(\overline{\mathrm{enc}}(c)).$

We call such encoding the witnessing encoding. In other words, the following
diagram commutes.

${\mathbf{2}}^{3}$$({\mathbf{2}}^{k})^{3}$${\mathbf{2}}$${\mathbf{2}}^{k}$$\circ$$\overline{\mathrm{enc}}$$\mathrm{enc}$$f$$\widetilde{g}^{k}$
(1)

We say that $\mathrm{ca}_{1}$ can be emulated by $\mathrm{ca}_{2}$ if there
exists some $k$ for which $\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$, and we
write $\mathrm{ca}_{1}\leq\mathrm{ca}_{2}$. If
$\mathrm{ca}_{1}\leq_{1}\mathrm{ca}_{2}$ we say that $\mathrm{ca}_{2}$
emulates $\mathrm{ca}_{1}$ trivially. If $\mathrm{ca}\leq_{k}\mathrm{ca}$ for
$k>1$ we say that the $\mathrm{ca}$ is self-similar.

In the next section, we give some simple proofs of the key properties of the
$\leq$ relation. Namely, we prove that $\leq$ is a preorder and that whenever
$\mathrm{ca}_{1}\leq\mathrm{ca}_{2}$ and $\mathrm{ca}_{1}$ is Turing complete,
then $\mathrm{ca}_{2}$ is also Turing complete.

### Properties of the Emulation Relation

#### Preservation of Turing Completeness

######  Lemma 6.

Suppose that $({\mathbf{2}},f)\leq_{k}({\mathbf{2}},g)$ with a witnessing
encoding $\mathrm{enc}$. Then, for each configuration $c\in{\mathbf{2}}^{l}$,
$l\geq 3$, it holds that
$\overline{\mathrm{enc}}(\widetilde{f}(c))=\widetilde{g}^{k}(\overline{\mathrm{enc}}(c))$.

###### Proof.

Let $c\in{\mathbf{2}}^{l}$ for some $l\geq 3$. Then,
$\widetilde{f}(c)\in{\mathbf{2}}^{l-2}$, and
$\overline{\mathrm{enc}}(\widetilde{f}(c))$ is a sequence of length $l-2$,
each of its elements being a binary $k$-tuple. Below, we show that for each
$1\leq i\leq l-2$ it holds that
$\overline{\mathrm{enc}}(\widetilde{f}(c))_{i}=\widetilde{g}^{k}(\overline{\mathrm{enc}}(c))_{i}$.

$\displaystyle\overline{\mathrm{enc}}(\widetilde{f}(c))_{i}$
$\displaystyle=\mathrm{enc}(\widetilde{f}(c)_{i})$
$\displaystyle=\mathrm{enc}(f(c_{i-1},c_{i},c_{i+1}))$
$\displaystyle=\widetilde{g}^{k}(\mathrm{enc}(c_{i-1}),\mathrm{enc}(c_{i}),\mathrm{enc}(c_{i+1}))$
$\displaystyle=\widetilde{g}^{k}(\overline{\mathrm{enc}}(c))_{i}.$

∎

Using the previous result it can easily be shown that for each
$t\in\mathbb{N}$ and for each $c\in{\mathbf{2}}^{+}$ sufficiently long it
holds that:

$\overline{\mathrm{enc}}(\widetilde{f}^{t}(c))=\widetilde{g}^{kt}(\overline{\mathrm{enc}}(c))$
(2)

simply using induction on t.

Therefore, if $({\mathbf{2}},f)\leq_{k}({\mathbf{2}},g)$ then the diagram (1)
commutes for arbitrarily long configurations and for arbitrarily many
iterations of the functions $\widetilde{f}$ and $\widetilde{g}^{k}$.
Therefore, any space-time diagram produced by
$\mathrm{ca}_{1}=({\mathbf{2}},f)$ can be efficiently encoded by
$\mathrm{enc}$ to a space-time diagram of $\mathrm{ca}_{2}=({\mathbf{2}},g)$.
Hence, if $\mathrm{ca}_{1}$ is Turing complete, $\mathrm{ca}_{2}$ must be
also.

#### Emulation Relation Is a Preorder

Clearly, for each $\mathrm{ca}$ it holds that
$\mathrm{ca}\leq_{1}\mathrm{ca}$. Hence, $\leq$ is reflexive.

######  Lemma 7.

The relation $\leq$ is transitive.

###### Proof.

Suppose that $({\mathbf{2}},f)\leq_{k}({\mathbf{2}},g)$ with witnessing
encoding $\varphi$ and $({\mathbf{2}},g)\leq_{l}({\mathbf{2}},h)$ with
witnessing encoding $\psi$ for some $k,\,l\in\mathbb{N}$. We define
$\mathrm{enc}=\overline{\psi}\circ\varphi:{\mathbf{2}}\rightarrow{\mathbf{2}}^{kl}$.
Then, for any configuration $c\in{\mathbf{2}}^{3}$ we have that:

$\displaystyle\mathrm{enc}(f(c))$
$\displaystyle=\overline{\psi}(\varphi(f(c)))$
$\displaystyle=\overline{\psi}(\widetilde{g}^{k}(\overline{\varphi}(c)))$
$\displaystyle=\widetilde{h}^{kl}(\overline{\psi}\circ\overline{\varphi}(c))$
$\displaystyle=\widetilde{h}^{kl}(\overline{\mathrm{enc}}(c)).$

In the third equality we use the result from (2). ∎

Hence, $\leq$ is transitive and therefore, a preorder.

In contrast, Mazoyer and Rapaport, (1999) define that $({\mathbf{2}},f)$ can
be emulated by $({\mathbf{2}},g)$ if there exist $k,\,l\in\mathbb{N}$ such
that $({\mathbf{2}}^{k},\widetilde{f}^{k})$ is a subautomaton of
$({\mathbf{2}}^{l},\widetilde{g}^{l})$. They present much deeper theoretical
results about the notion in their paper. In contrast, we concentrate on
explicitly computing which ECA can emulate which to form their computational
hierarchy and discuss the results.

#### Emulations as Subalgebras

For each ECA $({\mathbf{2}},g)$ we have a sequence of algebras:

$({\mathbf{2}},g),\,({\mathbf{2}}^{2},\widetilde{g}^{2}),\,({\mathbf{2}}^{3},\widetilde{g}^{3}),\,({\mathbf{2}}^{4},\widetilde{g}^{4}),\,\ldots$

It holds that $({\mathbf{2}},f)\leq_{k}({\mathbf{2}},g)$ if and only if
$({\mathbf{2}}^{k},\widetilde{g}^{k})$ contains a two-element subalgebra
isomorphic to $({\mathbf{2}},f)$; the witnessing encoding being the
corresponding algebra homomorphism.

Israeli and Goldenfeld, (2006) have defined a different notion of CA
simulation called the CA coarse-graining. It is interesting to notice that
their notion is dual to the CA emulation we have defined in this paper. More
specifically, the CA coarse-graining directly corresponds to the congruences
of algebras. Namely, $({\mathbf{2}},g)$ can be coarse-grained into
$({\mathbf{2}},f)$ if and only if there exists some $k\in\mathbb{N}$ such that
$({\mathbf{2}},f)$ is isomorphic to some quotient algebra of
$({\mathbf{2}}^{k},\widetilde{g}^{k})$.

The interpretation of CA emulation as subalgebras comes in handy when
developing an effective algorithm for computing the $\leq_{k}$ for a given
supercell size $k$.

### Emulation Computing Algorithm

Checking whether $\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$ gets infeasible for
large $k$. In this section, we present an algorithm computing the $\leq_{k}$
relation for a given $k$ and compare it to the “naive” algorithm in terms of
efficiency.

Input : $\mathrm{ca}_{1}=({\mathbf{2}},f),\,\mathrm{ca}_{2}=({\mathbf{2}},g)$,
supercell size $k$

Output : a witnessing $\mathrm{enc}$ if
$\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$;

”cannot emulate” message otherwise

_1_ for  _$\mathrm{enc}(0),\,\mathrm{enc}(1)\in{\mathbf{2}}^{k}$_

2 do

3 encoding is valid;

4 for  _$i,j,k\in{\mathbf{2}}$_ do

5 if
_$\mathrm{enc}(f(i,j,k))\neq\widetilde{g}^{k}(\overline{\mathrm{enc}}(i,j,k))$_
then

6 encoding is not valid;

7 break

8 end if

9

10 end for

11 if  _encoding is valid_ then

12 return $\mathrm{enc}$

13 end if

14

15 end for

16 return ”cannot emulate”;

Algorithm 1 Naive algorithm checking whether
$\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$.

Input : $\mathrm{ca}_{2}=({\mathbf{2}},g)$, supercell size $k$

Output : all $\mathrm{ca}_{1}=({\mathbf{2}},f)$ together with a witnessing
$\mathrm{enc}$ such that $\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$

_1_ for  _$u,v\in{\mathbf{2}}^{k}$_

2 do

3 $\\{u,v\\}$ is a valid sublagebra of $({\mathbf{2}}^{k},\widetilde{g}^{k})$;

4 for  _$i,j,k\in\\{u,v\\}$_ do

5 if  _$\widetilde{g}^{k}(i,j,k)\not\in\\{u,v\\}$_ then

6 $\\{u,v\\}$ is not a valid subalgebra;

7 break

8 end if

9

10 end for

11 if  _$\\{u,v\\}$ is a valid sublagebra_ then

12 put $\mathrm{enc}(0)=u,\,\mathrm{enc}(1)=v$;

13 determine $f$ for which
$f=\mathrm{enc}^{-1}\circ\widetilde{g}^{k}\circ\overline{\mathrm{enc}}$;

14 save ($f,\,\mathrm{enc}$)

15 end if

16

17 end for

18

Algorithm 2 Subalgebra algorithm computing all $\mathrm{ca}_{1}$ for which
$\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$.

The naive algorithm simply goes through all possible encodings to determine
whether $\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$. The subalgebra algorithm
computes all two-element subalgebras of $({\mathbf{2}}^{k},\widetilde{g}^{k})$
and thus, determines all $\mathrm{ca}_{1}$ for which
$\mathrm{ca}_{1}\leq\mathrm{ca}_{2}$.

The time complexity with respect to the supercell size $k$ is in
$\mathcal{O}(k\cdot 2^{2k})$ for both algorithms. However, Algorithm 1 needs
to iterate over all ECA to compute the same result as Algorithm 2.
Experimentally, we have indeed observed that Algorithm 1 takes approximately
250 times longer than Algorithm 2 to compute the emulated automata.

## Emulation Hierarchy of ECA

In this section, we present a computational hierarchy based on the emulation
relation. We were able to compute $\leq_{k}$ for supercell size $k$ ranging
from 1 to 11 due to the computational limitations.

In Example 1 we have seen that for $\mathrm{ca}=({\mathbf{2}},f)$ and its dual
$\mathrm{ca}^{\prime}=({\mathbf{2}},f^{\prime})$ it holds that
$\mathrm{ca}\leq_{1}\mathrm{ca}^{\prime}$ and
$\mathrm{ca}^{\prime}\leq_{1}\mathrm{ca}$. Thus, $\leq_{1}$ is an equivalence
relation on the set of ECA, each class containing exactly $f$ and its dual
$f^{\prime}$ (those might coincide for some rules). From the transitivity of
$\leq$ we know that $\mathrm{ca}$ can emulate (and be emulated by) exactly the
same rules as $\mathrm{ca}^{\prime}$. Using a simple program, we obtained 136
different ECA equivalence classes given by $\leq_{1}$. In the following text,
we will identify each equivalence class with the ECA it contains having the
smaller Wolfram number. We show the hierarchy for such representatives and for
supercells of size $k$ ranging from 2 to 11. Specifically, whenever we found
that $\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$ for some
$k\in\\{2,3,\ldots,11\\}$ we represent it by the following diagram.

$\mathrm{ca}_{1}$ $\mathrm{ca}_{2}$

Many ECA are capable of emulating simple rules, such as rule 0 or the identity
rule 204. Adding so many edges to the diagram would make it unreadable. Hence,
we present the hierarchy in three parts.

$184$ Main Part of the Hierarchy:pt $57$ pt $6$ pt $20$ pt $56$ pt $98$ pt
$148$ pt $134$ pt $97$ pt $41$ pt $176$ pt $162$ $34$ pt 2, 10, 66, 42, 46,
172, 130, 138 pt 14 pt 188 pt $48$ pt 16, 24, 80, 112, 116, 144, 208, 216 pt
84 pt 52, 88 pt 152 pt $15$ pt 82, 43, 180 53, 142, 11 pt $85$ pt 154, 113,
212, 81, 26, 27 pt 38, 74 pt $192$ pt 28 pt 94 pt 9, 13, 7, 224 pt $136$ pt 70
pt 156 pt 21, 65, 69, 168 pt $200$ pt 33 pt 1, 5, 29, 37 pt 132 pt $30$ pt
$45$ pt $86$ pt $89$ pt $102$ pt $54$ pt $50$ pt $44$ pt $12$ pt $100$ pt $68$
pt $4$ pt $36$ pt $104$ pt $72$ pt $108$ pt $76$ Figure 3: Emulation Hierarchy
of ECA computed for supercell sizes ranging from 2 to 11; main part of the
diagram. Some edges are depicted in light gray purely for better
understandability. A looped arrow marks self-similar rules. Frequently
Emulated Rules: $170$ pt 2, 6, 10, 11, 14, 20, 25, 26, 27, 34, 35, 37, 38, 40,
41, 42, 43, 46, 54, 56, 57, 58, 65, 66, 74, 81, 84, 85, 97, 98, 106, 113, 130,
134, 138, 142, 148, 154, 162, 168, 170, 172, 184, 188, 212 shift of 1 bit to
the leftpt $240$ pt 6, 9, 11, 14, 15, 16, 20, 24, 37, 41, 43, 48, 49, 52, 53,
54, 56, 57, 61, 80, 81, 82, 84, 88, 96, 97, 98, 112, 113, 114, 116, 120, 134,
142, 144, 148, 152, 176, 180, 184, 208, 212, 216, 224, 240 shift of 1 bit to
the rightpt $204$ pt 1, 4, 5, 12, 13, 18, 19, 22, 23, 28, 29, 33, 36, 37, 44,
50, 51, 54, 62, 68, 69, 70, 72, 73, 76, 77, 78, 92, 94, 100, 104, 108, 110,
118, 122, 124, 126, 132, 140, 146, 156, 164, 172, 178, 196, 200, 204, 216, 232
identitypt $51$ pt 5, 19, 23, 28, 29, 50, 51, 54, 70, 73, 94, 108, 156, 178
negationpt $0$ pt 15, 30, 45, 51, 60, 85, 86, 89, 90, 102, 105, 106, 120, 150,
154, 170, 180, 204, 240 costant 0 rules not capable of emulating 0 with
supercell sizes 2 to 11 are depicted pt $128$ pt 6, 11, 14, 20, 23, 32, 33,
37, 40, 41, 43, 50, 54, 56, 57, 58, 77, 81, 84, 96, 97, 98, 104, 113, 114,
122, 128, 130, 132, 134, 142, 144, 148, 160, 162, 164, 168, 176, 178, 184,
212, 224, 232 ternary AND Figure 4: Emulation Hierarchy of ECA computed for
supercell sizes ranging from 2 to 11, this part shows the most frequently
emulated rules. Emulated Linear Rules: $90$ pt $146$ pt $22$ pt 18, 26, 82,
94, 122, 126, 154, 164, 180 pt $150$ pt $105$ pt $60$ Figure 5: Emulation
Hierarchy of ECA computed for supercell sizes ranging from 2 to 11, this part
shows rules emulating particular linear ECA.

#### Main Part

For compactness, we often show multiple rules in the same node in Figure 3.
However, that does not imply that such rules can emulate one another. Rules in
one node are emulated and can emulate exactly the same rules contained in the
main part in Figure 3. However, they do not necessarily emulate the same
trivial or linear rules in Figures 4 and 5. Therefore, we note that rules in
the same node do not necessarily have an identical computational capacity.

A task frequently studied in the CA environment is to determine the majority
of 0’s and 1’s in the input configuration. The strict version of the task
requires the CA to enter a homogenous state of all 0’s or all 1’s to indicate
the result. It has been shown that no ECA can solve this strict version. If we
relax the requirements on the form of the output, Capcarrere et al., (1996)
have shown that rule 184 solves the majority task exactly. From the main part
of the hierarchy, we can see that by encoding the input configuration in a
simple way, at least nine more ECA and their dual rules can solve the task.

Figure 6: Space-time diagram of ECA rule 184 on the left. The emulated
computation by rule 148 showed on the right. The emulation uses supercells of
size 2 and every second time step is depicted.

#### Frequently Emulated Rules

The CA emulation offers a natural definition of a CA supporting memory: we say
that a CA has a capability of perfect memory if it can emulate the identity
rule 204 (alternatively, we could also tolerate the emulation of shifting
rules or the negation). Hence, we have a practical criterion for checking
whether a given CA can support memory effectively, i.e., with a small
supercell size. Rules, which cannot emulate either of the rules 51, 204, 170,
and 240 with the examined supercell sizes are: 0, 3, 7, 8, 17, 21, 22, 30, 32,
45, 60, 64, 86, 89, 90, 102, 105, 128, 136, 150, 160, 192. We can notice that
they include trivial rules (0, 3, 8, 128), linear rules (60, 90, 105, 150), as
well as the chaotic rules (22, 30, 45, 86, 89).

#### Emulating Linear Rules

We say an ECA is linear, if its local rule $f$ is a linear Boolean function;
i. e., it satisfies $f(x+y)=f(x)+f(y)$ for all $x,y\in{\mathbf{2}}^{3}$. Such
ECA can be studied algebraically which lead to exact results describing e.g.,
their transients or the structure of attractors (Martin et al., (1984), Jen,
(1988)).

Jen, (1999) has shown that the nonlinear rules 18, 126, and 146 can be mapped
onto the linear rule 90. Figure 5 agrees with the results and further shows
other nonlinear rules with this property which makes them interesting
candidates for further algebraic studies.

#### Bottom Part of the Hierarchy

Top and Bottom of the Emulation Hierarchy
---
ECA rules | Number of rules they can emulate
41, 97 |  9 |  Top
6, 20, 54, 57, 134, 148 |  7
14, 37, 56, 84, 94, 98, 156 |  6
0, 3, 8, 17, 60, 64, 90, 102, 105, 106, 120, 150, 170, 204, 240 |  1 |  Bottom
30, 45, 86, 89 |  0
Table 1: ECA Rules at the Top and Bottom of the Emulation Hierarchy computed
for supercells of size 1 to 11.

There are exactly four rules (and their duals) that seem to be unable to
emulate any ECA non-trivially (i.e., with a supercell size larger than one).
Rule 86 is obtained from rule 30 by changing the role of the “left and right
neighbor”. Rule 89 is obtained in the same way from the dual of rule 45. All
the four rules seem to belong to the most disordered ECA according to metrics
such as the compression size of space-time diagrams (Zenil, (2009)) or the
transient classification (Hudcová and Mikolov, (2020)). Due to the seemingly
random patterns it produces, rule 30 was implemented as a pseudorandom number
generator in Mathematica.

It is interesting that the most chaotic ECA seem to be exactly those unable to
emulate any ECA non-trivially. To further explore this, we asked whether the
four chaotic CA can emulate any 1D nearest-neighbors CA non-trivially, not
just the ECA.

Let $(S,f)$ be a CA. Since $S$ is finite, there must exist some $s\in S$ and
$k\in\mathbb{N}$ such that $\widetilde{f}^{k}(s,s,s,\ldots,s)=s$. Hence, for
any CA with one state $(\\{q\\},g)$ the encoding
$\mathrm{enc}:q\mapsto(\underbrace{s,s,\ldots,s}_{k-\text{times}})$ witnesses
that $(S,f)$ can emulate $(\\{q\\},g)$. Therefore, any CA can emulate every CA
with one state. The question is whether the four CA can emulate any CA with
more than one state non-trivially.

We have examined the sequence of algebras

$({\mathbf{2}},f),\,({\mathbf{2}}^{2},\widetilde{f}^{2}),\,({\mathbf{2}}^{3},\widetilde{f}^{3}),\,({\mathbf{2}}^{4},\widetilde{f}^{4}),\,\ldots,({\mathbf{2}}^{11},\widetilde{f}^{11})$

up to supercell size 11 for rules 30, 45, 86, and 89, and computed all their
subalgebras, not only the two-element ones. We found out that all of the four
CA only contain subalgebras with one element, hence, they cannot emulate any
CA with more than one state for supercell sizes 2 to 11.

It remains an open problem whether this would hold for supercells of arbitrary
size. However, we can conclude that if any of the four chaotic ECA can emulate
any non-trivial automaton, it cannot do so effectively, i.e., with a supercell
of small size.

We note that analogous experiments were conducted for the coarse-graining
relation. (Dzwinel and Magiera, (2015)) showed that the rules 30 and 45,
together with their negations, seem to be exactly those that cannot be coarse-
grained non-trivially. These results motivate us to study chaos from a
computational perspective in the next section.

## Chaos In Cellular Automata

Intuitively, a chaotic CA has unpredictable dynamics and seemingly random
space-time diagrams with no apparent higher-order structures forming. Below,
we discuss examples of formal definitions of chaos suggested in the past and
propose a new one based on the CA emulation notion.

For discrete systems, chaos is classically studied through topological
dynamics, and it is typically connected to the system’s sensitivity to initial
conditions (Devaney, (1989), Cattaneo et al., (1999)). Other studies of chaos
in CA examine e.g., their input entropy (Wuensche, (1999)) or properties of
the local rule itself (Wuensche, (2009)); a comprehensive study of ECA,
including their chaotic behavior, was introduced by Chua et al., (2007).

From a great review by Martinez, (2013) one can see that some definitions of
chaos admit either the shift CA (Devaney, (1989)) or linear CA (Cattaneo et
al., (1999)) to be chaotic. In the first case, shifting a configuration by one
bit does not intuitively feel very unpredictable. For linear automata, it is
known that they can be simulated on a computer significantly faster than
general CA (Robinson, (1987)). Hence, their dynamics can be predicted quite
efficiently. Below, we propose a new, much stricter definition of chaotic
behavior.

######  Definition 8.

We call a 1D nearest-neighbors CA computationally chaotic if it cannot emulate
any 1D nearest-neighbors CA with more than one state non-trivially.

Proving that a particular CA is computationally chaotic would, in principle,
require verifying an infinite amount of conditions and might require deep
theoretical insight into the dynamics of the CA. However, it is useful as a
new theoretical notion, which formalizes the unpredictability of a CA. Indeed,
if $\mathrm{ca}_{1}\leq_{k}\mathrm{ca}_{2}$ then for a subset of initial
conditions, the result of $\mathrm{ca}_{2}$ run for $k$ time-steps is
equivalent to simply running $\mathrm{ca}_{1}$ for 1 time-step. Hence, at
least some part of its dynamics can be predicted more effectively. In
contrast, computationally chaotic CA cannot contain any simpler CA as a
substructure in its space-time dynamics. This suggests that no part of its
dynamics can be predicted efficiently. Hence, the definition seems to agree
with the intuitive understanding of unpredictability.

In CA with trivial dynamics, structures tend to die out quite quickly. This
enables such rules to emulate either the rule 0 non-trivially (this can be
seen from Figure 4) or to be self-similar (e.g., rules 0, 240, or the shift
rule). Thus, they are not computationally chaotic. As linear CA are self-
similar, it follows that they are not computationally chaotic either.

From the results we presented, it follows that the only ECA that could be
computationally chaotic are rules 30 and 45, together with their symmetrical
rules. Hence, we can form the following hypothesis.

######  Hypothesis 9.

ECA rules 30 and 45 are computationally chaotic.

Proving this result would imply that such rules cannot compute any task in the
sense of CA emulation.

We note that the definition can be extended to CA in arbitrary dimensions with
various neighborhood shapes in a straightforward way. We also note that the
definition does not depend on whether the CA operates on a finite or infinite
grid and is purely a property of its local rule.

Most definitions of chaos in discrete systems are not related to their
computational capacity. Hence, there is an interesting ongoing debate whether
chaotic CA can compute any non-trivial tasks. In contrast, the notion of
computational chaos is strongly connected to the computational capacity of the
CA. If a CA is computationally chaotic, it cannot compute the dynamics of any
other CA. Hence, we might conclude that being computationally chaotic directly
implies the inability to compute a non-trivial task.

## Conclusion

We studied the notion of CA emulation as a method of determining the
computational capacity of CA. We showed that the CA emulation relation forms a
preorder on the ECA class and presented an approximation of the emulation
hierarchy produced by the preorder. We did notice that the most chaotic CA
seem to be the minimal elements of the hierarchy. This inspired us to
introduce a new definition of chaos in the CA environment. Our notion of chaos
is novel because it is connected to the computational capacity of a system. In
contrast to previous concepts of chaos, our definition does not regard linear
and shifting automata as chaotic. This agrees with the results that the
dynamics of such CA can be predicted more efficiently than for general CA.

The emulation relation can be defined for CA in any dimension and with an
arbitrary neighborhood. Though its computation requires verifying infinitely
many conditions, we can compute it just for supercells of small size.
Verifying such small supercell values helps us determine whether a CA can
emulate any others effectively.

A pair of ECA obtained by changing the role of their “left and right neighbor”
has equivalent dynamics; however, the emulation relation we presented is
unable to discover such an equivalence. For this reason, it is meaningful to
study other possible definitions of a CA subautomaton. It would be very
interesting to examine whether the different variants of the definition would
change the computational hierarchy results significantly.

As a possible application, we can use the method of CA emulation to construct
Turing-complete CA without having to give elaborate proofs of this fact. We
would simply embed a well-known Turing-complete CA into a newly constructed CA
with possibly much richer dynamics. It could be interesting to design a CA
capable of emulating many different CA and study the dynamics of such a rich
system.

## Acknowledgements

Our work was supported by Grant Schemes at CU, reg. no.
CZ.02.2.69/0.0/0.0/19_073/0016935, the Ministry of Education, Youth and Sports
within the dedicated program ERC CZ under the project POSTMAN with reference
LL1902, by the Czech project AI$\&$Reasoning CZ.02.1.01/0.0/0.0/15_003/0000466
and the European Regional Development Fund, by SVV-2020-260589, and is part of
the RICAIP project that has received funding from the European Union’s Horizon
2020 research and innovation programme under grant agreement No. 857306.

## References

  * Capcarrere et al., (1996) Capcarrere, M. S., Sipper, M., and Tomassini, M. (1996).  Two-state, $\mathit{r}\phantom{\rule{0.0pt}{0.0pt}}=\phantom{\rule{0.0pt}{0.0pt}}1$ cellular automaton that classifies density.  Phys. Rev. Lett., 77:4969–4971.
  * Cattaneo et al., (1999) Cattaneo, G., Formenti, E., Margara, L., and Mauri, G. (1999).  On the dynamical behavior of chaotic cellular automata.  Theoretical Computer Science, 217(1):31–51.
  * Chua et al., (2007) Chua, L. O., Guan, J., Sbitnev, V. I., and Shin, J. (2007).  A nonlinear dynamics perspective of wolfram’s new kind of science. part vii: Isles of eden.  International Journal of Bifurcation and Chaos, 17(09):2839–3012.
  * Cook, (2004) Cook, M. (2004).  Universality in elementary cellular automata.  Complex Systems, 15.
  * Culik II and Albert, (1987) Culik II, K. and Albert, J. (1987).  A simple universal cellular automaton and its one-way and totalistic version.  Complex Systems, 1:1–16.
  * Devaney, (1989) Devaney, R. L. (1989).  An Introduction To Chaotic Dynamical Systems, Second Edition.  Addison-Wesley advanced book program. Avalon Publishing.
  * Dzwinel and Magiera, (2015) Dzwinel, W. and Magiera, K. (2015).  Irreducible elementary cellular automata found.  Journal of Computational Science, 11:300–308.
  * Gutowitz, (1990) Gutowitz, H. A. (1990).  A hierarchical classification of cellular automata.  Physica D: Nonlinear Phenomena, 45(1):136–156.
  * Hudcová and Mikolov, (2020) Hudcová, B. and Mikolov, T. (2020).  Classification of complex systems based on transients.  Artificial Life Conference Proceedings, (32):367–375.
  * Israeli and Goldenfeld, (2006) Israeli, N. and Goldenfeld, N. (2006).  Coarse-graining of cellular automata, emergence, and the predictability of complex systems.  Physical Review E, 73(2).
  * Jen, (1988) Jen, E. (1988).  Linear cellular automata and recurring sequences in finite fields.  Communications in Mathematical Physics, 119:13–28.
  * Jen, (1999) Jen, E. (1999).  Exact solvability and quasiperiodicity of one-dimensional cellular automata.  Nonlinearity, 4:251.
  * Kari, (2005) Kari, J. (2005).  Theory of cellular automata: A survey.  Theoretical Computer Science, 334(1):3 – 33.
  * Kůrka, (2009) Kůrka, P. (2009).  Topological dynamics of cellular automata.  Cellular Automata. Encyclopedia of Complexity and Systems Science Series.
  * Martin et al., (1984) Martin, O., Odlyzko, A., and Wolfram, S. (1984).  Algebraic properties of cellular automata.  Communications in Mathematical Physics, 93.
  * Martinez, (2013) Martinez, G. J. (2013).  A note on elementary cellular automata classification.  Journal of cellular automata, 8.
  * Martinez et al., (2013) Martinez, G. J., Seck Tuoh Mora, J., and Zenil, H. (2013).  Computation and universality: Class iv versus class iii cellular automata.  Journal of cellular automata.
  * Mazoyer, (1986) Mazoyer, J. (1986).  An overview of the firing squad synchronization problem.  Automata Networks, LITP Spring School on Theoretical Computer Science Proceedings, 316.
  * Mazoyer and Rapaport, (1999) Mazoyer, J. and Rapaport, I. (1999).  Inducing an order on cellular automata by a grouping operation.  Discrete Applied Mathematics, 91(1):177–196.
  * Mitchell, (1998) Mitchell, M. (1998).  Computation in cellular automata: A selected review.  Non‐Standard Computation, pages 95–140.
  * Mitchell et al., (2000) Mitchell, M., Crutchfield, J., and Das, R. (2000).  Evolving cellular automata with genetic algorithms: A review of recent work.  First Int. Conf. on Evolutionary Computation and Its Applications, 1.
  * Ollinger, (2001) Ollinger, N. (2001).  Toward an algorithmic classification of cellular automata.  LIP 222001-10, (2).
  * Robinson, (1987) Robinson, A. (1987).  Fast computation of additive cellular automata.  Complex Systems, 1(1):211–216.
  * Toffoli, (1977) Toffoli, T. (1977).  Computation and construction universality of reversible cellular automata.  Journal of Computer and System Sciences, 15(2):213 – 231.
  * Wolfram, (1984) Wolfram, S. (1984).  Universality and complexity in cellular automata.  Physica D: Nonlinear Phenomena, 10(1):1–35.
  * Wolfram, (2002) Wolfram, S. (2002).  A New Kind of Science.  Wolfram Media, Champaign, USA.
  * Wuensche, (1999) Wuensche, A. (1999).  Classifying cellular automata automatically: Finding gliders, filtering, and relating space-time patterns, attractor basins, and thez parameter.  Complexity, 4:47–66.
  * Wuensche, (2009) Wuensche, A. (2009).  Cellular automata encryption: the reverse algorithm, z-parameter and chain-rules.  Parallel Processing Letters, 19:283–297.
  * Wuensche and Lesser, (2001) Wuensche, A. and Lesser, M. (2001).  The global dynamics of celullar automata: An atlas of basin of attraction fields of one-dimensional cellular automata.  J. Artificial Societies and Social Simulation, 4.
  * Zenil, (2009) Zenil, H. (2009).  Compression-based investigation of the dynamical properties of cellular automata and other systems.  Computing Research Repository - CORR, 19.