<figure><img src="image_1.png" /></figure>

HAL
open science

# Unsupervised Learning in Complex Systems

Hugo Cisneros

► To cite this version:

Hugo Cisneros. Unsupervised Learning in Complex Systems. Artificial Intelligence [cs.AI]. École normale supérieure, 2023. English. NNT: tel-04159511

HAL Id: tel-04159511

https://inria.hal.science/tel-04159511

Submitted on 27 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Unsupervised Learning in Complex Systems

Soutenue par

Hugo Cisneros

Le 15 Mai 2023

École doctorale n°386

Sciences Mathématiques de
Paris Centre

Spécialité

Informatique

## Composition du jury :

Lenka Zdeborová

École Polytechnique Fédérale de
Lausanne
Président du jury

Stefano Nichele

Østfold University College
Rapporteur

Hiroki Sayama

Binghamton University
Rapporteur

Josef Sivic

École Normale Supérieure, Inria
CIIRC, Czech Technical University
Directeur de thèse

Tomas Mikolov

CIIRC, Czech Technical University
Co-encadrant de thèse

<figure><img src="image_3.png" /></figure>
<figure><img src="image_4.png" /></figure>
<figure><img src="image_5.png" /></figure>

# Résumé

Dans cette thèse, nous explorons l'utilisation de systèmes complexes pour étudier l'apprentissage et l'adaptation dans les systèmes naturels et artificiels. L'objectif est de développer des systèmes autonomes capables d'apprendre sans supervision, de se développer de manière autonome et de devenir de plus en plus complexes avec le temps. Les systèmes complexes apparaissent comme un cadre adapté pour comprendre ces phénomènes en raison de leur capacité à croître en complexité. Être en mesure de construire des algorithmes d'apprentissage qui nécessitent peu ou pas de supervision permettrait une plus grande flexibilité et adaptabilité dans diverses applications. En tentant de comprendre les principes fondamentaux de l'apprentissage dans les systèmes complexes, nous espérons améliorer notre capacité à concevoir et à mettre en œuvre des algorithmes d'apprentissage à l'avenir. Cette thèse apporte les contributions suivantes : le développement d'une métrique de complexité générale que nous appliquons pour rechercher des systèmes complexes qui présentent une croissance de complexité, l'introduction d'une méthode pour étudier les calculs dans des systèmes complexes à grande échelle, et le développement d'une métrique pour l'efficacité de l'apprentissage ainsi que l'introduction d'un jeu de données de référence pour évaluer la vitesse des algorithmes d'apprentissage. Nos résultats contribuent significativement à notre compréhension de l'apprentissage et de l'adaptation dans les systèmes naturels et artificiels. De plus, notre approche contribue à une nouvelle direction prometteuse pour la recherche dans ce domaine. Nous espérons que ces résultats inspireront le développement d'algorithmes d'apprentissage plus efficaces et plus performants à l'avenir.

Apprentissage automatique ★ Intelligence artificielle ★ Systèmes complexes ★ Complexité
<figure><img src="image_6.png" /></figure>

Abstract

In this thesis, we explore the use of complex systems to study learning and adaptation
in natural and artificial systems. The goal is to develop autonomous systems that can
learn without supervision, develop on their own, and become increasingly complex
over time. Complex systems are identified as a suitable framework for understanding
these phenomena due to their ability to exhibit growth of complexity. Being able to
build learning algorithms that require limited to no supervision would enable greater
flexibility and adaptability in various applications. By understanding the fundamental
principles of learning in complex systems, we hope to advance our ability to design and
implement practical learning algorithms in the future. This thesis makes the following key
contributions: the development of a general complexity metric that we apply to search for
complex systems that exhibit growth of complexity, the introduction of a coarse-graining
method to study computations in large-scale complex systems, and the development of
a metric for learning efficiency as well as a benchmark dataset for evaluating the speed
of learning algorithms. Our findings add substantially to our understanding of learning
and adaptation in natural and artificial systems. Moreover, our approach contributes to
a promising new direction for research in this area. We hope these findings will inspire
the development of more effective and efficient learning algorithms in the future.
<figure><img src="image_7.png" /></figure>

# Acknowledgments

I express my sincere gratitude to Professor Josef Sivic and Tomas Mikolov for their invaluable guidance, enthusiasm, and endless encouragement, without which this thesis would not have been possible.

I am very grateful to a number of people who have directly or indirectly influenced and helped me in my work, both as colleagues and friends: Jelle, Barbora, Kateryna, Teven, and David. Special thanks to Barbora Hudcová and Josef Sivic for their meticulous proofreading of parts of this thesis.

I would also like to thank everyone at the Foundational AI Lab and the Automated Reasoning Lab of the Czech Institute of Informatics, Robotics and Cybernetics (CIIRC) for creating a welcoming and supportive work environment.

Most importantly, I would like to express my deep appreciation to my parents, family, and friends for their support and understanding throughout my academic journey.

Lastly, I would like to express my heartfelt thanks to Brune for her unwavering patience, help, support, and constant inspiration.
<figure><img src="image_8.png" /></figure>

# Contents

Résumé i

Abstract ii

Acknowledgments iii

1 Introduction 1

1.1 Goal 1

1.2 Motivation 4

1.3 Challenges 6

1.3.1 Design of a complex system 8

1.3.2 Complex systems control 8

1.3.3 Open-ended evolution without objectives 9

1.3.4 Complex systems inputs and outputs 10

1.4 Contributions 10

1.5 Thesis overview 11

1.6 Publications and software 12

2 Background 13

2.1 Cellular automata 13

2.1.1 Definition 13

2.1.2 Classification of cellular automata 16

2.1.3 Cellular automata variants 18

2.1.4 Parametrizing and sampling CA rules 23

2.2 Cellular automata and Recurrent neural networks 26

2.2.1 RNN formalism 27
<table>
  <tbody>
    <tr>
      <td>2.2.2</td>
      <td>Transition rule as a set of convolutions</td>
      <td>28</td>
    </tr>
    <tr>
      <td>2.2.3</td>
      <td>Adding inputs and outputs</td>
      <td>30</td>
    </tr>
    <tr>
      <td>2.2.4</td>
      <td>Consequences</td>
      <td>34</td>
    </tr>
    <tr>
      <td>2.2.5</td>
      <td>Beyond the naive rule representation</td>
      <td>35</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Reservoir computing</td>
      <td>35</td>
    </tr>
    <tr>
      <td>2.3.1</td>
      <td>Reservoir computing applications</td>
      <td>38</td>
    </tr>
    <tr>
      <td>2.3.2</td>
      <td>Echo-state network</td>
      <td>38</td>
    </tr>
    <tr>
      <td>2.3.3</td>
      <td>Reservoir cellular automata</td>
      <td>40</td>
    </tr>
    <tr>
      <td><b>3</b></td>
      <td><b>Literature Review</b></td>
      <td><b>43</b></td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Measuring complexity</td>
      <td>43</td>
    </tr>
    <tr>
      <td>3.1.1</td>
      <td>Information content</td>
      <td>45</td>
    </tr>
    <tr>
      <td>3.1.2</td>
      <td>Shannon Entropy</td>
      <td>45</td>
    </tr>
    <tr>
      <td>3.1.3</td>
      <td>Rényi Entropy</td>
      <td>46</td>
    </tr>
    <tr>
      <td>3.1.4</td>
      <td>Mutual information</td>
      <td>46</td>
    </tr>
    <tr>
      <td>3.1.5</td>
      <td>Computational complexity</td>
      <td>47</td>
    </tr>
    <tr>
      <td>3.1.6</td>
      <td>Solomonoff-Kolmogorov-Chaitin complexity (algorithmic complexity)</td>
      <td>47</td>
    </tr>
    <tr>
      <td>3.1.7</td>
      <td>Lempel-Ziv complexity</td>
      <td>48</td>
    </tr>
    <tr>
      <td>3.1.8</td>
      <td>Logical Depth</td>
      <td>49</td>
    </tr>
    <tr>
      <td>3.1.9</td>
      <td>Thermodynamic Depth</td>
      <td>49</td>
    </tr>
    <tr>
      <td>3.1.10</td>
      <td>Epsilon-machines</td>
      <td>50</td>
    </tr>
    <tr>
      <td>3.1.11</td>
      <td>Sophistication</td>
      <td>51</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Emergence</td>
      <td>52</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Evolutionary algorithms</td>
      <td>53</td>
    </tr>
    <tr>
      <td>3.3.1</td>
      <td>Genetic programming</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3.3.2</td>
      <td>Novelty search</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Open-ended evolution</td>
      <td>57</td>
    </tr>
    <tr>
      <td>3.4.1</td>
      <td>Defining open-endedness</td>
      <td>58</td>
    </tr>
    <tr>
      <td>3.4.2</td>
      <td>Open-endedness in cellular automata</td>
      <td>60</td>
    </tr>
    <tr>
      <td>3.4.3</td>
      <td>Artificial chemistries</td>
      <td>60</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Computing with complex systems</td>
      <td>62</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>3.5.1 Computing in cellular automata</td>
      <td>62</td>
    </tr>
    <tr>
      <td>3.5.2 Amorphous computing</td>
      <td>69</td>
    </tr>
    <tr>
      <td><strong>4 Measuring complexity in evolving complex systems</strong></td>
      <td><strong>71</strong></td>
    </tr>
    <tr>
      <td>4.1 Introduction</td>
      <td>71</td>
    </tr>
    <tr>
      <td>4.2 Related work</td>
      <td>72</td>
    </tr>
    <tr>
      <td>4.2.1 Artificial life and open-ended evolution</td>
      <td>72</td>
    </tr>
    <tr>
      <td>4.2.2 Cellular automata</td>
      <td>73</td>
    </tr>
    <tr>
      <td>4.2.3 Compression and complexity</td>
      <td>74</td>
    </tr>
    <tr>
      <td>4.3 Compression-based metric</td>
      <td>75</td>
    </tr>
    <tr>
      <td>4.4 Predictor-based metric</td>
      <td>77</td>
    </tr>
    <tr>
      <td>4.4.1 Joint compression</td>
      <td>78</td>
    </tr>
    <tr>
      <td>4.4.2 Count-based predictor</td>
      <td>79</td>
    </tr>
    <tr>
      <td>4.4.3 Neural network based predictor</td>
      <td>83</td>
    </tr>
    <tr>
      <td>4.5 Experiments</td>
      <td>85</td>
    </tr>
    <tr>
      <td>4.6 Discussion</td>
      <td>86</td>
    </tr>
    <tr>
      <td>4.7 Conclusion</td>
      <td>89</td>
    </tr>
    <tr>
      <td><strong>5 Visualization of computations in large-scale complex systems</strong></td>
      <td><strong>92</strong></td>
    </tr>
    <tr>
      <td>5.1 Introduction</td>
      <td>92</td>
    </tr>
    <tr>
      <td>5.2 Related work</td>
      <td>94</td>
    </tr>
    <tr>
      <td>5.2.1 Coarse-graining in cellular automata</td>
      <td>94</td>
    </tr>
    <tr>
      <td>5.2.2 Filtering</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.2.3 Scaling-up cellular automata</td>
      <td>95</td>
    </tr>
    <tr>
      <td>5.3 Proposed coarse-graining of cellular automata</td>
      <td>96</td>
    </tr>
    <tr>
      <td>5.3.1 Frequency histogram coarse-graining</td>
      <td>97</td>
    </tr>
    <tr>
      <td>5.3.2 Clustering</td>
      <td>100</td>
    </tr>
    <tr>
      <td>5.3.3 Autoencoders for coarse-graining</td>
      <td>100</td>
    </tr>
    <tr>
      <td>5.4 Results</td>
      <td>101</td>
    </tr>
    <tr>
      <td>5.4.1 Domains and filtering</td>
      <td>102</td>
    </tr>
    <tr>
      <td>5.4.2 Results on elementary cellular automata</td>
      <td>103</td>
    </tr>
    <tr>
      <td>5.4.3 Complexity metrics and coarse-graining</td>
      <td>105</td>
    </tr>
  </tbody>
</table>
<table><tbody><tr><td>5.4.4</td><td>Discussion</td><td>107</td></tr><tr><td>5.5</td><td>Conclusion</td><td>108</td></tr><tr><td>6</td><td>Learning efficiency in deep reservoir computing</td><td>113</td></tr><tr><td>6.1</td><td>Introduction</td><td>114</td></tr><tr><td>6.2</td><td>Related work</td><td>115</td></tr><tr><td>6.3</td><td>A benchmark for reservoir computing</td><td>117</td></tr><tr><td>6.3.1</td><td>Performance metric</td><td>119</td></tr><tr><td>6.3.2</td><td>Description of tasks in the benchmark</td><td>120</td></tr><tr><td>6.4</td><td>Standard language classification task</td><td>124</td></tr><tr><td>6.5</td><td>Reservoir computing and reservoir cellular automata</td><td>126</td></tr><tr><td>6.6</td><td>New benchmark: compared methods</td><td>129</td></tr><tr><td>6.6.1</td><td>Fully trained sequential models (RNN, LSTM and Transformer)</td><td>130</td></tr><tr><td>6.6.2</td><td>Echo-state networks and reservoir cellular automata</td><td>130</td></tr><tr><td>6.6.3</td><td>Experimental set-up</td><td>131</td></tr><tr><td>6.7</td><td>Human performance evaluation</td><td>131</td></tr><tr><td>6.8</td><td>Experimental parameters</td><td>133</td></tr><tr><td>6.8.1</td><td>Task generation parameters</td><td>133</td></tr><tr><td>6.9</td><td>Results</td><td>134</td></tr><tr><td>6.10</td><td>Reservoir cellular automata (ReCA) results</td><td>139</td></tr><tr><td>6.10.1</td><td>Binary sequences</td><td>139</td></tr><tr><td>6.10.2</td><td>Symbol counting</td><td>140</td></tr><tr><td>6.11</td><td>Conclusions</td><td>141</td></tr><tr><td>7</td><td>Conclusion</td><td>143</td></tr><tr><td>7.1</td><td>Contributions</td><td>143</td></tr><tr><td>8</td><td>Future Work</td><td>145</td></tr><tr><td>8.1</td><td>Evolving cellular automata for complexity</td><td>145</td></tr><tr><td>8.2</td><td>Feedback reservoir computing</td><td>146</td></tr><tr><td>8.3</td><td>Learning in dynamical systems</td><td>147</td></tr><tr><td colspan="2">Bibliography</td><td>151</td></tr></tbody></table>
## 1.1 Goal

Complex systems are systems composed of many interconnected elements that interact with each other, often in non-linear ways, giving rise to emergent properties that cannot be explained by the properties of the individual elements alone. In many naturally occurring complex systems, learning is a key process that allows the system to adapt and evolve over time. In this thesis, we explore complex systems as a framework for studying learning and adaptation in natural and artificial systems. The aim of this thesis is to develop methods for studying and using computations that take place in complex dynamical systems to eventually create learning algorithms that require limited to no supervision. Complex systems often exhibit features such as self-organization, adaptation, feedback loops, and the ability to undergo phase transitions, all of which can make them challenging to model and predict. Examples of complex systems that occur naturally include ecosystems, weather patterns, financial markets, social networks, and biological organisms. The objective of this thesis is broken down into the following subgoals:

1. The first subgoal is to identify complex dynamical systems that have the potential to display emergent open-ended growth, which refers to the ability of a system to spontaneously generate new levels of complexity over time without reaching a state of equilibrium. This growth is often associated with evolutionary-like properties, such as variation, selection, and inheritance. There are many ways to define complex systems, and as illustrated in Figure 1.1, some may exhibit more interesting and promising behaviors than others. A cellular automaton (CA) is an instance of a complex system that can be described as a grid of cells that can take on a finite set of states and change their state over time according to
a set of rules that depend on the states of their neighboring cells. Interesting Cellular automata (CAs), such as the one shown in Figure 1.1b may be hard to find depending on how the search space is defined. The work presented in Chapter 4 aims to construct a metric of complexity that can help to identify interestingly behaving complex dynamical systems.

<figure><img src="image_9.png" /><figcaption>**Figure 1.1** – Two examples of complex systems (cellular automata) with different behavior types. 1.1b and 1.1a show a single state of a randomly initialized 2D CA simulated for a fixed number of steps. Some CAs appear more promising than others for the design of unsupervised learning systems because of their emergent complex structures (visible in 1.1b), whereas the CA in 1.1a seems to behave randomly.</figcaption></figure>

2. The second subgoal is to measure the fraction of systems that have the most complex and rapidly evolving behavior. Defining these notions is also part of the goal. This sub-goal is distinct from the first one, which focuses on complex system design in general. Instead, this sub-goal involves establishing approaches to evaluate growth in complexity. We believe that systems with the most complex and rapidly evolving behavior are promising for further use since they may exhibit open-ended complexity growth. In Chapters 4 and 5, we present different methods to measure the evolving complexity and understand when the complexity increases over time. Chapter 5 of this thesis is dedicated to investigating the importance of multiscale analysis for the complexity of cellular automata, which is the process of studying a system at different levels of detail or resolution. Additionally, we identify complex systems with behavior that changes when we manipulate the
scale of the system, either by increasing or decreasing its size.

<figure><img src="image_10.png" /><figcaption>**Figure 1.2** – Filtering the behavior of elementary cellular automaton rule 18. This uncovers a highly structured behavior with the propagation of area boundaries within the apparent randomness. (a) shows 300 timesteps of a randomly initialized rule 18 simulation. Notice the complex structures made visible in (b) with a filtering method. This Figure is discussed in more detail in Chapter 5.</figcaption></figure>

3. Apply promising systems to challenging learning tasks where classical machine learning models may fail or become less efficient. The goal is to define various ways to apply evolving complex dynamical systems to some standard learning tasks and to find out if it improves the performance or efficiency of the learning algorithm. An example application that we explore in Chapter 6 is illustrated in Figure 1.3, where a CA is used to implement a language model, a probabilistic model that is used to generate text to complete sentences. In Chapter 2, we explore the similarities between CAs networks and highlight the numerous potential applications of CAs.

In this thesis, we do not focus our attention on what the machine learning community commonly refers to as “unsupervised learning”, that is, learning underlying structures from unlabeled data (Hinton and Sejnowski 1999). In machine learning, supervision refers to a label, a number, or a collection of numbers that represents an expected outcome associated with some input data. Our goal is to construct models that can develop autonomously without necessarily needing any external data or interactions. Therefore, we first seek systems that behave that way on their own and postpone the issue of learning to optimize a particular objective function to the end of the thesis
<figure><img src="image_11.png" /><figcaption>**Figure 1.3** – An example of using a cellular automaton and reservoir computing to implement a language model. Tokens are encoded within the cellular automaton internal state. The cellular automaton update rule is then applied, and an output value is decoded from the final state. Each vertical bar represents two consecutive internal CA states, which are enriched with the encoded input.</figcaption></figure>

(Chapter 6). We are particularly interested in systems that can develop through internal evolution rules without any input, which is the case for many complex systems.

Throughout the thesis, we work toward achieving these goals, focusing on one particular complex system: the cellular automaton (Von Neumann and Burks 1966). This model has been extensively studied due to its simple definition and its ability to simulate a wide range of complex behaviors. We provide a detailed description of cellular automata in Section 2.1. This thesis provides new insights into the role of learning in complex systems and open-ended evolution and demonstrates the potential of cellular automata as a model for studying these phenomena.

## 1.2 Motivation

It is possible that some form of evolutionary mechanism may be necessary in order to achieve advanced forms of artificial intelligence (AI), particularly if the goal is to create intelligent systems that can adapt and learn in complex and changing environments. Complex dynamical systems could be the key to overcoming the problems with existing learning algorithms, such as difficulties in generalization, robustness, or the ability to learn continuously (Parisi et al. 2019). The natural intelligence of biological systems
seems to depend on emerging properties selected through evolution. For example, bio-
logical life exhibits a pattern of major evolutionary transitions in which autonomously
replicating entities at the lower level merged to form a single more complex corporate
body (Lorenz, Jeng, and Deem 2011). This is believed to have occurred for the emergence
of eukaryotic cells, as well as for multicellular life (Hammerschmidt et al. 2014). All
complex and diverse biological entities have presumably emerged from a single common
ancestor and, even before, from inorganic components present on the surface of the Earth
(J. M. Smith and Szathmary 2000; Woese 1998).

So far, it remains uncertain which algorithmic characteristics might enable an artificial
system to exhibit a path comparable to the natural evolutionary process within its
state space, which involves the transition from basic constituents to intricate entities.
Producing emerging phenomena similar to those of nature *in silico* is a long-standing
challenge. Many algorithms try to mimic evolutionary properties to solve particular tasks,
called evolutionary algorithms (Bäck and Schwefel 1993; Fogel, Owens, and Walsh 1966;
Miller, Todd, and Hegde 1989). However, most of these methods focus on searching the
space of solutions using high-level evolutionary mechanisms, such as genetic mutations
and crossovers. Although effective in solving precise tasks, these methods often obscure
another crucial component of natural evolution by using an explicit fitness function and
hard-coded primitives.

Most existing machine learning algorithms rely on the choice of an objective function: a clearly defined mapping from the current state and parameters of a model to a real value, which indicates the performance of that model. The function depends on the objective of the model. For a supervised learning problem, we may count the number of misclassified objects or the distance between the predictions and the expected results. Even in unsupervised learning, the family of algorithms used for learning from unlabeled data, objectives are still central. For example, the well-known K-means clustering algorithm minimizes the sum of square distances of data points to cluster centers.

This reliance on objective functions creates two main issues: (i) The objective is
not always clearly defined or can be too broad for general-purpose applications. For
example, a possible objective function of a walking robot could be “not fall when stepping
through its surrounding environment”. This function is impractical to define and will vary
greatly depending on the parameters of the environment (e.g. depending on the terrain,
malfunctioning or missing limbs of the robots). Indirect rewards such as “head not
touching the ground” may help, but it does not cover all possible ways the model could
fail. This problem of defining an objective is also present within our goal of designing a
general-purpose, autonomously developing learning algorithm. Furthermore, (ii) using
predefined functions as goals can be counterproductive because, as many examples in
nature demonstrate, robust paths to complex objectives are often deceptive. They involve
developing in unexpected directions that may initially seem to be against the original
goal (Stanley and Lehman 2015).

In this thesis, the term *unsupervised* refers to a form of learning with no predefined
objective. Like in natural evolution, we expect true unsupervised algorithms to develop
new features autonomously and become progressively more complex over time. Such
algorithms would regularly learn to solve problems on their own without the need to be
explicitly guided towards the problem solutions, thereby discovering robust and diverse
solutions to deceptive problems.

1.3 Challenges

The study of complex systems presents a range of difficult challenges by itself (San
Miguel et al. 2012). The complexity of a system is an emergent property that arises from
various factors, including its intricate structure, the number of elements it contains, how
it functions, and how it responds to different types of external influences. These factors
contribute to the overall complexity of the system, which could be measured in various
ways. For many complex systems, their emergent mechanisms are poorly understood.
We identify four main challenges associated with the study of complex systems in the
context of this thesis: (i) the questions of the design choices in defining and sampling a
complex system (section 1.3.1), (ii) which will in turn define their potential to support
a form of open-ended evolution (section 1.3.3) and (iii) how we can expect to build an
interface to communicate with it (section 1.3.4), (iv) which is essential to achieve some
form of control of that system (section 1.3.2).
<figure><img src="image_12.png" /><figcaption>**Figure 1.4** – The challenges of working with complex systems encountered across this thesis can be broken down into several categories. (1) The choice and design of the system, (2) the control of the system which involves understanding the complex mapping between a parameter space and the possibly unpredictable behavior of a system, (3) the search for open-ended evolution properties similar to the ones found in nature, (4) the issue of sending inputs and reading outputs from the internal state of a complex system. This illustration uses *Lenia icon4.png* by Bert Wang-Chak Chan, licensed under CC BY 4.0.</figcaption></figure>
1.3.1 Design of a complex system

The first challenge is to construct a suitable complex system. The definition of complex systems is broad, and several systems with interesting dynamics have been studied under that name. For example, abstract models such as L-Systems (formal grammars used to model the growth and development of organic structures), Random Boolean networks (dynamical systems consisting of a fixed number of binary nodes that are randomly connected to each other and updated based on the state of its neighboring nodes), Turing machines (abstract deterministic computing machines consisting of a tape and a read-write head) and Cellular automata (CAs). In this thesis, we focus mainly on the last listed element: the cellular automaton. There are multiple benefits to working with this model. It is very simple to define, and its high parallelism makes its implementation straightforward. Furthermore, CAs have shown the ability to simulate a wide range of complex behaviors (Wolfram 2002).

Choosing the right complex system architecture is essential because it defines the search space over which interesting and useful systems can be found. Correctly parameterizing that space can also be challenging because too many degrees of freedom make it difficult to search for and find good systems, while too few might indicate a lack of expressivity. An example parameterization of CAs with little expressivity is Langton's lambda parameter (Langton 1990). At the other end of the spectrum, the CA rule is a simple parametrization with many degrees of freedom, which makes it impractical.

Even when we limit ourselves to the study of CAs, many variants can be considered (we list some of them in Chapter 2). They can have continuous or discrete states and operate in continuous or discrete space. For discrete state and space CAs, the number of states and the window of the update function can vary, as well as the topology of the simulation space. We tackle this issue from various angles throughout the thesis, exploring the question of rule definition in Chapter 4, the scale in Chapter 5, and the connectivity pattern in Chapter 6.

1.3.2 Complex systems control

Even if we want complex systems that evolve in unexpected directions and grow in an
open-ended way without supervision, it may still be useful to steer them locally toward
<figure><img src="image_13.png" /><figcaption>**Figure 1.5** – Comparison of the same 1 dimensional CA ran from two initial conditions differing only by one cell. Each row is the state of the CA at one time step, with time increasing from top to bottom. The effect of the small perturbation of initial conditions grows rapidly as the CA evolves in time.</figcaption></figure>

specific targets. A major challenge that follows from this goal is that many complex systems are unpredictable. Deterministic dynamical systems can exhibit behavior that is extremely sensitive to their initial conditions. This sensitivity is sometimes referred to as chaos, where even tiny perturbations in the initial conditions can result in significantly different outcomes or trajectories. This is illustrated with a simple CA rule in Figure 1.5. Because of this property, it is often very hard to predict how a given system will evolve over time, and therefore hard to steer its evolution towards a particular final state or along a chosen trajectory. Although it can be challenging to apply these systems to certain basic tasks, their ability to generate unexpected solutions to difficult problems is also an advantage for the development of open-ended systems. Such systems need to continuously produce new and diverse solutions without a predefined endpoint and their sensitivity to initial conditions can be helpful in achieving this goal.

This problem is also related to the issue of sending a control signal to a complex system that depends on the choice of encoding for that signal. This encoding can significantly impact the system's ability to respond to the control signal and generate the desired outcome, which we describe in Section 1.3.4.

### 1.3.3 Open-ended evolution without objectives

Another challenge is posed by the lack of a clear objective in the design of fully unsupervised learning systems. Even without an objective, it is still essential to understand which complex systems to choose from all available options or how to tune their
parameters so that they behave in interesting ways. For example, CAs such as the one shown in Figure 1.1a is unlikely to be useful for building systems that demonstrate a growth of complexity due to their disordered nature in both space and time. Because of this disordered behavior, they cannot preserve information about the past.

Our goal is to create a system that has the property of evolving in an open-ended way, and we need metrics that can help select systems. However, these metrics should not be used as an additional objective function, as doing so would result in the problems described in Section 1.2. We address this challenge in particular in Chapter 4, by designing a complexity metric that can help select interesting complex systems without relying on a specific task-based performance score.

### 1.3.4 Complex systems inputs and outputs

Some of the complex systems we study in this thesis are closed systems. They do not expect inputs or have well-defined outputs. For example, CA and random boolean network (RBN) do not have the notion of inputs and outputs built into the model. They are standalone objects that evolve according to a set of internal rules.

A common challenge for such systems is to define the inputs and outputs in a way that preserves their internal dynamics. This is also related to the control problem of Section 1.3.2, because controlling a complex system implies being able to send a control signal and read the current state of the system.

Throughout this work, we use a framework called *reservoir computing* (RC) for this purpose. Reservoir computing allows one to harvest the internal computations of complex systems, or *read* information from its internal state by learning a linear regression that maps that internal state to desired outputs (see Section 2.3).

## 1.4 Contributions

The main contributions of this thesis are as follows.

1. Our first contribution is a review of the literature (Chapter 3) describing the connections between cellular automata and other complex systems, open-ended evolution, and neural networks. Studying these fields from a single point of view
is a relatively novel endeavor, and we consider a review necessary to place the
rest of our work in its context.

2. Our second contribution is the development of a general complexity metric that can help identify complex systems with interesting behavior (Chapter 4). Our findings were published in (Cisneros, Sivic, and Mikolov 2019).

3. The third contribution is the development of a coarse-graining method to visualize computations in cellular automata and other discrete systems with local interactions (Chapter 5). Our findings were published in (Cisneros, Sivic, and Mikolov 2020).

4. Our fourth contribution is the introduction of a learning efficiency metric for learning algorithms and a benchmark dataset of progressively harder language tasks (Chapter 6). Our findings were published in (Cisneros, Mikolov, and Sivic 2022).

**1.5 Thesis overview**

Chapter 2 presents some background notions about complex systems, cellular automata, and reservoir computing.

In Chapter 3, we review relevant methods and tools for measuring the complexity of complex systems and use their computations for various tasks.

Chapter 4 introduces a complexity metric that allows us to select complex systems with interesting behavior. The metric measures the “novelty” of the temporal states of a system compared to a reference system. We built a dataset of interesting cellular automata to validate the quality of the metric.

Chapter 5 addresses the question of large-scale complex systems and the applicability of complexity metrics at multiple scales. We propose three algorithms for the coarse-graining of cellular automata. This allows us to reduce the size of large-scale systems while retaining the interesting parts of the behavior.

Chapter 6 presents a learning efficiency metric and a dataset to measure the learning speed of various systems. We show that reservoir computing-based systems using cellular automata can be more efficient than usual machine learning algorithms in constrained
data and computation settings.

Chapter 7 summarizes our contributions to this thesis.

Chapter 8 outlines potential avenues for further research in the field. It provides suggestions on how current work could be expanded or improved, highlighting areas that require additional research to advance understanding of the subject. We propose some experiments that could extend the potential applications of cellular automata, as well as a theoretical framework for learning in dynamical systems.

## 1.6 Publications and software

The thesis has led to the following publications.

- Hugo Cisneros, Josef Sivic, and Tomas Mikolov (Dec. 2019). « Evolving Structures in Complex Systems ». In: *2019 IEEE Symposium Series on Computational Intelligence (SSCI)*. 2019 IEEE Symposium Series on Computational Intelligence (SSCI). Xiamen, China: IEEE, pp. 230–237. ISBN: 978-1-72812-485-8. DOI: 10.1109/SSCI44817.2019.9002840

- Hugo Cisneros, Josef Sivic, and Tomas Mikolov (July 1, 2020). « Visualizing Computation in Large-Scale Cellular Automata ». In: *Artificial Life Conference Proceedings*. ALife 2020. Vol. 32. MIT Press, pp. 239–247. DOI: 10.1162/isal_a_00277

- Hugo Cisneros, Tomas Mikolov, and Josef Sivic (Nov. 28, 2022). « Benchmarking Learning Efficiency in Deep Reservoir Computing ». In: *Proceedings of The 1st Conference on Lifelong Learning Agents*. Conference on Lifelong Learning Agents (CoLLAs 2022). Vol. 199. PMLR, pp. 532–547

The code to reproduce the experiments of all three publications is available on GitHub: https://github.com/hugcis/evolving-structures-in-complex-systems, https://github.com/hugcis/benchmark-learning-efficiency. We also published a dataset that we used for our benchmark in our last publication. It is available at https://github.com/hugcis/incremental_tasks/.
# Background

In this chapter, we provide background on some of the topics that we study in the thesis. We begin by describing the cellular automaton (CA) model, which is the complex system we focus on in most of the thesis (Section 2.1). Next, we describe the similarities between CAs and recurrent neural networks (RNNs), which is useful for understanding the characteristics and applications these two models share (Section 2.2). CAs can be completely reformulated as a special case of RNNs. This reformulation makes it natural to replace recurrent RNNs with CAs in several use cases. We describe the reservoir computing (RC) model (Section 2.3), which enables sending inputs to and harvesting the computations of complex systems to perform machine learning tasks without requiring explicit knowledge of the underlying system dynamics or the need to solve complex optimization problems. We apply reservoir computing to sequential language-like tasks in Chapter 6. This model was developed with RNNs, but we describe how CAs and other complex systems can be a suitable replacement, and the benefits of this approach.

## 2.1 Cellular automata

Stanislaw Ulam and John von Neumann proposed the CA model in the 1940s as a means of modeling crystal growth and creating an autonomous self-replicating system (Von Neumann and Burks 1966).

### 2.1.1 Definition

It is usually defined in a regular lattice in one or two dimensions. Each of its components is called a cell and can be in a state $k \in S$. S is the space of available states for cells, usually chosen to be {0, 1} for binary CAs or {1, ..., n} for CAs with n states.
A neighborhood function $N$ is defined that associates each cell with the indexes of its neighbors in the grid. In general, a CA can be constructed on any space $\mathcal{L}$ where this function can be defined. The space $\mathcal{L}$ specifies an index and the relationship between the cells. In practice, regular finite or infinite grids are chosen, $\mathcal{L} \subset \mathbb{Z}$ or $\mathcal{L} \subset \mathbb{Z}^2$. For example, the grid could be a one-dimensional torus with 10 cells, that is, $\mathcal{L}_{T_{10}} = \{1, 2, \dots, 10\}$. The neighborhood function has the following general form:

$$ \mathcal{N}_{\mathcal{L}}: \mathcal{L} \to \mathcal{L}^s \\ i \mapsto [c_j]_{j \in \mathcal{N}_{c_i}}, \qquad (2.1) $$

where $\mathcal{N}_{c_i}$ is the neighborhood of cell $c_i$, $s$ is the number of cells in the neighborhood, and the returned value is a finite set of cells: the neighbors of cell $c_i$. For the torus $\mathcal{L}_{T_{10}}$ above, we can define the neighbors to be the cell itself and the two immediately adjacent cells. This type of 1D neighborhood is usually associated with a *radius* parameter $r = 1$. The neighborhood $r = 1$ is illustrated in Figure 2.1a. It corresponds to the following neighborhood function:

$$ \mathcal{N}_{\mathcal{L}_{T_{10}}}: \mathcal{L} \to \mathcal{L}^3 \\ \mathcal{N}_{\mathcal{L}_{T_{10}}}(i) = \begin{cases} [c_{i-1}, c_i, c_{i+1}], & \text{if } i \in \{2, \dots, 9\} \\ [c_{10}, c_1, c_2], & \text{if } i = 1 \\ [c_9, c_{10}, c_1], & \text{if } i = 10. \end{cases} \qquad (2.2) $$

A larger radius $r$ would correspond to including more cells in the neighborhood, going in each direction from the initial cell. On two-dimensional grids, there are multiple ways to define the neighborhood. Some common examples are the Moore neighborhood (see Figure 2.1b) and the von Neumann neighborhood (see Figure 2.1c). In the rest, we omit the subscript on the neighborhood function $N$ as we almost always work with regular grids on a torus, for which it is simpler to just mention the size and dimension.

A CA evolves in discrete time steps. An update rule $\Phi: \mathcal{S}^s \to \mathcal{S}$ defines the new state of a cell as a function of its local neighborhood at the current time step. The local group of cells is said to be in a particular configuration. The function $\Phi$ is applied in parallel to all cells. For a CA in its initial state at time step 0 — i.e., a set of cells
<figure><img src="image_14.png" /><figcaption>**Figure 2.1** – Illustration of commonly used neighborhoods for 1D and 2D CA.</figcaption></figure>

$(c_i^{(0)})_{i \in \mathcal{L}} \in \mathcal{S}^{|\mathcal{L}|}$, and a neighborhood function $\mathcal{N}$, we have the following update rule:

$$ \forall i \in \mathcal{L}, \quad c_i^{(t+1)} = \Phi \left( \left( c_j^{(t)} \right)_{j \in \mathcal{N}(c_i)} \right). \qquad (2.3) $$

The details of an update step in a 1-dimensional CA is shown in Figure 2.2 for example. The neighborhood of the cell $c_i$ is $c_i$ itself, as well as its two immediately adjacent cells.

<figure><img src="image_15.png" /><figcaption>**Figure 2.2** – Illustration of a cellular automaton update rule in 1 dimension. For each cell, we look up the neighboring cells and update their state according to the current state of the neighbors. This operation is applied in parallel to the entire grid.</figcaption></figure>

The function $\Phi$ is sometimes called a *rule table* because it associates an output state for each possible combination of input neighbor states. In an implementation of the CA model, these output states can be looked up from a table or array containing all possible transitions from input to output.
There are many ways to define CAs, and we gave a broad definition that includes many rarely used CA variants. For another more straightforward definition of CAs, see, for example, (Kari 2012).

**Rule representation.** A useful rule representation can be obtained by listing all output states corresponding to input neighborhood configurations in a predetermined order. This results in a list of values [$o_1, \dots, o_s|S$], with $\forall i, o_i \in S$, where $s$ is the number of cells in a neighborhood, $S$ is the space of available states, and $|S|$ is the number of available states per cell. Using $o_i$ as the digits of a base-$|S|$ number, each rule is uniquely represented by a number. For example, as explained in more detail in Section 2.1.3, a one-dimensional cellular automaton with a neighborhood of size three has eight possible configurations of neighboring cells (000, 001, 010, 011, 100, 101, 110, and 111), each of which can map to two states (0 or 1). To obtain the rule number, we concatenate the new states for each of the input configurations, resulting in an 8-bit number. The 256 possible binary rules in one dimension with neighborhood size 3 can be numbered from 0 to 255 with this standard method and are referred to by their number in the literature.

**Boundary conditions.** The grid of a CA can be finite or infinite. In the infinite case, the grid is assumed to be initialized to a uniform state, except for a few cells set to other states. The simulation is then run on these few cells, while the rest of the infinite grid does not have to be simulated from the start. For a finite grid, an exhaustive simulation can be run, but one needs to define boundary conditions. The boundaries can be set to wrap around the other side of the grid, forming a torus. An example of the corresponding indexing is given in equation 2.2. Other choices of boundary conditions consist of adding virtual padding cells outside of the main grid. They can be set to a fixed state, a randomly chosen state, or they can mirror the cells inside the grid. Each of these choices affects the evolution and properties of the CA, but the importance of these boundaries decreases for very large grids.

## 2.1.2 Classification of cellular automata

Stephen Wolfram approached CAs as discrete, spatially extended dynamical systems (Wolfram 1984). The analogy is only superficial since many concepts from dynamical
systems theory, such as “chaos”, “attractors” and “sensitivity to initial conditions” only
admit a rigorous definition in the continuous state and continuous time models. Wolfram
proposed a qualitative classification of CA update rule behaviors roughly analogous to
classifications in dynamical systems theory, with four classes defined as follows.

**Class 1** All initial configurations converge to a single fixed configuration, such as a configuration consisting entirely of 1s.

**Class 2** All initial configurations relax after a transient period to some fixed point or some temporally periodic cycle of configurations, but which one depends on the initial configuration. (CA defined on finite lattices always have periodic behavior because there is only a finite number of grid configurations. Class 2 does not refer to this type of periodic behavior but rather to cycles with periods much shorter than the total number of possible states).

**Class 3** All initial configurations exhibit chaotic behavior after a short transient period. (The term “chaotic” here refers to apparently unpredictable space-time behavior, not the standard notion of chaos in dynamical system theory.)

**Class 4** Some initial configurations result in complex localized structures that can persist for a long time.

This classification being qualitative and not rigorous, there is not even a clear con-
sensus on which of the Elementary cellular automata (ECAs) rules (one dimensional
binary CAs with neighborhood size 3) belong to class 4 or class 3. Class 4 CA rules
are speculated to be capable of universal computation (Wolfram 1984), which refers to
the ability of a system to simulate any other system or computation, given enough time
and memory. For example, Li and Norman H. Packard (1990) claimed that ECA rule
110 has class 4 behavior and was eventually proven to be universal (Cook 2004), but no
other ECA has been proven universal since.

Zenil (2010) studied the compression size of the space-time diagrams of all ECA for
a fixed-length simulation. Using a simple k-means clustering technique, he obtained two
groups that roughly match Wolfram's classes 1 and 2 and classes 3 and 4. We repro-
duce these results in Figure 4.1b in Chapter 4. They offer an interesting nonqualitative
confirmation of the results of Wolfram's classification. However, these results vary signif-
icantly if we modify the initial conditions, as well as the grid size, data representation, or compression algorithm (Hudcová and Mikolov 2020).

Wuensche and Lesser (1992) studied the behavior of ECAs when the simulations are reversed and computed the preimages of each configuration. They introduced the Z-parameter, which is the probability that a partial preimage can be extended by one symbol for each CA. The authors speculate that class 4 occurs at $Z \approx 0.75$. This is not a classification but a class 4 membership test that can be computed from the rule directly, making it practical compared to alternatives. However, when tested in practice, Wuensche's hypothesis that $Z \approx 0.75$ corresponds to class 4 is rarely verified.

Hudcová defined a classification of CA that is based on the form of the asymptotic growth of its transients, obtained by repeated simulation of the rules of CA starting from various initial conditions and using the best fitting function for the asymptotic growth of its transients (Hudcová and Mikolov 2020, 2022). For a detailed review of the current methods and challenges of CAs classifications, see (Vispoel, Daly, and Baetens 2022).

### 2.1.3 Cellular automata variants

Several CA variants have been proposed, modifying or restricting various parts of the definition in Section 2.1.1. Here, we list some common ones.

**Elementary cellular automata.** ECAs are 1 dimensional CAs with two states per cell and neighborhood size 3 — the cell and its two direct neighbors on the 1D grid. There are 8 possible configurations of a neighborhood with 3 cells and 2 states per cell, which corresponds to 256 possible ways to define an ECA — two possible outputs for each of these 8 possible configurations; hence $2^8 = 256$ possible CA. This relatively small number of rules allows for an exhaustive exploration of the rule space and mapping of the properties of ECA, which would not be possible for general CAs.

<figure><img src="image_16.png" /><figcaption>**Figure 2.3** – Illustration of the rule of ECA number 150.</figcaption></figure>

These CAs have been extensively studied and offer an interesting combination of
trivial definition and implementation, and complex and unpredictable properties. One
of the fundamental problems of CA research is to classify the 256 rules into well-defined
behavior types and order them by complexity, which was attempted in several previous
works (H. A. Gutowitz 1991; Hudcová and Mikolov 2020, 2021; Wolfram 2002; Wuensche
1999; Wuensche and Lesser 1992; Zenil 2010). We discuss these classifications in more
detail in Section 2.1.2.

<figure><img src="image_17.png" /><figcaption>**Figure 2.4** — Evolution of ECA number 150 simulated on a grid of size 100 for 50 steps. 2.4a shows a simulation starting from a blank grid with only one cell set to 1. 2.4b shows a simulation starting from a random initialization. Time flows from top to bottom, and each row of cells represents a single internal CA state at a single time step.</figcaption></figure>

ECA rules are easily visualized because there are only eight possible neighborhood
configurations that need to be shown. For example, Figure 2.3 shows the transition table
for ECA Rule 150. If we assign 1 to the black state and 0 to the white state, the possible
neighborhood configurations are ordered in their binary order from right to left, with
the leftmost bit being the most significant. This is how the rule number is computed,
using the output states (last row in Figure 2.3) as a binary representation. Example
simulations of ECA rule 150 starting from various initial states are shown in Figure 2.4.
In these representations, time flows from top to bottom, and each row of cells represents
a single internal CA state at a single time step. This representation of an ECA rule is
often called a spacetime diagram.

There are several well-known ECA rules with particularly interesting properties. For
example, a universal computer has been constructed in rule 110 (Cook 2004), making it
the first (and only to date, although ECA Rule 54 is postulated to be as well) Turing
complete ECA.
**Totalistic cellular automata.** Totalistic CA are a subset of CA whose rules can be expressed as a function of the sum of neighboring cell values. These CA were introduced by Stephen Wolfram (Wolfram 1983). ECA 150 is an example of a totalistic CA, and its rule and typical space-time diagram are depicted in Figures 2.3 and 2.4. Its rule could be summarized in natural language as "if exactly two states are 1 or all states are 0, the next state is 0. Otherwise, if one or three states are 1, the next state is 1". Conway's Game of Life, presented in the following paragraph, is an example of a totalistic CA in two dimensions.

**Game of Life.** The game of life is one of the most famous CA. It was proposed by the mathematician John Conway in 1970 (Gardner 1970). It is a two-dimensional binary CA. Its two states are often called “alive” and “dead”. Its rule can be summarized in three sentences as follows.

- Any live cell with two or three live neighbors survives.

- Any dead cell with three living neighbors becomes a live cell.

- Other live cells and already dead cells are dead in the next generation.

This CA has a particularly active community dedicated to finding interesting patterns with particular properties, such as long cycling periods or particular speeds of movement through the grid. Some of these patterns have been used as memory registers and communication channels to build a Universal Computer (Igblan - Life Universal Computer 2022). This construction proved that the game of life is Turing complete. This property is thought to be also shared by other CAs — it is proven for at least ECA rule 110 (Cook 2004). This makes CAs models theoretically appealing for the design of a learning algorithm because they can simulate any algorithm.

**Asynchronous cellular automata.** A cellular automaton is said to be asynchronous when its cells are not updated in parallel at each time step. Various cell update schemes can be chosen. Asynchronous CAs have to define an order of update of the cells, which is not easily defined on an infinite grid. Therefore, most implementations are done on finite grids. This update order can be random, follow a predefined order, or be decided by a global controller. Groups of cells can be updated simultaneously, or updates can be done cell by cell. The many effects of asynchronous updates can be difficult to predict
<figure><img src="image_18.png" /><figcaption>**Figure 2.5** – Some game of life patterns with various properties. The moving oscillator 2.5a moves 1 cell along the bottom right diagonal every 4 steps. The fixed pattern 2.5b never changes except if it interacts with others. The oscillator 2.5c loops through three different configurations.</figcaption></figure>

and can make these types of CAs relatively difficult to work with.

Asynchronous CAs are an attractive model when using large systems for which parallel update is prohibitively expensive. The rule could adapt itself and preferentially update cells in active or useful parts of the CA while doing slower updates to the less active parts of the CA. An asynchronous CA with evolutionary properties was constructed in (Nehaniv 2003). The authors argue that asynchronicity is a strong advantage in the construction of systems that can evolve within CAs.

**Stochastic cellular automata.** In a stochastic cellular automaton, the update function $\Phi$ is stochastic. This means that the next state of all cells in the grid is sampled from a probability distribution that depends on the current neighborhood configuration. We can write

$$ \forall i \in \mathcal{L}, s \in \mathcal{S} \quad p(c_i^{(t+1)} = s) = \Phi(N(c_i^{(t)}))(s). \qquad (2.4) $$

These CAs have had successful applications in physics (Ottavi and Parodi 1989; Vichniac 1984) or biology (Boas et al. 2018), because their stochasticity allows modeling complex physical phenomena.

**Continuous cellular automata.** Another family of CAs uses real-valued states, often restricted to the range [0, 1]. The update function is then a real multivariate function,
which we can write

$$ \Phi : [0, 1]^s \to [0, 1]. \qquad (2.5) $$

Examples of continuous CAs include Lenia, which uses convolutional operators followed by thresholding for the update rule (Chan 2019), or neural cellular automata (NCAs) (Mordvintsev et al. 2020) which represent CAs as neural networks (more details in Sections 2.2 and 3.5.1). Garzon and Botelho (1993) explores the condition for a real-valued CA to be able to compute real-valued functions.

**Higher dimensional cellular automata.** CAs of dimension greater than 2 have been comparatively less studied for several reasons, including the limits that arise when simulating and visualizing systems in more than 2 dimensions on a computer screen. Another problem is with the increasing number of possible rules in higher dimensions. The number of neighbors per cell grows exponentially with dimension, and the number of possible rules has a doubly exponential growth rate. This makes the convenient representation of CAs rules as tables infeasible since the size of that table quickly exceeds the available memory of most computers. For example, there are $2^9 = 512$ possible $3 \times 3$ configurations of a binary Moore neighborhood (see Figure 2.1b) in 2D that can each lead to 2 states, which means that there exists $2^{512} \approx 10^{154}$ distinct 2D rules with 2 states and a Moore neighborhood. This is already an incomprehensibly large number, but there are $2^{134,217,728} \approx 10^{40,403,562}$ such rules in 3D. In general, there are $2^{3n}$ possible binary CA rules with a Moore-like neighborhood in dimension $n$.

Despite these limitations, several works have studied higher dimensions CAs, although mostly in 3 dimensions (Sudhakaran et al. 2021; Tsalides, Hicks, and York 1989). For example, there are several examples of successful 3D CA simulations applied to material sciences (Arata et al. 1999; Di Caprio et al. 2016; Gandin and Rappaz 1997; Pan, Feng, and Hudson 2009).

**Hexagonal cellular automata.** Hexagonal CAs are defined on grids tiled with hexagons, where each cell has 6 direct neighbors.

One advantage of this type of CAs over regular grids is the absence of preferred direction since they are all equivalent. In a square grid, diagonally adjacent cells span
<figure><img src="image_19.png" /><figcaption>**Figure 2.6** – Example definition of the neighborhood for hexagonal cellular automata.</figcaption></figure>

more distance than directly adjacent cells, whereas on a hexagonal grid, all neighbors are equivalent (Moree and Riele 2004).

### 2.1.4 Parametrizing and sampling CA rules

As explained in Section 1.3, one of the main challenges of working with complex systems is the parametrization and control of the behavior of that system. CAs are no exception, and several parametrization schemes have been proposed, with their respective benefits and drawbacks. There are infinitely many ways to define rules for CAs, and for each of these definitions, a very large number of rules are available, making it impossible to sample a significant portion of the rule space. A good parametrization allows one to scan a wide range of behaviors of CAs by traversing the space in an interesting way.

**Langton's lambda.** Christopher Langton proposed one of the first parameters that could be used to control the behavior of CAs to a certain extent (Langton 1986, 1990). For a CA with $K = |S|$ states and a state arbitrarily chosen to be the *quiescent* state (this corresponds to the “dead” state in the game of life, for example), if there are $n$ transitions to the quiescent state within the rule table and $K^s - n$ remaining transitions that do not lead to the quiescent state, where $s$ is the number of cells in the neighborhood, we have

$$ \lambda = \frac{K^s - n}{K^s}. \qquad (2.6) $$
For Langton, the purpose of $\lambda$ is to search the space of CAs in an ordered manner by varying the value of the parameter. His procedure starts with a rule with all transitions leading to the quiescent state. The value of $\lambda$ is increased in discrete steps to 1 – $1/K^s$ by randomly changing the transitions of the rule to lead to a different state. An illustration of the behavior change of a CA rule under this procedure is shown in Figure 2.7. Langton collected various measures of the dynamical behavior of CA and studied them as a function of $\lambda$. He observes a phase transition as $\lambda$ approaches its maximal value, where CAs seem to behave in the most complex way, with long transients and large-scale propagating structures. He calls this phase transition “edge of chaos” because it corresponds to $\lambda$ values just before the generated CAs become chaotic.

<figure><img src="image_20.png" /><figcaption>**Figure 2.7** – The effect of varying $\lambda$ for a 1D binary CA rule. A CA rule is progressively modified along a trajectory of increasing $\lambda$. The rule starts with a fixed behavior, becomes periodic, and then complex when hitting the “edge of chaos”. Note the localized structures that spread over time (down along the y axis) and interact in the fourth figure. When $\lambda$ increases further the rule becomes chaotic.</figcaption></figure>

**Dirichlet sampling.** For an arbitrary number of states and an arbitrary neighborhood size, it can be challenging to sample CAs while scanning the wide range of behavior these models are capable of simulating. This is because the size of the space of CA rules is very large. The naive method of uniformly sampling each transition output is not useful in large rule spaces with multiple states, large neighborhoods, or grid dimensions larger than 1.

Rules with equal proportions of transitions leading to all states tend to be the most chaotic. This is also what Langton observed when studying the $\lambda$ parameter. This is explained by the fact that sampling each transition output uniformly is equivalent to sampling the rules from a multinomial distribution with the number of possible states ($K$), number of possible transitions ($n = K^{K^s}$, with $s$ the number of cells in a neighborhood), and the uniform probabilities $p_1 = \dots = p_K$ as parameters. A random
variable $X = (X_1, \dots, X_K)$ indicates the number of times each result is observed. The probability mass function of this distribution is

$$P(X_1 = x_1, \dots, X_K = x_K) = \frac{n!}{x_1! \cdots x_K!} p_1^{x_1} \cdots p_K^{x_K} \quad (2.7)$$

This quantity becomes vanishingly small for large $n$ (that is, rules with many states or large neighborhoods) and rules with an output transition skewed towards a specific state. For example, the binary Game of Life rule has 372 transitions leading to state 0 and 140 leading to state 1. The probability of sampling a rule with these transition proportions is $8.24e-26$ (compared to a 17% probability of sampling a rule with between 254 and 257 transitions to 0). Using a uniform transition sampling method makes it highly improbable to obtain a rule resembling the Game of Life, let alone the rule itself.

We propose a Dirichlet sampling method as an alternative way to sample rules with fixed ratios of transitions leading to each output state. For a rule with $K$ states, we sample a $K$-uple from a Dirichlet distribution of order $K$ with parameters $(\alpha_k)_{k \in [1, K]}$ where $\alpha_0 = \alpha_1 = \dots = \alpha_K = \alpha < 1$. The result is a quantile $K$-uple $(q_1, \dots, q_K)$. For $\alpha < 1$, the distribution is concentrated around the corners of a simplex of dimension $K$, which means that it preferably samples a tuple with one of its values dominating the others.

Rules are sampled so that the number of transitions to each output state matches the quantile generated from the Dirichlet distribution. The samples will be more likely to have a dominant quantile, which can be associated with the quiescent state in Langton's $\lambda$ calculation. The resulting sampling of cellular automata is much better for scanning the entire space of rules and generating CA rules that would be impossible to reach with the naive sampling method, as illustrated in Figure 2.8.

**Smooth sampling with the recurrent convolutional neural network analogy.**

The recurrent convolutional network analogy of CAs (see Section 2.2 for details) formulates a CA as a neural network. In this paradigm, the rule of a CA is completely determined by the parameters of that neural network. This allows us to modify the CA rules as we would change the weights of a neural network. For example, as long as all simulation steps are differentiable, the backpropagation algorithm can be applied and
<figure><img src="image_21.png" /><figcaption>**Figure 2.8** – An illustration of the difference between naive uniform rule sampling and Dirichlet sampling on binary rules. 10000 2D binary CA rules were sampled with each method, the plot displays a histogram of their average output state. Uniform sampling leads to more rules with an average output state close to 0.5, meaning approximately as many transitions toward both states. Dirichlet sampling can control, through the choice of parameter α, how many rules with more skewed transitions should be sampled, which is shown by the spikes on the blue histogram close to 0 and 1.</figcaption></figure>

used to update a CA rule by updating the weights of the neural network with gradient descent. Dimensionality reduction methods applied to the parameter space of the neural network can help to understand the structure of that space and the CAs behaviors produced by different sets of parameters.

## 2.2 Cellular automata and Recurrent neural networks

The purpose of this section is to show that cellular automata and recurrent convolutional neural networks have very strong connections and to draw this parallel as clearly as possible. We express the CA update function as a set of convolutional operations that can be derived from any CA rule.

This connection yields interesting consequences for both the theoretical properties of these models and the potential applications of CAs and recurrent networks. It shows, for instance, that emergent properties with increasing complexity and perhaps open-ended development that we expect from complex systems such as CAs are possible within the hidden state of a recurrent neural network.

The conceptual similarity between CAs and neural networks is not surprising. Fundamentally, CAs represent the paradigm of emergence, in which the origins of complex
behaviors are searched for in an assemblage of possibly very simple parts rather than
being viewed as a sum of complex building blocks. It is often argued that there is no
better-known example of a truly emergent phenomenon than that of the emergence of
consciousness out of the large network of functionally simple (and certainly unconscious)
neuronal components that make up the human brain. A biological neural network con-
sists of a large space of interconnected nodes with dynamic behavior that is a local
function of the other nodes to which it is connected, which is strikingly similar to the
principle of a CA. Artificial neural networks can be thought of as being a set of bio-
logically inspired CA rules (Ilachinski 2001). The more precise parallel between cellular
automata and a form of recurrent/convolutional network has also been drawn by several
other researchers (Gilpin 2018; Mordvintsev et al. 2020; Wulff and Hertz 1993).

In this section, we formally describe the relation between CAs and RNNs and discuss
the extensions this implies for the CA model, as well as some of the consequences of
that relation.

2.2.1 RNN formalism

<figure><img src="image_22.png" /><figcaption>**Figure 2.9** – Standard RNN architecture (top) and Game of Life seen as a RNN with no inputs and outputs (bottom). For the CA, the “hidden” state $h_t$ is the current state of the grid at timestep $t$. The operator $\Phi$ is the CA update rule which is equivalent to a convolutional neural network (CNN) (see Section 2.2.2). This illustration is based on *Recurrent neural network unfold* by fdeloche, licensed under *CC BY 4.0*.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>**Figure 2.10** – The CA update rule is local and can be represented by a set of CNNs: linear local uniform transformations followed by the application of a non-linear function.</figcaption></figure>

We write the definition of a CA with a formalism and notation inspired by RNNs. This parallel is illustrated in Figures 2.9 and 2.10. The state of the grid at time *t* is denoted *h*<sub>*t*</sub> and corresponds to the hidden state in a RNN. In the case of classical CAs, it is a 1 or 2D vector of discrete values (see Section 2.1.1), but (Mordvintsev et al. 2020) and other CAs extensions use continuous values, much like usual RNNs.

The mapping $\Phi$ operates only on the hidden state. Because it is local, it can be shown to be equivalent to a convolutional layer, as explained below in Section 2.2.2. The inputs and outputs $(\mathbf{x}, \mathbf{o})$ in Figure 2.9 are not included in the classical definition of CAs, but are an easy-to-implement extension, as discussed in Section 2.2.3.

2.2.2 Transition rule as a set of convolutions

Each cell $c_i$ from the CA definition can be in one of the $k$ states. We represent the cells as vectors of size $k$, where $k$ is the number of states, with only zeros excepted for a one in the $k$-th position, which can be interpreted as a one-hot encoding of the cell state. A neighborhood of size 3 in a 1D CA can be represented as a $3 \times k$ vector $\mathbf{u}_i = [u_{i-1}, u_i, u_{i+1}]$. Each $u_i$ is a vector of size $k$ with a 1 in position $s_i$. This is illustrated in the 1D (resp. 2D) case on the left (resp. right) of Figure 2.11. For a CA with only two states, it is redundant to have a $3 \times 2$ vector, but the “one-hot” encoding becomes useful when working with more states.

With this representation, we can express the transition rule $\Phi$ as a simple convolutional neural network. This network would be made up of two layers, which are shown in Figure 2.12 for the 2D case. We describe a simple construction of a CNN representing any binary CA rule below:
<figure><img src="image_24.png" /><figcaption>**Figure 2.11** – A CA neighborhood vector representation example with 2 states. A 3 × 3 square of cells with two states can be represented by a 3 × 3 × 2 tensor. Left: 1D 3-neighbors representation. Right: 2D 3 × 3 neighbors representation. The last dimension of the vector encodes a “one-hot” representation of the state of the cell.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(a) A representation of any binary 2D CA rule operating on 3 × 3 neighborhoods as a 2-layer CNN. There are 512 filters receptive to each 3 × 3 neighborhood states in the first layer. The second layer determines the output value for each neighborhood state.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) A representation of any binary 1D CA rule operating on neighborhoods of size 3 as a 2-layer CNN. There are 8 filters of size 3 × 3 in the first convolutional layer, and 8 filters of size 1 × 1 in the second one.</figcaption></figure>

<caption>**Figure 2.12** – Two examples of CA rules represented as a simple CNN. Padding ensures the grid obtained after the CNN step is of the same size as the original grid.</caption>
**The first convolutional layer.** It is receptive to each possible neighborhood configuration. In the 1D case with neighborhood size 3, this is all eight configurations 000,001,010, ..., 111. For a radius $r$ in 1D, this layer is composed of $k^{2r+1}$ filters of dimension $k \times (2r + 1)$ with only ones and zeros, which are the mirror of a possible configuration of size $(2r + 1)$ of the neighborhood.

The product of each filter with an input neighborhood in the grid will be an integer between 0 and $2r + 1$. Applying a ReLU nonlinearity preceded by a constant vector of value $2r$ as bias, we obtain the input to the second layer, a vector of size $k^{2r+1}$ for each cell containing all zeros except for one corresponding to the detected neighborhood. This vector plays the role of an indicator of the current input configuration.

**A second convolutional layer.** It has $k^{2r+1}$ filters of size 1 that are 0s or 1s depending on the desired output value of that transition. When applied to the indicator vector above, the output obtained will be the output value of the rule corresponding to the detected neighborhood state.

With some circular padding that wraps around the edges of the grid or applies constant values outside of the main grid, the size of the output is the same as the input. This means that the CNN step can be applied multiple times, simulating a single CA step every time.

### 2.2.3 Adding inputs and outputs

Figure 2.9 shows the locations where inputs and outputs can be incorporated into the Game of Life or any other CA. We provide a variety of examples of how input and output can be integrated and how they can interact with the rule of the CA.

#### Inputs

There are many ways to add inputs to the recurrent CNN model presented above. We divide our proposed methods into two classes: inputs that directly modify the hidden state during the CA evolution and inputs that augment the hidden state without changing it directly and affect the rule.
**Hidden-state augmentation and rule modulation.** Inputs are usually added to the hidden state just before the application of the non-linearity in standard RNNs.

In the case of cellular automata, it can also be done in the following ways:

<figure><img src="image_27.png" /><figcaption>**Figure 2.13** — Inputs $x$ and outputs $y$ can be directly encoded within the hidden state. $U$ and $V$ are the operators mapping $x$ to a hidden state tensor and decoding a hidden state tensor to output $y$.</figcaption></figure>

- **Vector state.** A straightforward way to add inputs to a CA is to consider the grid as having an additional (possibly read-only) dimension. For example, a 1D automaton of size $N$ would be represented by a vector of dimension $N \times 2$ where the first vector of size $N$ is the state of the grid and the other is the input. The update rule would be changed to take into account this new component. One can see this model as using multiple CA rules at the same time, with the input state conditioning the rule chosen for a given update step.

- **Variable update rule.** Inputs can also influence a CA's evolution by changing the update rule. With this configuration, the update rule $\Phi$ is now a function of the input $x$. We now have $\Phi_x = G(x)$, where

$$G: \mathcal{X} \to (\{1, \dots, k\}^{2r+1} \to \{1, \dots, k\}),$$

and $\mathcal{X}$ is the input space.

A similar approach was used in (Adams et al. 2017). It showed that conditioning a CA update rule on another CA's evolution could enable a higher diversity of behavior.

- **Concatenation.** One way to incorporate input is by concatenating it with the hidden state, such as at the boundaries of the grid, before applying the update rule. This creates a fixed read-only memory space, as shown in Figure 2.14. However, a potential drawback of this method is that information must propagate through the grid in order to be processed in locations far from the read-only memory's position.
This approach is more reminiscent of computer architecture, where designated areas are reserved for specific data functions.

<figure><img src="image_28.png" /><figcaption>**Figure 2.14** — Input $x_t$ is projected and concatenated to hidden state $h_t$, affecting the boundary conditions of the CA. For the CA, this vector is just like another part of the hidden state except that it is not writable</figcaption></figure>

**Hidden-state manipulation.** Another approach is to directly modify the hidden state to “communicate” information to the system. These methods are used in (Mordvintsev et al. 2020; Randazzo et al. 2020) to make interactive demos and allow users to directly modify that hidden state. A user can interact with the system by drawing with its mouse. This sets parts of the internal state to a fixed cell state.

— **Masking.** Input data can serve as a mask on top of the current grid state, forcing some cells into states that depend on the input values. If we view the state activations in the CA as some neural pattern, this approach can be seen as a form of *neuromodulation*, which has previously been used in machine learning (Beaulieu et al. 2020; Ishiguro, Fujii, and Hotz 2003; Soltoggio et al. 2008).

A mask can be constructed from an input vector $x$ by linearly transforming $x$ and applying an activation function (e.g., step function). It controls which channels are activated and where information can flow or not. This kind of mask can be applied with element-wise multiplication or addition but also more elaborate operations. It is illustrated in Figure 2.15.

## Outputs

Similarly, inputs can be extracted from the hidden state in several different ways. One might read outputs from the boundaries of the grid, from a transformation of the
<figure><img src="image_29.png" /><figcaption>**Figure 2.15** — Input is converted into a binary mask for the hidden state. The mask is applied through element-wise multiplication here.</figcaption></figure>

grid using a neural network, etc. The output could also be stored within an additional dimension of the “extended” grid state presented above.

## Cellular automata and computations

From a computational point of view, the hidden state can be seen as a working tape on which some kind of *parallel Turing machine* (the cellular automaton) performs computations. In this framework, the input is encoded as the initial state of the tape and decoded from the last state of the tape (as illustrated in Figure 2.13). The recurrent convolutional neural network (RCNN), or CA, then plays the role of a fixed computer program that is executed. This is the setting adopted in previous work on constructing computations with cellular automata: a task is chosen, and rules that can execute that task on input/output pairs are searched (Melanie Mitchell 2005).

Another interesting point of view can be taken when we consider the hidden state (or tape) as encoding both some data and a computer program, as is the case with a Turing machine. We can then expect a CA (or RNN) to not only compute the result of a particular function but to be equivalent to a general-purpose computer capable of computing the result of any chosen algorithm without any need for optimization. Making a universal CA compute some program amounts to encoding the program in the right language and communicating it by embedding it within the hidden state. In practice, figuring out the right encoding is often the biggest challenge. This is also the reason that the few Turing complete CAs are not so useful in practice since we do not know of
any efficient encoding for programming them.

### 2.2.4 Consequences

Viewing cellular automata as recurrent convolutional neural networks, as described above, has several interesting consequences, of which we list a few here.

#### Turing-completeness of the system

Because we can simulate rule 110 ECA in the above RCNN system, it follows from the Turing completeness of this CA rule that the RCNN is Turing complete. This is an interesting result, although not very significant since RNNs have already been proven to be Turing complete (Siegelmann and Sontag 1992), with very few practical implications. Our model is relatively far from a real-world CNN with a fixed number of layers independent from one another — compared to a variable number of steps and shared layers for the automaton-RCNN. Thus, this property does not seem to bear any consequence on CNNs.

#### Differentiable cellular automata

Each step of the computations involved in computing a step of a cellular automaton represented as a RCNN is differentiable. Therefore, we can theoretically couple this framework with backpropagation to create *learnable* cellular automata that can adapt their rules to minimize a target loss function.

This is the direction taken by Mordvintsev et al. (2020). The authors use supervised learning on a CA and train it to obtain a stable self-repairing target shape. However, because supervised CA can only do as much as they have been trained to do, it may defeat the purpose of working with a model capable of spontaneous complex emergent behavior. An open-ended complexity increase could very unlikely be achieved through pure supervision.

Gilpin (2018) takes the reverse approach and tries to train RCNNs to simulate a fixed CA rule, using the statistics of that training process as a way to help understand the structure of the CA rule space.
### 2.2.5 Beyond the naive rule representation

The representation of CA rules that we built in Section 2.2.2 is a one-to-one mapping. Each rule has its RCNN counterpart, and vice versa. However, there are several ways to make that representation of CA rules more efficient, using a neural network with fewer parameters. For example, we can make the first layer of filters receptive to both a given configuration and its inverse (e.g., (1, 0, 1) and (0, 1, 0) in an elementary cellular automaton) by using negative values in the filter (use a filter (1, -1, 1) instead of (1, 0, 1) and (0, 1, 0) separately).

For example, the rule of the game of life does not need more than two convolutional filters on the first layer to be represented by a CNN because it is totalistic. The new state of a cell depends only on its number of neighbors in the alive state and the current state of the cell.

## 2.3 Reservoir computing

Reservoir computing (RC) is a computational framework that aims to exploit the states of a complex dynamical system to perform some target task (G. Tanaka et al. 2019). It relies on a reservoir of computations, that is, a dynamical system performing some computations on its internal state. An input signal is fed into this reservoir that behaves as a black box from the point of view of the algorithm. Input values can be sequential or nonsequential. They are projected and combined with the dynamical system using a suitable mechanism. This projection can be learned, set randomly, or chosen. The internal state of the dynamical system evolves according to its update rule, which is defined by the chosen system. Finally, a decoder model (usually a linear regression) is trained in a supervised manner to extract necessary information from the internal state of the dynamical systems to predict the right output. A general diagram describing this process is shown in figure 2.16.

Some early known formulations of the RC idea were done by Kirby (1991) and Schomaker (1990, 1991, 1992). In another early work, Buonomano and Merzenich (1995) used a random spiking neural network with excitatory and inhibitory elements. The probabilities of neuron connection was inspired by the real rat brain data (Mason, Nicoll,
<figure><img src="image_30.png" /><figcaption>**Figure 2.16** – Diagram of the general functioning of a reservoir computing (RC) system. Input values are projected in the state or evolution function of a dynamical system which is run according to its internal rule. Output values are decoded from the internal state of the dynamical systems with a trainable layer.</figcaption></figure>

and Stratford 1991). The authors observed an interval-sensitive response of these neurons when subjected to pulses spaced by varying time intervals. This indicates an ability to encode temporal information. To measure this phenomenon, they trained an output linear layer to recognize specific patterns in the activation of the last layer of neurons. This idea of using a randomly initialized RNN with fixed weights and training a simple linear layer to decode outputs was reinvented and popularized later under the names “Echo-state networks” (Jaeger 2001) and “Liquid state machines” (Maass, Natschläger, and Markram 2002).

Another advantage of reservoir computing is that reservoirs can be implemented using a variety of physical systems, substrates, and devices. These types of reservoir are sometimes called “exotic” (Lukoševičius and Jaeger 2009) in contrast to RNN-based reservoirs. Physical RC is a promising area of research that may help us develop cheap and efficient machine-learning hardware and can inform us about the nature and behavior of the underlying physical systems. G. Tanaka et al. (2019) identifies four characteristics that make a reservoir suitable for solving a computational task.

1. The inputs must be mapped to a high-dimensional space, which corresponds to the internal state of the reservoir. This high-dimensionality ensures that the inputs are separated with high probability.

2. The reservoir should be non-linear. This will allow inputs that are not linearly separable to be separated as they are projected within the reservoir.

3. The reservoir should have the echo-state property (Jaeger 2001; Yildiz, Jaeger, and Kiebel 2012), or fading memory (Boyd and Chua 1985; Maass, Natschläger,
and Markram 2002, 2004). This property relates the asymptotic properties of
the excited reservoir dynamics to the driving signal. A reservoir with the echo-
state property will asymptotically forget all input values. This makes some
reservoirs unable to solve tasks where unbounded memory is possibly required,
such as context-free grammar parsing (Schmidhuber, Wierstra, et al. 2007). Later,
reservoir models that address this issue were developed (Pascanu and Jaeger 2011).

4. A reservoir should separate responses from different signals into different parts of the space while being insensitive to small perturbations of the input. For a given system, this trade-off is often obtained at the transition between chaotic and non-chaotic regimes (Bertschinger and Natschläger 2004; Legenstein and Maass 2007).

We list common types of RC that are based on various families of reservoirs below:

**Dynamical systems.** This refers to a broad category of RC models based on other dynamical systems than the random RNNs. This includes reservoir computing with Cellular automata, which we give more detail about in section 2.3.3.

**Electronic RC.** These models are implemented on electronic circuits, such as artificial neural networks implemented on physical circuits or other neuromorphic circuits. For example, RC models can be built on FPGAs (Alomar et al. 2014; Antonik 2018; Antonik, Smerieri, et al. 2015; David Verstraeten, Benjamin Schrauwen, and Stroobandt 2005) or memristive circuits (Donahue et al. 2015; Merkel et al. 2014; Yang, Chen, and F. Z. Wang 2016).

**Mechanical RC.** Soft mechanical bodies have complex nonlinear dynamics that can be used as a reservoir to build a RC model (Pfeifer, Bongard, and Grand 2007). For example, a reservoir using a mass-spring network was proposed by Hauser et al. (2011).

**Biological RC.** There have been multiple hypotheses on whether the parts of the brain behave like a RC system (Yamazaki and S. Tanaka 2007). The goal of this area of research is to study the RC-like properties of biological systems rather than build reservoir computers. For example, a RC model was proposed to understand the mechanism of context-dependent eye movement (P. Dominey,
Arbib, and Joseph 1995; P. F. Dominey 1995).

**Photonic RC.** The principle of photonic RC is to use optical reservoirs to generate the computations. A first example using semiconductor optical amplifiers (SOAs) was proposed by Vandoorne and colleague (Vandoorne, Dambre, et al. 2011; Vandoorne, Dierckx, et al. 2008) and subsequently built (Vandoorne, Mechet, et al. 2014). The photonic implementation of reservoir computing is based on the idea that photonic accelerators (Kitayama et al. 2019) can realize fast information processing with low learning costs (Antonik, Marsal, et al. 2019; Martinenghi et al. 2012; Paquot et al. 2012; Sugano, Kanno, and Uchida 2020). Some of the implementations use the scattering of light in complex media (Dong et al. 2020; Rafayelyan et al. 2020).

**Spintronics RC.** Some RC works used nanoscale electronics involving the charge and spin of electrons called *spintronics* (Wolf et al. 2001). Spintronics has been used in multiple ways to build RC models, using spin torque oscillators (STO) (Torrejon et al. 2017; Williame et al. 2019), or spin waves (Nakane, G. Tanaka, and Hirose 2018) for example.

### 2.3.1 Reservoir computing applications

RC has had successful applications in various fields. It has been applied for pattern classification; in audio with spoken digit recognition (D. Verstraeten et al. 2005) and waveform classification (Paquot et al. 2012); in computer vision with written digits recognition (Jalalvand, Van Wallendael, and Van De Walle 2015), and human motion classification (Antonik, Marsal, et al. 2019; Soh and Demiris 2012). Another popular application of RC is time-series forecasting (Jaeger 2001, 2002; Wyffels and B. Schrauwen 2010). RC has also been used in reinforcement learning to control agents in artificial environments (Kanno and Uchida 2022) or learn radio spectrum access (Chang et al. 2019).

### 2.3.2 Echo-state network

The echo-state network (ESN) (and simultaneously developed model *Liquid state machine* (LSM)) is the most well-known implementation of the RC principle (G. Tanaka
et al. 2019). The main idea behind this model is to drive a RNN with some input signal, which will excite the neurons within the reservoir and induce several nonlinear response signals. This is combined with a trainable linear regression to extract a desired output from a combination of these response signals. The ESN model is illustrated in Figure 2.17.

In practice, the input projection is often fixed, but it can also be trainable. The output can also be fed back to the reservoir as an input, and the structure of the RNN reservoir can be adapted to various applications.

<figure><img src="image_31.png" /><figcaption>**Figure 2.17** – Diagram of an echo-state network. The RNN is represented “flattened” in time. The circles represent the units of the hidden state of the RNN. An arrow between unit i and j corresponds to a non-zero entry in the matrix, $W_{h,ij} \neq 0$. Non-linear operations are not represented.</figcaption></figure>

We define an input sequence as $(x_t)^T \in \mathbb{R}^{n \times T}$, where n is the input dimension and T is the number of input time steps. The initial state of the reservoir at t = 0 is $r_0$. The echo-state network is based on the following update equation:

$$r_{t+1} = (1 - \beta)r_t + \beta \tanh(W_h r_t + W_{in} x_{t+1}), \quad (2.8)$$

where $r_t$ is the K-dimensional state vector corresponding to hidden neurons — at time t, $\beta$ is the leak rate, $W_h \in \mathbb{R}^{K \times K}$ is a sparsely connected random hidden layer matrix and $W_{in} \in \mathbb{R}^{L \times K}$ is the input projection matrix. The matrices $W_h$ and $W_{in}$ have random values.

For the random initialization of $W_h$ and $W_{in}$, Jaeger (2012) recommends: $W_h$ should have an average of 10 non-zeros entries per row, all sampled uniformly in [-1, 1]. The matrix is then scaled to achieve a set spectral radius $\rho$, which is optimized for each experiment in (Jaeger 2012). $W_{in}$ has its entries uniformly sampled in [-1, 1]. The columns of the matrix are then individually scaled with factors $\sigma_1, \dots, \sigma_L$ specific to
each experiment. In this formulation, there are multiple parameters to optimize:

- The size of the reservoir $K$

- The spectral radius of $W_h$, $\rho$

- The input weight scaling parameters $\sigma_1, \dots, \sigma_L$

- The leaking rate $\beta$.

Jaeger (2012) explores three optimization schemes for these parameters:

**Blind.** The input weight scaling parameters are set to a single value $\sigma$. The parameters $K$, $\rho$, $\sigma$, and $\beta$ are optimized. This corresponds to a search space of dimension 4.

**Basic.** Each $\sigma_i$ is optimized individually or in groups. The other parameters are also optimized.

**Smart.** All parameters are optimized as in the **Basic** scheme, and the weights of $W_h$ are also designed specifically for the target task.

The $L$-dimensional outputs are computed at times $t > 0$ as

$$ \tilde{x}_{t+1} = D(r_t), \qquad (2.9) $$

where $D: \mathbb{R}^K \to \mathbb{R}^L$ is a (trained) decoding function. For echo state networks, the decoder is often a linear transformation $D(r_t) = W_{\text{out}}r_t$ where $W_{\text{out}}$ is a matrix of dimensions $K \times L$. In our experiments in Chapter 6, we set $\beta = 0$, which was empirically observed to yield the best results on our tasks. Parameter $\beta$, as well as the randomly sampled weight matrix, are sometimes adjusted for each task (Jaeger 2012). We only use default values to obtain a task-independent setup with the least possible assumptions and to make the methods comparable.

### 2.3.3 Reservoir cellular automata

Cellular automata can also be used as a dynamical system in reservoir computing. Yilmaz originally proposed this in 2014 and later was named ReCA (Reservoir Cellular Automata) (Margem and Yilmaz 2017; Yilmaz 2014).

Several implementations of ReCA systems have been proposed and evaluated on various tasks (Babson and Teuscher 2019; Kleyko, Frady, and Sommer 2020; Margem and
Gedik 2018; McDonald 2017; Morán, Frasser, and Rosselló 2018; Nichele and Gundersen 2017; Nichele and Molund 2017a,b; Yilmaz 2014). In this section, we lay out a general description of reservoir computing with cellular automata based on these previous works and present our proposed extensions. We apply reservoir cellular automata to language-like tasks in Chapter 6.

**Encoding.** Input data is assumed to be categorical and sequential,

$$X = [X_1, \dots, X_t, \dots], \quad t \in \mathbb{N}, \quad \forall t \ X_t \in \mathcal{X} \subset \mathbb{N}. \qquad (2.10)$$

A CA is a discrete system not specifically designed to handle input. The purpose of the encoding step in ReCA is to convert the input data points to vectors that can be embedded in the state space of the cellular automaton. Such encoding should translate — or embed — the input space $\mathcal{X}$ into a new one that can be combined with the current state of the cellular automaton. The goal is to make a cellular automaton “react” to its input and leverage the resulting computation to solve a problem.

**Input projection.** First, each categorical input vector is one-hot encoded into a vector of size $L_{in}$, where $L_{in} = |\mathcal{X}|$ is the number of input categories. We have

$$\forall t, \quad \mathbf{x}_t \in \{0, 1\}^{L_{in}}, \qquad (2.11)$$

with $\sum_{i=1}^{L_{in}} (\mathbf{x}_t)_i = 1$.

Next, the vectors $\mathbf{x}_t$ are projected onto the vectors $\mathbf{p}_t = P(\mathbf{x}_t) \in \mathcal{P} = \{0, 1\}^{L_d}$ of fixed size $L_d$, where $L_d$ is the size of the internal state of CA. Usually, we have $L_d > L_{in}$, which means that we embed the input in a space of higher dimension. In our work, we compute this projection in one of three ways:

**One-to-one** each input bit is assigned to a single index in $\mathbf{p}$. This is the projection function first proposed in Yilmaz (2014).

**Building the reservoir** The final state of the reservoir $\mathbf{r}_{t+1}$ is obtained by stacking $I$ consecutive CA states obtained from a combination of inputs of a single state $\mathbf{s}'$ by applying $\Phi$ again. The parameter $R$ can be understood as a space redundancy parameter,
while *I* is a time redundancy parameter. If the size of the state is *n*, the resulting reservoir vector *r* has dimension *K* = *I* × *n*, where *I* is defined above. Similar to the ESN, the *L*-dimensional output tokens are calculated at times *t* > 0 as

$$ \tilde{x}_{t+1} = D(r_t), \qquad (2.12) $$

where *D*: $\mathbb{R}^K \to \mathbb{R}^L$ is the (trained) decoding function. Here, we predict $\tilde{x}_{t+1}$, which corresponds to a language-modeling-like task of predicting the next input. It would also be possible to predict another sequence $(y_t)_t$. Usually a linear decoding function is used such that $D(r_t) = W_{\text{out}}r_t$. The complete pipeline of a RC model with CA is summarized in Figure 2.18.

<figure><img src="image_32.png" /><figcaption>**Figure 2.18** – Illustration of the 1D ReCA model with the XOR-based input and state combination method. The highlighted boxes correspond to *I* stacked CA states that make up a single reservoir state.</figcaption></figure>
# Literature Review

In the previous chapter, we gave a detailed definition of cellular automata, neural networks, and reservoir computing, and described the interplay between them. The current chapter aims to provide an overview of several fields associated with learning and complex systems that have informed and influenced our work in this thesis. We begin by reviewing related work on complexity measures for dynamical and complex systems (Section 3.1). We also review works on defining and understanding emergence (Section 3.2), evolutionary methods for searching solutions in complex spaces (Section 3.3), methods for constructing open-ended evolving systems (Section 3.4), and methods for extracting and using the computations happening within complex systems (Section 3.5).

## 3.1 Measuring complexity

<figure><img src="image_33.png" /><figcaption>**Figure 3.1** — *Ideal* “complexity” curve as a function of the “randomness” of a system. The three CAs from figure 3.2 are displayed at their approximate expected location on that curve. This curve is just for illustrative purposes and does not correspond to any true function linking randomness and complexity.</figcaption></figure>

Measuring the complexity of a system is a fundamentally difficult task. Many complex systems exhibit what Peter Grassberger calls *self-generated complexity* (Grassberger 1986). This means that the formulation of the system is translationally invariant and the observed structure arises from a spontaneous breakdown of translational invariance. For
example, CA has a uniform update rule, but the complexity arises locally in some of them. Unfortunately, there is no universally accepted and formalized notion of “complexity”, although most intuitively agree that it exists. For example, Figure 3.2 shows three examples of behaviors generated by 1D ECAs. Most people would consider the leftmost figure 3.2a to be not complex. However, depending on one’s definition of complexity, the last picture could be labeled as complex.

<figure><img src="image_34.png" /><figcaption>**Figure 3.2** – Three 1D Cellular automata with qualitatively different behavior. The left (3.2a) and right (3.2c) CAs are usually not defined as “complex”, whereas the middle CA (3.2b) appears to be more complex than the others. The Figures show CAs space-time diagrams. Each row represents the state of the CA at a single time step, with the time increasing from top to bottom.</figcaption></figure>

No clear observable quantities and protocols have been proposed that would give a quantitative notion of complexity. Grassberger (1989a) argued that no single quantity is sufficient to measure complexity since it depends on exactly how the meaning is assigned to this term. Even when some meaning has been fixed by the definition of an observable quantity, the statistics of the measurement of this observable are also crucial to its interpretation (H. Gutowitz 1995).

Despite these drawbacks, we seek a measure of complexity that matches intuition as well as possible. Such a measure would assign low complexity to both simple objects and random objects and high complexity to objects in between. For illustration, we place the three ECAs from Figure 3.2 on this ideal complexity curve, shown in Figure 3.1. For example, in Chapter 4, we develop a complexity metric that matches the human intuition of complexity on a dataset of cellular automata.

One of the most widely known concepts of the complexity of symbol sequences, “algorithmic complexity,” should also be called a measure of information. This is the position held by one of its co-creators (Gregory J. Chaitin 1990). In the literature, information measures are often used as complexity metrics. These quantities are better measures of “randomness” than complexity, but they are nonetheless important for the
study of complex systems. In this section, we list several complexity measures, including information measures that have influenced our research.

### 3.1.1 Information content

For an event $E$ with probability $P$, the information content of this event is defined as the negative logarithm of its probability, that is,

$$I(E) := -\log(P). \qquad (3.1)$$

This metric quantifies how unlikely an event is and depends on how the probability $P$ was estimated in practice. This notion is at the foundation of information theory. The information content of a finite set of $N$ events is maximal when all events are equally probable, that is $I(E) = \log(N)$.

### 3.1.2 Shannon Entropy

Shannon entropy is defined as the expected information content of the input (Shannon and Weaver 1975). We have

$$H(X) := \mathbb{E}[-\log(P(x))]_{x \in X}.$$

This measure is a lower bound on the number of bits that the input could be compressed to.

In the case of a 1D ECA, the Shannon entropy could be computed cell-wise over the time distribution of the states, producing an entropy per cell score that can be averaged over an entire automaton. This measure is, for example, one of the measures used by Wolfram in (Wolfram 1983) and by Langton in (Langton 1990) to study how the parameter $\lambda$ affects the behavior of the automaton. There are several ways to compute it when dealing with a CA, depending on which part of the CA is considered as the main random variable. If the CA is finite, the state at timestep $t$ can be seen as a random variable that can take one of $2^N$ possible values (with $N$ the width of the automaton state). In that case, the probability of a state can easily be estimated by counting
its number of occurrences during evolution. For too large state spaces, this becomes challenging, as it is very unlikely that a given state will be seen again.

### 3.1.3 Rényi Entropy

The Rényi entropy is a generalization of Shannon entropy that gives different weights to events of various probabilities (Rényi 1961). It is formally defined for an $\alpha \ge 0, \alpha \ne 1$ and a random variable $X$ with possible outcomes 0, 1, ..., $n$ as

$$H_{\alpha}(X) = \frac{1}{1-\alpha} \log \left( \sum_{i=0}^{n} p_i^{\alpha} \right).$$

In the limit $\alpha \to 1$ it is equal to the Shannon entropy, where the probabilities of each event are equally important. $\alpha \to \infty$ yields the Min-entropy, and we have $H_{\infty}(X) = \min_i(\log(p_i))$. With $\alpha \to 0$ Rényi entropy is the same as Max-entropy, or topological entropy. Rényi entropy was used by Wolfram (1983) and Lindgren and G. Nordahl (1988) to estimate the complexity in infinite CAs.

### 3.1.4 Mutual information

The mutual information is a measure of the dependence between two variables. Two independent variables will have mutual information of zero. If the two are strongly dependent, for example, if one is a function of another, the mutual information between them will be larger.

Gregory J. Chaitin (1987) proposes to split a system into fixed blocks and calculate mutual information among its components. He uses the maximum value of the mutual information in all possible partitions to define "life" in a mathematical way. Shaw (1984) and Grassberger (1986) also use the mutual information measure between two semi-infinite blocks in a sequence to define complexity. Mutual information has also been used to study the complexity of CAs as the parameter $\lambda$ is changed (for details on $\lambda$ in the context of CA, see section 2.1.4) (H. A. Gutowitz and Langton 1988; Li, Norman H. Packard, and Langton 1990).
### 3.1.5 Computational complexity

Algorithms can be described by their space and time complexity, which respectively correspond to the amount of memory storage and CPU time necessary to run them (Hopcroft, Motwani, and Ullman 2007; Packel and Henryk Woźniakowski 1987; Traub, Wasilkowski, and H. Woźniakowski 1983). A problem that depend on a parameter *N* will be said to be NP-hard if the time needed to find a solution to them increases exponentially with the parameter *N*. NP-hard problems can be said to be complex as they require a supposedly irreducible amount of computation to compute their output. In practice, this description only applies to computations that were generated by a hand-designed algorithm and not the kind of computations self-generated by a complex system that we are interested in. For example, the computations generated by one CA cannot be said to be computationally more complex than those of another CA because the algorithm that implements them is fundamentally the same and computationally simple.

### 3.1.6 Solomonoff–Kolmogorov-Chaitin complexity (algorithmic complexity)

Introduced by Solomonoff (R. J. Solomonoff 1960), Kolmogorov (Kolmogorov 1968) and Chaitin (Gregory J Chaitin 1977, 1969, 1990), this complexity measure is defined for a string *s* of characters and a universal description language (e.g., a programming language) as the length of the shortest program that can generate the string *s*.

This number $K(s)$ is called the minimum description length of *s*. We note that Chaitin prefers to call his field "algorithmic information theory." His position is that "algorithmic complexity" is a measure of randomness rather than a measure of complexity (Gregory J. Chaitin 1990).

From the invariance theorem (Gregory J. Chaitin 1969; Kolmogorov 1968; R. Solomonoff 1964a,b), the difference in algorithmic complexity of the same string *s* in two different description languages is bounded, although this bound might be very large in practice.

The algorithmic complexity is uncomputable (Kolmogorov 1968; R. Solomonoff 1964a), and there exists strings of arbitrarily large complexity, which makes it hopeless to use this exact measure in practice. It can be approximated from above, but no accuracy
guarantee can be given, and the runtime of the approximators can grow arbitrarily large.
Moreover, algorithmic complexity tells us how much information is required to encode a
number, but does not tell us how difficult it is to recreate the number from that code
(Gell-Mann 1988).

A range of resource-bounded approximations of algorithmic complexity were devel-
oped to obtain computable complexity measures (R. P. Daley 1973; R. Daley 1977; Feder,
Merhav, and Gutman 1992; Ko 1986; Schmidhuber 2002). One rather straightforward
way of approaching the algorithmic complexity of an arbitrary string is to use a compres-
sion algorithm and use the length of the decompression program plus the length of the
compressed string as an upper bound to the algorithmic complexity. This is similar to
the optimal symbol compression ratio (OSCR) algorithm proposed by (Evans, Saulnier,
and Bush 2003). Zenil (2010) used the compression-based method to classify 1D ECA.
This method, which often uses the popular LZ algorithm, can be seen as an independent
complexity measure, also closely related to Lempel-Ziv complexity, described in more
detail in the next section 3.1.7.

However, it is worth noting that when a complex system is described by an algorithm
(such as for a CA), the algorithmic complexity is also easily upper bounded by a constant
value entirely defined by the algorithm that simulated the complex system. For example,
for a CA, the algorithmic complexity is lower than the size of the rule table of the
automaton, its characteristics (size, boundary conditions, etc.), its initial state, and the
number of simulation steps. This approximation does not inform us about the difference
in complexity between the two CAs.

3.1.7 Lempel-Ziv complexity

The Lempel-Ziv complexity, as defined in (Lempel and Ziv 1976) is the number of
steps in the LZ algorithm, directly related to the number of repeated substrings in the
input string. The main idea of this algorithm is to scan the input string while trying
to find some repetition of the previous input in the incoming data. This builds over
time a set of basic components called the exhaustive history of the string, from which
the complete string can be constructed. The number of components in that set is the
Lempel-Ziv complexity of that string.
The compressed length method for measuring complexity uses a compression algorithm to reduce the size of the input. An input with regularities and repetitive patterns will produce a small output, whereas a completely random input will be irreducible to a simpler, shorter string.

### 3.1.8 Logical Depth

Bennett's logical depth is a measure of complexity based on the algorithmic complexity (C. H. Bennett 1988, 1995). The main difference is that it takes into account the computation time (or number of steps) along with the length of the program used to generate the sequence. It is a combination of algorithmic complexity and computational complexity. For a universal computer $U$, the logical depth of a string $x$ at a significance level $s$ is defined as

$$ \min\{T(p) : (|p| - |p^*| < s) \wedge (U(p) = x)\}, \qquad (3.2) $$

where $T$ is the number of steps that the universal computer runs the program $p$ for. This quantity corresponds to the least time required to compute the string $x$ by a $s$-incompressible program, that is a program $p$ whose length is within $s$ symbols from the optimal program $p^*$ of the algorithmic complexity metric.

Luis Antunes et al. (2006) considered logical depth to be one instance of a more general concept: computational depth; and proposed several other variants.

### 3.1.9 Thermodynamic Depth

The thermodynamic depth, developed by Lloyd and Pagels (1988) is another measure of complexity based on the intuitive notion that complex systems lie somewhere in the continuum between order and chaos (Ceccatto and Huberman 1988; Gregory J. Chaitin 1990; Deutsch 1985). Similar to computational complexity and logical depth, this metric is designed to be a measure of how the system came to be in its final state. The complexity corresponds to how difficult it is to put the system together in that state.

By the definition of thermodynamic depth, the average complexity of a state must be proportional to the Shannon entropy (Shannon and Weaver 1975) of the set of trajectories
that the experiment determines can lead to that state. If we write $p_i$ the probability of reaching the target state with the i-th possible trajectory, we have

$$S = - \left( \sum_i p_i \log p_i \right). \qquad (3.3)$$

The measure of complexity of a macroscopic state $s$ of a system that has reached that state by the i-th possible trajectory is $-k(\log p_i)$, where $p_i$ is the probability that the system has reached $s$ by the i-th trajectory, and $k$ is an arbitrary positive constant. The thermodynamic depth of the state $s$ is defined as

$$\mathcal{D}(s) = -k(\log p_i). \qquad (3.4)$$

The quantity $\mathcal{D}$ can be thought of as the amount of information required to specify the trajectory that the system has followed to its present state. The thermodynamic depth of the whole system is then

$$\mathcal{D} = \sum_s \mathbb{P}(s) \mathcal{D}(s), \qquad (3.5)$$

where $\mathbb{P}(s)$ is the probability that the system followed the trajectory $s$.

James P. Crutchfield and C. R. Shalizi (1999) argued that the thermodynamic depth is a fundamentally flawed structural complexity measure because it relies on a set of chosen macroscopic states for the system, which are difficult to choose and need to be defined separately for each system. Lloyd and Pagels (1988) does not mention how these states are supposed to be chosen, nor do any follow-up work, making this complexity metric difficult to use in practice.

### 3.1.10 Epsilon-machines

The $\epsilon$-machine is an approach to complexity that seeks to construct a metric more suitable for physical systems that also addresses some issues of other existing complexity metrics (James P. Crutchfield 2012). It is the result of the field of computational mechanics, an extension of statistical mechanics that describes not only the statistical properties of a system but also how it stores information and how it computes (James P. Crutchfield 1994, 2012; James P. Crutchfield and Young 1989; Feldman and James P
Crutchfield 1998). Computational mechanics algorithms take as input the time series being analyzed and output a minimal, optimized model that can reproduce a time series that is statistically equivalent to the input time series. The size of the model produces a metric known as statistical complexity. The models are known as $\epsilon$-machines. Three optimality theorems say that $\epsilon$-machines capture all the properties of a process (James P. Crutchfield and C. R. Shalizi 1999; James P. Crutchfield and Young 1989; C. R. Shalizi and James P. Crutchfield 2001): prediction: the $\epsilon$-machine is its optimal predictor; minimality: compared to all other optimal predictors, the $\epsilon$-machine of a process is its minimal representation; uniqueness: any minimal optimal predictor is equivalent to the $\epsilon$-machine.

This model has many attractive properties and was successfully applied to study the symbolic dynamics of chaotic systems (James P. Crutchfield 1994), molecular dynamics (Ryabov and Nerukh 2011), single-molecule microscopy (Kelly et al. 2012), and the spatiotemporal complexity of CAs (James P. Crutchfield and Hanson 1993; Hanson and James P. Crutchfield 1997; C. R. Shalizi, K. L. Shalizi, and Haslinger 2004). The main drawback is the construction of the $\epsilon$-machine itself, for which there is no general-purpose approach.

### 3.1.11 Sophistication

Intuitively, sophistication is the complexity in a set of strings of which the string is a “typical” member. Mota et al. (2013) defines sophistication based on the original definition of Koppel (1988) and Koppel and Atlan (1991). It measures the amount of structural information contained in a string. It also uses a follow-up result by (Vitanyi 2006) that shows that Koppel’s definition is equivalent to measuring the complexity of a good model for a string up to low order terms. The definition uses an intermediate quantity called discrepancy that measures how far a set $S$ is from being a good model of a string $x$.

Formally, the discrepancy $\Delta$ is defined for a string $x$ and a set $S$ that contains $x$ as

$$ \Delta(x|S) := \log |S| - K(x) + K(S), \qquad (3.6) $$
where K is the Solomonoff–Kolmogorov-Chaitin complexity function.

The sophistication of $x$ is defined as the complexity of the simplest model of $x$ with a limited discrepancy:

$$ \text{soph}_c(x) := \min_S \{K(S) : \Delta(x|S) \le c\} \qquad (3.7) $$

where the significance level $c$ tells us how much discrepancy the set $S$ is allowed to have.

Koppel's definition of sophistication, Definition 3.1, may not be stable. Small changes in $c$ could cause large changes in $\text{soph}_c(x)$. For this reason, Luís Antunes and Fortnow (2009) introduces a new notion of *coarse sophistication* that incorporates the "constant" $c$ as a penalty in the formula to obtain a more robust measure. In the formalism of Mota et al. (2013) the definition is

$$ \text{csoph}(x) = \min_S \{K(S) + \Delta(x|S)\} \qquad (3.8) $$

which is also equivalent to

$$ \text{csoph}(x) = \min_c \{\text{soph}_c(x) + c\}. \qquad (3.9) $$

Like algorithmic complexity, sophistication is a useful theoretical notion to model ideas of entropy and complexity, but it cannot be directly applied in numerical simulations and can only be approximated because it is also uncomputable.

## 3.2 Emergence

The notion of emergence is central to the study of complex systems. It is a broad concept that could be phrased as *the properties a composite entity acquires that its composing parts did not possess*. Emergence is often described as being similar to self-organization. Although these two concepts are closely related, they are fundamentally different ideas (De Wolf and Holvoet 2005). Both properties can exist independently within a dynamical system. Some examples of emergent behavior include:

- The spontaneous formation of clusters of randomly distributed objects - a behavior common in ant colonies forming bridges out of individual ants or birds flocking into
"murmurations" - that naturally emerges out of a simple set of autonomous actions having nothing to do with that clustering. This was demonstrated by Beckers, O. E. Holland, and Deneubourg (2000) in the context of exploring collective robotics). The macroscopic behavior in each of these examples is unexpected even though the details of the microscopic dynamics are well-defined.

— The spirals of the Belousov- Zhabotinsky chemical reaction (Tyson 2013).

— The Navier-Stokes-like macroscopic behavior of a lattice gas that consists, of simple unit-bit billiards moving back and forth between discrete nodes along discrete links at the microscopic scale (Hasslacher 1987).

There seems to be qualitatively different types of emergence, which de Haan (2006) describes as Discovery, Mechanistic emergence, and Reflective emergence. These correspond to different levels, or “strengths”, of emergence, from fractal patterns to complex social systems and ecosystems.

A subset of emergence research has focused on complex adaptive systems. In such systems, emergence is explicitly used to refer to macro-level patterns arising from interacting agents (J. H. Holland 2000; S. Kauffman 1995; Langton 1986).

## 3.3 Evolutionary algorithms

Evolutionary algorithms are the class of algorithms inspired by natural biological processes, in particular evolution. There are many types of evolutionary algorithms, from genetic algorithms to quality diversity and novelty search (Lehman and Stanley 2011a,b). Evolutionary algorithms are used as search methods, in particular when the search space is hard to parametrize or very large (Poli and Logan 1996). Evolutionary algorithms facilitate the search for complex optimization spaces by using intermediate representations, as illustrated in Figure 3.3. CA rules are a good example of such search spaces. The CA rule space is very large, and it is difficult to predict how perturbations of the rule translate to changes in the behavior of a CA. Genetic algorithms were successfully used to evolve rules to perform complex computations (Melanie Mitchell, James P Crutchfield, and Das 1996).
<figure><img src="image_35.png" /><figcaption>**Figure 3.3** – Illustration of the general principle common to several evolutionary algorithms. There are three virtual spaces: (G) the genotype space where the candidate solutions are generated through some encoding. This is the space where the search process is happening. (P) the phenotype space, which corresponds to the decoded solution from a genotype x. It can be the policy of an agent, the rule of a CA, etc. (B), the behavior space, which is the space from which candidate solutions are evaluated. The smaller subspace G is the goal space, corresponding to the set of behaviors that would be considered successful. A small perturbation in the genotype will have potentially large consequences on the behavior of the candidate solution. In the case of a CA rule, the genotype could be a vector (x) that encodes a rule (p) that will translate into the CA behavior b.</figcaption></figure>
### 3.3.1 Genetic programming

Many of the problems that machine learning and artificial intelligence are attempting to solve require the discovery of a computer program that produces some desired mapping between inputs and outputs. Solving the problems therefore amounts to searching the space of computer programs in order to find a suitable individual with high enough fitness (Banzhaf et al. 1998; Booker, Goldberg, and J. H. Holland 1989; J. R. Koza 1994; Langdon and Poli 2002).

In genetic programming, candidate solutions to a problem are encoded as *chromosomes*, which is a structured data representation that can be modified incrementally. In practice, it is often chosen to be a string of bits. Every generation, a population of solutions is evaluated. The best solutions are kept and combined to form the candidates for the next generation. Random mutations may also be applied to introduce randomness in the search.

### 3.3.2 Novelty search

The idea behind novelty search (NS) is to drive a search algorithm only by the novelty of produced behavior (Lehman and Stanley 2011a). The most counter-intuitive feature of this algorithm is that the actual objective of the task is not taken into account at all during the search process. Generated agents are evaluated solely on the novelty of their behavior. Despite this, it appears to be at least as efficient as other search processes that are focused on goals, such as maze navigation and biped locomotion (Lehman and Stanley 2011a), swarm robotics (Gomes, Urbano, and Christensen 2013), or neural network design (Risi, Hughes, and Stanley 2010). NS can lead to the discovery of innovative and creative solutions to complex and deceptive objectives that would be difficult or impossible to achieve using traditional optimization methods. We illustrate this property in Figure 3.4.

The principle of guiding search by novelty alone is closely related to the goal of this thesis, and especially with the challenge of parametrizing and searching the space of available complex systems without using any explicit goal function. The purpose of NS is to perform a search without a goal function, and we draw inspiration from this algorithm throughout our work.
<figure><img src="image_36.png" /><figcaption>**Figure 3.4** – Illustration of the advantage of novelty search (NS) on a maze environment with deceptive dead-ends. Dotted trajectories illustrate the behavior of a goal-based search, which will start by minimizing the distance to the goal. Blue trajectories are from a novelty search algorithm exploring the space more efficiently. This type of maze environment was proposed by Lehman and Stanley (2011a).</figcaption></figure>

This method has some limits, which are caused by a relative lack of theoretical understanding of the effects of the parameters of NS. The choice of behavior descriptors for determining the novelty of agents is also not straightforward. Several empirical studies have investigated the impact of a range of parameters on the performance of NS (Gomes, Mariano, and Christensen 2015; Kistemaker and Whiteson 2011), while other works have begun to shed light on the theory (Doncieux, Laflaquière, and Coninx 2019). In practice NS algorithms are expected to cover the exploration space well (Cully and Demiris 2017; Pugh, Soros, and Stanley 2016), in some cases uniformly (Gomes, Mariano, and Christensen 2015). NS algorithms are related to genetic algorithms since they maintain a population, that is, a set of individuals used to measure the novelty of current solutions. There is also an archive of the behavior of previous individuals to ensure that the novelty is measured against all previously generated behaviors. There exists multiple strategies for maintaining this archive (Gomes, Mariano, and Christensen 2015). Some use individuals whose novelty was above a threshold when first evaluated (Lehman and Stanley 2011a), the most novel individuals of each generation (Liapis, Yannakakis, and Togelius 2015), random individuals from the past generations (Lehman and Stanley 2010), or none at all (Mouret and Doncieux 2012). An alternative to
explicitly measuring novelty in NS is *fitness sharing*. The goal of fitness sharing is to encourage the diversity of a population by making them share virtual “resources” (Goldberg 1987; J. H. Holland 1992). A particular implementation of NS using this idea, with rewards physically spread across the environment, was done by Herel et al. (2022).

Typically, novelty search algorithms use one of a few types of maps from genotypes to phenotypes, that is, a map from the space of encoded solution to agent and policy implementation (Mouret and Doncieux 2012). Some common ones are:

- A RNN (of type Elman 1990 or Jordan 1997), for which the weights are directly searched. For input size $n_i$, output size $n_o$, and hidden size $n_h$, $n_i n_h + n_h n_o + n_h^2$ weights must be evolved.

- Neuro-evolution of augmenting topology (NEAT, Stanley and Miikkulainen 2002), which is an encoding of both neural network weights and architectures that has a built-in mechanism to sustain diversity.

## 3.4 Open-ended evolution

The field of Artificial Life research has worked to figure out what the fundamental conditions for the emergence of living systems are, and how to create a process that can display analogous levels of creativity and complexity as natural evolution (Dyson 1999; Eigen and Schuster 1979; Langton 1989; N. Packard et al. 2019; Soros, Stanley, and Lehman 2017; Stanley 2019). It builds upon the data and understanding we collected about the process of life, but abstracts from any specific living process and attempts to integrate various approaches into one unified research to extract the first principles of life. A major assumption underpinning this research is that this natural evolutionary process can be implemented equally well in different media (Dennett 1996).

A system that behaves like natural evolution, producing a seemingly endless amount of novelty and complexity starting from elementary building blocks is called *open-ended*. The main challenge of open-ended evolution (OEE) research is that there is no single simple test for the phenomenon, but instead different kinds of open-ended evolution exist. Systems can exhibit more than one kind at a time. In the report from a workshop on OEE at York (T. Taylor et al. 2016), the authors summarized the different types of OEE, which were further refined in a follow-up work (N. Packard et al. 2019). They are:
1. Interesting new kinds of entities and interactions

2. Evolution of evolvability

3. Major transitions

4. Semantic evolution.

The first category describes the ability of a system to construct new entities with different properties, behaviors, or interactions with other entities. For example, the Tierra simulation sees entities emerge that exploit the computing power of others, acting as parasites (Ray 1991). The second type is related to how open-ended evolution itself can be evolved through the emergence of multiple “stepping stones” that allow individuals to use higher-level evolutionary units or through interactions (Pattee and Sayama 2019). “Major transitions” refer to the emergence of new levels of hierarchy in evolving populations. These transitions are characterized by groups of reproducing entities that interact tightly and form a new population of higher-level reproducing entities. The process can repeat, with certain groups in the new population forming even higher-level entities. This process of emergence of major transitions was explored in (Moreno and Ofria 2019; Sayama 2019) for example. Lastly, semantic evolution is related to the work of Ikegami, Hashimoto, and Oka (2019), who propose a new category of open-ended evolution called semantic evolution, which refers to the evolution of semantic relationships within a system. This can be observed in the evolution of web services, where the meaning of tags changes as new combinations of tags are created. It is not biological in nature and is also present in the analysis of technological evolution (Bedau, Gigliotti, et al. 2019), where the meaning and importance of keywords shift over time.

### 3.4.1 Defining open-endedness

Several necessary conditions and requirements have been identified in the OEE literature, forming an overlapping set of potential research directions for developing open-ended systems. We list a few of these conditions here. First, Maley (1999) identifies four requirements:

1. An open-ended evolutionary system must demonstrate unbounded diversity during its growth phase.
2. An open-ended evolutionary system must embody selection.

3. An open-ended evolutionary system must exhibit continuing new adaptive activity.

4. An open-ended evolutionary system must have an endogenous implementation of niches.

Next, Soros and Stanley (2014) identified four more necessary conditions, which are

1. A rule should be enforced that individuals must meet some minimal criterion (MC) before they can reproduce, and that criterion must be nontrivial.

2. The evolution of new individuals should create novel opportunities for satisfying the MC.

3. Decisions about how and where individuals interact with the world should be made by the individuals themselves.

4. The potential size and complexity of the individuals' phenotypes should be (in principle) unbounded.

Then, T. Taylor (2015) also stated five requirements for an open-ended system:

1. Robustly reproductive individuals.

2. A medium allowing the possible existence of a practically unlimited diversity of individuals and interactions, at various levels of complexity.

3. Individuals capable of producing more complex offspring.

4. An evolutionary search space that typically offers mutational pathways from one viable individual to another viable (and potentially fitter) individuals.

5. Drive for continued evolution.

Taylor also states a more general condition that should be sufficient for creating an open-ended system as “evolutionary dynamics in which new, surprising, and sometimes more complex organisms continue to appear” (T. Taylor 2015; T. Taylor et al. 2016).

All these conditions illustrate a major challenge of OEE research: the lack of a clear definition or notion of what exact conditions make a system open-ended. We seem to agree about what is not open-ended, but whenever a constraint or requirement for OEE is identified, subsequent evidence forces us to refine them later. This is related to the problem of complexity, for which no single definition exists (Johnson 2009).
A common process to produce systems that behave analogously to natural evolution is to start from evolvable units or building blocks (Channon 2003; Ofria and Wilke 2004; Ray 1991; Sims 1994; Soros and Stanley 2014; Spector, Klein, and Feinstein 2007; Yaeger 1994). The reason is that starting from higher-level primitive units whose emergence would be hard to characterize may be easier and faster than starting from lower-level components. The results obtained from these systems are often surprising, as they bear some key similarities to natural evolutionary processes. For example, we note the emergence of parasitic entities within the Tierra simulation (Ray 1991). The process of emergence of these building blocks from simple rules and components has also been investigated significantly (Bagley and Farmer 1991; Flamm et al. 2010; Hutton 2007; Sayama 2011). It appears harder to bridge the gap and create high-level evolutionary-like processes and behaviors from elementary rules and substrates.

### 3.4.2 Open-endedness in cellular automata

One class of systems that has rich interactions between each of its components, as well as no predefined evolvable units or assumptions about individuality is the CA. One of the very first automata, Von Neumann's self-reproducing machine, was designed with goals that align with open-ended evolution, which is to build a machine with no central controller and limited local interaction that can self-reproduce as a whole (Pesavento 1995; Von Neumann and Burks 1966). Later, other self-replicating structures such as Langton's loop (Langton 1984) and evoloops (Salzberg and Sayama 2004; Sayama 1999) showed that more properties of open-ended systems can be included in a CA. A potential limitation of CAs is the absence of the notion of conservation of "matter". For example, the game of life can start from a configuration with very few live cells and create many more at no cost during its evolution. Some authors believe that this conservation property is essential to the construction of an open-ended evolving system (T. Taylor 2002).

### 3.4.3 Artificial chemistries

Artificial chemistries, are computational models that aim to simulate the behavior of chemical systems, which are typically used to study the emergence of complex behav-
ior from simple interactions between individual molecules. Most of the time, artificial chemistries (AC) do not try to accurately model chemical processes (as in (Bedau, McCaskill, et al. 2000; Buliga and L. H. Kauffman 2014; Ostrovsky et al. 2001; Sayama 2011) for example). Instead, they build models of the dynamics of complex molecular processes that lead to evolutionary behavior (Dittrich, Ziegler, and Banzhaf 2001). By abstracting away from natural molecular processes, AC tries to uncover fundamental conditions for the emergence of organization, self-maintenance, or self-construction with basic building blocks. There are various approaches for building these artificial chemistries, of which we review a few here.

**Rewriting systems.** Rewriting systems are a class of computational models that use a set of rewriting rules to transform and manipulate strings of symbols. They are composed of entities or symbols that get modified according to their set of syntactic rules. Patterns of symbols or entities are replaced according to these rules. In AlChemy (Fontana and Buss 1994), molecules are represented as expressions in the λ-calculus. The λ-calculus is a mathematical formalism that, like Turing machines, is capable of describing any computable function. In AlChemy, pairs of randomly selected expressions are joined using function application, evaluated, and the resulting value is added back to the population of molecules. This process is repeated to simulate the behavior of a chemical system. Kruszewski and Mikolov (2022) hypothesized that self-reproducing metabolisms emerge as a subset of autocatalyzed reactions within a Turing-complete set. They introduced combinatory chemistry, an artificial chemistry designed to capture the core computational properties of living systems that is based on the rewriting rules of combinatory logic (Curry 1958; Schönfinkel 1924). The generated structures exhibit behaviors that are similar to natural metabolisms. The system is able to represent patterns of arbitrary complexity, and it is able to produce diverse self-reproducing patterns.

**Cellular automata.** Cellular automata (CAs) can be seen as a particular case of lattice molecular systems. For example, the autopoietic system (referring to the ability of an organism to continuously sustain itself through the exchange of matter and energy with its environment, while maintaining its own boundaries and identity) introduced by
F. G. Varela, Maturana, and Uribe (1991) is a 2D square grid with sites that be occupied by *atoms* which is similar to a CA with 4 states. Each of these states is analogous to a chemical component of the system: (Ø) empty site, (S) substrate, (K) catalyst, and (L) monomer. Basic rules are applied asynchronously and define how neighboring atoms interact with each other. Remarkably, stable self-repairing cells spontaneously arise from these basic rules. Their membrane is composed of a chain of monomers, which is maintained by the substrate and catalyst reacting according to the rules. Some key components of that model were investigated in other works, showing that they are crucial for this emergence autopoeisis to be possible (McMullin and F. J. Varela 1997; Zeleny 1977).

## 3.5 Computing with complex systems

The problem of computing within complex systems is closely related to the question of decentralized parallel computation, in general, for which there exists abundant literature. Different paradigms exist for controlling and harvesting the computations within complex systems. Several other names have been used for closely related topics, such as organic computing (Müller-Schloer 2011) which is the study of systems with life-like properties such as self-organization or the ability to adapt to a dynamically changing environment. Agent-based computing (Jennings 1999) focuses on computing systems composed of several relatively simple autonomous agents. Amorphous computing (Abelson et al. 2000; Nagpal 2002, 2008) is about making large amounts of individual computing elements work and ensuring “the cooperation of large numbers of unreliable parts interconnected in unknown, irregular, and time-varying ways”.

### 3.5.1 Computing in cellular automata

Cellular automata (CAs) are decentralized parallel systems with many identical components with local connectivity. Because of these properties, they have the potential to carry out robust and efficient computations. CA-based computing machines could recover from perturbations or carry out computations in stochastic environments. Moreover, they are also interesting for modeling the behavior of natural complex systems.

For more details on CAs, see Section 2.1.
**Von Neumann's self-reproducing CA.** The early developments of CAs were guided by Von Neumann's question "What kind of logical organization is sufficient for an automaton to be able to reproduce itself?" Since this question may admit very simple solutions, additional constraints were added so that the problem does not admit trivial solutions. For example, Von Neumann required that the automaton in question be equivalent in power to a universal Turing machine while having minimal complexity (Von Neumann and Burks 1966). The final model he proposed was an intricate CA composed of a tape containing the instructions to construct the next automaton and a construction unit that would progressively build that new automaton cell by cell. This was an early example of complex computations being carried out by a relatively simple decentralized system such as a CA.

**Universal Computation in Cellular Automata.** It is not difficult to see that a CA can be capable of universal computation. The basic approach is to show that it can simulate a Turing machine, which we assume has an infinite tape. A CA rule can be constructed by reproducing all the steps of a Turing machine's behavior while adding a state to each cell, indicating whether it is active or not, and enforcing that only one cell is active at a time. This turns the parallel CA into a sequential object that simulates a Turing machine. A similar construction was carried out by (A. R. Smith 1971), making a CA with one-to-one correspondence between its time steps and that of a target Turing machine.

Another approach was used to show the universality of the game of life (see details in Section 2.1.3). Instead of simulating a Turing machine, basic logical functions are built from gliders (see figure 2.5a) produced by a game of life structure called *glider gun*. Two gliders that collide within the simulation will either annihilate or keep moving depending on the exact position in which they collide. This allows the construction of NOT, AND and OR gates. For example, see Figure 3.5 for a construction of a NOT gate in the game of life.

Universal computation is an intriguing yet impractical property, and proving that a system is a universal computers often has limited usability, but it serves as a proof of principle that the studied system is theoretically as powerful as any computer. In
<figure><img src="image_37.png" /><figcaption>**Figure 3.5** – Construction of a NOT gate in game of life. The left image represents an input equal to 1, while the second shows an input equal to 0. Output streams behave like a NOT gate. A glider gun and a modified half-frequency glider gun are highlighted. They are necessary for the construction of the gate. Images are adapted from a blog post (Carlini 2020).</figcaption></figure>

the case of the game of life, this was illustrated by the implementation of the OTCA metapixel which is a tileable square cell object that can be used to simulate another instance of game of life within the game. It is important to note that such embedded computers and Turing machines are often very slow and inefficient, and the game of life has brittle structures which would be destroyed by the slightest perturbation, making these constructions unusable for anything else than demonstration.

**Evolving cellular automata with genetic algorithms.** An important advancement beyond the clever yet challenging manual design of computational rules in CAs is the utilization of learning algorithms for automated rule design (Melanie Mitchell, James P Crutchfield, and Das 1996). Genetic algorithms, which are search methods inspired by biological evolution, are one such type of learning algorithm (Booker, Goldberg, and J. H. Holland 1989). See Section 3.3.1 for more details about genetic algorithms.

Some early work on evolving CAs was done by Norman Packard and colleagues (Norman H. Packard 1988; Richards, Meyer, and Norman H. Packard 1990). Richards, Meyer, and Norman H. Packard (1990) used genetic algorithms to learn the rules of CA to match experimental data for the patterns created by dendritic solidification of NH₄Br. In (Norman H. Packard 1988), the author evolves CAs to perform a specific computation,
and observes that the evolution process tends to select CAs with a $λ$ parameter close to the critical value observed by Langton that corresponds to the edge of chaos (see Section 2.1.4 for details about the parameter). However, other authors were unable to replicate this experiment later (Melanie Mitchell, P. Hraber, and James P. Crutchfield 1993). J. Koza (1992) applied genetic algorithms to CAs to generate random numbers.

The problem of density classification was thoroughly explored since it fits the CA paradigm particularly well (Andre, F. H. Bennett, and J. R. Koza 1996; J. P. Crutchfield and M. Mitchell 1995; Das, Melanie Mitchell, and James P. Crutchfield 1994; Melanie Mitchell, James P. Crutchfield, and P. T. Hraber 1994; Melanie Mitchell, P. Hraber, and James P. Crutchfield 1993; Sipper 1996). The goal of this task is to have CAs figure out what the most abundant cell state in an initial binary configuration in one dimension is. This task is trivial for a central controller that can read the state of all cells in the initial grid, but it is much more complex for decentralized systems like CAs that have a small radius of interaction with neighbors. The evolved CAs developed various strategies to solve the task. Some are unsophisticated, such as expanding the 1s or 0s, which does not work on all validation examples. We show one of the more sophisticated solutions discovered during some runs of the genetic algorithm in figure 3.6. That solution uses the propagation of structures within the CA state space to communicate the local state to the rest of the grid, which usually ends with the right answer from the CA.

<figure><img src="image_38.png" /><figcaption>**Figure 3.6** – The particule-based solution to the density classification problem (Melanie Mitchell, James P Crutchfield, and Das 1996). When the density is above 0.5, the CA quickly reaches a state with only black cells, whereas when the density is lower, it goes to a state with white cells only.</figcaption></figure>
**Reliable computation in cellular automata.** The problem of carrying out reliable computations in CAs is studied early on in its history because of the potential of that model for a real hardware implementation. The first studies used probabilistic approaches for solving the error detection and correction problems of automata — which includes other models than CAs (Harao and Noguchi 1973; Neumann 1956; Tsertsvadze 1964; Winograd and Cowan 1963). Other methods make use of the special structure of CAs, but rely on strong assumptions about the likelihood of errors (Harao and Noguchi 1975; Nishio and Kobuchi 1975). Peter Gács has shown how to construct a 1D CA which can perform arbitrarily large reliable computations, assuming that each cell gets an error with a positive probability (Gács 1986). The fault model he describes is important from the point of view of ergodic theory because Gács' result invalidates the “positive probability conjecture” in statistical physics, which states that all one dimensional infinite particle systems with positive transition probabilities are ergodic, which means that it will visit all parts of its state space eventually. For more recent work on fault tolerance by Gács, see (Gács 2001).

**Neural cellular automata.** The principle of NCA is based on the analogy between CAs and two types of neural networks: RNNs and CNNs (see Section 2.2 for more details about the analogy). Following this analogy, CAs can be formulated as a continuous recurrent and convolutional neural network system called a neural cellular automaton (NCA). NCAs can be manipulated like neural networks, using automatic differentiation, backpropagation, and optimization algorithms to make these models learn some target task. The first example used this neural network representation to learn a NCA version of existing CA rules from recorded examples of their evolution (Gilpin 2018). The structure of the training process and the final error of the resulting model are used by the author as a tool to understand the complexity and properties of the original CA.

NCAs were also used to learn to produce stable self-repairing patterns that resemble a target image (Mordvintsev et al. 2020). The training process is structured in such a way that the NCA needs to learn to maintain the pattern stable across time and space, but also recover from random perturbations that are introduced at arbitrary points in time. This results in an interesting demo where a pattern is generated and maintained
in real-time by a NCA while a user can directly apply perturbations*.

<figure><img src="image_39.png" /><figcaption>**Figure 3.7** – Resistance of patterns to perturbations with NCAs. NCA can robustly maintain patterns, even under relatively strong perturbations (10-20% of the pattern’s size). This image was generated with an online notebook based on code from Alexander Mordvintsev†.</figcaption></figure>

The same authors also used NCAs to learn to classify digits from the MNIST database of handwritten digits (Randazzo et al. 2020). The particularity of this model is its ability to perform the classification in a decentralized way, by having the NCA change the color of the digits directly on the CA grid. Cells can only communicate with their direct neighbors, so the model has to reach a consensus by spreading information across the grid. This is done in parallel, so multiple digits can be classified simultaneously through this progressive consensus process‡. Another work used the same model to learn to continuously generate images with the same texture as a target (Niklasson et al. 2021). Again, an advantage of this method is the ability to recover from perturbations.

<figure><img src="image_40.png" /><figcaption>**Figure 3.8** – Examples of stable hybrids obtained with NCAs trained with the method of (Cisneros 2021).</figcaption></figure>

We used NCAs as part of our submission to the Minecraft open-endedness challenge

<footnote>* https://distill.pub/2020/growing-ca/</footnote>

<footnote>† https://colab.research.google.com/github/google-research/self-organising-systems/blob/master/notebooks/growing_ca.ipynb</footnote>

<footnote>‡ https://distill.pub/2020/selforg/mnist/</footnote>
2021§, a competition organized at the Gecco 2021 conference competition track¹ (Cisneros 2021). The goal was to build an open-ended system within the computer game Minecraft. This means building a systems that fulfill the following requirements:

- Divergence: Open-ended algorithms are not expected to slow down or converge but rather keep expanding and generating more complex outputs over time.

- Diversity: Does the algorithm produce entities with strong phenotypic diversity?

- Complexity: Can the algorithm produce complex entities or entities interactions that give rise to complex systems? Are hierarchical and modular structures present?

- Ecological interactions: Do the created entities interact with each other?

- Life-Like properties: Inspiration may be taken from other attributes of living systems.

In each of the examples in Figure 3.7, a separate CA was trained on a pattern starting from a single black pixel. We extend this setup by training a single NCA on multiple patterns from different seeds. Instead of a single black cell, a small number of black cells are arranged in simple shapes. We call these patterns *seeds*. NCAs are neural networks with many redundant parameters. This allowed us to learn (and encode) more than one pattern (or larger patterns) without adding parameters in the architecture. The resulting patterns are slightly lower quality but still stable under perturbations and capable of growing from their corresponding seed pattern.

We turned the process of combining seed patterns into a “game” by allowing user interaction with the evolving CA through the Minecraft video game, using an API to communicate between the program and the game. Our system is not open-ended by itself, but serves as an exploration space and uses human interactions to play with the shapes and seed patterns. This is in the spirit of open-ended game-like systems like Picbreeder (Secretan et al. 2011; Woolley and Stanley 2011) where players choose which generated image to mutate or evolve from, collectively constructing surprisingly complex images. The algorithms can generate endless novelty with the help of human interactions that provide the missing step for making the system potentially open-ended.

<footnote>§. https://evocraft.life/</footnote>

<footnote>¶. https://gecco-2021.sigevo.org/Competitions</footnote>
### 3.5.2 Amorphous computing

The goal of amorphous computing is to study and define the problem of computing with multiple interconnected components that are unreliable, irregular, time-varying, and with limited computational capacity (Abelson et al. 2000). The motivation for this area of research is the development of biological substrates that can compute, function as sensors and actuators, and robustly self-organize to compute with little cost. The precise manufacturing of chemical or biological substrates with the ability to communicate locally through chemical or physical interactions is within reach, making their use as computing vessels particularly attractive. We can embed millions of chips with sensors (Abelson et al. 2000) or program biological cells to serve as logic gates (R. Weiss, G. Homsy, and Nagpal 1998; R. Weiss, G. E. Homsy, and Knight 2002) while relying on cheaper, decentralized parts (Butera 2002).

<figure><img src="image_41.png" /><figcaption>**Figure 3.9** – An example of amorphous computing on a surface. Multiple computing particles are arranged randomly. A wave (represented by the color of the cells) propagates within the medium. They can communicate with neighbors within a radius *r*.</figcaption></figure>

In amorphous computing systems, a medium comprises “computational particles” spread on a surface or mixed in a volume. The particles have limited computing power, with no knowledge of their position or orientation. The particles can communicate with neighbors via some specified mechanism, which may be analogous to biological ones. The two main components of this research are *synthetic biology* and the choice of computing paradigm.

**Synthetic biology.** The design of a substrate from chemical processes or genetically engineered biological cells (R. Weiss, G. E. Homsy, and Knight 2002). It should behave correctly and follow predefined rules.

**Computing paradigm.** How to program individual agents to follow a predefined global goal (Nagpal 2001)?
The growing point language (GPL) is an example of a programming language that enables programmers to specify complex patterns for computational particles, which are internally represented in the particles as state machines (Coore 1999).

Nagpal (2002) developed another programming language inspired by biological cell differentiation (Lawrence 1992; Wolpert 1969) that can be compiled into individual agent programs to follow some global specifications.
# Measuring complexity in evolving complex systems

In this chapter, we propose an approach to measure the growth of the complexity of emerging patterns in complex systems such as Cellular automata (CAs). We discuss several ways in which a metric to measure complexity growth can be defined. This includes approaches based on compression algorithms and artificial neural networks. We believe that such a metric can be useful for designing systems that could exhibit open-ended evolution, which itself might be a prerequisite for the development of general artificial intelligence. We conduct experiments on 1D and 2D grid worlds and demonstrate that using the proposed metric, we can automatically construct computational models with emerging properties similar to those found in Conway's Game of Life, as well as many other emergent phenomena. Interestingly, some of the patterns we observe resemble forms of artificial life (Langton 1993). We test our metric of structural complexity growth on cellular automata, but it can also be applied to a wide range of complex systems.

## 4.1 Introduction

Recent advances in machine learning and deep learning have had success in reproducing some very complex feats traditionally thought to be only achievable by living beings. However, making these systems adaptable and capable of developing and evolving on their own remains a challenge that could be crucial to eventually developing AI with general learning capabilities (for example, as is further discussed in (Mikolov, Joulin, and Baroni 2018)). Building systems that mimic some key aspects of the behavior of existing intelligent organisms (such as the ability to evolve, improve, adapt, etc.) might represent a promising path. Intelligent organisms — for example, human beings but
also most living organisms, if we consider a broad definition of intelligence — are a form
of spontaneously occurring, ever-evolving complex systems that exhibit these kinds of
properties (Booker 2004). The ability to sustain open-ended evolution appears to be a
requirement in order to enable the emergence of arbitrarily complex adaptive systems.

Although a rigorous attempt at defining intelligence or life is beyond the scope of
this work, we assume that a system we might identify as evolving, with the potential
of developing intelligence, should have the property of self-preservation and the ability
to grow in complexity over time. These properties can be observed in living organisms
 Booker 2004) and should also be part of computational models that aim to mimic them.

To recognize self-preservation and growth in complexity, one should be able to detect
emerging macro-structures composed of smaller elementary components. For the purpose
of obtaining computational models that grow in complexity over time, one should also
be able to determine the amount of complexity these systems contain. In this chapter,
we propose and discuss several ways to estimate complexity and detect the presence of
emerging and stable patterns in complex systems such as cellular automata. We show
that such metrics are useful when searching the space of cellular automata with the
objective of finding those that seem to evolve over time.

4.2 Related work

4.2.1 Artificial life and open-ended evolution

Several works have attempted to artificially create open-ended evolution. A nonex-
haustive list of some well-known systems include Tierra (Ray 1991), Sims (Sims 1994),
Avida (Ofria and Wilke 2004), Polyworld (Yaeger 1994), Geb (Channon 2003), Division
Blocks (Spector, Klein, and Feinstein 2007) and Chromaria (Soros and Stanley 2014).
Designs that focus on an objective and make use of reinforcement learning methods
to drive evolution are also being studied, e.g. in (Pathak et al. 2019). Most of these
simulated “worlds” have had some success reproducing key aspects of evolving artificial
life, enabling the emergence of complex behavior from simple organisms. However, they
still work within constrained simulated environments and usually consider organisms
composed of elementary building blocks, while they do not work outside of this usually
very constrained framework. The divergent and creative evolutionary process could be
happening at a much lower conceptual level with fewer assumptions. For this reason,
we consider cellular automata in the rest of the chapter because they rely on very few
assumptions while offering a very large expressive power and a potentially wide range
of behaviors that can be discovered. However, the metrics defined in this work have the
potential to be applied to other types of complex system, as discussed in Section 4.7.

4.2.2 Cellular automata

Cellular automata are very simple systems, usually defined in one or two dimensions,
composed of cells that can be in a set of states. The cells are updated in discrete
time steps using a transition table that defines the next state of a cell given the states
of its neighbors. They were originally proposed by Stanislaw Ulam and studied by
Von Neumann (Von Neumann and Burks 1966), who was interested in designing a
computational system that can self-reproduce itself in a non-trivial way. The trivial
self-reproducing patterns were then those that do not have the potential to evolve, for
example, the growth of crystals.

Stephen Wolfram later took a more bottom-up approach, beginning with the study
of simple 1D binary cellular automata (CA), and identifying four qualitative classes of
cellular automaton behavior (Wolfram 1984):

**Class 1** evolves to a homogeneous state.

**Class 2** evolves to simple periodic patterns.

**Class 3** produces aperiodic disordered patterns.

**Class 4** produces complex aperiodic and localized structures, including propagating structures.

Wolfram and his colleagues also studied 2D CA using tools from information theory
and dynamical systems theory, describing the global properties of these systems in terms
of entropies and Lyapunov exponents (Norman H Packard and Wolfram 1985).

Christopher Langton and colleagues also studied CA dynamics (Li, Norman H.
Packard, and Langton 1990) — e.g. using the $\lambda$ parameter (Langton 1990) — and
designed a self-replicating pattern much simpler than Von Neumann's (Langton 1984),
now known as Langton loops. The main issues with his system and Von Neumann's
universal replicator is the fact that they are very fragile and based on a large amount
of human design. As a consequence, although they self-replicate, they cannot increase
in complexity and are not robust to perturbations or unexpected interactions with the
environment.

A genetic algorithm-based search for spontaneously occurring self-replicating patterns was carried out in 2D cellular automata with several states was undertaken in (Bilotta and Pantano 2011) using entropy measures of the frequency distribution of 3 × 3 patterns.

4.2.3 Compression and complexity

Compression has often been used as a measure of complexity. Lempel and Ziv have introduced in (Lempel and Ziv 1976) the now widely used Lempel-Ziv (LZ) algorithm as a method for measuring the complexity of a sequence. By constructing back-references to previous parts of a string, the LZ algorithm is capable of taking advantage of duplicate patterns in the input to reduce its size. The DEFLATE algorithm that we use in the following section combines LZ with Huffman coding for an efficient representation of the symbols obtained after the first step. It is one of the most widespread compression algorithms and is, for instance, used in gzip and PNG file compression standards.

The PAQ compression algorithm series (Mahoney 2000) is an ensemble of compression algorithms initially developed by Matt Mahoney with state-of-the-art compression ratio on several compression benchmarks. Better compression of an input means a better approximation of the minimum description length and implicit understanding of more of the underlying patterns in input data. The usefulness of a better compressor is that it can detect much more complex and intricate patterns that aren't simple repetitions of previous patterns.

In (Zenil 2010), H. Zenil investigates the effects of a compression-based metric to classify cellular automata and observes that it results in a partitioning of the space of 1D CA into several clusters that match Wolfram's classes of automata. He also used this approach in the output of simple Turing machines and some 1D automata with more than two states and larger neighborhoods. Extensions of this work include an asymptotic sensitivity analysis of compressed length for input configurations of increasing complexity,
as introduced in (Zenil 2014; Zenil and Villarreal-Zapata 2013).

Additionally, the image decompression time as an approximation of Bennet's logical depth (C. H. Bennett 1995; Zenil, Delahaye, and Gaucherel 2012) and the output distribution of simple Turing machines combined with block decomposition of CA to approximate their algorithmic complexity have also been investigated (Soler-Toscano et al. 2014; Zenil, Soler-Toscano, et al. 2015). However, the possible extent to which such measures of complexity could be applied to more complex automata and other complex systems has not yet been extensively studied. For a review of several measures of complexity and their applications, see (Grassberger 1989b).

## 4.3 Compression-based metric

A cellular automaton of size *n* in 1D can be represented at time *t* by its grid state $S^{(t)} = \{c_1^{(t)}, ..., c_n^{(t)}\}$ where each $c_i$ (also called a cell) can take one of *k* possible values (representing the possible states) and a transition rule $\phi$. The transition rule is defined with respect to a neighborhood radius *r* with the mapping $\phi(c_{i-r}^{(t)}, ..., c_i^{(t)}, ..., c_{i+r}^{(t)}) = c_i^{(t+1)}$ that maps $\{1, ..., k\}^{2r+1}$ to $\{1, ..., k\}$. The quantity $2r + 1$ corresponds to the number of cells taken into account to calculate the next state of a cell, namely the cell itself and the *r* neighboring cells in both directions.

Simulating a CA amounts to the recursive application of this mapping $\phi$ to an initial state $S^{(0)} = \{c_1^{(0)}, ..., c_n^{(0)}\}$.

In the rest of the chapter, we consider cyclic boundary conditions for the automata, meaning that the indices $i-r, ..., i+r$ above are taken modulo *n* the size of the automaton in 1D. Boundary conditions can have some effect on a CA's evolution, but cyclic boundaries have been empirically observed to have limited effect on the complexity of automata in 1D (LuValle 2019).

The definition given in the equation above can be extended to higher dimensional automata by modifying the neighborhood and the definition of $\phi$. A 2D neighborhood of radius 1 can be defined as the 3 by 3 square around the center cell — also called the Moore neighborhood — or by only considering the four immediate horizontal and vertical neighbors of the center cell — the Von Neumann neighborhood.

Elementary cellular automata (ECA) are 1D CA with $k = 2$ and $r = 1$. There are
$2^3$ elements in $\{1, ..., k\}^{2r+1}$ and $2^{2^3} = 256$ possible different set transition rules that are often compactly represented as a binary string with 8 bits. The relatively low number of rules of this type makes it possible to appreciate the performance of a metric and compare it with others.

We define the compressed length $C$ of a 1D cellular automaton at time $t$ as

$$C(S^T) = \text{length}(\text{comp}(c_1 || c_2 || ... c_n)) \qquad (4.1)$$

where $||$ denotes the string concatenation operator, and cells $c_i$ are implicitly converted into string characters (with one symbol per unique state). comp is a compression algorithm that takes a string as input and outputs a compressed string, and length is the operator that returns the length of an input string.

Similarly to (Kowaliw 2008; Zenil 2010), we use zlib's C implementation of DEFLATE to compress the final state of the automaton. If we apply the above metric to the 256 ECA run for 512 timesteps and initialized with one activated cell in the middle, we obtain the plot of Figure 4.1b. This example is re-used in the rest of the chapter as a way to easily visualize and check that the defined complexity measures are coherent with one another. The colors in Figure 4.1b were obtained with a KMeans clustering algorithm applied on the compressed length of the automata states.

As visible on Figure 4.1b, rules are clearly separated into several clusters that turn out to match Wolfram's classification of ECA. Class 3 behavior can be found at the top of the plot (highest compressed length, orange and blue clusters), Class 1 and 2 are clearly separated at the bottom part (not detailed here) and Class 4 rules (colored in green) lie in between the other types of behavior. The 6 highest scoring rules are shown on Figure 4.1a and correspond to Class 3 behavior in Wolfram's classification. Among the classes of behavior, some sub-clusters can be found that correspond to similarly behaving rules.

Ultimately, the theoretical goal of using compression algorithms is to approach the theoretical minimum description length of the input (Grunwald 2007). For very regular inputs, this length should be relatively small and inversely for random inputs. However, gzip and PAQ are crude approximations of the minimum description length and may
<figure><img src="image_42.png" /><figcaption>**Figure 4.1** – Compression-based metric on 1D ECA. 4.1a represents the 1D ECA evolution with each line being the state of the automaton at a given timestep, starting from a single cell set to 1. Cells which are in state 1 are represented in black and cells in state 0 are represented in white. Time increases downwards. Figure 4.1b represents the compressed length of the 256 ECA rules, with different marker and colors corresponding to the obtained KMeans clusters.</figcaption></figure>

only approach it in a given context. As an example, compressing text data (a task often
performed with gzip in practice) is much more efficient with a language model that can
assign a very low probability to non grammatically correct sentences. The Kolmogorov
complexity (Kolmogorov 1968) of a cellular automaton is upper bounded by a value that
is independent of the chosen rule, as it is entirely determined by the transition table,
the grid size, initial configuration and number of steps.

**4.4 Predictor-based metric**

One obvious limit of using compression length as a proxy for complexity is the
fact that interesting systems mostly have intermediate compressed length. Compressed
length increases with the amount of disorder in the string being compressed. Therefore,
extreme lengths correspond either to systems that do not increase in complexity on
the lower end of the spectrum, or systems that produce a maximal amount of disorder
on the higher end. Neither of them have the potential to create interesting behavior
and increase in complexity. Intermediate values of compressed length are also hard to
interpret, since average lengths might either correspond to interesting rules or slowly
growing disordered systems.
To cope with this limitation, one should also take into account the dynamics of complexity, that is, how the system builds on its complexity at a given time as it keeps evolving, while retaining some of the structures it had acquired earlier. Compression leverages the amount of repetitions in inputs to further compress, and this may also be used as an estimate of structure overlap, as explained in the following section.

### 4.4.1 Joint compression

As a way to both measure the complexity and the amount of overlap between two automata states apart in time, we define a joint compressed length metric for a delay $τ$ as

$$C'(S^{(T+\tau)}, S^{(T)}) = C(S^{(T)} || S^{(T+\tau)}) \qquad (4.2)$$

where $||$ represents the concatenation operator. This quantity is simply the compressed length of a pair of global states — defined at the beginning of 4.3, represented by the letter $S$ — at two timesteps separated by a delay $τ$. In 1D, concatenation means chaining the two string representations before compressing, and in 2D we can chain two flattened representations of the 2D grid. This introduces several issues which we discuss in Section 4.4.2.

To quantify the amount of overlap between the two global states, we can compute the ratio of this joint compressed length with the sum of the two compressed lengths $C(S^t)$ and $C(S^{t-τ})$, thereby forming the joint compression score

$$\mu = \frac{C(S^t) + C(S^{t-\tau})}{C'(S^t, S^{t-\tau})} \qquad (4.3)$$

defined for an automaton $S$, time $t$ and delay $τ$.

This metric is based on the intuition that if patterns occur at step $T - τ$ of the automaton's evolution and are still present at step $T$, the joint compressed length will be lower than the sum of the two compressed lengths. The idea is illustrated in Figure 4.2, where it is pointed out that a stable moving structure (sometimes called *glider* or *spaceships* in Game of Life) will yield lower joint compressed lengths. This also applies to structures that self-replicate, grow from a stable seed or maintain the presence of
some sub-structures. Bigger structures yield a higher compression gain.

<figure><img src="image_43.png" /><figcaption>**Figure 4.2** — Joint compression method illustration. If a structure persists through time, this will decrease the joint compressed length compared to the sum of compressed lengths. A persistent structure is circled in red.</figcaption></figure>

Joint compression alone is not sufficient since it only selects rules that either behave like identity or shift input because they maximize the conservation of structures through time — as illustrated in Figure 4.3a. However, one may combine the joint compression score with another complexity measure to only select rules that exhibit some disorder, or growth in complexity — as Figure 4.3b shows (the condition here was a threshold on the difference of compressed length between initial and final states).

<figure><img src="image_44.png" /><figcaption>(a) Highest joint compression score among (b) With condition on compressed length the 256 ECA. increase.</figcaption></figure>

<caption>**Figure 4.3** — Comparison of the raw joint compression score and the addition of a complexity increase condition. The high overlap in structures is not enough to get interesting rules a shown in 4.3a, but the addition of a complexity threshold allows to retrieve rules with complex but still structured behavior, as shown in 4.3b. Figures are from the same slice of 60 cells over 30 timesteps taken from larger automata with random initial states. The top row corresponds to *t* = 0 and time increases downwards.</caption>

### 4.4.2 Count-based predictor

A major issue with the joint compression metric is the fact that it is designed to compress a linear stream of data. This is not ideal when considering higher dimensional automata. Larger sets of transformations have to be considered such as translations, rotations, flips, etc. Theoretically this should not be a problem for a good enough linear compression algorithm, but hardware and software limitations make it impractical to
work with existing algorithms on higher dimensional structures — with e.g DEFLATE's upper limit on dictionary size.

These higher dimensional automata might be better at generating complex dynamics, and the large size of their rule spaces makes it a challenge to explore. There has been at least one attempt to deal with these higher dimensional systems (Zenil, Soler-Toscano, et al. 2015) that lacks the scalability to work with large inputs.

An alternative to the linear compression-based method presented above would be to use compressors optimized for n-dimensional data (e.g. PNG compression for 2D automata) to take advantage of spatial correlation for compressing. However, these compressors are rare for higher dimensional data, and are usually optimized for one type of input — e.g. images with PNG.

Another way to tackle the problem is to use a prediction based approach to compression. Similarly to methods described in (Schmidhuber and Heil 1996) and one of the first steps of the PAQ compression algorithm (Mahoney 2000), we learn a statistical model of input data to predict the content of a cell given its immediate predecessors. For compression, this is often followed by an encoding step — using Huffman or arithmetic coding — that encodes data which contains the least information (least “surprising” data) with the most compact representation. This approach can also be related to the texture synthesis method described in (Efros and Leung 1999), where the authors learn a non parametric model to predict the next pixel of a texture given a previously synthesized neighborhood. Additionally, because we do not need the operation to be reversible as in regular compression, it is not necessary to limit the prediction model to making prediction with predecessors only.

For a global state $S = (c_1, ..., c_i, ..., c_n)$, the neighborhood of cell $i$ with radius $r$, denoted $n_{r,i}$ is defined as the tuple $n_{r,i} = (c_{i-r}, ..., c_{i-1}, c_{i+1}, ..., c_{i+r})$ — without the middle cell. The goal of this method is to estimate the conditional probability distribution $p(s|n_r)$ of the middle states at timestep $T$ given a neighborhood of radius $r$. Assuming cell states given their neighborhood can be modeled by mutually independent random variables, the log-probability of global state $S^{(T)}$ is written

$$ \log p(S^{(T)}) = \log \prod_{i=1}^{N} p(c_i|n_{r,i}) = \sum_{i=1}^{N} \log p(c_i|n_{r,i}) \quad (4.4) $$
If the automaton has a very ordered behavior, a model will predict with high confidence the state of the middle cell given a particular neighborhood. On the other hand, in the presence of maximal disorder, the middle cell will have an equal probability of being in every state no matter the neighborhood. In the latter case, a predictive model minimizing $-\log p(S^{(T)})$ would yield a high negative log-probability.

A simple possible predictor for such purpose is a large lookup table that maps all visited neighborhoods to a probability distribution over the states that the middle cell can be in. State distributions for each neighborhood are obtained by measuring the frequency of cell states given some observed neighborhoods. We denote by $\Lambda$ this lookup table, defined for a window of radius $r$, which maps all possible neighborhoods of size $2r + 1$ (ignoring the middle cell) to a set of probabilities $p$ over the possible states $\{s_1, ..., s_n\}$, and $p$ can be written $[p_{s_1}, p_{s_2}, ..., p_{s_n}]$. $\Lambda$ is defined by

$$ \Lambda: \{s_1, ..., s_n\}^{2r} \rightarrow \Delta_n \qquad (4.5) \\ n_{r,i} \mapsto p $$

where $\Delta_n$ denotes the probability simplex in dimension $n$.

To measure the uncertainty of that predictor, we can compute the cross-entropy loss between the data distribution it was trained on and its output. We compute the log probability of the observed data given the model, or the quantity

$$ L = -\frac{1}{N} \sum_{i=1}^{N} \sum_{k=1}^{n} 1_{\{s_k\}}(c_i) \log \Lambda(n_{r,i})_{s_k} \qquad (4.6) $$

where $1_{\{s_k\}}$ denotes the indicator function of the singleton set $\{s_k\}$. An illustration of the counting process is represented in Figure 4.4. The quantity $L$ is minimal when the $\Lambda(n_{r,i})_{s_k}$ is always equal to one, which means that the state of every cell is entirely determined by its neighborhood.

We apply this metric to all 256 ECA, with a window radius of size 3 (the 6 closest neighbors are used for prediction), and the same settings as for Figure 4.1b. Cross-entropy loss of the lookup table gives the results of Figure 4.5a. Colors are the same as in Figure 4.1b for comparison purposes.

We note the similarity between this plot and the one from Figure 4.1b, with a roughly
<figure><img src="image_45.png" /><figcaption>**Figure 4.4** – Count-based predictor method for a radius *r* = 1. A frequency lookup table is computed from the global state at time *T* by considering all neighborhoods with radius *r* = 1 (3 consecutive cells but ignoring the middle cell). Cross-entropy with the automaton at time *T* quantifies the overall complexity. This can be compared to the cross-entropy at time *T* + *t* for the amount of overlap.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Figure 4.5** – Average cross entropy loss for the two predictor-based methods on all 256 ECA. Rules are separated in several clusters. The count-based predictor (left plot) and neural network-based predictor (right plot) were applied with a neighborhood radius *r* = 1 and 3.</figcaption></figure>
equivalent resulting classification of ECA rules, with the exception of rules with low score. Rules that produce highly disordered patterns are on top of the plot whereas the very simply behaving rules are at the bottom. This indicates coherence between the two metrics.

### 4.4.3 Neural network based predictor

<figure><img src="image_47.png" /><figcaption>**Figure 4.6** – Neural network architecture for predicting a central cell given its neighbors. Output probabilities are defined for all possible states of the central cell.</figcaption></figure>

The frequency-based predictor described above still has limitations:

- It does not take into account any redundancy in the input which may lead to suboptimal predictions (in a CA, very similar positions might have similar center cell state distribution, e.g. a glider in Game of Life should be recognized by the model no matter the rest of the neighborhood).

- For the same reasons, when considering large window sizes, the number of possible neighborhood configurations grows much larger than the number of observed ones, leading to an input sparsity problem.

More sophisticated models can cope with above limitations by dealing with high-dimensional inputs without sparsity problems, and taking into account redundancy of inputs and potential interactions between states for prediction.

We measure the cross-entropy loss of this simple model on the training set after a standard learning procedure, which is the same for all rules. The procedure is applied to a one-hidden-layer neural network with a fixed hidden-layer size. We use a ReLU nonlinearity for the hidden layer and a softmax to obtain the output probabilities.

For *n* possible states $s_1, ..., s_n$, a cell in state $s_k$ is represented as a vector of 0s of
size $n$ with a 1 in position $k$. The input to the network is the concatenation of these cell vectors for all cells in the neighborhood. The output of the network is a vector of size $n$ with the output probability for each state.

Gradient updates are computed during training to minimize the cross-entropy loss between outputs and target examples. For a timestep $T$, we use the training procedure in order to minimize with respect to $\theta$ the following quantity,

$$L_{\theta}^{(T)} = -\frac{1}{N} \sum_{i=1}^{N} \sum_{k=1}^{n} \mathbb{1}_{\{s_k\}} (c_{i}^{(T)}) \log [f_{\theta}(n_{r,i}^{(T)})] \quad (4.7)$$

where the neural network depending on parameter $\theta$ is denoted $f_{\theta}$, and $n_{r,i}^{(T)}$, the neighborhood of cell $i$ with radius $r$ at time $T$ is defined in the same way as in eq. equation 4.6. Loss is calculated with respect to the testing set at time $T + \tau$ by computing the same quantity at this subsequent timestep.

The training procedure is selected to achieve reasonable convergence of the loss for the tested examples. It must be well defined and stay the same to allow for comparison of the results across several rules. The score in timestep $T$ for a delay $\tau$ is computed with the following formula

$$\mu_{\tau} = \frac{L^{(T)}}{L^{(T+\tau)}} \qquad (4.8)$$

where $L^{(T+\tau)}$ is the log probability of the automaton state at timestep $T + \tau$ (defined in eq. equation 4.7) according to a model with parameters learned during training at timestep $T$ and $L^{(T)}$ is the same as in eq. equation 4.7. The value $\mu_{\tau}$ will be lower for easily “learnable” global states that do not translate well in future steps — they create more complexity or disorder — thereby discarding slowly growing very disordered structures. Higher values of $\mu_{\tau}$ correspond to automata that have a disordered global state at time $T$ that can be transposed to future timesteps relatively well. Those rules will tend to have interesting spatial properties — not trivially simple but not completely disordered because the model transposes well — as well as a large amount of overlap between a given step and the future ones, indicating persistence of the spatial properties from one state to another. We also selected the metric among other quantities computed from $L^{(T)}$ and $L^{(T+\tau)}$ because it gave the best score in our experimental datasets.
## 4.5 Experiments

We carried out several experiments on a dataset of 500 randomly generated 3 states ($n = 3$) rules with radius $r = 1$. Those rules were manually annotated for interestingness, defined as the presence of visually detectable non-trivial structures. The annotation was carried out by five persons, using a majority rule to assign a binary label to each rule. The dataset contains 46 rules labeled as interesting and 454 uninteresting rules. Ranking those rules with the metrics introduced above allows us to study the influence of parameters and the adequacy between interestingness as we understand it and what the metric measures.

The task of finding interesting rules can be framed as either a classification problem or a ranking problem with respect to the score we compute on the dataset. The performance of our metric can be measured with the usual evaluation metrics used on these problems, and, notably, the average precision (AP) of the resulting classifier.

The average precision scores for the neural network and the count-based methods for time windows of 5, 50, and 300 timesteps are given in Table 4.1. Scores were calculated in automata of size 256 × 256 cells, ran for 1000 timesteps ($T + \tau = 1000$). Scores were computed for radii ranging from 1 cell (8 nearest neighbors) to 5 cells (120 neighbors), with a one-layer neural network containing 10 hidden units trained for 30 epochs with batches of 8 examples. The best AP for each time window is shown in bold. Results for the frequency lookup table predictor are only shown for $r = 1, 2$ because of the sparsity issues with the lookup table from $r = 2$ and above, making it unpractical to use the table — 3²⁴ possible entries for the lookup table with $r = 2$ against only 256² observed states.

From these experiments, larger radii appear to perform slightly better, although not in a radical way. Since the number of neighbors scales with the square of the radius, reasonably small radii might be a good trade-off between performance and computational cost of the metric.

We also study the performance of our metrics — lookup table and neural network-based — as inputs of a binary classifier against two simple baselines on a random 70/30 split of our dataset. The first baseline classifies all examples as negative. The second baseline is based on compressed length as defined in (Zenil 2010) and computed by
<table><caption>**Table 4.1** – Experimental results - AP scores</caption><thead><tr><th>Neural network</th><th>r = 1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>5 steps</td><td>0.387</td><td>0.448</td><td><strong>0.541</strong></td><td>0.525</td><td>0.534</td></tr><tr><td>50 steps</td><td>0.377</td><td>0.433</td><td>0.517</td><td>0.491</td><td><strong>0.542</strong></td></tr><tr><td>300 steps</td><td>0.358</td><td>0.454</td><td>0.488</td><td><strong>0.527</strong></td><td>0.525</td></tr></tbody><tbody><tr><th>Lookup table</th><th>r = 1</th><th>2</th><th></th><th></th><th></th></tr><tr><td>5 steps</td><td>0.092</td><td>0.070</td><td></td><td></td><td></td></tr><tr><td>50 steps</td><td>0.102</td><td>0.070</td><td></td><td></td><td></td></tr><tr><td>300 steps</td><td>0.093</td><td>0.069</td><td></td><td></td><td></td></tr></tbody></table>

This table shows the average precision (AP) scores obtained on the dataset of section 4.5 with the neural network-based and lookup table-based methods. Results are shown for delays $\tau = 5, 50, 300$ and several radii values $r$.

<table><caption>**Table 4.2** – Experimental results - Accuracy of each metric of complexity when used to classify which automatons do evolve interestingly, compared against the trivial all-negative baseline and the compressed length metric (Zenil 2010).</caption><thead><tr><th>Metric</th><th>Baseline</th><th>Compressed length (Zenil 2010)</th><th>Lookup Table</th><th>Neural network</th></tr></thead><tbody><tr><td>Accuracy</td><td>0.90</td><td>0.913</td><td>0.926</td><td><strong>0.953</strong></td></tr></tbody></table>

choosing a pair of thresholds that minimize mean square error when classifying examples in between as positive — this is based on the observation made in Section 4.3 that interesting rules have intermediate compressed lengths. Results are in Table 4.2 where only the best radius is shown. The lookup table performs better than the baselines, but the neural network gives the best score.

Above experiments demonstrate the capability of our proposed metric to match a subjective notion of interestingness of our labeling. For instance, the top 5 and top 10 scoring rules of the best performing configuration ($r = 3, \tau = 5$) are all labeled as interesting, and top 100 scores contain 41 of the 46 rules labeled as interesting.

## 4.6 Discussion

In this section, we discuss the results obtained by using the metric of equation 4.8 and the way they can be interpreted.
**One dimensional cellular automata** By applying the metric on the same example as before, we again obtain a plot with a rule classification that matches a visual appreciation of complexity of 1D CA. Results are shown on Figure 4.5b. Similarly to the previous cases, rules we might label as interesting are unlikely to be either at the top or bottom of the plot.

**Two dimensional cellular automata** Simulations conducted with 2D CA used grids of size 256×256. Automata were ran for 1000 steps (the metric is measured with respect to the reference time $T = 1000$). Rules are defined with a table of transitions from all possible neighborhood configurations with radius $r = 1$ (3×3 squares) to a new state for the central cell. Unbiased sampling of rules, obtained by uniformly sampling the resulting state for each transition independently, overwhelmingly produces rules with a similar amount of transitions towards each state and fails to produce rules without completely disordered behavior more than 99% of the time.

Therefore, we adopt a biased sampling strategy of the rules, selecting the proportion of transitions towards each state uniformly on the simplex — e.g for 3 states we might get the triple (0.1, 0.5, 0.4) and sample transitions according to these proportions. This parametrization can be related to Langton's lambda parameter (Langton 1990) that takes into account the proportion of transitions towards a transient (inactive) state and all the other states. We obtain approximately 10% interesting rules with this sampling as the proportions of our experimental dataset show.

Using the neural network-based complexity metric, we were able to find rules with interesting behavior among a very large set through random sampling. Some of these rules are shown here. Figure 4.8 displays three 2D rules that were selected manually upon visual inspection among the 20 highest scoring for metric $\mu_{50}$ (defined in eq. equation 4.8) of a sample of 1700 randomly generated 2-states 3 by 3 neighborhood rules. For comparison, Conway's Game of Life rule (GoL) ranks in the top 1% of the 2500 rules mentioned above for runs that do not end in a static global state. We observe that spontaneous glider formation events appear to be captured by our metric. Although gliders in cellular automata are a simple process that can be created manually, detection of their spontaneous emergence within a random search setting is a first step towards
<figure><img src="image_48.png" /><figcaption>**Figure 4.7** – Rules with 3 states that have spontaneously occurring glider structures. The gliders are the small structures that are outside of the center disordered zone. Some of them move along the diagonals while some others follow horizontal or vertical paths. Note that some repeating patterns occur also in the more disordered center zone.</figcaption></figure>

finding more complex macro structures that can emerge out of simple components. Rules with low scores are overwhelmingly of the disordered kind.

Figures 4.7, 4.9 and 4.10 show some three states rules that were selected through random sampling on the simplex with the neural network-based metric. They were selected among the 30 highest scoring rules out of 2500 randomly selected 3 states rules. All of their behaviors involve the growth and interaction of some small structures made of elementary cells.

All automata were initialized with a random disordered square of 20 by 20 cells in the center. In the Figures mentioned above, colors were normalized with the most common state set to blue. Figure 4.7 shows rules that spontaneously emit gliders that go through space in a direction until they interact with some other active part of the automaton. Figure 4.9 shows rules that generate small structures of between four and thirty cells that are relatively stable and interact with each other. These elementary components could be a basis for the spontaneous construction of more complex and bigger components. Figure 4.10 shows some other rules from this set of high ranking automata. They highlight the wide range of behaviors that can be obtained with these systems. Interesting rules from our search process can be found, along with other examples, in the form of animated GIFs *.

For some of these rules interesting patterns appear less frequently in smaller grids,

<footnote>* https://bit.ly/interesting_automata</footnote>
<figure><img src="image_49.png" /><figcaption>**Figure 4.8** – Spontaneous glider formation and evolution is observed for some high scoring 2 states rules. Each row corresponds to a rule, with a 50 timesteps difference between the two columns. Gliders are marked with a gray square. Runs were initialized with a small 20 by 20 disordered square (uniformly sampled among possible configuration) in the center simulated for up to 400 steps.</figcaption></figure>

indicating that the size of the space might impact the ability to generate complex macro-
structures. Increasing the size of the state space to very large grids might therefore make
it easier to generate very complex patterns.

**4.7 Conclusion**

In this chapter, we have proposed compression-inspired metrics for measuring a type
of complexity occurring in complex systems. We demonstrated its usefulness for selecting
CA rules that generate interesting emergent structures from very large sets of possible
rules. In higher dimensions where linear compression, as in gzip, is not sufficient to find
complex patterns, our metric is also useful.

We study 2 and 3 states automata here and we plan to investigate the effects of
<figure><img src="image_50.png" /><figcaption>**Figure 4.9** – Rules with 3 states that generate cell-like interacting structures. These patterns are either static or moving and can interact with one another to generate copies of themselves and other patterns. Note the very similar micro-structures that are repeated at several places in the space.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>**Figure 4.10** – Rules with surprising behaviors that are highly structured but complex. Those rules were selected among high-ranking rules for the neural-network based complexity metric. They all exhibit structurally non trivial behavior.</figcaption></figure>

additional states or larger neighborhoods on the ability to evolve more structures and obtain more interesting behaviors.

The dataset and code to reproduce our experiments and improvement on the reported results are published on GitHub†. The metrics we introduce in this chapter could be used to design organized systems of artificial developing organisms that grow in complexity through an evolutionary mechanism. A possible path toward such systems could start by creating an environment where computational resource allocation favors the fraction of subsystems that perform the best according to our measure of complexity.

The proposed metric is theoretically applicable to any complex system where a

<footnote>†. https://github.com/hugcis/evolving-structures-in-complex-systems</footnote>
notion of state of an elementary component and locality can be defined. With these
requirements satisfied, we can build a similar prediction model that uses information
about local neighbors to predict the state of a component, and thereby assess the
structural complexity of an input.

We believe that the capability of creating evolving systems out of such elementary
components and with few assumptions could be a step towards AGI. By devising ways
to guide this evolution in a direction that we find useful, we would be able to find an
efficient solution to hard problems while retaining the adaptability of the system. It
might be suitable to avoid over-specialization that can happen in systems designed to
solve a particular task — e.g. reinforcement learning algorithms that can play games, and
supervised learning — by staying away from any sort of objective function to optimize
and by leaving room for open-ended evolution.
# Visualization of computations in large-scale complex systems

Emergent processes in complex systems such as cellular automata can perform computations of increasing complexity, and could possibly lead to artificial evolution. Such a feat would require scaling up current simulation sizes to allow for enough computational capacity. Understanding complex computations happening in cellular automata and other systems capable of emergence poses many challenges, especially in large-scale systems. We propose methods for coarse-graining cellular automata based on frequency analysis of cell states, clustering and autoencoders. These innovative techniques facilitate the discovery of large-scale structure formation and complexity analysis in those systems. They emphasize interesting behaviors in elementary cellular automata while filtering out background patterns. Moreover, our methods reduce large 2D automata to smaller sizes and enable identifying systems that behave interestingly at multiple scales.

## 5.1 Introduction

Cellular automata (CA) have been extensively studied since the 1960s. Originally designed and studied to create artificial evolution from self-replication (Langton 1984; Von Neumann and Burks 1966), previously studied cellular automata simulations were often of relatively modest sizes. Only specific rules with repetitive or predictable dynamics such as John Conway's Game of Life (Gardner 1970) have been scaled up to larger grid sizes ($10^4 \times 10^4$ or more cells).

For complex phenomena such as artificial evolution to exist and be open-ended within
those simulated worlds, there needs to be sufficient “capacity” — a large enough state-space. In nature, complex and significantly different dynamics often arise from uniform laws at a smaller scale (Anderson 1972). It seems unlikely that such complex processes, like artificial evolution, could happen in too small CAs because higher order dynamics do not have enough capacity to emerge. However, several issues arise when scaling CAs to large sizes:

- Time complexity rapidly becomes a bottleneck. Updating a large number of cells is costly. Tricks such as caching of some of the computations can help, but do not always improve performance significantly (Gosper 1984).

- Memory complexity can also become an issue when dealing with numerous states, and especially grids in 3 dimensions and more. In that case, even the underlying rule of the system cannot be stored within reasonable memory capacity.

- Visual inspection of these large grids is infeasible. Studying CA complexity is rendered difficult by the highly variable nature of emergent processes. It is especially the case for large-scale systems.

When working with such large systems, it is less relevant to focus on the local behaviors at the single cell level. This is similar to other complex systems like the weather, in which behaviors of individual atoms in a cloud are irrelevant to large-scale air mass movements. Much richer behaviors can be observed from studying large patterns' formation and their evolution. This should also hold true for CAs; we further discuss this question in **Conclusion**.

In this chapter, we investigate techniques that can help us visualize large space-time diagrams of CAs. We demonstrate that simple clustering and coarse-graining techniques can be used in order to perceive structures that cannot emerge on smaller grids. This is also useful for disordered cellular automata with hidden structures as is the case for the elementary cellular automaton rule 18, illustrated in Figure 5.1 — more details in **Results**.

Reducing large grids to smaller sizes while preserving interesting behaviors such as pattern formation is essential to apply to these CAs complexity metrics designed to work on modestly sized grids (Grassberger 1986; Soler-Toscano et al. 2014; Zenil 2010; Zenil, Soler-Toscano, et al. 2015). Common metrics of complexity are often limited by the
<figure><img src="image_52.png" /><figcaption>**Figure 5.1** – Hidden structures in rule 18 are uncovered by filtering the space-time diagram with our frequency histogram-based method. (a) shows 300 timesteps of a randomly initialized rule 18 simulation. Notice the complex structures made visible in (b) with our method.</figcaption></figure>

number of components in the systems (number of cells in a CA grid, timesteps, etc.) or
may not be effective when small-scale patterns are less relevant than large-scale ones.

**5.2 Related work**

Previous work on coarse-graining cellular automata focused either on conserving the main computational properties of CA rules through exact coarse-graining or on filtering interesting behaviors without reducing the amount of computations. Our work both highlights interesting behaviors and compresses the representation, which we argue are necessary to study complexity in large cellular automata.

**5.2.1 Coarse-graining in cellular automata**

Coarse-graining is an approximation procedure used to speed up computations in
systems made of many components. It originated in (Levitt and Warshel 1975) and
is now widely used in physics to model complex systems at various granularity levels
and is successful at modeling bio-molecules (Ingólfsson et al. 2014; Kmiecik et al. 2016;
Potoyan, Savelyev, and Papoian 2013).

The exact coarse-graining of elementary cellular automata (ECA) has been investi-
gated extensively in (Israeli and Goldenfeld 2004, 2006). Authors found ways of rewriting