Finite-size analysis in neural network classification of critical
  phenomena
We analyze the problem of supervised learning of ferromagnetic phase
transitions from the statistical physics perspective. We consider two systems
in two universality classes, the two-dimensional Ising model and
two-dimensional Baxter-Wu model, and perform careful finite-size analysis of
the results of the supervised learning of the phases of each model. We find
that the variance of the neural network (NN) output function (VOF) as a
function of temperature has a peak in the critical region. Qualitatively, the
VOF is related to the classification rate of the NN. We find that the width of
the VOF peak displays the finite-size scaling governed by the correlation
length exponent, $\nu$, of the universality class of the model. We check this
conclusion using several NN architectures -- a fully connected NN, a
convolutional NN and several members of the ResNet family -- and discuss the
accuracy of the extracted critical exponents $\nu$.
Vladislav Chertenkov, Evgeni Burovski, Lev Shchur
# Finite-size analysis in neural network classification of critical phenomena

###### Abstract

We analyze the problem of supervised learning of ferromagnetic phase transitions from the statistical physics perspective. We consider two systems in two universality classes, the two-dimensional Ising model and two-dimensional Baxter-Wu model, and perform careful finite-size analysis of the results of the supervised learning of the phases of each model. We find that the variance of the neural network (NN) output function (VOF) as a function of temperature has a peak in the critical region. Qualitatively, the VOF is related to the classification rate of the NN. We find that the width of the VOF peak displays the finite-size scaling governed by the correlation length exponent, \(\nu\), of the universality class of the model. We check this conclusion using several NN architectures--a fully connected NN, a convolutional NN and several members of the ResNet family--and discuss the accuracy of the extracted critical exponents \(\nu\).

_Introduction.--_ Deep learning is since recently emerging as a promising tool for studying phase transitions and critical phenomena. The pioneering observation of Ref. [1] is that training a neural network (NN) to perform a binary classification of microscopic spin states of a two-dimensional (2D) Ising model reproduces the critical temperature of the ferromagnetic phase transition, known from the exact solution [2]. Following the seminal work, a variety of approaches are being explored to test deep learning techniques in application to several models, including the Ising and \(q\)-state Potts models, percolation, XY- and clock models [3; 4; 5; 6; 7; 8; 9].

It is becoming clear that a neural network (NN) trained on an equilibrium ensemble of microscopic states can learn and predict phase transitions between macroscopic states, _in many situations_. This gives rise to a series of fundamental questions: How to interpret NN results from the physics perspective--specifically, does a NN learn the critical behavior of a universality class of a transition? What are relevant NN observables? How general is the NN approach and what are its failure modes? What limits the reliability and accuracy of these predictions? What is the role of the NN architecture?

In this Letter, we address these questions by considering two exactly solvable models in 2D, the Ising model [2] and the Baxter-Wu (BW) model [10; 11]. We train NNs to perform binary classification of microscopic spin configurations, and perform a careful finite-size scaling analysis of the classification results. We show that the second moment of the NN output displays finite-size scaling governed by the correlation length exponent, \(\nu\), of the universality class of the model. We compare predictions of several network architectures--fully connected networks (FCNN), shallow convolutional networks (CNN) and several members of the ResNet family.

We note that using the BW model turns out to be essential to be able to distinguish between the critical scaling, \(\sim 1/L^{\nu}\), from regular, analytic corrections, \(\sim 1/L\), to thermodynamic limit behavior of systems with finite linear size \(L\). While for the Ising model the correlation length exponent \(\nu=1\), the BW model belongs to the 4-state Potts universality class with \(\nu=2/3\)--thus making the critical scaling clearly distinguishable from analytic corrections. We note that the BW model, unlike other models in the same universality class, does not show any logarithmic corrections [10], which allows us to simplify the finite-size analysis.

_Models and methods.--_ We consider two classical, exactly solved models, formulated in terms of Ising spins, \(\sigma_{i}=\pm 1\) on an \(L\times L\) lattices. The Ising model [2] is defined by the Hamiltonian \(H_{\rm Is}=-J\sum_{\langle ij\rangle}\sigma_{i}\sigma_{i}\), where \(J\) is the coupling constant, and the summation runs over the pairs of nearest neighbors of the _square_ lattice with periodic boundary conditions. The BW model [10; 11] is defined on a triangular lattice, and contains three-spin interactions \(H_{\rm BW}=-J\sum_{\langle ijk\rangle}\sigma_{i}\sigma_{j}\sigma_{k}\), where the summation runs over triplets of spins which form triangular plaquettes of a triangular lattice with periodic boundary conditions.

To generate data sets for NN training and validation, we use the standard Monte Carlo (MC) simulations with Metropolis single spin flip updates [12]. We perform simulations for system sizes with \(L=48,72,96,114,\) and \(216\) for the Ising model, and \(L=48,72,96,114,\) and \(243\) for the BW model. For each system size, we perform simulations for \(N_{t}=114\) values of the temperatures between \([T_{c}-0.4;T_{c}+0.4]\), using the value of the critical temperature \(T_{c}\) from the exact solution of the corresponding model.

For each system size and for each value of the temperature, we collect \(N_{s}=1500\) "snapshots" of spin configurations generated by the MC process (here by a "snapshot" we mean a collection of \(L^{2}\) spin values, \(\pm 1\)). To make sure that snapshots are uncorrelated, we skip at least \(2\,\tau_{\rm corr}\) Monte Carlo steps between snapshots, where \(\tau_{\rm corr}\) is the integrated autocorrelation time for the magnetization [13]. For each simulation, we allow at least \(20\,\tau_{\rm corr}\) MC steps for equilibration (see Ref. [14] for a detailed discussion of our MC simulations).

_NN training.--_ We train a NN to perform binary clas
sification of snapshots for a given system size \(L\) into two classes, ferromagnetic, FM, (\(T<T_{c}\)) or paramagnetic, PM, (\(T>T_{c}\)) separately for the Ising model and the BW model.

A NN takes as input a "snapshot" of size \(L\times L\), and outputs the class scores for the FM and PM classes. We interpret the class scores as probabilities, since their sum equal unity.

We use three different network architectures: Convolutional neural network (CNN) [15], Fully-connected neural network (FCNN) [16], and Deep convolutional residual networks (ResNet) [17]. In the ResNet family we use networks with 10, 18, 34, and 50 layers. The detailed parameters of the networks and our training protocol can be found in Supplementary Materials.

_Analysis of NN outputs.--_ Once a NN is trained, we feed it with \(N\) snapshots from the testing dataset to perform the classification. In what follows we denote by \(f_{i}^{T}\) the FM class prediction for the \(i\)-th snapshot at temperature \(T\).

Averaging over the training dataset, we define the average prediction, \(F^{T}\),

\[F^{T}=\frac{1}{N}\sum_{i=1}^{N}f_{i}^{T} \tag{1}\]

and its variance, \(V^{T}\),

\[V^{T}=\frac{1}{N}\sum_{i=1}^{N}\left(f_{i}^{T}\right)^{2}-\left(\frac{1}{N} \sum_{i=1}^{N}f_{i}^{T}\right)^{2}. \tag{2}\]

Fig. 1 shows the dependence of the FM class prediction of (left image) the Ising model and (right image) the BW model with the CNN architecture. Other NN architectures give similar results. Here we only show the FM class prediction, because the PM class prediction is given by \(1-F^{T}\).

The network output, \(F^{T}\), for both models, is clearly similar to the observation of Ref. [1]: for low temperatures, \(F^{T}\approx 1\), for high temperatures, \(F^{T}\approx 0\), and the transition region clearly shrinks on increasing the system size \(L\), thus developing a step function for \(L\gg 1\). This behavior is qualitatively similar for all network architectures we considered.

According to Ref. [1]--for the Ising model, the FM prediction, \(F^{T}\), approaches the value of 0.5 for all values of the system size \(L\) at the exact value of the critical temperature, \(T_{c}=2/\ln(1+\sqrt{2})\)[2]. Since the PM prediction is simply \(1-F^{T}\), a straightforward interpretation would be that at \(T=T_{c}\), NNs are equally likely to classify a snapshot as either ferromagnetic or paramagnetic _for finite system sizes, \(L\)_.

However, our simulations of the Ising model and the BW model, Fig. 1, show that this interpretation is not entirely correct. For some lattice sizes for Ising model and for the BW model, the "equal prediction" point, \(F^{T}=1/2\) is shifted away from the value of \(T_{c}=2/\ln(1+\sqrt{2})\) known from the exact solution [10]. For FCNN architecture, the point \(F^{T}=1/2\) is shifted away to the paramagnetic phase for all lattice sizes both for the Ising and the BW models (see Fig. 1 of Supplementary Materials). Non-systematic shifts can be observed for the Ising and the BW models for different system sizes in the networks of the ResNet family. For the ResNet-50 (Fig. 5 of Supplementary Materials) for the Ising model large system sizes (96, 144, 216) are shifted to the ferromagnetic phase, while small sizes (48, 72) are shifted to the opposite side, to the paramagnetic phase. We thus conclude that \(F^{T}=1/2\) is not a reliable finite-size estimate of the critical temperature \(T_{c}\).

For the Ising model, Ref. [1], considered system sizes of up to \(L=60\) and observed that the \(F^{T}\) curves display data collapse with respect to the "scaling variable", \(tL^{1/\nu}\), where the reduced temperature \(t=(T-T_{c})/T\) is scaled by the critical exponent \(\nu\). The data collapse estimate of Ref. [1] for \(L\) up to \(L=60\), produces the values \(T_{c}=2.266\pm 0.002\) and \(\nu=1.0\pm 0.2\), consistent with the exact values of the critical temperature and the correlation length exponent for the 2D Ising universality class, \(\nu\)=1, [10]. Our numerical experiments show that data collapse is visually observed in a wide range of values of the critical exponent \(\nu\in[0.75,1.5]\), depending on the network architecture (see Figs 13-24 of Supplementary Materials for details). We stress that simply including larger system sizes does not improve correlation length exponent and critical temperature estimates due to increasing errorbars of the NN output in the critical region, cf Fig.1.

We note however, that the increase of the errorbars of \(F^{T}\), Eq. (1)--equivalently, the variance \(V^{T}\), Eq. (2)--around \(T=T_{c}\) is similar to the expected behavior of thermodynamic functions in the critical region, where second moments of observables are related to temperature derivatives of corresponding thermodynamic functions. In this spirit, we consider the second moment of the NN prediction of the FM class, Eq. (2), and hypothesize that the variance of the NN output, Eq. (2) is singular in the thermodynamic limit. This way, the observed increase of the errorbars of \(F^{T}\) around \(T\approx T_{c}\) is in fact nothing but a finite-size rounding of this divergence, governed by the correlation length exponent \(\nu\). Fig. 2 displays the temperature dependence of \(V^{T}\), which indeed shows a drastic increase around \(T=T_{c}\), and a characteristic Gaussian-like bell shape for both Ising and BW models and all network architectures. Furthermore, the widths of the bell-shaped curves decrease with increasing the system size, which is consistent with scaling behavior.

To test this hypothesis, we study the \(L\) dependence of the width of the peak of \(V^{T}\), Eq. (2). Specifically, for each value of \(L\), we fit \(V^{T}\) vs T with an unnormalized Gaussian-like Ansatz, \(V^{T}\sim\exp\left(-(T-T_{*})^{2}/2\sigma^{2}\right)\)
with \(\sigma\) and \(T_{*}\) being fit parameters, and extract the dependence of the width of \(\sigma\) on \(L\). Since there is no _a priori_ requirement that the profile is strictly Gaussian, we also perform a separate single-parameter fits the left-hand (\(T<T_{*}\)) and the right-hand (\(T>T_{*}\)) parts of the \(V^{T}\) curves. In this procedure, \(T_{*}\) is simply the location of the maximum of \(V^{T}\), and \(\sigma\) is the (only) fit parameter. For both fitting protocols, we then fit the resulting widths, \(\sigma(L)\), to a power-law Ansatz, \(\sigma(L)\sim 1/L^{1/\nu_{\sigma}}\). We perform this procedure for the Ising and the BW models and for all network architectures, and results are summarized in Tables 1 and 2.

For the Ising model, Table 1, the first observation is that the resulting values of the scaling exponent (both one-sided \(1/\nu_{\sigma}\) and two-sided \(1/\nu_{\sigma^{\pm}}\)) are consistent with the correlation length exponent for the Ising universality

Figure 1: \(F^{T}\) ferromagnetic phase predictions for the Ising model (left) and the Baxter-Wu model (right) with FCNN for various lattice sizes. The error bars correspond to the variance \(V^{T}\) of the NN prediction. The black vertical dashed line is the position of the critical temperature \(T_{c}=2/\ln(1+\sqrt{2})\), which by chance the same for both two models.

Figure 2: \(V^{T}\) variance for the Ising model (left) and the Baxter-Wu model (right) with FCNN for different lattice sizes. The black vertical dahed line is the position of the critical temperature. The solid lines are limited to the area where the Gaussian approximation was applied to extract the width \(\sigma\) for each lattice size \(L\).
class, \(\nu=1\). One notable exception is the ResNet 10- and 34-layer architectures, which shows vastly different values for exponents \(1/\nu_{\sigma^{\pm}}\) and \(1/\nu_{\sigma}\), and the resulting values are barely within the 4 standard deviations from the exact result, \(\nu=1\).

For the BW model, Table 2, the striking observation is that the scaling exponents, \(1/\nu_{\sigma}\), estimated from the width of \(V^{T}\), are consistent with the exact value of the correlation length exponent for universality class of the BW model, \(\nu=2/3\). The accuracy of fit results, Table 2, allows to conclusively distinguish this value from regular, non-singular corrections, \(\sim L^{-1}\). This is the major advantage of considering the BW model in addition to the Ising model where \(\nu=1\).

We also note that the shape of \(V^{T}\) is in fact not symmetric around the maximum--for both Ising and BW models. Allowing for different widths, \(\sigma^{+}\) and \(\sigma^{-}\) for \(T>T_{*}\) and \(T<T_{*}\), respectively, produces closer fits of \(V^{T}\). Moreover, scaling exponents, \(1/\nu_{\sigma^{+}}\) and \(1/\nu_{\sigma^{-}}\) are different--the low-temperature exponent, \(1/\nu_{\sigma^{-}}\), is consistently larger than the high-temperature exponent, \(1/\nu_{\sigma^{+}}\) -- again, for both Ising and BW models.

It is clear from Tables 1 and 2 that the values of the critical exponents, extracted from NN data are largely independent of the NN architecture, and that increasing the depth of an NN does not bring drastic improvements in exponent accuracy estimation. For networks of the ResNet family, both for the Ising model and for the BW model, some of the scaling exponents have larger errors than similar ones for simpler architectures FCNN and CNN.

We thus conclude that the width of the \(V^{T}\) peak displays finite-size scaling consistent with the universality class of a model, and that simple convolutional networks, CNN, or fully-connected, FCNN, are more appropriate for studying this class of problems, and that increasing the network depth does not automatically translate into better reliability or accuracy of the estimates--this is consistent with the conclusion of Ref. [3].

Given that the width of the \(V^{T}\) peak displays finite size scaling with the correlation length exponent, it is natural to study \(L\)-dependence of other properties of the peak: its maximum value, \(V^{T}_{\rm max}\), and the shift of the maximum from the thermodynamic limit value of \(T_{c}\). Our numerical experiments show that both maximum height and the peak shift are NN architecture dependent and do not display meaningful convergence with \(L\to\infty\).

This behavior must be contrasted with the behavior of more traditional thermodynamic observables. It is well known [18] that the position of the specific heat maximum \(T^{*}\) shifts from the critical point \(T_{C}\) with the correlation length index \(T^{*}-T_{C}\propto 1/L^{1/\nu}\), and the same behavior is found for other thermodynamic quantities due to fluctuation cutoff, when the correlation length becomes comparable with the dimensions of the system, similar to i.e. to the rounding of the magnetic susceptibility at a temperature close to the critical one [19].

We tested the deviation of the maximum VOT for both models and six networks, and the results are placed in Table 3 for the Ising model and Table 4 for the Baxter-Wu model. Note that the critical temperature values are coincidentally the same for the two models, but \(1/\nu\) is different - it is 1 for the Ising model and 1.5 for the Baxter-Wu model, and we use these values when analyzing the VOT data. A demonstration of fits can be found in Supplemental Materials. The results of the fitting are in most cases consistent within no more than five standard deviations and follow the assumption that the shift of the VOT function follows the Ferdinand-Fischer law with an exactly known exponent. The testing of the Ising model with the ResNet-50 network is the worst, and at the same time the \(T*\) estimates for the largest systems are very close to the critical temperature \(T_{C}\), as can be seen from the Fig. 23 in Supplemental materials. Surprisingly, the values of \(T*\) change more regularly with \(L^{-1/\nu}\) for the Baxter-Wu model than for the Ising model. This may be due to weaker corrections to scaling for the Baxter-Wu model (see, for discussion Ref. [20]).

_Conclusion.--_ The main result of the presented analysis is that the most reliable information on the classification of snapshots of the spin configuration of statistical mechanics systems experiencing phase transitions of the second kind is contained in the output variation (VOT) of neural networks. Namely, VOT contains information about the critical temperature and the correlation length

\begin{table}
\begin{tabular}{|l|r|r|r|} \hline \hline NN & \(1/\nu_{\sigma}\) & \(1/\nu_{\sigma^{-}}\) & \(1/\nu_{\sigma^{+}}\) \\ \hline FCNN & 1.01(1) & 1.02(13) & 0.98(4) \\ CNN & 1.06(3) & 1.11(5) & 1.07(2) \\ ResNet-10 & 1.25(3) & 1.24(7) & 1.24(3) \\ ResNet-18 & 1.17(11) & 1.41(6) & 1.08(10) \\ ResNet-34 & 1.15(16) & 1.26(7) & 1.12(24) \\ ResNet-50 & 1.20(5) & 1.21(5) & 1.31(6) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Peak widths for the Ising model. Here \(\nu_{\sigma}\) is the estimate from fitting the Gaussian profile to the \(V^{T}\). \(\nu_{\sigma}^{+}\) and \(\nu_{\sigma}^{-}\) are similar estimates where we only fit the right-hand side (resp., left-hand-side) of the \(V^{T}\) curves. See the text for discussion.
exponent. We present a VOT analysis method and extract estimates for the critical temperature and correlation length exponent of two systems in two universality classes. The results are stable when using three different architectures in the NN deep pool - CNN, FCNN and Resnet with four configurations.

We do not have theory for the network output function as the thermodynamic function in the same ensemble as the statistical mechanics model which we tested with the neural network. At the same time we found evidence that the VOT width scales with the critical length exponent \(\nu\) and demonstrated that clearly for two universality classes. This means that output function \(F(T)\) somehow connected to the fluctuation of the physical quantities of the model although the clear connection is not directly found [24].

We find no evidence that the network output function \(F(T)\) should be equal to \(1/2\) at the critical point, as stated in the pioneering work [1] -- our claim is based on careful analysis using different network architecture. Instead, we show that the variation bias of the VOT output function does not contradict the Ferdinand-Fischer picture and can be used to estimate the critical temperature. This estimate is still not under control of the desired accuracy and more work needs to be done on a sound methodology.

We would like to emphasize that the width dependence of VOT on the system size is a good candidate for extracting the exponent \(\nu\) of the critical length and gives better accuracy than the approach proposed in Ref. [1] using \(F(T)\) collapse data. We should note again that more research is needed to find a reliable way to estimate \(\nu\) from the VOT width, since not all network architectures produce \(\nu\) with the desired precision.

Research supported by the grant 22-11-00259 of the Russian Science Foundation. The simulations were done using the computational resources of HPC facilities at HSE University [21].
