Unsupervised learning of topological phase transitions using the Calinski-Harabaz index

Here, we explain the DM method of Rodriguez-Nieva and Scheurer nature . Assume that we have $M$ configurations $\{x_{l}\}$, where $l=1,\dots M$. The connectivity between $x_{l}$ and $x_{l^{{}^{\prime}}}$ is denoted by the elementary Gaussian kernel

The normalization of $K_{\varepsilon}(x_{l},x_{l^{{}^{\prime}}})$ can be done by performing the sum over $l^{{}^{\prime}}$,

We can also perform the normalization along the direction of $l$ (i.e., the sum over $l$). The eigenvalue equation of the diffusion matrix $P_{l,l^{{}^{\prime}}}$ is $P\psi_{k}=\lambda_{k}\psi_{k}$, where $\psi_{k}$’s are the right eigenvectors, with the corresponding eigenvalues $\lambda_{k}\leq 1$, for $k=0,1,\dots,m-1$.

To find the phase transition point, Ref. nature and its earlier preprint version prep use two different ways, respectively:
(a) Intersection of the mean distance of cluster centers $\overline{D}/2n$, where $n$ is the number of clusters, and the dispersion $\overline{\sigma}$ of the data points in each cluster. The quantities $\overline{D}/2n$ and $\overline{\sigma}$ are obtained from the scatter plot, where $n$ is the number of topological sectors.

After obtaining the eigenvectors of the $P$ matrix given by

the authors project them unto a two-dimensional space, namely, a two-column matrix, then $\overline{D}/2n$ and the dispersion $\overline{\sigma}$ can be obtained from the scatter plot of the two column vectors or their modified version. The detailed application to the one-dimensional XY model is reproduced in Appendix A.
(b) Intersection of the mean fluctuation $\sigma_{\lambda}$ and gap of eigenvalues $\Delta\lambda$, where

and the gap of eigenvalues between the topological sectors $n$ and $n-1$,

We propose to use indices instead of intersections. Based on the first two leading eigenvectors $\Psi_{0}$ and $\Psi_{1}$ of the PCA, kernel-PCA (kPCA) or the second and third vectors $\Psi_{1}$ and $\Psi_{2}$ of the DM method, for a manually chosen cluster numbers $k$, the $ch$ index is given in terms of the following ratio:

where $B_{k}$ is the between-group dispersion matrix and $W_{k}$ is the within-cluster dispersion matrix and they are defined as follows,

where $N$ is the number of data points, $C_{q}$ is the set of points in cluster $q$, $c_{q}$ is the center of cluster $q$, and $c$ denotes the average center of all cluster centers $\{c_{q}\}$, and $n_{q}$ the number of points in cluster $q$. The $sc$ index of the $i$-th sample is

where $a(i)$ is the mean distance between sample $i$ and all other data points in the same cluster, $b(i)$ is the mean dissimilarity of point $i$ to some cluster $C$ expressed as the mean of the distance from $i$ to all points in $C$ (where $C\neq C_{i}$). The mean $sc=\sum sc(i)/N$ over all points of a cluster is a measure of how tightly grouped all the points in the cluster are. Sometimes $sc_{a}=\sum a(i)/N$ or $sc_{b}=\sum b(i)/N$ are also very useful 11index .

The performance of a total of 11 indices is shown in Appendix C and the results of indices $ch$, $dn$, $pbm$ and $Ii$ applied to the $q=8$ generalized XY model are all reasonable choices. Here, we choose two representative indices $ch$ and $sc$ as examples.

Fig. 1 draws the flowchart of the RNS method and the place of our modifications. Depending on whether or not we are considering topology, we perform dimensional reduction by DM or PCA (k-PCA) respectively, to get the features of the data, i.e., the eigenvectors of the diffusion matrix $P$ or the covariance matrix. Then we apply k-means to cluster the two-dimensional or higher dimensional scattering points. If we do not choose “ch”, then intersections of $\Delta\lambda$ and $\sigma_{\lambda}$, or intersections of $\bar{D}/2n$ and $\bar{\sigma}$ nature ; prep can be used instead.

To test whether the clustering is good or bad, if “$ch$” is chosen, then the peaks of $ch$ or its components $ch_{t}$ and $ch_{b}$ will be used directly. For the topological phase transition, the index $ch$ or its component $ch_{t}$ and $ch_{b}$ can give signatures of phase transition when it is not easy to determine the intersection of the RNS method or when there is no intersection.

The index $ch$ can also be applied to non-topological phase transitions, such as the order-disorder phase transition of the Ising model. This does not require any prior knowledge except for the configurations generated by e.g. Monte Carlo methods or from real experiments.

For topological phase transitions, such as the XY and generalized XY models, although this type of unsupervised learning is not similar to supervised learning (such as fully-connected layers, or convolutional neural networks), it still needs labels of configurations, i.e., the topological winding number and the number of possible phases. However, the label of the topological winding number does not mean the label of phases, and essentially, this method is still an unsupervised learning method.

Now we analyze the first example, i.e., the two-dimensional XY model on the square lattice. Firstly, we generate the configurations with five fixed winding number pairs at $\nu=(\nu_{x},\nu_{y})=(0,0)$, $(0,1)$, $(1,0)$, $(0,-1)$ and $(-1,0)$. For each fixed winding number pair, it should be noted that, for the two-dimensional geometry, the winding number component $\nu_{x}=1$ means that the spins in each row form a winding number of $1$ rather than just the spins on one row randomly selected. Cluster simulation algorithms wf ; SW are not suitable to update the spins because the global flips break the topological winding number easily and therefore the Metropolis algorithm is used while trying to rotate the spin vector with a very small step each time so as to preserve the winding number.

For each topological sector, we generate $m=500$ configurations. Combining all configurations $\{x_{l},l=1,\cdots,5m\}$ from all five sectors, we construct the kernel $K_{\varepsilon}$. The elements between $K_{\varepsilon}(x_{l},x_{l^{{}^{\prime}}})$ is defined in Eq. (1) and the normalized matrix $P_{\varepsilon}(x_{l},x_{l^{{}^{\prime}}})$ is obtained, obeying its eigenvalue equation $P\Psi_{k}=\lambda_{k}\Psi_{k}$. Using the scatter plot of the second and third leading eigenvectors $\Psi_{1}$ and $\Psi_{2}$ of $P_{\varepsilon}$, the $ch$ index are obtained and displayed in Figs. 3 (a)- (d).

There is a parameter $\varepsilon$ in the above matrices and in order to deduce consistent results, we need to make sure the results are converging for a finite range of $\varepsilon/\varepsilon_{0}$. We see that this is indeed the case for different values of $\varepsilon/\varepsilon_{0}$ in Figs. 3(a)-(d). For $\varepsilon/\varepsilon_{0}=2.5$, $3$, $3.5$ & $4$, $T_{c}$ is less than $0.9$, and then becomes $0.93$ when $\varepsilon/\varepsilon_{0}=4.5$, $5$, $5.5$, $6$, $6.5$ and $7$. Finally, the peaks move left and deviate from $0.93$ again when $\varepsilon/\varepsilon_{0}=8$ or larger. It should be noted that here, $10,000$ or more Monte Carlo steps, are used in order to reach the equilibrium of systems.

The estimated points are labeled by the circles in Fig. 3 (e) and they all distribute nearby or on the red lines representing the latest result critical temperature $T_{c}=0.8935$kao The intersections of $\bar{D}$ and $\bar{\sigma}$ are labeled by the gray regions in the critical regimes $T_{c}=0.9\pm 0.1$ from Ref. nature . It appears that it is easier to use $ch$ to locate the phase transition as we only need to identify the peak location. The results from Ref. nature have greater uncertainty than those by using the $ch$ index.

The histogram of our estimated $T_{c}$ is shown in Fig. 3 (f), which helps to determine the hyperparameter $\varepsilon/\varepsilon_{0}$ (see Sec. VI).

In Fig. 3 (g)-(i), the finite size effect of $T_{c}$ is checked using the $ch$ with two smaller sizes $L=8$ and $L=16$ with $\varepsilon/\varepsilon_{0}=5.5$, $6$ & $6.5$. Combining the estimated $T_{c}(L)$ with $L=8,16,32$, we use two different ways of (linear, exponential) extrapolation to get $T_{c}(\infty)$ in the thermodynamic limit. The results are $0.88(5)$ and $0.89(5)$, respectively.

Here, we study the second example , i.e., the pure XY model on the honeycomb lattice, as the critical point is known exactly at $T_{c}=1/\sqrt{2}\approx 0.707$ 0.7 .

The geometry of the honeycomb lattice is equivalent to the $8\times 8$ brick-wall lattice shown in Fig. 4, where every spin has three nearest-neighbor spins. To initialize the configurations with a fixed winding number ($\nu_{x}$, $\nu_{y}$)=(0, 1), the spins connected by solid gray lines are defined as forming $\nu_{y}=1$ in the vertical direction. Specifically, if we start a spin at position (2,1), and then go left to $\rightarrow$ (1,1) $\rightarrow$ (1,2) $\cdots$ (2,8) and go back to (2,1) through the red dashed lines, connecting the sites at the boundaries for periodic boundary conditions, the spins sweep an angle of $2\pi$ counter-clockwise.

It should be noted that two spins, such as (1,1) and (2,1), connected by the horizontal gray lines will contribute to the winding number $\nu_{y}$, and they also contribute to $\nu_{x}$. This poses a problem when fixing $\nu_{x}$ and $\nu_{y}$ independently. Fortunately, this problem can be solved. Specifically, in the first row, labeled by 1 in the vertical ($y$) direction, the relation of angles obeys $\theta_{(2,1)}-\theta_{(1,1)}=-(\theta_{(1,1)}-\theta_{(2,1)})\approx-(\theta_{(3,1)}-\theta_{(2,1)})$. During the simulation, small fluctuations are allowed if they do not break the winding number.

In Fig. 5, using configurations constrained in five topological sectors on the $32\times 32$ lattices, we find that the peaks are located at the exact value $T_{c}=\frac{1}{\sqrt{2}}\approx 0.7$ for several different values of $\varepsilon/\varepsilon_{0}$ in the interval [2,6] with intervals of 0.5. We also calculate the intersections by $\Delta\lambda$ and $\sigma_{\lambda}$ at $\varepsilon/\varepsilon_{0}=$ 5, 6, 7. The intersection can also arrive at $0.7$ when $\varepsilon/\varepsilon_{0}=7$, but not 5 and 6. It indicates that the $ch$ performs better than the intersection as the intersection method may not give a high confidence of the transition.

The $q=2$ GXY model has a richer phase diagram than the XY model and has an additional new phase mlpxy . Thus, it is a good candidate model to test our method away from the pure XY limit. The phase diagram is illustrated in Fig. 6 (a) on the $\Delta-T$ plane, and we also show the results from our unsupervised method and those from PCA as comparison. The symbols $N$, $F$, $P$ represent nematic, ferromagnetic and paramagnetic phases, respectively. The dashed lines are data from the MC simulations q2mc mainly of $L=50$ up to $L=300$. The color indicates the value of $ch$ index obtained from simulation with the system size $L=32$. We now discuss the $F$-$P$, $N$-$F$ and $N$-$P$ transitions as follows.

(i) The $F$-$P$ phase transition: Interestingly, the positions of the peaks by $ch$ are consistent with the dashed line of the phase boundary in the whole region $\Delta>0.4$. The index $ch$ performs very well where $\Delta$ is away from the pure XY model limit. The essential nature of the $F$-$P$ phase transition is still BKT.

(ii) The $N$-$P$ phase transition: In the regime $\Delta=0$, the $ch$ peaks around $T=0.7$, which is $0.2$ less than $0.9$. This discrepancy is likely due to the nature of half-vorticity in the $N$ phase. We can improve the result by limiting the configurations to have the half-winding number $\nu_{x(y)}=1/2$ as the topological constraint in our Monte Carlo simulation. The half-vortex physics has been discussed in Ref. q2mc .

Thus, we only consider ($\nu_{x}$,$\nu_{y}$) in the four types of combinations ($\pm 1/2$, 0) and (0, $\pm 1/2$). To form ($\nu_{x}=1/2$,0), the difference between a pair of spins located at the left most and right most boundaries is fixed as $\pi$ and we assign each spin using the Eq. (14). With the half vortex constraint, our results illustrated by the red symbols move closer to the dashed lines of the $N-P$ transition in Fig. 6 (a) than the results using the integer-value constraint.

Figs. 6 (b) and (c) show the details at $\Delta=0$ and $\Delta=0.1$, respectively. For $L=64$, we generate two groups of data and the peaks almost converge in the interval $[0.82,0.85]$, and the convergence is closer to $0.89$ than the result from the integer vortex constraint. For $\Delta=0.1$, the half vortex constraint gives better results at $T_{c}=0.8$.

(iii) The $N$-$F$ phase transition: When we realize that the topological constraint makes the peaks (features) of the distribution for the spin directions implicit in the $F$ and $P$ phases, we use the configurations without any constraint in this case. The results are labeled by purple symbols with legend ‘$pca$, $L32$’.

The behavior of $ch$ depends on the structure of the sample points, and sometimes $ch$ emerges as a sharp jump at the phase transition point such as in the Ising model discussed in the main text. Sometimes it is a local minimum value at the phase transition point. Take the $q=2$ GXY model at $T=0.5$ as an example, in Fig. 7, both $\ln(ch_{t})$ and $\ln(ch_{b})$ increase as a function of $\Delta$. However the difference $\alpha=\ln(ch_{t})-\ln(ch_{b})=-\mbox{gap}$ decreases first and then increases around $\Delta_{2}=0.22$, where the gap has a positive value.

According to the following equation:

clearly, $min(ch)\Rightarrow\max(\mbox{gap})$ and the local minimum of the $ch$ is at the location of $\Delta_{2}$.

For the $N$-$F$ phase transition in Fig. 6 (a), the PCA is stronger than the DM method. Here, we would like to give a physical discussion about this because such a difference is related to the advantage of the proposed method, and thus a more detailed discussion should be made.

The reasons for PCA being stronger for the $N-F$ transition can be explained as follows:

(i) The $N-F$ phase transition is not a topological phase transitionq8mc . The DM method designed here is to determine the topological phase transition. It is still interesting to see the $ch$ of DM by inputting Ising configurations without any topological constraint. In Fig. 9, using k-PCA and the DM method with complete same configurations, the signatures of phase transition emerge and the results are consistent at $T_{c}=2.3$. However, the signature of k-PCA method is clearer.

(ii) From the view of the data, it is also understandable that PCA is stronger. Without any topological constraint, in the nematic ($N$) phase, the spins prefer two dominant directions and their histograms obey a double peaks structure. In the Ferromagnetic ($F$) phase, one main direction remains and one peak emerges for the histograms. Fig. 8 (b1)-(b3) show the typical histograms in three phases. The main feature difference between the phases is obvious.

However, with a topological constraint, the spin directions are mainly distributed according to the winding numbers. For example, the spin angles in each row are $\theta_{x}=2\pi x/L$ with additional fluctuations, where $x$ and $L$ are the number index and total number, respectively, in each row. Therefore, the distribution sits almost in the range $[0,2\pi]$ as shown in Figs. 8(a1)-(a3). The feature difference of spin angles disappears when using the constraint. Therefore, the DM method cannot get transition points using the configurations with constraints.

For the $q=8$ GXY model q8mc , Fig. 10 (a) shows the global phase diagram using the values of the $ch$ index. The phases $N$, $F_{2}$, $F$, $P$ are obtained and the distributions are also shown.

In the new $F_{2}$ phase, the distribution of the spin vectors has 8 peaks, but is dominated by $4$ possible angles. The ‘X’ shape dashed lines are from Refs. q2mc ; q8mc . The orange color represents the values of $ch$ by the DM method.

Let us first discuss the $F$($F_{2}$)-$P$ transition. Clearly, for $F$-$P$ transition in the range $0.5<\Delta<1$, the peak positions of $ch$ align well with the dashed line. In the center of ’X’, the $F_{2}$ and $P$ transition are also consistent with MC result i.e., $T_{c}=0.5$.

However, we could not use the intersection of the cluster average distance $\bar{D}$ and within-cluster dispersion $\bar{\sigma}$ as described in Ref. nature because $\bar{D}$ does not vary too much and there is no intersection as shown in Figs. 10 (e)-(g). The five different colors therein represent the five typical topological sectors. However, fortunately, when zooming in the figures in Figs. 10 (g)-(i), we find that, near the transition region, the shape of a fixed cluster shrinks because the data points gather closer together, and hence the within-cluster dispersion $\bar{\sigma}$ is smaller. The index $ch=\frac{ch_{t}}{ch_{b}}$ thus develops a peak around $T_{c}=0.5$ as shown in Fig. 10 (c), with $ch_{t}$ and $ch_{b}$ also displayed in Fig. 10 (d). In contrast, Fig. 10 (b) shows that $sc$ cannot be used to signal the transition temperature $T_{c}$ because it evaluates the difference, $sc_{b}-sc_{a}$, but $sc_{a}\ll sc_{b}$ in spite of the fact that $sc_{a}$ has a local minimum.

It should be mentioned that for the $N$-$P$ transition, the use of either the integer or half vortex constraint is not suitable. Instead the $\nu_{x(y)}=1/8$ vortex constraint is needed in generating the configurations. Moreover, after comparison, we find that kPCA works better if we use $\theta_{i}$ as the input into the kPCA (using PCA does not work well). The kernel used for kPCA defined as a radial basis functional kernel $\exp(-\gamma||x/L-y/L||^{2})$, where $L=32$ is the system size, $x$ and $y$ are the configurations $\{\theta_{i}\}$ and $\gamma=1$ is the default value gp .

The other details of the $F$-$P$, $F_{2}$-$F$, $N$-$F_{2}$ and $N$-$P$ transitions will be discussed as follows.

For the $q=8$ GXY model, Figs. 11 show the distributions of spin vectors in the (a1) $N$, (a2) $F_{2}$ and (a3) $F$ phases with the integer-vortex constraint. The distributions of spin vectors in the $N$ and $F_{2}$ phase have no obvious differences. To distinguish the $N$ and $F_{2}$ phases, the constraints are canceled and the distributions are shown in Figs. 11 (b1)-(b3) respectively.

Figs. 12 (a)-(d) show the detail of phase transition $F$-$P$, $F_{2}$-$F$, $N$-$F_{2}$, respectively. For the $F$-$P$ transition, in Fig. 12 (a), fixing $\Delta$ in the interval $[0.5,1]$ at steps of $0.1$, the peaks of $ch$ are located at $0.5$ and $1$ respectively, completely matching the dashed lines in the global phase diagram in Fig.10.

In Fig. 12 (b), for the $F_{2}$-$F$ transition, by fixing $\Delta=0.8$, the peaks of $ch$ are located at $0.2$ with $L=8$, $16$, $32$, $48$ and $64$. The other values of $\Delta$ are not shown. In Fig. 12 (c), fixing $T$=$0.2$, $0.3$ and $0.4$, the positions $\Delta$ of local minimum of $ch$ located at $0.2$, $0.3$ and $0.4$ respectively.

In Fig.12 (d), for $\Delta=0$, using the $1/8$ vortex constraint results in the most accurate critical points. The peak position is at $0.9$ more closely to $1$ than $0.7$ by the integer vortex constraint. The error bars are obtained by the bootstrap method using 400 randomly chosen configurations between the total 2000 configurations and total of 20 bins.