Data-driven Reconstruction of Nonlinear Dynamics from Sparse Observation

TL;DR

This paper introduces a fixed-point echo state network method to reconstruct complex nonlinear dynamical systems from extremely sparse time series data, demonstrating high accuracy even with minimal observations.

Contribution

It proposes a novel fixed-point ESN approach for reconstructing nonlinear dynamics from sparse data, leveraging the universal approximation capability of ESNs.

Findings

Reconstructs chaotic systems from up to 95% missing data

Successfully reconstructs simple systems from only 1-2% data

Outperforms traditional methods in data-scarce scenarios

Abstract

We present a data-driven model to reconstruct nonlinear dynamics from a very sparse times series data, which relies on the strength of the echo state network (ESN) in learning nonlinear representation of data. With an assumption of the universal function approximation capability of ESN, it is shown that the reconstruction problem can be formulated as a fixed-point problem, in which the trajectory of the dynamical system is a fixed point of the ESN. An under-relaxed fixed-point iteration is proposed to reconstruct the nonlinear dynamics from a sparse observation. The proposed fixed-point ESN is tested against both univariate and multivariate chaotic dynamical systems by randomly removing up to 95% of the data. It is shown that the fixed-point ESN is able to reconstruct the complex dynamics from only 5 ~ 10% of the data. For a relatively simple non-chaotic dynamical system, the numerical experiments on a forced van der Pol oscillator show that it is possible to reconstruct the nonlinear dynamics from only 1~2% of the data.

Tables

Table 1. Table 1: Free parameters of ESN.

$\lambda$	Temporal relaxation parameter
$\beta$	$L_2$-regularization parameter
$\nu$	Fraction of nonzero elements of $𝑨$
${\rho }_{max}$	Spectral radius of $𝑨$
${\xi }_A$	$A_{ij}\sim 𝒰(-{\xi }_A,{\xi }_A)$
${\xi }_B$	$B_{ij}\sim 𝒰(-{\xi }_B,{\xi }_B)$

Table 2. Table 2: Relaxation ( α 𝛼 \alpha ) and regularization ( β 𝛽 \beta ) parameters of the fixed-point ESN

$\omega$	0.60	0.80	0.90	0.95
$\alpha$	0.2	0.2	0.2	0.4
$\beta$	${10}^{-9}$	${10}^{-9}$	${10}^{-9}$	${10}^{-8}$

Table 3. Table 3: Relaxation ( α 𝛼 \alpha ) and regularization ( β 𝛽 \beta ) parameters and NRMSE for the partial observation of the Lorenz-63 time series.

$\omega$	$\alpha$	$\beta$	${\sigma }_{lin}$
0.90	0.2	${10}^{-7}$	0.06
0.95	0.4	${10}^{-6}$	0.37

Table 4. Table 4: Normalized RMSE of the multivariate reconstruction of the Lorenz time series.

$\omega$	0.90	0.95
${\sigma }_{lin}$	0.18	0.51
${\sigma }_{csp}$	0.17	0.36

Equations

$${\mathbf{s}}_{t+1}=\lambda {\mathbf{s}}_t+(1-\lambda )\Psi ({\mathbf{s}}_t,{\mathbf{x}}_t).$$

$$\Psi ({\mathbf{s}}_t,{\mathbf{x}}_t)=\tanh (\mathbf{A}{\mathbf{s}}_t+\mathbf{B}{\mathbf{x}}_t).$$

$$\frac{{\mathbf{s}}_{t+1}-{\mathbf{s}}_t}{\delta t}=\frac{1-\lambda }{\delta t}\left(-{\mathbf{s}}_t+{\mathbf{f}}_t\right),$$

$${\gamma }^{-1}\frac{d\mathbf{s}}{dt}=-\mathbf{s}(t)+\mathbf{f}(t).$$

$${\mathbf{y}}_{t+1}={\mathbf{\theta }}^T{\mathbf{S}}_{t+1},$$

$$\underset{\mathbf{\theta }}{\min }\sum _{i=1}^N\frac{1}{2}\Vert {\mathbf{y}}_i-{\mathbf{\theta }}^T{\mathbf{S}}_i{\Vert }_2^2+\frac{\beta }{2}\Vert \mathbf{\theta }{\Vert }_F^2,$$

$$\mathbf{\theta }={\left(\sum _{i=1}^N{\mathbf{S}}_i{\mathbf{S}}_i^T+\beta \mathbf{I}\right)}^{-1}\left(\sum _{i=1}^N{\mathbf{S}}_i{\mathbf{y}}_i^T\right).$$

$${\mathbf{s}}_{t+1}$$

$${\mathbf{y}}_{t+1}$$

$$\frac{d\mathbf{y}}{dt}=\mathcal{F}(\mathbf{y},\mathbf{u}),t\in [0,T\delta t],$$

$$y_{l_j}^{\ast }=y_j(l\delta t),\text{for}l=0,\cdots ,T,\text{and}j=1,\cdots ,N_y.$$

$$y_{l_j}^O=\{\begin{array}{ll}{y}_{l_j}^{\ast }& \text{if observed,}\\ \eta & \text{otherwise.}\end{array}$$

$$m_{l_j}=\{\begin{array}{ll}0& \text{if}{y}_{l_j}^O=\eta ,\\ 1& \text{otherwise}.\end{array}$$

$$\omega =\frac{1}{T\times N_y}\sum _{t=1}^T\sum _{j=1}^{N_y}(1-m_{t_j}),$$

$${\mathbf{y}}_t^{\ast }={{\mathbf{\theta }}^{\ast }}^T{\mathbf{S}}_t^{\ast },\text{for}t=1,\dots ,T$$

$${\mathbf{s}}_t^{\ast }=\lambda {\mathbf{s}}_{t-1}^{\ast }+(1-\lambda )\tanh (\mathbf{A}{\mathbf{s}}_{t-1}^{\ast }+{\mathbf{B}}_y{\mathbf{y}}_{t-1}^{\ast }+{\mathbf{B}}_u{\mathbf{u}}_{t-1}),$$

$${\mathbf{Y}}^{\ast }=\mathcal{E}({\mathbf{Y}}^{\ast },\mathbf{U}),$$

$${\mathcal{E}}_0({\mathbf{Y}}^{\ast },\mathbf{U})$$

$${\mathcal{E}}_l({\mathbf{Y}}^{\ast },\mathbf{U})$$

$${\mathbf{\theta }}_j^k={\left(\sum _{t=T_s}^T{m}_{t_j}{\mathbf{S}}^k{{\mathbf{S}}^k}^T+\beta \mathbf{I}\right)}^{-1}\left(\sum _{t=T_s}^T{m}_{t_j}{y}_{t_j}^O{\mathbf{S}}_t^k\right),\text{for}j=1,\cdots ,N_y.$$

$${\mathbf{Y}}^{R,k+1}=\alpha {\mathbf{Y}}^{R,k}+(1-\alpha )\mathcal{E}({\mathbf{Y}}^{R,k},\mathbf{U})$$

$$e^k={\left[\frac{1}{N_{y}T}\sum _{j=1}^T\Vert {\mathbf{y}}_j^{R,k+1}-{\mathbf{y}}_j^{R,k}{\Vert }_2^2\right]}^{1/2}$$

$${\mathbf{y}}_t^R=(1-{\mathbf{m}}_t)\circ {\hat{\mathbf{y}}}_t+{\mathbf{m}}_t\circ {\mathbf{y}}_t^O,$$

$${\mathbf{Y}}^R=\mathcal{E}({\mathbf{Y}}^R,\mathbf{U}).$$

$${\mathbf{Y}}^{R,k+1}=\mathcal{G}({\mathbf{Y}}^{R,k},\mathbf{U})=\alpha {\mathbf{Y}}^{R,k}+(1-\alpha )\mathcal{E}({\mathbf{Y}}^{R,k},\mathbf{U}),\text{for}k\ge 0,$$

$${\mathbf{\theta }}^k=\underset{\mathbf{\theta }}{\min }\sum _{t=T_s}^T\frac{1}{2}\Vert {\mathbf{m}}_t\circ ({\mathbf{y}}_t^O-{\mathbf{\theta }}^T{\mathbf{S}}_t^k){\Vert }_2^2+\frac{\beta }{2}\Vert {\mathbf{\theta }}^k{\Vert }_F^2.$$

$${\mathbf{\theta }}_j^k={\left(\sum _{t=T_s}^T{m}_{t_j}{\mathbf{S}}^k{{\mathbf{S}}^k}^T+\beta \mathbf{I}\right)}^{-1}\left(\sum _{t=T_s}^T{m}_{t_j}{y}_{t_j}^O{\mathbf{S}}_t^k\right),\text{for}j=1,\cdots ,N_y.$$

$$\lambda =0.6,\nu =0.01,{\rho }_{max}=0.9,{\xi }_A=1,{\xi }_B=0.4,\text{and}{T}_s=200.$$

$$\frac{dy(t)}{dt}=\frac{{\gamma }_{1}y(t-\tau )}{1+y^{{\gamma }_2}(t-\tau )}-{\gamma }_{3}y(t),$$

$$\text{NRMSE}={\left[\frac{{\sum }_i(y_i^{R,k}-y_i^{\ast }{)}^2}{{\sum }_i{y_i^{\ast }}^2}\right]}^{1/2}.$$

$$e^k={\left[\frac{1}{T-T_s}\sum _{j=T_s+1}^T(y_j^{R,k+1}-y_j^{R,k}{)}^2\right]}^{1/2}.$$

Full text

1 Introduction

Over the last decade, artificial neural network (ANN) has been extensively studied due to its strength in identifying complex nonlinear structure in data [1, 2, 3]. Recent studies on the application of ANN on the data-driven modeling of physical systems have shown promising results [4, 5, 6, 7]. In particular, in [8, 9], it is shown that ANN is capable of predicting the dynamics of a chaotic attractor. When training an ANN for the modeling of a dynamical system, it is typical to use a sufficiently long time series data, so that the ANN fully explores the phase space. Given the huge amount of training data, it is unclear whether the long-term forecast capability of ANN is due to a high dimensional interpolation in a phase space, similar to the delay-coordinate embedding [10, 11], or because ANN really learns the nonlinear dynamics. As a first step to answer this question, we first consider the problem, where only a partial observation of a time series data is available. We investigate if ANN can reconstruct the nonlinear dynamics from the incomplete information.

Reconstructing dynamics from an incomplete observation has practical applications in many physical, biological, and engineering problems, where the state of the system is accessible only through a sensor network. For example, in geophysics, it is not uncommon to find time series data, which contains a large amount of missing information or is irregularly sampled, due to sensor malfunction, atmospheric conditions, or physical limitations of the sensors [12, 13, 14]. When a priori knowledge on the physical system is available, one of the standard approaches to reconstruct the nonlinear dynamics from an incomplete data set is to design a statistical model that incorporates the physical knowledge [15, 16, 17]. In [18], it is shown that a high-resolution temporal dynamics can be reconstructed from a low-resolution time series data, which has a 10 times lower time resolution, by using only a subset of the governing equations. Incorporating the physics knowledge on the training of ANN, [6] proposed “physics informed neural network”, which can recover the complex nonlinear dynamics from only a small fraction of the data.

It is challenging to reconstruct nonlinear dynamics without having any prior knowledge on the underlying process. There are three major approaches to reconstruct the missing information; (a) identifying / modeling (auto-)correlation structures [19, 20, 21], (b) low-dimensional modeling [22, 23, 24], and (c) pattern matching [25]. Although there is a significant progress in imputing the missing data, most of the classical approaches are based on strong assumptions, such as linearization, existence of a low-dimensional structure, or a priori knowledge on the correlation structures. The reconstruction accuracy of these models degrade quickly when the fraction of missing data is large [26], e.g., more than 50%, and the underlying dynamical system is nonlinear.

Recently, a few ANN-based methods have been proposed for the data-driven reconstruction of the missing information by using recurrent neural networks (RNN). [27] showed that simply filling in a missing period with the last known observation and using a missing data indicator as one of the input data to RNN can increase the prediction accuracy in a clinical binary classification problem. In [28], a Bidirectional-RNN based method is developed for a probabilistic reconstruction of the missing information for a binary time series data. [29] proposed an unidirectional RNN model, which assumes a relaxation of the unobserved dynamics towards an empirical mean of the data. Based on the relaxation model of [29], [30] and [31] developed Bi-directional RNN and Generative Adversarial Network models, respectively. Those ANN models showed better performance compared to the conventional statistical imputation approaches for the applications in the health care domain. However, it should be noted that the relaxation towards an empirical mean is essentially based on the assumption of a stationary Markov process, which does not consider temporal structures or dynamics of the underlying process.

In this study, we develop a novel RNN model for the data-driven reconstruction of nonlinear dynamics from temporally sparse observations. For the modeling of nonlinear dynamics, the echo state network (ESN) [32] is employed. The reconstruction problem is formulated as a fixed-point problem. The behavior of the proposed fixed-point ESN is thoroughly studied by using chaotic time series data. It is shown that ESN can reconstruct the chaotic, nonlinear dynamics from as little as 5% of the data. This paper is organized as follows; the echo state network is reviewed in section 2.1. In 2.2, a fixed-point iteration method is proposed. Numerical experiments of the fixed-point ESN model is presented in section 3. Finally, conclusions are given in section 4.

2 Echo State Network

2.1 Review of Echo State Network

The echo state network (ESN) has been studied extensively due to its strength in learning complex dynamics from time series data [5, 8, 33, 34]. ESN is a nonlinear state-space model, of which latent state consists of a set of randomly generated dynamical systems. Following [32, 35], the evolution equation of the latent state is

[TABLE]

Here, $\bm{s}_{t}\in\mathbb{R}^{N_{s}}$ is the latent state of ESN, $\bm{x}_{t}\in\mathbb{R}^{N_{x}}$ is a vector of the input signal to ESN, $\lambda\in[0,1)$ is a temporal relaxation coefficient, and $\Psi(\cdot)$ is a nonlinear function. The subscript $t$ denotes a time stamp, $\phi_{t}=\phi(t\delta t)$, where $\delta t$ denotes a time step size, or data sampling interval. The input signal may consist of the observation of the physical system, $\bm{y}_{t}\in\mathbb{R}^{N_{y}}$, and an exogenous forcing or auxiliary input, $\bm{u}_{t}\in\mathbb{R}^{N_{u}}$, such that $\bm{x}_{t}=(\bm{y}_{t},\bm{u}_{t})$. For the nonlinear function, usually a hyperbolic tangent function is used;

[TABLE]

Here, $\bm{A}\in\mathbb{R}^{N_{s}\times N_{s}}$ and $\bm{B}\in\mathbb{R}^{N_{s}\times N_{x}}$ are randomly generated weight matrices. Independent uniform random variables are used to generate $\bm{A}$ and $\bm{B}$ [32], e.g., $A_{ij},\,B_{ij}\sim\mathcal{U}(-\xi,\xi)$, in which $\mathcal{U}(a,b)$ denotes a uniform random distribution in $(a,b)$. It is worthwhile to note that the evolution equation of the standard recurrent neural network (RNN) has a very similar structure with (1–2). One of the most significant differences between RNN and ESN is that, in ESN, the weight matrices in (2) are generated randomly, while they are computed by a maximum likelihood method in RNN. The weight matrix, $\bm{A}$, decides the connectivity between the latent state, $\bm{s}$. It is suggested that making the connection sparse results in a richer internal dynamics [32]. Typically, $\bm{A}$ has only 1 $\sim$ 2% of nonzero elements. A more detailed explanation about generating $\bm{A}$ is given in section 3.

The echo state network solves an initial value problem. While it is straightforward to specify the initial condition of the physical variables, e.g., $\bm{y}_{0}$ and $\bm{u}_{0}$, it is challenging to find the correct initial condition of the latent state, $\bm{s}_{0}$. In [32], it is proposed that an ESN must have the “echo state property”, which is similar to a “fading memory effect” [35, 36]. For $\lambda>0$, it is straightforward to show the fading memory effect. Observe that (1) can be rewritten as

[TABLE]

in which $\bm{f}_{t}=\Psi(\bm{s}_{t},\bm{x}_{t})$. If we choose $\lambda=1-\gamma\delta t$ for $0<\gamma<1/\delta t$, it is clear that the evolution equation of the latent state, (1), is the forward Euler approximation to the system of relaxation equations,

[TABLE]

Since $|\bm{f}(t)|<\infty$, (4) implies that the effects of the misspecification of $\bm{s}_{0}$ will vanish after some initial transient period. It is also clear that ESN relies on a set of randomly generated relaxation processes to model the nonlinear behaviors in the time series data. When $\lambda=0$, the “echo state property” can be guaranteed by making the spectral radius of $\bm{A}$ smaller than unity, i.e., $\rho(\bm{A})<1$ [35].

Once the latent state of ESN is updated by (1), the physical process at the next time step ($t+\delta t$) is computed by a linear projection of the latent state onto the physical space,

[TABLE]

in which $\bm{S}_{t+1}\in\mathbb{R}^{N_{s}+1}$ is an augmented latent state, $\bm{S}_{t}=(1,\bm{s}^{T}_{t})^{T}$ and $\bm{\theta}\in\mathbb{R}^{(N_{s}+1)\times N_{y}}$ is a linear map. Conceptually, ESN is similar to a “kernel method” in the machine learning literature [37], where a nonlinear projection is used to map the data into a high dimensional space with a goal of finding a linear pattern in the feature space. It should be noted that, although this class of methods has been successfully applied to solve complex problems in practice, it lacks theoretical analysis, such as the optimality of the projection and the convergence.

In the standard ESN, “model training” is to find the linear map, $\bm{\theta}$, from the data. Let $\bm{Y}=(\bm{y}_{0},\cdots,\bm{y}_{N})$ be a time series data of length $N+1$. The linear map is computed by solving the following regularized optimization problem;

[TABLE]

in which $\|\cdot\|_{F}$ denotes the Frobenius norm, and $\beta$ is a regularization parameter. The analytical solution of (6) is

[TABLE]

The second term in (6) corresponds to the standard $L_{2}$ regularization, often called the Tikhonov regularization. In ESN, first the physical variable, $\bm{x}\in\mathbb{R}^{N_{x}}$, is projected onto a high dimensional ($\mathbb{R}^{N_{s}}$) latent space by a nonlinear map, $\Psi$, where $N_{s}\gg N_{x}$. Because $\Psi$ is generated randomly, there is no guarantee that $\bm{s}_{t}$ is linearly independent. Hence, the $L_{2}$ regularization is usually required to deal with the possible multicollinearity.

In summary, an ESN consists of the following update equations;

[TABLE]

The free parameters of an ESN are listed in Table 1.

2.2 Fixed-point ESN for Sparse Observation

Here, we consider the time series data generated by a nonlinear dynamical system,

[TABLE]

in which $\mathcal{F}(\cdot)$ is an unknown time marching operator and $\delta t$ is the sampling interval of the observation. Let $\bm{Y}^{*}=(\bm{y}^{*}_{0},\cdots,\bm{y}^{*}_{T})$ be the full observation obtained by an equidistance sampling of $\bm{y}(t)$, i.e.,

[TABLE]

We assume that the full observation, $\bm{Y}^{*}$, is not accessible, and only a partial observation of $\bm{Y}^{*}$, which is denoted by $\bm{Y}^{O}$, is known. The missing (unobserved) elements of $\bm{Y}^{O}=(\bm{y}^{O}_{0},\cdots,\bm{y}^{O}_{T})$ are marked by an arbitrarily large number, $\eta$, i.e.,

[TABLE]

Based on $\bm{Y}^{O}$, we can define an indicator set of the observation, $\bm{M}=(\bm{m}_{0},\cdots,\bm{m}_{T})$, in which

[TABLE]

Here, we are interested in reconstructing $\bm{Y}^{*}$ from a sparse observation, $\bm{Y}^{O}$, where the fraction of the missing data,

[TABLE]

is larger than 0.5. We assume that the exogenous forcing, $\bm{U}=(\bm{u}_{0},\cdots,\bm{u}_{T})$, is fully known. Figure 1 shows an example of such a sparse observation.

The standard method of training an ANN, so called the “teacher forcing” method, consists of three steps; 1) provide an input data $\bm{x}_{t-1}$, 2) make a prediction, $\hat{\bm{y}}_{t}$, and 3) compare the difference between $\hat{\bm{y}}_{t}$ and $\bm{y}^{*}_{t}$ at every time step, i.e., for $t=1,\cdots,T$. However, when there are missing entries in the data set, such sequential process is not possible. One of the conventional approaches to circumvent the problem is to fill the missing entries by using a statistical imputation method and to train an RNN against the imputed data set, assuming that the imputed data is close to $\bm{Y}^{*}$ and has no systemic bias. When the missing fraction is small, e.g., $\omega<0.2\sim 0.3$, the imputation-and-training method works fairly well. For a very high missing fraction, e.g., $\omega>0.8$, the standard imputation methods become increasingly unreliable.

To develop an ESN for the reconstruction of complex dynamics from a sparse observation, we exploit the capability of ESN as a universal function approximator [38, 39, 40]. We assume that, for a large enough $N_{s}$, there is an ESN, which can accurately computes the time series data, $\bm{Y}^{*}$, such that

[TABLE]

and

[TABLE]

in which $\bm{B}_{y}\in\mathbb{R}^{N_{s}\times N_{y}}$ and $\bm{B}_{u}\in\mathbb{R}^{N_{s}\times N_{u}}$ are the submatrices of $\bm{B}$, which correspond to $\bm{y}$ and $\bm{u}$, respectively. We can rewrite (10 – 11) as

[TABLE]

where

[TABLE]

Note that $\bm{S}^{*}_{l}$ is a function of $(\bm{y}^{*}_{0},\cdots,\bm{y}^{*}_{l-1},\bm{u}_{0},\cdots,\bm{u}_{l-1})$ and $\bm{\theta}^{*}$ depends on $(\bm{Y}^{*},\bm{U})$ through (6). It is clear from (12) that $\bm{Y}^{*}$ is a fixed point of the nonlinear operator, $\bm{\mathcal{E}}(\bm{\cdot},\bm{U})$.

Now, we define a “reconstructed” observation vector,

[TABLE]

in which $\circ$ is the Hadamard product, $\bm{1}\in\mathbb{R}^{N_{y}}$ is a vector of ones, and $\hat{\bm{y}}_{t}$ is a model prediction, or reconstruction. Simply, $\bm{y}^{R}_{t}$ takes the observation if available, otherwise the missing entry is filled by the output of the ESN. From (12), the goal of the nonlinear reconstruction of the dynamics is equivalent to finding $\bm{Y}^{R}=(\bm{y}^{R}_{0},\cdots,\bm{y}^{R}_{T})$ and $\bm{\theta}$, which satisfy

[TABLE]

Hence, we propose the following fixed-point iteration method,

[TABLE]

in which $0<\alpha<1$ is an under-relaxation parameter of the fixed point iteration, and the superscript $k$ is the iteration count. The under relaxation is required to guarantee that the Jacobian of $\bm{\mathcal{G}}$ exists (see 37). It is trivial to show that $\bm{Y}^{*}=\bm{Y}^{R}$ upon convergence.

The linear map at $k$-th iteration is computed by solving the following optimization problem,

[TABLE]

In other words, $\bm{\theta}^{k}$ is computed by comparing only with the available observation data. Note that the time index starts from $T_{s}$. As discussed in section 2.1, we cannot specify the correct initial condition of the latent state, $\bm{s}_{0}$. Hence, we disregard the ESN outputs in the first $T_{s}-1$ time steps to remove the effects of the initial condition. The analytical solution of (18) is

[TABLE]

Here, $\bm{\theta}_{j}$ is the $j$-th column of $\bm{\theta}$. Note that after a proper rescaling, the summations in the original ESN formulation (7) are time averaging, while the summations in (19) correspond to an ensemble averaging, or a Monte Carlo sampling. Hence, when the dynamical system is ergodic, the time series data is missing in random, and the number of observation, $\omega T$, is large enough, the linear map of the fixed-point method, $\bm{\theta}^{k}$, approaches to $\bm{\theta}^{*}$ as the fixed-point iteration converges.

The fixed point iteration procedure is outlined in Algorithm 1. Computing $\bm{Y}^{R,\,k+1}$ involves updating the ESN twice over the entire time series ($t=0\sim T$); first to compute $\bm{\theta}^{k}$ and then to update $\bm{Y}^{R,\,k+1}$. It is possible to reduce the computation cost by storing $\mathcal{S}^{k}=(\bm{S}^{k}_{1},\cdots,\bm{S}^{k}_{T})$ in the first passage and updating $\bm{Y}^{R}$ by computing only the projection ${\bm{\theta}^{k}}^{T}\mathcal{S}^{k}$. However, this approach requires a huge memory space. For example, when $T\sim O(10^{4})$ and $N_{s}\sim O(10^{3})$, storing $\mathcal{S}^{k}$ with double precision requires a few gigabytes of memory.

One of the advantages of a fixed point iteration method is that it does not require derivative information of the nonlinear operator, $\mathcal{G}$. Computing the derivatives requires an additional storage space of $O(N_{s}\times T)$ and has computational complexity of $O(N_{s}^{2}\times T^{2})$. However, it is challenging to analyze the convergence of a fixed point iteration and it usually exhibits a linear convergence. A theoretical analysis of the proposed fixed point iteration is provided in A.

3 Numerical experiments

In this section, a series of numerical experiments on the reconstruction of nonlinear dynamical systems is presented. The same ESN parameters are used across all the numerical experiments, except for $\beta$;

[TABLE]

The size of the ESN is also fixed at $N_{s}=10^{3}$. To generate the sparse matrix $\bm{A}$, fist, $\nu\times N_{s}^{2}=10^{4}$ elements of $\bm{A}$ are randomly selected and filled by sampling from the uniform distribution, $\mathcal{U}(-\xi_{A},\xi_{A})$. Then, $\bm{A}$ is scaled by multiplying a constant to make $\rho(\bm{A})=\rho_{max}$. As shown in [8], the accuracy of an ESN depends on the choice of those free parameters. However, we do not make an attempt to fine tune those parameters, because, without having a priori knowledge on the dynamical system, hand-tuning those parameters with a sparse observation is impractical.

3.1 Mackey–Glass Time Series

First, the behavior of the proposed fixed-point ESN is investigated by using a sparse observation of a univariate time series. The ground truth is generated by a discrete sampling of the Mackey-Glass delay-time dynamical system with a sampling interval of $\delta t=1$. The Mackey-Glass equation is [41],

[TABLE]

where we use $\gamma_{1}=0.2$, $\gamma_{2}=10$, $\gamma_{3}=0.1$, and $\tau=17$. For this set of parameters, the Mackey-Glass time series becomes chaotic [42]. The length of the time series is $T=5\times 10^{4}$. The sparse observation is generated by randomly removing $\omega$ fraction of the data. Figure 2 (a) shows the sparse observation when 95% of the time series data is removed. In general, it is very challenging to reconstruct the missing information of a univariate time series [25]. Without a prior knowledge, the usual practice is to use an interpolation. However, as shown in figure 2 (b), a standard cubic spline interpolation results in a poor estimation.

Figure 3 shows the improvement of $\bm{Y}^{R,\,k}$ over the iterations for $\omega=0.95$. The initial condition, $\bm{Y}^{R,\,0}$, is given as a linear interpolation (fig. 3 a). It is clearly shown that, over the iterations, the fixed-point ESN learns the complex nonlinear, time-delay dynamics closer and closer to the ground truth. By the 200-th iteration, the reconstructed trajectory, $\bm{Y}^{R,200}$, becomes almost indistinguishable from $\bm{Y}^{*}$.

Figure 4 shows the reconstructed chaotic attractor of the Mackey-Glass time series. Note that only the linearly interpolated time series (fig. 4 a) together with the locations of missing data is provided to the fixed-point ESN. After 200 iterations, the fixed-point ESN is able to closely reconstruct the chaotic attractor. This result suggests that the capability of ESN is beyond a high-dimensional interpolation in the phase space.

In figure 5, a normalized root mean-square error of $\bm{Y}^{R}$ is shown as a function of the iteration number. The normalized root mean-square error (NRMSE) is defined as,

[TABLE]

There is a rapid decrease in NRMSE in the first 100 iterations. NRMSE is decreased from 0.93 of the initial condition to 0.15 at $k=100$. After $k=280$, $\bm{Y}^{R}$ seems to be converged and there is no further improvement in NRMSE.

In figure 6, the behavior of the fixed-point ESN is investigated for a range of the missing rate, $0.6\leq\omega\leq 0.95$. The relaxation and regularization parameters are listed in table 2. It is found that, when $\omega$ is large, both $\alpha$ and $\beta$ need to be increased to make the fixed-point iteration stable. The data sets have the equal length, $T=5\times 10^{4}$. The number of observations changes from $2\times 10^{4}$ at $\omega=0.6$ to $2.5\times 10^{3}$ at $\omega=0.95$.

Figure 6 (a) shows the convergence of the fixed-point ESN for the range of $\omega$. During the fixed point iteration, the convergence is checked by the $l_{2}$-improvement,

[TABLE]

It is shown that, when $\omega$ is modest, the fixed point iteration converges very rapidly. For $\omega=0.6$, i.e., 40% of the data is available, the fixed-point ESN converges under 5 iterations, which increases to 50 iterations for $\omega=0.90$.

In figure 6 (b), the accuracy of the fixed-point ESN is compared with linear and cubic spline interpolations. The normalized RMSE is computed as,

[TABLE]

Here, $y^{lin}$ and $y^{csp}$, respectively, denote the linear and cubic spline interpolations. The fixed-point ESN outperforms both the linear and cubic spline interpolations even at low missing rate, $\omega=0.6$. The root mean-square error of the fixed-point ESN is about 30% of the linear interpolation and 80% of the cubic spline interpolation. The strength of the fixed-point ESN is more pronounced at higher $\omega$. When $\omega=0.95$, RMSE of the fixed-point ESN is less than 8% of the linear and cubic spline interpolations.

The effects of the free parameters and the random generation of ESN are shown in figure 7. In the experiments, four free parameters considered, $\lambda$, $\xi_{B}$, $\nu$, and $N_{s}$. In each experiment, all of the free parameters are fixed at the reference values (20) except for the one being tested. The normalized RMSEs with respect to the linear interpolation ($\sigma_{lin}$) are computed from 20 randomly generated fixed-point ESNs. It is shown that the expectation of $\sigma_{lin}$, for the reference parameter set, e.g., $\lambda=0.6$ in figure 7 (a), is much lower than $\sigma_{lin}$ from a single realization in figure 6 (b); $E[\sigma_{lin}]\simeq 0.05$ versus $\sigma_{lin}\simeq 0.10$. Note that, for the range of parameters tested, the variations of $E[\sigma_{lin}]$ is less than 0.04, indicating that the fixed-point ESN is not so sensitive to the changes in the free parameters.

The effects of the temporal relaxation parameter, $\lambda$, are shown in figure 7 (a). In general, the error decreases at smaller $\lambda$ for the Mackey-Glass problem. However, it is found that, when $\lambda\leq 0.2$, the fixed-point ESN does not converge. The temporal relaxation parameter, $0\leq\lambda<1$, works effectively as a regularization. Increasing $\lambda$ enforces a longer autocorrelation, or higher inertia, which results in a smoother trajectory. On the other hand, $\lambda=0$ does not impose such constraint in the temporal dynamics. When the missing fraction, $\omega$, is high, a regularization in the temporal structure is required to guarantee temporal smoothness in the solution, which explains why the fixed-point ESN does not converge when $\lambda\leq 0.2$.

Figure 7 (b) shows the changes in $\sigma_{lin}$ as a function of the scale of the weight matrix, $B_{ij}\sim\mathcal{U}(-\xi_{B},\xi_{B})$. It is shown that, while $E[\sigma_{lin}]$ is not so sensitive to $\xi_{B}$, the standard deviation of $\sigma_{lin}$ is an increasing function of $\xi_{B}$. The scale of $\bm{B}$ determines the relative contribution between the input variable, $\bm{x}_{t}$, and the latent state, $\bm{s}_{t}$, to the time evolution of the latent state, $\bm{s}_{t+1}$. Because the time evolution of $\bm{s}$ becomes more autonomous when $\xi_{B}$ is small, the standard deviation of $\sigma_{lin}$ also becomes smaller [43].

It is argued that making the connection of the latent state sparse leads to a richer dynamics [32]. In figure 7 (c), $\sigma_{lin}$ is computed as a function of the sparsity, $\nu$. It is shown that, in the range of $\nu$ explored, the fixed-point ESN is not sensitive to the changes in $\nu$.

In figure 7 (d), $E[\sigma_{lin}]$ is shown to be an increasing function of $N_{s}$ when $N_{s}>1000$. In general, it is expected that increasing $N_{s}$ makes the nonlinear representation capability of an ESN stronger. However, when the number of data is fixed, using a large $N_{s}$ makes it difficult to train an ESN because the number of unknown parameters increases linearly proportional to $N_{s}$. For $\omega=0.9$, the number of observations is $(1-\omega)T=5\times 10^{3}$, while the number of unknown parameters is $\text{dim}(\bm{\theta})=N_{s}+1$. The increase in the reconstruction error for $N_{s}>1000$ seems to be related with the decrease in the ratio between the number of unknown parameters to the number of data.

3.2 Lorenz-63 System

In this section, we use the Lorenz-63 system [44] to validate the fixed-point ESN. The Lorenz-63 system is given by the following set of equations,

[TABLE]

We used the coefficients from [44], $\gamma_{1}=10$, $\gamma_{2}=8/3$, and $\gamma_{3}=28$. The time series data is generated by sampling at every, $\delta t=0.02$. The length of the time series is $T=5\times 10^{4}$.

First, we consider a partial observation of the Lorenz-63 system [8]. The data set consists of only ($x_{t},z_{t}$). We assume that $z_{t}$ is completely known, while $\omega$ fraction of the data is randomly missing in $x_{t}$ . Figure 8 shows an example of the partial observation with missing data. Here, we aim to recover the missing data in $x_{t}$ from ($x^{O}_{t},z_{t}$).

Figure 9 shows the reconstruction, $\bm{Y}^{R}$, at four different iteration counts. The fixed-point ESN starts from a linear interpolation (Fig. 9 a). Similar to the Mackey-Glass time series (Fig. 5), the reconstruction error quickly reduces in the first 100 iterations (Fig. 9 b-c). By $k=100$, $\bm{Y}^{R}$ is already very close to the ground truth, $\bm{Y}^{*}$. At $k=300$, $\bm{Y}^{R}$ becomes almost indistinguishable from $\bm{Y}^{*}$. The RMSE normalized by the linear interpolation is $\sigma_{lin}=0.07$ for the time window shown in figure 9 (d).

The normalized RMSEs ($\sigma_{lin}$) for $\omega=0.9$ and 0.95 are shown in table 3, together with the parameters ($\alpha$, $\beta$). $\sigma_{lin}$ is computed over the entire trajectory, $T=5\times 10^{4}$. It is shown that, when $\omega=0.9$, RMSE of the fixed-point ESN is only about 6% of the linear interpolation, which changes to 37% at $\omega=0.95$. Figure 10 explains why there is a significant increase in RMSE at $\omega=0.95$. Because of the random sampling, in the time window shown in figure 10, there is only one observation at $t=193\delta t$ for $t\in[96\delta t,285\delta t]$. In other words, two extremely large missing intervals ($97\delta t$ and $92\delta t$) exist consecutively. This is a rare case, because the average interval of the missing observation is $20\delta t$ for $\omega=0.95$. As discussed in [8], it is not possible to correctly predict $x(t)$ given only $z(t)$, because, when $y(t)$ is not observed, the symmetry of the Lorenz-63 system implies that $x(t)$ and $-x(t)$ are equally possible. Hence, the reconstructed trajectory, $\bm{Y}^{R}$, predicts a positive peak at $t/\delta t=165$, while the ground truth, $\bm{Y}^{*}$, reaches a negative peak at $t/\delta t=180$. Then, for $193\leq t/\delta t\leq 285$, the fixed-point ESN cannot adjust its dynamics due to the lack of the data. On one hand, this observation demonstrates the limitation of the fixed-point ESN, which fails to capture the dynamics when the missing interval is too long. One the other hand, the Lorenz system is chaotic, of which long-term predictability is quite limited.

In the next experiment, the sparse observation of all three variables, $(x_{t},y_{t},z_{t})$, are given, and we aim to reconstruct the dynamics of the three-dimensional time series simultaneously. Figure 11 (a–c) shows an example of the sparse observation of the Lorenz-63 time series. All three variables have the same missing rate, i.e., $\omega_{x}=\omega_{y}=\omega_{z}=\omega$. Figure 11 (d) shows the availability of the data. Since the time series data is removed at random, the probability that all three variables are observed at the same time frame is 0.001 for $\omega=0.9$.

Figure 12 shows the reconstructed Lorenz-63 time series from the 90% randomly missing data after 90 iterations. The relaxation and regularization parameters are, respectively, $\alpha=0.2$ and $\beta=10^{-6}$. It is shown that the fixed-point ESN is capable of accurately reconstruct the nonlinear dynamics of the Lorenz-63 system.

The normalized RMSEs with respect to $\bm{Y}^{*}$ are shown in table 4 for $\omega=0.90$ and 0.95. The normalized RMSE is defined as,

[TABLE]

For $\omega=0.90$, the fixed-point ESN makes a much more accurate reconstruction compared to the linear and cubic spline interpolations. However, the accuracy of the fixed-point ESN significantly decreases at $\omega=0.95$. The fixed-point ESN still outperforms those interpolation methods. But, there is about three-fold increase in $\sigma_{lin}$ as $\omega$ changes from 0.9 to 0.95. At $\omega=0.95$, the cubic spline interpolation results in a large error due to the overshoot, which makes $\sigma_{csp}$ smaller than $\sigma_{lin}$.

3.3 Lorenz-96 System

In this experiment, the fixed-point ESN is tested against a 6-node Lorenz-96 system [45];

[TABLE]

The Lorenz-96 system is periodic, i.e., $y^{(i-N_{y})}=y^{(i+N_{y})}=y^{(i)}$. Here, we consider $N_{y}=6$ and $F=8$. The time series data is generated by sampling at every $\delta t=0.02$ and the length of the time series data is $T=5\times 10^{4}$.

Figure 13 (a–c) show a subset of the 6-node Lorenz-96 system for $\omega=0.9$. Figure 13 (d) shows an example of the data availability. Because the probability that every $y^{(i)}$ is observed at the same time frame is $10^{-6}$ for $\omega=0.9$, the data set does not contain a time frame at which every $y^{(i)}$ is observed.

Figure 14 shows the reconstructed 6-node Lorenz-96 system from the sparse observation of $\omega=0.9$. The parameters of the fixed-point ESN are $\alpha=0.2$ and $\beta=10^{-7}$. It is again shown that $\bm{Y}^{R}$ is almost indistinguishable from $\bm{Y}^{*}$. The normalized RMSEs are $\sigma_{lin}=0.11$ and $\sigma_{csp}=0.10$. When the fixed-point ESN is used to recover the dynamics from a sparser data set of $\omega=0.95$, the normalized RMSEs are significantly increased to $\sigma_{lin}=0.57$ and $\sigma_{csp}=0.46$.

3.4 Forced van der Pol oscillator

Now, we consider a forced van der Pol oscillator, which is given by the following equations,

[TABLE]

The exogenous forcing, $u(t)$, is given by an Ornstein-Uhlenbeck process as

[TABLE]

in which $W$ is the Wiener process. The parameters used in this simulation are, $\gamma_{1}=2$, $\gamma_{2}=0.2$, and $\gamma_{3}=5\sqrt{2\gamma_{2}}$. Similar to the partial observation of the Lorenz-63 system, $y_{2}$ is not given and the data set consists of $(y_{1}(t),u(t))$, which is sampled at every $\delta t=0.5$. Hereafter, we use $y_{t}$ to denote $y_{1}(t)$. The length of the time series is $T=5\times 10^{4}$.

In sections 3.1 – 3.3, we consider the reconstruction of chaotic time series. As shown in figure 10, the accuracy of the reconstruction becomes lower as the interval of the missing information increases. The decrease in the accuracy seems to be related with the chaotic nature of the time series. On the other hand, the unforced van der Pol oscillator has a stable limit cycle. The dynamics of the forced van der Pol oscillator in (27–28) depends on the restoring force to the limit cycle and excursion due to the exogenous forcing, $u_{t}$. Hence, in theory, if we know the governing equations and the full history of $u_{t}$, we can find $\bm{Y}^{*}$ no matter how long the interval between the observations.

The sparse observation of the van der Pol oscillator is shown in figure 15. Here, we consider an extremely sparse data, in which 99% of the observations are randomly removed. Since $T=5\times 10^{4}$, the number of observations is only 500. Although the average interval of the missing data is $100\delta t$, due to the random removal, it is not uncommon to find a missing interval of $200\delta t\sim 300\delta t$. The exogenous forcing is shown in figure 15 (b). While the Ornstein-Uhlenbeck process is computed with $dt=2.5\times 10^{-3}$ to numerically integrate (27–28), $u_{t}$ is downsampled with the sampling interval of $\delta t=0.5$. In a sense, we introduce an uncertainty, in which the sub-timescale information of $u_{t}$ is unknown. Hence, the time series becomes a Markov process, $y_{t+1}\sim p(y_{t+1}|y_{t},u_{t})$. Although the subscale fluctuations of $u_{t}$ are not significant, it breaks the assumption of the fixed-point ESN that there exists a deterministic mapping, $y_{t+1}=f(y_{t},u_{t})$, which guarantees the existence of a fixed point.

Figure 16 shows the $l_{2}$-improvement in terms of the iteration count. In both cases, $\alpha=0.2$ and $\beta=10^{-6}$ are used. It is shown that, for $\omega=0.98$, the fixed-point ESN quickly converges to a solution. The $l_{2}$-improvement becomes $\simeq 10^{-6}$ by the 100-th iteration. When only 1% of the data is available ($\omega=0.99$), the convergence becomes much slower. The $l_{2}$-improvement reaches $e^{k}\simeq 0.008$ at the iteration number $k=160$, and there is no further reduction of $e^{k}$ afterward.

Figure 17 shows the reconstructed dynamics for $\omega=0.99$. Even with the extremely sparse data, the fixed-point ESN is able to reconstruct the nonlinear dynamics. The normalized RMSE is $\sigma_{lin}=0.25$ for $\omega=0.99$ and $\sigma_{lin}=0.21$ for $\omega=0.98$. It is observed that, although the reconstructed trajectory, $\bm{Y}^{R}$, closely follows the dynamics, $\bm{Y}^{R}$ is much noisier compared to the previous experiments on the deterministic time series, which seems to be related with the uncertainty in $u_{t}$.

4 Concluding remarks

In this study, a novel model-free method is developed to reconstruct nonlinear dynamics from a temporally sparse observation, in which a large fraction of the time series data is randomly missing, or unobserved. Based on the assumptions of the noise-free observation and universal function approximation capability of a recurrent neural network, we show that the reconstruction problem can be solved by a fixed-point problem. The fixed-point method consists of two major components; outer-iteration to find the fixed-point solution and inner-loop to compute the parameters of the recurrent neural network. In theory, any recurrent neural network architecture can be used for the inner-loop, but the computational cost of solving the inner-loop with a stochastic gradient method will make it impractical. Hence, we employ the echo state network for the recurrent neural network, of which parameter can be computed by solving a simple ordinary least square problem. The fixed-point method is simple to implement and can be solved with $O(N_{s}^{2}\times T)$ floating point operations, because it does not require to compute the derivatives of the nonlinear operator. Although it may be possible to develop a gradient-based optimization algorithm to solve the problem, computing the derivatives of the echo state network will require additional memory of $O(N_{s}\times T)$ and floating operation of $O(N_{s}^{2}\times T^{2})$, which makes it difficult to use in practice.

The proposed method, referred as the fixed-point ESN, is tested against time series data generated by chaotic dynamical systems. For the Mackey-Glass time series, it is shown that the fixed-point ESN can accurately reconstruct the nonlinear dynamics even when the 95% of the data are randomly removed. For the missing fraction $\omega=0.95$, the average time interval between two consecutive observations is $\Delta T=20$, while the characteristic period of the Mackey-Glass time series is $T_{c}\simeq 50$ [46]. The Lorenz-63 and 6-node Lorenz-96 systems are used to demonstrate the capability of the fixed-point ESN for the reconstruction of nonlinear dynamics of multivariate time series. In both cases, the fixed-point ESN provides very good approximation of the complex nonlinear dynamics up to $\omega=0.9$. Finally, we use a forced van der Pol oscillator, of which dynamics is largely determined by the exogenous forcing, to investigate the behavior of the fixed-point ESN. It is shown that the fixed-point ESN can learn the relation between the exogenous forcing and the dynamical system with only 1% of the time series data.

It is demonstrated that the echo state network is capable of learning the dynamics from only a partial observation of a complex nonlinear dynamical system. In this study, we limit our interest only to relatively low dimensional dynamical systems. It is a subject of further investigation how the nonlinear reconstruction method can be extended to a high dimensional dynamical system with complex correlation structures, such as a spatio-temporal process. In many practical applications, observations are typically corrupted by a sensor noise, which makes the observation a stochastic process. Although, it is shown that the fixed-point ENS still can learn the nonlinear dynamics for a small noise, the accuracy of the fixed-point ESN is not guaranteed under a small signal-to-noise ratio. Finally, although it is empirically shown that the fixed-point ESN is very effective in the nonlinear reconstruction, our understanding on the theoretical understanding, e.g., accuracy and convergence, is far from complete.

Appendix A Convergence analysis of fixed-point ESN

It is challenging to theoretically show the convergence of a fixed point iteration of a nonlinear map. One of the widely used method to analyze the convergence is based on the Banach fixed-point theorem, which provides a sufficient condition. The Banach fixed point theorem requires a contraction mapping;

[TABLE]

in which $0\leq L<1$ is a Lipschitz constant. In [9], it is assumed that ESN provides a contraction mapping. However, it is a too strong assumption and in general such a globally uniform contraction is not expected for a nonlinear map. Instead, here, we explore a local convergence around a fixed point, $\bm{Y}^{*}$.

Let $\mathcal{J}_{\mathcal{G}}(\bm{Y})$ denote the Jacobian of $\bm{\mathcal{G}}(\bm{Y},\bm{U})$ in (17). From the contraction mapping theorem, it can be shown that, if $\|\mathcal{J}_{\mathcal{G}}(\bm{Y}^{*})\|<1$ for a subordinate matrix norm $\|\cdot\|$, there is a neighborhood around the fixed point, $\bm{Y}^{*}$, in which the fixed-point iteration converges to $\bm{Y}^{*}$.

For simplicity, here we assume that $y_{t}$ is univariate, i.e., $N_{y}=1$. But, extending the analysis to a multivariate time series is straightforward. The Jacobian of the fixed point iteration is given as

[TABLE]

in which $\delta_{ij}$ is the Kronecker delta.

In the second term on the right hand side (RHS) of (31),

[TABLE]

Here, $\bm{M}_{1}$ is a set of nonzero entries of $\bm{M}$, i.e., the summation is over the available observations, and $\bm{D}=(\sum\bm{S}\bm{S}^{T}+\beta\bm{I})^{-1}$. By definition,

[TABLE]

and, from (10), $y^{O}_{l}={\bm{\theta}^{*}}^{T}\bm{S}^{*}_{l}$. Then, (32) can be written as

[TABLE]

At a fixed point, $\bm{Y}=\bm{Y}^{*}$, we have $\bm{S}_{l}=\bm{S}^{*}_{l}$ and $\bm{\theta}=\bm{\theta}^{*}$, which makes $\partial\bm{\theta}/\partial y_{j}=0$. Therefore, the second term on the RHS of (31) vanishes.

In the last term on the RHS of (31), because $\bm{S}_{i}=(1,\bm{s}_{i}^{T})^{T}$, it is sufficient to investigate the behavior of $\partial\bm{s}_{i}/\partial y_{j}$. From (1), $\bm{s}_{i}$ depends only on the past trajectory, $\bm{s}_{i}=f(y_{0},\cdots,y_{i-1},\bm{u}_{0},\cdots\bm{u}_{i-1})$. Thus,

[TABLE]

For $i=j+1$,

[TABLE]

in which $\bm{Z}_{j}$ is a diagonal matrix,

[TABLE]

Similarly, when $i>j+1$,

[TABLE]

Here,

[TABLE]

To sum up, the Jacobian, $\mathcal{J}_{\mathcal{G}}(\bm{Y}^{*})$, is a lower triangular matrix, of which element is

[TABLE]

Here, $\widetilde{\bm{\theta}}=(\theta_{2},\cdots,\theta_{N_{s}+1})$ is a submatrix of $\bm{\theta}$, which corresponds to $\bm{s}_{t}$. The $l_{1}$-norm of $\mathcal{J}_{\mathcal{G}}(\bm{Y}^{*})$ is defined as

[TABLE]

Let $J$ be the column index of the maximum absolute summation. The sufficient condition for a local convergence, $\|\mathcal{J}_{\mathcal{G}}(\bm{Y}^{*})\|_{1}<1$, yields

[TABLE]

However, evaluating (39) a priori is not possible, because of its dependence on $\bm{\theta}^{*}$ and $\bm{Z}^{*}$.

The sufficient condition (39) involves a summation over the length of the time series data, $T$, which is usually very long, e.g., $T\sim O(10^{4})$. For (39) to hold, $\bm{K}^{i-1}_{J+1}$, should be a monotonically decreasing sequence, which requires

[TABLE]

For (39) to hold, it requires that $|\bm{K}^{i-1}_{J+1}|$, which dictates the sensitivity of the latent state in the future ($t>J$) to the changes in the latent state at $t=J$, decays rapidly in time. It is challenging to evaluate $\bm{K}^{i-1}_{J+1}$ for a general case. Since the temporal decay of $\bm{K}^{i-1}_{J+1}$ depends on the diagonal matrix $\bm{Z}_{J}$, where $0\leq\text{diag}(\bm{Z}_{J})_{i}\leq 1$ for $i=1,\cdots,N_{s}$, let assume $\bm{Z}=\bm{I}$. Also, note that, since $\bm{A}$ and $\bm{B}$ are generated by independent uniform distributions,

[TABLE]

where $\zeta_{A}\sim\mathcal{U}(-\xi_{A},\xi_{A})$ and $\zeta_{B}\sim\mathcal{U}(-\xi_{B},\xi_{B})$. Then, the left hand side of (39) is

[TABLE]

Now, the sufficient condition for the convergence is simply

[TABLE]

Note that $\widetilde{\bm{\theta}^{*}}^{T}\bm{B}_{y}$ corresponds to a weighted average of a zero-mean independent uniform random variables, which behaves $\widetilde{\bm{\theta}^{*}}^{T}\bm{B}_{y}\simeq 0$ as $N_{s}\rightarrow\infty$. Although the convergence analysis is performed in a limited scope, it shows the effects of the temporal relaxation parameter, $\lambda$, on the convergence. It is consistent with the finding in figure 7 (a) that the fixed-point ESN fails to converge when $\lambda$ becomes small.