We develop a theory characterizing the fundamental capability of deep neural networks to learn, from evolution traces, the logical rules governing the behavior
of cellular automata (CA). This is accomplished by first establishing a novel connection between CA and Łukasiewicz propositional logic. While binary CA have been known for decades to essentially perform operations in Boolean logic, no such relationship exists for general CA. We demonstrate that many-valued (MV) logic, specifically Łukasiewicz propositional logic, constitutes a suitable language for characterizing general CA as logical machines.
This is done by interpolating CA transition functions to continuous piecewise linear functions, which, by virtue of the McNaughton theorem, yield formulae in MV logic characterizing the CA. Recognizing that deep rectified linear unit (ReLU) networks realize continuous piecewise linear functions, it follows that these formulae are naturally extracted from CA evolution traces by deep ReLU networks. A corresponding algorithm together with a software implementation is provided.

Finally, we show that the dynamical behavior of CA can be

realized by recurrent neural networks.

§ INTRODUCTION

Neural networks were originally introduced as a computational model whose functionality imitates that of the human brain <cit.>. Machine learning based on neural networks has achieved state-of-the-art results in numerous applications, such as pattern recognition <cit.>, game intelligence <cit.>, or protein structure prediction <cit.>. Informally speaking, a neural network consists of layers of nodes that are connected by weighted edges. In practice, the network topology and weights are learned through training on data either in a supervised <cit.> or an unsupervised <cit.> manner.

An ability ingrained in the human brain is to draw inferences from data it is presented with and apply what has been learned to new problems. In machine learning parlance, this aspect is usually referred to as generalization or extrapolation. For example, given the geometric sequence of integers $1,2,4,8,\ldots$, one would conclude that the rule generating the sequence is: The next number is obtained by multiplying the present number by $2$. Based on this extracted logical rule, the entire sequence can now be calculated. Such logical reasoning tasks, which appear natural to humans, may be challenging to learn for neural networks. Moreover, even if the network had learned, e.g., the mechanism generating the geometric sequence above, it would have the corresponding logical rule encoded in its topology and weights and it is unclear how to extract the rule and present it in a manner accessible to humans.

In this paper, we report an attempt at building a theory characterizing the fundamental capability of neural networks to learn logical rules from data. This immediately leads to the question, “What is the logical structure behind data?” Here, we shall make a first step towards developing the corresponding research program. Concretely, we consider data generated by cellular automata[We will abbreviate both the singular “cellular automaton” and the plural form “cellular automata” as “CA”.] (CA). Abstractly speaking, a CA is a discrete dynamical system evolving on a regular lattice whose points take values in a finite set. Starting from an initial configuration, all lattice points change their states at synchronous discrete time steps according to a transition function that takes the present values of the lattice point under consideration and its neighbors as inputs. The evolution of the lattice point states over time furnishes the data sequence we shall be interested in. Now, assume that the underlying CA is binary, i.e., the state set is given by $\{0,1\}$. In this case the transition function governing the CA evolution is a Boolean function <cit.>, which can be viewed

as determining the logical structure behind the data sequence generated by the CA.

The question we shall ask is whether one can extract the logical rule governing a general (i.e., not only binary) CA by training a neural network on data generated by the CA. While the connection between binary CA and Boolean logic has been known for decades <cit.>, <cit.>, to the best of our knowledge no such relationship has been reported in the literature for CA with arbitrary state sets. One of the main conceptual contributions of this paper is to demonstrate that many-valued[With slight abuse of terminology, we shall use the term MV logic to refer to propositional logic.] (MV) logic is a suitable language for interpreting general CA as logical machines. We then show that

all possible transition functions can be realized by deep ReLU networks[ReLU stands for the Rectified Linear Unit nonlinearity, defined as $x\mapsto \max\{0,x\}$.], which, in turn, are found to naturally express statements in MV logic.

The dynamical system component of CA is demonstrated to be realizable by recurrent neural networks (RNNs). Finally, we propose and analyze a procedure for extracting the logical formulae behind CA transition functions from neural networks trained on corresponding CA evolution data.

Notation: $\mathbbm{1}_{ \{\cdot \}}$ denotes the truth function which takes the value $1$ if the statement inside $\{\cdot\}$ is true and equals $0$ otherwise. $0_N$ stands for the $N$-dimensional column vector with all entries equal to $0$. $\mathbb{N}_0:= \mathbb{N} \cup \{0\}$ and $|A|$ stands for the cardinality of the set $A$. $||\cdot||_1$ is the $\ell_1$-norm.

§.§ Cellular automata

CA were invented in the 1940s by von Neumann <cit.> and Ulam <cit.> in an effort to build models that are capable of universal computation and self-reproduction. Von Neumann's conceptualization emphasized the aspect of self-reproduction, while Ulam suggested the use of finite state machines on two-dimensional lattices. A widely known CA is the two-dimensional Game of Life devised by Conway <cit.>.
Despite the simplicity of its rules, the Game of Life exhibits remarkable behavioral complexity and has therefore attracted widespread and long-standing interest. We begin by briefly reviewing the Game of Life.

Consider an infinite two-dimensional grid of square cells centered on the points of an underlying lattice (symbolized by dashed lines in fig:neighborhood). Initially, each cell (or equivalently lattice point) is in one of two possible states, namely “live” or “dead”. Each cell has eight neighbors, taking into account horizontal, vertical, and diagonal directions (see fig:neighborhood). For a given cell, we refer to the set made up of the cell itself and its neighbors as the neighborhood of the cell. The cells change their states synchronously at discrete time steps, all following the same rule given by:

* Every live cell with two or three live neighbors stays live; other live cells turn into dead cells.

* Every dead cell with three live neighbors becomes live; other dead cells stay dead.

[x=1cm, y=1cm, >=stealth, scale=1, transform shape]

at (0.75, 0.75) $z$;

[dashed] (-0.25,0.25) – (1.75,0.25);

[dashed] (-0.25,0.75) – (1.75,0.75);

[dashed] (-0.25,1.25) – (1.75,1.25);

[dashed] (0.25, -0.25) – (0.25, 1.75);
[dashed] (0.75, -0.25) – (0.75, 1.75);
[dashed] (1.25, -0.25) – (1.25, 1.75);

[line width=0.5mm] (0,0) – (1.5,0);
[line width=0.5mm] (0,0.5) – (1.5,0.5);
[line width=0.5mm] (1.5,0) – (1.5,1.5);
[line width=0.5mm] (0,1) – (1.5,1);
[line width=0.5mm] (1.5,1.5) – (0,1.5);
[line width=0.5mm] (0,1.5) – (0,0);
[line width=0.5mm] (0.5,1.5) – (0.5,0);
[line width=0.5mm] (1,1.5) – (1,0);

[x=1cm, y=1cm, >=stealth, scale=1, transform shape]

at (0.75, 0.75) $z$;

[dashed] (-0.25,0.25) – (1.75,0.25);

[dashed] (-0.25,0.75) – (1.75,0.75);

[dashed] (-0.25,1.25) – (1.75,1.25);

[dashed] (0.25, -0.25) – (0.25, 1.75);
[dashed] (0.75, -0.25) – (0.75, 1.75);
[dashed] (1.25, -0.25) – (1.25, 1.75);

[line width=0.5mm] (0,0) – (1.5,0);
[line width=0.5mm] (0,0.5) – (1.5,0.5);
[line width=0.5mm] (1.5,0) – (1.5,1.5);
[line width=0.5mm] (0,1) – (1.5,1);
[line width=0.5mm] (1.5,1.5) – (0,1.5);
[line width=0.5mm] (0,1.5) – (0,0);
[line width=0.5mm] (0.5,1.5) – (0.5,0);
[line width=0.5mm] (1,1.5) – (1,0);

Cell $z$ and its neighborhood.

This process repeats at discrete time steps, forming an evolution of the system. We refer to fig:evolution for an example illustration. As already mentioned, despite the simplicity of the transition rule, the evolutions induced by different initial configurations exhibit rich and complex long-term behavior. It is shown in <cit.>, <cit.>, <cit.>, for example, that given an initial configuration with a finite number of live cells, it is an undecidable problem whether all these live cells would eventually die out. We refer to <cit.> for an in-depth discussion of the behavioral patterns of the Game of Life and proceed to formally define CA.

[x=1cm, y=1cm, >=stealth, scale=1, transform shape]

[] (0,1) – (0,5);   [dashed] (0.5,0) – (0.5,6);
[] (1,1) – (1,5);[dashed] (1.5,0) – (1.5,6);
[] (2,1) – (2,5); [dashed] (2.5,0) – (2.5,6);
[] (3,1) – (3,5); [dashed] (3.5,0) – (3.5,6);
[] (4,1) – (4,5); [dashed] (4.5,0) – (4.5,6);
[] (5,1) – (5,5);
[] (0,1) – (5,1);   [dashed] (-1,1.5) – (6,1.5);
[] (0,2) – (5,2);  [dashed] (-1,2.5) – (6,2.5);
[] (0,3) – (5,3);  [dashed] (-1,3.5) – (6,3.5);
[] (0,4) – (5,4); [dashed] (-1,4.5) – (6,4.5);
[] (0,5) – (5,5);
[fill=black]  (1,4) – (1,5) – (2,5) – (2,4) – cycle;
[fill=black]  (2,3) – (2,4) – (3,4) – (3,3) – cycle;
[fill=black]  (0,2) – (0,3) – (3,3) – (3,2) – cycle;
[-stealth, line width=2mm] (7,3) – (8,3) node[anchor=west] ;

[] (10,1) – (10,5);[dashed] (10.5,0) – (10.5,6);
[] (11,1) – (11,5); [dashed] (11.5,0) – (11.5,6);
[] (12,1) – (12,5); [dashed] (12.5,0) – (12.5,6);
[] (13,1) – (13,5); [dashed] (13.5,0) – (13.5,6);
[] (14,1) – (14,5); [dashed] (14.5,0) – (14.5,6);
[] (15,1) – (15,5);
[] (10,1) – (15,1);   [dashed] (9,1.5) – (16,1.5);
[] (10,2) – (15,2);  [dashed] (9,2.5) – (16,2.5);
[] (10,3) – (15,3); [dashed] (9,3.5) – (16,3.5);
[] (10,4) – (15,4);[dashed] (9,4.5) – (16,4.5);
[] (10,5) – (15,5);
[fill=black]  (12,3) – (12,4) – (13,4) – (13,3) – cycle;
[fill=black]  (10,3) – (10,4) – (11,4) – (11,3) – cycle;
[fill=black]  (11,2) – (11,3) – (13,3) – (13,2) – cycle;
[fill=black]  (11,1) – (11,2) – (12,2) – (12,1) – cycle;
[-stealth, line width=2mm] (17,3) – (18,3) node[anchor=west] ;

[] (20,1) – (20,5); [dashed] (20.5,0) – (20.5,6);
[] (21,1) – (21,5); [dashed] (21.5,0) – (21.5,6);
[] (22,1) – (22,5); [dashed] (22.5,0) – (22.5,6);
[] (23,1) – (23,5);[dashed] (23.5,0) – (23.5,6);
[] (24,1) – (24,5); [dashed] (24.5,0) – (24.5,6);
[] (25,1) – (25,5);
[] (20,1) – (25,1);  [dashed] (19,1.5) – (26,1.5);
[] (20,2) – (25,2);   [dashed] (19,2.5) – (26,2.5);
[] (20,3) – (25,3);   [dashed] (19,3.5) – (26,3.5);
[] (20,4) – (25,4); [dashed] (19,4.5) – (26,4.5);
[] (20,5) – (25,5);
[fill=black]  (20,2) – (20,3) – (21,3) – (21,2) – cycle;
[fill=black]  (21,1) – (21,2) – (22,2) – (22,1) – cycle;
[fill=black]  (22,1) – (22,4) – (23,4) – (23,1) – cycle;
[-stealth, line width=2mm] (27,3) – (28,3) node[anchor=west] ;

[] (30,1) – (30,5);[dashed] (30.5,0) – (30.5,6);
[] (31,1) – (31,5); [dashed] (31.5,0) – (31.5,6);
[] (32,1) – (32,5); [dashed] (32.5,0) – (32.5,6);
[] (33,1) – (33,5); [dashed] (33.5,0) – (33.5,6);
[] (34,1) – (34,5); [dashed] (34.5,0) – (34.5,6);
[] (35,1) – (35,5);
[] (30,1) – (35,1);  [dashed] (29,1.5) – (36,1.5);
[] (30,2) – (35,2);  [dashed] (29,2.5) – (36,2.5);
[] (30,3) – (35,3); [dashed] (29,3.5) – (36,3.5);
[] (30,4) – (35,4);[dashed] (29,4.5) – (36,4.5);
[] (30,5) – (35,5);
[fill=black]  (31,3) – (31,4) – (32,4) – (32,3) – cycle;
[fill=black]  (31,1) – (31,2) – (33,2) – (33,1) – cycle;
[fill=black]  (32,2) – (32,3) – (34,3) – (34,2) – cycle;
[-stealth, line width=2mm] (37,3) – (38,3) node[anchor=west] ;

[] (40,1) – (40,5); [dashed] (40.5,0) – (40.5,6);
[] (41,1) – (41,5);[dashed] (41.5,0) – (41.5,6);
[] (42,1) – (42,5); [dashed] (42.5,0) – (42.5,6);
[] (43,1) – (43,5);[dashed] (43.5,0) – (43.5,6);
[] (44,1) – (44,5);[dashed] (44.5,0) – (44.5,6);
[] (45,1) – (45,5);
[] (40,1) – (45,1);  [dashed] (39,1.5) – (46,1.5);
[] (40,2) – (45,2);  [dashed] (39,2.5) – (46,2.5);
[] (40,3) – (45,3); [dashed] (39,3.5) – (46,3.5);
[] (40,4) – (45,4);[dashed] (39,4.5) – (46,4.5);
[] (40,5) – (45,5);

[fill=black]  (42,3) – (42,4) – (43,4) – (43,3) – cycle;
[fill=black]  (43,2) – (43,3) – (44,3) – (44,2) – cycle;
[fill=black]  (41,1) – (41,2) – (44,2) – (44,1) – cycle;

[x=1cm, y=1cm, >=stealth, scale=1, transform shape]

[] (0,1) – (0,5);   [dashed] (0.5,0) – (0.5,6);
[] (1,1) – (1,5);[dashed] (1.5,0) – (1.5,6);
[] (2,1) – (2,5); [dashed] (2.5,0) – (2.5,6);
[] (3,1) – (3,5); [dashed] (3.5,0) – (3.5,6);
[] (4,1) – (4,5); [dashed] (4.5,0) – (4.5,6);
[] (5,1) – (5,5);
[] (0,1) – (5,1);   [dashed] (-1,1.5) – (6,1.5);
[] (0,2) – (5,2);  [dashed] (-1,2.5) – (6,2.5);
[] (0,3) – (5,3);  [dashed] (-1,3.5) – (6,3.5);
[] (0,4) – (5,4); [dashed] (-1,4.5) – (6,4.5);
[] (0,5) – (5,5);
[fill=black]  (1,4) – (1,5) – (2,5) – (2,4) – cycle;
[fill=black]  (2,3) – (2,4) – (3,4) – (3,3) – cycle;
[fill=black]  (0,2) – (0,3) – (3,3) – (3,2) – cycle;
[-stealth, line width=2mm] (7,3) – (8,3) node[anchor=west] ;

[] (10,1) – (10,5);[dashed] (10.5,0) – (10.5,6);
[] (11,1) – (11,5); [dashed] (11.5,0) – (11.5,6);
[] (12,1) – (12,5); [dashed] (12.5,0) – (12.5,6);
[] (13,1) – (13,5); [dashed] (13.5,0) – (13.5,6);
[] (14,1) – (14,5); [dashed] (14.5,0) – (14.5,6);
[] (15,1) – (15,5);
[] (10,1) – (15,1);   [dashed] (9,1.5) – (16,1.5);
[] (10,2) – (15,2);  [dashed] (9,2.5) – (16,2.5);
[] (10,3) – (15,3); [dashed] (9,3.5) – (16,3.5);
[] (10,4) – (15,4);[dashed] (9,4.5) – (16,4.5);
[] (10,5) – (15,5);
[fill=black]  (12,3) – (12,4) – (13,4) – (13,3) – cycle;
[fill=black]  (10,3) – (10,4) – (11,4) – (11,3) – cycle;
[fill=black]  (11,2) – (11,3) – (13,3) – (13,2) – cycle;
[fill=black]  (11,1) – (11,2) – (12,2) – (12,1) – cycle;
[-stealth, line width=2mm] (17,3) – (18,3) node[anchor=west] ;

[] (20,1) – (20,5); [dashed] (20.5,0) – (20.5,6);
[] (21,1) – (21,5); [dashed] (21.5,0) – (21.5,6);
[] (22,1) – (22,5); [dashed] (22.5,0) – (22.5,6);
[] (23,1) – (23,5);[dashed] (23.5,0) – (23.5,6);
[] (24,1) – (24,5); [dashed] (24.5,0) – (24.5,6);
[] (25,1) – (25,5);
[] (20,1) – (25,1);  [dashed] (19,1.5) – (26,1.5);
[] (20,2) – (25,2);   [dashed] (19,2.5) – (26,2.5);
[] (20,3) – (25,3);   [dashed] (19,3.5) – (26,3.5);
[] (20,4) – (25,4); [dashed] (19,4.5) – (26,4.5);
[] (20,5) – (25,5);
[fill=black]  (20,2) – (20,3) – (21,3) – (21,2) – cycle;
[fill=black]  (21,1) – (21,2) – (22,2) – (22,1) – cycle;
[fill=black]  (22,1) – (22,4) – (23,4) – (23,1) – cycle;
[-stealth, line width=2mm] (27,3) – (28,3) node[anchor=west] ;

[] (30,1) – (30,5);[dashed] (30.5,0) – (30.5,6);
[] (31,1) – (31,5); [dashed] (31.5,0) – (31.5,6);
[] (32,1) – (32,5); [dashed] (32.5,0) – (32.5,6);
[] (33,1) – (33,5); [dashed] (33.5,0) – (33.5,6);
[] (34,1) – (34,5); [dashed] (34.5,0) – (34.5,6);
[] (35,1) – (35,5);
[] (30,1) – (35,1);  [dashed] (29,1.5) – (36,1.5);
[] (30,2) – (35,2);  [dashed] (29,2.5) – (36,2.5);
[] (30,3) – (35,3); [dashed] (29,3.5) – (36,3.5);
[] (30,4) – (35,4);[dashed] (29,4.5) – (36,4.5);
[] (30,5) – (35,5);
[fill=black]  (31,3) – (31,4) – (32,4) – (32,3) – cycle;
[fill=black]  (31,1) – (31,2) – (33,2) – (33,1) – cycle;
[fill=black]  (32,2) – (32,3) – (34,3) – (34,2) – cycle;
[-stealth, line width=2mm] (37,3) – (38,3) node[anchor=west] ;

[] (40,1) – (40,5); [dashed] (40.5,0) – (40.5,6);
[] (41,1) – (41,5);[dashed] (41.5,0) – (41.5,6);
[] (42,1) – (42,5); [dashed] (42.5,0) – (42.5,6);
[] (43,1) – (43,5);[dashed] (43.5,0) – (43.5,6);
[] (44,1) – (44,5);[dashed] (44.5,0) – (44.5,6);
[] (45,1) – (45,5);
[] (40,1) – (45,1);  [dashed] (39,1.5) – (46,1.5);
[] (40,2) – (45,2);  [dashed] (39,2.5) – (46,2.5);
[] (40,3) – (45,3); [dashed] (39,3.5) – (46,3.5);
[] (40,4) – (45,4);[dashed] (39,4.5) – (46,4.5);
[] (40,5) – (45,5);

[fill=black]  (42,3) – (42,4) – (43,4) – (43,3) – cycle;
[fill=black]  (43,2) – (43,3) – (44,3) – (44,2) – cycle;
[fill=black]  (41,1) – (41,2) – (44,2) – (44,1) – cycle;

Steps of the evolution of the Game of Life. Black cells are live; white cells are dead.

Let $d, n,k \in \mathbb{N}, k\geq 2$. A cellular automaton is an ordered quadruple $(\mathbb{Z}^d, K, \mathcal{E}, f)$, where

* the $d$-dimensional lattice $\mathbb{Z}^d$ is referred to as the cellular space,

* $K = \{0,\frac{1}{k-1},\ldots \hspace{-0.02cm}, \frac{k-2}{k-1}, 1\}$ is the state set with $k$ states,

* $\mathcal{E}  = \{0_d, z_1,\ldots \hspace{-0.02cm}, z_{n-1}\}\subset \mathbb{Z}^d$ is the neighborhood set,

* $f: K^n \rightarrow K $ is the transition function.

The lattice dimension $d$ is called the cellular space dimension. Cells are centered on the lattice points of $\mathbb{Z}^d$. A configuration, denoted by $c$, over the cellular space assigns a state in the state set to each cell (or equivalently lattice point). For $z \in \mathbb{Z}^d$, denote by $c[z]$ its state under the configuration $c$. Let $\Ccal$ be the set of all possible configurations over the cellular space. The CA map $F: \Ccal \rightarrow \Ccal$ effecting the configuration evolution from one time step to the next is determined by the transition function $f$ according to
\begin{equation}\label{eq:globalmap}
F(c)[z] = f(c[z], c[z+z_1],\ldots \hspace{-0.02cm}, c[z+z_{n-1}]), \quad \text{ for all } z\in \mathbb{Z}^d.
\end{equation}

The CA map $F$ is often referred to as the “global mapping function”, whereas the transition function $f$ is called the “local mapping function”.

Note that our choice of the state set $K$ in defn:CA is without loss of generality, as for every general state set $K^*$, there is a bijection $h: K^* \rightarrow K$, with $|K^*| = |K|$. For example, for the state set $K^* = \{1,2,3\}$, $h: K^* \rightarrow \{0,1/2,1\}, h(x) = (x-1)/2$ is a bijection.

We continue with examples illustrating defn:CA.

Consider a CA with cellular space dimension $d =1$, neighborhood set $\mathcal{E}=\{0,-1\}\subset \mathbb{Z}$, state set $K = \{0,1\}$, and the transition function $f:K^2\rightarrow K$ as provided in tbl:example1. The CA map $F$ is given by

\begin{equation*}
F(c)[z] = f( c[z], c[z-1]), \quad \text{ for all } z\in \mathbb{Z}.
\end{equation*}

fig:example1 depicts one step of the CA evolution.

$x_0\,x_{-1}$    00     10     01     11
$f(x_0, x_{-1})$     0     0     0     1

tableTransition function $f$.

(-4,1) – (4,1) node[anchor= east] ;

(-4,0) – (4,0) node[anchor= east] ;

(-4,1) – (-4,0) node[anchor= east] ;
(-3,1) – (-3,0) node[anchor= east] ;
(-2,1) – (-2,0) node[anchor= east] ;
(-1,1) – (-1,0) node[anchor= east] ;
(0,1) – (0,0) node[anchor= east] ;
(1,1) – (1,0) node[anchor= east] ;
(2,1) – (2,0) node[anchor= east] ;
(3,1) – (3,0) node[anchor= east] ;
(4,1) – (4,0) node[anchor= east] ;
[fill=black]  (-4,0) – (-4,1) – (-3,1) – (-3,0) – cycle;
[fill=black]  (-3,0) – (-3,1) – (-2,1) – (-2,0) – cycle;
[fill=black]  (0,0) – (0,1) – (1,1) – (1,0) – cycle;
[fill=black]  (2,0) – (2,1) – (3,1) – (3,0) – cycle;
at (-4.5,0.5) $\cdots$;
at (4.5,0.5) $\cdots$;

[-stealth, line width=1mm] (0,-0.4) – (0,-1) node[anchor= east] ;

(-4,-1.4) – (4,-1.4) node[anchor= east] ;

(-4,-2.4) – (4,-2.4) node[anchor= east] ;

(-4,-1.4) – (-4,-2.4) node[anchor= east] ;
(-3,-1.4) – (-3,-2.4) node[anchor= east] ;
(-2,-1.4) – (-2,-2.4) node[anchor= east] ;
(-1,-1.4) – (-1,-2.4) node[anchor= east] ;
(0,-1.4) – (0,-2.4) node[anchor= east] ;
(1,-1.4) – (1,-2.4) node[anchor= east] ;
(2,-1.4) – (2,-2.4) node[anchor= east] ;
(3,-1.4) – (3,-2.4) node[anchor= east] ;
(4,-1.4) – (4,-2.4) node[anchor= east] ;

[fill=black]  (-3,-2.4) – (-3,-1.4) – (-2,-1.4) – (-2,-2.4) – cycle;

at (-4.5,-1.9) $\cdots$;
at (4.5,-1.9) $\cdots$;

[dashed] (-2.5,-0.3) – (-2.5,1.3) node[anchor= east] ;
[dashed] (-3.5,-0.3) – (-3.5,1.3) node[anchor= east] ;
[dashed] (-1.5,-0.3) – (-1.5,1.3) node[anchor= east] ;
[dashed] (-0.5,-0.3) – (-0.5,1.3) node[anchor= east] ;
[dashed] (0.5,-0.3) – (0.5,1.3) node[anchor= east] ;
[dashed] (1.5,-0.3) – (1.5,1.3) node[anchor= east] ;
[dashed] (2.5,-0.3) – (2.5,1.3) node[anchor= east] ;
[dashed] (3.5,-0.3) – (3.5,1.3) node[anchor= east] ;

[dashed] (3.5,-2.7) – (3.5,-1.1) node[anchor= east] ;
[dashed] (2.5,-2.7) – (2.5,-1.1) node[anchor= east] ;
[dashed] (1.5,-2.7) – (1.5,-1.1) node[anchor= east] ;
[dashed] (0.5,-2.7) – (0.5,-1.1) node[anchor= east] ;
[dashed] (-0.5,-2.7) – (-0.5,-1.1) node[anchor= east] ;
[dashed] (-1.5,-2.7) – (-1.5,-1.1) node[anchor= east] ;
[dashed] (-2.5,-2.7) – (-2.5,-1.1) node[anchor= east] ;
[dashed] (-3.5,-2.7) – (-3.5,-1.1) node[anchor= east] ;

(-4,1) – (4,1) node[anchor= east] ;

(-4,0) – (4,0) node[anchor= east] ;

(-4,1) – (-4,0) node[anchor= east] ;
(-3,1) – (-3,0) node[anchor= east] ;
(-2,1) – (-2,0) node[anchor= east] ;
(-1,1) – (-1,0) node[anchor= east] ;
(0,1) – (0,0) node[anchor= east] ;
(1,1) – (1,0) node[anchor= east] ;
(2,1) – (2,0) node[anchor= east] ;
(3,1) – (3,0) node[anchor= east] ;
(4,1) – (4,0) node[anchor= east] ;
[fill=black]  (-4,0) – (-4,1) – (-3,1) – (-3,0) – cycle;
[fill=black]  (-3,0) – (-3,1) – (-2,1) – (-2,0) – cycle;
[fill=black]  (0,0) – (0,1) – (1,1) – (1,0) – cycle;
[fill=black]  (2,0) – (2,1) – (3,1) – (3,0) – cycle;
at (-4.5,0.5) $\cdots$;
at (4.5,0.5) $\cdots$;

[-stealth, line width=1mm] (0,-0.4) – (0,-1) node[anchor= east] ;

(-4,-1.4) – (4,-1.4) node[anchor= east] ;

(-4,-2.4) – (4,-2.4) node[anchor= east] ;

(-4,-1.4) – (-4,-2.4) node[anchor= east] ;
(-3,-1.4) – (-3,-2.4) node[anchor= east] ;
(-2,-1.4) – (-2,-2.4) node[anchor= east] ;
(-1,-1.4) – (-1,-2.4) node[anchor= east] ;
(0,-1.4) – (0,-2.4) node[anchor= east] ;
(1,-1.4) – (1,-2.4) node[anchor= east] ;
(2,-1.4) – (2,-2.4) node[anchor= east] ;
(3,-1.4) – (3,-2.4) node[anchor= east] ;
(4,-1.4) – (4,-2.4) node[anchor= east] ;

[fill=black]  (-3,-2.4) – (-3,-1.4) – (-2,-1.4) – (-2,-2.4) – cycle;

at (-4.5,-1.9) $\cdots$;
at (4.5,-1.9) $\cdots$;

[dashed] (-2.5,-0.3) – (-2.5,1.3) node[anchor= east] ;
[dashed] (-3.5,-0.3) – (-3.5,1.3) node[anchor= east] ;
[dashed] (-1.5,-0.3) – (-1.5,1.3) node[anchor= east] ;
[dashed] (-0.5,-0.3) – (-0.5,1.3) node[anchor= east] ;
[dashed] (0.5,-0.3) – (0.5,1.3) node[anchor= east] ;
[dashed] (1.5,-0.3) – (1.5,1.3) node[anchor= east] ;
[dashed] (2.5,-0.3) – (2.5,1.3) node[anchor= east] ;
[dashed] (3.5,-0.3) – (3.5,1.3) node[anchor= east] ;

[dashed] (3.5,-2.7) – (3.5,-1.1) node[anchor= east] ;
[dashed] (2.5,-2.7) – (2.5,-1.1) node[anchor= east] ;
[dashed] (1.5,-2.7) – (1.5,-1.1) node[anchor= east] ;
[dashed] (0.5,-2.7) – (0.5,-1.1) node[anchor= east] ;
[dashed] (-0.5,-2.7) – (-0.5,-1.1) node[anchor= east] ;
[dashed] (-1.5,-2.7) – (-1.5,-1.1) node[anchor= east] ;
[dashed] (-2.5,-2.7) – (-2.5,-1.1) node[anchor= east] ;
[dashed] (-3.5,-2.7) – (-3.5,-1.1) node[anchor= east] ;

The state $1$ is indicated by black cells, the state $0$ by white cells.

[Game of Life]
Consider a CA with cellular space dimension $d =2$, state set $K = \{0,1\}$, and neighborhood set
$\mathcal{E} = \{0_2,z_1,\ldots \hspace{-0.02cm},\\z_8\} \subset \mathbb{Z}^2$, where

\begin{equation*}
\begin{aligned}
z_1 &= (-1,1), \, z_2 = (0,1), \, z_3 = (1,1), \, z_4=(-1,0), \\[1mm]
z_5 &= (1,0), \, z_6 = (-1,-1), \, z_7 = (0,-1), \, z_8=(1,-1).
\end{aligned}
\end{equation*}

The transition function $f:K^9\rightarrow K$ is given by

\begin{equation*}
f(x_0,\ldots \hspace{-0.02cm},x_8) =\begin{cases}
\mathbbm{1}_{x_1+\cdots+x_8=3}   , &\text{if } x_0=0,\\
\mathbbm{1}_{x_1+\cdots+x_8=2}+\mathbbm{1}_{x_1+\cdots+x_8=3},    & \text{if } x_0=1.
\end{cases}
\end{equation*}

The resulting CA is the Game of Life with dead cells corresponding to state $0$ and live cells to state $1$.

One of the main goals of this paper is to show that deep ReLU networks are capable of learning CA transition functions from CA evolution sequences. This result also establishes that every CA transition function can be realized by a deep ReLU network. Concretely, we build this universal CA transition function representation theorem in two steps. First, we establish a correspondence between CA transition functions and MV algebras <cit.>. Then, we show that
statements in MV logic are naturally expressed by deep ReLU networks. En route it is shown that MV logic is a suitable language for describing CA with general state sets as logical machines.

§.§ One-dimensional binary CA

In order to develop intuition on the connection between CA, MV logic, and ReLU networks, we consider a simple example setup. Specifically, we investigate CA with cellular space $\mathbb{Z}$, state set $K = \{0,1 \}$, and neighborhood set $\mathcal{E} = \{-1,0,1\}$. This class of CA is referred to as elementary CA <cit.>. As the transition function depends on the states of the cell itself and its left and right neighbor, there are $2^3$ possible states for a given cell's neighborhood and hence a total of $2^{2^3} = 256$ transition functions. Wolfram proposed a scheme <cit.> to index the associated CA using an $8$-digit binary string. For example, tbl:truth30 specifies the transition function of the elementary CA of index 30 ($30_{10} = 0001 1110_{2}$), and tbl:truth110 contains that corresponding to elementary CA 110 ($110_{10} = 0110 1110_{2}$). It is remarkable that elementary CA, despite their simplicity, are capable of universal computation. Specifically,

Cook <cit.> proved that elementary CA 110 is Turing complete.

$x_{-1}\,x_0\,x_1$    111     110     101     100     011     010     001     000
$f_{30}(x_{-1}, x_0, x_1)$     0     0     0     1     1     1     1     0

tableTransition function of elementary CA 30.

$x_{-1}\,x_0\,x_1$    111     110     101     100     011     010     001     000
$f_{110}(x_{-1}, x_0, x_1)$     0     1     1     0     1     1     1     0

tableTransition function of elementary CA 110.

We now bring Boolean logic into the picture, concretely by interpreting the state set elements 0 and 1 as truth values in Boolean algebra, and assuming that $x_{-1}, x_0, x_1$ are Boolean propositional variables.
The transition function $f_{110}$, for example, can in effect be regarded as a Boolean function, namely $f_{110} = \texttt{OR}(\texttt{XOR}(x_0,x_1), \texttt{AND}(\texttt{NOT}(x_{-1}),\texttt{OR}(x_0,x_1)))$.

In the same manner, we get $f_{30}(x_{-1}, x_0, x_1) = \texttt{XOR}(x_{-1}, \texttt{OR}(x_0,x_1))$. As shown in <cit.>, <cit.>, the transition function of every elementary CA can be expressed as a formula in Boolean logic, particularly in a canonical form called disjunctive normal form (DNF), which we define next.

In Boolean logic, a finite string of symbols is called a conjunction if it is given by  connections of a finite number of propositional variables (possibly negated), and a disjunction if the connections are effected through the  operation. A logical formula is said to be in disjunctive normal form if it is given by the  connection of a finite number of conjunctions.

Let the $\texttt{OR}$ operation be denoted by $\oplus$, $\texttt{AND}$ by $\odot$, and $\texttt{NOT}$ by $\lnot$. Following the procedure in <cit.>, we can express the transition function in tbl:truth30 in DNF according to

\begin{equation}\label{eq:rule30logic}
f_{30} = (x_{-1}\odot  \lnot x_0 \odot \lnot x_1) \oplus (\lnot x_{-1}\odot   x_1 ) \oplus (\lnot x_{-1}\odot   x_0).
\end{equation}

Having established the connection between elementary CA and Boolean logic, by expressing the transition function in terms of a Boolean formula, we are ready to have ReLU networks enter the story. Specifically, this shall be effected by realizing the Boolean logical operations $\oplus, \odot,$ and $\lnot$ through ReLU networks (see defn:ReLU below). With the ReLU function $\rho: \mathbb{R} \rightarrow \mathbb{R}, \hspace{0.2cm} \rho (x):=\max\{0,x\}$, the binary and $n$-ary conjunction can be written as

\begin{align}
x_1 \odot x_2 &= \rho\hspace{0.02cm}(x_1 + x_2 - 1), \\\label{eq:ReLUandgate}
x_1 \odot \cdots \odot x_n &=\rho\hspace{0.02cm}(\Sigma_{i = 1}^n x_i - (n-1)).
\end{align}

Likewise, the binary and $n$-ary disjunction can be realized according to

\begin{align}
x_1 \oplus x_2 &=\rho\hspace{0.02cm}(x_1+x_2) -\rho\hspace{0.02cm}(x_1 + x_2 - 1), \\\label{eq:ReLUorgate}
x_1 \oplus \cdots \oplus x_n &=\rho\hspace{0.02cm}(\Sigma_{i = 1}^n x_i) -\rho(\Sigma_{i = 1}^n x_i-1).
\end{align}

Finally, the Boolean operation $\lnot$ can be expressed as

\[
\lnot x = 1-x.
\]

We now combine these results to realize the logical formula $f_{30}$ in (<ref>) in terms of the ReLU network  $\Phi^{f_{30}}:\mathbb{R}^3 \rightarrow \mathbb{R} $ with three layers according to

\begin{equation}\label{eq:NNf30}
\Phi^{f_{30}}:= W_3^{f_{30}} \circ \rho \circ W_2^{f_{30}} \circ \rho \circ  W_1^{f_{30}}
\end{equation}

\begin{equation*}
W_1^{f_{30}}(x_{-1},x_0,x_1) = \begin{pmatrix}
-1 & -1 & 1 \\
1 & 0 & -1 \\
0 & 1 & -1
\end{pmatrix} \begin{pmatrix}
x_1 \\x_0 \\x_{-1}
\end{pmatrix} \text{, }
\end{equation*}

\begin{equation*}
W_2^{f_{30}}(x) = \begin{pmatrix}
1 & 1 & 1  \\
1 & 1 & 1
\end{pmatrix} x + \begin{pmatrix}
0 \\-1
\end{pmatrix}\text{, } \quad \text{ for } x\in \mathbb{R}^3,
\end{equation*}

\begin{equation*}
W_3^{f_{30}}(x) = \begin{pmatrix}1 & -1 \end{pmatrix} x,  \quad \text{ for } x\in \mathbb{R}^2.
\end{equation*}

It can now be verified directly that

\begin{equation*}
\Phi^{f_{30}}(x_{-1},x_0,x_1) = f_{30}(x_{-1},x_0,x_1), \text{ for } (x_{-1},x_0,x_1 ) \in \{0,1\}^3.
\end{equation*}

We hasten to add that ReLU network realizations of Boolean formulae are not unique. For example, upon noting that $\Phi^{\text{Id}}(x) = \rho(x) - \rho(-x)=x, x\in \mathbb{R},$ is a valid ReLU network, the map $\Phi^{f_{30}}$ in (<ref>) can be augmented according to

\begin{equation*}
\Phi^{\text{Id}} \circ \cdots \circ \Phi^{\text{Id}} \circ \Phi^{f_{30}}
\end{equation*}

without changing the input-output function it realizes. This results in an infinite collection of different ReLU networks, all realizing the Boolean formula $f_{30}$.

In addition, the Boolean logical formula corresponding to a given binary truth table is not unique either <cit.>. These issues of nonuniqueness will be addressed in more detail in sec:Extract.

We can summarize what was done in this section as follows. In an effort to implement the transition function of a one-dimensional binary CA by a ReLU neural network, we first expressed the transition function, originally specified in the form of a truth table, by a formula in Boolean logic. This formula was then realized by a ReLU network. The bridge via Boolean algebra effected in this manner serves to illustrate that CA are actually “logical machines”, an aspect elaborated on in the literature before  <cit.>, <cit.>, <cit.>, but notably only in the binary case.

Inspired by this connection, we shall show that the transition functions of general CA with arbitrary state sets, in particular of arbitrary cardinality, can be expressed in terms of formulae in MV logic <cit.>, which, in turn, are found to be naturally realized by ReLU networks. Embedding these (feedforward) networks into a recurrent neural network structure <cit.>, we get a dynamical system realization of general CA, entirely in terms of neural networks.

§.§ Previous work

Despite the fact that neural networks and CA stand for parallel efforts in building computational models of the human brain <cit.>, <cit.>, a profound theoretical understanding of the connections between the two structures seems to be lacking. Wulff and Hertz <cit.> consider shallow neural networks with the threshold activation function to learn the transition functions of one- and two-dimensional binary CA. While no details on the training algorithm are provided, the experiments in <cit.> lead to the conclusion that the approach employed is capable of learning only a small subset of the CA considered. More recently, Gilpin <cit.> designed a convolutional neural network simulating two-dimensional CA with arbitrary state sets and neighborhood size.

The connection between general CA and MV logic we report here seems completely new. For binary CA and Boolean logic, besides the classical results in <cit.>, <cit.>, more recent work studying CA evolution through the lens of formal logic is reported in <cit.>, <cit.>. Specifically, Sukanta and Chakraborty <cit.> show that binary CA on finite one-dimensional cellular spaces are models of certain logical languages in the domain of binary strings. Inokuch et al. <cit.> use the interpretation of binary CA transition functions as Boolean logical formulae to demonstrate that the multiplication of logical formulae, defined by monoid action, corresponds to the composition of global CA maps associated with the logical formulae under consideration.

As for the connection between MV logic and neural networks, Amato et al. <cit.> establish a correspondence between formulae in
Łukasiewicz logic and neural networks with rational weights and the clipped ReLU (CReLU) nonlinearity $\sigma(x) = \min\{1, \max\{0,x\}\}$. In <cit.> this result is further extended to a correspondence between Riesz MV algebras and CReLU networks with real weights; this is accomplished by adding scalar multiplication over the real numbers to propositional logic.

§.§ Outline of the paper

sec:MVLogic establishes that CA are in essence machines performing MV logic.

In sec:NetworkConstruction, we show how formulae in MV logic can be realized by (deep) ReLU networks.

The problem of extracting CA transition functions, along with underlying formulae in MV logic, from ReLU networks that have been trained on corresponding CA evolution data is addressed
in sec:Extract.

§ MANY-VALUED LOGIC AND CELLULAR AUTOMATA

We now turn to generalizing the correspondence between one-dimensional binary CA and Boolean logic. Specifically, we show that CA with arbitrary cellular space dimension and arbitrary state set (cardinality) can be seen as machines carrying out operations in propositional logic <cit.>, a many-valued extension of Boolean logic. The corresponding algebraic counterpart is known as Chang's many-valued algebras <cit.>. The connection we uncover is built on a fundamental result in MV logic,  which states that every piecewise linear function, whose linear pieces have integer coefficients, corresponds to a formula in MV logic and vice versa. Known in the literature as McNaughton theorem <cit.>, this result constitutes the equivalent of translating between binary truth tables and Boolean logical formulae.

§.§ Many-valued algebras

We start with a brief review, following <cit.>, of the basic elements in the theory of MV algebras.

A many-valued algebra is a structure $\mathcal{A} = \langle A, \oplus, \lnot,0\rangle$ consisting of a nonempty set $A$, a constant $0 \in A$, a binary operation  $\oplus$, and a unary operation $\lnot$ satisfying the following axioms:

\begin{align}
\label{eq:associativity}    x \oplus (y \oplus z) &= (x \oplus y) \oplus z\\
\label{commutivity}   x\oplus y  &= y\oplus x\\
\label{eq:NeutralZero}    x \oplus 0  &=x\\
\label{eq:NOTNOT}  \lnot \lnot x  &= x\\
\label{eq:NeutralOne}    x \oplus \lnot 0  &= \lnot 0\\
\label{eq:wedge operations}   \lnot (\lnot x \oplus y)\oplus y  &= \lnot (\lnot y \oplus x)\oplus x.
\end{align}

Specifically, (<ref>)-(<ref>) state that the structure $\langle A, \oplus, 0 \rangle$ is an abelian monoid. An MV algebra $\langle A, \oplus, \lnot,0\rangle$ is said to be nontrivial iff $|A| > 1$. On each MV algebra we can define a constant $1$ and a binary operation $\odot$ as follows:
\begin{align}
\label{eq:1}    1 &:=\lnot 0  \\
\label{eq:odot}    x \odot y&:= \lnot (\lnot x \oplus \lnot y).
\end{align}
The ensuing identities are then direct consequences of defn:MV algebra:

\begin{align}
\label{eq:associativity for times} x \odot (y \odot z) &= (x \odot y) \odot z\\
x\odot y  &= y\odot x\\
x \odot 1  &=x\\
\label{eq:0 for times}   x \odot  0  &=  0.
\end{align}

We will frequently use the notions of MV term and term function formalized as follows.

Let $n\in \mathbb{N}$ and $S_n = \{(, ), 0,\lnot, \oplus, x_1,\ldots \hspace{-0.02cm},$ $x_n\}$. An MV term in the variables $x_1, \ldots \hspace{-0.02cm}, x_n$ is a string over $S_n$ arising from a finite number of applications of the operations $\lnot$ and $\oplus$ as follows. The elements $0$ and $x_i$, for $i=1,\ldots \hspace{-0.02cm},n$, considered as one-element strings, are MV terms.

* If the string $\tau$ is an MV term, then $\lnot \tau $ is also an MV term.

* If the strings $\tau$ and $\gamma $ are MV terms, then $(\tau \oplus \gamma) $ is also an MV term.

In the remainder of the paper, we write $\tau (x_1, \ldots \hspace{-0.02cm}, x_n)$ to emphasize that $\tau$ is an MV term in the variables $x_1, \ldots \hspace{-0.02cm},x_n$.

For instance, the following strings over $S_2 = \{(, ), 0,\lnot, \oplus, x_1,x_2\}$ are MV terms in the variables $x_1$ and $x_2$:

\[
0, x_1, x_2,\lnot 0, \lnot x_2, x_1\oplus \lnot x_2.
\]

MV terms are syntactic expressions without semantic meaning. To endow them with semantics, an underlying MV algebra must be specified.

The resulting functions will be referred to as term functions.

Let $\tau(x_1, \ldots \hspace{-0.02cm}, x_n)$  be an MV term and

$\mathcal{A} = \langle A, \oplus, \lnot,0\rangle$ an MV algebra. The term function $\tau^\mathcal{A}:A^n\rightarrow A $ associated with $\tau$
under $\mathcal{A}$ is obtained by substituting, for $i=1,\ldots \hspace{-0.02cm}, n$, $a_i \in A$ for all occurrences of $x_i$ in $\tau$ and
interpreting the symbols $\oplus$ and $\lnot$ according to how they are specified in $\mathcal{A}$.

defined as follows. For every input $(a_1,\ldots \hspace{-0.02cm}, a_n) \in A^n$, first substitute $a_i$ for all occurrences of the variable $x_i$ in $\tau$:

* $x_i^\Acal: (a_1,\ldots \hspace{-0.02cm}, a_n)\mapsto a_i$,  for $i=1,\ldots \hspace{-0.02cm}, n$.

Then, proceed by induction on the number of operation symbols, i.e., $\oplus$ and $\lnot$, occuring in $\tau$ by applying the following rules:

* $(\lnot \tau)^\Acal = \lnot \tau^\Acal$,  for MV term $\tau$,

* $ (\tau \oplus \gamma )^\Acal = \tau ^\Acal \oplus \gamma ^\Acal$,   for MV terms $\tau$ and $\gamma$.

We next demonstrate that Boolean algebra is an MV algebra.

Consider the set $B = \{0,1\}.$ Define the binary operation $\oplus$ and the unary operation $\lnot$ on $B$ according to

\begin{equation*}
\begin{aligned}[c]
0\oplus 0 &=0\\
0\oplus 1 &=1\\
1\oplus 0 &=1 \\
1\oplus 1 &=1
\end{aligned}
\qquad \text{and} \qquad
\begin{aligned}[c]
\lnot 0 &=1\\
\lnot 1 &=0.\\
\end{aligned}
\end{equation*}
It is readily verified that the resulting structure $ \mathcal{B} = \langle B, \oplus, \lnot  ,0\rangle$ satisfies the axioms in defn:MV algebra and is hence an MV algebra. In fact, by further defining the binary operation $\odot$ on $B$ through

\begin{equation*}
x  \odot y := \lnot(\lnot x \oplus \lnot y)\text{,}
\end{equation*}
which yields
\begin{equation*}
\begin{aligned}
0 \odot 0 &= 0 \\
0 \odot 1 &= 0 \\
1 \odot 0 &= 0 \\
1 \odot 1 &= 1,
\end{aligned}
\end{equation*}
we immediately see that the operations $\odot, \oplus$, and $ \lnot$ correspond to the , , and  operations, respectively, in Boolean logic with $0$ designating  and $1$ standing for .

Recall the informal discussion in Section <ref> explaining how the transition functions of binary one-dimensional CA can be expressed in terms of Boolean logic. We now formalize this observation (for arbitrary cellular space dimension) by casting it into MV algebras.

Consider a CA with cellular space dimension $d \in \mathbb{N}$, neighborhood size $n \in \mathbb{N}$, state set $K = \{0,1\}$, and transition function $f:K^n \rightarrow K$ specified in the form of a (Boolean) truth table. There exists an MV term $\tau(x_1,\ldots \hspace{-0.02cm}, x_n)$ with associated term under the Boolean algebra $\mathcal{B} =  \langle  \{0,1\}, \oplus, \lnot  ,0\rangle$ satisfying
\begin{equation*}
\tau^{\mathcal{B}}(x_1,\ldots \hspace{-0.02cm}, x_n) = f(x_1,\ldots \hspace{-0.02cm}, x_n), \quad \text{ for all } (x_1,\ldots \hspace{-0.02cm}, x_n) \in \{0,1\}^n.
\end{equation*}

We first transform the truth table specifying the transition function $f$

into a Boolean formula <cit.>.

Identifying, in the resulting Boolean algebraic expression,

the Boolean operations
$\texttt{AND}$, $\texttt{OR}$, and $\texttt{NOT}$ with the operations
$\odot$, $\oplus$, and $\lnot$ in a general MV algebra yields a valid MV term
$\tau$. Converting $\tau$ into its associated term function $\tau^{\mathcal{B}}$ according to defn:TermFunction establishes the result.

In order to extend lemma2.6 to CA with state sets of arbitrary cardinality, we need a logic that can cope with more than two truth values. We shall see that the theory of MV logic <cit.> provides a suitable framework and start by developing the underlying algebraic structure, namely that of MV algebras.

§.§ MV algebras and CA

We begin with a simple example of an MV algebra.

For $k\in\mathbb{N}$, $k \geq 2$, consider the CA state set
\begin{equation}\label{eq:KStateSet}
K = \Bigl\{  0, \frac{1}{k-1}, \ldots \hspace{-0.02cm}, \frac{k-2}{k-1},1 \Bigr\}
\end{equation}
of cardinality $k$. Define the binary operation $\oplus$ and the unary operation $\lnot$ on $K$ according to

\begin{align*}
x \oplus y &= \min \{1, x+y\} \\
\lnot x &= 1-x.
\end{align*}
It is readily seen that the structure $\mathcal{A}_k:=\langle K, \oplus, \lnot,0 \rangle$ satisfies the axioms in defn:MV algebra and hence constitutes an MV algebra. For $k=2$, $\mathcal{A}_k$ reduces to the Boolean algebra in defn:boolean.

For the MV algebra $\mathcal{A}_k$ in exmp:MVAlgebraAk, given an MV term $\tau$ in $n$ variables, the associated term function $\tau^{\mathcal{A}_k}: K^n \rightarrow K $ can be interpreted as a CA transition function.
We illustrate this observation through a simple example, but first note that $\tau^{\mathcal{A}_k}$ is guaranteed to map
$K^n$ to $K$ owing to the specific structure of the state set $K$ and the operations $\oplus$ and $\lnot$.

Consider the MV algebra $\mathcal{A}_3$ according to Example <ref> with the associated set $K = \{0, 1/2,1\}$.

The term function $\tau^{ \mathcal{A}_3}: K^3 \rightarrow K$ corresponding to the
MV term $\tau = x_{-1}\oplus x_0 \oplus x_1$ equals the transition function of the $3$-color totalistic CA <cit.> with cellular space $\mathbb{Z}$, neighborhood set $\mathcal{E}=\{-1, 0, 1\}$, state set $K$, and transition function $f$ as specified in tbl:totalistic (note that $f$ depends only on the sum of the neighborhood state values).

$x_{-1}+x_0+x_1$    3     5/2     2     3/2     1     1/2     0
$f(x_{-1}, x_0, x_1)$     1     1    1     1     1     1/2     0

tableTransition function $f$ of the totalistic CA.

Inspired by this example, we now ask whether every CA transition function has an underlying MV term with associated term function under a suitable MV algebra equal to the transition function. Formally speaking, we seek a generalization of Lemma <ref> to the case of CA with arbitrary state set cardinality $k$. If answered in the affirmative, this would show that CA are essentially machines that perform operations in a suitable logic. It turns out that this is, indeed, the case, but we need to consider an MV algebra that accommodates state sets $K$ of arbitrary cardinality $k$.

For a CA with the state set $K$ in (<ref>), neighborhood size $n\in \mathbb{N}$, and arbitrary cellular space dimension $d\in \mathbb{N}$, the associated transition rule is a function $f: K^n \rightarrow K$. Conversely, as illustrated by exmp:2.1 and exmp:2.2, for the MV algebra $\mathcal{A}_k$ in exmp:MVAlgebraAk, given an MV term $\tau$ in $n$ variables, the associated term function $\tau^{\mathcal{A}_k}: K^n \rightarrow K $ can be interpreted as a CA transition rule. Based on this simple observation, we will now develop a systematic correspondence between CA and MV term functions. This correspondence endows CA transition rules with logical meaning. In order to elucidate the fundamental structure of the so-identified connection, we will enlarge $\mathcal{A}_k$ to encompass all possible CA state sets. This will require more general underlying MV algebras, which, as it turns out, are given by the unit interval on $\mathbb{R}$ endowed with truncated addition and negation.

Consider the unit interval $[0, 1]$ on $\mathbb{R}$, and define $x \oplus y = \min \{1, x+y\}$ and $\lnot x = 1-x$, for  $x, y \in [0,1]$. It can be verified that the structure $\mathcal{I} = \langle [0,1], \oplus, \lnot,0 \rangle $ is an MV algebra. In particular, $\mathcal{I}$ constitutes the algebraic counterpart of propositional logic <cit.>. We further define the operation $x  \odot y := \lnot(\lnot x \oplus \lnot y) = \max\{0,x+y-1\}$.

The MV algebra $\mathcal{I}$ in defn:MV01 is the so-called standard MV algebra.

As shown in <cit.>, <cit.>, an equation holds in every MV algebra iff it holds in the standard MV algebra $\mathcal{I}$, endowing $\mathcal{I}$ with universality.
Furthermore, the unit interval $[0,1]$ accommodates all possible CA state set cardinality. Let $\tau(x_1,\ldots \hspace{-0.02cm},x_n)$ be an MV term and $\tau^{\mathcal{I}}:[0,1]^n \rightarrow [0,1]$ the associated term function under $\mathcal{I}$. For every state set cardinality $k \in \mathbb{N}, k\geq 2$ in (<ref>),  the function $f: K^n \rightarrow K$ given by

\[
f(x_1,\ldots \hspace{-0.02cm},x_n) := \tau^{\mathcal{I}}(x_1,\ldots \hspace{-0.02cm},x_n), \quad \text{ for all } (x_1,\ldots \hspace{-0.02cm},x_n) \in K^n,
\]

is the transition function of some CA with state set $K$ and neighborhood size $n$. In other words, CA can be interpreted as machines carrying out MV logical operations. Conversely, given a CA transition function, we wish to construct an associated logical formula to endow the underlying CA with logical operational meaning.
Next, recall that for the binary state set $\{0,1\}$, in the proof of Lemma <ref> we started by converting the binary truth table specifying
the transition function $f$ into a formula in Boolean logic. In MV logic truth tables are given by mappings from $[0,1]^n$ to $[0,1]$.
While in the Boolean case,

every truth table can be cast into a formula in Boolean logic <cit.>, this does not
hold for the MV algebra $\mathcal{I}$.

The McNaughton theorem, stated next, explicitly characterizes the class of functions from $[0,1]^{n}$ to $[0,1]$ that have underlying MV terms $\tau$. Here, “underlying” means that the term function associated with $\tau$ under $\mathcal{I}$ equals the function under consideration.

[McNaughton theorem <cit.>]
Consider the MV algebra $\mathcal{I} = \langle [0,1], \oplus, \lnot,0 \rangle $ in defn:MV01. Let $n\in \mathbb{N}$. For a function $f_{\text{c}}: [0,1]^n \rightarrow [0,1]$ to have an associated MV term $\tau$ such that $\tau ^{\mathcal{I}} = f_{\text{c}}$ on $[0,1]^n$, it is necessary and sufficient that

* $f_{\text{c}}$ is continuous with respect to the natural topology on $[0,1]^n$,

* there exist linear functions $p_1, \ldots \hspace{-0.02cm}, p_\ell$ with integer coefficients, i.e.,

\begin{equation*}
p_j(x_1, \ldots \hspace{-0.02cm}, x_{n}) = m_{j1}x_1+\cdots+m_{jn}x_{n}+b_j, \quad  j = 1,\ldots \hspace{-0.02cm},\ell,
\end{equation*}
where $ m_{j1},\ldots \hspace{-0.02cm}, m_{jn},b_j \in \mathbb{Z}$, for $j=1,\ldots \hspace{-0.02cm}, \ell$, such that for every $x \in [0,1]^n$, there is a $j \in \{1,\ldots \hspace{-0.02cm}, \ell\}$ with $ f_{\text{c}}(x) = p_j(x)$.
Functions satisfying these conditions are called McNaughton functions.

We now demonstrate how the McNaughton theorem can be employed to decide whether a given CA transition function has an underlying MV term. Since term functions under $\mathcal{I}$ are continuous piecewise linear functions according to them:McNaughtonTheorem, we start by linearly interpolating CA transition functions $f:K^n\rightarrow K$ to continuous piecewise linear functions $f_{\text{c}}:[0,1]^n\rightarrow [0,1]$. The following simple example illustrates the approach we pursue.

Consider a one-dimensional CA with cellular space $\mathbb{Z}$, neighborhood set $\mathcal{E} = \{ 0\}$, state set $K = \{0, 1/3, 2/3,1\}$, and transition function $f$ according to

$x$    $0$     $1/3$     $2/3$     $1$
$f(x)$     $0$     $1$     $1$     $0$

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]
[->] (0,0) – (0,3.2) node[anchor=west] $f(x)$;
[->] (0,0) – (3.2,0) node[anchor= west] $x$;
(0.07,1) – (-0.07,1) node[anchor= east] $1/3$;
(0.07,2) – (-0.07,2) node[anchor= east] $2/3$;
(0.07,3) – (-0.07,3) node[anchor= east] $1$;

(2,-0.02) – (2,0.02) node[anchor= north] $2/3$;
(1,-0.02) – (1,0.02) node[anchor= north] $1/3$;
(3,-0.02) – (3,0.02) node[anchor= north] ;
at (3,-0.2) $1$;
at (0.2,-0.3 cm) $0$;
[black,fill=black] (0,0) circle (.5ex);
[black,fill=black] (1,3) circle (.5ex);
[black,fill=black] (2,3) circle (.5ex);
[black,fill=black] (3,0) circle (.5ex);

[->] (5,0) – (5,3.2) node[anchor=west] $f_{\text{c}}(x)$;
[->] (5,0) – (8.2,0) node[anchor= west] $x$;
(5.07,1) – (4.93,1) node[anchor= east] $1/3$;
(5.07,2) – (4.93,2) node[anchor= east] $2/3$;
(5.07,3) – (4.93,3) node[anchor= east] $1$;
(7,-0.02) – (7,0.02) node[anchor= north] $2/3$;
(6,-0.02) – (6,0.02) node[anchor= north] $1/3$;
(8,-0.02) – (8,0.02) node[anchor= north] ;
at (8,-0.2) $1$;
at (5.2,-0.3 cm) $0$;
[black,fill=black] (5,0) circle (.5ex);
[black,fill=black] (6,3) circle (.5ex);
[black,fill=black] (7,3) circle (.5ex);
[black,fill=black] (8,0) circle (.5ex);
(5,0) – (6,3) node[anchor= north] ;
(7,3) – (6,3) node[anchor= north] ;
(7,3) – (8,0) node[anchor= north] ;

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]
[->] (0,0) – (0,3.2) node[anchor=west] $f(x)$;
[->] (0,0) – (3.2,0) node[anchor= west] $x$;
(0.07,1) – (-0.07,1) node[anchor= east] $1/3$;
(0.07,2) – (-0.07,2) node[anchor= east] $2/3$;
(0.07,3) – (-0.07,3) node[anchor= east] $1$;

(2,-0.02) – (2,0.02) node[anchor= north] $2/3$;
(1,-0.02) – (1,0.02) node[anchor= north] $1/3$;
(3,-0.02) – (3,0.02) node[anchor= north] ;
at (3,-0.2) $1$;
at (0.2,-0.3 cm) $0$;
[black,fill=black] (0,0) circle (.5ex);
[black,fill=black] (1,3) circle (.5ex);
[black,fill=black] (2,3) circle (.5ex);
[black,fill=black] (3,0) circle (.5ex);

[->] (5,0) – (5,3.2) node[anchor=west] $f_{\text{c}}(x)$;
[->] (5,0) – (8.2,0) node[anchor= west] $x$;
(5.07,1) – (4.93,1) node[anchor= east] $1/3$;
(5.07,2) – (4.93,2) node[anchor= east] $2/3$;
(5.07,3) – (4.93,3) node[anchor= east] $1$;
(7,-0.02) – (7,0.02) node[anchor= north] $2/3$;
(6,-0.02) – (6,0.02) node[anchor= north] $1/3$;
(8,-0.02) – (8,0.02) node[anchor= north] ;
at (8,-0.2) $1$;
at (5.2,-0.3 cm) $0$;
[black,fill=black] (5,0) circle (.5ex);
[black,fill=black] (6,3) circle (.5ex);
[black,fill=black] (7,3) circle (.5ex);
[black,fill=black] (8,0) circle (.5ex);
(5,0) – (6,3) node[anchor= north] ;
(7,3) – (6,3) node[anchor= north] ;
(7,3) – (8,0) node[anchor= north] ;

The transition function $f$ and the associated function $ f_{\text{c}}$ in example McNaughton is not enough.

The continuous piecewise linear function $ f_{\text{c}}:[0,1] \rightarrow [0,1]$ obtained by interpolating $f$ is given by

\begin{equation}\label{eq:egf}
\begin{cases}
3x, &0\leq x< 1/3 \\
1,&1/3\leq x< 2/3 \\
-3x+3,&2/3\leq x\leq 1 \\
\end{cases}.
\end{equation}

As $f_{\text{c}}$ in (<ref>) has exclusively integer coefficients and is piecewise linear and continuous by construction, it
satisfies Conditions 1 and 2 in them:McNaughtonTheorem. Hence, there is an underlying MV term,
shown in Appendix <ref> to be given by
\begin{equation*}
\tau=  (x\oplus x\oplus x) \wedge  \lnot 0 \wedge (x\odot x\odot x),
\end{equation*}

where, for brevity, we write $x \wedge y = (x\oplus \lnot y)\odot y$.

However, for the same cellular space and state set, a slightly different transition function $g$, given by

$x$    $0$     $1/3$     $2/3$     $1$
$g(x)$     $0$     $1/3$     $1$     $0$

yields the associated piecewise linear function

\begin{equation*}
\begin{cases}
x, &0\leq x< 1/3 \\
2x-1/3,&1/3\leq x< 2/3 \\
-3x+3,&2/3\leq x\leq 1 \\
\end{cases},
\end{equation*}

which fails to satisfy Condition 2 in them:McNaughtonTheorem, as one of its linear pieces exhibits a noninteger coefficient. As the McNaughton theorem is iff, we can conclude that $g_{\text{c}}$ does not have an underlying MV term $\tau$ such that $\tau^{\mathcal{I}} =  g_{\text{c}}$.

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]

[->] (0,-5) – (0,-1.8) node[anchor=west] $g(x)$;
[->] (0,-5) – (3.2,-5) node[anchor= west] $x$;
(0.07,-4) – (-0.07,-4) node[anchor= east] $1/3$;
(0.07,-3) – (-0.07,-3) node[anchor= east] $2/3$;
(0.07,-2) – (-0.07,-2) node[anchor= east] $1$;

(2,-5.02) – (2,-4.98) node[anchor= north] $2/3$;
(1,-5.02) – (1,-4.98) node[anchor= north] $1/3$;
(3,-5.02) – (3,-4.98) node[anchor= north] ;
at (3,-5.2) $1$;
at (0.2,-5.2) $0$;
[black,fill=black] (0,-5) circle (.5ex);
[black,fill=black] (1,-4) circle (.5ex);
[black,fill=black] (2,-2) circle (.5ex);
[black,fill=black] (3,-5) circle (.5ex);

[->] (5,-5) – (5,-1.8) node[anchor=west] $g_{\text{c}}(x)$;
[->] (5,-5) – (8.2,-5) node[anchor= west] $x$;
(5.07,-4) – (4.93,-4) node[anchor= east] $1/3$;
(5.07,-3) – (4.93,-3) node[anchor= east] $2/3$;
(5.07,-2) – (4.93,-2) node[anchor= east] $1$;
(7,-5.02) – (7,-4.98) node[anchor= north] $2/3$;
(6,-5.02) – (6,-4.98) node[anchor= north] $1/3$;
(8,-5.02) – (8,-4.98) node[anchor= north] ;
at (8,-5.2) $1$;
at (5.2,-5.2) $0$;
[black,fill=black] (5,-5) circle (.5ex);
[black,fill=black] (6,-4) circle (.5ex);
[black,fill=black] (7,-2) circle (.5ex);
[black,fill=black] (8,-5) circle (.5ex);
(5,-5) – (6,-4) node[anchor= north] ;
(7,-2) – (6,-4) node[anchor= north] ;
(7,-2) – (8,-5) node[anchor= north] ;

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]

[->] (0,-5) – (0,-1.8) node[anchor=west] $g(x)$;
[->] (0,-5) – (3.2,-5) node[anchor= west] $x$;
(0.07,-4) – (-0.07,-4) node[anchor= east] $1/3$;
(0.07,-3) – (-0.07,-3) node[anchor= east] $2/3$;
(0.07,-2) – (-0.07,-2) node[anchor= east] $1$;

(2,-5.02) – (2,-4.98) node[anchor= north] $2/3$;
(1,-5.02) – (1,-4.98) node[anchor= north] $1/3$;
(3,-5.02) – (3,-4.98) node[anchor= north] ;
at (3,-5.2) $1$;
at (0.2,-5.2) $0$;
[black,fill=black] (0,-5) circle (.5ex);
[black,fill=black] (1,-4) circle (.5ex);
[black,fill=black] (2,-2) circle (.5ex);
[black,fill=black] (3,-5) circle (.5ex);

[->] (5,-5) – (5,-1.8) node[anchor=west] $g_{\text{c}}(x)$;
[->] (5,-5) – (8.2,-5) node[anchor= west] $x$;
(5.07,-4) – (4.93,-4) node[anchor= east] $1/3$;
(5.07,-3) – (4.93,-3) node[anchor= east] $2/3$;
(5.07,-2) – (4.93,-2) node[anchor= east] $1$;
(7,-5.02) – (7,-4.98) node[anchor= north] $2/3$;
(6,-5.02) – (6,-4.98) node[anchor= north] $1/3$;
(8,-5.02) – (8,-4.98) node[anchor= north] ;
at (8,-5.2) $1$;
at (5.2,-5.2) $0$;
[black,fill=black] (5,-5) circle (.5ex);
[black,fill=black] (6,-4) circle (.5ex);
[black,fill=black] (7,-2) circle (.5ex);
[black,fill=black] (8,-5) circle (.5ex);
(5,-5) – (6,-4) node[anchor= north] ;
(7,-2) – (6,-4) node[anchor= north] ;
(7,-2) – (8,-5) node[anchor= north] ;

The transition function $g$ and the continuous function $ g_{\text{c}}$ in example McNaughton is not enough.

Example <ref> manifests missing structure for the MV algebra $\mathcal{I}$ to allow universal conversion from CA transition functions to MV terms. Concretely, inspection of the state set $K$ shows that we would need a McNaughton theorem that allows rational coefficients in Condition 2 of them:McNaughtonTheorem. It turns out that so-called divisible many-valued (DMV) algebras resolve this issue.

Consider the MV algebra $\mathcal{I} = \langle [0,1], \oplus, \lnot,0 \rangle $ in defn:MV01. Define the family of unary operations $\{\delta_i:[0,1]\rightarrow [0,1] \}_{i \in \mathbb{N}}$ according to
\begin{equation*}
\delta_i x = \frac{1}{i}x,\quad x\in[0,1], \quad \text{ for all } i\in \mathbb{N}.
\end{equation*}

It is easily verified that the structure $\mathcal{I}_{\text{d}}=\langle [0,1], \oplus, \lnot,0,  \{\delta_i \}_{i \in \mathbb{N}}\rangle$ is a DMV algebra <cit.>. In particular, it constitutes the algebraic counterpart of Rational logic <cit.>.

The following result, often referred to as rational McNaughton theorem, explicitly characterizes the class of term functions under the DMV algebra $\mathcal{I}_{\dt}$.

[Rational McNaughton theorem <cit.>]
Consider the DMV algebra $\Ical_{\dt}=\langle [0,1], \oplus, \lnot,0,  \{\delta_i \}_{i \in \mathbb{N}}\rangle$ in defn:DivisibleMV01. Let $n \in \mathbb{N}$. For a function $f_{\text{c}}: [0,1]^n \rightarrow [0,1]$ to have an associated DMV term $\tau$ such that $\tau^{\Ical_{\dt}} = f_{\text{c}}$ on $[0,1]^n$, it is necessary and sufficient that

* $f_{\text{c}}$ is continuous with respect to the natural topology on $[0,1]^n$,

* there exist linear functions $p_1, \ldots \hspace{-0.02cm}, p_\ell$ with rational coefficients, i.e.,

\begin{equation*}
p_j(x_1, \ldots \hspace{-0.02cm}, x_{n}) = m_{j1}x_1+\cdots+m_{jn}x_{n}+b_j, \quad j=1,\ldots \hspace{-0.02cm},\ell,
\end{equation*}
where $b_j, m_{j1},\ldots \hspace{-0.02cm}, m_{jn} \in \mathbb{Q}$, for $j=1,\ldots \hspace{-0.02cm}, \ell$, such that for every $x \in [0,1]^n$, there is a $j \in \{1,\ldots \hspace{-0.02cm}, \ell\}$ with $f_{\text{c}}(x) = p_j(x)$.

The central result we have been working towards in this section now follows readily.

Consider a CA with cellular space dimension $d$, neighborhood size $n $, state set $K = \{0, 1/(k-1),\ldots \hspace{-0.02cm}, (k-2)/(k-1), 1\}$ of cardinality $k\in \mathbb{N}, k\geq 2$, and transition function $f: K^{n} \rightarrow K$. There exists a DMV term $\tau(x_1,\ldots \hspace{-0.02cm}, x_n)$ with associated term function under the DMV algebra $\mathcal{I}_{\text{d}}$ satisfying
\begin{equation*}
\tau^{\Ical_{\dt}}(x_1,\ldots \hspace{-0.02cm}, x_n) = f(x_1,\ldots \hspace{-0.02cm}, x_n), \quad \text{ for all } (x_1,\ldots \hspace{-0.02cm}, x_n) \in K^n.
\end{equation*}

We start by linearly interpolating the CA transition funtion $f: K^n \rightarrow K$ to a continuous function $f_{\text{c}}: [0,1]^n\rightarrow [0,1]$. As all the points
\begin{equation*}
\left(x_1,\ldots \hspace{-0.02cm}, x_n, f(x_1,\ldots \hspace{-0.02cm}, x_n)\right)\in K^{n+1}
\end{equation*}
have exclusively rational coordinates, the linear pieces of the interpolated function $f_{\text{c}}$ have rational coefficients as well. Application of rationalMcNaughtonTheorem to $f_{\text{c}}$ then yields a DMV term $\tau$ satisfying $\tau^{\Ical_{\dt}} = f_{\text{c}}$ on $[0,1]^n$. This implies $\tau^{\Ical_{\dt}} = f_{\text{c}} = f$ on $K^n$, as desired.

Note that we extend the MV algebra $\mathcal{I}$ to the DMV algebra $\mathcal{I}_{\dt}$ with the goal of constructing logical formulae associated with CA transition functions. The choice of linear interpolation is attributed to the McNaughton theorem. Whilst all term functions under $\Ical$ can work as the transition function of some CA, as argued above, not all term functions under $\mathcal{I}_{\dt}$ can be viewed as CA transition rules. Take a simple example, namely $\tau = \delta_2 x$, its associated term function $\tau^{\Ical_{\dt}}: x \mapsto  x/2$, cannot induce a valid CA transition rule. In other words, the class of linearly interpolated CA transition functions is a proper subset of term functions under $\mathcal{I}_{\dt}$.

Consider a CA with cellular space $\mathbb{Z}$, the neighborhood set $\mathcal{E}=\{-1,0,1\}\subset \mathbb{Z}$, and state set $K$ in (<ref>) with $k\geq 3$. The transition rule $f$ is given by

\begin{equation*}\label{eq:kcolorRule}
f(x_{-1},x_0,x_1) =\begin{cases}
(x_0 + \frac{1}{k-1}) \mathbbm{1}_{ \{ x_{-1} = x_0 +\frac{1}{k-1} \text{ or } x_1 = x_0+\frac{1}{k-1}\}} &\quad \\
\hspace{0.7cm}  +\; x_0\mathbbm{1}_{\{ x_{-1} \neq x_0+\frac{1}{k-1}\} }\mathbbm{1}_{ \{ x_{1} \neq x_0+\frac{1}{k-1}\} }, & \text{ if } x_0 \in K \text{\textbackslash} \{1\}, \\
x_0\mathbbm{1}_{\{ x_{-1} \neq 0\} }\mathbbm{1}_{ \{ x_{1} \neq 0\} }, & \text{ if } x_0=1.
\end{cases}
\end{equation*}

which has a corresponding MV term under the MV algebra $\mathcal{A}_k$ in exmp:MVAlgebraAk. The resulting CA is the $k$-colored cyclic CA created by Bramson and Griffeath <cit.> to model particle systems.

§ RECURRENT NEURAL NETWORKS REALIZE CA DYNAMICS

In this section, we show that recurrent neural networks can realize the overall dynamical behavior of CA. This will be accomplished in two steps. First, we demonstrate that ReLU  networks naturally realize operations in DMV algebras, which, in turn, leads to a universal realization theorem for CA transition functions. The ReLU network realizing the CA transition function is then embedded into an RNN that emulates the dynamics of the CA.

§.§ DMV algebras and ReLU neural networks

We start by formally defining ReLU neural networks <cit.>.

Let $L \in \mathbb{N}$ and $N_0,$ $N_1, \ldots \hspace{-0.02cm},$ $N_L \in \mathbb{N}$. A ReLU neural network is a map $\Phi: \mathbb{R}^{N_0} \rightarrow \mathbb{R}^{N_L} $ given by

\[ \Phi= \begin{cases}
W_1, & L=1 \\
W_2 \circ \rho \circ W_1, & L=2 \\
W_L \circ \rho \circ W_{L-1} \circ \rho \circ \cdots \circ \rho  \circ W_1, & L \geq 3
\end{cases} \text{,}
\]
where, for $\ell \in \{ 1,2,\ldots \hspace{-0.02cm}, L\}$, $W_\ell: \mathbb{R}^{N_{\ell-1}}\rightarrow \mathbb{R}^{N_\ell}, W_\ell(x):= A_\ell x+b_\ell $ are affine transformations with weight matrices $A_\ell = \mathbb{R}^{N_\ell\times N_{\ell-1}}$ and bias vectors $b_\ell \in \mathbb{R}^{N_\ell}$, and the ReLU activation function $\rho: \mathbb{R}\rightarrow \mathbb{R},\rho\hspace{0.02cm}(x):=\max\{0,x\}$ acts component-wise. We denote by $\mathcal{N}_{d, d'}$ the set of ReLU neural networks of input dimension $N_0=d$ and output dimension $N_L = d'$. The number of layers of the network $\Phi$, denoted by $\mathcal{L}(\Phi)$, is defined to equal $L$.

It follows immediately from defn:ReLU that ReLU networks realize continuous piecewise linear functions. We shall frequently make use of basic ReLU network constructions, namely compositions <cit.>, augmentations <cit.>, and parallelizations <cit.>, collected here for completeness. The proofs of the corresponding results lem:ReLUConcatenating-<ref> provided in <cit.> present explicit network constructions that we will occasionally refer to.

Let $d_1, d_2, d_3 \in \mathbb{N}$, $\Phi_1 \in \mathcal{N}_{d_1,d_2}$, and $\Phi_2 \in \mathcal{N}_{d_2,d_3}.$ There exists a network $\Psi \in \mathcal{N}_{d_1,d_3}$ with $\mathcal{L}(\Psi) = \mathcal{L}(\Phi_1)+\mathcal{L}(\Phi_2)$, satisfying
\begin{equation*}
\Psi(x) = (\Phi_2\circ \Phi_1)(x), \quad \text{ for all } x \in \mathbb{R}^{d_1}.
\end{equation*}

Let $d_1, d_2, L \in \mathbb{N}$, and $\Phi \in \mathcal{N}_{d_1,d_2}$ with $\mathcal{L}(\Phi)<L$. There exists a network $\Psi \in \mathcal{N}_{d_1,d_2}$ with $\mathcal{L}(\Psi)=L$, satisfying $\Psi(x) = \Phi(x)$, for all $x\in \mathbb{R}^{d_1}$.

Let $n, d, L \in \mathbb{N}$ and, for $i \in \{1,\ldots \hspace{-0.02cm}, n\}$, let $d_i'\in \mathbb{N}$ and $\Phi_i \in \Ncal_{d,d_i'}$ with $\mathcal{L}(\Phi_i) = L$. There exists a network $\Psi \in \Ncal_{d, \sum_{i=1}^nd_i'}$ with $\mathcal{L}(\Psi) = L$, satisfying
\begin{equation*}
\Psi(x) = (\Phi_1(x),\ldots \hspace{-0.02cm}, \Phi_n(x)) \in \mathbb{R}^{\sum_{i=1}^nd_i'}, \quad \text{ for all } x\in \mathbb{R}^{d}.
\end{equation*}

We proceed with the  ReLU network constructions realizing operations in the DMV algebra $\mathcal{I}_{\text{d}}=\langle [0,1], \oplus, \lnot,0, \{\delta_i \}_{i \in \mathbb{N}}\rangle$. To this end, we start by noting that the operation $\lnot x = 1-x$ is trivially implemented by a ReLU network with one layer according to

\begin{equation*}
\lnot x = W_1(x) = 1-x, \quad \text{ for } x \in [0,1].
\end{equation*}

By the same token, as the operations $\delta_i x = \frac{1}{i}x$, for $i\in \mathbb{N}$, are affine mappings, they can be realized by single-layer ReLU networks. The following lemma details the ReLU network constructions realizing the operations $x\oplus y =\min\{1,x+y\}$ and $x\odot y = \max\{0,x+y-1\}$ in $\mathcal{I}_{\text{d}}$.

There exist ReLU networks $\Phi^{\oplus}\in \mathcal{N}_{2,1}$ and $\Phi^{\odot } \in \mathcal{N}_{2,1}$ satisfying
\begin{align*}
\label{eq:networkAND}    \Phi^{\oplus }(x,y) &= \min\{1,x+y\} \\
\Phi^{\odot }(x,y) &= \max\{0,x+y-1\},
\end{align*}
for all $x,y\in [0,1].$

First, to realize the operation $x\oplus y =\min\{1,x+y\} $, we note that addition can be implemented by a single-layer ReLU network according to

\begin{equation*}
x+y =  \begin{pmatrix} 1 & 1 \end{pmatrix}  \begin{pmatrix} x \\y \end{pmatrix}.
\end{equation*}

For the “$\min$” operation, we observe that

\begin{equation*}
\min\{1,x\} = 1-\rho(1-x) = (W_2 \circ \rho \circ W_1)(x), \quad x \in [0,1],
\end{equation*}

\begin{equation*}
W_1(x) = -x+1, \quad \quad  W_2(x) = -x+1.
\end{equation*}

Now, applying lem:ReLUConcatenating

to concatenate the networks  $\Phi_1(x,y) = \begin{pmatrix} 1 & 1 \end{pmatrix}  \begin{pmatrix} x \\y \end{pmatrix} $ and $\Phi_2(x) = (W_2 \circ \rho \circ W_1) (x)$ yields the desired ReLU network realization of $x\oplus y$ according to

\begin{equation*}
x \oplus y = (W_2^{\oplus} \circ \rho \circ W_1^{\oplus})(x,y), \quad x,y\in [0,1],
\end{equation*}

\begin{equation*}
W_1^{\oplus}(x,y)= \begin{pmatrix}
-1 & -1
\end{pmatrix} \begin{pmatrix}
x \\y
\end{pmatrix}+1, \quad
W_2^{\oplus}(x) = -x+1.
\end{equation*}

To realize the operation $x\odot y = \max\{0,x+y-1\}$, we directly note that
\begin{equation*}
\max\{0,x+y-1\} =  \rho \left( \begin{pmatrix}
1 & 1
\end{pmatrix}   \begin{pmatrix}
x \\ y
\end{pmatrix} -1 \right) = (W_2^{\odot} \circ \rho \circ W_1^{\odot})(x,y),
\end{equation*}
for $x,y\in[0,1]$, where
\begin{equation*}
W_1^{\odot} (x,y) = \begin{pmatrix}
1 & 1
\end{pmatrix}   \begin{pmatrix}
x \\ y
\end{pmatrix} -1,\quad \quad  W_2^{\odot}(x) = x.
\end{equation*}

Equipped with the ReLU network realizations of the logical operations underlying $\mathcal{I}_{\text{d}}$,

we are now ready to state the following universal representation result.

Consider the DMV algebra $\mathcal{I}_{\text{d}}$ in defn:DivisibleMV01. Let $n \in \mathbb{N}$. For each DMV term $\tau(x_1,\ldots \hspace{-0.02cm}, x_n)$ and its associated term function $\tau^{\mathcal{I}_{\text{d}}}:[0,1]^n \rightarrow [0,1]$, there exists a ReLU network $\Phi \in \mathcal{N}_{n,1}$ such that
\begin{equation*}
\Phi(x_1,\ldots \hspace{-0.02cm}, x_n) = \tau^{\Ical_{\dt}}(x_1,\ldots \hspace{-0.02cm}, x_n), \text{ for all } (x_1,\ldots \hspace{-0.02cm}, x_n)\in [0,1]^n.
\end{equation*}

The proof follows by realizing the logical operations appearing in the term function $\tau^{\mathcal{I}_{\text{d}}}$ through corresponding concatenations, according to lem:ReLUConcatenating, of ReLU networks implementing the operations $\oplus$ and $\odot$ as by lma:MinMax and noting that $\lnot x= 1-x, \delta_ix = \frac{1}{i}x$ are trivally ReLU networks with one layer.

We note that, in general, the ReLU network $\Phi$ in prop:ReLU realize divisible will be a properly deep network as it is obtained by concatenating the networks realizing the basic logical operations $\lnot$, $\oplus, \odot$, and $\{\delta_i \}_{i \in \mathbb{N}}$.

Now the ground has been prepared for the central result in this section, namely a universal ReLU network realization theorem for CA transition functions. Specifically, this will be effected by combining the connection between CA and DMV algebras established in sec:MVLogic with the ReLU network realizations of the logical operations in DMV algebras presented above.

Consider a CA with cellular space dimension $d \in \mathbb{N}$, neighborhood size $n \in \mathbb{N}$, state set $K = \{0,\frac{1}{k-1},\ldots \hspace{-0.02cm}, \frac{k-2}{k-1}, 1\}$ of cardinality $k \in \mathbb{N}, k\geq 2$, and transition function $f: K^{n} \rightarrow K$. There exists a ReLU network $\Phi \in \mathcal{N}_{n,1}$ satisfying

\begin{equation*}
\Phi(x_1,\ldots \hspace{-0.02cm}, x_n) = f(x_1,\ldots \hspace{-0.02cm}, x_n),\quad \text{ for all } (x_1,\ldots \hspace{-0.02cm},x_n)\in K^n.
\end{equation*}

By them:rationalMcnaughton, there exists a  DMV term $\tau(x_1,\ldots \hspace{-0.02cm}, x_n)$ whose corresponding term function
under the DMV algebra $\Ical_{\dt}$ satisfies

\begin{equation}\label{eq:15}
\tau^{\Ical_{\dt}}(x_1,\ldots \hspace{-0.02cm},x_n) = f(x_1,\ldots \hspace{-0.02cm},x_n),\quad \text{ for all } (x_1,\ldots \hspace{-0.02cm}, x_n) \in K^n.
\end{equation}

Application of prop:ReLU realize divisible then yields a ReLU network $\Phi \in \mathcal{N}_{n,1}$ such that
\begin{equation}\label{eq:16}
\Phi(x_1,\ldots \hspace{-0.02cm},x_n) = \tau^{\Ical_{\dt}}(x_1,\ldots \hspace{-0.02cm},x_n), \quad \text{ for all } (x_1,\ldots \hspace{-0.02cm},x_n)\in [0,1]^n.
\end{equation}
The proof is finalized by combining (<ref>) and (<ref>) to get

\begin{equation*}
\Phi(x_1,\ldots \hspace{-0.02cm},x_n) = \tau^{\Ical_{\dt}}(x_1,\ldots \hspace{-0.02cm},x_n) = f(x_1,\ldots \hspace{-0.02cm},x_n),
\end{equation*}

for all $(x_1,\ldots \hspace{-0.02cm},x_n) \in K^n$.

Consider a CA with cellular space $\mathbb{Z}$, the neighborhood set $\mathcal{E}=\{-1,0,1\}\subset \mathbb{Z}$, and state set $K$ in (<ref>) with $k\geq 3$. The transition rule $f$ is given by

\begin{equation*}\label{eq:kcolorRule}
f(x_{-1},x_0,x_1) =\begin{cases}
(x_0 + \frac{1}{k-1}) \mathbbm{1}_{ \{ x_{-1} = x_0 +\frac{1}{k-1} \text{ or } x_1 = x_0+\frac{1}{k-1}\}} &\quad \\
\hspace{0.7cm}  +\; x_0\mathbbm{1}_{\{ x_{-1} \neq x_0+\frac{1}{k-1}\} }\mathbbm{1}_{ \{ x_{1} \neq x_0+\frac{1}{k-1}\} }, & \text{ if } x_0 \in K \text{\textbackslash} \{1\}, \\
x_0\mathbbm{1}_{\{ x_{-1} \neq 0\} }\mathbbm{1}_{ \{ x_{1} \neq 0\} }, & \text{ if } x_0=1.
\end{cases}
\end{equation*}

which has a corresponding MV term under the MV algebra $\mathcal{A}_k$ in exmp:MVAlgebraAk. The resulting CA is the $k$-colored cyclic CA created by Bramson and Griffeath <cit.> to model particle systems.

§.§ Realizing the dynamical behavior of CA

We now turn our attention to the dynamical systems aspects of CA. Concretely, we show how CA evolution can be realized through RNNs. The basic idea underlying the construction we present is to suitably impose a recurrent structure on top of the ReLU network realizing the transition function of the CA under consideration. Concretely, we build on the RNN construction techniques developed in <cit.> and, for simplicity of exposition, consider the case of one-dimensional CA. An extension to the multi-dimensional case is quite readily obtainable using tools from multi-dimensional signal processing <cit.> and multi-dimensional RNNs <cit.>, <cit.>. For the sake of simplicity of exposition, however, we do not provide these extensions here.

An RNN is a discrete dynamical system mapping an input sequence to an output sequence, both possibly of infinite length, through—at each discrete time step—application of a feedforward neural network that updates a hidden-state vector and computes the next output signal sample <cit.>. The formal definition of an RNN is as follows.

For hidden-state vectordimension $m \in \mathbb{N}$, let $\Phi \in \mathcal{N}_{m+1,m+1}$ be a ReLU neural network. The recurrent neural network associated with $\Phi$ is the operator $\mathcal{R}_{\Phi}$ mapping input sequences $(x[z])_{z \in \mathbb{N}_0}$ in $\mathbb{R}$ to output sequences $(y[z])_{z \in \mathbb{N}_0}$ in $\mathbb{R}$ according to
\begin{equation}\label{eq:RNNdefinition}
\begin{pmatrix} y[z] \\ h[z] \end{pmatrix} = \Phi  \left( \begin{pmatrix}
x[z] \\ h[z-1]
\end{pmatrix}  \right), \quad z \in \mathbb{N}_0,
\end{equation}
where $h[z] \in \mathbb{R}^{m}$ is the hidden-state vector with initial value $h[-1] = 0_{m}$.

The key to emulating the evolution of a CA with an RNN is  the construction of an appropriate hidden-state vector in combination with a suitable ReLU network that encodes the CA transition function. Before stating the corresponding result, by way of preparation, we introduce a decomposition of
$\Phi$ in (<ref>) according to

\begin{equation}\label{eq:C}
\Phi = \begin{pmatrix}
\Phi^f \\
\Phi^h
\end{pmatrix},
\end{equation}

where $\Phi^f \in \Ncal_{m+1,1}$ is responsible for the computation of the output sample according to

\begin{equation}\label{eq:ypart}
y[z] = \Phi^f  \left( \begin{pmatrix}
x[z] \\ h[z-1]
\end{pmatrix}  \right), \quad z \in \mathbb{N}_0,
\end{equation}

and $\Phi^h \in  \Ncal_{m+1,m}$ effects the evolution of the hidden-state vector such that

\begin{equation}\label{eq:hpart}
h[z] = \Phi^h  \left( \begin{pmatrix}
x[z] \\ h[z-1]
\end{pmatrix}  \right), \quad z \in \mathbb{N}_0.
\end{equation}

The following theorem states the announced universal realization theorem for one-dimensional CA by RNNs.

Consider a CA with cellular space dimension $d =1$, neighborhood set $\mathcal{E}=\{0,-1,\ldots \hspace{-0.02cm}, -n+1\}$ of size $n \in \mathbb{N}, n\geq 2$, state set $K = \{0,\frac{1}{k-1},\ldots \hspace{-0.02cm}, \frac{k-2}{k-1}, 1\}$ of cardinality $k \in \mathbb{N}, k\geq 2$, and transition function $f: K^{n} \rightarrow K$ with associated CA map $F$. There exists an RNN that maps every configuration $(c[z])_{z\in \mathbb{Z}}$ over the cellular space $\mathbb{Z}$ to the next configuration $\left(F(c)[z]\right)_{z\in \mathbb{Z}}$ according to

\begin{equation}\label{eq:GlobalMap2}
F(c)[z] = f(c[z], c[z-1], \ldots \hspace{-0.02cm}, c[z-n+1]), \quad  \forall z\in \mathbb{Z}.
\end{equation}

Note that here we do not work with the general neighborhood set $\mathcal{E}=\{0, z_1,\ldots \hspace{-0.02cm}, z_{n-1}\} \subset \mathbb{Z}$, but rather consider the concrete neighborhood set $\mathcal{E} = \{0, -1, \ldots \hspace{-0.02cm}, -n+1\}$. This does not result in a loss of generality as the following argument shows. Take a general neighborhood set $\{0,z_1,\ldots \hspace{-0.02cm},z_{n-1}\} \subset \mathbb{Z}$ and expand it into a set of the form

\begin{equation}\label{eq:neighborhood}
\{\ell_1, \ell_1-1,\ldots \hspace{-0.02cm}, \ell_2+1,\ell_2\},
\end{equation}

where $\ell_1:=\max\{0,z_1,\ldots \hspace{-0.02cm},z_{n-1}\}$ and $\ell_2 := \min\{0,z_1,\ldots \hspace{-0.02cm},z_{n-1}\}$ by, when needed, adding void neighbors that do not affect the CA transition function. For example, consider the CA with neighborhood set $\mathcal{E}_*=\{1,0,-2\}$, transition function $f_*$, and corresponding global map

\begin{equation*}\label{eq:ExampleF}
F_*(c_*)[z] = f_*(c_*[z+1], c_*[z], c_*[z-2]), \quad  \text{ for } (c_*[z])_{z\in \mathbb{Z}}.
\end{equation*}

Then, expand $\mathcal{E}_*$ to $\mathcal{E}=\{1,0,-1,-2\}$ and make $f_*$ formally depend on the void neighbor $z=-1$ according to

\begin{multline*}
f(c_*[z+1], c_*[z],  c_*[z-1], c_*[z-2]) =\\f_*(c_*[z+1], c_*[z], c_*[z-2]),  \text{ for } (c_*[z])_{z\in \mathbb{Z}}.
\end{multline*}

These modifications lead to the global map

\begin{equation}\label{eq:examplec}
F_*(c_*)[z] =  f( c_*[z+1], c_*[z], c_*[z-1], c_*[z-2]), \quad \text{ for } (c_*[z])_{z\in \mathbb{Z}}.
\end{equation}

The general form (<ref>) is finally obtained upon making the substitution

\begin{equation*}
c[z]=c_*[z+1], \text{ for all }z\in \mathbb{Z},
\end{equation*}
to yield

\begin{equation*}
F(c)[z] =  f(c[z], c[z-1], c[z-2], c[z-3]), \quad \text{ for } (c[z])_{z\in \mathbb{Z}}.
\end{equation*}

For ease of exposition, we provide the proof for one-sided infinite sequences $(c[z])_{z\in \mathbb{N}_0}$ only. The general case follows mutatis mutandis, but is based on exactly the same ideas. The proof is constructive and will be effected by explicitly specifying the ReLU networks $\Phi^f$ and $\Phi^h$ in (<ref>) underlying the desired RNN. By (<ref>), the output sample $F(c)[z]$, which we identify with $y[z]$, is a function of the states of the cells $\{z, z-1, \ldots \hspace{-0.02cm}, z-n+1\}$. Therefore, equating the current input sample $x[z]$ in (<ref>) with $c[z]$, we choose the hidden-state vector $h[z-1]$ such that it stores the states of the neighbors $\{z-1,\ldots \hspace{-0.02cm}, z-n+1\}$ of the current input cell, i.e.,

\begin{equation}\label{eq:hiddenvec}
h[z-1] = \begin{pmatrix}
c[z-1] \\ c[z-2]\\ \vdots \\ c[z-n+1]
\end{pmatrix}.
\end{equation}

This leads to

\[
h[z] = \begin{pmatrix}
c[z] \\ c[z-1]\\ \vdots \\ c[z-n+2]
\end{pmatrix}
\]

and informs the choice of $\Phi^h$, which must satisfy

\begin{equation*}
\Phi^h \left(\begin{pmatrix}
c[z] \\ h[z-1]
\end{pmatrix} \right) = \Phi^h \left( \begin{pmatrix}
c[z] \\ c[z-1] \\\vdots \\c[z-n+2] \\ c[z-n+1]
\end{pmatrix}\right) =  \begin{pmatrix}
c[z] \\ c[z-1] \\\vdots \\ c[z-n+2]
\end{pmatrix} = h[z].
\end{equation*}

The evolution of the hidden-state vector hence proceeds by dropping the oldest value $c[z-n+1]$ and inserting the new value $c[z]$ at the top. Following the methodology developed in <cit.>, $\Phi^h$ can be realized by a two-layer ReLU network according to

\begin{equation*}
\Phi^h: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n-1}, \quad
\Phi^h(x)= (W^h_2 \circ \rho \circ W^h_1) (x),
\end{equation*}
with affine maps $W^h_1: \mathbb{R}^{n} \rightarrow \mathbb{R}^{2n}$, $W^h_2: \mathbb{R}^{2n} \rightarrow \mathbb{R}^{n-1} $ given by

\begin{equation}\label{eq:HiddenNetwork}
\begin{aligned}
W^h_1(x) &= \begin{pmatrix}
\mathbb{I}_{n} \\ - \mathbb{I}_{n}
\end{pmatrix} x, \quad \text{for } x \in \mathbb{R}^{n},\\
W^h_2(x) &= \begin{pmatrix}
\mathbb{I}_{n-1} \quad 0_{n-1} \quad -\mathbb{I}_{n-1} \quad 0_{n-1}
\end{pmatrix}x,\quad \text{for } x \in  \mathbb{R}^{2n}.
\end{aligned}
\end{equation}

We are left with the design of the network $\Phi^f$ which has to satisfy

\begin{equation*}
F(c)[z]   = \Phi^f  \left( \begin{pmatrix}
c[z] \\ c[z-1] \\ \vdots \\ c[z-n+1]
\end{pmatrix}  \right).
\end{equation*}

It follows by inspection of (<ref>) that this amounts to realizing the CA transition function $f$ through the ReLU network $\Phi^f$. Based on this insight, we can apply thm:CA-ReLU to conclude the existence of $\Phi^f$. In fact, the proof of thm:CA-ReLU spells out how $\Phi^f$ can be obtained explicitly by composing the logical operations in the DMV term associated with $f$. The proof is now completed by noting that the overall network $\Phi$ is obtained by applying lma:augmenting and lem:parallelizing to combine $\Phi^f$ and $\Phi^h$ according to (<ref>).

§ IDENTIFICATION OF CA LOGIC FROM TRAINED NEURAL NETWORKS

With suitably chosen transition functions and initial configurations, CA can simulate a plethora of dynamical behavior characteristics <cit.>. The inverse problem of deducing CA transition functions from observations of their evolution is, however, extremely difficult <cit.>, <cit.>.
Formally, this is known as the CA identification problem <cit.>: Given a finite sequence of consecutive configurations $\{c^1,c^2,\ldots \hspace{-0.02cm}, c^T\}$ collected during the evolution, construct a CA $(\mathbb{Z}^d, K, \mathcal{E}, f)$ whose associated CA map $F$ satisfies

\[
F(c^t) = c^{t+1}, \quad \text{ for } t=1, \ldots \hspace{-0.02cm}, T-1.
\]

Note that the initial configuration need not be recoverable from observed evolution traces <cit.>, <cit.>.

Whilst the cellular space $\mathbb{Z}^d$ and the state set $K$ can be read off directly from
evolution traces, the neighborhood set $\mathcal{E}$ can either be selected manually <cit.>, <cit.> or determined
through the application of specific criteria such as, e.g., mutual information <cit.>. The most challenging aspect of the CA identification problem resides in determining the transition function $f$.
Sakama et al. <cit.> construct decision trees to classify pairs of combinations of neighborhood states and new states collected from the evolution traces, which is then used to derive the transition function by genetic programming. In the case of binary state set $K = \{0,1\}$, Zhao et al. <cit.> exploit the fact that the transition function can be represented by a polynomial, which takes states of neighborhood cells as input variables, to reduce the problem from searching through all possible transition functions to computing the polynomial coefficients. Either write:

A detailed discussion of known CA identification methods can be found in <cit.>.

We next propose a novel approach to CA identification, namely reading out the DMV formula underlying the CA under consideration from an RNN
trained on its evolution traces. The focus will be on elucidating the fundamental principles of this idea.
Accordingly, we will not be concerned with the performance of specific RNN training algorithms, but will rather
assume that the feedforward network part $\Phi^{f}$ inside the RNN has been trained to achieve what is called in machine learning parlance “interpolation”, i.e.,

\begin{equation}\label{eq:interpolation}
\Phi^{f}(x_1,\ldots \hspace{-0.02cm},x_n) = f(x_1,\ldots \hspace{-0.02cm},x_n),
\end{equation}

for all $(x_1,\ldots \hspace{-0.02cm},x_n) \in K^n$.
This, of course, requires that the RNN being trained
have “seen” all possible combinations of neighborhood states and thereby the transition function $f$ on its entire domain $K^{n}$, a condition
met when the training configuration sequences are sufficiently long and their initial configurations exhibit sufficient richness <cit.>.

We will make this a standing assumption. In addition, the cellular space, the state set, and the neighborhood set are taken to be known a priori.

§.§ Interpolation

The procedure for extracting DMV terms underlying CA evolution data presented in Section <ref> below works off $\Phi^f$ as a (continuous) function mapping $[0,1]^{n}$ to $[0,1]$.
It turns out, however, that condition (<ref>) does not uniquely determine $\Phi^f$ on $[0,1]^n$ as there are—in general infinitely
many—different ways of interpolating $f(x_1,\ldots \hspace{-0.02cm},x_n)$ to a continuous piecewise linear function; recall that ReLU networks always realize continuous piecewise linear functions.

Since rationalMcNaughtonTheorem states that the truth functions associated with DMV terms under $\mathcal{I}_{\text{d}}$ are continuous piecewise linear functions with rational coefficients, we can constrain the weights of $\Phi^{f}$ to be rational. In fact, as the next two lemmata show, when interpolating CA transition functions, we can tighten this constraint even further,
namely to integer weights and rational biases. We first establish that CA transition functions can be interpolated to continuous piecewise linear functions with integer weights and rational biases[To be consistent with terminology used in the context of neural networks, for a linear function of the form $p(x_1, \ldots \hspace{-0.02cm}, x_{n}) = m_{1}x_1+\cdots+m_{n}x_{n}+b$, the coefficients $m_i$ will often be referred to as weights and $b$ as bias.].
Then, it is shown that such functions can be realized by ReLU networks with integer weights and rational biases.

as stated in lma:refinefc. Based on Condition  2 of rationalMcNaughtonTheorem, we can constrain that each linear piece of the interpolated functions has rational coefficients. In fact, when the underlying discrete function $f$ is a CA transition function, we can further tighten the constraint on coefficients, as stated in lma:refinefc.
In fact, generally, there exist infinitely many different $\Phi^f(x),\,x\in[0,1]^{n}$ satisfying (<ref>).
Before presenting our procedure for extracting DMV terms from trained networks, we discuss the interpolation condition (<ref>) and how it affects the extraction process.
With the goal of extracting DMV terms in mind, we only consider continuous piecewise linear interpolations of the discrete function $f$, because on the one hand, ReLU networks realize piecewise linear functions; on the other hand, the rational McNaughton theorem in rationalMcNaughtonTheorem states that the truth functions associated with DMV terms under $\mathcal{I}_{\text{d}}$ are continuous piecewise linear functions. Based on Condition  2 of rationalMcNaughtonTheorem, we can constrain that each linear piece of the interpolated functions has rational coefficients. In fact, when the underlying discrete function $f$ is a CA transition function, we can further tighten the constraint on coefficients, as stated in lma:refinefc.

Consider a CA with cellular space dimension $d \in \mathbb{N}$, neighborhood size $n \in \mathbb{N}$, state set $K = \{0,\frac{1}{k-1},\ldots , \frac{k-2}{k-1}, 1\}$ of cardinality $k \in \mathbb{N}, k\geq 2$, and transition function $f: K^{n} \rightarrow K$. There exists a function $f_{\text{c}}:[0,1]^n \rightarrow [0,1]$ with the following properties:

* $f_{\text{c}}(x_1, \ldots ,x_n) = f(x_1,\ldots \hspace{-0.02cm}, x_n)$,   for $(x_1, \ldots \hspace{-0.02cm},x_n)\in K^n$,

* $f_{\text{c}}$ is continuous with respect to the natural topology on $[0,1]^n$,

* there exist $\ell \in \mathbb{N}$ and linear functions

\begin{equation}%\label{eq:coeff}
p_j(x_1,\ldots \hspace{-0.02cm},x_n) = m_{j1}x_1+\cdots+m_{jn}x_n+\frac{b_j}{k-1}, \quad j=1,\ldots \hspace{-0.02cm},\ell,
\end{equation}

with $b_j,m_{j1},\ldots \hspace{-0.02cm},m_{jn} \in \mathbb{Z}$, for $j=1,\ldots \hspace{-0.02cm}, \ell$ such that, for every $x\in[0,1]^n$ there is an index $j\in \{1,\ldots \hspace{-0.02cm}, \ell\}$ with $f_{\text{c}}(x) = p_j(x)$.

We explicitly construct a function $f_{\text{c}}$ satisfying Properties 1-3. For $n=1$, one simply performs linear interpolation between
each pair of points $(\frac{i}{k-1}, f(\frac{i}{k-1}))$ and $(\frac{i+1}{k-1}, f(\frac{i+1}{k-1}))$, $i=0,\ldots \hspace{-0.02cm},k-2$, to obtain
an $f_{\text{c}}$ with the desired properties. While it is immediate that the resulting function $f_{\text{c}}$ satisfies Properties $1$ and $2$ and has the structure as demanded by Property $3$,
the condition $b_j,m_{j1},\ldots \hspace{-0.02cm},m_{jn} \in \mathbb{Z}$ will be verified summarily, for all $n$, at the end of the proof.

Before proceeding to the cases $n \geq 2$, we note that throughout the proof an $n$-simplex, $n\in \mathbb{N}$, with vertices $\{v_0,v_1,\ldots \hspace{-0.02cm},v_n\} \subset \mathbb{R}^n$ will be denoted by $\Delta^n = [v_0, v_1, \ldots \hspace{-0.02cm}, v_n]$.
We turn to the case $n=2$ and refer to fig:K2 for an illustration of the domain $K^2$ of $f$ for $k=3$. First, we divide $[0,1]^2$ into the squares $[\frac{i_1}{k-1}, \frac{i_1+1}{k-1}] \times [\frac{i_2}{k-1},\frac{i_2+1}{k-1}]$, $i_1,i_2 \in \{0,\ldots \hspace{-0.02cm},k-2\}$, followed by a subdivision of each of the resulting squares into $2$-simplices. There are precisely two possibilities for each of these subdivisions. To illustrate this, we consider the square $[0,\frac{1}{k-1}]\times [0,\frac{1}{k-1}]$ and note that
it can be split into two $2$-simplices according to either

\begin{align*}
\Delta^2_1 &= \left[(0,0), \left(0,\frac{1}{k-1}\right), \left(\frac{1}{k-1}, \frac{1}{k-1}\right)\right] \\
\Delta^2_2 &= \left[(0,0), \left(\frac{1}{k-1},0\right), \left(\frac{1}{k-1}, \frac{1}{k-1}\right)\right]
\end{align*}

\begin{align*}
\Delta^2_3 &= \left[(0,0), \left(0,\frac{1}{k-1}\right), \left(\frac{1}{k-1}, 0\right)\right] \\
\Delta^2_4 &= \left[\left(\frac{1}{k-1},\frac{1}{k-1}\right), \left(0,\frac{1}{k-1}\right), \left(\frac{1}{k-1}, 0\right)\right],
\end{align*}

as depicted in fig:subdide2.

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]

[->] (0,-5) – (0,-1.8) node[anchor=west] ;
[->] (0,-5) – (3.2,-5) node[anchor= west] ;
(0.07,-3.5) – (-0.07,-3.5) node[anchor= east] $1/2$;
(0.07,-2) – (-0.07,-2) node[anchor= east] $1$;

(1.5,-5.02) – (1.5,-4.98) node[anchor= north] $1/2$;
(3,-5.02) – (3,-4.98) node[anchor= north] ;
at (3,-5.2) $1$;
at (0.2,-5.2) $0$;
[black,fill=black] (0,-5) circle (.5ex);
[black,fill=black] (1.5,-5) circle (.5ex);
[black,fill=black] (3,-5) circle (.5ex);
[black,fill=black] (0,-3.5) circle (.5ex);
[black,fill=black] (1.5,-3.5) circle (.5ex);
[black,fill=black] (3,-3.5) circle (.5ex);
[black,fill=black] (0,-2) circle (.5ex);
[black,fill=black] (1.5,-2) circle (.5ex);
[black,fill=black] (3,-2) circle (.5ex);

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]

[->] (0,-5) – (0,-1.8) node[anchor=west] ;
[->] (0,-5) – (3.2,-5) node[anchor= west] ;
(0.07,-3.5) – (-0.07,-3.5) node[anchor= east] $1/2$;
(0.07,-2) – (-0.07,-2) node[anchor= east] $1$;

(1.5,-5.02) – (1.5,-4.98) node[anchor= north] $1/2$;
(3,-5.02) – (3,-4.98) node[anchor= north] ;
at (3,-5.2) $1$;
at (0.2,-5.2) $0$;
[black,fill=black] (0,-5) circle (.5ex);
[black,fill=black] (1.5,-5) circle (.5ex);
[black,fill=black] (3,-5) circle (.5ex);
[black,fill=black] (0,-3.5) circle (.5ex);
[black,fill=black] (1.5,-3.5) circle (.5ex);
[black,fill=black] (3,-3.5) circle (.5ex);
[black,fill=black] (0,-2) circle (.5ex);
[black,fill=black] (1.5,-2) circle (.5ex);
[black,fill=black] (3,-2) circle (.5ex);

The set $K^2$ for $k=3$.

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]

(0,-5) – (0,-2) node[anchor=west] ;
(0,-2) – (3,-2) node[anchor= west] ;
(0,-5) – (3,-5) node[anchor= west] ;
(3,-5) – (3,-2) node[anchor= west] ;
(0,-5) – (3,-2) node[anchor= west] ;

at (3,-5.3) $(\frac{1}{k-1},0)$;
at (0.0,-5.3) $(0,0)$;
at (3.0,-1.7) $(\frac{1}{k-1},\frac{1}{k-1})$;
at (0.0,-1.7) $(0, \frac{1}{k-1})$;
[black,fill=black] (0,-5) circle (.5ex);

[black,fill=black] (3,-5) circle (.5ex);
[black,fill=black] (0,-2) circle (.5ex);
[black,fill=black] (3,-2) circle (.5ex);

(5,-5) – (5,-2) node[anchor=west] ;
(5,-2) – (8,-2) node[anchor= west] ;
(5,-5) – (8,-5) node[anchor= west] ;
(8,-5) – (8,-2) node[anchor= west] ;
(5,-2) – (8,-5) node[anchor= west] ;

at (8,-5.3) $(\frac{1}{k-1},0)$;
at (5.0,-5.3) $(0,0)$;
at (8.0,-1.7) $(\frac{1}{k-1},\frac{1}{k-1})$;
at (5.0,-1.7) $(0, \frac{1}{k-1})$;
[black,fill=black] (5,-5) circle (.5ex);

[black,fill=black] (8,-5) circle (.5ex);
[black,fill=black] (5,-2) circle (.5ex);
[black,fill=black] (8,-2) circle (.5ex);

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]

(0,-5) – (0,-2) node[anchor=west] ;
(0,-2) – (3,-2) node[anchor= west] ;
(0,-5) – (3,-5) node[anchor= west] ;
(3,-5) – (3,-2) node[anchor= west] ;
(0,-5) – (3,-2) node[anchor= west] ;

at (3,-5.3) $(\frac{1}{k-1},0)$;
at (0.0,-5.3) $(0,0)$;
at (3.0,-1.7) $(\frac{1}{k-1},\frac{1}{k-1})$;
at (0.0,-1.7) $(0, \frac{1}{k-1})$;
[black,fill=black] (0,-5) circle (.5ex);

[black,fill=black] (3,-5) circle (.5ex);
[black,fill=black] (0,-2) circle (.5ex);
[black,fill=black] (3,-2) circle (.5ex);

(5,-5) – (5,-2) node[anchor=west] ;
(5,-2) – (8,-2) node[anchor= west] ;
(5,-5) – (8,-5) node[anchor= west] ;
(8,-5) – (8,-2) node[anchor= west] ;
(5,-2) – (8,-5) node[anchor= west] ;

at (8,-5.3) $(\frac{1}{k-1},0)$;
at (5.0,-5.3) $(0,0)$;
at (8.0,-1.7) $(\frac{1}{k-1},\frac{1}{k-1})$;
at (5.0,-1.7) $(0, \frac{1}{k-1})$;
[black,fill=black] (5,-5) circle (.5ex);

[black,fill=black] (8,-5) circle (.5ex);
[black,fill=black] (5,-2) circle (.5ex);
[black,fill=black] (8,-2) circle (.5ex);

Two different subdivisions.

Fixing the subdivisions of all squares, on each $2$-simplex $\Delta^2 =[v_0,v_1,v_2]$, we define a linear function $p_{\Delta^2}$ that interpolates $f(x)$, i.e., $p_{\Delta^2}(x) = f(x)$, for $x\in \{v_0,v_1,v_2\}$. As $\{v_0,v_1,v_2\}$ is affinely independent, $p_{\Delta^2}$ is uniquely determined <cit.>. We then stitch the resulting linear pieces together as follows.

For each point $x\in [0,1]^2$, if it falls into the interior of some $2$-simplex $\Delta^2$, we set $f_{\text{c}}(x) =p_{\Delta^2}(x)$; if $x$ resides in the intersection of two or more $2$-simplices, then $f_{\text{c}}(x)$ is taken to have the (shared) value of the interpolating functions associated with the intersecting simplices evaluated at $x$. The resulting function $f_{\text{c}}$ satisfies Properties $1$ and $2$ and exhibits the structure as demanded by Property $3$.

For $n\geq 3$, we first divide the unit cube $[0,1]^n$ into the $n$-dimensional cubes $[\frac{i_1}{k-1}, \frac{i_1+1}{k-1}] \times \cdots \times [\frac{i_n}{k-1}, \frac{i_n+1}{k-1}], i_1,\ldots \hspace{-0.02cm},i_n \in \{0,\ldots \hspace{-0.02cm},k-2\}$.

Each of the resulting smaller cubes is then subdivided following the procedure described in <cit.>, which
divides, e.g.,

$\Delta^{n-1} \times [0,\frac{1}{k-1}] $ into $n$-simplices as follows. Let $\Delta^{n-1} \times \{0\} = [v_0,\ldots \hspace{-0.02cm},v_{n-1}]$ and $\Delta^{n-1} \times \{\frac{1}{k-1}\} = [w_0,\ldots \hspace{-0.02cm},w_{n-1}]$ and note that the coordinates of $v_j \in \mathbb{R}^{n}$ and $w_j \in \mathbb{R}^{n}$ coincide in the first $n-1$ indices.
Then, $\Delta^{n-1} \times [0,\frac{1}{k-1}]$ is given by the union of the $n$-simplices $[v_0,\ldots \hspace{-0.02cm},v_j,w_j,\ldots \hspace{-0.02cm},w_{n-1}]$, $j=0,\ldots \hspace{-0.02cm},n-1$, each intersecting the next one in an $(n-1)$-simplex face. The result of this procedure yields a
subdivision of $[0,1]^n$ into $n$-simplices with vertices $K^n$.
As in the case $n=2$, an interpolating linear function is uniquely determined on each of these $n$-simplices <cit.>; we stitch the resulting linear pieces together to obtain a function $f_{\text{c}}$ satisfying Property $1$ and exhibiting the structure as demanded by Property $3$.

As $f_{\text{c}}$ constitutes a simplex interpolation of the discrete function $f$ on the cube $[0,1]^n$, it is continuous <cit.>, <cit.>, thereby satisfying Property $2$.

We hasten to add that, for $n\geq 2$, the set of $n$-simplices resulting from the subdivision procedure employed here is not unique, see fig:subdide2 for an illustration in the case $n=2$. But, each fixed subdivision into $n$-simplices uniquely determines a continuous piecewise linear function $f_{\text{c}}$ that interpolates $f$.

It remains to prove that the coefficients of the constituent linear polynomials

\begin{equation}%\label{eq:coeff}
p_j(x_1,\ldots \hspace{-0.02cm},x_n) = m_{j1}x_1+\cdots+m_{jn}x_n+\frac{b_j}{k-1}, \quad j=1,\ldots \hspace{-0.02cm},\ell,
\end{equation}

of $f_{\text{c}}$, indeed, satisfy $b_j,m_{j1},\ldots \hspace{-0.02cm},m_{jn} \in \mathbb{Z}$, for $j=1,\ldots \hspace{-0.02cm}, \ell$.

As already mentioned, the coefficients of each $p_j$ are uniquely determined and obtained by solving a system of linear equations induced
by the vertex coordinates of the $n$-simplex $p_j$ resides on. We now show that the vertices $v_0,v_1,\ldots \hspace{-0.02cm},v_n\subset K^n$ of
every $n$-simplex in $[0,1]^n$ resulting from the procedure described above satisfy, after reordering if necessary,

\begin{equation}\label{eq:differ}
||v_{i}-v_{i+1}||_1 = \frac{1}{k-1}, \quad \text{for } i=0,\ldots \hspace{-0.02cm},n-1.
\end{equation}

In other words, each of the (reordered) vertices differs from the next one in only one position, specifically by $\pm \frac{1}{k-1}$. In the case
$n=1$ the condition (<ref>) holds without reordering. For $n=2$, no reordering is needed in $\Delta^2_1$ and $\Delta^2_2$, while in each of
$\Delta^2_3$ and $\Delta^2_4$ we only need to swap the second vertex with the first one. For general $n \geq 3$,

we make use of the fact that,

for each $n$-simplex $\Delta^{n}$ of $[0,1]^{n}$ constructed according to the procedure described above, there exists an $(n-1)$-simplex $\Delta^{n-1} = [v_0,v_1,\ldots \hspace{-0.02cm},v_{n-1}]$ such that

\[
\Delta^{n} =  \begin{bmatrix}
\begin{pmatrix}
v_0 \\ \frac{j}{k-1}
\end{pmatrix}, \ldots \hspace{-0.02cm}, \begin{pmatrix}
v_i \\ \frac{j}{k-1}
\end{pmatrix}, \begin{pmatrix}
v_i \\ \frac{j+1}{k-1}
\end{pmatrix},\ldots \hspace{-0.02cm},\begin{pmatrix}
v_{n-1} \\ \frac{j+1}{k-1}
\end{pmatrix}
\end{bmatrix},
\]

for some $j\in \{0,\ldots \hspace{-0.02cm},k-2\}$ and with $i \in \{0,\ldots \hspace{-0.02cm},n-1\}$. If $\Delta^{n-1}$ satisfies (<ref>), it follows directly that
$\Delta^{n}$ satisfies (<ref>). Now, noting that for $n=3$, $\Delta^{2}$ satisfies (<ref>) as already established, and
proceeding by induction across $n$ establishes the desired statement.

Finally, assume that the linear piece

\begin{equation*}
m_1x_1+m_2x_2+\cdots +m_nx_n +m_0
\end{equation*}

of $f_{\text{c}}$ is obtained, w.l.o.g. thanks to (<ref>), by interpolating $n+1$ points at

\begin{equation*}
\begin{aligned}
\biggl( &\frac{a_1}{k-1}, \frac{a_2}{k-1}, \ldots \hspace{-0.02cm}, \frac{a_n}{k-1}, \frac{b_0}{k-1} \biggl)  \\
\biggl( &\frac{a_1+1}{k-1}, \frac{a_2}{k-1}, \ldots \hspace{-0.02cm}, \frac{a_n}{k-1}, \frac{b_1}{k-1} \biggl)\\
\biggl( &\frac{a_1+1}{k-1}, \frac{a_2+1}{k-1}, \ldots \hspace{-0.02cm}, \frac{a_n}{k-1}, \frac{b_2}{k-1} \biggl)\\
&\quad \quad \quad \quad \quad \quad \vdots \\
\biggl( &\frac{a_1+1}{k-1}, \frac{a_2+1}{k-1}, \ldots \hspace{-0.02cm}, \frac{a_n+1}{k-1}, \frac{b_n}{k-1} \biggl),
\end{aligned}
\end{equation*}

with $a_1,\ldots \hspace{-0.02cm}, a_n,b_0,\ldots  \hspace{-0.02cm},b_n\in \{0,1,\ldots \hspace{-0.02cm}, k-1\}$.
The coefficients $m_0, m_1,\ldots \hspace{-0.02cm}, m_n$ hence have to satisfy the following system of linear equations

\begin{equation}\label{eq:Gauss}
\begin{aligned}
\frac{a_1}{k-1}m_1+ \frac{a_2}{k-1}m_2\, + &\cdots+ \frac{a_n}{k-1}m_n + m_0 = \frac{b_0}{k-1} \\
\frac{a_1+1}{k-1}m_1+ \frac{a_2}{k-1}m_2\, + &\cdots+ \frac{a_n}{k-1}m_n + m_0 = \frac{b_1}{k-1} \\
\frac{a_1+1}{k-1}m_1+ \frac{a_2+1}{k-1}m_2\, + &\cdots+ \frac{a_n}{k-1}m_n + m_0 = \frac{b_2}{k-1} \\
&\;\; \vdots \\
\frac{a_1+1}{k-1}m_1+ \frac{a_2+1}{k-1}m_2\, + &\cdots+ \frac{a_n+1}{k-1}m_n + m_0 = \frac{b_n}{k-1}.
\end{aligned}
\end{equation}

Subtracting the first row in (<ref>) from the second and multiplying the result by $k-1$, we obtain

\begin{equation*}
m_1 = b_1-b_0 \in \mathbb{Z}.
\end{equation*}

Continuing likewise by subtracting the first row from the third yields

\begin{equation*}
m_2  = b_2-b_0 - m_1 \in  \mathbb{Z}
%      & \;\; \vdots \\
%m_n &= b_{n}-b_0,
\end{equation*}

and so forth, establishing that $m_i \in \mathbb{Z}$, for all $i=1,2,\ldots \hspace{-0.02cm},n$.

Substituting back into any of the equations in (<ref>) establishes that $m_0$ is of the desired form as well and thereby finalizes the proof.

We next show that functions $f_{\text{c}}$ satisfying Properties $1$-$3$ in lma:refinefc can, indeed, be realized by ReLU networks with integer
weights and biases in $\mathbb{Q}_{k}:=\{\frac{b}{k-1}:b\in \mathbb{Z}\}$.

Let $n\in \mathbb{N}$ and let $f_{\text{c}}:[0,1]^n\rightarrow [0,1]$ be a continuous piecewise linear function with linear pieces $p_1, \ldots \hspace{-0.02cm}, p_\ell, \ell\in \mathbb{N}$, of the form

\begin{equation}\label{eq:linearpj1}
p_j(x_1, \ldots \hspace{-0.02cm}, x_{n}) =m_{j1}x_1+\cdots+m_{jn}x_{n}+ \frac{b_j}{k-1}, \quad  j = 1,\ldots \hspace{-0.02cm},\ell,
\end{equation}
where $k \geq 2$ and $m_{j1},\ldots \hspace{-0.02cm}, m_{jn}, b_j \in \mathbb{Z}$, for $j=1,\ldots \hspace{-0.02cm}, \ell$. There exists a ReLU network $\Phi \in \mathcal{N}_{n,1}$ with integer weights and biases in $\mathbb{Q}_k$, satisfying $\Phi(x) = f_{\text{c}}(x)$, for all $x\in [0,1]^n$.

By lem:piecewise, the function $f_{\text{c}}$ can be written in terms of $p_1,\dots,p_\ell$ through “$\min$” and “$\max$” operations in the form

\begin{equation}\label{eq:gamma}
f_{\text{c}} = \max_I \min_{J} p_{j},
\end{equation}

where $I, J \subset \{1,\ldots \hspace{-0.02cm},\ell\}$ are index sets. The linear pieces $p_j$, for $j=1,\ldots \hspace{-0.02cm},\ell$, in (<ref>) can straightforwardly be realized by single-layer ReLU networks with integer weights $m_{ji}$ and biases $b_{j}/(k-1)$. Next, we note that the “$\min$” operation can be implemented by a ReLU network with integer weights and biases equal to $0$ according to

\begin{equation*}
\min\{x_1,x_2\} = \rho(x_1) - \rho(-x_1) - \rho(x_1-x_2), \text{ for } x_1,x_2\in \mathbb{R}.
\end{equation*}

\begin{equation*}
\min\{x_1,x_2,x_3\} = \min\{x_1, \min\{x_2,x_3\}\}, \text{ for } x_1,x_2,x_3 \in \mathbb{R},
\end{equation*}

it follows that the minimum over any number of variables can be realized by nesting minimum operations. Hence, inspection of the proof of
<cit.> reveals that there exists a ReLU network with integer weights and biases equal to $0$ computing the minimum over any number of variables. Likewise, we have

\begin{equation*}
\max\{x_1,x_2\} = \rho(x_1) - \rho(-x_1) +\rho(x_2-x_1), \text{ for } x_1,x_2\in \mathbb{R},
\end{equation*}

\begin{equation*}
\max\{x_1,x_2,x_3\} = \max\{x_1, \max\{x_2,x_3\}\}, \text{ for } x_1,x_2,x_3 \in \mathbb{R},
\end{equation*}

which analogously establishes that the operation taking the maximum over any number of variables can be realized by a ReLU network with integer weights and biases equal to $0$. Based on (<ref>), $f_{\text{c}}$ can be written as the composition of ReLU networks with integer weights and biases either equal to $0$ (for the min and max operations) or equal to $b_{j}/(k-1)$ (for the $p_{j}$). Inspection of the proof of <cit.> then shows that the resulting overall ReLU network will have integer weights and biases in $\mathbb{Q}_k$, as desired.

The proof of lma:refinefc, combined with lma:intReLU, reveals that even if we constrain the network $\Phi^{f}$ interpolating the transition function $f$ to have integer weights and biases in $\mathbb{Q}_k$, the network will not be unique. This is a consequence of the
$n$-simplices underlying the linear pieces $p_{j}$ constructed in the proof of lma:refinefc not being unique, see e.g. fig:subdide2 for an illustration of this phenomenon in the case $n=2$.
Correspondingly the DMV formulae extracted from the networks realizing the functions $f_{\text{c}}$ based on different choices of $n$-simplices,
following the procedure described in Section <ref> below, will be functionally different. Beyond this nonuniqueness, we can even get interpolating ReLU networks with integer weights and biases in $\mathbb{Q}_{k}$ that do not effect simplex interpolation.
This aspect is illustrated in the next example.

Consider a CA with cellular space dimension $d=1$, neighborhood size $n=2$, state set $\{0,1/2,1\}$ with $k=3$, and transition function $f :\{0,1/2,1\}^2 \rightarrow \{0,1/2,1\}$ specified as follows:

$x_{-1}\,x_0$    (0, 0)     (0, 1/2)     (0, 1)     (1/2, 0)     (1/2, 1/2)
$f(x_{-1},\,x_0)$     0     0     0     0     1

$x_{-1}\,x_0$     (1/2, 1)     (1, 0)     (1, 1/2)     (1, 1)
$f(x_{-1},\,x_0)$      0     0    0     0

fig:eg2d shows two different ReLU networks $\Phi^f_1,\Phi^f_2:[0,1]^2 \rightarrow [0,1]$, both with integer weights and biases in $\mathbb{Q}_{k}$, interpolating the transition function. We remark that $\Phi^f_1$ effects simplex interpolation, while $\Phi^f_2$ does not.

Now, applying the DMV formula extraction algorithm introduced in sec:reading-out below yields the DMV term associated with $\Phi^f_1$ as

\[
\tau_1 = (x_{-1}\oplus x_{-1}) \wedge (\lnot x_{-1} \oplus \lnot x_{-1}) \wedge (x_0 \oplus x_0) \wedge ( \lnot x_0 \oplus \lnot x_0),
\]

and that associated with $\Phi^f_2$ according to

\[
\tau_2=(x_{-1} \oplus x_{-1}) \wedge (\lnot x_{-1} \oplus \lnot x_{-1}) \wedge (x_0 \oplus x_0 \oplus x_0) \wedge ( \lnot x_0 \oplus \lnot x_0 \oplus \lnot x_0 ).
\]

These two DMV terms are algebraically and functionally different, but the associated term functions under $\mathcal{I}_{\text{d}}$ coincide on $K^2$, i.e.,

\[
\tau_1^{\mathcal{I}_{\text{d}}}(x_{-1},x_0)=\tau_2^{\mathcal{I}_{\text{d}}}(x_{-1},x_0) =f(x_{-1},x_0),
\]

for $(x_{-1},x_0) \in \{0,1/2,1\}^2$.

We shall consider all term functions that coincide with a given CA transition function on $K^n$ to constitute
an equivalence class in the sense of faithfully describing the logical behavior of the underlying CA.

Different ReLU networks interpolating $f$.

Consider a CA with cellular space dimension $d=1$, neighborhood size $n=2$, state set $\{0,1/2,1\}$ with $k=3$, and transition rule $f :\{0,1/2,1\}^2 \rightarrow \{0,1/2,1\}$ specified as follows:

$x_{-1}\,x_0$    (0, 0)     (0, 1/2)     (0, 1)     (1/2, 0)     (1/2, 1/2)
$f(x_{-1},\,x_0)$     0     0     0     0     1

$x_{-1}\,x_0$     (1/2, 1)     (1, 0)     (1, 1/2)     (1, 1)
$f(x_{-1},\,x_0)$      0     0    0     0

There exist different ReLU networks $\Phi^f:[0,1]^2 \rightarrow [0,1]$, with integer weights and biases in $\mathbb{Q}_{k}$, satisfying $\Phi^f(x_{-1},x_0)  = f(x_{-1},x_0)$. See fig:eg2d for two examples, namely $\Phi^f_1$ and $\Phi^f_2$. The DMV term associated with $\Phi^f_1$ is

\[
\tau_1 = (x_{-1}\oplus x_{-1}) \wedge (\lnot x_{-1} \oplus \lnot x_{-1}) \wedge (x_0 \oplus x_0) \wedge ( \lnot x_0 \oplus \lnot x_0),
\]

which is different from the DMV term associated with $\Phi^f_2$ given by

\[
\tau_2=(x_{-1} \oplus x_{-1}) \wedge (\lnot x_{-1} \oplus \lnot x_{-1}) \wedge (x_0 \oplus x_0 \oplus x_0) \wedge ( \lnot x_0 \oplus \lnot x_0 \oplus \lnot x_0 ).
\]

But the associated term functions under $\mathcal{I}_{\text{d}}$ satisfy

\[
\tau_1^{\mathcal{I}_{\text{d}}}(x_{-1},x_0)=\tau_2^{\mathcal{I}_{\text{d}}}(x_{-1},x_0) =f(x_{-1},x_0),
\]

for $(x_{-1},x_0) \in \{0,1/2,1\}^2$.

Two different interpolations of $f$.

There are two further sources of potential nonuniqueness discussed next.

rationalMcNaughtonTheorem states that the truth functions associated with DMV terms under $\mathcal{I}_{\text{d}}$ are continuous piecewise linear functions with the individual linear pieces exhibiting rational coefficients, we constrain the coefficients of $\Phi^{f}$ to be rational. In fact, as the next result shows, when interpolating CA transition functions, we can tighten the constraint on the coefficients even further, namely to integer weights and
rational biases.

prescribing the piecewise linear function $\Phi^{f}$ on $K^n$ only does not lead to a unique interpolated
function $f_{\text{c}}$.

Correspondingly the DMV formulae extracted from these networks, following the procedure described in Section <ref> below, will be functionally different. We shall consider all these formulae to be logically equivalent in the sense of faithfully characterizing the logical behavior underlying the CA we want to learn and leave the study of the implications of this observation for further work.

lma:refinefc states that for every CA transition function, there always exists a interpolation satisfying Conditions 1$\sim$3 in lma:refinefc, and, in particular, have integer coefficients and rational biases. These conditions do not lead to unique interpolations though. One reason is that, as already mentioned, the subdivisions of the unit cube are not unique. But we emphasize that all the DMV terms extracted from interpolations satisfying Conditions 1$\sim$3 in lma:refinefc are logically equivalent under the resolution $k$, in that their associated truth functions under $\mathcal{I}_{\text{d}}$ coincide with the underlying CA transition function at all possible inputs $K^n$. We would come back to this point with illustrating examples after introducing our extraction procedure.

While the aspects discussed so far exhibit nonuniqueness that can either easily be mitigated—by enforcing a specific structure on the RNN to be trained—or lead to structural differences in the ReLU network $\Phi^{f}$ while preserving functional equivalence, there is a last aspect that points to proper nonuniqueness. Specifically, recalling that the network $\Phi^{f}$ is assumed to be trained to match the transition function $f$ on $K^{n}$ according to (<ref>), we note that

written by Yani ends

Before we proceed to the development of the algorithm for reading out the DMV formula corresponding to the trained network $\Phi^{f}$, we address the question of whether this is possible, in general, in a unique fashion.

For all resolutions $k\in \mathbb{N}, k\geq 2$, the transition table follows readily from the DMV term simply by direct evaluation. Take a simple example, namely $\tau = x_{-1} \oplus x_0$. Consider a CA with cellular space dimension $d=1$ and neighborhood set $\mathcal{E} = \{-1,0\}$. We first take the state set to be $K = \{0,1\}$. The tabular CA transition rule corresponding to $\tau$ is as follows

$x_{-1}\,x_0$    11     10     01     00
0pt15pt $\tau^{\mathcal{I}}(x_{-1}, x_0)$     1     1     1     0

For the state set $K = \{0,1/2,1\}$, the resulting CA transition rule in tabular form reads:

$x_{-1}\,x_0$    00     (0, 1/2)     (0,1)     (1/2, 0)     (1/2,1/2)
0pt15pt  $\tau^{\mathcal{I}}(x_{-1}, x_0)$     0     1/2     1     1/2     1

If we further increase the resolution of the state set to $K = \{0,\frac{1}{k-1},\ldots \hspace{-0.02cm},$ $\frac{k-2}{k-1},1\}$, for $k=4,5,\ldots$, we generate bigger CA transition tables (of sizes $k^2$). In other words, we can find infinitely many different CA whose transition rules have the same DMV term, and, correspondingly, can be realized by the same neural network according to prop:ReLU realize divisible. Conversely, starting from a neural network $\Phi$ that has learned some CA transition function, the network realizes a mapping $\Phi: [0,1]^n \rightarrow [0,1], n \in \mathbb{N}$, which is the term function of some DMV term. Hence, all DMV terms $\tau$ extracted from the network $\Phi$, such that $\tau^\mathcal{I}_{\text{d}} = \Phi$, are equivalent in the sense of having the same term function. In other words, different values of $k$ reflects different levels of resolution of discretizing the continuous term function. The two examples above with $k=2$ and $k=3$ are associated with the same DMV term $\tau = x_{-1} \oplus x_0$ and correspondingly the same term function $ \tau^{\mathcal{I}_{\text{d}}}(x_{-1}, x_0) = \min\{1,x_{-1}+x_0\}$, with $k=3$ resulting in a higher resolution of discretization and a larger transition table.

On the other hand, neural networks trained on traces of CA evolutions are likely to be over-resolved compared to the CA resolution $k$. Specifically, the trained ReLU network, which realizes a piecewise linear mapping, might not be the linear interpolation of the underlying CA transition function. For example, a network trained on sample data generated by the discrete function $f$ in  fig:contrast(a) might end up realizing the function $f^*$ in fig:contrast(b), when we actually want the network to realize the function $f^{**}$ in fig:contrast(c) instead. Intuitively, the network realizing $f^*$ in (b) should be more complex than the one realizing  $f^{**}$ in (c), and, furthermore,  the DMV term extracted from (b) are likely to have higher complexity than the one extracted from (c). How the complexities of neural networks and DMV terms interplay with the CA resolution $k$, and how to control the network complexity during training such that it is not over-resolved compared to the CA resolution $k$, will be our future research direction.

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]
[->] (0,0) – (0,3.2) node[anchor=west] $f$;
[->] (0,0) – (3.2,0) node[anchor= west] $x$;
(0.07,1.5) – (-0.07,1.5) node[anchor= east] $1/2$;

(0.07,3) – (-0.07,3) node[anchor= east] $1$;

(1.5,-0.02) – (1.5,0.02) node[anchor= north] $1/2$;
(3,-0.02) – (3,0.02) node[anchor= north] ;
at (3,-0.2) $1$;
at (0.2,-0.3 cm) $0$;
[black,fill=black] (0,0) circle (.5ex);
[black,fill=black] (1.5,3) circle (.5ex);
[black,fill=black] (3,0) circle (.5ex);

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]

[->] (0,0) – (0,3.2) node[anchor=west] $f^*$;
[->] (0,0) – (3.2,0) node[anchor= west] $x$;
[black,fill=black] (0,0) circle (.5ex);

[black,fill=black] (1.5,3) circle (.5ex);

[black,fill=black] (3,0) circle (.5ex);
(0,0) –(0.5,2);
(0.5,2) –(1,1);
(1,1) – (1.5,3);
(1.5,3) – (2,2.5);
(2,2.5) –  (2.5,1.5);
(3,0) –  (2.5,1.5);

(0.07,1.5) – (-0.07,1.5) node[anchor= east] $1/2$;

(0.07,3) – (-0.07,3) node[anchor= east] $1$;

(1.5,-0.02) – (1.5,0.02) node[anchor= north] $1/2$;
(3,-0.02) – (3,0.02) node[anchor= north] ;
at (3,-0.2) $1$;
at (0.2,-0.3 cm) $0$;

[x=1.5cm, y=1.5cm, >=stealth, scale=1, transform shape]
[->] (0,0) – (0,3.2) node[anchor=west] $f^{**}$;
[->] (0,0) – (3.2,0) node[anchor= west] $x$;
(0.07,1.5) – (-0.07,1.5) node[anchor= east] $1/2$;

(0.07,3) – (-0.07,3) node[anchor= east] $1$;

(1.5,-0.02) – (1.5,0.02) node[anchor= north] $1/2$;
(3,-0.02) – (3,0.02) node[anchor= north] ;
at (3,-0.2) $1$;
at (0.2,-0.3 cm) $0$;
[black,fill=black] (0,0) circle (.5ex);
[black,fill=black] (1.5,3) circle (.5ex);
[black,fill=black] (3,0) circle (.5ex);

(0,0) – (1.5,3);
(1.5,3) –  (3,0);

(a) discrete function $f$ (b) over-resolved function $f^{**}$ (c) linearly interpolated function $f^*$

Linear interpolation of the two CA transition tables above results in the continuous function $f_{\text{c}}=\min \{1,x_{-1}+x_0\}$. In general, continuous functions obtained by linearly interpolating CA transition functions $f:K^n \rightarrow K$, for all $k \in \mathbb{N}$, are rational McNaughton functions. For the construction of DMV terms from rational McNaughton functions, three different approaches have been proposed in the literature. The method reported in <cit.> is effected by multiplying the rational McNaughton function by a large enough integer to make it a McNaughton function. The other two methods are proposed by Gerla et al. <cit.>, where the DMV term is constructed from Schauder hats <cit.> by subdividing the unit cube into simplices, or from the decomposition the rational McNaughton function into linear polynomials connected by the “min” and “max” operators. Whereas it is straightforward to specify a rational McNaughton function on the interval $[0,1]$, it is not easy to effectively characterize it on the multi-dimensional domain $[0,1]^n, n\in \mathbb{N}$. Compared to <cit.>, <cit.>, the proposed method works with, instead of a rational McNaughton function, a network that has learned the function. Logical operations are reflected through the networks's affine transformations and compositions, and the proposed procedure shall effect by translating these affine transformations and compositions into DMV terms.

§.§ Uniqueness properties of DMV formula extraction

We first comment on a structural property pertaining to the RNN emulating the CA evolution. Recall the decomposition of $\Phi$ in (<ref>). Our construction built on the idea of having the hidden-state vector store the neighbors of the current input cell, leading to the specific form of $\Phi^{h}$, and making the

subnetwork $\Phi^f$ be responsible for the computation of the output samples. But this separation need not be the only way the RNN can realize the CA evolution. The resulting ambiguity is, however, easily eliminated by enforcing the split of $\Phi$ according to (<ref>) on the RNN to be trained. In practice, this means that
the weights of $\Phi^h$ are fixed according to (<ref>) and only the weights of $\Phi^f$ are trained. In the light of lma:refinefc and lma:intReLU, the training algorithm has to enforce integer weights and biases in $\mathbb{Q}_k$ on $\Phi^f$.
This can be done using techniques such as those described in e.g. <cit.>, <cit.>.

At the end of the training process, the DMV formula underlying the CA can be read out from $\Phi^f$ using the algorithm presented in Section <ref> below. Recall that given what was said after lma:intReLU, this formula will not be unique in general.

The second aspect we wish to dwell on here revolves around the fact that algebraically different DMV formulae can be functionally equivalent, just like what one has in classical algebra.

To demonstrate the matter, we analyze a simple example in Boolean algebra, namely the elementary CA of rule 30 with transition function specified in DNF in (<ref>). A functionally equivalent, but algebraically different, expression for this transition rule is given by

\begin{multline}\label{eq:rule30logic2}
f_{30} = (x_{-1}\odot  \lnot x_0 \odot \lnot x_1) \oplus  (\lnot x_{-1}\odot   x_0 \odot  x_1) \oplus \\(\lnot x_{-1}\odot   x_0 \odot  \lnot x_1) \oplus (\lnot x_{-1}\odot   \lnot x_0 \odot  x_1).
\end{multline}

This shows that even if we restrict ourselves to DNF, the algebraic expression will not be unique in general.

Moreover, we can also express rule 30 in conjunctive normal form (CNF), i.e., as a concatenation of a finite number of clauses linked by the Boolean $\odot$ operation, where each clause consists of a finite number of variables or negated variables connected by the Boolean $\oplus$ operation, e.g.,

\begin{equation}\label{eq:rule30logic3}
f_{30} = ( x_{-1}\oplus   x_0 \oplus  x_1) \odot ( \lnot x_{-1}\oplus  \lnot x_1 ) \odot (\lnot x_{-1}\oplus  \lnot x_0),
\end{equation}

\begin{equation}\label{eq:rule30logic4}
f_{30} = ( x_{-1}\oplus   x_0 \oplus  x_1) \odot ( \lnot x_{-1}\oplus \lnot   x_1 ) \odot ( \lnot x_{-1}\oplus \lnot  x_0 \oplus  x_1).
\end{equation}

If we do not limit ourselves to normal forms, starting from an expression $f_{30}$ (with $x_1$ as one of its variables) for the transition function of CA 30, equivalent expressions for $f_{30}$ are given by

\begin{equation}\label{eq:rule30logic5}
f_{30} \oplus (x_1\odot \lnot x_1),
\end{equation}

\begin{equation}\label{eq:rule30logic6}
f_{30} \odot (x_1\oplus \lnot x_1).
\end{equation}

Now, noting that the ReLU network realization of $f_{30}$ as expressed in (<ref>), is given by

\begin{equation*}
\hat{ \Phi}^{f_{30}}:= \hat{W}_3^{f_{30}} \circ \rho \circ \hat{W}_2^{f_{30}} \circ \rho \circ  \hat{W}_1^{f_{30}},
\end{equation*}

\begin{equation*}
\hat{W_1}^{f_{30}}(x_{-1},x_0,x_1) = \begin{pmatrix}
-1 & -1 & 1 \\
1 & 1 & -1 \\
-1 & 1 & -1 \\
1 & -1 & -1
\end{pmatrix} \begin{pmatrix}
x_{1} \\x_0 \\x_{-1}
\end{pmatrix}  + \begin{pmatrix}
0 \\ -1 \\ 0 \\ 0
\end{pmatrix}\text{,}
\end{equation*}

\begin{equation*}
\hat{ W_2}^{f_{30}}(x) = \begin{pmatrix}
1 & 1 & 1 & 1 \\
1 & 1 & 1 & 1
\end{pmatrix} x + \begin{pmatrix}
0 \\-1
\end{pmatrix}\text{, } \quad x\in \mathbb{R}^4
\end{equation*}

\begin{equation*}
\hat{W_3}^{f_{30}}(x) = \begin{pmatrix}1 & -1 \end{pmatrix} x, \quad x\in \mathbb{R}^2,
\end{equation*}

we can conclude, by comparing to the network (<ref>) built based on (<ref>), that different algebraic expressions for $f_{30}$ lead to different ReLU network realizations. Yet, these two networks must exhibit identical input-output relations.
The augmentation of the Boolean formula in (<ref>), by $x_1 \odot \lnot x_1 = \rho(x_1 + 1-x_1-1)=0$, amounts to composing the network $\Phi^{f_{30}}$ with a null neuron through an extra layer for the operation $\oplus$ according to (<ref>):

\[
\rho(\Phi^{f_{30}} + 0 ) - \rho(\Phi^{f_{30}} + 0-1 );
\]

and in (<ref>), by $x_1 \oplus \lnot x_1 = \rho(x_1 + 1-x_1) -\rho(x_1 + 1-x_1-1) =\rho(1) - \rho(0)=1$, amounts to composing with a constant neuron through an extra layer for the operation $\odot$ according to (<ref>):

\[
\rho(
\Phi^{f_{30}} +1-(2-1)
\]

Conversely, the observation just made shows that ReLU networks can be modified without changing their input-output relations.

Concretely, suppose we are given a function $f:\mathbb{R}^n \rightarrow \mathbb{R}$ and an associated ReLU network $\Phi^f $ realizing $f$. Upon noting that $x = \rho(x) -\rho(-x)$ is a ReLU network and the composition of ReLU networks again yields a ReLU network <cit.>, it follows that composing the network $\Phi^f$ with the identity network
\begin{equation*}
\Phi^{\text{Id}} := \rho(x) - \rho(-x)
\end{equation*}
leads to infinitely many different ReLU network realizations of $f$ as follows:

\begin{equation}\label{eq:augId}
\Phi^{\text{Id}} \circ \cdots \circ \Phi^{\text{Id}} \circ \Phi^f.
\end{equation}

A more sophisticated network transformation that preserves the input-output relation is obtained by starting from a so-called
null-network <cit.> that maps all input values $x$ to $0$, such as, e.g.,

\begin{equation}\label{eq:nullnet}
\rho(x) - \frac{1}{2}\rho(2x) = 0, \quad \text{ for all } x \in \mathbb{R}.
\end{equation}

Rewriting (<ref>) according to $\rho(x) = \frac{1}{2}\rho(2x)$ and applying this relationship recursively yields

\begin{equation}\label{eq:nullnet-replacement}
\rho(x) = 2^{-n} \rho(2^n x), \quad \text{ for all } x \in \mathbb{R},
\end{equation}

for all $n \in \mathbb{N}$.

We can hence generate infinitely many functionally equivalent ReLU networks by replacing individual neurons in a given network according to
(<ref>). The DMV formulae extracted from these networks—following the procedure detailed in Section <ref> below—will be algebraically different in general, yet functionally equivalent.

We have thereby established a way of manipulating DMV formulae

by effecting functional equivalence preserving transformations in the ReLU network domain.

This approach to manipulating formulae in MV logic either in terms of their algebraic expression or in their ReLU network representation can be seen as a generalization of Shannon's results on circuits realizing Boolean functions <cit.>. Boolean formula of a given circuit can be written down and manipulated according to Boolean algebra. In this way, many circuit can be found equivalent in operating characteristics to the given circuit. The synthesis of a circuit to represent a given Boolean formula, in turn, depends upon other considerations from the implementation perspective such as economy of switching elements. We refer to (ICML, 2014) for an in-depth development of this theory.

§.§ Extracting DMV formulae from trained networks

We now discuss how DMV formulae can be read out from the network $\Phi^f$. Recall that the idea is that $\Phi^{f}$ was trained on CA evolution data and the extracted DMV formula describes the logical behavior underlying this CA. In the Boolean case, with $K=\{0,1\}$, the truth table representing the CA transition function can be obtained by passing all possible neighborhood state combinations through the trained network $\Phi^f$ and recording the corresponding outputs.
Following, e.g., the procedure in <cit.>, this truth table can then be turned into a Boolean formula.
For state sets $K$ of cardinality larger than two, we are not aware of systematic procedures that construct DMV terms from truth tables.

The approach we develop here applies to state sets of arbitrary cardinality and does not work off truth tables, but rather starts from

a (trained) ReLU network $\Phi^{f}$ that, by virtue of satisfying the interpolation property (<ref>), realizes a linear interpolation of the truth table.

We start with a result on the fine structure of piecewise linear functions obtained by linearly interpolating CA transition rules.
where Lemma 4.1 was

where Proof of Lemma 4.1 was

The structural property identified in lma:refinefc immediately leads to constraints on ReLU networks realizing interpolated CA transition functions.

We start by noting that $\Phi^{f}$ encodes its underlying DMV formula both through the network topology and the network weights and biases.
We proceed to discuss a procedure for extracting DMV terms from ReLU networks with integer weights and rational biases. We start by noting that a trained ReLU network stores the DMV terms not only in its weights and biases, but also in its architecture. Every piecewise linear function can be realized by a single-hidden-layer ReLU network and, conversely, a deep ReLU network can always be converted to an equivalent single-hidden-layer ReLU network without changing its input-output relation.
Different network architectures lead to different algebraic expressions extracted (see an example in Section 4.3).
Honoring how the network architecture, i.e., layer-by-layer compositions of affine mappings and the ReLU nonlinearity, affects the extracted DMV term, we aim to proceed in a decompositional manner and on a node-by-node basis. It turns out, however, that the individual neurons, in general, will not correspond to DMV terms. To see this, recall that DMV term functions map $[0,1]^n$ to $[0,1]$ and note that the range of the function $\rho$ is not in $[0,1]$ in general, e.g., $\rho(3x): [0,1] \rightarrow [0,3]$. We will deal with this matter by transforming the $\rho$-neurons in $\Phi^{f}$ into $\sigma$-neurons with a suitably chosen $\sigma:\mathbb{R} \rightarrow [0,1]$. This will be done in a manner that preserves the network's input-output relation and is reversible in a sense to be made precise below.

We start by defining $\sigma$-networks.

Let $L \in \mathbb{N}$ and $N_0, N_1, \ldots \hspace{-0.02cm}, N_L \in \mathbb{N}$. A $\sigma$-network $\Psi$ is a map $\Psi: \mathbb{R}^{N_0} \rightarrow \mathbb{R}^{N_L} $ given by

\begin{equation}\label{eq:defSigma}
\Psi= \begin{cases}
W_1, & L=1 \\
W_2 \circ \sigma \circ W_1, & L=2 \\
W_L \circ \sigma \circ W_{L-1} \circ \sigma \circ \cdots \circ \sigma  \circ W_1, & L \geq 3
\end{cases} \text{,}
%    \Psi=
%    W_L \circ \sigma \circ W_{L-1} \circ \cdots \circ \sigma  \circ W_1
%  \text{,}
\end{equation}

where, for $\ell \in \{ 1,2,\ldots \hspace{-0.02cm}, L\}$, $W_\ell: \mathbb{R}^{N_{\ell-1}}\rightarrow \mathbb{R}^{N_\ell}, W_\ell(x):= A_\ell x+b_\ell $ are affine transformations with weight matrices $A_\ell = \mathbb{R}^{N_\ell\times N_{\ell-1}}$ and bias vectors $b_\ell \in \mathbb{R}^{N_\ell}$, and the activation function $\sigma: \mathbb{R}\rightarrow [0,1], \sigma(x):= \min\{1,\max\{0,x\}\}$ acts component-wise. Moreover, we define the depth of the network $\Psi$ as $\mathcal{L}(\Psi):=L$.

We proceed to formally establish the result on $\rho$-network to $\sigma$-network transformation announced above.

For $n \in \mathbb{N}$, let $\Phi:[0,1]^n \rightarrow [0,1]$ be a ReLU network with integer weights and biases in $\mathbb{Q}_k$, with $k \geq 2$. There exists a $\sigma$-network $\Psi: [0,1]^n \rightarrow [0,1]$, with $\mathcal{L}(\Psi) = \mathcal{L}(\Phi)$, integer weights, and biases in $\mathbb{Q}_k$, satisfying
\begin{equation*}
\Psi(x) = \Phi(x), \quad \text{ for all } x\in [0,1]^n.
\end{equation*}

By defn:ReLU, we can write[To keep the exposition simple, we consider the form of $\Phi$ for $L \geq 3$, the cases $L=1,2$ are trivially contained in the discussion.]

\begin{equation*}
\Phi =  W_L \circ \rho \circ W_{L-1} \circ \rho \circ \cdots \circ \rho  \circ W_1,
\end{equation*}
with $N_L = 1, N_0=n, W_\ell(x) = A_\ell x+b_\ell, A_\ell\in \mathbb{Z}^{N_\ell\times N_{\ell-1}}$, and $b_\ell \in \mathbb{Q}_{k}^{N_\ell}$,

all for $\ell\in\{1,\ldots \hspace{-0.02cm}, L\}$. We start with the conversion of neurons in the first hidden layer of $\Phi$. These neurons are all of the form

\begin{equation*}
\rho(a_1x+b_1), \quad \text{ with } a_1 \in \mathbb{Z}^{1\times n}, b_1\in \mathbb{Q}_{k}.
\end{equation*}
As the domain of $\Phi$ is the unit cube $[0,1]^n$, $a_1x+b_1$ resides in a bounded interval in $\mathbb{R}$. Denote this interval by $[-u, u]\subset \mathbb{R}, u\geq 0$. If $u\leq 1,$ no conversion is needed as $ \rho(x')  = \sigma(x'),$ for $x'\in [-u,u]$. For $u>1$, we note that

\begin{equation}\label{eq:relu2sigma}
\rho(x)=  \sigma(x) + \sigma(x-1)+\cdots +\sigma(x-m), \quad x\in [-u,u],
\end{equation}

with $m\in \mathbb{N}, m\geq u-1$. In summary, we can replace each $\rho$-neuron in the first hidden layer by one or several $\sigma$-neurons according to

\begin{equation*}
\rho(a_1x+b_1) = \sigma(a_1x+b_1)
\end{equation*}
whenever $a_1x+b_1 \leq 1$, for all $x\in [0,1]^n$, or

\begin{equation*}
\rho(a_1x+b_1) =\sigma(a_1x+b_1) + \cdots +\sigma(a_1x+b_1-m)
\end{equation*}
when $a_1x+b_1 \leq m+1,$ for all $x\in [0,1]^n$, for some $m>1, m\in \mathbb{N}$. The resulting network $\Phi^{(1)}:[0,1]^n\rightarrow[0,1]$ is given by

\begin{equation*}
\Phi^{(1)} =  W_L \circ \rho \circ \cdots \circ W_3 \circ \rho \circ W_2^{(1)} \circ \sigma  \circ W_1^{(1)},
\end{equation*}
with $W_1^{(1)}(x) = A_1^{(1)}x+b_1^{(1)}, A_1^{(1)}\in \mathbb{Z}^{N_1^{(1)}\times n}, b_1^{(1)} \in \mathbb{Q}_{k}^{N_1^{(1)}},$  and $W_2^{(1)}(x) = A_2^{(1)}x+b_2, A_2^{(1)}\in \mathbb{Z}^{N_2\times N_1^{(1)}}$, where $ N_1^{(1)}$ is the number of $\sigma$-neurons in the first hidden layer of $ \Phi^{(1)}$, and satisfies

\begin{equation*}
\Phi^{(1)}(x) = \Phi(x), \text{ for } x\in [0,1]^n.
\end{equation*}