Computational capabilities at the edge of chaos for one dimensional systems undergoing continuous transitions

The transition from regular towards chaotic behavior is one of the defining characteristics of complex non-linear dynamical systems. There has been a keen interest in the behavior of a dynamical system as it undergoes such transition, which comes from the fact that it is precisely in such region where enhanced or improved computational capabilities of the system are claimed to be found crutchfield90 ; langton90 ; kauffman92 , this is known as the edge of chaos (EOC) hypothesis. In this context, computation in a dynamical system can have different meanings crutchfield93 . The most common interpretation of computational capacity is related to the ability of the system to perform some “meaningful” task. Meaningful can be understood, given some initial condition, as the system transforming the input in some prescribed way into some output data. More broadly, it can be understood related to the system capability to behave as a universal Turing machine given some proper input and “programming”. A related but different point of view is to associate computation to the generic production, storage, transmission and manipulation of information. It is in this last sense that any dynamical system can be seen as some device amenable to be characterized by its intrinsic computational abilities crutchfield12 . Measures of the computational capabilities of a system have been developed that aim at quantifying the efficiency, in terms of resources, while characterizing the balance between irreducible randomness and pattern production crutchfield12 .

In favor of the EOC hypothesis several systems modeling natural and artificial phenomena exhibit interesting, if not some type of critical behavior, precisely at the edge of chaos su89 ; sole96 ; bertschinger04 ; beggs08 ; boedecker12 . It is not clear, for some one dimensional systems, if the EOC criticality plays a role in the enhancement of computational capabilities, for example in cellular automata. Cellular automata (CA) are discrete space, dynamical systems that act over a spatially extended grid taking values over a finite alphabet, resulting in the evolution of the spatial configuration in discrete time steps. There has been a discussion if meaningful computation abilities in CA are found and should be sought at EOC regions.

Langton langton90 has used the density of quiescent states in the cellular automaton rule, usually denoted by the letter $\lambda$, as the control parameter to explore such questions. The quiescent state is chosen as some arbitrary value in the CA alphabet, and the $\lambda$ parameter is the ratio between the number of times the quiescent state does not appear in the CA rule table and the number of entries in the same table. For a binary alphabet CA with nearest neighbor rules (elementary cellular automata or ECA), $\lambda=1-n\times 2^{-3}=1-n\times 8$, where $n$ is the number of times the quiescent state (say value $1$) appears in the rule table. Monte Carlo simulations were performed over the CA space to conclude that, on average, as $\lambda$, increases from $0$, the rules behavior changes from a fixed point regime, to periodic, to complex and finally to chaotic, each regime corresponding to one of Wolfram classification scheme wolfram02 . Packard packard88 further developed these ideas to conclude that CA capable of meaningful computations happens precisely at those critical values of $\lambda$ where a transition to the chaotic behavior occurs. Packard results were questioned by Crutchfield et al. crutchfield93 who performed Monte-Carlo simulations similar to the ones made by Packard failing to find EOC enhanced computation capabilities.

In this paper, we show that improvement of computational capabilities indeed occurs at the edge of chaos for some continuously extended CA rules, but not following Langton $\lambda$ parameter and therefore, not in the sense of Packard simulations. We will not be seeking to discover useful computational of the CA but instead, improvement of intrinsic computation as measured by entropic based measures of structure and disorder. To support this claim we follow the behavior of a CA as one entry in its rule table is continuously varied between $0$ and $1$ and use it as control parameter the value of the continuously varying entry. After discussing CA, we turn our attention to a system of one dimensional non-linear locally coupled oscillators that have been found to exhibit dynamical behavior similar to those of CA alonso17 . The quest is to find if in this system with a continuous control variable space, enhanced computation at the EOC can also be found.

To start formalizing some of the expressed ideas, for our purposes it suffices to consider one dimensional CA as a tuple $\{\mathcal{A},r,\phi\}$ consisting of a set of states $\mathcal{A}$, a radius $r$ which defines a neighborhood structure from $-r$ to $r$, and $\phi$ a local rule $\phi:\mathcal{A}^{2r+1}\rightarrow\mathcal{A}$. The local rule acts simultaneously over a numerable ordered set $S$ of cells which can be bi-infinite, but in any case, each cell is properly labeled by natural numbers as in $\{\ldots,s_{-1},s_{0},s_{1},\ldots\}$. At a given time step $t$ each cell is in a state belonging to $\mathcal{A}$ and the whole cell array configuration is denoted by $S^{(t)}$. The evolution of the cell array from time $t$ to time $t+1$ is given by the global function induced by $\phi$: $\forall S\in\mathcal{A}^{\mathbb{Z}},\;S^{(t+1)}=\Phi(S^{t})$ such that $s_{i}^{(t+1)}=\phi(s_{i-r}^{(t)},\ldots,s_{i}^{(t)},\ldots,s_{i+r}^{(t)})$. $S^{(t)}$ is also called the spatial configuration of the cell array at time $t$. The set of $\{S^{(0)},S^{(1)},\ldots,S^{(t)}\}$ is called the spatiotemporal diagram or evolution up to time $t$. The spatiotemporal diagram of the CA can show a wide range of behaviors including the emergence of complex patterns, fractal nature, chaotic behavior, periodic, random evolution, or evolve to homogeneous spatial configurations. CA have been used as models for different natural phenomena as well as a model for computing devices. It is known at least of a one dimensional CA that behaves as a Universal Turing Machine, the so-called rule 110 cook04 ; wolfram02 .

Consider two state $\mathcal{A}=\{0,1\}$ cells with a rule defined over a one neighborhood radius $r=1$ such that, a cell $s_{i}$ will be updated in the next time step according to the current value of the cell and its two nearest neighbors: $s_{i}^{(t+1)}=\phi(s_{i-1}^{(t)},s_{i}^{(t)},s_{i+1}^{(t)})$. Such cellular automata are named elementary cellular automata or ECA wolfram02 . To specify a rule function $\phi$ a look up table can be given, for example as shown in Table 1

When ordered lexicographically the rule has the binary representation $00101110$ equivalent to the decimal number $46$ (the reader should notice that according to the definition given in the introducion, the $\lambda$ parameter for this rule is $0.5$). CA rules are named by their decimal representation. The space of cellular automata with binary alphabet and nearest neighbor interaction will then be a discrete space of eight dimensions or an eight dimensional hypercube where there are $2^{2^{3}}=2^{8}=256$ rules.

There is no apparent geometry in the CA discrete space, as there is no natural way to define a transition between CA rules in the hypercube allowing for the smooth change of one cellular automaton to a neighboring one while their behavior also changes gradually. One may think in ordering the CA rules as Gray codes savage97 , where successive CA differ by a change of only one entry in the rule table. However, Gray codes ordering does not lead to a corresponding smooth transition in the CA regimes. There is not even any apparent relation in behavior between adjacent rules. Consider rules $238$ and $110$ with binary representation $00101110$ and $01101110$, respectively. Both rules differ only in the seventh entry and yet, the former result always in a homogeneous state where each cell in the array has the same value, while the later, as already stated, can perform as a Universal Turing Machine.

The inability to define a natural ordering among the CA rules along which behavior changes gradually, has hindered the analysis of transitions between CA rules. To overcome this limitation, Pedersen devised a procedure to continuously deform a CA rule into another at the expense of sacrificing the discrete nature of the cells states petersen90 . Consider rule $46$ and rule $110$ , a continuous transition from the former to the later can be done if we define a CA rule as $0\xi 101110$, where $\xi$ goes from $0$ to $1$ in a continuous manner (Figure 1), such that for $\xi=0$ we get rule $46$ and for $\xi=1$ rule $110$ is recovered. In order for the transition to work, first a real valued function $\beta$ must be defined, for example $\beta(s_{i-1}s_{i}s_{i+1})=4s_{i-1}+2s_{i}+s_{i+1}$, next, an interpolating function $f(x)$ is used over the values of $\beta$, in our case: $f(\beta(s_{i-1}s_{i}s_{i+1})):[0,7]\rightarrow[0,1]$. The interpolating function is a real valued continuous function that must be consistent with the discrete case for integer values of its argument, but apart from that, it can be chosen in any way. Pedersen used three interpolating functions, one lineal, one quadratic and a third by trigonometric splines, finding that the results were qualitatively similar for every choice. We follow Pedersen and settle for the quadratic interpolating function

where $n(\equiv\lfloor x\rfloor)$ is the largest integer smaller than $x$. As explained, if $n$ is an integer number $f(n)$ is given by the discrete CA rule.

The price paid in Pedersen extension is that the states of the cells are no longer binary values $\{0,1\}$ but any real value in the interval $[0,1]$. In order to recover the discrete binary nature of the cell states, we observe the system via a binarizing measuring device: we use the mean value over the whole cell array for a given time step as a threshold value, any cell with a state value over the mean is projected to $1$, otherwise to $0$.

In our simulations, an array of $5\times 10^{3}$ cells were used and the initial random configuration (random binary digits were taken from www.random.org ) was left to evolve in $5\times 10^{3}$ time steps. Simulation results are then presented by averaging over ten simulations for each $\xi$ parameter value, each with different initial configurations. The last ten spatial configurations were taken for each run, so the final averaging was over $10^{2}$ spatial configurations. For testing the robustness of the result, calculations were occasionally performed for $10^{4}$, $5\times 10^{4}$, $10^{5}$ and $10^{6}$ cells, and correspondingly equal time steps, showing no change in the results. Cyclic boundary conditions are used. In all cases, and contrary to Packard Monte Carlo simulations, the variance of the mean values for our simulations is too small even to be represented as error bars.

Before we proceed, we need some measures that can capture the computational capabilities of a system and distinguish from the random portion of the output on the one hand, and data structuring on the other.

As a measure of unpredictability we use Shannon entropy rate $h_{\mu}$. When a bi-infinite sequence $S=\ldots s_{-3}s_{-2}s_{-1}s_{0}s_{1}s_{2}s_{3}\ldots$ of emitted symbols $s_{i}$ has been observed, $h_{\mu}$ is the amount of information produced per symbol. If we consider the bi-infinite sequence as a time series, then, it is the amount of new information in an observation of let say cell $s_{i}$, considering the state of all previous cells $s_{j}$, $j<i$ crutchfield03 . In terms of Shannon entropy

where $H(X|Y)$ denotes the Shannon conditional entropy of random variable $X$ given variable $Y$cover06 .

The block entropy of a bi-infinite string $S$ measures the average uncertainty of finding blocks of a given length within the string. Block entropy is calculated by the equation

where the sum goes over all sequences $S^{L}$ of length $L$ drawn from the alphabet $\mathcal{A}$, and $p(S^{L})$ is the probability of observing one particular string $S^{L}$ in the bi-infinite string $S$. Entropy density, also known as entropy rate, can be understood as the unpredictability that remains when all that can be seen has been observed, and can also be defined through a limiting process of the normalized block entropies

For a necessarily finite data, the entropy rate has to be estimated. There is a body of literature dealing with the estimation of entropy rate for finite sequences schurmann99 ; rapp01 ; lesne09 . Here $h_{\mu}$ will be estimated through Lempel-Ziv factorization procedure that has been extensively discussed in the literature and has been used before in the context of CA kaspar87 ; estevez15 . In short, Lempel-Ziv factorization performs a sequential partition of a finite length string adding a new factor if and only if, it has not been found in the prefix up to the current position of analysis lz76 . In order to emphasize that a finite size estimate of $h_{\mu}$ through the Lempel-Ziv procedure is being used, we change to the notation $h_{LZ}$ instead.

Formally, if $s(i,j)$ is the substring of consecutive symbols $s_{i}s_{i+1}\ldots s_{j}$ and $\pi$ the ”drop” operator defined as $s(i,j)\pi=s(i,j-1)$ and consequently, $\pi^{l}=s(i,j-l)$, where $l$ is a positive number. The Lempel-Ziv factorization is a partition $F(S^{N})$ of a finite string $S^{N}$ of length $N$

in $m$ factors such that each factor $s(l_{k-1}+1,l_{k})$ complies with

1. 1.

   $s(l_{k-1}+1,l_{k})\pi\subset s(1,l_{k})\pi^{2}$

2. 2.

   $s(l_{k-1}+1,l_{k})\not\subset s(1,l_{k})\pi$ except, perhaps, for the last factor $s(l_{m-1}+1,N)$.

The partition $F(S^{N})$ is unique for every string. For example the sequence $u=101101110001$ has a Lempel-Ziv factorization $F(s)=1.0.11.0111.00.01$, each factor is delimited by a dot. For a data stream of length $1\ll N<\infty$, $h_{LZ}$ is defined by

where $C_{LZ}(S)=|F(S^{N})|$ is the number of factors in the Lempel-Ziv factorization. The entropy rate for an ergodic source is then given byziv78

For a finite data stream of length $N$ where the above limit can not be strictly realized, the entropy rate is then estimated by $h_{LZ}$ where we drop the argument $N$ when no ambiguity arises. The error in the estimation for a given value of $N$ has been discussed beforeamigo06 .

Excess entropy $E$, introduced by Grassberger as effective complexity measuregrassberger86 , has been used in several contexts crutchfield12 including the analysis of CA as a measure of correlation along different scales in a process and related to the intrinsic memory of a system. For a bi-infinite string, the excess entropy measures the mutual information between two infinite halves of the sequence and is related to the correlation of the symbols at all possible scales. Excess entropy can be defined in several equivalent ways, for example as the convergence of the entropy rategrassberger86 ; crutchfield03

$I[X:Y]=H[X]+H[Y]-H[X,Y]$ is the mutual information between $X$ and $Y$, and $h_{L}(S)=H_{S}(L)-H_{S}(L-1)$, where $H(X)$ is the Shannon entropy of the random variable $X$. As already explained, mutual information is a measure of the amount of information one variable carries regarding the other, and it is a symmetric measure $I[X:Y]=I[Y:X]$. The information gain $\Delta H[L]=h_{L}(S)-h_{\mu}(S)$ measure how much the entropy density is being overestimated as a result of making measurements of the probability of only string blocks up to length $L$, in other words, how much information still has to be learned in order to assert the true entropy of the sourcecrutchfield03 . As entropy decreases with length $L$, if such decrease is fast enough, the sum (7) will converge. $E(S)$ is the intrinsic redundancy of the source or apparent memory of the sourcecrutchfield03 . It measures structure in the sense of redundancy being a measure of pattern production. Feldmancrutchfield03 has proven that equation (7) is equivalent, in the ergodic case, to the mutual information between two infinite halves of a bi-infinite string coming from the source output. Structure in this context then means how much information one halve carries about the other and vice-versa, a measure related again with pattern production and context preservation as redundancy.

Numerically, when dealing with finite data streams, excess entropy has to be effectively estimated. Here we use a random shuffle based excess entropy estimate using the Lempel-Ziv estimation of entropy rate and given by melchert15

$S_{(M)}$ is a surrogate string obtained by partitioning the string $S$ in non-overlapping blocks of length $M$ and performing a random shuffling of the blocks. The shuffling for a given block length $M$ destroys all correlations between symbols for lengths larger than $M$ while keeping the same symbol frequency. $M_{max}$ is chosen appropriately given the sequence length as to avoid fluctuations. In spite that $E_{LZ}$ is not strictly equivalent to the excess entropy as given by equation (7), it is expected to behave in a similar manner melchert15 .

Finally, information distance $d(s,p)$ will be used. Information distance is derived from Kolmogorov randomness theory. The Kolmogorov randomness $s^{*}=K(s)$ of a sequence $s$, is the length of the shortest program ($s^{*}=|K(s)|$), that can reproduce the sequence $s$kolmogorov65 , accordingly, $K(s|p^{*})$, known as Kolmogorov conditional randomness, is the length of the shortest program able to reproduce $s$ if the program $p^{*}$, reproducing the string $p$ is given. The information distance is defined byli04

$d(s,p)$ is an information based distance measuring, from an algorithmic perspective, how correlated are two sequences: if two sequences can be derived one from the other by a small-sized algorithm, then the corresponding $d(s,p)$ is small.

Following Estevez et. al.estevez15 , we estimate $d(s,p)$ via Lempel-Ziv by

which will have the same interpretation than $d(s,p)$ as much as the normalized (5) estimates the entropy density.

We finally look into the complexity-entropy map crutchfield12 of both systems. The entropy map for cellular automaton rule $78$ is shown in figure 5a. Rule $46$ has a similar diagram. The map was calculated by continuously varying the sixth entry in the rule table. Feldman et al. feldman08 discussed the interpretation of complexity-entropy diagrams; it measures the system intrinsic computation capabilities in a parameterless way. As already discussed, we also found enhanced computation at the EOC for rule $78$ (See supplementary material). The complexity-entropy map of figure 5a is typical of systems transforming from non-chaotic to a chaotic regime, the reader may compare our plot with that reported by Crutchfield and Young for the logistic map exhibiting the well-known period-doubling cascading crutchfield89 . For values of entropy rate near zero, the system is rather simple, and the $E_{LZ}$ values are small. It is only with the introduction of some unpredictability that the system can accommodate more complex behavior. It is precisely at those value around $h_{LZ}=0.2$ where the $E_{LZ}$ attains a maximum that corresponds to the EOC. However, it is only to a certain value that a system can accommodate randomness while attaining higher intrinsic computation capability.

Further increase of entropy rate results in the decrease of the system intrinsic computation capability and at $h_{LZ}=1$, $E_{LZ}$ is again at zero. A similar interpretation can be made of the complexity-entropy diagram for the NLCO the main difference is that the curve is smoother pointing to a continuous first derivative. The striking similarity between our plots and those reported by Crutchfield and Young for well-studied systems, further points to the fact that the studied systems undergo a transition to chaotic behavior as the parameter $\xi$ is varied.

In conclusion, we have argued for the existence of enhanced computation at the edge of chaos for continuously deformed cellular automata, different from the reports of Langton and Packard. By using a continuously varying parameter that allows transforming from one CA rule to another, geometry in the CA space was recovered. Our approach does not carry the limitations pointed out against Langton studies. The evidence heavily points to the existence, for specific CA rules, of a transition towards a chaotic regime where enhanced intrinsic computation can be found at the EOC. We directly measured such enhancement by looking into a measure of structural complexity while following the behavior of the entropy rate of the system. $k$-balanced analysis seems to be a useful tool to identify other rules where the same phenomena could be expected. The analysis was extended to a system of nonlinear coupled oscillators whose behavior has a strong resemblance to CA spatiotemporal patterns. For this system also improved intrinsic computation at the EOC was identified. The complexity-entropy diagrams show a striking resemblance to the ones reported for well studied non-linear systems were chaotic transitions are known to happen.

This work was partially financed by FAPEMIG under the project BPV-00047-13. EER which to thank PVE/CAPES for financial support under the grant 1149-14-8. Infrastructure support was given under project FAPEMIG APQ-02256-12.