Physical extrapolation of quantum observables by generalization with Gaussian Processes

TL;DR

This paper demonstrates that Gaussian processes can be used to accurately predict quantum properties outside the training data range, enabling physical extrapolation across quantum phase transitions with models of varying complexity.

Contribution

The work introduces a method for physical extrapolation using Gaussian processes, including kernel complexity control with BIC, to predict quantum properties beyond training data.

Findings

GP models can predict quantum phase transitions.

Simple kernels suffice for analytic, simple physical functions.

Complex kernels improve predictions for intricate physical systems.

Abstract

For applications in chemistry and physics, machine learning (ML) is generally used to solve one of three problems: interpolation, classification or clustering. These problems use information about physical systems in a certain range of parameters or variables in order to make predictions at unknown values of these variables within the same range. The present work illustrates the application of ML to prediction of physical properties outside the range of the training parameters. We define `physical extrapolation' to refer to accurate predictions $y(\mathbf{x^\ast})$ of a given physical property at a point $\mathbf{x^\ast} = [ x^\ast_1, ..., x^\ast_{\cal D} ]$ in the $\cal D$-dimensional space, if, at least, one of the variables $x^\ast_i \in \left [ x^\ast_1, ..., x^\ast_{\cal D} \right ]$ is {\it outside} of the range covering the training data. We show that Gaussian processes (GPs) can be used to build ML models capable of physical extrapolation of quantum properties of complex systems across quantum phase transitions. The approach is based on training GP models of variable complexity by the evolution of the physical functions. We show that, as the complexity of the models increases, they become capable of predicting new transitions. We also show that, where the evolution of the physical functions is analytic and relatively simple, GP models with simple kernels already yield accurate generalization results, allowing for accurate predictions of quantum properties in a different quantum phase. For more complex problems, it is necessary to build models with complex kernels. The complexity of the kernels is increased using the Bayesian Information Criterion (BIC). We illustrate the importance of the BIC by comparing the results with random kernels of various complexity and illustrate a method to obtain meaningful extrapolation results without direct validation in the extrapolated region.

Tables

Table 1. Table 1: The numerical values of the BIC ( 14 ) for the models with different simple kernels ( 9 ) - ( 13 ) used for the predictions of the order parameter depicted in Figure 10 .

Kernel type	BIC
RQ	8667.10
RBF	8657.13
MAT	7635.20
LIN	- 104437128213.0
LIN$\times$LIN	-10397873744.9

Equations

$$\mathcal{H}=\sum _k{ϵ}_k{c}_k^{†}{c}_k+\sum _q{\omega }_q{b}_q^{†}{b}_q+V_{\mathrm{e}-\mathrm{ph}},$$

$$V_{\mathrm{e}-\mathrm{ph}}=\alpha H_1+\beta H_2,$$

$$H_1=\sum _{k,q}\frac{2i}{\sqrt{N}}\left[\sin (k+q)-\sin (k)\right]c_{k+q}^{†}{c}_k\left(b_{-q}^{†}+b_q\right)$$

$$H_2=\sum _{k,q}\frac{2i}{\sqrt{N}}\sin (q)c_{k+q}^{†}{c}_k\left(b_{-q}^{†}+b_q\right)$$

$$H=-\frac{J}{2}\sum _{⟨i,j⟩}{\vec{S}}_i\cdot {\vec{S}}_j.$$

$$f(T,m)\approx \frac{1}{2}\left(1-\frac{T_c}{T}\right)m^2+\frac{1}{12}{\left(\frac{T_c}{T}\right)}^3{m}^4,$$

$$\mu ({\mathbf{x}}_{\ast })$$

$$\sigma ({\mathbf{x}}_{\ast })$$

$${\mathbf{x}}_i\Rightarrow \{\mathrm{polaron}\mathrm{momentum}K,\mathrm{Hamiltonian}\mathrm{parameter}\alpha ,$$

$$\mathrm{Hamiltonian}\mathrm{parameter}\beta ,\mathrm{phonon}\mathrm{frequency}\omega \}.$$

$${\mathbf{x}}_i\Rightarrow \{\mathrm{Temperature}T,\mathrm{magnetization}m\}$$

$$k_{\mathrm{LIN}}({\mathbf{x}}_i,{\mathbf{x}}_j)$$

$$k_{\mathrm{RBF}}({\mathbf{x}}_i,{\mathbf{x}}_j)$$

$$k_{\mathrm{MAT}}({\mathbf{x}}_i,{\mathbf{x}}_j)$$

$$k_{\mathrm{RQ}}({\mathbf{x}}_i,{\mathbf{x}}_j)$$

$$\log p(\mathbf{y}|\mathbf{X},\mathbf{\theta })=-\frac{1}{2}{\mathbf{y}}^{\top }{K}^{-1}\mathbf{y}-\frac{1}{2}\log |K|-\frac{n}{2}\log (2\pi ),$$

$$\text{BIC}({\mathcal{M}}_i)=\log p(\mathbf{y}|\mathbf{x},\hat{\theta },{\mathcal{M}}_i)-\frac{1}{2}|{\mathcal{M}}_i|\log n$$

$$c_0{k}_0+c_i{k}_i$$

$$c_i\times k_0\times k_i,$$

$$\mathcal{H}=-\frac{J}{2}\sum _{⟨i,j⟩}{\bar{S}}_i\cdot {\bar{S}}_j,$$

$$f(T,m)\approx \frac{1}{2}\left(1-\frac{T_c}{T}\right)m^2+\frac{1}{12}{\left(\frac{T_c}{T}\right)}^3{m}^4,$$

Full text

I Introduction

As described throughout this book, machine learning has in recent years become a powerful tool for physics research. A large number of machine learning applications in physics can be classified as supervised learning, which aims to build a model ${\cal F}(\cdot)$ of the $\bm{x}\mapsto y$ relation, given a finite number of $\bm{x}_{i}\mapsto y_{i}$ pairs. Here, $\bm{x}$ is a vector of (generally multiple) parameters determining the physical problem of interest and $y$ is a physics result of relevance. For example, $\bm{x}$ could be a vector of coordinates specifying the positions of atoms in a polyatomic molecule and $y$ the potential energy of the molecule calculated by means of a quantum chemistry method fitting-1 ; fitting-2 ; fitting-3 ; fitting-4 ; fitting-5 ; sergei-nn ; JiePRL ; JieJPhysB ; VargasBO ; GPvsNN ; GPper . In this case, ${\cal F}(\bm{x})$ is a model of the potential energy surface constructed based on $n$ energy calculations $\bm{y}=\left(y_{1},...,y_{n}\right)^{\top}$ at $n$ points $\bm{x}_{i}$ in the configuration space of the molecule. To give another example, $\bm{x}$ could represent the parameters entering the Hamiltonian of a complex quantum system (e.g., the tunnelling amplitude, the on-site interaction strength and/or the inter-site interaction strength of an extended Hubbard model) and $y$ some observable such as the energy of the system. Trained by a series of calculations of the observable at different values of the Hamiltonian parameters, ${\cal F}(\bm{x})$ models the dependence of the observable on the Hamiltonian parameters WangPRB2016 ; CarasquillaNatPhys2017 ; NieuwenburgNatPhys2017 ; BroeckerarXiv2017 ; WetzelScherzerarXiv2017 ; Wetzel2017 ; LiuarXiv2017 ; ChangPRX2017 ; BroeckerSciRep2017 ; SchindlerPRB2017 ; OhtsukiJPSJapan2016 ; ArsenaultPRB2014 ; ArsenaultarXiv2015 ; BeacharXiv2017 ; RafaelarXiv2017 ; YoshiokaarXiv2017 ; VenderleyarXiv2017 ; RBM_Troyer ; SchmittarXiv2017 ; CaiPRB ; HuangarXiv ; DengPRB ; NomuraarXiv2017 ; DengPRX ; GaoNatComm ; TorlaiarXiv2017 ; GPdeepNN_2015 ; DanielyNIPS2016 ; GPdeepNN_2017 ; MLcrystal1 ; MLcrystal2 ; MLcrystal3 ; MLcrystal4 , which could be used to map out the phase diagram of the corresponding system.

The ability of a machine learning model to predict previously unseen data is referred to as ‘generalization’. These previously unseen data must usually come from the same distribution as the training data, but may also come from a different distribution. Most of the applications of machine learning in physics aim to make predictions within the range of training data. In the present work, we discuss a method for building machine learning models suitable for physical extrapolation. We define ‘physical extrapolation’ to refer to accurate predictions $y(\bm{x}^{\ast})$ of a given physical property at a point $\bm{x}^{\ast}=\left[x^{\ast}_{1},...,x^{\ast}_{\cal D}\right]$ in the $\cal D$-dimensional input space, if, at least, one of the variables $x^{\ast}_{i}\in\left[x^{\ast}_{1},...,x^{\ast}_{\cal D}\right]$ is outside of the range covering the training data. Thus, in the present work, the training data and test data distributions are necessarily separated in input space. We will refer to the predictions of machine learning models as generalization and the physical problems considered here as extrapolation.

Our particular goal is to extrapolate complex physical behaviour without a priori knowledge of the physical laws governing the evolution of the system. For this purpose, we consider a rather challenging problem: prediction of physical properties of complex quantum systems with multiple phases based on training data entirely in one phase. The main goal of this work is schematically illustrated in Figure 1. We aim to construct the machine learning models that, when trained by the calculations or experimental measurements within one of the Hamiltonian phases (encircled region in Figure 1), are capable of predicting the physical properties of the system in the other phases. Of particular interest is the prediction of the phase transitions, which are often challenging to find with rigorous quantum calculations.

This problem is challenging because the wave functions of the quantum systems – as well as the physical observables characterizing the phases – undergo sharp changes at the phase transitions. Most of the machine learning models used for interpolation/fitting are, however, smooth functions of $\bm{x}$. So, how can one construct a machine learning model that would capture the sharp and/or discontinuous variation of the physical properties? The method discussed here is based on the idea put forward in our earlier work extrapolation-paper .

We assume that the properties of a physical system within a given phase contain information about multiple transitions and that, when a system approaches a phase transition, the properties must change in a way that is affected by the presence of the transition as well as the properties in the other phase(s). In addition, a physical system is often characterized by some properties that vary smoothly through the transition. The goal is then to build a machine learning model that could be trained by such properties within a given phase, make a prediction of these properties in a different phase and predict the properties that change abruptly at the transition from the extrapolated properties. We will use Gaussian processes to build such models.

I.1 Organization of this chapter

The remainder of this chapter is organized as follows. The next section describes the quantum problems considered here. Understanding the physics of these problems is not essential for understanding the contents of this chapter. The main purpose of Section II is to introduce the notation for the physical problems discussed here. These problems are used merely as examples. Section III briefly discusses the application of Gaussian process regression for interpolation in multi-dimensional spaces, mainly to set the stage and define the notation for the subsequent discussion. Section IV describes the extension of Gaussian process models to the extrapolation problem. Section V presents the results and Section VI concludes the present chapter. We will abbreviate ‘Gaussian process’ as GP, ‘Artificial Neural Network’ as NN and ‘machine learning’ as ML throughout this chapter.

II Quantum systems

In this section, we describe the quantum systems considered in the present work. In general, we consider a system described by the Hamiltonian $\hat{H}=\hat{H}(\Gamma)$ that depends on a finite number of free parameters $\Gamma=\{\alpha,\beta,\dots\}$. The observables depend on these Hamiltonian parameters as well as the intrinsic variables ${V}=\{v_{1},v_{2},\dots\}$ such as the total linear momentum for few-body systems or thermodynamic variables for systems with a large number of particles. The set $\Gamma+{V}$ comprises the independent variables of the problems considered here. The ML models $\cal F$ will be functions of $\Gamma+{V}$.

More specifically, we will illustrate the extrapolation method using two completely different quantum models: the lattice polaron model and the mean-field Heisenberg model.

II.1 Lattice polarons

The lattice polaron model describes low-energy excitations of a quantum particle hopping on a lattice coupled to the bosonic field provided by lattice phonons. We consider a quantum particle (often referred to as the ‘bare’ particle) in a one-dimensional lattice with $N\rightarrow\infty$ sites coupled to a phonon field:

[TABLE]

where $c_{k}$ and $b_{q}$ are the annihilation operators for the bare particle with momentum $k$ and phonons with momentum $q$, $\epsilon_{k}=2t\cos(k)$ is the energy of the bare particle and $\omega_{q}=\omega={\rm const}$ is the phonon frequency. The particle-phonon coupling is chosen to represent a combination of two qualitatively different polaron models:

[TABLE]

where

[TABLE]

describes the Su-Schrieffer-Heeger (SSH) ssh particle-phonon coupling, and

[TABLE]

is the breathing-mode model breathing-mode . We will focus on two specific properties of the polaron in the ground state: the polaron momentum and the polaron effective mass. The ground state band of the model (1) represents polarons whose effective mass and ground-state momentum are known to exhibit two sharp transitions as the ratio $\alpha/\beta$ increases from zero to large values HerreraPRL . At $\alpha=0$, the model (1) describes breathing-mode polarons, which have no sharp transitions no-transition . At $\beta=0$, the model (1) describes SSH polarons, whose effective mass and ground-state momentum exhibit one sharp transition in the polaron phase diagram ssh . At these transitions, the ground state momentum and the effective mass of the polaron change abruptly.

II.2 The Heisenberg model

The second model we consider here is the Heisenberg model

[TABLE]

This model describes a lattice of interacting quantum spins $S_{i}$, which – depending on the strength of the interaction $J$ – can be either aligned in the same direction (ferromagnetic phase) or oriented randomly leading to zero net magnetization (paramagnetic phase). The parameter $J$ is the amplitude of the interaction and the $\langle..\rangle$ brackets indicate that the interaction is non-zero only between nearest neighbour spins.

Within a mean-field description, this many-body quantum system has free energy density Chaikin ; sachdev

[TABLE]

where $m$ is the magnetization, $T$ is the temperature and $T_{c}$ is the critical temperature of the phase transition. At temperatures ${T}>T_{c}$, the model yields the paramagnetic phase, while ${T}<T_{c}$ corresponds to the ferromagnetic phase. The main property of interest here will be the order parameter. This property undergoes a sharp change at the critical temperature ovf the paramagnetic - ferromagnetic phase transitions.

III Gaussian process regression for interpolation

The purpose of GP regression is to make a prediction of some quantity $y$ at an arbitrary point $\bm{x}\in\left[\bm{x}_{\rm min},\bm{x}_{\rm max}\right]$ of a $\cal D$-dimensional space, given a finite number of values $\bm{y}=\left(y_{1},...,y_{n}\right)^{\top}$, where $y_{i}$ is the value of $y$ at $\bm{x}_{i}$. Here, $\bm{x}_{i}$ is a $\cal D$-dimensional vector specifying a particular position in the input space and it is assumed that the values $\bm{x}_{i}$ sample the entire range $\left[\bm{x}_{\rm min},\bm{x}_{\rm max}\right]$. If the training data are noiseless (as often will be the case for data coming from the solutions of physical equations), it is assumed that $y$ is represented by a continuous function $f$ that passes through the points $y_{i}$, so the vector of given results is $\bm{y}=\left(f(\bm{x}_{1}),...,f(\bm{x}_{n})\right)^{\top}$. The goal is thus to infer the function $f(\bm{x})$ that interpolates the points $y_{i}\equiv f(x_{i})$. The values $y_{i}$ in the vector $\bm{y}$ represent the ‘training data’.

GPs infer a distribution over functions $p(f|{\bm{y}})$ given the training data, as illustrated in Figure 2. The left panel of Figure 2 shows an example of the GP prior, i.e. the GP before the training. The right panel shows the GP conditioned by the training data (red dots). The GP is characterized by a mean function ${\bm{\mu}}({\bm{x}})$ and covariance $\Sigma({\bm{x}})$. The matrix elements of the covariance are defined as $\Sigma_{ij}=k({\bm{x}}_{i},{\bm{x}}_{j})$, where $k(\cdot,\cdot)$ is a positively defined kernel function.

It is possible to derive the closed-form equations for the conditional mean and variance of a GP gpbook , yielding

[TABLE]

where $\bm{x}_{\ast}$ is a point in the input space where the prediction $\bm{y}_{\ast}$ is to be made; $K(\bm{x},\bm{x})$ is the $n\times n$ square matrix with the elements $K_{i,j}=k(\bm{x}_{i},\bm{x}_{j})$ representing the covariances between $y(\bm{x}_{i})$ and $y(\bm{x}_{j})$. The elements $k(\bm{x}_{i},\bm{x}_{j})$ are represented by the kernel function. Eq. (7) can then be used to make the prediction of the quantity $y$ at point $\bm{x}^{\ast}$, while Eq. (8) can be used to define the error of the prediction.

In this work, the GP models are trained by the results of quantum mechanical calculations. For the case of the polaron models considered here,

[TABLE]

For the case of the Heisenberg model considered here,

[TABLE]

As already mentioned, $\bm{y}\Rightarrow f(\bm{x})$ is a vector of quantum mechanics results at the values of the parameters specified by $\bm{x}_{i}$. For the case of the polaron models considered here, $\bm{y}\Rightarrow{\rm polaron~{}energy~{}}E$. For the case of the Heisenberg model considered here, $\bm{y}\Rightarrow{\rm free~{}energy~{}density}$.

To train a GP model, it is necessary to assume some analytic form for the kernel function $k(\cdot,\cdot)$. In the present work, we will use the following analytic forms for the kernel functions:

[TABLE]

where $r^{2}({\bm{x}}_{i},{\bm{x}}_{j})=({\bm{x}}_{i}-{\bm{x}}_{j})^{\top}\times{M}\times({\bm{x}}_{i}-{\bm{x}}_{j})$ and ${M}$ is a diagonal matrix with different length-scales $\ell_{d}$ for each dimension of ${\bm{x}}_{i}$. The unknown parameters of these functions are found by maximizing the log marginal likelihood function,

[TABLE]

where $\bm{\theta}$ denotes collectively the parameters of the analytical function for $k(\cdot,\cdot)$ and $|K|$ is the determinant of the matrix $K$. Given the kernel functions thus found, Eq. (7) is a GP model, which can be used to make a prediction by interpolation.

III.1 Model selection criteria

As Eq. (7) clearly shows, the GP models with different kernel functions will generally have a different predictive power. In principle, one could use the marginal likelihood as a metric to compare models with different kernels. However, different kernels have different numbers of free parameters and the second term of Eq. (13) directly depends on the number of parameters in the kernel. This makes the log marginal likelihood undesirable to compare kernels of different complexity.

As shown in Ref. bic , a better metric could be the Bayesian information criterion (BIC) defined as

[TABLE]

where $|{\cal M}_{i}|$ is the number of kernel parameters of the kernel ${\cal M}_{i}$. In this equation, $p({\bm{y}}|{\bm{x}},\hat{\mathbf{\theta}},{\cal M}_{i})$ is the marginal likelihood for the optimized kernel $\hat{\bm{\theta}}$ which maximizes the logarithmic part. The assumption – one that will be tested in the present work for physics applications – is that more physical models have a larger BIC. The last term in Eq. (14) penalizes kernels with a larger number of parameters. The optimal BIC will thus correspond to the kernel yielding the largest value of the log marginal likelihood function with the fewest number of free parameters.

IV Physical extrapolation by generalization with Gaussian Processes

As shown Refs. kernel_comb ; gp-ss , one can use the BIC to increase the generalization power of GP models. The algorithm proposed in Refs. kernel_comb ; gp-ss aims to build up the complexity of kernels, starting from the simple kernels (9) - (13), in a greedy search algorithm guided by the values of the BIC. Here, we employ this algorithm to extrapolate the quantum properties embodied in lattice models across phase transitions.

IV.1 Learning with kernel combinations

The approach adopted here starts with the simple kernels (9) - (13). For each of the kernels, a GP model is constructed and the BIC is calculated. The kernel corresponding to the highest BIC is then selected as the best kernel. We will refer to such kernel as the ‘base’ kernel and denote it by $k_{0}$. The base kernel is then combined with each of the kernels (9) - (13). The new ‘combined’ kernels are chosen to be either of the sum form

[TABLE]

or of the product form

[TABLE]

where $c_{0}$ and $c_{i}$ are treated as independent constants to be found by the maximization of the log marginal likelihood. The GP models with each of the new kernels are constructed and the BIC values are calculated. The kernel of the model with the largest BIC is then chosen as $k_{0}$ and the process is iterated. We thus have an ‘optimal policy’ algorithm RL that selects the kernel assumed optimal based on the BIC at every step in the search.

We note that a similar procedure could be used to improve the accuracy of GP models for the interpolation problems. We have done this in one of our recent articles jun-paper , where GP models were used to construct a six-dimensional potential energy surface for a chemically reactive complex with a very small number of training points. In the case of interpolation problems, it may also be possible to use cross-validation for kernel selection arthur-paper . Cross-validation could also be applied to kernel selection for the extrapolation problems. We have not attempted to do this in the present work. We will compare the relative performance of the validation error and the BIC as the kernel selection metric in a future work.

V Extrapolation of quantum properties

In this section, we present the results illustrating the performance of the algorithm described above for the prediction of the quantum properties of complex systems outside the range of the training data. Our particular focus is on predicting properties that undergo a sharp variation or discontinuity at certain values of the Hamiltonian parameters. Such properties cannot be directly modelled by GPs because the mean of a GP is a smooth, differentiable function.

The main idea extrapolation-paper is to train a GP model with functions (obtained from the solutions to the Schrödinger equation) that vary smoothly across phase transitions and derive the properties undergoing sharp changes from such smoothly varying function. We thus expect this procedure to be generally applicable to extrapolation across second-order phase transitions. Here, we present two examples to illustrate this. The particular focus of the discussion presented below is on how the method converges to the accurate predictions as the complexity of the kernels increases.

V.1 Extrapolation across sharp polaron transitions

As discussed in Section II, the Hamiltonian describing a quantum particle coupled to optical phonons through a combination of two couplings defined by Eq. (2) yields polarons with unusual properties. In particular, it was previously shown HerreraPRL that the ground-state momentum of such polarons undergoes two sharp transitions as the ratio $\alpha/\beta$ in Eq. (2) as well as the parameter $\lambda=2\alpha^{2}/t\hbar\omega$ are varied. The dimensionless parameter $\lambda$ is defined in terms of the bare particle hopping amplitude $t$ and the phonon frequency $\omega$. It quantifies the strength of coupling between the bare particle and the phonons. One can thus calculate the ground-state momentum or the effective mass of the polaron as a function of $\lambda$ and $\alpha/\beta$. The values of $\lambda$ and $\alpha/\beta$, where the polaron momentum and effective mass undergo sharp changes, separate the ‘phases’ of the Hamiltonian (1).

The GP models are trained by the polaron energy dispersions (i.e. the full curves of the dependence of the polaron energy on the polaron momentum) at different values of $\lambda,\alpha$ and $\beta$. These models are then used to generalize the full energy dispersions to values of $\lambda,\alpha$ and $\beta$ outside the range of the training data and the momentum of the polaron with the lowest energy is calculated from these dispersion curves. The results are shown in Figure 3. Each of the white dots in the phase diagrams depicted specifies the values of $\alpha,\beta$ and $\lambda$, for which the polaron dispersions were calculated and used as the training data. One can thus view the resulting GP models as four-dimensional, i.e. depending on $\alpha$, $\beta$, $\lambda$ and the polaron momentum.

Figure 3 illustrates two remarkable results:

• $\circ$

The left panel illustrates that the GP models are capable of predicting multiple new phase transitions by using the training data entirely in one single phase. This proves our conjecture extrapolation-paper that the evolution of physical properties with the Hamiltonian parameters in a single phase contains information about multiple phases and multiple phase transitions.

• $\circ$

While perhaps less surprising, the right panel illustrates that the accuracy of the predictions increases significantly and the predictions of the phase transitions become quantitative if the models are trained by data in two phases. The model illustrated in this panel extrapolates the polaron properties from high values of $\lambda$ to low values of $\lambda$. Thus, the extrapolation becomes much more accurate if the models are trained by data in multiple phases.

In the following section we analyze how the kernel selection algorithm described in Ref. extrapolation-paper and briefly above arrives at the models used for the predictions in Figure 3.

V.2 Effect of kernel complexity

Figure 4 illustrates the performance of the models with kernels represented by a simple addition of two simple kernels, when trained by the data in two phases, as in the rigth panel of Figure 3. The examination of this figure shows that the generalization accuracy, including the prediction of the number of the phase transitions, is sensitive to the kernel combination. For example, the models with the combination of the RBF and LIN kernels do not predict any phase transitions. Most of the other kernel combinations predict only one of the two transitions. Remarkably, the combination of two RBF kernels already leads to the appearance of the second phase transition, and allows the model to predict the location of the first transition quite accurately. The combination of Figures 3 and 4 thus illustrates that the BIC is a meaningful metric to guide the kernel selection algorithm, as it rules out many of the kernels leading to incorrect phase diagrams shown in Figure 4. The results in Figure 4 also raise the question, how many combinations are required for kernels to allow quantitative predictions?

To answer this question, we show in Figures 5 and 6 the convergence of the phase diagrams to the results in Figure 3 with the number of iterations in the kernel selection algorithm. We use the following notation to label the figure panels: GPL-$X$, where $X$ is the number of iteration. Thus, $X=0$ corresponds to level zero of the kernel selection algorithm, i.e. GPL-0 is the phase diagram predicted by the model with the single simple kernel leading to the highest BIC. Level $X=1$ corresponds to kernels constructed as the simple combinations (15) or (16). Level $X=2$ corresponds to kernels of the form (15) or (16), where $k_{i}$ is a combination of two kernels. As can be seen from Figures 5 and 6, level $X=2$ and $X=3$ produce kernels with sufficient complexity for accurate predictions.

It must be noted that increasing the complexity of the kernels further (by increasing $X$) often decreases the accuracy of the predictions. This is illustrated in Figure 7. We assume that this happens either due to overfitting or because the kernels become so complex that it is difficult to optimize them and the maximization of the log marginal likelihood gets stuck in a local maximum. To overcome this problem, one needs to optimize kernels multiple times starting from different conditions (either different sets of training data or different initial kernel parameters) and stop increasing the complexity of kernels when the optimization produces widely different results. Alternatively, the models could be validated by a part of the training data and the complexity of the kernels must be stopped at level $X$ that corresponds to the minimal validation error, as often done to prevent overfitting with NNs overfitting_NN .

V.3 Extrapolation across paramagnetic - ferromagnetic transition

In this section, we discuss the Heisenberg spin model described by the lattice Hamiltonian

[TABLE]

where $\langle i,j\rangle$ only account for nearest-neighbour interactions between different spins $\bar{S}_{i}$. The free energy of the system can be calculated within the mean-field approximation to yield

[TABLE]

where $m$ is the magnetization and $T_{c}=1.25$ is the critical temperature of the phase transition between the paramagnetic (${T}>T_{c}$) and ferromagnetic (${T}<T_{c}$) phase.

We train GP models by the entire free-energy curves at temperatures far above $T_{c}$. The free energy curves are then predicted by the extrapolation models at temperatures decreasing to the other side of the transition. The oder parameter $m_{0}$ – defined as the value of magnetization that minimizes free energy – is then computed from the extrapolated predictions. The results are shown in Figure 8.

As evident from Eq. (18), the free-energy curves have an analytic dependence on temperature $T$ so this is a particularly interesting case for testing the generalization models. Can the kernel selection algorithm adopted here converge to a model that will describe accurately the analytic dependence of the free energy (18) as well as the order parameter derived from it? We find that the temperature dependence of Eq. (18) can be rather well captured and accurately generalized by a model already with one simple kernel! However, this kernel must be carefully selected. As Figure 9 illustrates, the accuracy of the free-energy prediction varies widely with the kernel. This translates directly into the accuracy of the order-parameter prediction illustrated by Figure 10. Figure 10 illustrates that the RBF and RQ kernels capture the evolution of the order parameter quantitatively, while the LIN, MAT and quadratic kernels produce incorrect results.

Table 1 lists the BIC values for the models used to obtained the results depicted in Figure 10, clearly demonstrating that the higher value of the BIC corresponds to the model with the better prediction power.

V.4 Validation of extrapolation

Validation of the generalization predictions in the extrapolated region presents a major problem. By definition, there are no data in the extrapolated region. One can, of course, divide the given data into a training set and a validation set outside of the training data range. The validation set can then be used to verify the accuracy of the extrapolation. This is what is done throughout this work. However, this does not guarantee the accuracy of the predictions beyond the range of the validation data. Finding a proper method to validate the extrapolation predictions is particularly important for applications of the present approach to making predictions of observables at physical parameters, where no theoretical or experimental results are available.

A possible way to verify the accuracy of the extrapolation predictions without using data in the extrapolated region is to examine the sensitivity of the predictions to the positions and number of training points. If the predictions are stable to variations of the training data, one might argue that the predictions are valid. To illustrate this, we rebuild the models of the phase diagram depicted in Figure 1 with a variable number of training points. Figure 11 shows the results obtained with models trained by the quantum calculations at different values of $\lambda$ and $\alpha/\beta$. The figure illustrates the following:

• $\circ$

The generalization models capture both transitions even when trained by the quantum calculations far removed from the transition line and with a random distribution of training points.

• $\circ$

The predictions of the transitions become more accurate as the distribution of the training points approaches the first transition line.

One may thus conclude that the predictions of the sharp transitions are physical. If possible, this can be further validated by training generalization models with data in a completely different part of the phase diagram. This is done in Figure 12 that shows the same phase diagram obtained by the generalization models trained with quantum results entirely in the middle phase. Given the remarkable agreement of the results in the last panel of Figure 11 and in the last panel of Figure 12, one can argue that the predicted phase diagram is accurate.

Based on these results, we suggest the following algorithm to make stable predictions of unknown phase transitions by extrapolation:

• $(1)$

Sample the phase diagram with a cluster of training points at random.

• $(2)$

Identify the phase transitions by extrapolation in all directions.

• $(3)$

Move the cluster of the training points towards any predicted transition.

• $(4)$

Repeat the calculations until the predictions do not change with the change of the training point distributions.

• $(5)$

If possible, rebuild the models with training points in a completely different part of the phase diagram.

While step (5) is not necessary, the agreement of the results in steps (4) and (5) can be used as an independent verification of the extrapolation. The comparison of the results in steps (4) and (5) may also point to the part of then phase diagram, where the predictions are least reliable.

V.5 Power of the Bayesian Information Criterion

As explained in Section IV.1, the generalization models used here are obtaining by gradually increasing the complexity of kernels using the BIC (14) as a kernel selection criterion. The algorithm starts with a simple kernel that leads to a model with the largest BIC. This kernel is them combined with multiple simple kernels leading to multiple models. The kernel of the model with the largest BIC is selected as a new kernel, which is again combined with multiple simple kernels. The procedure is iterated to increase the kernel complexity, one simple kernel at a time. Since the BIC (14) is closely related to the log marginal likelihood and the kernel parameters are optimized for each step by maximizing the log marginal likelihood, why not to simply select some complex kernel function at random and maximize the log marginal likelihood of the model with this kernel?

To illustrate the power of the BIC in the greedy search algorithm, we repeat the calculations of Figure 12 with kernels of various complexity selected at random. We mimic the iterative process used above, but, instead of using the BIC as a kernel selection criterion, we select a new kernel at each complexity level at random. Every model is optimized by maximizing the log marginal likelihood as usual. The results are depicted in Figure 13. The different panels of Figure 13 are obtained with models using kernels of different complexity. The models are not physical and there is no evidence of model improvement with increasing kernel complexity. We thus conclude that the BIC is essential as the kernel selection criterion to build GP models for applications targeting physical extrapolation.

VI Conclusion

The present article presents clear evidence that Gaussian process models can be designed to predict the physical properties of complex quantum systems outside the range of the training data. As argued in Ref. BML , the generalization power of GP models in the extrapolated region is likely a consequence of the Bayesian approach to machine learning that underlies GP regression. For this reason, the authors believe that Bayesian machine learning has much potential for applications in physics and chemistry BML . As illustrated here, it can be used as a new discovery tool of physical properties, potentially under conditions, where neither theory nor experiment are feasible.

The generalization models discussed here can also be used to guide rigorous theory in search of specific phenomena (such as phase transitions) and/or particular properties of complex systems. Generating the phase diagram, such as the one depicted in Figure 1, presents no computational difficulty (taking essentially minutes of CPU time). One can thus envision the following efficient approach for the generation of the full phase diagrams based on a combination of the GP models with rigorous calculations or experiments:

• $(1)$

Start with a small number of rigorous calculations or experimental measurements.

• $(2)$

Generate the full phase diagram with the GP models with complex kernels. This diagram is likely to be inaccurate at the system parameters far away from the initial training points.

• $(3)$

Use the rigorous calculations or experiments to add training points in the parts of the parameter space, where (a) the system exhibits desired properties of interest; and (b) where the system properties undergo the most rapid change.

• $(4)$

Repeat the calculations until the predictions in the extrapolated region do not change with the change of the training point distributions.

With this approach, one can envision generating complete ${\cal D}$-dimensional phase diagrams with about $10\times{\cal D}$ rigorous calculations or experimental measurements. Training the models and making the predictions in step (2) will generally take a negligibly small fraction of the total computation time.

It should be pointed out that the results presented in this work suggest algorithms to construct complex GP models capable of meaningful predictions in the extrapolated region without direct validation. To do this, one can examine the sensitivity of the predictions to the distribution of the training points for models with the same level of kernel complexity as well as models with different complexity. Increase of the sensitivity to the training points with the kernel complexity would suggest overfitting or insufficient optimization of the kernel parameters. In such cases, the iterative process building up the kernel complexity should be stopped or the process of optimizing the kernel parameters revised. Constructing algorithms for physical extrapolation without the need for validation should be the ultimate goal of the effort aimed at designing ML models for physics and chemistry. Such models could then use all available chemistry and physics information to make meaningful discoveries.

Acknowledgments

We thank Mona Berciu for the quantum results used for training and verifying the ML models for the polaron problem. We thank John Sous and Mona Berciu for the ideas that have led to work published in Ref. extrapolation-paper and for enlightening discussions.