Machine Learning Percolation Model

In this section, we will focus on the two-dimensional percolation model and the dataset obtained by the Monte Carlo simulation. This section will discuss several machine learning approaches, including CNNs, VAE and cVAE, and PCA, to deduce physical characteristics.

Let us first introduce CNNs. CNNs, supervised learning methods, are particularly useful in solving realistic problem for many disciplines, such as physics[3], chemistry [19], medicine [20], economics [21], biology [22], and geophysics [23, 24], ect. In the applications of phase transitions and critical phenomena, many predecessors utilize CNNs to detect physical features, especially order parameter [1, 25, 26, 27, 10, 5, 2, 6]. In this study, the four kinds of CNNs are not only used to detect the two order parameters ($\varPi(\bm{p},L)$ and $P(\bm{p},L)$), but also to detect the permeability $\bm{p}$.

Next, we demonstrate the architecture of the CNNs (see Fig. 2). The structure of the CNNs has four layers, including two convolution layers and two fully connected layers. The percolating configurations $\{\bm{X}_{1},\bm{X}_{2},\ldots,\bm{X}_{M}\}$ are taken as inputs. Consequently, the CNNs receive the corresponding order parameters ($\varPi(\bm{p},L)$ and $P(\bm{p},L)$) or the permeability $\bm{p}$ as outputs.

In the former section (see Sec. III.1), the raw configuration $\bm{X}$ at different permeability $\bm{p}$ is generated by Monte Carlo simulation. However, what if configurations are not in the raw configurations, can we still identify the physical features as well? Here we consider to use the VAE and the cVAE to generate new configuration $\hat{\bm{X}}_{\text{VAE}}$ and $\hat{\bm{X}}_{\text{cVAE}}$, respectively.

VAE (see Fig. 3), a generative network, bases on the variational Bayes inference proposed by [28]. Contracted with traditional AE (see Fig. S. 1), the VAE describes latent variables with probability. From that point, the VAE shows great values in data generation. Just like AE, the VAE is composed of an encoder and a decoder. The VAE uses two different CNNs as two probability density distributions. The encoder in the VAE, called the inference network $p_{\text{encoder}}(\bm{Z}|\bm{X})$, can generate the latent variables $\bm{Z}$. And the decoder in the VAE, called the generating network $p_{\text{decoder}}(\hat{\bm{X}}|\bm{Z})$, reconstructs the raw configuration $\bm{X}$. Unlike AE, the encoder and the decoder in VAE are constrained by the two probability density distributions.

To better reconstruct the raw two-dimensional configuration $\bm{X}$, the fully connected layer in the VAE is replaced by the convolution layer. Now we have got the cVAE. The architecture of the cVAE is shown in Fig. 3. Usually, the performance of the cVAE is better than the VAE due to the configuration $\bm{X}$ with spatial attribute. In our work, the VAE and the cVAE are both used for generating THE new configuration $\hat{\bm{X}}$.

The Sec. III.1 and Sec. III.2 focus on supervised learning, which hypothesize that the labels exist for the raw configuration $\bm{X}$ on the percolation model. However, though we can detect the permeability values $\{0.41,0.42,\ldots,0.80\}$ and the two order parameters ($\varPi(\bm{p},L)$ and $P(\bm{p},L)$) by supervised learning, label dearth often occurs. Thus, it is imperative to identify the labels, such as $\varPi(\bm{p},L)$, $P(\bm{p},L)$, and $\bm{p}$. Some recent studies have shown that the first principal component obtained by PCA can be regarded as typical physical quantities [15, 8, 14, 17]. Base on these studies, we explore the meaning of the first principal component on the percolation model.

As it is well-known, PCA can reduce the dimension of the matrix $\bm{X}$. First, we compute the mean value $\bm{X}_{\text{mean}}=1/M\sum_{i=1}^{M}a_{ij}(i=1,2,\ldots,M;j=1,2,\ldots,N)$ for each column $\bm{C}_{j}(j=1,2,\ldots,N)$ in the matrix $\bm{X}$. Then we get the centered matrix $\bm{X}_{\text{centered}}$ that is expressed as $\bm{X}_{\text{centered}}=\bm{X}-\bm{X}_{\text{mean}}$. After obtaining $\bm{X}_{\text{centered}}$, by an orthogonal linear transformation expressed as $\bm{X}_{\text{centered}}^{T}\bm{X}_{\text{centered}}\bm{W}=\bm{\lambda}\bm{W}$, we extract the eigenvectors $\bm{W}$ and the eigenvalues $\bm{\lambda}$. The eigenvectors $\bm{W}$ are composed with $\bm{w}_{1},\bm{w}_{2},\ldots,\bm{w}_{N}$. The eigenvalues are sorted in the descending order, i.e., $\bm{\lambda}_{1}\geq\bm{\lambda}_{2}\geq\ldots\geq\bm{\lambda}_{N}\geq 0$. The normalized eigenvalues $\tilde{\bm{\lambda}}_{j}(j=1,2,\ldots,N)$ are expressed as $\bm{\lambda}_{j}/\sum_{j=1}^{N}{\bm{\lambda}_{j}}$. The row $\bm{R}_{i}$ in matrix $\bm{X}$ can be transformed into $\bm{X}^{{}^{\prime}}_{i}=\bm{R}_{i}\bm{W}$. Eq. 13 represents the statistic average every 40 intervals for each permeability $p$. This process is quite similar to the process of calculating the two order parameters, i.e., $\varPi(\bm{p},L)$ and $P(\bm{p},L)$. Table 1 shows the procedure of PCA algorithm.

In this section, we consider to use the approches in Sec. III to capture the physic features. First we make use of TensorFlow 2.2 library to perform the CNNs with four layers. To predict the two order parameters ($\varPi(\bm{p},L)$ and $P(\bm{p},L)$), two kinds of CNNs (CNNs-I and CNNs-II) are constructed. The first two layers of the former two kinds of CNNs are composed of two convolution layers (“Con1” and “Con2”), both of which possessed “filter1”=32 and “filter2”=64 filters with the size of $3\times 3$, and a stride of 1. Each convolution layer is followed by a max-pooling layer with the size of $2\times 2$. The final convolution layer “Con2” is strongly interlinked to a fully connected layer “FC” with 128 variables. The output layer “Output”, following by “FC”, is a fully connected layer. For the two convolution layers and the fully connected layer “FC”, a rectified linear unit (ReLU) $\bm{a}=\text{max}(0,\bm{x})$ [29] is chosen as activation function due to its reliability and validity. However, the output layer has no activation function.

After determining the framework of CNNs-I and CNNs-II, here we mention how to train CNNs-I and CNNs-II for deducing $\varPi(\bm{p},L)$ and $P(\bm{p},L)$. First, we carry out an Adam algorithm [30] as an optimizer to update parameters, i.e., weights and biases. Then, a mini batch size of 256 and a learning rate of $10^{-4}$ are selected for its timesaving. Following this treatment, CNNs-I and CNNs-II are trained on 1000 epochs for 40,000 uncorrelated and shuffled configurations, respectively. Before training, we split $\{\bm{X}_{1},\bm{X}_{2},\ldots,\bm{X}_{M}\}$, $\varPi(\bm{p},L)$ and $P(\bm{p},L)$ into 32,000 training set and 8,000 testing set. While training, we monitor three indicators, including the loss function (i.e., mean squared error (MSE, see Eq. 14)), mean average error (MAE, see Eq. 15), and root mean squared error (RMSE, see Eq. 16), for training and testing set [24]. In Eq. 14-16, ${y_{i}}^{\text{raw}}$ and ${y_{i}}^{\text{pred}}$ refer to $\varPi(\bm{p},L)$/$P(\bm{p},L)$ and its predictions. If the loss function in testing set reaches the minimum, then the optimal CNNs-I and CNNs-II will be obtained. As can be seen from Fig. S. 2, these indicators gradually decrease. In Table. S. 3, the errors of the optimal CNNs-I and CNNs-II are very small. What stands out in Fig. S. 2 is that CNNs-I and CNNs-II have high stability, consistency, and faster convergence rate.

Before assess CNNs-I and CNNs-II, we have to explain what is meant by statistic average. Statistic average can be defined as the averages of CNNs-I’s or CNNs-II’s outputs for each permeability $p$. As shown in Fig. 5, there is a clear trend of phase transition between the permeability values $\{0.41,0.42,\ldots,0.80\}$ and $\varPi(\bm{p},L)$/$P(\bm{p},L)$. The two grey lines in Fig. 5, which are the same as the two blue lines in Fig. 1, represent the relationship between the permeability values $\{0.41,0.42,\ldots,0.80\}$ and the raw $\varPi(\bm{p},L)$ or $P(\bm{p},L)$. Likewise, the two blue lines in Fig. 5, represent the relationship between the permeability values $\{0.41,0.42,\ldots,0.80\}$ and the statistic average from the outputs of CNNs-I and CNNs-II. The overlapping of the two kinds of lines shows that CNNs-I and CNNs-II can be used to deduce the two order parameters and the process of phase transition.

To overcome the difficulty associated with the percolation model being near the critical transition point, we truncate the dataset. Specifically, we remove the data near the critical transition point, and only retain the data far away from the critical transition point. Here we take the simulation of $\varPi(\bm{p},L)$ as an example. The retained data with the raw $\varPi(\bm{p},L)$ rangs from 0 to 0.1, and 0.9 to 1. As shown in Fig. S. 3 and the middle red points in Fig. 6, we find that CNNs-I can extrapolate $\varPi(\bm{p},L)$ to missing data by learning the retained data.

This section demonstrate that CNNs-I and CNNs-II can be two effective tools for detecting $\varPi(\bm{p},L)$ and $P(\bm{p},L)$, respectively. Additional test should be made to verify that whether or not CNNs-I and CNNs-II are robust against noise. To address this issue, we deliberately invert a proportion, i.e., 5%, 10%, and 20%, of the labels for the raw $\varPi(\bm{p},L)$ and $P(\bm{p},L)$ and verify that whether or not the “artificial” noises can affect the predicted $\varPi(\bm{p},L)$ and $P(\bm{p},L)$. Fig. 7 and Fig. S. 4 demonstrate that CNNs-I and CNNs-II are robust against noise. As the labeling error rates increase, the same trend is evident in the outputs of CNNs-I and CNNs-II within a relatively small difference (see Fig. 7). Therefore, we draw the conclusion that noises have little effect on detecting $\varPi(\bm{p},L)$ and $P(\bm{p},L)$.

Just like we use CNNs-I and CNNs-II in Sec. IV.1, here we use the same structure for CNNs-III and strategies for training the permeability $\bm{p}$. The only distinction between CNNs-III and CNNs-I/CNNs-II is that the outputs for CNNs-III are the permeability $\bm{p}$ instead of the two order parameters. Fig. S. 5 shows the performance of CNNs-III. With successive increases in epochs, the MSE, MAE, and RMSE continue to decrease until no longer dropping.

Another measure of CNNs-III’s performance is concerned with the difference between the raw permeability $\bm{p}$ and its prediction $\hat{\bm{p}}$. The blue circles in Fig. 8 show that there is a strong positive correlation between the raw permeability $\bm{p}$ and its prediction $\hat{\bm{p}}$. Further statistical tests reveal that most of the gap between the raw permeability $\bm{p}$ and its prediction $\hat{\bm{p}}$ is less than 0.1. The result proves that CNNs-III has the advantage of convenient use and high precision.

On the other hand, extrapolation ability can also reflect the performance of CNNs-III. To do this, we use Monte Carlo simulation to generate new dataset. The new permeability $\bm{p}$ not only ranges from 0.01 to 0.40, but also from 0.81 to 1.00, with an interval of 0.01. Just like Sec. II, we perform 1050 Monte Carlo steps and keep the last 1000 steps. As a consequence, 60,000 configurations are generated. The red circles in Fig. 8 exhibit the extrapolation ability of CNNs-III. Our results show that no significant correlation between the extrapolated permeability $\bm{p}_{\text{ex}}$ ranging from 0.81 to 1.00 and their prediction $\hat{\bm{p}}_{\text{extrapolated}}$. However, in Fig. 8, the results indicate that CNNs-III has good extrapolation ability for the extrapolated permeability $\bm{p}_{\text{ex}}$ from 0.01 to 0.40.

Though CNNs-I, CNNs-II, and CNNs-III are valid when detecting the two order parameters ($\varPi(\bm{p},L)$ and $P(\bm{p},L)$) and the permeability $\bm{p}$, the validity of these three CNNs is unkown for percolating configurations outside of the dataset. As is shown in Fig. 3 and Fig. 4, we use the same network structures for the VAE and the cVAE to generate new configurations ($\hat{\bm{X}}_{\text{vae}}$ and $\hat{\bm{X}}_{\text{cvae}}$). Actually, $\hat{\bm{X}}_{\text{vae}}$ and $\hat{\bm{X}}_{\text{cvae}}$ can be regarded as adding some noise into the raw configuration $\bm{X}$.

Let us first consider VAE. Just like AE (see Fig. S. 1), the VAE is also composed of an encoder and a decoder. The encoder of the VAE owns two fully connected layers, both of which follows with a ReLU activation function. The first layer of the encoder possess “size1”=512 neurons. Another 512 neurons, including 256 mean “$\bm{\mu}$” and 256 variance “$\bm{\sigma}$”, are taken into account for the second layer of the encoder. By resampling from the Gaussian distribution with the mean “$\bm{\mu}$” and the variance “$\bm{\sigma}$”, we obtain 256 latent variables “$\bm{Z}$” which are the inputs of the decoder. For the decoder of the VAE, two fully connected layers follow with the outputs of the encoder. For symmetry, the first layer of the decoder also contains “size1”=512 neurons and follows with a ReLU activation function. And the output layer of the decoder contains 784 neurons which are used to reconstruct the raw configuration $\bm{X}$. Thus, the neurons in the output layer are the same as that in the input layer.

Moving on now to consider cVAE. The encoder of the cVAE is composed of one input layer, two hidden convolution layers with “filter1”=32 and “filter2”=64 filters with the size of 3 and a stride of 2, and a ReLU activation function. The output layer in the encoder is a fully connected flatten layer with 800 neurons (400 mean “$\bm{\mu}$” and 400 variance “$\bm{\sigma}$”) without activation function. By resampling from the Gaussian distribution with the mean “$\bm{\mu}$” and the variance “$\bm{\sigma}$”, we obtain 400 latent variables “$\bm{Z}$”. The reason why latent variables in the cVAE is more than the VAE is that the cVAE needs to consider more complex spatial characteristic. The decoder of the cVAE is composed of an input layer with 400 latent variables “$\bm{Z}$”. A fully connected layer “FC” with 1,568 neurons is followed by “$\bm{Z}$”. After reshaping the outputs of “FC” into three dimension, we feed the data into two transposed convolution layers (“Conv3” and “Conv4”) and one output layer “Output”. The filters in these deconvolution layers are 64, 32 and 1 with the size of 3 and the stride of 2. After excluding the output layer “Output” without activation functions, there exist two ReLU activation functions followed by “Conv3” and “Conv4”, respectively.

We train the VAE and the cVAE over $10^{3}$ epochs using the Adam optimizer, a learning rate of $10^{-3}$, and a mini batch size of 256. To train the VAE and the cVAE, we use the sum of binary cross-entropy (see Eq. 17) and the Kullback-Leibler (KL) divergence (see Eq. 18) as the loss function [18]. In Eq. 17-18, $\bm{x}_{i}^{\text{raw}}$ and $\bm{x}_{i}^{\text{pred}}$ represent each raw configuration with one/two dimension and its prediction. As shown in Fig. S. 6, the loss function, the binary cross-entropy, and the KL divergence vary with epochs for the VAE and the cVAE. Here we focus on the minimum value of loss function. From Fig. S. 6, the optimal cVAE performs better than the optimal VAE.

For a more visual comparison, we show the snapshots of the raw configurations $\{\bm{X}_{1},\bm{X}_{2},\ldots,\bm{X}_{M}\}$, and compare them to the VAE-generated and cVAE-generated configurations in Fig. 9. As we can see from Fig. 9, the configurations from the Monte Carlo simulation, the VAE and the cVAE are very close to each other.

After reconstructing the configurations ($\hat{\bm{X}}_{\text{vae}}$ and $\hat{\bm{X}}_{\text{cvae}}$) through the VAE and the cVAE and pouring $\hat{\bm{X}}_{\text{vae}}$ and $\hat{\bm{X}}_{\text{cvae}}$ into CNNs-I, CNNs-II, and CNNs-III, we can detect $\varPi(\bm{p},L)$, $P(\bm{p},L)$, and the permeability $\bm{p}$. Fig. 10 shows the relationships between the permeability values $\{0.41,0.42,\ldots,0.80\}$ and the statistic average from the outputs in CNNs-I and CNNs-II, respectively. From the red and purple lines in Fig. 10, by using $\hat{\bm{X}}_{\text{vae}}$ and $\hat{\bm{X}}_{\text{cvae}}$, we can obtain the two order parameters as well. From the Fig. 11, the raw permeability $\bm{p}$ is remarkably correlated linearly with its prediction $\hat{\bm{p}}$ through the VAE/cVAE and CNNs-III. Thus, subtle change in the raw configurations does not effect the catch of physical features.

This section will discuss how to capture the physic characteristics without labels. Various studies suggest to use PCA to identify order parameters [15, 8, 14, 17]. Therefore, we try to verify the feasibility of PCA in capturing order parameter on the percolation model.

First, we perform the PCA to the raw configuration $\bm{X}$. Fig. 12 exhibit the $N$ normalized eigenvalues $\tilde{\bm{\lambda}}_{n}=\bm{\lambda}_{n}/\sum_{n=1}^{N}\bm{\lambda}_{n}$. $\tilde{\bm{\lambda}}_{n}$ is also called as the explained variance ratios. The most noteworthy information in Fig. 12 is that there is one dominant principal component $\tilde{\bm{\lambda}}_{1}$, whcih is the largest one among $\tilde{\bm{\lambda}}_{n}$ and much larger than other explained variance ratios. Thus, $\tilde{\bm{\lambda}}_{1}$ plays a key role when dealing with dimension reduction. Based on $\tilde{\bm{\lambda}}_{1}$, the raw configuration $\bm{X}$ are mapped to another matrix $\bm{Y}=\bm{X}\tilde{\bm{\lambda}}_{1}$.

In Fig 13, we construct the matrix $\bm{Y}^{{}^{\prime}}=\{{\bm{X}\tilde{\bm{\lambda}}_{1},\bm{X}\tilde{\bm{\lambda}}_{2}}\}$ by the first two eigenvalues and their eigenvectors. We use 40,000 blue scatter points to plot the the relationship between $\bm{X}\tilde{\bm{\lambda}}_{1}$ and $\bm{X}\tilde{\bm{\lambda}}_{2}$ on 40 permeability values ranging from 0.41 to 0.8. Just like in Fig. 12, there is only one dominant representation on the percolation model due to the first principal component is much more important than the second principal component.

Having analysed the importance of the first principal component, we now move on to discuss the meaning of the first principal component. In Fig. 14, we focus on the quantified first principal component as a function of the permeability $\bm{p}$ and the two order parameters, i.e., $\varPi(\bm{p},L)$ and $P(\bm{p},L)$. From Fig. 14, we can see that there is a strong linear correlation between the quantified first principal component and the permeability $\bm{p}$. And the relationships between the quantified first principal component and two order parameters ($\varPi(\bm{p},L)$ and $P(\bm{p},L)$) are similar to the relationships between the permeability values $\{0.41,0.42,\ldots,0.80\}$ and the two order parameters. Our results are significant different from former study which demonstrates that the quantified first principal component can be taken as order parameter by data preprocessing [17]. A possible explanation may be that [17] evaluates the first principal component by removing the raw dataset with certain attributes on the percolation model. Therefore, we assume that the quantified first principal component obtained by PCA may not be taken as order parameter for various physical models. Another possible explanation is that, under certain conditions, the quantified first principal component can be regarded as order parameter.

From Sec. IV.4, no significant corresponding is found between the normalized quantified first principal component and the two order parameters. Here we eagerly wonder how to identify order parameter. And another physical characteristics we desire to explore is the critical transition point. So we try to find a way to capture order parameter from the raw configurations $\{\bm{X}_{1},\bm{X}_{2},\ldots,\bm{X}_{M}\}$ and their permeability $\bm{p}$ on the percolation model.

As it is well-known, for the two-dimensional percolating configurations, the critical transition point is equal to 0.593 in theoretical calculation. Though the critical transition point is already known, we wonder that whether or not the critical transition point can be found by machine. To do this, we first use a cluster analysis algorithm named $k$-means [31] to separate the raw configurations $\{\bm{X}_{1},\bm{X}_{2},\ldots,\bm{X}_{M}\}$ into two categories. The minimum and maximum value of their permeability $\bm{p}$ are 0.55 and 0.80 in the first cluster, and 0.41 and 0.65 in the second cluster. Note that there is an overlapping interval between 0.55 and 0.65 for the two categories. According to the overlapping interval, we hypothesize that the critical threshold $p_{c}$ is set to be 11 values, i.e., 0.55, 0.56, $\ldots$, and 0.65. For each raw configuration, if the permeability is smaller than $p_{c}$, there will exist no percolating cluster and its label will be marked as 0; otherwise, there will exist at least one percolating cluster and its label will be marked as 1.

To detect physic characteristics, we use the fourth kinds of CNNs (CNNs-IV) with the structure in Fig. 2 for 11 preset critical thresholds from 0.55 to 0.65. Note that the output layer use a sigmoid activation function expressed as $\bm{a}=1/(1+e^{-\bm{x}})$, to make sure the outputs are between 0 and 1. Another critical thing to pay attention is that the He normal distribution initializer [32] and L2 regularization [33] are used in the layer of “Conv1”, “Conv2”, and “FC” on CNNs-IV. To avoid overfiting, in addition to L2 regularization, we also use a dropout layer with a dropout rate of 0.5 on “FC”. A Mini batch size of 512 and a learning rate of $10^{-4}$ are chosen while training CNNs-IV. The binary cross-entropy (see Eq. 17) is taken as the loss function on CNNs-IV. Another metric, used to measure the performance of CNNs-IV, is the binary accuracy (see Eq. 19). The other hyper-parameters are the same as the CNNs-I, CNNs-II, and CNNs-III.

Turning now to the experimental evidence on the inference ability of capturing relevant physic features. After obtaining the well-trained CNNs-IV with high accuracy (see Fig. S. 7), we obtain the outputs by pouring the raw 40,000 configurations $\{\bm{X}_{1},\bm{X}_{2},\ldots,\bm{X}_{M}\}$ into CNNs-IV. The statistical average of the outputs is calculated according to the 40 independent permeability values $\{0.41,0.42,\ldots,0.80\}$. The results of the correlational analysis are shown in Fig. 15. We set the horizontal dashed line as the threshold value 0.5. Hence, each curve is divided into two parts by the horizontal dashed line. The lower part indicates that the percolation system is not penetrated; while the upper part implies that the percolation system is penetrated. The crosspoint, where the horizontal dashed line and the red curve are intersected, has a permeability value of 0.594, which is very close to the theoretical value of 0.593 that is marked by the vertical dashed line in Fig. 15. Remarkably, the critical transition point can be calculated by CNNs-IV with the preset value of 0.60 for the sampling interval of 0.01. Therefore, CNNs-IV with the preset threshold value of 0.60 is the most effective model. In further studies, the preset threshold value may need to be enhanced by smaller sampling intervals for higher precision.