Active Learning Algorithm for Computational Physics

TL;DR

This paper introduces an active learning algorithm using neural networks and query by committee to efficiently compute multi-dimensional physical functions, significantly reducing data requirements while maintaining accuracy.

Contribution

The work presents a novel active learning protocol for fitting multi-dimensional functions in physics, demonstrating efficiency and accuracy improvements over uniform sampling methods.

Findings

Achieves few percent error with less than 10% data compared to uniform sampling.

Effectively identifies regions with rapid function variation for targeted data acquisition.

Applicable to diverse physics problems like atomic and condensed matter physics.

Abstract

In large-scale computation of physics problems, one often encounters the problem of determining a multi-dimensional function, which can be time-consuming when computing each point in this multi-dimensional space is already time-demanding. In the work, we propose that the active learning algorithm can speed up such calculations. The basic idea is to fit a multi-dimensional function by neural networks, and the key point is to make the query of labeled data economically by using a stratagem called "query by committee". We present the general protocol of this fitting scheme, as well as the procedure of how to further compute physical observables with the fitted functions. We show that this method can work well with two examples, which are quantum three-body problem in atomic physics and the anomalous Hall conductivity in condensed matter physics, respectively. In these examples, we show that one reaches an accuracy of few percent error for computing physical observables with less than $10\%$ of total data points compared with uniform sampling. With these two examples, we also visualize that by using the active learning algorithm, the required data are added mostly in the regime where the function varies most rapidly, which explains the mechanism for the efficiency of the algorithm. We expect broad applications of our method on various kind of computational physics problems.

Figures

Equations

$$1:$$

$$2:$$

$$3:$$

$$4:$$

$$5:$$

$$\hat{H}=\sum _{i=1}^3-\frac{{\nabla }_i^2}{2m}+\sum _{i<j,=1}^{3}V(|{\mathbf{r}}_i-{\mathbf{r}}_j|),$$

$$\int \frac{d{k}^{\prime }{k}^{\prime 2}}{2{\pi }^2}\mathcal{U}(k,k^{\prime })D_s\left(k^{\prime },E\right)\mathcal{A}(k^{\prime })-\mathcal{A}(k)=\mathcal{U}(k,0),$$

$$D_s^{-1}\left(k^{\prime },E\right)=\frac{m}{2\pi a_{\text{s}}}-\frac{m}{2\pi }\sqrt{\frac{3k^{\prime 2}}{4}-mE},$$

$$a_{\mathrm{ad}}=-\frac{8}{3m{a}_s}\mathcal{A}(k=0).$$

$$\hat{H}=\sum _{i,j}{t}_{ij}{c}_i^{+}{c}_j,$$

$${\sigma }_{xz}(\mu )$$

$${\Omega }_{xz}^n(k)$$

Full text

I Background

Neural network (NN) based supervised learning methods has nowadays found broad applications in studying quantum physics of condensed matter materials and atomic, molecular and optical systems review1 ; review2 . On the theoretical side, applications include, for example, finding orders and topological invariants in quantum phases Phase1 ; Phase2 ; Phase3 ; Phase4 ; Phase5 ; Phase6 ; Phase7 , generating variational wave functions for quantum many-body states VariationalWF ; WFRBM ; WF ; WFFewBody ; TheoryTomography and speeding up quantum Monte Carlo sampling MC1 ; MC2 . On the experimental side, these methods can help both optimizing experimental protocols State ; Tomography and analyzing experimental data STM ; ExpPhase ; ExpFHM . Usually the supervised learning scheme requires a huge set of labelled data. However, in many physics applications, labelling data can be quite expensive, for instance, performing computation or experiments repeatedly can be time- and resources-demanding. Therefore, in many cases labelled data are not abundant, which is a challenge that have prevented many applications.

The active learning is a scheme to solve this problem ALbook . It starts from training a NN with a small initial dataset, and then actively queries the labelled data based on the prediction of the NN and iteratively improves the performance of the NN until the goal of the task is reached. With this approach, sampling the large parameter space can be made more efficiently, and the request of labelled data is usually much more economical than normal supervised learning methods. Recently a few works have applied the active learning algorithm to determine the inter-atomic potentials in quantum materials PotentialAL1 ; PotentialAL2 ; PotentialAL3 and to optimize control in quantum experiments ExpControl1 ; ExpControl2 , where labelled data have to be obtained either by ab initio calculation or by repeating experiments, which are both time consuming.

In this work we focus on a class of general and common task in computational physics that is to numerically determining a multi-dimensional function, say, $\mathcal{F}(\alpha_{1},\alpha_{2},\dots,\alpha_{n})$, where $\alpha_{i}$ are parameters. Considering a uniform discretization, suppose we discretize each parameter into $L$ points, there are totally $L^{n}$ number of data points that need to be calculated. In many cases, calculation of each point already takes quite some time, and thus the total computational cost is massive. Nevertheless, for most functions, there are regimes where the function varies smoothly and regimes where the function varies rapidly. Ideally, one should sample more points in the steep regimes and less points in the smooth regimes in order to obtain a good fitting in the entire parameter space efficiently. However, it seems to be paradox because one does not know which regimes the function varies more rapidly prior to computing the function.

Here we show that this goal can actually be achieved by using the active learning algorithm and the “query by committee” stratagem Query_by_committee to add data points iteratively. Below we will first introduce the general protocol and then demonstrate the algorithm in two concrete problems. One is the quantum three-body problem and the other is the anomalous Hall conductivity problem. These two are very representative examples in atomic and condensed matter physics, respectively. Using these two examples, we will illustrate how the active learning and “query by committee” guide adding data to steep regimes of the fitting function. We will also discuss how to compute physical observables from the fitted functions.

II General Protocol.

The main procedure is summarized in Fig. 1 and explained as follows:

1. We start with an initial dataset with the number of data $\mathcal{S}_{0}\ll L^{n}$, whose values of $\mathcal{F}$ have been computed exactly. We use this dataset to train a number of different NNs.

2. We ask all NNs to make predictions of $\mathcal{F}$ on all $L^{n}$ number of data points, and for each point we compute the variance among predications made by different NNs.

3. We select $\mathcal{S}_{t}$ number of data point with the largest variance, again with $\mathcal{S}_{t}\ll L^{n}$, and we query the accurate value of $\mathcal{F}$ of these points by numerical calculation.

4. We add the $\mathcal{S}_{t}$ number of new data into the training set to train all NNs again, and then repeat from step 2). We repeat the procedure for $m$-epochs until the results from all NNs converge.

5. We use NN to calculate the value of $\mathcal{F}$ on all $L^{n}$ data points.

In this protocol, one only needs to query the value of $\mathcal{F}$ on $\mathcal{S}_{0}+m\mathcal{S}_{t}$ number of data points, and we should keep $\mathcal{S}_{0}+m\mathcal{S}_{t}\ll L^{n}$. The trade-off is that we need to train a number of different NNs and keep using all NN to make predications on all data points. It can save computational cost if the computational cost for training NN and making predication with NN is much less than the computational cost of $\mathcal{F}$, which is often the case in many applications.

In this protocol, the key idea is “query by committee” Query_by_committee . Here the committee is made of different NNs, which contain different number of layers, different number of nodes at each layer and different activation functions. Thanks to the great expressive power of NN, we do not need to assume a specific form of the fitting functions. In the regime where the function varies smoothly, different NNs can quickly reach a consensus and the variance will be small. On the other hand, in the regime where the function varies rapidly, it is hard for different NNs to converge and the variance will be large. Hence, we can use this variance to guide sampling data point. In fact, as we will see in the examples below, the data points added in later epoch are all added in the regime where the function changes rapidly.

III Three-body Problem

In the first example, we consider a quantum three-boson problem. These bosons interact with a two-body pairwise potential described by an $s$-wave scattering length $a_{\text{s}}$ and a high-energy cut-off $\Lambda$. When $a_{\text{s}}$ is positive, there exists a two-body bound dimer state, and it is important to compute the atom-dimer scattering length $a_{\text{ad}}$ by solving the three-body problem. Here a key quantity is the scattering kernel $\mathcal{U}({\bf k},{\bf k^{\prime}})$. Focusing on the $s$-wave scattering only, $\mathcal{U}$ only depends on the amplitude of momentum $k=|{\bf k}|$ and $k^{\prime}=|{\bf k}^{\prime}|$, and we can simplify $\mathcal{U}$ as $\mathcal{U}(k,k^{\prime})$. Once we know $\mathcal{U}(k,k^{\prime})$, one can calculate $a_{\text{ad}}$ through the Skorniakov-Ter-Martirosion (STM) equation Braaten . The details of how to compute $\mathcal{U}(k,k^{\prime})$ and how to obtain $a_{\text{ad}}$ from $\mathcal{U}(k,k^{\prime})$ are summarized in Appendix A.

In the conventional method, we uniformly discretize both $k$ and $k^{\prime}$ into $L=100$ points between zero and $\Lambda$, and we need to compute $\mathcal{U}$ for all $10^{4}$ data points. We then solve the STM equation with these $U(k,k^{\prime})$ to obtain $a^{\rm exact}_{\text{ad}}$, and we have checked that such a discretization can ensure reaching the convergence, which is referred as the exact result and is shown by the dash line as reference in the Fig. 3. Here we will follow the general procedure described above to fit $\mathcal{U}(k,k^{\prime})$ using NNs. We will show that i) at most only $300$ number of data points are needed in order to obtain a reliable fitting, and ii) most of the queried data occupy the parameter regime where the function $\mathcal{U}(k,k^{\prime})$ is steep, and iii) we use the trained NN to generate $\mathcal{U}(k,k^{\prime})$ on all $10^{4}$ data points to solve STM equation, and we find that the error is only few percent compared with the exact result.

To be more concrete, we start with an initial dataset with uniformly sampled $S_{0}=100$ points. Here we design five different fully connected NNs, whose structures are schematically shown in Fig. 2. The input of all NN are two numbers $k$ and $k^{\prime}$, and the output is $\mathcal{U}$. Each layer of a NN is characterized by the number of nodes $N_{\alpha}$ and an activation function $f_{\alpha}$, and we describe each NN with by $(N_{1},f_{1};N_{2},f_{2};\dots)$. The five different NN used in this work are

[TABLE]

where ${\rm LS}$ denotes the logistic sigmoid activation function with $f(a)=1/(1+e^{-a})$. It is important to keep these NNs different but their detailed structures are not important for final results. For each NN, the training results also depend on the initialization, which is particularly so when the number of data points are not enough. Therefore, for each NN, we also consider $20$ different initialization.

Therefore, at each iteration, we totally have $5\times 20$ trained NN to predict $\mathcal{U}$ on all $\{k,k^{\prime}\}$ points in the set $\mathcal{M}$. For each point $i\equiv\{k,k^{\prime}\}$, we can compute variance $\sigma_{i}$ for all $100$ predications. We will select $\mathcal{S}_{\text{t}}=10$ points with the largest variance to compute $\mathcal{U}$ at these points, and add them into the training set. At each iteration, we can compute the mean variance $\sigma_{\rm mean}=1/L^{2}\sum_{i\subset\mathcal{M}}\sigma_{i}$. When the mean variance $\sigma_{\rm mean}$ is small enough, we start to evaluate the physical quantity $a_{\rm ad}$.

Here it comes to another important consideration of our method. On one hand, we have utilized the discrepancy between NNs to guide the query processes more efficiently, but on the other hand, when we start to compute observables, we need to properly average out these variances between different NNs to obtain a converged result. Hence, we compute $a_{\text{ad}}$ as follows:

1. For each NN, since there are $20$ copies depending on different initializations, we first average the predications over these $20$ copies to obtain a mean value $\bar{\mathcal{U}}_{\rm t}({k,k^{\prime}})$ for each point. We then choose one of $\mathcal{U}_{\rm t}^{j}({k,k^{\prime}})$ that is the closest to $\bar{\mathcal{U}}_{\rm t}({k,k^{\prime}})$, and we use this $\mathcal{U}_{\rm t}({k,k^{\prime}})$ (ignoring the upper index) as the representative of the $t$-th NN.

2. We solve the STM equation with $\mathcal{U}_{\rm t}({k,k^{\prime}})$ and obtain $a^{t}_{\text{ad}}$. Among all five $a^{t}_{\text{ad}}$, we discard the largest and the smallest one, and take an average over the rest three, which is taken as $a_{\text{ad}}$ predicted by the active learning approach. The results are shown in Fig. 3 with blue circles. One can see that the results fluctuate around the exact value, and the fluctuation decreases as $\mathcal{S}$ increases.

3. To further suppress the fluctuation, we can further take a self-average of $a_{\text{ad}}$, that is to average over $a_{\text{ad}}$ for five successive iterations. This result is denoted by $\bar{a}_{\text{ad}}$ and shown in Fig. 3 with yellow solid dots. Indeed, one can see that the fluctuation is already suppressed strongly at $\mathcal{S}\sim 200$.

When the result does not change with training epoch, we take a self-average over the last five iterations to obtain $\bar{a}_{\text{ad}}$, and we compare this result with $a^{\rm exact}_{\rm ad}$. For instance, when $a_{s}\Lambda=10$, the relative error $(\bar{a}_{\rm ad}-a^{\rm exact}_{\rm ad})/a^{\rm exact}_{\rm ad}$ is about $1\%$.

We apply this method to scan different values of $a_{\text{s}}$. We stop the iteration at $\mathcal{S}=300$, a number much smaller than uniformly sampled $10^{4}$ data points, and we take a self-average over the last five iterations to obtain $\bar{a}_{\text{ad}}$. The results are shown in Fig. 4. Our calculation can obtain the right location for the atom-dimer resonances, and our active learning method can well reproduce the Efimov scaling law. By compared with the exact result, we define a relative error as $(\bar{a}_{\rm ad}-a^{\rm exact}_{\rm ad})/a^{\rm exact}_{\rm ad}\equiv\epsilon$‰. The $\epsilon$ around the first resonance is shown in the inset of Fig. 4, which shows that the relative error does increases nearby the Efimov resonance, but overall it is kept within a few thousandths.

To further reveal the mechanism of how our method works, we plot in Fig. 5(a) the function $\mathcal{U}(k,k^{\prime})$ generated by the uniform sampling of all $10^{4}$ points. All the data points added during iteration (with $100<\mathcal{S}\leqslant 300$) are shown in Fig. 3(b). It is very clear that nearly all points are added in the regime around $k\sim 0$ and $k^{\prime}\sim 0$, where $U(k,k^{\prime})$ is the steepest as one can see from Fig. 5(a) and the equal number contour in Fig. 5(b). In Fig. 5(c) and (d), we compare the function $\mathcal{U}(k,k^{\prime})$ generated by NNs and the function generated by uniform sampling. Here the NN results are averaged over five different NNs and different initializations of different NNs. It shows that the fitting is perfect when the calculation of $a_{\text{ad}}$ converges.

IV Anomalous Hall Conductivity Problem

In the second example, we consider the anomalous Hall conductivity of a magnetic Weyl semimetal Mn3Ge Kubler2014 , which has been extensively studied recently Nakatsuji2015 ; Nayak2016 . The anomalous Hall conductivity ($\sigma_{\text{Hall}}$) is an intrinsic quantity induced by the Berry curvature of the band structure. We have computed $\sigma_{\text{H}}=\sum_{\mathbf{k}}\sigma(k_{x},k_{y},k_{z})$ based on the $ab~{}initio$ band structure calculations. The model of this material and the details of computing the Berry curvature are shown in Appendix B. For each $k_{x}$ and $k_{z}$ point, we can obtain a value $\sigma(k_{x},k_{z})$ by performing an integration over $k_{y}$. Here without the active learning method, when we discretize both $k_{x}$ and $k_{z}$ into $100$ points, there are totally $10^{4}$ data points to be computed. The total Hall conductivity is obtained by summing over $k_{x}$ and $k_{z}$ as $\sigma_{\text{H}}=\sum_{k_{x},k_{z}}\sigma(k_{x},k_{z})$.

As one can see from Fig. 6(a), the function $\sigma(k_{x},k_{z})$ is much more singular than $\mathcal{U}(k,k^{\prime})$ in the previous case. There are four broad peaks and four narrow peaks located around $(0.5\pm 0.1,0.5\pm 0.4)$ and $(0.5\pm 0.4,0.5\pm 0.15)$. We follow the same active learning procedure discussed above. In this case we take five different initializations for each NN. Compared with the first example, computing physical observable is easier because we only need to average over all five different initializations for each NN. Similarly, to suppress the fluctuation, at each iteration, we also average over three different predictions of $\sigma_{\text{H}}$ by discarding the largest and the smallest results. We also take a self-average over five successive iterations. As we show in Fig. 7, the prediction of the active learning method approaches the exact value when $\mathcal{S}_{\rm max}>700$, which is less than $10\%$ of the total number of uniformly sampled data, and the relative error is about $0.8$‰. As we shown in Fig. 6(b), the “query by committee” stratagem guides most of the queried date distributed in the areas of four broad peaks. In Fig. 6(c) and (d), we show the function $\sigma(k_{x},k_{z})$ generated by NNs. For $\mathcal{S}=1200$ shown in Fig. 6(c), the fitting has already captured the four main peaks, and since the contribution to the Hall conductance mainly comes from the four main peaks, the results of $\sigma_{\text{H}}$ already converges very well. Nevertheless, when we continue to add labelled data until $\mathcal{S}=4000$, as shown in Fig. 6(d), one can see that these extra data are mainly added in the four narrow peaks such that these four small peaks can also be generated very well.

V Conclusion and Outlook

In summary, we have developed a neural network based machine learning method to determine a multi-dimensional function efficiently. The method combines the great expressivity of NN and the advantage of the active learning scheme to reduce the demand for labeled data to a minimum. There are two key ingredients of this method. On the one hand, we make use of the variances between NNs to guide the queried data to be sampled into the regime where the function varies rapidly, which makes the calculation more efficiently. On the other hand, we need to properly average out these variances when calculating physical observables using the fitted results. With two examples, we show our method can work remarkably well. Compared with uniform sampling, our method can achieve an accuracy of less than $1\%$ error with only less than $10\%$ of total data points, even when the function has multiple sharp peaks. We believe our method can find broad applications in many areas of computational physics.

VI Acknowledge

We thank colleagues in Microsoft Research Asia for discussing active learning. This work is supported Beijing Outstanding Young Scientist Program (HZ), MOST under Grant No. 2016YFA0301600 (HZ) and NSFC Grant No. 11734010 (HZ), the Willner Family Leadership Institute for the Weizmann Institute of Science (BY), the Benoziyo Endowment Fund for the Advancement of Science (BY), and Ruth and Herman Albert Scholars Program for New Scientists (BY).

Appendix A Three-boson problem

The Hamiltonian of the three-boson system is given by

[TABLE]

where ${\bf r}_{i}$ ($i=1,2,3$) is the spatial coordinate of the $i$-th particle, and $V(r)$ is an isotropic short-range potential. Normally we focus on the $s$-wave scattering for ultracold low-energy atoms, and $V(r)$ can be described by the $s$-wave scattering length $a_{\text{s}}$ and a high-momentum cut-off $\Lambda$. Most features of this quantum three-body problem only depend on the dimensionless parameter $a_{\text{s}}\Lambda$ Braaten . We consider the case that $a_{\text{s}}$ is large and positive, where $V(r)$ supports a low-energy bound state called a dimer. A key process to this three-body problem is the atom-dimer scattering illustrated in Fig. 8(a), and the quantity to describe this process is a scattering kernel, which is usually denoted by $\mathcal{U}({\bf k},{\bf k}^{\prime})$, and it can be computed diagrammatically with the Feynman diagram shown in Fig. 8(b) Braaten . When focusing on the $s$-wave scattering only, $\mathcal{U}$ only depends on the amplitude of $k=|{\bf k}|$ and $k^{\prime}=|{\bf k}^{\prime}|$, and we can simplify $\mathcal{U}$ as $\mathcal{U}(k,k^{\prime})$.

Knowing $\mathcal{U}(k,k^{\prime})$, one can calculate the few-body quantity, such as atom-dimer scattering length $a_{\text{ad}}$, through the celebrated Skorniakov-Ter-Martirosion (STM) equation, which can be written as

[TABLE]

where

[TABLE]

is the full propagator for $s$-wave dimer. Solving Eq. 3 with $E$ fixed at the dimer energy $-\frac{1}{ma_{\text{s}}^{2}}$, we can obtain atom-dimer scattering amplitude $\mathcal{A}(k)$, and $a_{\text{ad}}$ is given by Braaten

[TABLE]

In practices, suppose that we discretize both $k$ and $k^{\prime}$ into $L=100$ points between zero and $\Lambda$, Eq. 3 becomes a matrix equation after the discretization, which can be solved by inverting the matrix. Let $\mathcal{M}$ denotes the set of total data points, and the number of data points in $\mathcal{M}$ is $10^{4}$. With uniform sampling, we need to compute $\mathcal{U}$ for all $10^{4}$ points, with which we solve Eq. 3 to obtain $a^{\rm exact}_{\text{ad}}$. This is referred as the exact results in the main text.

Appendix B Anomalous Hall conductivity

The band structure of Mn3Ge calculated with the density-functional theory from Vienna ab-initio simulation package kresse1996 in the framework of the generalized-gradient approximation. The calculated plane wave basis result were projected to atomic-orbital-like Wannier functions Mostofi2008 to get the tight-binding parameters ($t_{ij}$). Based on the tight-binding Hamiltonian,

[TABLE]

we calculated the Berry curvature and the anomalous Hall conductivity in the clean limit. For the sufficient data set, uniform K point grid as $100\times 100\times 100$ employed to produce $10^{6}$ numbers of data points. We have evaluated the AHC ($\sigma_{zx}$) and Berry curvature ($\Omega_{zx}$) by the Kubo-formula approach in the linear response scheme Xiao2010 ,

[TABLE]

Here $n,m$ are the band index, $\epsilon_{n}$ is the eigenvalue of the $|n\rangle$ eigenstate of $H_{k}$ (Fourier transform of $\hat{H}$), and $\hat{v}_{i}=dH_{k}/(\hbar dk_{i})$ ($i=x,z$) is the velocity operator, $\mu$ is Fermi level of the system. The calculated anomalous Hall conductivity value is 296 S/cm where S denoting unit Siemens.