{'fraction_non_alphanumeric': 0.03310713551305779, 'fraction_numerical': 0.017684556857906642, 'mean_word_length': 4.903793626707132, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Network-Initialized Monte Carlo Based on Generative Neural Networks

We design generative neural networks that generate Monte Carlo configurations with complete absence of autocorrelation from which only short Markov chains are needed before making measurements for physical observables, irrespective of the system locating at the classical critical point, fermionic Mott insulator, Dirac semimetal, or quantum critical point. We further propose a network-initialized Monte Carlo scheme based on such neural networks, which provides independent samplings and can accelerate the Monte Carlo simulations by significantly reducing the thermalization process. We demonstrate the performance of our approach on the two-dimensional Ising and fermion Hubbard models, and expect it can systematically speed up the Monte Carlo simulations especially for the very challenging many-electron problems.

Introduction -Monte Carlo (MC) method is a widely used numerical method to investigate problems in statistical and quantum many-body physics. Although being generally polynomial when the sign problem is absent, the complexity of MC simulation can still be prohibitively high for certain challenging problems. In particular, quantum MC simulations of generic interacting fermion systems [1][2][3][4][5], such as Hubbard model for correlated electrons [2,3], Holstein model for electron-phonon interaction [6][7][8][9], the spin-fermion model of non-Fermi liquid [5,[10][11][12][13], and the momentum space QMC for quantum moiré systems [14][15][16][17], to name a few, have very high complexity that scales at least to ∼ 3 where = 1/ is the inverse temperature and the = the total site number for a lattice model in -dimension. Such type of MC approach for the many-electron problems is termed as determinant quantum Monte Carlo (DQMC) where the complexity comes from the matrix operation on the fermion determinant. The situation worsens when dealing with critical points like the novel quantum critical points ubiquitously present in the metal-insulator transition, pseudogap and non-Fermi liquids, and dynamics and transport properties in high-temperature superconductivity, etc. This can largely be ascribed to that the widely-used Markov-chain Monte Carlo (MCMC) method with local updates (or the Metropolis algorithm) suffers from very long Markov-chain autocorrelation that scales with the system size to a high power due to the critical slowing down.
Although there exist powerful and highly efficient cluster/loop update schemes for classical problems such as the Ising model [18,19] and quantum spin/boson systems such as Heisenberg and Bose-Hubbard models [20][21][22][23], these ingeniously designed algorithms do not exist for generic quantum many-body problems, especially for the interacting fermion systems, for example in DQMC. Therefore, the efforts for more efficient algorithms that could overcome these difficulties and boost the (quantum) MC simulations for the aforementioned systems are called for.
Inspired by these recent developments, we design generative neural networks to approximate the distribution of MC configurations, and find the network can provide configurations with complete absence of autocorrelation, thus serving as a good starting point for MCMC with significantly shortened Markov chain. We demonstrate the powerful fitting ability of the neural networks by successfully generating configurations from which the correct physical observables can be directly measured after a short Markov chain, as showcased in both the classical = 2 Ising model and quantum critical point of = 2 fermion Hubbard model. In particular, the latter cannot be mapped into a classical model with short-range interactions due to the existence of gapless Dirac fermions, but can still be successfully dealt with by the network. This suggests that our approach can be applied to quantum many-body systems, in particular the interacting fermion models.
Moreover, we design a network-initialized Monte Carlo (NIMC) scheme with the assistance of such neural networks,

MCMC Simulation
Observables Ⅰ which "heals" the initial bias due to the unthermalized training set of the neural networks and thus provides very accurate results. We find the NIMC scheme with independent samplings can accelerate simulations by reducing the long thermalization time and measurement processes, which is especially important for simulating systems with long autocorrelation time in traditional MCMC.
Models -In this work, we aim to train neural networks to generate MC configurations that approach the distribution of the original Hamiltonian using convolutional and transposed convolutional architectures. Different from the reference work such as autoregressive models and referring to the idea of the Generator in Generative Adversarial Nets (GAN) [59], we use random configurations (each element of the configurations is randomly set to be ±1 with equal probability) as input, which, after operations inside the networks, turn into the output of MC configurations distributed according to the physical Hamiltonian. The difference of our approach from GAN is that we do not need to build any discriminator, but directly compare the physical observables, like the internal energy and magnetization (or magnetic structure factor) measured from a batch of configurations generated by the neural network with those obtained from MC simulations to optimize the parameters inside neural networks, as demonstrated in the flow diagram for training in Fig.1.
For the 2d square lattice Ising model with we study the system at the critical point ( ≈ 2.269 ) with system sizes of 16×16 and 32×32. We design the loss function directly related to the physical observables, i.e., where is the -th generated configuration, and refer to the magnetization 1 | | and the energy density 1 measured from the corresponding -th MCMC configuration, and 1 , 2 , and 3 are constants to balance each part of the loss to be in the same magnitude.
For the 2d Hubbard model with we consider the half-filled honeycomb lattices where it experiences a chiral Heisenberg Gross-Neveu quantum critical point at / ≈ 3.83 between the Dirac semimetal and the antiferromagnetic Mott insulator [60][61][62][63][64]. Here the configuration is made of the auxiliary fields used in the discrete Hubbard-Stratonovich transformation to decouple the fermion interaction terms in DQMC, so the configurational space is of the size × × compared with the × for the Ising case (see detailed explaination in Sec. I of Supplemental Materials (SM) [65]). We then design the loss function as: where is the -th generated configuration, refers to the magnetic structure factor (so ( ) and ( ) ( ) refer to the structure factor measured from the -th training sample and from the -th generated configuration) and (Γ = (0, 0)) is chosen for honeycomb lattice as the Γ = (0, 0) is the ordered wave vector for the antiferromagnetic long-range order, and is the kinetic energy 1 , , , + ℎ. .) , with constants 1 and 2 .
The design of the loss functions comes from physical intuition, and we choose the key physical observables such as the energy and magnetic structure factor, which reflect the phase transition of our interests in the loss function. Similar considerations can be generalized to other problems. In practice, our networks are all based on Tensorflow [66], and we complement with more details of the Ising and Hubbard models in Sec.I and those of neural networks in Sec.II of SM [65].
Autocorrelation analysis -In MCMC, the autocorrelation time, associated with the particular update scheme, can be extremely long in the physically interesting parameter regime like, e.g., near the classical and quantum critical points. And it is even worse when one performs the finite size analysis as the autocorrelation time usually increase with system size with a high power [8,18,41]. It therefore benefits a lot that the independent and direct samplings from the neural networks are completely free of autocorrelation. We demonstrate these behaviours quantitatively by measuring the autocorrelation function (defined in Sec.I.3 of SM [65]) of magnetization for 2d Ising model at = 2.269 in Fig. 2(a), and of (0, 0) for the honeycomb-lattice Hubbard model at the Gross-Neveu QCP ( / = 3.83) in Fig. 2(b), between the results generated by neural networks and MCMC results. It can be observed that, for the MCMC the autocorrelation functions decay slowly in Monte Carlo steps and the autocorrelation time for classical (quantum) critical points increases with the system size, but in results from neural networks, such autocorrelation has been completely eliminated.
NIMC algorithm based on neural networks -We now introduce the NIMC algorithm. In order to demonstrate its ability of correcting the error in the training sample, we no longer train the networks to approximate well-thermalized MC results, but take unthermalized MC configurations generated from short Markov chains. After training the neural networks and generate independent configurations, we add MCMC updates starting from these generated configurations (one Markov chain for each configuration) [67] and find the NIMC results quickly converge to the well-thermalized results. The whole  Numerical experiments on the Ising model -In the standard MCMC, the Markov chains must first be well thermalized such that it converges to the correct distribution, and then generate enough bins of statistically uncorrelated configurations to yield expectation values of observables with satisfactory statistical error, which requires the thermalization and each bin to be longer than the system-size-dependent autocorrelation time. To demonstrate that NIMC largely reduces the thermalization process, we compare the performance of NIMC  . We run 50 iterations for all three sets of data with each set consisted of 1000 configurations (for each iteration, we re-run Markov chains from all the configurations), and the results are shown in Fig. 3. For training samples, almost about 150 sweeps are needed to reach a convergence, and for random configurations the number is more than 200. Remarkably, the NIMC scheme greatly lowers the thermalization time and arrive at the convergence after about only 30 sweeps.
To further quantify the comparison, we take the training part for neural networks into consideration, and compute on the same CPU (Intel Xeon E5-2680) that for the 150 epochs of network training that are more than sufficient to converge the parameters in the neural network. It turns out that the training process takes around 125 seconds (could be made even faster on GPU), and each MCMC sweep for 1000 configurations takes around 2.5 seconds. That is, NIMC takes 125 + 2.5 × 30 = 200 seconds for the 1000 thermalized and independent configurations here from training samples, while the 1000 samples themselves need up to 2.5 × 150 = 375 seconds for thermalization. Therefore, we claim that NIMC can overall speed up the MC simulations. Besides, the independent NIMC samples also save time in the process of measurements as there is completely no autocorrelation and we can directly measure observables once the thermalization is completed.
To better compare the accuracy, we then show the observables by taking 1000 configurations from the beginning of a rather short Markov-chain as traing samples to run NIMC and the numerical results are shown in Tab. I. Besides, we also plot the corresponding histograms of distribution of magnetization in Fig. 4 after 50 iterations of NIMC (1000 configurations for and corresponds to the distributions to be compared and we divide them into 20 bins. Again we count the NIMC samples from 50 iterations for statistical expectation with standard error. The result is % = 94.78 (18), which reconfirms that the distribution of NIMC with just 1000 configurations has a large overlap with that of Good MC. Hubbard model on the honeycomb lattice -Beyond the classical Ising model, here we take the 2d Hubbard model on a honeycomb lattice at its unique Gross-Neveu quantum criticality between the Dirac semimetal and antiferromagnetic Mott insulator ( / = 3.83, = 6) as an example of in-teracting fermion systems. We run NIMC with 15 DQMC sweeps from 1000 generated configurations and the mean values with errorbars of physical observables are listed in Tab. II. The NIMC results again present almost identical expectation values of physical observables as the Good DQMC (containing 10000 configurations in total that are divided into 1000 bins). Furthermore, in Fig. 5 we show the magnetic structure (Q) along the high-symmetry path of the Brilliouin zone, which are observables not used in training the networks. We find again excellent agreements between Good MC and NIMC only after 15 sweeps. The peak at the Γ point highlight the quantum critical fluctuations towards the magnetic order that will gap out the Dirac cones. These results confirm the neural networks with optimizers like ADAM [68] are capable of fixing the large-scale determinant computations. We have to clarify that in NIMC we do not necessarily need up to 1000 configurations (as we use for example here), which can further reduce the computational cost.
Discussion -By developing a NIMC method based on the generative neural networks, our numerical results support firmly that the latter as independent and direct sampling approach can approximately capture the weight distributions of MC configurations and completely get rid of autocorrelation in the study of classical and quantum many-body systems. NIMC therefore provides a scheme for speeding up the MC simulations in reducing the thermalization time and saving time for direct measurement of observables from those independent NIMC samples. Admittedly, our NIMC still needs a few steps of the MCMC simulation, therefore it does not change the overall scaling of the computational cost but can nevertheless lead to a significant factor reduction.
From the perspective of neural networks, we verify that large-scale simulations, such as DQMC calculations, can be implemented through the evaluation of the neural network loss functions. The next step from here is to use generative neural networks to provide efficient and accurate system-size extrapolations of MC configurations and in this way many more complicated yet fundamental problems, such as the quantum moiré materials that the DQMC has just been shown to be able to solve [14][15][16][17][69][70][71] but with heavy computational costs, could be improved in NIMC.

Determinant Quantum Monte Carlo (DQMC)
The partition function = Tr{ − } is expressed as a path integral by discretizing the inverse temperature into slices of length Δ . Then, after Trotter-Suzuki decomposition, the Hamiltonian are separated in each time slice and can be written as To treat the quartic interaction term, the discrete Hubbard-Stratonovich transformation [1][2][3] can be applied and where the scalar is defined by cosh = Δ 2 . Putting kinetic part and interaction part together, the term in partition function becomes where is rewritten with operator = − , † , , + ℎ. .. According to the feature of fermion operator [1], a fermion operator (ˆfor example) with a quadratic form likeˆ= , † ( ) satisfies Then the partition function could be finally written in determinant form as in which where = {1, −1} in calculation corresponding mark {↑, ↓}, and is the matrix corresponding to the operator . Thus, we have introduced a sum over the field of auxiliary variables , in a ( + 1)-dimension space ( for spatial denoted by and 1 for imaginary time denoted by ) as shown in Fig. 6 (b), and the fermionic degrees of freedom in the quadratic form have been integrated out analytically. Note is an × matrix that depends on the auxiliary configurations.

Autocorrelation Function
The autocorrelation function for an observable is defined as: For MCMC, and denotes the simulation time, normally in units of the MC sweeps (one sweep means doing flipping attempts over all the spins of the configuration), and the averages are over the reference time . For Neural Network results, when computing the formal autocorrelation function, and simply refers to the serial numbers of the configurations. For 2d Ising model on square lattice (16 × 16 for example), the network structure is shown in Fig. 6 (a). For the random input, we fix its shape to be the same as that of the Monte Carlo configurations. Inside the network, inspired by the network structure in Ref. [30], we use three 2d convolutional layers, with 64 filters of size 8 × 8 for the first layer, 32 filters of size 1 for the second, and 1 filter of size 6 × 6 for the third layer. We choose rectified linear function ReLU(x) = max(0, x) in the first and second convolutional layers and sigmoid function ( ) = 1/(1 + − ) in the third. We also apply periodic boundary condition (PBC) layers each time before the convolutional layers if the kernel dimension is larger than one, which provides the configuration tensors with paddings of boundary elements instead of 0s (the padding size is decided by the size of filters in the corresponding convolutional layer, so as to ensure the shape invariance of the configuration tensors after convolutional operations).
In the case of 2d Hubbard model on the honeycomb lattice (6 × 6 × 2 for example), we switch from classical Ising configurations to the auxiliary field configurations of the DQMC [3,33], with the configuration space of × × × 2 where = 1/ is the inverse temperature. We set = 6 and the Trotter discretization Δ = 0.1, so the DQMC auxiliary field configurations are of shape 60 × 12 × 6 (note the 2 site per unit cell for the honeycomb lattice), as schematically shown in Fig. 6 (b). As shown in Eq.(4), we have to carry out large-scale determinant computation to measure the observables from generated configurations in order to compute the loss, and thus the complexity for the optimization is much larger than the case of classical model. Therefore, instead of the previous PBC-based network structure, we try another architecture of transposed convolutional layers and feed smaller random configurations as input. We find it is still capable of fitting the results we want and optimize faster. As shown in Fig. 6 (c) for the example model, we build two 3d transposed convolutional layers, with all filters using valid padding and sigmoid function as activation. There are 20 filters of size 16 × 6 × 3 in the first layer and 1 filter of size 36 × 6 × 3 in the second and the random input configurations are of shape 10 × 2 × 2 while the predicted output are of shape 60 × 12 × 6.

Optimization
In order to optimize such neural networks, for both classical and quantum models, we prepare 1000 sets of observables measured from MC simulations and take 1000 input random configurations which will be processed into generated configurations. The comparison in loss function is randomly distributed without any grouping. For the choice of observables in the defined loss functions, we simply pick some from the ones we usually focus on. In the Ising case, we take the batch size to be 5 and epoch number to be up to 150 as the computation is quite easy. While in the Hubbard case, we run at most 15 epochs with a batch size of 3. We optimize the network parameters using Adam [63] with conventional learning rate 10 −3 , 1 = 0.9, and 2 = 0.999.