Chaos

1. Aperiodicity:

We try to visualize periodicity in our 4 functions in the Fourier domain. Below shows the Discrete Fast Fourier Transform (DFFT) of the first 1000 samples of each function. Sine function as expected has a sharp peak at its frequency = 0.04 i.e its periodic with T = 25 which can be observed in Fig 1 above. Meanwhile one can notice multiple frequencies in other functions. For Eq 2, it looks as if some frequencies are affecting the function more while others negligibly. For Eq 3 and Eq 4, there is somewhat even distribution of frequencies with large fft coefficients. Even though the periodicity is clear for the sine function, from Fourier transforms only, we are unable to ascertain the periodicity of other functions.

• 

• 

• 

• 

Below is the autocorrelation plots for the 4 functions. Autocorrelation or self-convolution is done in the Laplace domain (here for the amazing Fourier and Laplace connection). From autocorrelation plots we can clearly see that Eq 2 is self similar (autocorrelation approximately 1) for t somewhere near 70. In fact, we find that f(t) = f(t + 76) for Eq 2. Hence, Eq 2 is not chaotic and has periodicity T =76. For Eq 3 and Eq 4, no such autocorrelation is observed (10000 samples is the largest Python allows for acf).

• 

• 

• 

• 

2. Sensitivity to Initial Conditions:

We are at the heart of the matter. Newtonian determinism implied that perfectly calibrated inputs would give rise to perfectly certain outputs. Such determinism led to a model of clockwork universe. However, Poincaré was earliest to point out a flaw in the model when he added the Moon in the Earth-Sun dynamics where the Newtons’ Law of Gravitation failed, the problem of scientific ménage à trois.

Sensitivity to initial conditions is crucial to chaos and is often termed as the butterfly effect. It essentially means that no matter how close two points start out with, their pathways eventually diverge. Therefore, one cannot naively treat tiny differences in inputs as noise, which stabs in the heart of reductive approaches such as linear and non-linear regressions. Below shows the plots showing sensitivity to initial conditions for Eq 3. Plot 1 & 3 show the diverging paths due to slight differences in starting points. Plot 2 & 4 show the absolute differences in the paths in which one can observe tiny differences initially that become large after a certain time.

• 

• 

• 

• 

Such sensitivity to initial conditions is measured in terms of Lyapunov exponents. The idea behind Lyapunov exponents is similar to the idea of a derivative, we measure the change in dependent variable with respect to small change in the independent variable. Here the independent variable is the initial condition and dependent variable is the state of the function at a later time. Strictly, for any two initial conditions (x_{00}, x_{10}), let the state of the system at time t be (x_{0t}, x_{1t}) = (x_{00}+dx_{0t}, x_{10}+dx_{1t}), we define \delta_t = |x_{0t}-x_{1t}| as the divergence between the two trajectories at time t. The Lapunov exponent is then defined⁴ long term logarithmic rate of divergence of two trajectories: \sigma = \lim_{t\to\infty}\frac{1}{t}ln(\frac{\delta_t}{\delta_0}) = \lim_{t\to\infty}\frac{1}{t}\sum_{i=0}^{t-1}ln|f'(x_i)| \hspace{10mm} \text{(5)} Note that there will be n Lyapunov exponents for n-dimensional maps. Its also enlightening to rewrite Eq. 5 as \delta_t^{\varepsilon} = \delta_0e^{\sigma t}\hspace{10mm} \text{(6)} where \delta_t^{\varepsilon} = |x_{0t}-x_{1t}|\le \varepsilon for some error tolerance \varepsilon \in \mathbf{R^+} from which one can compute the Lyapunov time t_{\varepsilon} such that \delta_t \le {\varepsilon}. Using Eq 5, we compute the Lyapunov exponent for Eq 3 to get \sigma \approx 0.693 . This value plugged in to Eq 6 provides t_{\varepsilon} \approx 7 and t_{\varepsilon} \approx 37 for \varepsilon = 0.1 \text{ and } \delta_0 = \text{1e-3}, \delta_0 = \text{1e-12} , respectively, which aligns with the observation in Fig 5.

3. Deterministic Inputs Indeterminate Outputs:

To me, this is the most exciting concept in the ideas of chaos. We noted earlier that aperiodicity implied lack of rigid temporal structures. From Fig 1 one can note the y-axis values for Eq 3 seem to be random. However, Chaos and Random aren’t the same but they’re close. For example, below are the 20 iterations of the Eq 3:

[0.2,0.64,0.9216,0.289,0.8219,0.5854,0.9708,0.1133,0.402,0.9616,0.1478,0.5039,0.9999,0.0002,0.001,0.0039,0.0157,0.0617,0.2317,0.7121]

For each value above, if its greater than 0.5 we set it to ‘Heads’ else ‘Tails’:

['T', 'H', 'H', 'T', 'H', 'H', 'H', 'T', 'T', 'H', 'T', 'H', 'H', 'T', 'T', 'T', 'T', 'T', 'T', 'H']

What we just did above is Symbolic coding of a chaotic function. We then claim that the series resulting from such coding of our Eq 3 indistinguishable from random⁵. This is an absurd claim. However for all practical purpose its true. So true that its the basis of random number generation in a deterministic system such as a computer, and are called Psuedo-Random Numbers. I suggest you try it out, compare the values with a real random coin toss sequence from Random.org (here of an amazing blog on tests of randomness).

Now you’d have to say, wait a minute, each sample in a coin toss is independent of the previous or next or any other tosses. This is clearly not true for Eq 3 since x_{n+1} depends on x. We can see this in (x_{n},x_{n+1}) Poincaré plots below which reveals persistent structures for our four equations. We note the inverted parabola for Eq 3 while a dense discrete network for Eq 4. This should give us the estimate of the complexity of the functions we are tying to learn. Eq 1 and Eq 3 seem easier compared to Eq 1 and especially Eq 4 for function learning due to discontinuities involved in the latter. Finally, since there exist functions that generate pseudo-random numbers with persistent structures before being masked by symbolic codings, could we learn this hidden structure from the series of psuedo-random numbers and win bets on average with computers. We try this out in our experiments.

• 

• 

• 

• 

Learnability

Some notes on learnability before we move to experiments. We know that learning is more than memorization. A perfect memorization fails to generalize to unseen examples and is called overfitting. Learning is more of an encoding. Perfect encoding such as formulae, eg. Newton’s laws, are elegant and generalize very well to unseen data. The idea of reduction is to follow occam’s razor i.e. the simplest formula that describes the phenomenon well enough is chosen. Most of the time, such simple formula describes trends and not individual data points and especially not anomalies, which are usually treated as noise. In our case, we will treat our signal as void of noise. We want to be able to encode such a way that we describe individual data points well. Hence, we will try to fit individual training points while maintaining testing accuracy and avoiding the overfitting. We are not interested in seasons, but the weather, esp given the weather now and in the past, if its going to rain in an hour, figuratively.

Based on Information Theory, one can also think of encoding as the number of bits it takes to describe something⁶. Lets take a number 2, it takes two bit to describe it perfectly, while a number 1/3 takes around 5 bits (1 for 1, 2 for 3 and 2 for ‘/’ if its chosen from a set of (+,-,* , / ), etc. Other representations of 1/3 are the string ‘1/3’, '0.\bar{3}', etc. Hence, we can ignore how strings or ints are actually represented in a computer but think of it in terms of information representation⁶. All rational numbers have finite representation. Almost all irrational numbers have no finite representation and thus no encoding. However, there are some which do, for example, \sqrt{2}. It took me two letters to write it down and communicate perfectly what it is. The same \sqrt{2} written as ‘1.4142135623..’ is not a finite representation. Hence, while the data may seem to be infinite in one representation, it may still have a finite encoding in another representation. So the question is if, from the time series data generated by an irrational number \sqrt{2}, we could learn the function itself, i.e. sqrt().

Btw, any solution to analytic function, if it exists, has finite representation, the representation being the function itself. I’ve chosen \pi for Eq 4. \pi‘s finite representations include the symbol \pi itself, the formula \frac{Perimeter}{2r}, its power series expansions, continued fractions, and others. The idea is that each digit of \pi is not random for the generator function exists, BBP for example, and we want to learn it using function approximators such as neural nets. Finally, you maybe wondering why we are using an irrational number in the analysis of chaotic functions. The symbolic coding such as one we did earlier generates a unique irrational number for a given initial condition, for example 011011100… is binary for some irrational. Hence there is a deep connection between the two and symbolic dynamics are usually enough to understand the dynamics of chaotic systems⁵.

Experiments

1. Discrete Fast Fourier Transform (DFFT):

From Fig 3 and Fig 4, we were able to ascertain periodicity of Eq 1 and Eq 2. With finite periods, DFFT is the go-to method for signal reconstruction. Here we initially map our time series in Frequency Domain, filter high order frequencies (noise) if needed and remap the signal back to the time domain. Below shows the results of DFFT for the 4 equations. With minimal training samples of 1 period, we were able to reconstruct the signals quite perfectly for Eq 1 and 2. However, one can notice this method failing for Eq 3 and 4. Due to lack of periodicity, we are forced to choose finite frequencies for reconstruction. This implies that the DFFT method will reconstruct only a subset of the signal which will repeat. Hence, the method fails to generalize right away to unseen data (frequencies).

• 

• 

• 

• 

2. Kalman Filters:

A discrete-time Kalman Filter (KF) is equivalent to a Bayesian belief propagation technique⁷ and is widely used in engineering. A Bayesian method of inference uses Bayes’ rule to update probability for a hypothesis or estimates as more information becomes available. A belief underlies uncertainties to the hypothesis. Hence, KF is two step algorithm. In the prediction step, estimates of next state of a variable with uncertainties are computed. The state is then measured, with noise, and the prediction estimates are updated using weighted averages. The underlying model is HMM (here for an amazing lecture on Bayesian networks).

• 

• 

• 

• 

Above shows the predictions of the KF for Eq 3 & 4. We’ve also plotted the performance of Extended Kalman filters (EKF). KF is linear in state space and doesn’t capture non-linear dynamics well. One way to get around this is to linearize the dynamics around the weighted mean which leads to EKF. However in our case, we measure at each iteration and the dynamics is linear between two measurements, hence EKF & KF performs similar. Due to lack of analytical understanding for Eq 4, I’ve ad hoc used derivative of the sigmoid as the non-linear Jacobian.

3. AutoRegressive Integrated Moving Average (ARIMA):

KF worked well for single step predictions. However, we would like to peek further into the future. ARIMA uses KF along with other tools: auto-regression (AR), differencing (I) , and moving averages (MA) to model time series in a more robust manner. Auto-regression builds linear regression model by lagging the time-series with itself and is defined by parameter ‘p’. Differencing is a technique where you subtract the timeseries by its previous iteration. The parameter ‘d’ defines the number of time you difference the series to get a stationary model. Moving averaging is where you update the current model with new measurements with decay parameters. Parameter ‘q’ sets the order if the MA. Below is the results of ARIMA on Eq 3 and 4 with params (p,d,q). Note that the SCIPY.stats implementation only allows d < 2 and thus we self create stationary version by differencing ‘l’ times on plots 3 & 4. We note that in both naive and stationary models, the prediction is good to at most 2 steps in the future.

• 

• 

• 

• 

5. Long Short Term Memory (LSTM):

LSTM are a class of Recurrent Neural Networks (RNNs), which are the neural nets with feedback loops. The feedback loop from the RNN helps take the output of previous iteration as a input feature for the current iteration helping the neural network learn time series data such as language, video, stock prices etc. While in principle the network should learn from past experience, RNN fails to learn long term context due to vanishing gradients. The LSTM improves upon RNN’s by utilizing CEC’s that help counteract the vanishing gradients making long term contextual learning feasible⁸. We will be utilizing the network for chaos prediction.

5.1. The Logistic map

5.1.1. Forecasting & Generalization:

The LSTM regressor architecture⁹ comprised of 2 layers of Bidirectional LSTM with 4 & 2 units respectively followed by two dense hidden layers of size 8 and 4 units. The activation functions were all Relus and the GD optimized with adam. The cost function used was RMS loss. No dropout were used since the goal was to minimize the training error as much as possible. The dataset used was first 10K samples of the iterations of Eq 3 with x0 = 0.2 with 67% of the data used in training and rest used in testing. A 1 step look-back gave the best results. This makes sense since Eq 3 only depends on its previous recurrence.

Below shows the results of LSTM on Eq 3. I was able to achieve a training error of 1e-3 on around 1000 epochs. We can see from Fig 10 below that the LSTM fits the function fairly well for the training data. For the test data, we make multi step prediction by taking the last training sample and making the prediction, then taking that prediction and making another one and so on. We see that LSTM does great in multi time step forecasting beyond the training data for about 7 iterations into the future.

• 

• 

Fig 11 a) shows the maximum iterations before the test error \gt 0.1 for each test sample with average ~ 7.68 iterations. This is consistent with the Lyapunov time we computed for 1e-3 difference in the initial condition signalling the possibility that had we been able to reach a lower training error, the prediction window would increase as well. Fig 11 b) shows a correction scheme where after 6 iterations we correct for the residual error in our predicted path similar to the KF update, which allows the LSTM to predict infinitely into the future.

• 

• 

Fig 12 below shows generalization of the learned model to Eq 3 with different initial conditions. Similar prediction window of around 8 iterations is observed before the chaos prevails.

• 

• 

• 

• 

5.1.2. Forecasting with Noise:

All real life dynamics include noise (sensor noise, precision, etc). We would like to see the performance of LSTM with noisy signal. Here we have added small noise n to Eq 3 such that n is gaussian N(0,0.01) if 0 \le | Eq 3 + n| \le 1 else n = 0. We used the same architecture in 5.1.1 above except that small dropouts were now added to our hidden layers (5% and 1% respectively). Fig 13 shows the prediction performance of LSTM. The minimum training error achieved was around 0.034. Fig 14 a) shows the reduction in the prediction horizon to an average of 4 steps with Fig 14 b) showing consistent reconstruction still possible with correction updates after each t=3 prediction.

• 

• 

• 

• 

5.1.3. Predicting Symbolic Dynamics:

This experiment converted the dataset to binary by thresholding at 0.5 similar to when we discussed the symbolic coding. LSTM classifier was trained with architecture similar to 5.1.1. The max training accuracy achieved was ~0.52 i.e. 52%. While this is great and slightly better than a random guess, the test accuracy was 0.5 and wasn’t significantly different from random for any prediction intervals beyond training. Hence, the PRNG scheme still looks safe with this regard.

5.2. Digits of Pi

The LSTM regressor architecture comprised of 2 layers of Bidirectional LSTM with 16 & 4 units respectively followed by three dense hidden layers of size 8, 4 and 4 units. The activation functions were all Relus and the GD optimized with adam. The cost function used was RMS loss. No dropout were used. The dataset used was first 1500 samples of digits of Pi scaled to 0-1. 67% of the data was used in training. While the look-back in range 64 to 128 were tested with larger LSTM and hidden units, they failed to converge. A 6 step look-back trained on 4K epochs gave the training error of 0.09 but fails to generalize. Below we note the overfitting on the training set which massive overshoots in the test data. This is most likely Runge’s phenomenon as the LSTM tries to fit the training data will higher order polynomials as it fails to find a low order representation.

• 

• 

As with 5.1.3, I’ve tried to run a LSTM classifier with 10 output nodes predicting the next digit of Pi from the previous values. I was able to get training accuracy of ~0.24 with the testing accuracy was no better than 0.1 for any prediction window and hence, the classifier was no better than random. I’d also noticed the erratic nature of the loss possibly due to the non-continuous nature of the cost function, which did make me appreciate classification as much of a harder problem than regression.