Prediction in a driven-dissipative system displaying a continuous phase transition

TL;DR

This study investigates the predictability of event sizes in a critical system using neural networks, revealing that predictability diminishes near criticality and is limited to large, non-scaling events, implying challenges in forecasting earthquakes.

Contribution

It demonstrates the limitations of neural network-based prediction in critical systems, especially near phase transitions, and highlights the difficulty of forecasting earthquakes with Gutenberg-Richter scaling.

Findings

Prediction accuracy decreases near criticality.

Large, non-scaling events are more predictable.

Earthquake faults with Gutenberg-Richter scaling are hard to forecast.

Abstract

Prediction in complex systems at criticality is believed to be very difficult, if not impossible. Of particular interest is whether earthquakes, whose distribution follows a power law (Gutenberg-Richter) distribution, are in principle unpredictable. We study the predictability of event sizes in the Olmai-Feder-Christensen model at different proximities to criticality using a convolutional neural network. The distribution of event sizes satisfies a power law with a cutoff for large events. We find that prediction decreases as criticality is approached and that prediction is possible only for large, non-scaling events. Our results suggest that earthquake faults that satisfy Gutenberg-Richter scaling are difficult to forecast.

Equations

$$n_s\sim s^{-\tau }\exp (-(s/s_c{)}^{\sigma })$$

$$s_c=(\eta -{\eta }_c{)}^{-1/\sigma },$$

Full text

I Introduction

A subclass of driven-dissipative systems modeled by a two-dimensional cellular automaton has been proposed to understand the existence of power laws in many complex systems. Examples include the Bak-Tang-Wiesenfeld sandpile model BTW_1987 , the Rundle-Jackson model Rundle_1977 and the Olami-Feder-Christensen (OFC) model OFC_1992 ; the latter two models have been used to gain insight into the nature of earthquakes. Many earthquake fault systems display a power-law event size distribution spanning many orders of magnitude. Such a power-law distribution is known as Gutenberg-Richter scaling in the seismology literature GR_1944 . For example, the Gutenberg-Richter scaling in Southern California (1984–2000) spans about six orders of magnitude Bak_2002 .

There has been a substantial interest in forecasting or predicting earthquakes. However, it has been conjectured that systems at criticality are inherently unpredictable BTW_1987 . That is, events of different sizes that satisfy a scale-free distribution are due to the same physical mechanism, and thus there are no distinct precursors to distinguish one event from another Bak_1989 . The idea of unpredictability at criticality has been challenged over the years. One school of thought Laherrere_1999 has proposed that very large events are due to inherently different mechanisms such as self-reinforcement, synchronization Sornette_2009 and nucleation Klein_2000 , and are thus in principle distinguishable from smaller events. Some support for this proposal is the use of a technique called the log-periodic-power-law fitting procedure; it has been shown to successfully predict large, non-power law events, such as ruptures in materials Anifrani_1995 ; Johansen_2000 and the end of financial bubbles Feigenbaum_1996 ; Sornette_1996 .

In this paper we address the question of predictability near and at criticality by applying machine learning to the OFC model. Previous work Pepke_1994 has shown that predictability in the OFC model decreases as the conservative limit (a critical point) is approached. We find results consistent with Ref. Pepke_1994, and investigate the predictability of events near another critical point in the OFC model: the recently observed noise transition critical point Matin_draft . By using a convolutional neural network (CNN), we find that the event sizes are more difficult to forecast as the critical point is approached and that only large events that do not satisfy power law scaling can be successfully predicted.

II Criticality in the OFC model

The Olami-Feder-Christensen (OFC) model OFC_1992 is a modified version of the spring-block model first proposed by Rundle and Jackson Rundle_1977 , which is a simplification of the Burridge-Knopoff model Burridge_1967 . The nearest-neighbor OFC model that we will consider consists of a two-dimensional lattice of linear dimension $L$ with each site initially assigned a random value of stress $\sigma$ between the mean residual $\sigma_{R}$ and the failure threshold $\sigma_{F}$. We denote the stress on each site, the stress grid, by the vector $\vec{\sigma}=(\sigma_{1},\ldots,\sigma_{N=L\times L})$. The system is then driven so that one site reaches $\sigma_{F}$, a procedure known as the zero velocity limit Rundle_1977 . This site is said to fail, and its stress is reduced to $\sigma_{R}\pm\eta r$, where $\eta$ is the magnitude of the noise and $r$ is a uniform random variable in the range $[-1,1]$. A failing site with stress $\sigma$ distributes stress $(1-\alpha)(\sigma-\sigma_{R}\mp\eta r)/4$ to its four nearest-neighbor sites, where $\alpha$ is the dissipation parameter. The failure of one site can trigger other sites to fail, thus creating an avalanche. The avalanche or event stops when the stress of all sites is less than $\sigma_{F}$. We denote the number of failing sites, or the size of the event, by $s$. The system is then driven again using the zero velocity limit.

The OFC model is believed to approach criticality in the conservative limit $\alpha\to 0$ Christensen_1991 ; Serino_2001 . Recently, it has been found that even for $\alpha>0$, there exists a phase transition at a critical value of the noise $\eta_{c}\approx 0.07$ Matin_draft . This phase transition is characterized by the event size distribution $n_{s}$ of the form

[TABLE]

with

[TABLE]

and $\tau=1.04\pm 0.14$ and $\sigma=0.43\pm 0.03$. The mean cluster size $\chi$ diverges as $(\eta-\eta_{c})_{+}^{-\gamma}$ and the connectedness length diverges as $(\eta-\eta_{c})_{+}^{-\nu}$ Stauffer_1979 with $\gamma=2.01\pm 0.14$ and $\nu=1.20\pm 0.13$, consistent with the scaling relations $\gamma=(3-\tau)/\sigma$ and $\nu=(\tau-1)/d\sigma$, where $d=2$ is the spatial dimension. Note that $n_{s}$ satisfies power law scaling for $s\lesssim s_{c}$ Matin_draft .

III Supervised machine learning

Our goal is to predict the event size (the number of failed sites) given the stress grid before stress has been added using the zero velocity limit and before the onset of an event. Figure 1 shows that the event size $s$ is strongly correlated with the average stress for $\eta<\eta_{c}$ and weakly correlated for $\eta\geq\eta_{c}$. The event size is weakly correlated with the variance of the stress for all values of $\eta$.

To force the CNN to learn higher order features, we first remove the correlations between the event size $s$ and the first and second moments of the stress grid. We normalize each stress grid $\vec{\sigma}$ by its average stress and the spatial variance. That is, we rescale the stress $\sigma_{i,\mu}$ at site $i$ for sample $\mu$ to $\tilde{\sigma}_{i,\mu}\equiv(\sigma_{i,\mu}-\langle{\sigma_{\mu}}\rangle)/\sqrt{\mbox{var}_{\mu}}$, where $\langle\sigma_{\mu}\rangle=1/N\sum_{i=1}^{N}\sigma_{i,\mu}$ is the mean stress per site of sample $\mu$ and $\mbox{var}_{\mu}\equiv\sum_{i=1}^{N}(\sigma_{i,\mu}-\langle\sigma_{\mu}\rangle)^{2}/N$ is the spatial variance of the stress. In the following all references to the stress will be to the rescaled stress and we will omit the tilde symbol. We will train the CNN regressor using the rescaled stress grid $\vec{\sigma}$, sampled with quasi–uniform event sizes (see the appendix).

To assess the performance of the machine, we show in Fig. 2 the predicted event size $\hat{s}$ versus the true event size $s$. The top row shows the event size distribution $n_{s}$ at different values of $\eta$. The bottom row shows the predicted event sizes versus the true event sizes. We see that for $\eta<\eta_{c}$, the machine performs impressively at predicting events that are larger than $s_{c}$ [see Fig. 2(a)]. For $\eta>\eta_{c}$, the machine performs less impressively [see Fig. 2(c)]. At $\eta=\eta_{c}$, the machine fails at predicting events of all sizes [Fig. 2(b)].

In Fig. 3 we plot the testing error $\text{Err}(\log\hat{s},\,\log s)\equiv 1/N\sqrt{\sum_{i=1}^{N}(\log s_{i}-\log\hat{s}_{i})^{2}}$ as a function of $\eta$. The reason for using $\log s$ instead of $s$ in the error function is because of the larger fluctuations in $(\hat{s}-s)$ for larger $s$ and because we are interested in the relative error rather than the absolute error. The peak at $\eta_{c}$ indicates that prediction is not possible in the OFC model at criticality using only the stress grid and the CNN.

We next look at how the values of the dissipation parameter $\alpha$ affect the predictability of the system as $\alpha\to 0$. No scaling function has been found to fit the dependence of $n_{s}$ on $\alpha$ in the nearest-neighbor OFC model. Nevertheless, we can determine the cutoff $s_{c,\alpha}$ from $n_{s}$ as the value of $s$ for which $n_{s}$ deviates from a power law (see Fig. 4). As $\alpha$ decreases, the cutoff $s_{c,\alpha}$ increases. We observe that the onset of predictability is close to $s_{c,\alpha}$ and the trend persists for different values of $\alpha$.

IV Visualizing the CNN

We next explore the features that the machine has learned which allow it to successfully forecast the size of the non-scaling events and discuss why the critical events are difficult to forecast. We will use occlusion sensitivity analysis to identify the regions of importance of the images that are used by the CNN Zeiler_2014 . For example, the face of a dog is expected to contain the most relevant features in determining the type of animal. Hence, blocking the face of the dog should increase the classification error of the CNN. We implement a similar analysis by defining an occluded region in the stress grid and sweeping the occluded region across the entire image to create a map that shows the regions that are the most sensitive to the occlusion. In this way we associate the region that gives the largest change in the predicted event size $\hat{s}$ with the region that is most useful in determining the size of the event.

In Fig. 5 we visualize three randomly chosen samples for which $s>s_{c}$ for (a) $\eta<\eta_{c}$, (b) $\eta=\eta_{c}$, and (c) $\eta>\eta_{c}$. In the top row we show the failure maps, which correspond to the number of times that a site has failed. In the second row we show the sensitivity maps from the occlusion sensitivity analysis. Because we chose events of size $s>s_{c}$, an occlusion that yields a decrease in the predicted event size $\hat{s}$ implies a worse prediction. For $\eta\neq\eta_{c}$, the region that gives the largest decrease in $\hat{s}$ if occluded coincides with the failure region. We call the region with the largest increase in $\hat{s}$ if occluded the sensitive region. In the third row, we plot the local average stress map. To determine this map, the local average stress of a site is computed by averaging the stress of sites within a square of linear dimension $b=10$, centered at that site. Note that the region of high local stress overlaps with the failure and most sensitive regions. This consistency is reasonable because regions with high local stress have a greater probability of initiating and sustaining a large event.

Among the 32 channels in the third layer of the CNN we chose the channel that is visually the most similar to the structure of the failure region. We interpret the channel as the high level feature learned by the CNN. We plot the channels in the bottom row of Fig. 5. From these channels, we see that the machine has learned the connection between the high local stress and failure regions.

To understand why prediction is difficult at $\eta=\eta_{c}$, we look at the failure maps in the top row of Fig. 5. We see that the failure regions become more diffuse at $\eta_{c}$ compared to the more compact failure regions away from $\eta_{c}$. Although the local average stress map and the failure map remain qualitatively similar, the stress gradient between the high local stress region and the surrounding background is much smaller at $\eta=\eta_{c}$.

More quantitatively, we define a high stress region as a collection of nearest-neighbor connected sites whose local stress is above the cutoff $\overline{\sigma}_{c}$ (the bar denotes the stress is the local average stress). We measure the radius of gyration $R_{g}$ of the largest high stress region in each sample and define the density of the high stress region $\tilde{\phi}$ as the sum of the local stress within the area of radius $R_{g}$ divided by $\pi R^{2}_{g}$. The density of the high stress region $\tilde{\phi}$ decreases as $\eta_{c}$ is approached (see Fig. 6). The smaller density difference makes it more difficult for the machine to obtain the appropriate cutoff for the high stress region, thus making prediction more difficult. Multiple failures are more prominent for very large events for $\eta\neq\eta_{c}$, which is why the machine underestimates the event sizes of the very large events [Figs. 2(a) and (c)].

V Discussion

Since we have normalized the stress grid by the average stress before training the machine, the first and second moments of the stress grids do not contain information that can be used by the CNN to predict event sizes. The fact that the machine learns the association between the region of high local average stress and the event size means that the machine has learned the optimal cutoff that separates the high stress regions from the low stress regions. This task becomes increasingly difficult as the critical point is approached because the high stress region becomes less distinguishable from its background.

We have found evidence that events whose size distribution satisfies a power law lack distinguishable features which allow the machine to predict their size. This lack of distinguishable features is related to the difficulty of distinguishing between the fluctuations and the background at critical points Coniglio_1980 . For the large non-scaling events, there exists features that allow the machine to successfully predict the event sizes. Similar conclusions are found for both the noise Matin_draft and dissipation Christensen_1991 transitions. Our results suggest that large non-scaling events are qualitatively different from the smaller scaling events. This conclusion agrees with the conjecture Bak_1989 that prediction is not possible at a true critical point ($L\to\infty$), where there is no deviation from a power law for large events.

It is known that small, large, and very large events in the long-range Rundle-Jackson model Rundle_1977 are due to different mechanisms, namely, fluctuations about the spinodal critical point, and failed nucleation and arrested nucleation events, respectively Klein_2000 ; Anghel_2000 ; Xia . These different mechanisms suggest that very large events are in principle distinguishable from other events. The caveat is that all three types of events follow a power law, albeit with different exponents. It will be interesting to see if a machine can learn the difference between the different scaling events. It is important to note that the failed and arrested nucleation events, despite the fact that they satisfy a power-law distribution, do not exhibit the same diffusive nature as the smaller events (spinodal fluctuations) on the Gutenburg-Richter scaling plot. This difference appears to be what the CNN picks up.

Acknowledgement

We would like to thank Huang Shan and Pankaj Mehta of the Department of Physics at Boston University for useful discussions.

Appendix A Sampling method and CNN architecture

After discarding the transient ($10^{6}$ plate updates), we run for an additional $10^{7}$ plate updates and record the event sizes and the random number seed. We then construct the event size distribution and randomly chose five samples from each bin of the distribution (or the number of samples in that bin if there are less than five samples) and record the time of events in each bin. We then re-run the simulation using the same random number seed and save the stress grids at the recorded times. This procedure ensures that the number of samples in each bin remains the same for different values of the noise.

The architecture of the CNN consists of 8 alternating layers of convolutional layers and maxpooling layers (4 layers each). The depths of the convolutional layers are 8, 16, 32 and 64, each with a filter of size 5 $\times 5$. We used zero padding on the boundaries to ensure the same size after each convolution. The output of the last maxpooling layer is connected to a fully connected neural network with one hidden layer of 25 nodes. All layers use relu (rectified linear unit) as the activation function except the last layer, which uses a linear activation function. Dropout Dropout_2014 with $\mbox{dropout rate}=0.1$ is applied to the layer immediate before the fully connected layer.