Permutations uniquely identify states and unknown external forces in non-autonomous dynamical systems

In addition to the theoretical proofs provided above, we tested our hypothesis numerically. For testing our idea, we used the logistic map³⁴34. R. May, Nature 261, 459 (1976). https://doi.org/10.1038/261459a0 and the Hénon map³⁵35. M. Hénon, Commun. Math. Phys. 50, 69 (1976). https://doi.org/10.1007/BF01608556 subject to dynamical noise. The logistic map we used is defined as follows:

$$x_{t+1}=(3.7+{ϵ}_t)x_t(1-x_t),$$

(6)

where ${ϵ}_t$ is a source of independent uniform noise on $[-0.1,0.1]$. We chose the initial condition $x_0$ from a uniform distribution on $[0,1]$ and observed the time evolution of the variable $x_t$. Similarly, the Hénon map we used is defined as follows:

$$\begin{array}{r}\begin{array}{rl}{x}_{t+1}& =1-(1.2+{\eta }_t){x_t}^2+0.3y_t,\\ y_{t+1}& =x_t,\end{array}\end{array}$$

(7)

where ${\eta }_t$ independently follows a uniform distribution on $[-0.05,0.05]$. In addition, we choose the initial conditions $x_0$ and $y_0$ from the uniform distribution on $[0,1]$. Therefore, here we are considering dynamical systems subject to dynamical noise and discuss whether we can specify a state $x_t$ by a permutation.

Although our theorem is restricted to one-dimensional interval dynamics, we use the Hénon map to observe whether our claim holds for higher dimensional dynamical systems. We adopt $x_t$ as the observed time series. Our numerical approach aims to represent a time series using a series of permutations (see Fig. 1 for the intuitive illustration). Thus, by specifying a permutation, we can eventually specify both series of states as well as stochastic inputs, simultaneously. Since a permutation eventually specifies a state in the limit for the length of permutation approaching infinity, representing points sharing the same permutation with a point becomes a reasonable approximation. We generated two time series $x$ and $x^{\prime }$ of the length $N$ from the same system. Here, we set $N=1000000$. For each time series $x$ and $x^{\prime }$, we also obtained a series of permutations $\{{\pi }_t(n)\}$ and $\{{{\pi }^{\prime }}_t(n)\}$ by using permutations of length $n$. Then, we use the first time series $x$ and its series of permutations $\{{\pi }_t(n)\}$ to compute the mean state

$$M_{\pi }^x(\kappa )=\frac{1}{|\mathrm{\#}\{{\pi }_t=\pi |t=1,2,\dots ,N\}|}\sum _{t:{\pi }_t=\pi }{x}_{t+\kappa }$$

(8)

for the $\kappa$th point of each appearing permutation $\pi$. This step is similar to a step of the k-means algorithm.³⁶36. N. Gershenfeld, The Nature of Mathematical Modeling (Cambridge University Press, Cambridge, 1998). Especially, we define the middle point $K=\lfloor n/2\rfloor$ of the permutation. These mean states become our estimates for states corresponding to a permutation $\pi$ since each permutation corresponds to a single initial condition as well as a single series of the external noise in the limit of the size of the permutation approaching to infinity. These means enable us to represent the second time series $x^{\prime }$ by replacing each permutation ${{\pi }^{\prime }}_t(n)$ with the corresponding mean $M_{{{\pi }^{\prime }}_t(n)}^x(\kappa )$ state of the $\kappa$th point of permutation ${{\pi }^{\prime }}_t(n)$ obtained from the first time series $x$. (If $\kappa \ge n$, then such an estimation becomes a time series forecast.) Finally, we evaluate the estimation error

$${\epsilon }^x(\kappa -n+1)=E_t\left[|{x^{\prime }}_{t+\kappa }-M_{{{\pi }^{\prime }}_t(n)}^x(\kappa )|\right]$$

(9)

by comparing the second time series $x^{\prime }$ and its representation $\{M_{{{\pi }^{\prime }}_t(n)}(\kappa )|t=1,2,\dots ,N-n+1\}$ constructed by permutations. This approach is called the mean representation.⁶6. Y. Hirata, K. Judd, and D. Kilminster, Phys. Rev. E 70, 016215 (2004). https://doi.org/10.1103/PhysRevE.70.016215

We first estimated states using the mean representation method (Fig. 2). We found that states corresponding to the time period of permutations were estimated more accurately than the cases where we just used the mean states for all over the points of the entire time series.

Moreover, we found a general tendency in the estimation error convergence to $0$ with increasing length of permutations [Figs. 3(a) and 3(b), for the logistic map³⁴34. R. May, Nature 261, 459 (1976). https://doi.org/10.1038/261459a0 and the Hénon map,³⁵35. M. Hénon, Commun. Math. Phys. 50, 69 (1976). https://doi.org/10.1007/BF01608556 respectively]. When we rigorously compared the model of exponential decrease converging to a constant in the limit of $n\to \mathrm{\infty }$ with the model of exponential decrease converging to zero using the Akaike information criterion,³⁷37. H. Akaike, IEEE Trans. Autom. Control 19, 716 (1974). https://doi.org/10.1109/TAC.1974.1100705 the model of exponential decrease converging to zero was selected for both cases (see Fig. 4). Thus, these results imply that a permutation corresponds to an initial condition, or a state, for these models.

Furthermore, the mean representation method was extended so that we considered the means for $q$ steps ahead by defining $\kappa =n-1+q$ in Eqs. (8) and (9) to make them “forecasts.” Then, we found that we could forecast short-term prediction horizons up to $10$ and four steps ahead better, in the noisy logistic map [Fig. 2(a)] and the noisy Hénon map [Fig. 2(b)], respectively, than the method of control where we considered the simple means over all the points of the attractor. The accuracy of these forecasts was achieved because the permutations could specify the past states and noise realization, even though uncertainty was generated due to the current and future parts of dynamical noise as well as the sensitive dependence on the initial conditions.

Theorem 1 implies that we can estimate the realization of the external force when there is a mathematical model for the dynamical system. However, our results imply that even if such a model $f$ is not available, we can estimate the realization of external force in the following two cases:

1.	Dynamics of the external force is gradual compared with the driven intrinsic dynamics.
2.	We observe multiple time series whose intrinsic dynamics is governed by identical dynamics and are subject to a common external force.

Therefore, Theorem 1 provides a way of obtaining the distribution of external noise as well.

To evaluate this possibility, we have estimated the dynamical noise by using the mean representation method and recurrence plot^32,3332. J.-P. Eckmann, S. O. Kamphorst, and D. Ruelle, Europhys. Lett. 4, 973 (1987). https://doi.org/10.1209/0295-5075/4/9/00433. N. Marwan, M. C. Romano, M. Thiel, and J. Kurths, Phys. Rep. 438, 237–329 (2007). https://doi.org/10.1016/j.physrep.2006.11.001 (see the Appendix). We found that the mean representation method achieved lower estimation errors than merely using the same means over all the permutations when the time steps are those corresponding to the positions of permutations (Fig. 5). In addition, the estimation errors for the noise realization also tend to decrease as the length of permutations increases (Fig. 6). Although we assumed here that a time series of noise realization for modeling is available, this result implies that we could narrow down the possible realization of dynamical noise using permutations. This direction enables us to construct a random dynamical system model.

When we used recurrence plots, we tested with three systems: the logistic maps with the common dynamical noise ${\eta }_t$. For an additive noise case, we use

$$x_{i,t+1}=3.7x_{i,t}(1-x_{i,t})+{\eta }_t$$

(10)

for ${\eta }_t\in [0,0.05]$; for a multiplicative noise case, we use

$$x_{i,t+1}=(3.7+{\eta }_t)x_{i,t}(1-x_{i,t})$$

(11)

for ${\eta }_t\in [-0.2,0.2]$.

In the Hénon maps, we use

$$\begin{array}{rl}{y}_{i,t+1}& =1-1.2y_{i,t}^2+0.3z_{i,t}+{\zeta }_t,\end{array}$$

(12)

$$\begin{array}{rl}{z}_{i,t+1}& =y_{i,t}\end{array}$$

(13)

for the common additive dynamical noise ${\zeta }_t\in [-0.1,0.1]$ or

$$\begin{array}{rl}{y}_{i,t+1}& =1-(1.2+{\zeta }_t)y_{i,t}^2+0.3z_{i,t},\end{array}$$

(14)

$$\begin{array}{rl}{z}_{i,t+1}& =y_{i,t}\end{array}$$

(15)

for the common multiplicative dynamical noise ${\zeta }_t\in [-0.1,0.1]$.

Here, we also use a chaotic neuron model³⁸38. K. Aihara, T. Takabe, and M. Toyoda, Phys. Lett. A 144, 333–340 (1990). https://doi.org/10.1016/0375-9601(90)90136-C to examine a more complicated situation. A chaotic neuron model is an extension of the Nagumo and Sato’s neuron model³⁹39. J. Nagumo and S. Sato, Kybernetik 10, 155–164 (1972). https://doi.org/10.1007/BF00290514 by replacing the Heaviside function with the sigmoid function, which defines whether or not a neuron fires. In a chaotic neuron model,³⁸38. K. Aihara, T. Takabe, and M. Toyoda, Phys. Lett. A 144, 333–340 (1990). https://doi.org/10.1016/0375-9601(90)90136-C we use

$$w_{i,t+1}=0.5w_{i,t}-\frac{1}{1+e^{-w_{i,t}/0.04}}+0.24+0.02{\sigma }_t$$

(16)

for the common additive dynamical noise ${\sigma }_t\in [-0.02,0.02]$ or

$$w_{i,t+1}=0.5w_{i,t}-\frac{1}{1+e^{-w_{i,t}/(0.04+0.01{\sigma }_t)}}+0.24$$

(17)

for the common multiplicative dynamical noise ${\sigma }_t\in [-0.01,0.01]$. Namely, states $w_{i,t}$ for multiple neurons at time $t$ are forced by the common dynamical noise ${\sigma }_t$. Here, we assume that each of ${\eta }_t$, ${\zeta }_t$, and ${\sigma }_t$ follows the independent and identical uniform noise, respectively. The number $L$ of the maps were decided as the minimum number in the form of $10\times 2^n$ with whose corresponding network, each time point is connected with all the other time points within ten steps. We assigned each of the initial conditions $x_{i,0},y_{i,0},z_{i,0},w_{i,0}$ by the uniform distribution between $0$ and $1$ and generated a time series of length $500$ each. We repeated this simulation 30 times to examine the robustness of our findings. We identified the tendency for longer permutations to perform more effectively in estimating the underlying dynamical noise [see panels (a), (c), and (e) of Figs. 7 and 8] for the additive dynamical noise and multiplicative dynamical noise, respectively. Examples shown in panels (b), (d), and (f) demonstrate that the dynamical noise reconstructed via recurrence plots agreed well with the true dynamical noise.

Overall, these numerical simulations confirmed that a permutation can specify a realization of the external forces, especially dynamical noise, which can be fast in its time scale.