Machine learning technique to find quantum many-body ground states of bosons on a lattice

We develop a variational method to obtain many-body ground states of the Bose-Hubbard model using feedforward artificial neural networks. A fully-connected network with a single hidden layer works better than a fully-connected network with multiple hidden layers, and a multi-layer convolutional network is more efficient than a fully-connected network. AdaGrad and Adam are optimization methods that work well. Moreover, we show that many-body ground states with different numbers of atoms can be generated by a single network.


Introduction
Recent developments in the field of artificial neural networks (ANNs), in combination with high-performance computers, have dramatically increased the ability of artificial intelligence. Techniques used in ANNs and machine learning have been applied to a wide variety of fields not only in engineering, but also in science, including physics research.
Pattern recognition is an important application of machine learning. By training an ANN with a large amount of data, e.g., many sample pictures, some features are extracted from the image data, and the ANN becomes able to classify the pictures. In physics, this ability of ANNs can be used to classify numerically or experimentally obtained data, which are sometimes complicated and cannot be identified by humans. Trained ANNs can discriminate different phases of numerically or analytically obtained many-body states. [1][2][3][4][5][6][7][8][9][10][11][12] A similar approach was developed for the analysis of experimental data 13) and telescope images. 14) In the above example of the classification of pictures, the memory size of the ANN is much smaller than the total size of the image data used for training. Nevertheless, the ANN acquires features of sample pictures and can even reproduce similar pictures. This implies that ANNs can efficiently encode and store the features of large amount of data, which is also applicable to physics problems. Many-body quantum states in a large Hilbert space can be efficiently stored in ANNs. [15][16][17][18][19][20][21][22] Thermal fluctuations can also be learned by ANNs; i.e., ANNs trained by Monte Carlo samples at finite temperature can reproduce thermodynamic properties. 23,24) Such trained ANNs can be used for efficient Monte Carlo updates. 25,26) Complicated functions of many variables can be stored in ANNs, which has been used for efficient simulation of molecular dynamics. 27) Quantum error corrections are also possible using ANNs. 28) Recently, a method to solve quantum many-body problems using ANNs was proposed. 15) It was demonstrated that the ground states and time evolutions of the quantum Ising and Heisenberg models can be obtained using ANNs. In this method, a quantum state is represented by a restricted Boltzmann machine, which consists of input and hidden units. When a spin configuration ↑↓ · · · is set to the input units, the corresponding wave function ψ(↑↓ · · · ) is obtained from the ANN. The internal parameters of the network are optimized in such a way that the wave functions produced by the ANN satisfy the desired properties, e.g., energy minimization.
Motivated by Ref. 15, a method to treat the Bose-Hubbard model was proposed in Ref. 29, where the feedforward network was used instead of the restricted Boltzmann machine. It was shown that the many-body ground state obtained by the method of ANN agrees very well with that obtained by exact diagonalization, even when the number of bases in the Hilbert space is much larger than the number of network parameters. This implies that the method proposed in Ref. 15 is applicable to a broad class of quantum many-body problems and that ANNs with machine learning are powerful tools to explore quantum many-body physics.
The present paper provides the extended results of the Letter in Ref. 29, in which a fully-connected network with a single hidden layer was only used with a simple steepestdescent method to optimize the network. In the present paper, we examine fully-connected networks and convolutional networks with multiple hidden layers. We show that a fullyconnected network having a single hidden layer with sufficient units yields better ground states than that with multiple hidden layers. On the other hand, a convolutional network with multiple hidden layers is more efficient than a convolutional network with a single hidden layer or a fully-connected network. With respect to methods for optimizing ANNs, Ada-Grad 30) and Adam 31) allow faster convergence than the simple steepest-descent method. We also show that many-body ground states with different numbers of atoms can be generated by a single ANN that has been multiply optimized for these numbers of atoms. Even when an ANN is optimized for a specific number of atoms, ground states with other numbers of atoms can be extrapolated approximately.
The remainder of the present paper is organized as follows.

Input
Hidden Output

Re logψ
Im logψ Fig. 1. Schematic diagram of the feedforward network. The input and output units are expressed as u (0) and u (L) , respectively, and u (1) , · · · , u (L−1) are hidden units. Given the values in the input units, the feedforward propagation generates the values in the output units, from which the wave function is cal- . Some or all of the units in the (n − 1)th layer are connected to those in the nth layer through the weights W (n) . Section 2 explains the method to obtain the quantum manybody ground state using an ANN. Section 3 provides the results of numerical calculations. Section 4 presents the conclusions of the present study.

Network architectures
We use the feedforward network illustrated in Fig. 1. The ANN consists of input, hidden, and output layers. The units in the input and output layers are denoted by u (0) and u (L) , respectively, and those in the hidden layers are denoted by u (1) , u (2) , · · · , u (L−1) . The number of units in the nth layer is written as N n , and the number of units in the output layer is fixed as N L = 2. Some or all of the units in u (n−1) are connected to those in u (n) through the weights W (n) .
In the present paper, we examine two types of feedforward neural networks. The first is a fully-connected neural network, in which each unit in the (n − 1)th layer is connected to all of the units in the nth layer. The network parameters are the weights W (n) and biases b (n) , which are real N n ×N n−1 matrices and N n -component vectors, respectively. First, we set values to the input units u (0) i , which are transferred to the next layer as In the hidden layers, an activation function f is applied as where f must be a nonlinear function, and here we adopt In the final layer, the bias is absent, b (L) = 0, because the wave function is only multiplied by an overall factor (see Eq. (12)). The total number of network parameters for the fully-connected network is then The second one is the convolutional neural network, in which the weights W (n) act as local filters and the units u (n) consist of multiple channels. The values set in the input units u (0) are transferred to the next layer as where the index k denotes the channel, and F 1 is the size of the filter W (1) pk for each channel. In Eq. (5), u (1) j,k is generated only from u (0) j , u (0) j+1 , · · · , and u (0) j+F 1 −1 , and thus the local feature in the input units is captured by each filter and transferred into each channel. The number of units in u (1) is N 1 = N 0 C 1 , where C n is the number of channels in the nth layer. The subsequent convolutional layers have the form, where all of the C n−1 channels in u (n−1) are filtered and summed to generate each output channel. The numbers of units in u n−1 and u n are thus N 0 C n−1 and N 0 C n , respectively. Finally, the units u (L−1) produced by the convolution layers are fully-connected to the output units as where j = 1, 2. The total number of network parameters is thus In the terminology of ANN, this network consists of multiple convolution layers with a unit stride and no pooling layers, followed by a fully-connected layer.

Quantum many-body states
A quantum many-body state of bosons on a lattice is expressed by the feedforward ANN as follows. An arbitrary state can be expanded by the Fock states as where n = (n 1 , n 2 , · · · , n M ) represents a particle distribution on the lattice sites, with M being the number of sites. For the total number of particles N, the number of Fock state bases, i.e., the number of n satisfying M i=1 n i = N is The number of input units is taken to be N 0 = M. We set the input units as and the feedforward propagation in Eqs. (1) and (2) [or Eqs. (5) and (6)] is then performed. The wave function ψ(n) is calculated from the output units u (L) as Although the ground-state wave function of bosons can be taken to be real and positive, we include the phase u (L) 2 in Eq. (12) for future use. By calculating ψ(n) for all possible n, we can construct the many-body state in Eq. (9). The information of the many-body quantum state is therefore stored in the network parameters W (n) and b (n) . Our aim is to optimize the network parameters so that the corresponding many-body quantum state is as close to the ground state as possible.
Expectation values of quantities are calculated by the Monte Carlo method with Metropolis sampling. If we adopt the trial n 1 → n 2 with the probability min[1, |ψ(n 2 )/ψ(n 1 )| 2 ], where |ψ(n 2 )/ψ(n 1 )| 2 can be obtained from the network by the above procedure, the sampling probability distribution of n becomes |ψ(n)| 2 / n ′ |ψ(n ′ )| 2 ≡ P(n). Using this sampling of n, we can approximate as if the number of samples N s is sufficient. Using Eq. (13), the expectation value of a quantityÂ is calculated as

Network optimization
The Hamiltonian for the system is given bŷ where J > 0 is the hopping coefficient, and U is the on-site interaction energy. The operatorâ i annihilates a particle in the ith site, andn i =â † iâ i is the number operator, where the Bose commutation relation [â i ,â † j ] = δ i j is satisfied. In order to optimize the network parameters W (n) and b (n) , we need to calculate the derivative of the expectation value of the Hamiltonian Ĥ with respect to these network param-eters. Since the wave function ψ(n) depends on the network parameters, we have where w is one of the network parameters and The derivative ∂ψ(n)/∂w can be calculated systematically using the method of back propagation. 32) Using the stochastic approximation in Eqs. (13) and (14), Eq. (16) is obtained as There are various ways to update the network parameters to reduce the expectation value of the Hamiltonian. In the steepest-descent method, the ith network parameter w i is updated as where α < 1 controls the magnitude of change in each update.
More efficient methods to update the network parameters have been developed in the field of machine learning. In the Ada-Grad method, the network parameters are updated as 30) v where ǫ ≪ 1 avoids division by zero, and γ < 1. The Adam method is given by 31) for the ℓth update, where β 1 = 0.9 and β 2 = 0.999 are usually used, and δ < 1. The initial values of v i and m i in Eqs. (20) and (21) are zero, and the parameters α, γ, and δ are chosen so that the optimization works efficiently. It is known that if a valley exists in the energy landscape, the steepest-descent method may cause oscillation between the two sides of the valley, while AdaGrad and Adam can avoid oscillation and  efficiently decrease the energy. The stochastic reconfiguration method 33) used in Ref. 15 may also be efficient, but it has a higher computational cost for each update. For the steepest descent, AdaGrad, and Adam, the cost for each update is O(N w ), 34) where N w is the number of network parameters, whereas the cost for the stochastic reconfiguration method is O(N 3 w ). The procedure for obtaining the approximate many-body ground state using the ANN is as follows. First, we initialize the network parameters W (n) and b (n) with random numbers obeying normal distributions. The standard deviations of the normal distributions are taken to be 1/ √ N n−1 for W (n) and b (n) . 35

Numerical Results
We consider a one-dimensional (1D) system of N bosons on M sites with the periodic boundary condition. The ground state of the Hamiltonian in Eq. (15) We first compare three optimization methods: steepest descent, AdaGrad, and Adam. Figure 2 shows the energy H M as a function of the number of updates using the three methods. We can see that the energy decreases quickly for the initial ∼ 1000 updates and then gradually converges to the final value. The convergence is faster for AdaGrad and Adam  than for the steepest-descent method. The update rates used in Fig. 2 are α = 0.3, γ = 0.3, and δ = 0.005, which are chosen so that the energy convergence becomes the most efficient. When these parameters are smaller, the convergence becomes slower, and when the parameters are too large, the calculation becomes unstable. In the following calculations, we adopt Adam with δ = 0.005 for the network updates. Figure 3 shows the U-dependence of the energy convergence with respect to the network updates. For the noninteracting case, the deviation from the exact energy rapidly converges to (E − E exact )/J 0.002, whereas for a large on-site interaction U, the convergence is slow and (E − E exact )/J ∼ 0.02 even after 20000 updates. Slow convergence occurs for U/J 5, which is the Mott insulator regime. Increasing the number of updates to 60000, the energy difference becomes (E − E exact )/J ≃ 0.01 for U/J = 10.
In order to improve the slow convergence for large U, we change the size and depth of the fully-connected network. Figure 4(a) shows the dependence of the energy convergence on the number of hidden units N 1 for a single hidden layer. The energy convergence is improved as N 1 is increased, and saturates at N 1 100. Figure 4(b) shows the case of a fullyconnected network with two hidden layers. Although the convergence for N 1 = N 2 = 20 in Fig. 4(b) is better than that of the single hidden layer with N 1 = 20 in Fig. 4(a), the convergence for N 1 = N 2 = 40 becomes worse compared with that of the single hidden layer with N 1 = 40. This implies that the increase in the depth of the fully-connected network may improve the energy convergence due to the increase in the capability of the network, but it can also be counterproductive due to the increase in the complexity of the network.
In Fig. 4, the value of U is linearly increased from 0 to 10J in the initial 1000 updates. Such a gradual ramp of U prevents the network from becoming trapped in a local minimum of the energy. Furthermore, the initial ramp of U accelerates the energy convergence (compare the N 1 = 40 line in Fig. 4(a) with the U/J = 10 line in Fig. 3). In the following calculations, we linearly ramp the value of U in the initial 1000 updates. We next consider the convolutional neural network in Eqs. (5)- (7). The input layer is connected to the first convolution layer through filters of size F 1 = F, which produces C 1 = C channels; i.e., N 0 = M input units are connected to N 1 = CM hidden units. When L = 2 (a single convolution layer), N 1 = CM units produced by the convolution layer are fully connected to the output units. When L > 2 (L−1 consecutive convolution layers), CM output units of the convolution layer are input into the next convolution layer, which also produces CM output units, and finally the (L − 1)th layer is fully connected to the output layer. Com-  paring (L, C, F) = (2,4,8), (3,4,8), and (4,4,8) in Fig. 5, we find that a network with two convolution layers (L = 3) exhibits much better convergence than a network with a single convolution layer (L = 2). However, a further increase in the convolution layer (L = 4) results in no improvement, and therefore the improvement saturates at L = 3. The decrease in the number of channels (C = 1) or the size of the filters (F = 4) makes the convergence worse, while no improvement is obtained for C > 4 and F > 8. Thus, in the present case, the convolutional ANN with (L, C, F) = (3, 4, 8) yields the fastest convergence and the smallest energy, which is much better than the fully-connected ANNs. We evaluate the accuracy of the ground-state wave function obtained by the optimized network. The fidelity of the wave function is defined as where ψ is the wave function generated by the ANN, and ψ exact is that obtained by the exact diagonalization of the Hamiltonian. The sum in Eq. (22) is taken for all possible n. When the wave function ψ is the exact ground state, the fidelity in Eq. (22) is f = 1. Figure 6 shows the error 1 − f as a function of U/J. In order to obtain each value in Fig. 6, 10000 updates of the network parameters are performed, and the results of five runs with different initial network parameters produced by random numbers are averaged, where error bars represent standard deviation. We find that the errors 1 − f are smaller for the convolutional neural network than for the fully-connected network, which is consistent with the results On the other hand, the number of network parameters for the fully-connected network used in Fig. 6 is N FC = 680 for N 1 = 40 (Eq. (4)), and that for the convolutional network is N conv = 280 for F = 8 and C = 4 (Eq. (8)), which are much smaller than N Fock . Therefore, the information of the many-body wave function is compressed and stored in the ANN very efficiently. It is remarkable that such compressibility of the wave function is automatically achieved in the optimization process of the network.
To see how the quantum many-body states are stored in the ANNs, the network weights W (n) after 10000 updates are visualized in Fig. 7. The fluctuations in the weights are larger for U/J = 10 than for U/J = 0. For the fully-connected network in Fig. 7(a), no regularity or meaningful pattern appears in the weights. Although the many-body wave function produced by the network has translational symmetry, the weights have no apparent translational symmetry with respect to the site index i for the fully-connected network in Fig. 7(a). By contrast, for the convolutional network in Fig. 7(b), the weights W (3) in the final fully-connected layer are independent of the site index i. This indicates that the network has translational symmetry, since the convolution layers have translational symmetry by the definition (i.e., the filters are shared by all of the site indices j in Eqs. (5) and (6)).
The translational symmetry of the convolutional network may be related to the faster convergence and smaller energy achieved in the optimization process, as compared with the fully-connected ANN. Let us consider the capability of both ANNs. For example, the convolutional ANN with (L, C, F) = (3,4,8) in Fig. 5 is a subset of the fully-connected ANN having two hidden layers with N 1 = N 2 = MC = 56, namely, the former is realized by the latter with constraints on W (n) and b (n) . Therefore, potentially, the latter has the ability to represent the many-body wave function more accurately than the former. Nevertheless, the fully-connected ANN is worse than the convolutional ANN, because the large degree of freedom of the fully-connected ANN makes its optimization inefficient. The convolutional ANN, on the other hand, takes into account the translational symmetry of the system, and a smaller number of network parameters makes the optimization efficient. The filters may also be suitable for capturing the local correlations produced by the on-site local interaction. Thus, presumably, the convolutional ANN is quite compatible with physical systems with local interaction and translational symmetry. We next consider larger systems for which exact diagonalization is difficult or almost impossible. Instead of exact diagonalization, we adopt the method of the density matrix renormalization group (DMRG), 36,37) which is known to be very accurate for 1D many-body problems. Figure 8 shows the difference between the energy per particle obtained by the present ANN method, E/N, and that by the DMRG, E DMRG /N, as a function of the number of sites M. For the fully-connected neural network, the error ∆ε appears to increase exponentially with M as ∆ε ∼ exp(κM) with κ ≃ 0.12. When the number of hidden units is increased from N 1 = 40 to 80, the accuracy is improved by an order of magnitude. The convolutional neural network yields energies with much better precision. The error ∆ε also appears to increase with M as ∆ε ∼ exp(κM), where κ ≃ 0.05 is smaller than that of the fully-connected neural network. The results for larger particle density N/M = 1.5 are similar to those for N/M = 1.
We have thus far considered the case in which a single many-body ground state is stored in the ANN. We next try to obtain the ground states for different numbers of atoms N by a single optimized ANN. In order to optimize the network in such a manner, in each update step described in Sec. 2.3, we choose N randomly, and the Metropolis samplings are per- formed for n with N particles. In Fig. 9, we randomly choose N from 5 to 15 in each update step, and a total of 10000 updates are performed. The fidelity is then calculated for each N using the optimized network. Figure 9 shows that 1 − f is always smaller than 0.01 for U/J = 1, which indicates that the multiple many-body ground states are stored in the single ANN. However, the precision of each state is worse than in the case in which the ANN is optimized for a specific N (see Fig. 6). For U/J = 10, the error is prominent at N = 10, which is the Mott insulator state with unit filling. The convolutional network is better than the fully-connected network also in this case. Figure 10 shows similar results, but the network is optimized only for a specific N (= 7, 10, or 13). For example, when the network is optimized for N = 7, the fidelity is best at N = 7, as expected. Note that the fidelity is good not only for N = 7, but also for N = 6 and 8. For the convolutional network optimized for N = 10, the fidelity is below 0.01 in the range of N shown in Fig. 10(b). These results imply that the present method may also be used for extrapolating (or interpolating) the quantum many-body states, i.e., the ANN optimized for certain parameters may generate approximate many-body states for other parameters.

Conclusions
In conclusion, we have developed a method to obtain the quantum many-body ground state of the Bose-Hubbard model using a feedforward artificial neural network. Although the simple steepest-descent method was only employed in the previous Letter, 29) we examined AdaGrad and Adam in the present paper and found that the convergence is better (Fig 2).  The accuracy of the present method becomes worse as the on-site interaction U is increased (Fig. 3). In order to increase the accuracy, we investigated a deep (multi-layer) fully-connected network. However, a single hidden layer with a sufficient number of hidden units is found to be better than multiple hidden layers (Fig. 4). We then investigated the deep convolutional network and found it to be much more efficient than the fully-connected network (Figs. 5 and 6). We found that the convolutional network has translational symmetry (Fig. 7). The convolutional network is also promising for studying large systems (Fig. 8). Multiple quantum manybody states can be stored in a single ANN (Figs. 9 and 10). At present, it is unclear whether the present method of obtaining quantum many-body states can surpass other existing methods in terms of precision and computational resources. For 1D cases, the DMRG seems better, whereas for 2D and 3D cases, the present method may have advantages. At the very least, the ANN is a very versatile scheme for representing quantum many-body states and may provide an initial choice to tackle quantum many-body problems for which effective solution methods are unknown. In order to confirm this possibility, we must confirm that the method works for various other quantum many-body problems.