Reconstructing dynamics of complex systems from noisy time series with hidden variables

The Lorenz system is a famous deterministic dissipative system [61], which is written as

$$\begin{align} \begin{aligned} \dot{x}&=a(y-x), \\ \dot{y}&=r x-y-x z, \\ \dot{z}&=x y-b z \,. \end{aligned} \end{align} \tag{ 10 }$$

With the often-chosen parameters $a = 10,b = 2.6,r = 28$, this system displays very typical chaotic behaviour with the maximum Lyapunov exponent close to 1. If the time series of $x,y,z$ are all known, it is not hard to recover the coefficients of the polynomials on the right hand side of equation (10) [20, 22, 35, 42, 43]. Instead, if only one variable, say x, is observed, the recovery becomes much harder [59]. Here, the new formulation is applied with two hidden variables $y,z$, three unknown parameters $a,r,b$ and two auxilliary parameters $e_2,e_3$. Thus, according to the above discussion, the total number of equations is $n_t = 2\times(1+3+2)+(3+2)\times 2 = 22$. In the following simulations, $[a,r,b] = [1,1,1], [y,z] = [1,1]$ are conveniently used as the initial conditions for the corresponding variables. With the above algorithm, the three unknown parameters $a,r,b$ converge well as shown in figures 2(A)–(C)and 3(C), (D). For different hyper-parameters, the convergence seems to follow similar patterns. Initially, the unknown parameters change steadily but at some point start to rush to the true values and then remain there throughout the simulation. Nevertheless, the convergence is different in disparate conditions and could even go awry if improper hyper-parameters are chosen as displayed in figures 2(D) and 3(A), (B).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this example, the integration time step is $t_{\mbox{unit}} = 0.01$. As can be seen from the figures, a total number of $N_{step}\sim 10^5$ steps are needed for the convergence, which is determined by the two hyper-parameters $\alpha, \gamma$. Currently, we take α = 3, indicating a memory of $1/\alpha = 0.33$ which is about a third of the Lyapunov time of the system. The learning rate is chosen to be γ = 0.01, about one hundredth of the Lyapunov exponent. γ is used to control the update of the parameters in the evolution equation and should keep small compared with the time rate of the system state (y, z) in order not to disrupt the generic dynamics of the original system. As shown in figures 2(A)–(C), a reduction of γ extends the learning time and requires more data for convergence. Nevertheless, if the available time series has a limited length, a technique of data recycling may be employed for a successful reconstruction as shown in appendix 'Recycling of data'. On the other hand, if γ is overly large, the parameters change rapidly and the original dynamics may be altered so that no convergence results, as displayed in figure 2(D). So during the reconstruction, the learning rate should be chosen properly to ensure both the stability and efficiency.

The decay rate α should not be too small in order to bound the integrals in equation (6). On the other hand, if α is too large, the memory becomes really short and these integrals highly depend on the current state which varies rapidly and easily overruns the slow parameter changes. In figure 3, the influence of the decay rate α on the convergence is displayed. As α increases, the learning becomes slow and more data is needed, just like reducing γ. When α is relatively small (figure 3(A)), the effective integration time becomes long and the integral will be very sensitive to parameter changes, causing instability in the system, unless the learning rate γ is really small and averaging out the fast oscillations. Therefore, for a system with a large Lyapunov exponent, when the learning rate is fixed, it is necessary to apply a proper decay rate to reduce the sensitivity of the parameters (appendix figures 21 and 22). In general, it is important to adjust the learning and the decay rate to keep a balance between the two so as to achieve a stable and fast convergence.

From the above discussions, it seems that the current reconstruction scheme is quite fragile. Actually, it is not and able to work quite well even in the presence of noise. Next, we investigate the impact of noise on the reconstruction since it inevitably exists in every experiment or observation. For this purpose, we add a noise term to each of the evolution equation (10), which describes white noise with the following characteristics

$$\begin{align}<\eta_{i}(t)> &= 0,\kern22.6pt \end{align} \tag{ 11a }$$

$$\begin{align} <\eta_{i}(t)\eta_{j}(t^{^{\prime}})> &= 2D\delta_{ij}(t-t^{^{\prime}}), \end{align} \tag{ 11b }$$

where $\langle \cdot \rangle$ denotes an ensemble average and the delta function is a Kronecker one for the discrete index $i, j$, and a Dirac one for the continuous time variables $t-t^{^{\prime}}$. D is the noise intensity and when D = 0 we recover the deterministic autonomous system.

Even though a noisy orbit is generated with equation (10) with noise (D ≠ 0), we still use the same deterministic formulation in the previous section. As an example, suppose that only the noisy $x-$orbit is given and we still aim to deduce the unknown parameters and the hidden variables $y(t),z(t)$. If the noise is not very strong, x(t) may be directly substituted into the 22 reconstruction evolution equations mentioned above, together with the $\dot{x}(t)$ computed with any fair smoothening technique [43]. As a result, only the average orbit is reconstructed with the unknown parameters. Figure 4 shows the reconstruction results at different noise intensities, i.e. D = 0.1 in figures 4(A) and (B) and D = 1 in (C) and (D). The convergence pattern of the parameters displayed in figures 4(A) and (C) looks very similar to that in the deterministic case but small jittering exists due to noise perturbation, which is more clearly seen in (C) for D = 1. After the small oscillations are averaged out, we get a fairly good estimate of the parameters, which is listed in table 1. In fact, the results of extra tries are also displayed in the table up to D = 5. In general, the estimation gets worse for greater noise intensity as expected but the errors do not appear monotone. The reconstructed orbits $y(t),z(t)$ are displayed in figures 4(B) and (D), which seem almost overlapping with the true orbits, just as in the deterministic case.

Zoom In Zoom Out Reset image size

From the above analysis, it can be seen that the new construction scheme is quite robust against noise as long as it is not too strong. Next, the scheme will be applied to several typical examples in nonlinear dynamics with high dimensions to show its validity.

First, we consider the inference problem in a network with N nodes each of which is occupied by a Lorenz system. The connection between nodes marks their interaction. More specifically, the ith node is governed by the equation below

$$\begin{align} \dot x_{i} &=a_{i}(y_{i}-x_{i})+\eta_{i}(t), \notag\\ \dot y_{i} &=r_{i}x_{i}-y_{i}-x_{i}z_{i}+\sum_{j=1}^N A_{ij}(y_{j}-y_{i})+\eta_{i}(t), \notag \\ \dot z_{i} &=x_{i}y_{i}-b_{i}z_{i}+\eta_{i}(t)\, \end{align} \tag{ 12 }$$

where $\eta_i(t)$ is the above-mentioned Gaussian white noise and A_ij is the interaction matrix. Here, without less of generality, the interaction is assumed to take place between $y-$components. The local parameters are randomly selected from intervals on the positive real axis: $a_{i}\in [10,12], r_{i}\in [28,30], b_{i}\in [2.6,2.8]$. The noise intensity D = 0.01 and the sampling interval $t_ {\text{unit}} = 0.01$. Suppose that only the $x-$components are measured. Our goal is to reconstruct the coupling matrix A_ij and all the local parameters $\{a_i,b_i,r_i\}_{i = 1,2,\ldots,N}$ as well as the time-dependent hidden variables.

With the above scheme, for an Erdős-Rényi (ER) network with N = 16, the coupling matrix and the node parameters are correctly deduced as displayed in figure 5. For simplicity, we take the initial conditions of interaction matrix A_ij to be zero unless otherwise stated, and the initial values of local parameters are $[a_{i},r_{i},b_{i}] = [1,1,1], [y_{i},z_{i}] = [1,1]$. After the reconstruction, the absolute errors of the computation in the coupling matrix are of the order of 10⁻³, and around 10⁻² in the node parameters. The time series of all hidden variables are also reconstructed at the same time, which is not shown here for brevity. Please check figure 17 in the appendix for one example.

Zoom In Zoom Out Reset image size

When N is small, we may start with a full matrix $[A_{i,j}]$ containing N² elements to be determined. The network structure can be drawn from the non-zero elements of $[A_{i,j}]$ after the reconstruction and we do not need any prior information about the network topology. Nevertheless, this type of analysis is not scalable since the quadratic increase of the number of the unknown matrix elements with the network size prevents their efficient estimation based only on the linearly increasing information. In fact, the current brute-force method works up to $N \sim 20$. For larger networks, extra information about the network topology helps a lot, especially those that cut down the number of unknown matrix elements. For example, if the structure of the network is approximately known (which may still contain some spurious links), our computation remains valid as long as the number of links increases not that fast.

Here, we tried the reconstruction on a small-world network with N = 100 and an average degree K = 4 as portrayed in figure 6(A). For an easy comparison, the non-zero coupling is allowed taking only four values $0.1,0.3,0.5,0.8$. The results are displayed in figures 6(B) and (C) which agree well with the original values. The convergence of the local node parameters looks quite similar as in precious cases (see figure 6(B)). However, the evolution of the couplings starts more erratically. The fluctuations grow larger and larger and culminate just before a sudden convergence to the true values (see figure 6(C)). The whole process closely resembles what happens in chaos control: only when the system state gets close to the stable manifold of the true dynamics and the true parameters does the gradient descent scheme in equation (8) start to work and drive the wandering point to the right solution [62].

Zoom In Zoom Out Reset image size

In the current example, we have 700 unknown system parameters, including 300 node parameters and 400 coupling strengths. In addition, 200 auxiliary parameters e are needed to drive $y_i(t)\,,z_i(t)$ to the right trajectories. Altogether, we have 900 hundred parameters, which will generate over 180 000 equations in the new formulation, if we exactly follow the procedures in the theory part. However, a convenient approximation could be adopted based on the given network structure and the observation that the influence of a node on other nodes decays quickly in an asynchronous state. According to the equation of motion, each of the above-mentioned parameters could be assigned to a node or a connection between nodes. To the lowest order, it is reasonable to assume that the dynamics at a node mainly depends on the parameters of neighbouring nodes or edges. With this assumption, the partial derivatives with parameters are taken to be zero if they locate outside the immediate neighbourhood. As a result, the number of equations is greatly cut down, which in fact grows only linearly with the size of the network with a fixed average degree.

As the true values of the parameters or state variables are obtained through a co-evolution with the observed time series, the proposed method can be applied online to monitor changes in network topology or node dynamics if the change does not happen so frequently. As an example, the above formulation is directly applied to the network with 8 nodes displayed in figure 17(A) in the appendix. Starting from a somewhat arbitrary initial guess value 1, the coupling converges to the true value 0.1 after a quite long transient. Later on, as the coupling is switched to 0.5 in the network, the reconstruction algorithm detects and quickly follows this change. It is interesting to see that the reaction is really fast this time and the transient is almost unnoticeable as shown in figure 7(A). The node parameters that change at the same time can be similarly observed, as shown in figure 7(B). One reason for the quick response for the parameter change is that the state variables $\{y_i(t),z_i(t)\}_{i = 1,2,\ldots,N}$ remain close to the true values in this process so that the convergence seems monotone and exponential!

Zoom In Zoom Out Reset image size

A salient feature of the FitzHugh–Nagumo (FHN) system is the excitability—once the potential exceeds a threshold, the system state immediately shoots to a very high value and then relaxes back to the base. Therefore, the temporal pattern is quite different from that of a Lorentz system. Is it still possible to apply our method to this new situation? Below, we check the continued validity of the current scheme.

The FHN equation is a second-order nonlinear system used to simulate the dynamics of a neuron. To model neural networks, these neurons connect to each other and interact with their neighbours. The equations of motion can be written as

$$\begin{align} \begin{aligned} &\dot{v}_{i}=v_{i}-v_{i}^{3}-w_{i}-\sum_{j=1}^{N} A_{i j}\left(v_{j}-v_{i}\right)+\Gamma_{1,i} \\ &\dot{w}_{i}=a_{i}+b_{i}+c_{i} w_{i}+\Gamma_{2,i} \end{aligned} \end{align} \tag{ 13 }$$

where $i = 1,2,\ldots,N$ denotes different neurons, v_i is the rate of change of the membrane potential, w_i is the slowly changing recovery variable of the ith neuron and N is the number of neurons. In this model, the dynamical behaviour of the system largely hinges on the coefficients $a_{i}\in[0.3,0.5],b_{i}\in[0.5,0.7],c_{i}\in[-0.04,-0.03]$ and the interaction weight $A_{i j}\in[0,1]$. $\Gamma_{1,i}, \Gamma_{2,i}$ are noise terms of the two components, which are usually taken to be white satisfying equation (11), with D = 0.01.

In practice, only the potential $v_{i}(t)$ of each neuron is subject to a direct measurement. Is it possible to recover $w_{i}(t)$ and all the parameters in the equation with this information? figure 8(A) displays a ER network with N = 25 from which the measured data was obtained. As usual, we apply the above scheme by initially setting all unknown coupling strength A_ij to 1 and the local system parameters to 0. With this network size, all the parameters could be obtained directly without other assumptions. The obtained coupling strengths are depicted in figure 8(B) with the absolute error of the order of 10⁻², while the local parameters are portrayed in figure 8(C) with the absolute error of the order of 10⁻³. Figure 9(A) shows the convergence of local parameters for a node from the initial values 0 to the correct values, where the oscillatory character at the beginning is still clearly seen. Figure 9(B) displays the time series of $v_{i}(t)$ where the profile is contaminated partly by visible noise. However, the reconstructed time series of w_i seems smooth and overlaps with the true one almost perfectly.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For small-world FHN networks with known structures, for example, the one in figure 6(A), how does the reconstruction accuracy depend on the network size? If the average degree is fixed to 4, this dependence is shown in figure 10(A). It is easy to see that the inference accuracy can be maintained above $95\%$ after averaging over 20 realizations. The accuracy is defined as [20]

$$\begin{align} Ac_{\rho}=\frac{1}{K \times N} \sum_{i, j} H\left((1-\rho)-\Delta A_{i j}\right) \,, \end{align} \tag{ 14 }$$

where ρ is the required accuracy and we take ρ = 0.1 here. H is the Heaviside function and the relative error between the reconstructed ($\hat{A}_{i j}$) and the real coupling ($A_{i j}$) is:

$$\begin{align} \Delta A_{i j}=\left|\hat{A}_{i j}-A_{i j}\right| /\left(A_{i j}\right) \,. \end{align} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

The main source of error is the synchronization among part of the nodes, which makes it hard to distinguish them and thus leads to incorrect inferences. More discussion will be made below. On the other hand, if we keep the network size and increase the average degree, the inference will become harder since more couplings need deducing. Nevertheless, according to figure 10(B), the accuracy still maintains near the value $95\%$ throughout.

The ability to deduce a large number of parameters from the time series may be applied to the reconstruction of networks with hidden nodes, about which little information is known. Even their existence has to be deduced from the measurement on neighbouring nodes as done in previous works to determine the possible connections of hidden nodes to the known ones in the network [51–54]. Nevertheless, this type of deduction is able to confirm the existence of hidden nodes but the precise connections may not be inferred very accurately. In our scheme, all possible links may be included in the deduction as long as the node's existence is expected. The spurious ones will give a weight near zero and the weights will be recovered for the true ones. Also, the dynamics and internal parameters of the hidden nodes could be deduced as well.

Here, Kuramoto oscillator model [63] defined on networks will be used as an example for this computation. Each oscillator in the model is described with a phase, so D = 1. The coupling between them induces synchronization at various levels, which has been studied comprehensively [64, 65]. Consider a system with N oscillators

$$\begin{align} \dot{\theta}_{i}=\omega_{i}+\frac{K}{N} \sum_{j=1}^{N} a_{i j} \sin \left(\theta_{j}-\theta_{i}\right),(i=1, \ldots, N) \,, \end{align} \tag{ 16 }$$

where θ_i and ω_i are the phase and the natural frequency of the ith oscillator, K > 0 is the coupling strength and a_ij is the adjacency matrix of the network. When node i and node j are connected, $a_{i j} = 1$; otherwise $a_{i j} = 0$.

To characterize synchronization, in the literature [65], an order parameter R is defined as follows

$$\begin{align} R e^{i\psi}=\frac{1}{N}\sum_j e^{i \theta_j} \,, \end{align} \tag{ 17 }$$

where $R,\psi \in \mathbb{R}$. If $K\ll 1$, each oscillator rotates near its natural frequency and $R\sim 0$. With increasing K, the average frequencies of connected oscillators approach the average natural frequency in general and some of them may start to rotate together, reaching a partial synchronization such that R > 0. For sufficiently large K, all the oscillators become synchronous and rotate with a common frequency $\sum_i \omega_i/N$, so that $R \sim 1$. During the course of dynamics or parameter deduction, synchronization is harmful as seen previously and will be explained below.

In the current example, 32 oscillators are connected to form an ER network as shown in figure 11(A) and ω_i 's are picked up from a uniform distribution in the interval $[0,20]$. The coupling strength K = 1 is not very large to prevent synchronization. We assume that four hidden nodes exist in the network. Their dynamics and natural frequencies are not revealed to us, while the time series of the remaining nodes are known. First, the very existence could be deduced together with their rough connections with existing methods in the literature [51]. Next, the whole network structure and the natural frequencies of these hidden nodes are reconstructed with the formalism introduced in this manuscript. In the following, the hidden nodes are randomly selected from the network with no neighbouring pair.

Zoom In Zoom Out Reset image size

In this example, as the number of nodes is not large, as long as the very existence of hidden four nodes are known, we may start with a full connection matrix in the reconstruction. That is, the initial values of all entries of $[a_{ij}]$ are assumed to be 0.5,and the natural frequency ω_i 's start with the value 1. In figure 11(B), a comparison is made between the simulation results and true values of the natural frequencies. The error seems small, not exceeding a few percent. The convergence of the natural frequencies of the hidden nodes is displayed in figure 11(C). The familiar oscillatory profiles appear again but they decay much faster than in previous examples. In figure 12, it is easy to see that the accuracy of the reconstruction (see equation (14)) decreases with the increase of hidden nodes. The larger the network, the faster the decay. Nevertheless, from the figure we see that when the hidden nodes are less than $15\%$, no matter how large the network size is, the accuracy is higher than $90\%$.

Zoom In Zoom Out Reset image size

In fact, our method is even able to infer the structure of a cluster of hidden nodes. In figure 13, there are 20 hidden nodes (green dots) forming a small cluster. The structure as well as all the parameters is unknown and none of them is observed directly.We can only measure their neighbouring nodes (black ones). From the measured time series, the natural frequencies of nodes in the hidden cluster and the coupling strengths among nodes can be inferred quite accurately, which is clearly displayed in figure 14(A). A full connection between hidden nodes is assumed at the start and the algorithm will computed the actual strengths of the links: about 0 for the non-existing ones and 0.5 for the true connections. The convergence seems fast (see figure 14(B)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

If the network is fully synchronized, all nodes behave in essentially the same way. The time series look similar apart from a phase or a scaling factor, so it is hard to carry out the inference (see appendix figure 23). Generally speaking, the more synchronized the network, the harder it is to infer the structure or coupling. Taking the Kuramoto oscillator network as an example, the natural frequency ω of each node follows the uniform distribution $W(0,\mu)$, and a larger µ indicates that the dynamics has higher heterogeneity. As shown in the figure 15, for $\mu\in [1,10]$, with the increase of µ, the inference accuracy increases significantly. When µ reaches 10, the accuracy is over $95\%$, and fluctuates around this value upon further increment of µ.

Zoom In Zoom Out Reset image size

In order to address in more detail the effect of synchronization on inference, we study a network of six oscillators coupled in a specific configuration as shown in figure 16(A), two of which are assumed to be hidden nodes. The natural frequencies of the oscillators conform to a uniform distribution in $[0,1]$ and the coupling strength is set to 3. In figure 16(B), the coupling strengths indeed converge, but to completely wrong values, since in this case the oscillators fully synchronize. Nevertheless, if the distribution of the natural frequencies is expanded to a larger interval $[0,30]$ so that synchronization is not reached, our algorithm is able to dig out the coupling strengths as shown in figure 16(C). To quantify the performance of the algorithm at different levels of synchronization, we employ two statistical indicators—the order parameter $0\leqslant R \leqslant 1$ defined in equation (17) and the G-P dimension [66] for the oscillator network. The order parameter reflects the degree of synchronization, which is between 0 (fully random) and 1 (fully synchronous). The fractal dimension gives the amount of entropy contained in the time series. The larger the value, the more entropy the time series contains.

Zoom In Zoom Out Reset image size

When the natural frequency is distributed in $[0,1]$, $\lt R\gt = 0.97$ as given in figure 16(D), which is close to a full synchronization. Accordingly, the G-P dimension is only 1.45, indicating that the time series contains very little variability besides the nearly periodic oscillation, which makes it hard to distinguish different nodes. However, the entropy contained in the time series considerably increases if the frequency distribution is extended to the interval $[0,30]$, since the order parameter and the G-P dimension turn 0.38 and 5.2 respectively as plotted in figure 16(E), signalling little synchronization in the network. To see clearly the influence of heterogeneity induced by the frequency distribution $[0,\mu]$, we scan different values of µ and depict the results in figure 16(F). The order parameter of the system decreases while the normalized dimensionality gradually increases, which indicates the waning of the synchronization. As expected, the accuracy of the reconstruction is gradually increasing and approaching 100%.