Neural Network Solver for Small Quantum Clusters

Abstract

Machine learning approaches have recently been applied to the study of various problems in physics. Most of these studies are focused on interpreting the data generated by conventional numerical methods or the data on an existing experimental database. An interesting question is whether it is possible to use a machine learning approach, in particular a neural network, for solving the many-body problem. In this paper, we present a neural network solver for the single impurity Anderson model, the paradigm of an interacting quantum problem in small clusters. We demonstrate that the neural-network-based solver provides quantitative accurate results for the spectral function as compared to the exact diagonalization method. This opens the possibility of utilizing the neural network approach as an impurity solver for other many-body numerical approaches, such as the dynamical mean field theory.

1. Introduction

A single quantum impurity is the simplest quantum many-body problem for which interaction plays a crucial role [1,2]. It was designed as a model to describe diluted magnetic impurities in a non-magnetic metallic host. In the 1960s, it was showed that the perturbation series in coupling strength diverges even with an infinitesimal anti-ferromagnetic coupling value [1]. This early unexpected result motivated the development of innovative approaches to model the strongly correlated systems [1,2,3].

While the physics of a single impurity problem has been rather well studied, interest in the quantum impurity problem was revived during the 1990s. This was partly due to the interest in mapping lattice models onto impurity models [4,5,6,7,8], since at infinite dimensions, the lattice models are equivalent to single impurity models in a mean-field represented by the density of states of the host. This approximated mapping is known as the dynamical mean-field theory. It has been further generalized to cluster impurity models which include some of the effects present in finite dimensional systems [9,10,11].

These mappings provide a systematic tractable approximation for the lattice models and have become a major approach in the field of strongly correlated systems [11]. Combined with density functional theory, they provide one of the best available methods for the study of the properties of materials in which strong interactions are important [12].

The mean-field where the single impurity is embedded, i.e., the density of the bath, can be rather complicated. Therefore, there is in general no analytic method for an accurate solution. Many different methods for solving the effective impurity problem have been proposed. They can be broadly divided in two categories: semi-analytic and numeric.

Between the semi-analytic methods, the most widely used one is the iterative perturbation theory. It interpolates the self-energy at both the weak and strong coupling limits and incorporates some exact constraints, such as the Luttinger theorem [13]. Another widely used semi-analytic method is the local moment approximation. It considers the perturbation on top of the strong coupling limit represented by the unrestricted Hatree Fock solution [14].

Numerical methods can also be divided into two main classes, diagonalization-based and Monte-Carlo-based. Diagonalization methods usually require discretizing the host by a finite number of so-called bath sites. The Hamiltonian which includes the bath sites and the impurity site are diagonalized exactly [15]. Another digonalization-based method is the numerical renormalization group in which the bath sites are mapped onto a one-dimensional chain of sites. The hopping amplitude decreases rapidly down the chain. The model is then diagonalized iteratively as more sites are included [16]. Density matrix renormalization group and coupled cluster theory have also been used as impurity solvers [17,18,19,20].

Hirsh and Fye were the first to propose the quantum Monte Carlo method for solving impurity problems [21]. The time axis of the simulations is divided up using the Trotter–Suzuki approximation. The interaction in each time segment is handled by the Hubbard–Stratonovich approximation [22]. The Monte Carlo method is then used to sample the Hubbard–Stratonovich fields. Continuous time quantum Monte Carlo methods have seen a lot of development over the last decade. They directly sample the partition function without using the Trotter–Suzuki approximation [23], similar to the Stochastic Series Expansion in the simulation of quantum spin models. Notably, the continuous time quantum Monte Carlo using the expansion with respect to the strong coupling limit has been proposed and complicated coupling functions beyond simple Hubbard local density-density coupling terms can now be studied [24].

On the other hand, the past few years have seen tremendous development of machine learning (ML), both algorithms and their implementation [25,26,27]. Many of the ML approaches in physics are designed to detect phase transitions or accelerate Monte Carlo simulations. It is a tantalizing proposal to utilize ML approaches to build a solver for quantum systems.

To build a quantum solver based on the ML approach, we need to identify the feature vector (input data) and the label (output data) for the problem. Then a large pool of data must be generated to train the model, specifically a neural network model. The Anderson impurity problem is a good test bed for the validity of such a solver. We note that similar ideas have been explored using machine learning approaches [28]. This paper is focused on using the kernel polynomial expansion and supervised ML, specifically a neural network, as the building blocks for a quantum impurity solver.

While it is relatively inexpensive to solve a single impurity problem using the above methods in modern computing facilities, current interest in the effects of disorder warrants the new requirement of solving a large set of single or few impurities problems in order to calculate the random disorder average [29,30,31,32,33,34,35]. Our hope is that a fast neural-network-based numerical solver in real frequency can expand the range of applicability of the recently developed typical medium theory to interacting strongly correlated systems, such as the Anderson–Hubbard model [36,37,38,39].

This paper is organized as follows. In the next section, we map the continuous Green function into a finite cluster as has been performed in many dynamical mean field theory calculations. In Section 3, we discuss the expansion of the spectral function in terms of Chebyshev polynomials. In Section 4, we discuss how to use the results from Section 2 and Section 3 as the feature vectors and labels of the neural network. In Section 5, we present the spectral function calculated by the neural network approach. We conclude and discuss future work in the last section.

2. Representing the Host by a Finite Number of Bath Sites

We first identify the input and the output data of a single impurity Anderson model. The input data includes the bare density of states, the chemical potential and the Hubbard interaction on the impurity site. For a system in the thermodynamic limit, the density of states is represented by a continuous function. Since inputting a continuous function to the neural network presents a problem, we describe the continuous bath by a finite number of poles as it is performed within the exact diagonalization approach [15,40,41,42,43]. We approximate the host Green function by a cluster of bath sites,

$$\begin{array}{c}\hfill G_0\left(i{\omega }_n\right)\approx G_0^{cl}\left(i{\omega }_n\right).\end{array}$$

(1)

In the exact diagonalization method, the continuum bath is discretized and represented by a finite number of so-called bath sites, see Figure 1. Assuming that there are $N_b$ bath sites, each bath site is characterized by a local energy (${ϵ}_i$) and a hopping ($t_i$) term to the impurity site. Two additional variables, one the local Hubbard interaction (U) and the other the chemical potential (${ϵ}_f$), are required to describe the impurity site. Therefore, there is a total of $2+2N_b$ variables representing the impurity problem.

The full Hamiltonian in the discretized form (Figure 1) is given as

$$\begin{array}{c}H={\displaystyle \sum _{i,\sigma }}{t}_i(c_{i,\sigma }^{†}{c}_{0,\sigma }+H.c.)+{\displaystyle \sum _{i,\sigma }}{ϵ}_i{c}_{i,\sigma }^{†}{c}_{i,\sigma }+\\ U(c_{0,\uparrow }^{†}{c}_{0,\uparrow }-1/2)(c_{0,\downarrow }^{†}{c}_{0,\downarrow }-1/2)-{ϵ}_f{\displaystyle \sum _{\sigma }}{c}_{0,\sigma }^{†}{c}_{0,\sigma },\end{array}$$

(2)

$c_{i,\sigma }^{†}$ and $c_{i,\sigma }$ are the creation and annihilation operators for site i with spin $\sigma$, respectively. The impurity site is denoted as the 0-th site. The sum of the bath sites goes from 1 to $N_b$, and the spin $\sigma =\pm 1/2$.

The host Green function represented in such a finite cluster can be written exactly as,

$$\begin{array}{c}\hfill G_0^{cl}\left(i{\omega }_n\right)={(i{\omega }_n+{ϵ}_f-\sum _{k=1}^{N_b}\frac{t_k{t}_k^{*}}{i{\omega }_n-{ϵ}_k})}^{-1}.\end{array}$$

(3)

Many different prescriptions for the parameterization of the host Green function have been investigated in detail for exact diagonalization solvers [44]. Conceptually, practical applications of the numerical renormalization group method also require the approximated mapping of the problem onto a finite cluster chain. Unlike the exact diagonalization method, the cluster chain can be rather large; therefore much higher accuracy can be attained in general.

These approaches do not mimic the continuum in the time dimension as it is done by continuous time Quantum Monte Carlo methods, and the mapping onto the finite cluster may represent a nuisance. Nonetheless, this is a necessity for any diagonalization-based method. In our case, the mapping presents an opportunity to naturally adapt the method to a machine learning approach in which a finite discretized set of variables is required.

Under the above approximation, the finite set of variables, $\left\{t_i\right\},\left\{{ϵ}_i\right\},U,{ϵ}_f$, can be treated as the input feature vector for the machine learning algorithm. The next question is what is the desired output or label. We will focus on the spectral function in this study. For this purpose, the next step is to represent the spectral function with a finite number of variables instead of a continuous function. The kernel polynomial method fulfills this goal [45].

3. Expanding the Impurity Green Function by Chebyshev Polynomials

In this section, we briefly discuss the kernel polynomial method for the calculation of the spectral function of a quantum interacting model. Once the host parameters, the impurity interaction and chemical potential are fixed, the ground state of the cluster is obtained using a diagonalization method for sparse matrices. We use the Lanczos approach in the present study [46,47]. Once the ground state is found, the spectral function can be calculated by applying the resolvent operator, $1/(\omega -H)$, to the ground state. A popular method is the continuous fraction expansion [47]. The challenge is that the continuous fraction tends to be under-damped and produce spurious peaks [46,47]. A more recent method is to use an orthogonal polynomial expansion. We will argue that for the application of the ML method, the polynomial expansion method tends to produce better results as we will explain later [45].

The zero temperature single particle retarded Green function corresponding to a generic many-body Hamiltonian is defined as

$$\begin{array}{c}\hfill G\left(\omega \right)=\langle GS\left|c\frac{1}{\omega +i{0}^{+}-H}{c}^{†}\right|GS\rangle ,\end{array}$$

(4)

where $|GS\rangle$ is the ground state of H, and $c^{†}$ and c are the creation and annihilation operators, respectively [45,48]. The spectral function is given as $A\left(\omega \right)=-(1/\pi )Im\left(G\right(\omega \left)\right)$.

Consider the Chebyshev polynomials of the first kind defined as $T_n\left(x\right)=cos\left(narccos\left(x\right)\right)$. Two important properties of the Chebyshev polynomial are their orthogonality and their recurrence relation. The product of two Chebyshev polynomials integrated over $x=[-1,1]$ and weighted by the function $w_n=\frac{1}{\pi \sqrt{1-x^2}}$ is given as

$$\begin{array}{c}\hfill \int dx{w}_n\left(x\right)T_n\left(x\right)T_m\left(x\right)=\frac{1+{\delta }_{n,0}}{2}{\delta }_{n,m}.\end{array}$$

(5)

The recurrence relation is given as

$$\begin{array}{c}\hfill T_n\left(x\right)=2x{T}_{n-1}\left(x\right)-T_{n-2}\left(x\right).\end{array}$$

(6)

The Chebyshev polynomials expansion method is based on the fact that the set of Chebyshev polynomials form an orthonormal basis as defined in Equation (5). Thus a function, $f\left(x\right)$, defined within the range of $x=[-1,1]$ can be expanded as

$$\begin{array}{c}\hfill f\left(x\right)={\mu }_0+2\sum _{n=1}^{\infty }{\mu }_n\left(x\right)T_n\left(x\right),\end{array}$$

(7)

and the expansion coefficient can be obtained by the inner product of the function $f\left(x\right)$ and the Chebyshev polynomials as follow

$$\begin{array}{c}\hfill {\mu }_n={\int }_{-1}^{1}dxf\left(x\right)T_n\left(x\right)w_n\left(x\right).\end{array}$$

(8)

Practical calculations involve truncation at a finite order. The truncation is found to be problematic when the function to expand is not smooth. In our problem, the function is the spectral function of a finite size cluster, which is a linear combination of a set of delta functions. For this reason, a direct application of the above formula will not provide a smooth function. This is analogue to the Gibbs oscillations in the Fourier expansion. The remedy is to introduce a damping factor (kernel) in each coefficient of the expansion [45,49,50,51,52]. Here we use the Jackson kernel given as

$$\begin{array}{c}\hfill f\left(x\right)\approx \sum _{n=0}^N{g}_n{\mu }_n\left(x\right)T_n\left(x\right),\mathrm{where}\end{array}$$

(9)

$$\begin{array}{c}\hfill g_n=\frac{(N-n+1)cos\left(\frac{\pi n}{N+1}\right)+sin\left(\frac{\pi n}{N+1}\right)cot\left(\frac{\pi }{N+1}\right)}{N+1}.\end{array}$$

(10)

We refer the choice of the damping factor to the review by Weiße et al. [45].

We list the steps for calculating the coefficients as follows.

• The input bare Green function is approximated by the bare Green function of a finite size cluster. The set of parameters $\mu ,\left\{t_i\right\},\left\{{ϵ}_i\right\}$ are obtained by minimizing the difference between the left-hand side and the right-hand side of Equation (1) according to some prescriptions [15,40,41,42,43].
• The ground state $\left(\right|GS\rangle )$ and the corresponding energy $\left(E_{GS}\right)$ are obtained by the Lanczos algorithm.
• The spectrum of the Hamiltonian is scaled within the range [−1,1] as required by the Chebyshev expansion. $H\Rightarrow (H-E_{GS})/a$, where a is a real positive constant. The units of energy are also scaled in terms of a.
• The expansion coefficients are given by the inner product between the spectral function and the Chebyshev polynomials.

  $$\begin{array}{c}\hfill {\mu }_n=<{\alpha }_0|{\alpha }_n>,\end{array}$$

  (11)

  where $|{\alpha }_0\rangle =c^{†}|GS>$ and $|{\alpha }_n\rangle =T_n\left(H\right)|{\alpha }_0\rangle$. With the $|{\alpha }_0\rangle$ and the $|{\alpha }_1\rangle =H|{\alpha }_0\rangle$ ready, all the higher order coefficients can be obtained via the recurrence relation

  $$\begin{array}{c}\hfill |{\alpha }_n\rangle =2H|{\alpha }_{n-1}\rangle -|{\alpha }_{n-2}\rangle .\end{array}$$

  (12)
• The spectral function is obtained by feeding the coefficients into Equation (7).

All the coefficients can be obtained by repeated use of Equation (12) which involves matrix vector multiplication. The matrix for an interacting system is usually very sparse, and the computational complexity of the matrix vector multiplication is linear with respect to the vector length, which grows as $4^{N_b+1}$ assuming no reduction by symmetry is applied.

Since practical calculations are limited to a finite order, N, the impurity Green function can be represented by N coefficients of the Chebyshev polynomials expansion. It is worthwhile to mention that the expansion of the Green function in terms of orthogonal polynomials is independent of the method for obtaining the coefficients of the polynomials. Instead of employing exact diagonalization, a recent study shows that one can obtain the expansion coefficients by representing the quantum states by matrix products [48].

4. Feature Vectors and the Labels for the Machine Learning Algorithm

Our strategy is to train a neural network with a large set of variables for the host, i.e., the bath sites, the impurity interaction and the impurity chemical potential. The impurity solver calculates the impurity Green function for a given bath Green function, local impurity interaction and chemical potential that is a total of $2+2N_b$ variables for the input.

The output is the impurity Green function which can be represented by N coefficients of the Chebyshev polynomials expansion. Using the above method, the spectral function is effectively represented in terms of N coefficients. It allows us to naturally employ the supervised learning method by identifying the $2+2N_b$ variables as the input feature vectors and the N variables as the output labels.

While the kernel polynomial method grows exponentially with the number of sites, the end result is represented by a finite number of coefficients which presumably does not scale exponentially with the number of sites. Once the neural network is properly trained, we can use it to predict the impurity Green function without involving a calculation which scales exponentially.

5. Results

We generated 5000 samples for training the neural networks with randomly chosen parameters which are drawn uniformly from the range listed as follows:

$$\begin{array}{c}\hfill t_{i,\uparrow }=t_{i,\downarrow }=[0,1.5],\\ \hfill {ϵ}_{i,\uparrow }={ϵ}_{i,\downarrow }=[-5,5],\\ \hfill U=[0,10],\\ \hfill {ϵ}_f=[-2.5,2.5].\end{array}$$

(13)

We assume that the electron bath has a symmetric density of states. That is, $t_i=t_{i+N_b/2}$ and ${ϵ}_i=-{ϵ}_{i+N_b/2}$ for $i=1$ to $N_b/2$ and $N_b$ even. This further reduces the number of variables in the feature vector to $N_b+2$.

Before embarking on training the neural network, we would like to have some understanding of the range of Chebyshev coefficients. For this purpose, we randomly pick 32 samples and plot the coefficients in Figure 2. There are two prominent features of the coefficients: 1. There are clear oscillations and the coefficients do not decrease monotonically. 2. For all cases shown here, the coefficients essentially vanish for orders around 200 and higher. Due to these two reasons, we decide to train the neural network for the coefficients between 0 and 255.

With the above approximations, the task of solving the Anderson impurity model boils down to mapping a vector containing $N_b+2$ variables to a vector containing N coefficients. For the particular case, we study we choose $N_b=6$ and $N=256$. Machine learning algorithms can thus be naturally applied to this mapping.

We set up an independent dense neural network for each coefficient. The neural network has 14 layers. The input layer contains $N_{in}=N_b+2$ units, and the output layer contains the expansion coefficient for one specific order. The 12 hidden layers have the following number of units: $2N_{in},2N_{in},4N_{in},4N_{in},8N_{in},8N_{in},$$8N_{in},8N_{in},4N_{in},4N_{in},2N_{in},2N_{in}$.

As we consider a total of 256 orders, we train 256 independent neural networks. Considering the coefficients at different orders separately may lose some information contained in the correlations between them. While it is possible to predict a few coefficients by only one neural network, we do not obtain a good prediction using a single neural network for all 256 coefficients without an elaborated fine tuning. Therefore, here instead of searching for a optimal number of coefficients to be predicted for one neural network, we consider each coefficient independently.

In Figure 3, we show the spectral functions corresponding to the same 32 samples used in Figure 2. Both the results from the direct numerical calculation based on the Lanczos method and recurrence relations and those predicted by the neural network are plotted. They basically overlap each other. There is a slight difference for the range of energies where the spectral function is nearly zero. This is perhaps due to the incomplete cancellation among the expansion terms at different orders. An improvement may be attainable if we consider the correlations of the coefficients for different orders. The input parameters of each of the 32 samples are plotted in Figure 4.

Evidence of the capability of the neural network approach can be seen in Figure 5 where we plot the comparison of the first 32 expansion coefficients obtained by the direct numerical calculation and the neural network prediction. We find that two methods give very close results for the 1000 testing samples we consider. Please be reminded that 5000 samples are used for training. All the 1000 testing samples exhibit a linear trend. This clearly shows that a neural network is capable of providing a good prediction. There were no exceptional outliers among the 1000 testing samples we tested.

6. Conclusions

We demonstrate that the supervised machine learning approach, specifically the neural network method, can be utilized as a solver for small quantum interacting clusters. This could be potentially useful for the statistical DMFT or the typical medium theory for which a large number of impurity problems have to be solved to perform the disorder averaging [29,30,31,32,33,36]. The main strategy is to devise a finite number of variables as the feature vector and the label for the supervised machine learning method. In line with the exact diagonalization method for the single impurity Anderson model, the feature vector is represented by the hoppings and the energies of the lattice model. The output we consider, the spectral function, is represented in terms of Chebyshev polynomials with a damping kernel. The labels are then the coefficients of the expansion. By comparing the coefficients directly calculated by the Lanczos method and the recurrence relation with the ones by the neural network, we find the agreement between these two methods is very good. Notably, among the 1000 samples we tested, there is no exceptional outlier. They all have a good agreement with the results of the numerical method.

For a simple impurity problem, the present method may not have an obvious benefit, as a rather large pool of samples have to be generated for training purposes. The situation is completely different for the study of disorder models, such as those being studied by the typical medium theory, where the present method has a clear advantage. Once the neural network is trained, the calculations are computationally inexpensive. For systems in which disorder averaging is required, this method can beat most if not all other numerical methods in term of efficiency. Moreover, the present approach is rather easy to be generalized for more complicated models, such as the few impurities models required in the dynamical cluster approximation. In addition, this method can be easily adapted to the matrix product basis proposed for the kernel polynomial expansion [48,53].

The range of parameters we choose in the present study covers the metal insulator transition of the Hubbard model within the single impurity dynamical mean field theory. An interesting question is the validity of the trained neural networks for the parameters well outside the range of the training data. We expect that the results could deteriorate without additional training data.

The size of the bath we choose in the present study is $N_b=6$. This choice is somewhat arbitrary, for more accurate results for the single impurity problem a larger number of bath sites is preferred. The computational time required to generate the training set scales exponentially with the number of bath sites; however the computational time required to train the neural networks only scales as the third power with the number of bath sites. A larger number of bath sites should pose no problem for the present approach.

The ideas presented in this paper are rather generic. They can be generalized to the solutions by different solvers. For example, this approach can be adapted to Quantum Monte Carlo results as long as they can be represented in some kind of series expansion [54,55]. Our approach can also be adapted to predict the the coefficients of the coupled cluster theory [19,20].