Accelerating lattice quantum Monte Carlo simulation using artificial neural networks: an application to the Holstein model

TL;DR

This paper introduces neural network-enhanced Monte Carlo methods that significantly speed up simulations of the Holstein model by learning effective models, enabling faster and scalable lattice quantum Monte Carlo computations.

Contribution

The authors develop neural network-based local and global move strategies for lattice Monte Carlo simulations, achieving substantial speedups and broad applicability to quantum and classical models.

Findings

Neural networks accurately reproduce MC configuration weights.

Achieved an order of magnitude speedup over traditional DQMC.

Applicable to various classical and quantum lattice MC algorithms.

Abstract

Monte Carlo (MC) simulations are essential computational approaches with widespread use throughout all areas of science. We present a method for accelerating lattice MC simulations using fully connected and convolutional artificial neural networks that are trained to perform local and global moves in configuration space, respectively. Both networks take local spacetime MC configurations as input features and can, therefore, be trained using samples generated by conventional MC runs on smaller lattices before being utilized for simulations on larger systems. This new approach is benchmarked for the case of determinant quantum Monte Carlo (DQMC) studies of the two-dimensional Holstein model. We find that both artificial neural networks are capable of learning an unspecified effective model that accurately reproduces the MC configuration weights of the original Hamiltonian and achieve an order of magnitude speedup over the conventional DQMC algorithm. Our approach is broadly applicable to many classical and quantum lattice MC algorithms.

Equations

$$p_c=\frac{W(\{X^{\prime }\})}{W(\{X\})}\frac{\exp (-\beta E_{\text{NN}}[\{X\}])}{\exp (-\beta E_{\text{NN}}[\{X^{\prime }\}])}.$$

Full text

Introduction — As their full potential becomes apparent, machine learning algorithms are assuming more prominent roles in the process of scientific discovery. Meanwhile, the boundary lines between industry applications of machine learning, data and computer science, and other disciplines have blurred. Applications ranging from the high-quality feature extraction from astrophysical images of galaxies Schawinski et al. (2017) to helping with the real-time data analysis of particle accelerators Edelen et al. (2016); Radovic et al. (2018); Scheinker et al. (2018) at the Fermilab Edelen et al. (2016) and the Large Hadron Collider Radovic et al. (2018) to discovering phases of matter Carrasquilla and Melko (2017); Wang (2016); Zhang et al. (2018) have emerged.

A series of early studies have underscored the potential for machine learning in the context of condensed matter physics by using artificial neural networks and dimension-reduction techniques to locate phase transitions Carrasquilla and Melko (2017); Wang (2016), or represent ground states of quantum many-body systems Carleo and Troyer (2017). Machine learning algorithms have also been employed to help gain more insight into classical and quantum systems Torlai and Melko (2016); Ch’ng et al. (2017); Deng et al. (2017); van Nieuwenburg et al. (2017); Zhang and Kim (2017); Wetzel (2017); Hu et al. (2017); Torlai et al. (2018); Schindler et al. (2017) as well as accelerate specific numerical algorithms Huang and Wang (2017); Liu et al. (2017a); Li and Wang (2018); Fournier et al. (2018); Nelson et al. (2018); Suwa et al. (2018). These applications are not only helping to automate and streamline scientific processes that could take many years to accomplish with more conventional computational approaches, but they are also uncovering previously inaccessible phenomena.

One machine learning application that has attracted significant attention is in accelerating MC simulations Huang and Wang (2017); Liu et al. (2017b, a); Xu et al. (2017); Nagai et al. (2017); Shen et al. (2018); Chen et al. (2018); Nagai et al. (2018); Inack et al. (2018); Song and Lee (2019). For example, in the so-called self-learning Monte Carlo (SLMC) method Liu et al. (2017b), an effective bosonic model is trained to mimic the statistics of the original Hamiltonian. Once trained, the effective model is then used to perform the same simulations much more efficiently. The primary advantage of this approach is that the action of the effective model is often much easier to compute than the action for the full fermion model, thus granting access to larger system sizes. This approach has also been extended to include correlations in both the real space and imaginary time domains Xu et al. (2017); Nagai et al. (2017). Despite their power, however, the SLMC methods require that the form of the effective model be known a priori. This limitation can be significant as different effective models may be required for the same fermionic model as the model parameters, system size, or simulation temperature changes, and the overall effectiveness of these approaches is severely limited if the wrong effective model is chosen. To overcome this problem, several groups have used artificial neural networks to learn the form of the model in some instances Shen et al. (2018); Nagai et al. (2018); Inack et al. (2018); Song and Lee (2019) (e.g. QMC simulations of the Anderson impurity model); however, generalizing this approach to lattice QMC problems has not yet been achieved. One reason for this is the fact that such problems typically involve thousands of auxiliary spacetime fields and any neural network using that many input features will often generalize poorly.

Here, we show how to design artificial neural networks that can be trained to represent an effective bosonic model for lattice QMC simulations. Inspired by applications of the traveling cluster approximation to spin-fermion models Sen et al. (1984); Mukherjee et al. (2015), we design fully-connected and convolutional neural networks that only require information from surrounding auxiliary fields (see. Fig. 1) to perform both local and global moves of the MC configurations. This method does not suffer from the scaling issues restricting other self-learning methods and can be easily generalized across models and parameter regimes without changes in the underlying algorithm, provided the neural networks are versatile enough to learn the effective models. As such, this approach can be integrated directly into existing MC codes. To demonstrate the efficiency of this approach, here we apply it to determinant quantum Monte Carlo (DQMC) simulations of the two-dimensional Holstein model. This problem is particularly challenging owing to long autocorrelation times Hohenadler and Lang (2008), the need for both local and global MC moves to ensure ergodicity Scalettar et al. (1991); Johnston et al. (2013a), and competition between charge-density-wave (CDW) and superconducting ground states Dee et al. (2019) that may require different effective boson models. Reproducing known results, we obtain an order of magnitude of speedup with our algorithm.

Model — The single-band Holstein Hamiltonian Holstein (1959) is $H=H_{0}+H_{\text{lat}}+H_{\text{e-ph}}$, where $H_{0}=-t\sum_{\langle i,j\rangle,\sigma}c_{i,\sigma}^{\dagger}c_{j,\sigma}^{\phantom{\dagger}}-\mu\sum_{i,\sigma}\hat{n}_{i,\sigma}$, $H_{\text{lat}}=\sum_{i}\left(\frac{1}{2M}\hat{P}_{i}^{2}+\frac{M\Omega^{2}}{2}\hat{X}_{i}^{2}\right)$ , and $H_{\text{e-ph}}=g\sum_{i,\sigma}\hat{X}_{i}\hat{n}_{i,\sigma}$. Here, $\langle\cdots\rangle$ denotes a summation over nearest neighbors; $c_{i,\sigma}^{\dagger}$ ($c_{i,\sigma}^{\phantom{\dagger}}$) creates (annihilates) an electron with spin $\sigma$ on site $i$; $\hat{n}_{i,\sigma}=c_{i,\sigma}^{\dagger}c^{\phantom{\dagger}}_{i,\sigma}$ is the particle number operator; $t$ is the nearest-neighbor hoping integral; $M$ is the ion mass and $\Omega$ is the phonon frequency; $\hat{X}_{i}$ and $\hat{P}_{i}$ are the lattice position and momentum operators, respectively; and $g$ is the strength of the e-ph coupling. Throughout, we set $M=t=1$ as the unit of mass and energy, and we study this Hamiltonian on an $N=N_{x}^{2}$ square lattice, where $N_{x}$ is the linear size of the cluster. To facilitate a direct comparison with a recent state-of-the-art simulation Chen et al. (2018), we focus on $\Omega=t/2$ and dimensionless e-ph coupling strength $\lambda=\frac{g^{2}}{8t\Omega^{2}}=0.5$.

Determinant quantum Monte Carlo — DQMC is an auxiliary field, imaginary time technique that computes expectation values of an observable within the grand canonical ensemble. In a DQMC simulation, the imaginary time interval $\tau\in[0,\beta]$ is evenly divided into $L$ discrete slices of length $\Delta\tau=\beta/L$ ($=0.1$ in this work). Using the Trotter approximation, the partition function is then given by $Z=\mathrm{Tr}\left(e^{-\Delta\tau LH}\right)\approx\mathrm{Tr}\left(e^{-\Delta\tau H_{\text{{\it e}-ph}}}e^{-\Delta\tau(H_{0}+H_{\mathrm{lat}})}\right)^{L}$. After integrating out the electronic degrees of freedom, the partition function can be reduced to $Z=\int W\left(\{X\}\right)\mathrm{d}X$, where the configuration weight is $W\left(\{X\}\right)=e^{-S_{\mathrm{ph}}\Delta\tau}\mathrm{det}M^{\uparrow}\mathrm{det}M^{\downarrow}$. Here, $\int\mathrm{d}X$ is shorthand for integrating over all of continuous displacements $X_{i,l}$ defined at each spacetime point $(i,l)$, the matrices $M^{\sigma}$ are defined as $M^{\sigma}=I+B_{L}^{\sigma}B_{L-1}^{\sigma}\cdots B_{1}^{\sigma}$, where $I$ is an $N\times N$ identity matrix, and $B_{l}^{\sigma}=e^{-\Delta\tau H_{\text{{\it e}-ph}}}e^{-\Delta\tau H_{0}}$, and $S_{\mathrm{ph}}=\frac{M}{2\Delta\tau^{2}}\sum_{i,l}\left(X_{i,l+1}-X_{i,l}\right)^{2}+\frac{M\Omega^{2}}{2}\sum_{i,l}X_{i,l}^{2}$ is lattice’s contribution to the total action. Note that $B^{\sigma}_{l}$ matrices for the Holstein model do not depend on spin but are implicitly dependent on the fields $X_{i,l}$ through $H_{\text{{\it e}-ph}}$. For more details, we refer the reader to Refs. White et al., 1989; Scalettar et al., 1989; Johnston et al., 2013b.

As mentioned, two types of MC updates are needed in the simulation. The first are local updates of the type $X_{i,l}\rightarrow X^{\prime}_{i,l}=X_{i,l}+\Delta X_{i,l}$, which are made at each spacetime point. The second are global or block updates, where the field for a given site are updated simultaneously at all timeslices $X_{i,l}\rightarrow X_{i,l}+\Delta X_{i},~{}\forall l\in[0,L]$. Such block updates are needed to help move phonon configurations out of local minima at low temperatures and large couplings  Scalettar et al. (1991); Johnston et al. (2013a). DQMC accepts both kinds of moves with a probability $p=W\left(\{X^{\prime}\}\right)/W\left(\{X\}\right)\equiv e^{-\beta\Delta E}$, which requires the costly evaluation of matrix determinants. Moreover, since the matrices $M^{\sigma}$ depend on the fields, these must also be updated after every accepted change in the phonon fields. While an efficient update algorithm exists for performing local updates White et al. (1989), no such algorithm is known for block updates. The computational cost for performing a full sweep of (fast) local updates and block updates is $\mathcal{O}(N^{3}L)$ and $\mathcal{O}(N^{4}L)$ Lee et al. (2010), respectively.

To reduce this cost, we train our networks to predict $\beta\Delta E$ appearing in the definition of the configuration weight given only changes in, and local information of, the phonon fields and their expected behavior at large displacements as input. This reduces the total computational complexity of determining whether both kinds of updates will be accepted to the time needed to evaluate the networks, which is $\mathcal{O}(1)$ for the case of the simply connected network and $\mathcal{O}(L)$ for the convolutional neural network. As with other SLMC methods, we then use the neural networks to propose many MC updates that are ultimately accepted or rejected based on the configuration weights of the original model. While determining this final acceptance probability requires the evaluation of the matrix determinants, this task can be done infrequently enough that a considerable speedup is achieved. Another advantage of our approach is that the networks can be trained using data generated by the conventional DQMC algorithm on inexpensive small clusters before being generalized to larger systems. In this way, our method combines the flexibility of neural networks with the inexpensive training costs seen in SLMC approaches making use of largely local effective models.

Local Updates — Local updates are performed using a fully-connected network with two hidden layers [Fig. 1(a)]. Assuming that the update is proposed at spacetime site $(i,l)$, the learning objective is to predict $\beta\Delta E$ given only $\Delta X_{i,l}$ and the field values at the surrounding spacetime points as input features. Here, we include nearest- and next-nearest neighbor phonon fields in both space and imaginary time, and neglect long-range correlations. While there is justification for a short-range effective interaction in proximity to the CDW phase at half-filling Dee et al. (2019), this approximation can also be systematically improved by taking more input features. We have found, however, that next-nearest-neighbor inputs are sufficient. We also supply an additional neuron in the input layer that enforces known behavior at large displacements Chen et al. (2018); Sup (2019).

Global Updates — Global updates are performed using a convolutional neural network (CNN) with four layers [Fig. 1(b)], where the objective again is to predict $\beta\Delta E$ given only local information about the phonon fields. Assuming the update occurs at site $i$, the input layer has three columns of input features: the first contains fields $X_{i,l}$ across all imaginary time slices; the second and third columns contain averages $\bar{X}_{i,l}^{(1)}=\frac{1}{4}\sum_{\langle j\rangle}X_{j,l}$ and $\bar{X}^{(2)}_{i,l}=\frac{1}{4}\sum_{\langle\langle j\rangle\rangle}X_{j,l}$, respectively, at all time slices, where $\langle j\rangle$ and $\langle\langle j\rangle\rangle$ denote nearest- and next-nearest-neighbor sums around site $i$. The use of $\bar{X}_{i,l}^{(1)}$ and $\bar{X}^{(2)}_{i,l}$ enforces $C_{4}$ rotational symmetry and reduces the cost of training the CNN. The convolution operation from the input layer to the first hidden layer is standard Sup (2019).

For each set of (fixed) input parameters $\{\beta,\mu,\Omega,g\}$ we train both networks using training examples generated with the conventional DQMC algorithm on a $6\times 6$ cluster. Throughout, we generated $8\times 10^{4}$ samples, which were randomly partitioned into $6\times 10^{4}$ training and $2\times 10^{4}$ test samples. We first show results for their performance; the insets of Figs. 1(a) and 1(b) compare the predicted $\beta\Delta E_{\text{NN}}$ against the exact $\beta\Delta E$ values obtained from our test data sets for the local and global updates, respectively. This simulation was performed close to the CDW transition for the model [Fig. 3]. Both networks accurately predict the MC configure weights but the fully-connected neural network is slightly more accurate. While the accuracy can be systematically improved by taking more input features, we find that the knowledge learned by both networks can be transferred to larger clusters remarkably well based on Fig. 2 (a).

Once our networks have been trained and tested, we then define a full MC sweep as consisting of $N_{\text{u}}$ complete sweeps of local updates performed at each spacetime point ($i,l$) using the fully-connected neural network, followed by $N_{\text{u}}$ sweeps of global updates performed at every lattice site $i$ using the CNN. (This sampling procedure differs from the conventional one Johnston et al. (2013a), where global updates are performed on a subset of sites to minimize the total computational cost.) After performing these sweeps, the original field configuration $\{X\}$ is replaced with a newly proposed one $\{X^{\prime}\}$ in a cumulative update Liu et al. (2017a, b); Xu et al. (2017) with a probability $\mathrm{min}\left[1,p_{c}\right]$, where

[TABLE]

Benchmarks — To benchmark the NNMC, we performed direct comparisons with the conventional DQMC algorithm for the half-filled model $\langle\hat{n}\rangle=1$ at $\beta=4.1/t$, which is close to the CDW transition temperature for this parameter set. We emphasize that both the DQMC and NNMC simulations used the same sampling protocol with $N_{\text{u}}=1$. Figure 2(a) plots the CDW structure factor $S({\bf q})$ at ${\bf q}=(\pi,\pi)/a$ Sup (2019) as a function of the linear cluster size $N_{x}$, and demonstrates that the NNMC algorithm accurately reproduces the results of the conventional DQMC algorithm for the accessible lattice sizes. Figure 2(b) compares the autocorrelation time $\tau_{L}$ of $S_{\text{CDW}}({\bf q})$ for both techniques, which again yields similar results. We note, however, that the autocorrelation time can be reduced significantly by increasing the number of update sweeps $N_{\text{u}}$ that are performed before computing the cumulative update acceptance probability, as shown in the inset of Fig. 2(b).

To address how NNMC reduces the computational cost, we compare the time to solution for both algorithms in Fig. 2(c). Fitting a power law $t_{\text{CPU}}=\alpha N^{z}$ to the data yields $z=3.41$ and $2.35$ for DQMC and NNMC, respectively, a significant reduction in the scaling. We note that a similar speedup was obtained using SLMC Chen et al. (2018); however, the NNMC does not require the functional form of the effective model to be specified a priori. Moreover, the NNMC method is more efficient at generating accepted MC moves, particularly as it is generalized to larger system sizes. We highlight this aspect in Fig 2(d), which shows the cumulative acceptance ratio $p_{c}$ obtained using NNMC and compares it with SLMC. As the methods are generalized to larger cluster sizes, $p_{c}$ decreases for the SLMC method while the NNMC method proposes cumulative moves that are almost always accepted, and becomes more accurate on larger cluster sizes. The decrease of $p_{c}$ in SLMC is due to the poor performance of the regression model for predicting global updates, which requires a more sophisticated effective model.

We now demonstrate that the NNMC approach can also be used to study the finite-size scaling of the CDW structure factor and obtain the transition temperature in the thermodynamic limit. Fig. 3 presents a similar analysis carried out in the same temperature region $\beta=3.8/t\sim 4.5/t$ using the NNMC approach. At the critical point, the finite-size scaling behavior has the form $S_{\text{CDW}}(\pi,\pi)/N_{x}^{2}=N_{x}^{-2\gamma/\nu}f\left(N_{x}^{1/\nu}\frac{T-T_{\text{c}}}{T_{\text{c}}}\right)$, where $\gamma=\frac{1}{8}$ and $\nu=1$ are the 2D Ising critical exponents. The critical temperature $T_{\text{c}}/t\approx 0.244$ ($\beta_{c}=4.1/t$) is determined by the common intersection point of the curves. The inset of Fig. 3 replots $S_{\text{CDW}}N_{x}^{-7/4}$ against $N_{x}^{1/\nu}\frac{T-T_{\text{c}}}{T_{\text{c}}}$, showing the expected data collapsing to a single curve, consistent with Ref. Chen et al. (2018).

Summary and Conclusions — We have extended the use of artificial neural networks in self-learning Monte Carlo methods to lattice Monte Carlo simulations. Our approach overcomes many of the scaling issues associated with other SLMC implementations and can be widely applied to classical and quantum Monte Carlo simulations on extended lattices. We then applied this methodology to DQMC studies of the Holstein model. In doing so, we designed fully-connected and convolutional neural networks capable of performing accurate local and, for the first time, global moves in configuration space. Using this method we are able to reproduce results of the charge-density-wave transition in this model but using an approach that does not require the form of the effective model to be specified in advance.

Our approach constitutes a generalizable method for performing machine-learning accelerated lattice Monte Carlo simulations (even the Fermi-Hubbard model), provided that the neural networks are sophisticated enough to learn the underlying effective model for the relevant parameter ranges. In this context, we note that the specific network architectures depicted in Fig. 1 accurately predict the CDW transition; however, they are not well-suited to predict the metal-to-superconducting transition away from half-filling. This is also true for the effective model used in the SLMC approach as applied to the Holstein model in Ref. Chen et al. (2018). To learn and predict $\beta\Delta E$ near the superconducting transition, different network architectures that can capture the physics in this region should be considered. Finally, we stress that our approach allows to compute physical quantities that can be compared with experiments. Previous quantum Monte Carlo + machine learning approaches used to study the Hubbard model on extended lattices Ch’ng et al. (2017) have generally focused on classifying phases in the parameter space, while physical measurable quantities (such as one- and two-particle functions) are inaccessible.

Acknowledgements — S. L., P. D., and S. J. are supported by the Scientific Discovery through Advanced Computing (SciDAC) program funded by the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research and Basic Energy Sciences, Division of Materials Sciences and Engineering. S. L. is also supported by the National Science Foundation (NSF) under Grant No. DMR-1606348. E. K. is supported by the National Science Foundation Grant No. DMR-1609560. An award of computer time was provided by the INCITE program. This research also used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725, and the Spartan high-performance facility at San Jose State University, which is supported by the NSF under Grant No. OAC-1626645.