Discovering Phases, Phase Transitions and Crossovers through Unsupervised Machine Learning: A critical examination

TL;DR

This paper critically examines how unsupervised machine learning, especially PCA and autoencoders, can identify phases, phase transitions, and critical points in various classical spin models, revealing both capabilities and limitations.

Contribution

It demonstrates that PCA can effectively explore phases and symmetry-breaking, interpret physical features, and locate critical points, while also highlighting PCA's limitations in capturing certain correlations.

Findings

PCA identifies phases and symmetry-breaking in classical spin models.

PCA can distinguish phase transition types and locate critical points.

Autoencoder methods can also detect phase transitions and critical points.

Abstract

We apply unsupervised machine learning techniques, mainly principal component analysis (PCA), to compare and contrast the phase behavior and phase transitions in several classical spin models - the square and triangular-lattice Ising models, the Blume-Capel model, a highly degenerate biquadratic-exchange spin-one Ising (BSI) model, and the 2D XY model, and examine critically what machine learning is teaching us. We find that quantified principal components from PCA not only allow exploration of different phases and symmetry-breaking, but can distinguish phase transition types and locate critical points. We show that the corresponding weight vectors have a clear physical interpretation, which is particularly interesting in the frustrated models such as the triangular antiferromagnet, where they can point to incipient orders. Unlike the other well-studied models, the properties of the BSI model are less well known. Using both PCA and conventional Monte Carlo analysis, we demonstrate that the BSI model shows an absence of phase transition and macroscopic ground-state degeneracy. The failure to capture the `charge' correlations (vorticity) in the BSI model (XY model) from raw spin configurations points to some of the limitations of PCA. Finally, we employ a nonlinear unsupervised machine learning procedure, the `antoencoder method', and demonstrate that it too can be trained to capture phase transitions and critical points.

Equations

$$H$$

$$\displaystyle H=-J\sum_{\langle ij\rangle}{\rm cos}\,\big{(}\,\theta_{i}-\theta_{j}\,\big{)}$$

$${\mathbf{w}}_1=\underset{\Vert \mathbf{w}\Vert =1}{\mathrm{argmax}}\left\{\sum _i{\left({\mathbf{x}}_i\cdot \mathbf{w}\right)}^2\right\}$$

$$p_{ij}={\mathbf{S}}_i\cdot {\mathbf{w}}_j$$

$$⟨|p_j|⟩=\frac{1}{n}\sum _i|p_{ij}|$$

$${\mathbf{X}}^T\mathbf{X}{\mathbf{w}}_n={\lambda }_n{\mathbf{w}}_n$$

$${\mathbf{w}}_2^{\prime }=\frac{1}{L}[\cos ({\mathbf{r}}_1{\mathbf{k}}_1),\dots ,\cos ({\mathbf{r}}_N{\mathbf{k}}_1)]$$

$$+\frac{1}{L}[\cos ({\mathbf{r}}_1{\mathbf{k}}_2),\dots ,\cos ({\mathbf{r}}_N{\mathbf{k}}_2)]$$

$$⟨E⟩$$

$$⟨C⟩$$

$$⟨F⟩$$

$$S(T_b)-S(T_a)$$

$$=\frac{E(T_b)}{T_b}-\frac{E(T_a)}{T_a}+{\int }_{T_a}^{T_b}\frac{E(T^{\prime })}{T^{\prime 2}}d{T}^{\prime }$$

$$⟨S(T)⟩$$

$$-\frac{⟨E(T_b)⟩}{T_b}+\frac{⟨E(T_a)⟩}{T_a}-{\int }_{T_a}^{T_b}\frac{⟨E(T^{\prime })⟩}{T^{\prime 2}}d{T}^{\prime }$$

$$⟨|M|⟩$$

$$⟨\chi ⟩$$

$$⟨B⟩$$

$${\oint }_C\Delta \theta dl=2\pi kk=\pm 1,\pm 2,\pm 3\dots$$

Full text

I Introduction

In the age of big data, machine learning has become an indispensable tool whose utility transcends scientific and academic boundariesWitten . From social networkingNowell ; Kyumin to object and image recognitionBishop ; Rosten , from advertising to financePang , from engineering to medicineKing , from biological physicsMcKinney ; Schafer to astrophysicsVanderPlas , wherever there is preponderance of information and real data, machine learning is helping to find and quantify patterns and even discover basic lawsCrutchfield .

Condensed Matter Physics is a relatively late arrival in this field. Although isolated applications of machine-learning have been made over many years, it is only recently that a concerted effort towards the use of these methods for addressing problems in many-body physics has started to emergeCarleo ; Broecker ; Torlai ; Carrasquilla ; Wang ; Schoenholz ; Wei ; Torlai2 ; Nieuwenburg ; Khatami . On the one hand, machine learning algorithms, such as deep learningHinton2 ; Krizhevsky , have profound connections to the foundations of statistical physicsEngel ; Mehta . Scaling and renormalizationCardy are core principles that underlie our ability to make sense of macroscopic phenomena, in particular, their simplicity and universality despite incredible microscopic complexity. It is perhaps not surprising that the way forward for machines to learn from large data sets would incorporate similar principles.

On the other hand, machine learning is also being used, mostly in conjunction with Monte CarloWolff ; Robert ; Katzgraber data, to recognize phases and phase-transitionsCarrasquilla ; Broecker ; Wang ; Nieuwenburg ; Khatami ; Wetzel , to maximize accuracy of variational Monte Carlo methods in quantum many-body systemsCarleo , to guide choices of Monte Carlo moves Walter ; Graham , to explore overcoming the famous sign-problem bottleneckBroecker , and to infer spectral functions by performing analytic continuation from imaginary-time to real-time data Arsenault . Both supervised learningCarrasquilla ; Liu_Junwei , which depends on first training the system with trial data-sets and unsupervised learningWang , which works with no prior training, have been found useful.

In this paper, we focus further on Monte Carlo simulations and examine the extent to which unsupervised machine learning can succeed in deciphering and distinguishing different physics. We introduce a simple biquadratic-exchange spin-one Ising model that has a macroscopic ground-state degeneracy which can be lifted by a small bilinear-exchange term. We contrast its behavior with more well understood and studied cases– the Ising ferromagnet on a square latticeBinder_Ising and Ising antiferromagnet on a triangular latticeWannier ; Stephenson1 ; Stephenson2 , the Blume-Capel modelBlume ; Capel ; BEG ; Jain ; Liu (which has both first and second order transitions to phases with true long range order), and the two-dimensional XY modelBarouch ; Kosterlitz , which has a low temperature power-law correlated phase.

Our primary results are: (a) Simple machine learning methods are good at recognizing order and symmetry breaking and the temperature region of sharpest change. Thus, it is possible to locate the transition temperature with reasonable accuracy. (b) Strong first-order transitions can be easily distinguished from second-order transitions by the evolution of the principal component distribution. (c) A crossover can be distinguished from a phase-transition. (d) Principal components are related to Fourier modes and order parameters and their size dependence can be used to extract critical exponents. (e) The principal components are particularly interesting in the fully frustrated triangular antiferromagnet, where they point to incipient orders upon perturbation which are not evident at first glance. (f) Finally, the analysis is sensitive to the nature of the data provided to the machine-learning algorithms. The use of bare spin-configurations versus bond, triangle or square-plaquette spin, or local vorticity configurations provide complementary inferences and are not readily discovered from the other.

In the conclusion, we will discuss why automated machine learning can be an important part of our toolbox for studying phase transitions using large scale numerical simulations.

II Models

In this paper, we investigate several classical models of phase transitions. The simplest is the two-dimensional square-lattice Ising model:

[TABLE]

where the ‘spin’ $S_{i}=\pm 1$. Positive $J$ corresponds to ferromagnetism, where the energy favors aligned spins, while for negative $J$ antiferromagnetic configurations have lower energy. On a bipartite lattice with nearest neighbor interactions, the thermodynamics of $J>0$ and $J<0$ are identical and, for the square lattice, a phase transition to a magnetically ordered state occurs as the temperature $T$ is reduced below the critical value $T_{c}/J=2/\ln(1+\sqrt{2})\approx 2.269$Onsager . We also study the antiferromagnetic triangular lattice Ising model (TLIM) as an example of a fully frustrated model with no finite temperature phase transition and macroscopic ground state degeneracy.

The Blume-Capel model (BCM) is a generalization of Eq. 1 which allows three values $S_{i}=\pm 1,0$ with associated energy,

[TABLE]

Again we choose the coupling $J=+1$ between nearest neighbor spins on a square lattice with periodic boundary conditions. Qualitatively speaking, the BCM allows ‘site vacancies’ ($S_{i}=0$) in addition to non-zero magnetic moments ($S_{i}=\pm 1$). The parameter $\Delta$ of Eq. 2 controls the vacancy density. For $\Delta\rightarrow-\infty$, $S_{i}=0$ is energetically unfavorable, and the BCM reduces to the Ising model. As $\Delta$ increases, the second order Ising-like transition from ferromagnet to paramagnet becomes first order above a tricritical point $(T/J,\Delta/J)=(0.609(4),1.965(5))$ Plascak98 and occurs with a discontinuous jump in magnetization and vacancy density. The BCM provides a description of systems ranging from metamagnets and ternary alloys, to multicomponent spins and the dynamics of rough surfaces BEG ; MukamelBlume74 ; Lajzerowicz75 ; Newman83 ; Zahraouy04 ; Brito07 . It has also been studied for nonuniform $\Delta_{i}$ Pittman08 to understand situations where the spin density varies across the lattice.

The BCM was further generalized by Blume, Emery, and Griffiths to

[TABLE]

to incorporate an additional biquadratic interaction $K\sum_{\left\langle i,j\right\rangle}S_{i}^{2}S_{j}^{2}$, and used to study 3He-4He mixtures BEG . A variant of Eq. 3, which we call the biquadratic-exchange spin-one Ising (BSI) model,

[TABLE]

will be studied here. As with the BCM and BEG model, $S_{i}$ takes values: $+1$, [math] or $-1$. Our primary interest is in the purely biquadratic limit with $J=0$. Because it has been less well studied, we discuss in this section some of the basic physics of the BSI model revealed by conventional Monte Carlo before turning to the results of machine learning approaches.

Since $K>0$, it is energetically unfavorable to have adjacent sites of the BSI model both ‘occupied’ with $S_{i}=\pm 1$. Therefore, as $T$ is lowered one expects all occupied sites to be surrounded by vacancies. As long as $J=0$, there is no preference for one $S_{i}=\pm 1$ orientation over the other. If one sublattice is occupied with $S_{i}=0$, the other sublattice is free to take any spin value. Thus the $J=0$ BSI model has macroscopic ground-state entropy. A question arises: Is there a phase transition to a checkerboard ‘charge-ordered’ phase where the $S_{i}=\pm 1$ variables preferentially occupy one sublattice? Thus, along with the triangular-antiferromagnet, this model provides machine learning with the challenge of dealing with high degeneracy and deciding between the presence and absence of a phase transition.

Simulations of this $J=0$ case Fig. 1 panels (a,b,c) can be contrasted with a second-neighbor bilinear exchange interaction $J=0.1$, where conventional Monte Carlo analysis reveals a sharp second order phase transition at temperature $T/K=0.163$. Fig. 1 panels (d,e,f) show spin configuration snapshots that go from the ordered ferromagnetic phase at $T=0.1$ to an intermediate degenerate regime at $T=0.3$, to a high temperature phase at $T=2.0$.

Finally, the XY model in two dimensions,

[TABLE]

in which sites are occupied by a continuous spin $\theta_{i}\in(-\pi,\pi]$, offers a situation in which the performance of machine learning for a Kosterlitz-Thouless (KT) Barouch ; Kosterlitz phase transition between exponentially decaying spin correlations and a line of critical points with power law decaying correlations can be evaluated when provided with various sorts of real space snapshots - the spin directions themselves or measures of local vorticity.

III Methods of Machine Learning

Broadly speaking, machine learning can be divided into ‘unsupervised’ and ‘supervised’ approaches. In the latter case, a neural network, for example, can be ‘trained’ by providing a set of patterns along with a desired ‘answer’ (eg in the case of statistical mechanics whether the pattern is ordered or not). Thus the trained (supervised) network is able to ‘generalize’ and distinguish order from disorder on patterns to which it had not been previously exposed. In the former case, patterns are given to the analysis procedure, but no ‘answers’ are provided. In this paper, we utilize two methods for unsupervised machine learning, ‘Principal Component Analysis’ and the ‘Autoencoder’ approach. These two approaches have been described in Pearson ; Wikipedia ; Jolliffe ; Bourlard ; Hinton ; Ruslan . We provide a brief review here, with further detail in the Appendix.

III.1 Principal Component Analysis

Principal Component Analysis (PCA) Pearson ; Wikipedia ; Jolliffe is a strikingly simple method. In its implementation for our statistical models, one constructs a matrix $\mathbf{S}$, each row of which is a snapshot of the instantaneous values of the degrees of freedom from a Monte Carlo simulation, e.g. a listing of spins on each lattice site. We denote by $M$ the total number of such configurations, and therefore the number of rows of $\mathbf{S}$. It is convenient to compute the mean value $m_{j}=(1/M)\sum_{i}S_{ij}$, of each column, and subtract that value from the entries in the column to obtain the ‘centered data matrix’ $\mathbf{X}$ with ‘column-wise zero empirical mean’. To study phase transitions using PCA, we choose $t$ evenly separated temperatures (or some other parameter) and for each, we generate $n$ uncorrelated spin configurations using a Monte Carlo simulation. The dimension of $\mathbf{X}$ is $M=nt$.

PCA then extracts features from the data in $\mathbf{X}$ by performing an orthogonal transformation. The result is the conversion of a set of configurations of possibly correlated variables into a set of values of linearly uncorrelated variables, named principal components, such that the greatest variance of the data comes to lie on the first principal component, the second greatest variance on the second principal component, and so onWikipedia . Each succeeding component in turn has the largest variance possible under the constraint that it is orthogonal to the preceding components and the resulting vectors are an uncorrelated orthogonal basis set.

To put this more formally, given the centered data matrix $\mathbf{X}$, the transformation is defined by iteratively constructing a set of $N$-dimensional ‘weight vectors’. The first is found by

[TABLE]

where the arguments of the maxima (abbreviated to arg max) are the points of the domain of some function at which the function values are maximized. Succeeding weight vectors are obtained by subtracting the already calculated principal components from $\mathbf{X}$ and repeating Eq. 6. The full principal components decomposition of $\mathbf{S}$ can therefore be given as $\mathbf{P}=\mathbf{S}\mathbf{W}$, where $\mathbf{W}=[{\mathbf{w}}_{1},\dots,{\mathbf{w}}_{N}]$.

The inner products of the spin configuration $\mathbf{S}_{i}$, with the weight vectors $\mathbf{w}_{j}$ are termed its components $p_{ij}$,

[TABLE]

The ‘quantified principal components’ are defined as the averages:

[TABLE]

where the average is taken over some appropriate subset of all samples. In our statistical mechanics application of PCA, it is natural to do the quantification by averaging over all $n$ spin configurations with the same temperature. The central power of PCA is in cases where the nature of the system is well described by a small number of ‘principal components’, e.g. $p_{i1}$ and $p_{i2}$ only, and their associated averages.

The principal components transformation can, equivalently, be thought of in terms of another matrix factorization method, the singular value decomposition (SVD). It can be shown that the PCA weight vectors are eigenvectors of the $N\times N$ real symmetric matrix $\mathbf{X}^{T}\mathbf{X}$, determined byJolliffe ; White :

[TABLE]

Eigenvalues returned from Eq. 9, sorted in a descending order $\lambda_{1}\geq\lambda_{2}\dots\geq\lambda_{N}\geq 0$, represent variances of the input matrix $\mathbf{X}$ along corresponding weight vectors. It is convenient to denote the normalized eigenvalues $\tilde{\lambda}_{n}=\lambda_{n}/\sum_{i=1}^{N}\lambda_{i}$ as relative variances.

III.2 The Autoencoder Approach

PCA has been shown to be effective in extracting aspects of phase transitionsWang ; Bradde . As discussed above, however, PCA is based on linear transformations of the input data. Thus, it is natural to explore other, non-linear transformation, methods.

In this spirit, we have also examined the use of an ‘autoencoder’Bourlard ; Hinton ; Ruslan , i.e. an artificial neural network used for unsupervised learning of efficient encodings by exploiting the discovery of compressed representations. The basic idea is to pass information, in our case spin configurations, through an intermediate layer whose number of ‘neurons’ is considerably smaller than the number of bits needed to encode the original, full spin configuration. The goal is to reproduce the original spin configuration as best possible despite the constraint imposed by the reduced information storage capacity of the hidden layer. By attempting this accurate reconstruction, the network may discover important structure/patterns in the input data.

A schematic autoencoder structure with two hidden neurons is shown in Fig. 2. Although the autoencoder and PCA both use the idea of dimensional reduction (in the case of the PCA by focussing only on the one or two largest (‘principal’) components of the decomposition in Eq. 9) the autoencoder approach is generally more powerful and more flexible than PCA. Rather than being restricted to linear transformations of the input data into principal components, the autoencoder can incorporate nonlinear representations. Given a set of spin configurations $\{\mathbf{S}_{i}\}$, in an autoencoder the weights in the neural network can be trained through, for example, backpropagationWitten ; Rumelhart , to return target values equal to the inputs. In this work, we choose convolutional neural networks (CNNs)Krizhevsky ; Lawrence ; Schmidhuber as the encoder and the decoder networks. As an alternative, artificial neural network based variational autoencoders can also be applied to tackle the problem, as reported in Ref. Wetzel .

IV Results and Analysis

IV.1 PCA results of the Ising model

By feeding a set of raw Ising spin configurations $\{\mathbf{S}_{i}\}$ into PCA, we generate results for the square-lattice Ising model, shown in Fig. 3. We have chosen $t=40$ temperatures in the range from $2.0$ to $2.8$ with $\Delta T=0.02$. These bracket the critical temperature $T_{c}\approx 2.269$. For each temperature, we generate $n=10,000$ uncorrelated spin configurations, so that $M=400,000$ configurations are used for calculating weight vectors and relative variances. In Fig. 3 (a), a single dominant principal component is observed; the largest $\lambda_{k}$, obtained through Eq. 9, is more than an order of magnitude greater than the next largest. As we shall further argue below, the physical implication is that there is a single dominant spin pattern for the Ising model, corresponding to alignment of all spins (i.e. ferromagnetism).

Projecting raw spin configurations into the 2D plane spanned by the first two principal components $p_{i1}$ and $p_{i2}$, we obtain Fig. 3 (b). For clarity of presentation, only $100$ scatter points for each temperature are shown. We see that above the critical temperature the points lie in a single blob at the origin, but that below $T_{c}$ a separation arises into two distinct regions.

In Fig. 3 (c), we plot the quantified first leading component $\left\langle|p_{1}|\right\rangle/L$ versus temperature, normalized to the system size $L$. This is seen to mimic the behavior associated with a measurement of the (absolute value of) magnetization from a conventional Monte Carlo analysis. Figure 3 (d) displays the quantified second leading component $\left\langle|p_{2}|\right\rangle$ versus temperature. Its development is reminiscent of that of the susceptibility $\chi$.

The weight vector corresponding to the first leading component is shown in Fig. 4 (a). Ignoring statistical/numerical variations, this weight vector is seen to be roughly constant: $\mathbf{w}_{1}=\frac{1}{L}[1,\dots,1]$. Since the first principal component $p_{i1}$ is the inner product of $\mathbf{w}_{1}$ with $\mathbf{S}_{i}$ (Eq. 6) the connection to magnetization is clear. The ability of an automatic machine learning algorithm to recognize this order parameter is not surprising as the order parameter of the model is so evident. We have checked that studying the variation of this with size at the critical temperature $T_{c}$ can be related to known critical exponents for this 2D Ising system. This is also not a surprise given that this principal component is just the order parameter of the system.

Similarly, the weight vector corresponding to the second leading component is shown in Fig. 4 (b). To understand the physical meaning of this weight vector, we compare it in Fig. 4 (c) with

[TABLE]

where $\mathbf{r}_{i}$ are the lattice sites and $\mathbf{k}_{1}=(0,\frac{2\pi}{L})$ and $\mathbf{k}_{2}=(\frac{2\pi}{L},0)$, are the two Fourier wave vectors closest to the origin $\mathbf{k}_{0}=(0,0)$. There is clear similarity between panels (b,c).

We conclude that, for the ferromagnetic Ising model, PCA is building up weight vectors corresponding to the Fourier modes of the spin configuration. In its ordered phase, the physics of the Ising model is dominated by the single point $\mathbf{k}_{0}=(0,0)$, and hence PCA reveals a single dominant eigenvalue. Subleading eigenvalues are associated with the next most ordered arrangements, i.e. with a single horizontal or vertical domain wall. When supplied with configurations of the ferromagnetic Ising model in the zero magnetization ensemble, PCA generates four large eigenvalues corresponding to horizontal or vertical domain wallsWang .

Finally, Fig. 4 (d), shows the peaks ($T^{*}$) returned from the quantified second leading component, versus the inverse of the lattice dimension ($1/L$). A linear finite size scaling fit to these peaks yields $T_{c}\sim 2.278\pm 0.015$, which agrees reasonably well with the exact result $T_{c}/J\approx 2.269$.

IV.2 PCA results of the Blume-Capel model

Following the same procedure, we feed a set of raw Blume-Capel spin configurations $\{\mathbf{S}_{i}\}$ into the PCA to generate the results shown in Figs. 5, 6. For both plots, we choose $t=40$ values of $\Delta$ and $n=10,000$ uncorrelated spin configurations for each $\Delta$. As noted earlier, the phase diagram of the BCM differs from that of the Ising model in one important respect. It possesses both first and second order transitions, separated by a tricritical point at $(T/J,\Delta/J)=(0.609(4),1.965(5))$. In Fig. 5, the PCA was supplied with a set of spin configurations in a sweep of $\Delta$ in the range $0.5<\Delta<2.5$ at fixed $T=1.0$, thus crossing the phase boundary in a second order transition. Similarly, in Fig. 6 the sweep covers the range $1.0<\Delta<3.0$ at fixed $T=0.4$, and crosses the phase boundary in a first order transition.

The panels of Fig. 5 are very similar to that of the Ising model, Fig. 3. Using the peaks returned from Fig. 5 (d) and finite size scaling, we can locate $\Delta_{c}\approx 1.70\pm 0.01$ in the thermodynamic limit. This is close to the value $\Delta_{c}\approx 1.63$ at $T=1.0$ reported in Pittman08 . Of course, to locate the critical point more accurately, one could refine the $\Delta$ range and add more samples at each temperature, much as improvements in critical values are obtained with conventional analysis.

At fixed $T=0.4$, the BCM undergoes a first order transition at $\Delta_{c}\approx 1.996$Kwak . The corresponding PCA data are shown in Fig. 6. These first order phase transition results can be easily distinguished from those at the second order transition, which were seen in Figs. 3 (b), 5 (b) to spread out uniformly from the origin. In stark contrast, the scatter points in Fig. 6 (b) spread out mostly in the second principal component, and there are few intermediate points between two phases. The latter observation is of course expected from a first order transition. We see that, PCA clearly behaves differently across first and second order transitions.

IV.3 Monte Carlo and PCA results of the BSI model

To explore thermodynamic properties of the BSI model beyond the preliminary results presented in Fig. 1 we first study its energy, specific heat and the ‘checkerboard’ structure factor, which are shown in Fig. 7 (a), (b) and (c). Their definitions are as follows:

[TABLE]

where $\sigma_{i}$ is the spin configuration from the Monte Carlo simulation, $E_{\sigma_{i}}$ is the total energy of the spin configuration $\sigma_{i}$, $N^{\prime}$ is the total number of spin configurations used for calculation and $\mathbf{Q}$ is the ordering wavevector $(\pi,\pi)$

With only the biquadratic term, the specific heat of the BSI model has a peak at $T\approx 0.8$, however, the intensity of the peak doesn’t increase with system size, indicative of only short-range order. The structure factor is shown in Fig. 7 (c). It shows a sharp increase around the same temperature as the specific heat peak. However, once again, the peak doesn’t grow with system size, indicating an absence of a phase transition for the $J=0$ BSI model. If we turn on the next nearest neighbor bilinear interaction ($J=0.1$), we can see a clear sharp peak in the specific heat at $T\approx 0.17$ which grows with system size, indicating a true phase transition. This peak is very sharp because the strength of the bilinear exchange interaction $J=0.1$ is small compared with the leading energy scale $K=1$.

For frustrated systems, an important quantity is the entropy of the system. In Fig. 7 (d), we show the entropy, calculated from:

[TABLE]

Replacing $E$ with the Monte Carlo average $\left\langle E\right\rangle$ and approximating $S(T_{b})$ with $S(+\infty)-c/T^{2}_{b}$, where $c$ is a fitting constant, we can estimate the low temperature entropy with:

[TABLE]

With only the biquadratic term, the system has a non-zero ground-state entropy of approximately $0.584$, which clearly exceeds $\ln{3}/2\approx 0.5493$ obtained from fixing one sublattice to have $S_{i}=0$ and allowing spins to take arbitrary values on the other sublattice. With a non-zero $J=0.1$, below $T_{c}$ the entropy quickly drops to zero. For the Ising antiferromagnet on a triangular lattice, we have checked that we obtain the correct ground state entropyWannier ; Stephenson1 ; Stephenson2 , with the value $0.323$, as reported in Ref. Singh, .

To locate the critical temperature accurately in the latter case, we calculate the absolute magnetization $\left\langle|M|\right\rangle$, the susceptibility $\left\langle\chi\right\rangle$ and the Binder ratioBinder $\left\langle B\right\rangle$, which are shown in Fig. 8 and defined as:

[TABLE]

where $N^{\prime}$ is the total number of spin configurations used for the calculation. Similar to $\left\langle C\right\rangle$, both $\left\langle|M|\right\rangle$ and $\left\langle\chi\right\rangle$ have small finite size effects and their transition signals agree well with the peak of $\left\langle C\right\rangle$. In Fig. 8 (c), the Binder ratio plot shows a crossing at $T_{c}=0.163$, which gives the critical point for the $J=0.1$ BSI model.

Next, we apply the PCA method to the BSI model and generate results as shown in Fig. 9. In the top row of Fig. 9, we display PCA results for the $J=0$ case. We first feed raw spin configurations $\{\mathbf{S}_{i}\}$ into PCA, however, no dominant principal components are observed (see red symbols). Projecting raw spin configurations onto the leading principal components plane, we observe a complete mixture of points of all temperatures. This is one of the failures of an automated machine learning algorithm, which we will see again in the study of the XY model. The physics of the problem relates to how the $S_{i}=\pm 1$ spin variables arrange with respect to the $S_{i}=0$ variables, not with respect to each other. But, PCA averages these $\pm 1$ variables out to zero and fails to recognize any emergent order.

However, if we feed the squared spin configurations $\{\mathbf{S}^{2}_{i}\}$ into PCA, we can get dominant principal components (blue symbols). If we further project squared spin configurations into two leading principal components, as shown in Fig. 9 (c), we are able to distinguish the low temperature phase with the high temperature phase. Actually, the first leading principal component is just the local ‘charge’ order parameter, where we take $s_{i}=+1$ and $s_{i}=-1$ as ‘charge’ one and $s_{i}=0$ as ‘charge’ zero. The charge density changes from $2/3$ at high temperatures to below $0.5$ at low temperatures. PCA analysis shows that this as a greater discriminant of different temperatures than the true emergent charge order in the checkerboard pattern. The checkerboard order forms only the second principal component. Although it has a larger relative change with temperature, in magnitude it is smaller than the extensive change in the total number of charges present.

The scatter points in Fig. 9 (c) are different from the case of first order phase transition in Fig. 6 (b) or the second order phase transition in Fig. 5 (b). The most evident difference is the absence of symmetry breaking, which in Fig. 6 (b) and Fig. 5 (b) shows up by the bifurcation of scatter points for the low temperature data at left and right ends. This is clear signature for a gradual crossover to a low temperature behavior and an absence of a phase transition.

Turning on the bilinear exchange interaction, by feeding either the raw spin configurations $\{\mathbf{S}_{i}\}$ or the squared spin configurations $\{\mathbf{S}^{2}_{i}\}$ into PCA, we observe dominant principal components, as shown in Fig. 9 (d). There is only one dominant relative variance returned from $\{\mathbf{S}^{2}_{i}\}$ but two equally weighted dominant relative variances returned from $\{\mathbf{S}_{i}\}$. In Fig. 9 (e), projecting raw spin configurations onto the plane of two leading principal components $p_{i1}$ and $p_{i2}$, we see four distinct branches extending from the center, and interestingly, all of them correspond to temperatures lower than $T_{c}$. In this model, there exist four distinct ground-states, which are shown in Fig. 10, and each distinct ground-state is associated with a branch in Fig. 9 (e). Thus we see that the number of dominant relative variances $N_{v}$ are related to the number of ground-states $N_{g}$ by $N_{g}=2^{N_{v}}$.

On the other hand, by feeding $\{\mathbf{S}^{2}_{i}\}$ into PCA and projecting data onto two leading principal components, we generate Fig. 9 (f). The original four branches become two, because replacing $s_{i}$ with $s_{i}^{2}$, will only leave us with two distinct ground-states corresponding to the charge order. In this case, it is evident that $s_{i}^{2}$ can not capture the full phase behavior of the model.

IV.4 PCA results of the TLIM

We now turn to the PCA analysis of the fully frustrated antiferromagnetic Triangular-Lattice Ising Model (TLIM). This is a particularly interesting case as there is no long-range order in the model all the way down to $T=0$. Instead, power-law spin correlations develop as temperature is lowered Stephenson1 ; Stephenson2 . We feed raw spin configurations $\{\mathbf{S}_{i}\}$ into PCA. We observe two equally weighted relative variances, as shown in Fig. 11 (a). If we project them into leading two principal components, we see a separation of high temperature scatter points, forming a central blob, from the low temperature scatter points, forming the outside periphery of the circle, as shown in Fig. 11 (b). Figures (c) and (d) show the growth of the projected principal components as the temperature is lowered. Figures (e) and (f) show the spatial patterns in the projected components.

It is interesting to observe that the two principal components correspond to patterns of ordering in which the three sub-lattices of the triangular-lattice form roughly an (m,0,-m) and an (m,-m/2,-m/2) pattern respectively. PCA analysis shows an emergent continuous XY type symmetry between these patterns. It is known moessner ; grest that these are the two incipient orders of the triangular antiferromagnet when either a weak transverse-field is added to the system or if triangular-Ising layers are stacked on top of each other with a weak interlayer coupling. These two patterns are nearly degenerate. Within the Ginzburg-Landau approach, the degeneracy is only lifted by an irrelevent 6th order anisotropy term. It is remarkable that machine learning can automatically point to such incipient order in the Monte Carlo data. Independently recognizing such incipient patterns may be one of the strengths of machine learning.

IV.5 PCA results of the XY model

We now focus on the PCA analysis of the 2D XY model by feeding various sorts of real space snapshots, the spin directions themselves and measures of local vorticity, into PCA. To express spin directions, we define an angle configuration as $[\cos(\theta_{1}),\sin(\theta_{1}),\dots,\cos(\theta_{N}),\sin(\theta_{N})]$, and a collection of angle configurations simply as $\{\cos(\theta_{i}),\sin(\theta_{i})\}$. By providing these projections of the XY spins along x axis and y axis to PCA, we observe two equally weighted relative variances in Fig. 12 (a) and concentric circles in Fig. 12 (b), which essentially recovers the rotational symmetry of XY spins in the spin space. In Fig. 12 (b), we can see that in the high temperature phase, scatter points aggregate in the center, and are separated from the low temperature scatter points, which spread out in the periphery of the circle. In figure (c) we show a plot of $\left\langle\sqrt{p^{2}_{1}+p^{2}_{2}}\right\rangle/L$, which is the order parameter for the system. Also shown in figure (c) by a dashed vertical line is the known KT transition temperature $T_{c}=0.892$Olsson . The data show an abrupt change in the temperature dependence around $T_{c}$. Differentiating power-law correlations from long-range order is possible in PCA but will require careful finite-size study of the projected principal components, which will again take us back to conventional analysis.

The KT transition is known to be driven by unbinding of vortex-antivortex pairs. Since, PCA is a linear transformation method, it can not deduce vortex-related information from the spin configurations themselves. However, we can first preprocess raw spin configurations to generate local vortex and anti-vortex configurations by the definition:

[TABLE]

where $k$ is the winding number and the integral is done for each square plaquette. Explicitly, each spin difference along a plaquette bond is rescaled into the range $(-\pi,\pi]$. The spin differences are then summed over the plaquette to give the ‘integral’ of Eq. 15. The winding number $k$ is non-zero only when the integral is a multiple of $2\pi$. Thus each square plaquette is expressed using its winding number, which is either $+1$, indicating a vortex, [math], indicating no vorticity, or $-1$, indicating an anti-vortex. However, when we feed the raw vorticity configurations into PCA, we do not observe any dominant principal components, as shown in Fig. 12 (d) red triangles. This again shows the limitation of the PCA method. If equally weighted positive and negative vorticities exist in one configuration, they will cancel out each other, and PCA fails to capture their proliferation with temperature. This is analogous to the failure of PCA in capturing the ‘charge’ density in the BSI model when feeding $\{\mathbf{S}_{i}\}$.

On the other hand, if we feed absolute vorticity configurations (defined to be +1 if there is a vortex or antivortex present and 0 otherwise) into PCA, a clear dominant principal component is observed. See blue triangles in Fig. 12 (d). The scatter points of the leading two principal components are shown in figure (e), which appears more like a crossover from previous examples than a true phase transition. We believe, seeing the vortex-antivortex unbinding transition in the 2D XY model may be difficult in PCA or any machine learning technique. The change of projected principal component with temperature is shown in Fig. 12 (f) and it shows the proliferation of vortices with increase in temperature. The dashed vertical line again shows the known transition temperature and is consistent with the temperature of onset for proliferation of vortices.

A clear message from these studies is that PCA with spatial configurations of local variables is primarily a study of relevant patterns, their symmetry, and symmetry-breaking in the model. The discrete symmetry breaking shows up in the 2 or 4-fold patterns of separation of principal component data. Continuous rotational symmetry in the spin or order parameter space is reflected in the circular pattern of scatter points.

IV.6 Autoencoder results of the Ising model

The autoencoder method and architecture are explained in Sec. IIIB and the Appendix. In Fig. 13 (a), we display reconstruction results using $200$ hidden neurons. The top row is the raw Ising configurations at several temperatures and the bottom row is the reconstructed Ising configurations. Although the reconstructed configurations lose detailed spin structures, the network succeeds to capture essential informations of the input data, such as the domain area, the domain position and also the overall magnetization.

Next, we limit the autoencoder network to only two hidden neurons and results have been shown in Fig. 13 (b). Interestingly, we obtain similar results with the PCA method, as shown in Fig. 3, and the ferromagnetic phase is reasonably distinguished from the high temperature phase.

Furthermore, we limit the autoencoder network with a single hidden neuron and we still try to reconstruct original inputs, with results revealed in Fig. 13 (c). From the plot, we can clearly observe that the transition temperature is around $2.3$, above which, the autoencoder network learns a trivial mapping function by projecting all high temperature configurations into an almost constant. However, when the temperature is below the transition temperature, the trained network is able to capture two distinct branches, which actually corresponds to two distinct ground-states of the Ising model. Thus, even with a single hidden neuron, the autoencoder network can be trained to learn key phase transition informations and help locate the critical point.

V Conclusions

In this work, we have employed unsupervised machine learning techniques, such as PCA and autoencoder, to study phase transitions in the Ising, Blume-Capel, BSI and XY models. One of our key goals was to examine critically what these methods are capturing from the Monte Carlo data and as a result where they succeed and where they fail. In the variety of models we studied here, some had a clearly discernible order parameter. In others it was a more subtle one. Some had no finite temperature phase transition just a gradual crossover from uncorrelated high temperature behavior to strongly correlated low temperature behavior, some had a power-law ordered phase, and others a long-range ordered one. Some had a first order phase transition, others a second order one. Some had continuous symmetry, others a discrete symmetry. The number of ordered phases varied between models. It is clear from our study that even a straightforward implementation of PCA or the autoencoder-network machine-learning methods can be useful tools for studying phase transitions and can distinguish between many different scenarios.

Much of our focus was on the Principal Component Analysis (PCA) method. It is clear that the PCA method, when one feeds in spin configurations from Monte Carlo data over a range of parameters, is primarily about recognizing spatial patterns of order and symmetry breaking. This is what discriminates one temperature from another. When we have a rather obvious order parameter, such as in the square-lattice Ising model, the dominant principal component must correspond to this order parameter. The subleading principal components are related to small-q behavior, and hence capture the domain-walls and low energy fluctuations in the system. Careful analysis of temperature and size dependence of the projections on to these principal components can be used to obtain the transition temperature and the critical exponents, much like in conventional analysis.

In a fully frustrated system, such as the antiferromagnetic triangular-lattice Ising model (TLIM), there is no obvious order-parameter, only subtle incipient ones. PCA automatically brings them out. This may be one of the real strengths of automated machine-learning methods. The TLIM has power-law correlations at T=0. Furthermore, it has been shownmoessner ; grest that a stack of triangular-lattice Ising models with arbitrarily weak interlayer coupling, or TLIM in a weak transverse quantum-field is susceptible to order in a 3-sublattice pattern in which the three sublattices have (m,0,-m) or (m,-m/2,-m/2) magnetization. These two patterns of order are nearly degenerate, leading to an emergent XY symmetry which is only broken by a higher order anisotropy term. PCA picks out these as the principal components and the emergent XY symmetry between them. To fully understand these one would clearly need further analysis. But, PCA could serve as a useful starting point to the analysis of Monte Carlo data.

PCA can distinguish between first and second order transitions and between crossovers and phase transitions. The most clear distinction comes in the way the projections on to the lowest few principal components spread out. The degeneracy or number of PCA components with large comparable variances and the groups in which the scatter points of data arrange themselves very clearly and nicely illustrate the number of ordered phases and the symmetries of the model.

One glaring weakness of PCA is the inability to recognize an order in $s_{i}^{2}$ when data on $s_{i}$ are fed in. We saw that problem in the BSI model, where short-range order is associated with checkerboard arrangement of charge ($s_{i}^{2}$) not the original spin ($s_{i}$) variables. We encountered the same problem in studying vortex-antivortex pairs in the 2D XY model. The PCA fails to discriminate between configurations where opposite sign charges equal in number are changing substantially with temperature. Only when the squared variables are fed into the PCA algorithms it recognizes charge order or proliferation of vortices. It suggests that it may generally be useful to apply PCA with different powers of the local spin-variable to get a more complete picture.

The 2D XY model is particularly challenging for any numerical study. It would be too much to expect an automated machine-learning algorithm to learn about unbinding of vortex-antivortex pairs from the data on raw spin configurations. In the current form, the algorithms can not even recognize the existence of vortices. Similarly, while we believe the power-law and exponentially decaying phases can be distinguished by studying the size dependence of the projected principal components, accurately locating the phase-transition between them would be difficult. Designing machine learning algorithms that can deal with such a situation will be an important achievement.

An interesting question is whether automated machine-learning approaches bear any parallels with the study of entanglement entropies in quantum systems entanglement ; Sarma . A major theme of condensed matter physics over the last decade has been the exploration of phases with unusual, or perhaps even non-existent, local patterns of order. The entanglement entropyOsborne ; Calabrese ; Cincio of quantum lattice models is touted as a quantity that, in principle, catches them all – conventional order, conformal invariance, Goldstone modes, Fermi-surfaces and topological order– although, in practice, it takes substantial effort to address each of these individually. It is clear from our study that the simple PCA, without any pre-processing of data, is more closely tied to conventional order. One can imagine defining a ‘discriminant-entropy’ in terms of the relative variances of the PCA method as a general characteristic of the complexity of the model. It remains to be seen if such a properly defined quantity can be helpful in studies of statistical models in an order-parameter independent manner. This may be an interesting avenue for future research.

Finally, we note that because they directly analyze the real space configurations of the degrees of freedom, machine learning algorithms might make useful connections to experimental situations where the order varies over the sample. For example, in studies of cold atomic gases, the confining potential leads to a variation of density and energy scales as one moves outward from the cloud center. As a consequence, metallic, magnetic, and Mott insulating phases can coexist in different regions of a single sample. By allowing the machine learning algorithm to be exposed to different spatial windows on real space data obtained from quantum microscopyHaller_Hudson ; Bakr_Gillen ; Sherson_Weitenberg , useful information about this coexistence and domain boundaries may be generated.

Acknowledgements.

The work of WH and RTS was supported by DE-NA0002908. The work of RRPS is supported by US National Science Foundation DMR-130608.

Appendix A Autoencoder Architecture

In the preceding discussion, we have demonstrated detailed results about the autoencoder method using KerasChollet . In Fig. 14, we show the detailed autoencoder architecture. In the encoder, which is surrounded by the left blue dashed rectangle, the convolutional layer 1, abbreviated as Conv1, uses $16$ learnable filters, each of which is an array of weights, of size $3\times 3$ with stride $1$ and zero padding of $1$. The filter is spatially small, but extends through the full depth of the input volume. During the forward pass, each filter can move across the width and height of the input volume to compute dot products between the filter and the input. Following the convention, stride $1$ stands for moving the filters one pixel at a time and zero padding indicates padding the input volume with zeros around the border.

Similarly, Conv2, Conv3 and Conv4 use $8$ filters of size $3\times 3$ with stride $1$ and zero padding of $1$. All pooling layers are using max pooling, which returns the maximum of the sub-region covered by the filter, with filters of size $2\times 2$ and stride $2$. In the decoder, surrounded by the right blue dashed rectangle, Conv5, Conv6 and Conv7 use $8$ filters of size $3\times 3$ with stride $1$ and zero padding of $1$. Conv8 (Conv9) uses $16$ (1) filters of size $3\times 3$ with stride $1$ and zero padding of $1$. Instead of using pooling layers, in the decoder we apply upsampling layers, which upsamples low dimensional data to high dimensional data through deconvolutional filters that can be learned just like normal convolutional kernels, with filters of size $2\times 2$. The rectified linear unit (ReLU) function is activated after every convolutional layer and the loss function used here is the binary cross entropy loss. More details about the convolutional neural network can be found in Ref. Wikipedia_CNNs, . The network is trained for $100$ epochs in each case (Fig. 13 (a), (b) or (c)) and the training data is the raw spin configurations $\{\mathbf{S}_{i}\}$ of the Ising model on a square lattice, with fixed size $N=40\times 40$.