Stability of implicit neural networks for long-term forecasting in
  dynamical systems
Forecasting physical signals in long time range is among the most challenging
tasks in Partial Differential Equations (PDEs) research. To circumvent
limitations of traditional solvers, many different Deep Learning methods have
been proposed. They are all based on auto-regressive methods and exhibit
stability issues. Drawing inspiration from the stability property of implicit
numerical schemes, we introduce a stable auto-regressive implicit neural
network. We develop a theory based on the stability definition of schemes to
ensure the stability in forecasting of this network. It leads us to introduce
hard constraints on its weights and propagate the dynamics in the latent space.
Our experimental results validate our stability property, and show improved
results at long-term forecasting for two transports PDEs.
Leon Migus, Julien Salomon, Patrick Gallinari
# Stability of Implicit Neural Networks for Long-term Forecasting in Dynamical Systems

###### Abstract

Forecasting physical signals in long time range is among the most challenging tasks in Partial Differential Equations (PDEs) research. To circumvent limitations of traditional solvers, many different Deep Learning methods have been proposed. They are all based on auto-regressive methods and exhibit stability issues. Drawing inspiration from the stability property of implicit numerical schemes, we introduce a stable auto-regressive implicit neural network. We develop a theory based on the stability definition of schemes to ensure the stability in forecasting of this network. It leads us to introduce hard constraints on its weights and propagate the dynamics in the latent space. Our experimental results validate our stability property, and show improved results at long-term forecasting for two transports PDEs.

## 1 Introduction and motivation

Numerical simulations are one of the main tools to study systems described by PDEs, which are essential for many applications including, e.g., fluid dynamics and climate science. However, solving these systems and even more using them to predict long term phenomenon remains a complex challenge, mainly due to the accumulation of errors over time. To overcome the limitations of traditional solvers and to exploit the available data, many different deep learning methods have been proposed.

For the task of forecasting spatio-temporal dynamics, Ayed et al. (2019) used a standard residual network with convolutions and Sorteberg et al. (2019); Lino et al. (2020); Fotiadis et al. (2020) used Long short-term memory (LSTM) and Convolutional neural network (CNN) for the wave equation. In Wiewel et al. (2019); Kim et al. (2019), a good performances is obtained by predicting within the latent spaces of neural networks. More recently, Brandstetter et al. (2022) used graph neural networks with several tricks and showed great results for forecasting PDEs solutions behavior. Importantly, these methods all solve the PDE iteratively, meaning that they are auto-regressive, the output of the model is used as the input for the next time step. Another line of recent methods that have greatly improved the learning of PDE dynamics are Neural Operators (Li et al., 2020). These methods can be used as operators or as auto-regressive methods to forecast. However, when used as operators, they do not generalize well beyond the times seen during training. Crucially, these auto-regressive methods tend to accumulate errors over time with no possible control, and respond quite poorly in case of change in the distribution of the data. This leads to stability problems, especially over long periods of time beyond the training horizon.

In the numerical analysis community, the stability issue has been well studied and is usually dealt with implicit schemes. By definition, they imply to solve an equation to go from a time step to the next one but they are generally more stable than explicit schemes. This can be seen on the test equation \(\frac{\mathrm{d}y}{\mathrm{d}t}=\lambda y\), where Euler implicit schemes are always stable while Euler explicit schemes are not. Interpreting residual neural networks as numerical schemes, one can apply such schemes and gain theoretical insights on the properties of neural networks. This has already been done in various forms in Haber and Ruthotto (2017); Chen et al. (2018), but not applied to forecasting. Moreover,
these works use either the stability of the underlying continuous equation or the stability of the numerical scheme on the test equation and its derivatives, which is not the stability of the numerical scheme on the studied equation. Since a network is discrete, the latter is the most relevant.

We therefore use the stability in norm of schemes, as defined in 2.1. In deep learning (DL), this definition has only been applied to image classification problems (Zhang and Schaeffer, 2020). To the best of our knowledge, this work is the first attempt to forecast PDEs with neural networks using stability as studied in the numerical analysis community.

Using implicit schemes in DL has already been done in different contexts. The earliest line of works tackles image problems, with Haber et al. (2019) designing semi-implicit ResNets and Li et al. (2020); Shen et al. (2020); Reshniak and Webster (2021) designing different implicit ResNets. The second line of works focuses on dynamical problems. In this way, Nonnenmacher and Greenberg (2021) designed linear implicit layers, which learn and solve linear systems, and Horie and Mitsume (2022) used an implicit scheme as part of an improvement of graph neural network solvers to improve forecasting generalization with different boundary condition shapes. Tackling different types of problems, none of these methods guarantees the forecast stability. For our analysis, we restrict ourselves to ResNet-type networks, i.e., networks with residual connections. We introduce hard constraints on the weights of the network and predict within the latent space of our network. Hence, by modifying the classic implicit ResNet architecture, our method can forecast dynamical system at long range without diverging. We apply these theoretical constraints in our architecture to two 1D transport problems: the _Advection equation_ and _Burgers' equation_.

To better investigate networks stability, we perform our experiments under the following challenging setting: for training, we only give to the networks the data from \(t=0\) to a small time \(t=\Delta t\), and consider the results in forecasting in the long time range at \(t=N\cdot\Delta t\), with \(N\gg 1\). Note that our setting is harder that the conventional setting presented for e.g. in Brandstetter et al. (2022). Indeed, we only use changes between a single time step for training.

## 2 Method

To guarantee structural forecasting stability, we take inspiration from implicit schemes. We focus our study on an implicit ResNet with a ReLU activation function. In our approach, an equation is solved at each layer, namely \(x_{n+1}:=x_{n}+R_{n}(x_{n+1})\) with \(x\) in \(\mathbb{R}^{M}\) and \(n\) in \(\mathbb{N}\) and \(R_{n}(x)=\text{ReLU}(W_{n}x+b_{n})\) where \(W_{n}\) is upper triangular. The latter constraint is motivated below.

### Implicit ResNet stability analysis

To ensure that our proposed architecture is well-defined, we need to solve \(x=x_{n}+R_{n}(x)\). This equation has a solution, as proven in El Ghaoui et al. (2019) and detailed in Appendix 5.1. We can then study its stability. The recursion defining \((x_{n})_{n\in\mathbb{N}}\) reads as an implicit Euler scheme with a step size of 1. As described in the introduction, an implicit scheme is usually more stable than an explicit one. We first recall the definition of stability for schemes. This property ensures that our architecture has an auto-regressive stability.

**Definition 2.1** (Stability in norm).: _A scheme \((x_{n})_{n\in\mathbb{N}}\) of dimension \(M\) is stable in norm \(L^{p}\) if there exists \(C(T)\) independent of the time discretization step \(\Delta t\) such that:_

\[\forall\;x_{0}\in\mathbb{R}^{M},\;n\geq 0;\;n\Delta t\leq T,\;\|x_{n}\|_{p} \leq C(T)\|x_{0}\|_{p}\,.\]

This general definition of stability in norm ensures that a scheme does not amplify errors. This definition is equivalent to several others in the numerical analysis community.

Suppose that the spectrum of \(W_{n}\) is contained in \([-1,0[\) for every integer \(n\), we can assert that \((x_{n})_{n\in\mathbb{N}}\) is well-defined, using theorem 5.2. The proof of the stability of our Implicit ResNet network is then by induction on the dimension and is given in appendix 5.4:

**Theorem 2.1** (Stability theorem).: _If the elements of the diagonal of \(W_{n}\) are in \([-1,0)\) for every integer \(n\), then \((x_{n})_{n\in\mathbb{N}}\) is stable in norm \(L^{p}\)._

This theorem leads to hard constraints on the weights of the network.
### Implementation

To validate practically our theoretical results, we choose a setting that highlights stability issues. We then test our implementation of an implicit ResNet. In order to respect the assumptions of theorem 2.1, we forecast the dynamics in the latent space, as detailed below.

SettingWe first learn the trajectory at a given small time step \(\Delta t\). We only give data \(t=0\) to \(t=\Delta t\) for the training. We then forecast in long-term, at \(N\cdot\Delta t\) with \(N\gg 1\). This very restricted setting allows us to see how the network will react in forecasting with changes in the distribution and error accumulation. Usually neural network forecasting methods use data from \(t=0\) to \(T=L\cdot\Delta t\) for the training which allows to use different tricks to stabilize the network, such as predicting multiple time steps at the same time. However, in this work, we want to analyze how the network behaves without any trick that can slow down divergence. Indeed, the tricks used in the other settings do not actually guarantee stability of the network. The training is performed with a mean-squared error (MSE) loss.

Implicit neural network architectureTo implement a neural network from Theorem 2.1, we use the following architecture; \(z_{\Delta t}=f_{dec}\circ f^{K}_{res}\circ...\circ f^{1}_{res}\circ f_{enc}(z_ {0})\), with \(f^{k}_{res}(x)=x+\text{ReLU}(W_{k-1}\cdot f^{k}_{res}(x)+b_{k-1})\). The encoder and decoder are linear, and the encoder projects the data to a smaller dimension \(M\). The full architecture is illustrated in Figure 2. The specificity of our architecture is that the residual blocks are connected with an implicit Euler scheme iteration. To do so, we use a differentiable root-finding algorithm (Kasim and Vinko, 2020). More precisely, we use the first Broyden method van de Rotten and Lunel (2005) to solve the root-finding problem. It is a quasi-Newton method. This helped getting better results compared to other algorithms.

Latent space forecastingAs in Wiewel et al. (2019); Kim et al. (2019), the forecast can be done within the latent space of the network; \(z_{N\cdot\Delta t}=f_{dec}\circ f_{block}\circ...\circ f_{block}\circ f_{enc}(z_ {0})\), with \(f_{block}=f^{K}_{res}\circ...\circ f^{1}_{res}\). To predict at time \(N\cdot\Delta t\) from time \(t=0\), we apply \(N\) times the residual blocks. It is illustrated in Figure 4. This propagation allows our network to respect the assumptions of theorem 2.1 and thus be stable.

## 3 Experiments

We evaluate the performance of our method on the _Advection equation_ and _Burgers' equation_.

Datasets_Advection equation_. We construct our first dataset with the following 1-D linear PDE, \(\frac{\partial\psi}{\partial t}=-\frac{1}{2}\frac{\partial\psi}{\partial z},z \in(0,2\pi)\), \(t\in\mathbb{R}^{+}\) and \(\psi(z,0)=f_{0}(z),z\in(0,2\pi)\).

Burgers' equation.In the second dataset, we consider the following 1-D non-linear PDE, \(\frac{\partial\psi}{\partial t}=-\frac{1}{2}\frac{\partial\psi^{2}}{\partial z }+\frac{\partial^{2}\psi^{2}}{\partial z},z\in(0,2\pi),\)\(t\in(0,1]\) and \(\psi(z,0)=f_{0}(z),z\in(0,2\pi)\).

Both equations have periodic boundary conditions. We approximate the mapping from the initial condition \(f_{0}\) to the solution at a given discretization time \(\Delta t\), i.e. \(u_{0}\mapsto u(\cdot,\Delta t)\). We then forecast to longer time ranges. \(\Delta t\) is set to 1 for the _Advection equation_ with a grid of 100 points and 0.0005 for _Burgers' equation_ with a grid of 128 points.

Baseline methodsWe compare our Implicit ResNet with respect to an Explicit ResNet with ReLU activation function and a Fourier Neural Operator (FNO) (Li et al., 2020). We have also implemented two Explicit ResNet, with a tanh activation function and with batch normalization. We design the Explicit ResNet in the same way as our implicit version, with \(K\) layers of residual blocks that are linked by \(x_{n+1}=x_{n}+R_{n}(x_{n})\). Traditional methods forecast by using the output of the network at time \(t\) to predict the dynamics at time \(t+\Delta t\). So to predict at time \(N\cdot\Delta t\) from time \(t=0\), the baseline networks are iteratively applied \(N\) times, as illustrated in Figure 3.

ResultsPrediction errors are reported in Table 1 for the _Advection equation_ and _Burgers' equation_. We also show the error according to the forecast time in Figure 1.

We found that traditional deep learning methods diverge with the forecast time. They reach a MSE of more than \(10^{8}\) for the _Advection equation_ and go to infinity for _Burgers' equation_, respectively at
time \(400\) and \(0.15\). Moreover, we see in Figure 1 that their curve in time is convex, so the increase in error is accelerating over time. We also found that our proposed Implicit ResNet presents better results by several orders of magnitude for both datasets. Moreover, we can see in Figure 1 that its curve in time is reaching a stable plateau, as expected from our theorem 2.1.

As for the training, traditional deep learning methods manage to learn very well the dynamics at \(t=\Delta t\), with two orders of magnitude better than our Implicit ResNet. However, the latter still manages to learn well the dynamics with a MSE of \(10^{-2}\) for the _Advection equation_ and \(10^{-3}\) for _Burgers' equation_. This difference in training can mainly be explained by the longer training time of the Implicit ResNet, which made us take a smaller number of epochs for this network (1250 against 2500).

DiscussionFigure 1 demonstrates the main benefits of our constrained implicit neural network. Our network is stable whereas the other methods diverge in time. However, although being stable and far better than the baselines, it does not manage to forecast accurately the long-term dynamics. This is further confirmed by Table 4, which shows high relative errors. Said otherwise, when stability is guaranteed, convergence is not. We can also note that constraining the weights makes our network harder to train, but guarantees structural forecasting stability.

## 4 Conclusion

In this work, we studied the challenging task of long-term forecasting in dynamical systems. To do so, we developed a theoretical framework to analyze the stability of deep learning methods for forecasting dynamical systems. We then designed a constrained implicit neural network out of this framework. To the best of our knowledge, this is the first work that proposes to study deep learning architectures for forecasting dynamical systems from a numerical schema standpoint. We showed improved results with respect to deep learning baselines for two transport PDEs.
This work opens new perspectives to study neural networks forecasting stability from a numerical schema standpoint, thus offering more robust architectures. However, this analysis still needs improvements. Even though it ensures forecasting stability of our proposed network, it does not guarantee good convergence properties. We believe that developing this line of research could help overcome these challenges, and provide more robust architectures for forecasting dynamical systems in long time range.
## Acknowledgments

We thank Yuan Yin and Etienne Le Naour for helpful insights and discussions on this project.
