0.0.1 Introduction↩︎

0.0.1.1 Background

In recent years, artificial intelligence (AI) has become increasingly integrated into various aspects of our daily lives, from healthcare to autonomous vehicles. This widespread adoption has raised concerns regarding AI system reliability and raised the need for robust quality assurance measures. In response to these trends, multiple organizations and regulatory bodies have published guidelines for ensuring AI safety and reliability, such as the European Union’s AI Act [1] and AI HLEG guidelines [2]. Although these guidelines have contributed to the gradual establishment of AI quality assurance processes, the development of specific evaluation metrics and implementation methodologies remains in its nascent stages. The development of these metrics and implementation methodologies is a critical issue, especially for control systems incorporating AI, because such systems interact directly with the physical world, potentially impacting human lives in extreme scenarios.

In this paper, we focus on the control of time-series models, particularly long short-term memory (LSTM) networks, having real-time capability, which is an important property in control systems. Although typical time-series models can extract nonlinear features from time-series data, store them as internal states, and utilize them to make inferences at each time step, especially, LSTM surpasses other time-series models in the accuracy of prediction [3]–[5]. This ability is well-suited for predicting future outputs of nonlinear systems in model predictive control (MPC), and therefore LSTM’s capacity to track target values is better than traditional control methods [6]–[8]. Moreover, LSTM can be implemented on Field Programmable Gate Array (FPGA), which are suitable for machine embedding [9]. However, the state-space nature of time-series models also presents significant risks: If a model receives temporarily anomalous input values, it may retain inaccurate state memories, potentially leading to persistent errors in prediction that could significantly compromise the system’s safety and reliability.

The concerns arising from the use of time-series models for control systems can be mitigated by ensuring stability, robustness, and resilience in these models [10]. Stability refers to consistent performance across similar inputs, which forms the foundation for robustness and resilience. Robustness refers to patience against perturbations, while resilience refers to the ability to return to a normal state after transitioning to an abnormal mode due to anomalous inputs. See also Section 0.0.2.1 for related work on these concepts. Time-series models with these properties should exhibit consistent behavior under normal conditions. However, these properties cannot be adequately evaluated using conventional metrics focused on inference accuracy, such as precision, recall, and F1-score [11]. Therefore, it is necessary to formulate new evaluation metrics and develop techniques to achieve the required metrics. This issue might be addressed by comprehensively collecting and learning from data across various scenarios for evaluation, but there are difficulties in defining comprehensiveness and collecting data at scale. Given this context, we mathematically formulate the resilience of LSTM, which is different from data-driven approaches. For theoretical and experimental approaches in assurance of LSTM stability and quality, see also Section 0.0.2.2 and Section 0.0.2.3.

0.0.1.2 Objective

Our research focuses on formulating evaluation metrics and implementation methodologies to assess and ensure the resilience of LSTMs. Specifically, our goals are as follows:

• To introduce a novel mathematical formulation of the resilience, recovery time, defined as the duration required for a system to return to its normal state after transitioning to an abnormal mode due to anomalous inputs (see Figure 1).

• To theoretically evaluate the upper bound of recovery time and to develop a practical data-independent method for estimating it.

• To present parameter adjustment techniques for controlling the recovery time of the model during the training process.

The key tool to achieve our goals is the incremental input-to-state stability (\(\delta\)ISS) [12], [13], because \(\delta\)ISS guarantees the resilience of LSTM. Specifically, even if temporary perturbations in the input sequence disrupt inference, the state variables can converge to the same inference values that would have been obtained without perturbations. Sufficient conditions for the LSTM to exhibit the \(\delta\)ISS properties have already been derived as inequalities of parameters [14], [15]. However, it is not clarified in these studies how long it takes for state variables to converge. We will refine and improve upon these previous studies to achieve the above three goals.

0.0.1.3 Contribution

The main contributions of this study are as follows:

0.0.1.3.0.1 Theoretical Contributions

• Refinement of an evaluation of LSTM’s invariant sets, which form the foundation for stability analysis (Section 0.0.4).

• Improvement of sufficient conditions for LSTM to be \(\delta\)ISS (Section 0.0.4).

• Definition of recovery time that quantifies the resilience of LSTM and derivation of a data-independent upper bound for it (Section 0.0.5).

0.0.1.3.0.2 Practical Contributions

• Introduction of a quantitative metric for quality assurance through the recovery time (Section 0.0.6).

• Development of a method to control the trade-off between recovery time and inference accuracy for LSTMs (Section 0.0.6).

0.0.2 Related Work↩︎

0.0.2.1 Control System and Stability Theory

In control theory, stability and resilience are important properties. Furthermore, it is well-established that the stability of controlled objects fundamentally determines the overall stability and resilience of control systems. For instance, systems possessing input-to-state stability (ISS) properties achieve global asymptotic stability under closed-loop control [16]. Moreover, systems with \(\delta\)ISS characteristics acquire robustness when implemented with model predictive control (MPC) [13].

Therefore, in control methods such as MPC where controlled objects are modeled as dynamical systems, stability requirements for the dynamical systems are essential. While stability can be comprehensively guaranteed by general theory for linear systems, practical applications require specific and detailed engineering for a nonlinear system. This research focuses on LSTM networks as a specific case of nonlinear systems.

0.0.2.2 Theoretical LSTM Stability Assurance

Theoretical approaches have investigated sufficient conditions for LSTM networks to satisfy various stability properties. For single-layer LSTMs, conditions for global asymptotic stability [17], ISS [18], and \(\delta\)ISS [14] have been successively established. Research on sufficient conditions for \(\delta\)ISS has been extended to multilayer LSTM architectures [15].

However, existing studies face a significant limitation: while they can guarantee recovery, they cannot provide estimates of time to recover. Our research addresses these challenges by advancing the stability theory of LSTM to enable both quantitative evaluation and adjustment of LSTM resilience.

0.0.2.3 Experimental LSTM Quality Assurance

This paper discusses theoretical LSTM quality assurance; however, there are also some experimental methods. For example,  [19] evaluated the robustness of LSTM inference using noise-added datasets.  [20] proposed creating robust LSTMs by expanding training datasets to include abnormal scenes. [21] proposed a data-driven stability assurance method.

These data-driven verification methods, like the ones mentioned above, are comprehensive when sufficient data are available, yet face challenges in determining adequate data quantities and provide only sparse guarantees for nonlinear models like LSTMs. Theoretical verification offers limited but exhaustive guarantees within its scope, and is therefore complementary to, rather than in conflict with, a data-driven approach.

0.0.3 Preliminary↩︎

In this section, we introduce essential mathematical preliminaries, an important stability concept: \(\delta\)ISS, and the LSTM.

0.0.3.1 Notation

For a vector \(v \in \mathbb{R}^n\), \(\|v\|=\|v\|_2\) represents the Euclidean norm of \(v\), and \(\|v\|_{\infty}=\max_{i \in [n]} |v_{(i)}|\) represents the maximum norm, where \(v_{(i)}\) is the \(i\)-th component of \(v\) and \([n]=\{1,2,\ldots,n\}\). For a discrete-time sequence of vector \(v(k), k \in \mathbb{Z}_{\geq 0}\), we denote \(v(k_1 : k_2)=[v(k_1),\ldots, v(k_2)]\), and \(\|v(k_1:k_2)\|_{2,\infty}=\|(\|v(k_1)\|_2,\ldots,\|v(k_2)\|_2)^\top\|_\infty\).

For a matrix \(A \in \mathbb{R}^{m \times n}\), \(\|A\|=\|A\|_2\) represents the induced 2-norm of \(A\), and \(\|A\|_{\infty}\) does the induced \(\infty\)-norm, that is, \(\|A\|_{\infty}=\max_i \sum_j A_{(i,j)}\). In addition, \(\left\lvert{A}\right\rvert\) denotes a matrix whose components are the absolute values of \(A\), that is, \(\left\lvert{A}\right\rvert \mathrel{\vcenter{:}}= (|A_{(i,j)}|)\). For a matrix \(A \in \mathbb{R}^{n \times n}\), \(\rho(A)\) denotes the spectral radius of \(A\), that is, \(\rho(A) \mathrel{\vcenter{:}}= \max_{i \in [n]} |\lambda_i|\), where \(\lambda_{i}\) are the eigenvalues of \(A\).

We denote the sigmoid function as \(\sigma(w)=\frac{1}{1+e^{-w}}\) and the hyperbolic tangent as \(\phi(w)=\tanh{w}=\frac{e^w-e^{-w}}{e^w+e^{-w}}\).

0.0.3.2 Stability

Now we consider a (autonomous) discrete-time nonlinear system \[\label{NS} s(t+1)=f(s(t),x(t))\tag{1}\] with \(f(0,0)=0\), for initial state \(s(0)\) and input \(x(t), t \geq 0\). In the following, we denote \(s_i(t)\) as the system whose initial state is \(s_i(0)\) and input is \(x_i(t)\).

0.0.3.3 LSTM Network

The LSTM is described by the following relations: \[\begin{align} c^{(l)}(t+1) &=\sigma(W^{(l)}_f x^{(l)}(t) + U^{(l)}_f h^{(l)}(t) + b^{(l)}_f) \odot c^{(l)}(t) \notag\\ &\quad + \sigma(W^{(l)}_i x^{(l)}(t) + U^{(l)}_i h^{(l)}(t) + b^{(l)}_i) \notag\\ &\qquad \odot \phi(W^{(l)}_c x^{(l)}(t) + U^{(l)}_c h^{(l)}(t) + b^{(l)}_c), \\\tag{3}\\ h^{(l)}(t+1) &= \sigma(W^{(l)}_o x^{(l)}(t) + U^{(l)}_o h^{(l)}(t) + b^{(l)}_o)\notag \\ &\quad \odot \phi(c^{(l)}(t+1)), \tag{4}\\ y(t)&= U_y h^{(L)}(t) + b_y, \tag{5}\noeqref{y94L-LSTM}\\ x^{(l)}(t)&= \begin{cases} x & \text{for } l=1,\\h^{(l-1)}(t+1) & \text{for } l=2,\ldots L,\end{cases} \tag{6} \end{align}\] where \(x \in [-x_\mathrm {max},x_\mathrm {max}]^{n_x}\), \(c^{(l)} \in \mathbb{R}^{n_c}\), \(h^{(l)} \in (-1,1)^{n_c}\), \(y \in \mathbb{R}^{n_y}\) are the input, the cell state, the hidden state, and the output, respectively, for each time \(t \in \mathbb{Z}_{\geq 0}\), each layer \(l \in [L]\), and \(x_\mathrm {max}>0\). Moreover, \(W^{(l)}_{\ast} \in \mathbb{R}^{n_c \times n_x}, U^{(l)}_{\ast} \in \mathbb{R}^{n_c \times n_c}, U_y \in \mathbb{R}^{n_y \times n_c}\) are the weight matrices, \(b^{(l)}_{\ast} \in \mathbb{R}^{n_c}\) and \(b_y \in \mathbb{R}^{n_y}\) are the biases, where \(\ast\) denotes \(f, i, c\), or \(o\).

In this paper, we regard the LSTM as a discrete-time nonlinear system 1 . We define the state \(s(t)\) and the input \(x(t)\) of 1 as \(s^{(1:L)}(t)=[s^{(1)}(t),\dots, s^{(L)}(t)]\) and \((x^{(l)}(t)^\top,b_c^{{(l)}\top})^\top\), respectively, where \(s^{(l)}(t)=(c^{(l)}(t)^\top, h^{(l)}(t)^\top)^\top\).

0.0.4 \(\delta\)ISS of LSTM↩︎

In this section, we relax the existing \(\delta\)ISS condition of the LSTM [14] by refining the invariant set of the LSTM dynamics. This relaxation contributes to an accurate evaluation of the recovery time introduced in Section 0.0.5.

0.0.4.1 Invariant Set in LSTM

Define \(G^{(l)}_*, G^{(l)}_c:\mathbb{R}_{\ge 0}\to \mathbb{R}_{\ge 0}\) (\(*\) denotes \(f, i\), or \(o\)) by \[\begin{align} G^{(l)}_*(\eta)&= \left\|{\left(x_\mathrm {max}|W^{(l)}_*| \boldsymbol{1}_{n_x}+ \eta |U^{(l)}_*| \boldsymbol{1}_{n_c} + b^{(l)}_*\right)_+}\right\|_{\infty}, \\ G^{(l)}_c(\eta)&=\left\|{x_\mathrm {max}|W^{(l)}_c| \boldsymbol{1}_{n_x}+ \eta |U^{(l)}_c| \boldsymbol{1}_{n_c} + |b^{(l)}_c|}\right\|_{\infty}, \end{align}\] where \(\boldsymbol{1}_n\) is a \(n\)-dimensional vector of all components being one, and \(v_+\) denotes a vector \((\max\{v_{(i)},0\})_{i \in [n_c]}\) for a vector \(v\). Note that \(x_\mathrm {max}=1\) if \(n\ge 2\). In addition, we define the sequences \(\{\overline{\sigma}^{(l)}_*(k)\}_{k=0}^{\infty}\), \(\{\overline{\phi}^{(l)}_c(k)\}_{k=0}^{\infty}\), \(\{\eta^{(l)}(k)\}_{k=-1}^{\infty}\), and \(\{\overline{c}^{(l)}(k)\}_{k=0}^{\infty}\) using following recursive formulas: \[\begin{align} \overline{\sigma}^{(l)}_*(k)&=\sigma(G^{(l)}_*(\eta^{(l)}(k-1))), \tag{7}\\ \overline{\phi}^{(l)}_c(k)&=\phi(G^{(l)}_c(\eta^{(l)}(k-1))),\tag{8}\\ \eta^{(l)}(k)&=\phi(\overline{c}^{(l)}(k)) \overline{\sigma}^{(l)}_o(k),\tag{9} \\ \overline{c}^{(l)}(k) &=\frac{\overline{\sigma}^{(l)}_i(k) \overline{\phi}^{(l)}_c(k)}{1-\overline{\sigma}^{(l)}_f(k)}, \end{align}\] for \(k \in \mathbb{Z}_{\geq 0}\), with initial term \(\eta^{(l)}(-1)=1\), where \(\ast\) denotes \(f,i\) or \(o\). Moreover, we define that \[\begin{align} \mathcal{C}^{(l)}(k) &\mathrel{\vcenter{:}}= \left\{c^{(l)}\in \mathbb{R}^{n_c}:\|c^{(l)}\|_{\infty} \leq \overline{c}^{(l)}(k) \right\},\\ \mathcal{H}^{(l)}(k) &\mathrel{\vcenter{:}}= \left\{h^{(l)}\in \mathbb{R}^{n_c}:\|h^{(l)}\|_{\infty} \leq \phi(\overline{c}^{(l)}(k))\overline{\sigma}^{(l)}_o(k) \right\}. \end{align}\]

The exact proof is described in Appendix 0.0.10.1. We note that the invariant set \(\mathcal{S}^{(l)}(0)\) is considered in the literature [14], [15].

0.0.4.2 \(\delta\)ISS of LSTM

Using Proposition 1, we provide the following sufficient condition of \(\delta\)ISS:

The rigorous proof is described in Appendix 0.0.10.3, and where we also see that a similar relaxation can be obtained for ISS.

0.0.5 Recovery Time↩︎

In this section, we define a novel mathematical formulation, the recovery time of LSTM, which accurately expresses the time required for the system to recover from disturbances and return to its normal state. This definition makes it difficult to evaluate in practical applications due to data dependence, but we overcome this limitation by deriving a practical and computable upper bound for the recovery time. This upper bound constitutes the primary novel contribution of our research, offering a fresh perspective on LSTM performance evaluation.

0.0.5.1 Definition

Evaluating the recovery time \(T_R\) is difficult because it depends on the input data, then we utilize the function \(\beta\) of \(\delta\)ISS in Definition 2.

0.0.5.2 Upper Bound of Recovery Time

Theorem 2 indicates that there exists a \(\mathcal{KL}\) function \(\beta\) satisfying \[\begin{align} &\| s(t, s_1(0), x(0:t))- s(t, s_2(0), x(0:t)) \| \\ &\le \beta(\|s_1(0) - s_2(0) \|,t) \label{beta} \end{align}\tag{10}\] for any initial states \(s_1(0), s_2(0)\in\prod^L_{l=1}\mathcal{S}^{(l)}(k)\) and input sequence \(x(0:t) \in [-x_\mathrm {max},x_\mathrm {max}]^{n_x \times (t+1)}\) of the LSTM and \(t>0\). The proof of Theorem 2 also leads to a more precise and explicit evaluation of 10 :

Detailed proof is stated in Appendix 10.

0.0.6 Proposal↩︎

We provide a comprehensive strategy for translating theoretical foundations into practical applications.

0.0.6.1 Estimate

0.0.6.1.0.1 Aim

Quantitatively assess the resilience of LSTMs.

0.0.6.1.0.2 Method

Evaluate the obtained model’s resilience using the following index: \[\begin{align} &\bar{T}_{R}(e;k) \\ &=\min \{t_1 \ge 0 \;: \;\tilde{\beta} (2\bar{s}(k),t; k) \le e/\|U_y\| \;\text{for any} \;t \ge t_1\} \end{align}\] for some \(k\). Since \(\bar{T}_{R}\) depends solely on the model’s weight parameters, it can be evaluated independently of the data. According to Theorem 5, a smaller \(\bar{T}_{R}\) indicates a stronger resilience capacity of the model.

0.0.6.1.0.3 Use case

We propose to establish this metric as a formal requirement for assessing LSTM quality, contextualized within the framework of control periods.

0.0.6.2 Training

0.0.6.2.0.1 Aim

Build an LSTM with a specified degree of resilience.

0.0.6.2.0.2 Method

Control the model’s resilience (i.e., \(T_R\)) during training by the following loss function: \[\begin{align} Loss &= \mathcal{L} + \lambda \Phi(\theta,\varepsilon), \tag{11}\\ \Phi(\varepsilon) &= \sum^{L}_{l=1} \max \{\rho(A^{(l)}_s(k))-1+\varepsilon, 0 \} \tag{12}, \end{align}\] where \(\mathcal{L}\) denotes the existing loss function, \(\Phi\) denotes the penalty term for the resilience, \(\lambda>0\) denotes the penalty intensity, and \(\varepsilon \in [0, 0.5]\) denotes the control parameter for the resilience. The resilience of LSTM can be ensured through the design of the loss function, without modifications to the model architecture. In addition, the degree of resilience can be precisely calibrated by adjusting the parameter \(\varepsilon\). In fact, as training progresses, \(\Phi(\theta,\varepsilon)\) is expected to approach 0, and after training is completed, the model is expected to satisfy \(\rho(A^{(l)}_s(k)) \thickapprox 1-\varepsilon\). Based on Theorems 2 and 5, \(T_R\) is expected to decrease as \(\rho(A^{(l)}_s(k))\) decreases. Consequently, increasing the value of \(\varepsilon\) enhances the resilience of the model. This amplification of resilience concurrently imposes stronger constraints on the model, potentially leading to a trade-off between resilience and model accuracy.

0.0.6.2.0.3 Use case

We tune \(\varepsilon\) to balance resilience and accuracy requirements, observing the trade-off and selecting an optimal model.

0.0.7 Experiments↩︎

In this section, we verify the following two properties of our proposed method:

• Increasing \(k\) in \(\bar{T}_R(k)\) leads to more accurate recovery time estimates.

• Adjusting \(\varepsilon\) enables the tuning of \(T_R\), and there is a trade-off between recovery time and inference loss.

Detailed experimental settings are given in Appendix 0.0.11.

We also compare our resilience training method and adversarial training in Appendix 0.0.12.

0.0.7.1 Experiment I: Simplified Model

In this section, we employ a simplified model to validate (A).

0.0.7.1.1 Experimental Setup

We employed a simplified model using an LSTM network with randomly initialized weight parameters. We applied step input to this model, introduced instantaneous noise, and observed the time required for inference recovery (for a visual representation of the recovery time measurement, see Figure 1).

0.0.7.1.2 Results and Analysis

We calculated \(T_R\) and \(\bar{T}_R(k)\) \((k=0,1,\cdots,20)\) for randomly initialized \(20\) models, deriving estimated average test errors and correlation coefficients. As shown in Figure 2, the estimation test error consistently decreases as \(k\) increases, eventually converging at around \(k=20\) steps. Furthermore, the correlation coefficient exhibits a generally increasing (except for a temporary decrease) trend as \(k\) increases, which indicates a strong relationship between the true and our estimated values.

Figure 3 compares the estimation accuracy for \(k=0\) and \(k=20\), demonstrating improved accuracy with a larger \(k\). The left side of the red line is the area where \(\delta\)ISS is guaranteed by Theorem 2. Increasing \(k\) from \(0\) to \(20\) decreases \(\rho(A(k))\) by 32% on average, then all data points are included in the \(\delta\)ISS guaranteed area.

To illustrate this improvement in a specific case, we focus on one of these models. As evident in Figure 4, the \(\bar{T}_R(k)\) value for \(k=20\) is smaller than that for \(k=0\). Furthermore, this value is closer to the true \(T_R\). This indicates that using a larger \(k\) value enables a more accurate assessment.

These results offer empirical support for our theoretical statements in Theorem 5 and Remark 6: increasing \(k\) not only relaxes the sufficient condition of \(\delta\)ISS but also enhances the accuracy of recovery time estimation.

0.0.7.2 Experiment \(\rm{I}\rm{I}\): Two-Tank System

Here, we demonstrate the effectiveness of our approach from the perspective of (B) in a more practical scenario using the two-tank system, a recognized benchmark in control experiments [22], [23].

0.0.7.2.1 Experimental Setup

0.0.7.2.1.1 System

The two-tank system is described by the following differential equations [23]: \[\begin{align} \frac{dh_1}{dt}&=-p_1\sqrt{2gh_1}+p_2u,\\ \frac{dh_2}{dt}&=p_3\sqrt{2gh_1}-p_4\sqrt{2gh_2}, \end{align}\] where \(u\) is the control input, \(h_1\) and \(h_2\) are the tank levels, \(g\) is the gravitational acceleration and \(p_1\), \(p_2\), \(p_3\), and \(p_4\) are hyperparameters depending on the system properties. The LSTM model is designed to predict the water level at the next time step based on the current water level and control input. The current water level is obtained from sensors, which are assumed to be subject to sensor noise. Formally, we define the input-output relationship of the LSTM as \(x(t) =(u(t), h(t)+w(t))\) and \(y(t) = h(t+1)\), where \(w \sim \mathcal{N}(0, {0.01}^2)\).

0.0.7.2.1.2 Baseline

To evaluate our proposed recovery time estimation method and adjustment technique, we prepared a baseline model: a standard model trained without considering resilience. As shown in Figure 5, the baseline model closely approximates the target system for inference.

0.0.7.2.2 Results and Analysis

We will compare the models trained using our proposed loss function 11 against the baseline model.

First, regarding (A), although the result in Figure 6 is not contradict Theorem 5 in the sense of weakly decreasing, the magnitude of this decay was notably smaller than the results in Section 0.0.5.1.

Next, concerning (B), we observed a trade-off between recovery times and MAE loss when comparing the baseline model and our models, and when varying \(\varepsilon\) of our model, as illustrated in Figure 6. Specifically, as \(\varepsilon\) increases from \(0.1\) to \(0.5\), we observe a trend of increasing MAE loss and decreasing recovery time. This trend aligns with our theoretical predictions in Section 0.0.6.2.

The experimental results demonstrate that LSTM resilience can be tuned via \(\varepsilon\), validating our proposed method 12 that introduces recovery time as a quality assurance metric and balances inference accuracy with recovery time.

0.0.8 Discussion↩︎

• Theoretical and data-independent guarantees may not always reflect actual data distributions. Therefore, we strongly recommend complementing our approach with data-driven quality assurance methods for practical implementation. Also, be aware that applying our method to an unstable system can be counterproductive.

• As tasks become increasingly dependent on long-term memory, the decline in the accuracy of recovery time estimation and the accuracy degradation caused by the constraints for resilience are likely to emerge as significant challenges. Relaxing the \(\delta\)ISS sufficient conditions could be a key approach to address these challenges. This is because, as demonstrated in Theorem 5, the \(\delta\)ISS conditions form the foundation for the recovery time estimation, and relaxing these conditions would contribute to weakening model constraints and mitigating inference accuracy degradation.

• Although our current research focused on the resilience of LSTM modules in control systems, comprehensive quality assurance for control systems remains a challenging issue. Refining studies that have demonstrated the relationship between LSTM stability and control system stability [14] could be a key approach. By more precisely analyzing the calculations used to prove control system stability in these studies, it may be possible to quantitatively determine the convergence time of the state of the system to the target value.

Acknowledgment↩︎

We are grateful to Professor Tomoyuki Tanaka of Yokohama National University for helpful discussions and comments on the manuscript.

Yoshihara Sota was financially supported by JST SPRING, Grant Number JPMJSP2125. He would like to take this opportunity to thank the “THERS Make New Standards Program for the Next Generation Researchers”.

Appendix↩︎

0.0.9 Other Related Work↩︎

In this paper, we define recovery time as the duration from when different input sequences align until the output difference falls below a threshold. We prove that LSTMs with \(\delta\)ISS properties can effectively estimate this time.

Finite time stability (FTS) [24] rigorously discusses convergence time, characterized by a settling-time function. Finite-time input-to-state stability (FTISS), combining FTS and ISS, enables robust control against disturbances [25]. Recent studies [26], [27] propose Lyapunov functions satisfying FTISS, eliminating the need to directly consider the settling-time function.

However, deriving sufficient FTS conditions involves strict model constraints, established only for specific models like Hopfield networks [28] and Clifford-valued RNNs [29], but not for conventional RNNs or LSTMs. While achieving FTS for LSTMs would theoretically allow for precise recovery time evaluation, we argue that from a practical standpoint, considering numerical errors and rounding, it suffices for output differences to fall below a threshold rather than strictly reaching equilibrium. Thus, we opted for \(\delta\)ISS instead of FTS.

0.0.10 Proofs↩︎

0.0.10.1 Proof of Proposition 1

We show the following lemma before proving Proposition 1.

Here we restate and prove Proposition 1:

0.0.10.2 Approximation of \(\eta^{(l)}(\infty)\)

Denote the limits of sequences \(\{\overline{\sigma}^{(l)}_*(k)\}_{k=0}^{\infty}\) (\(*\) denotes \(f,i\), or \(o\)), \(\{\overline{\phi}^{(l)}_c(k)\}_{k=0}^{\infty}\), and \(\{\eta^{(l)}(k)\}_{k=-1}^{\infty}\) as \(\overline{\sigma}^{(l)}_*(\infty)\), \(\overline{\phi}^{(l)}_c(\infty)\) and \(\eta^{(l)}(\infty)\), respectively. Since \(\overline{\sigma}^{(l)}_*(\infty)=\sigma(G^{(l)}_*(\eta^{(l)}(\infty)))\) and \(\overline{\phi}^{(l)}_c(\infty)=\phi(G^{(l)}_c(\eta^{(l)}(\infty)))\), we have \[\begin{align} \eta^{(l)}(\infty)&=\phi\left(\frac{\sigma(G^{(l)}_i(\eta^{(l)}(\infty)))\phi(G^{(l)}_c(\eta^{(l)}(\infty)))}{1-\sigma(G^{(l)}_f(\eta^{(l)}(\infty)))}\right)\sigma(G^{(l)}_o(\eta^{(l)}(\infty))). \label{eq:eta-infty} \end{align}\tag{15}\]

If the equation 15 could be solved, that is, \(\overline{\sigma}^{(l)}_*(\infty)\) and \(\overline{\phi}^{(l)}_c(\infty)\) could be obtained, we would find the loosest sufficient condition for the LSTM to be \(\delta\)ISS, and the optimal estimate of recovery time. However, it is difficult to solve it analytically, so we approximate it using properties of convex functions.

Let \(w \geq 0\). Since \(\sigma(w)\) is convex upward and \(\frac{d\sigma(w)}{dw} = \sigma(w) (1 - \sigma(w))\), then for \(w\) and \(w_0 \ge 0\), \[\sigma(w) \leq \sigma(w_0) + \sigma(w_0) \left( 1 - \sigma(w_0) \right) (w - w_0) \eqqcolon \widecheck{\sigma}(w; w_0).\] As \(\frac{d\phi(w)}{dw} = 1 - \phi^2(w)\), it holds that \[\phi(w) \leq \phi(w_0) + (1 - \phi^2(w_0)) (w - w_0) \eqqcolon \widecheck{\phi}(w; w_0).\] Therefore, for \(w\) and \(w_0\ge 0\), we have \[\begin{align} \sigma(G^{(l)}_*(w)) &\leq \widecheck{\sigma}(G^{(l)}_*(w); G^{(l)}_*(w_0)),\\ \phi(G^{(l)}_c(w)) &\leq \widecheck{\phi}(G^{(l)}_c(w); G^{(l)}_c(w_0)). \end{align}\] From these inequalities, the following inequality holds: \[\begin{align} &\sigma(G^{(l)}_o(w)) \phi\left( \frac{ \sigma(G^{(l)}_i(w)) \phi(G^{(l)}_c(w)) }{1 - \sigma(G^{(l)}_f(w))} \right) \eqqcolon \overline{g}^{(l)}(w) \\ &\leq \widecheck{\sigma}(G^{(l)}_o(w); G^{(l)}_o(w_0)) \, \widecheck{\phi}\left( \frac{ \widecheck{\sigma}(G^{(l)}_i(w); G^{(l)}_i(w_0)) \widecheck{\phi}(G^{(l)}_c(w); G^{(l)}_c(w_0)) }{1 - \sigma(G^{(l)}_f(w))} ; \frac{ \sigma(G^{(l)}_i(w_0)) \phi(G^{(l)}_c(w_0)) }{1 - \sigma(G^{(l)}_f(w_0))} \right) \\ &\eqqcolon {\widecheck{g}}^{(l)}(w; w_0). \end{align}\] Here, \(G_*(w)\), \(\widecheck{\sigma}(w, w_0)\), \(\widecheck{\phi}(w, w_0)\), and \((1 - \sigma(w))^{-1}\) are non-negative, monotonically increasing, and convex downward functions. Since \(\widecheck{g}^{(l)}(w; w_0)\) is represented by the composition and product of these functions, it is also convex downward. Hence, for \(0 \leq \kappa^{(l)}\leq \eta^{(l)}\), \(0 \leq \lambda^{(l)}\leq 1\), and \(\eta^{{(l)}+} \mathrel{\vcenter{:}}= (1 - \lambda) \kappa^{(l)}+ \lambda \eta^{(l)}\), it follows that \[\overline{g}^{(l)}(\eta^{{(l)}+}) \leq \widecheck{g}^{(l)}(\eta^{{(l)}+}; \eta^{(l)}) \leq (1 - \lambda^{(l)}) \widecheck{g}^{(l)}(\kappa^{(l)}; \eta^{(l)}) + \lambda^{(l)}\overline{g}^{(l)}(\eta^{(l)}).\] In particular, if \(\eta^{(l)}(\infty)\) holds \(\overline{g}^{(l)}(\eta^{(l)}(\infty)) = \eta^{(l)}(\infty)\) and \(\kappa^{(l)}\leq \eta^{(l)}(\infty) \leq \eta^{(l)}\), there exists \(\lambda^{(l)}(\infty)\) such that \(\eta^{(l)}(\infty) = (1 - \lambda^{(l)}(\infty)) \kappa^{(l)}+ \lambda^{(l)}(\infty) \eta^{(l)}\). In this case, we have the following: \[\begin{gather} (1 - \lambda^{(l)}(\infty)) \kappa^{(l)}+ \lambda^{(l)}(\infty) \eta^{(l)} \leq (1 - \lambda^{(l)}(\infty)) \widecheck{g}^{(l)}(\kappa^{(l)}; \eta^{(l)}) + \lambda^{(l)}(\infty) \overline{g}^{(l)}(\eta^{(l)}), \\ \lambda^{(l)}(\infty) \leq \frac{ \widecheck{g}^{(l)}(\kappa^{(l)}; \eta^{(l)}) - \kappa^{(l)} }{ \eta^{(l)}- \kappa^{(l)}- \overline{g}^{(l)}(\eta^{(l)}) + \widecheck{g}^{(l)}(\kappa^{(l)}; \eta^{(l)}) }, \\ \eta^{(l)}(\infty) \leq \frac{ \eta^{(l)}\, \widecheck{g}^{(l)}(\kappa^{(l)}; \eta^{(l)}) - \kappa^{(l)}\, \overline{g}^{(l)}(\eta^{(l)}) }{ \eta^{(l)}- \kappa^{(l)}- \overline{g}^{(l)}(\eta^{(l)}) + \widecheck{g}^{(l)}(\kappa^{(l)}; \eta^{(l)}) } = \eta^{(l)}- \frac{ (\eta^{(l)}- \kappa^{(l)}) (\eta^{(l)}- \overline{g}^{(l)}(\eta^{(l)})) }{ \eta^{(l)}- \kappa^{(l)}- \overline{g}^{(l)}(\eta^{(l)}) + \widecheck{g}^{(l)}(\kappa^{(l)}; \eta^{(l)}) }. \end{gather}\]

Therefore, we can approximate \(\eta^{(l)}(\infty)\) by following three steps:

1. Take \(\kappa^{(l)}\) and \(\eta^{(l)}\) such that \(\kappa^{(l)}\leq \eta^{(l)}(\infty) \leq \eta^{(l)}\), for example, \(\kappa^{(l)}=\overline{g}^{(l)}(0)\) and \(\eta^{(l)}=\overline{g}^{(l)}(1)\).

2. Calculate \(\dfrac{ \eta^{(l)}\, \widecheck{g}(\kappa^{(l)}; \eta^{(l)}) - \kappa^{(l)}\, \overline{g}^{(l)}(\eta^{(l)})}{ \eta^{(l)}- \kappa^{(l)}- \overline{g}^{(l)}(\eta^{(l)}) + \widecheck{g}^{(l)}(\kappa^{(l)}; \eta^{(l)})}\).

3. Substitute the result of Step 2 for \(\eta^{(l)}\). Return to Step 2.

0.0.10.3 Proof of Theorem 2

In this section, we prove the input-to-state stability and the incremental input-to-state stability of LSTM.

First, we introduce the definition of input-to-state stability.

Remarking the invariant set in Proposition 1, the sufficient condition for ISS of LSTM is relaxed (in the case of \(k \geq 1\)) than [14] as follows.

The next is the same as Theorem 2.

0.0.10.4 Proof of Theorem 5

0.0.11 Details of Experimential Setup↩︎

This section provides a detailed description of the experimental settings used in the experiments discussed in the main text. The experiments were conducted using the PyTorch framework.

0.0.11.1 Simplified Model

0.0.11.1.1 LSTM Setup

The model dimensions are provided in Table 1. Twenty models were generated with weights sampled uniformly from the intervals \(W_* \in [0,0.1]\), \(U_* \in [0, 1.0]\), and \(C \in [0,1.0]\).

0.0.11.1.2 Evaluation

We define the nominal input data \(x\) and the perturbed input data \(\widehat{x}\) that satisfy the assumptions of Theorem 4.3 with \(t_0=20\) as follows. \[\begin{align} x(t) = 0.5 , \;\widehat{x}(t) = \begin{cases} 1.0 & \text{if } t = 20, \\ 0.5 & \text{else}, \end{cases} \end{align}\] where \(t=0,1,2,\ldots,99\). \(T_R(x,\widehat{x})\) and \(\bar{T}_R(k)\) were calculated according to Algorithm 7 and Algorithm 8 where \(\mathcal{T} =100\). It should be noted that in cases where the model exhibits instability or \(\rho(A_s(k)) \ge 1\), this algorithm computes \(T_R = \infty\) or \(\bar{T}_R(k) = \infty\). Calculating statistical measures such as estimation errors or correlation coefficients using infinite values would result in arbitrary and meaningless numbers. Indeed, for some values of \(k\) and specific models, our experiments showed that \(\bar{T}_R(k) = \infty\) and \(\rho(A_s(k)) \ge 1\) occurred. This result is theoretically correct and does not indicate any deficiency in the theory. However, using infinite values to calculate statistical measures such as estimation errors or correlation coefficients would result in arbitrary and meaningless numbers. To address this issue, we replaced infinity values with finite proxies: \(\bar{T}_R(k)=\mathcal{T}=100\) when computing these statistical measures.

0.0.11.2 Two Tank System

0.0.11.2.1 LSTM Setup

The model dimensions are provided in Table 2.

0.0.11.2.2 Data Generation

The training and test data were randomly generated through the following process [22], [23]. First, the continuous two-tank system \[\begin{align} \frac{dh_1}{dt}&= -p_1\sqrt{2gh_2}+p_2u,\\ \frac{dh_2}{dt}&= p_3\sqrt{2gh_1}-p_4\sqrt{2gh_2}, \end{align}\] was discretized using the Runge-Kutta method with \(dt = 0.01\). The system parameters are shown in Table 3. We randomly generated a sequence of control inputs \(\{u(t)\}_t\), set the initial state of the system to \((h_1(0), h_2(0)) = (0,0)\), and obtained the system states \(\{(h_1(t),h_2(t))\}_t\). Here, the generation of control inputs was set to a time series length of \(100,000\), with values randomly switched every \(4,000\) step following a uniform distribution \(\mathcal{U}(0.5, 3.0)\). Using this method, we created two independent control input sequences, one for training and one for testing.

As the final stage of the data generation process, we created training and test datasets. In this step, we simulated more realistic scenarios by adding noise to the previously generated state variables. Specifically, we constructed the input data as follows: \[\begin{align} x(t) = (u(t), h_1(t) + w_1(t), h_2(t)+w_2(t)). \end{align}\] This addition of noise simulates real-world factors such as sensor measurement errors and environmental uncertainties. The noise distribution was set to \(w_1,w_2 \sim \mathcal{N}(0,{0.1}^2)\) for the training data and to \(w_1,w_2 \sim \mathcal{N}(0,{0.01}^2)\) for the test data. This is because while noise with a standard deviation of about \(0.01\) is expected under operational conditions, intentionally introducing larger noise during the training phase aims to improve the model’s robustness and accuracy. On the other hand, the inference target data was set using the state of the next time step as follows: \[\begin{align} y(t) = (h_1(t+1), h_2(t+1)). \end{align}\] This setup requires the model to predict the state of the next time step from the current control input and the state including noise. This simulates the important task of one-step-ahead state prediction in actual system control.

The time series data are normalized to the interval \([-1, 1]\) for both input \(x(t)=(u(t),h_1(t)+w(t), h_2(t)+w(t))\) and output \(y(t)=(h_1(t+1),h_2(t+1))\) of the LSTM before conducting inference. The equation for normalization is given below: \[\begin{align} x &= 2 (x-x_\mathrm {min})/(x_\mathrm {max}-x_\mathrm {min}) -1, \\ y &= 2 (y-y_\mathrm {min})/(y_\mathrm {max}-y_\mathrm {min}) -1, \end{align}\] where \(x_\mathrm {min}= (0.5, 0,0)\), \(x_\mathrm {max}= (3.0, 10,10)\), \(y_\mathrm {min}= (0,0)\), and \(y_\mathrm {max}= (10,10)\).

0.0.11.2.3 Train

0.0.11.2.3.1 Data set

The training dataset is split into training and validation sets, with the first \(50\%\) of the sequences allocated for training and the remaining \(50\%\) used for validation purposes. To reduce computational load and augment data, the training dataset is created by segmenting one length time series \(50,000\) into individual sequences \(49,000\), each with a length of \(1,000\), using a sliding window of step size \(1\). On the other hand, as data augmentation is not required for the validation dataset, we create the dataset using both a window length and a step size of 1,000, solely considering computational efficiency.

0.0.11.2.3.2 Loss

We design the loss function for backward as follows: \[\begin{align} \mathcal{L} = \text{MAE}(y_{\textrm{true}}, y_{\textrm{pred}}) + \lambda \Phi(\theta,\varepsilon). \end{align}\] When calculating the MAE on the batched data, we discard the first 10 steps of the time series predictions. This adjustment accounts for the inherent instability of LSTM predictions during the initial time steps, which is a characteristic behavior of LSTM networks.

0.0.11.2.3.3 Hyperparameters

We used the following hyperparameters. The weights were initialized with uniform distribution \(\mathcal{U}(-1/\sqrt{n_c},1/\sqrt{n_c})\) and optimized using the Adam (\(\beta =(0.9, 0.999)\), \(\varepsilon=10^{-8}\)) optimizer. The mini-batch size was \(64\), the epoch size was \(200\), the learning rate was \(0.001\) with ReduceLROnPlateau (\(\text{mode}=\text{min}\), \(\text{factor}=0.1\), \(\text{patience}=10\)) scheduler. The baseline model was constructed with \(\lambda = 0\), while the stability-guaranteed model was implemented with \(\lambda = 1.0\). The recovery time control parameter \(\varepsilon\) was selected from the set \(\{0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45, 0.5 \}\).

0.0.11.2.4 Evaluation

The inference accuracy was evaluated using MAE as the metric for the test data. For the time series, we calculated the MAE, ignoring the first \(10\) steps and using non-normalized values for the evaluation.

Separate from the data used for evaluating inference accuracy, we generate the nominal input data \(x\) and the perturbed input data \(\widehat{x}\) for the estimation of recovery time \(T_R\) using the following steps. First, we set the time series length to \(3500\) steps, the control input to \(u = 1.0\), and the initial water level of the system to \((h_1(0),h_2(0)) = (0,0)\). Then, following the data generation process, we created \((h_1(t), h_2(t))\) and set \(x(t) = (u(t),h_1(t), h_2(t))\). Then, we define \(\widehat{x}\) that satisfies the assumptions of Theorem 5 with \(t_0=1010\) as \[\begin{align} \widehat{x}(t) = \begin{cases} (u(t), h_1(t)+ 1.0, h_2(t) +1.0) & \text{if } 1000 \le t \le 1010, \\ x(t) & \text{else} .\end{cases} \end{align}\] This definition models the phenomenon where observational noise temporarily has a magnitude larger than expected. Based on these settings, we calculate \(T_R\) and \(\bar{T}_R\). The algorithm is similar to the one used with the toy model, but there is one point to note due to the effect of data normalization. As the LSTM output is normalized, \(T_R\) is calculated to ensure the error remains within e on the original scale. Consequently, \(\bar{T}_R\) must also be adjusted to the original scale. In the calculation of \(\bar{T}_R\), we apply the following scale transformation to \(e\) in Theorem 5: \[\begin{align} e \gets e/(h_\mathrm {max}- h_\mathrm {min}) = e/(10-0) =e/10. \end{align}\]

0.0.12 Theoretical Aproaches Compare to Data-Driven Aproaches↩︎

We conducted comparative experiments between our proposed theoretical method and conventional data-driven approaches. In the data-driven approach, training was performed by randomly adding pulses along the time series direction to the training data.

0.0.12.1 Setting

0.0.12.1.1 TRAIN

0.0.12.1.1.1 Data set

For the adversarial learning, training and validation data were generated by adding pulse signals to the input data generated in Section 0.0.11.2, excluding control inputs. For data with time series length \(T=50,000\), we generated randomly occurring sustained pulses. Each pulse occurrence follows a Poisson process with an average rate of \(\lambda = 0.001\) per unit time. The total number of pulses over the entire period \(T\) is determined according to the Poisson distribution Poisson\((\lambda T)\).Each pulse has the following characteristics:

• The start time is randomly selected from a uniform distribution in the range \([0, T-d]\), where \(d=10\) is the duration of each pulse.

• The amplitude is sampled from a uniform distribution within \([0, M]\), where \(M\) is the maximum intensity of the pulse.

• Each pulse is added to the input signal at a constant value for \(d\) steps from the selected start time.

We created training and validation data with variations of \(M = 1, 3, 10\) using the method described above. Models trained with each of these datasets are referred to as "DD pulse \(M\) model" (where DD stands for data-driven), while our model was trained using data without added pulses.

0.0.12.1.1.2 Loss

We used the same loss function as defined in Section 0.0.11.2.

0.0.12.1.1.3 Hyperparameters

With the exception of the following two points, we used the same hyperparameters as in Section 0.0.11.2:

• For our model, \(\varepsilon=0.05\).

• For the DD pulse \(M\) models, \(\lambda=0\).

0.0.12.1.2 EVALUATION

We evaluated each model using the following evaluation data, under the same settings as in Section 0.0.11.2.

We define \(\widehat{x}_p\) that satisfies the assumptions of Theorem 5 with \(t_0=1010\) as \[\begin{align} \widehat{x}_p(t) = \begin{cases} (u(t), h_1(t)+ p, h_2(t) + p) & \text{if } 1000 \le t \le 1010, \\ x(t) & \text{else}, \end{cases} \end{align}\] where \(p\) was selected from the set \(\{1.0, 5.0, 9.0\}\), and let \(\widehat{y}_p\) represent the output of the LSTM when \(\widehat{x}_p\) is used as the input.

0.0.12.2 Results and Analysis

As shown in Table 4, regarding \(T_R\), our proposed method demonstrated superiority. Our model achieved faster recovery times compared to models obtained using data-driven methods. Although models trained with data having stronger noise pulses showed improvement in recovery time, they could not reach the performance of our theoretically guaranteed model. These results indicate that our proposed method is effective in enhancing LSTM’s recovery capability.

In achieving stability criteria, our model with theoretical guarantees also showed clear advantages. Models trained with these data-driven techniques failed to meet the theoretical stability criterion of \(\rho(A_s)<1\). This result suggests that approaches relying solely on data are insufficient for theoretical guarantees of system resilience.

On the other hand, in terms of noise resistance, models trained with data-driven approaches outperformed our model in certain cases. This is likely because our proposed method did not directly aim to maximize noise resistance. Noise resistance can be theoretically bounded using \(\gamma\) from Definition 2. As a future challenge, by explicitly deriving \(\gamma\), there is potential to improve noise resistance metrics.

In general, our theoretical approach has been proven to provide a better solution than data-driven methods to ensure resilience under any noise conditions.

References↩︎