Deep Learning-Enhanced Variational Monte Carlo Method for Quantum Many-Body Physics

TL;DR

This paper introduces an efficient deep neural network approach with a new optimization algorithm to accurately compute ground states in quantum many-body systems, demonstrating significant speed and accuracy improvements.

Contribution

The paper presents a novel importance sampling gradient optimization algorithm and an efficient convolutional DNN architecture for quantum many-body calculations.

Findings

Achieved ground-state energies matching exact solutions for 1D SU(N) spin chains.

Demonstrated improved training speed for deep neural networks in VMC.

Successfully computed loop correlation functions using the trained wave functions.

Abstract

Artificial neural networks have been successfully incorporated into variational Monte Carlo method (VMC) to study quantum many-body systems. However, there have been few systematic studies of exploring quantum many-body physics using deep neural networks (DNNs), despite of the tremendous success enjoyed by DNNs in many other areas in recent years. One main challenge of implementing DNN in VMC is the inefficiency of optimizing such networks with large number of parameters. We introduce an importance sampling gradient optimization (ISGO) algorithm, which significantly improves the computational speed of training DNN in VMC. We design an efficient convolutional DNN architecture to compute the ground state of a one-dimensional (1D) SU($N$) spin chain. Our numerical results of the ground-state energies with up to 16 layers of DNN show excellent agreement with the Bethe-Ansatz exact solution. Furthermore, we also calculate the loop correlation function using the wave function obtained. Our work demonstrates the feasibility and advantages of applying DNNs to numerical quantum many-body calculations.

Tables

Table 1. Table 1: Comparison of the variational energies of a one-hidden layer CNN for N site = 60 subscript 𝑁 site 60 N_{\text{site}}=60 SU( N 𝑁 N ) spin chains with exact Bethe Ansatz (BA) results. The kernel size is K = 19 𝐾 19 K=19 . Both results from the one-hot and the value encoding are included. The incertitudes on the VMC data are smaller than 10 − 4 superscript 10 4 10^{-4} .

$N$	BA	one-hot encoding	error [%]	value encoding	error [%]
2	-0.3868	-0.3865	0.06	-0.3866	0.05
3	-0.7038	-0.7027	0.16	-0.6885	2.18
4	-0.8258	-0.8237	0.25	-0.8020	2.88
5	-0.8854	-0.8827	0.32	-0.8586	3.03

Equations

$$H=\sum _{i=1}^{N_{\mathrm{site}}}{P}_{i,i+1},$$

$$|\Psi ⟩=\sum _{\{s_i\}}\Psi (s_1,s_2,...,s_{N_{\text{site}}})|s_1,s_2,...,s_{N_{\text{site}}}⟩,$$

$$A_{i,f^{\prime }}^{[1]}=\sigma (\sum _{k=1}^K\sum _{f=1}^{N_{\text{in}}}{W}_{k,f,f^{\prime }}^{[1]}{S}_{i+k,f}+b_{f^{\prime }}^{[1]}),$$

$$A_{i,f^{\prime }}^{[l]}=\sigma (\sum _{k=1}^K\sum _{f=1}^F{W}_{k,f,f^{\prime }}^{[l]}{A}_{i+k,f}^{[l-1]}+b_{f^{\prime }}^{[l]}),l=2,3,...,L,$$

$$\text{ln}\Psi (S)=\sum _{f=1}^F{a}_f\sum _{i=1}^{N_{\text{site}}}{A}_{i,f}^{[L]},$$

$$\begin{array}{r}{\partial }_w{E}_v\approx \sum _x{I}^0{E}_x^0{\partial }_w\text{ln}{\Psi }_x^0-\sum _x{I}^0{E}_x^0\sum _x{I}^0{\partial }_w\text{ln}{\Psi }_x^0,\end{array}$$

$$\begin{array}{r}{\partial }_w{E}_v\approx \sum _x{I}_x{E}_x{\partial }_w\text{ln}{\Psi }_x-\sum _x{I}_x{E}_x\sum _x{I}_x{\partial }_w\text{ln}{\Psi }_x,\end{array}$$

$$S_{m,n}=(-1{)}^{m-n}⟨(m\cdots n)⟩,$$

$$\begin{array}{rl}& N\theta (2\alpha )-\sum _{{\alpha }^{\prime }}\theta (\alpha -{\alpha }^{\prime })+\sum _{\beta }\theta (2\alpha -2\beta )+2\pi J_{\alpha }=0,\\ & \sum _{\alpha }\theta (2\beta -2\alpha )-\sum _{{\beta }^{\prime }}\theta (\beta -{\beta }^{\prime })+\sum _{\gamma }\theta (2\beta -2\gamma )+2\pi J_{\beta }=0,\\ & ⋮\\ & \sum _{\delta }\theta (2\zeta -2\delta )-\sum _{{\zeta }^{\prime }}\theta (\zeta -{\zeta }^{\prime })+2\pi J_{\zeta }=0\end{array}$$

$$N\theta (2\alpha )-\sum _{{\alpha }^{\prime }}\theta (\alpha -{\alpha }^{\prime })+2\pi J_{\alpha }=0,$$

$$J_j=-\frac{M_j-1}{2},-\frac{M_j-1}{2}+1,...,\frac{M_j-1}{2}.$$

$$E=N-2(N-1)M-2\sum _{l=1}^{(N-1)M}\text{cos}{K}_l.$$

$$E=N-\sum _{l=1}^{(N-1)M}\frac{1}{1/4+{\alpha }_l^2}.$$

$$\tilde{P}|\alpha ,\beta ⟩=\{\begin{array}{ll}|\beta ,\alpha ⟩& \alpha =\beta \\ -|\beta ,\alpha ⟩& \alpha \ne \beta \end{array}$$

Full text

Introduction — Over the past few years, artificial neural networks have been introduced to study quantum many-body systems Carleo and Troyer (2017); Nomura et al. (2017); Saito and Kato (2018); Cai and Liu (2018); Carleo et al. (2018); Liang et al. (2018); Choo et al. (2018a); McBrian et al. (2019). In their seminal paper Carleo and Troyer (2017), Carleo and Troyer proposed to represent the quantum many-body states by a Restricted Boltzmann Machine (RBM), which contains one visible and one hidden layer. The many-body wave function is represented by the visible layer after integrating out the hidden layer. The parameters in the RBM are trained in the VMC method. Following this work, the RBM and a few other networks have been applied to study several quantum many-body systems with good accuracy Carleo and Troyer (2017); Nomura et al. (2017); Saito and Kato (2018); Cai and Liu (2018); Carleo et al. (2018); Liang et al. (2018); Choo et al. (2018a); McBrian et al. (2019); Torlai et al. (2018); Glasser et al. (2018); Kaubruegger et al. (2018). So far, the networks that have been implemented in quantum physics studies are not deep and hence not powerful enough to represent more complicated many-body states. As a result, there has been no clear evidence that their performance far exceeds the more traditional state-of-the-art numerical algorithms sch as quantum Monte Carlo, density matrix reformalization group, or tensor networks, etc. To overcome this problem, deep neural networks (DNNs) have been suggested. Theoretical studies have shown that the DNNs can efficiently represent any tensor network states and most quantum many-body states, and possess distinct advantages over shallow networks Gao and Duan (2017); Clark (2018); Levine et al. (2019). In fact, the DNN-based deep learning has become the most successful model of many machine learning tasks and has dominated the field since 2012. The DNNs have been demonstrated to have comparable or superior performance in various tasks when compared to human experts, such as playing Atari games Mnih et al. (2015), Go Silver et al. (2016, 2017), and manipulating robots Gu et al. (2016); Levine et al. (2016), etc, and have led to rapid advances in artificial intelligence.

Despite of great interest, there has been relatively few works in applying DNNs to quantum many-body computations Sharir et al. (2019); Choo et al. (2019). This perhaps is due to the fact that applying DNNs to represent quantum many-body states faces two main challenges: inefficient optimization and insufficient information for the proper choice of DNN architectures. The former arises because a DNN typically contains a large number of parameters to train, while a proper choice of the architecture often requires physical insights about the nature of the quantum systems.

In this Letter, we propose an efficient convolutional DNN architecture to represent the ground state of quantum many-body systems. Most of the quantum systems consist of particles interacting with each other through finite range. Such a local interacting character can be ideally captured by the convolutional neural network (CNN). We have developed an importance sampling gradient optimization (ISGO) algorithm within VMC method, which significantly improves the optimization speed and hence enables us to utilize DNN architectures. Our method can take advantage of the automatic differentiation, which automatically computes gradient update via backward-propagation algorithm  Hecht-Nielsen (1992). We show that our method can be parallelized and take full advantage of the acceleration provided by graphics processing units (GPUs). The ISGO method achieves at least one order of magnitude speed up when trained on GPUs sup .

For benchmark purpose, we construct DNNs to represent the ground-state wave function of the 1D SU($N$) spin chain, which has an exact solution under the Bethe ansatz. We systematically test different DNN architectures with ISGO for systems with different complexity. Our numerical results for the ground-state energies of 1D SU($N$) spin chain show excellent agreement with the exact solutions. Furthermore, we are able to compute correlation functions which are extremely difficult to obtain by Bethe Ansatz. The convolutional DNN architecture we constructed for this work can be readily generalized to represent the ground states of other quantum many-body systems. The ISGO method can also be used to accelerate the computation based on general VMC methods.

Network Architectures — We consider a homogeneous 1D SU($N$) spin chain with $N_{\text{site}}$ spins, which is the simplest prototypical model with SU($N$) symmetry, governed by the Hamiltonian:

[TABLE]

where $P_{i,i+1}$ is the spin exchange operator which exchanges two neighboring spins: $P_{i,i+1}|a_{i},b_{i+1}\rangle=|b_{i},a_{i+1}\rangle$. This model can describe the behavior of 1D strongly interacting quantum spinor gases Deuretzbacher et al. (2014); Volosniev et al. (2014); Yang et al. (2015); Yang and Pu (2016), and has attracted significant attentions both experimentally and theoretically Wu et al. (2003); Honerkamp and Hofstetter (2004); Gorshkov et al. (2010); Taie et al. (2012); Pagano et al. (2014); Scazza et al. (2014); Zhang et al. (2014); Hofrichter et al. (2016). Here we will use the DNN to represent the ground-state wave function of this model.

A general state takes the form

[TABLE]

where each $s_{i}$ represents one of the $N$ spin states for the SU($N$) model. The goal is to build a network that takes the basis state $|\{s_{i}\}\rangle$ as input and compute the ground state wave function $\Psi(\{s_{i}\})$ such that the energy functional $\langle\Psi|H|\Psi\rangle/\langle\Psi|\Psi\rangle$ is minimized. The first step is to encode the input basis state into a 2D tensor $S_{j,\beta}$, where the first and the second indices $j$ and $\beta$ represent the spatial site and the local spin state, respectively. In this work, we consider two kinds of state encodings: value encoding which encodes each spin state into a number and one-hot encoding which encodes the spin state into a one-hot boolean vector. The tensor $S$ is fed into the DNN as an input. The output of the first hidden layer, which follows immediately after the input layer, is given by:

[TABLE]

where $A^{[1]}$ is the activation (or feature map) of the first hidden layer, $\sigma(x)=\text{max}(x,0)$ is the Rectified Linear Units (ReLU) activation function, which has been demonstrated to outperform traditional sigmoid activation function for DNNs Glorot et al. (2011), $W^{[1]}$ is a 3D tensor of shape ($K$, $N_{\text{in}}$, $F$), where $K$ is the convolution kernel size, $F$ the number of channels of the hidden layer, and $N_{\text{in}}$ the number of the channels of the input layer which is 1 for value encoding and $N$ for one-hot encoding. In this work, we use the same $K$ and $F$, which determines the width of the network, for every hidden layer. $b^{[1]}$ is a bias vector of size $F$. The output from the remaining hidden layers are given by:

[TABLE]

where $L$ is the total number of hidden layers that determines the depth of the network, $W^{[l]}$ a 3D tensor of shape $[K,F,F]$, and $b^{[l]}$ a bias vector of size $F$. After the last hidden layer, its output is summed along the spatial dimension, and followed by a single fully connected layer to give the final output of the network

[TABLE]

where $a$ is a weight vector of size $F$. The full structure of the network is illustrated in Fig. 1. Each magenta rectangular object corresponds to a convolutional filter. We use periodic padding for each convolutional layer to enforce periodic boundary condition. This network is fully convolutional Long et al. (2014), which means the network architecture is compatible with different system sizes, and we can easily do transfer learning. Here, $W^{[l]}$, $b^{[l]}$, and $a$ are the network parameters that need to be optimized. The total number of parameter is roughly $KF^{2}L$.

Importance Sampling Gradient Optimization — Before introducing the ISGO method, we first revisit the conventional gradient optimization method in VMC. The wave function $\Psi(\{w\})$ is encoded by the set of network parameters $w\in\{W^{[l]},b^{[l]},a\}$. In every iteration step, $N_{\text{sample}}$ quantum states following distribution $P^{0}_{x}\propto|\Psi^{0}_{x}|^{2}$ are sampled using a Markov chain. Here $\Psi^{0}$ is the input wave function from the previous step, and $x$ indexes the sampled states. The variational energy functional $E_{v}(\{w\})=\frac{\bra{\Psi^{0}(\{w\})}H\ket{\Psi^{0}(\{w\})}}{\braket{\Psi^{0}(\{w\})}{\Psi^{0}(\{w\})}}$ is then computed. To minimize $E_{v}(\{w\})$, the network parameters are updated as $w\leftarrow w-\alpha\partial_{w}E_{v}$, where the “learning rate” $\alpha$ is a small parameter opt . In our work, we use the Adam Kingma and Ba (2014), a variant of the Stochastic Gradient Descent algorithm. Here, $\partial_{w}E_{v}$ is approximated by the variational wave function $\Psi^{0}$ with the $N_{\text{sample}}$ samples

[TABLE]

where $E^{0}_{x}=\sum_{x^{\prime}}H_{x,x^{\prime}}\Psi^{0}_{x^{\prime}}/\Psi^{0}_{x}$ is the local energy under $\Psi^{0}$ and $I^{0}=2/N_{\text{sample}}$. After the parameters $w$ are updated, the wave function changes from $\Psi^{0}$ to $\Psi$ which serves as the input for the next iteration step, where a new set of states are sampled based on $\Psi$ and the previously sampled states based on $\Psi^{0}$ are discarded.

Inspired by the off-policy policy gradient method in Reinforcement Learning (RL) Meuleau et al. (2001); Sutton and Barto (2018), we develop an efficient importance sampling gradient optimization (ISGO) method that utilizes the mismatched samples, as shown in Fig. 2(a). The key is to renormalize the distribution of those mismatched samples to $|\Psi_{x}|^{2}$ by multiplying the local energies and derivatives in Eq. (5) with importance sampling factors:

[TABLE]

where $E_{x}$ is the local energy under the last updated wave function $\Psi$, and $\frac{I_{x}}{I^{0}}=\frac{P_{x}}{P^{0}_{x}}={\cal C}\frac{|\Psi_{x}|^{2}}{|\Psi_{x}^{0}|^{2}}$ with ${\cal C}$ the normalization factor which can also be approximated using $\sum_{x}I_{x}/I^{0}=1$. The key difference, in comparison to the conventional method, is that, within each iteration step, the network parameters $w$ (and hence the wave functions) are updated multiple times. This enables us to use the $N_{\rm sample}$ sampled states much more efficiently. Furthermore, the update procedure can be efficiently parallelized and run on GPUs.

We plot the variational energies versus iteration step and wall time for a 60 sites SU(2) chain with a 16 layers CNN in Fig. 2(b) and (c) sup . As can be seen, the ISGO method converges with much fewer samples and much faster on GPU than the conventional method. We also implement the ISGO method using Tensorflow with auto differentiation sto , which allows us to try new network architectures much more easily. Our code for both RBM and CNN can be found in Ref. cod . We emphasize that, although we choose a particular gradient optimization algorithm Adam in our work, the concept of ISGO is general and can be implemented with any other optimization methods.

Numerical Results for 1D $SU(N)$ Spin Chain — We test our DNN on the Sutherland model, the 1D homogeneous SU($N$) spin chain governed by Hamiltonian (1). We pick this model for two main reasons. First, the ground state energy of this model can be exactly solved by Bethe Ansatz Sutherland (1975), which allows us to benchmark our results sup . Second, the number of spin states $N$ controls the complexity of the system, which allows us to systematically study the efficiency and accuracy of the DNN as the complexity of model grows. Numerical details can be found in sup .

Figure 3 shows our main results of the ground-state energies for various $N$ on an ${N}_{\text{site}}=60$ chain with a DNN with varying depth (i.e., number of layers $L$) and width (i.e., kernel size $K$). We tested two encoding methods for the input state. Fig. 3(a,b) correspond to the value encoding, while (c,d) correspond to the one-hot encoding. The value encoding imposes ordinality, i.e., different spin states are encoded into a number with the average value to be zero. For one-hot encoding, different spin states are encoded into a vector orthogonal to each other, and thus are not ordinal. For example, for an $N=3$ system, the three spin states are encoded into values of $-1/2$, 0, $1/2$ in value encoding, and into vectors $(1,0,0)$, $(0,1,0)$, $(0,0,1)$ in one-hot encoding. The one-hot encoding requires more computational resources (in terms of both memory and computational time) than the value encoding, and scales lineally with respect to $N$. However, in general it yields better accuracy than the value encoding method. This could be due to the fact that one-hot encoding encodes each spin state into a vector which effectively enlarges the dimension of the parameter space. Optimization in such an artificially enlarged space helps to prevent the system stuck in metastable states Pickard (2019).

Figures 3(a,c) display the results from a single-layer network with varying width. As one can see, increasing the kernel size $K$ helps bring the ground-state energy closer to the exact solution represented by horizontal dashed lines, which indicates that, for such a shallow network, large kernel size is necessary for capturing long-range effects mediated by nearest-neighbor interactions. Here, one-hot encoding performs significantly better than value encoding especially for SU($N>2$), where, no matter how large the kernel size is, the energies computed via value encoding do not converge to exact solutions. In Fig. 3 (b, d), we fix the width of the network, but vary its depth by adjusting the number of layers $L$. Even for a relatively small kernel size $K=3$, increasing $L$ helps to bring the computed ground-state energy closer to the exact result. Therefore, a DNN can capture long-range effect even with a small kernel size. For the SU(5) model (the largest $N$ we used in the calculation), the energy does not converge to exact solution using value encoding with kernel size 3. Simply by increasing the kernel size to 5, we can reduce the computed ground-state energy and match it with the exact solution as the black squares shown in Fig. 3(b). We vary the number of channels $F$ and find that the energy results are not sensitive to $F$. More details can be found in the Supplemental Materials sup .

The Bethe Ansatz method can yield energy spectrum and, in principle, the many-body wave function for exactly solvable models such as the one considered here. However, due to the complexity of a general many-body wave function, it remains a tremendous challenge to compute other useful quantities such as the correlation functions. Often, advanced numerical techniques are needed for such tasks Dufour et al. (2015); Nataf and Mila (2018). To further demonstrate the power of DNN, here we show our results for the loop correlation function

[TABLE]

where, the expectation value is taken with respect to the ground state, and $(m\cdots n)$ is the loop permutation operator that permutes the spatial indices in the wave function by $m\rightarrow m+1,m+1\rightarrow m+2,...,n-1\rightarrow n,n\rightarrow m$. Physically, this operator puts the spin in the original $n^{\rm th}$ position to the $m^{\rm th}$ position and correspondingly move the spins at the original $i^{\rm th}$ (with $i=m,...,n-1$) positions to their neighboring position on the right. The loop correlation function appears in the definition of the one-body density matrix of 1D strongly interacting quantum spinor gases whose ground state can be represented by a strong coupling ansatz wave function Levinsen et al. (2015); Yang et al. (2015); Yang and Pu (2016, 2017); Deuretzbacher et al. (2014) due to the fact that such wave functions must obey the permutation symmetry rule originated from quantum indistinguishability.

For the homogeneous system we considered here, $S_{m,n}=S_{r}$ with $r\equiv n-m$. We plot $S_{r}$ for an SU($N$) spin chain with $N_{\text{site}}=60$ spins in Fig. 4(a), and its discrete Fourier transform $S_{k}=|\sum_{r}S_{r}e^{ikr}|$ in Fig. 4(b). $S_{r}$ and $S_{k}$ characterize the correlation in the real and the momentum space, respectively. By taking the Jordan-Wigner transformation, an SU($N$) spin chain with each spin component having $N_{\text{site}}/N$ spins can be mapped to a nearest-neighbor interacting $N$-component fermionic system with each component having $N_{\text{site}}/N$ fermions Ogata and Shiba (1990); Yang and Pu (2017). The peaks of $S_{k}$ at $k=\pm\pi/N$, that can be clearly seen in Fig. 4(b), correspond to the Fermi points of those fermions. These peaks lead to the singularities of momentum distribution of strongly interacting spinor Fermi gases at the same momentum point Ogata and Shiba (1990); Yang and Pu (2017).

Conclusion and Outlook — We have constructed a DNN, combined with VMC, to study the ground state of the 1D SU($N$) spin chain. The key in our work is the development of the ISGO algorithm, which can be straightforwardly applied to any type of variational wave functions, for the optimization procedure. This algorithm allows us to efficiently train the network, and is particularly suitable for training DNNs which typically contains a large number of parameters. Note that the VMC with the ISGO algorithm may be interpreted as an RL process if we identify the Markov-chain state trajectories in the former as the state transitions/policies in the latter. We tested the network to solve the 1D SU($N$) spin chain model and systematically investigated the performance of the network by varying its depth and width. We have found that, when using value state encoding, as the complexity of the model increases by increasing $N$, it is not sufficient just to increase the width (i.e., kernel size) of the network, one needs to add more depth to capture the long-range correlation of the quantum state. We only show numerical results computed by the DNN up to 16 layers. We do not observe any significant benefit by using much deeper networks up to 100 layers on this model. This could be due to a potential problem of vanishing gradients in very deep networks (see Ref. Srivastava et al. (2015); He et al. (2016) and references therein), which may be alleviated via using other network architectures such as ResNet He et al. (2016), which we leave for future studies. Finally, we note that another key finding from our work, which has not been discussed in previous works, is the importance of input states encoding. We find that one-hot encoding, although requires more computational resource, in general leads to much more accurate results than value encoding.

In conclusion, our study clearly demonstrates that it is feasible to use DNNs to represent quantum many-body wave functions and to significantly enhance the efficiency of numerical quantum many-body computation. Applying machine learning techniques to quantum many-body physics is still a young and emerging field with many open questions. We believe that such investigation will not only benefit quantum physics, but may also help us to gain deeper insights into the neural networks.

I Bethe Ansatz Solution for Ground-State Energy of SU($N$) Spin Chain

The 1D homogeneous SU($N$) spin-chain model (Sutherland model) can be solved using Bethe-Ansatz as reported in Ref. Sutherland (1975). In this section, we provide details of how to numerically calculate the ground-state energy for a finite size SU($N$) spin-chain model with periodic boundary condition. Consider the case with all $N$ spin components equally populated with M spins. There are a total number of $N_{\text{site}}=MN$ spins. Following Eq. (63)($\alpha$)($\beta$)($\zeta$) in Ref. Sutherland (1975), the Bethe Ansatz equations consist of the following $(N-1)$ set of coupled equations:

[TABLE]

where $\theta(x)=-2\text{arctan}(x)$. $\alpha,\beta,\gamma,...$ are the $(N-1)$ sets of rapidities to be solved. The $j^{\rm th}$ set of equations (rapidities) has size $M_{j}=(N-j)M$, where $j=1,2,...,N-1$. For example, there are $(N-1)M$ different equations for $(N-1)M$ different $\alpha$’s, $(N-2)M$ different equations for $(N-2)M$ different $\beta$’s, and so on. So there are totally $(N-1)NM/2$ different equations. $J_{\alpha}$’s are integers or half odd integers and serve as quantum numbers Sutherland (1975). For $N=2$, there is only a single set of equations with size $M$:

[TABLE]

which are just the Bethe Ansatz Equations for 1D Heisenberg model.

For every set of equations with size $M_{j}$, the quantum numbers $J$’s are concentrated around the origin.

[TABLE]

And the ground state energy can be calculated by Eq. (56) in Ref. Sutherland (1975):

[TABLE]

Using the the relation of momentum and rapidity: $K_{l}=-\theta(2\alpha_{l})$, we have:

[TABLE]

Note that the rapidities in Eq. (LABEL:BAE) can be numerically solved for up to hundreds of spins using standard zero-finding libraries such as SciPy. Therefore, the exact ground state energy can be efficiently computed via Eq. (LABEL:BAE) and Eq. (12). We provide the python code for solving the ground state energy for arbitrary $N$ and $N_{\text{site}}$ in cod .

II The Unitary Transformation

In general, the ground-state wave function of the SU($N$) spin model Eq. (1) in the main text is not always positive. However, current DNNs only support positive wave functions. The reason is that at the final output, first $\text{ln}\Psi$ is calculated then an exponential is used to output $\Psi$, which is a positive value. In order to overcome this constraint, one needs to either add sign in the DNN by hand or do a unitary transformation for the Hamiltonian. We choose the latter by replacing the spin exchange operator $P$ in Eq. (1) of the main text by

[TABLE]

The new Hamiltonian has exactly the same eigenspectrum, but its ground-state wave function is positive. In this paper, we focus on studying the new model with the spin exchange operator $\tilde{P}$ under periodic boundary conditions. And we also need to change the loop correlation operators in Eq. (7) of the main text by adding corresponding signs.

III More Numerical Details and Results

We typically use 1500 iteration steps, the initial learning rate for Adam is set to $10^{-3}$ for the first 100 steps and then changed to $10^{-4}$, $N_{\text{sample}}=20000$ except for the first 10 steps to avoid out-of-memory, and $N_{\text{optimize}}=100$ for the first 1000 steps and then set to 10. We average the network parameters over the last 100 steps when calculating physical quantities.

Besides the kernel size $K$ and number of layers $L$, we have also studied the effects of changing the number of filters $F$. We fix the kernel size $K=11$ with one hidden layer, and the results are shown in Fig. 5. We found that for the value encoding in Fig. 5(a), the energies are not sensitive for SU(2) and SU(3) models. But for SU(4) and SU(5) models, the variational energies are still far away from the exact solutions. While for the one-hot encoding in Fig. 5(b), all energies are not sensitive to number of filters.

In Table I, we list the ground state energies obtained from our numerical calculation with a one-layer CNN and compare them with the Bethe Ansatz (BA) exact results. One can clearly see that the one-hot encoding yields much more accurate results than the value encoding. The one-hot results are comparable to those given in Ref. Dufour et al. (2015).

IV Comparison with Restricted Boltzmann Machines

So far, Restricted Boltzmann Machines (RBMs) remain as the most widely used machine learning tool in solving quantum many-body physics. In representing wave functions, an RBM can be regarded as a 1-hidden layer neural network Choo et al. (2018b), with activation function being $\ln(2\cosh(x))$. Also, the RBM wave function in Ref. Carleo and Troyer (2017) uses translational symmetry, which is equivalent to a CNN with kernel size $K=N_{\text{site}}$ and periodic boundary condition. The number of channels $\alpha$ in Ref. Carleo and Troyer (2017) plays the same role as the number of filters $F$ in our work. We plot the ground state variational energies of a convolutional RBM ($K\leq N_{\text{site}}$) with different $K$ in Fig. 6, which should be compared with Fig. 3(a)(c) in the main text, where we plot the energies as a function of convolutional kernal size $K$ in a one-layer CNN. It can be seen that the RBM results are very similar to those of the one-layer CNN although they use different activation functions. In addition, Fig. 6 also shows that, just as in our work, one-hot encoding with RBM has significant advantage over value encoding for $N>2$.

V Comparison of GPU and CPU performance

To compare the computational speed of the ISGO method on GPUs and CPUs, the wall time for the calculations of the SU(2) model on a $N_{\text{site}}=60$ size cluster with model architecture (L, F, K) = (16, 8, 3) is plotted in Fig. 2(c) in the main text. We have performed simulations on a single K80 GPU with 11GB memory of Google Colaboratory and two 2.30GHz Intel(R) Xeon(R) CPUs. The ISGO method has two steps: Monte Carlo sampling and importance sampling gradient optimization. The Monte Carlo sampling step is always calculated on CPUs. The simulation with $N_{\text{optimize}}=1$ corresponds to conventional GO method, whose performance on GPU (blue line) is similar to the one of the ISGO method on CPUs (green line), as in main text Fig. 2(c). Remarkably, the simulation with the ISGO method on GPUs (red line) arrives at a good variational energy at about 0.5h, while the conventional GO method (blue line) requires about 5h to obtain the similar accuracy for variational energy. The simulation with the conventional GO method on CPUs (magenta line) is the slowest one due to python overhead. The results shows that the simulations with the ISGO method on GPU are at least one order of magnitude faster than the conventional GO method on any hardware. We also observe the speedup is even greater for more complex networks.