Using machine learning to predict extreme events in complex systems.

Extreme events and the related anomalous statistics are ubiquitously observed in many natural systems, and the development of efficient methods to understand and accurately predict such representative features remains a grand challenge. Here, we investigate the skill of deep learning strategies in the prediction of extreme events in complex turbulent dynamical systems. Deep neural networks have been successfully applied to many imaging processing problems involving big data, and have recently shown potential for the study of dynamical systems. We propose to use a densely connected mixed-scale network model to capture the extreme events appearing in a truncated Korteweg-de Vries (tKdV) statistical framework, which creates anomalous skewed distributions consistent with recent laboratory experiments for shallow water waves across an abrupt depth change, where a remarkable statistical phase transition is generated by varying the inverse temperature parameter in the corresponding Gibbs invariant measures. The neural network is trained using data without knowing the explicit model dynamics, and the training data are only drawn from the near-Gaussian regime of the tKdV model solutions without the occurrence of large extreme values. A relative entropy loss function, together with empirical partition functions, is proposed for measuring the accuracy of the network output where the dominant structures in the turbulent field are emphasized. The optimized network is shown to gain uniformly high skill in accurately predicting the solutions in a wide variety of statistical regimes, including highly skewed extreme events. The technique is promising to be further applied to other complicated high-dimensional systems.

Extreme events and their anomalous statistics are ubiquitous in various complex turbulent systems such as the climate, material science, and neuroscience, as well as engineering design (1–4). Understanding and accurate prediction of such phenomena remain a grand challenge, and have become an active contemporary topic in applied mathematics (5–8). Extreme events can be isolated rare events (2, 9, 10), or they can often be intermittent and even frequent in space and time (6, 8, 11, 12). The curse of dimension forms one important obstacle for the accurate prediction of extreme events in large complex systems (3, 4, 6, 13), where both novel models and efficient numerical algorithms are required. A typical example can be found in recent laboratory experiments for turbulent surface water waves going through an abrupt depth change revealing a remarkable transition to anomalous extreme events from near-Gaussian incoming flows (1).

A statistical dynamical model is then proposed in refs. 14 and 15 that successfully predicts the anomalous extreme behaviors observed in the shallow water wave experiments. The truncated Korteweg–de Vries (tKdV) equation is proposed as the governing equation modeling the flow surface displacement. Gibbs invariant measures are induced based on the Hamiltonian form of the tKdV equation to describe the probability distributions at equilibrium. A statistical transition from symmetric near-Gaussian statistics to a highly skewed probability density function (PDF) is achieved by simply controlling the “inverse temperature” parameter in the Gibbs measure (15).

In recent years, machine learning strategies, particularly the deep neural networks, have been extensively applied to a wide variety of problems involving big data, such as image classification and identification (16–19). On the other hand, the construction of proper deep learning strategies for the study of complex turbulent flows still remains an actively growing topic. The deep neural network tools developed for imaging processing have been suggested to be applied for data-driven predictions of chaotic dynamical systems (20–22), climate and weather forecasts (23, 24), and parameterization of unresolved processes (25–27). In the statistical prediction of extreme events, the available data for training are often restricted in limited regimes. A successful neural network is required to maintain adaptive skill in wider statistical regimes with vastly distinct statistics away from the training dataset. Besides, a working prediction model for turbulent systems would also require a prediction time scale longer than the decorrelation time that characterizes the mixing rate of the state variables.

In this paper, we investigate the extent of skill of the deep neural networks in predicting statistical solutions of complex turbulent systems, especially involving highly skewed PDFs. The statistical tKdV equations serve as a difficult first test model for extreme event prediction with simple trackable dynamics but a rich variety of statistical regimes from near-Gaussian to highly skewed PDFs showing extreme events. The important questions to ask are whether the deep networks can be trained to learn the complex hidden structures in the highly nonlinear dynamics purely from restricted data, and what are the essential structures required in the network to gain the ability to capture extreme events.

Our major goal here is to get accurate statistical prediction for the extreme events in time intervals significantly longer than the decorrelation time of the complex turbulent system. To achieve this, a convolutional neural network (mixed-scale dense neural network [MS-D Net]) which exploits multiscale connections and densely connected structures (28) is adopted to provide the basic network architecture to be trained using the model data from the tKdV equation solutions. This network enjoys the benefits of simpler model implementation and a smaller number of tunable training hyperparameters. Thus it becomes much easier to train, requiring less computational cost and technical tuning of the hyperparameters.

The key structures for the neural network to successfully capture extreme events include: 1) the use of a relative entropy (Kullback–Leibler divergence) loss function to calibrate the closeness between the target and the network output as distribution functions, so that the crucial shape of the model solutions is captured instead of a pointwise fitting in the turbulent output field values; and 2) calibrating the output data under a combination of empirical partition functions emphasizing the large positive and negative values in the model prediction so that the main features in the solutions are further emphasized. This convolutional neural network model enjoys the following major advantages for the prediction of extreme events: 1) The simple basic network architecture makes the model easier to train and efficient to predict the solutions among different statistical scenarios. 2) The network structure approximates the original model dynamics with the designed model loss function and treatment of output data, so it gives a better approximation to the complex structure of the system dynamics. 3) The temporal and spatial correlations in different scales are modeled automatically from the design of the network with convolution kernels representing different scales in different layers. 4) The method shows robust performance with different model hyperparameters, and can be generalized for the prediction of more complicated turbulent systems.

Direct numerical tests show high skill of the neural network in successfully capturing the extreme values in the solutions with the model parameters only learned from the near-Gaussian regime of vastly different statistics. The model also displays accurate prediction in much longer time beyond the decorrelation time scale of the state, proving the robustness of the methods. The successful prediction in the tKdV equation implies the potential of future applications of the network to more complicated high-dimensional systems.

Background for Extreme Events and Neural Network Structures

The tKdV Equations with Extreme Events.

The tKdV model provides a desirable set of equations capable of capturing many complex features in surface water wave turbulence with simple trackable dynamics. Through a high wavenumber cutoff at (with grid points), the Galerkin projected state induces stronger turbulent dynamics than the original continuous one (15). The tKdV equation has been adapted to describe the sudden phase transition in statistics (14) where highly skewed extreme events are generated from near-Gaussian statistics for waves propagating across an abrupt depth change. The tKdV model is formulated on a periodic domain aswhere the state variable represents the surface wave displacement to be learned directly using the deep neural network. The model is nondimensionalized using the characteristic scales as the total energy, as the length scale, and as the water depth. The steady-state distribution of the tKdV solution can be described by the invariant Gibbs measure derived from the equilibrium statistical mechanics (29),with a competition between the cubic term and the quadratic term from the Hamiltonian. The last term in Eq. is the delta function constraining the total energy conservation, . The only parameter as the inverse temperature determines the skewness of the PDF of (14). The Gibbs measures [] with different values of can be used to provide initial samples for the direct simulations of the model [. Different final equilibrium statistics (with various skewness) can be obtained based on the initial configuration of the ensemble set (that is, from picking different inverse temperature values ). A detailed description of the statistical tKdV model together with the simulation setup is provided in .

Training and prediction data from the same model with distinct statistics.

The basic idea in training the deep neural network is to use a training dataset with solutions sampled from Eq. using near-Gaussian statistics. The tKdV model dynamics can be learned from the training set without explicitly knowing the model dynamics. Then the question is, What is the range of skill in the trained neural network to predict the highly skewed non-Gaussian distributions among different sets of data? The training and prediction datasets are proposed from the ensemble solutions of the tKdV model [ based on the following strategy:

In the training dataset, we generate solutions from an ensemble simulation starting from a near-Gaussian PDF using a small absolute value of the inverse temperature (as shown in the first row and the near-Gaussian PDF in Fig. 1). On one hand, the model dynamics is represented by the group of solutions to be learned through the deep neural network. On the other hand, only near-Gaussian statistics is obtained in this training dataset, so the neural network cannot know about the skewed rare events appearing in other statistical regimes directly from the training process.

Fig. 1.

Solutions and statistics of the tKdV equation in three typical parameter regimes with different statistics. The first three rows plot solution trajectories in the three regimes with near-Gaussian (first row), mildly skewed (second row), and highly skewed (third row) statistics. The corresponding equilibrium PDFs of the three cases are shown in the fourth row. The autocorrelation functions and decorrelation time of each Fourier mode of the model state are compared in the last row.

For the prediction dataset, we test the model skill using the data generated from various different initial inverse temperatures (as shown in the second and third rows as well as the skewed PDFs in Fig. 1). It provides an interesting test bed to check the scope of skill in the optimized neural network for capturing the distinctive statistics and extreme events.

The choices of the training and prediction datasets are illustrated in Fig. 1, which first shows, in the first three rows, realizations of the tKdV model solutions from different inverse temperatures . A smaller amplitude of gives near-Gaussian statistics in the model state , while larger amplitudes of give trajectories with skewed PDFs. The corresponding equilibrium PDFs of the state from direct ensemble dynamical solutions are compared next to illustrate explicitly the transition in statistics. Turbulent dynamics with multiscale structures are observed in all of the tKdV solutions. The autocorrelation function and the decorrelation time are plotted in the last row of Fig. 1, confirming the rapid mixing in the solutions.

Data structures for the deep neural networks.

The deep convolutional networks can be viewed as a function , mapping the input signal with rows, columns, and channels to the output data with rows, columns, and channels. We consider an ensemble simulation of trajectories of the tKdV equation evaluated at the grid points and measured in a time interval . Thus the input data for the network comes from the ensemble solutions at the first time measurements at ,as a tensor with rows for the spatial grid points, columns for the discretized time evaluations, and only one channel for each of the input samples . An ensemble of total independent solutions from the Monte Carlo simulation is divided into minibatches to feed in the network in the training process. For simplicity, the output data are designed in the same shape as the predicted states evaluated at a later time starting from the previous initial data,The forwarding time controls how long we would like the network to push forward the states in one time update. For one effective neural network for the complex system, the time scale is expected to be longer than the decorrelation time, . The above construction is supposed to feed both the time and spatial correlations of the original dynamical model into the neural network to be learned in the approximation map .

Deep Convolutional Neural Network Architecture.

The basic structures of the convolutional neural network include the operators in each single convolution layer, and the connections between multiple layers. We would like to first keep the neural network in its simplest standard setup, so that we are able to concentrate on the improvement in key structures without risking getting lost in manipulating the various complicated ad hoc hyperparameters. More detailed convolutional network construction is described in , following the general neural network architecture as in refs. 19 and 28.

Basic convolutional neural network unit.

In each single convolution layer, the input data from the previous layer output are updated in the general formA convolutional operator is first applied on the input data in a small symmetric window with size , where the first dimension controls the correlation in the spatial direction and the second one controls the temporal correlation. A bias is added to the convolved data before applying a final nonlinear operator using the common choice of rectified linear unit function. The convolution kernel starts with a small size, (that is, using only the two nearest neighbor points in space and time), which enables fast computation and easy control. Naturally, periodic boundary condition is applied on the spatial dimension, and replicate boundary is added in time before and after for the boundary padding. No additional structures are implemented in the convolution layer unit, to keep the basic standard architecture used in imaging processing (19).

A densely connected and mixed-scale structure.

Next, we need to propose the connection between different layers. The common feedforward deep neutral network feeds the input data in the th layer only to the next th layer. The feedforward network often requires a larger number of layers to work; thus it is expensive to train and difficult to handle. Proper downscaling and upscaling steps going through the layers may also be required, while these downscaling and upscaling operations may not be a feasible approach for simulating the dynamical model time integration steps.

As an alternative approach, an MS-D Net is introduced, in ref. 28, by mixing different scales within each layer using a dilated convolution, and densely connecting all of the feature maps among all of the layers. First, to learn the multiscale structures, the convolution kernels in different layers are dilated differently by adding zeros between the values in the original kernel . The dilated convolutions become especially appealing for capturing the multiscale structures in the turbulent dynamics. Different spatial and temporal scales are included adaptively with different convolution length scales. Second, the dense network connection includes all of the previous layer information to update the output data in the next layer. In the implementation, all of the previous layer outputs are piled together as input channels for the next layer. Together with the multiscale convolution kernels used in different layers, the output in the next layer combines the information in different scales and produces a balanced update in the next step.

The MS-D Net requires fewer feature maps and trainable parameters, so it is easier to handle compared with the direct feedforward network. It provides a desirable setup for the prediction in dynamical systems by feeding in all of the data in previous layers decomposed into different-scale structures. Then dynamical structures at different scales communicate with each other through the dense network connection. Intuitively, this is a reasonable structure for the turbulent solutions, since all of the history information is useful for the prediction in the next steps. The densely connected network structure is also comparable to the time integration scheme, where all of the history information is used to update the state at the next time step without using any downscaling and upscaling steps.

A Learning Strategy for Extreme Event Prediction

In this section, we construct the specific network structures designed for learning turbulent system dynamics and then the prediction of extreme events. In this case, with data from the turbulent models, small fluctuations in the solutions may introduce large errors in optimizing the loss function. We aim at capturing the dominant emerging features such as extreme events and are more interested in the statistical prediction than in the exact locations of the extreme values.

Loss Function Calibrating the Main Dynamical Features.

In the training process, the model parameters are achieved through the optimization for the proper loss function (or cost) proposed depending on the target to be captured using the network. The two popular standard choices for the loss functions are using the error and the mean-square error: The error loss criterion measures the mean absolute error between each element in the output data and target through the distance,where is the training data size in one cycle, and is one member of output data and target data from the training minibatch. error loss criterion measures the mean-squared error between each element in the output and target through the mean-square norm,among all of the numbers of training samples .

The above two loss functions offer pointwise measurements of the errors in space-time for each predicted sample from the network. This may cause problems, especially when the system for prediction is highly turbulent, with internal instability, and is fast mixing. Small error perturbation in the input data may lead to vastly different solutions shortly after the mixing decorrelation time. The pointwise measurements focus on the accuracy at each value of the solutions; thus the small fluctuation errors might be accumulated and amplified under such metrics, and unnecessary large weights are added to correct errors in the moderate-amplitude fluctuation parts.

On the other hand, we are most interested in the prediction of statistical features in the extreme events rather than the exact trajectory solutions of the system. The small shifts in the extreme value locations should be tolerated in the loss function. Therefore, a more useful choice could be the relative entropy loss function that measures the Kullback–Leibler divergence in the predicted density functions: The relative entropy loss function computes the distance between two distribution functions,where the superscript represents the minibatch members to be measured in the relative entropy metric, and the subscript goes through all of the dimensions of the normalized variables to be described next.

Under the relative entropy loss function in Eq. , the input data and are treated as distribution functions. The “shapes” between the output data and the target are compared, rather than the pointwise details, so that the loss function guides the network to focus on the main model dynamical features instead of the turbulent fluctuations that are impossible to fit accurately. The additional difficulty in training the network using the relative entropy loss function is the constraint on the form of the output data to measure. The input of the relative entropy is required to be in the form of a density distribution function.

Scaling the Output Data with Empirical Partition Functions.

In this case of using the relative entropy loss function, we need to propose proper preprocessing of the output and target data to fit the required structure as a probability distribution. One direct way to do so is by taking the softmax function for the output data from the neural network normalized by a partition function,before measuring the error in the loss function. In this way, the data to put into the relative entropy loss function are normalized inside the range [0,1] with summation 1. This agrees with the definition in the relative entropy inputs. More importantly, this normalization emphasizes the large positive values of the data. Thus it offers a better calibration for the occurrence of positive major flow structures to be captured in the solutions.

Furthermore, a better choice for balancing both the positive and negative dominant values in the training data is to introduce scales with “temperatures.” We use the following two empirical partition functions with both positive and negative coefficients to rescale the output data:where are the positive and negative temperatures weighing the importance of dominant large-amplitude features in the scaled measure. Accordingly, the loss function to minimize under the relative entropy metric becomes a combination with the two empirical partition functions,where we use as a further balance between the positive and negative temperature components. In this combined empirical partition function metric using Eqs. and , the major flow structures in the turbulent field represented by the dominant extreme values are better characterized from both the positive and negative sides in the statistics of the model. In the following computational experiments, we always pick and for simplicity. More discussion of role of the balance weight is provided in .

Model Performance in Training and Prediction Stages

Using the previous model construction, we show the training and prediction performance using the MS-D Net combined with the relative entropy loss function and rescaled output data using empirical partition functions applied on the tKdV equations. In the numerical tests, we first consider the optimal prediction skill using the deep neural network within one updating cycle. Especially, we are more interested in capturing the non-Gaussian statistics from the network rather than the exact recovery of the single time series which should give large difference with small perturbations due to its turbulent nature. The basic strategy is to train the model using data from the near-Gaussian solutions with a small inverse temperature in the Gibbs measure [; then the optimized neural network is used to predict the more skewed model statistics for regimes with larger absolute values of .

For simplicity, we set the input and output data of the network in the same shape. Then the model output from one single iteration of the network is compared with the target data from the true model solution. In the structure of the neural network, the standard setup is adopted for the tKdV solutions. In each layer of the network, a kernel with minimum size is taken for the convolution update. The dilation size for mixed scales changes from to repeatedly as the network grows in depth. Different network depths of layers are tested for the model performance, while it is found that a moderate choice is enough to produce desirable training and prediction results.

In calibrating the errors from model predictions, we propose the normalized square error between the network output and the true target ,where the subscript represents the th component in the prediction/target set . The function acting on each component of can be used to extract the useful features to be calibrated. In the following tests, we use to compare the original output of the data, and use the exponential scaling to check the prediction in positive extreme values.

Training the Neural Network Using Near-Gaussian Data.

In the training process for turbulent system statistics, we include the temporal and spatial correlations together in the input data by considering a short time series of the solution. The training data are drawn from the model solutions of Eq. only among the near-Gaussian regime statistics. In summary, we use the training dataset in the following structure:

The input data are from the ensemble solutions of the tKdV equation. It is organized in a tensor of the shape , where is the total ensemble size, is the spatial discretization size, and is the sampled time instants with the time step . Thus the initial data for training are given as the tKdV solutions in the time window .

The target data for the training result are the solutions of the tKdV solution. The data are organized into the tensor with the same size as the input data. We can consider different starting times for the target data by changing the time length . The prediction time scale here is taken as , that is, to consider the prediction in time window .

The input data (with total ensemble size 10,000) are divided into 100 minibatches with size 100 for each training group in one epoch. In total, we use 1,000 epochs in the entire training process. The minibatches are randomly selected, with the batch indices resampled from random numbers in each step.

Notice that the above prediction time length used in the experiments is much longer than the decorrelation time of the tKdV states. From the direct numerical results in Fig. 1, the autocorrelation function decays to near zero at , and the longest decorrelation time among the spectral modes is below 0.1.

First, we compare the training performance using the standard and loss functions in Eqs. and with the same MS-D Net structure. Fig. 2 shows the evolution of training loss functions and the mean relative square errors [ among the training samples during the stochastic gradient descent iterations. According to the loss functions under the and distances, the training appears to be effective, and the error quickly dropped to smaller values in the first few steps. However, if we compare more carefully regarding the relative errors in the results, both cases get saturated quickly at a high error level near 1. The errors then become difficult to improve by training with a larger number of iterations and applying deeper layers. This is because, under both metrics, the model tries hard to fit the small-scale turbulent fluctuations in small values while missing the most important large-scale events in the solutions. It can be seen more clearly in Fig. 3 for typical training output snapshots compared with the truth. No desirable prediction can be reached.

Fig. 2.

Training loss function and the mean relative square error in the data (A) using and error loss functions and (B) using the relative entropy loss function with rescaled output data during the training iterations. Both networks with and loss are set to have densely connected layers; the networks with the relative entropy loss are compared using layers.

Fig. 3.

One snapshot of the final training results with three different loss functions. (Left) the input data, (Middle) the true target to fit, and (Right) the output results from the trained networks. All networks contain layers in the tests.

In contrast, significant improvement is achieved by switching to the relative entropy loss function and adopting empirical partition functions to normalize the output data. In this case as illustrated in Fig. 2, both the relative entropy loss function and the relative error drop to very small values in the final steps of the training iterations, implying high skill of the network in producing accurate predictions in the prediction time range (although it appears that the accuracy could be further improved by applying more iteration steps in the training, the results are already good enough after about 200 iterations). Note that we use logarithmic coordinates so the small values are emphasized. According to the last row of Fig. 3, for one typical training result snapshot, both the extreme values and the small-amplitude structures are captured in the model.

As a final remark, we check the proper depth needed for accurate predictions in the neural network. Fig. 2 also compares the same network under relative entropy loss but using different numbers of layers. A deeper network clearly can further improve the prediction skill and push the final optimized error to an even lower value, with the cost of a larger computational requirement. Still, from the comparison, it shows that a moderate number of layers (such as ) is sufficient to produce accurate results with relatively low cost. By pushing the network to deeper layers with , the improvement in error just becomes small and may not be necessary with the additional computational cost. The last row of Fig. 3 shows already the quite accurate recovery of the solution field purely learned from data.

Predicting Extreme Events Using Deep Neural Network.

In checking the prediction skill of the optimized network, we pick the neural network with densely connected layers as the standard model to test its predictions among different statistical regimes. It has been shown, in the training process, to have a high skill in recovering the original flow structures. Next, we should confirm that the neural network has really learned the dynamical structure of the original model, instead of merely overfitting the data.

Three statistical regimes ranging from near-Gaussian, mildly skewed, and highly skewed PDFs as shown in Fig. 1 are taken, for testing the range of prediction skill in the neural network model. An ensemble of 500 new trajectories from the tKdV solutions in different statistical regimes is used to show the robustness of the method. The mean and variance of the relative square errors [ among the samples for the state and the errors under the exponential function are listed in Table 1 for different statistical regimes. Uniform high accuracy in the mean with tiny variance is achieved among the vastly different regimes with distinct statistical features.

Table 1.

Mean and variance of the relative square errors among a test with 500 samples for the state and the scaled state using the same trained network for different statistical regimes

Error	Near-Gaussian	Mildly skewed	Highly skewed
u
Mean	0.2682	0.2556	0.2690
Variance	0.0039	0.0048	0.0087
exp (u)
Mean	0.0733	0.0764	0.0985
Variance	0.0005	0.0011	0.0060

Especially, we are interested in the case with highly non-Gaussian statistics representing the frequent occurrence of extreme events. In Fig. 4, the network is used to predict the flow solutions in the regime with highly skewed statistics (results for the other two cases can be found in ). The extreme values are represented by high peaks of a dominant wave moving along the field. This feature is not shown at all in the training data, where only near-Gaussian statistics is presented. As shown in the results, the trained network displays uniform skill among all of the tested samples in capturing the exact dynamical solutions in the extreme event domain unknown from the training data. By looking at the errors in the scaled data using the exponential function, the error amplitude even becomes smaller, confirming the accurate characterization of large extreme values through the network. In the typical snapshot of one sample, both the extreme values in the transporting waves and nonextreme detailed turbulent fluctuating structures are captured by the model.

Fig. 4.

Prediction in the regime with highly skewed statistics using the trained neural network with layers. (Top and Middle) The relative square errors for the state and the scaled state among the 500 tested samples. (Bottom) One typical snapshot of the prediction compared with the truth.

Furthermore, the prediction errors as the absolute difference between the truth and the model prediction in the three statistical regimes are displayed in Fig. 5. The network predictions maintain accuracy with small errors in a much longer time scale than the decorrelation time of the states. This gives final confirmation of the general high skill in the deep neural network for capturing key representative features in complex dynamical systems once the essential structures are learned from the training procedure.

Fig. 5.

The prediction error in the absolute difference between the truth and model output in the three tested regimes with different statistics.

Concluding Remarks

A strategy using a densely connected multiscale deep neural network with relative entropy loss function for calibrating rescaled output data is proposed for the prediction of extreme events and anomalous features from data. It needs to be noticed that the extreme events are often represented by highly skewed PDFs and have frequent occurrence in the turbulent field (3, 8, 12), in contrast to the other situation of isolated rare events which can be studied with machine learning models (9). The prediction skill of the optimized deep neural network is tested on the tKdV equation, where different Gibbs states create a wide range of statistics from near-Gaussian to highly skewed distributions. By adopting the densely connected and multiscale structures, the deep neural network is easy to train with a standard model setup and fewer model hyperparameters. Many further investigations about the optimized structure of the network, such as the minimum required number of layers as well as a cost function as a combination of the or loss function together with the relative entropy loss, for efficient and accurate prediction in general turbulent systems can be considered next, based on the basic framework.

Using training data only drawn from the near-Gaussian regime of the dynamical model, the deep neural network displays high skill in learning the essential dynamical structures of the complex system and provides uniformly accurate prediction among a wide range of different regimes with distinct statistics. The network also shows robustness among tests in a large ensemble of samples. The robust performance in the test model implies the potential of more general applications using the neural network framework for the prediction of extreme events and important statistical features in a wider group of more realistic high-dimensional turbulent dynamical systems.