Learning to find order in disorder

The Edwards–Anderson Ising spin glass model [36] is described by the Hamiltonian

$$\begin{equation}H=-\sum _{\left\langle ij\right\rangle }{J}_{ij}{s}_{i}{s}_{j}-\sum _{i}{h}_{i}{s}_{i},\end{equation} \tag{ 1 }$$

where each J_ij and h_i is a random variable drawn from a given symmetric probability distributions. We refer to each such sample as an instance of the Ising spin glass model, while each microstate of spin configurations is called a sample of the instance. In zero-field models, h_i = 0  ∀i. In this paper we primarily consider the couplings J_ij to be drawn from a Gaussian distribution of unit variance. In nonzero field models, each field h_i is drawn from a Gaussian distribution standard deviation h.

Carrasquilla and Melko [4] implemented a multi-layer perceptron (MLP)—a neural network with a single hidden layer—to distinguish the ferromagnetic and paramagnetic states on the two-dimensional Ising model and detecting its phase transition regime by finding the point where the classification probabilities cross. Using arrays of spin configurations sampled from Monte Carlo simulations as inputs, their neural network was trained to classify single configuration samples as being in a ferromagnetic or paramagnetic state. With the resulting classification probabilities, one can identify a transition temperature T_c from where the probabilities cross. The neural network successfully discriminates between the ferromagnetic (FM) and paramagnetic (PM) phases for linear system sizes between L = 10 and L = 60, with a natural classification error occurring in the vicinity of the phase transition temperature T_c ≈ 2.269 [37]. They find that the representation of the ferromagnetic phase in the MLP is tied to learning the average magnetization m by comparing it with a toy-model MLP in which the hidden layer is simplified to directly measure m, while leaving a single free parameter that is learned during training.

Following the initial results of reference [4], Ch'ng et al [5] implemented a three-dimensional convolutional neural network (CNN) to distinguish the Néel transition in the three-dimensional Fermi–Hubbard model at half-filling, and extrapolated its magnetic phase diagram as the chemical potential is varied, which is problematic to study as the sign problem of quantum Monte Carlo becomes significant. This paper has the analogous objective of identifying the spin-glass state in the Edwards–Anderson model and establishing its phase boundary, with the benefit that past Monte Carlo studies are available to compare with.

An outline of the design framework for a phase classifying neural network on spin glasses is as follows:

• (a)  
  The spin configuration samples of the spin-glass model at low temperatures are visually indistinguishable from paramagnetic ones. To overcome this problem, as inputs for learning we use replica overlap samples

  $$\begin{equation}{q}_{i}={s}_{i}^{\left(1\right)}{s}_{i}^{\left(2\right)}\end{equation} \tag{ 2 }$$

  taken from two independent spin samples (1) and (2). The usual overlap parameter q [11] is the mean of the q_i with i ∈ {1, ..., N} and N = L³ the number of spins.
• (b)  
  The networks in reference [4] infer a phase using a single Monte Carlo configuration sample at a time as the input to the neural network. To ensure reliable phase classification, we implement a CNN that averages over multiple configuration samples from a Monte Carlo simulation.
• (c)  
  To learn a reasonable representation for the spin-glass state, a neural network should consider configurations across multiple instances simulated at the same temperature T and with the same field strength h, as the critical properties for a spin-glass transition are considered from a quenched average. Thus, we propose that an averaging step over a sample of instances (instance averaging) is necessary for a NN to classify the spin-glass state faithfully.
• (d)  
  It is possible that a supervised binary classifier for the spin glass (SG) and PM phases may learn to represent the SG state in a way that is only sensitive to the magnitude of the overlap parameter q. While such a classifier would probably be successful in a zero-field model, one could question whether such a classifier has truly "learned" the spin-glass state. However, teaching a classifier to also distinguish between the SG and FM phases helps create a better knowledge representation of the SG state that is not too directly tied to |q|. For both a typical spin-glass replica and for the Ising ferromagnet, there is a large probability of q ≈ ±1 below the critical point. However, intermediate values of q between −1 and 1 are likely only in a spin-glass model due to nontrivial ergodicity breaking. Thus, this motivates a design for a phase classifier that is taught three categories, i.e., a ternary classifier, of the FM, PM, and SG states.
• (e)  
  When generalizing to the nonzero field case, the mean of each q_i is biased because the spins become more likely to point in the direction of their local field. Thus, each q_i should be centered around its thermal mean $\left\langle {q}_{i}\right\rangle$.

While our framework sets the conditions to make the machine learning of the spin-glass state reachable, it is not without likely limitations in applicability and explanatory power. It establishes the machine learning task as that of one where a neural network learns the typical pattern of spatial correlations in thermally sampled spin overlap configurations from typical instances of spin glasses. In the non-zero field case of spin glasses, atypical values of q during Monte Carlo simulations can cause strong fluctuations that may lead to an unsuccessful finite size scaling analysis [38, 39]. Likewise, one could expect that the dominance of rare events in q may reduce the capacity for a neural network to learn and identify a phase in the thermodynamic limit.

Suppressing rare-event fluctuations and constructing a processed data set that is more representative of the thermodynamic limit would likely improve the performance of machine learning methods by reducing the amount of noise in the data. While we do not apply additional processing on our data set here, such a method could be based on simply constraining the sampled value of q to a fixed value and training on the filtered data set [40].

While we have proposed the generalization task of classifying a phase at non-zero field based on training at zero-field, there is substantial theoretical and numerical evidence that the spin-glass state is inherently different at non-zero field [13, 29]. As such, simply training on FM, PM, and SG states at zero field may result in a learned representation that is not accurate in the thermodynamic limit at non-zero field.

In spite of theoretical concerns for three dimensional systems, our proposed framework is not necessarily restricted to 3D spin glasses. One could in principle implement analogous higher-dimensional convolutions to probe the physics of higher-dimensional spin glasses closer to the critical dimension. The four-dimensional case would be a natural modification, given its more numerically conclusive AT line [41]. However, this is not without some additional software design effort, as high N-d convolutions are neither ubiquitous in the machine learning field nor natively supported in the Tensorflow API (see https://github.com/tensorflow/tensorflow/issues/9513).

The ternary state classifying neural network (hereafter referred to as the classifier) is a three-dimensional convolutional neural network with two convolutional layers that extract spatial information from the overlap samples from each instance, followed by two fully-connected layers with intermediate averaging steps. The classifier is illustrated in figure 1 and detailed in algorithm 1. The predictions of the CNN for a single example manifest in line 10 as a vector y_i of three real numbers, which are not necessarily bounded. The last step applies the 'softmax' function to create probabilities

$$\begin{equation}{p}_{i}=\frac{{\mathrm{e}}^{{y}_{i}}}{{\sum }_{j}{\mathrm{e}}^{{y}_{j}}}\end{equation} \tag{ 3 }$$

that can be interpreted as the degree of belief that a given example belongs to in each class (state) i. We refer to these probabilities as softmax probabilities. However, the results presented here utilize a classification probability of a category i, which is the rate that i is the most likely class (largest p_i) for the instances at a given temperature. This is the more practical measure for testing a NN.

Zoom In Zoom Out Reset image size

During training we wish to align the prediction of the CNN as close as possible to the actual label of the example by adjusting the parameters of each layer of the CNN. In other words, we would like to match the softmax probabilities p_i as close as possible with the indicator probabilities d_i, where d_i = 1 if i is the index of the correct state, and 0 otherwise. This can be cast as a numerical optimization problem of a cost function C(p_i, d_i). Here, the cost function is the negative cross entropy of the probabilities,

$$\begin{equation}C\left({p}_{i},{d}_{i}\right)=-\sum _{i}\left[{d}_{i} \mathrm{ln}{p}_{i}+\left(1-{d}_{i}\right) \mathrm{ln}\left(1-{p}_{i}\right)\right].\end{equation} \tag{ 4 }$$

The general strategy in machine learning is to calculate the numerical gradient of the cost function for a batch of examples and adjust the NN parameters in the negative gradient direction until a minimum for C is reached. The simplest gradient strategy is stochastic gradient descent (SGD), where each iteration of optimization uses a random batch of examples. However, in this paper we use the ADAM [42] method, an adaptive variant of SGD, for training the parameters of the neural network along with an ℓ² penalty on the cost function. Further details on NN training can be found, for example, in references [1, 2].

In algorithm 1, the hyperparameters B, Q, and R specify the grouping sizes of the training examples for every training step. B is the number of examples in each step (the batch size), Q is the number of spin-glass instances randomly selected in each example, and R is number of overlap configuration samples for each instance. Per the second and third point of the design framework in section 2.2, a single example is a 2-tensor of three-dimensional configurations, with dimensions Q × R × L³, which is progressively reduced through the layers of the CNN into the vector p_i of softmax probabilities. In every step of training, a batch of these examples of size B is evaluated into softmax probabilities, and the cost function is averaged across the batch to evaluate the gradient of the cost function.

The convolutional layers are 'pool-less', i.e., the output tensor sizes are instead controlled through the use of kernel strides and periodic padding of the overlap samples. The primary activation function used is the rectified linear unit (ReLU) function, which simply applies max(x, 0) on each element x of the outgoing tensor. The first convolutional layer uses the concatenated ReLU (CReLU) variant activation function, which uses twice the number of output channels as kernels to additionally utilize the negatives of learned features. This avoids the neural network having to duplicate learning effort on closely related kernels [43].

The first layer is capable of capturing correlations within a three spin radius. As the kernels have a stride of two, information redundancy is minimized while not completely decimating the spins into blocks of side length three. In the second layer, dilation [44] introduces the capacity for the CNN to represent longer-range correlations. In our model, each individual unit in the second convolutional layer is capable of capturing correlations within a region with a side length of 7 spins. Thus, if a large number of units in a particular channel capture the correlations they are trained to find, the spatial average is capable of expressing any long-range correlations that may be present. While deep neural networks excel at integrating such correlations for inference task, it should be advised for our particular task that they may not necessarily make use of features that are directly physically relevant after training. However, in order to successfully distinguish the SG from the FM and PM states, the neural network should at minimum be able to recognize nontrivial wave-vectors of the spin overlap, hence our inclusion of all three states in our training procedure, and our choice of basing the neural network around two convolutional layers, each with a large receptive field.

The systematic testing of models with additional layers and parameters would substantially increase the computational time and complexity of our intended objectives. Thus, the design we present here is only one possible example with the desired learning capacities, as discussed above. Adding layers and fine-tuning the model would likely increase the test accuracy, but at diminishing returns for the mere task of classifying between three possibilities.

We use parallel tempering (PT) Monte Carlo [45] to generate spin configurations from spin-glass instances, as well as a few ferromagnet runs. Sample decorrelation is enhanced by two concurrent, independent runs of PT for each instance. Four concurrent runs are used in the case of nonzero field spin glasses. Our training procedure involves the shuffling and sampling of spin states obtained during Monte Carlo sampling. In combination with PT and training in batches of multiple instances, the trained neural network is expectedly robust against time correlations.

We generate instances for the training data as outlined in table 1. The system size of L = 8 is small enough to thermalize easily and thus to simulate a large number of instances, yet large enough to fit with larger systems in a finite-size scaling analysis, e.g., as done in in reference [46]. A large amount of L = 8 samples can thus provide the classifier a 'big picture' view of the difference between the spin-glass state and the paramagnetic state with much more uncorrelated data. A larger system size (L = 12) is then useful for improving the classifier's characterization near the critical temperature. For our implementation, the size parameters of algorithm 1, the learn rate parameter η, and the ℓ² penalty weight parameter λ, are specified in table 2. Error bars in classification probabilities and performance are derived from independently training the neural network 16 times on the training instances.

Because the classifier relies on supervised training, a critical temperature needs to be presumed before labeling the data. While this could be estimated with Monte Carlo statistics or quoted from previous results, we instead utilize the confusion method introduced in Ref. [47] to find the critical temperature. We repeat the entire training procedure for various choices of T'_c within the critical region, and select the choice that maximizes the accuracy of the classifier. We choose as our accuracy measure:

$$\begin{align}\hfill \mathcal{P}& ={\int }_{{T}_{\mathrm{min}}}^{{T}_{\mathrm{max}}}\mathrm{d}T \left[{p}_{\text{SG}}\left(T\right)\theta \left({T}_{\mathrm{c}}^{\prime }-T\right)\right.\hfill \\ \hfill & \quad \left.+ {p}_{\text{PM}}\left(T\right)\theta \left(T-{T}_{\mathrm{c}}^{\prime }\right)\right],\hfill \end{align} \tag{ 5 }$$

where p_SG and p_PM are the classification probability densities of the SG and PM states between T_min and T_max, respectively. Naturally, because data are only collected for a discrete set of temperatures, $\mathcal{P}$ should be calculated by a trapezoid-rule summation of the classification probabilities divided by the temperature interval measure T_min − T_max. It is important to weigh the classification probabilities as a numerical integral–and not just sum them directly–to take the uneven spacing of the simulation temperatures into account.

To simplify training, the critical temperature for classifying the FM state is simply quoted as T_c = 4.51 [48]. We do not attempt to optimize this transition temperature with the confusion method.

The confusion method is carried out with a validation set of 100 spin-glass instances with L = 12, which were not used for training the classifier. The average accuracy as a function of T'_c shown in figure 2. The accuracy peak suggests a transition temperature of approximately T_c = 0.958(8). This is consistent with the critical temperature of the Edwards–Anderson model with Gaussian disorder as found, for example, in reference [46]. The precise temperatures used as the boundary between the PM and SG states that are within the range of both T_c values are 0.953 and 0.958. In the figures that follow, we use the classifier model trained at T'_c = 0.953, however, the results from either model are practically identical.

Zoom In Zoom Out Reset image size

The configuration samples for testing the classifier are obtained with parallel tempering Monte Carlo simulations, with parameters listed in table 3. For spin glasses in a field, we use the same temperature range and at least 10 times the number of sweeps as reference [22] used up to L = 8 to ensure termalization. The samples are collected after a thermalization period of 75% of the sweeps, for a time between samples of 2¹⁶ sweeps.

Figure 3 shows the classification probabilities of the classifier evaluated on the sample of spin-glass test instances. In all figures, vertical dashed lines represent transition temperatures and horizontal dashed lines the 50% probability line. The classification probabilities follow a smooth transition near the critical temperature. The classifier performs well, except for an anomalous jump in FM classification at the lowest temperature. This is likely due to the possibility of sampling instances with simple energy landscapes when dealing with finite-size spin glasses. This is supported by observing the finite size behavior as L increases to 12, for which the FM anomaly vanishes.

Zoom In Zoom Out Reset image size

Figure 4 shows the classification probabilities of the classifier on a test set of three-dimensional Ising ferromagnet samples, note that in the CNN algorithm, the instance averaging step is kept for code consistency. An 'instance' of the Ising ferromagnet is simply another set of uncorrelated configuration samples. The distinction that the classifier learns is affected by finite-size effects, i.e., the labeling that the classifier learns for L = 8 (T'_c = 4.51) shifts for L = 16 (T'_c = 4.46) to a lower value. Observables that would diverge at the phase transition, such as heat capacity, in the thermodynamic limit are instead critical at higher temperatures for small system sizes. The shift away from T_c at system sizes larger the training set is expected since the larger system sizes become critical at slight lower temperature than the smallest system size. This causes the predicted transition of a larger system to seem lower than the learned representation.

Zoom In Zoom Out Reset image size

The classifier's sharp transition region is likely a consequence of the number of available training parameters in the model, compared to the complexity of the learning task of distinguishing between the PM and FM states on a given system size, as well as the enhancement furnished by averaging over multiple configuration samples.

In principle, the characteristics of these classification probabilities may be useful as a NN-based finite-size scaling technique. However, neither learning the three-dimensional Ising ferromagnetic transition temperature to high precision nor inferring its critical properties from the NN were initial objectives in this work.

Zoom In Zoom Out Reset image size

As a first generalization test, we generate a test set of spin-glass instances with bimodal instead of Gaussian disorder. In this case, the critical temperature shifts to T_c ≈ 1.12 [46]. As seen in figure 5, the crossover point of the bimodal test set easily shifts to this new temperature. We take this as very good evidence that the classifier learns a good representation of the spin-glass state, and can make accurate predictions for different spin-glass models with different disorder distributions.

Now that we have demonstrated that our NN can detect a spin-glass transition at zero field and even detect correct transition temperatures when the disorder is changed, we evaluate its prediction for a spin-glass state in a field. The classification probabilities for nonzero fields are shown in figure 6. It can be observed that at h = 0.05, the spin-glass signal for the classifier becomes weak. At a stronger field, h = 0.10, the PM phase becomes dominant with larger linear system sizes L. Figure 7 illustrates a break down of any crossing temperature beyond system sizes of L = 6.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We emphasize again that the model is trained under the limitations explained for the framework, and the system sizes examined here may not be conclusive in the thermodynamic limit. However, at least for the system sizes examined here, we can observe that the transition from the spin-glass state to the PM state with increasing h becomes sharper as the system size is increased. From this trend, one would conclude from the results of the classifier that if there is a spin-glass state in a field, it occurs for fields h ≲ 0.05. This is consistent with the correlation length-based analysis of reference [22] for similar system sizes.