Mutual information of spin systems from autoregressive neural networks
We describe a new direct method to estimate bipartite mutual information of a
classical spin system based on Monte Carlo sampling enhanced by autoregressive
neural networks. It allows studying arbitrary geometries of subsystems and can
be generalized to classical field theories. We demonstrate it on the Ising
model for four partitionings, including a multiply-connected even-odd division.
We show that the area law is satisfied for temperatures away from the critical
temperature: the constant term is universal, whereas the proportionality
coefficient is different for the even-odd partitioning.
Piotr Białas, Piotr Korcyl, Tomasz Stebel
# Mutual information of spin systems from autoregressive neural networks

###### Abstract

We describe a direct approach to estimate bipartite mutual information of a classical spin system based on Monte Carlo sampling enhanced by autoregressive neural networks. It allows studying arbitrary geometries of subsystems and can be generalized to classical field theories. We demonstrate it on the Ising model for four partitionings, including a multiply-connected even-odd division. We show that the area law is satisfied for temperatures away from the critical temperature: the constant term is universal, whereas the proportionality coefficient is different for the even-odd partitioning.

Mutual Information, Hierarchical Autoregressive Neural Networks, Monte Carlo simulations, Ising model _Introduction_: The discovery of topological order in quantum many-body systems [1] initiated a very fruitful exchange of ideas between solid-state physics and information theory. Many new theoretical tools developed to quantitatively describe the flow of information or the amount of information shared by different parts of the total system have been employed in the studies of physical systems [2]. Among these tools, quantum entanglement entropy or mutual information and their various alternatives were found to be particularly useful [3]. With their help, it was understood that many new phases of matter respect the same set of symmetries but differ in long-range correlations quantified by bipartite, tripartite, or higher information-theoretic measures such as mutual information. Turning this fact around, it is expected that calculating mutual information can provide hints about the topological phase of the system, playing a role similar to order parameters in the usual Landau picture of phase transitions (see for example [4]). Indeed, such quantities are not only useful for theoretical understanding but are also measurable observables in experiments. For example, in Ref. [5] the mutual information in a quantum spin chain was measured demonstrating the area law [6] governing the scaling of mutual information with the volume of the bipartite partition. The law [7] claims that the discussed quantities, the mutual information being one of them, scale with the boundary area separating the two parts of the system, instead of the expected scaling with the volume. The increased interest comes from the fact that the black hole entropy was shown a long time ago to follow similar law: the black hole entropy depends only on its surface and not on the interior. To a large surprise, this connection was recently made explicit with the realization that a certain quantum spin model in \(0\) spatial dimensions called the SYK model [8] is dual to a black hole via the gravity/gauge duality [9]. From this perspective, quantum entanglement entropy or quantum mutual information which can be defined and calculated on both sides of this correspondence became even more attractive from the theoretical point of view. Therefore there is a strong pressure to develop computational tools able to estimate such quantities, which in itself is, however, very difficult in the general case.

Although an important part of Condense Matter physics [6; 10; 11; 12] aims at studying entanglement entropy and quantum mutual information in quantum spin systems, already in classical spin systems the Shannon and Renyi entropies were found to follow the area law [13]. With the interface boundary sizes up to \(64\) spins the coefficients of the area law were precisely determined and their universality was verified using different lattice shapes. Later, the mutual information of two halves of the classical Ising model on an infinitely long cylinder was calculated using the transfer matrix approach [14]. Both quantities were discussed with more precision in [15] using a different method, called the bond propagation algorithm.

In this work, we propose a new approach to directly estimate bipartite mutual information. It is based on the incorporation of machine learning techniques into Monte Carlo simulation algorithms. The calculations are stochastic and yield provably exact results within their statistical uncertainties. The main advantage of the method is its flexibility, allowing its application: to any geometry of the partitioning, to any statistical system with a finite number of degrees of freedom, in an arbitrary number of space dimensions.

We consider a classical system of spins which we divide into two arbitrary parts \(A,B\). In this case, the Shannon mutual information is defined as

\[I=\sum_{\mathbf{a},\mathbf{b}}p(\mathbf{a},\mathbf{b})\log\frac{p(\mathbf{a},\mathbf{b})}{p(\mathbf{a})p(\mathbf{b})}, \tag{1}\]

where a particular configuration \(\mathbf{s}\) of the full model has parts \(\mathbf{a}\) and \(\mathbf{b}\), and where the Boltzmann probability dis
tribution of states, depending on inverse temperature \(\beta\), is given by (we omit the explicit dependence of \(Z\) on \(\beta\)),

\[p(\mathbf{a},\mathbf{b})=\frac{1}{Z}e^{-\beta E(\mathbf{a},\mathbf{b})},\quad Z =\sum_{\mathbf{a},\mathbf{b}}e^{-\beta E(\mathbf{a},\mathbf{b})} \tag{2}\]

and

\[p(\mathbf{a})=\sum_{\mathbf{b}}p(\mathbf{a},\mathbf{b}),\ \ p(\mathbf{b})=\sum_{ \mathbf{a}}p(\mathbf{a},\mathbf{b}) \tag{3}\]

are probability distributions of subsystems. We shall use the same symbol \(p\) for all probability distributions defined on different state spaces distinguishing them by the arguments. In the above expressions, the summation was performed over all configurations of subsystems \(A\) or \(B\). Inserting Eq. (2) and Eq. (3) into Eq. (1) we obtain

\[I=\log Z-\sum_{\mathbf{a},\mathbf{b}}p(\mathbf{a},\mathbf{b}) \Big{[}\beta E(\mathbf{a},\mathbf{b})+\\ +\log Z(\mathbf{a})+\log Z(\mathbf{b})\Big{]}, \tag{4}\]

where \(Z(\mathbf{a})=\sum_{\mathbf{b}}e^{-\beta E(\mathbf{a},\mathbf{b})}\) and \(Z(\mathbf{b})=\sum_{\mathbf{a}}e^{-\beta E(\mathbf{a},\mathbf{b})}\). Please note that in typical Monte Carlo approaches the partition functions \(Z\), \(Z(\mathbf{a})\), and \(Z(\mathbf{b})\) are not available.

_Method:_ Below we argue that \(I\) can be obtained from a Monte Carlo simulation enhanced with autoregressive neural networks (ANNs). It was recently shown that ANNs can be used to approximate the Boltzmann probability distribution \(p(\mathbf{a},\mathbf{b})\) for spin systems and provide a mean to simulate from this approximate distribution [16; 17; 18; 19]. Let us call this approximate distribution \(q_{\theta}(\mathbf{a},\mathbf{b})\), \(\theta\) stands here for the parameters of the neural network that are to be tuned so that \(q_{\theta}\) is as close to \(p\) as possible under an appropriate measure, typically the backward Kullback-Leibler divergence

\[D_{\mathrm{KL}}(q_{\theta}|p)=\sum_{\mathbf{a},\mathbf{b}}q_{\theta}(\mathbf{a },\mathbf{b})\,\log\left(\frac{q_{\theta}(\mathbf{a},\mathbf{b})}{p(\mathbf{a },\mathbf{b})}\right). \tag{5}\]

The formula Eq. (4) can be rewritten in terms of averages with respect to the distribution \(q_{\theta}\)

\[I=\log Z-\frac{1}{Z}\beta\left\langle\hat{w}(\mathbf{a},\mathbf{ b})E(\mathbf{a},\mathbf{b})\right\rangle_{q_{\theta}(\mathbf{a},\mathbf{b})}+\\ +\frac{1}{Z}\left\langle\hat{w}(\mathbf{a},\mathbf{b})\log Z( \mathbf{a})\right\rangle_{q_{\theta}(\mathbf{a},\mathbf{b})}+\\ +\frac{1}{Z}\left\langle\hat{w}(\mathbf{a},\mathbf{b})\log Z( \mathbf{b})\right\rangle_{q_{\theta}(\mathbf{a},\mathbf{b})}, \tag{6}\]

where the importance ratios are defined as

\[\hat{w}(\mathbf{a},\mathbf{b})=\frac{e^{-\beta E(\mathbf{a},\mathbf{b})}}{q_{ \theta}(\mathbf{a},\mathbf{b})}. \tag{7}\]

The crucial feature of ANN-enhanced Monte Carlo is that contrary to standard Monte Carlo, we can estimate directly the partition functions \(Z\), \(Z(\mathbf{a})\) and \(Z(\mathbf{b})\) (see Appendix A for details). In this way, we have expressed the mutual information \(I\) only through averages with respect to the distribution \(q_{\theta}\). It can be now estimated by sampling from this distribution. The procedure of sampling configurations from approximate probability distribution provided by neural network together with reweighting observables with importance ratios was proposed in Ref. [20] and named Neural Importance Sampling (NIS). This paper reports on the first application of this technique to information theory observables.

Autoregressive neural networks rely on the product rule _i.e._ factorization of \(q_{\theta}\) into the product of conditional probabilities

\[q_{\theta}(\mathbf{a},\mathbf{b})=\prod_{i=1}^{L^{2}}q_{\theta}(s^{i}|s^{1},s^ {2},\ldots,s^{i-1}). \tag{8}\]

Due to the fact that the labeling of spins in Eq. (8) is arbitrary, we can choose it in such a way that we first enumerate all spins from part \(A\), \(\mathbf{a}=(s^{1},s^{2},\ldots,s^{n_{A}})\) and only afterward all spins from part \(B\), \(\mathbf{b}=(s^{n_{A}+1},\ s^{n_{A}+2},\ldots,s^{n_{A}+n_{B}})\). We then obtain

\[q_{\theta}(\mathbf{s})\equiv q_{\theta}(\mathbf{a},\mathbf{b})=q_{\theta}( \mathbf{a})q_{\theta}(\mathbf{b}|\mathbf{a}), \tag{9}\]

with

\[q_{\theta}(\mathbf{a})=\prod_{i=1}^{n_{A}}q_{\theta}(s^{i}|s^{1},s^{2},\ldots, s^{i-1}), \tag{10}\]

and

\[q_{\theta}(\mathbf{b}|\mathbf{a})=\prod_{i=1}^{n_{B}}q_{\theta}(s^{n_{A}+i}|s^ {n_{A}+1},s^{n_{A}+2},\ldots,s^{n_{A}+i-1},\mathbf{a}). \tag{11}\]

This feature of ANN stands behind the fact that we can readily estimate \(\log Z(\mathbf{a})\) and \(\log Z(\mathbf{b})\) required by Eq. (6) directly. In this respect, the ANN approach differs from _normalizing flows_ used to approximate continuous probability distributions and employed recently in the context of the lattice field theory [21; 22; 23; 24; 25; 26; 27]. There, the conditional probability \(q_{\theta}(\mathbf{b}|\mathbf{a})\) as well the marginal distribution \(q_{\theta}(\mathbf{a})\) are not so easily available.

_Results:_ We leave the technical details of the evaluation of the terms in Eq. (6) to the Appendix A and now provide an example of its application. We demonstrate it on the Ising model on a periodic \(L\times L\) lattice, with ferromagnetic, nearest-neighbor interactions, defined by the Hamiltonian

\[E(\mathbf{a},\mathbf{b})=-\sum_{\langle i,j\rangle}s^{i}s^{j}, \tag{12}\]

where \(s^{i}\in\{1,-1\}\). We consider the following divisions and respective mutual information observables: "strip" geometry: the system is divided into equal rectangular subsystems; "square" geometry: subsystem \(A\) is the
square of size \(\frac{1}{2}L\times\frac{1}{2}L\); "quarter" geometry: \(A\) is rectangular of size \(\frac{1}{4}L\times L\). We show them schematically in the three sketches on the left in Fig. 1. We note that all of them have the same length of the border between the subsystems, i.e. \(2L\) (as we used periodic boundary conditions) although may they differ in their volumes. In the discussion below, we refer to those three partitionings as "block" partitionings. In addition, we also consider a division, called "chessboard", where the system is divided using even-odd labeling of spins. In this case, the boundary between spins is \(L^{2}\), i.e. every spin of the system is at the boundary between parts \(A\) and \(B\). Note that in the chessboard partitioning the subsystems are not simply connected, as is usually considered in the Literature.

We investigate a wide range of system sizes, reaching \(L=66\) or \(130\) (for some \(\beta\)). We combine results obtained by the Variational Autoregressive Network (VAN) approach of Ref. [16] and our recently proposed modification called Hierarchical Autoregressive Network (HAN) algorithm [18]. Both approaches are applied to several divisions, as summarized in Table 1. We describe the details of the VAN and HAN architectures and the quality of their training in Appendix B. In Appendix C we discuss the details of mutual information calculation for the chessboard partitioning. Simulations are performed at 13 values of the inverse temperature: from 0.1 to 0.4 with a step 0.05, 0.44, and from 0.5 to 0.9 with a step 0.1. For each \(\beta\) we collect the statistics of at least \(10^{6}\) configurations and estimated the statistical uncertainties using the jackknife resampling method. In all cases, the relative total uncertainty of \(I\) combining the systematic and statistical uncertainties is of the order of 1 %, or smaller.

In Fig. 2 we show \(I\) as a function of \(L\) for the three block partitionings and for three representative inverse temperatures \(\beta\): 0.1 (high-temperature regime), 0.44 (very close to critical temperature), and 0.9 (low-temperature regime). We clearly observe a linear dependence on the system size with a strongly \(\beta\)-dependent slope and intercept. Similar behavior can be observed for the chessboard geometry when one plots \(I\) as the function of \(L^{2}\) - see Fig. 3. The emerging picture, supporting the conjectured area law, is that \(I\) can be written in a compact form as

\[I_{geom}(\beta,L)=\alpha_{geom}(\beta)B(L)+r_{geom}(\beta), \tag{13}\]

where the variable \(B(L)\) corresponds to the length of the boundary between the two parts \(A\) and \(B\): \(B(L)=L^{2}\) for the chessboard partitioning and \(B(L)=2L\) for other partitionings considered in this work. Eq. (13) has two

\begin{table}
\begin{tabular}{|c|c|c|} \hline Partitioning & VAN & HAN \\ \hline strip & 8,12,16, & 10,18,34, \\  & 20,24,28 & 66,130(for \(\beta\leq 0.3\)) \\ \hline square & & 10,18,34, \\  & 66,130(for \(\beta\leq 0.3\)) \\ \hline quarter & 8,12,16, & \\  & 20,24,28 & \\ \hline chessboard & 8,12,16,20, & \\  & 24,28,32 & \\ \hline \end{tabular}
\end{table}
Table 1: The system sizes \(L\) which were considered for given partitioning type. Note that HAN requires total system size \(L=2^{n}+2\), where \(n\) would be the number of levels in the hierarchy.

Figure 1: Four considered partitioning geometries. Periodic boundary conditions are applied. The small blocks in the chessboard partitioning represent single spins.

Figure 3: General dependence of \(I\) on \(L^{2}\) for the chessboard partitioning at three inverse temperatures: \(\beta=0.10\), \(\beta=0.44\) and \(\beta=0.90\). Lines correspond to attempted fits using the area law Ansatz. Statistical uncertainties are shown but are much smaller than the symbol size.

Figure 2: General dependence of \(I\) on \(L\) for various geometries and three representative inverse temperatures: \(\beta=0.10\) in the disordered phase, \(\beta=0.44\) close to the phase transition and \(\beta=0.90\) in the ordered phase. Lines correspond to attempted fits using the area law Ansatz.
parameters, \(\alpha\) and \(r\), which, as we denoted, depend on \(\beta\) and can differ for different partitioning geometries.

We expect that Eq. (13) is valid as long as finite volume effects in \(I\) are small. In general, they may depend on three length scales present in the system: the size of the whole system \(L\), the size of the smallest subsystem, and correlation length \(\xi\) determined by the temperature of the system. Finite volume effects should vanish when \(\xi\) is smaller than other scales and they may depend on the geometry of the partitioning. A closer look at the data indeed reveals additional contributions to mutual information which spoil the area law (13). To see this we plot in Fig. 4 the \(\sqrt{\chi^{2}/{\rm DOF}}\) of the fits, which were performed assuming relation Eq. (13). Three regions of inverse temperatures can be easily defined. At small \(\beta\) the Ansatz Eq. (13) describes all system sizes since \(\sqrt{\chi^{2}/{\rm DOF}}\approx 1\), suggesting that finite volume effects are smaller than the statistical uncertainties. A similar situation occurs for large \(\beta\), where again the fit involved all available system sizes. Contrary, in the region close to the phase transition, where the correlation length \(\xi\) is the largest, the fits have clearly bad quality, \(\sqrt{\chi^{2}/{\rm DOF}}\gg 1\). This picture is further confirmed by the fact that the quality of the fit improves when we discard smaller system sizes, \(L<L_{min}\), as shown for the strip partitioning in the inset.

The values of \(\sqrt{\chi^{2}/{\rm DOF}}\) close to critical temperature are particularly large for the chessboard partitioning. However, since the errors of the individual points are much smaller than for other partitionings (due to the simplified calculation of \(Z({\bf a})\) - see Appendix C), we refrain from drawing conclusions about the size of finite-size effects in this geometry compared to block partitionings.

_Discussion of \(\alpha(\beta)\) and \(r(\beta)\) coefficients:_ For \(\beta\) far from critical inverse temperature we can reliably describe our data with the area law Ansatz Eq. (13). This allows us to extract the coefficients \(\alpha\) and \(r\) for the four different partitionings and compare them. We show the results in Fig. 5. In the main figures, we show the dependence on \(\beta\) whereas in the insets we show the difference between the "strip" and "square" partitionings (the differences between other partitionings look similar). We have shown only the values which were obtained from fits with \(\sqrt{\chi^{2}/{\rm DOF}}\lesssim 1\), as to be sure that the postulated dependence Eq. (13) is indeed reflected in the data. This means \(\beta\leq 0.35\) and \(\beta\geq 0.5\) for block partitionings, \(\beta\leq 0.25\) and \(\bar{\beta}\geq 0.5\) for chessboard.

In the left panel, we show the coefficient \(\alpha\) as a function of \(\beta\). Qualitative behavior of \(\alpha(\beta)\) seems to be universal: it goes to \(0\) at \(\beta=0\) and \(\beta\to\infty\) and rises around the critical temperature. The data clearly show that the chessboard partitioning yields different values than the three other possibilities which seem to be compatible with each other. The differences shown in the inset are indeed compatible with zero within their uncertainties. Therefore, we conclude that in this range of \(\beta\) the partitioning geometry does not influence the mutual information (when block partitionings are considered). In the right panel, we show the value of \(r\). In that case, all four partitionings give the same result: \(0\) for \(\beta<\beta_{c}\) and \(\log 2\) for \(\beta>\beta_{c}\). This seems to be quite a universal behavior [13].

_Conclusions and prospects:_ In this work, we have provided a numerical demonstration that the Shannon bipartite mutual information can be readily obtained from the Neural Importance Sampling algorithm for the classical Ising model on a square lattice. Our approach allows for studies of different partitionings and we provided comparisons of the mutual information estimated for four geometries. We successfully exploited the hierarchical algorithm (HAN) to reach larger system sizes than achievable using standard Variational Autoregressive Networks (VAN). We discussed the validity and universality of the area law for different partitionings. We found that at low and high temperatures the area law is satisfied whereas such an Ansatz does not describe our data in the vicinity of the phase transition, supposedly due to inherent finite volume effects.

Our proposal can be applied to many other spin systems with short- and long-ranged interactions, such as various classes of spin-glass systems. Long-ranged interactions prohibit the application of the HAN algorithm limiting, at the moment, the available system sizes down to \(L\sim 30\).

We believe that by exploiting the Feynman path integral quantization prescription one may use the approach based on autoregressive networks to estimate also the entanglement entropy in quantum spin systems. In such a picture, a \(D\) dimensional quantum system is described by a \(D+1\) statistical system, where the machine learning enhanced Monte Carlo is applicable. In particular, thanks to the replica method [28], one can directly express the Renyi entropies in terms of partition functions [29] readily obtainable in the NIS. In this approach sys

Figure 4: Quality of fits involving all available data points and assuming the area law functional form of the mutual information for the four geometries shown in Figure 1: \(\sqrt{\chi^{2}/{\rm DOF}}\) plotted as a function of the inverse temperature. In the inset, we show how the fit improves as we restrict the system sizes included in the fit to be bigger or equal to \(L_{\rm min}\).
tems in any space-time dimension can be studied, with the obvious limitation that performance is limited by the total number of sites in the system, hence reducing the available volumes in higher dimensions. Still, due to its straightforwardness, it should be seen as a valuable alternative to studying information-theoretic properties of one-, two-, or three-dimensional quantum systems.

One should also observe that our proposal in principle can work unaltered for systems with continuous degrees of freedom. Neural generative networks have been already discussed in the context of \(\phi^{4}\) classical field theory [30] as well as \(U(1)\), \(SU(N)\)[31], and the Schwinger model gauge theories [32]. In all these cases, one would introduce a conditional normalizing flow or another proposal that would give access to conditional probabilities (see for example Ref. [33]). For instance, a hierarchical construction similar to the HAN algorithm employing normalizing flows could be used to simulate the \(\phi^{4}\) classical field theory. In this way, conditional probabilities would be naturally introduced and could be used to calculate the mutual information for some specific partitioning. The combination of the proposed method of measuring information-theoretic quantities with the recent advancements [27] in the machine learning enhanced algorithms for simulating four-dimensional Lattice Quantum Chromodynamics (LQCD) can open new ways of investigating this phenomenologically important theory, for instance by studying quantum correlations and entanglement of the QCD vacuum.

###### Acknowledgements.

 Computer time allocation 'plpng' on the Prometheus and ARES supercomputers hosted by AGH Cyfronet in Krakow, Poland was used through the Polish PLGRID consortium. T.S. kindly acknowledges the support of the Polish National Science Center (NCN) Grants No. 2019/32/C/ST2/00202 and 2021/43/D/ST2/03375 and support of the Faculty of Physics, Astronomy and Applied Computer Science, Jagiellonian University Grant No. LM/23/ST. P.K. acknowledge that this research was partially funded by the Priority Research Area Digiworld under the program Excellence Initiative - Research University at the Jagiellonian University in Krakow. P.K. and T.S thank Alberto Ramos and the University of Valencia for their hospitality during the stay when part of this work was performed and discussed. We also acknowledge very fruitful discussions with Leszek Hadasz on entanglement entropy.
