{'fraction_non_alphanumeric': 0.045653021442495126, 'fraction_numerical': 0.012787524366471735, 'mean_word_length': 4.421898118790954, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Restricted-Boltzmann-Machine Learning for Solving Strongly Correlated Quantum Systems

We develop a machine learning method to construct accurate ground-state wave functions of strongly interacting and entangled quantum spin as well as fermionic models on lattices. A restricted Boltzmann machine algorithm in the form of an artificial neural network is combined with a conventional variational Monte Carlo method with pair product (geminal) wave functions and quantum number projections. The combination allows an application of the machine learning scheme to interacting fermionic systems. The combined method substantially improves the accuracy beyond that ever achieved by each method separately, in the Heisenberg as well as Hubbard models on square lattices, thus proving its power as a highly accurate quantum many-body solver.


I. INTRODUCTION
Obtaining accurate ground-state wave functions of many-body quantum Hamiltonians is one of the grand challenges in condensed matter physics. Great successes so far are, for example, Bardeen-Cooper-Schrieffer (BCS) wave functions for conventional superconductivity [1], Bethe-ansatz wave function for one-dimensional interacting systems [2], and Laughlin wave functions for fractional quantum Hall effect [3].
However, in order to construct the ground-state wave functions of many-body interacting systems and grasp the essential physics encoded in them, we often need to resort to numerical estimates. Currently, many numerical techniques are available such as the variational Monte Carlo (VMC) method [4][5][6][7][8][9], the density matrix renormalization group [10,11], tensor network methods [12,13], and the path-integral renormalization group [14]. Among them, the VMC method offers an accurate ground-state wave function for quantum spins as well as fermions on various lattices.
Recently, alternative approaches, based on machine learning, have attracted growing attention in many-body physics . In particular, Carleo and Troyer [49] have proposed a machine-learning algorithm, which uses a restricted Boltzmann machine (RBM) as a variational wave function |Ψ for representing the ground states of quantum spin systems. In this scheme, hidden artificial neurons are introduced on top of the physical degrees of freedom (quantum spins), to mediate entanglement in the state. The RBM variational wave functions are selfoptimized through machine learning.
We can express a general quantum state |Ψ by using the Fock space basis {|x } in the form of a variational function as with a correlation factor F (x) [50,51] and a reference state |φ ref . The RBM wave function in Ref. [49] is obtained by employing an RBM for F (x) and a product state for |φ ref so that x|φ ref = 1 is satisfied for any orthonormalized complete set x. The product state |φ ref is not able to describe nonlocal quantum entanglement, although it is essential in strongly correlated systems.
Then the entanglement has to be represented solely by the RBM factor F (x). However, alternative choices of |φ ref may already incorporate typical quantum correlations and can potentially allow |Ψ to more efficiently capture ground state entanglement. In fact, in the manyvariable VMC (mVMC) method [9], a pair-product (PP) wave function (or equivalently geminal wave function in quantum chemistry [52][53][54]) is chosen as |φ ref , which can efficiently capture a substantial part of the non-local entanglement in strongly correlated quantum systems by using many variational parameters.
In this paper, we propose a variational wave function for studying strongly correlated quantum systems called RBM+PP which combines flexible and nonempirical correlation factor given by RBM and entangled reference state given by the PP wave function to inherit the advantages of both. When applied to the two-dimensional (2D) Heisenberg model on a square lattice, we show that our method significantly outperforms the original RBM method [49], which itself outperforms existing numerical techniques for finite lattices based on tensor networks.
The PP wave function can also flexibly incorporate non-local correlations in fermionic systems and account for the fermionic sign, allowing RBM+PP to be applied to interacting systems of fermions. When applied to the Hubbard model we show that the combined method achieves greater accuracy than either method applied separately. To the best of our knowledge, this is the first application of the RBM-based wave functions to interacting fermions. The RBM+PP method thus provides a powerful tool not only for quantum spins but also for highly entangled quantum states such as strongly correlated itinerant fermions.
The structure of the paper is as follows. In Sec. II, we explain the RBM+PP method after introducing individual RBM and PP wave functions. In Sec. IV, we apply the RBM+PP scheme to Heisenberg and Hubbard models and show significant improvement from RBM and mVMC results. The representability of the RBM+PP wave function is discussed in detail in Sec. V. Finally, we give a summary and present future perspectives in Sec. VI.

II. RBM+PP METHOD
In this section, we will define RBM and PP states and explain how these are combined in RBM+PP method.

A. RBM wave function
The RBM state in Ref. [49] for spin Hamiltonians is given by setting F (x) = N (x) with a neural-network correlation factor N (x) and |φ ref to be the product state |φ product ( x|φ product = 1) in Eq.(1), which from now we refer to as product-basis RBM (P-RBM) ( Fig. 1(a)). N (x) is defined by an artificial neural network (ANN) as where x = (σ 1 , σ 2 , . . . , σ N visible ) is a real space configuration of N visible physical variables and σ i is the i th discrete-valued physical variable (visible-layer spin variable). In the S = 1 2 Heisenberg model, we take σ i = 2S z i = ±1 with S z i being the z-component of the S = 1 2 spin at site i. Here N visible is equal to the number of sites N site . The auxiliary pseudo spin variables h k = ±1 are for the hidden neurons and {a i , W ik , b k } is a set of variational parameters. In this study, we take variational parameters to be real. Importantly, as there are no weights connecting hidden neurons, the sum over hidden variables can be evaluated exactly and Eq.(2) can be reduced to the form N (x) ≡ k 2 cosh b k + i W ik σ i ×e i aiσi , which can be computed efficiently for each x.
In fermionic models we define a different RBM state which we refer to as F-RBM. The correlation factor F (x) = N (x) is taken similarly to, but by slightly modifying the spin case. For |φ ref , our most primitive choice is a Fermi-sea state rather than the product state. This is because the product state is too poor in representing fermionic entanglement. The Fermi sea is the ground state of the noninteracting fermion lattice models, much like product states are ground states of noninteracting spin models, and is thus able to account for the most primitive part of the fermionic entanglement and signs. We can use the same form of N (x) as Eq. (2) by doubling the number of the visible-layer variables (N visible = 2N site ) and mapping fermionic modes to these spins as (σ 2i , σ 2i−1 ) = (2n i↑ −1, 2n i↓ −1) where n iσ is the number operator for the fermions at site i with spin σ.
As discussed in Ref. [49], the accuracy of the wave function can be controlled by the "hidden variable density" α, which is defined as N hidden /N visible with the number of neurons N hidden in the hidden layer and the number of physical variables N visible in the visible layer.

B. PP wave function
Whereas the RBM states in the first step use the product and Fermi-sea states for |φ ref , more sophisticated choices of |φ ref are able to incorporate more involved entanglement directly into the state. In this work we will use the pair-product (PP) state |φ pair as |φ ref . The PP wave function is given by where N e is the number of electrons, f σσ ′ ij are variational parameters, and c † iσ is the operator creating a σ-spin electron at site i. For a given real space configuration x, φ pair (x) = x|φ pair can be expressed as the Pfaffian of a matrix, and can be calculated efficiently, much like N (x). Note that an accurate description of the node position is crucially important for fermionic wave function while the simple product state with the positive definite coefficients does not describe the node. In contrast, the Pfaffian wave function is able to optimize the nodal structure within the framework of the Pfaffian wave function. Therefore, the PP wave function can account for typical non-local entanglement not only in non-frustrated spin systems but also in frustrated spin and fermionic systems.
In spin models, to prohibit the double occupation, PP wave function is supplemented by the Gutzwiller factor P ∞ G = i (1 − n i↑ n i↓ ). This form of reference function |φ ref = P ∞ G |φ pair is able to represent resonating valence bond (RVB) wave functions [55,56].

C. RBM+PP wave function
In this paper, we will study the combined wave function RBM+PP ( Fig. 1 (1). We remark that the RBM+PP wave function is similar to that used in the mVMC method [9], except that the mVMC method uses an empirical form of F (x) instead of more flexible and unbiased neuralnetwork factor N (x). Specifically, in the mVMC method, F (x) is given by F (x) = x|P G P J |x with Gutzwiller P G (controlling double occupancy) [51] and Jastrow P J (for long-ranged charge-charge correlation) [50] factors. In Appendix A, we show that the neural-network factor N (x) is indeed more flexible than the empirical factors by showing that N (x) can represent both the Gutzwiller and Jastrow factors. The N (x) factor can also represent many-body (more than two-body) correlations [57] at the same time.
Various symmetries can be imposed on the wave function to improve accuracy and reduce computational cost [58].
In this study, as in Ref. [49], we impose translational symmetry in the variational parameters {a i , W ik , b k } in the RBM. In the antiferromagnetic Heisenberg and half-filled Hubbard models, because i σ i = 0 holds in the ground state, the translationallyinvariant bias term a i = a becomes irrelevant. Therefore, we neglect it. Furthermore, as the ground states of the Hamiltonians considered have total spin S = 0 and momentum K = 0, we apply the projections onto these subspaces, respectively, L S=0 and L K=0 , to |φ pair or P ∞ G |φ pair to improve accuracy. If we apply both L S=0 and L K=0 , the reference state becomes

D. Machine learning of variational parameters
The form of the wavefunction in Eq. (1) allows calculating physical quantities, and derivatives with respect to variational parameters to be approximated efficiently using Markov chain Monte Carlo sampling over the probability distribution p(x) = Ψ|x x|Ψ / Ψ|Ψ . We use a machine learning method (called stochastic reconfiguration in Ref. [8] and natural gradient in Refs. [59,60]) to optimize the variational parameters in the wave function with respect to the energy. The computational cost of the optimization scales as O(N 3 site ) for RBM+PP, compared to O(αN 2 site ) for P-RBM [49]. Thus, the improved accuracy of RBM+PP over P-RBM comes at some additional computational cost. Details and comparisons of the variational wave functions we introduced in this section are listed in Appendix B.

III. MODELS
We apply the RBM+PP scheme to calculate the ground states of 2D S = 1 2 antiferromagnetic (AFM) Heisenberg and 2D Hubbard models on the square lattice. Their Hamiltonians are defined as follows: The sum over sites i, j is restricted to nearest-neighbor pairs. We take the exchange constant J and hopping t as an energy unit in each case. The onsite repulsion U controls the strength of correlation in the Hubbard model. In our calculations, fully periodic boundary conditions are imposed for the Heisenberg model, while periodic (in x-direction) and anti-periodic (in y-direction) (P-AP) boundary conditions are imposed for the Hubbard model. Further details of the computation conditions are available in Appendix C.

IV. RESULTS
A. Heisenberg model Figure 2 shows the RBM+PP ground-state energy for 8 × 8 2D Heisenberg model compared to quantum Monte Carlo calculations using the stochastic series expansion (SSE-QMC) at sufficiently low temperature T of 1/T = 64 [61], which gives practically exact ground state energy. For comparison, the mVMC results [9] and the P-RBM wave function employed in Ref. [49] are shown. Here, mVMC is equivalent to α = 0 of RBM+PP. This is because the Gutzwiller factor P G is fixed to freeze charge degrees of freedom and Jastrow factor P J becomes irrelevant in the absence of charge degrees of freedom.
As discussed in Refs. [64,65], the resonating valence bond (RVB) wave function is known to provide a highly accurate description of the 2D Heisenberg model. The relative error in the result of the mVMC function, which can represent the RVB wave function, is indeed less than 0.2 percent. We see that the non-empirical P-RBM wave function is also powerful, giving a comparable accuracy to the mVMC results. The RBM+PP wave functions, which take advantages of the above two, substantially improve the accuracy of independent mVMC and P-RBM schemes.
It is interesting to note that each of the three curves in Fig. 2 corresponds to the results with the very same form of correlation factor N (x) but with different reference wave functions |φ ref . The lowest energy is obtained G |φ pair , with the product state |φ ref = |φ product having the highest energy. Therefore, improving reference function helps the RBM to learn the ground state more efficiently.
In Fig. 2(b), we plot the total energy as a function of its variance ∆ var = H 2 − H 2 / H 2 . The variance is zero in the case of an exact ground state (or more generally, an exact eigenstate of Hamiltonian). By the linear fit of the energy as a function of the variance and extrapolating to ∆ var = 0, we can obtain a more accurate estimate of the ground state energy [8,14,[66][67][68]. The variance extrapolation works better for the RBM+PP than the P-RBM wave function because the variance is TABLE I. P-RBM and RBM+PP results for spin structure factor S(π, π) in 2D Heisenberg model. For comparison, S(π, π)/Nsite = 0.05986(3) in mVMC and 0.059280(3) in SSE-QMC [61] results.
RBM+PP can also be used to accurately calculate other physical quantities, besides the energy. For instance, we also measure spin structure factor S(q) = 1 3Nsite i,j α=x,y,z S α i S α j e iq·(ri−rj ) . The result for S(q peak ) = S(π, π) is shown in Table I In both cases (U/t = 4 and U/t = 8), the F-RBM has an error of several percent. The RBM+PP method, in which the variational parameters in both the N and |φ pair are optimized, achieves significantly higher accuracy. We see the same trend as the Heisenberg model; namely, the accuracy is improved by choosing a better reference wave function |φ ref .
The RBM+PP wave function also surpasses the accuracy of the mVMC wave function, which indicates superiority of more unbiased neural-network factor N (x) to the empirical Gutzwiller-Jastrow factors P(x) = x|P G P J |x . We expect that the advantage of a selfoptimized neural network will be more substantial for more complicated Hamiltonians than the single-band Hubbard model. In more complex Hamiltonians, more flexible forms for correlation factors are likely to be necessary. Another advantage of the RBM+PP to the mVMC methods is that the accuracy improves as U increases (see Appendix D and Fig. 4), while the mVMC results show the opposite trend [9].
At U/t = 4 ( Fig. 3(a)), the best RBM+PP accuracy is comparable to the accuracy of TNVMC (tensor net- work combined with mVMC), where the relative error is ∼ 0.25 percent as obtained at available maximum tensor bond dimension D = 16, which is the computationally practical upper limit [70]. In the TNVMC method, the mVMC wave function is supplemented by the tensor network enabling one of the most accurate schemes among existing numerical methods [70]. Compared to the TNVMC method, the RBM+PP has an advantage, because it can be applied flexibly and easily to any kind of lattice and does not require an involved contraction procedure [13] in contrast to the TNVMC method.
The spin structure factor S(π, π) is listed in Table II. At both U/t = 4 and U/t = 8, with increasing α, the value becomes closer to the exact AF-QMC value.

V. DISCUSSION
While the physical properties of the RBM have only recently started being discussed in condensed matter physics [72][73][74][75][76][77][78][79][80], more general discussion of representability can be traced back to earlier studies [81][82][83], which show that the RBM is able to describe any smooth function, if arbitrarily large α is allowed. In the present 2D Heisenberg model, a gauge transformation can make probability amplitude of the exact ground-state wave function |Ψ GS positive ( x|Ψ GS > 0 for all x). Therefore, the exact ground state can be represented by realvariable RBM with infinite α. Accordingly, the relative error should go to zero as 1/α → 0. Indeed, in Fig. 2(a), it looks that the RBM+PP energy curves (red and green) start bending toward 0 as 1/α decreases. It is likely that a better reference state makes x|Ψ GS / x|φ ref (to be represented by N (x)) smoother and helps to reach faster convergence at small α.
For fermionic problems, nodal structure of wave functions is crucial [84], which is beyond the representability of the real-variable RBM giving positive N (x) [79]. In the RBM+PP, |φ pair is expected to accurately reproduce the nodal structure. Then, x|Ψ GS / x|φ pair may become smooth enough so that with moderate α, N (x) can represent a quick convergence to the exact value at α → ∞. However, rigorously speaking, the nodal structure of |φ pair is likely to be different from the exact one even when the variational parameters contained in |φ pair are ideally optimized. Therefore, introduction of complex variational parameters in the RBM part may be useful to adjust the nodes to the exact positions beyond the framework of the Pfaffian wave function. Although it is an interesting open question whether the energy curve as a function of 1/α in the Hubbard model (Fig. 3) goes to 0 as 1/α → 0, in practical computations, the optimization at larger α becomes more and more difficult, which might hamper the expected convergence.

VI. SUMMARY AND PERSPECTIVES
In this work, we have proposed a new variational ansatz for studying the ground states of many-body interacting quantum systems. Our variational wave function, which we call RBM+PP, combines the RBM based on the ANN and the mVMC methods. We have shown that, in both 2D Heisenberg and Hubbard models, the RBM+PP results show a dramatic improvement of accuracy over simple neural network wave functions (the P-RBM and F-RBM wave functions). We also see the superiority of the RBM+PP to the mVMC method. Since the RBM+PP method can be flexibly applied not only to bosonic (or spin) systems but also to fermionic problems, the RBM+PP method offers a wide range of applications with high accuracy and a reasonable computational cost.
As a future perspective, it would be interesting to go beyond the RBM structure and introduce second hidden layer (deep Boltzmann machine (DBM)). DBM is argued to have more efficient representation of certain manybody wave-functions than RBM [75]. In DBM, the spin variables of the neurons distributed in more than one hidden layer cannot be traced out analytically, thus we need to introduce additional Monte Carlo samplings for hidden spins.

ACKNOWLEDGMENTS
We acknowledge useful discussions with Giuseppe Carleo, Takahiro Ohgoe, Hui-Hai Zhao, and Kota Ido. We also thank Shiwei Zhang for providing us with the auxiliary-field quantum Monte Carlo result of energy and spin structure factor in the 2D Hubbard model. The implementation of the RBM+PP scheme is done based on the mVMC package [85]. The computation was mainly done at Supercomputer Center, Institute for Solid State Physics, University of Tokyo. This work was financially supported by Grant-in-Aids for Scientific Research (JSPS KAKENHI) (No. 16H06345 and No. 17K14336) from Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan. This work was also supported in part by MEXT as a social and scientific priority issue (Creation of new functional devices and highperformance materials to support next-generation industries (CDMSI)) to be tackled by using post-K computer. We also thank the support by the RIKEN Advanced Institute for Computational Science through the HPCI System Research project (hp160201,hp170263) supported by MEXT. Y. Y. was also supported by PRESTO, JST (JP-MJPR15NF).
Appendix A: Ability of RBM to represent Gutzwiller and Jastrow factors The Gutzwiller factor P G = exp(−gn i↑ n i↓ ) at site i can be recast as (except for trivial constant factor and one-body potential) where (σ 2i , σ 2i−1 ) = (2n i↑ −1, 2n i↓ −1) are physical variables in the visible layer defined in Sec. II A. This interaction between physical variables can be mediated by adding one hidden neuron h as (except for trivial constant factor) with W 1 = 1 2 arcosh (exp (|g|/2)) and W 2 = −sgn g × W 1 . This form is consistent with the neural-network factor N defined in Eq. (2). In the very same way, we can show that the neural-network correlation factor N can represent the Jastrow factor P J = exp(− 1 In the present study, the parameters are optimized by the stochastic reconfiguration (SR) method [8]. The same optimization scheme is called natural gradient in Refs. [59,60]. This optimization is equivalent to the imaginary-time evolution e −τ H |Ψ of the wave function |Ψ to reach the ground state for sufficiently large imaginary time τ in the truncated Hilbert space which is spanned by the variational wave function. In the SR optimization, the variational parameters γ m (m = 1, 2, . . . , N v ) at the p-th iteration are updated as where the difference in the update ∆γ with a small imaginary time step ∆τ . Here, S is a positive definite matrix given by with a normalized variational wave function |Ψ = |Ψ / Ψ|Ψ and |∂ γmΨ = ∂ ∂γm |Ψ . The g vector is the gradient of energy with respect to the γ parameters: To stabilize the optimization, we introduce the stabilization factor to the diagonal elements of the S matrix as Here, ǫ 1 scales the diagonal elements and ǫ 2 gives a constant shift to the S matrix. Although ǫ 2 strongly stabilize the optimization, it sometimes makes convergence slower, so we typically put a small number to ǫ 2 : ǫ (p=0) 2 ∼ 10 −3 and we gradually decrease to ǫ ∞ 2 ∼ 10 −7 -10 −6 through the first several hundred iterations of optimization. As for ǫ 1 , we found that smaller ǫ 1 (∼ 10 −5 -10 −4 ) factor for neural-network-related variational parameters {b k , W ik } sometimes helps to lower the energy. Meanwhile, we also found that a small ǫ 1 ( 10 −3 ) makes the optimization of f ↑↓ ij parameters in the pair-product (geminal) wave parameters in the PP wave function are always set to be a pair of ↑ and ↓ spins. We impose 2 × 2 sublattice structure in f ↑↓ ij parameters in the PP (geminal) wave function. In this case, the number of independent f ↑↓ ij parameters becomes 2 × 2 × Nsite = 4Nsite, and the other f ↑↓ ij parameters are determined by using spatial translation operations. The RBM part ({b k , W ik }) is taken to be fully translationally invariant (1 × 1 sublattice structure) [49]. When the double occupancy is completely prohibited by the Gutzwiller factor P ∞ G = i (1 − n i↑ n i↓ ), the onsite f ↑↓ ii parameters become completely irrelevant, i.e., the wave function and the energy do not depend on f ↑↓ ii at all. Thus the number of f ↑↓ ij parameters in the table are reduced from 4Nsite to 4(Nsite −1).
As for the initial {b k , W ik } parameters, we use small random numbers. We run several calculations with different seeds for generating random numbers and adopt the wave function that has the lowest energy.

Particle-hole transformation of Hubbard model
When we analyze many-body Hamiltonians by theoretical or numerical solvers, we can utilize transformations of the Hamiltonians to find a representation suitable for the solver in hand [86]. In the conventional variational Monte Carlo method, many-body interacting models defined in the real-space basis such as the Hubbard model is analyzed as they are, because the empirical form of the correlation factors is also defined in the real-space basis. On the other hand, in cases where we employ machine learning technique, because the neural network will find a way to lower the energy even when the form of the Hamiltonian is complicated, we could think of the "best" transformation of the Hamiltonian such that the neural network can lower the energy with small number of parameters.
Though it would not be the best of the best, we find that performing a staggered particle-hole transformation (c i↓ → (−1) i c † i↓ and c † i↓ → (−1) i c i↓ ) and mapping onto the attractive Hubbard model help the RBM+PP wave function to lower the energy of the Hubbard model. This transformation is also helpful for mVMC calculations to get a better energy [87]. Thus, we have solved the Hubbard model with this transformation. Finding a better transformation (by again employing machine learning technique) is an interesting future problem.
Appendix D: U dependence of energy in RBM+PP Figure 4 shows the U dependence of RBM+PP (N L K=0 |φ pair with α = 32) energy of the 2D 8 × 8 Hubbard model. The result shows that the error decreases with increasing U/t.