Accelerated Continuous time quantum Monte Carlo method with Machine Learning

TL;DR

This paper introduces a machine learning-enhanced continuous time quantum Monte Carlo method that significantly speeds up computations while accurately predicting impurity Green's functions and double occupancy at low temperatures.

Contribution

It combines machine learning with CTQMC to eliminate matrix multiplications, reducing computational time and maintaining accuracy for strongly correlated materials.

Findings

Predicts impurity Green's function accurately at low temperatures

Reduces computational time compared to conventional CTQMC

Maintains physical properties at high Matsubara frequencies

Abstract

An acceleration of continuous time quantum Monte Carlo (CTQMC) methods is a potentially interesting branch of work as they are matchless as impurity solvers of a density functional theory in combination with a dynamical mean field theory approach for the description of electronic structures of strongly correlated materials. The inversion of the $k \times k$ matrix given by the diagram expansion order $k$ in the CTQMC update and the multiplication of the $k \times k$ matrix and the non-interacting Green's function to measure the impurity Green's function are computationally time-consuming. Here, we propose the CTQMC method in combination with a machine learning technique, which would eliminate the need for multiplication of the matrix with the non-interacting Green's function. This method predicts the accurate impurity Green's function and double occupancy at low temperature, and also considers the physical properties of high Matsubara frequency in a much shorter computational time than the conventional CTQMC method.

Equations

$$Z=\int D[d^{†},d]\exp (-S_0-S_U),$$

$$\frac{Z}{Z_0}=\sum _{C_k}(-U{)}^k<n_{1\uparrow }^{{}^{\prime }}({\tau }_1)n_{1\downarrow }^{{}^{\prime }}({\tau }_1)\dots n_{k\uparrow }^{{}^{\prime }}({\tau }_k)n_{k\downarrow }^{{}^{\prime }}({\tau }_k){>}_0,$$

$$G_{\sigma }(\tau )=\frac{<T{c}_{\sigma }^{†}(\tau )c_{\sigma }(0)c_{\sigma }^{†}({\tau }_1)c_{\sigma }({\tau }_1)\dots c_{\sigma }^{†}({\tau }_k)c_{\sigma }({\tau }_k){>}_0}{\det D_{\sigma }}.$$

$$G_{\sigma }(i{\omega }_n)=G_{\sigma }^0(i{\omega }_n)-G_{\sigma }^0(i{\omega }_n)[\sum _{i,j}{M}_{i,j,\sigma }{e}^{i{\omega }_n({\tau }_i-{\tau }_j)}]G_{\sigma }^0(i{\omega }_n),$$

$$G_{\sigma }^0(i{\omega }_n)=\frac{2}{i{\omega }_n+\sqrt{{\omega }_n^2+1}+2\mu }.$$

$$D=\frac{2<k>-U{T}^{-1}}{2U{T}^{-1}},$$

Full text

Introduction.- Gaining an understanding of unconventional behaviors observed in strongly correlated electronic materials is an interesting and challenging research subject in theoretical condensed matter physics and material science communities, because exact diagonalization of a Hamiltonian with an electronic correlation $U$ requires extremely huge numerical work with increasing lattice size Imada1998 . Dynamical mean field theory (DMFT) approximation is a useful tool in exploring such a Hamiltonian, although it ignores the spatial correlations to avoid excessive numerical burden Georges1996 . DMFT captures unconventional behaviors such as a Mott insulator and a non-Fermi-liquid state beyond the mean field approximations, as well as supplies the reliable spectral functions in combination with the first-principles density functional theory method, which can be directly compared with angle-resolved photoemission spectroscopy experimental ones of real correlated materials such as iron-based superconductors and transition metal oxides Kotliar2006 ; Han2018 ; Kim2018 ; Brito2016 ; Medici2014 ; Karolak2015 ; Sim2018 ; Ruff2013 ; Lee2012 .

The most challenging part of DMFT approximation is to solve a quantum impurity problem that describes how the electrons on an impurity site interact with ones in a thermodynamic bath. In this sense, developments of the exact continuous-time quantum Monte Carlo (CTQMC) approaches for the quantum impurity problem have made remarkable progress in DMFT as well as in strongly correlated electronic material communities Rubtsov2005 ; Werner2006 ; Gull2011 . The unconventional natures of many strongly correlated materials could be explained by the density functional theory plus DMFT (DFT+DMFT) tools in combination with CTQMC methods Kotliar2006 ; Han2018 ; Kim2018 ; Brito2016 ; Medici2014 ; Karolak2015 ; Sim2018 ; Ruff2013 ; Lee2012 . On the other hand, the CTQMC tools still retain heavy computational works in the multi-orbital systems with low temperature $T$ and high Matsubara frequencies Rubtsov2005 ; Werner2006 ; Gull2011 . Therefore, accelerating the CTQMC algorithm for those systems would be highly beneficial.

Computationally, the CTQMC methods have two time-consuming parts: (i) One is the inversion of the $k\times k$ matrix given by the diagram expansion order $k$ in the CTQMC update. (ii) Another one is the multiplication of the $k\times k$ matrix and the non-interacting Green’s function $G_{\sigma}^{0}(i\omega_{n})$ in order to measure the impurity Green’s function $G_{\sigma}(i\omega_{n})$, where $\sigma$ and $\omega_{n}$ are the spin index and the Matsubara frequency, respectively. Computational time in the second one increases proportionally with increasing values of $\omega_{n}$ and $k$.

Recently, modern machine learning (ML) approaches, which conjecture the results from massive data, have been extensively employed in various fields including that of theoretical condensed matter physics Carrasquilla2017 ; Arsenault2014 ; Carleo12017 ; Yoon2018 ; Huang2017 ; Liu2017 ; Otsuki2017 ; Amin2018 . Here, we propose the CTQMC method in combination with the ML (CTQMC+ML) technique, which eliminates the time-consuming step of of the matrix. The CTQMC+ML algorithm is as follows. The raw datasets of $G_{\sigma}(i\omega_{n})$ obtained from the data of the imaginary time $\tau_{k}$ with an array of $k$ order are constructed by the CTQMC simulation. $\tau_{k}$ and $G_{\sigma}(i\omega_{n})$ are employed for input data and output data, respectively. The ML kernel is constructed through a one dimension convolutional autoencoder (CAE) and a stochastic gradient descent based optimizer implemented in Tensorflow-gpu tsflow ; LeCun2015 ; Krizhevsky2012 ; Bengio2014 . After finishing the training of ML, $G_{\sigma}(i\omega_{n})$ is directly predicted by the ML process, where the multiplication of the $k\times k$ matrix and the non-interacting Green’s function $G_{\sigma}^{0}(i\omega_{n})$ with the number of $\omega_{n}$ disappears in the calculation. Finally, we find that our ML results show accurate $G_{\sigma}(i\omega_{n})$ values with high Matsubara frequency and double occupancy in a much shorter computational time than the CTQMC approach.

Continuous time quantum Monte Carlo method in combination with Machine learning.- The Hamiltonian $H$ in the DMFT approximation includes $H_{\text{Local}}$, $H_{\text{Bath}}$, and $H_{\text{Hyb}}$ for an impurity part, a thermodynamic bath part, and a coupling part between the impurity and the thermodynamic bath, respectively. The CTQMC approaches have two versions, one where the partition functions are expanded on the interaction or hybridization. Although the expansion quantities differ rather significantly, the simulation procedures of both algorithms are quite similar. Therefore, we will explain the interaction-expansion CTQMC method in combination with ML, and discuss the results of $G_{\sigma}(i\omega_{n})$ and double occupancy.

Here, we discuss the formalism of the interaction-expansion CTQMC tool briefly Rubtsov2005 . The partition function of an impurity problem is written as

[TABLE]

where $S_{0}$ and $S_{U}$ are the non-interacting part including the hybridization and the electronic interaction part given as $S_{U}=U\int_{0}^{T^{-1}}d\tau(n_{\uparrow}(\tau)-\frac{1}{2})(n_{\downarrow}(\tau)-\frac{1}{2})$, respectively. Here, for simplicity of the notation to remove the Fermionic sign problem in the cases of the repulsive interaction $(U>0)$ we introduce $n_{\uparrow}^{{}^{\prime}}n_{\downarrow}^{{}^{\prime}}$ operator given as $n_{\uparrow}^{{}^{\prime}}n_{\downarrow}^{{}^{\prime}}=U(n_{\uparrow}-\frac{1}{2})(n_{\downarrow}-\frac{1}{2})=\frac{U}{2}n_{\uparrow}(n_{\downarrow}-1)+\frac{U}{2}(n_{\uparrow}-1)n_{\downarrow}$. The expansion of the partition function in the powers of the modified interaction part is given as

[TABLE]

where $\sum_{C_{k}}=\sum_{k=0}^{\infty}\int_{0}^{T^{-1}}d\tau_{1}\dots\int_{0}^{\tau_{k-1}}d\tau_{k}$ and $Z_{0}=\int D[d^{\dagger},d]e^{-S_{0}}$ is the single-site case. The CTQMC update is performed on the diagram expansion order $k$ and the continuous time $\tau$ occurred between [math] and $T^{-1}$. $<n_{1\uparrow}^{{}^{\prime}}(\tau_{1})n_{1\downarrow}^{{}^{\prime}}(\tau_{1})\dots n_{k\uparrow}^{{}^{\prime}}(\tau_{k})n_{k\downarrow}^{{}^{\prime}}(\tau_{k})>_{0}$ is a determinant of a matrix obtained from the Wick’s theorem, and the matrix can be separated into spin up matrix $D_{\uparrow}$ and spin down matrix $D_{\downarrow}$. The impurity Green’s function $G_{\sigma}(i\omega_{n})$ can be expressed as

[TABLE]

In the end, via the Fourier transformation and the Wick’s theorem, $G_{\sigma}(i\omega_{n})$ as a function of the Matsubara frequency is computed with

[TABLE]

where the size $k\times k$ matrix $M_{i,j,\sigma}$ is the inversion of $D_{\sigma}/T$. The average matrix size $<k>$ is proportional to $T^{-1}$. Note that most computational times in the CTQMC algorithm are used to calculate the inversion of $k\times k$ matrix $D_{\sigma}/T$ and the multiplication of the matrix in Eq. (4).

Here, we employ the ML technique to substitute for the calculation of Eq. (4). To construct the ML kernel, we adapt the CAE, which is customized as a one dimensional input layer with values of the $k$ size imaginary time $\tau_{k}$ and deep layered structure with three convolution and deconvolution layers. Figure 1 shows the schematic CAE architecture from input data $\tau_{k}$ to output data $G_{\sigma}(i\omega_{n})$ via deep layered structureBengio2014 ; Kingma2014 . During the convolution process, we apply max pooling and keep the probability larger than $0.7$ with a filter of $1\times 5\times m$, where $m$ denotes a feature map $1-16-32-64-32-16-1$ of ongoing progress. After the deconvolution layer, we add one fully connected layer to match the desired number of output nodes. The softmax activation and stochastic gradient descent based optimizer are applied for the one-dimensional CAE training implemented with the TensorFlow-gpu versiontsflow . We prepare $8\times 10^{4}$ different combinations of $\tau_{k}$ as the input dataset and use batch learning with a few thousand parallel procedures. In all cases, the training times of ML process take approximately one hour with NVIDIA GeForce-1080Ti. The results of the CTQMC+ML approach are compared with those of conventional CTQMC measurements of Eq. (4) with $6\times 10^{5}$.

Results.- We consider the single impurity system on the Bethe lattice with the semi-circular density of states. We set the bandwidth to $W=2.0$ and the temperature $T=0.01$ in all calculations. In these cases the non-interacting Green’s function $G_{\sigma}^{0}(i\omega_{n})$ with the chemical potential $\mu$ is given as

[TABLE]

We study the cases of half-filling with $\mu=0$ in Eq. (5). Figure 2 shows the imaginary part of the impurity Green’s function $\text{Im}[G_{\sigma}(i\omega_{n})]$ obtained from exact Eq. (4) and the CTQMC+ML approach with $8\times 10^{4}$ different trained datasets for $U=2.0$ and $3.0$. The trainings of the ML for $U=2.0$ and $3.0$ are independently performed to obtain each $\text{Im}[G_{\sigma}(i\omega_{n})]$. We approximately spend one hour for the training of the ML process. (i) We find that the computational time in the measurement step of $G_{\sigma}(i\omega_{n})$ of the CTQMC+ML method is almost independent, regardless of the numbers of $k$ order and $\omega_{n}$. It means that the computational time of the CTQMC+ML method is much less than one of the conventional CTQMC approaches with proportional increase, as the numbers of $\omega_{n}$ and $k$ increase. (ii) Although we do not show the CTQMC+ML results obtained from $2\times 10^{4}$ different trained datasets in Figure 2, $\text{Im}[G_{\sigma}(i\omega_{n})]$ of both methods for $U=2.0$ are already in good agreement as results show in Figure 2, while they are quite different for $U=3.0$. This means that the numbers of the trained datasets affect the accuracy of the CTQMC+ML results more strongly than other parameters such as the batch size, the training time, and the filter length.

Next, we perform the training of the ML with input and output datasets of both $U=2.0$ and $3.0$, and predict the double occupancy $D$ and $\text{Im}[G_{\sigma}(i\omega_{n})]$ between $U=2.0$ and $3.0$. $D$ is given as

[TABLE]

where $<k>$ is an average diagram expansion order. Although we already know $<k>$ in Eq. (6) while performing the CTQMC update, we estimate $<k>$ from ML approach again to compute $D$. The results of the imaginary part of impurity Green’s function $\text{Im}[G_{\sigma}(i\omega_{0})]$ at the first Matsubara frequency $\omega_{0}$ and the $D$ as a function of $U$ are shown in Figure 3. $D$ in all $U$ and $\text{Im}[G_{\sigma}(i\omega_{0})]$ in the weak interaction are in a very good agreement, while $\text{Im}[G_{\sigma}(i\omega_{0})]$ for $U=2.8$ and $3.0$ are barely in agreement within error the bar, owing to strong numerical fluctuations and lack of trained datasets. We believe that the results with strong numerical uncertainty will be improved with increasing the measurements in Eq. (4) and the trained dataset in the CTQMC+ML method.

Finally, we confirm the cases of away half-filling for $U=2.0$, where real and imaginary parts of $G_{\sigma}(i\omega_{n})$ as a function of $\omega_{n}$ are shown in Figure 4 (a) and (b), respectively. We employ $\mu=0.4$ and $0.8$ as dopants. The results obtained from Eq. (4) and CTQMC+ML approach are in a very good agreement, as in the cases of half-filling.

Conclusion.- The first-principles DFT+DMFT approaches have proven to be successful in describing physical properties of strongly correlated materials, while the numerically time-consuming impurity step is presented in the DMFT part of these calculations. Therefore, the acceleration of the exact continuous time quantum Monte Carlo (CTQMC) approach sustaining accuracy as the impurity solver of the DFT+DMFT method is a step forward for the multi-orbital systems with low $T$ and high $\omega_{n}$. Here, we propose the CTQMC method in combination with the ML technique based on a one dimensional CAE and a stochastic gradient descent based optimizer in Tensorflow-gpu. We find that our CTQMC+ML can accurately predict $G_{\sigma}(i\omega_{n})$ and $D$, but also considers the physical properties of high Matsubara frequency in much less computational time than one of the conventional CTQMC approaches.

We only confirm the properties of the two-point correlators like $G_{\sigma}(i\omega_{n})$ in the CTQMC+ML approach. One future work is immediately raising. The dynamical vortex and dual fermion approximations, which can capture the nonlocal correlations beyond DMFT approximation, require the four-point correlators in the CTQMC part, where the computational times are proportional to square of the numbers of $k$ and $\omega_{n}$ Toschi2007 ; Rubtsov2008 ; Rohringer2018 . We think that our CTQMC+ML idea will reduce the computational time of the four-point correlators more dramatically than one of the two-point correlators.

I Acknowledgements

This work was supported by the Kangwon National University and Ministry of Science through NRF-2018R1D1A1B07048139 (Hunpyo Lee), and the Ministry of Science, ICT and Future Planning through NRF-2017R1D1A1B03034600 (Taegeun Song).