Unsupervised learning of phase transitions: from principal component analysis to variational autoencoders

TL;DR

This paper demonstrates that unsupervised machine learning methods, including PCA and variational autoencoders, can identify phase transitions and order parameters in statistical physics models without prior knowledge.

Contribution

It introduces the use of neural network-based autoencoders to detect phase transitions and order parameters in Ising and XY models, extending previous PCA-based approaches.

Findings

Latent parameters match known order parameters.

States form clusters indicating phases.

Reconstruction loss signals phase transitions.

Abstract

We employ unsupervised machine learning techniques to learn latent parameters which best describe states of the two-dimensional Ising model and the three-dimensional XY model. These methods range from principal component analysis to artificial neural network based variational autoencoders. The states are sampled using a Monte-Carlo simulation above and below the critical temperature. We find that the predicted latent parameters correspond to the known order parameters. The latent representation of the states of the models in question are clustered, which makes it possible to identify phases without prior knowledge of their existence or the underlying Hamiltonian. Furthermore, we find that the reconstruction loss function can be used as a universal identifier for phase transitions.

Equations

$$H=H(\mathbf{S})=-J\sum _{<i,j{>}_{nn}}{s}_i{s}_j,$$

$$\underline{S}\widehat{=}{\left(\begin{array}{ccccc}\uparrow & \downarrow & \uparrow & \dots & \uparrow \\ \downarrow & \uparrow & \uparrow & \dots & \uparrow \\ ⋮& ⋮& ⋮& & ⋮\\ \downarrow & \downarrow & \uparrow & \dots & \downarrow \end{array}\right)}_{L\times L}.$$

$$M(\mathbf{S})=\frac{1}{N}\sum _i{s}_i.$$

$$Z=\sum _{\mathbf{S}\in \Lambda }\exp (-H(\mathbf{S})/T)$$

$$⟨\Vert M(T)\Vert ⟩=\frac{1}{Z}\sum _{\mathbf{S}\in \Lambda }\Vert M(\mathbf{S})\Vert \exp (-H(\mathbf{S})/T).$$

$$M_{st}=M(\underline{S}\odot A).$$

$$H(\mathbf{S})=-J\sum _{<i,j{>}_{nn}}{\text{s}}_i\cdot {\text{s}}_j,$$

$$\underset{\Vert \mathbf{w}\Vert =1}{\mathrm{argmax}}\left[\sum _j{\left(({\mathbf{S}}_{(j)}-\mu )\cdot \mathbf{w}\right)}^2\right],$$

Full text

I Introduction

Inferring macroscopic properties of physical systems from their microscopic description is an ongoing work in many disciplines of physics, like condensed matter, ultra cold atoms or quantum chromo dynamics. The most drastic changes in the macroscopic properties of a physical system occur at phase transitions, which often involve a symmetry breaking process. The theory of such phase transitions was formulated by Landau as a phenomenological model Landau (1937) and later devised from microscopic principles using the renormalization group Kadanoff (1966); Wilson (1975). One can identify phases by knowledge of an order parameter which is zero in the disordered phase and nonzero in the ordered phase.

Whereas in many known models the order parameter can be determined by symmetry considerations of the underlying Hamiltonian, there are states of matter where such a parameter can only be defined in a complicated non-local way gang Wen (2004). These systems include topological states like topological insulators, quantum spin hall states Kane and Mele (2005) or quantum spin liquids Anderson (1973). Therefore, we need to develop new methods to identify parameters capable of describing phase transitions.

Such methods might be borrowed from machine learning. Since the 1990s this field has undergone major changes with the development of more powerful computers and artificial neural networks. It has been shown that such neural networks can approximate every function under mild assumptions Cybenko (1989); Hornik (1991). They quickly found applications in image classification, speech recognition, natural language understanding and predicting from high-dimensional data. Furthermore, they began to outperform other algorithms on these tasks Krizhevsky et al. (2012).

In the last years physicists started to employ machine learning techniques. Most of the tasks were tackled by supervised learning algorithms or with the help of reinforcement learning Curtarolo et al. (2003); Rupp et al. (2012); Li et al. (2015); LeDell et al. (2012); Pilania et al. (2015); Saad et al. (2012); Ovchinnikov et al. (2009); Arsenault et al. (2014); Snyder et al. (2012); Hautier et al. (2010); Carrasquilla and Melko (2016); Kasieczka et al. (2017); Carleo and Troyer (2017). Supervised learning means one is given labeled training data from which the algorithm learns to assign labels to data points. After successful training it can then predict the labels of previously unseen data with high accuracy.

In addition to supervised learning, there are unsupervised learning algorithms which can find structure in unlabeled data. They can also classify data into clusters, which are however unlabelled. It is already possible to employ unsupervised learning techniques to reproduce Monte-Carlo-sampled states of the Ising model Torlai and Melko (2016). Phase transitions were found in an unsupervised manner using principal component analysis Wang (2016); van Nieuwenburg et al. (2017). We employ more powerful machine learning algorithms and transition to methods that can handle nonlinear data. A first nonlinear extension is kernel principal component analysis Scholkopf et al. (1999).

The first versions of autoencoders have been around for decades Bourlard and Kamp (1988); Hinton and Zemel (1993) and were primarily used for dimensional reduction of data before feeding it to a machine learning algorithm. They are created from an encoding artificial neural network, which outputs a latent representation of the input data, and a decoding neural network that tries to accurately reconstruct the input data from its latent representation. Very shallow versions of autoencoders can reproduce the results of principal component analysis Baldi and Hornik (1989).

In 2013, variational autoencoders have been developed as one of the most successful unsupervised learning algorithms Kingma and Welling (2013). In contrast to traditional autoencoders, variational autoencoders impose restrictions on the distribution of latent variables. They have shown promising results in encoding and reconstructing data in the field of computer vision.

In this work we use unsupervised learning to determine phase transitions without any information about the microscopic theory or the order parameter. We transition from principal component analysis to variational autoencoders, and finally test how the latter handles different physical models. Our algorithms are able to find a low dimensional latent representation of the physical system which coincides with the correct order parameter. The decoder network reconstructs the encoded configuration from its latent representation. We find that the reconstruction is more accurate in the ordered phase, which suggests the use of the reconstruction error as a universal identifier for phase transitions.

Whereas for physicists this work is a promising way to find order parameters of systems where they are hard to identify, computer scientists and machine learning researchers might find an interpretation of the latent parameters.

II Models

II.1 Ising Model in 2d

The Ising model is one of the most-studied and well-understood models in physics. Whereas the one-dimensional Ising model does not possess a phase transition, the two-dimensional model does. The Hamiltonian of the Ising model on the square lattice with vanishing external magnetic $h$ field reads

[TABLE]

with uniform interaction strength $J$ and discrete spins $s_{i}\in\{+1=\ \uparrow,\ -1=\ \downarrow\}$ on each site $i=1,...,N$. The notation $<i,j>_{nn}$ indicates a summation over nearest neighbors. A spin configuration $\mathbf{S}=(s_{1},...,s_{N})$ is a fixed assignment of a spin to each lattice site, $\Lambda$ denotes the set of all possible configurations $\mathbf{S}$. We set the Boltzmann constant $k_{B}=1$ and the interaction strength $J=1$ for the ferromagnetic case and $J=-1$ for the antiferromagnetic case. A spin configuration $\mathbf{S}$ can be expressed in matrix form as

[TABLE]

Lars Onsager solved the two dimensional Ising model in 1944 Onsager (1944). He showed that the critical temperature is $T_{c}=2/\ln(1+\sqrt{2})=2.269$.

For the purpose of this work, we assume a square lattice with length $L=28$ such that $L\times L=N=784$, and periodic boundary conditions. We sample the Ising model using a Monte-Carlo algorithm Metropolis and Ulam (1949) at temperatures $T\in[0,5]$ to generate $50\,000$ samples in the ferromagnetic case and $10\,000$ samples in the antiferromagnetic case. The Ising model obeys a discrete $\mathbb{Z}_{2}$-symmetry, which is spontaneously broken below $T_{c}$. The magnetization of a spin sample is defined as

[TABLE]

The partition function

[TABLE]

allows us to define the corresponding order parameter. It is the expectation value of the absolute value of the magnetization at fixed temperature

[TABLE]

Similarly, with the help of the matrix $A_{ij}=(-1)^{i+j}$, we define the order parameter, as the expectation value of the staggered magnetization. The latter is calculated from an element-wise product with a matrix form of the spin configurations

[TABLE]

II.2 XY Model in 3d

The Mermin-Wagner-Hohenberg theorem Mermin and Wagner (1966); Hohenberg (1967) prohibits continuous phase transitions in $d\leq 2$ dimensions at finite temperature when all interactions are sufficiently short-ranged. Hence, we choose the XY model in three dimensions as a model to probe the ability of a variational autoencoder to classify phases of models with continuous symmetries. The Hamiltonian of the XY model reads

[TABLE]

with spins on the one-sphere $\textbf{s}_{i}\in\mathbb{R}^{2},\left\lVert\textbf{s}_{i}\right\rVert=1$. Employing $J=1$, the transition temperature of this model is $T_{c}=2.2017$ Gottlob and Hasenbusch (1993) Using a cubic lattice with $L=14$, such that $N=L^{3}=2744$, we perform Monte-Carlo simulations to create 10 000 independent sample spin configurations in the temperature range of $T\in[0,5]$. The order parameter is defined analogously to the Ising model magnetization (5), but with the $L^{2}$-norm of a magnetization consisting of two components.

III Methods

Principal component analysis Pearson (1901) is an orthogonal linear transformation of the data to an ordered set of variables, sorted by their variance. The first variable, which has the largest variance, is called the first principal component, and so on. The linear function $\langle\cdot,\mathbf{w}\rangle$, which maps a collection of spin samples $(\mathbf{S}_{(1)},...,\mathbf{S}_{(n)})$ to its first principal component, is defined as

[TABLE]

where $\mathbf{\mu}$ is the vector of mean values of each spin averaged over the whole dataset. Further principal components are obtained by subtracting the already calculated principal components and repeating (8).

Kernel principal component analysis Scholkopf et al. (1999) projects the data into a kernel space in which the principal component analysis is then performed. In this work the nonlinearity is induced by a radial basis functions kernel.

Traditional neural network-based autoencoders consist of two artificial neural networks stacked on top of each other. The encoder network is responsible for encoding the input data into some latent variables. The decoder network is used to decode these parameters in order to return an accurate recreation of the input data, shown in Fig. 1. The parameters of this algorithm are trained by performing gradient descent updates in order to minimize the reconstruction loss (reconstruction error) between input data and output data.

Variational autoencoders are a modern version of autoencoders which impose additional constraints on the encoded representations, see latent variables in Fig. 1. These constraints transform the autoencoder to an algorithm that learns a latent variable model for its input data. Whereas the neural networks of traditional autoencoders learn an arbitrary function to encode and decode the input data, variational autoencoders learn the parameters of a probability distribution modeling the data. After learning the probability distribution, one can sample parameters from it and then let the encoder network generate samples closely resembling the training data.

To achieve this, variational autoencoders employ the assumption that one can sample the input data from a unit Gaussian distribution of latent parameters. The weights of the model are trained by simultaneously optimizing two loss functions, a reconstruction loss and the Kullback-Leibler divergence between the learned latent distribution and a prior unit Gaussian.

In this work we use autoencoders and variational autoencoders Chollet (2014) with one fully connected hidden layer in the encoder as well as one fully connected hidden layer in the decoder, each consisting of 256 neurons. The number of latent variables is chosen to match the model from which we sample the input data. The activation functions of the intermediate layers are rectified linear units. The activation function of the final layer is a sigmoid in order to predict probabilities of spin $\uparrow$ or $\downarrow$ in the Ising model, or tanh for predicting continuous values of spin components in the XY model. We do not employ any $L^{1},L^{2}$ or Dropout regularization. However, we tune the relative weight of the two loss functions of the variational autoencoder to fit the problem at hand. The Kullback-Leibler divergence of the variational autoencoder can be regarded as reguarization of the traditional autoencoder. In our autoencoder the reconstruction loss is the cross-entropy loss between the input and output probability of discrete spins, as in the Ising model. The reconstruction loss is the mean-squared-error between the input and the output data of continuous spin variables in the XY model.

To understand why a variational autoencoder can be a suitable choice for the task of classifying phases, we recall what happens during training. The weights of the autoencoder learn two things: on the one hand, they learn to encode the similarities of all samples to allow for an efficient reconstruction. On the other hand, they learn a latent distribution of the parameters which encode the most information possible to distinguish between different input samples.

Let us translate these considerations to the physics of phase transitions. If all the training samples are in the unordered phase, the autoencoder learns the common structure of all samples. The autoencoder fails to learn any random entropy fluctuations, which are averaged out over all data points. However, in the ordered phase there exists a common order in samples belonging into the same phase. This common order translates to a nonzero latent parameter, which encodes correlations on each input sample. It turns out that in our cases this parameter is the order parameter corresponding to the broken symmetry. It is not necessary to find a perfect linear transformation between the order parameter and the latent parameter as it is the case in Fig. 3. A one-to-one correspondence is sufficient, such that one is able to define a function that maps these parameters onto each other and captures all discontinuities of the derivatives of the order parameter.

We point out similarities between principal component analysis and autoencoders. Although both methods seem very different, they both share common characteristics. Principal component analysis is a dimensionality reduction method which finds the linear projections of the data that maximizes the variance. Reconstructing the input data from its principal components minimizes the mean squared reconstruction error. Although the principal components do not need to follow a Gaussian distribution, principal components have the highest mutual agreement with the data if it emerges from a Gaussian prior. Moreover, a single layer autoencoder with linear activation functions closely resembles principal component analysis Baldi and Hornik (1989). principal component analysis is much easier to apply and in general uses less parameters than autoencoders. However, it scales very badly to a large dataset. Autoencoders based on convolutional layers can reduce the number of parameters. In extreme cases this number can be even less than the parameters of principal component analysis. Furthermore, such autoencoders can promote locality of features in the data.

IV Results

IV.1 Ising Model

The four different algorithms can be applied to the Ising model to determine the role of the first principal components or the latent parameters. Fig. 3 shows a clear correlation between these parameters and the magnetization for all four methods. However, the traditional autoencoder is inaccurate; this fact leads us to enhancing traditional autoencoders to variational autoencoders. The principal component methods show the most accurate results, slightly better than the variational autoencoder. This is to be expected, since the former are modeled by fewer parameters.

In the following results section, we concentrate on the variational autoencoder as the most advanced algorithm for unsupervised learning.

To begin with, we choose the number of latent parameters in the variational autoencoder to be one. After training for 50 epochs and a saturation of the training loss, we visualize the results in Fig. 2. On the left, we see a close linear correlation between the latent parameter and the magnetization. In the middle we see a histogram of encoded spin configurations into their latent parameter. The model learned to classify the configurations into three clusters. Having identified the latent parameter to be a close approximation to the magnetization $M(S)$ allows us to interpret the properties of the clusters. The right and left clusters in the middle image correspond to an average magnetization of $M(S)\approx\pm 1$, while the middle cluster corresponds to the magnetization $M(S)\approx 0$. Employing a different viewpoint, from Fig. 2 we conclude that the parameter which holds the most information on how to distinguish Ising spin samples is the order parameter. On the right panel, the average of the magnetization, the latent parameter and the reconstruction loss are shown as a function of the temperature. A sudden change in the magnetization at $T_{c}\approx 2.269$ defines the phase transition between paramegnetism and ferromagnetism. Even without knowing this order parameter, we can now use the results of the autoencoder to infer the position of the phase transition. As an approximate order parameter, the average absolute value of latent parameter also shows a steep change at $T_{c}$. The averaged reconstruction loss also changes drastically at $T_{c}$ during a phase transition. While the latent parameter is different for each physical model, the reconstruction loss can be used as a universal parameter to identify phase transitions. To summarize, without any knowledge of the Ising model and its order parameter, but sample configurations, we can find a good estimation for its order parameter and the occurrence of a phase transition.

It is a priori not clear how to determine the number of latent neurons in the creation of the neural network of the autoencoder. Due to the lack of theoretical groundwork, we find the optimal number by experimenting. If we expand the number of latent dimensions by one, see Fig. 5, the results of our analysis only change slightly. The second parameter contains a lot less information compared to the first, since it stays very close to zero. Hence, for the Ising model, one parameter is sufficient to store most of the information of the latent representation.

While the ferromagnetic Ising model serves as an ideal starting ground, in the next step we are interested in models where different sites in the samples contribute in a different manner to the order parameter. We do this in order to show that our model is even sensitive to structure on the smallest scales. For the magnetization in the ferromagnetic Ising model, all spins contribute with the same weight. In contrast, in the antiferromagnetic Ising model, neighboring spins contribute with opposite weight to the order parameter (6).

Again the variational autoencoder manages to capture the traditional order parameter. The staggered magnetization is strongly correlated with the latent parameter, see Fig. 4. The three clusters in the latent representation make it possible to interpret different phases. Furthermore, we see that all three averaged quantities - the magnetization, the latent parameter and the reconstruction loss - can serve as indicators of a phase transition.

Fig. 6 demonstrates the reconstruction from the latent parameter. In the first row we see the reconstruction from samples of the ferromagnetic Ising model, the latent parameter encodes the whole spin order in the ordered phase. Reconstructions from the antiferromagnetic Ising model are shown in the second and third row. Since the reconstructions clearly show an antiferromagnetic phase, we infer that the autoencoder encodes the spin samples even to the most microscopic level.

IV.2 XY Model

In the XY model we examine the capabilities of a variational autoencoder to encode models with continuous symmetries. In models like the Ising model, where discrete symmetries are present, the autoencoder only needs to learn a discrete set, which is often finite, of possible representations of the symmetry broken phase. If a continuous symmetry is broken, there are infinitely many possibilities of how the ordered phase can be realized. Hence, in this section we test the ability of the autoencoder to embed all these different realizations into latent variables.

The variational autoencoder handles this model equally well as the Ising model. We find that two latent parameters model the phase transition best. The latent representation in the middle of Fig. 7 shows the distribution of various states around a central cluster. The radial symmetry in this distribution leads to the assumption that a sensible order parameter is constructed from the $L^{2}$-norm of the latent parameter vector. In Fig. 7, one sees the correlation between the magnetization and the absolute value of the latent parameter vector. Averaging the samples for the same temperature hints to the facts that the latent parameter and the reconstruction loss can serve as an indicator for the phase transition.

V Conclusion

We have shown that it is possible to observe phase transitions using unsupervised learning. We compared different unsupervised learning algorithms ranging from principal component analysis to variational autoencoders and thereby motivated the need for the upgrade of the traditional autoencoder to a variational autoencoder. The weights and latent parameters of the variational autoencoder approach are able to store information about microscopic and macroscopic properties of the underlying systems. The most distinguished latent parameters coincide with the known order parameters. Furthermore, we have established the reconstruction loss as a new universal indicator for phase transitions. We have expanded the toolbox of unsupervised learning algorithms in physics by powerful methods, most notably the variational autoencoder, which can handle nonlinear features in the data and scale very well to huge datasets. Using these techniques, we expect to predict unseen phases or uncover unknown order parameters, e.g. in quantum spin liquids. We hope to develop deep convolutional autoencoders which have a reduced number of parameters compared to fully connected autoencoders and can also promote locality in feature selection. Furthermore, since there exists a connection between deep neural networks and renormalization group Mehta and Schwab (2014), it may be helpful to employ deep convolutional autoencoders to further expose this connection.

Acknowledgments We would like to thank Timo Milbich, Björn Ommer, Michael Scherer, Manuel Scherzer and Christof Wetterich for useful discussions. We thank Shirin Nkongolo for proofreading the manuscript. S.W. acknowledges support by the Heidelberg Graduate School of Fundamental Physics.